Function ionization_fraction

pub fn ionization_fraction(z: f64, cosmo: &Cosmology) -> f64

Ionization fraction X_e(z) with Peebles TLA correction.

Returns the total free electron fraction (hydrogen + helium) per hydrogen atom.

Regimes

• z > 8000: Fully ionized H; He from Saha equations.
• z_switch < z ≤ 8000: Saha equilibrium for H + He.
• z ≤ z_switch: Peebles three-level atom ODE for H, plus Saha He.

Saha-subtracted ODE

The raw Peebles ODE has catastrophic cancellation: α_B n_H X_e² and β_B (1−X_e) are both ~10² s⁻¹ but their difference is ~10⁻⁴. Using the Saha relation β_B = α_B X_S² n_H / (1−X_S), we rewrite:

dX_e/dz = C × α_B × n_H / (H(1+z)) × [X_e² − X_S² (1−X_e)/(1−X_S)]

This is O(X_e − X_S) near equilibrium, eliminating the cancellation.

References:

• Peebles (1968) — Three-level atom
• Seager, Sasselov & Scott (1999) — RECFAST
• Liu et al. (2020) — DarkHistory implementation