API index

adiabatic_cooling_rate(c, rec, z)

Fractional heating (dQ/dz)/ργ from baryon cooling – a negative source. The baryons cool adiabatically as (1+z)² and, while Compton coupling holds, drain heat from the photons: the one guaranteed source of the opposite sign to Silk damping.

(dQ/dz)/ργ = -(3/2) N_tot k_B T_γ / ργ · 1/(1+z),   N_tot = nₑ + n_H + n_He

damping_scale(c, rec, a)

Photon diffusion (Silk) wavenumber kD, in CLASS's kD = 2π/rd normalisation (envelope e^{-2k²/kD²}):

k_D⁻² = ∫₀^a da'/(a'³ H nₑ σ_T) · [16/15 + R²/(1+R)]/[6(1+R)]

Checked against CLASS's r_d (indexthr_d): agrees to <0.5% across 10³ < z < 10⁶.

distortion_amplitudes(c, rec; silk = true, cooling = true, nuclear = true)

Compute all spectral-distortion amplitudes at once by convolving the heating history with the exact branching ratios:

X = ∫ dz  J_X(z)  (dQ/dz)/ργ

Returns a NamedTuple (g, y, μ, residual) where residual is the vector of PCA residual amplitudes μ_k. The integral runs over the redshift range of the branching-ratio table (z ≈ 10³–5×10⁶), which brackets the μ and y eras.

distortion_spectrum(c, rec, ν_GHz; kws...)

The spectral distortion ΔI(ν) in 10⁻²⁶ W/(m² Hz sr), summing every component:

ΔI(ν) = g·G_T(ν) + y·Y_SZ(ν) + μ·M_μ(ν) + Σ_k μ_k S_k(ν)

The residual sum Σ μk Sk is the r-type distortion – the part of the signal that is neither a pure μ, y, nor temperature shift, and which distinguishes energy- release histories rather than just totals. ν_GHz may be a scalar or array.

firas_chi2(c, rec; kws...)

χ² of the predicted distortion against the COBE/FIRAS null result (Fixsen et al. 1996): Σi [ΔI(νi)/σ_i]². A ΛCDM distortion is ~4 orders of magnitude below the FIRAS noise, so this is ≪ 1 – FIRAS does not constrain standard cosmology, which is exactly why a next-generation spectrometer (PIXIE) is the goal.

Temperature-shift (g-type) amplitude. A monopole shift, absorbed into T_CMB.

heating_rate(c, rec, z; silk = true, cooling = true, nuclear = true)

Total fractional heating rate (dQ/dz)/ρ_γ, summing the guaranteed sources.

silk_heating_rate(c, rec, z)

Fractional heating rate (dQ/dz)/ργ from Silk damping of acoustic modes, reproducing CLASS's `noninjectionrateacousticdiss` integrand exactly (Chluba, Khatri & Sunyaev 2012):

(dQ/dz)/ργ = ∫ dk  4 A²  k  Δ²_ℛ(k)  e^{-2k²/k_D²}  · 𝒟(z)

A    = 1 / (1 + 4f_ν/15)                        neutrino anisotropic-stress factor
Δ²_ℛ = A_s (k/k_pivot)^(n_s-1)                  dimensionless primordial spectrum
𝒟(z) = 1/(H κ̇) · [16/15 + R²/(1+R)] / [6(1+R)]   how fast the damping scale sweeps

kD grows as the diffusion length, so the heating at redshift z probes P(k) at k ≈ kD(z). In the μ era that is k ~ 10²–10⁴ /Mpc – scales no other observable reaches, which is why a μ measurement constrains inflation.

y (Compton) distortion from the pre-recombination sources. FIRAS |y| < 1.5e-5.

μ chemical-potential distortion. Planck-cosmology ΛCDM ≈ 2e-8; FIRAS |μ| < 9e-5.

H_dc(x)

Effective double-Compton correction factor (Chluba & Sunyaev 2012, eq. 13):

H_dc(x) ≈ e^{-2x} [ 1 + (3/2)x + (29/24)x² + (11/16)x³ + (5/12)x⁴ ]

H_dc(0) = 1; it carries the entire frequency dependence of the DC emissivity beyond the soft-photon limit, and falls as x⁴e^{-2x} at large x.

Kept as the published reference point; the solver itself uses g_dc_exact, which replaces this fit (and its 1/(1+14.16θ) companion) with the full tree-level QED result.

K_emission(c, rec, z, x)

Photon-production coefficient K(x) = KDC + KBR (Chluba & Sunyaev 2012, §2.2.2), which enters the photon Boltzmann equation as

dn/dτ |_{DC+BR} = (1/x³) [ 1 - n(e^{x_e} - 1) ] · K(x)

with τ the Compton optical depth (dτ = σ_T nₑ c dt). Dimensionless.

Double Compton (their eqs. 10, 11 give the structure; the factor itself is exact):

K_DC = (4α/3π) θ² · g_dc_exact(x, θ)

where gdcexact is the tabulated exact tree-level QED result (see above) whose soft limit is I₄^pl = ∫x⁴ nPl(nPl+1)dx = 4π⁴/15 ≈ 25.976; CS2012's [I₄^pl/(1+14.16θ)]·H_dc(x) is its fitted stand-in, kept only as a reference.

That integral is the seed-photon phase space: double Compton emission is stimulated by the ambient blackbody, so the rate carries a factor of ~26. Reading KDC = (4α/3π)θ² and setting gdc = 1 – which is what the headline formula invites – underestimates photon production 26-fold and throws the thermalization redshift badly off. The 14.16θ term is the leading relativistic correction and only matters for z ≳ few×10⁶.

Bremsstrahlung (their eq. 14), fully ionized H+He:

K_BR = [ α λ_e³ θ^{-7/2} e^{-x} / (2π √(6π)) ] · Σᵢ Zᵢ² Nᵢ · g_ff

with λe = h/(me c) the unreduced Compton wavelength. DC dominates for z ≳ few×10⁵; BR takes over below that and is what keeps the low-frequency spectrum a blackbody down to recombination.

The ⁷Li/⁷Be trace ions are not folded into ΣZ²N: even taking every such nucleus at its largest ionic charge gives ΣLi,Be Z²N/(ΣH,He Z²N) < 5×10⁻⁹ for the network abundances. Their capture energy is retained by nuclear_deposition_rate; their bounded free–free opacity is numerically irrelevant at the stated precision.

chluba_J_y(z), chluba_J_mu(z), chluba_J_vis(z)

Chluba (2013), MNRAS 434, 352, eqs. (5)-(6): fits to the exact thermalization Green's function, and the reference this module is validated against.

J_y(z)  = [ 1 + ((1+z)/6.0e4)^2.58 ]⁻¹
J_μ(z)  = 1 - exp[ -((1+z)/5.8e4)^1.88 ]
𝒥(z)    = exp[ -(z/z_μ)^{5/2} ]            distortion visibility, z_μ ≈ 1.98e6

G_th(ν,z) ≈ 1.401 J_μ(z) 𝒥(z) M(ν) + (J_y/4) Y_SZ(ν) + ((1-𝒥)/4) G(ν)

so the coefficient of M – which is exactly what [distortion_from_spectrum] returns as μ – is 1.401·Jμ·𝒥, and the coefficient of YSZ is J_y/4. The two cross at z ≈ 5.3e4, which is what makes 5e4 the conventional μ/y boundary.

Why not validate against CLASS's FIRAS_branching_ratios.dat? Because those are not the physical branching ratios. They are a PCA decomposition in a basis built from FIRAS's 43 frequency channels, and Chluba & Jeong (2014) say outright that "these ratios are not unique but depend on the experimental settings". Their J_μ reaches 1.896 at z = 1e5, which is above the photon-number ceiling of 1.401 and is a projection coefficient, not an energy branching. Chluba's own physical value there is 1.312; this solver gets 1.354.

distortion_from_spectrum(grid, n_final)

Project the evolved distortion Δn = nfinal - nPlanck onto the y and μ spectral shapes to extract the amplitudes, using the photon-number and energy weights that define them:

y = ⟨Δn, Y⟩ / ⟨Y, Y⟩,   μ = ⟨Δn, M⟩ / ⟨M, M⟩

with the same G/M/Y basis as the Green's-function method. Returns (μ, y, Δρ_ρ).

g_dc_exact(x, θ)

Exact double-Compton emission factor gdc(x, θ) (see the table notes above); gdc_exact(x→0, θ→0) → I₄^pl = 4π⁴/15.

solve_thermalization(grid; z_start = 5e6, z_end = 200, sources...)

Evolve the photon distortion from z_start (blackbody; thermalization complete) down to z_end, injecting the guaranteed heating, and return the solution.

Stiff (the Compton rate dwarfs H at high z), so an implicit solver is used. The Jacobian is tridiagonal plus a rank-1 term (the electron temperature is an integral over the whole spectrum), so it is taken dense.

thermalization_distortions(c, rec; kws...)

Convenience: build the grid, solve the thermalization equation, and return the distortion amplitudes (μ, y, Δρ_ρ) derived from the evolved spectrum.

thermalization_grid(c, rec; xmin = 1e-5, xmax = 50, nx = 320)

The frequency grid.

xmin must sit below the critical frequency xc = √(Λ/κ) at which photon production balances Compton scattering. xc ≈ 4e-3 - 9e-3 across the μ era, so the default sits an order of magnitude under it. Below xc the emission term physically pins the spectrum to a blackbody; if the grid ends near xc instead, the boundary does the pinning by fiat and acts as a huge artificial photon sink that over-thermalizes μ toward zero.

Pushing xmin far lower does not help: the emission rate grows as 1/x², so xmin = 1e-4 makes the system stiff enough (~1e7 × H) to break the implicit solver on the initial transient, while the answer is already converged by xmin = 1e-3.

Thermal detailed-balance residual; should vanish up to quadrature error.

Kinematic output-frequency interval for one electron momentum p0.

Minimum incident electron momentum able to scatter omega0 to omega.

single_electron_compton_kernel(omega0, omega, p0)

Exact redistribution density P(omega0 -> omega, p0) per unit omega and per unit Thomson optical depth for an isotropic monoenergetic electron population.

thermal_compton_kernel(omega0, omega, theta; rtol=1e-10)

Exact Klein-Nishina kernel folded with the normalized relativistic Maxwell-Juttner electron distribution at theta=kT_e/(m_e*c^2).

Cosmic

A cosmology library in pure Julia, following Baumann's Cosmology and Dodelson & Schmidt's Modern Cosmology.

The package is built in layers, each resting on the one below:

1. constants.jl – CODATA constants, internal Mpc/c=1 unit system
2. species.jl – energy components: γ, ν (massless and massive), CDM, baryons, curvature, dark energy with arbitrary w(a)
3. background.jl – E(a), distances, times, conformal time
4. thermodynamics.jl – recombination, x_e(z), optical depth, visibility
5. perturbations.jl – the Boltzmann hierarchy
6. observables.jl – P(k), σ₈, growth, CMB C_ℓ

The organising idea is that a cosmology is a list of species, each of which knows only its own ρ(a). Everything above the species layer is generic in the composition, so a new component costs one method rather than a new struct in every file.