# Recombination and spectral-distortion materials (AC75)

These files reproduce every number in the Sandbox Physics Lab
`/labs/recombination-spectral-distortions/` and the tutorial
`/stories/cosmic-transparency/compare-recombination-codes/`.
They use only the Python 3 standard library. The saved outputs were made with
Python 3.9.6.

## Run

    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

`compare_class_files.py` reads four files from CLASS v3.3.4
(commit e85808324f51fc694d12e3ed7439552a3c3f9540) and stops if a SHA256 differs:

| File in the CLASS tree | SHA256 |
| --- | --- |
| external/HyRec2020/output_xe.dat | 1ac2ddad6aa0e3c4080883c71479e2e398669468e4813075ec61bca86b7afc40 |
| external/distortions/Greens_data.dat | 03972c42a530315a768c973a45c343ee744a8d64da93493c9c68057f92327bed |
| external/distortions/FIRAS_branching_ratios.dat | 966daa8fa99d7e15de6961f33e0fac8800ae5d64464c21d3526c52cd8fa42864 |
| external/distortions/FIRAS_nu_delta_I.dat | 92f583021393c543b512c038a9a829b766027ab18938a2b527a1d8fdac889aa7 |

Raw copies at that commit:
https://raw.githubusercontent.com/lesgourg/class_public/e85808324f51fc694d12e3ed7439552a3c3f9540/<path above>

`make_hyrec_extract.py` rebuilds `hyrec-reference-extract.json` (60 rows)
from output_xe.dat.

CLASS, HyRec and CosmoTherm were NOT executed. Every reference value is read
from a file shipped with the CLASS v3.3.4 release.

## Files

| File | Role |
| --- | --- |
| ac75_model.py | Reference model mirrored by the Lab's TypeScript code |
| inputs.json | Background, solver settings and the five Lab releases |
| run_recombination.py → recombination-output.json | Saha, three-level atom, HyRec-2 extract, visibility, step check, failure cases, releases |
| check_distortion_kernels.py → kernel-checks.json | Independent checks of the G, Y and M shapes |
| compare_class_files.py → class-comparison.json | Full HyRec-2 file, Green's-function energy, FIRAS branching table |
| hyrec-reference-extract.json | 60-row attributed sample of output_xe.dat |
| THIRD-PARTY-NOTICES.txt | Attribution and licence status of the extract |

## Units and conventions

- Redshift z; the solver integrates in s = ln(1+z) with backward Euler and a
  step of 5e-4, from Saha at z = 1800 down to z = 50.
- SI units inside the scripts: m^-3, s^-1, J. Spectra in kJy/sr, where
  1 kJy/sr = 1e-23 W m^-2 Hz^-1 sr^-1. x = h nu / (k T); x = 1 is 56.79 GHz today.
- Background: h = 0.6735837, T0 = 2.7255 K, Omega_b = 0.0494143,
  Omega_cb = 0.3124208, Y_He = 0.245, N_eff = 3.046, matching
  external/HyRec2020/input.dat. The massive neutrino is treated as massless.
- Hydrogen only. Matter temperature equals photon temperature throughout.
  Reionization is not included in the optical depth.
- Energy release: Gaussian in ln(1+z) of width 0.15; amplitude is the total
  Delta rho / rho_gamma.
- Deposition shares (Chen and Kamionkowski 2004, eq. 14): ionization (1-x)/3,
  excitation (1-x)/3, heat (1+2x)/3.
- Distortions: mu = 1.401 * integral of J * J_mu * heat; y = 1/4 * integral of
  J_y * heat; Delta T / T = 1/4 * integral of (1 - J) * heat; J = exp[-(z/z_th)^2.5]
  with z_th = 1.973e6 from the scaling in CLASS source/distortions.c; J_mu and
  J_y from Chluba (2013) eq. 5.

## The excitation term: why 1 - C

C is the probability that an atom in n = 2 reaches the ground state before it
is photo-ionized (the Peebles factor). An extra excitation to n = 2 therefore
produces a net ionization with probability 1 - C. The extra ionization rate per
hydrogen nucleus from excitation energy deposited at rate eps_exc per nucleus is

    (1 - C) * eps_exc / E_alpha,  E_alpha = 10.2 eV,

which adds to the ionization rate and so makes x larger. Chen and Kamionkowski
(2004) print this factor as (C - 1) in their eq. 21; read with their eq. 15
exactly as printed (dx/dz = [R_s - I_s - I_X]/((1+z)H)), that sign would make
extra Lyman-alpha photons lower x. The model uses 1 - C from the probability
argument above. This is our derivation; we did not contact the authors.

## Known errors and failure cases

- Three-level atom vs HyRec-2: at most 3.9 % for sampled z <= 1500;
  8.8 % at z = 50 over the full 50-1600 range (most likely the equal-temperature
  assumption; not isolated).
- Optical depth reaches 1 at z = 1092.0 (three-level) vs 1089.8 (HyRec-2 file).
- Step halving: Delta z(tau = 1) = 0.078; relative change of x(1100) = -5.2e-4.
- Explicit Euler with step 5e-3 gives nan by z = 1500 (stiff equation).
- The textbook root x = (-S + sqrt(S^2 + 4S))/2 loses digits of 1 - x when S
  is large (about 0.7 % at z = 2400). The scripts use the equivalent form
  x = 2 / (1 + sqrt(1 + 4/S)).
- Starting from Saha at z = 1300 gives x(1100) = 0.105 instead of 0.142.
- The soft branching fits J_mu + J_y reach 1.15 near z = 5e4, so mu and y can
  overstate distortion energy by up to 15 % there. They also differ strongly
  from the FIRAS Gram-Schmidt table between z = 1e4 and 1e5 (different
  definition of the split; see the tutorial).
- Green's-function energy check uses only the tabulated G_th column; it closes
  to 0.1 % for z <= 2.6e5 and drifts by up to 16 % above.

## Licence and attribution

Scripts and outputs written for Sandbox Physics. The HyRec-2 values come from
the frozen output file distributed inside CLASS v3.3.4 (HyRec-2: Y. Ali-Haimoud
and C. M. Hirata, with N. Lee; CLASS: J. Lesgourgues and collaborators). No
licence file was found in the CLASS v3.3.4 tree, so only a 60-row extract is
redistributed; see THIRD-PARTY-NOTICES.txt. Download the original files from
the CLASS release.