Nonlinear Structure Formation and N-body Simulations

Linear perturbation theory is a triumph of cosmology, but it is a theory of the infinitesimal: it holds only while the density contrast $\delta=\delta\rho/\bar\rho$ stays far below one. The moment a region grows to $\delta\sim1$ that description fails. Fourier modes that evolved independently begin to couple, cold streams of matter cross and fold, and violent relaxation scrambles orbits into a bound object. Everything that makes the Universe interesting to look at — galaxies, clusters, and the filamentary cosmic web threaded between the voids — lives in this nonlinear regime, and for a cold, pressureless, collisionless fluid there is no general analytic solution to reach it. So cosmologists run the experiment numerically. The dark matter is a smooth phase-space distribution $f(\mathbf{x},\mathbf{v},t)$ obeying the collisionless Boltzmann (Vlasov) equation coupled self-consistently to Poisson's equation. A six-dimensional grid is hopeless, so one Monte-Carlo samples $f$ with $N$ tracer particles, each coasting along a characteristic of the Vlasov flow — nothing more than Newton's law in an expanding box. Each particle is a parcel of dark matter, not a physical WIMP; its mass fixes the resolution and a softening length tames the short-range force. The whole art is then computing the mutual gravity of $N$ bodies fast enough. Direct summation is exact but scales as $N^2$ ; tree codes group distant matter into multipoles for $N\log N$ ; particle-mesh methods solve Poisson on an FFT grid; hybrids marry the two. The outputs — halo catalogues, profiles, mass functions — are checked against spherical collapse, Press–Schechter theory, and the NFW profile, and against galaxy surveys, lensing, and the baryon acoustic peak. This is the machinery that turns $\Lambda$ CDM from a parameter set into a universe.

Key takeaways

Linear theory dies at $\delta\sim1$ ; beyond it, mode coupling, shell crossing, and virialization make structure formation intrinsically nonlinear, with no general closed-form solution.
The governing physics is the Vlasov (collisionless Boltzmann) equation coupled to Poisson; the practical solution samples the phase-space distribution with $N$ particles following its characteristics — Newtonian gravity in comoving coordinates.
Gravity is evaluated by direct summation ( $O(N^2)$ ), tree codes ( $O(N\log N)$ ), particle-mesh FFT solvers, or TreePM/P $^3$ M hybrids, and time is advanced by a symplectic leapfrog.
Spherical collapse supplies the analytic anchors every code is calibrated against: turnaround overdensity $9\pi^2/16\approx5.55$ , linear threshold $\delta_c=\tfrac{3}{20}(12\pi)^{2/3}\approx1.686$ , and virial overdensity $\Delta_\mathrm{vir}=18\pi^2\approx178$ .
Halos emerge with a near-universal NFW profile and an abundance described by Press–Schechter and Sheth–Tormen theory, assembled into the cosmic web of sheets, filaments, and knots.
Key idea: nonlinear collapse is followed by sampling the collisionless phase-space distribution with $N$ particles moving along the characteristics of the Vlasov equation in an expanding background, with self-gravity solved by tree, PM, or TreePM algorithms; spherical collapse and Press–Schechter theory are the analytic scaffolding.

The core insight

Nonlinear collapse is followed by sampling the collisionless phase-space distribution with $N$ particles that move along the characteristics of the Vlasov equation in an expanding background, with self-gravity solved via tree/PM/TreePM algorithms; the spherical collapse model and Press–Schechter theory provide the analytic scaffolding.

Where this lives today

Cosmological $N$ -body simulation is now an industrial-scale enterprise at the centre of precision cosmology. Flagship dark-matter runs — Millennium, Bolshoi, MultiDark, the exascale AbacusSummit and Uchuu suites following trillions of particles — provide the templates against which surveys are read. When DESI measures the baryon acoustic peak, or Euclid and the Vera Rubin Observatory's LSST map weak-lensing shear across billions of galaxies, the interpretation runs through mock catalogues that predict the nonlinear power spectrum, the halo mass function, and their covariance to percent accuracy. Because a full simulation per trial cosmology is impossibly expensive, the field has turned to emulators — machine-learning surrogates that interpolate observables across parameter space in milliseconds. Meanwhile the small-scale predictions of collisionless CDM are where the model is most sharply tested: the "core–cusp" tension and the "missing satellites" and "too-big-to-fail" problems drive both baryonic-feedback modelling and simulations of warm, fuzzy, and self-interacting dark matter. The largest hydrodynamic runs — IllustrisTNG, EAGLE, SIMBA — fold gas, stars, black holes, and feedback onto the dark-matter skeleton, closing the loop to observed galaxies.

The physical picture

Picture the dark matter as a smooth, cold fluid with no thermal pressure to hold it up. Wherever it is slightly overdense, its extra self-gravity slows the local expansion; the region lags behind the Hubble flow, turns around, and falls back on itself. In linear theory each Fourier mode of $\delta$ grows by the same factor $D(a)$ and never notices its neighbours; nonlinearity is the moment the modes begin to talk to each other and a patch detaches to become a self-bound halo. Because the fluid is collisionless it passes straight through itself, so the first structures are not blobs but sheets — matter falls onto a plane, overshoots, and forms a caustic "pancake." Pancakes drain into filaments, filaments into knots, and the knots are the halos where galaxies light up: the cosmic web that galaxy redshift surveys map. The single most important quantitative anchor is the spherical collapse model, and the first figure is its whole story. Treat one overdense sphere as a closed sub-universe: it expands, turns around, and collapses along a cycloid. Real matter never reaches the singular point; it virializes at half the turnaround radius. Following the numbers through gives three constants every simulation is measured against — turnaround overdensity $9\pi^2/16\approx5.55$ , virial overdensity $18\pi^2\approx178$ , and linear collapse threshold $\delta_c\approx1.686$ .

Figure

Spherical collapse: the cycloid R(t) expands to turnaround at θ = π where 1 + δ = 9π²/16 ≈ 5.55, then collapses toward θ = 2π; virializing at half the turnaround radius gives mean overdensity Δ_vir ≈ 18π² ≈ 178, while linear theory extrapolated to collapse reaches only δ_c ≈ 1.686

An $N$ -body simulation realizes exactly this collapse for millions of overlapping spheres at once. Lay down $N$ particles displaced from a uniform grid by the Zel'dovich approximation, seed the displacements from the linear power spectrum $P(k)$ , and integrate Newtonian gravity in an expanding box from $z\sim50$ to today. The only real obstacle is speed: summing the mutual attraction of every tracer at every timestep by brute force costs $N^2$ operations — ruinous for the $N\sim10^9$ to $10^{12}$ particles of a modern run. The second figure is why the field exists in its present form: hierarchical and mesh algorithms replace $N^2$ with $N\log N$ , turning centuries of computing into days.

Figure

Why N-body gravity needs clever algorithms: direct summation costs O(N²) per step while tree and TreePM cost O(N log N); for N = 10⁹ that is 10¹⁸ versus ≈ 3×10¹⁰ operations — a speedup of about 3×10⁷ separating the feasible from the impossible

The core equations

The fundamental object is not a particle but a distribution. Dark matter fills a six-dimensional phase space, and its evolution under gravity alone is governed by the collisionless Boltzmann — or Vlasov — equation, which states that the phase-space density is conserved along every trajectory.

\frac{\partial f}{\partial t}+\frac{\mathbf{p}}{a^2 m}\cdot\frac{\partial f}{\partial\mathbf{x}}-m\,\nabla\Phi\cdot\frac{\partial f}{\partial\mathbf{p}}=0 \qquad (1)

where: $f(\mathbf{x},\mathbf{p},t)$ — phase-space distribution of the dark matter; $\mathbf{x}$ — comoving position; $\mathbf{p}=a^2m\dot{\mathbf{x}}$ — comoving canonical momentum; $a$ — scale factor; $m$ — particle mass; $\Phi$ — peculiar potential. Reading it: $f$ is conserved in phase space, only sheared and folded by gravity — the basis for solving it with tracer particles. The potential is sourced by the matter itself through Poisson's equation. In comoving coordinates only departures from the mean gravitate, so the source is the density contrast $\delta$ , not the full density.

\nabla^2\Phi=4\pi G a^2\,\bar\rho\,\delta=\frac{3}{2}\,\Omega_m H^2 a^2\,\delta \qquad (2)

where: $\Phi$ — peculiar potential; $\delta=\delta\rho/\bar\rho$ — density contrast; $\bar\rho$ — mean matter density; $\Omega_m$ — matter density parameter; $H$ — Hubble parameter. Reading it: a uniform universe produces no peculiar force ( $\delta=0\Rightarrow\nabla^2\Phi=0$ ) — structure gravitates only where matter clumps. Each tracer obeys the characteristic of (1) — Newton's second law dressed for an expanding universe, with a "Hubble drag" term that damps peculiar velocities as the box stretches.

\ddot{\mathbf{x}}+2H\dot{\mathbf{x}}=-\frac{1}{a^2}\nabla\Phi \qquad (3)

where: $\mathbf{x}$ — comoving coordinate; $2H\dot{\mathbf{x}}$ — Hubble drag; $a^{-2}\nabla\Phi$ — peculiar acceleration. Reading it: with $\Phi=0$ peculiar velocities redshift away as $\dot{\mathbf x}\propto a^{-2}$ ; with $H=0$ it is ordinary Newtonian dynamics. Evaluating $\nabla\Phi$ particle-by-particle gives the softened direct-summation force. The softening length $\epsilon$ flattens the $1/r^2$ divergence at short range so close-passing tracers do not scatter violently — an artifact of finite $N$ , not real physics.

\mathbf{F}_i=-Gm_i\sum_{j\neq i}m_j\,\frac{\mathbf{r}_{ij}}{(r_{ij}^2+\epsilon^2)^{3/2}} \qquad (4)

where: $\mathbf{r}_{ij}$ — separation between particles $i$ and $j$ ; $\epsilon$ — softening length; $m_i,m_j$ — tracer masses. Reading it: for $r_{ij}\gg\epsilon$ this is exact Newtonian gravity; for $r_{ij}\lesssim\epsilon$ the force saturates instead of diverging. The sum in (4) is what costs $O(N^2)$ . Barnes–Hut tree codes avoid it by lumping a distant clump into one multipole, opening a cell only when it subtends too large an angle.

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

where: $s$ — node (cell) size; $d$ — distance from the node's centre of mass to the target; $\theta$ — opening angle, typically $0.5$ to $0.7$ . Reading it: a small, faraway cell ( $s/d<\theta$ ) is replaced by one multipole, so each particle sees only $O(\log N)$ nodes — the scaling of the second figure. Spherical collapse supplies the calibration constants. The first is the linear collapse threshold: extrapolate the linear growth of a perturbation to the moment its exact nonlinear counterpart collapses, and in an Einstein–de Sitter universe it has reached this universal value.

\delta_c=\frac{3}{20}(12\pi)^{2/3}\approx1.686 \qquad (6)

where: $\delta_c$ — linearly extrapolated overdensity at collapse. Reading it: a bookkeeping number, not a real density — the barrier a linear field must cross for the true field to have collapsed, feeding the mass function below. The second constant is the actual mean overdensity of the settled halo relative to the background — the number that defines a halo's edge.

\Delta_\mathrm{vir}=\frac{\bar\rho_\mathrm{halo}}{\bar\rho}\approx18\pi^2\approx178 \qquad (7)

where: $\Delta_\mathrm{vir}$ — mean halo overdensity relative to background at virialization. Reading it: halos are defined as spheres enclosing $178$ (often rounded to $200$ ) times the mean or critical density, giving $R_{200}$ and $M_{200}$ . Inside those halos, simulations find a near-universal shape — the Navarro–Frenk–White profile — independent of mass and cosmology, with a $1/r$ cusp at the centre and a $1/r^3$ falloff outside the scale radius.

\rho(r)=\frac{\rho_s}{(r/r_s)\,(1+r/r_s)^2} \qquad (8)

where: $\rho_s$ — characteristic density; $r_s$ — scale radius; concentration $c=R_{200}/r_s$ . Reading it: two numbers — or equivalently mass and concentration — fix the entire density run of a relaxed halo. Halo abundance versus mass is captured by the Press–Schechter mass function, the statistical bridge from the initial Gaussian field to the collapsed-object count.

\frac{\mathrm{d}n}{\mathrm{d}M}=-\sqrt{\frac{2}{\pi}}\,\frac{\bar\rho}{M}\,\frac{\delta_c}{\sigma^2}\,\frac{\mathrm{d}\sigma}{\mathrm{d}M}\,\exp\!\left(-\frac{\delta_c^2}{2\sigma^2}\right) \qquad (9)

where: $n$ — comoving halo number density; $M$ — halo mass; $\sigma(M)$ — rms fluctuation smoothed on scale $M$ ; $\delta_c\approx1.686$ . Reading it: the exponential cutoff makes cluster abundance exquisitely sensitive to $\sigma_8$ — the lever that turns halo counts into a cosmological measurement. The motion is advanced by the symplectic leapfrog integrator in kick–drift–kick form ( $\mathbf{v}_{n+1/2}=\mathbf{v}_n+\mathbf{a}_n\Delta t/2$ , then $\mathbf{x}_{n+1}=\mathbf{x}_n+\mathbf{v}_{n+1/2}\Delta t$ , then a second half-kick). Because it conserves phase-space volume and is time-reversible, it shows no secular energy drift over the billions of steps of a run — which is why it beats a high-order Runge–Kutta scheme.

Derivation ASpherical collapse — turnaround, virialization, and the number 178

The goal is to derive, from Newton's law alone, the three constants that calibrate every cosmological simulation: the turnaround overdensity $9\pi^2/16$ , the virial overdensity $\Delta_\mathrm{vir}=18\pi^2\approx178$ , and the linear collapse threshold $\delta_c=\tfrac{3}{20}(12\pi)^{2/3}\approx1.686$ . The trick is that a uniform overdense sphere in an Einstein–de Sitter background is itself a closed Friedmann universe, so its exact evolution is the familiar cycloid — and comparing that cycloid to the background at matched times converts a one-line orbit into the numbers cosmologists quote.

Step 1. Take a uniform sphere of mass $M$ and radius $R$ . By Newton's shell theorem the boundary obeys $\ddot R=-GM/R^2$ ; first-integrating gives $\tfrac12\dot R^2-GM/R=E$ , and for an overdensity $E<0$ , so the orbit is bound and must turn around. Why: stronger self-gravity than the background makes the sphere bound — collapse is guaranteed.

Step 2. The bound solution is the cycloid $R=A(1-\cos\theta)$ , $t=B(\theta-\sin\theta)$ with $A^3=GMB^2$ . The angle $\theta$ runs from $0$ at the Big Bang to $\pi$ at turnaround ( $R_\mathrm{ta}=2A$ ) and on to $2\pi$ where $R\to0$ . Why: the same parametric solution as a closed Friedmann universe, so the sphere re-collapses like a mini big crunch.

Step 3. Compare densities at turnaround. The sphere has $\rho_\mathrm{ta}=3M/(4\pi R_\mathrm{ta}^3)$ ; the EdS background at $t_\mathrm{ta}=B\pi$ has $\bar\rho=1/(6\pi G t^2)$ . Substituting $A^3=GMB^2$ , $R_\mathrm{ta}=2A$ gives $1+\delta_\mathrm{ta}=9\pi^2/16\approx5.55$ . Why: the region is already five times the mean density when it stops expanding — the first hard number.

Step 4. A collisionless sphere cannot reach $R=0$ ; it virializes. The virial theorem $2K+U=0$ with energy conservation $E=U_\mathrm{ta}$ (all potential at turnaround, where $\dot R=0$ ) forces $R_\mathrm{vir}=\tfrac12 R_\mathrm{ta}$ , raising the density by $2^3=8$ . Why: orbits phase-mix and settle at half the turnaround radius — this is where a real halo lives.

Step 5. The background keeps diluting between turnaround ( $\theta=\pi$ ) and formal collapse ( $\theta=2\pi$ ), a factor $(t_\mathrm{coll}/t_\mathrm{ta})^2=4$ . Multiplying: $\Delta_\mathrm{vir}=8\times(9\pi^2/16)\times4=18\pi^2\approx178$ . Why: the halo is measured against a thinned-out background, boosting the contrast to the canonical $178$ .

Step 6. Separately, expand $1+\delta$ to lowest order in $\theta$ for the linear growth $\delta_\mathrm{lin}=\tfrac{3}{20}(6\pi)^{2/3}(t/t_\mathrm{coll})^{2/3}$ , and evaluate at $t=t_\mathrm{coll}$ : $\delta_c=\tfrac{3}{20}(12\pi)^{2/3}\approx1.686$ . Why: the value linear theory reports at collapse — the Press–Schechter threshold, cleanly separated from the true overdensity $178$ . *The exact EdS values are $9\pi^2/16=5.5517$ , $18\pi^2=177.65$ , and $\tfrac{3}{20}(12\pi)^{2/3}=1.6865$ . In $\Lambda$ CDM the values shift slightly: $\delta_c\approx1.675$ and $\Delta_\mathrm{vir}\approx100$ relative to critical density.*

Go deeper · Derivation B · 1 minThe Press–Schechter halo mass function

Derivation BThe Press–Schechter halo mass function

The goal is to predict how many collapsed halos of each mass exist, starting only from a Gaussian initial density field and the threshold $\delta_c$ derived above. Press and Schechter's idea sidesteps the intractable nonlinear dynamics entirely: smooth the linear field on a mass scale $M$ , declare any region above $\delta_c$ collapsed into a halo of at least that mass, and count. It works better than it has any right to, and repairing its one embarrassing factor of two runs from a 1974 fudge to modern excursion-set theory.

Step 1. Smooth the linear density field with a filter enclosing mass $M$ ; because the field is Gaussian with variance $\sigma^2(M)$ , the fraction of space above the collapse barrier is $P(>\delta_c)=\tfrac12\,\mathrm{erfc}\!\left(\delta_c/\sqrt2\,\sigma(M)\right)$ . Why: identifying "collapsed on scale $M$ " with "linear field exceeds $\delta_c$ when smoothed on $M$ " converts a dynamics problem into a Gaussian tail integral.

Step 2. As $\sigma\to\infty$ (no smoothing) this fraction tends to $\tfrac12$ , not $1$ — half of all mass appears never to collapse, which is absurd. Press and Schechter multiply by $2$ , arguing underdense material is swept into larger collapsed structures. Why: the "cloud-in-cloud" problem; the factor of two, a 1974 heuristic, was later justified by the excursion-set formalism.

Step 3. Differentiate the corrected fraction with respect to mass and multiply by $\bar\rho/M$ to get a number density per unit mass — equation (9), compactly $\tfrac{\mathrm{d}n}{\mathrm{d}M}\,\mathrm{d}M=\sqrt{\tfrac{2}{\pi}}\,\tfrac{\bar\rho}{M}\,\nu\,e^{-\nu^2/2}\,\tfrac{\mathrm{d}\nu}{\nu}\,\tfrac{\mathrm{d}M}{M}$ , $\nu=\delta_c/\sigma(M)$ . Why: the exponential $e^{-\nu^2/2}$ makes massive clusters (large $\nu$ ) exponentially rare and so cosmology-sensitive.

Step 4. Calibrating against simulations, PS overpredicts low-mass and underpredicts high-mass halos; Sheth and Tormen fix this with ellipsoidal collapse, a moving barrier $\nu'^2=a\nu^2$ with $a\approx0.707$ , $p\approx0.3$ . Why: real peaks collapse along their shortest axis first, so the effective threshold depends on scale — the correction that makes the mass function percent-accurate. *The form reduces to standard PS with $\delta_c=1.686$ , and by construction $\int (M/\bar\rho)\,\mathrm{d}n=1$ , so all mass is accounted for.*

Go deeper · Derivation C · 1 minFrom the Vlasov equation to N bodies

Derivation CFrom the Vlasov equation to N bodies

The goal is to justify the central sleight of hand of the whole field: that throwing down $N$ point particles and integrating Newtonian gravity actually solves the continuous Vlasov–Poisson system (1) and (2). It is not obvious — a smooth six-dimensional distribution is replaced by a cloud of discrete points — and the justification rests on the method of characteristics, which turns the PDE for $f$ into ordinary differential equations for trajectories.

Step 1. The Vlasov equation (1) says $\mathrm{d}f/\mathrm{d}t=0$ : the convective derivative vanishes, so phase-space density is conserved following a fluid element — Liouville's theorem for a collisionless system. Why: with no collisions to scatter particles between cells, the fluid density is frozen along its streamlines.

Step 2. Those streamlines — the characteristics — are exactly Hamilton's equations $\dot{\mathbf{x}}=\mathbf{p}/(a^2m)$ , $\dot{\mathbf{p}}=-m\nabla\Phi$ , the single-particle equations of motion (3). Why: solving the PDE reduces to following trajectories that are just the orbit of one particle in the mean field.

Step 3. Sample the initial $f$ with $N$ points, advance each along its characteristic with leapfrog, and recompute $\Phi$ from their density via Poisson (2) each step. The evolving cloud is a Monte-Carlo estimate of $f$ . Why: because $f$ is conserved along characteristics, the tracer density tracks the true $f$ — the $N$ -body method, whose error falls as sampling improves.

Step 4. The loop closes: particles source the potential, the potential accelerates the particles. The only unphysical ingredient is discreteness, controlled by softening $\epsilon$ and by using $N$ large enough that the relaxation time exceeds the age of the universe. Why: exact as $N\to\infty$ , $\epsilon\to0$ ; finite- $N$ noise is the price of tractability and must be quantified. *The characteristics of (1) reproduce the equations of motion (3), and the collisionless limit is recovered as $N\to\infty$ .*

Limiting cases

Nonlinear structure formation is bracketed by regimes that are individually well understood; a correct simulation must reproduce each in its corner of parameter space.

Case	Result	Meaning
Linear regime $\delta\ll1$	$\delta(\mathbf{x},a)=D(a)\,\delta_0(\mathbf{x})$	Every mode grows by the same factor $D(a)$ , independently — a code must reproduce this before nonlinearity sets in, the primary validation test.
Quasi-linear (Zel'dovich)	$\mathbf{x}(t)=\mathbf{q}+D(t)\,\boldsymbol{\psi}(\mathbf{q})$	First-order Lagrangian theory: particles stream along their initial displacement field until trajectories cross and form caustic pancakes. Used to seed initial conditions.
Fully nonlinear	$\delta\gtrsim1,\;\Delta_\mathrm{vir}\approx178$	Virialized halos, tractable only numerically. The overdensity threshold defines the halo boundary ( $R_{200}$ ) and enclosed mass.
Collisionless limit ( $N\to\infty$ )	$t_\mathrm{relax}\to\infty$	The physical target: dark matter never two-body relaxes, retaining its phase-space structure (streams, caustics) indefinitely.
Discrete artifact (finite $N$ )	$t_\mathrm{relax}\sim\dfrac{N}{8\ln N}\,t_\mathrm{cross}$	Finite $N$ injects spurious relaxation and heating; kept negligible by large $N$ and softening $\epsilon$ .

Common trapThe N particles are real dark-matter particles, and δ_c = 1.686 is a real overdensity

Both halves are wrong, and both are common. First, the $N$ bodies are not physical dark-matter particles — a single one stands in for $10^{60}$ or more actual particles. They are Monte-Carlo tracers of a smooth collisionless fluid, which is why two may pass right through each other (real dark matter streams multiply, it does not collide) and why any close-encounter scattering is a numerical artifact to suppress with softening, not physics. Second, $\delta_c=1.686$ is never an actual density anywhere. It is the value linear theory — deliberately extrapolated past its own validity — would assign to a perturbation at the instant the true nonlinear object collapses. The genuine mean overdensity of the halo is $\Delta_\mathrm{vir}\approx178$ , roughly a hundred times larger. Confusing the two is the single most frequent error in this subject.

Worked examples

1 · Why a billion-body run needs a tree (fully worked)

A simulation follows $N=10^9$ particles. Compare the cost per timestep of direct summation against a TreePM solver, and translate the difference into wall-clock time on a petascale machine sustaining $10^{15}$ floating-point operations per second.

Step 1. Direct summation touches every pair: $\sim N^2=(10^9)^2=\mathbf{10^{18}}$ operations per force evaluation. On a $10^{15}$ flop/s machine one step takes $10^{18}/10^{15}=10^3$ s, about 17 minutes — for a single timestep, of thousands.

Step 2. A tree/TreePM solver scales as $N\log_2 N$ . With $\log_2(10^9)=29.9$ , that is $10^9\times29.9\approx\mathbf{3.0\times10^{10}}$ operations — tens of microseconds of raw arithmetic per step (real codes are slower from memory traffic, but the scaling is the point).

Step 3. The ratio is $10^{18}/3.0\times10^{10}\approx\mathbf{3.3\times10^{7}}$ : the tree does the same physics with thirty million times fewer operations. Over $\sim10^3$ steps, direct summation would burn $\sim10^{21}$ operations (months of supercomputer time) while the tree finishes in seconds. This ~3×10⁷ factor is why cosmological N-body simulation is an O(N log N) discipline, not O(N²).

2 · Is a simulated halo collisionless? (fully worked)

A dark-matter halo is resolved with $N=10^6$ particles and has a crossing time $t_\mathrm{cross}\approx1$ Gyr. Estimate its two-body relaxation time and compare to the age of the universe, $t_H\approx1.4\times10^{10}$ yr.

Step 1. The standard estimate is $t_\mathrm{relax}\approx\dfrac{N}{8\ln N}\,t_\mathrm{cross}$ . With $N=10^6$ , $\ln N=13.8$ , so the prefactor is $10^6/(8\times13.8)=10^6/110\approx\mathbf{9.0\times10^3}$ .

Step 2. Multiply by the crossing time: $t_\mathrm{relax}\approx9.0\times10^3\times1\ \mathrm{Gyr}=9.0\times10^3\ \mathrm{Gyr}=\mathbf{9\times10^{12}\ yr}$ .

Step 3. Compare: $t_\mathrm{relax}/t_H=9\times10^{12}/1.4\times10^{10}\approx\mathbf{650}$ . The halo relaxes $650$ times more slowly than the universe is old, so two-body relaxation is negligible and the halo faithfully samples a collisionless system. The scaling matters: going to $N=10^7$ pushes $t_\mathrm{relax}$ up by roughly $8.6$ , buying deeper collisionless fidelity. Large N is not vanity; it is what keeps the numerics collisionless.

Physics of the portrait

Nonlinear structure formation was built by many hands, but no single career embodies its arc — from the first crude simulations to the industrial runs that define $\Lambda$ CDM — better than Simon White's.

Origins · 1 minSimon White and how this was discovered

Simon White

born 1951 · cosmological structure formation and N-body simulations

Simon White is a British-German cosmologist who, more than almost anyone, made numerical simulation the arbiter of how cosmic structure forms. In 1978, with Martin Rees, he laid out the two-stage picture in which gas cools inside dark-matter halos that grow by hierarchical merging — the skeleton of modern galaxy formation. His defining contribution came in 1985, when Marc Davis, George Efstathiou, Carlos Frenk, and White published the "DEFW" simulations of a cold dark matter universe. Tiny by today's standards, they showed for the first time that CDM plus gravity naturally produces the observed web of filaments and voids, effectively selecting CDM over hot dark matter by direct comparison with galaxy clustering. White then chased the halos' internal structure: with Julio Navarro and Frenk he found, in 1996 and 1997, that relaxed halos of every mass share one universal density profile — the NFW form of equation (8), so robust it became a diagnostic of the dark matter's cold, collisionless nature. As a director of the Max Planck Institute for Astrophysics he oversaw the 2005 Millennium simulation, which followed ten billion particles and seeded a generation of mock catalogues. Threaded through it all is one conviction: a cosmological model is only as good as the nonlinear, observable predictions you can extract from it numerically — the simulation, not the pencil, is where theory meets the sky.

History

1970

— Peebles runs early N-body simulations of galaxy clustering

1974

— Press and Schechter predict the halo mass function from a Gaussian field

1985

— Davis, Efstathiou, Frenk and White simulate the cold dark matter cosmogony

1986

— Barnes and Hut publish the tree algorithm for O(N log N) gravity

2005

— Springel's Millennium simulation follows over ten billion particles

Going further

Texts. Carroll & Ostlie, An Introduction to Modern Astrophysics, 2nd ed., Ch. 30 (structure formation) · Ryden, Introduction to Cosmology, 2nd ed., Ch. 12 (nonlinear collapse) · Dodelson, Modern Cosmology, Ch. 7 and Ch. 9 (nonlinear structure) · Mo, van den Bosch & White, Galaxy Formation and Evolution, Ch. 5 to 7 (spherical collapse, mass functions, halos) · Binney & Tremaine, Galactic Dynamics, 2nd ed., Ch. 4 (collisionless systems). Spherical collapse yields $9\pi^2/16=5.5517$ , $\delta_c=\tfrac{3}{20}(12\pi)^{2/3}=1.6865$ , and $\Delta_\mathrm{vir}=18\pi^2=177.65$ . The Sheth–Tormen parameters are $a\approx0.707$ , $p\approx0.3$ . The cost ratio for $N=10^9$ is $10^{18}/(10^9\log_2 10^9)\approx3\times10^7$ , and the relaxation-to-crossing-time ratio is $t_\mathrm{relax}/t_\mathrm{cross}=N/(8\ln N)\approx9\times10^3$ at $N=10^6$ .

Prerequisites. Linear Growth of Structure · Density Perturbations and the Density Contrast · Jeans Instability (Cosmological) · Transfer Function and the Matter Power Spectrum · Cold Dark Matter. Leads to. The ΛCDM Concordance Model · Baryon Acoustic Oscillations · Dark Energy and the Cosmological Constant.

Ask the tutor: Want to walk through how the spherical-collapse cycloid turns into the numbers 1.686 and 178, or how a tree code actually opens cells to hit O(N log N)?

Are the simulation particles real dark-matter particles?

No. A single simulation particle has a mass of $10^6$ to $10^{10}$ solar masses — it stands in for a parcel of $10^{60}$ or more actual particles. The $N$ bodies are Monte-Carlo samples of a smooth phase-space distribution, chosen to trace its density and velocity field. That is why two can stream through one another: real cold dark matter is collisionless and multi-streaming, not a gas of colliding balls. Any genuine two-body scattering between them is numerical noise, suppressed by softening.

What does the softening length do, and how big should it be?

The softening length $\epsilon$ replaces the singular $1/r^2$ force at short range with a smooth, finite one, as in equation (4), preventing artificial large-angle scattering when two tracers pass close — spurious heat that would mimic a collisional system. Too small reintroduces two-body noise and tiny timesteps; too large washes out substructure and flattens halo centres. In practice $\epsilon$ is a few percent of the mean interparticle spacing — a key convergence parameter to test.

Why do simulations start at high redshift like z ~ 50 rather than at the Big Bang?

At very high redshift $\delta$ is tiny and known analytically — linear theory and the Zel'dovich approximation give the exact displacement field from $P(k)$ ; starting earlier only buys noise, as displacements shrink below the grid discreteness. Starting later (say $z=5$ ) is wrong: the largest scales have already gone mildly nonlinear, so linear initial conditions no longer apply. The choice $z\sim50$ to $100$ is the sweet spot where $\delta\ll1$ everywhere yet displacements are cleanly resolved.

What is the practical difference between tree, PM, and TreePM force solvers?

A pure tree code (Barnes–Hut) multipoles distant matter — accurate at all scales, $O(N\log N)$ , but expensive per particle at short range. A particle-mesh (PM) code deposits mass on a grid and solves Poisson by FFT — very fast, but its resolution is capped at one grid cell, so it cannot resolve halo interiors. TreePM marries them: the FFT mesh handles the smooth long-range force while a local tree supplies the sharp short-range force. Almost every modern large simulation — GADGET, for instance — uses TreePM.

Why does the cosmic web form sheets and filaments rather than just round clumps?

Because gravitational collapse is generically anisotropic. A random overdense region has three unequal tidal axes, so it collapses along its shortest axis first into a two-dimensional sheet (a Zel'dovich pancake), then along the second into a one-dimensional filament, and only where filaments intersect does the third axis collapse into a knot. This ordered sequence — sheet, filament, knot — is what the Zel'dovich approximation predicts and what surveys display as the cosmic web. Spherical collapse is an idealization for thresholds, not a map of where matter ends up.

Do baryons matter, or is a dark-matter-only simulation enough?

For the large-scale skeleton, dark-matter-only runs are excellent: dark matter dominates the mass, so the cosmic web, the mass function, and the BAO scale are set by gravity on the collisionless component. But baryons matter decisively inside halos — gas cools, forms stars, and is blasted out by supernovae and AGN, feedback that can turn a cusp into a core and shift the power spectrum at the percent level next-generation lensing measures. So dark-matter-only is the workhorse for large-scale structure; full hydrodynamic runs (IllustrisTNG, EAGLE, SIMBA) are needed for galaxies themselves.