Part II: Thermal History | Chapter 1

Big Bang & Thermal History

From the Planck epoch through nucleosynthesis and recombination: a rigorous treatment of the thermal evolution of the early universe and the origin of light elements

Historical Context

The idea that the universe began in a hot, dense state was first proposed by Georges Lemaître in 1931 as the “primeval atom” hypothesis. George Gamow, Ralph Alpher, and Robert Herman developed the hot Big Bang model in the late 1940s, predicting both the synthesis of light elements and a residual thermal radiation field. The detection of the cosmic microwave background (CMB) by Penzias and Wilson in 1965 confirmed these predictions and established the thermal history of the universe as a cornerstone of modern cosmology.

Today, Big Bang Nucleosynthesis (BBN) and the CMB provide two independent precision tests of the standard cosmological model, constraining the baryon density, the number of neutrino species, and physics beyond the Standard Model down to times as early as $t \sim 1$ second after the Big Bang.

1. Introduction: The Thermal Timeline

The thermal history of the universe describes the sequence of phase transitions, particle freeze-outs, and nucleosynthetic processes that occurred as the universe expanded and cooled from an initial state of extraordinary temperature and density. In the radiation-dominated era, the Friedmann equation connects the Hubble rate directly to temperature:

Friedmann Equation in Terms of Temperature

$$H^2 = \frac{8\pi G}{3}\,\rho_{\text{rad}} = \frac{8\pi G}{3}\cdot\frac{\pi^2}{30}\,g_*(T)\,T^4$$

where we use natural units $\hbar = c = k_B = 1$ and $g_*(T)$ counts the effective number of relativistic degrees of freedom.

This gives the time-temperature relation in the radiation era:

$$t \approx 0.301\,\frac{M_{\text{Pl}}}{\sqrt{g_*}\,T^2} \approx 2.42\,\frac{1}{\sqrt{g_*}}\left(\frac{1\;\text{MeV}}{T}\right)^2\;\text{s}$$

where $M_{\text{Pl}} = (8\pi G)^{-1/2} = 2.435 \times 10^{18}$ GeV is the reduced Planck mass.

Effective Relativistic Degrees of Freedom

The total energy density of a relativistic gas at temperature T is determined by the number of spin states weighted by species type:

$$g_*(T) = \sum_{\text{bosons}} g_i\!\left(\frac{T_i}{T}\right)^{\!4} + \frac{7}{8}\sum_{\text{fermions}} g_i\!\left(\frac{T_i}{T}\right)^{\!4}$$

The factor 7/8 arises from the Fermi-Dirac vs. Bose-Einstein statistics. At $T \gg 100$ GeV with the full Standard Model particle content, $g_* = 106.75$. As the universe cools and species become non-relativistic ($T \ll m_i$), they decouple from the thermal bath and$g_*$ decreases stepwise.

2. The Planck Epoch ($t < 10^{-43}$ s)

The very earliest moments of the universe lie beyond the reach of known physics. When the cosmic temperature approaches the Planck scale, quantum gravitational effects become dominant and general relativity breaks down as an effective theory.

Planck Units

$$T_{\text{Pl}} = \sqrt{\frac{\hbar c^5}{G k_B^2}} \approx 1.42 \times 10^{32}\;\text{K} \approx 1.22 \times 10^{19}\;\text{GeV}$$

$$t_{\text{Pl}} = \sqrt{\frac{\hbar G}{c^5}} \approx 5.39 \times 10^{-44}\;\text{s}$$

$$\ell_{\text{Pl}} = \sqrt{\frac{\hbar G}{c^3}} \approx 1.62 \times 10^{-35}\;\text{m}$$

$$\rho_{\text{Pl}} = \frac{c^5}{\hbar G^2} \approx 5.16 \times 10^{93}\;\text{g/cm}^3$$

At the Planck density, the Schwarzschild radius of a region equals its Compton wavelength. The de Broglie wavelength of typical particles is comparable to $\ell_{\text{Pl}}$, meaning that spacetime itself undergoes quantum fluctuations. A complete description requires a theory of quantum gravity — whether string theory, loop quantum gravity, or some as-yet-unknown framework. The classical Friedmann equations cannot be trusted in this regime.

Key criterion: Quantum gravity effects become important when the curvature scale $R \sim H^2 \sim \ell_{\text{Pl}}^{-2}$. Below the Planck temperature, gravitational interactions are weak enough that a semiclassical description of spacetime is valid, and we can treat matter quantum mechanically on a classical curved background.

3. Grand Unification and Symmetry Breaking

3.1 The GUT Scale

Grand Unified Theories (GUTs) hypothesize that the strong, weak, and electromagnetic interactions are unified into a single gauge group at energies above the GUT scale:

$$E_{\text{GUT}} \sim 10^{16}\;\text{GeV}, \qquad T_{\text{GUT}} \sim 10^{29}\;\text{K}, \qquad t_{\text{GUT}} \sim 10^{-36}\;\text{s}$$

The running of the Standard Model gauge couplings $\alpha_1, \alpha_2, \alpha_3$ suggests approximate unification at this scale, particularly in supersymmetric extensions (MSSM).

As the universe cools below $T_{\text{GUT}}$, the unified symmetry group G undergoesspontaneous symmetry breaking via the Higgs mechanism:

$$G \;\xrightarrow{T_{\text{GUT}}}\; SU(3)_C \times SU(2)_L \times U(1)_Y$$

Popular candidates for G include SU(5) (Georgi-Glashow), SO(10), and $E_6$.

3.2 Baryogenesis

The GUT phase transition is one candidate epoch for baryogenesis — the generation of the observed matter-antimatter asymmetry. Sakharov (1967) identified three necessary conditions:

1. Baryon Number Violation

Processes that change $B$ must exist. In GUTs, leptoquark gauge bosons X and Y mediate proton decay and $\Delta B \neq 0$ interactions.

2. C and CP Violation

To produce an excess of baryons over antibaryons, both charge conjugation (C) and the combined CP symmetry must be violated in baryon-number-changing processes.

3. Departure from Equilibrium

In thermal equilibrium, CPT invariance guarantees equal baryon and antibaryon densities. A departure from equilibrium (e.g., out-of-equilibrium decay of heavy X bosons) is required.

The observed baryon asymmetry is quantified by the baryon-to-photon ratio:

$$\eta \equiv \frac{n_b - n_{\bar{b}}}{n_\gamma} = \frac{n_b}{n_\gamma} \approx 6.1 \times 10^{-10}$$

Measured independently from BBN light element abundances and from CMB acoustic peak ratios (Planck 2018).

4. The Electroweak Phase Transition

At a temperature $T_{\text{EW}} \sim 100$ GeV ($t \sim 10^{-12}$ s), the electroweak symmetry $SU(2)_L \times U(1)_Y$ is spontaneously broken to the electromagnetic $U(1)_{\text{em}}$ by the Higgs mechanism:

$$SU(2)_L \times U(1)_Y \;\xrightarrow{\langle\phi\rangle = v/\sqrt{2}}\; U(1)_{\text{em}}$$

The Higgs field acquires its vacuum expectation value $v = 246$ GeV, giving mass to the $W^\pm$ and $Z^0$ bosons while the photon remains massless.

The nature of this transition has profound implications. In the Standard Model with a Higgs mass of 125 GeV, lattice calculations show the EW transition is a smooth crossover, not a first-order phase transition. However, extensions of the Standard Model (e.g., additional scalar fields, MSSM) can make it strongly first-order, which would satisfy Sakharov's third condition forelectroweak baryogenesis.

Sphaleron Processes

Above the EW transition temperature, sphaleron processes — non-perturbative field configurations that interpolate between different vacuum states of$SU(2)_L$ — are unsuppressed and violate both baryon number B and lepton number L while conserving $B - L$. Their rate per unit volume is$\Gamma_{\text{sph}} \sim \alpha_W^5 T^4$. Below $T_{\text{EW}}$, sphalerons are exponentially suppressed by the Boltzmann factor $e^{-E_{\text{sph}}/T}$ where$E_{\text{sph}} \approx 8\pi v/g \sim 10$ TeV.

5. The QCD Phase Transition

At $T_{\text{QCD}} \approx 150\text{--}170$ MeV ($t \sim 10^{-5}$ s), the universe undergoes the quark-hadron transition. Above this temperature, quarks and gluons exist as a deconfined quark-gluon plasma (QGP). Below it, the strong interaction confines quarks into color-neutral hadrons (protons, neutrons, pions, etc.).

Key Features of the QCD Transition

Lattice QCD calculations indicate a smooth crossover for physical quark masses, occurring at$T_c \approx 155$ MeV (HotQCD Collaboration, 2019).
The effective degrees of freedom drop dramatically: from $g_* \approx 61.75$ (quarks, gluons, leptons, photons) above $T_c$ to $g_* \approx 17.25$ (pions, muons, electrons, neutrinos, photons) just below.
Chiral symmetry breaking accompanies confinement: the quark condensate$\langle\bar{q}q\rangle \neq 0$ generates constituent quark masses of$\sim 300$ MeV, much larger than the bare current quark masses.
Most hadrons produced are pions, which are light and decay or annihilate quickly. Protons and neutrons freeze out of pair annihilation at $T \sim 20$ MeV, setting the initial baryon inventory for nucleosynthesis.

6. Neutrino Decoupling

Neutrinos interact with the thermal plasma via weak interactions. The weak interaction rate for processes such as $\nu_e + e^- \leftrightarrow \nu_e + e^-$ scales as:

$$\Gamma_{\text{weak}} \sim G_F^2\,T^5$$

where $G_F = 1.166 \times 10^{-5}\;\text{GeV}^{-2}$ is the Fermi constant.

Comparing this to the Hubble expansion rate $H \sim \sqrt{G}\,T^2 \sim T^2/M_{\text{Pl}}$:

$$\frac{\Gamma_{\text{weak}}}{H} \sim G_F^2\,M_{\text{Pl}}\,T^3 \sim \left(\frac{T}{1\;\text{MeV}}\right)^3$$

When $\Gamma_{\text{weak}} < H$, neutrinos can no longer maintain thermal equilibrium and they decouple (freeze out) from the plasma.

Neutrino decoupling occurs at $T_{\nu,\text{dec}} \approx 1$ MeV ($t \approx 1$ s). After decoupling, neutrinos free-stream and their temperature redshifts as $T_\nu \propto a^{-1}$, independently of the photon temperature. This distinction becomes important after $e^+e^-$ annihilation.

Neutrino Temperature After $e^+e^-$ Annihilation

$$\boxed{T_\nu = \left(\frac{4}{11}\right)^{1/3} T_\gamma \approx 0.714\,T_\gamma}$$

Today: $T_\gamma = 2.725$ K and $T_\nu = 1.945$ K (the cosmic neutrino background).

7. Big Bang Nucleosynthesis (BBN)

Big Bang Nucleosynthesis is one of the three pillars of the hot Big Bang model. It predicts the primordial abundances of light elements synthesized in the first few minutes of the universe. BBN is sensitive to a single free parameter — the baryon-to-photon ratio $\eta$ — and the number of light neutrino species $N_\nu$.

7.1 Neutron-to-Proton Freeze-Out

At temperatures $T \gg 1$ MeV, neutrons and protons are maintained in chemical equilibrium by weak interactions:

$$n + \nu_e \leftrightarrow p + e^-$$

$$n + e^+ \leftrightarrow p + \bar{\nu}_e$$

$$n \leftrightarrow p + e^- + \bar{\nu}_e$$

In thermal equilibrium, the neutron-to-proton ratio is determined by the Boltzmann factor:

$$\frac{n_n}{n_p} = \exp\!\left(-\frac{\Delta m\,c^2}{k_B T}\right) = \exp\!\left(-\frac{Q}{T}\right)$$

where $Q = (m_n - m_p)c^2 = 1.293$ MeV is the neutron-proton mass difference.

The weak interaction rates for these processes are:

$$\lambda_{n \to p} = \frac{1+3g_A^2}{2\pi^3}\,G_F^2\!\int_0^\infty \frac{E_e^2\,(Q+E_e)^2}{1+e^{E_e/T}}\,\frac{1}{1+e^{-(Q+E_e)/T_\nu}}\,dE_e$$

where $g_A \approx 1.2723$ is the axial-vector coupling constant.

The total weak rate $\Gamma_{n \leftrightarrow p}$ falls below the Hubble rate at thefreeze-out temperature:

$$T_f \approx 0.7\;\text{MeV} \qquad \Rightarrow \qquad \left.\frac{n}{p}\right|_f = e^{-Q/T_f} \approx e^{-1.293/0.7} \approx \frac{1}{6}$$

After freeze-out, the neutron-to-proton ratio evolves only through free neutron decay ($\tau_n = 879.4 \pm 0.6$ s):

$$\frac{n}{p}(t) = \frac{n}{p}\bigg|_f \cdot e^{-t/\tau_n}$$

By the time nucleosynthesis begins at $t_{\text{BBN}} \approx 180$ s, the ratio has decreased to $n/p \approx 1/7$.

7.2 The Deuterium Bottleneck

Nuclear reactions cannot proceed efficiently until deuterium (D) can survive photodissociation. The relevant reaction is:

$$p + n \leftrightarrow D + \gamma \qquad (B_D = 2.224\;\text{MeV})$$

Although the deuterium binding energy is $B_D = 2.224$ MeV, the enormous photon-to-baryon ratio $\eta^{-1} \sim 1.6 \times 10^9$ means that the high-energy tail of the Planck distribution contains enough photons to dissociate deuterium even at temperatures well below $B_D$. Deuterium becomes stable only when:

$$T_{\text{BBN}} \approx 0.07\;\text{MeV} \approx 8 \times 10^8\;\text{K} \qquad (t \approx 180\;\text{s} \approx 3\;\text{min})$$

The factor $\sim 30$ between $B_D/T_{\text{BBN}}$ and unity reflects the logarithmic sensitivity to the enormous photon-to-baryon ratio:$T_{\text{BBN}} \sim B_D / \ln(\eta^{-1})$.

7.3 Light Element Synthesis

Once deuterium survives, a rapid chain of nuclear reactions builds up heavier elements. The key reactions are:

$$D + p \to {}^3\text{He} + \gamma$$

$$D + D \to {}^3\text{He} + n \qquad D + D \to T + p$$

$$\,{}^3\text{He} + n \to T + p \qquad T + D \to {}^4\text{He} + n$$

$$\,{}^3\text{He} + D \to {}^4\text{He} + p$$

The chain terminates predominantly at ${}^4\text{He}$ because there are no stable nuclei with mass number A = 5 or A = 8 (the “mass-5 and mass-8 gaps”).

7.4 The Helium-4 Mass Fraction

Since virtually all neutrons end up in ${}^4\text{He}$ (the most tightly bound light nucleus), the primordial helium mass fraction can be estimated analytically:

Primordial Helium Mass Fraction

$$Y_p = \frac{4 \cdot n_{{}^4\text{He}}}{n_b} = \frac{4 \cdot (n_n/2)}{n_n + n_p} = \frac{2(n/p)}{1 + (n/p)}$$

$$\boxed{Y_p = \frac{2(n/p)}{1+(n/p)} \approx \frac{2 \times 1/7}{1 + 1/7} = \frac{2/7}{8/7} = \frac{1}{4} = 0.25}$$

Detailed numerical calculations give $Y_p = 0.2470 \pm 0.0002$ for the Planck best-fit baryon density. Observations of metal-poor H II regions yield $Y_p = 0.245 \pm 0.003$, in excellent agreement.

7.5 Trace Element Abundances

Predicted Primordial Abundances (by number relative to H)

Element	BBN Prediction	Observation	Sensitivity
${}^4\text{He}$ (mass fraction)	$0.2470 \pm 0.0002$	$0.245 \pm 0.003$	Sensitive to $N_\nu$, $\tau_n$
D/H	$(2.57 \pm 0.13) \times 10^{-5}$	$(2.55 \pm 0.03) \times 10^{-5}$	Most sensitive to $\eta$
${}^3\text{He}/\text{H}$	$\sim 10^{-5}$	$\sim 1.1 \times 10^{-5}$	Complicated by stellar production
${}^7\text{Li}/\text{H}$	$(4.7 \pm 0.5) \times 10^{-10}$	$(1.6 \pm 0.3) \times 10^{-10}$	The “lithium problem”

The Cosmological Lithium Problem

The predicted primordial ${}^7\text{Li}$ abundance exceeds the observed value in metal-poor halo stars (the “Spite plateau”) by a factor of ~3. This longstanding discrepancy may indicate stellar depletion mechanisms, systematic errors in nuclear cross sections (particularly the ${}^3\text{He}(\alpha,\gamma){}^7\text{Be}$ rate), or new physics beyond the Standard Model.

7.6 BBN as a Precision Cosmological Probe

The Boltzmann equation governing the evolution of nuclear abundances takes the form:

$$\frac{dY_i}{dt} = \sum_{j,k,l} N_i\!\left[\frac{Y_j^{N_j} Y_k^{N_k}}{N_j!\,N_k!}\,\langle\sigma v\rangle_{jk\to il}\,\rho_b^{N_j+N_k-1} - \frac{Y_i^{N_i} Y_l^{N_l}}{N_i!\,N_l!}\,\langle\sigma v\rangle_{il\to jk}\,\rho_b^{N_i+N_l-1}\right]$$

where $Y_i = n_i/n_b$ is the abundance by number and $\langle\sigma v\rangle$ denotes thermally averaged cross sections.

BBN constrains fundamental physics in several important ways:

Number of neutrino species: Additional relativistic species increase $g_*$, raising $H$, causing earlier freeze-out at higher $n/p$, thus increasing $Y_p$. BBN constrains$N_{\text{eff}} = 2.89 \pm 0.28$, consistent with $N_\nu = 3$.
Baryon density: D/H is the best “baryometer,” yielding $\Omega_b h^2 = 0.0224 \pm 0.0001$, in remarkable agreement with the CMB determination.
Neutron lifetime: $Y_p$ depends sensitively on $\tau_n$ through the freeze-out temperature and subsequent decay.

8. Electron-Positron Annihilation and Photon Reheating

When the temperature drops below $T \sim m_e c^2 \approx 0.511$ MeV, electron-positron pairs annihilate: $e^+ + e^- \to \gamma + \gamma$. The entropy released heats the photon bath but not the already-decoupled neutrinos.

Entropy Conservation Argument

The entropy density of relativistic species is $s = (2\pi^2/45)\,g_{*s}\,T^3$. Since entropy is conserved in a comoving volume ($sa^3 = \text{const}$):

$$g_{*s}^{\text{before}}\,T_{\text{before}}^3\,a_{\text{before}}^3 = g_{*s}^{\text{after}}\,T_{\text{after}}^3\,a_{\text{after}}^3$$

Before: $g_{*s} = 2 + \frac{7}{8} \times 4 = \frac{11}{2}$ (photons + $e^\pm$). After: $g_{*s} = 2$ (photons only).

$$\frac{T_\gamma^{\text{after}}}{T_\gamma^{\text{before}}} = \left(\frac{g_{*s}^{\text{before}}}{g_{*s}^{\text{after}}}\right)^{1/3} = \left(\frac{11/2}{2}\right)^{1/3} = \left(\frac{11}{4}\right)^{1/3} \approx 1.401$$

Since neutrinos have already decoupled and simply redshift as $T_\nu \propto a^{-1}$, while photons receive additional entropy from $e^\pm$ annihilation:

$$\boxed{\frac{T_\nu}{T_\gamma} = \left(\frac{4}{11}\right)^{1/3} \approx 0.7138}$$

This ratio is maintained for all subsequent evolution. The effective number of neutrino species contributing to the radiation density today is $N_{\text{eff}} = 3.044$ (slightly above 3 due to partial $e^\pm$ heating of neutrinos before full decoupling).

9. Recombination

As the universe cools to $T \sim 0.3$ eV ($\sim 3000$ K), free electrons combine with protons to form neutral hydrogen. This process — historically called “recombination” despite the fact that electrons and protons had never been combined before — is a pivotal event: it makes the universe transparent to photons.

9.1 The Saha Equation

In thermal equilibrium, the ionization fraction is governed by the Saha equation, derived from chemical equilibrium ($\mu_p + \mu_e = \mu_H$) for the reaction$p + e^- \leftrightarrow H + \gamma$:

The Saha Equation

$$\boxed{\frac{n_p\,n_e}{n_H} = \left(\frac{m_e\,T}{2\pi}\right)^{3/2}\exp\!\left(-\frac{B}{T}\right)}$$

where $B = 13.6$ eV is the hydrogen ionization energy (natural units: $\hbar = c = k_B = 1$).

Defining the free electron fraction $x_e \equiv n_e/n_b = n_p/n_b$ and using$n_H = n_b(1-x_e)$, $n_p = n_b x_e$:

$$\frac{x_e^2}{1-x_e} = \frac{1}{n_b}\left(\frac{m_e T}{2\pi}\right)^{3/2} e^{-B/T}$$

Using $n_b = \eta\,n_\gamma = \eta \cdot \frac{2\zeta(3)}{\pi^2}\,T^3$:

$$\frac{x_e^2}{1-x_e} = \frac{\pi^2}{2\zeta(3)\,\eta}\left(\frac{m_e}{2\pi T}\right)^{-3/2}\left(\frac{T}{m_e}\right)^{3/2} e^{-B/T} = \frac{1}{\eta}\cdot\frac{\sqrt{\pi}}{4\zeta(3)}\left(\frac{2T}{m_e}\right)^{3/2}e^{-B/T}$$

The Saha equation predicts that $x_e = 0.5$ (50% ionized) at:

$$T_{\text{rec}} \approx 0.32\;\text{eV} \approx 3740\;\text{K}, \qquad z_{\text{rec}} \approx 1370$$

Note that $T_{\text{rec}} \ll B = 13.6$ eV because of the enormous entropy ($\eta^{-1} \sim 10^9$ photons per baryon): even though $T \ll B$, the Boltzmann tail still contains enough high-energy photons to ionize hydrogen.

9.2 Beyond Saha: Peebles' Three-Level Atom

The Saha equation assumes perfect thermal equilibrium, but recombination is actually anon-equilibrium process. Recombination directly to the ground state produces a Lyman-continuum photon that immediately ionizes another atom. The effective recombination proceeds through excited states, and Lyman-alpha photons can be trapped.

Peebles (1968) developed a simplified three-level model (ground state, first excited state, continuum) leading to an effective recombination equation:

$$\frac{dx_e}{dt} = -C\!\left[\alpha^{(2)}(T)\,n_b\,x_e^2 - \beta^{(2)}(T)\,(1-x_e)\,e^{-B_{2s}/T}\right]$$

where $\alpha^{(2)}$ is the case-B recombination coefficient (excluding captures to the ground state), $\beta^{(2)}$ is the corresponding photoionization rate from n=2,$B_{2s} = B/4 = 3.4$ eV, and C is the Peebles correction factor:

$$C = \frac{\Lambda_{2s\to 1s} + \Lambda_\alpha}{\Lambda_{2s\to 1s} + \Lambda_\alpha + \beta^{(2)}(T)}$$

Here $\Lambda_{2s\to 1s} = 8.22\;\text{s}^{-1}$ is the two-photon decay rate from the metastable 2s state, and $\Lambda_\alpha$ is the net Lyman-alpha escape rate (cosmological redshifting removes photons from the Ly$\alpha$ line).

The Peebles model shows that recombination proceeds more slowly than the Saha prediction, with the ionization fraction “freezing out” at a residual level:

$$x_e^{\text{freeze}} \approx 2.7 \times 10^{-5}\;\left(\frac{\Omega_b h^2}{0.022}\right)^{-1/2}$$

Modern codes (RECFAST, CosmoRec, HyRec) achieve sub-percent accuracy by including multi-level hydrogen, helium recombination, and radiative transfer effects.

10. Photon Decoupling and the Last Scattering Surface

Photons decouple from matter when the mean free path for Thomson scattering exceeds the Hubble radius. The Thomson scattering rate is:

$$\Gamma_T = n_e\,\sigma_T\,c = x_e\,n_b\,\sigma_T\,c$$

where $\sigma_T = 6.65 \times 10^{-25}\;\text{cm}^2$ is the Thomson cross section.

The optical depth looking back from today to redshift z is:

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

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

$$\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 sharply peaked at $z_* \approx 1090$ ($T_* \approx 2970$ K) with a width $\Delta z \approx 80$. This defines the last scattering surfacefrom which the CMB photons we observe today were released.

Physical Interpretation

The last scattering surface is not infinitely thin but has a finite thickness ($\Delta z \approx 80$), corresponding to a comoving distance of about 15 Mpc. This finite thickness causes small-scale damping of CMB anisotropies (Silk damping) and sets a minimum angular scale for primary CMB fluctuations. The optical depth to the last scattering surface is $\tau(z_*) \approx 1$ by definition, while $\tau \gg 1$ for$z \gg z_*$ (the universe is opaque) and $\tau \ll 1$ for $z \ll z_*$(the universe is transparent).

11. The Dark Ages and Cosmic Dawn

After recombination at $z \approx 1100$, the universe enters thecosmic dark ages — a period with no luminous sources. The matter content is predominantly neutral hydrogen and helium, slowly cooling adiabatically as$T_{\text{gas}} \propto (1+z)^2$ (faster than the CMB temperature$T_\gamma \propto (1+z)$ because the gas decouples thermally from the radiation).

Timeline of the Post-Recombination Universe

$z \sim 1100\text{--}200$: Dark ages proper. No stars or galaxies. Small density perturbations grow linearly under gravity:$\delta \propto a$ (matter domination).
$z \sim 200\text{--}30$: First dark matter halos form via gravitational collapse. The 21-cm absorption signal from neutral hydrogen against the CMB becomes potentially observable.
$z \sim 30\text{--}15$: “Cosmic dawn” — the first Population III stars ignite in minihalos of mass$M \sim 10^{5}\text{--}10^{6}\;M_\odot$, powered by primordial (zero-metallicity) gas.
$z \sim 15\text{--}6$: The epoch of reionization. UV photons from early galaxies and quasars gradually re-ionize the intergalactic medium. Observations of Gunn-Peterson troughs in high-z quasar spectra indicate reionization was largely complete by $z \approx 6$. The CMB polarization constrains the optical depth to reionization: $\tau_{\text{reion}} = 0.054 \pm 0.007$ (Planck 2018).

12. Numerical Simulation: Recombination History

The following Python code computes the recombination history of the universe by solving the Saha equation and the Peebles equation for the ionization fraction $x_e(z)$. It also computes the optical depth and visibility function.

Python

recombination_history.py218 lines

#!/usr/bin/env python3
"""
recombination_history.py
Compute the ionization fraction x_e(z), optical depth tau(z),
and visibility function g(z) for the recombination epoch.
Implements both Saha and Peebles three-level atom equations.
(Following Dodelson & Schmidt, "Modern Cosmology", eq. 4.36)
"""
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ============================
# Physical Constants (CGS)
# ============================
k_B   = 1.3806e-16       # Boltzmann constant [erg/K]
m_e   = 9.1094e-28        # Electron mass [g]
m_p   = 1.6726e-24        # Proton mass [g]
hbar  = 1.0546e-27        # Reduced Planck constant [erg s]
c     = 2.9979e10          # Speed of light [cm/s]
sigma_T = 6.6524e-25      # Thomson cross section [cm^2]
eV    = 1.6022e-12         # 1 eV in erg
G     = 6.6743e-8          # Gravitational constant [cm^3/g/s^2]
Mpc   = 3.0857e24          # 1 Mpc in cm

# ============================
# Cosmological Parameters (Planck 2018)
# ============================
H0_km = 67.4              # [km/s/Mpc]
H0    = H0_km * 1e5 / Mpc  # [1/s]
T_cmb0 = 2.7255           # CMB temperature today [K]
Omega_b = 0.0493
Omega_m = 0.315
Omega_r = 9.15e-5         # Radiation density
Omega_L = 0.685
Y_He = 0.245              # Helium mass fraction

# Hydrogen energy levels
B_H   = 13.6 * eV         # Ground-state binding energy [erg]
B_2s  = B_H / 4.0         # n=2 binding energy: 3.4 eV [erg]
E_Lya = B_H - B_2s        # Lyman-alpha energy: 10.2 eV [erg]

# Two-photon decay rate from 2s state
Lambda_2s = 8.2245         # [s^-1]

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

def n_H(z):
    """Hydrogen number density at redshift z [cm^-3]."""
    rho_crit = 3 * H0**2 / (8 * np.pi * G)
    n_H0 = Omega_b * rho_crit * (1 - Y_He) / m_p
    return n_H0 * (1 + z)**3

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

def alpha_B(T):
    """Case-B recombination coefficient [cm^3/s]. T in Kelvin."""
    t = T / 1e4
    return 4.309e-13 * t**(-0.6166) / (1 + 0.6703 * t**0.5300)

def beta_B(T):
    """Photoionization rate from n=2 [s^-1]. T in Kelvin."""
    return alpha_B(T) * (m_e * k_B * T / (2*np.pi*hbar**2))**1.5 \
           * np.exp(-B_2s / (k_B * T))

# ============================
# Saha Equation
# ============================
def saha_xe(z):
    """Ionization fraction from the Saha equation."""
    T_erg = T_cmb(z) * k_B
    nb = n_H(z)
    rhs = (1.0/nb) * (m_e * T_erg / (2*np.pi*hbar**2))**1.5 \
          * np.exp(-B_H / T_erg)
    disc = rhs**2 + 4*rhs
    xe = (-rhs + np.sqrt(max(disc, 0.0))) / 2.0
    return min(max(xe, 0.0), 1.0)

# ============================
# Peebles Equation (ODE)
# ============================
def peebles_rhs(z, xe_arr):
    """RHS of dx_e/dz for the Peebles three-level atom."""
    xe = float(np.ravel(xe_arr)[0])
    xe = max(min(xe, 1.0 - 1e-15), 1e-30)

T_K  = T_cmb(z)
    nb   = n_H(z)
    H_val = Hz(z)

aB = alpha_B(T_K)
    bB = beta_B(T_K)

# Lyman-alpha escape rate
    n_1s = max((1 - xe) * nb, 1e-30)
    lambda_alpha = H_val * (3*B_H)**3 / ((8*np.pi)**2 * n_1s * (hbar*c)**3)

# Peebles C factor
    C = (Lambda_2s + lambda_alpha) / (Lambda_2s + lambda_alpha + bB)

# Dodelson eq. (4.36): Lyman-alpha exponential converts n=2 to ground state
    dxe_dt = C * (bB * (1 - xe) * np.exp(-E_Lya / (k_B * T_K))
                  - aB * nb * xe**2)

# Convert dx_e/dt -> dx_e/dz
    dxe_dz = dxe_dt / (-(1+z) * H_val)
    return [dxe_dz]

# ============================
# Solve the System
# ============================
z_start = 2000
z_end   = 200
z_arr = np.linspace(z_start, z_end, 5000)

# Saha solution over full range
xe_saha = np.array([saha_xe(z) for z in z_arr])

# Switch from Saha to Peebles when x_e drops below 0.99
z_switch = z_start
for z in z_arr:
    if saha_xe(z) < 0.99:
        z_switch = z
        break

# Build combined solution: Saha at high z, Peebles at lower z
xe_peebles = np.zeros_like(z_arr)
switch_idx = 0
for i in range(len(z_arr)):
    if z_arr[i] >= z_switch:
        xe_peebles[i] = saha_xe(z_arr[i])
        switch_idx = i
    else:
        break

# Solve Peebles ODE from switch point using stiff solver (Radau)
z_ode = z_arr[switch_idx:]
sol = solve_ivp(peebles_rhs,
    [z_ode[0], z_ode[-1]], [saha_xe(z_switch)],
    t_eval=z_ode, method='Radau',
    rtol=1e-8, atol=1e-12, max_step=1.0)
xe_peebles[switch_idx:switch_idx+len(sol.y[0])] = sol.y[0]
xe_peebles = np.clip(xe_peebles, 1e-30, 1.0)

# ============================
# Optical Depth and Visibility Function
# ============================
# tau(z) = optical depth from observer (z=0) to redshift z
tau_arr = np.zeros_like(z_arr)
for i in range(len(z_arr)-2, -1, -1):
    dz = z_arr[i] - z_arr[i+1]
    nb = n_H(z_arr[i])
    dtau = xe_peebles[i] * nb * sigma_T * c / (Hz(z_arr[i]) * (1+z_arr[i]))
    tau_arr[i] = tau_arr[i+1] + dtau * dz

# Visibility function: g(z) = |dtau/dz| * exp(-tau)
dtau_vals = np.gradient(tau_arr, z_arr)
g_arr = np.abs(dtau_vals) * np.exp(-tau_arr)
norm = np.trapezoid(g_arr, z_arr)
if abs(norm) > 0:
    g_arr /= abs(norm)

# ============================
# Print Key Results
# ============================
idx_half = np.argmin(np.abs(xe_peebles - 0.5))
z_rec = z_arr[idx_half]
T_rec = T_cmb(z_rec)
idx_peak = np.argmax(g_arr)
print("=" * 60)
print("  RECOMBINATION HISTORY (Three-Level Atom)")
print("=" * 60)
print(f"  Saha-Peebles switch:        z = {z_switch:.0f}")
print(f"  Recombination (x_e = 0.5):  z_rec = {z_rec:.0f}")
print(f"  Recombination temperature:  T_rec = {T_rec:.0f} K")
print(f"  Residual ionization (z=200):x_e = {xe_peebles[-1]:.2e}")
print(f"  Visibility function peak:   z_* = {z_arr[idx_peak]:.0f}")
print(f"  (Planck 2018 value:         z_* ~ 1089)")
print("=" * 60)

# ============================
# Plot Results
# ============================
fig, axes = plt.subplots(3, 1, figsize=(8, 12), sharex=True)

axes[0].semilogy(z_arr, np.maximum(xe_saha, 1e-30),
                 'b--', label='Saha equation', alpha=0.7)
axes[0].semilogy(z_arr, np.maximum(xe_peebles, 1e-30),
                 'r-', label='Peebles equation', lw=2)
axes[0].set_ylabel(r'$x_e$', fontsize=14)
axes[0].set_ylim(1e-4, 2)
axes[0].legend(fontsize=12)
axes[0].set_title('Recombination History', fontsize=14)
axes[0].invert_xaxis()

axes[1].semilogy(z_arr, np.maximum(tau_arr, 1e-30), 'g-', lw=2)
axes[1].axhline(1, color='gray', ls='--', alpha=0.5, label=r'$\tau = 1$')
axes[1].set_ylabel(r'$\tau(z)$', fontsize=14)
axes[1].legend(fontsize=12)

axes[2].plot(z_arr, g_arr, 'm-', lw=2)
axes[2].set_xlabel('Redshift z', fontsize=14)
axes[2].set_ylabel(r'$g(z)$', fontsize=14)
axes[2].set_title('Visibility Function', fontsize=14)

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

Click Run to execute the Python code

Code will be executed with Python 3 on the server

13. Summary: The Thermal History Timeline

Complete Thermal History Timeline

Epoch	Time	Temperature	Key Physics
Planck epoch	$< 10^{-43}$ s	$> 10^{19}$ GeV	Quantum gravity regime
GUT transition	$\sim 10^{-36}$ s	$\sim 10^{16}$ GeV	Symmetry breaking, baryogenesis (?)
EW transition	$\sim 10^{-12}$ s	$\sim 100$ GeV	Higgs mechanism, EW baryogenesis (?)
QCD transition	$\sim 10^{-5}$ s	$\sim 150$ MeV	Quark confinement into hadrons
$\nu$ decoupling	$\sim 1$ s	$\sim 1$ MeV	Neutrino freeze-out
n/p freeze-out	$\sim 2$ s	$\sim 0.7$ MeV	$n/p \approx 1/6$
$e^+e^-$ annihilation	$\sim 10$ s	$\sim 0.5$ MeV	Photon reheating, $T_\nu/T_\gamma = (4/11)^{1/3}$
BBN	3 min	$\sim 0.07$ MeV	${}^4\text{He}$, D, ${}^7\text{Li}$ synthesis
Matter-radiation eq.	$\sim 47$ kyr	$\sim 0.75$ eV	$z_{\text{eq}} \approx 3400$
Recombination	$\sim 370$ kyr	$\sim 0.26$ eV	$z_* \approx 1090$, CMB released
Reionization	$\sim 0.1\text{--}1$ Gyr	$\sim$ meV	$z \sim 6\text{--}15$, first stars/galaxies
Today	13.8 Gyr	$T_\gamma = 2.725$ K	Dark energy domination

Key Equations of Thermal History

Friedmann Equation (Radiation Era)

$$H^2 = \frac{8\pi G}{3}\cdot\frac{\pi^2}{30}\,g_*(T)\,T^4$$

Neutron-Proton Freeze-Out

$$\frac{n_n}{n_p} = \exp\!\left(-\frac{Q}{T_f}\right), \qquad T_f \approx 0.7\;\text{MeV}$$

BBN Helium Mass Fraction

$$Y_p = \frac{2(n/p)}{1+(n/p)} \approx 0.245$$

Neutrino Temperature

$$T_\nu = \left(\frac{4}{11}\right)^{1/3} T_\gamma$$

Saha Equation (Recombination)

$$\frac{n_p\,n_e}{n_H} = \left(\frac{m_e T}{2\pi}\right)^{3/2}\exp\!\left(-\frac{B}{T}\right)$$

Visibility Function

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

Logical Chain of the Thermal History

1. Friedmann equation + radiation EOS $\Rightarrow$ $T \propto a^{-1}$, $t \propto T^{-2}$

2. Interaction rate vs. Hubble rate $\Gamma/H$ determines freeze-out of each species

3. Neutrinos decouple at $T \sim 1$ MeV; $e^\pm$ annihilation reheats photons by $(11/4)^{1/3}$

4. n/p freeze-out at $T_f \sim 0.7$ MeV sets initial conditions for BBN

5. Deuterium bottleneck delays nucleosynthesis until $T \sim 0.07$ MeV; neutron decay reduces n/p to 1/7

6. Rapid nuclear burning produces $Y_p \approx 0.25$; trace D, ${}^3\text{He}$, ${}^7\text{Li}$

7. Saha equation governs recombination at $T \sim 0.3$ eV; Peebles equation for non-equilibrium corrections

8. Photon decoupling at $z_* \approx 1090$ releases the CMB; visibility function peaked at last scattering

9. Dark ages $\to$ cosmic dawn $\to$ reionization ($z \sim 6\text{--}15$)

10. Concordance: BBN and CMB independently constrain $\Omega_b h^2$ and $N_\nu$ with remarkable agreement

Bibliography

Textbooks & Monographs

Kolb, E.W. & Turner, M.S. (1990). The Early Universe. Addison-Wesley. — The definitive reference for thermal history, Big Bang nucleosynthesis, and particle cosmology.
Weinberg, S. (2008). Cosmology. Oxford University Press. — Rigorous treatment of recombination, nucleosynthesis, and the thermal history of the universe.
Dodelson, S. & Schmidt, F. (2020). Modern Cosmology, 2nd ed. Academic Press. — Detailed derivation of the Peebles equation and Saha equilibrium.
Mukhanov, V. (2005). Physical Foundations of Cosmology. Cambridge University Press. — Comprehensive treatment of thermodynamics, nucleosynthesis, and recombination physics.
Ryden, B. (2017). Introduction to Cosmology, 2nd ed. Cambridge University Press. — Accessible introduction to BBN, recombination, and the CMB formation.
Bernstein, J. (1988). Kinetic Theory in the Expanding Universe. Cambridge University Press. — Boltzmann equation approach to freeze-out and relic abundances.

Key Papers

Gamow, G. (1948). “The Evolution of the Universe,” Nature 162, 680–682. — Prediction of a hot early universe and relic radiation.
Alpher, R.A., Bethe, H. & Gamow, G. (1948). “The Origin of Chemical Elements,” Physical Review 73, 803–804. — The celebrated αβγ paper on primordial nucleosynthesis.
Peebles, P.J.E. (1968). “Recombination of the Primeval Plasma,” Astrophysical Journal 153, 1–11. — The three-level atom model for hydrogen recombination (Peebles equation).
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. — Independent derivation of the recombination equation.
Penzias, A.A. & Wilson, R.W. (1965). “A Measurement of Excess Antenna Temperature at 4080 Mc/s,” Astrophysical Journal 142, 419–421. — Discovery of the cosmic microwave background radiation.
Fixsen, D.J. (2009). “The Temperature of the Cosmic Microwave Background,” Astrophysical Journal 707, 916–920. — Precise measurement of T = 2.7255 ± 0.0006 K.
Planck Collaboration (2020). “Planck 2018 results. VI. Cosmological parameters,” Astronomy & Astrophysics 641, A6. arXiv:1807.06209. — Precision constraints on recombination redshift and optical depth.
Fields, B.D., Olive, K.A., Yeh, T.-H. & Young, C. (2020). “Big-Bang Nucleosynthesis after Planck,” Journal of Cosmology and Astroparticle Physics 2020(03), 010. arXiv:1912.01132. — Modern review of BBN predictions and observational constraints.

← Friedmann Equations Cosmic Inflation →

Element	BBN Prediction	Observation	Sensitivity
\({}^4\text{He}\) (mass fraction)	\(0.2470 \pm 0.0002\)	\(0.245 \pm 0.003\)	Sensitive to \(N_\nu\), \(\tau_n\)
D/H	\((2.57 \pm 0.13) \times 10^{-5}\)	\((2.55 \pm 0.03) \times 10^{-5}\)	Most sensitive to \(\eta\)
\({}^3\text{He}/\text{H}\)	\(\sim 10^{-5}\)	\(\sim 1.1 \times 10^{-5}\)	Complicated by stellar production
\({}^7\text{Li}/\text{H}\)	\((4.7 \pm 0.5) \times 10^{-10}\)	\((1.6 \pm 0.3) \times 10^{-10}\)	The “lithium problem”

Epoch	Time	Temperature	Key Physics
Planck epoch	\(< 10^{-43}\) s	\(> 10^{19}\) GeV	Quantum gravity regime
GUT transition	\(\sim 10^{-36}\) s	\(\sim 10^{16}\) GeV	Symmetry breaking, baryogenesis (?)
EW transition	\(\sim 10^{-12}\) s	\(\sim 100\) GeV	Higgs mechanism, EW baryogenesis (?)
QCD transition	\(\sim 10^{-5}\) s	\(\sim 150\) MeV	Quark confinement into hadrons
\(\nu\) decoupling	\(\sim 1\) s	\(\sim 1\) MeV	Neutrino freeze-out
n/p freeze-out	\(\sim 2\) s	\(\sim 0.7\) MeV	\(n/p \approx 1/6\)
\(e^+e^-\) annihilation	\(\sim 10\) s	\(\sim 0.5\) MeV	Photon reheating, \(T_\nu/T_\gamma = (4/11)^{1/3}\)
BBN	3 min	\(\sim 0.07\) MeV	\({}^4\text{He}\), D, \({}^7\text{Li}\) synthesis
Matter-radiation eq.	\(\sim 47\) kyr	\(\sim 0.75\) eV	\(z_{\text{eq}} \approx 3400\)
Recombination	\(\sim 370\) kyr	\(\sim 0.26\) eV	\(z_* \approx 1090\), CMB released
Reionization	\(\sim 0.1\text{--}1\) Gyr	\(\sim\) meV	\(z \sim 6\text{--}15\), first stars/galaxies
Today	13.8 Gyr	\(T_\gamma = 2.725\) K	Dark energy domination

Cosmology