spectroxide

Numerical solver for cosmic microwave background (CMB) spectral distortions from energy and photon injection in the early Universe.

Energy released into the photon-baryon plasma at redshifts $z \sim 10^3 - 5 \times 10^6$ creates deviations from a perfect blackbody spectrum. These spectral distortions — $\mu$-type (chemical potential) and $y$-type (Compton) — encode information about early-universe physics, from the dissipation of primordial acoustic waves to exotic particle decays and dark sector interactions.

spectroxide solves the coupled photon-electron Boltzmann equation including Compton scattering (Kompaneets equation), double Compton emission, bremsstrahlung, and Hubble expansion. It provides both a full partial differential equation (PDE) solver (Rust) and a fast Green's function approximation (Rust and Python).

spectroxide provides:

• Full PDE solver in Rust: implicit Kompaneets and coupled double Compton (DC) and bremsstrahlung (BR) with adaptive stepping
• Green's function mode for fast approximate calculations (pure Python, no compilation needed)
• 9 built-in injection scenarios: single burst, decaying particles (heat or photon channel), dark matter (DM) annihilation (s-wave or p-wave), dark photon oscillation, monochromatic photon injection, and tabulated sources (plus custom heating through the Rust API)
• Comprehensive test suite: 470+ unit, integration, and doc-tests
• Zero production dependencies in Rust (pure std library)

The install script handles everything — Rust toolchain, compilation, and Python package:

git clone https://github.com/bakerem/spectroxide.git
cd spectroxide

# Into a new conda environment (recommended)
./install.sh --conda spectroxide --extras notebook

# Or into your current Python environment
./install.sh

The --extras flag controls which Python dependencies are installed:

Every install pulls in the base dependencies numpy and scipy; the extras add on top of those:

Run ./install.sh --help to see all options, such as the flags that skip steps or print verbose output.

The solve function runs the Rust solver for one injection scenario and returns the distortion parameters and the spectrum. The run_sweep function runs many injection redshifts in one call:

from spectroxide import solve, run_sweep
import numpy as np

# Single-burst injection at z = 2e5 (mu-y transition era)
result = solve(
    injection={"type": "single_burst", "z_h": 2e5, "sigma_z": 5000},
    delta_rho=1e-5,
    z_start=3e5, z_end=1e3,
)
print(f"mu = {result.mu:.3e}, y = {result.y:.3e}")

# Full spectrum is in result.x, result.delta_n
x, delta_n = result.x, result.delta_n

# Decaying particle (gamma_x is the decay rate Gamma_X in 1/s;
# gamma_x = 1e-10 => lifetime ~1e10 s, decay peaks near z ~ 5e4)
result = solve(
    injection={"type": "decaying_particle", "f_x": 5e5, "gamma_x": 1e-10},
    z_start=5e6, z_end=1e3,
)

# DM annihilation (s-wave)
result = solve(
    injection={"type": "annihilating_dm", "f_ann": 1e-22},
    z_start=5e6, z_end=1e3,
)

# Dark photon oscillation (resonant conversion, narrow-width approximation)
result = solve(
    injection={"type": "dark_photon_resonance", "epsilon": 1e-7, "m_ev": 1e-5},
    z_end=1e3,  # z_start auto-set to z_res
)

# Multi-redshift single-burst sweep (one Rust call, parallelised internally)
sweep = run_sweep(z_injections=np.geomspace(2e3, 5e5, 15).tolist(), delta_rho=1e-5)
r = sweep["results"][0]
print(f"z_h = {r['z_h']:.1e}: mu = {r['pde_mu']:.3e}, y = {r['pde_y']:.3e}")

The Green's function mode is pure Python. It gives $\mu$ and $y$ in milliseconds, at the cost of accuracy in the $\mu$-$y$ transition era:

from spectroxide import run_single, mu_from_heating, y_from_heating
import numpy as np

# Quick estimate via Green's function
result = run_single(z_h=2e5, delta_rho=1e-5)
print(f"mu = {result['mu']:.3e}, y = {result['y']:.3e}")

# mu and y from an arbitrary heating history dQ/dz(z)
dq_dz = lambda z: 1e-5 * np.exp(-((z - 2e5)**2) / (2 * 5000**2)) / (5000 * np.sqrt(2*np.pi))
mu = mu_from_heating(dq_dz, 1e3, 5e6)
y = y_from_heating(dq_dz, 1e3, 5e6)

To call the solver from Rust, import the prelude, which re-exports the main types:

use spectroxide::prelude::*;

// Green's function (fast, approximate)
let g = greens::greens_function(3.0, 2e5);

// Full PDE solver
let cosmo = Cosmology::default();
let mut solver = ThermalizationSolver::new(cosmo, GridConfig::default());
solver.set_injection(InjectionScenario::SingleBurst {
    z_h: 2e5, delta_rho_over_rho: 1e-5, sigma_z: 5000.0,
}).unwrap();
solver.set_config(SolverConfig {
    z_start: 5e5, z_end: 1e3, ..SolverConfig::default()
});
let snapshots = solver.run_with_snapshots(&[1e3]);

The spectroxide binary exposes the same solver through subcommands:

# Show help and available subcommands
cargo run --release --bin spectroxide

# Specific injection scenarios (solve subcommand)
cargo run --release --bin spectroxide -- solve decaying-particle --f-x 5e5 --gamma-x 1e-10
cargo run --release --bin spectroxide -- solve annihilating-dm --f-ann 1e-22
cargo run --release --bin spectroxide -- solve dark-photon-resonance --epsilon 1e-7 --m-ev 1e-5
# Sweep over injection redshifts
cargo run --release --bin spectroxide -- sweep --delta-rho 1e-5

# Green's function mode (no PDE)
cargo run --release --bin spectroxide -- greens --z-h 2e5 --delta-rho 1e-5

The CLI writes JSON to stdout by default. Use --format csv or --format table for other formats, and --output <path> to write to a file.

The tutorial notebooks in notebooks/tutorials/ are numbered in the suggested reading order:

Additional notebooks cover physics topics (adiabatic cooling, photon injection, injection-width resolution) and observational constraints (FIRAS photon injection limits, dark-photon limit conventions). Development and validation notebooks are in the development notebooks directory.

The Python solve(injection={"type": ...}) dict uses the snake_case type key in the middle column (not the Rust CamelCase name in the first column). Pass arbitrary sources as callables instead of an injection dict.

Units: $f_X$ [eV] is the energy released per hydrogen nucleus ($\Gamma_X$ [1/s] the decay rate); $f_{\rm ann}$ [eV/s] is energy per hydrogen nucleus per second (both rates multiply $n_{\rm H}$, not the baryon number density); $\Delta N/N$ is the fractional photon-number injection.

To recover the scalar $(\mu, y, \Delta T/T)$ decomposition from any distortion $\Delta n(x)$, use decompose_distortion(x, delta_n) (a solve result exposes these directly as result.mu and result.y).

The repository is laid out as follows:

src/
├── lib.rs                 # Library root + prelude
├── main.rs                # CLI binary entry
├── cli.rs                 # CLI argument parsing and dispatch
├── output.rs              # JSON serialization of solver results
├── kompaneets.rs          # Compton scattering (implicit-explicit Newton solver)
├── double_compton.rs      # DC emission (photon-number changing)
├── bremsstrahlung.rs      # BR emission (non-relativistic Gaunt factor)
├── solver.rs              # PDE integrator (coupled Kompaneets + DC/BR)
├── greens.rs              # Green's function approximation
├── energy_injection.rs    # Injection scenarios
├── dark_photon.rs         # Narrow-width-approximation helpers for γ↔A' (plasma freq, γ_con)
├── distortion.rs          # mu/y/DeltaT decomposition + intensity conversion
├── cosmology.rs           # Flat LCDM background
├── spectrum.rs            # Planck, Bose-Einstein, spectral shapes
├── grid.rs                # Non-uniform frequency grid
├── electron_temp.rs       # Electron temperature (quasi-stationary)
├── recombination.rs       # Peebles 3-level atom (RECFAST fudge + Lyman-α correction) + Saha; X_H evolved at T_e
├── constants.rs           # CODATA 2018 constants
└── bin/check_adiabatic.rs # Adiabatic cooling validation utility

python/spectroxide/
├── __init__.py            # Public API
├── cosmology.py           # Flat ΛCDM background (mirrors src/cosmology.rs)
├── greens.py              # Pure Python Green's function (NumPy)
├── greens_table.py        # Precomputed Green's function tables
├── solver.py              # Rust binary wrapper + run_single()
├── firas.py               # FIRAS data and constraint utilities
├── dark_photon.py         # Narrow-width-approximation helpers (γ_con, z_res) — Python port
├── cosmotherm.py          # CosmoTherm data loaders (submodule import)
├── plot_params.py         # Plot constants (submodule import)
├── style.py               # Matplotlib style helpers
└── _validation.py         # Input validation (errors + warnings)

tests/
├── common/mod.rs          # Shared test setup and cached PDE runs
├── pde_heat.rs            # PDE: heat injection (bursts, decay, annihilation)
├── pde_photon.rs          # PDE: photon injection
├── solver_numerics.rs     # PDE: stability, stepping, null tests
├── gf_visibility.rs       # Green's functions, visibility, decomposition
├── components.rs          # Cosmology, grid, rates, kernel, table I/O
├── dark_sector.rs         # Photon depletion and resonant conversion
├── heat_delivery.rs       # Post-recombination heat delivery
├── ionization_coupling.rs # Hydrogen ionization evolved with T_e
├── adversarial_inputs.rs  # Edge cases, bad inputs
├── cosmotherm_comparison.rs # PDE vs CosmoTherm reference data
├── greens_function_checks.rs # GF spectral shapes, limits, conservation
├── coverage_gaps.rs       # Additional coverage
├── convergence_order.rs   # Grid/timestep convergence
├── cli_integration.rs     # CLI end-to-end
├── science_suite.rs       # End-to-end physics validation
├── physics_identities.rs  # Closed-form/published identities (audit anchors)
├── mms_convergence.rs     # Manufactured solutions (verifies the scheme)
├── conservation_fuzz.rs   # Randomized energy/number-closure fuzzing
├── kompaneets_moments.rs  # Exact moment hierarchy vs published equation
├── compton_equilibrium_analytic.rs # Perturbative Δρ_eq vs mpmath quadrature
├── rate_coefficients_first_principles.rs # DC/BR magnitudes from literal CODATA
└── mu_photosphere_profile.rs # Fitted x_c(z) vs Chluba (2015) Eq. 25

dev/
├── scripts/               # Validation and diagnostic scripts
└── notebooks/             # Validation and dev notebooks (6 notebooks)

If you use spectroxide in your research, please cite the accompanying paper and Chluba & Sunyaev (2012). If you use the Green's function mode, please also cite Chluba (2013) and Chluba (2015):

• Baker, Liu & Mishra-Sharma (2026), "spectroxide: a code package for computing cosmic microwave background spectral distortions" (arXiv:2604.24838) — this paper
• Chluba & Sunyaev (2012), "The evolution of CMB spectral distortions in the early Universe", MNRAS 419, 1294 (arXiv:1109.6552, doi:10.1111/j.1365-2966.2011.19786.x) — CosmoTherm thermalization solver
• Chluba (2013), "Green's function of the cosmological thermalization problem", MNRAS 434, 352 (arXiv:1304.6120, doi:10.1093/mnras/stt1025) — Green's function I
• Chluba (2015), "Green's function of the cosmological thermalization problem — II. Effect of photon injection and constraints", MNRAS 454, 4182 (arXiv:1506.06582, doi:10.1093/mnras/stv2243) — Green's function II (photon injection)

A machine-readable CITATION.cff is also included in the repository root.

Contributions — new injection scenarios, improved physics, validation, bug fixes — are welcome. Most contributions to spectroxide (including the bulk of the original codebase) are written with LLM assistance, and the workflow is built around that: the human supplies the physics (analytic limits, paper references, dimensional checks), the LLM supplies the implementation, and the human is responsible for verifying that tests have independent targets rather than ones calibrated to the code's own output.

See the contributing guide for full details. If you are using an LLM, also drop the LLM context file into its system prompt — it encodes the numerical pitfalls and review rules that have caught real bugs during development.

MIT