Einstein–Boltzmann Evolution

Linear cosmological observables are obtained by evolving gravity, fluids, and particle phase-space distributions as one constrained system. The Einstein equations determine the metric potentials from stress–energy moments; the Boltzmann equation transports each distribution and generates those moments. This page derives the hierarchy, explains its controlled approximations, and gives a small reproducible closure test. It does not substitute that test for a production Boltzmann code. The displayed equations treat linear scalar modes on a spatially flat FLRW background; vector and tensor hierarchies, second-order sources, and nonlinear evolution lie outside the scope of this page.

Helpful background. The Einstein equation and relativistic matter introduce the constraints and stress–energy sources. Formulations of general relativity explains gauge and constraint monitoring; the cosmology overview locates the background history; and PDEs and Green functions provides Fourier and evolution methods. The coupled gravitational equations are derived in Linearized Einstein and Conservation Equations.

From phase space to moments

Use Newtonian gauge,

$$ds^2=a^2\left[-(1+2\Phi)d\eta^2 +(1-2\Psi)\delta_{ij}dx^i dx^j\right],$$

and the Fourier convention $e^{i\mathbf k\cdot\mathbf x}$. For a species with comoving momentum magnitude $q$, direction $\mathbf n$, and $\mu=\hat{\mathbf k}\cdot\mathbf n$, write $f=f_0(q)[1+\chi(\eta,\mathbf k,q,\mu)]$. The first-order kinetic equation is

$$\frac{Df}{d\eta} =\frac{\partial f}{\partial\eta} +\frac{dx^i}{d\eta}\frac{\partial f}{\partial x^i} +\frac{dq}{d\eta}\frac{\partial f}{\partial q} =\left(\frac{\partial f}{\partial\eta}\right)_{\!C},$$

where the geodesic equation supplies the gravitational redshift terms and the right side is the collision operator. The exact phase-space derivative also contains $(dn^i/d\eta)\partial f/\partial n^i$. It does not contribute at this order because $f_0$ is isotropic and the change of direction is already first order. For a massless, collisionless species, integrate over $q^3f_0\,dq$ and define a brightness perturbation $F(\mu)$ by

$$F(\mu) =\frac{\int_0^\infty q^3f_0(q)\chi(q,\mu)\,dq} {\int_0^\infty q^3f_0(q)\,dq}.$$

In this convention it obeys

$$F'+ik\mu F=4\left(\Psi'-ik\mu\Phi\right).$$

Expand its angular dependence as

$$F(\mu)=\sum_{\ell=0}^{\infty} (-i)^\ell(2\ell+1)F_\ell P_\ell(\mu).$$

Orthogonality of the Legendre polynomials and $(2\ell+1)\mu P_\ell=(\ell+1)P_{\ell+1}+\ell P_{\ell-1}$ then turn one transport equation into neighboring moment equations. With the Ma–Bertschinger normalization,

$$\delta=F_0,\qquad \theta=\frac{3k}{4}F_1,\qquad \sigma_{\rm MB}=\frac12F_2,$$

and

$$\begin{aligned} \delta'&=-\frac43\theta+4\Psi',\\ \theta'&=k^2\left(\frac14\delta-\sigma_{\rm MB}\right)+k^2\Phi,\\ F_\ell'&=\frac{k}{2\ell+1} \left[\ell F_{\ell-1}-(\ell+1)F_{\ell+1}\right], \qquad \ell\ge2. \end{aligned}$$

Here $\sigma_{\rm MB}$ is particle anisotropic shear, not the metric $\sigma=E'-B$ used earlier in the chapter. For a massless neutrino species, $F=F_\nu$; a massive species retains the $q$ dependence instead of closing on one momentum-integrated hierarchy. Ma and Bertschinger derive the phase-space equation in Ma and Bertschinger 1995, §5.1, printed pp. 14–16, eqs. (33)–(41) (Open PDF) and the massless-neutrino normalization and hierarchy in Ma and Bertschinger 1995, §5.3, printed pp. 17–18, eqs. (44)–(51) (Open PDF).

Photons, baryons, dark matter, and gravity

Photons have the same free-streaming left side, but Thomson scattering couples their intensity and polarization multipoles to electrons. Define the positive conformal opacity

$$\dot\kappa\equiv an_e\sigma_T\gt0.$$

Here $n_e$ is the physical free-electron number density, so $an_e\sigma_T$ is a conformal-time scattering rate. Suppressing the polarization hierarchy only to display the low moments, define $F_{\gamma\ell}=4\Theta_\ell$ for a linear blackbody temperature perturbation. Then

$$\delta_\gamma=4\Theta_0, \qquad \theta_\gamma=3k\Theta_1, \qquad \sigma_{\gamma,{\rm MB}}=2\Theta_2.$$

The Newtonian-gauge density and velocity equations are

$$\begin{aligned} \delta_\gamma'&=-\frac43\theta_\gamma+4\Psi',\\ \theta_\gamma'&= k^2\left(\frac14\delta_\gamma-\sigma_{\gamma,{\rm MB}}\right) +k^2\Phi+\dot\kappa(\theta_b-\theta_\gamma),\\ \delta_b'&=-\theta_b+3\Psi',\\ \theta_b'&=-\mathcal H\theta_b+c_b^2k^2\delta_b+k^2\Phi +\frac{4\rho_\gamma}{3\rho_b}\dot\kappa (\theta_\gamma-\theta_b),\\ \delta_c'&=-\theta_c+3\Psi',\\ \theta_c'&=-\mathcal H\theta_c+k^2\Phi. \end{aligned}$$

The photon quadrupole is not simply damped at the rate $\dot\kappa$: Thomson scattering mixes it with the linear-polarization brightness multipoles $G_{\gamma\ell}$. In the same normalization,

$$\begin{aligned} F_{\gamma2}'={}&\frac{8}{15}\theta_\gamma -\frac35 kF_{\gamma3} -\frac9{10}\dot\kappa F_{\gamma2} +\frac1{10}\dot\kappa(G_{\gamma0}+G_{\gamma2}),\\ F_{\gamma\ell}'={}&\frac{k}{2\ell+1} \left[\ell F_{\gamma,\ell-1}-(\ell+1)F_{\gamma,\ell+1}\right] -\dot\kappa F_{\gamma\ell}, \qquad \ell\ge3. \end{aligned}$$

The $G_{\gamma\ell}$ hierarchy must be evolved alongside these equations; the quadrupole source shows why an intensity-only hierarchy is not a closed model of CMB polarization.

The opposite photon and baryon drag terms conserve total momentum. To see the cancellation, multiply the photon velocity equation by $4\rho_\gamma/3$ and the baryon velocity equation by $\rho_b$: their collision terms are equal and opposite. A complete calculation also evolves photon polarization, free-streaming neutrinos, massive-species momentum bins, ionization and baryon temperature, and any declared interactions. The total density, momentum, pressure, and anisotropic stress moments feed the Einstein constraints and evolution equations from the preceding page, closing the system. Initial amplitudes must satisfy those constraints and correspond to a declared adiabatic or isocurvature basis. The displayed baryon closure treats baryons and electrons as one nonrelativistic fluid, retains the scalar pressure perturbation $\delta p_b=c_b^2\delta\rho_b$, and neglects baryon anisotropic stress.

The cold-dark-matter limit is derived in Ma and Bertschinger 1995, §5.2, printed p. 16, eqs. (42)–(43) (Open PDF).

The full photon–baryon hierarchy and collision terms are given in Ma and Bertschinger 1995, §§5.5–5.7, printed pp. 21–25, eqs. (61)–(75) (Open PDF).

Tight coupling, stiffness, and closure

Before recombination, $\dot\kappa\gg\mathcal H,k$. The direct equations are stiff: an explicit or otherwise nonstiff solver must resolve the short collision time for stability, even though the desired mode evolves much more slowly. An implicit stiff integrator can avoid that stability restriction, while a tight-coupling expansion instead evolves the common photon–baryon motion and expands the slip and photon quadrupole in $k/\dot\kappa$ and $\mathcal H/\dot\kappa$. It is an asymptotic regime, not an identity: switching to the explicit hierarchy must occur where both descriptions overlap, and observables must be stable when the switch threshold and approximation order are varied.

Free streaming creates an infinite angular hierarchy. At $\ell_{\max}$, simply setting $F_{\ell_{\max}+1}=0$ reflects power back toward low multipoles. For a source-free spherical-Bessel stream, $F_\ell=j_\ell(k\eta)$, the exact recurrence suggests

$$F_{\ell_{\max}+1} =\frac{2\ell_{\max}+1}{k\eta}F_{\ell_{\max}} -F_{\ell_{\max}-1}.$$

This is far better than a hard wall in that benchmark, but time-dependent metric sources and a small $\ell_{\max}$ make it only an approximation. Production solvers combine controlled tight-coupling, fluid, streaming, and stiff-integration regimes. Blas, Lesgourgues, and Tram give the full tight-coupling equations and switch comparisons in Blas, Lesgourgues, and Tram 2011, §2.1, printed pp. 4–5, and §§2.5–2.7, pp. 9–13 (Open PDF), the ultra-relativistic-fluid and radiation-streaming approximations in Blas, Lesgourgues, and Tram 2011, §3, pp. 13–19, and §4, pp. 19–25 (Open PDF), and the stiff integrator and hierarchy derivation in Blas, Lesgourgues, and Tram 2011, Appendices A–B, pp. 28–34 (Open PDF).

For current physical context, the Particle Data Group summarizes production CMB codes in Scott and Smoot 2026, §29.5, printed p. 6 (official PDF) and the driven photon–baryon acoustic picture in Scott and Smoot 2026, §29.5.2, printed pp. 7–8 (official PDF).

A reproducible closure benchmark

The repository benchmark uses deterministic binary64 RK4 under Node.js 24.18.0 on Darwin arm64. It uses the JavaScript standard library only, with no external packages or dependencies. Its first case takes the zero-baryon-loading tight-coupling limit, sets the metric forcing to zero, and uses $c_s^2=1/3$. It therefore solves the deliberately stripped acoustic equation

$$\Theta_0''+\omega^2\Theta_0=0, \qquad \omega=\frac{k}{\sqrt3},$$

with $k=4$, $\eta_i=0$, $\eta_f=2.3$, $\Theta_0(0)=1$, and $\Theta_0'(0)=0$. Since $\Theta_0'=-k\Theta_1$ in the source-free moment normalization, the initial dipole is $\Theta_1(0)=0$. The tabulated final-state error is

$$E_\Theta=\left\{ [\Theta_0^{\rm num}-\Theta_0^{\rm exact}]^2 +\left[\frac{(\Theta_0')^{\rm num}-(\Theta_0')^{\rm exact}}{\omega}\right]^2 \right\}^{1/2}.$$

Five refinements give:

Steps	State error	Maximum invariant error	Observed order
24	$1.061438\times10^{-4}$	$3.893132\times10^{-5}$	—
48	$6.636467\times10^{-6}$	$1.222247\times10^{-6}$	$3.9995$
96	$4.148171\times10^{-7}$	$3.823916\times10^{-8}$	$3.9999$
192	$2.592666\times10^{-8}$	$1.195316\times10^{-9}$	$4.0000$
384	$1.620425\times10^{-9}$	$3.735390\times10^{-11}$	$4.0000$

The state error converges at fourth order. The independently monitored invariant is

$$I=\Theta_0^2+\left(\frac{\Theta_0'}{\omega}\right)^2 =\Theta_0^2+3\Theta_1^2,$$

and the table reports $\max_n|I_n-1|$ over accepted time steps, including the finest-grid value.

The second case evolves the source-free hierarchy from $\eta_i=4$ to $\eta_f=8$ with $k=3$, $\ell_{\max}=6$, and exact $F_\ell(\eta_i)=j_\ell(k\eta_i)$. Here the reported norm is $E_{01}=[(\Delta F_0)^2+(\Delta F_1)^2]^{1/2}$, not a full-hierarchy norm:

Steps	Hard-cutoff $E_{01}$	Recurrence $E_{01}$	Recurrence order
32	$1.698445\times10^{-1}$	$7.861282\times10^{-5}$	—
64	$1.699686\times10^{-1}$	$4.867218\times10^{-6}$	$4.0136$
128	$1.699754\times10^{-1}$	$3.030333\times10^{-7}$	$4.0055$
256	$1.699758\times10^{-1}$	$1.890771\times10^{-8}$	$4.0024$
512	$1.699758\times10^{-1}$	$1.180808\times10^{-9}$	$4.0011$

The recurrence error again approaches fourth order. The hard closure instead plateaus: its finest monopole–dipole polar amplitude is $4.66274684$ times the exact value and its polar-angle error is $0.7787306$ rad. The recurrence gives amplitude ratio $1.00000000$ and angle error $-2.772989\times10^{-8}$ rad.

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 boltzmann

These are unit tests for integration order, normalization, an invariant, and closure reflection. They contain no recombination history, metric feedback, polarization, massive particles, or observable spectra. No CLASS executable or current source release was run here. CLASS I documents the module pipeline and its input and precision files in Lesgourgues 2011, §3.4, printed pp. 13–14 (Open PDF). The cited CLASS method papers support the discussion; any future numerical comparison must additionally pin and retain the exact code commit, input and precision files, gauge and output settings, compiler environment, and executable or source hash. Any larger executable system must retain that complete record together with its convergence and truncation checks.

The principal failure modes are inconsistent gauge normalization or initial constraints, stiff drag outside tight coupling, a switch made without an overlap test, too small $\ell_{\max}$, closure reflection, phase error, and mistaking this solution-matched boundary test for evidence about a sourced production hierarchy. Recombination, metric forcing, arbitrary angular initial data, and transfer functions require independent sensitivity and convergence tests.

Exercises

Starting from the source-free transport equation $F'+ik\mu F=0$, derive the moment hierarchy and show why $F_\ell=j_\ell(k\eta)$ solves it.

Solution

Insert the Legendre series and use

$$\mu P_\ell(\mu) =\frac{\ell+1}{2\ell+1}P_{\ell+1}(\mu) +\frac{\ell}{2\ell+1}P_{\ell-1}(\mu).$$

Projecting onto $P_\ell$ gives

$$F_\ell'=\frac{k}{2\ell+1} \left[\ell F_{\ell-1}-(\ell+1)F_{\ell+1}\right].$$

For $x=k\eta$, the spherical-Bessel derivative identity

$$(2\ell+1)\frac{dj_\ell}{dx} =\ell j_{\ell-1}-(\ell+1)j_{\ell+1}$$

makes the right side equal to $k\,dj_\ell/dx=dj_\ell/d\eta$. Thus the stream solves every member of the infinite hierarchy. Truncating it with $F_{\ell_{\max}+1}=0$ violates the recurrence unless that omitted multipole happens to vanish, explaining the reflected error.

References

• Blas, D., J. Lesgourgues, and T. Tram, “The Cosmic Linear Anisotropy Solving System (CLASS) II: Approximation Schemes,” Journal of Cosmology and Astroparticle Physics 2011(07), 034 (2011), DOI, Open PDF.
• Lesgourgues, J., “The Cosmic Linear Anisotropy Solving System (CLASS) I: Overview,” arXiv:1104.2932v2 (2011), Open PDF.
• Ma, C.-P., and E. Bertschinger, “Cosmological Perturbation Theory in the Synchronous and Conformal Newtonian Gauges,” The Astrophysical Journal 455, 7–25 (1995), DOI, Open PDF.
• Scott, D., and G. F. Smoot, “Cosmic Microwave Background,” in F. Takahashi et al. (Particle Data Group), Review of Particle Physics, International Journal of Modern Physics A 41, 2630011 (2026), Official PDF.