Treecode2: The Power of Pluralism. I. Static Tests

I describe an `oct-tree’ N-body code which randomly shifts, reorients, and resizes the root cell at each time step. Averaging over a plurality of root cell positions and orientations statistically restores translational and rotational invariance. The potentials and forces which result can be much more accurate than those obtained from a single force calculation. In this paper, the principle of averaging is tested on static configurations. The next paper will show how this technique can substantially improve global energy, momentum, and angular momentum conservation at a negligible computational cost.

💡 Research Summary

**
The paper introduces “Treecode2”, a novel variant of the classic oct‑tree N‑body gravity solver that seeks to eliminate the inherent translational and rotational bias of Eulerian tree codes. In traditional Barnes‑Hut (BH86) implementations the root cell is fixed in a particular position, orientation and size. Consequently, the computed gravitational field depends on the absolute placement and rotation of the simulated system, violating Galilean invariance and Newton’s third law. This manifests as non‑zero bulk forces that vary periodically with the system’s offset inside the root cell, especially pronounced for flattened configurations (see Figure 1).

Treecode2’s key idea is to randomize the root cell at every integration step. Three independent random variables are drawn: a rotation matrix (R_k) (uniform over SO(3)), a scale factor (S_k) (log‑uniform within (|\log S_k|\le \log S_{\max})), and a translation vector (T_k) (uniform inside a sphere of radius (T_{\max})). Particle positions (r_p) are transformed to tree coordinates (r’_p = R_k S_k (r_p - T_k)) before tree construction. The root cell is then centered at the origin of this transformed frame and sized to the smallest power‑of‑two cube that encloses all particles, guaranteeing that cell boundaries are exactly representable in binary floating‑point arithmetic.

The tree itself is built exactly as in standard oct‑tree codes: each cell stores total mass, centre‑of‑mass, and a traceless quadrupole tensor (Q_c) computed from the relative positions of its constituent particles. The potential of a cell is expanded to quadrupole order (Eq. 4) and softened using a refined formula (Eq. 6) that includes a correction term based on the scalar quadrupole (\tilde Q_c). This improves accuracy for interactions where the softening length (\epsilon) is comparable to the cell size.

Opening criteria are modified to avoid pathological “detonating satellite” failures. The classic geometric test (s < \ell_c/\theta) (where (s) is the distance from a particle to a cell centre and (\ell_c) the cell edge) is augmented by the offset (d_c) between the cell’s geometric centre and its centre‑of‑mass, yielding (s < \ell_c/\theta + d_c) (Eq. 9). This guarantees that any cell containing the target particle will be opened provided (\theta < 2/\sqrt{3}), a condition comfortably satisfied in typical simulations.

Because each timestep uses an independent random transformation, the systematic bias introduced by any single tree is cancelled by subsequent trees. Over many timesteps the ensemble of trees samples the full space of orientations, scales and translations, and the averaged forces converge to the translationally and rotationally invariant limit. The author demonstrates this with static tests: for a Gaussian sphere and a highly flattened Gaussian disc (both with (N=2^{18}) particles) the bulk force oscillations observed with a fixed root cell disappear when the forces are averaged over (10^4)–(10^5) random trees. Potential errors drop below 0.5 % and acceleration errors below 1 % on average, a substantial improvement over the single‑tree case. The disc, which previously exhibited bulk forces up to 50 times larger than the sphere, shows negligible net force after averaging.

These static results imply that, in dynamical simulations, energy, linear momentum and angular momentum conservation will be markedly better. The averaged forces have decorrelated errors from one step to the next, so the cumulative error growth is suppressed. Importantly, the randomization adds virtually no extra computational cost: generating the random rotation, scale and translation is O(1) per step, and the tree construction and force evaluation retain their (O(N\log N)) scaling.

The paper also contrasts Treecode2 with Dehnen’s (2000) symmetrized oct‑tree, which enforces exact pairwise force symmetry using multipole expansions but requires additional memory and more complex bookkeeping. Treecode2 achieves a comparable (or better) level of momentum conservation through a statistically symmetric approach that is easy to retrofit onto existing BH86‑style codes.

In summary, Treecode2 offers a conceptually simple yet powerful remedy to the long‑standing issue of coordinate‑system bias in oct‑tree gravity solvers. By randomizing the root cell and averaging over many realizations, the method restores Galilean invariance, reduces systematic force errors, and improves conservation properties without sacrificing performance. The author promises a follow‑up paper that will present dynamical tests (energy, momentum, angular momentum drift) and explore optimal choices for the randomization parameters ((\theta), (S_{\max}), (T_{\max})). This work opens a promising avenue for high‑precision, large‑scale N‑body simulations in astrophysics.