Part II: Thermal History | Chapter 2

Dark Ages & Reionization

From the last scattering surface through the cosmic dark ages to the epoch of reionization: the 21-cm hydrogen signal, first star formation, and the re-illumination of the universe

Historical Context

The concept of a “dark age” between recombination and the first luminous objects was recognized in the 1960s, but its observational study only became feasible with advances in radio astronomy and space-based telescopes. Gunn & Peterson (1965) predicted that neutral hydrogen in the intergalactic medium (IGM) would produce complete absorption of Lyman-α photons from distant quasars, providing a probe of the reionization epoch. Decades later, Fan et al. (2006) used SDSS quasars at $z > 6$ to reveal Gunn-Peterson troughs, confirming that reionization was completing around that epoch.

The EDGES collaboration (Bowman et al. 2018) reported a tentative detection of the global 21-cm absorption signal at $z \approx 17$, corresponding to the era when the first stars began coupling the hydrogen spin temperature to the cold gas kinetic temperature. This measurement, while still debated, opened a new window into the cosmic dawn. Most recently, JWST has revolutionized our understanding by discovering luminous galaxies at $z > 10$, challenging existing models of early galaxy formation and suggesting that reionization may have started earlier than expected.

The physics of the dark ages and reionization connects atomic physics (the 21-cm hyperfine transition), gravitational collapse (first structure formation), stellar astrophysics (Population III stars), and radiative transfer (ionization fronts), making it one of the richest frontiers in modern cosmology.

1. Introduction: The Cosmic Dark Ages and Reionization

After recombination at $z \approx 1100$, the universe entered a period known as the cosmic dark ages. The CMB photons decoupled from matter, and the universe was filled with neutral hydrogen and helium gas but no luminous sources. This era persisted from roughly $z \sim 1100$ down to $z \sim 30$, when the first gravitationally collapsed structures began to form.

The key physical processes during the dark ages are governed by the thermal evolution of the gas. After decoupling, the gas temperature $T_K$ evolved adiabatically as:

Gas Temperature Evolution

$$T_K(z) \approx T_{\text{CMB}}(z_{\text{dec}}) \left(\frac{1+z}{1+z_{\text{dec}}}\right)^2 \quad \text{(after Compton decoupling)}$$

The gas cools as $T_K \propto (1+z)^2$ (adiabatic cooling) while the CMB cools as $T_{\text{CMB}} \propto (1+z)$, so the gas becomes colder than the CMB radiation.

Compton scattering keeps $T_K$ coupled to $T_{\text{CMB}}$ down to redshift $z_{\text{dec}} \approx 150$. Below this redshift, the residual free electron fraction is too small to maintain thermal equilibrium, and the gas temperature drops below the CMB temperature. This thermal decoupling sets the stage for the 21-cm signal.

The epoch of reionization (EoR) refers to the period during which ultraviolet radiation from the first stars and galaxies re-ionized the neutral hydrogen in the IGM. Observations constrain the bulk of reionization to occur between $z \sim 6$ and $z \sim 12$, with the midpoint at $z \sim 7\text{--}8$. The process was highly inhomogeneous, with ionized bubbles growing around luminous sources and eventually percolating to fill all of space.

Timeline Summary

$z \approx 1100$: Recombination and CMB last scattering
$z \sim 200\text{--}150$: Compton decoupling; gas begins cooling below CMB
$z \sim 30\text{--}20$: First dark matter halos collapse; molecular hydrogen cooling enables first star formation
$z \sim 20\text{--}15$: Cosmic dawn — first stars (Pop III) emit Lyman-α and UV radiation
$z \sim 15\text{--}6$: Epoch of reionization; ionized bubbles grow and overlap
$z \sim 6$: Reionization essentially complete; Gunn-Peterson troughs disappear

2. The 21-cm Signal from Neutral Hydrogen

The hydrogen 21-cm line arises from the hyperfine splitting of the $1s$ ground state due to the interaction between the proton and electron magnetic moments. The spin-spin coupling energy is:

Hyperfine Transition Energy

$$\Delta E_{21} = \frac{4}{3}\,\alpha^2\,g_p\,\frac{m_e}{m_p}\,E_{1s} \approx 5.87 \times 10^{-6}\;\text{eV}$$

where $\alpha \approx 1/137$ is the fine structure constant, $g_p \approx 5.586$ is the proton g-factor, and $E_{1s} = 13.6$ eV is the hydrogen ground-state energy.

This energy corresponds to a photon with frequency and wavelength:

$$\nu_{21} = \frac{\Delta E_{21}}{h} = 1420.405\;\text{MHz}, \qquad \lambda_{21} = \frac{c}{\nu_{21}} = 21.106\;\text{cm}$$

The spin temperature $T_S$ is defined by the ratio of hydrogen atoms in the triplet (parallel spins, $F=1$) and singlet (antiparallel, $F=0$) states:

Spin Temperature Definition

$$\frac{n_1}{n_0} = 3\,\exp\!\left(-\frac{T_\star}{T_S}\right) \approx 3\left(1 - \frac{T_\star}{T_S}\right)$$

where $T_\star = h\nu_{21}/k_B = 0.068\;\text{K}$ and the factor of 3 reflects the spin degeneracy ratio $g_1/g_0 = 3$.

The spin temperature is determined by three competing processes: (1) absorption/emission of CMB photons driving $T_S \to T_{\text{CMB}}$, (2) collisional coupling driving $T_S \to T_K$, and (3) the Wouthuysen-Field (WF) effect from Lyman-α scattering, also driving $T_S \to T_K$:

$$T_S^{-1} = \frac{T_{\text{CMB}}^{-1} + x_c\,T_K^{-1} + x_\alpha\,T_\alpha^{-1}}{1 + x_c + x_\alpha}$$

where $x_c$ is the collisional coupling coefficient, $x_\alpha$ is the WF coupling coefficient, and $T_\alpha \approx T_K$ is the color temperature of the Lyman-α radiation field.

2.1 Brightness Temperature

The 21-cm signal is observed as a perturbation to the CMB spectrum. The optical depth for the 21-cm transition through a uniform IGM patch is:

21-cm Optical Depth

$$\tau_{21} = \frac{3\,c^3\,\hbar\,A_{10}\,n_{\text{HI}}}{16\,k_B\,\nu_{21}^2\,T_S\,H(z)}$$

where $A_{10} = 2.869 \times 10^{-15}\;\text{s}^{-1}$ is the Einstein A coefficient for the spin-flip transition, $n_{\text{HI}}$ is the neutral hydrogen number density, and $H(z)$ is the Hubble parameter.

The radiative transfer equation along a line of sight through the IGM gives the observed brightness temperature relative to the CMB:

21-cm Brightness Temperature

$$T_b = \frac{T_S - T_{\text{CMB}}}{1 + z}\left(1 - e^{-\tau_{21}}\right) \approx \frac{T_S - T_{\text{CMB}}}{1 + z}\,\tau_{21}$$

The approximation holds because $\tau_{21} \ll 1$ in the diffuse IGM. The sign of $T_b$ reveals the physics:

Absorption ($T_b < 0$): When $T_S < T_{\text{CMB}}$, neutral hydrogen absorbs CMB photons at 21 cm. This occurs during the dark ages and cosmic dawn when the gas is colder than the CMB.
Emission ($T_b > 0$): When $T_S > T_{\text{CMB}}$, the 21-cm line appears in emission. X-ray heating from the first sources eventually raises $T_K$ above $T_{\text{CMB}}$.
No signal ($T_b = 0$): When $T_S = T_{\text{CMB}}$ (CMB coupling dominates) or the hydrogen is fully ionized ($n_{\text{HI}} = 0$).

Substituting the optical depth expression and evaluating at mean density, the brightness temperature becomes:

$$T_b \approx 27\,x_{\text{HI}}\left(1 - \frac{T_{\text{CMB}}}{T_S}\right)\left(\frac{\Omega_b h^2}{0.023}\right)\left(\frac{0.15}{\Omega_m h^2}\,\frac{1+z}{10}\right)^{1/2}\;\text{mK}$$

This is the key observable for 21-cm cosmology experiments such as EDGES, HERA, and the SKA.

3. First Star Formation: Population III Stars

The formation of the first luminous objects requires dark matter halos massive enough to trap baryonic gas and allow it to cool and condense. In the absence of heavy elements (metals), the only available coolant in primordial gas is molecular hydrogen ($\text{H}_2$), which forms through the gas-phase reactions:

$$\text{H} + e^- \to \text{H}^- + \gamma, \qquad \text{H}^- + \text{H} \to \text{H}_2 + e^-$$

The $\text{H}^-$ channel is the dominant formation pathway for $\text{H}_2$ in primordial gas at temperatures $T \sim 200\text{--}2000\;\text{K}$.

3.1 Jeans Mass in the IGM

The minimum mass for gravitational collapse against gas pressure is given by the Jeans mass. In the post-recombination IGM with gas temperature $T_K$ and mean density $\bar{\rho}$:

Jeans Mass

$$M_J = \frac{\pi}{6}\,\bar{\rho}\,\lambda_J^3, \qquad \lambda_J = c_s\sqrt{\frac{\pi}{G\bar{\rho}}}$$

where $c_s = \sqrt{k_B T_K / (\mu m_p)}$ is the sound speed and $\mu \approx 1.22$ is the mean molecular weight of neutral primordial gas.

Evaluating at the mean baryon density of the universe, the cosmological Jeans mass is:

$$M_J \approx 5.7 \times 10^3\,M_\odot \left(\frac{\Omega_m h^2}{0.15}\right)^{-1/2}\left(\frac{\Omega_b h^2}{0.023}\right)^{-3/5}\left(\frac{1+z}{10}\right)^{3/2}$$

At $z \sim 20$, $M_J \sim 10^4\,M_\odot$, much smaller than the minimum halo mass for star formation because additional conditions (cooling, angular momentum) must be satisfied.

3.2 Minimum Halo Mass for H₂ Cooling

For gas to collapse and form stars, it must be able to radiate away its gravitational energy. The virial temperature of a dark matter halo of mass $M$ collapsing at redshift $z$ is:

Virial Temperature

$$T_{\text{vir}} = \frac{\mu\,m_p}{2\,k_B}\,\frac{GM}{r_{\text{vir}}} \approx 1800\;\text{K}\,\left(\frac{M}{10^6\,M_\odot}\right)^{2/3}\left(\frac{1+z}{21}\right)$$

where $r_{\text{vir}}$ is the virial radius defined by the overdensity $\Delta_c \approx 200$ relative to the critical density.

Molecular hydrogen cooling is effective for gas temperatures in the range $200\;\text{K} \lesssim T \lesssim 10^4\;\text{K}$. Requiring $T_{\text{vir}} \gtrsim 200\text{--}1000\;\text{K}$ for $\text{H}_2$ cooling to operate yields the minimum halo mass:

Minimum Halo Mass for First Stars

$$M_{\min} \sim 10^5\text{--}10^6\,M_\odot \quad \text{at } z \sim 20\text{--}30$$

These are the so-called “minihalos.” At lower redshifts, atomic cooling halos with $T_{\text{vir}} \gtrsim 10^4\;\text{K}$ and $M \gtrsim 10^8\,M_\odot$ can cool via Lyman-α emission and are more robust to feedback.

3.3 Population III Star Properties

In the absence of metals, the gas cannot fragment as efficiently as in present-day molecular clouds. The characteristic mass scale is set by the Jeans mass at the loitering point where $\text{H}_2$ cooling becomes inefficient ($T \sim 200\;\text{K}$, $n \sim 10^4\;\text{cm}^{-3}$):

$$M_{J,\text{frag}} \sim 10^3\,M_\odot \left(\frac{T}{200\;\text{K}}\right)^{3/2}\left(\frac{n}{10^4\;\text{cm}^{-3}}\right)^{-1/2}$$

Hydrodynamic simulations show that the protostellar accretion rate is very high, leading to massive Pop III stars with characteristic properties:

Population III Star Characteristics

Mass: $M_\star \sim 10\text{--}1000\,M_\odot$, with a typical mass $\sim 100\,M_\odot$. Some simulations suggest a bimodal IMF with lower-mass fragments.
Effective temperature: $T_{\text{eff}} \sim 10^5\;\text{K}$ for a $100\,M_\odot$ zero-metallicity star, much hotter than present-day massive stars ($T_{\text{eff}} \sim 4 \times 10^4\;\text{K}$).
Luminosity: $L \sim 10^6\text{--}10^7\,L_\odot$, radiating predominantly in the far-UV.
UV emission: The high surface temperature means a large fraction of photons have $E > 13.6\;\text{eV}$ (ionizing), with ionizing photon production rate $\dot{N}_\gamma \sim 10^{50}\;\text{s}^{-1}$ per star.
Lifetime: $\sim 2\text{--}3\;\text{Myr}$, ending as pair-instability supernovae ($140\text{--}260\,M_\odot$) or direct-collapse black holes ($M > 260\,M_\odot$).

The copious UV emission from Pop III stars has two critical effects: (1) producing a Lyman-α background that couples $T_S$ to $T_K$ via the Wouthuysen-Field effect, making the 21-cm signal visible, and (2) beginning to ionize the surrounding IGM, initiating the epoch of reionization.

4. Reionization Physics

Reionization is fundamentally a competition between the production of ionizing photons by luminous sources and the recombination of hydrogen ions in the IGM. We derive the key equations governing this process.

4.1 Ionization Balance Equation

The volume-averaged ionization fraction $Q_{\text{HII}}$ (the filling factor of ionized regions) evolves according to:

Reionization Equation

$$\frac{dQ_{\text{HII}}}{dt} = \frac{\dot{n}_{\text{ion}}}{\bar{n}_H} - \frac{Q_{\text{HII}}}{t_{\text{rec}}}$$

The first term represents ionization by sources; the second represents recombinations in the IGM.

The ionizing photon production rate per unit volume is:

$$\dot{n}_{\text{ion}} = f_{\text{esc}}\,\xi_{\text{ion}}\,\rho_{\text{SFR}}$$

where $f_{\text{esc}}$ is the escape fraction of ionizing photons, $\xi_{\text{ion}}$ is the ionizing photon production efficiency (photons per unit stellar mass formed), and $\rho_{\text{SFR}}$ is the star formation rate density.

4.2 Recombination Timescale

The recombination timescale in the clumpy IGM is:

Recombination Timescale

$$t_{\text{rec}} = \frac{1}{\alpha_B\,n_e\,C}$$

where $\alpha_B \approx 2.6 \times 10^{-13}\;\text{cm}^3\,\text{s}^{-1}$ is the case-B recombination coefficient at $T = 10^4\;\text{K}$, $n_e \approx \bar{n}_H(1+z)^3$ is the electron density, and $C = \langle n^2 \rangle / \langle n \rangle^2$ is the clumping factor.

The clumping factor accounts for the enhanced recombination rate in overdense regions. Simulations give $C \sim 2\text{--}5$ during reionization. Evaluating the recombination time:

$$t_{\text{rec}} \approx 0.93\;\text{Gyr}\;\left(\frac{C}{3}\right)^{-1}\left(\frac{1+z}{8}\right)^{-3}\left(\frac{T}{2\times 10^4\;\text{K}}\right)^{0.7}$$

4.3 Reionization Criterion

Reionization is complete when $Q_{\text{HII}} = 1$. A necessary condition is that the cumulative number of ionizing photons produced per hydrogen atom exceeds unity plus the number of recombinations:

Reionization Criterion

$$\frac{\dot{n}_{\text{ion}}}{\bar{n}_H} > \frac{1}{t_{\text{rec}}} = \alpha_B\,\bar{n}_H\,(1+z)^3\,C$$

The ionizing photon production rate per atom must exceed the recombination rate. This is often expressed as a critical star formation rate density needed to maintain reionization.

4.4 Thomson Optical Depth

The reionization history is constrained by the Thomson scattering optical depth of the CMB, which measures the integrated column of free electrons along the line of sight:

Thomson Optical Depth

$$\tau_e = \int_0^{z_{\text{re}}} n_e(z)\,\sigma_T\,c\,\frac{dt}{dz}\,dz$$

where $\sigma_T = 6.652 \times 10^{-25}\;\text{cm}^2$ is the Thomson cross-section and $dt/dz = -1/[(1+z)H(z)]$.

For a simple model of instantaneous reionization at redshift $z_{\text{re}}$, the electron density is $n_e = \bar{n}_H(1+z)^3$ for $z < z_{\text{re}}$. Evaluating the integral in the matter-dominated era ($\Omega_m = 1$ approximation):

$$\tau_e \approx \frac{2\,\bar{n}_{H,0}\,\sigma_T\,c}{3\,H_0\,\Omega_m^{1/2}}\left[(1+z_{\text{re}})^{3/2} - 1\right]$$

Planck (2018) measures $\tau_e = 0.054 \pm 0.007$, implying an effective instantaneous reionization redshift of $z_{\text{re}} \approx 7.7 \pm 0.7$.

The relatively low value of $\tau_e$ rules out very early reionization models and indicates that the process occurred primarily at $z \lesssim 10$. For a more realistic extended reionization history, $\tau_e$ is computed numerically by integrating over the ionization fraction history $x_e(z)$.

5. The Gunn-Peterson Trough

The Gunn-Peterson (GP) effect is the most direct observational probe of neutral hydrogen in the high-redshift IGM. Lyman-α photons ($\lambda_\alpha = 1216$ Å) from a background quasar are resonantly absorbed by neutral hydrogen along the line of sight.

5.1 Gunn-Peterson Optical Depth

The GP optical depth at observed wavelength $\lambda_{\text{obs}} = \lambda_\alpha(1+z)$ is:

Gunn-Peterson Optical Depth

$$\tau_{\text{GP}} = \frac{\pi\,e^2}{m_e\,c}\,\frac{f_\alpha\,n_{\text{HI}}}{H(z)}$$

where $f_\alpha = 0.4162$ is the Lyman-α oscillator strength and $n_{\text{HI}} = x_{\text{HI}}\,\bar{n}_H(1+z)^3$ is the neutral hydrogen density.

Evaluating this expression numerically:

$$\tau_{\text{GP}} \approx 4.9 \times 10^5\,x_{\text{HI}}\left(\frac{\Omega_b h^2}{0.023}\right)\left(\frac{\Omega_m h^2}{0.15}\right)^{-1/2}\left(\frac{1+z}{7}\right)^{3/2}$$

5.2 Implications for Reionization

The enormous prefactor ($\sim 5 \times 10^5$) means that even a tiny neutral fraction of $x_{\text{HI}} > 10^{-4}$ produces $\tau_{\text{GP}} \gg 1$, completely absorbing the quasar flux at the Lyman-α wavelength. This has profound consequences:

Complete GP absorption: For $x_{\text{HI}} \gtrsim 10^{-4}$, the transmitted flux $e^{-\tau_{\text{GP}}} \approx 0$. The GP trough is essentially a binary indicator: it signals $x_{\text{HI}} > 10^{-4}$ but cannot measure higher neutral fractions precisely.
Quasar observations: The Sloan Digital Sky Survey (SDSS) and subsequent surveys discovered quasars at $z > 6$ showing complete GP absorption (zero transmitted flux blueward of Lyman-α), while quasars at $z \lesssim 5.5$ show a transmitted Lyman-α forest. This transition marks the end of reionization.
Fan et al. (2006): Analysis of 19 SDSS quasars at $5.7 < z < 6.4$ showed a rapid increase in the effective GP optical depth from $\tau_{\text{eff}} \sim 2$ at $z = 5.5$ to $\tau_{\text{eff}} > 5$ at $z > 6$, indicating that the volume-averaged neutral fraction rises from $\sim 10^{-4}$ to $\gtrsim 10^{-3}$ over this interval.

The GP damping wing — the broad absorption feature from the natural Lorentzian profile of the Lyman-α line — provides sensitivity to higher neutral fractions ($x_{\text{HI}} \gtrsim 0.1$) and has been used to infer $x_{\text{HI}} \sim 0.1\text{--}0.6$ at $z \sim 7\text{--}8$ from quasar and gamma-ray burst spectra.

6. Observational Probes and Applications

6.1 EDGES 21-cm Detection

The Experiment to Detect the Global Epoch of Reionization Signature (EDGES) reported a detection of the sky-averaged 21-cm absorption signal centered at $\nu = 78 \pm 1\;\text{MHz}$, corresponding to $z \approx 17$. The measured absorption profile had an amplitude of approximately $T_b \approx -500^{+200}_{-500}\;\text{mK}$, which is roughly twice the depth predicted by standard models ($\sim -200\;\text{mK}$).

If confirmed, the excess absorption could indicate either that the gas was colder than expected (perhaps due to baryon-dark matter interactions as proposed by Barkana 2018) or that the background radiation temperature was enhanced beyond the CMB alone (excess radio background). The SARAS3 experiment has placed constraints that are in tension with the EDGES profile, and the result remains under active investigation.

6.2 JWST Early Galaxies

The James Webb Space Telescope has discovered an unexpectedly large number of luminous galaxies at $z > 10$, with spectroscopic confirmations reaching $z \approx 14$. These observations suggest that star formation efficiency in the early universe may have been higher than pre-JWST models predicted. The UV luminosity density at $z \sim 10\text{--}14$ appears sufficient to begin reionization early, consistent with extended reionization models.

JWST spectroscopy has also measured the ionizing photon production efficiency $\xi_{\text{ion}}$ and Lyman-α emission from galaxies within the EoR, providing direct constraints on the sources responsible for reionization.

6.3 Planck Thomson Optical Depth

The Planck satellite measured $\tau_e = 0.054 \pm 0.007$ from the large-scale polarization of the CMB (the E-mode signal at low multipoles). This value is sensitive to the integrated ionization history and constrains the midpoint of reionization to $z_{\text{mid}} = 7.7 \pm 0.7$ for a hyperbolic tangent reionization model.

The relatively low $\tau_e$ disfavors models with significant ionization at $z > 15$ and is consistent with reionization driven primarily by galaxies at $z \sim 6\text{--}10$.

6.4 21-cm Intensity Mapping

Future 21-cm experiments aim to map the three-dimensional distribution of neutral hydrogen during the EoR. The power spectrum of 21-cm brightness temperature fluctuations encodes information about the sizes and distribution of ionized bubbles:

$$P_{21}(k, z) = \bar{T}_b^2\left[\bar{x}_{\text{HI}}^2\,P_{\delta\delta}(k) + 2\bar{x}_{\text{HI}}\,P_{\delta x}(k) + P_{xx}(k)\right]$$

where $P_{\delta\delta}$ is the matter power spectrum, $P_{xx}$ is the ionization power spectrum, and $P_{\delta x}$ is the cross-spectrum. Experiments such as HERA (Hydrogen Epoch of Reionization Array) and the upcoming SKA-Low are designed to measure this signal.

7. Historical Context and Key References

Foundational Papers

Gunn, J.E. & Peterson, B.A. (1965). “On the Density of Neutral Hydrogen in Intergalactic Space,” Astrophysical Journal 142, 1633–1641. — Prediction of the Gunn-Peterson effect: complete Lyman-α absorption by neutral IGM hydrogen.
Field, G.B. (1958). “Excitation of the Hydrogen 21-cm Line,” Proceedings of the IRE 46, 240–250. — Theory of the 21-cm spin temperature and its coupling to the radiation field and gas kinetic temperature.
Wouthuysen, S.A. (1952). “On the Excitation Mechanism of the 21-cm Interstellar Hydrogen Emission Line,” Astronomical Journal 57, 31–32. — Prediction of the Lyman-α pumping mechanism (Wouthuysen-Field effect).
Madau, P., Meiksin, A. & Rees, M.J. (1997). “21 Centimeter Tomography of the Intergalactic Medium at High Redshift,” Astrophysical Journal 475, 429–444. — Foundational paper on 21-cm cosmology.

Observational Milestones

Fan, X. et al. (2006). “Constraining the Evolution of the Ionizing Background and the Epoch of Reionization with z ~ 6 Quasars. II. A Sample of 19 Quasars,” Astronomical Journal 132, 117–136. — SDSS quasar spectra revealing Gunn-Peterson troughs at $z > 6$.
Bowman, J.D. et al. (2018). “An Absorption Profile Centred at 78 Megahertz in the Sky-Averaged Spectrum,” Nature 555, 67–70. — EDGES detection of the global 21-cm absorption signal at $z \approx 17$.
Planck Collaboration (2020). “Planck 2018 results. VI. Cosmological parameters,” Astronomy & Astrophysics 641, A6. — Thomson optical depth measurement $\tau_e = 0.054 \pm 0.007$.
Robertson, B.E. et al. (2023). “Identification and properties of intense star-forming galaxies at redshifts z > 10,” Nature Astronomy 7, 611–621. — JWST discovery of luminous galaxies at $z > 10$ and implications for early reionization.
Abel, T., Bryan, G.L. & Norman, M.L. (2002). “The Formation of the First Star in the Universe,” Science 295, 93–98. — First ab initio simulation of Population III star formation in a cosmological context.

8. Python Simulation: 21-cm Signal and Reionization History

The following simulation computes the expected global 21-cm brightness temperature signal as a function of redshift, capturing the key features: the dark-age absorption trough, the cosmic dawn absorption dip (enhanced by Wouthuysen-Field coupling), X-ray heating turning absorption to emission, and the reionization cutoff. It also computes and plots the ionization fraction history $x_{\text{HII}}(z)$.

Python

script.py284 lines

#!/usr/bin/env python3
"""
dark_ages_reionization.py
Compute the global 21-cm brightness temperature signal and ionization
fraction history through the dark ages and epoch of reionization.

Uses only numpy. No scipy.
"""
import numpy as np

# ============================
# Physical Constants (CGS)
# ============================
c       = 2.998e10       # speed of light [cm/s]
k_B     = 1.381e-16      # Boltzmann constant [erg/K]
h_P     = 6.626e-27      # Planck constant [erg s]
m_p     = 1.673e-24      # proton mass [g]
m_e     = 9.109e-28      # electron mass [g]
sigma_T = 6.652e-25      # Thomson cross-section [cm^2]
G       = 6.674e-8       # gravitational constant [cgs]
eV_erg  = 1.602e-12      # 1 eV in erg
Mpc_cm  = 3.086e24       # 1 Mpc in cm
yr_s    = 3.156e7        # 1 year in seconds
Gyr_s   = yr_s * 1e9

# ============================
# Cosmological Parameters
# ============================
H0_km   = 67.4           # Hubble constant [km/s/Mpc]
H0      = H0_km * 1e5 / Mpc_cm  # [s^-1]
Omega_m = 0.315
Omega_b = 0.0493
Omega_L = 1.0 - Omega_m
Omega_r = 9.15e-5
T_CMB0  = 2.7255         # CMB temperature today [K]
Y_He    = 0.245          # helium mass fraction
n_H0    = (1 - Y_He) * Omega_b * 3 * H0**2 / (8 * np.pi * G * m_p)

# 21-cm constants
nu_21   = 1.4204e9       # 21-cm frequency [Hz]
T_star  = h_P * nu_21 / k_B   # 0.068 K
A_10    = 2.869e-15      # Einstein A coefficient [s^-1]

# ============================
# Helper Functions
# ============================
def T_CMB(z):
    return T_CMB0 * (1.0 + z)

def Hz(z):
    return H0 * np.sqrt(Omega_r*(1+z)**4 + Omega_m*(1+z)**3 + Omega_L)

def n_H(z):
    return n_H0 * (1.0 + z)**3

def dt_dz(z):
    return -1.0 / ((1.0 + z) * Hz(z))

# ============================
# Gas Temperature Model
# ============================
# Compton decoupling at z_dec ~ 150
# Above z_dec: T_K = T_CMB (Compton coupling)
# Below z_dec: T_K ~ T_CMB(z_dec) * ((1+z)/(1+z_dec))^2
z_dec = 150.0

def T_gas(z):
    if z > z_dec:
        return T_CMB(z)
    else:
        return T_CMB(z_dec) * ((1.0 + z) / (1.0 + z_dec))**2

# ============================
# Spin Temperature Model
# ============================
# Coupling coefficients model the physics of the spin temperature
# x_c: collisional coupling (important at high z in dark ages)
# x_alpha: Wouthuysen-Field coupling (turns on with first stars)

def x_collisional(z):
    """Collisional coupling coefficient (approximate)."""
    T_K = T_gas(z)
    nH = n_H(z)
    # Collisional de-excitation rate ~ 3.1e-11 T^0.357 exp(-32/T) cm^3/s
    # for H-H collisions (Zygelman 2005 approximation)
    kappa_HH = 3.1e-11 * T_K**0.357 * np.exp(-32.0 / T_K)
    x_c = (4.0/3.0) * kappa_HH * nH * T_star / (A_10 * T_CMB(z))
    return x_c

def x_wf(z):
    """Wouthuysen-Field coupling coefficient (phenomenological model).
    Turns on around z ~ 25 as first stars produce Lyman-alpha photons."""
    # Model: x_alpha grows exponentially as stars form, peaks around z~15
    if z > 30:
        return 0.0
    elif z > 6:
        # Smooth turn-on of Lyman-alpha coupling
        z_onset = 28.0
        z_peak  = 12.0
        if z > z_onset:
            return 0.0
        # Ramp up
        frac = (z_onset - z) / (z_onset - z_peak)
        frac = min(frac, 1.0)
        x_a = 100.0 * frac**3
        return x_a
    else:
        return 100.0  # strongly coupled during and after reionization

def T_spin(z):
    """Compute spin temperature from coupling coefficients."""
    Tcmb = T_CMB(z)
    Tk   = T_gas(z)
    xc   = x_collisional(z)
    xa   = x_wf(z)

# T_alpha ~ T_K for the WF effect
    T_alpha = Tk

numerator   = Tcmb**(-1) + xc * Tk**(-1) + xa * T_alpha**(-1)
    denominator = 1.0 + xc + xa

Ts_inv = numerator / denominator
    return 1.0 / Ts_inv

# ============================
# X-ray Heating Model
# ============================
# X-rays from first accreting black holes and HMXBs heat the IGM
# Model: heating turns on around z ~ 20, heats gas above T_CMB by z ~ 12

def T_gas_heated(z):
    """Gas temperature including X-ray heating."""
    T_adiabatic = T_gas(z)
    if z > 22:
        return T_adiabatic
    elif z > 8:
        # X-ray heating ramps up
        z_heat_start = 22.0
        z_heat_sat   = 10.0
        frac = (z_heat_start - z) / (z_heat_start - z_heat_sat)
        frac = max(0.0, min(frac, 1.0))
        # Heating raises T_K; saturates well above T_CMB
        T_heated = T_adiabatic + frac**2 * 5000.0  # X-ray heated to ~5000 K
        return T_heated
    else:
        return T_adiabatic + 5000.0

def T_spin_full(z):
    """Spin temperature with X-ray heating included."""
    Tcmb = T_CMB(z)
    Tk   = T_gas_heated(z)
    xc   = x_collisional(z)
    xa   = x_wf(z)
    T_alpha = Tk

numerator   = Tcmb**(-1) + xc * Tk**(-1) + xa * T_alpha**(-1)
    denominator = 1.0 + xc + xa
    Ts_inv = numerator / denominator
    return 1.0 / Ts_inv

# ============================
# Ionization Fraction Model
# ============================
# Simple tanh reionization model consistent with Planck tau_e
z_re   = 7.7     # midpoint of reionization
dz_re  = 1.5     # width parameter

def x_HII(z):
    """Volume-averaged ionization fraction (tanh model)."""
    if z > 35:
        return 1.2e-4  # residual from recombination
    arg = (z - z_re) / dz_re
    Q = 0.5 * (1.0 - np.tanh(arg))
    # Ensure residual ionization at high z
    return max(Q, 1.2e-4)

def x_HI(z):
    """Neutral fraction."""
    return 1.0 - x_HII(z)

# ============================
# 21-cm Brightness Temperature
# ============================
def T_brightness(z):
    """Global 21-cm brightness temperature [mK]."""
    Ts   = T_spin_full(z)
    Tcmb = T_CMB(z)
    xHI  = x_HI(z)

# Optical depth
    nHI  = n_H(z) * xHI
    tau_21 = (3.0 * c**3 * h_P * A_10 * nHI) / (16.0 * k_B * nu_21**2 * Ts * Hz(z))

# Brightness temperature in mK
    Tb = (Ts - Tcmb) / (1.0 + z) * (1.0 - np.exp(-tau_21)) * 1000.0
    return Tb

# ============================
# Thomson Optical Depth
# ============================
def compute_tau_e(z_arr):
    """Compute cumulative Thomson optical depth."""
    tau = np.zeros(len(z_arr))
    for i in range(1, len(z_arr)):
        z_mid = 0.5 * (z_arr[i-1] + z_arr[i])
        dz = z_arr[i] - z_arr[i-1]
        ne = x_HII(z_mid) * n_H(z_mid)
        dtdz = abs(dt_dz(z_mid))
        tau[i] = tau[i-1] + ne * sigma_T * c * dtdz * dz
    return tau

# ============================
# Compute and Print Results
# ============================
# Redshift array from z=5 to z=50
z_vals = np.linspace(5.0, 50.0, 500)

# Compute 21-cm brightness temperature
Tb_vals = np.array([T_brightness(z) for z in z_vals])

# Compute ionization fraction
xHII_vals = np.array([x_HII(z) for z in z_vals])

# Compute Thomson optical depth (integrate from z=0)
z_tau = np.linspace(0.0, 30.0, 3000)
tau_vals = compute_tau_e(z_tau)

# Print key results
print("=" * 60)
print("DARK AGES & REIONIZATION: GLOBAL 21-CM SIGNAL")
print("=" * 60)

# Find absorption minimum
idx_min = np.argmin(Tb_vals)
print(f"Absorption minimum: T_b = {Tb_vals[idx_min]:.1f} mK at z = {z_vals[idx_min]:.1f}")

# Find emission maximum
idx_max = np.argmax(Tb_vals)
print(f"Emission maximum:   T_b = {Tb_vals[idx_max]:.1f} mK at z = {z_vals[idx_max]:.1f}")

# Thomson optical depth
print(f"Thomson optical depth (z=0 to 30): tau_e = {tau_vals[-1]:.4f}")
print(f"  (Planck 2018 measurement: tau_e = 0.054 +/- 0.007)")

# Reionization milestones
for frac in [0.1, 0.5, 0.9, 0.99]:
    for i, z in enumerate(z_vals):
        if x_HII(z) >= frac and z_vals[min(i+1, len(z_vals)-1)] < z:
            continue
        if x_HII(z) <= frac:
            print(f"x_HII = {frac:.2f} at z ~ {z:.1f}")
            break

print()
print("=" * 60)
print("21-cm Brightness Temperature vs Redshift")
print("=" * 60)
print(f"{'z':>8s}  {'T_b [mK]':>12s}  {'x_HI':>8s}  {'T_S [K]':>10s}  {'T_CMB [K]':>10s}")
print("-" * 60)
for i in range(0, len(z_vals), 25):
    z = z_vals[i]
    print(f"{z:8.1f}  {Tb_vals[i]:12.2f}  {x_HI(z):8.4f}  {T_spin_full(z):10.2f}  {T_CMB(z):10.2f}")

print()
print("=" * 60)
print("Ionization Fraction History")
print("=" * 60)
print(f"{'z':>8s}  {'x_HII':>10s}")
print("-" * 25)
z_ion_arr = np.linspace(5.0, 20.0, 30)
for z in z_ion_arr:
    print(f"{z:8.2f}  {x_HII(z):10.6f}")

print()
print("=" * 60)
print("Thomson Optical Depth Buildup")
print("=" * 60)
print(f"{'z':>8s}  {'tau_e':>12s}")
print("-" * 25)
for i in range(0, len(z_tau), 200):
    print(f"{z_tau[i]:8.2f}  {tau_vals[i]:12.6f}")
print(f"{z_tau[-1]:8.2f}  {tau_vals[-1]:12.6f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Thermal History Cosmic Inflation →