N-Body Simulation Methods

An N-body simulation is a numerical experiment that follows $N$ self-gravitating masses as they pull on one another under Newton's law of gravity. The particles may be literal — the stars of a globular cluster — or coarse tracers sampling a smooth phase-space distribution, so a billion of them stand in for a galaxy halo's dark matter. Either way the reason the field exists is stark: the gravitational $N$ -body problem has no closed-form solution for $N\geq 3$ , so once you want to know how a star cluster relaxes, how a disk goes bar-unstable, or how the cosmic web grows from a nearly uniform early universe, direct numerical integration is the only route. Three problems then organize everything. The first is arithmetic: gravity couples every particle to every other, so the naive force sum costs $\mathcal{O}(N^2)$ per step — fine for the $10^5$ stars of a cluster, hopeless for the $10^{10}$ particles of a cosmological box; the subject's history is largely the history of cheaper force solvers (Barnes–Hut trees, particle-mesh grids, hybrid TreePM). The second is physical: a real galaxy is collisionless, its relaxation time far longer than its age, yet a coarsely sampled simulation exaggerates graininess, and softening — blurring each point mass over a length $\epsilon$ — is the fix. The third is geometric: gravity is Hamiltonian, so the integrator must preserve phase-space volume, and the symplectic leapfrog does exactly that, which is why energy errors oscillate instead of drifting over millions of steps.

Key takeaways

An N-body simulation integrates Newton's equations for $N$ self-gravitating particles, supplying numerical solutions to the gravitational $N$ -body problem where analytic methods fail for $N\geq 3$ .
Direct summation evaluates all pairwise forces at cost $\mathcal{O}(N^2)$ — accurate and appropriate for collisional star clusters, but ruinous for large $N$ .
Tree codes group distant particles into multipoles for $\mathcal{O}(N\log N)$ cost; particle-mesh solvers put mass on a grid and solve Poisson's equation by FFT at $\mathcal{O}(N_g\log N_g)$ .
Gravitational softening replaces the singular $1/r^2$ force with a smoothed kernel, capping the peak pairwise force at $\tfrac{2}{3\sqrt3}\,Gm/\epsilon^2$ and suppressing spurious two-body scattering.
Symplectic leapfrog integration conserves phase-space volume and holds the energy error bounded over long runs, unlike higher-order Runge–Kutta; method choice follows the physics — collisional systems want direct summation, collisionless systems want tree, PM, or TreePM with softening and large $N$ .
Key idea: every N-body method answers one question — how to evaluate $N$ gravitational forces cheaply and accurately — under two constraints: preserve the collisionless mean-field limit through softening, and preserve Hamiltonian structure through symplectic time-stepping.

The core insight

Every N-body method is an answer to one question — how to evaluate $N$ gravitational forces cheaply and accurately — subject to two constraints: preserve the collisionless (mean-field) limit through softening, and preserve Hamiltonian structure through symplectic time-stepping.

Where this lives today

N-body simulation is the workhorse of computational astrophysics, its frontier set by hardware as much as algorithms. Cosmological codes such as GADGET-4, PKDGRAV3, and AbacusSummit now evolve tens of billions to trillions of particles on GPU supercomputers, producing the mock universes that Euclid, DESI, and Rubin's LSST compare against to constrain dark energy and the neutrino mass through the matter power spectrum and halo mass function. On the collisional side, direct codes descended from Aarseth's NBODY6 — on GPUs as NBODY6++GPU and PeTar — model million-star clusters, tracking mass segregation, core collapse, and the assembly of the black-hole binaries whose mergers LIGO and Virgo detect. A stubborn open problem sits between the two traditions: the "final parsec problem" of how a pair of supermassive black holes at a merged galaxy's center sheds its last orbital energy to coalesce, demanding integration accurate enough to resolve individual stellar scatterings yet large enough to hold a whole nucleus.

The physical picture

Every timestep, each particle must answer "who pulls on me, and how hard?" The three great force solvers are three answers. Direct summation asks every other particle, one at a time — honest, exact, and quadratic in cost. A tree notices that a distant clump of a million particles pulls almost exactly like a single point at its center of mass, so it lumps the clump together whenever it subtends a small enough angle. Particle-mesh goes furthest: it forgets individual particles, smears their mass onto a grid, and solves for the potential everywhere at once with one Fourier transform. The first figure sets the three strategies side by side, with the cost each pays.

Figure

Three ways to answer who pulls on me. Direct summation touches every pair at cost N², a Barnes–Hut tree replaces a well-separated clump by its center of mass whenever the angle s/d is below θ giving N log N, and particle-mesh deposits mass on a grid and solves Poisson by FFT at cost N_g log N_g

Softening is a deliberate blurring, and the second figure makes precise what it buys. A true point mass gives an infinite force at zero separation, so two particles that pass close scatter violently, injecting artificial relaxation a real collisionless galaxy never suffers. Replace $1/r$ by the Plummer-softened $1/\sqrt{r^2+\epsilon^2}$ and the force no longer diverges — it rises, peaks, and turns over, so particles glide past like extended blobs rather than colliding like billiard balls. The peak is capped at the exact value $\tfrac{2}{3\sqrt3}\,Gm/\epsilon^2\approx0.385\,Gm/\epsilon^2$ at $r=\epsilon/\sqrt2$ , and beyond a few $\epsilon$ the softened and Newtonian curves are indistinguishable.

Figure

Softening caps the pairwise force. The Newtonian force Gm/r² diverges as r → 0 (red), while the Plummer-softened force Gm r /(r²+ε²) to the 3/2 power turns over and peaks at r = ε/√2 with F_max = 2/(3√3)·Gm/ε² ≈ 0.385 Gm/ε² (blue), so close encounters no longer scatter violently; beyond a few ε the two curves agree

The leapfrog integrator earns its name from positions and velocities that hop over each other, staggered half a step apart in time — a staggering that makes the scheme time-reversible and symplectic, so the energy error oscillates instead of drifting and an orbit survives millions of steps without spiraling in or out. The deepest tension is statistical: with $N_{\rm sim}\ll N_{\rm real}$ , each simulation particle carries an enormous mass and its graininess is exaggerated, so the craft of collisionless work is to push $N$ high and choose $\epsilon$ so discreteness relaxation stays negligible while the structures of interest remain resolved.

The core equations

The starting point is Newton's second law with the softened force: each particle is accelerated by the vector sum of the smoothed pulls of all the others. Setting $\epsilon\to0$ recovers the exact inverse-square law; a finite $\epsilon$ is what makes the sampled system behave as a collisionless fluid rather than a cloud of scattering point masses.

\ddot{\mathbf{r}}_i = -\,G\sum_{j\neq i}^{N} m_j\,\frac{\mathbf{r}_i-\mathbf{r}_j}{\left(|\mathbf{r}_i-\mathbf{r}_j|^2+\epsilon^2\right)^{3/2}} \qquad (1)

where: $\mathbf{r}_i$ — position of particle $i$ ; $m_j$ — mass of particle $j$ ; $G$ — gravitational constant; $\epsilon$ — softening length; $N$ — number of particles. Reading it: the denominator would be $|\mathbf{r}_i-\mathbf{r}_j|^3$ for true point masses; the added $\epsilon^2$ is the entire content of softening. The choice of force solver is dictated by how this sum is evaluated, and the three canonical strategies have sharply different asymptotic costs per timestep — the reason the field moved from direct summation to trees and grids as $N$ grew.

\text{cost}_{\rm direct}\propto N^2,\qquad \text{cost}_{\rm tree}\propto N\log N,\qquad \text{cost}_{\rm PM}\propto N_g\log N_g \qquad (2)

where: $N$ — number of particles; $N_g$ — number of grid cells. Reading it: direct summation touches all $\tfrac12 N(N-1)$ pairs, the tree lets each particle see only $\mathcal{O}(\log N)$ nodes, and PM cost is set by the FFT on the mesh. A tree decides when to approximate a whole cell by its center of mass with a single geometric test: a cell of width $s$ seen from distance $d$ is "well separated" if the angle it subtends is small enough.

\frac{s}{d} < \theta \qquad (3)

where: $s$ — width of a tree cell; $d$ — distance to the cell's center of mass; $\theta$ — opening-angle tolerance, typically $0.3$ – $0.7$ . Reading it: smaller $\theta$ opens more cells, buying accuracy at higher cost; $\theta\to0$ recovers direct summation. Particle-mesh replaces the pairwise sum by solving the field equation on a grid. Poisson's equation is differential in real space but becomes an algebraic multiplication in Fourier space — the shortcut the FFT exploits.

\nabla^2\Phi = 4\pi G\,\rho \quad\Longrightarrow\quad \tilde{\Phi}(\mathbf{k}) = -\frac{4\pi G}{k^2}\,\tilde{\rho}(\mathbf{k}) \qquad (4)

where: $\Phi$ — gravitational potential; $\rho$ — grid density; $\tilde{\Phi},\tilde{\rho}$ — Fourier transforms; $\mathbf{k}$ — wavevector of magnitude $k$ . Reading it: $\nabla^2$ becomes multiplication by $-k^2$ , so dividing mode-by-mode solves for $\Phi$ ; the $\mathbf{k}=0$ mode is zeroed to subtract the mean density. Time advance uses the kick–drift–kick leapfrog: a half-step velocity kick, a full-step position drift, and a closing half-step kick at the new position. That symmetric structure is what makes it second-order and symplectic.

\mathbf{v}_{n+1/2}=\mathbf{v}_n+\tfrac{1}{2}\mathbf{a}_n\,\Delta t,\ \ \mathbf{r}_{n+1}=\mathbf{r}_n+\mathbf{v}_{n+1/2}\,\Delta t,\ \ \mathbf{v}_{n+1}=\mathbf{v}_{n+1/2}+\tfrac{1}{2}\mathbf{a}_{n+1}\,\Delta t \qquad (5)

where: $\mathbf{r}_n,\mathbf{v}_n$ — position and velocity at step $n$ ; $\mathbf{a}_n$ — acceleration at $\mathbf{r}_n$ ; $\Delta t$ — timestep. Reading it: only one new force evaluation per step is needed (the closing kick is reused as the next opening kick), so leapfrog is as cheap as Euler. The step size is set adaptively from each particle's acceleration, so a particle plunging through a dense region is resolved more finely than one drifting through a void.

\Delta t_i = \eta\,\sqrt{\frac{\epsilon}{|\mathbf{a}_i|}} \qquad (6)

where: $\Delta t_i$ — timestep of particle $i$ ; $\eta$ — accuracy parameter, $\sim0.01$ – $0.05$ ; $\epsilon$ — softening length; $|\mathbf{a}_i|$ — acceleration magnitude. Reading it: $\sqrt{\epsilon/|\mathbf{a}|}$ has units of time and estimates how long the particle needs to cross a fraction of $\epsilon$ . Finally, discreteness sets a floor on how collisionless a simulation can be: finite $N$ induces an artificial two-body relaxation, the numerical echo of the physical relaxation time.

t_{\rm relax} \approx \frac{0.1\,N}{\ln N}\;t_{\rm cross} \qquad (7)

where: $t_{\rm relax}$ — two-body relaxation time; $t_{\rm cross}$ — crossing time; $N$ — number of particles. Reading it: a faithful collisionless run of duration $T$ needs $N$ large enough that $t_{\rm relax}\gg T$ — millions of particles for a Hubble-time integration.

Derivation AThe softened force from the Plummer kernel

The goal is to show how replacing each point mass by a Plummer sphere of scale $\epsilon$ turns the divergent Newtonian force into a smooth, bounded one that still reproduces the far field, and to locate exactly where the softened force peaks — the quantitative fingerprint of softening the second figure displays. We start from the softened potential, take its gradient, find the peak, and read off the near and far limits so the cap on close-encounter scattering is explicit.

Step 1. Replace the point-mass potential $-Gm/r$ by the Plummer-softened form $\Phi(r)=-Gm/\sqrt{r^2+\epsilon^2}$ , which is finite at the origin, equal to $-Gm/\epsilon$ . Why: the singularity that drives unphysical scattering lives in $\Phi$ at $r=0$ ; smoothing it there is the whole point.

Step 2. The acceleration is $\mathbf{a}=-\nabla\Phi=-Gm\,\nabla\big[(r^2+\epsilon^2)^{-1/2}\big]$ ; using $\nabla r=\mathbf{r}/r$ gives $\mathbf{a}(r)=-Gm\,\mathbf{r}/(r^2+\epsilon^2)^{3/2}$ . Why: this is the softened force of equation (1), differing from Newton only by the $\epsilon^2$ in the denominator.

Step 3. Set $\mathrm{d}F/\mathrm{d}r=0$ on the magnitude $F(r)=Gm\,r/(r^2+\epsilon^2)^{3/2}$ ; the derivative's numerator is proportional to $\epsilon^2-2r^2$ , so the force peaks at $r=\epsilon/\sqrt2$ , where $F_{\max}=\tfrac{2}{3\sqrt3}\,Gm/\epsilon^2\approx0.385\,Gm/\epsilon^2$ . Why: this finite ceiling is the mechanism — the maximum kick one particle can deliver is bounded, so no close pass runs away.

Step 4. Check the limits: for $r\gg\epsilon$ the $\epsilon^2$ term is negligible and $F\to Gm/r^2$ (exactly Newtonian); for $r\to0$ the force vanishes linearly, $F\to Gm\,r/\epsilon^3$ , rather than diverging. Why: softening must be invisible on resolved scales and act only below $\epsilon$ , which these limits confirm. *These two limits — recovering the inverse-square law as $\epsilon\to0$ and the peak $F_{\max}=2/(3\sqrt3)\,Gm/\epsilon^2$ at $r=\epsilon/\sqrt2$ — are the quantitative fingerprints of Plummer softening.*

Go deeper · Derivation B · 1 minWhy leapfrog conserves energy and Runge–Kutta does not

Derivation BWhy leapfrog conserves energy and Runge–Kutta does not

The goal is to explain the most important fact about time integration in gravitational dynamics: leapfrog, though only second order, keeps the energy error bounded over arbitrarily many orbits, while a higher-order Runge–Kutta scheme with smaller per-step error drifts secularly and corrupts the orbit. The reason is that leapfrog exactly integrates a nearby shadow Hamiltonian — exposed by splitting the evolution operator and applying the Baker–Campbell–Hausdorff expansion — so geometric structure, not formal order, is what matters.

Step 1. Write $H=T(\mathbf{p})+V(\mathbf{q})$ and the exact one-step evolution as $e^{\Delta t(\hat T+\hat V)}$ , with $\hat T,\hat V$ the Lie operators of the drift and kick. The KDK leapfrog is the symmetric product $e^{\frac{\Delta t}{2}\hat V}\,e^{\Delta t\hat T}\,e^{\frac{\Delta t}{2}\hat V}$ . Why: each factor is an exact partial flow, hence symplectic, and a composition of symplectic maps is symplectic — phase-space volume is preserved exactly.

Step 2. The Baker–Campbell–Hausdorff formula gives $e^{\Delta t\,\hat{\tilde H}}$ with $\tilde H = H + \tfrac{\Delta t^2}{12}\big(\{\{T,V\},T\}+\tfrac12\{\{T,V\},V\}\big)+\mathcal{O}(\Delta t^4)$ , the shadow Hamiltonian. Why: symmetric splitting cancels the odd-order terms, so $\tilde H$ differs from $H$ only at $\mathcal{O}(\Delta t^2)$ — the numerical orbit exactly solves a small, fixed deformation of the real problem.

Step 3. Because the map is symplectic, $\tilde H$ is conserved exactly along the numerical orbit, so the true energy $H=\tilde H-\mathcal{O}(\Delta t^2)$ oscillates within a band of width $\mathcal{O}(\Delta t^2)$ and never drifts. Why: this is the payoff — bounded energy error for all time, so a halo runs for a Hubble time without spurious heating.

Step 4. By contrast RK4 has a smaller local error, $\mathcal{O}(\Delta t^5)$ , but is not symplectic: with no conserved shadow Hamiltonian its energy error accumulates roughly linearly and eventually dominates. Why: for long runs the geometric property beats formal order, which is why production codes integrate with leapfrog, not RK. *The leading energy error scales as $\Delta t^2$ and stays bounded for all time — a result that can be confirmed directly by integrating a Kepler orbit over many periods and watching the energy oscillate within a fixed band rather than drift.*

Go deeper · Derivation C · 1 minThe particle-mesh Poisson solve by FFT

Derivation CThe particle-mesh Poisson solve by FFT

The goal is to trace the full particle-mesh cycle that turns a cloud of particles into the force on each one in only $\mathcal{O}(N_g\log N_g)$ operations, by recasting the global gravity problem as a single convolution the FFT performs cheaply. We deposit mass on a grid, transform, multiply by the Poisson Green's function, transform back, and interpolate to the particles — seeing why the mean-density mode is dropped and why deposit and interpolate must share a kernel.

Step 1. Assign mass to the grid: $\rho(\mathbf{x}_p)=\tfrac{1}{V_{\rm cell}}\sum_i m_i\,W(\mathbf{x}_p-\mathbf{r}_i)$ , spreading each particle over nearby grid points with a kernel $W$ such as Cloud-In-Cell. Why: the FFT acts on gridded fields, so masses must first become a density at the mesh points.

Step 2. Forward transform $\tilde\rho(\mathbf{k})=\mathrm{FFT}[\rho]$ at cost $\mathcal{O}(N_g\log N_g)$ , with periodic boundaries automatic in the discrete Fourier basis. Why: in Fourier space the Laplacian is diagonal — a differential operator becomes plain multiplication.

Step 3. Apply the Green's function mode-by-mode, $\tilde\Phi(\mathbf{k})=-\dfrac{4\pi G}{k^2}\,\tilde\rho(\mathbf{k})$ , and set $\tilde\Phi(\mathbf{0})=0$ . Why: the $\mathbf{k}=0$ mode is the uniform mean density, whose potential diverges as $1/k^2$ ; zeroing it is the Jeans swindle for a periodic box, leaving only the fluctuation potential that drives the dynamics.

Step 4. Inverse transform for $\Phi$ , take $\mathbf{g}=-\nabla\Phi$ (finite differences or $i\mathbf{k}$ in Fourier space), and interpolate back, $\mathbf{g}_i=\sum_p \mathbf{g}(\mathbf{x}_p)\,W(\mathbf{x}_p-\mathbf{r}_i)$ , with the same kernel $W$ as Step 1. Why: identical deposit and interpolate kernels make the grid forces antisymmetric, so total momentum is conserved exactly. *Momentum conservation follows from the shared deposit-and-interpolate kernel; the Green's function $-4\pi G/k^2$ follows directly from $\nabla^2\to-k^2$ in Fourier space; and the isolated-mass potential can be cross-checked against Ewald summation.*

Limiting cases

The right method, and the regime in which its output can be trusted, follows from the system's type and scale.

Case	Result	Meaning
Collisional (small $N$ )	$t_{\rm relax}\lesssim t_{\rm age}$ , $N\sim10^4$ – $10^6$	Star clusters: two-body relaxation is real physics. Use direct summation with adaptive timesteps and regularized close binaries; softening small or zero.
Collisionless (large $N$ )	$t_{\rm relax}\gg t_{\rm age}$ , $N\gtrsim10^7$	Galaxies and dark-matter halos sample the collisionless Boltzmann equation. Use tree, PM, or TreePM with finite softening to suppress artificial relaxation.
Cosmological volumes	$L\sim10^2$ – $10^3\ \mathrm{Mpc}$ , $N\sim10^9$ – $10^{12}$	Comoving coordinates with $a(t)$ in the equations of motion, periodic boundaries via FFT or Ewald summation; TreePM dominates.
High dynamic range	$\rho_{\max}/\rho_{\min}\gg10^3$	Fixed-grid PM under-resolves dense clumps; a tree, TreePM, or adaptive-mesh solve restores small-scale force resolution.

Common trapMore particles and smaller softening always mean a more accurate simulation.

Not so — this is the trap that ruins under-resolved collisionless runs. Larger $N$ does help, but $\epsilon$ and $N$ jointly set the resolution and must be chosen together, not by force accuracy alone. Take $\epsilon$ too small and the pairwise force is barely capped, so a system meant to be a smooth collisionless galaxy suffers artificial scattering: it relaxes, its core heats and expands, particles evaporate, spurious binaries form — none of which the real galaxy does. Take $\epsilon$ too large and you erase genuine small-scale structure, flattening cusps and dissolving subhalos. And too few particles shorten the numerical relaxation time $t_{\rm relax}\approx0.1N/\ln N\;t_{\rm cross}$ below the run duration, so discreteness noise, not physics, drives the evolution. Accuracy is a joint property of $(N,\epsilon,\Delta t)$ and the integrator, verified by convergence tests, not assumed from any one knob.

Worked examples

1 · The cost of a billion-particle step (fully worked)

A flagship run evolves $N=10^{10}$ dark-matter particles. Compare the per-step cost of direct summation, a Barnes–Hut tree, and a PM solver on a $2048^3$ grid.

Step 1. Direct summation evaluates $\tfrac12 N(N-1)\approx\tfrac12 N^2=5\times10^{19}$ pairwise terms per step. Even at $10^9$ interactions/s/core that is $5\times10^{10}$ core-seconds, roughly 1600 core-years for one timestep — and a run needs thousands of steps. Impossible at this $N$ .

Step 2. A tree scales as $N\log_2 N$ . With $\log_2(10^{10})=33.2$ the work is $3.3\times10^{11}$ node interactions per step, so the speedup over direct is $N/\log_2 N=10^{10}/33.2=\mathbf{3.0\times10^{8}}$ — about 300 million times cheaper per step, turning centuries into seconds.

Step 3. PM cost is the FFT on $N_g=2048^3=8.6\times10^9$ cells: $N_g\log_2 N_g=2.8\times10^{11}$ operations — comparable to the tree but with a smaller prefactor, though forces below one cell are lost, which is why production cosmology uses TreePM (PM long-range, tree short-range). Trees and PM cut $\sim10^{20}$ operations to $\sim10^{11}$ , a factor $\sim10^{9}$ — the entire reason billion-particle cosmology is feasible.

2 · How many particles keep a halo collisionless (fully worked)

A halo has crossing time $t_{\rm cross}\approx50\ \mathrm{Myr}$ and must stay faithful over a Hubble time, $t_{\rm age}\approx14\ \mathrm{Gyr}\approx280\,t_{\rm cross}$ . How large must $N$ be so that numerical relaxation is negligible, $t_{\rm relax}\gtrsim10\,t_{\rm age}$ ?

Step 1. From equation (7), $\dfrac{0.1N}{\ln N}\,t_{\rm cross}\geq 10\,t_{\rm age}=2800\,t_{\rm cross}$ ; the crossing times cancel, leaving $\dfrac{N}{\ln N}\geq2.8\times10^{4}$ .

Step 2. Solve by trial: $N=3\times10^5$ gives $N/\ln N=2.4\times10^4$ (too small); $N=4\times10^5$ gives $3.1\times10^4$ (enough). So $\mathbf{N\gtrsim4\times10^{5}}$ just to keep relaxation below the physics.

Step 3. That is only the discreteness floor. Resolving the halo's cusp, subhalos, and streams demands far more, $N\sim10^6$ – $10^9$ . Relaxation sets a hard minimum of $\sim4\times10^5$ , but structural resolution is what pushes real halo runs to $10^6$ – $10^9$ .

3 · Softened force ceiling and the timestep (fully worked)

A Milky-Way-scale run uses $m=10^{6}\,M_\odot$ and $\epsilon=100\ \mathrm{pc}$ . Find the maximum pairwise acceleration and estimate the timestep from $\Delta t=\eta\sqrt{\epsilon/|\mathbf{a}|}$ with $\eta=0.03$ .

Step 1. From Derivation A the acceleration peaks at $a_{\max}=\tfrac{2}{3\sqrt3}\,Gm/\epsilon^2$ . With $G=4.30\times10^{-3}\ \mathrm{pc}\,(\mathrm{km/s})^2/M_\odot$ , the scale $Gm/\epsilon^2=4.30\times10^{-3}\times10^6/100^2=0.430\ (\mathrm{km/s})^2/\mathrm{pc}$ , so $a_{\max}=0.385\times0.430=\mathbf{0.166\ (km/s)^2/pc}$ .

Step 2. Then $\sqrt{\epsilon/a_{\max}}=\sqrt{100/0.166}=24.6\ \mathrm{pc}/(\mathrm{km/s})$ ; with $1\ \mathrm{pc}/(\mathrm{km/s})=0.978\ \mathrm{Myr}$ and $\eta=0.03$ the step is $\Delta t=\mathbf{0.72\ Myr}$ — about 1.4% of the $\sim50\ \mathrm{Myr}$ crossing time, roughly seventy steps per crossing: fine enough to keep the leapfrog energy error small, coarse enough that a Hubble-time run stays affordable.

Physics of the portrait

N-body simulation has no single founding equation and no single founder, but it has an archetype: the patient builder of direct-summation codes who insisted that gravity be integrated honestly, star by star, for the collisional systems where every close encounter matters.

Origins · 1 minSverre Aarseth and how this was discovered

Sverre Aarseth

1934–2024 · direct N-body simulation of star clusters

Sverre Aarseth, born in Norway and for six decades based at the Institute of Astronomy in Cambridge, was the field's great constant. Where cosmologists chased ever-cheaper approximate forces, Aarseth devoted his career to doing the gravitational $N$ -body problem exactly — every pairwise force summed directly — because the collisional systems he cared about, open and globular clusters, live and die by the close encounters that approximate methods smear away. Beginning around 1963 he built and freely gave away the NBODY series, inventing much of the machinery that still makes direct integration tractable: individual and hierarchical timesteps so a tight binary does not throttle the whole cluster; the Ahmad–Cohen neighbour scheme splitting fast near forces from slow far forces; and Kustaanheimo–Stiefel regularization, which removes a close binary's $1/r$ singularity by a coordinate change rather than by softening. He rode every wave of hardware — mainframes, the GRAPE boards, GPUs in NBODY6++GPU — to push $N$ higher without abandoning exactness. His 2003 monograph and freely shared codes trained several generations, making him the embodiment of the collisional, direct-summation tradition that stands opposite the tree and mesh approximations.

History

1941

— Holmberg computes galaxy encounters with 74 movable lightbulbs, reading gravity off photocell currents

1960

— von Hoerner performs the first digital N-body integrations, following up to 16 bodies on an electronic computer

1963

— Aarseth launches the direct-summation NBODY code series that will run for six decades

1986

— Barnes and Hut introduce the hierarchical O(N log N) tree algorithm

2005

— the Millennium Run follows about ten billion particles to map the growth of the cosmic web

Going further

Texts. Binney and Tremaine, Galactic Dynamics, 2nd ed., Ch. 2 and Appendix C (equations of motion, softening, integrators) · Aarseth, Gravitational N-Body Simulations (2003) · Carroll and Ostlie, An Introduction to Modern Astrophysics, 2nd ed. · Dodelson, Modern Cosmology, 2nd ed. and Ryden, Introduction to Cosmology, 2nd ed. (collisionless Boltzmann context) · Hockney and Eastwood, Computer Simulation Using Particles (1988) · Barnes and Hut, Nature 324, 446 (1986). Key numerical benchmarks: the softened-force peak is $F_{\max}=2/(3\sqrt3)\,Gm/\epsilon^2$ at $r=\epsilon/\sqrt2$ , recovering the inverse-square law in the $\epsilon\to0$ limit; the PM Green's function $-4\pi G/k^2$ follows from $\nabla^2\to-k^2$ ; the leapfrog energy error is bounded via the BCH shadow-Hamiltonian argument; the tree-vs-direct cost ratio is $N/\log_2 N\approx3\times10^8$ at $N=10^{10}$ ; and the collisionless relaxation floor is $N\gtrsim4\times10^5$ .

Builds on: Gravitational N-Body Problem · Gravitational Two-Body Problem · Collisionless Boltzmann Equation · Two-Body Relaxation Time · Distribution Function and Phase Space. Leads to: Virial Theorem in Stellar Systems · Galactic Potential Models · Jeans Equations · Orbits in Axisymmetric and Triaxial Potentials.

Ask the tutor: ask me to walk through choosing $N$ , $\epsilon$ , and $\Delta t$ for a specific system — a globular cluster, a galaxy halo, or a cosmological box — and we will size the run together.

Why not just use direct summation for everything, since it is exact?

Because its cost grows as $N^2$ , and the systems astrophysicists most want to simulate are enormous — Worked Example 1 shows a single direct step for $10^{10}$ particles would take of order 1600 core-years, while a tree does the same physics 300 million times faster. Direct summation stays the choice only for collisional systems of modest size — clusters of $10^4$ to $10^6$ members — where two-body relaxation is real physics and softening would wrongly erase the close encounters. For galaxies and cosmological volumes, $\mathcal{O}(N\log N)$ solvers are the only viable option.

How do I choose the softening length ε?

It is a compromise, not a free parameter. Too small and the force is barely capped, so an under-sampled system scatters artificially and its core heats; too large and you smear out real structure below $\epsilon$ , flattening cusps and dissolving subhalos. A common rule ties $\epsilon$ to a small fraction of the mean interparticle separation in the densest resolved regions. Crucially it must be chosen jointly with $N$ : more particles without smaller $\epsilon$ wastes resolution, smaller $\epsilon$ without more particles reintroduces relaxation noise. The honest test is convergence — vary $\epsilon$ and confirm the result does not.

Isn't a run with 10⁹ particles still absurdly under-sampling a galaxy of 10¹¹ stars?

In number, yes — and in a dark-matter halo the mismatch is far larger still. But that is the whole point of the collisionless framework. The particles are not individual stars; they are Monte Carlo samples of a smooth phase-space distribution function obeying the collisionless Boltzmann equation. What must be faithful is the mean-field potential and its evolution, not the identity of any particle. Softening and large $N$ together keep the sampling noise — quantified by the numerical relaxation time of equation (7) — negligible over the run, so the coarse sampling tracks the smooth fluid it represents.

Why leapfrog rather than a high-order Runge–Kutta integrator?

Because for gravity the integrator's geometric structure matters more than its formal order. Leapfrog is only second order and RK4 has a much smaller per-step error, yet leapfrog is symplectic: it exactly integrates a shadow Hamiltonian differing from the true one by $\mathcal{O}(\Delta t^2)$ , so the energy error oscillates within a fixed band forever. RK4 has no conserved quantity; its energy drifts secularly until, over millions of steps, orbits spiral in or out. Add that leapfrog needs one force evaluation per step against RK4's four, and it is both cheaper and more trustworthy.

How do cosmic expansion and periodic boundaries enter the equations?

Cosmological runs use comoving coordinates: the background expansion is factored out, the scale factor $a(t)$ enters the equations of motion, and the particles track only the peculiar motions driven by density fluctuations. The Poisson equation is written for the fluctuation potential, which is why the $\mathbf{k}=0$ mode is dropped in the PM solve. Periodic boundaries model an infinite universe with a finite box — a particle leaving one face re-enters the opposite one, enforced by the FFT and by Ewald summation over image lattices in tree codes — so a $100$ – $1000\ \mathrm{Mpc}$ box is a fair, homogeneous sample of the cosmos.

How do you know a result is real physics and not a numerical artifact?

By convergence testing and conservation checks, never by trusting a single run. Standard health checks: total energy stays constant (bounded oscillation, no drift), and an equilibrium system satisfies the virial ratio $2T/|W|\to1$ — a Plummer or Hernquist realization should sit still, not relax or evaporate. Then vary $N$ , $\epsilon$ , $\Delta t$ , and $\theta$ and confirm the result does not change. A feature that weakens as $N$ grows or $\epsilon$ shrinks is discreteness noise, not physics; only what survives convergence can be believed.