diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 528b5ca..6448628 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -35,6 +35,13 @@ jobs:
run: python scripts/make_figures.py
env:
MPLBACKEND: Agg
+ - name: Run validation against CAMB
+ run: python scripts/validate.py
+ env:
+ MPLBACKEND: Agg
+ - name: Check committed results tables match the code
+ if: matrix.python-version == '3.13'
+ run: git diff --exit-code -- results/key_results.md results/validation.md
- name: Execute notebook top to bottom
run: jupyter nbconvert --to notebook --execute --output /tmp/executed.ipynb notebooks/recombination_and_decoupling.ipynb
@@ -54,3 +61,7 @@ jobs:
run: python scripts/make_figures.py
env:
MPLBACKEND: Agg
+ - name: Run validation against CAMB
+ run: python scripts/validate.py
+ env:
+ MPLBACKEND: Agg
diff --git a/CITATION.cff b/CITATION.cff
index 7979c3d..8655793 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -10,8 +10,8 @@ authors:
website: "https://sethuastro.github.io/"
repository-code: "https://github.com/sethuastro/CMB_Decoupling"
license: MIT
-version: 0.1.0
-date-released: 2026-10-22
+version: 1.0.0
+date-released: 2026-11-02
keywords:
- cosmology
- cosmic microwave background
diff --git a/README.md b/README.md
index 26f16f7..4b19593 100644
--- a/README.md
+++ b/README.md
@@ -3,104 +3,229 @@
[](https://github.com/sethuastro/CMB_Decoupling/actions/workflows/tests.yml)
[](LICENSE)
-When is the cosmic microwave background actually emitted, and how thick is the "surface" it
-comes from? This project is my attempt to answer that from first principles, by following the
-free-electron fraction of the Universe, $`x_e(z)`$, through hydrogen recombination.
-
-I started with the Saha equation, which is the textbook answer, and quickly found that it puts last
-scattering at $`z \approx 1291`$, about 200 in redshift too early. The reason is that recombination is
-a bottleneck rather than an equilibrium process. An electron captured straight into the ground state
-emits a photon that ionises the next atom, and Lyman-α photons are trapped by resonant scattering.
-So I implemented the Peebles three-level atom, in which hydrogen reaches the ground state only through
-two slow channels: Lyman-α photons redshifting out of the line, and the two-photon 2s→1s decay. From
-the resulting $`x_e(z)`$ I compute the Thomson scattering rate, the optical depth $\tau(z)$ and the
-visibility function $`g(z)`$, which is the probability distribution for when a CMB photon last scattered.
-
-Everything uses Planck 2018 parameters, and I check each step against CAMB, which uses the RECFAST
-recombination code.
+When is the cosmic microwave background actually emitted, and how thick is the "surface" it comes
+from? I built this code to learn recombination physics from first principles by following the
+free-electron fraction $`x_e(z)`$ through hydrogen recombination. I started with Saha equilibrium,
+which puts last scattering at $`z_\ast = 1291`$, far too early. I then implemented the Peebles
+three-level atom, and from its $`x_e(z)`$ I computed the Thomson optical depth $\tau(z)$ and the
+visibility function $`g(z)`$. I validated every step against CAMB, which uses the RECFAST
+recombination code, run on exactly the same Planck 2018 background. The plain Peebles model, with no
+fudge factor, reproduces CAMB's last-scattering redshift to 0.45 and the width of last scattering to
+0.67 in redshift. Its residual ionisation at $`z = 200`$ is 14.05% too high, and I explain why below.
+
+## Key results
- 
+ 
-
- 
-
-
-*Top: the free-electron fraction for each model, with the residual relative to CAMB underneath.
-Bottom: the visibility function (FWHM shaded) and the optical depth. Saha equilibrium recombines far
-too early; the Peebles model sits almost on top of CAMB through last scattering.*
-## Results
+*(a) Free-electron fraction from CAMB/RECFAST, my Peebles model and the Peebles model with the
+RECFAST fudge factor $`F = 1.14`$. (b) Fractional residual relative to CAMB; the grey band is the
+last-scattering window (the FWHM of CAMB's $`g(z)`$). (c) The visibility function.*
-Planck 2018 cosmology, no reionisation. "This work" is the plain Peebles model, with no fudge factor.
-I ran CAMB with exactly the same background, so the differences come from the recombination physics
-alone.
+Planck 2018 background, no reionisation. CAMB 2.0.5 is run with the same parameters (massless
+neutrinos, same $`Y_p`$ and $`T_0`$), so the differences come from the recombination physics alone.
+Nothing has been tuned to match CAMB.
-| | This work | CAMB 2.0.5 (RECFAST) | Difference |
+| | Peebles (this work) | Peebles, $`F = 1.14`$ (this work) | CAMB 2.0.5 (RECFAST) |
|---|---|---|---|
-| Last-scattering redshift $`z_\ast`$ ($`\tau = 1`$) | 1089.44 | 1089.89 | −0.45 |
-| Peak of $`g(z)`$, per unit redshift | 1078.80 | 1078.37 | +0.43 |
-| Peak of visibility, per unit conformal time | 1089.11 | 1088.78 | +0.32 |
-| Width of last scattering (FWHM of $`g`$) | 203.3 | 204.0 | −0.7 |
-| $`z_{\rm dec}`$, where $`\Gamma_T = H`$ | 903.2 | 903.6 | −0.4 |
-| Residual ionisation $`x_e(z = 200)`$ | 3.85 × 10⁻⁴ | 3.38 × 10⁻⁴ | +14% |
-| Residual ionisation $`x_e(z = 0)`$ | 2.67 × 10⁻⁴ | 1.74 × 10⁻⁴ | +53% |
-
-For comparison, Saha alone gives $`z_\ast = 1291`$.
-
-The last-scattering numbers agree with CAMB to better than one unit of redshift. That's because the
-Peebles model tracks RECFAST to within about 2% between $z \approx 1000$ and 1500, which is where
-$\tau$ builds up. The late-time residual ionisation is a different story: my model is 14% high at
-$`z = 200`$. I think this comes from two simplifications rather than a bug. First, a three-level atom
-recombines too slowly at late times, which is why RECFAST introduced a fudge factor; with their
-$F = 1.14$ my value drops to within 1% of CAMB. Second, I assume the gas stays at the radiation
-temperature, while in reality it cools faster below $`z \sim 300`$. Helium is the other visible gap.
-I treat it in Saha equilibrium, so at $`z \approx 1800\text{–}2800`$ my $`x_e`$ is up to 6% lower than CAMB's,
-where helium recombination is delayed.
-
-A detail worth pointing out: the familiar "$`z \approx 1090`$" is not where $g(z)$ peaks when you
-plot it per unit redshift. That peak is at 1079. The 1090 comes from the $\tau = 1$ definition, or
-from the visibility per unit conformal time. I report all three in the table, and
-[`docs/theory.md`](docs/theory.md#63-which-redshift-of-last-scattering) explains why they differ.
-
-To make sure my optical-depth integrator wasn't hiding errors, I also fed CAMB's own $`x_e(z)`$ into
-it. It recovers CAMB's internal $`z_\ast = 1089.892`$ to within 0.01.
-
-## How it's organised
+| $`z_\ast`$, where $\tau = 1$ | 1089.44 | 1090.28 | 1089.89 |
+| Peak of $`g(z)`$ per unit redshift | 1078.80 | 1079.05 | 1078.37 |
+| Peak of visibility per unit conformal time | 1089.11 | 1089.32 | 1088.78 |
+| Width of last scattering (FWHM of $`g`$) | 203.33 | 202.73 | 204.00 |
+| $`z_{\rm dec}`$, where $`\Gamma_T = H`$ | 903.18 | 904.44 | 903.55 |
+| $`x_e(z = 200)`$ | 3.849 × 10⁻⁴ | 3.404 × 10⁻⁴ | 3.375 × 10⁻⁴ |
+| $`x_e(z = 0)`$ | 2.665 × 10⁻⁴ | 2.351 × 10⁻⁴ | 1.738 × 10⁻⁴ |
+
+**How far $`x_e`$ is from CAMB.** Fractional residual $`\Delta x_e/x_e = x_e/x_e^{\rm CAMB} - 1`$
+on $`z = 200\text{–}2000`$, sampled every $\Delta z = 0.1$:
+
+| | Peebles (this work) | Peebles, $`F = 1.14`$ (this work) |
+|---|---|---|
+| (a) Largest $`\lvert\Delta x_e/x_e\rvert`$ in the last-scattering window, $`z = 966.5\text{–}1170.5`$ | 1.08% (at $`z = 1170.5`$) | 1.04% (at $`z = 1170.5`$) |
+| (b) Most negative residual | −3.58% (at $`z = 2000`$) | −3.58% (at $`z = 2000`$) |
+| (c) Most positive residual | +14.05% (at $`z = 200`$) | +1.46% (at $`z = 1264`$) |
+| Range over the hydrogen-only era, $`z = 200\text{–}1800`$ | −0.33% to +14.05% | −1.46% to +1.46% |
+
+**(a) Last scattering.** This is the part that matters for the CMB. Across the window that holds half
+the visibility, my $`x_e`$ is within 1.08% of CAMB's. That is why $`z_\ast`$ and the width agree to
+better than one unit of redshift. The largest value sits at the window's upper edge, because just
+above the window both curves in panel (b) rise to a bump, which reaches +1.46% at $`z = 1264`$ for
+the $`F = 1.14`$ run.
+CAMB's RECFAST has two corrections that my model lacks: its own fudge factor $`F = 1.125`$, and two
+Gaussian corrections to the recombination rate centred at $`z \approx 1452`$ and 867 (the
+`RECFAST_Hswitch` option, fitted to multi-level calculations). I think these are the likely origin of
+the bump, but I have not isolated it.
+
+**(b) Helium.** Above $`z = 1800`$ I treat helium in Saha equilibrium, while RECFAST evolves helium
+recombination with a rate equation. Real helium recombination lags Saha, so my $`x_e`$ is too low.
+The −3.58% minimum on this range sits at its upper edge, $`z = 2000`$. Beyond that range the gap
+reaches −5.81% at $`z = 2299`$. Helium has recombined well before last scattering, so this does not
+affect $`z_\ast`$.
+
+**(c) Late times.** Below the last-scattering window the plain Peebles residual grows steadily,
+reaching +14.05% at $`z = 200`$. The largest value over the whole range therefore sits at
+the lower end of the range I chose. There, $`C_r \to 1`$ and the freeze-out value of $`x_e`$ is set by
+the recombination rate $`\alpha_B`$. A three-level atom with the case-B rate recombines too slowly
+compared with a full multi-level atom. Seager, Sasselov & Scott (1999) introduced the fudge factor
+$`F = 1.14`$ on $`\alpha_B`$ to correct for exactly this. With $`F = 1.14`$, my residual stays within
+±1.46% over the whole hydrogen era, and the excess at $`z = 200`$ drops from +14.05% to +0.87%. So the
+late-time excess is the three-level approximation, not a numerical error.
+
+The comparison is not fully independent, because CAMB's RECFAST is itself a fudged three-level model.
+The fudge also moves $`z_\ast`$ slightly *away* from CAMB, from 1089.44 to 1090.28, which is why I
+keep the plain Peebles model as my headline result. The small positive residual that remains at the lowest
+redshifts with $`F = 1.14`$ (+0.87% at $`z = 200`$) is what I would expect from my second simplification. I keep the
+gas at the radiation temperature, while RECFAST lets it cool faster, and cooler gas recombines faster.
+By $`z = 0`$ this matters a lot: even with $`F = 1.14`$, my $`x_e(0)`$ is 2.351 × 10⁻⁴ against CAMB's
+1.738 × 10⁻⁴.
+
+**Two more checks.**
+- I fed CAMB's own $`x_e(z)`$ into my optical-depth integrator. It gives $`z_\ast = 1089.8901`$,
+ against CAMB's internal 1089.8923, so the integrator is not hiding anything.
+- The familiar "$`z \approx 1090`$" is $`z_\ast`$ or the peak per unit conformal time. The peak of
+ $`g(z)`$ per unit *redshift* is about 10 lower, at 1078.80.
+ [`docs/theory.md`](docs/theory.md#63-which-redshift-of-last-scattering) explains why.
+
+Every number above is printed by `python scripts/validate.py` (written to `results/validation.md`) or
+`python scripts/make_figures.py` (written to `results/key_results.md`).
+
+## Physics
+
+The full derivations, with each approximation stated where it is made, are in
+[`docs/theory.md`](docs/theory.md). These are the equations the code solves.
+
+**Background.** Flat ΛCDM with $`H(z) = H_0\sqrt{\Omega_m(1+z)^3 + \Omega_r(1+z)^4 + \Omega_\Lambda}`$,
+radiation temperature $`T = T_0(1+z)`$ and hydrogen density $`n_{\rm H} = (1 - Y_p)\,\rho_b/m_{\rm H}`$.
+Here $`Y_p`$ is the helium mass fraction, $`\rho_b`$ the baryon density and $`m_{\rm H}`$ the hydrogen
+atom mass. The parameters are Planck 2018 VI, Table 2 (TT,TE,EE+lowE+lensing).
+
+### Saha equilibrium, with helium
+
+If ionisation and recombination are in balance, the ionised fraction $`x_e = n_e/n_{\rm H}`$ of pure
+hydrogen obeys
+
+```math
+\frac{x_e^2}{1 - x_e} = \frac{1}{n_{\rm H}}\left(\frac{m_e k_B T}{2\pi\hbar^2}\right)^{3/2} e^{-E_{\rm ion}/k_BT},
+```
+
+where $`m_e`$ is the electron mass, $`k_B`$ Boltzmann's constant, $`\hbar`$ the reduced Planck constant
+and $`E_{\rm ion} = 13.598`$ eV the hydrogen ionisation energy. Helium adds two more equilibria:
+```math
+\frac{n_{\rm He\,II}\, n_e}{n_{\rm He\,I}} = 4\left(\frac{m_ek_BT}{2\pi\hbar^2}\right)^{3/2} e^{-\chi_{\rm He\,I}/k_BT},
+\qquad
+\frac{n_{\rm He\,III}\, n_e}{n_{\rm He\,II}} = \left(\frac{m_ek_BT}{2\pi\hbar^2}\right)^{3/2} e^{-\chi_{\rm He\,II}/k_BT}.
```
-src/cmb_decoupling/ the package
- cosmology.py Planck 2018 parameters (defined once), H(z), n_H(z)
- constants.py physical and atomic constants, with sources
- saha.py Saha equilibrium for hydrogen and helium
- peebles.py Peebles three-level atom (stiff ODE)
- visibility.py Thomson rate, tau(z), g(z), peak and width
- reference.py the CAMB comparison run
-notebooks/ a walk-through from Saha to last scattering
-scripts/ make_figures.py rebuilds every figure and the results table
-tests/ 38 tests: physical limits, convergence, comparison with CAMB
-docs/theory.md the derivations, step by step, with references
-archive/ my original scripts, kept as they were
+
+Here $`\chi_{\rm He\,I} = 24.587`$ eV and $`\chi_{\rm He\,II} = 54.418`$ eV are the helium ionisation
+energies, and 4 and 1 are statistical-weight factors. These are closed by charge neutrality:
+
+```math
+x_e = x_{\rm H\,II} + f_{\rm He}\,(y_{\rm He\,II} + 2y_{\rm He\,III}),
```
-The derivations in [`docs/theory.md`](docs/theory.md) go from the Boltzmann distribution to the Saha
-equation, then through the rate balance for the $n = 2$ level to the Peebles equation and the $`C_r`$
-factor, and finally to $\tau(z)$ and $`g(z)`$. I've tried to state each approximation at the point
-where it's made.
+with $`x_{\rm H\,II} = n_p/n_{\rm H}`$, $`y_i = n_i/n_{\rm He}`$ and $`f_{\rm He} = n_{\rm He}/n_{\rm H}`$.
+Saha fails during recombination for two reasons:
+- A direct capture to the ground state emits a photon that immediately ionises a neighbouring atom.
+- Lyman-α photons are trapped by resonant scattering.
+
+### The Peebles three-level atom
+
+Below $`z = 1800`$ the code integrates the Peebles (1968) equation. It treats hydrogen as a ground
+state, an $n = 2$ level and the continuum:
+
+```math
+\frac{dx_e}{dz} = \frac{C_r}{(1+z)H(z)}\left[\alpha_B\, n_{\rm H}\, x_e^2 - \beta_B\,(1 - x_e)\, e^{-E_{{\rm Ly}\alpha}/k_BT}\right],
+\qquad
+C_r = \frac{\Lambda_{2s1s} + \Lambda_\alpha}{\Lambda_{2s1s} + \Lambda_\alpha + \beta_B}.
+```
+
+- $`\alpha_B(T)`$ is the case-B recombination coefficient: recombinations to $n \ge 2$ only (Péquignot,
+ Petitjean & Boisson 1991 fit, multiplied by a fudge factor $F$, which is 1 by default).
+- $`\beta_B = \alpha_B\,(m_ek_BT/2\pi\hbar^2)^{3/2}\,e^{-E_{\rm ion}/4k_BT}`$ is the photo-ionisation rate
+ from $n = 2$, fixed by detailed balance.
+- $`E_{{\rm Ly}\alpha} = \tfrac34 E_{\rm ion}`$ is the Lyman-α energy.
+- $`\Lambda_{2s1s} = 8.22458\ {\rm s^{-1}}`$ is the two-photon 2s → 1s decay rate.
+- $`\Lambda_\alpha = 8\pi H/(\lambda_\alpha^3\, n_{1s})`$ is the rate at which Lyman-α photons redshift
+ out of the line, with $`\lambda_\alpha`$ the Lyman-α wavelength and
+ $`n_{1s} = (1 - x_e)\,n_{\rm H}`$ the ground-state density.
+
+$`C_r`$ is the probability that an atom in $n = 2$ reaches the ground state before being
+re-ionised. During recombination $`C_r \ll 1`$, and recombination is throttled by the two slow
+escape routes. The bracket vanishes exactly at the Saha solution, so the ODE starts from Saha at
+$`z = 1800`$. It is stiff, so it is integrated with the implicit Radau method.
+
+### Thomson optical depth
+
+Photons scatter off free electrons at the rate $`\Gamma_T = n_e\sigma_T c`$, where $`\sigma_T`$ is the
+Thomson cross-section and $c$ the speed of light. The optical depth back to redshift $z$ is
+
+```math
+\tau(z) = \int_0^z \frac{n_e(z')\,\sigma_T\, c}{(1+z')\,H(z')}\,dz', \qquad n_e = x_e\, n_{\rm H}.
+```
+
+I define the last-scattering redshift $`z_\ast`$ by $`\tau(z_\ast) = 1`$, as Planck does.
+
+### Visibility function
+
+$`e^{-\tau(z)}`$ is the probability that a photon seen today has not scattered since $z$. So the
+probability density that it last scattered at $z$ is
+
+```math
+g(z) = \frac{d\tau}{dz}\,e^{-\tau(z)}, \qquad \int_0^{z} g\,dz' = 1 - e^{-\tau(z)} \to 1 .
+```
+
+Its full width at half maximum is the thickness of the last-scattering surface.
+
+
+
+
+
+*Saha equilibrium (orange, pink) recombines far too early. The Peebles model sits almost on top of
+CAMB through last scattering. More figures: [`figures/visibility_function.png`](figures/visibility_function.png)
+(also shows $`\tau(z)`$) and [`figures/rates.png`](figures/rates.png) ($`\Gamma_T/H`$ and $`C_r`$).*
+
+## Approximations and limitations
+
+| Approximation | What it does to the result | What would fix it |
+|---|---|---|
+| Three-level atom with case-B $`\alpha_B`$, no fudge factor | $`x_e`$ too high after recombination: +14.05% at $`z = 200`$ (+0.87% with $F = 1.14$) | multi-level atom (HyRec, CosmoRec) |
+| Gas temperature equals radiation temperature | $`x_e`$ too high at late times; $`x_e(0)`$ is 2.351 × 10⁻⁴ vs 1.738 × 10⁻⁴ even with $F = 1.14$ | evolve the matter temperature with Compton heating |
+| Helium in Saha equilibrium; neutral below $`z = 1800`$ | $`x_e`$ too low in the helium era: −5.81% at $`z = 2299`$ | helium rate equations (as in RECFAST) |
+| No Lyman-series feedback, Raman scattering or two-photon corrections; none of RECFAST's extra rate corrections | ~1% wiggles in $`x_e`$ around last scattering (e.g. +1.46% at $`z = 1264`$ with $F = 1.14$) | HyRec / CosmoRec |
+| No reionisation | $\tau$ is the optical depth to recombination only; reionisation adds a further optical depth at low redshift | add a reionisation model |
+| Massless neutrinos | none in the comparison (CAMB uses the same) | massive-neutrino background |
-## Running it
+CAMB's RECFAST includes helium rate equations, matter-temperature evolution, $`F = 1.125`$ and further
+corrections. HyRec and CosmoRec are not in the pip-installed CAMB, so I could not compare with them.
+The code is a learning tool and a check of the physics. It is not a replacement for those codes in
+parameter estimation.
-You'll need Python 3.10 or newer.
+## Installation and quick start
+
+You'll need Python 3.10 or newer. With pip:
```bash
git clone https://github.com/sethuastro/CMB_Decoupling.git
cd CMB_Decoupling
python -m venv .venv
source .venv/bin/activate # on Windows: .venv\Scripts\activate
+pip install -e ".[dev]" # or, for the exact versions I used: pip install -r requirements.txt && pip install -e . --no-deps
+```
+
+With conda:
+
+```bash
+conda create -n cmb python=3.13 && conda activate cmb
pip install -e ".[dev]"
+```
-pytest # run the tests
-python scripts/make_figures.py # rebuild figures/ and results/
+Reproduce the validation figure and tables (well under a minute):
+
+```bash
+python scripts/validate.py # figures/validation.{png,pdf}, results/validation.{md,json}
+python scripts/make_figures.py # the other figures, results/key_results.{md,json}
```
Or, from Python:
@@ -108,52 +233,83 @@ Or, from Python:
```python
from cmb_decoupling import PLANCK2018, solve_peebles, last_scattering
-history = solve_peebles(PLANCK2018)
-ls = last_scattering(history.x_e, PLANCK2018)
-print(ls.z_tau_one(), ls.z_peak(), ls.fwhm()[0]) # 1089.44 1078.80 203.33
+history = solve_peebles(PLANCK2018) # Peebles x_e(z) on the Planck 2018 background
+ls = last_scattering(history.x_e, PLANCK2018) # tau(z) and g(z)
+print(f"{ls.z_tau_one():.2f} {ls.fwhm()[0]:.2f}") # 1089.44 203.33
+```
+
+**How it's organised**
+
+```
+src/cmb_decoupling/ the package
+ cosmology.py Planck 2018 parameters (defined once), H(z), n_H(z)
+ constants.py physical and atomic constants, with sources
+ saha.py Saha equilibrium for hydrogen and helium
+ peebles.py Peebles three-level atom (stiff ODE)
+ visibility.py Thomson rate, tau(z), g(z), peak and width
+ results.py summary numbers and residuals against a reference
+ reference.py the CAMB comparison run
+scripts/ validate.py (comparison with CAMB), make_figures.py (all other figures)
+tests/ pytest suite
+notebooks/ a walk-through from Saha to last scattering
+docs/theory.md the derivations, step by step, with references
+archive/ my original scripts, kept as they were
```
-## How I tested it
+## Tests and CI
+
+```bash
+pytest
+```
-The tests check things that physics says must be true:
-- $`x_e \to 1`$ at early times;
+There are 41 tests. Some check things that physics says must be true:
+- $`x_e \to 1`$ at high redshift in Saha;
- the Peebles rate vanishes exactly at the Saha solution (detailed balance);
- charge is conserved in the hydrogen–helium Saha solution;
-- $`\int g\,dz = 1 - e^{-\tau}`$.
+- $`\int g\,dz = 1 - e^{-\tau}`$ exactly, and $`\int g\,dz = 1`$ to $10^{-6}$ since $`\tau(3000) \gg 1`$.
-They also check that the answers don't change when I tighten the ODE tolerances, refine the redshift
-grid, or move the point where the ODE takes over from Saha. The comparison tests against CAMB use
-tolerances that describe the agreement I actually measured; I didn't adjust the model to pass them.
-GitHub Actions runs the tests, rebuilds the figures and runs the notebook on every push.
+Others check that the answers don't change when I tighten the ODE tolerances, refine the redshift
+grid or move the Saha-to-Peebles switch.
-## What it doesn't do (yet)
+The comparison tests against CAMB use tolerances that describe the agreement I measured, not targets
+I tuned the model to meet:
+- $`z_\ast`$, the peak and $`z_{\rm dec}`$ within 2;
+- the FWHM within 3;
+- $`\lvert\Delta x_e/x_e\rvert < 2\%`$ in the last-scattering window.
-- **A multi-level hydrogen atom.** This would fix most of the late-time excess. HyRec and CosmoRec do
- this properly.
-- **A proper helium recombination equation** instead of Saha.
-- **Evolving the matter temperature** separately from the radiation temperature.
-- **Reionisation.** The optical depth here stops at recombination; Planck measures an extra
- $\tau \approx 0.054$ from reionisation.
-- **Parameter sensitivity.** I'd like to look at how $`z_\ast`$ and the width of last scattering respond
- to $`\Omega_b h^2`$ and $`Y_p`$.
+GitHub Actions runs on every push (badge at the top):
+- lint (ruff);
+- tests on Linux (Python 3.11 and 3.13) and on Windows;
+- both scripts, rebuilt from scratch;
+- a check that the committed results tables match the code;
+- the notebook, run top to bottom.
-## References
+## Roadmap
-- Peebles, P. J. E. 1968, *ApJ* **153**, 1
-- Zel'dovich, Ya. B., Kurt, V. G. & Sunyaev, R. A. 1968, *Zh. Eksp. Teor. Fiz.* **55**, 278
-- Seager, S., Sasselov, D. D. & Scott, D. 1999, *ApJ* **523**, L1
-- Seager, S., Sasselov, D. D. & Scott, D. 2000, *ApJS* **128**, 407
-- Planck Collaboration 2020, *A&A* **641**, A6 (Planck 2018 results VI)
-- Lewis, A., Challinor, A. & Lasenby, A. 2000, *ApJ* **538**, 473 (CAMB)
-- Dodelson, S. & Schmidt, F. 2020, *Modern Cosmology*, 2nd ed.
-- Péquignot, D., Petitjean, P. & Boisson, C. 1991, *A&A* **251**, 680
+These are **planned, not done**:
-The full list is at the end of [`docs/theory.md`](docs/theory.md#9-references).
+1. **Sound horizon $`r_s(z_\ast)`$**: the comoving distance a sound wave in the photon–baryon fluid
+ travels before last scattering. This sets the angular scale of the acoustic peaks.
+2. **Angular power spectrum $`C_\ell`$**: starting from the visibility function computed here.
-## About me
+## Citation and licence
-I'm Sethulekshmi Beena, an MSc Astrophysics student at the University of Hertfordshire.
+If you use this code, please cite it with the metadata in [`CITATION.cff`](CITATION.cff); GitHub's
+"Cite this repository" button generates BibTeX from it. The code is released under the
+[MIT licence](LICENSE).
+
+**References:**
+- Peebles 1968, *ApJ* **153**, 1
+- Zel'dovich, Kurt & Sunyaev 1968, *Zh. Eksp. Teor. Fiz.* **55**, 278
+- Seager, Sasselov & Scott 1999, *ApJ* **523**, L1
+- Seager, Sasselov & Scott 2000, *ApJS* **128**, 407
+- Péquignot, Petitjean & Boisson 1991, *A&A* **251**, 680
+- Planck Collaboration 2020, *A&A* **641**, A6
+- Lewis, Challinor & Lasenby 2000, *ApJ* **538**, 473 (CAMB)
+- Dodelson & Schmidt 2020, *Modern Cosmology*, 2nd ed.
+
+The full list is in [`docs/theory.md`](docs/theory.md#9-references).
+
+**About me.** I'm Sethulekshmi Beena, MSc Astrophysics (Distinction), University of Hertfordshire.
ORCID [0009-0004-4767-2733](https://orcid.org/0009-0004-4767-2733) ·
[homepage](https://sethuastro.github.io/)
-
-If this is useful to you, there's a [`CITATION.cff`](CITATION.cff) file. The code is MIT-licensed.
diff --git a/archive/README.md b/archive/README.md
index 3161f3b..c61e63e 100644
--- a/archive/README.md
+++ b/archive/README.md
@@ -18,7 +18,10 @@ Known issues, fixed in the package:
- cosmological parameters and constants were duplicated across files without
sources.
-Using the Saha x_e, the visibility function peaks at z ≈ 1277, much earlier
+Using the Saha x_e, the visibility function peaks at z = 1283.6 (`python scripts/validate.py`), much earlier
than the z ≈ 1080 obtained once the Peebles model is used. This is the
physical failure of Saha equilibrium at late times, not a coding error; see
`docs/theory.md`.
+
+An earlier version of this note gave z ≈ 1277; that figure came from the archive scripts above
+(before the fixes listed), not from the package.
diff --git a/figures/validation.pdf b/figures/validation.pdf
new file mode 100644
index 0000000..6a677e9
Binary files /dev/null and b/figures/validation.pdf differ
diff --git a/figures/validation.png b/figures/validation.png
new file mode 100644
index 0000000..68cfee0
Binary files /dev/null and b/figures/validation.png differ
diff --git a/pyproject.toml b/pyproject.toml
index c5ded50..a777599 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
[project]
name = "cmb-decoupling"
-version = "0.1.0"
+version = "1.0.0"
description = "Hydrogen recombination (Saha, Peebles) and CMB photon decoupling, validated against CAMB/RECFAST"
readme = "README.md"
requires-python = ">=3.10"
diff --git a/results/validation.json b/results/validation.json
new file mode 100644
index 0000000..89ab6aa
--- /dev/null
+++ b/results/validation.json
@@ -0,0 +1,97 @@
+{
+ "camb_version": "2.0.5",
+ "camb_recombination": "Recfast",
+ "last_scattering": {
+ "camb": {
+ "z_star": 1089.8900710994267,
+ "z_peak": 1078.3679452200388,
+ "fwhm": 204.00059486771897
+ },
+ "peebles": {
+ "z_star": 1089.4370745392039,
+ "z_peak": 1078.7955467016986,
+ "fwhm": 203.33067774818244
+ },
+ "peebles_fudge": {
+ "z_star": 1090.2769961061613,
+ "z_peak": 1079.0473446893725,
+ "fwhm": 202.72657409606586
+ }
+ },
+ "camb_recfast_settings": {
+ "fudge": 1.125,
+ "Hswitch": true,
+ "gauss_centres_z": [
+ 1451.884439702605,
+ 867.4446628700608
+ ]
+ },
+ "residuals_vs_camb": {
+ "peebles": {
+ "window_lo": 966.507145886229,
+ "window_hi": 1170.507740753948,
+ "max_abs_window": 0.010813545422290627,
+ "z_max_abs_window": 1170.5,
+ "min_residual": -0.03579640434004827,
+ "z_min_residual": 2000.0,
+ "max_residual": 0.1404936065241551,
+ "z_max_residual": 200.0,
+ "z_min": 200.0,
+ "z_max": 2000.0,
+ "residual_at_freeze_out": 0.1404936065607214
+ },
+ "peebles_fudge": {
+ "window_lo": 966.507145886229,
+ "window_hi": 1170.507740753948,
+ "max_abs_window": 0.010390256944198484,
+ "z_max_abs_window": 1170.5,
+ "min_residual": -0.03579640434004827,
+ "z_min_residual": 2000.0,
+ "max_residual": 0.014606252660196661,
+ "z_max_residual": 1264.0,
+ "z_min": 200.0,
+ "z_max": 2000.0,
+ "residual_at_freeze_out": 0.008650244424860842
+ }
+ },
+ "residuals_hydrogen_era": {
+ "peebles": {
+ "window_lo": 966.507145886229,
+ "window_hi": 1170.507740753948,
+ "max_abs_window": 0.01081268868056795,
+ "z_max_abs_window": 1170.488888888889,
+ "min_residual": -0.003268465909092355,
+ "z_min_residual": 1800.0,
+ "max_residual": 0.14049360653415421,
+ "z_max_residual": 200.0,
+ "z_min": 200.0,
+ "z_max": 1800.0
+ },
+ "peebles_fudge": {
+ "window_lo": 966.507145886229,
+ "window_hi": 1170.507740753948,
+ "max_abs_window": 0.010389339233394024,
+ "z_max_abs_window": 1170.488888888889,
+ "min_residual": -0.014599229913939693,
+ "z_min_residual": 861.5111111111112,
+ "max_residual": 0.014606253069484154,
+ "z_max_residual": 1264.0,
+ "z_min": 200.0,
+ "z_max": 1800.0
+ }
+ },
+ "residuals_helium_era": {
+ "z_min": 1800.0,
+ "z_max": 3000.0,
+ "min_residual": -0.05810588572572306,
+ "z_min_residual": 2299.4
+ },
+ "tau_integrator_check": {
+ "z_star_from_camb_xe": 1089.8900710994267,
+ "z_star_camb_internal": 1089.8922558485872
+ },
+ "saha": {
+ "z_star": 1291.1323350280427,
+ "z_peak": 1283.606148602622
+ }
+}
diff --git a/results/validation.md b/results/validation.md
new file mode 100644
index 0000000..6bdb5fd
--- /dev/null
+++ b/results/validation.md
@@ -0,0 +1,23 @@
+| | Peebles (this work) | Peebles, F = 1.14 (this work) | CAMB 2.0.5 (RECFAST) |
+|---|---|---|---|
+| $z_\ast$ where $\tau = 1$ | 1089.44 | 1090.28 | 1089.89 |
+| Peak of $g(z)$ per unit $z$ | 1078.80 | 1079.05 | 1078.37 |
+| FWHM of $g(z)$ | 203.33 | 202.73 | 204.00 |
+
+Residuals $\Delta x_e/x_e = x_e/x_e^{\rm CAMB} - 1$ on $z = 200$–$2000$ (grid spacing $\Delta z = 0.1$). Last-scattering window = FWHM of CAMB's $g(z)$, $z = 966.5$–$1170.5$.
+
+| | Peebles (this work) | Peebles, F = 1.14 (this work) |
+|---|---|---|
+| (a) max $\lvert\Delta x_e/x_e\rvert$ in the last-scattering window | 1.08% at $z = 1170.5$ | 1.04% at $z = 1170.5$ |
+| (b) most negative residual | -3.58% at $z = 2000.0$ | -3.58% at $z = 2000.0$ |
+| (c) most positive residual | +14.05% at $z = 200.0$ | +1.46% at $z = 1264.0$ |
+| residual at $z = 200$ | +14.05% | +0.87% |
+| range on $z = 200$–$1800$ (hydrogen only) | -0.33% to +14.05% | -1.46% to +1.46% |
+
+Helium era, $z = 1800$–$3000$: most negative residual -5.81% at $z = 2299.4$ (both models; helium is in Saha).
+
+$\tau$-integrator check: CAMB's $x_e$ through this code gives $z_\ast = 1089.8901$; CAMB's internal value is 1089.8923 (difference -0.0022).
+
+CAMB RECFAST settings: $F = 1.125$, Hswitch = True (Gaussian corrections centred at $z \approx 1452$ and 867).
+
+Saha (H + He) for comparison: $z_\ast = 1291.13$, peak of $g(z)$ at 1283.61.
diff --git a/scripts/validate.py b/scripts/validate.py
new file mode 100644
index 0000000..e43a922
--- /dev/null
+++ b/scripts/validate.py
@@ -0,0 +1,265 @@
+"""Validate the Peebles ionisation history against CAMB (RECFAST) on a matched background.
+
+Usage
+-----
+ python scripts/validate.py
+
+Both codes use the same Planck 2018 background (massless neutrinos, same Y_p
+and T_cmb, no reionisation), so differences come from the recombination
+physics alone. Nothing here is fitted or tuned to the reference.
+
+Outputs
+-------
+figures/validation.{pdf,png}
+ (a) x_e(z), (b) fractional residual Delta x_e / x_e against CAMB over
+ z = 200-2000 for the pure Peebles model and the F = 1.14 literature variant,
+ (c) the visibility function g(z).
+results/validation.{json,md}
+ Last-scattering quantities and the residual summary, also printed to the console.
+"""
+
+from __future__ import annotations
+
+import json
+import sys
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from cmb_decoupling import PLANCK2018, last_scattering, solve_peebles, x_e_saha_full
+from cmb_decoupling.peebles import DEFAULT_Z_SWITCH, RECFAST_FUDGE
+from cmb_decoupling.plotting import GUIDE_COLOR, STYLES, save, use_style
+from cmb_decoupling.reference import camb_reference
+from cmb_decoupling.results import (
+ Z_FREEZE_OUT,
+ Z_VALIDATION_MAX,
+ Z_VALIDATION_MIN,
+ fractional_residual,
+ residual_summary,
+)
+
+ROOT = Path(__file__).resolve().parents[1]
+FIG_DIR = ROOT / "figures"
+RES_DIR = ROOT / "results"
+
+N_PLOT = 2000
+Z_TICKS = (200, 300, 500, 700, 1000, 1500, 2000)
+#: Redshift range shown in the visibility-function panel.
+G_PLOT_RANGE = (700.0, 1500.0)
+#: Upper redshift for the helium-era residual (He I and He II recombination in Saha).
+Z_HELIUM_MAX = 3000.0
+#: Grid points on [z_switch, Z_HELIUM_MAX] (spacing 0.1 in z, as for the main range).
+N_HELIUM = 12001
+
+MODELS = ("peebles", "peebles_fudge")
+LABELS = {"peebles": "Peebles (this work)", "peebles_fudge": f"Peebles, F = {RECFAST_FUDGE} (this work)"}
+
+
+def compute() -> dict[str, object]:
+ """Run both codes and collect every number reported by this script."""
+ cosmo = PLANCK2018
+ ref = camb_reference(cosmo)
+ x_e = {
+ "camb": ref.x_e,
+ "peebles": solve_peebles(cosmo).x_e,
+ "peebles_fudge": solve_peebles(cosmo, fudge=RECFAST_FUDGE).x_e,
+ }
+ ls = {k: last_scattering(f, cosmo) for k, f in x_e.items()}
+ last_sc = {k: {"z_star": v.z_tau_one(), "z_peak": v.z_peak(), "fwhm": v.fwhm()[0]} for k, v in ls.items()}
+ resid = {k: residual_summary(x_e[k], ref.x_e, cosmo).as_dict() for k in MODELS}
+ for k in MODELS:
+ resid[k]["residual_at_freeze_out"] = float(
+ fractional_residual(x_e[k], ref.x_e, np.array(Z_FREEZE_OUT))
+ )
+ # Hydrogen-only era (Peebles ODE region) and the helium era beyond the main range.
+ z_switch = DEFAULT_Z_SWITCH
+ hydrogen_era = {
+ k: residual_summary(x_e[k], ref.x_e, cosmo, z_min=Z_VALIDATION_MIN, z_max=z_switch).as_dict()
+ for k in MODELS
+ }
+ z_he = np.linspace(z_switch, Z_HELIUM_MAX, N_HELIUM)
+ r_he = fractional_residual(x_e["peebles"], ref.x_e, z_he)
+ i_he = int(np.argmin(r_he))
+ saha = last_scattering(lambda z: x_e_saha_full(z, cosmo), cosmo)
+ rec = ref.params.Recomb
+ recfast = {
+ "fudge": rec.RECFAST_fudge,
+ "Hswitch": bool(rec.RECFAST_Hswitch),
+ "gauss_centres_z": [float(np.exp(rec.zGauss1) - 1.0), float(np.exp(rec.zGauss2) - 1.0)],
+ }
+ return {
+ "x_e": x_e,
+ "ls": ls,
+ "numbers": {
+ "camb_version": __import__("camb").__version__,
+ "camb_recombination": ref.recombination_model,
+ "last_scattering": last_sc,
+ "camb_recfast_settings": recfast,
+ "residuals_vs_camb": resid,
+ "residuals_hydrogen_era": hydrogen_era,
+ "residuals_helium_era": {
+ "z_min": z_switch,
+ "z_max": Z_HELIUM_MAX,
+ "min_residual": float(r_he[i_he]),
+ "z_min_residual": float(z_he[i_he]),
+ },
+ "tau_integrator_check": {
+ "z_star_from_camb_xe": ls["camb"].z_tau_one(),
+ "z_star_camb_internal": ref.z_star,
+ },
+ "saha": {"z_star": saha.z_tau_one(), "z_peak": saha.z_peak()},
+ },
+ }
+
+
+def make_figure(data: dict[str, object]) -> None:
+ """Three-panel validation figure: x_e, fractional residual, visibility function."""
+ x_e, ls, num = data["x_e"], data["ls"], data["numbers"]
+ window = (
+ num["residuals_vs_camb"]["peebles"]["window_lo"],
+ num["residuals_vs_camb"]["peebles"]["window_hi"],
+ )
+ z = np.geomspace(Z_VALIDATION_MIN, Z_VALIDATION_MAX, N_PLOT)
+
+ fig = plt.figure(figsize=(6.4, 8.4))
+ # Row 2 is an empty spacer: panel (c) has its own redshift axis.
+ gs = fig.add_gridspec(4, 1, height_ratios=[2.0, 1.4, 0.55, 2.0], hspace=0.08)
+ ax_x = fig.add_subplot(gs[0])
+ ax_r = fig.add_subplot(gs[1], sharex=ax_x)
+ ax_g = fig.add_subplot(gs[3])
+
+ for k in ("camb", *MODELS):
+ ax_x.plot(z, x_e[k](z), **STYLES[k])
+ ax_x.set_yscale("log")
+ ax_x.set_ylabel(r"Free-electron fraction $x_e = n_e/n_{\rm H}$")
+ ax_x.legend(loc="lower left")
+ ax_x.tick_params(labelbottom=False)
+
+ for ax in (ax_x, ax_r):
+ ax.axvspan(*window, color=GUIDE_COLOR, alpha=0.15, lw=0)
+ ax_r.axhline(0.0, color=GUIDE_COLOR, lw=0.8)
+ for k in MODELS:
+ ax_r.plot(z, 100.0 * fractional_residual(x_e[k], x_e["camb"], z), **STYLES[k])
+ ax_r.set_ylabel(r"$\Delta x_e / x_e^{\rm CAMB}$ [%]")
+ ax_r.text(
+ np.sqrt(window[0] * window[1]), 0.95, "last-scattering\nwindow", transform=ax_r.get_xaxis_transform(),
+ ha="center", va="top", fontsize=8, color="#404040",
+ ) # fmt: skip
+ ax_r.set_xscale("log")
+ ax_r.set_xlim(Z_VALIDATION_MAX, Z_VALIDATION_MIN)
+ ax_r.set_xticks(Z_TICKS, [str(t) for t in Z_TICKS])
+ ax_r.xaxis.set_minor_locator(plt.NullLocator())
+ ax_r.set_xlabel(r"Redshift $z$")
+
+ ax_g.axvspan(*window, color=GUIDE_COLOR, alpha=0.15, lw=0, label="FWHM of CAMB $g(z)$")
+ for k in ("camb", *MODELS):
+ ax_g.plot(ls[k].z, ls[k].g * 1e3, **STYLES[k])
+ ax_g.set_xlim(*G_PLOT_RANGE[::-1])
+ ax_g.set_ylim(bottom=0)
+ ax_g.set_xlabel(r"Redshift $z$")
+ ax_g.set_ylabel(r"$g(z) = \frac{d\tau}{dz}\,e^{-\tau}$ [$10^{-3}$ per unit $z$]")
+ ax_g.legend(loc="upper left")
+
+ # Panel tags placed where no curve passes.
+ for ax, tag, y, va in ((ax_x, "a", 0.95, "top"), (ax_r, "b", 0.06, "bottom"), (ax_g, "c", 0.95, "top")):
+ ax.text(0.98, y, f"({tag})", transform=ax.transAxes, ha="right", va=va, fontsize=11)
+ save(fig, FIG_DIR, "validation")
+
+
+def write_tables(num: dict[str, object]) -> str:
+ """Write results/validation.{json,md} and return the Markdown."""
+ ls, res = num["last_scattering"], num["residuals_vs_camb"]
+ cols = ("peebles", "peebles_fudge", "camb")
+ head = f"| | {LABELS['peebles']} | {LABELS['peebles_fudge']} | CAMB {num['camb_version']} (RECFAST) |"
+ lines = [head, "|---|---|---|---|"]
+ for key, label in (
+ ("z_star", r"$z_\ast$ where $\tau = 1$"),
+ ("z_peak", "Peak of $g(z)$ per unit $z$"),
+ ("fwhm", "FWHM of $g(z)$"),
+ ):
+ lines.append(f"| {label} | " + " | ".join(f"{ls[c][key]:.2f}" for c in cols) + " |")
+
+ def pct(x: float) -> str:
+ return f"{100 * x:+.2f}%"
+
+ p = res["peebles"]
+ lines += [
+ "",
+ r"Residuals $\Delta x_e/x_e = x_e/x_e^{\rm CAMB} - 1$ "
+ f"on $z = {p['z_min']:.0f}$–${p['z_max']:.0f}$ (grid spacing $\\Delta z = 0.1$). "
+ "Last-scattering window = FWHM of CAMB's $g(z)$, "
+ f"$z = {p['window_lo']:.1f}$–${p['window_hi']:.1f}$.",
+ "",
+ f"| | {LABELS['peebles']} | {LABELS['peebles_fudge']} |",
+ "|---|---|---|",
+ ]
+ rows = (
+ (
+ r"(a) max $\lvert\Delta x_e/x_e\rvert$ in the last-scattering window",
+ "max_abs_window",
+ "z_max_abs_window",
+ True,
+ ),
+ ("(b) most negative residual", "min_residual", "z_min_residual", False),
+ ("(c) most positive residual", "max_residual", "z_max_residual", False),
+ )
+ for label, key, zkey, absolute in rows:
+ cells = [
+ (f"{100 * res[m][key]:.2f}%" if absolute else pct(res[m][key])) + f" at $z = {res[m][zkey]:.1f}$"
+ for m in MODELS
+ ]
+ lines.append(f"| {label} | " + " | ".join(cells) + " |")
+ lines.append(
+ f"| residual at $z = {Z_FREEZE_OUT:.0f}$ | "
+ + " | ".join(pct(res[m]["residual_at_freeze_out"]) for m in MODELS)
+ + " |"
+ )
+ h_era = num["residuals_hydrogen_era"]
+ he = num["residuals_helium_era"]
+ lines.append(
+ f"| range on $z = {Z_VALIDATION_MIN:.0f}$–${DEFAULT_Z_SWITCH:.0f}$ (hydrogen only) | "
+ + " | ".join(f"{pct(h_era[m]['min_residual'])} to {pct(h_era[m]['max_residual'])}" for m in MODELS)
+ + " |"
+ )
+ tc = num["tau_integrator_check"]
+ rf = num["camb_recfast_settings"]
+ lines += [
+ "",
+ f"Helium era, $z = {he['z_min']:.0f}$–${he['z_max']:.0f}$: most negative residual "
+ f"{pct(he['min_residual'])} at $z = {he['z_min_residual']:.1f}$ (both models; helium is in Saha).",
+ "",
+ "$\\tau$-integrator check: CAMB's $x_e$ through this code gives "
+ f"$z_\\ast = {tc['z_star_from_camb_xe']:.4f}$; "
+ f"CAMB's internal value is {tc['z_star_camb_internal']:.4f} "
+ f"(difference {tc['z_star_from_camb_xe'] - tc['z_star_camb_internal']:+.4f}).",
+ "",
+ f"CAMB RECFAST settings: $F = {rf['fudge']}$, Hswitch = {rf['Hswitch']} "
+ f"(Gaussian corrections centred at $z \\approx {rf['gauss_centres_z'][0]:.0f}$ "
+ f"and {rf['gauss_centres_z'][1]:.0f}).",
+ ]
+ saha = num["saha"]
+ lines += [
+ "",
+ f"Saha (H + He) for comparison: $z_\\ast = {saha['z_star']:.2f}$, "
+ f"peak of $g(z)$ at {saha['z_peak']:.2f}.",
+ ]
+ md = "\n".join(lines) + "\n"
+ RES_DIR.mkdir(exist_ok=True)
+ (RES_DIR / "validation.json").write_text(json.dumps(num, indent=2) + "\n", encoding="utf-8")
+ (RES_DIR / "validation.md").write_text(md, encoding="utf-8")
+ return md
+
+
+def main() -> None:
+ """Compute, plot and tabulate the validation against CAMB."""
+ if hasattr(sys.stdout, "reconfigure"):
+ sys.stdout.reconfigure(encoding="utf-8")
+ use_style()
+ data = compute()
+ make_figure(data)
+ print(write_tables(data["numbers"]))
+
+
+if __name__ == "__main__":
+ main()
diff --git a/src/cmb_decoupling/__init__.py b/src/cmb_decoupling/__init__.py
index 015a99b..b9c65f1 100644
--- a/src/cmb_decoupling/__init__.py
+++ b/src/cmb_decoupling/__init__.py
@@ -19,4 +19,4 @@
"x_e_saha_hydrogen",
]
-__version__ = "0.1.0"
+__version__ = "1.0.0"
diff --git a/src/cmb_decoupling/results.py b/src/cmb_decoupling/results.py
index e2fa736..58ea5d8 100644
--- a/src/cmb_decoupling/results.py
+++ b/src/cmb_decoupling/results.py
@@ -4,6 +4,9 @@
from dataclasses import asdict, dataclass
+import numpy as np
+from numpy.typing import NDArray
+
from .cosmology import Cosmology
from .visibility import XeFunction, decoupling_redshift, last_scattering
@@ -12,6 +15,12 @@
#: approximation of this package becomes poor (see docs/theory.md).
Z_FREEZE_OUT: float = 200.0
+#: Redshift range over which x_e is compared with the reference (scripts/validate.py).
+Z_VALIDATION_MIN: float = 200.0
+Z_VALIDATION_MAX: float = 2000.0
+#: Number of points of the uniform comparison grid (spacing 0.1 in z).
+N_VALIDATION: int = 18001
+
@dataclass(frozen=True)
class KeyResults:
@@ -69,3 +78,98 @@ def key_results(x_e: XeFunction, cosmo: Cosmology) -> KeyResults:
x_e_freeze_out=float(x_e(Z_FREEZE_OUT)),
x_e_today=float(x_e(0.0)),
)
+
+
+def fractional_residual(x_e: XeFunction, x_e_ref: XeFunction, z: NDArray[np.float64]) -> NDArray[np.float64]:
+ """Fractional residual Delta x_e / x_e = x_e(z) / x_e_ref(z) - 1 [dimensionless].
+
+ Parameters
+ ----------
+ x_e, x_e_ref : callable
+ Ionisation history under test and reference ionisation history, x_e(z) = n_e / n_H.
+ z : ndarray
+ Redshifts at which to compare.
+ """
+ return x_e(z) / x_e_ref(z) - 1.0
+
+
+@dataclass(frozen=True)
+class ResidualSummary:
+ """Where and by how much an ionisation history departs from a reference.
+
+ All residuals are fractional, Delta x_e / x_e = x_e / x_e_ref - 1 [dimensionless].
+
+ Attributes
+ ----------
+ window_lo, window_hi : float
+ Last-scattering window: the half-maximum redshifts of the *reference* g(z).
+ max_abs_window : float
+ Largest |Delta x_e / x_e| inside the last-scattering window.
+ z_max_abs_window : float
+ Redshift where ``max_abs_window`` occurs.
+ min_residual, z_min_residual : float
+ Most negative residual over the full comparison range, and its redshift.
+ max_residual, z_max_residual : float
+ Most positive residual over the full comparison range, and its redshift.
+ z_min, z_max : float
+ Comparison range.
+ """
+
+ window_lo: float
+ window_hi: float
+ max_abs_window: float
+ z_max_abs_window: float
+ min_residual: float
+ z_min_residual: float
+ max_residual: float
+ z_max_residual: float
+ z_min: float
+ z_max: float
+
+ def as_dict(self) -> dict[str, float]:
+ """Return the summary as a plain dictionary."""
+ return asdict(self)
+
+
+def residual_summary(
+ x_e: XeFunction,
+ x_e_ref: XeFunction,
+ cosmo: Cosmology,
+ z_min: float = Z_VALIDATION_MIN,
+ z_max: float = Z_VALIDATION_MAX,
+ n: int = N_VALIDATION,
+) -> ResidualSummary:
+ """Summarise the fractional x_e residual against a reference.
+
+ Extrema are taken on a uniform grid of ``n`` points over [``z_min``, ``z_max``]
+ (no interpolation between grid points).
+
+ Parameters
+ ----------
+ x_e, x_e_ref : callable
+ Ionisation history under test and reference, x_e(z) = n_e / n_H.
+ cosmo : Cosmology
+ Background cosmology (used to compute the reference visibility function).
+ z_min, z_max : float, optional
+ Comparison range in redshift.
+ n : int, optional
+ Number of grid points.
+ """
+ z = np.linspace(z_min, z_max, n)
+ r = fractional_residual(x_e, x_e_ref, z)
+ _, w_lo, w_hi = last_scattering(x_e_ref, cosmo).fwhm()
+ in_window = (z >= w_lo) & (z <= w_hi)
+ i_w = int(np.argmax(np.abs(r[in_window])))
+ i_min, i_max = int(np.argmin(r)), int(np.argmax(r))
+ return ResidualSummary(
+ window_lo=float(w_lo),
+ window_hi=float(w_hi),
+ max_abs_window=float(np.abs(r[in_window][i_w])),
+ z_max_abs_window=float(z[in_window][i_w]),
+ min_residual=float(r[i_min]),
+ z_min_residual=float(z[i_min]),
+ max_residual=float(r[i_max]),
+ z_max_residual=float(z[i_max]),
+ z_min=z_min,
+ z_max=z_max,
+ )
diff --git a/tests/test_regression.py b/tests/test_regression.py
index 2f1e53f..1a94c3d 100644
--- a/tests/test_regression.py
+++ b/tests/test_regression.py
@@ -9,7 +9,7 @@
import pytest
-from cmb_decoupling.results import key_results
+from cmb_decoupling.results import key_results, residual_summary
@pytest.fixture(scope="module")
@@ -34,6 +34,24 @@ def test_freeze_out_against_camb(ours, ref):
assert ours.x_e_freeze_out == pytest.approx(ref.x_e_freeze_out, rel=0.25)
+@pytest.fixture(scope="module")
+def residuals(cosmo, peebles, camb_ref):
+ return residual_summary(peebles.x_e, camb_ref.x_e, cosmo)
+
+
+def test_xe_residual_in_last_scattering_window(residuals):
+ # Measured: 1.08% (z = 966.5-1170.5, the FWHM of CAMB's g(z)); see scripts/validate.py.
+ assert residuals.max_abs_window < 0.02
+
+
+def test_xe_residual_full_range(residuals):
+ # Measured on z = 200-2000: -3.58% (helium in Saha, z = 2000) and +14.05%
+ # (three-level atom at late times, z = 200). Bounds describe this, not a target.
+ assert -0.05 < residuals.min_residual < 0.0
+ assert 0.0 < residuals.max_residual < 0.20
+ assert residuals.z_max_residual == pytest.approx(200.0)
+
+
def test_self_regression(ours):
# Values recorded on 2026-10-22 (results/key_results.json).
assert ours.z_tau_one == pytest.approx(1089.44, abs=0.05)
diff --git a/tests/test_visibility.py b/tests/test_visibility.py
index 957f0a6..95f39df 100644
--- a/tests/test_visibility.py
+++ b/tests/test_visibility.py
@@ -29,6 +29,13 @@ def test_visibility_normalisation(cosmo, peebles):
assert trapezoid(ls.g, ls.z) == pytest.approx(1 - np.exp(-ls.tau[-1]), rel=1e-6)
+def test_visibility_integrates_to_one(cosmo, peebles):
+ # tau(z_max = 3000) >> 1, so e^{-tau} is negligible and g(z) is a normalised density.
+ ls = last_scattering(peebles.x_e, cosmo)
+ assert ls.tau[-1] > 100.0
+ assert trapezoid(ls.g, ls.z) == pytest.approx(1.0, abs=1e-6)
+
+
def test_grid_convergence(cosmo, peebles):
coarse = last_scattering(peebles.x_e, cosmo, n=15001)
fine = last_scattering(peebles.x_e, cosmo, n=60001)