Cosmological Perturbation Theory

Why a linear theory is enough — for most of cosmic history

The observable universe averages out, on scales of a few hundred megaparsecs, to a remarkably uniform fluid: matter density, expansion rate, and even temperature are the same in every direction to about one part in a hundred thousand. The smooth Friedmann-Lemaître-Robertson-Walker (FLRW) solution captures this background. But the structures we actually see — galaxies, clusters, the cosmic web — required that perfect uniformity to be slightly broken. Cosmological perturbation theory is the mathematical machinery for those small departures from uniformity.

The setup is to write every cosmological quantity as background plus a small bump:

ρ(x,t) = ρ̄(t) · (1 + δ(x,t))           density contrast
g_μν(x,t) = ḡ_μν(t) + h_μν(x,t)          metric perturbation
v(x,t) = 0 + v_pec(x,t)                  peculiar velocity

and then keep only the first-order terms in δ, h, v_pec when the full Einstein-fluid system is expanded around the background. The result is a set of linear equations that evolve each independent Fourier mode δ_k(t) with constant coefficients depending only on time. This linearity is the key technical point: it lets us solve mode by mode, build a power spectrum, and then project it cleanly onto the CMB sky.

How long does the linear regime last? Until structures begin to turn around and collapse, which roughly tracks δ ≈ 1 on the relevant scale. For the dark-matter density field, that happens at z ~ 0 on scales below about 10 megaparsecs; for the CMB temperature anisotropy (δT/T ~ 10⁻⁵), it has never been violated. Linear theory therefore describes the CMB to sub-percent precision, and matches galaxy surveys on large scales (k ≲ 0.1 h/Mpc). Smaller scales require N-body simulations or higher-order perturbation theory.

Fourier modes — turning a PDE into infinitely many ODEs

Because the FLRW background is translation-invariant in space, linear perturbations on it decouple in Fourier space. Decomposing the density contrast as

δ(x,t) = ∫ d³k / (2π)³ · δ_k(t) · e^(ik·x)

turns the partial differential equation for δ(x,t) into an ordinary differential equation for each δ_k(t). Different k modes evolve completely independently in linear theory — there is no mode-mode coupling, no cascade, no turbulence. Each comoving wavenumber k corresponds to a fixed comoving wavelength λ = 2π/k that stretches with the expansion. The crucial geometrical fact is the comparison of this wavelength to the Hubble radius c/H:

• Super-horizon modes (k < aH) are larger than the causal horizon at that epoch. They are essentially frozen: gauge-invariant curvature perturbations stay constant.
• Sub-horizon modes (k > aH) are smaller than the horizon. Causal physics — gravity, pressure, diffusion — can act on them, so they evolve, oscillate, grow, or decay.

Inflation makes every mode start out super-horizon (because inflation blew up the universe so fast that every comoving scale was once outside the Hubble radius), then re-enter as the universe decelerates. Modes re-enter at different times — small scales first, large scales much later — so the imprint on each is set at different cosmic epochs. This is what gives the matter power spectrum its characteristic broken-power-law shape.

SVT decomposition — scalars, vectors, tensors

On a homogeneous, isotropic background, the perturbation modes can be classified by how they transform under spatial rotations. This Stewart-Walker (or scalar-vector-tensor) decomposition splits a general metric perturbation into three completely decoupled pieces in linear theory:

Type	Spin	Physical content	Behaviour in expanding universe	Observable signature
Scalar	0	Density, potential, peculiar velocity divergence	Grow into structure	Galaxies, CMB temperature, E-mode polarisation
Vector	1	Vorticity, frame-dragging	Decay as 1/a²	None at observable level
Tensor	2	Transverse-traceless metric perturbation = gravitational waves	Free-stream after horizon crossing	B-mode polarisation, stochastic GW background

The decoupling is exact at linear order: a scalar mode never sources a vector or tensor mode, and so on. At second and higher orders the modes mix (scalar-scalar terms can source tensor perturbations, for example — the dominant source of stochastic gravitational waves in the standard cosmology). But for everything observable in the CMB or galaxy surveys, the linear SVT split is the right description.

The vector modes are the boring case: they decay rapidly with the expansion and contain almost no information about the early universe. The interesting modes are scalars (which become galaxies) and tensors (which produce B-mode polarisation and a stochastic gravitational-wave background).

The growth equation

For a single Fourier mode δ_k of cold dark matter on sub-horizon scales, linearised general relativity (or, equivalently, the Newtonian limit) gives

δ̈_k + 2H δ̇_k − (3/2) H² Ω_m δ_k = 0

The three terms are: inertia, Hubble friction, gravitational self-attraction. The solutions split into a growing mode and a decaying mode. Only the growing mode matters at late times; it scales as the linear growth factor D(a). Limiting cases are:

Epoch	Dominant component	Growth factor D(a)	Note
Radiation era (a < a_eq ≈ 1/3400)	radiation	logarithmic in a	Mészáros suppression
Matter era (a_eq ≲ a ≲ a_Λ ≈ 0.7)	matter	D ∝ a	Pure linear growth
Λ era (a ≳ 0.7)	cosmological constant	D → const	Growth saturates

So the entire history of structure growth in our universe is: a brief logarithmic suppression in the radiation era, three orders of magnitude of clean linear growth in the matter era, and a saturation today as dark energy takes over. The factor by which a mode is amplified from CMB to today is about D(today) / D(z_dec) ≈ 1100 — the CMB anisotropies of 10⁻⁵ have grown into δ ~ 10⁻² on the scales we observe in galaxy clustering.

The Mészáros effect — why sub-horizon perturbations stall during radiation domination

During radiation domination, the universe expands as a ∝ t^(1/2). The Hubble rate is dominated by radiation, while sub-horizon dark-matter perturbations have only their own self-gravity to drive them — and the radiation pressure prevents the photon-baryon fluid from clustering at all. The dark-matter growth equation in this regime reduces to

δ̈_DM + 2H δ̇_DM ≈ 0   (gravitational source ≪ Hubble friction)

which has two solutions: a constant and a logarithm. Péter Mészáros showed in 1974 that the late-time solution scales as δ ∝ 1 + (3/2)(a/a_eq) — essentially flat during radiation domination, then transitioning smoothly to linear growth after equality. The kink at k_eq ≈ 0.01 h/Mpc is one of the most precisely measured features of the matter power spectrum and a powerful constraint on Ω_m h².

The power spectrum P(k) — the central two-point statistic

If the density field is statistically homogeneous, isotropic and Gaussian — all assumptions strongly supported by observation — then its full statistical content is captured by the two-point function. In Fourier space this is the power spectrum:

⟨δ_k δ*_k'⟩ = (2π)³ δ_D(k − k') · P(k)

where δ_D is the three-dimensional Dirac delta and the angle brackets average over realisations. The same two-point information lives in the position-space correlation function ξ(r) = ⟨δ(x) δ(x+r)⟩, related to P(k) by Fourier transform. Cosmologists work mostly with P(k) because the equations of motion are diagonal in k.

The shape of today's matter power spectrum is the product of two factors:

P(k, a₀) = T²(k) · P_primordial(k)

The primordial spectrum P_primordial(k) is set by inflation and is nearly scale-invariant: P_primordial(k) ∝ k^(n_s) with the Planck-measured spectral index n_s = 0.9649 ± 0.0042. The transfer function T(k) encodes the linear evolution from end-of-inflation to today; it is approximately constant for k < k_eq, then falls as 1/k² for k > k_eq because those modes spent the radiation era in the Mészáros suppression. The product gives the canonical broken-power-law shape: ∝ k^(n_s) at large scales, turnover at k_eq, ∝ k^(n_s − 4) at small scales. Galaxy surveys (BOSS, eBOSS, DESI) measure this curve directly and pin down Ω_m h² to the few-percent level.

Projecting P(k) onto the sky — the CMB angular power spectrum C_l

The CMB last-scattering surface at z ≈ 1100 is a two-dimensional spherical slice through the perturbation field at that epoch. We see its temperature anisotropies projected onto the sky:

δT/T (n̂) = Σ_l Σ_m  a_lm  Y_lm(n̂)

with the angular power spectrum defined as the variance of the spherical-harmonic coefficients: C_l = ⟨|a_lm|²⟩. The relationship to the underlying 3D mode spectrum is

C_l = (4π) ∫ dk/k · Δ²_primordial(k) · |T_l(k)|²

where Δ²_primordial(k) = k³ P_primordial(k) / (2π²) is the dimensionless primordial spectrum and T_l(k) is the transfer function that propagates a primordial mode k into a temperature pattern at multipole l. The transfer function carries all the physics: gravitational redshift (Sachs-Wolfe at large angles), baryon-photon acoustic oscillations, photon diffusion (Silk damping) at small angles, and an integrated Sachs-Wolfe contribution from the late-time dark-energy era. Modern Einstein-Boltzmann codes such as CAMB and CLASS compute T_l(k) numerically by integrating the linearised Einstein-Boltzmann hierarchy.

Acoustic peaks — what each bump on C_l means

Before recombination, the photon-baryon plasma is a tightly coupled fluid in which gravity compresses and radiation pressure restores. The result is standing acoustic waves in every Fourier mode that was inside the horizon at that epoch. At the moment of recombination (z ≈ 1090), each mode is captured at a particular phase of its oscillation. Modes that happen to be at maximum compression appear hot in the CMB; modes at maximum rarefaction appear cold; modes at zero amplitude are silent. The result is a series of acoustic peaks in C_l at the multipoles corresponding to the angular size of half-, full-, and multiple sound-horizon wavelengths at last scattering:

Peak	Multipole l	Physical interpretation	Parameter it constrains
1st	~220	Mode just compressed once before recombination	Angular sound horizon → spatial curvature Ω_k
2nd	~540	Mode rarefied once	Baryon-to-photon ratio Ω_b h² (suppressed by baryons)
3rd	~810	Mode compressed twice	Dark matter / radiation ratio (matter loading)
Damping tail	l ≳ 1500	Silk damping from photon mean free path	N_eff, primordial helium fraction Y_He

Planck has measured these peaks to about l ≈ 2500. Their precise positions and heights fix the six base ΛCDM parameters — Ω_b h², Ω_c h², θ_MC, τ, n_s, ln(10¹⁰ A_s) — to sub-percent precision, and from these every other cosmological quantity can be derived.

Gauge invariance and the Bardeen variables

In general relativity, a metric perturbation can be partly a coordinate artefact rather than physical. If you change time-slicing — for example, from a slicing in which the universe is uniformly expanding to one in which the density is uniform — the perturbation variables look completely different. Naïve perturbation theory can therefore produce "growth modes" that are just gauge transformations.

James Bardeen's 1980 paper introduced gauge-invariant scalar perturbations Ψ and Φ — the Bardeen potentials — that have unambiguous physical meaning regardless of slicing. Schematically, Ψ acts like the Newtonian potential in the metric, and Φ like a spatial curvature perturbation. In the absence of anisotropic stress they are equal: Ψ = Φ. A direct comparison of weak-lensing measurements (which see Ψ + Φ) and galaxy clustering (which sees Ψ) tests this relation and is sensitive to modifications of gravity. The two main practical gauges in use are:

• Synchronous gauge. Time coordinate is comoving with freely-falling observers; perturbations are entirely spatial. The native gauge of CAMB. Has a residual gauge mode that can confuse super-horizon results unless treated carefully.
• Conformal-Newtonian (longitudinal) gauge. Metric perturbations are diagonal: ds² = a²[−(1+2Ψ) dτ² + (1−2Φ) dx²]. The Bardeen Ψ and Φ appear directly. Most intuitive for sub-horizon Newtonian limits.
• Comoving gauge / curvature perturbation R. The curvature perturbation on comoving slices, R, is conserved on super-horizon scales for adiabatic modes — a powerful theoretical tool that lets inflationary predictions be carried unchanged across the unknown reheating epoch.

What inflation predicts

Inflation generates the primordial spectrum by stretching vacuum quantum fluctuations of the inflaton field to super-horizon scales, where they freeze into classical curvature perturbations. The leading predictions of single-field slow-roll inflation are:

• Adiabatic. All species fluctuate together — δ_DM/(1+w_DM) = δ_γ/(1+w_γ) = …, with no extra isocurvature degree of freedom. The CMB pins isocurvature contamination below a few percent.
• Gaussian. The inflaton fluctuations are very nearly free-field, so their statistics are Gaussian. Three-point function f_NL ≲ 1 in single-field; observed to be |f_NL^local| < 5 (Planck), consistent with the simplest models and ruling out many multi-field variants.
• Nearly scale-invariant. A perfectly de Sitter expansion would give n_s = 1; the slight tilt (n_s ≈ 0.965) reflects the slow rolling of the inflaton down its potential. Observed at the 8σ level, a non-trivial confirmation of the slow-roll picture.
• Tensor modes. Same vacuum mechanism generates gravitational waves, with amplitude r = 16 ε set by the slow-roll parameter ε. Current limit r < 0.036; next-generation experiments aim for r ~ 0.001.

The Lyth bound — what r tells us about the inflaton

David Lyth derived in 1997 a relationship between r and the field range Δφ traversed by the inflaton during the observable e-folds of inflation. Schematically

Δφ / M_Pl ≳ (r / 0.01)^(1/2)

So a detection of r > 0.01 would force the inflaton to traverse more than a reduced Planck mass, generally classed as a "large-field" model. Such models are theoretically delicate: in any effective field theory with cutoff at M_Pl, the inflaton potential receives generic Planck-suppressed corrections that look like (φ/M_Pl)ⁿ — exactly the corrections that should be ignorable in a controlled low-energy theory. Large-field inflation therefore needs a UV-complete protection mechanism (axion shift symmetry, monodromy, supergravity). Conversely, small-field models (Starobinsky's R + R², Higgs-like inflation) predict r ≲ 0.005 and are perfectly compatible with current limits. The forthcoming round of B-mode experiments will discriminate.

Where the theory has been tested

• COBE (1992). First detection of CMB anisotropy at δT/T ~ 10⁻⁵, vindicating the very existence of the perturbations linear theory was designed to describe.
• WMAP (2003–2012). Resolved the first two acoustic peaks, fixed n_s < 1, gave first robust ΛCDM parameter set.
• Planck (2009–2018). Six base parameters to sub-percent precision, |f_NL^local| < 5, n_s = 0.9649 ± 0.0042, r < 0.10 by itself.
• BICEP/Keck + Planck (2021). Combined B-mode analysis pushed r < 0.036 (95 percent).
• BOSS / eBOSS / DESI (2014–present). Map the matter power spectrum out to z ~ 2 across hundreds of cubic gigaparsecs; detect the BAO ruler at sub-percent precision.
• Cosmic shear (KiDS, DES, HSC, Euclid). Weak-lensing two-point statistics give an independent handle on σ_8, currently in mild tension with the CMB-inferred value.
• Simons Observatory / LiteBIRD (2024–). Target r ~ 0.001, the discriminating threshold for many inflationary models.

Where the linear theory breaks down — and what comes next

Linear theory cleanly predicts the universe down to scales where δ ~ 1. Beyond that, three families of techniques take over:

• Higher-order perturbation theory. Expand to second, third, … order in δ. Captures non-linear corrections analytically and predicts the leading non-Gaussian signal from gravitational mode coupling. Effective field theory of large-scale structure (EFTofLSS) is the modern, renormalised version that reliably extends linear results into the weakly non-linear regime k ~ 0.2 h/Mpc.
• N-body simulations. Direct numerical integration of self-gravitating cold dark matter particles. Millennium, IllustrisTNG, AbacusSummit and similar codes follow the full non-linear collapse, formation of haloes, mergers, and the cosmic web. Hydrodynamic extensions include gas physics for galaxy formation.
• Halo model and hybrid emulators. Phenomenological compositions of linear theory at large scales with empirical fits or machine-learned emulators at small scales. The workhorse for joint analyses of cosmic shear and clustering.

Common pitfalls

• Confusing the density contrast δ with the absolute density. δ can be negative (an under-density), and is dimensionless. Voids have δ ≈ −0.8; clusters have δ > 200 inside the virial radius.
• Ignoring gauge issues on super-horizon scales. Many spurious "growth modes" of the synchronous-gauge perturbations are pure gauge artefacts; they vanish in the conformal-Newtonian frame or in gauge-invariant variables.
• Forgetting that P(k) is the variance per mode, not per volume. Discrete tracers (galaxies) measure a binned estimator whose noise is set by sample variance plus shot noise — quantitatively very different from the ideal continuous-field result.
• Treating CMB anisotropies as a 3D snapshot. They are a 2D projection of the perturbation field at z_dec onto the sky; angular features encode physical scales mapped through the comoving angular-diameter distance to last scattering.
• Misreading n_s = 1 as "no inflation". Many simple inflationary models do predict n_s very close to 1; the Planck deviation to n_s = 0.9649 is a non-trivial selection, but n_s = 1 itself is not in itself a contradiction with the inflationary paradigm.
• Quoting r without specifying the pivot. The tensor-to-scalar ratio depends on the scale at which it is evaluated. Conventional pivots are k_pivot = 0.002 Mpc⁻¹ and k_pivot = 0.05 Mpc⁻¹; the implied bounds differ slightly.