Part II: Thermal History | Chapter 3

Recombination & the Cosmic Microwave Background

How the universe became transparent: the physics of hydrogen recombination, photon decoupling, and the origin of the oldest light in the cosmos

1. Introduction: When the Universe Became Transparent

For the first 380,000 years after the Big Bang, the universe was a hot, opaque plasma of protons, electrons, and photons. The photons were tightly coupled to the baryonic matter through Thomson scattering off free electrons, maintaining thermal equilibrium. As the universe expanded and cooled, a remarkable phase transition occurred: free protons and electrons combined to form neutral hydrogen atoms, a process cosmologists call recombination (despite the fact that electrons and protons had never been combined before).

Why “Recombination”?

The term is a historical misnomer. In laboratory physics, recombination refers to electrons rejoining ions from which they were previously separated. In cosmology, the electrons and protons were never previously combined — they are combining for the first time. Despite this, the name has stuck and is universally used in the literature.

Recombination is one of the most precisely understood epochs in cosmology. The physics involves atomic hydrogen, statistical mechanics, and radiative transfer — all well-tested in terrestrial laboratories. This makes theoretical predictions for the CMB exceptionally reliable.

The key physical process is the reaction:

$$p + e^- \longleftrightarrow \text{H} + \gamma$$

When the temperature drops below the binding energy of hydrogen ($B_1 = 13.6$ eV), the equilibrium shifts to the right and neutral atoms form. However, due to the enormous photon-to-baryon ratio ($\eta \sim 10^{-9}$), this transition occurs at a temperature far below $B_1$.

Once neutral atoms form, the photon mean free path grows rapidly and eventually exceeds the Hubble radius. The photons decouple from matter and begin to free-stream through the universe. These freely-streaming photons — redshifted by a factor of about 1090 — are what we observe today as the cosmic microwave background (CMB), a nearly perfect blackbody at$T_0 = 2.7255 \pm 0.0006$ K.

In this chapter, we derive the physics of recombination from first principles. We begin with the Saha equation for equilibrium ionization, then move to the Peebles three-level atom model that captures the non-equilibrium dynamics, compute the optical depth and visibility function that define the last scattering surface, and finally discuss photon decoupling and Silk damping.

2. The Saha Equation for Hydrogen Recombination

We begin with the equilibrium treatment of recombination. Consider a gas of protons, electrons, and hydrogen atoms in thermal equilibrium at temperature $T$. We define the ionization fraction:

$$x_e \equiv \frac{n_e}{n_b} = \frac{n_p}{n_H + n_p}$$

where $n_e$ is the electron number density, $n_b = n_H + n_p$ is the total baryon (hydrogen) number density, and we ignore helium for this derivation.

2.1 Derivation from Chemical Equilibrium

In thermal equilibrium, the chemical potentials satisfy:

$$\mu_p + \mu_e = \mu_H$$

For non-relativistic species in thermal equilibrium, the number density of species $i$ with mass $m_i$, spin degeneracy $g_i$, and chemical potential $\mu_i$ is given by the Maxwell-Boltzmann distribution:

$$n_i = g_i \left(\frac{m_i k_B T}{2\pi \hbar^2}\right)^{3/2} \exp\!\left(\frac{\mu_i - m_i c^2}{k_B T}\right)$$

We apply this to protons ($g_p = 2$), electrons ($g_e = 2$), and hydrogen atoms ($g_H = 4$, counting both spin states of the electron and proton in the ground state). Taking the ratio $n_p n_e / n_H$:

$$\frac{n_p\, n_e}{n_H} = \frac{g_p\, g_e}{g_H}\left(\frac{m_p m_e}{m_H}\right)^{3/2}\!\left(\frac{k_B T}{2\pi\hbar^2}\right)^{3/2} \exp\!\left(-\frac{(m_p + m_e - m_H)c^2}{k_B T}\right)$$

The mass difference in the exponential is precisely the binding energy of hydrogen:

$$(m_p + m_e - m_H)c^2 = B_1 = 13.6\;\text{eV}$$

Since $m_H \approx m_p \gg m_e$, the prefactor $(m_p m_e / m_H)^{3/2} \approx m_e^{3/2}$, and the spin degeneracy ratio $g_p g_e / g_H = 4/4 = 1$. Using $n_p = n_e = x_e n_b$ and$n_H = (1 - x_e) n_b$, we obtain the Saha equation:

The Saha Equation

$$\boxed{\frac{x_e^2}{1 - x_e} = \frac{1}{n_b}\left(\frac{m_e\, k_B T}{2\pi\hbar^2}\right)^{3/2} \exp\!\left(-\frac{B_1}{k_B T}\right)}$$

This is a quadratic equation in $x_e$ that can be solved analytically given $n_b$ and $T$.

2.2 The Recombination Temperature

A naïve estimate sets $k_B T_{\text{rec}} \sim B_1 = 13.6$ eV, giving$T_{\text{rec}} \sim 1.6 \times 10^5$ K. This is dramatically wrong! The actual recombination temperature is $T_{\text{rec}} \approx 3000$ K, about 50 times lower. The reason lies in the enormous photon-to-baryon ratio.

The baryon number density scales as $n_b = n_{b,0}(1+z)^3$, and the CMB temperature as $T = T_0(1+z)$. The right-hand side of the Saha equation contains the factor:

$$\frac{1}{n_b}\left(\frac{m_e k_B T}{2\pi\hbar^2}\right)^{3/2} \propto \frac{T^{3/2}}{n_b} \propto \frac{(1+z)^{3/2}}{(1+z)^3} = \frac{1}{(1+z)^{3/2}}$$

More precisely, this factor can be written as:

$$\frac{1}{n_b}\left(\frac{m_e k_B T}{2\pi\hbar^2}\right)^{3/2} = \frac{1}{\eta}\,\frac{1}{n_\gamma}\left(\frac{m_e k_B T}{2\pi\hbar^2}\right)^{3/2} \sim \frac{1}{\eta}\left(\frac{m_e}{T}\right)^{3/2} \sim \frac{10^{15}}{\eta}$$

where $\eta = n_b/n_\gamma \approx 6 \times 10^{-10}$ is the baryon-to-photon ratio. The factor $1/\eta \sim 10^{9}$ hugely favors ionization. To overcome this and drive$x_e \to 0$, the exponential $\exp(-B_1/k_BT)$ must become extremely small, requiring $k_B T \ll B_1$. Setting $x_e = 0.5$ in the Saha equation and solving numerically gives:

$$T_{\text{rec}} \approx 3740\;\text{K} \quad (z_{\text{rec}} \approx 1370)$$

The Saha prediction for $x_e = 0.5$. The more careful Peebles treatment gives$T_{\text{rec}} \approx 3000$ K ($z \approx 1100$), and the Planck satellite measures the last scattering redshift as $z_* = 1089.80 \pm 0.21$.

The Photon-to-Baryon Ratio Argument

The key insight is that for every baryon, there are $\sim 10^9$ photons. Even when the average photon energy $\langle E_\gamma \rangle \sim 3 k_B T$ is far below $B_1$, the high-energy tail of the Planck distribution still contains enough photons above 13.6 eV to keep hydrogen ionized. Only when $k_B T \lesssim B_1 / \ln(1/\eta) \approx 0.3$ eV does the exponential suppression finally win.

3. The Peebles Three-Level Atom Model

The Saha equation assumes thermal equilibrium, but recombination is fundamentally a non-equilibrium process. As the universe expands, the recombination rate eventually falls below the Hubble expansion rate, and $x_e$ freezes out at a value significantly above the Saha prediction. To capture this, we need the kinetic theory of recombination.

3.1 The Problem with Direct Recombination

Consider an electron recombining directly to the ground state of hydrogen:

$$e^- + p \to \text{H}(1s) + \gamma \quad (E_\gamma \geq 13.6\;\text{eV})$$

The emitted photon has energy $E_\gamma \geq B_1 = 13.6$ eV. In a universe with$\sim 10^9$ photons per baryon and many neutral hydrogen atoms, this photon is almost immediately absorbed by a neighboring hydrogen atom, re-ionizing it:

$$\gamma + \text{H}(1s) \to e^- + p$$

Therefore, direct recombination to the ground state produces no net recombination. Every recombination is immediately undone by a photoionization. This is the fundamental bottleneck that Peebles (1968) and Zel'dovich, Kurt & Sunyaev (1968) independently recognized.

3.2 Recombination Through Excited States

The way around this bottleneck is recombination through excited states. Consider the process:

$$e^- + p \to \text{H}(n \geq 2) + \gamma$$

Recombination to excited states ($n \geq 2$) produces photons with $E_\gamma < 13.6$ eV that cannot ionize ground-state hydrogen.

The atom in the excited state then cascades down to $n = 2$. From the $n = 2$ level, there are two paths to the ground state:

Two Escape Routes from n = 2

Route 1: Lyman-$\alpha$ emission ($2p \to 1s$)

The $2p \to 1s$ transition emits a Lyman-$\alpha$ photon with energy$E_{{\rm Ly}\alpha} = B_1(1 - 1/4) = 10.2$ eV. This photon is resonantly absorbed by neighboring ground-state hydrogen atoms. The only way Ly$\alpha$ photons escape is through cosmological redshifting: the Hubble expansion redshifts them out of the line profile before they are reabsorbed.

Route 2: Two-photon decay ($2s \to 1s + 2\gamma$)

The $2s$ state is metastable (single-photon decay to $1s$ is forbidden by selection rules). It decays via simultaneous emission of two photons with a combined energy of 10.2 eV. Each photon individually has $E < 10.2$ eV, so neither can excite a Lyman-$\alpha$ transition. The two-photon decay rate is$\Lambda_{2s \to 1s} = 8.227\;\text{s}^{-1}$.

3.3 The Peebles Equation

Peebles simplified the hydrogen atom to an effective three-level system: the ground state ($1s$), the first excited state ($n = 2$), and the continuum (free electrons and protons). The rate equation for the ionization fraction is:

The Peebles Recombination Equation

$$\boxed{\frac{dx_e}{dt} = \left[\beta\,(1 - x_e) - \alpha^{(2)}\, x_e^2\, n_H\right]\, C}$$

where $\alpha^{(2)}$ is the case-B recombination coefficient, $\beta$ is the photoionization rate from $n = 2$, and $C$ is the Peebles correction factor.

Let us define each quantity in this equation:

Case-B Recombination Coefficient

$$\alpha^{(2)} = \sum_{n=2}^{\infty} \alpha_n$$

The sum of recombination coefficients to all excited states ($n \geq 2$), excluding direct recombination to the ground state. A useful fitting formula (Pequignot et al. 1991):

$$\alpha^{(2)} \approx \frac{F \times 10^{-13}\;\text{cm}^3\text{s}^{-1}}{t^{0.6166}}\,\frac{1}{1 + 0.6703\, t^{0.5300}}$$

where $t = T / 10^4\;\text{K}$ and $F = 4.309$.

Photoionization Rate from n = 2

The photoionization rate $\beta$ is related to $\alpha^{(2)}$ by detailed balance:

$$\beta = \alpha^{(2)}\left(\frac{m_e\, k_B T}{2\pi\hbar^2}\right)^{3/2} \exp\!\left(-\frac{B_2}{k_B T}\right)$$

where $B_2 = B_1/4 = 3.4$ eV is the binding energy of the $n = 2$ level.

3.4 The Peebles C Factor

The crucial ingredient is the Peebles correction factor $C$, which accounts for the bottleneck in reaching the ground state:

The Peebles C Factor

$$\boxed{C = \frac{\Lambda_{2s \to 1s} + \Lambda_\alpha}{\Lambda_{2s \to 1s} + \Lambda_\alpha + \beta}}$$

Here $\Lambda_{2s \to 1s} = 8.227\;\text{s}^{-1}$ is the two-photon decay rate from $2s$, and $\Lambda_\alpha$ is the rate at which Lyman-$\alpha$ photons escape by cosmological redshifting:

$$\Lambda_\alpha = \frac{H\,(3 B_1)^3}{(8\pi)^2\, n_{1s}\,(\hbar c)^3}$$

where $n_{1s} = (1 - x_e) n_H$ is the number density of ground-state hydrogen atoms, and$H$ is the Hubble parameter. The factor $(3B_1)^3$ comes from the Lyman-$\alpha$ photon energy ($E_{{\rm Ly}\alpha} = \frac{3}{4}B_1$).

The physical meaning of $C$ is transparent:

The numerator $\Lambda_{2s} + \Lambda_\alpha$ is the total rate at which atoms successfully transition from $n = 2$ to the ground state (by two-photon decay or Ly$\alpha$ escape).
The denominator adds $\beta$, the rate at which $n = 2$ atoms are photoionized back to the continuum.
Thus $C$ is the probability that an atom reaching $n = 2$ actually makes it to the ground state rather than being re-ionized.
At early times, $\beta \gg \Lambda_{2s} + \Lambda_\alpha$ and $C \ll 1$, suppressing recombination. At late times, $\beta \to 0$ and $C \to 1$.

Converting to Redshift

Since redshift is a more convenient variable, we convert using$dt = -dz / [(1+z)H(z)]$:

$$\frac{dx_e}{dz} = \frac{1}{(1+z)\,H(z)}\left[\beta\,(1 - x_e) - \alpha^{(2)}\, x_e^2\, n_H\right]\, C$$

4. Optical Depth and the Last Scattering Surface

As neutral hydrogen forms and the free electron density drops, the universe becomes transparent to photons. We quantify this transition using the optical depth and the visibility function.

4.1 Thomson Optical Depth

A photon traveling through the universe encounters free electrons. The probability of scattering is governed by the Thomson cross section $\sigma_T = 6.652 \times 10^{-25}\;\text{cm}^2$. The optical depth from the observer (today, $z = 0$) looking back to redshift $z$ is:

Thomson Optical Depth

$$\boxed{\tau(z) = \int_0^z \frac{n_e(z')\,\sigma_T\,c}{H(z')\,(1 + z')}\,dz'}$$

This integral accumulates the scattering probability along the photon's path from $z = 0$ to redshift $z$.

The derivation proceeds from the definition of optical depth along a path element $dl$:

$$d\tau = n_e\,\sigma_T\,dl = n_e\,\sigma_T\,c\,dt$$

Converting from $dt$ to $dz$ using $dt = -dz/[(1+z)H]$ and integrating:

$$\tau(z) = \int_0^z n_e(z')\,\sigma_T\,c\,\frac{dz'}{H(z')\,(1+z')}$$

The free electron density depends on the ionization fraction and the hydrogen density:

$$n_e(z) = x_e(z)\,n_H(z) = x_e(z)\,n_{H,0}\,(1+z)^3$$

4.2 The Visibility Function

The visibility function $g(z)$ gives the probability density that a CMB photon last scattered at redshift $z$:

The Visibility Function

$$\boxed{g(z) = -\frac{d\tau}{dz}\,e^{-\tau(z)} = \frac{n_e(z)\,\sigma_T\,c}{H(z)\,(1+z)}\,e^{-\tau(z)}}$$

The visibility function is normalized: $\int_0^\infty g(z)\,dz = 1$.

The physical interpretation is elegant. The factor $-d\tau/dz$ is the scattering rate (probability of scattering per unit redshift), and $e^{-\tau}$ is the probability that the photon has not scattered between redshift $z$ and today. Their product gives the probability that the photon's last scattering occurred at $z$.

4.3 Properties of the Last Scattering Surface

The visibility function $g(z)$ is sharply peaked, resembling a Gaussian. Numerical computation reveals:

Quantity	Value	Description
$z_*$	$1089.80 \pm 0.21$	Peak of visibility function (Planck 2018)
$T_*$	$\approx 2970$ K	Temperature at last scattering
$\Delta z$	$\approx 80$	Width of last scattering surface
$t_*$	$\approx 370{,}000$ yr	Cosmic time at last scattering
$\tau(z_*)$	$\approx 1$	Optical depth at last scattering (by definition)

A Surface or a Shell?

The “last scattering surface” is not an infinitely thin surface but a shell of finite thickness $\Delta z \approx 80$, corresponding to a comoving depth of about 15 Mpc. The finite thickness has important observational consequences: it causes damping of CMB anisotropies on small angular scales (Silk damping) and means that the CMB polarization signal probes the velocity field at slightly different epochs across the thickness of the shell.

5. Photon Decoupling and Free Streaming

After recombination, the dramatic drop in free electron density causes the photon mean free path to grow rapidly. When it exceeds the Hubble radius, photons effectively decouple from matter and begin free streaming through the universe.

5.1 The Photon Mean Free Path

The mean free path of a photon scattering off free electrons is:

Photon Mean Free Path

$$\boxed{\lambda_{\text{mfp}} = \frac{1}{n_e\,\sigma_T} = \frac{1}{x_e\, n_H\,\sigma_T}}$$

Before recombination, when $x_e \approx 1$:

$$\lambda_{\text{mfp}} = \frac{1}{n_{H,0}\,(1+z)^3\,\sigma_T} \approx 2.5 \times 10^{19}\;\text{cm}\;\left(\frac{1100}{1+z}\right)^3$$

At $z = 1100$, $\lambda_{\text{mfp}} \approx 2.5 \times 10^{19}$ cm $\approx 8$ kpc (proper).

After recombination, as $x_e \to 0$, the mean free path diverges:$\lambda_{\text{mfp}} \to \infty$. Photons effectively travel without scattering — they are free-streaming. The photons we observe in the CMB have traveled freely since last scattering at $z_* \approx 1090$.

5.2 The Decoupling Condition

Photon decoupling occurs when the scattering rate drops below the expansion rate:

$$\Gamma_{\text{Thomson}} = n_e\,\sigma_T\,c = x_e\, n_H\,\sigma_T\,c < H(z)$$

Equivalently, the photon mean free path exceeds the Hubble radius: $\lambda_{\text{mfp}} > c/H$.

Note that decoupling and recombination are not exactly simultaneous. Recombination (defined as $x_e = 0.5$) occurs at $z_{\text{rec}} \approx 1270$ (Saha) or$\sim 1100$ (Peebles), while photon decoupling (the peak of the visibility function) occurs at $z_* \approx 1090$. The universe doesn't need to be fully neutral for photons to decouple; even a small residual ionization fraction is too low to keep photons coupled.

5.3 Silk Damping

During the epoch when photons are still partially coupled to baryons, they undergo a random walk with mean free path $\lambda_{\text{mfp}}$. This diffusion erases perturbations on small scales. The Silk damping scaleis the characteristic distance photons diffuse during the recombination epoch:

Silk Damping Scale

$$\lambda_S^2 = \int_0^{t_*} \frac{c^2}{n_e\,\sigma_T\,a^2}\,\frac{1}{6}\,\frac{R^2 + \frac{16}{15}(1+R)}{(1+R)^2}\,dt$$

where $R = 3\rho_b / (4\rho_\gamma)$ is the baryon-to-photon momentum density ratio. The factor of $1/6$ comes from the 3D random walk.

An approximate expression for the Silk damping scale is:

$$\lambda_S \sim \sqrt{\frac{c}{n_e\,\sigma_T\,H}} \sim \sqrt{\frac{\lambda_{\text{mfp}}\,c}{H}}$$

This is the geometric mean of the mean free path and the Hubble radius. Numerically,$\lambda_S \sim 3$ Mpc (comoving) at $z_*$, corresponding to an angular scale of about $\ell \sim 1500$ in the CMB power spectrum.

Silk damping produces an exponential suppression of the CMB power spectrum at high multipoles:

$$C_\ell \propto e^{-2(\ell/\ell_S)^2}$$

where $\ell_S \sim 1500$ is the Silk damping multipole. This exponential damping envelope is clearly observed in the measured CMB power spectrum and provides an independent constraint on$\Omega_b h^2$ and $\Omega_m h^2$.

6. Applications: The CMB and Beyond

6.1 The CMB Blackbody Spectrum (COBE/FIRAS)

The tight coupling between photons and baryons before recombination ensures that the CMB has an almost perfect blackbody spectrum. The COBE/FIRAS instrument (Mather et al. 1994) measured the CMB spectrum with extraordinary precision, finding:

$$T_0 = 2.7255 \pm 0.0006\;\text{K}$$

with deviations from a perfect Planck spectrum constrained to be less than 50 parts per million. This is the most perfect blackbody ever measured in nature.

The blackbody spectrum is preserved during the expansion of the universe because the Planck distribution is invariant under the cosmological redshift: a blackbody at temperature $T$redshifted by a factor $(1 + z)$ remains a perfect blackbody at temperature$T/(1+z)$. This can be shown explicitly: the photon occupation number$n(\nu) = 1/[\exp(h\nu/k_BT) - 1]$ depends only on the ratio $\nu/T$, and both $\nu$ and $T$ scale as $(1+z)$.

6.2 Spectral Distortions

While the CMB is remarkably close to a blackbody, small departures are expected from energy injection processes in the early universe. These spectral distortionsfall into two categories depending on when the energy was injected:

$\mu$-type distortions ($z \gtrsim 5 \times 10^4$)

At high redshifts, Compton scattering is efficient at redistributing photon energies but cannot create or destroy photons (number-changing processes like bremsstrahlung and double Compton scattering have frozen out). The spectrum relaxes to a Bose-Einstein distribution with a chemical potential $\mu$:

$$n(\nu) = \frac{1}{e^{h\nu/(k_BT) + \mu} - 1}$$

$y$-type distortions ($z \lesssim 5 \times 10^4$)

At lower redshifts, even Compton scattering is inefficient. Energy injection produces a$y$-type distortion characterized by the Compton $y$-parameter. This is the same physics as the thermal Sunyaev-Zel'dovich effect in galaxy clusters:

$$y = \int \frac{k_B(T_e - T_\gamma)}{m_e c^2}\,n_e\,\sigma_T\,c\,dt$$

The COBE/FIRAS upper limits are $|\mu| < 9 \times 10^{-5}$ and$|y| < 1.5 \times 10^{-5}$ at 95% confidence. Future experiments like PIXIE and PRISTINE aim to improve these limits by factors of 100–1000, potentially detecting the$\mu$-distortion from Silk damping of primordial perturbations, which is predicted to be $\mu \sim 2 \times 10^{-8}$.

6.3 21-cm Cosmology

After recombination, neutral hydrogen fills the universe. The hyperfine transition of hydrogen (the spin-flip transition between the $F = 1$ and $F = 0$ states of the $1s$ level) produces radiation at a rest wavelength of 21 cm ($\nu_0 = 1420$ MHz). This provides a unique probe of the cosmic dark ages and the epoch of reionization.

21-cm Brightness Temperature

$$\delta T_b \approx 27\, x_{\text{HI}}\,(1 + \delta_b)\,\left(\frac{1 + z}{10}\right)^{1/2}\!\left(1 - \frac{T_\gamma}{T_S}\right)\;\text{mK}$$

where $x_{\text{HI}}$ is the neutral fraction, $\delta_b$ is the baryon overdensity, and $T_S$ is the spin temperature that characterizes the relative populations of the hyperfine levels. The 21-cm signal is observed in absorption when $T_S < T_\gamma$and in emission when $T_S > T_\gamma$.

Experiments like EDGES, SARAS, HERA, and the SKA aim to detect this signal from the cosmic dark ages ($30 \lesssim z \lesssim 200$) and the epoch of reionization ($6 \lesssim z \lesssim 30$), opening a new window on the universe between recombination and the formation of the first galaxies.

7. Historical Context

The Prediction and Discovery of the CMB

1948 — Gamow, Alpher & Herman

George Gamow proposed that the early universe was hot and dense enough for nuclear reactions. His students Ralph Alpher and Robert Herman predicted that the relic radiation from this hot early phase should still be present today as a thermal background with a temperature of about 5 K. This prediction was remarkably close to the actual value of 2.725 K but was largely forgotten by the community.

1964–1965 — Dicke, Peebles, Roll & Wilkinson

At Princeton, Robert Dicke and P.J.E. Peebles independently re-derived the prediction of a thermal background radiation. David Roll and David Wilkinson began building a radiometer to detect it.

1965 — Penzias & Wilson

Arno Penzias and Robert Wilson at Bell Labs, while calibrating a horn antenna for satellite communication, discovered an excess noise corresponding to a temperature of about 3.5 K that was isotropic and persistent regardless of the antenna's orientation. After learning of the Princeton group's work, they realized they had detected the CMB. This discovery earned them the 1978 Nobel Prize in Physics.

1968 — Peebles; Zel'dovich, Kurt & Sunyaev

P.J.E. Peebles in the United States and Ya.B. Zel'dovich, V.G. Kurt, and R.A. Sunyaev in the Soviet Union independently developed the theory of hydrogen recombination, recognizing the bottleneck of Lyman-$\alpha$ trapping and the importance of the two-photon$2s \to 1s$ decay. Peebles' three-level atom model remains the foundation of modern recombination calculations.

1992 — COBE

The COBE satellite's FIRAS instrument confirmed the perfect blackbody spectrum of the CMB, and the DMR instrument detected anisotropies at the level of $\Delta T/T \sim 10^{-5}$. John Mather and George Smoot received the 2006 Nobel Prize for this work.

1999–2000 — Seager, Sasselov & Scott

Sally Seager, Dimitar Sasselov, and Douglas Scott developed RECFAST, a modern multi-level atom recombination code that tracks hydrogen and helium recombination with percent-level accuracy. This code became the standard for CMB data analysis.

2003–2018 — WMAP & Planck

The WMAP and Planck satellites measured the CMB with exquisite precision, determining cosmological parameters to percent-level accuracy. Planck's measurement of $z_*$ to 0.02% precision is one of the most accurately determined quantities in all of cosmology.

8. Numerical Simulation: Recombination History

We now solve the Saha equation numerically and compute the ionization fraction $x_e(z)$and the visibility function $g(z)$. The Saha equation provides the equilibrium approximation, which is accurate at early times but deviates from the full Peebles solution at later times due to the freeze-out of recombination.

The simulation below solves the Saha equation across a range of redshifts, computes the Thomson optical depth by numerical integration, and evaluates the visibility function. All calculations use numpy only — no scipy required.

Python

recombination_saha.py197 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ============================
# Physical Constants (CGS)
# ============================
c       = 2.99792458e10      # speed of light [cm/s]
k_B     = 1.380649e-16       # Boltzmann constant [erg/K]
hbar    = 1.054571817e-27    # reduced Planck constant [erg s]
h_P     = 2.0 * np.pi * hbar  # Planck constant [erg s]
m_e     = 9.1093837015e-28   # electron mass [g]
m_p     = 1.67262192e-24     # proton mass [g]
sigma_T = 6.6524587e-25      # Thomson cross section [cm^2]
G       = 6.67430e-8         # gravitational constant [cm^3 g^-1 s^-2]
eV      = 1.602176634e-12    # 1 eV in erg

# ============================
# Cosmological Parameters (Planck 2018)
# ============================
T_cmb0   = 2.7255            # CMB temperature today [K]
H0_km    = 67.36             # Hubble constant [km/s/Mpc]
Mpc      = 3.0857e24         # 1 Mpc in cm
H0       = H0_km * 1e5 / Mpc  # H0 in s^-1
Omega_m  = 0.3153            # matter density parameter
Omega_r  = 9.14e-5           # radiation density parameter
Omega_L  = 1.0 - Omega_m - Omega_r  # dark energy (flat universe)
Omega_b  = 0.0493            # baryon density parameter
Y_p      = 0.245             # primordial helium mass fraction
n_H0     = (1.0 - Y_p) * 3.0 * H0**2 * Omega_b / (8.0 * np.pi * G * m_p)

# Hydrogen binding energy
B_H  = 13.6 * eV             # ground state [erg]
B_2  = B_H / 4.0             # n=2 binding energy [erg]

# ============================
# Helper Functions
# ============================
def T_cmb(z):
    """CMB temperature at redshift z [K]."""
    return T_cmb0 * (1.0 + z)

def n_H(z):
    """Hydrogen number density at redshift z [cm^-3]."""
    return n_H0 * (1.0 + z)**3

def Hz(z):
    """Hubble parameter at redshift z [s^-1]."""
    return H0 * np.sqrt(Omega_r * (1.0 + z)**4
                        + Omega_m * (1.0 + z)**3
                        + Omega_L)

# ============================
# Saha Equation
# ============================
def saha_xe(z):
    """Solve the Saha equation for the ionization fraction at redshift z."""
    T = T_cmb(z)
    T_erg = T * k_B
    nb = n_H(z)

# RHS of Saha equation: x_e^2/(1-x_e) = S
    S = (1.0 / nb) * (m_e * T_erg / (2.0 * np.pi * hbar**2))**1.5 \
        * np.exp(-B_H / T_erg)

# Solve quadratic: x_e^2 + S*x_e - S = 0
    # x_e = [-S + sqrt(S^2 + 4S)] / 2
    disc = S**2 + 4.0 * S
    xe = (-S + np.sqrt(max(disc, 0.0))) / 2.0
    return min(max(xe, 0.0), 1.0)

# ============================
# Compute Ionization History
# ============================
z_start = 2000
z_end   = 200
N_pts   = 5000
z_arr   = np.linspace(z_start, z_end, N_pts)

# Saha ionization fraction
xe_saha = np.array([saha_xe(z) for z in z_arr])

# ============================
# Optical Depth (numerical integration, trapezoid rule)
# ============================
# tau(z) = integral from 0 to z of n_e * sigma_T * c / [H(z') * (1+z')] dz'
# We integrate from low z to high z (right to left in our array since z_arr is descending...
# actually z_arr goes from 2000 to 200, so we integrate from index -1 to 0)

# Build integrand: n_e * sigma_T * c / [H * (1+z)]
ne_arr = xe_saha * n_H0 * (1.0 + z_arr)**3
H_arr  = np.array([Hz(z) for z in z_arr])
integrand = ne_arr * sigma_T * c / (H_arr * (1.0 + z_arr))

# tau(z) accumulated from z=0 (we approximate tau from z_end to z for each point)
# Since our z range doesn't go to 0, tau values are relative
# We integrate from the lowest z upward
tau_arr = np.zeros(N_pts)
# z_arr goes from high z to low z, so reverse for integration from low z to high z
for i in range(N_pts - 2, -1, -1):
    dz = z_arr[i] - z_arr[i + 1]  # positive since z_arr[i] > z_arr[i+1]
    tau_arr[i] = tau_arr[i + 1] + 0.5 * (integrand[i] + integrand[i + 1]) * dz

# ============================
# Visibility Function
# ============================
# g(z) = |dtau/dz| * exp(-tau)
dtau_dz = np.gradient(tau_arr, z_arr)
g_arr   = np.abs(dtau_dz) * np.exp(-tau_arr)

# Normalize
norm = np.trapezoid(g_arr, z_arr)
if abs(norm) > 0:
    g_arr = g_arr / abs(norm)

# ============================
# Find Key Quantities
# ============================
# Recombination redshift (x_e = 0.5)
idx_rec = np.argmin(np.abs(xe_saha - 0.5))
z_rec   = z_arr[idx_rec]
T_rec   = T_cmb(z_rec)

# Visibility function peak
idx_peak = np.argmax(g_arr)
z_star   = z_arr[idx_peak]
T_star   = T_cmb(z_star)

# Width of visibility function (FWHM)
g_half = g_arr[idx_peak] / 2.0
above_half = np.where(g_arr > g_half)[0]
if len(above_half) > 1:
    dz_width = abs(z_arr[above_half[0]] - z_arr[above_half[-1]])
else:
    dz_width = 0.0

# ============================
# Print Results
# ============================
print("=" * 60)
print("  RECOMBINATION HISTORY (Saha Equation)")
print("=" * 60)
print(f"  Saha recombination (x_e=0.5):  z_rec = {z_rec:.0f}")
print(f"  Recombination temperature:     T_rec = {T_rec:.0f} K")
print(f"  Visibility function peak:      z_*   = {z_star:.0f}")
print(f"  Temperature at last scatter:   T_*   = {T_star:.0f} K")
print(f"  Visibility FWHM:               dz    = {dz_width:.0f}")
print(f"  Planck 2018 value:             z_*   ~ 1090")
print(f"  Residual x_e at z=200:         x_e   = {xe_saha[-1]:.2e}")
print("=" * 60)
print()
print("  Note: The Saha equation overestimates the recombination")
print("  redshift because it assumes thermal equilibrium. The full")
print("  Peebles treatment gives z_rec ~ 1100 and a larger residual")
print("  ionization fraction (freeze-out at x_e ~ 10^-3).")
print("=" * 60)

# ============================
# Plot Results
# ============================
fig, axes = plt.subplots(3, 1, figsize=(8, 12), sharex=True)
fig.suptitle('Hydrogen Recombination History (Saha Equation)', fontsize=15, y=0.98)

# Panel 1: Ionization fraction
axes[0].semilogy(z_arr, np.maximum(xe_saha, 1e-30), 'b-', lw=2, label='Saha equation')
axes[0].axhline(0.5, color='gray', ls='--', alpha=0.5, label='$x_e = 0.5$')
axes[0].axvline(z_rec, color='r', ls=':', alpha=0.6, label=f'$z_{{rec}} = {z_rec:.0f}$')
axes[0].set_ylabel('Ionization fraction $x_e$', fontsize=13)
axes[0].set_ylim(1e-15, 2)
axes[0].legend(fontsize=11, loc='lower left')
axes[0].set_title('Ionization Fraction $x_e(z)$', fontsize=13)
axes[0].invert_xaxis()
axes[0].grid(True, alpha=0.2)

# Panel 2: Optical depth
axes[1].semilogy(z_arr, np.maximum(tau_arr, 1e-30), 'g-', lw=2)
axes[1].axhline(1.0, color='gray', ls='--', alpha=0.5, label='$\\tau = 1$')
axes[1].set_ylabel('Optical depth $\\tau(z)$', fontsize=13)
axes[1].legend(fontsize=11)
axes[1].set_title('Thomson Optical Depth $\\tau(z)$', fontsize=13)
axes[1].grid(True, alpha=0.2)

# Panel 3: Visibility function
axes[2].plot(z_arr, g_arr, 'm-', lw=2)
axes[2].axvline(z_star, color='r', ls=':', alpha=0.6, label=f'$z_* = {z_star:.0f}$')
axes[2].fill_between(z_arr, g_arr, alpha=0.15, color='purple')
axes[2].set_xlabel('Redshift $z$', fontsize=13)
axes[2].set_ylabel('Visibility $g(z)$', fontsize=13)
axes[2].set_title('Visibility Function $g(z)$', fontsize=13)
axes[2].legend(fontsize=11)
axes[2].grid(True, alpha=0.2)

plt.tight_layout()
plt.savefig('recombination_saha.png', dpi=150, bbox_inches='tight')
print("  Plot saved to recombination_saha.png")
print("=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

9. Summary

Key Results of Recombination Physics

1. The Saha equation describes equilibrium recombination: $x_e^2/(1-x_e) = (1/n_b)(m_e k_B T/2\pi\hbar^2)^{3/2} e^{-B_1/k_BT}$. The high photon-to-baryon ratio delays recombination to $T \sim 3000$ K, far below $B_1 = 13.6$ eV.

2. The Peebles three-level atom model captures the non-equilibrium dynamics. Direct recombination to the ground state is ineffective due to Lyman-$\alpha$ trapping. Net recombination proceeds via the $2s \to 1s$ two-photon decay and cosmological redshifting of Ly$\alpha$ photons.

3. The visibility function $g(z) = |d\tau/dz|\,e^{-\tau}$ peaks sharply at $z_* \approx 1090$ with width $\Delta z \approx 80$, defining the last scattering surface.

4. After recombination, the photon mean free path $\lambda_{\text{mfp}} = 1/(n_e \sigma_T) \to \infty$ as $x_e \to 0$, and photons free-stream. Diffusion before full decoupling causes Silk damping of small-scale CMB anisotropies.

5. The CMB preserves a near-perfect blackbody spectrum ($T_0 = 2.7255$ K, confirmed by COBE/FIRAS). Spectral distortions ($\mu$-type and $y$-type) encode information about energy injection in the early universe.

6. The 21-cm hyperfine transition of neutral hydrogen offers a probe of the post-recombination universe — the dark ages, cosmic dawn, and reionization.

Bibliography

Textbooks & Monographs

Dodelson, S. & Schmidt, F. (2020). Modern Cosmology, 2nd ed. Academic Press. — Detailed derivation of the Saha equation, Peebles equation, and visibility function (Chapter 4).
Weinberg, S. (2008). Cosmology. Oxford University Press. — Rigorous treatment of recombination and the CMB in Chapters 2 and 7.
Mukhanov, V. (2005). Physical Foundations of Cosmology. Cambridge University Press. — Comprehensive treatment of recombination physics and spectral distortions.
Kolb, E.W. & Turner, M.S. (1990). The Early Universe. Addison-Wesley. — Classic reference for the thermal history and particle cosmology of the early universe.
Ryden, B. (2017). Introduction to Cosmology, 2nd ed. Cambridge University Press. — Accessible introduction to recombination and the CMB.

Key Papers

Gamow, G. (1948). “The Evolution of the Universe,” Nature 162, 680–682.
Alpher, R.A. & Herman, R. (1948). “Evolution of the Universe,” Nature 162, 774–775. — Prediction of relic radiation at ~5 K.
Penzias, A.A. & Wilson, R.W. (1965). “A Measurement of Excess Antenna Temperature at 4080 Mc/s,” Astrophysical Journal 142, 419–421.
Peebles, P.J.E. (1968). “Recombination of the Primeval Plasma,” Astrophysical Journal 153, 1–11. — The three-level atom model.
Zel'dovich, Ya.B., Kurt, V.G. & Sunyaev, R.A. (1968). “Recombination of Hydrogen in the Hot Model of the Universe,” Soviet Physics JETP 28, 146–150.
Seager, S., Sasselov, D.D. & Scott, D. (1999). “A New Calculation of the Recombination Epoch,” Astrophysical Journal Letters 523, L1–L5. — The RECFAST code.
Mather, J.C. et al. (1994). “Measurement of the Cosmic Microwave Background Spectrum by the COBE FIRAS Instrument,” Astrophysical Journal 420, 439–444.
Planck Collaboration (2020). “Planck 2018 results. VI. Cosmological parameters,” Astronomy & Astrophysics 641, A6. arXiv:1807.06209.

← Thermal History Cosmic Inflation →

Quantity	Value	Description
\(z_*\)	\(1089.80 \pm 0.21\)	Peak of visibility function (Planck 2018)
\(T_*\)	\(\approx 2970\) K	Temperature at last scattering
\(\Delta z\)	\(\approx 80\)	Width of last scattering surface
\(t_*\)	\(\approx 370{,}000\) yr	Cosmic time at last scattering
\(\tau(z_*)\)	\(\approx 1\)	Optical depth at last scattering (by definition)