Comparing Saha, a three-level atom and a frozen CLASS reference

The reference is the CLASS repository at tag v3.3.4, commit e85808324f51fc694d12e3ed7439552a3c3f9540. CLASS v3.3.4 on GitHub We read four files from it. We did not compile or run CLASS, HyRec or CosmoTherm for this page. Every "reference" number below is a value printed in a file that ships with the release.

Files, versions and units

All scripts use only the Python 3 standard library and were run with Python 3.9.6. The materials are in this collection's materials folder:

• ac75_model.py, the reference model that the Lab's TypeScript code mirrors;
• inputs.json, the background, solver settings and five energy releases;
• run_recombination.py and its saved recombination-output.json;
• check_distortion_kernels.py and kernel-checks.json;
• compare_class_files.py and class-comparison.json;
• hyrec-reference-extract.json, a 60-row sample of the HyRec-2 file, and the README.

Redshift is dimensionless and the scripts integrate in $s=\ln(1+z)$. Densities are in $\mathrm{m^{-3}}$, rates in $\mathrm{s^{-1}}$, energies in joules internally and electronvolts in the text. Spectra are in $\mathrm{kJy\,sr^{-1}}$, where $1\,\mathrm{kJy\,sr^{-1}}=10^{-23}\,\mathrm{W\,m^{-2}\,Hz^{-1}\,sr^{-1}}$. The background takes $h=0.6736$, $T_0=2.7255\,\mathrm K$, $\Omega_{\mathrm b}=0.04941$, $\Omega_{\mathrm{cb}}=0.3124$, helium mass fraction $0.245$ and $N_{\mathrm{eff}}=3.046$. These match external/HyRec2020/input.dat in CLASS, except that the one massive neutrino is treated as massless radiation.

mkdir ac75-check && cd ac75-check
base=https://sandboxphysics.com/stories/cosmic-transparency/materials
for name in ac75_model.py inputs.json run_recombination.py check_distortion_kernels.py compare_class_files.py hyrec-reference-extract.json; do curl -fsSLO "$base/$name"; done
python3 -I run_recombination.py inputs.json hyrec-reference-extract.json my-recombination-output.json
python3 -I check_distortion_kernels.py my-kernel-checks.json
git clone --depth 1 --branch v3.3.4 https://github.com/lesgourg/class_public.git
python3 -I compare_class_files.py class_public my-class-comparison.json

The comparison script refuses to run if any of the four CLASS files has a different SHA256 from the one written in the script. Compare your three new JSON files with the saved ones; on another machine the last digits may differ.

Step 1: three histories at the same redshifts

The script solves the three-level equation from Saha at redshift 1800 down to 50, with a backward-Euler step of $5\times10^{-4}$ in $s$. Backward Euler evaluates the rate at the end of each step and solves for it with Newton iterations, which keeps the stiff early phase stable.

Where to look: the last column. The largest size for redshifts at or below 1500 is $0.039$, at redshift 100. Common misreading: the factor 1.14 barely matters at redshift 1100 but changes the free fraction at redshift 200 by 13 per cent, because it enters the recombination coefficient that controls the freeze-out.

The visibility summaries follow. The three-level history peaks at 1081.0 and reaches optical depth one at 1092.0. Integrating the HyRec-2 file below redshift 1600, where helium no longer matters, gives 1080 and 1089.8. Planck 2018 lists $z_\ast=1089.92\pm0.25$. Planck Collaboration, 2018, Table 2 That paper's definition of $z_\ast$ was not checked here.

Step 2: numerical checks and failure cases

Halving the step moves the optical-depth-one redshift by $0.078$ and changes $x$ at redshift 1100 by $-5.2\times10^{-4}$ relative. Both are far below the gap to HyRec-2, so the gap is physics, not step size.

Two failures are saved on purpose. An explicit Euler step of $5\times10^{-3}$ overshoots the Saha balance at the start and produces nan by redshift 1500. Starting the three-level solution from Saha at redshift 1300 instead of 1800 gives $x=0.105$ at redshift 1100, against $0.142$. By 1300 the true history has already left Saha, so the starting value is already wrong, and that error carries through to redshift 1100.

Step 3: the full HyRec-2 file

The comparison script reads all 8000 rows of external/HyRec2020/output_xe.dat. Over redshifts 50 to 1600 the largest relative difference is $+8.8$ per cent, at redshift 50. The model lets gas share the photon temperature down to redshift 50. Real gas cools faster late in the dark ages, which changes the recombination rate, so this assumption is the most likely cause; we have not isolated it. At redshift 8000 the file gives $x_{\mathrm e}=1.16341$, against $1+2f_{\mathrm{He}}=1.16343$ for fully ionized helium. This confirms that the file counts helium electrons, so comparisons above redshift 1600 must allow for them.

Step 4: the distortion shapes without reusing their formulas

check_distortion_kernels.py tests the closed forms $G$, $Y$ and $M$ against quantities they do not contain:

1. Simpson integrals of energy and number against exact Planck integrals: $G$ carries 4 times the blackbody energy and 3 times its number; $Y$ carries 4 and $4\times10^{-9}$; $M$ carries $0.71395$ and $4\times10^{-5}$. The inverse of $0.71395$ is $1.4007$, the coefficient used for $\mu$.
2. A finite-difference Kompaneets operator applied to the Planck occupation number reproduces $Y$ to $1.7\times10^{-7}$ of its peak.
3. A Bose–Einstein spectrum with $\mu=10^{-6}$, its temperature refitted to keep the photon number, reproduces $M$ to $6.6\times10^{-5}$ of its peak.

Zero crossings come out at $217.51\,\mathrm{GHz}$ for $Y$ and $124.5\,\mathrm{GHz}$ for $M$.

Step 5: energy in the CosmoTherm Green's functions

Greens_data.dat in CLASS stores the spectral response to a unit energy release at 120 redshifts. The CLASS README says it is adapted from the public CosmoTherm 1.0.3 by Jens Chluba; the method is described by Chluba, 2013, section 2. Integrating each tabulated response over the file's frequency grid and dividing by the blackbody energy gives the share of the release that the tabulated column carries.

The energy closes to 0.1 per cent up to $2.6\times10^{5}$. Above that it drifts by up to 16 per cent. The README says the file also stores a final blackbody temperature shift and the energy it carries, in separate rows, and this integral leaves those rows out. Most early energy ends up in that shift, so it is the likely reason, but we have not added it back to confirm. Use the column alone for energy checks only below $2.6\times10^{5}$.

Step 6: two different ways to split a distortion

FIRAS_branching_ratios.dat lists, for each release redshift, the $\mu$, $y$ and temperature shift per unit $\Delta\rho/\rho_\gamma$. The comparison script repeats the steps of external/distortions/generate_PCA_files.py on the 43 FIRAS frequencies from 68 to 639 GHz. It builds orthonormal vectors from $Y$, then $M$, then $G$ by Gram–Schmidt and projects each Green's function onto them, the procedure written as equations 3.10 to 3.15 by Lucca and colleagues. Lucca and colleagues, 2020, section 3.2.2 The repeat matches the table: at $z=5\times10^{4}$ the table gives $\mu$ per unit energy of $1.501$, and our projection at the nearest tabulated redshift, $4.96\times10^{4}$, gives $1.488$.

The Lab uses the soft branching fits instead, and these differ strongly in the transition era. At $z=5\times10^{4}$ they give $0.743$ against the table's $1.501$; at $z=2\times10^{4}$, $0.177$ against $0.341$. Above $5\times10^{5}$ their $\mu$ values agree within 3 per cent, and at $2\times10^{5}$ within 9 per cent. Below $5\times10^{3}$ their $y$ values agree within 0.3 per cent.

Common misreading: treating one of these as the correct $\mu$. A projection depends on the frequencies and the order of the basis. Lucca and colleagues note in footnote 13 that the ordering changes the relative sizes of the $y$, $\mu$ and temperature parts, and that the separation between $y$ and $\mu$ is artificial. What both methods must agree on is the total spectrum and its energy. That is why the Lab shows $\mu$ and $y$ as the amplitudes of a stated approximation and draws the FIRAS limits as reference shapes.

What these comparisons do and do not establish

They establish that the hydrogen recombination history and visibility used by the Lab follow a detailed published calculation to a few per cent, that the step size is converged, and that the distortion shapes and their energy factors are right. They do not test helium, the gas temperature after redshift 200, Lyman-alpha radiative transfer, or how the injected energy cascades down to heat. CLASS was not executed, so nothing here checks CLASS's own handling of a given energy release.

Reusable checklist for the next code comparison:

1. Pin the reference version and record file checksums.
2. Match every background parameter, and list the ones you could not match.
3. Compare at fixed redshifts and in derived summaries such as the visibility peak.
4. Halve the step and confirm the gap to the reference is larger than the change.
5. Save at least one failure case, so a reader can see what a wrong setup looks like.
6. Check every shape against a quantity it does not contain: a moment, a limit or an independent operator.
7. Before comparing split parameters, confirm both sides use the same definition of the split.

The HyRec-2 extract and the three CLASS files are read under attribution to their authors. No licence file was found in the CLASS v3.3.4 tree, so Sandbox Physics redistributes only the 60-row extract and asks readers to download the original files from the release.