An interactive teaching resource on how neutral hydrogen between a distant quasar and us becomes the Lyman-α forest: atoms meet light, light stretches, gas maps to velocity, optical depth adds, transmitted light multiplies, and the forest can be read backward only up to a point. These notes list, for every scientific claim the scenes make, its sources, the assumptions under which it holds, the simplifications made for teaching, and how the browser's numbers were checked.
No hydrodynamics. The gas is a toy Zel'dovich volume The toy Zel'dovich volume (synthetic teaching data), not an evolved simulation: temperature follows a power-law temperature–density relation and the velocities follow linear theory; nothing is computed from a thermal or dynamical history.
Neutral hydrogen is set by photoionisation equilibrium in a uniform ultraviolet background Photoionisation equilibrium and the neutral-amount control. There is no transfer of ionising radiation, no self-shielding and no collisional ionisation; the scene that changes the neutral amount holds density and temperature fixed.
Only Lyman-α is drawn, except in the last scene, where Lyman-β and a heavier-ion line appear as extra information. There is no continuum fitting, no quasar proximity effect, and the fine structure of Lyman-α is below every width shown Lyα (H I 1s → 2p) transition data.
Independent calculations check the browser's physics: SciPy's Faddeeva function for the Voigt profile, exact-Voigt optical-depth sums, numerical integration of distances, and the fake_spectra code for the real-space method. Their results are listed under Accuracy.
Each claim states what the app teaches, its sources at their exact location, its assumptions and validity, the simplifications made for teaching, and how its numbers are checked. Read in order, the chapters follow the story of the app: one atom, one line, the expanding universe, many contributions, the gas itself, and the inverse problem.
1 · Lyα: one atom and the colour it absorbs
What a single hydrogen atom does to the quasar's light: the Lyα transition, resonance in the atom's own frame, and resonant scattering out of the beam.
Neutral hydrogen in its ground state resonantly absorbs Lyα photons at a rest (vacuum) wavelength of 1215.67 Å, with absorption oscillator strength f ≈ 0.4164; the 2p level decays back to 1s with Einstein coefficient A21 ≈ 6.26×10⁸ s⁻¹, i.e. a lifetime ≈ 1.6 ns.
λ0 = 1215.6701 Å is NIST's unresolved 1s–2 multiplet value (its upper level includes 2s); the gf-weighted 2p centre is 1215.67004 Å (0.015 km/s different). The ~1.3 km/s fine-structure splitting (1215.6736 / 1215.6682 Å) is below every width the resource shows and is not taught.
Assumptions
vacuum wavelengths; fine structure unresolved; ground-state population only
Valid for
all IGM conditions considered (A₂₁ far exceeds the collisional and photo-excitation rates at IGM densities, so n = 2 is essentially empty)
¶ Resonance in the atom's frame (Doppler condition)
An atom moving with line-of-sight velocity v sees incoming light Doppler-shifted; it is resonant with light whose observed (lab-frame) wavelength is λ0(1 + v/c) to first order in v/c. An atom moving away from the observer (v > 0, toward the quasar) therefore absorbs a redder colour.
¶ Natural broadening from the finite lifetime of 2p
Because the 2p state has a finite lifetime (1/A21 ≈ 1.6 ns), even a stationary atom's response is a Lorentzian with half-width γ = Γλ0/4π ≈ 0.0061 km/s, Γ = A21 (2p decays only to 1s). This is 1765× narrower (FWHM) than the thermal width at 10⁴ K, but falls off only as 1/Δv².
γ=4πΓλ0=6.06×10−3kms−1,a=γ/b
Convention
Γ is the total decay rate out of the upper level (the lower level is the ground state, width 0)
Assumptions
radiative damping only; collisional (pressure) broadening negligible at IGM densities
Valid for
IGM densities (n ≲ 10⁻² cm⁻³)
Simplifications
The damped oscillation (microscope only, labelled) is the classical Lorentz-oscillator analogue; the physical statement is the finite lifetime.
¶ Photoionisation equilibrium and the neutral-amount control
In the highly ionised IGM, photoionisation balances recombination: n_HI Γ_HI = α(T) n_e n_p, so n_HI ≈ α(T) n_e n_H / Γ_HI ∝ Δ² T^(−0.7) / Γ_HI. Changing the UV background rescales n_HI (and τ) by 1/Γ_HI; it does not move lines or change their τ-shape.
2 · The line profile: thermal core and natural wings
How motion and the excited state's short life shape one absorption line — the Gaussian census, the Lorentzian wings, their Voigt combination and how the browser computes it.
¶ Maxwellian velocities; one line-of-sight component is Gaussian
In gas at temperature T each Cartesian velocity component of the atoms follows a Gaussian with variance kT/m. The beam is sensitive only to the component along it.
f(v∥)=2πkTme−mv∥2/2kT
Convention
v_∥ along the sightline, positive away from the observer
Assumptions
thermal equilibrium of the atoms' motions; no turbulence
Valid for
IGM gas (collisional thermalisation of velocities is fast compared with dynamical times)
¶ Doppler parameter and the thermal (Gaussian) profile
The census of line-of-sight speeds is the line's thermal profile: a Gaussian in Δv of width b = √(2kT/m), normalised so that heating spreads the same absorption over more velocities without adding any. b(10⁴ K) = 12.85 km/s for hydrogen.
b=2kT/mH+bturb2,ϕG(Δv)=πbe−Δv2/b2
Convention
m = mass of the ¹H atom (1.00782503 u = 1.6735328×10⁻²⁴ g; an earlier value, 1.6735575×10⁻²⁴ g, was corrected — a 1.5×10⁻⁵ effect); b_turb = 0 unless stated. (fake_spectra's mass convention is recorded in fake_spectra conventions (the validation reference).)
¶ The Voigt profile — normalisation, definition and asymptotics
The line profile is the convolution of the thermal Gaussian with the natural Lorentzian (each atom's Lorentzian at each atom's Doppler shift). In velocity units φ(Δv) = H(a, Δv/b)/(√π b), with H(a,x) = Re w(x + ia) and ∫H dx = √π. For a ≪ 1 the core is ≈ e^(−x²) and the far wings ≈ a/(√π x²): core from motion, wings from lifetime.
The app evaluates H(a,x) with Humlicek's (1982) W4 rational approximation at all (a, x). Against the Faddeeva function (SciPy's wofz) it is accurate to < 10⁻⁴ (relative) over a ∈ [10⁻⁶, 1], |x| ≤ 2×10⁴ (measured values in Voigt function accuracy).
Valid for
a ∈ [1e-6, 1], |x| ≤ 2e4
Sources
Humlíček 1982 — algorithm W4 (regions I–IV): all 29 coefficients and 4 border constants are identical, digit by digit, in five independent open-source transcriptions that cite Humlicek 1982 (RFM, ExoCross, SKIRT 9, SOCRATES, monoRTM). The original paper was not consulted; the accuracy authority is SciPy's Faddeeva function (wofz), against which Voigt function accuracy tests every region down to a = 10⁻⁶; the regime where Zaghloul (2015) reports W4 failing (a < 10⁻⁷) lies below every damping parameter used here
Zaghloul 2015 — Table 1 (original borders of the regions), p. 3, arXiv:1505.05596v1 — matches the app's region borders(secondary)
Once τ0 ≳ 3 the line core transmits almost nothing and the line's depth stops measuring N; its equivalent width then grows only logarithmically (∝ b√ln τ0) until the damping wings pass τ = 1, after which W ∝ √(N γ). In τ the profile shape never changes; only in F does the line change character.
W∝N(τ0≪1);W∝blnτ0(flat);W∝Nγ(damped)
Assumptions
single Voigt line
Valid for
isolated lines
Sources
Draine 2011 — ch. 9 Absorption Lines: The Curve of Growth (p. 75 ff.)
From the Lorentzian asymptote, τ(Δv) ≈ Nσ_vγ/(πΔv²), so the wing reaches τ = 1 at Δv ≈ √(Nσ_vγ/π) ≈ 160 km/s (N/10¹⁹ cm⁻²)^½, independent of b (this algebra is ours, from the far-wing profile in Meiksin 2009). For forest lines (N ≲ 10¹⁵) this lies inside the Doppler core: their wings are invisible. The wings become detectable only in rare absorbers of very high column density (high N_HI).
Δvτ=1≈Nσvγ/π
Assumptions
a ≪ 1
Valid for
damping wings become detectable from N ≳ 10^18.5 (≈ 3×10¹⁸) cm⁻² and set the line's width at τ = 1 above about 10¹⁹ cm⁻² (the sources below); the asymptote assumes a ≪ 1
Sources
Meiksin 2009 — §II.A, Eqs. (6), (12)–(13) and the text before Eq. (26) (arXiv v3 p. 9): far from line centre φ_V ≃ a/(πx²); §II.B.1 (p. 9): forest systems are well fit by Doppler profiles, DLAs show the damping wings
Rauch 1998 — §3.4 (arXiv v1 p. 15): 'From N ≳ 10^18.5 on, the damping wings of Lyα become detectable'
Wolfe, Gawiser & Prochaska 2005 — §1.1 (arXiv v2 p. 4): unit optical depth in the damping wings at Δν_τ=1 ∝ [A21 f21 N(H I)]^1/2 when N(H I) > 10^19 cm⁻²(secondary)
How it is checked
Optical depth: 164 vs 161 km/s (10¹⁹), 720.7 vs 720.0 km/s (2×10²⁰)
3 · Cosmological coordinates: redshift, distance and velocity space
How the expanding universe stretches the light, how real-space distance and peculiar motion map to observed velocity and wavelength, and the conventions used for both.
¶ Cosmological redshift stretches every wavelength by the same factor
Light absorbed by gas at redshift z_abs at rest-frame wavelength λ_rest (1215.67 Å for Lyα, in the gas's own frame) reaches us at λ_obs = λ_rest (1 + z_abs), where 1 + z_abs = a(t_obs)/a(t_abs). Every wavelength, and every wavelength interval, is multiplied by the same factor. The reverse operation, λ_rest = λ_obs/(1 + z_abs), re-expresses the observed spectrum in the absorber's rest frame: there the absorber's own line is at exactly 1215.67 Å. (In the quasar's rest frame, λ_obs/(1 + z_q), the same absorption lies blueward of 1215.67 Å — The quasar's frame — why the forest lies blueward of the quasar's own Lyα.)
Three frames are kept distinct: the absorber's rest frame (where the atom sees 1215.67 Å), the observer's frame (us), and the quasar's rest frame (The quasar's frame — why the forest lies blueward of the quasar's own Lyα). "Measured at a point along the path" means by an observer at rest with the cosmic expansion there (comoving), for whom λ ∝ a.
Assumptions
FLRW expansion; comoving observers; the gas's own small motion folded into z_abs here
Valid for
any z; the beat uses 0 ≤ z_abs < z_q = 3.2 with z = 3 as the default (the toy's redshift)
Simplifications
Redshift is presented as the stretching of the light's wavelength with the expansion (λ ∝ a); the Doppler-vs-expansion interpretation debate is not raised.
Sources
Hogg 1999 — §3 'Redshift', Eq. (8): z = λ_o/λ_e − 1; Eq. (12): 1 + z = a(t_o)/a(t_e)
How it is checked
frame conversions computed directly from 1 + z = λ_obs/λ_rest (science/oracle/redshift_oracle.py)
¶ The quasar's frame — why the forest lies blueward of the quasar's own Lyα
The quasar at z_q emits its own Lyα at 1215.67 Å in its frame, observed at 1215.67 (1 + z_q) Å. Gas between us and the quasar (z_abs < z_q) absorbs at 1215.67 (1 + z_abs) Å observed, i.e. at 1215.67 (1 + z_abs)/(1 + z_q) < 1215.67 Å in the quasar's rest frame: the Lyα forest lies blueward of the quasar's Lyα emission.
¶ On a logarithmic ruler a stretch is a slide; small slides are velocities
Because the stretch multiplies every wavelength by (1 + z), it adds the same amount, ln(1 + z), to every ln λ. Small differences in ln λ are velocities, Δu = c Δln λ; this is the velocity coordinate u used from Beat 6 on, with λ_obs = λ0 (1 + z0) e^(u/c) (Velocity coordinate ↔ observed wavelength).
¶ Comoving distance along the light's path (Planck 2018)
The line-of-sight comoving distance to redshift z is χ(z) = (c/H0) ∫0^z dz'/E(z'), E = (Ωm(1+z)³ + ΩΛ)^½ for flat ΛCDM; with Planck 2018 parameters χ(3) ≈ 4384 Mpc/h. It places the gas and the travelling light along the drawn path in Beat 5, and gives the distance to us printed beside the volume (Beats 0, 1, 11, 12), where the beam is broken because the 20 Mpc/h sample is about 0.5 % of that distance.
χ(z)=H0c∫0zE(z′)dz′
Convention
Planck 2018 TT,TE,EE+lowE+lensing (Ωm = 0.3153, ΩΛ = 0.6847, h = 0.6736); radiation neglected (≲ 10⁻³ of E at z ≲ 4)
Assumptions
flat ΛCDM
Valid for
0 ≤ z ≤ 6 tested
Simplifications
The path is drawn to scale in comoving distance, so equal stretches do not cover equal lengths of path.
As the quasar's light travels, its wavelength measured by local comoving observers grows, λ_here = λ_emit (1 + z_q)/(1 + z_here). Each intervening hydrogen structure at z_abs absorbs the part of the spectrum that is at 1215.67 Å in its own frame — the part emitted at λ_emit = 1215.67 (1 + z_abs)/(1 + z_q) — and that shadow then stretches on with the light, reaching us at 1215.67 (1 + z_abs). Structures at different redshifts therefore cast shadows at different observed wavelengths, ordered by redshift, all blueward of the quasar's own Lyα at 1215.67 (1 + z_q): the Lyα forest.
Lyβ and higher Lyman lines, the Lyman limit, and the atmospheric cut-off are not drawn.
Sources
Meiksin 2009 — §II.A 'Resonance absorption lines', text preceding Eq. (15) (arXiv v3 p. 7): absorption appears at observed wavelengths between λ_lu and λ_lu a(t)/a(t₀)
Gas at comoving position x along the sightline, with line-of-sight peculiar velocity v_pec, absorbs at velocity coordinate u = aH(z) x_com + v_pec. Across a box of tens of Mpc/h the expansion rate is treated as constant.
Positions are comoving (Mpc/h); a column density uses the proper path length Δx_prop = a·Δx_com/h (in cm); velocities are physical km/s; densities are proper cm⁻³. The velocity per comoving Mpc/h, aH(z)/h = 100E(z)/(1+z), does not depend on h.
v_pec > 0 means the gas moves away from the observer (toward the quasar, +x); it lands at larger u, i.e. at a longer observed wavelength (redder). This is the same sign as the resonance condition Resonance in the atom's frame (Doppler condition).
vpec>0⇒u↑⇒λobs↑
Convention
as stated; used identically by the app, the independent checks and the toy generator
The velocity ruler is converted to observed wavelength with λ_obs = λ0(1+z0)e^(u/c) (the log-λ convention); to first order this equals λ0(1+z0)(1+u/c). Over 20 Mpc/h at z = 3 the two differ by 2.9e-05 in 1+z.
λobs=λ0(1+z0)eu/c
Convention
z0 is the redshift of the u = 0 point (the near end of the sightline)
Valid for
|u| ≪ c
Sources
Karaçaylı et al. 2024 — §3.2 'Quadratic estimator', Eq. (14) (arXiv v3 p. 6): v = c ln(λ/λ_Lyα), 'following the convention in cosmology' — the same convention, with this resource's zero point at λ0(1+z0)
Treating aH as constant across the box deviates from the exact distance–redshift relation by 1e-04 per Mpc/h (0.2% over 20 Mpc/h, 3.9 km/s at z = 3; 1.0% at 100 Mpc/h).
¶ Velocity-gradient (Hubble) broadening of an extended absorber
Gas that is extended along the line of sight spreads its absorption in velocity: the Hubble flow across its depth, and any internal peculiar-velocity gradient, give its parts different line-of-sight velocities, and this bulk spread adds to the thermal motion of its atoms. For a Gaussian clump of width σ_x with a linear velocity field, the spread is a Gaussian of dispersion |aH + dv/dx|·σ_x, i.e. a Doppler parameter b_size = √2·|aH + dv/dx|·σ_x, which adds in quadrature to the thermal b; this expression is the resource's own algebra (a Gaussian velocity spread convolved with the thermal Gaussian), checked analytically and by One parcel, four separate variables — the sources give the physics, not the formula. In CDM collapse models the low-column-density forest clouds are cool but large, so this term can dominate their widths.
Rauch 1998 — §3.3 (arXiv v1 p. 14): in CDM collapse models, low-column clouds are cool but large, 'so that the Hubble expansion dominates the line broadening'
Meiksin 2009 — §II.A, Eqs. (10)–(11): b = (2kT/m_a)^½; a kinematic component adds in quadrature
How it is checked
a Gaussian velocity spread convolved with a Gaussian thermal profile is a Gaussian whose b² is the sum; One parcel, four separate variables checks the drawn width against it
Why optical depths add while transmitted light multiplies, what saturation does, and what an instrument changes — the record, not the gas.
¶ Lyα "absorption" in the IGM is resonant scattering out of the beam
After excitation to 2p the atom decays radiatively (A21 ≫ collisional and photo-excitation rates at IGM densities), re-emitting a Lyα photon in an essentially random direction. For a narrow beam from a point source this removes the photon from the line of sight; scattering of other photons into the beam is negligible. The beam's extinction therefore obeys the same e^−τ law as true absorption.
1s+γLyα→2p→1s+γ(random direction)
Convention
Say "scattered out of our beam" in Beat 2 and again in Beat 5, where "shadow" is introduced; "removed"/"absorbed" elsewhere as shorthand.
Assumptions
point source; beam solid angle ≪ 4π; negligible re-scattering into the beam
Valid for
quasar sightlines through the diffuse IGM
Simplifications
Later beats say "absorbed" for brevity.
Sources
Meiksin 2009 — §I footnote 7; §II.A opening (features are scattering of resonance photons out of the line of sight)
An absorber of neutral column N_HI has τ(Δv) = N σ_v φ(Δv) with σ_v = (πe²/m_e c) f λ0 = 1.3435×10⁻¹² cm² km/s, so ∫τ dv = N σ_v (independent of temperature) and the Doppler line-centre value is τ0 = 0.758 (N/10¹³)(10 km/s / b).
τ(Δv)=NHIσvϕ(Δv),σv=mecπe2fλ0,τ0=πbNσv
Convention
velocity-space cross-section (cm² km/s); φ normalised to 1 in km/s
Assumptions
ground-state population only; stimulated emission negligible
¶ Line-of-sight accumulation of optical depth in velocity space
The optical depth at velocity u is the sum (integral) over all neutral hydrogen along the sightline of each element's profile centred on its own u(x): τ(u) = ∫ n_HI σ_v φ(u − u(x); b(x)) dx_prop. Discretised on cells it is Σ_j N_j σ_v φ(u − u_j; b_j).
¶ Discretisation limits — pixel size, pixel averaging and the periodic box
The cell sum is accurate only when (i) real-space pixels are finer than the narrowest line (aH·Δx ≤ b_min/2); (ii) one states whether a pixel's value is its centre τ, its mean τ or its mean F, since these differ for lines narrower than a pixel; and (iii) damping wings of strong absorbers are not cut by the periodic box.
aHΔx≤bmin/2,⟨e−τ⟩pix=e−⟨τ⟩pix
Convention
the app's default = point sampling of τ at pixel centres on a uniform u grid; nearest periodic image per cell
Valid for
Closure against fake_spectra: lines narrower than a coarse pixel are wrong by up to 0.46 in F (EW up to +93%); oversampling to aH·Δx ≤ b/2 with exact pixel averaging brings this to 0.019. The pixel-mean flux and the flux of the pixel-mean τ differ by 7e-03 to 0.31 at 10 km/s pixels (the larger for narrow saturated lines). Both codes truncate wings at half a box: for an LLS/DLA in a 5 Mpc/h box the optical depth dropped against a true periodic sum reaches 8.8. None of these occur in the toy at its sampling.
Optical depth is the path integral of the opacity; the opacity of independent absorbers (here: different gas elements at the same observed velocity) adds linearly because each atom's absorption is independent of the others and of the radiation field (no stimulated emission, no saturation of the level populations). Flux decrements do not add.
τ(u)=i∑τi(u)
Assumptions
linear (optically passive) medium; ground-state populations unaffected by the beam
With no emission into the beam along the path, the transmitted fraction of the quasar's light is F = e^(−τ). F is the observed spectrum divided by the quasar's unabsorbed (continuum) spectrum.
F(u)=e−τ(u)
Convention
continuum-normalised flux, 0 ≤ F ≤ 1
Assumptions
no source term in the beam; continuum known
Valid for
quasar sightlines
Simplifications
The quasar continuum is divided out and never drawn (continuum fitting is out of scope).
A spectrograph records the transmitted flux convolved with its line-spread function, approximated as a Gaussian (we quote its FWHM in km/s; P1D analyses often quote the rms width, e.g. Palanque-Delabrouille et al. 2013 Eq. 6). Resolution changes neither the gas nor its optical depth — only the recorded spectrum.
Fobs(u)=∫F(u′)LSF(u−u′)du′(+noise)
Convention
Gaussian LSF; FWHM in km/s
Assumptions
Gaussian
velocity-independent LSF
Valid for
illustrative; real LSFs vary with wavelength and fibre
The toy noise is white Gaussian with constant σ = 1/(S/N) per detector pixel in continuum-normalised flux; detector pixels are max(2.5 km/s, FWHM/3) wide. Any claim of "indistinguishable" is made only at a stated S/N and pixel size.
σF=(S/N)−1per pixel
Convention
S/N quoted per 2.5 km/s pixel unless stated
Assumptions
flux-independent noise; no sky
no continuum error
Valid for
illustrative
Simplifications
constant σ; Gaussian
Sources
No literature source: an idealised model chosen by this resource and declared here.
How it is checked
a declared, idealised noise model of this resource (white, Gaussian, flux-independent); the degeneracy check Same-shadow degeneracy uses it
5 · The gas and how it is represented: a continuous sightline
The intergalactic gas itself — density, ionisation, temperature — and the teaching representations (volume, atoms, parcels, cells) that stand for it, with the toy data behind them.
¶ One velocity per sampling element (unresolved velocity structure is not represented)
A real-space skewer gives each sampling element (pixel) one peculiar velocity and one temperature. If the gas inside one element has several velocities — in particle codes, several SPH particles or kernel contributions moving differently that overlap the same pixel — that sub-element velocity structure is not represented, and its absorption is placed at one velocity instead of being spread over several.
τ(u)=components∑∣du/dx∣nHIσv(cold limit)vs one u per sampling element
one mass-weighted velocity and temperature per real-space sampling element
Valid for
Measured in Closure against fake_spectra on synthetic configurations (largest difference in F): gas-like 1D forest 4.9e-03 (inside the tolerance); two velocity components at one position (C3c) 0.089; a velocity fold (C4c) 0.021; a 1D Zel'dovich field evolved past shell crossing (C10, a deliberately exaggerated stress test, not a claim about the IGM) up to 0.2. Any reduced sightline representation of this kind loses the velocity and density structure inside one sampling element.
Simplifications
The toy (The toy Zel'dovich volume (synthetic teaching data)) stores one mass-weighted, smoothed velocity per pixel (the Zel'dovich field is just past first shell crossing in places, but the deposit averages over it), so the app and the toy data agree by construction; the limitation matters only when comparing with particle-based codes.
the quoted numbers hold at z ≈ 3 and mean density for those inputs; x_HI scales as Δ (1+z)³ T^−0.7 / Γ_HI (the equation), so denser or cooler gas is more neutral
Simplifications
The teaching copy says "about one in a hundred thousand" (computed: one in 135 000–160 000).
A skewer is delivered as real-space pixels along the line of sight: x (comoving Mpc/h, uniform Δx), n_HI (proper cm⁻³), T (K), v_pec (km/s, + away from observer), with metadata z, cosmology (Ωm, ΩΛ, h), Δx, and optionally the producing code's own τ(u) for validation. The app computes N_j = n_HI,j Δx_prop and τ by Line-of-sight accumulation of optical depth in velocity space.
Interface is dataset-independent (no field names or units from any one code). Optional metadata `mass_convention` and per-pixel `b_turb` let the generating code's own τ be reproduced (the app's import options for mass and b_turb; the validation mode of Closure against fake_spectra uses fake_spectra's H mass, 1.00794 m_p, and so its thermal widths).
Sources
No literature source: checked numerically against an independent code.
¶ fake_spectra conventions (the validation reference)
fake_spectra 2.2.4 conventions for units, velocity mapping, thermal broadening mass, profile evaluation, SPH kernel integration and line data are documented with source citations in science/notes/FAKE_SPECTRA_CONVENTIONS.md, and their effect on our real-space method is quantified in Closure against fake_spectra.
Convention
documented in the source repository (science/notes/FAKE_SPECTRA_CONVENTIONS.md)
Valid for
fake_spectra 2.2.4 (PyPI sdist sha256 4cc80e6c…eb84), MIT licence
Simplifications
Differences that matter for closure: thermal width uses 1.00794 m_p (b 0.37% small; up to |ΔF| 5.5e-03 in a forest); stored τ is the pixel-averaged τ when the pixel is at least b/2 wide, the pixel-centre value otherwise (Discretisation limits — pixel size, pixel averaging and the periodic box); wings truncated at half a box; tail floor τ < 10⁻⁷. Its total difference from an exact calculation is up to 1.8e-04 in F for isolated lines, almost all from its pixel quadrature (1.8e-04); its 7-node kernel quadrature alone contributes up to 2.2e-07 there, and dominates for broad cold kernels (1.5e-02 in F).
fake_spectra (Bird) — v2.2.4 sdist (PyPI 2025-08-11): H mass 1.00794 m_p, line_data.py:21 and absorption.h:4; pixel-centre value or pixel average (the b/2 rule), singleabs.h:104-126; wings stopped at half a box, absorption.cpp:248-278; tail floor τ < 10⁻⁷, spectra.py:135 and absorption.cpp:250-278; 7-node kernel quadrature, singleabs.h:8, 143-167; Lyα data, atom.dat:7. Full notes in the source repository (science/notes/FAKE_SPECTRA_CONVENTIONS.md)
¶ Exporting skewers from fake_spectra 2.2.4 — a velocity-convention compatibility note
During compatibility testing, the velocity returned through the tested export path — fake_spectra 2.2.4, Spectra.get_velocity on a Gadget-format snapshot — differed from the convention this resource requires (the proper peculiar velocity in km/s, Comoving vs physical coordinates and h-units). On a mock snapshot at a = 0.25 with known particle velocities, the returned value was 0.5 of the true peculiar velocity in the pixel of every moving particle, equal to √a — consistent with the √a unit conversion being applied once more than this convention needs. The resource therefore never takes skewer velocities from that path: it constructs the line-of-sight velocity independently from the stored quantities (the toy from its displacement field). One other observation on the same version (2.2.4), for anyone exporting skewers under NumPy 2: get_velocity fails unless its arrays, promoted to float64, are cast back to float32, and HDF5Snapshot.get_npart overflows (2³² × uint32) on headers that store particle counts as uint32, as Gadget does — the mock snapshots store them as int64.
vreturned/vpec=aon the tested path (mock snapshot);used here: vpeccomputed independently
fake_spectra 2.2.4 (PyPI sdist), the tested path only; the √a observation also holds on the upstream master branch as inspected on 2026-10-02. If an upstream issue is filed, this note will link it.
not used by a scene · status: verified · id SCI-SIM-003
¶ The toy Zel'dovich volume (synthetic teaching data)
The 3D gas, sightline and parcels in this resource come from a declared toy, not a simulation: a 256³ Zel'dovich displacement of a Gaussian field over 20 Mpc/h (P(k) ∝ k⁻² e^(−k²R²), R = 0.075 Mpc/h, ψ_rms = 0.36 Mpc/h), linear peculiar velocities v = f aHψ with f = Ω_m(z)^0.55 (no free velocity scaling), T = T0 Δ^(γ−1) with T0 = 1.3×10⁴ K, γ = 1.5, Planck 2018 cosmology, and n_HI from Photoionisation equilibrium and the neutral-amount control. Γ_HI is tuned so ⟨F⟩ = 0.68 on the chosen sightline (Mean transmitted flux at z ≈ 3); the tuned value (2.41×10⁻¹² s⁻¹) is 3.0× the measured 0.8×10⁻¹² s⁻¹, so the toy's absolute ionisation is not a physical claim. The skewer file also carries the generating code's own τ (exact Voigt) for the reference view.
Rescaling τ to an observed mean flux is standard practice in simulation analysis (fake_spectra ships a mean-flux rescaler); here it is done through Γ_HI.
Sources
Becker et al. 2013 — Table 2: ⟨F⟩ = 0.6966 (z = 2.95), 0.6670 (z = 3.05)
n̄_H(z) = (1−Y) Ω_b ρ_crit,0 (1+z)³ / m_H ≈ 1.21×10⁻⁵ cm⁻³ at z = 3 for Planck 2018 (Ω_b h² = 0.02237, Y_P = 0.2454) — about a dozen hydrogen nuclei per cubic metre; ~80 per m³ in a filament at 6.6× the mean.
nˉH=(1−Y)Ωbρcrit,0(1+z)3/mH
Convention
ρ_crit,0 = 1.87834×10⁻²⁹ h² g cm⁻³
Valid for
the number holds at z = 3; it scales as (1+z)³ (the formula holds at any z)
At z = 3 one comoving Mpc/h corresponds to ≈ 114.2 km/s (Planck 2018), and peculiar velocities of tens of km/s (toy rms 40 km/s per component) shift gas by a few tenths of a Mpc/h in the spectrum — up to ~1 Mpc/h in infall regions.
Lyα-forest lines have N_HI ≲ 10¹⁷ cm⁻²; Lyman-limit systems N ≳ 1.59×10¹⁷ (τ at 912 Å ≥ 1, σ_912 = 6.304×10⁻¹⁸ cm²; surveys also use 10¹⁷ or 10^17.2); damped Lyα systems N ≥ 2×10²⁰ cm⁻².
The atoms drawn when zooming into a parcel are a few hundred representative neutral atoms, not a count and not simulation particles; their motions are straight-line cartoons slowed ~10¹³×, with velocity components drawn from the correct Maxwellian.
Valid for
always declared on screen
Simplifications
only neutral atoms drawn; ~10⁵ ionised nuclei per neutral atom omitted
Sources
A teaching representation chosen by this resource — described here so it is not mistaken for physics.
¶ The 3D slab is a teaching volume; only gas on the beam absorbs
The slab is a 20 × 8 × 8 Mpc/h toy volume (The toy Zel'dovich volume (synthetic teaching data); 7.97 Mpc/h across the beam) drawn as stacked translucent slices; gas near the beam is darkened (focus rule) because only gas on the line of sight contributes to τ.
Valid for
always
Simplifications
slices are downsampled 2× from the toy grid; darkness ∝ log density
the drawing fades toward the sample's edges (6 % of its length, 20 % across): a drawing convention — the sample ends, the gas does not
Sources
A teaching representation chosen by this resource — described here so it is not mistaken for physics.
¶ Parcels are a teaching decomposition of a continuous sightline
The three parcels a, b, c summarise three n_HI peak regions of the real toy sightline (HI-weighted position, width, column, temperature, velocity and gradient). They are not discrete clouds; Beat 8 dissolves them back into cells. A parcel has a size along the beam, so the velocity spread across it widens its line too (the physics is Velocity-gradient (Hubble) broadening of an extended absorber). Beat 6 isolates the declared terms by setting the parcel's own peculiar-velocity gradient to zero (dv/dx = 0 there), so its size term is the Hubble flow alone, b_size = √2·aH·σ_x; for the parcel it starts from, b_size = 23.8 km/s against a thermal 21.4 km/s, so the drawn line is 32.0 km/s wide (b) — Beat 6's figure names both parts, and "show the physics" prints each and the drawn line's own width. The same toy parcel with its own gradient (-70 km/s per Mpc/h: infall, partly cancelling the Hubble flow), as Beats 7–8 use it, has b_size = 9.1 km/s and b = 23.3 km/s.
b2≃bthermal2+bsize2,bsize=2aH+dxdvσx
Convention
computed once in science/data/toy/zeldovich3d.py (decompose)
Assumptions
Gaussian thermal and spatial profiles (the damping wings aside)
Valid for
always declared
Simplifications
Gaussian parcel shape
Sources
A teaching representation chosen by this resource — described here so it is not mistaken for physics.
How it is checked
Beat 8 reports the largest flux difference between the 3-parcel and full calculations; One parcel, four separate variables checks the width Beat 6 prints
The natural (Lorentzian) width is 0.006 km/s; wherever it must be seen it is magnified and labelled (×500 resonance window in Beat 2, ×1200 inset in Beat 4, ×1–1000 control in the micro view).
Valid for
always labelled
Sources
A teaching representation chosen by this resource — described here so it is not mistaken for physics.
¶ Ink is optical depth, and the ink rendering performs Beer–Lambert
"Ink" has one meaning throughout: a contribution to optical depth τ, laid in velocity space (Beats 6–13); neutral atoms are drawn as dark dots and motions as arrows, never called ink. The ribbon is rendered by compositing ink layers of opacity 1 − e^(−τᵢ) over a backlight C₀, which yields the transmitted fraction (C − C_ink)/(C₀ − C_ink) = Π e^(−τᵢ) = e^(−Στᵢ): the picture performs the rule it teaches (checked pixel by pixel, The drawing performs Beer–Lambert).
¶ Four representations of the same gas, never mixed silently
The resource draws the same gas in distinct representations, each labelled where it first appears: (1) the continuous gas field (density, Beats 0, 7–8, 11–12); (2) representative atoms, not to scale (Beats 1–4, and inside Beat 6's parcel); (3) teaching parcels — a summary of one stretch of the sightline (Beats 6–8; Beat 13's "clouds" are drawn the same way); (4) cells — the sightline sampled finely (Beats 8–12). Simulation particles are never drawn.
Valid for
representation rule
Sources
A teaching representation chosen by this resource — described here so it is not mistaken for physics.
¶ The stretching ribbon and the written shadows (Beat 5 representation)
In Beat 5 the quasar's light is a ribbon carrying its own labels (the wavelength each part left with), laid against a ruler of wavelength measured where the light is now, starting at 0 Å, so the stretch by (1 + z_q)/(1 + z_here) is drawn exactly (anchored at zero, every interval scaled). Lyα of the gas being passed is a fixed mark at 1215.67 Å. The gas along the line of sight is a schematic continuous field, not a set of separate clouds: a seeded Gaussian random field in ln Δ (correlation length ~7 Mpc/h comoving, σ = 0.9; far coarser than real structure, chosen to be visible on a 4.5 Gpc/h path), quieter near the quasar, with three density peaks set exactly for the narration (the selected one, placed where the toy stretch's deepest absorption falls on the forest's own wavelength scale, so its shadow has the same wavelength in Beats 5, 11 and 12; and z = 2.90, 2.81). Every cell writes at its own Lyα with strength ∝ Δ² (photoionisation equilibrium, Photoionisation equilibrium and the neutral-amount control) × ((1+z)/4)^1.6, normalised so ⟨F⟩ = 0.74 over z = 2–3; the forest therefore thins toward low z. Density peaks (Δ > 2.2, with their surroundings Δ > 1.3) are the "structures" the narration names and the reader can select. What each cell writes is drawn ~60× wider than a real forest line (6 Å Doppler width when written, versus ~0.1 Å for b = 25 km/s); shadows combine by Beer–Lambert. The close-up of the mark is magnified ~20×. The surface spectrum is continuum-normalised; the quasar's Lyα emission (a schematic bump) is drawn only with "show the physics". The light packet's drawn wavelength is proportional to its true stretch. Path positions are comoving distances (Comoving distance along the light's path (Planck 2018)).
No colours are painted on the ribbon; ultraviolet and visible light are not distinguished.
The close-up keeps a pencil trail of the drift of the quasar's own Lyα (and a cinnabar one for the selected shadow) from the Lyα mark — a record of motion, drawn only while it fits the close-up.
Sources
A teaching representation chosen by this resource — described here so it is not mistaken for physics.
A constructed degeneracy example, not a model of typical forest absorbers: (i) one unusually hot, idealised compact gas clump (T = 4.5×10⁴ K, Hubble width 2.23 km/s — about 7 kpc proper, far more compact than typical forest gas — N = 2.5×10¹⁴, at rest); (ii) two cold clouds (T = 10⁴ K, Hubble width 12.8 km/s each, 78 km/s apart in Hubble flow, 9.75×10¹³ each, approaching at ±24.9 km/s); (iii) one cold cloud stretched along the line (T = 10⁴ K, Hubble width 11.15 km/s, N = 2.5×10¹⁴) whose smooth extra velocity across 1σ (5.988 km/s) is computed from (i), not tuned: it is the value that gives (iii) the same total Doppler parameter as (i), b_tot² = 2kT/m_H + 2σ_u² = (27.431 km/s)², with the same column. At S/N = 20 per 2.5 km/s pixel with perfect resolution the three Lyα lines are not significantly distinguishable (Δχ² ≈ 2.8 for the blend over 121 pixels); the blend shows by S/N = 50 (Δχ² ≈ 17) and grows as (S/N)², while (iii) stays at Δχ² ≤ 0.002 at every S/N up to 200, in Lyα and in Lyβ. In this constructed model (i) and (iii) are degenerate in every H I Lyman-series line: for a Gaussian clump with a smooth (linear) velocity field the H I profile is a Voigt profile fixed by b_tot and N alone. These are two different kinds of degeneracy. The blend (ii) is measurement-limited: higher S/N, resolution or wavelength coverage can break it. Heat against smooth motion, (i) against (iii), is information-limited for hydrogen: every H I line depends only on the line-of-sight distribution of velocity and optical depth, in which thermal and smooth bulk broadening of one species are the same Gaussian, so another H I line gives no independent distinction — a line of a different mass does (What breaks the degeneracy). The blend hides because its two components overlap inside the saturated Lyα core, where the transmitted flux is near zero however the optical depth is split (Saturation and the curve of growth); Lyβ, with a smaller f λ, is not saturated there and separates it (What breaks the degeneracy). What remains of (iii)'s difference is numerical: the clumps are cut at ±5σ (exact-Voigt path, max |ΔF| 4e-06), and the app adds its Voigt approximation's own error (max |ΔF| 4e-05; Numerical Voigt implementation in the app).
noise-free Δχ² against configuration (i), σ_F = 1/(S/N) per 2.5 km/s pixel; configurations defined in velocity space (single source science/validation/degeneracy/configs.json, written by science/tools/same_shadow_configs.py), so they do not depend on the cosmology
(ii) against (i) — only at the stated S/N and resolution; (i) against (iii) — any S/N and resolution, in H I absorption alone (every Lyman-series line), for this constructed model's Gaussian clumps with smooth (linear) velocity fields
Simplifications
(i) and (ii) are chosen by hand; (iii) is computed from (i) by the matching condition
idealised Gaussian clumps, not typical forest absorbers
Sources
Rauch et al. 1996 — §3, Eq. (1) (arXiv p. 5): b_i² = b_b² + 2kT/m_i — a bulk width common to all species and a thermal width scaling as m^(−½); the paper separates temperature from bulk motion with ions of different mass — with lines of one mass, Eq. (1) cannot separate T from b_b (our inference)
Meiksin 2009 — §II.A, Eqs. (10)–(11): b = (2kT/m_a)^½; a kinematic component adds in quadrature
How it is checked
independent recomputation with wofz (science/oracle/degeneracy_oracle.py)
Thermal width scales as m^(−½) while bulk motion broadens all species alike, so a heavier ion's line in the same gas separates heat from motion (Si IV: m ≈ 28 m_H). Lyβ (f λ 6.24× smaller) saturates differently and can, in some cases, expose columns hidden in saturated Lyα.
the gas holds silicon in this ion — the tie-breaker applies only if it does
the metal line traces the same gas (single phase)
a common non-thermal (bulk) width shared by both species — an assumption that need not hold everywhere (Rauch et al. 1996), and a stronger one for H I with a metal ion than for two metal ions
illustrative metal column
Valid for
Illustrative. Si IV separates heat from motion only where silicon is present in the same gas. Lyβ helps only in some cases: in real spectra it lies inside the Lyα forest of lower-redshift gas, so blending and contamination must be modelled; its smaller f λ makes it useful mainly where Lyα saturates; and it needs wavelength coverage down to 1026 Å in the absorber's frame. Like every H I line, Lyβ cannot separate heat from smooth motion (The "same shadow" configurations).
Simplifications
the Si IV column is chosen for visibility
Lyβ is drawn as a clean line, without the lower-redshift Lyα forest it lies in
Sources
Rauch et al. 1996 — §3, Eq. (1) (arXiv p. 5): b_i² = b_b² + 2kT/m_i — two ions of different mass (C IV and Si IV) separate temperature from bulk motion
NIST ASD (Kramida et al. 2024) — Si IV lines output 1390–1405 Å (ASD v5.12, accessed 2026-10-04): 3s–3p ²P°3/2 at 1393.76 Å (vacuum), f = 0.513
NIST ASD (Kramida et al. 2024) — H I lines output 1025–1026 Å (ASD v5.12, accessed 2026-10-03): 1s–3p 2P°3/2 and 2P°1/2, Ritz λ 1025.72182505 / 1025.72296565 Å, f 0.052761 / 0.026380 (multiplet value 0.079142; f-weighted centre 1025.72221 Å)
Wiese & Fuhr 2009 — Table 6, entry No. 2 (1s–3p), p. 578: f = 0.079142
Rauch 1998 — §2.4 (arXiv v1 p. 7): 'For stronger Lyα lines the higher order Lyman lines can provide additional constraints when fitted simultaneously'; §3.4 (p. 15): for saturated, blended lines 'simultaneous fits to the higher order Lyman lines (which are less saturated owing to their lower oscillator strengths) can only help to some degree'
Griesmann & Kling 2000 — Table 1 'Resonance transitions' (arXiv p. 5): Si IV 3s ²S1/2 – 3p ²P3/2 at 1393.76018 Å (vacuum, laboratory Fourier-transform spectroscopy)(secondary)
Each check is an automated test in the source repository; its measured values are written by the test, never typed. Status: verified — claim, equation and source location checked; the numerical check passes where one applies; verified-numeric — the numerical check passes; it is an identity or a calculation, so no literature location applies.
The app's Voigt–Hjerting function H(a, x) is compared with SciPy's Faddeeva function (wofz) on a 42 × 321 grid covering a from 10⁻⁶ to 1 and x from 0 to 2×10⁴.
max relative error 8.5e-05 over 13482 grid points (physical a (1e-4 ≤ a ≤ 6.1e-3), core |x| ≤ 4: 5.0e-05; physical a, wings 4 < |x| ≤ 2e4: 4.9e-05; tiny a (1e-6 ≤ a < 1e-4), all x: 5.0e-05; teaching exaggeration (6.1e-3 < a ≤ 1), all x: 8.5e-05)
Criterion: max relative error < 1e-3 against wofz in every tested region of (a, x)
The map u = aH(z)·x + v_pec is checked against an independent formula, for independence of h, for its sign convention, and against the exact distance–redshift relation (how well a constant H holds across the box); the velocity-to-wavelength convention is checked too.
5/5 checks; linear map vs exact over 20 Mpc/h: 2.0e-03
Criterion: exact identities to 1e-12; linear map within 1 % over ≤ 20 Mpc/h
Optical depth of single, blended, saturated, damped, under-resolved and periodic lines is compared with an independent exact-Voigt calculation; the line-centre optical depth τ₀ and the integral ∫τ du are checked against their analytic values.
8/8 checks; point sampling max rel τ 5.0e-05, max |ΔF| 1.8e-05
Criterion: point sampling: rel τ < 1e-3 and |ΔF| < 1e-4
¶ Toy sightline: line kernels, neutral fraction, and the mean-flux calibration pass
On the toy sightline, the spectrum as drawn (Gaussian thermal kernels) is compared with the same sum using full Voigt kernels. The neutral fraction at mean density is computed from the measured ultraviolet background and checked to be about one in 10⁵. The mean transmitted flux is a calibration, not a test: the toy's ultraviolet background is tuned so that its mean flux matches the observed value, which is why the toy's ultraviolet is about three times the measured one.
3/3 checks; Gauss vs Voigt max |ΔF| 5.7e-04; neutral fraction x_HI = 6.16e-06 at mean density and T0 = 13000 K, for the measured Γ_HI = 0.8e-12 s⁻¹; ⟨F⟩ = 0.680, a calibration (Γ_HI tuned to 3.0× measured to meet it)
Criterion: max |ΔF| < 1e-3 between Gaussian and Voigt kernels; x_HI at mean density between 5×10⁻⁶ and 8×10⁻⁶; ⟨F⟩ = 0.68 ± 0.01 (a calibration: the toy's Γ_HI is tuned to meet it)
Each calculation this resource implements has a short, independent Python reproduction (NumPy and SciPy only) of the equation it publishes. Every reproduction is run on its own: its printed output must equal the expected output formatted from the app's own values at the same inputs, and every full-precision value must agree with the app's within the tolerance it states.
10/10 reproductions: each snippet, run alone, prints its expected output (formatted from the app's own values); 74/74 full-precision values agree with the app within their stated tolerances (largest relative deviation 1e-04)
Criterion: every snippet's output equals its expected output; every script exits 0; every value within its stated tolerance
The third “same shadow” configuration is checked to follow from its matching condition (the same total Doppler parameter and column as the first); the Δχ² separating the configurations is recomputed independently with exact Voigt profiles, in Lyα and Lyβ, at signal-to-noise 20, 50, 100 and 200.
13/13 checks; (iii) constructed by the matching rule (dv = 5.9883 km/s, b_tot = 27.4306 km/s for both); Lyα Δχ²(ii) at S/N 20 / 50 / 100 / 200: 2.79 / 17.4 / 70 / 279, Lyβ 103 at S/N 20; Δχ²(iii) ≤ 9e-06 exact-Voigt, ≤ 2e-03 as the app draws it (S/N ≤ 200, Lyα and Lyβ); max |ΔF| vs an independent calculation 3.6e-05
Criterion: (iii) equals the matching rule to 1e-9 km/s and b_tot(iii) = b_tot(i) to 1e-12; Δχ²(ii) browser vs oracle within 5 % at every S/N in Lyα and Lyβ; Δχ²(iii) < 1e-4 on the exact-Voigt path; the app's max |F(iii) − F(i)| < 1e-4 (its Voigt approximation)
Comoving distance χ(z) and its inverse are compared with SciPy integration (Planck 2018 parameters); the frame conversions λ_obs = λ_rest(1 + z) and their reverses, the quasar frame, and the correspondence of log-wavelength and velocity are checked.
6/6 checks; χ(z) vs scipy max rel 8e-11; z(χ) inverse 2e-09; frame conversions exact to 1e-16; Lyα at z = 3 → 4862.68 Å
¶ Toy spectrum against the code that generated it pass
The app's toy spectrum, as drawn, is compared with the exact-Voigt optical depth computed by the independent code that produced the toy data.
the app (Gaussian kernels, as drawn) vs the generating code's exact-Voigt τ: max |ΔF| 5.8e-04; the app with Voigt kernels 3.2e-05; ⟨F⟩ 0.6804 vs 0.6804
Criterion: as drawn max |ΔF| < 1e-3; Voigt kernels < 1e-4; ⟨F⟩ to 1e-3
In “Who ate this colour?”, the per-cell contributions traced back from a colour are checked to sum to the spectrum's own optical depth at every pixel; the share the display leaves untraced is recorded.
Σ of traced contributions vs the spectrum's τ over 2048 pixels: max relative difference 5e-15; τ in cells under the 1% display cut (kept in the ledger's “rest”): median 3.2%, max 24.8% of a pixel's τ
Criterion: max relative |Σ τ_j − τ| < 1e-12 at every pixel; thresholds read from the beat manifest
In “Where does one parcel land?”, place, own motion, temperature and neutral amount are each changed alone across their ranges and checked to change only what they should — where the optical depth lands, how wide it is, or how deep — and the comparison the scene prints is checked against a fine-grid calculation. The width and peak it prints are checked to be those of the line it draws: heat alone and the parcel's size (the Hubble flow across it) add in quadrature.
91/91 checks over 18 one-variable changes (cells and τ, fine grid and the scene's grid): place and motion move the centre by aH·Δx and Δv (max error 4e-07 km/s) at fixed FWHM and area; temperature changes FWHM² by 4 ln2 Δ(b²) (max rel. error 1e-04) at fixed centre and area; amount scales τ exactly (max rel. 2e-15) at fixed centre and FWHM. The width the scene prints is the drawn line's: heat alone b = 21.4 km/s and the parcel's size 23.8 km/s add in quadrature to the drawn 32.0 km/s (max rel. error 2e-04; peak 3e-04). Flux is not linear (recorded): N_HI ×2 gives τ area ×2 but equivalent width ×1.19 (τ₀ = 4.4); T ×2 at fixed N_HI changes it ×1.13
Criterion: cells: only the intended property changes (relative 1e-12); τ on the fine grid: centre within 1e-3 km/s of aH·Δx or Δv, FWHM and area within 1e-3 (relative) when they should not change, FWHM² increment within 1e-2 of 4 ln2 Δ(b²), τ ratio = N ratio within 1e-9 where τ > 10⁻⁴ of its peak; the drawn line's b within 2e-3 of √(b_thermal² + b_size²) and its peak within 2e-3 of N σ_v/(√π b); the scene's printed comparison, width and peak correct to their printed precision
In “Light multiplies”, each pixel of the drawn light table is read back from the canvas: its transmitted fraction equals e^−τ at that colour, within 8-bit quantisation.
Beat 10 ribbon read back from the canvas: max |T − e^(−τ)| 2.1e-03 over 1000 columns (8-bit bound 2.1e-03); mean 9.0e-04
Criterion: max |T − e^(−τ)| ≤ the 8-bit quantisation bound 0.5/(light − ink)
¶ Closure against fake_spectra pass (declared domain)
The app's real-space method is compared with fake_spectra 2.2.4 on controlled skewers; differences are decomposed by source, and the domain in which the app's method applies is declared.
in domain (one velocity per sampling element, b ≥ 2Δu; 25 cases): max |ΔF| 4.9e-03, |Δ⟨F⟩| 2.9e-04, |ΔEW/EW| 2.8e-03, forest P1D 0.17 % (k ≤ 0.05 s/km), 0.32 % (k ≤ 0.1); C8b compared with the exact reference (fake_spectra's own quadrature error exceeds half the tolerance). Outside the domain, stated not tolerated (max |ΔF| per test case, by case name): several velocities inside one real-space sampling element (unresolved velocity structure; C10 is a deliberately exaggerated 1D shell-crossing stress test): C10s0 0.14, C10s1 0.059, C10s2 0.13, C10s3 0.2, C3c 0.089, C4c 0.021; lines narrower than two velocity pixels (b < 2Δu); fixable by oversampling the skewer: C6a 0.45, C6b 0.21, C6c 0.049, C6d 0.46. The app's code vs its independent Python port: 6e-15 in τ
Criterion: declared domain (one velocity per sampling element, b ≥ 2Δu): max |ΔF| ≤ 0.01, |Δ⟨F⟩| ≤ 1e-3, |ΔEW/EW| ≤ 1 %, P1D within 1 % (k ≤ 0.05 s/km) and 2 % (k ≤ 0.1); where fake_spectra's own error exceeds half the tolerance, compare with the exact reference; out-of-domain cases listed with measured values
Every scientific number, curve and profile the app shows has a computation record: who computed it, from which relation, with which inputs, constants, conventions and assumptions, by which numerical method, where the code is, how it was checked independently, and — for the calculations this resource implements — a short Python reproduction you can run yourself. Four kinds of provenance are kept apart: the literature (why a physical relation is accepted: the claims above), the computation (how this app produced the displayed result: here), the validation (how the implementation was tested: section 7), and the teaching choices (what was simplified or exaggerated: the simplifications and the teaching visuals below). A record is computed by lya_app (code written for this resource), by an external package, by both (lya_app's result, compared on the same physical inputs with an independent external forward-model code — a lya_app computation whose functions are checked against a mathematical library such as SciPy stays lya_app), or is a teaching visualization that computes nothing. The same records, as a table: science/COMPUTATION_MATRIX.md in the source repository.
An independent reproduction in Python, with NumPy and SciPy only. Copy it into your own Python, a notebook or Colab; nothing runs on this page.
import numpy as np
k_B = 1.380649e-16 # Boltzmann constant, erg/K (exact, SI 2019)
m_H = 1.6735328e-24 # mass of the 1H atom, g (1.00782503223 u, NIST)
def b_thermal(T):
"""Doppler parameter of the thermal (Gaussian) profile, km/s: b = sqrt(2 k_B T / m_H)"""
return np.sqrt(2 * k_B * T / m_H) / 1e5
for T in (4.0e3, 1.0e4, 4.2e4):
print(f"T = {T:7.0f} K -> b = {b_thermal(T):.3f} km/s")
Expected output — the app's own values, checked automatically (Python reproductions):
T = 4000 K -> b = 8.124 km/s
T = 10000 K -> b = 12.845 km/s
T = 42000 K -> b = 26.325 km/s
natural broadening from the 2p lifetime only — collisional broadening is negligible at IGM densities
Approximations
w(z) is evaluated with Humlíček's (1982) W4 rational approximation, accurate to < 10⁻⁴ (relative) where tested
Numerical method
Humlíček W4: four regions chosen by |x| + a and by a against 0.195|x| − 0.176, each with its own rational or asymptotic form (region borders as in Zaghloul 2015, Table 1)
An independent reproduction in Python, with NumPy and SciPy only. Copy it into your own Python, a notebook or Colab; nothing runs on this page.
import numpy as np
from scipy.integrate import quad
c = 299792.458 # km/s
lambda0 = 1215.6701 # Lya rest wavelength, Angstrom (vacuum)
Om, OL = 0.3153, 0.6847 # Planck 2018, flat; distances in Mpc/h do not depend on h
z = 3.0
E = lambda zz: np.sqrt(Om * (1 + zz)**3 + OL)
lambda_obs = lambda0 * (1 + z) # every wavelength is stretched by 1 + z
chi = c / 100 * quad(lambda zz: 1 / E(zz), 0, z)[0] # comoving distance, Mpc/h
aH_h = 100 * E(z) / (1 + z) # Hubble velocity per comoving Mpc/h at z, km/s
u = aH_h * 4.4 # across a 4.4 Mpc/h stretch of gas, km/s
dlambda = lambda_obs * (np.exp(u / c) - 1) # its width in observed wavelength (log-lambda convention)
print(f"Lya absorbed at z = 3 is observed at {lambda_obs:.2f} Angstrom")
print(f"comoving distance to z = 3: {chi:.1f} Mpc/h")
print(f"aH/h at z = 3: {aH_h:.2f} km/s per Mpc/h")
print(f"4.4 Mpc/h of Hubble flow: {u:.1f} km/s = {dlambda:.2f} Angstrom observed")
Expected output — the app's own values, checked automatically (Python reproductions):
Lya absorbed at z = 3 is observed at 4862.68 Angstrom
comoving distance to z = 3: 4383.8 Mpc/h
aH/h at z = 3: 114.19 km/s per Mpc/h
4.4 Mpc/h of Hubble flow: 502.4 km/s = 8.16 Angstrom observed
An independent reproduction in Python, with NumPy and SciPy only. Copy it into your own Python, a notebook or Colab; nothing runs on this page.
import numpy as np
Om, OL, z = 0.3153, 0.6847, 3.0
aH_h = 100 * np.sqrt(Om * (1 + z)**3 + OL) / (1 + z) # km/s per comoving Mpc/h (independent of h)
x = np.array([1.9303074084405747, 3.175154331970289, 3.8278357071063325]) # parcels a, b, c: comoving Mpc/h
v_pec = np.array([62.0, 87.0, 20.0]) # their own motion, km/s
u_hubble = aH_h * x # expansion alone: the order in velocity follows the order in space
u = aH_h * x + v_pec # with their own motion: b and c land close together
for name, uh, uu in zip("abc", u_hubble, u):
print(f"parcel {name}: expansion alone {uh:.1f} km/s, with its own motion {uu:.1f} km/s")
Expected output — the app's own values, checked automatically (Python reproductions):
parcel a: expansion alone 220.4 km/s, with its own motion 282.4 km/s
parcel b: expansion alone 362.6 km/s, with its own motion 449.6 km/s
parcel c: expansion alone 437.1 km/s, with its own motion 457.1 km/s
a Gaussian clump of width σ_x with a linear velocity field has a Gaussian spread of velocities that adds in quadrature to the thermal width; τ₀ is the line-centre optical depth of the Doppler profile
what the scene draws — the clump cut into sub-cells s_k ∈ [−5, 5], each a thermal line at its own velocity (its sum reproduces the quadrature law above)
Variables
σ_x — the clump's 1σ width along the beam [comoving Mpc/h]
dv/dx — the parcel's own peculiar-velocity gradient (Beat 6: set to zero) [km/s per Mpc/h]
N_HI — neutral hydrogen column [cm⁻²]
σ_v — integrated cross-section per unit velocity, (πe²/m_e c) f λ₀ [cm² km/s]
x comoving from the near edge of the stretch; the parcel's column N is the total over the clump
Assumptions
Gaussian density profile along the beam; one temperature; a linear velocity field
Beat 6 sets the parcel's own peculiar-velocity gradient to zero, isolating heat and the Hubble flow across its size; with its own gradient the same parcel has b = 23.3 km/s
Approximations
sub-cell sampling: 161 sub-cells over ±5σ in Beats 6–10, 201 in Beat 13; the damping wings are kept (Voigt kernels)
Numerical method
parcelToCells (sub-cells with Gaussian weights) → tauFromCells (direct summation); parcelWidths gives the analytic widths
Validated for
Beat 6's control ranges (place, motion, temperature, amount)
An independent reproduction in Python, with NumPy and SciPy only. Copy it into your own Python, a notebook or Colab; nothing runs on this page.
import numpy as np
k_B, m_H = 1.380649e-16, 1.6735328e-24 # erg/K; g (the 1H atom)
sigma_v = 2.6540e-2 * 0.4164 * 1215.6701e-8 / 1e5 # (pi e^2 / m_e c) f lambda0: cm^2 km/s
Om, OL, z = 0.3153, 0.6847, 3.0
aH_h = 100 * np.sqrt(Om * (1 + z)**3 + OL) / (1 + z) # km/s per comoving Mpc/h
# the toy's parcel b, as Beat 6 starts it
T, sigma_x, N = 27744.0, 0.14721893675599934, 10**14.27 # K; comoving Mpc/h (1 sigma); cm^-2
def widths(dvdx):
"""thermal b, size b (Gaussian clump with a linear velocity field), and their quadrature sum, km/s"""
b_th = np.sqrt(2 * k_B * T / m_H) / 1e5
b_size = np.sqrt(2) * abs(aH_h + dvdx) * sigma_x
return b_th, b_size, np.hypot(b_th, b_size)
b_th, b_size, b = widths(0.0) # this scene: the parcel's own gradient switched off
tau0 = N * sigma_v / (np.sqrt(np.pi) * b) # peak optical depth of the Doppler profile
print(f"thermal b = {b_th:.2f} km/s, size b = {b_size:.2f} km/s, together b = {b:.2f} km/s")
print(f"peak optical depth tau0 = {tau0:.2f}")
_, b_size_own, b_own = widths(-70.41475616076316) # the same parcel with its own (infalling) gradient
print(f"with its own gradient: size b = {b_size_own:.2f} km/s, b = {b_own:.2f} km/s")
Expected output — the app's own values, checked automatically (Python reproductions):
thermal b = 21.40 km/s, size b = 23.77 km/s, together b = 31.98 km/s
peak optical depth tau0 = 4.41
with its own gradient: size b = 9.11 km/s, b = 23.26 km/s
each contribution's column times the cross-section per unit velocity times its normalised profile; along a sightline each real-space cell contributes its own column
profile sampled at pixel centres (point sampling); each contribution's kernel extends to 8 b and its damping wings until τ < 10⁻⁶; periodic boxes use the nearest image
Assumptions
optical depths of independent absorbers add (a linear
optically passive medium); one velocity per sampling element
Approximations
the toy forest is drawn with Gaussian kernels: max |ΔF| 6e-04 against Voigt kernels
Numerical method
direct summation over contributions on a uniform velocity grid
An independent reproduction in Python, with NumPy and SciPy only. Copy it into your own Python, a notebook or Colab; nothing runs on this page.
import numpy as np
sigma_v = 2.6540e-2 * 0.4164 * 1215.6701e-8 / 1e5 # (pi e^2 / m_e c) f lambda0: cm^2 km/s
for t in (1.73, 11.9, 0.37):
print(f"tau = {t:5.2f} -> F = exp(-tau) = {np.exp(-t):.3e}")
N, b = 1e15, 25.0 # cm^-2, km/s
tau0 = N * sigma_v / (np.sqrt(np.pi) * b) # line-centre optical depth of a Doppler line
print(f"N = 1e15 cm^-2, b = 25 km/s: tau0 = {tau0:.2f}, F at line centre = {np.exp(-tau0):.1e} (saturated)")
Expected output — the app's own values, checked automatically (Python reproductions):
tau = 1.73 -> F = exp(-tau) = 1.773e-01
tau = 11.90 -> F = exp(-tau) = 6.790e-06
tau = 0.37 -> F = exp(-tau) = 6.907e-01
N = 1e15 cm^-2, b = 25 km/s: tau0 = 30.32, F at line centre = 6.8e-14 (saturated)
the spectrograph convolves the transmitted flux (never τ) with its line-spread function; the detector adds white noise per pixel
Variables
FWHM — the resolution [km/s]
S/N — signal-to-noise per pixel (continuum-normalised) [dimensionless]
Conventions
order: intrinsic τ → intrinsic F → line-spread function → pixels → noise; the gas and τ are untouched
pixels max(2.5 km/s, FWHM/3) wide
Assumptions
a Gaussian line-spread function; white
flux-independent noise; no sky
no continuum error
Approximations
the kernel is cut at 5σ
Numerical method
discrete convolution with a normalised Gaussian kernel (edges extended; the toy forest is periodic); noise from a seeded generator (Box–Muller), so a given seed draws the same noise
Validated for
FWHM 0–200 km/s and S/N 5, 20, 50 or no noise (Beat 11); S/N 20 (Beat 13)
Validation
Python reproductions an independent Python convolution (the noise draw is not reproduced: it depends on the generator) tolerance: relative 1e-4
An independent reproduction in Python, with NumPy and SciPy only. Copy it into your own Python, a notebook or Colab; nothing runs on this page.
import numpy as np
# Planck 2018 (TT,TE,EE+lowE+lensing): baryon density and helium mass fraction; the redshift of the toy
Obh2, Yp, z = 0.02237, 0.2454, 3.0
G = 6.67430e-8 # cm^3 g^-1 s^-2 (CODATA 2022)
m_H = 1.6735328e-24 # g, the 1H atom
Mpc = 3.0856775814913673e24 # cm
H100 = 100e5 / Mpc # 100 km/s/Mpc in 1/s
rho_crit_h2 = 3 * H100**2 / (8 * np.pi * G) # critical density today per h^2, g/cm^3
nH_bar = Obh2 * rho_crit_h2 * (1 - Yp) / m_H * (1 + z)**3 # mean hydrogen density at z, cm^-3 (proper)
ne_over_nH = 1 + 2 * Yp / (4 * (1 - Yp)) # electrons per H nucleus, H and He fully ionised
def x_HI(T, Gamma12, Delta=1.0):
"""neutral fraction in photoionisation equilibrium (case A): alpha_A(T) n_e / Gamma_HI"""
alpha_A = 4.2e-13 * (T / 1e4) ** -0.7 # cm^3/s
return alpha_A * ne_over_nH * nH_bar * Delta / (Gamma12 * 1e-12)
print(f"mean hydrogen density at z = 3: {nH_bar:.4e} cm^-3")
for T in (1.0e4, 1.3e4):
print(f"x_HI at mean density, T = {T:.1e} K, Gamma_HI = 0.8e-12 /s: {x_HI(T, 0.8):.3e}")
Expected output — the app's own values, checked automatically (Python reproductions):
mean hydrogen density at z = 3: 1.2126e-05 cm^-3
x_HI at mean density, T = 1.0e+04 K, Gamma_HI = 0.8e-12 /s: 7.401e-06
x_HI at mean density, T = 1.3e+04 K, Gamma_HI = 0.8e-12 /s: 6.159e-06
An independent reproduction in Python, with NumPy and SciPy only. Copy it into your own Python, a notebook or Colab; nothing runs on this page. Note: this reproduces the sum the trace relies on — contributions add, and the transmission is their product — on a small example with Voigt kernels; the trace itself uses the Gaussian kernels above, and Tracing a colour back to the gas checks that its contributions sum to τ.
import numpy as np
from scipy.special import wofz
sigma_v = 2.6540e-2 * 0.4164 * 1215.6701e-8 / 1e5 # (pi e^2 / m_e c) f lambda0: cm^2 km/s
gamma = 6.2649e8 * 1215.6701e-8 / (4 * np.pi) / 1e5 # natural half width, km/s
def tau(u, u0, N, b):
"""one contribution: N sigma_v phi(u - u0), phi = H(a, x) / (sqrt(pi) b) with x = (u - u0)/b, a = gamma/b"""
return N * sigma_v * wofz((u - u0) / b + 1j * gamma / b).real / (np.sqrt(np.pi) * b)
cells = [(-20.0, 5e13, 15.0), (0.0, 1e14, 20.0), (25.0, 3e13, 12.0)] # (u0 km/s, N cm^-2, b km/s)
for u in (0.0, 20.0):
each = [tau(u, *c) for c in cells]
total = sum(each)
print(f"u = {u:4.0f} km/s: tau = " + " + ".join(f"{t:.4f}" for t in each) + f" = {total:.4f}; F = exp(-tau) = {np.exp(-total):.5f} = product of exp(-tau_j) = {np.prod(np.exp(-np.array(each))):.5f}")
Expected output — the app's own values, checked automatically (Python reproductions):
u = 0 km/s: tau = 0.4273 + 3.7886 + 0.0249 = 4.2408; F = exp(-tau) = 0.01440 = product of exp(-tau_j) = 0.01440
u = 20 km/s: tau = 0.0022 + 1.3943 + 1.5922 = 2.9887; F = exp(-tau) = 0.05036 = product of exp(-tau_j) = 0.05036
clumps cut into 201 sub-cells over ±5σ; Lyβ keeps Lyα's damping constant (its wings are immaterial at these columns)
Numerical method
construction by the matching rule (science/tools/same_shadow_configs.py); lines by parcelToCells → tauFromCells; Δχ² by direct summation
Validated for
S/N 20 to 200
Validation
Same-shadow degeneracy · Python reproductions an exact-Voigt sub-cell oracle (Same-shadow degeneracy) and a Python sub-cell sum with SciPy's Faddeeva function — two implementations, neither the app's tolerance: Δχ²(ii) within 5 % (oracle) and 1e-3 (Python); Δχ²(iii) < 1e-4 on the exact-Voigt paths (the oracle and Python), < 1e-2 as the app draws it with Humlíček's W4
An independent reproduction in Python, with NumPy and SciPy only. Copy it into your own Python, a notebook or Colab; nothing runs on this page. Note: the matching condition is closed form; the lines are summed over each clump's sub-cells, with no closed form, so that (iii) matches (i) is computed, not assumed.
import numpy as np
from scipy.special import wofz
k_B, m_H = 1.380649e-16, 1.6735328e-24
sigma_v = 2.6540e-2 * 0.4164 * 1215.6701e-8 / 1e5 # Lya: (pi e^2 / m_e c) f lambda0, cm^2 km/s
gamma = 6.2649e8 * 1215.6701e-8 / (4 * np.pi) / 1e5 # km/s
lyb = 0.079142 * 1025.72221 / (0.4164 * 1215.6701) # Lyb optical depth relative to Lya: (f lambda) ratio
b_T = lambda T: np.sqrt(2 * k_B * T / m_H) / 1e5 # thermal Doppler parameter, km/s
# gases in velocity space: centre u, Hubble-flow width su (1 sigma), bulk velocity v, extra smooth velocity dv across 1 sigma
I = [dict(u=0.0, su=2.23, N=2.5e14, T=45000.0, v=0.0, dv=0.0)] # one hot, compact clump
II = [dict(u=-39.0, su=12.8, N=9.75e13, T=1e4, v=24.9, dv=0.0), # two cold clumps, converging
dict(u=39.0, su=12.8, N=9.75e13, T=1e4, v=-24.9, dv=0.0)]
dv = np.sqrt(2.23**2 + (b_T(45000.0)**2 - b_T(1e4)**2) / 2) - 11.15 # matching: b_T^2 + 2 sigma_u^2 equal
III = [dict(u=0.0, su=11.15, N=2.5e14, T=1e4, v=0.0, dv=dv)] # one cold clump, expanding
u = np.linspace(-150, 150, 121) # 2.5 km/s pixels
s = np.linspace(-5, 5, 401); w = np.exp(-s**2 / 2); w /= w.sum() # a clump in sub-cells over +-5 sigma
def tau(gas, scale=1.0): # no closed form here: every sub-cell is a thermal line at its own velocity
t = np.zeros_like(u)
for c in gas:
b = b_T(c['T']); x = (u[:, None] - c['u'] - c['v'] - (c['su'] + c['dv']) * s) / b
t += scale * c['N'] * sigma_v * (w * wofz(x + 1j * gamma / b).real).sum(1) / (np.sqrt(np.pi) * b)
return t
sig3 = lambda x: f"{float(f'{x:.3g}'):g}"
print(f"(iii)'s smooth extra velocity from the matching condition: dv = {dv:.4f} km/s")
for line, scale in (("Lya", 1.0), ("Lyb", lyb)):
F = {k: np.exp(-tau(g, scale)) for k, g in (("i", I), ("ii", II), ("iii", III))}
for snr in (20, 50, 100, 200):
chi2 = {k: np.sum((F[k] - F["i"]) ** 2) * snr**2 for k in ("ii", "iii")} # noise-free, against (i)
print(f"{line} S/N {snr:3d}: dchi2(ii) = {sig3(chi2['ii']):>6} dchi2(iii) = {chi2['iii']:.2f}")
Expected output — the app's own values, checked automatically (Python reproductions):
An independent reproduction in Python, with NumPy and SciPy only. Copy it into your own Python, a notebook or Colab; nothing runs on this page. Note: the same accumulation on a small example; the full sightline is checked against the two independent codes.
import numpy as np
from scipy.special import wofz
sigma_v = 2.6540e-2 * 0.4164 * 1215.6701e-8 / 1e5 # (pi e^2 / m_e c) f lambda0: cm^2 km/s
gamma = 6.2649e8 * 1215.6701e-8 / (4 * np.pi) / 1e5 # natural half width, km/s
def tau(u, u0, N, b):
"""one contribution: N sigma_v phi(u - u0), phi = H(a, x) / (sqrt(pi) b) with x = (u - u0)/b, a = gamma/b"""
return N * sigma_v * wofz((u - u0) / b + 1j * gamma / b).real / (np.sqrt(np.pi) * b)
cells = [(-20.0, 5e13, 15.0), (0.0, 1e14, 20.0), (25.0, 3e13, 12.0)] # (u0 km/s, N cm^-2, b km/s)
for u in (0.0, 20.0):
each = [tau(u, *c) for c in cells]
total = sum(each)
print(f"u = {u:4.0f} km/s: tau = " + " + ".join(f"{t:.4f}" for t in each) + f" = {total:.4f}; F = exp(-tau) = {np.exp(-total):.5f} = product of exp(-tau_j) = {np.prod(np.exp(-np.array(each))):.5f}")
Expected output — the app's own values, checked automatically (Python reproductions):
u = 0 km/s: tau = 0.4273 + 3.7886 + 0.0249 = 4.2408; F = exp(-tau) = 0.01440 = product of exp(-tau_j) = 0.01440
u = 20 km/s: tau = 0.0022 + 1.3943 + 1.5922 = 2.9887; F = exp(-tau) = 0.05036 = product of exp(-tau_j) = 0.05036
fake_spectra.spectra.Spectra on mock Gadget-format snapshots (science/oracle/fakespectra_mocksnap.py writes them; science/oracle/fakespectra_run_closure.py runs Spectra and reads τ)
Inputs
particle positions, velocities, temperatures, smoothing lengths and neutral hydrogen of synthetic test cases (science/validation/closure/cases)
Conversions by this resource
this project writes its synthetic skewers as mock snapshots in Gadget units (comoving kpc/h; velocities with the √a factor) and compares τ and F on the same pixel grid
skewer velocities are never taken from get_velocity
Validation
Closure against fake_spectra compared with this project's real-space method and an exact reference on the same cases tolerance: declared-domain criteria (max |ΔF| ≤ 0.01; |Δ⟨F⟩| ≤ 1e-3)
three Gaussian parcels summarising peak regions of the continuous sightline, and Beat 13's idealised clumps — a teaching decomposition, not separate clouds
Becker & Bolton 2013. Becker, G. D. & Bolton, J. S. 2013, MNRAS 436, 1023–1039, 'New measurements of the ionizing ultraviolet background over 2 < z < 5 and implications for hydrogen reionization' DOI · arXiv
Becker et al. 2011. Becker, G. D., Bolton, J. S., Haehnelt, M. G. & Sargent, W. L. W. 2011, MNRAS 410, 1096–1112, 'Detection of extended He II reionization in the temperature evolution of the intergalactic medium' DOI · ADS · arXiv
Becker et al. 2013. Becker, G. D., Hewett, P. C., Worseck, G. & Prochaska, J. X. 2013, MNRAS 430, 2067–2081, 'A refined measurement of the mean transmitted flux in the Lyα forest over 2 < z < 5 using composite quasar spectra' DOI · arXiv
Bird et al. 2015. Bird, S., Haehnelt, M., Neeleman, M., Genel, S., Vogelsberger, M. & Hernquist, L. 2015, MNRAS 447, 1834–1846, 'Reproducing the kinematics of damped Lyman α systems' DOI · ADS · arXiv
CODATA 2022 (Mohr et al. 2025). Mohr, P. J., Newell, D. B., Taylor, B. N. & Tiesinga, E. 2025, Rev. Mod. Phys. 97, 025002, 'CODATA recommended values of the fundamental physical constants: 2022' (also J. Phys. Chem. Ref. Data 54, 033105); values as served by NIST CUU DOI · link
Dijkstra 2014. Dijkstra, M. 2014, PASA 31, e040, 'Lyα Emitting Galaxies as a Probe of Reionisation' DOI · ADS · arXiv
Draine 2011. Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton University Press)
fake_spectra (Bird). fake_spectra (S. Bird), MIT license, https://github.com/sbird/fake_spectra; PyPI package 'fake_spectra' v2.2.4 link · PyPI
Griesmann & Kling 2000. Griesmann, U. & Kling, R. 2000, ApJ 536, L113, 'Interferometric Measurement of Resonance Transition Wavelengths in C IV, Si IV, Al III, Al II, and Si II' DOI · arXiv
Hogg 1999. Hogg, D. W. 1999, 'Distance measures in cosmology', arXiv:astro-ph/9905116 (v4, 16 Dec 2000); pedagogical note, not published in a journal DOI · arXiv
Hui & Gnedin 1997. Hui, L. & Gnedin, N. Y. 1997, MNRAS 292, 27–42, 'Equation of state of the photoionized intergalactic medium' DOI · ADS · arXiv
Hui, Gnedin & Zhang 1997. Hui, L., Gnedin, N. Y. & Zhang, Y. 1997, ApJ 486, 599–622, 'The Statistics of Density Peaks and the Column Density Distribution of the Lyα Forest' DOI · ADS · arXiv
Humlíček 1982. Humlíček, J. 1982, J. Quant. Spectrosc. Radiat. Transfer 27, 437–444, 'Optimized computation of the Voigt and complex probability functions' DOI · ADS
Karaçaylı et al. 2024. Karaçaylı, N. G. et al. 2024, MNRAS 528, 3941, 'Optimal 1D Lyα forest power spectrum estimation — III. DESI early data' DOI · ADS · arXiv
McQuinn 2016. McQuinn, M. 2016, ARA&A 54, 313–362, 'The Evolution of the Intergalactic Medium' DOI · ADS · arXiv
Meiksin 2009. Meiksin, A. A. 2009, Rev. Mod. Phys. 81, 1405–1469, 'The physics of the intergalactic medium' DOI · ADS · arXiv
NIST ASD (Kramida et al. 2024). Kramida, A., Ralchenko, Yu., Reader, J., and NIST ASD Team (2024). NIST Atomic Spectra Database (ver. 5.12), https://physics.nist.gov/asd [accessed 2026-10-02]. National Institute of Standards and Technology, Gaithersburg, MD. link
NIST isotopic compositions. NIST Atomic Weights and Isotopic Compositions with Relative Atomic Masses, https://physics.nist.gov/comp — relative atomic masses from 'The AME2012 Atomic Mass Evaluation' (Wang, M. et al. 2012, Chinese Physics C 36(12); DOI 10.1088/1674-1137/36/12/003 as given on the NIST page) link
NIST isotopic compositions (silicon). NIST Atomic Weights and Isotopic Compositions with Relative Atomic Masses, https://physics.nist.gov/comp — silicon entry; relative atomic masses from 'The AME2012 Atomic Mass Evaluation' (Wang, M. et al. 2012, Chinese Physics C 36(12); DOI 10.1088/1674-1137/36/12/003 as given on the NIST page) link
Palanque-Delabrouille et al. 2013. Palanque-Delabrouille, N., Yèche, Ch., Borde, A., Le Goff, J.-M., Rossi, G., et al. 2013, A&A 559, A85, 'The one-dimensional Lyα forest power spectrum from BOSS' DOI · ADS · arXiv
Planck Collaboration 2020. Planck Collaboration (Aghanim, N. et al.) 2020, A&A 641, A6, 'Planck 2018 results. VI. Cosmological parameters'; Corrigendum A&A 652, C4 (2021) DOI · ADS · arXiv
Rauch 1998. Rauch, M. 1998, ARA&A 36, 267–316, 'The Lyman Alpha Forest in the Spectra of Quasistellar Objects' DOI · ADS · arXiv
Rauch et al. 1996. Rauch, M., Sargent, W. L. W., Womble, D. S. & Barlow, T. A. 1996, ApJ 467, L5, 'Temperature and Kinematics of C IV Absorption Systems' DOI · arXiv
Rybicki & Lightman 1979. Rybicki, G. B. & Lightman, A. P. 1979, Radiative Processes in Astrophysics (New York: Wiley); Wiley-VCH reprint 2004 ISBN 978-0-471-82759-7
SciPy scipy.special.wofz. SciPy scipy.special.wofz documentation (v1.15.2 docs page, the version used here; the v1.18.0 page is the same): w(z) = exp(−z²) erfc(−iz); Reference [1] Steven G. Johnson, Faddeeva W function implementation, http://ab-initio.mit.edu/Faddeeva link
Springel 2005 (GADGET-2 user guide). Springel, V. 2005, 'User guide for GADGET-2' (dated May 1, 2005), https://wwwmpa.mpa-garching.mpg.de/gadget/users-guide.pdf link
Wiese & Fuhr 2009. Wiese, W. L. & Fuhr, J. R. 2009, J. Phys. Chem. Ref. Data 38, 565–720, 'Accurate Atomic Transition Probabilities for Hydrogen, Helium, and Lithium'; Publisher's Note: J. Phys. Chem. Ref. Data 38, 1129 (2009) DOI · ADS
Wolfe, Gawiser & Prochaska 2005. Wolfe, A. M., Gawiser, E. & Prochaska, J. X. 2005, ARA&A 43, 861–918, 'Damped Lyα Systems' DOI · ADS · arXiv
Zaghloul 2015. Zaghloul, M. R. 2015, arXiv:1505.05596, 'A simple reform for treating the loss of accuracy of Humlicek's W4 algorithm near the real axis' arXiv
Cite and version
Lyα: from gas to forest — science notes, version 0.1.1. How to cite: see CITATION.cff in the source repository. Fonts and KaTeX are redistributed under their own licences: third-party licences.