Nonlinear Gravitational Clustering

Linear growth preserves independent Fourier modes. Nonlinear gravity does not: continuity and advection couple wavevectors, overdensities collapse anisotropically, and a cold single stream eventually folds into multiple streams. This page develops that transition through Eulerian and Lagrangian descriptions, then uses a one-dimensional Zel’dovich map to locate the first shell crossing reproducibly. The map tests geometry and convergence; it is not an $N$-body simulation.

Helpful background. The Einstein Equation supplies Einstein Equations and their weak-field limit; Relativistic Matter and Compact Objects supplies Relativistic Matter and Compact Objects kinetic and fluid closures; Relativistic Cosmology supplies the Relativistic Cosmology background; Computational Relativity supplies convergence and benchmark practice; and Probability and Statistical Inference supplies random-field statistics and sampling uncertainty.

From Vlasov flow to mode coupling

On subhorizon scales, in the weak-field and nonrelativistic limit of general relativity, cold collisionless matter is described by the Vlasov–Poisson system. Taking momentum moments gives density, bulk velocity, and a velocity-dispersion tensor. Before shell crossing, a negligibly warm population has one velocity at each position, so the dispersion may be set to zero. In conformal time $\eta$, let $\mathbf u=d\mathbf x/d\eta$ be the physical peculiar velocity, $\mathcal H=a'/a$, and $\nabla$ differentiate comoving position. Then

$$\delta'+\nabla\cdot[(1+\delta)\mathbf u]=0,$$ $$\mathbf u'+\mathcal H\mathbf u +(\mathbf u\cdot\nabla)\mathbf u=-\nabla\Phi, \qquad \nabla^2\Phi=4\pi Ga^2\bar\rho_m\delta.$$

Linearization removes the products and leaves independent Fourier modes. In a general background, nonlinear perturbations instead have time-dependent mode-coupling kernels,

$$\delta(\mathbf k,\eta) =\sum_{n=1}^{\infty}\delta^{(n)}(\mathbf k,\eta),$$

where each $\delta^{(n)}$ is a convolution of $n$ initial modes whose wavevectors sum to $\mathbf k$, with a kernel $F_n(\eta;\mathbf q_1,\ldots,\mathbf q_n)$. For the fastest-growing solution in Einstein–de Sitter, the time dependence factorizes exactly:

$$\delta^{(n)}(\mathbf k,\eta) =D_+^n(\eta)\,\delta_{\rm EdS}^{(n)}(\mathbf k).$$

Using these EdS kernels with the exact $D_+$ is often an accurate $\Lambda$CDM approximation, but it is not an identity in an arbitrary cosmology. In either form, nonlinear evolution transfers power and phase information between scales and generates non-Gaussian higher-point correlations even from Gaussian initial conditions.

The moment equations, single-stream assumption, linear modes, and nonlinear Fourier kernels are developed in Bernardeau et al. 2002, §§2.1–2.4, printed pp. 13–27, eqs. (9)–(87) (Open PDF). Pressureless perturbation theory is controlled while fluctuations and neglected stress remain small. It does not acquire a physical continuation through multistreaming merely by adding more formal orders.

The Lagrangian map and shell crossing

Label a fluid element by its initial coordinate $\mathbf q$ and write

$$\mathbf x(\mathbf q,D) =\mathbf q+D\,\mathbf s(\mathbf q),$$

where $D$ is the linear growing amplitude. Mass conservation before the map folds gives

$$1+\delta(\mathbf x,D) =\frac{1}{J(\mathbf q,D)}, \qquad J=\det\left( \delta^i{}_j+D\frac{\partial s^i}{\partial q^j} \right).$$

The Zel’dovich approximation keeps the linear displacement in this nonlinear coordinate map. It is exact for plane-parallel collapse until the first crossing. At $J=0$, the density formula diverges and distinct labels can reach the same Eulerian position. Afterward, $1/J$ is not a single-valued physical density and the pressureless fluid closure has generated velocity dispersion. Bernardeau et al. derive this map, its density Jacobian, and its shell-crossing boundary in Bernardeau et al. 2002, §§2.5–2.6, printed pp. 28–32, eqs. (88)–(107) (Open PDF).

Spherical collapse as a complementary benchmark

The Zel’dovich map resolves where a cold flow first becomes multistreaming; it does not assign the resulting object a halo boundary. A different idealized benchmark follows a uniform spherical overdensity in an Einstein–de Sitter background. Birkhoff’s theorem, or its Newtonian spherical analogue, lets the patch evolve as a closed dust universe. Its radius can be written

$$r(\vartheta)=A(1-\cos\vartheta), \qquad t(\vartheta)=B(\vartheta-\sin\vartheta).$$

The patch turns around at $\vartheta=\pi$ and formally collapses at $\vartheta=2\pi$. Extrapolating the early linear growing mode to that collapse time gives the often-used Einstein–de Sitter threshold

$$\delta_c=\frac{3}{20}(12\pi)^{2/3}\simeq1.68647.$$

If energy conservation and the virial theorem are imposed after turnaround, the idealized virial radius is half the turnaround radius and the mean enclosed density is $18\pi^2\simeq177.65$ times the background density at the formal dust-collapse time. This epoch convention matters: the turnaround density ratio is $(3\pi/4)^2$, halving the radius raises the patch density by eight, and the EdS background density falls by four between turnaround and formal collapse because $t_{\rm coll}=2t_{\rm ta}$. Their product is $32(3\pi/4)^2=18\pi^2$. Virialization is an added physical prescription; the pressureless parametric solution itself collapses to a point. Tong derives the parametric solution and both numbers in Tong 2019, §§3.3.1–3.3.2, printed pp. 164–168, eqs. (3.63)–(3.71) (Official PDF).

These constants are calibration checks, not universal halo definitions. They change with the background cosmology and collapse convention, while real collapse is triaxial, environmentally coupled, and followed by mixing rather than an instantaneous top-hat virialization. The next page therefore treats mass functions as model-dependent mappings between a smoothing filter, collapse barrier, and chosen halo mass.

A resolved crossing benchmark

Use a periodic one-dimensional domain $0\le q\lt 2\pi$ and

$$x(q,D)=q-DA\sin(q-q_*), \qquad J(q,D)=1-DA\cos(q-q_*).$$

For $A=0.4$ and $q_*=2\pi/3$, the exact first crossing occurs at $q=q_*$ and $D_*=A^{-1}=2.5$. A grid estimates the largest sampled compression and hence

$$D_*^{(N)} =\frac{1}{A\max_i\cos(q_i-q_*)}.$$

The deterministic benchmark uses Node.js 24.18.0 on Darwin arm64, IEEE-754 binary64 arithmetic, and the JavaScript standard library only; it has no external dependencies. Its results are:

Grid points $N$	Estimated $D_*$	Absolute error	Minimum sampled $J(D_*)$	Observed order
32	$2.5053641770$	$5.364177\times10^{-3}$	$2.141077\times10^{-3}$	—
64	$2.5013392484$	$1.339248\times10^{-3}$	$5.354125\times10^{-4}$	$2.0019$
128	$2.5003347000$	$3.347000\times10^{-4}$	$1.338621\times10^{-4}$	$2.0005$
256	$2.5000836680$	$8.366801\times10^{-5}$	$3.346608\times10^{-5}$	$2.0001$
512	$2.5000209166$	$2.091656\times10^{-5}$	$8.366556\times10^{-6}$	$2.0000$

At $D=0.8D_*$, every grid also gives $M/\bar M=1.000000000000$ when integrating $\rho/\bar\rho=1/J$ with $dx=J\,dq$. This is an exact Jacobian identity for the analytic map, useful for catching a density or volume-factor error; it is not a demanding test of a force solver.

Download the MIT-licensed standard-library benchmark, then reproduce the calculation with the standalone download-and-run recipe:

curl -fL https://blackhole.org/code/cosmology/cosmology-benchmarks.mjs -o cosmology-benchmarks.mjs
node cosmology-benchmarks.mjs --case zeldovich

Beyond the one-dimensional map

An $N$-body calculation samples the distribution with macro-particles and approximates gravitational forces by a mesh, tree, multipole, or hybrid algorithm. A credible result must vary particle number, box size, force resolution or softening, time step, starting redshift, and initial-condition order; it must also compare mass, momentum, large-scale linear growth, and code-independent benchmarks. A single higher-resolution run is not a convergence argument.

Principal failure modes include discreteness, force anisotropy, aliasing, finite-volume sample variance, transient decaying modes, inaccurate time integration, and interpreting the Zel’dovich map after shell crossing. Baryonic cooling and feedback are explicitly excluded here, as are theory-specific gravity tests handled in dedicated volumes. Full simulations that include either require their own equations, calibration data, code version, halo or observable definition, and uncertainty budget. Horizon-scale relativistic projection effects are treated on the final page of this chapter.

Exercises

Explain the observed second-order convergence of $D_*^{(N)}$. For the grids above, show that the nearest point to $q_*$ is displaced by $\epsilon_N=2\pi/(3N)$.

Solution

Because powers of two are not divisible by three, $q_*/\Delta q=N/3$ lies one third of a grid interval from the nearest point. Hence

$$\epsilon_N=\frac{\Delta q}{3} =\frac{2\pi}{3N}.$$

Near the compression maximum,

$$\cos\epsilon_N =1-\frac{\epsilon_N^2}{2}+O(\epsilon_N^4).$$

Therefore

$$D_*^{(N)} =\frac{D_*}{\cos\epsilon_N} =D_*\left[ 1+\frac{\epsilon_N^2}{2} +O(\epsilon_N^4) \right].$$

The leading localization error is proportional to $N^{-2}$, so doubling $N$ reduces it by four and gives observed order two. This rate diagnoses sampling of a smooth maximum; it says nothing about the post-crossing dynamics.

References

• Bernardeau, F., S. Colombi, E. Gaztañaga, and R. Scoccimarro, “Large-Scale Structure of the Universe and Cosmological Perturbation Theory,” Physics Reports 367, 1–248 (2002), DOI, Open PDF.
• Tong, D., Cosmology, University of Cambridge Part II Mathematical Tripos lecture notes, Michaelmas 2019, §§3.3.1–3.3.2, printed pp. 164–168, Official PDF.