🌌 Barnes–Hut N-body

Gravitational softening replaces the 1/r² force law with 1/(r² + ε²), which prevents the force from diverging to infinity when two bodies come very close together. Without softening, close encounters would require an infinitesimally small time step to integrate accurately, making the simulation crash or produce nonsensical "slingshot" velocities. The softening length ε is typically set to a fraction of the mean inter-particle separation; in real cosmological simulations it is around 1 kpc for dark-matter particles.

The quadtree has O(n) nodes and depth O(log n). For each of the n bodies the tree traversal visits O(log n) nodes on average before enough cells satisfy the θ criterion. Multiplying gives O(n log n) force evaluations per step. Building the tree itself also costs O(n log n). In practice the constant factors are small, and for n = 10,000 the speed-up over O(n²) is around 100–1,000-fold.

When two galaxies pass through each other, tidal forces strip stars from the outskirts and fling them into long "tidal tails" and bridges that stretch across millions of light-years. The process is governed almost entirely by gravity; the stellar discs pass through each other largely unimpeded because stars are so widely spaced. The Milky Way and Andromeda galaxies are on a collision course and are expected to merge in about 4.5 billion years, producing a similar tidal tail structure visible in this preset.

The simulation displays the Barnes-Hut evaluation count alongside the naïve n(n–1)/2 comparisons needed for direct summation. For n = 500 and θ = 0.7, Barnes-Hut typically evaluates around 3,000–5,000 force pairs per body per step instead of 499, giving a 50–100× speed-up. As n grows the advantage compounds: at n = 2,000 the naïve method needs nearly 2 million pair evaluations per step while Barnes-Hut needs only around 20,000.

Kinetic energy KE = Σ ½mᵢvᵢ² and gravitational potential energy PE = –Σᵢ≠ⱼ Gmᵢmⱼ/rᵢⱼ (summed over all pairs). In a true conservative simulation total energy E = KE + PE should be constant, but numerical integration introduces errors. The leapfrog (Verlet) integrator used here is symplectic, meaning it conserves a modified Hamiltonian and keeps long-term energy drift very small. Fluctuations you see are partly genuine physical exchanges (e.g. a body in an eccentric orbit swaps KE and PE) and partly numerical truncation error from finite dt.

When you tick "Highlight θ-descent for picked body", clicking a body shows which quadtree cells it approximates vs recurses. The same quadtree cell may be treated differently for bodies at different positions, illustrating that Barnes-Hut is an adaptive algorithm. For chaos in the actual N-body trajectories, use the two-trajectory comparison in related simulations: even a Δ = 10⁻⁵ perturbation to initial conditions leads to exponentially diverging paths, demonstrating that the N-body problem is in general chaotic.

Yes — the GADGET code family (Springel 2001, 2005) uses a hierarchical tree very similar to Barnes-Hut and has been used to simulate the Millennium Simulation (10 billion particles) and the IllustrisTNG galaxy formation runs. Modern codes often combine Barnes-Hut with fast-multipole methods (FMM), particle-mesh (PM) techniques, or GPU parallelism to simulate up to 10¹² particles. The algorithm is also used outside astrophysics, in molecular dynamics, fluid simulation, and rendering of dense particle systems.

Each internal node of the quadtree stores the total mass M = Σmᵢ and the mass-weighted average position r_cm = Σmᵢrᵢ / M of all bodies within that cell. When a distant cell is approximated, the gravitational force on body i is computed as if all that cell's mass were concentrated at r_cm. This single-pole approximation is accurate to first order in s/d; higher-order multipole expansions (quadrupole, octopole) can be included for better accuracy at the cost of more storage and computation.