Thermal history: recombination and reionization

No reionization at all. Useful for isolating the recombination history.

Injection(; p_ann = 0, Γ_dec = 0, f_dec = 0, f_eff = 1, f_halo = 0, z_halo = 30)

Exotic energy injection into the pre-reionization gas: dark-matter annihilation and/or decay, in the same parameterisation as CLASS's injection module.

• p_ann — annihilation parameter ⟨σv⟩/mDM in m³ s⁻¹ kg⁻¹. The injected power is dE/dV dt = ρcdm² c² pann, so it scales as (1+z)⁶ and matters earliest. (CLASS's `DMannihilation_efficiency` is this divided by c².)
• f_halo, z_halo — the erfc((1+z)/(1+z_halo))/(1+z)³ boost once structure forms.
• Γ_dec, f_dec — decay rate (1/s) and the decaying mass fraction of the DM. The power is ρcdm c² fdec Γ e^{−Γt}, scaling as (1+z)³ and mattering latest.
• f_eff — deposition efficiency: the fraction of injected energy deposited into the gas at the same redshift ("on the spot"). The energy that does deposit is split between heating, hydrogen ionization and Lyman-α excitation with the Chen–Kamionkowski fractions χheat = (1+2xe)/3, χion = χLyα = (1−x_e)/3.

Injection deposits into three different equations at once — Tm directly, xH through both ionization and excitation — which is why it lives here and not in a bolt-on: the CMB constrains pann precisely because extra xe at z ~ 600 widens the visibility function.

RecombinationSolution

Interpolated thermal history. Call the accessors rather than touching fields: x_e, T_matter, n_e.

Reionization(; z_reio = 7.67, Δz = 0.5, ...)

The late-time reionization of the intergalactic medium by the first stars and quasars, in the tanh parameterisation used by CAMB and CLASS.

Recombination leaves the universe neutral with a residual xe ≈ 2e-4. It does not stay that way: by z ≈ 6 the IGM is ionized again, and those free electrons scatter CMB photons. About 5% of them do, which suppresses the CMB anisotropies by e^{-2τ} – roughly 10% in Cℓ – so this is not a detail that can be skipped before computing a power spectrum.

x_e(z) = (1 + f_He - x_rec) · [1 + tanh(y)]/2 + x_rec + (helium term)
y      = [(1 + z_reio)^1.5 - (1+z)^1.5] / [1.5 (1+z_reio)^0.5 Δz]

The tanh is taken in (1+z)^1.5 rather than in z because that makes the resulting optical depth nearly independent of the width Δz, so zreio and τ are close to interchangeable. The asymptote is 1 + fHe: hydrogen fully ionized and helium singly ionized. Helium's second ionization needs harder photons and happens later, around z ≈ 3.5, and is given its own tanh.

The default zreio = 7.67 reproduces Planck 2018's τ = 0.0544. Use [`zreiofromτ`](@ref) to solve for it given a target τ instead.

R_baryon(c, a)

Baryon-to-photon momentum density ratio, R = 3ρb/(4ργ). It is what makes the photon-baryon fluid heavy: the sound speed is c/√(3(1+R)), so R suppresses the oscillation frequency and is directly responsible for the alternating peak heights in the CMB.

T_matter(r, z)

Baryon temperature in K. Locked to T_cmb while Compton scattering is efficient, then falls as (1+z)² once the gas decouples thermally.

z_drag(r)

The baryon drag epoch, where the baryons stop being dragged by the photons. Defined by τdrag(zd) = 1 with

τ_drag(z) = ∫₀^z  (dτ/dz') / R  dz'

The 1/R weighting is the whole point, and it is why zdrag ≈ 1060 differs from z* ≈ 1090: what decouples a photon is one scattering, but to change a baryon's momentum takes R scatterings, since the baryons carry R times the momentum. In our universe R ≈ 0.6 at recombination, so the baryons are released slightly after the photons.

dτ_dz(r, z)

-dτ/dz, the differential scattering probability per unit redshift. Positive.

f_He(c)

Helium-to-hydrogen number ratio, nHe/nH = Yp / (3.9715 (1 - Yp)).

Yp is a mass fraction, so converting it to a number ratio divides by the actual He/H mass ratio – 3.9715 (CLASS's _not4_), not 4. The α particle's binding energy is 0.7% of its mass, and that 0.7% lands directly in the helium abundance recombination sees.

Injected power in J m⁻³ s⁻¹ at redshift z. t_s is cosmic time in seconds, needed only for the decay exponential (pass 0.0 when Γ_dec = 0).

n_H_of_z(c, z)

Number density of hydrogen nuclei (ionized or not), m⁻³.

Free-electron number density in m⁻³.

optical_depth(r, z)

Thomson optical depth back to redshift z,

τ(z) = ∫₀^z  n_e σ_T c  dz'/[(1+z') H(z')]

Dimensionless; τ(z_*) = 1 defines the surface of last scattering. Precomputed at construction, so this is an interpolation lookup.

Sound horizon at the drag epoch, in Mpc – the BAO standard ruler.

Sound horizon at last scattering, in Mpc.

recombination(c; z_start = 5000, z_end = 0, ...)

Integrate the recombination history and return a RecombinationSolution.

The state vector is (x_H, x_HeII, log T_m):

• x_H = nHII / nH, hydrogen ionization fraction
• x_HeII = nHeII / nHe, singly-ionized helium fraction
• T_m = matter temperature in K

with xe = xH + fHe · xHeII. Doubly-ionized helium is handled by Saha above z_start (it recombines by z ≈ 5000 and never leaves equilibrium).

Accuracy

This is RECFAST 1.5 (Wong, Moss & Scott 2008; Scott & Moss 2009), whose correction terms are calibrated against the full multi-level-atom codes – HyRec (Ali-Haïmoud & Hirata 2010) and CosmoRec – and against the Switzer & Hirata (2008) helium papers. Beyond the textbook Peebles treatment it carries:

• fudge_H = 1.125, retuned in 2012 to match HyRec/CosmoRec. A plain three-level hydrogen atom recombines measurably too slowly; this absorbs the leading effect of the ~100 bound levels a real atom has.
• gaussian_H, a double-Gaussian correction to the Lyman-α escape rate K. This is how the two-photon and radiative-transfer physics of the multi-level codes is folded into a two-level equation.
• Helium via Sobolev escape probability rather than a bare λ³/8πH, so the optically thin limit is right and not just the thick one.
• H I continuum opacity during helium recombination (Switzer & Hirata 2008 I; Kholupenko et al. 2007 fit). Neutral hydrogen is already present while helium is still recombining, and it absorbs the He I line photons that would otherwise have re-ionized helium. This accelerates He I recombination substantially and is the single largest helium correction.
• The 2³P₁ → 1¹S₀ intercombination line (triplets), a second, forbidden recombination pathway. Its Einstein coefficient is 177 s⁻¹ against the singlet's 1.8e9 s⁻¹ – ten million times weaker – but the singlet channel is so throttled by its own resonant opacity that the forbidden line still competes.

Set fudge_H = 1.14, gaussian_H = false, sobolev_He = false, triplets = false to recover the original RECFAST 1.0 of Seager, Sasselov & Scott.

The system is stiff: the reaction rates dwarf H at early times, which is precisely why Saha holds there. It needs an implicit solver.

recombination_optical_depth(r, z)

Optical depth back to z counting only the electrons left free by recombination, excluding those liberated again by reionization.

reionization_x_e(reio, c, z, x_rec)

Free-electron fraction including reionization, smoothly matched onto the residual x_rec left by recombination. At high z the tanh saturates at -1 and this returns x_rec exactly, so it can be applied at every redshift without a switch.

saha_x_e(c, z)

Free-electron fraction xe = ne/n_H from full Saha equilibrium, including singly and doubly ionized helium.

Three ionization balances are imposed at once,

n_e n_HII  / n_HI   = S_H
n_e n_HeII / n_HeI  = S_He1
n_e n_HeIII/ n_HeII = S_He2

each of the form S = g (2π me kT/h²)^{3/2} e^{-E/kT}, with statistical weights g = 1, 4, 1 respectively. They share the same ne, so the system is coupled and solved by root-finding on x_e rather than species by species.

Valid while the rates beat the expansion, i.e. above z ≈ 5000. Below that, real recombination lags equilibrium badly (Saha would have recombination finish far too early) and you want recombination.

Hydrogen ionization fraction from full Saha equilibrium.

Free-electron fraction from full Saha equilibrium (H, He I, He II).

sound_horizon(c, z)

Comoving distance a sound wave has travelled by redshift z, in Mpc:

r_s(z) = ∫_z^∞  c_s(z') dz' / H(z')  =  ∫₀^a  c_s da'/(a'² H)

Evaluated at the drag epoch this is the standard ruler: it sets the BAO scale and the CMB acoustic peak spacing.

Sound speed of the photon-baryon fluid, in units of c.

visibility(r, z)

The visibility function g(z) = (-dτ/dz)·e^{-τ}, normalised so ∫g dz = 1.

This is the probability density that a CMB photon we see today last scattered at redshift z. Its peak is the last-scattering surface, and its finite width is why that surface has a thickness rather than being a sharp shell – which in turn is what damps the high-ℓ acoustic peaks.

x_e(r, z)

Free-electron fraction xe = ne/n_H.

Above the ODE start redshift this falls back to full Saha equilibrium, which is exact there: the reaction rates beat the expansion by many orders of magnitude, so the gas has no choice but to sit in equilibrium. Extrapolating the interpolant instead would be both slower and wrong.

z_reio_from_τ(c, τ_target; kws...)

Invert τ(zreio) to find the reionization redshift giving a target optical depth. CAMB and CLASS both offer this, because τ is what the CMB measures and zreio is what the model is parameterised by.

z_star(r)

Redshift of last scattering, defined by τ(z_*) = 1.

Deliberately measured against the recombination optical depth, not the total. Reionization contributes τ ≈ 0.054 of its own, and counting it would push the τ = 1 crossing down to z ≈ 1085 – but those low-redshift scatterings are not what "last scattering" means. The surface that stamps the acoustic peaks onto the CMB is the recombination one, and this is the convention Planck's 1089.92 uses.

α_H(T)

Case-B recombination coefficient for hydrogen, m³/s. Pequignot, Petitjean & Boisson (1991) fit, as used by RECFAST.

"Case B" means recombinations to the ground state are excluded: those produce an ionizing photon that is immediately reabsorbed in the optically thick medium, so they achieve nothing.

α_He(T)

Case-B recombination coefficient for neutral helium, m³/s. Verner & Ferland (1996) fit, as used by RECFAST.

θ_star(r)

Angular size of the sound horizon at last scattering, θ* = rs(z*)/DM(z*). The single best-measured number in cosmology (Planck: 100θ* = 1.04109 ± 0.00030).

τ_dot(r, a)

dτ/dη = -a ne σT, in 1/Mpc. Negative because τ is measured backward from the observer. This is the quantity the Boltzmann hierarchy actually consumes: it is the photon scattering rate, and its magnitude relative to k sets whether the photons are tightly coupled to the baryons or free-streaming.

τ_reio(r)

Optical depth to reionization: the Thomson depth accumulated between us and z = 50, which is well above the reionization transition but far below recombination.

This is one of the six parameters of ΛCDM, and it is almost entirely degenerate with the amplitude As in the CMB – both scale Cℓ, and only the large-scale polarization breaks the degeneracy. Planck 2018: τ = 0.0544 ± 0.0073.

nuclear_deposition_rate(c, rec, z; history=nothing)

Electromagnetic power from primordial ³H and ⁷Be in J m⁻³ s⁻¹. Neutrino and antineutrino energy is absent by construction.

Fractional positive heating magnitude (dQ/dz)/ρ_γ from trace nuclei.

trace_nuclear_history(c, rec)

Solve the post-BBN ³H and ⁷Be histories. Returns nothing when c was built with a manual helium abundance and therefore carries no trace-nuclide solution.