Enum InjectionScenario

#[non_exhaustive]
pub enum InjectionScenario {
    SingleBurst {
        z_h: f64,
        delta_rho_over_rho: f64,
        sigma_z: f64,
    },
    DecayingParticle {
        f_x: f64,
        gamma_x: f64,
    },
    AnnihilatingDM {
        f_ann: f64,
    },
    AnnihilatingDMPWave {
        f_ann: f64,
    },
    MonochromaticPhotonInjection {
        x_inj: f64,
        delta_n_over_n: f64,
        z_h: f64,
        sigma_z: f64,
        sigma_x: f64,
    },
    DecayingParticlePhoton {
        x_inj_0: f64,
        f_inj: f64,
        gamma_x: f64,
    },
    DarkPhotonResonance {
        epsilon: f64,
        m_ev: f64,
    },
    TabulatedHeating {
        z_table: Vec<f64>,
        rate_table: Vec<f64>,
    },
    TabulatedPhotonSource {
        z_table: Vec<f64>,
        x_grid: Vec<f64>,
        source_2d: Vec<Vec<f64>>,
    },
    Custom(Box<dyn Fn(f64, &Cosmology) -> f64 + Send + Sync>),
}

Energy injection scenario specification.

Variants (Non-exhaustive)

Non-exhaustive enums could have additional variants added in future. Therefore, when matching against variants of non-exhaustive enums, an extra wildcard arm must be added to account for any future variants.

SingleBurst

Single burst at redshift z_h with fractional energy Δρ/ρ.

Fields

z_h: f64

Central injection redshift.

delta_rho_over_rho: f64

Fractional energy injected into the photon bath Δρ/ρ.

sigma_z: f64

Gaussian width in redshift of the injection profile.

DecayingParticle

Decaying particle with lifetime 1/Γ_X

Fields

f_x: f64

f*_X: energy per baryon released, in eV.

gamma_x: f64

Γ_X: decay rate, in 1/s.

AnnihilatingDM

Annihilating dark matter (s-wave, <σv> = const)

Fields

f_ann: f64

f_ann: energy injection rate parameter [eV/s]. Defined as f_eff × ⟨σv⟩ × m_χ × n_χ,0² / n_H,0 (paper convention). Rate: dE/(dt dV) = f_ann × n_H(z) × (1+z)³. This matches the CosmoTherm convention (e.g. f_ann = 1e-22 eV/s for s-wave).

AnnihilatingDMPWave

Annihilating dark matter (p-wave, <σv> ∝ v² ∝ T ∝ (1+z))

Rate: dE/(dt dV) = f_ann × n_H(z) × (1+z)⁴, an extra (1+z) relative to s-wave capturing the velocity-dependent cross section <σv> ∝ v² ∝ T ∝ (1+z).

Fields

f_ann: f64

f_ann: energy injection rate parameter [eV/s]. Same definition as AnnihilatingDM but includes the present-day value of the velocity-dependent cross section.

MonochromaticPhotonInjection

Monochromatic photon injection/removal at a specific frequency.

Injects ΔN/N photons as a Gaussian in both frequency and redshift, approximating a delta-function injection at (x_inj, z_h). This creates both energy AND photon number perturbations, producing qualitatively different spectral distortions from pure energy injection.

Key physics: injection at x < x₀ ≈ 3.60 produces negative μ.

The heating rate method returns the energy injection rate from the photon injection: d(Δρ/ρ)/dt = (α_ρ × x_inj) × d(ΔN/N)/dt. The frequency-dependent source is applied separately via photon_source_rate.

References: Chluba (2015), arXiv:1506.06582 Arsenadze et al. (2025), arXiv:2409.12940, Appendix C+D

Fields

x_inj: f64

Injection frequency (dimensionless x = hν/(kT_z))

delta_n_over_n: f64

Total fractional photon number to inject: ΔN/N

z_h: f64

Central injection redshift

sigma_z: f64

Gaussian width in redshift

sigma_x: f64

Gaussian width in frequency (should match grid resolution)

DecayingParticlePhoton

Decaying particle X → γγ (spontaneous vacuum decay).

Each decay produces two photons at x_inj(z) = x_inj_0 / (1+z), where x_inj_0 = E_γ/(kT_0) = m_X c²/(2kT_0).

The photon source term follows Bolliet & Chluba (2021), Eq. 3: dn/dt|_inj = G₂ × f_inj × Γ_X × S(z) × G(x, x_inj, σ_x) / x²

where S(z) = exp(-Γ_X × t(z)) is the vacuum survival fraction.

NOTE: Stimulated emission (Bose enhancement) is NOT included. Because X → γγ deposits both final-state photons into the same mode at x_inj, the physical rate carries a factor (1 + n(x_inj))² (one (1+n) per emitted photon, squared because both land in the same occupied mode). This is significant at x_inj ≪ 1 where n_pl ≈ 1/x_inj ≫ 1, and is sub-percent for x_inj ≳ few. See commit 0b6cfa4 for a prototype implementation.

All injected photons are routed through photon_source_rate() and are evolved self-consistently by the PDE solver (Compton + DC/BR).

Reference: Bolliet & Chluba (2021), MNRAS 507, 3148 [arXiv:2012.07292]

Fields

x_inj_0: f64

x_inj,0 = E_γ/(kT_0) = m_X c²/(2kT_0)

f_inj: f64

Dimensionless injection amplitude (Eq. 5 of B&C 2021). f_inj ≈ 2(Ω_cdm/Ω_γ) × f_dm × G₃/(G₂ × x_inj,0)

gamma_x: f64

Γ_X: vacuum decay rate [1/s]

DarkPhotonResonance

Dark photon (γ ↔ A’) resonant conversion in the narrow-width approximation.

Applied as an initial condition at the resonance redshift: Δn(x) = -[1 - exp(-γ_con/x)] × n_pl(x) at z_start = z_res, where γ_con = π ε² m² / (|d ln ω_pl²/d ln a|_{z_res} × T_γ(z_res) × H(z_res)). The solver then evolves this IC with Kompaneets + DC/BR.

The 1/x factor in the conversion probability captures the frequency dependence P(x) ∝ 1/ω for ultrarelativistic photons.

References: Mirizzi, Redondo & Sigl (2009), JCAP 0903, 026 Chluba, Cyr & Johnson (2024), MNRAS 535, 1874 Arsenadze et al. (2025), JHEP 03, 018

Fields

epsilon: f64

Kinetic mixing parameter ε

m_ev: f64

Dark photon mass m_{A’}, in eV.

TabulatedHeating

Tabulated heating rate loaded from a file.

The z_table and rate_table are sorted ascending in z. heating_rate_per_redshift(z) interpolates linearly in log(z), returning 0 outside the table bounds.

This enables Python (or any external tool) to define arbitrary heating rates by writing a CSV table and passing it to the Rust PDE solver.

Fields

z_table: Vec<f64>

Redshift grid (ascending)

rate_table: Vec<f64>

d(Δρ/ρ)/dz at each z (positive = heating)

TabulatedPhotonSource

Tabulated frequency-dependent photon source loaded from a file.

The source function is defined on a 2D grid (z, x), with bilinear interpolation. Returns 0 outside the table bounds.

Fields

z_table: Vec<f64>

Redshift grid (ascending)

x_grid: Vec<f64>

Frequency grid (ascending)

source_2d: Vec<Vec<f64>>

Source rate: source_2d[iz][ix] = d(Δn)/dz at (z_table[iz], x_grid[ix]).

Custom(Box<dyn Fn(f64, &Cosmology) -> f64 + Send + Sync>)

Custom heating function

Implementations

impl InjectionScenario

pub fn name(&self) -> &str

CLI-friendly name for this injection scenario.

pub fn validate(&self) -> Result<(), String>

Validate parameters, returning an error message if invalid.

pub fn heating_rate(&self, z: f64, cosmo: &Cosmology) -> f64

Compute the heating rate d(Δρ_γ/ρ_γ)/dt at redshift z.

Returns the rate in units of [1/s].

pub fn photon_source_rate(&self, x: f64, z: f64, cosmo: &Cosmology) -> f64

Frequency-dependent photon injection/removal rate.

Returns d(Δn)/dt at frequency x, in units of [1/s]. This is the direct modification to the photon occupation number at each frequency, separate from the bulk heating captured by heating_rate.

For MonochromaticPhotonInjection, returns a Gaussian profile in both x and z centered at (x_inj, z_h).

Returns 0.0 for scenarios without frequency-dependent photon injection.

pub fn has_photon_source(&self) -> bool

Whether this scenario has frequency-dependent photon injection/depletion.

pub fn refinement_zones(&self) -> Vec<RefinementZone>

Return refinement zones for adaptive grid resolution near injection features.

Photon injection scenarios need extra grid points near the injection frequency to resolve the narrow Gaussian source profile and the DC/BR absorption/re-emission dynamics at low x.

pub fn characteristic_redshift(&self) -> Option<(f64, f64)>

Return the characteristic injection redshift(s) for this scenario.

For burst-like scenarios, returns Some((z_center, z_upper)) where z_upper is the highest redshift at which injection is active (typically z_h + 5σ_z for Gaussians). z_center is the peak.

For continuous scenarios (decaying particles, annihilation), returns None — injection happens at all z.

pub fn dark_photon_params(&self, cosmo: &Cosmology) -> Option<(f64, f64)>

Dark-photon NWA parameters (γ_con, z_res), if applicable.

Returns Some((γ_con, z_res)) for DarkPhotonResonance, None for all other scenarios. Returns None if the resonance falls outside the supported redshift range.

pub fn initial_delta_n( &self, x_grid: &[f64], cosmo: &Cosmology, ) -> Option<Vec<f64>>

Initial-condition perturbation Δn(x) to be installed at z_start.

Scenarios that deposit their distortion as an impulsive event (notably DarkPhotonResonance) return Some(Δn_init) here; the solver applies it at z_start before beginning redshift evolution. All other scenarios return None and are evolved from Δn = 0.

pub fn suggested_x_min(&self) -> Option<f64>

Suggest a lower x_min for the frequency grid when needed.

Low-frequency photon injection can be artificially absorbed by the Dirichlet boundary at x_min if the source support extends below the existing grid floor.

pub fn warn_strong_distortion(&self) -> Vec<String>

Check for strong distortion regime and return warnings.

The Kompaneets equation solver uses a linearized perturbation approach (Δn = n - n_pl) that breaks down for |Δρ/ρ| > ~0.01. Returns a list of warning messages for scenarios that may exceed this limit.

pub fn warn_dark_photon_range(&self, cosmo: &Cosmology) -> Vec<String>

Warn when the dark-photon NWA resonance falls outside the validated redshift range (roughly z ∈ [50, 3×10⁶] per CLAUDE.md).

Returns Some(Err(msg)) when no resonance exists at all in the supported band (hard error), Some(Ok(warnings)) with a regime warning when z_res lands outside the validated window, and None for non-dark-photon scenarios or when everything is fine.

pub fn warn_tabulated_coverage(&self, z_start: f64, z_end: f64) -> Vec<String>

Warn when a tabulated-source table doesn’t cover the solver’s integration range [z_end, z_start].

interp_log_z / interp_2d return 0.0 outside the table — a silent extrapolation that produces zero injection at redshifts the user may have intended to cover. Surface the mismatch as a warning at build time so it can’t go unnoticed.

pub fn warn_stimulated_emission(&self) -> Vec<String>

Warn if stimulated emission (Bose enhancement) is missing for photon decay.

DecayingParticlePhoton uses the vacuum decay rate only. For the canonical X → γγ channel, the physical rate carries a factor (1 + n(x_inj))² (one (1+n) per emitted photon, both into the same mode at x_inj); for single-photon channels (e.g. X → γ X’) the factor is (1 + n(x_inj)). The squared form is significant at x_inj ≪ 1 where n_pl ≈ 1/x_inj ≫ 1.

pub fn heating_rate_per_redshift(&self, z: f64, cosmo: &Cosmology) -> f64

Compute the physical d(Δρ/ρ)/dz.

Sign convention: For positive energy injection (heating_rate > 0), this returns a NEGATIVE value because energy enters as z decreases (dt > 0, dz < 0). This is the correct physical sign.

WARNING: The Green’s function routines (mu_from_heating, etc.) expect a POSITIVE dq/dz for heating. Use heating_rate_per_redshift().abs() or pass a positive Gaussian directly when calling Green’s function methods.