Cosmology Module

The hod_mod.core sub-package computes the fundamental cosmological quantities that drive the forward model: the matter power spectrum, the halo mass function, halo density profiles, and geometric distances.

—

Distances and Volumes

(hod_mod.core.distances)

The comoving distance to redshift \(z\) is

\[\chi(z) = \frac{c}{H_0} \int_0^z \frac{dz'}{E(z')}\]

where \(E(z) = H(z)/H_0 = [\Omega_m(1+z)^3 + \Omega_\Lambda]^{1/2}\) for a flat \(\Lambda`CDM cosmology. The angular diameter distance is :math:`D_A = \chi/(1+z)\) and the luminosity distance is \(D_L = \chi (1+z)\).

The comoving volume element per steradian is

\[\frac{dV}{dz\,d\Omega} = \frac{c}{H_0} \frac{\chi^2(z)}{E(z)}.\]

Cosmological distances and volumes: JAX-differentiable flat w0waCDM.

All integrals are computed via 256-point Gauss-Legendre quadrature for accuracy and JAX-JIT compatibility. Physical conventions:

• Distances in Mpc (not Mpc/h unless noted)

• \(c = 299\,792.458\) km/s

• Flat geometry (\(\Omega_k = 0\), so \(\Omega_\mathrm{DE} = 1 - \Omega_m\))

Dark energy equation of state — Chevallier-Polarski-Linder (CPL):

\[w(a) = w_0 + w_a (1 - a) = w_0 + w_a \frac{z}{1+z}\]

The dark energy density factor integrates to (Linder 2003):

\[f_\mathrm{DE}(z) = (1+z)^{3(1+w_0+w_a)} \exp\!\left(\frac{-3\,w_a\,z}{1+z}\right)\]

For \(\Lambda\mathrm{CDM}\): \(w_0 = -1,\;w_a = 0 \Rightarrow f_\mathrm{DE} = 1\).

The Hubble function is

\[E(z) = \frac{H(z)}{H_0} = \sqrt{\Omega_m (1+z)^3 + (1-\Omega_m)\,f_\mathrm{DE}(z)}\]

The comoving distance is

\[\chi(z) = \frac{c}{H_0} \int_0^z \frac{dz'}{E(z')}\]

Derived distances (Hogg 2000):

\[D_A(z) = \frac{\chi(z)}{1+z}, \qquad D_L(z) = (1+z)\,\chi(z)\]

For flat geometry the angular diameter distance between \(z_1\) and \(z_2\) is (Hogg 2000, Eq. 19):

\[D_A(z_1, z_2) = \frac{\chi(z_2) - \chi(z_1)}{1 + z_2}\]

Comoving volume element and total:

\[\frac{dV_c}{dz\,d\Omega} = \frac{c}{H_0}\,\frac{\chi^2(z)}{E(z)}, \qquad V_c(<z) = 4\pi \int_0^z \frac{dV_c}{dz'\,d\Omega}\, dz'\]

hod_mod.core.distances.age_of_universe(h: float, omega_m: float, w0: float = -1.0, wa: float = 0.0, z_max: float = 1000.0) → Array

Age of the Universe at redshift 0 [Gyr].

\[t_0 = \frac{1}{H_0} \int_0^\infty \frac{dz}{(1+z)\,E(z)}\]

Integrated to z_max (default 1000) — contribution beyond is negligible.

Parameters:

• h (float) – Dimensionless Hubble constant.

• omega_m (float) – Total matter density.

• w0 (float) – Dark energy equation-of-state today. Default \(-1\).

• wa (float) – CPL time-variation. Default \(0\).

• z_max (float) – Upper integration limit (default 1000).

Returns:

t0 (jnp.ndarray) – Age of the Universe [Gyr].

hod_mod.core.distances.angular_diameter_distance(z: Array, h: float, omega_m: float, w0: float = -1.0, wa: float = 0.0) → Array

Angular diameter distance \(D_A(z) = \chi(z)/(1+z)\) [Mpc].

Parameters:

• z (jnp.ndarray) – Redshift.

• h (float) – Dimensionless Hubble constant.

• omega_m (float) – Total matter density.

• w0 (float) – Dark energy equation-of-state today. Default \(-1\).

• wa (float) – CPL time-variation. Default \(0\).

Returns:

• D_A (jnp.ndarray) – Angular diameter distance [Mpc].

• Accuracy

• ——–

• Exact identity D_A = χ/(1+z) verified to < 1e-6 relative error pointwise

• for z ∈ [0.01, 3] (N=100, 2026-04-23).

• Timing

• ——

• ~ 212 µs / call (JIT-compiled, N=100 redshifts, CPU x86-64, 2026-04-23).

hod_mod.core.distances.angular_diameter_distance_z1z2(z1: Array, z2: Array, h: float, omega_m: float, w0: float = -1.0, wa: float = 0.0) → Array

Angular diameter distance between \(z_1\) and \(z_2\) [Mpc].

For flat geometry (Hogg 2000, Eq. 19):

\[D_A(z_1, z_2) = \frac{\chi(z_2) - \chi(z_1)}{1 + z_2}\]

Parameters:

• z1 (jnp.ndarray) – Near redshift.

• z2 (jnp.ndarray) – Far redshift (\(z_2 > z_1\)).

• h (float) – Dimensionless Hubble constant.

• omega_m (float) – Total matter density.

• w0 (float) – Dark energy equation-of-state today. Default \(-1\).

• wa (float) – CPL time-variation. Default \(0\).

Returns:

D_A12 (jnp.ndarray) – Angular diameter distance between \(z_1\) and \(z_2\) [Mpc].

hod_mod.core.distances.comoving_distance(z: Array, h: float, omega_m: float, w0: float = -1.0, wa: float = 0.0) → Array

Comoving distance \(\chi(z)\) [Mpc].

Evaluated via 256-point Gauss-Legendre quadrature on \([0, z]\):

\[\chi(z) = \frac{c}{H_0} \int_0^z \frac{dz'}{E(z')}\]

Parameters:

• z (jnp.ndarray) – Redshift array.

• h (float) – Dimensionless Hubble constant (\(H_0 = 100\,h\) km/s/Mpc).

• omega_m (float) – Total matter density \(\Omega_m\).

• w0 (float) – Dark energy equation-of-state today. Default \(-1\).

• wa (float) – CPL time-variation. Default \(0\).

Returns:

• chi (jnp.ndarray) – Comoving distance [Mpc].

• Accuracy

• ——–

• χ(z=0) = 0 exactly. χ(z=1) ≈ 3395 Mpc (Planck 2018) agrees with the Pen

• (1999) fitting formula to < 0.2%. 256-point GL quadrature gives < 0.01%

• error vs 4096-point reference (2026-04-23).

• Timing

• ——

• ~ 421 µs / call (JIT-compiled, N=100 redshifts, CPU x86-64, 2026-04-23).

hod_mod.core.distances.comoving_distance_z1z2(z1: Array, z2: Array, h: float, omega_m: float, w0: float = -1.0, wa: float = 0.0) → Array

Line-of-sight comoving distance between redshifts \(z_1\) and \(z_2\) [Mpc].

For flat geometry (Hogg 2000, §2):

\[D_C(z_1, z_2) = \chi(z_2) - \chi(z_1)\]

Parameters:

• z1 (jnp.ndarray) – Near redshift.

• z2 (jnp.ndarray) – Far redshift (\(z_2 > z_1\)).

• h (float) – Dimensionless Hubble constant.

• omega_m (float) – Total matter density.

• w0 (float) – Dark energy equation-of-state today. Default \(-1\).

• wa (float) – CPL time-variation. Default \(0\).

Returns:

D_C12 (jnp.ndarray) – Comoving distance between the two redshifts [Mpc].

hod_mod.core.distances.comoving_volume(z: Array, h: float, omega_m: float, w0: float = -1.0, wa: float = 0.0) → Array

Total comoving volume within redshift \(z\) [Mpc^3].

\[V_c(<z) = 4\pi \int_0^z \frac{c}{H_0}\,\frac{\chi^2(z')}{E(z')}\,dz'\]

Parameters:

• z (jnp.ndarray) – Redshift.

• h (float) – Dimensionless Hubble constant.

• omega_m (float) – Total matter density.

• w0 (float) – Dark energy equation-of-state today. Default \(-1\).

• wa (float) – CPL time-variation. Default \(0\).

Returns:

• Vc (jnp.ndarray) – Comoving volume [\({\rm Mpc}^3\)].

• Accuracy

• ——–

• V_c(z=1) ≈ 1.6 × 10¹¹ Mpc³ (Planck 2018); monotonically increasing.

• Agrees with comoving_volume_element numerical integration to < 0.1%

• for z ∈ [0.1, 3] (256-point GL, 2026-04-23).

• Timing

• ——

• ~ 4.7 ms / call (JIT-compiled, N=100 redshifts, CPU x86-64, 2026-04-23).

hod_mod.core.distances.comoving_volume_element(z: Array, h: float, omega_m: float, w0: float = -1.0, wa: float = 0.0) → Array

Comoving volume element per steradian per unit redshift [Mpc^3 / sr].

Hogg (2000) Eq. 28:

\[\frac{dV_c}{dz\,d\Omega} = \frac{c}{H_0}\,\frac{(1+z)^2 D_A^2(z)}{E(z)}\]

For flat geometry \(D_A = \chi/(1+z)\), so this reduces to:

\[\frac{dV_c}{dz\,d\Omega} = \frac{c}{H_0}\,\frac{\chi^2(z)}{E(z)}\]

Parameters:

• z (jnp.ndarray) – Redshift.

• h (float) – Dimensionless Hubble constant.

• omega_m (float) – Total matter density.

• w0 (float) – Dark energy equation-of-state today. Default \(-1\).

• wa (float) – CPL time-variation. Default \(0\).

Returns:

dVdzdOmega (jnp.ndarray) – Comoving volume element [\({\rm Mpc}^3/{\rm sr}\)].

hod_mod.core.distances.distance_modulus(z: Array, h: float, omega_m: float, w0: float = -1.0, wa: float = 0.0) → Array

Distance modulus \(\mu(z) = 5\log_{10}[D_L/{\rm Mpc}] + 25\) [mag].

Parameters:

• z (jnp.ndarray) – Redshift.

• h (float) – Dimensionless Hubble constant.

• omega_m (float) – Total matter density.

• w0 (float) – Dark energy equation-of-state today. Default \(-1\).

• wa (float) – CPL time-variation. Default \(0\).

Returns:

• mu (jnp.ndarray) – Distance modulus [mag].

• Accuracy

• ——–

• Identity μ = 5 log₁₀(D_L / Mpc) + 25 verified to < 1e-4 mag for

• z ∈ [0.01, 3] against luminosity_distance (N=100, 2026-04-23).

• At z=0.1 (Planck 2018) (μ ≈ 38.3 mag (Hubble diagram anchor).)

• Timing

• ——

• ~ 586 µs / call (JIT-compiled, N=100 redshifts, CPU x86-64, 2026-04-23).

hod_mod.core.distances.hubble_e(z: Array, omega_m: float, w0: float = -1.0, wa: float = 0.0) → Array

Dimensionless Hubble function \(E(z) = H(z)/H_0\).

Flat wCDM with CPL dark energy (Chevallier & Polarski 2001; Linder 2003):

\[E^2(z) = \Omega_m(1+z)^3 + (1-\Omega_m)\,(1+z)^{3(1+w_0+w_a)} \exp\!\left(\frac{-3\,w_a\,z}{1+z}\right)\]

Parameters:

• z (jnp.ndarray) – Redshift.

• omega_m (float) – Total matter density parameter \(\Omega_m\).

• w0 (float) – Dark energy equation-of-state today. Default \(-1\) (ΛCDM).

• wa (float) – CPL time-variation parameter. Default \(0\) (ΛCDM).

Returns:

• E (jnp.ndarray) – \(H(z)/H_0\).

• Accuracy

• ——–

• E(z=0) = 1.0 exactly (flat ΛCDM (Ω_m + Ω_Λ = 1, f_DE(0) = 1).)

• Timing

• ——

• ~ 19 µs / call (JIT-compiled, N=100 redshifts, CPU x86-64, 2026-04-23).

hod_mod.core.distances.lookback_time(z: Array, h: float, omega_m: float, w0: float = -1.0, wa: float = 0.0) → Array

Lookback time \(t_L(z)\) [Gyr].

\[t_L(z) = \frac{1}{H_0} \int_0^z \frac{dz'}{(1+z')\,E(z')}\]

Parameters:

• z (jnp.ndarray) – Redshift.

• h (float) – Dimensionless Hubble constant.

• omega_m (float) – Total matter density.

• w0 (float) – Dark energy equation-of-state today. Default \(-1\).

• wa (float) – CPL time-variation. Default \(0\).

Returns:

t_L (jnp.ndarray) – Lookback time [Gyr].

hod_mod.core.distances.luminosity_distance(z: Array, h: float, omega_m: float, w0: float = -1.0, wa: float = 0.0) → Array

Luminosity distance \(D_L(z) = (1+z)\,\chi(z)\) [Mpc].

Parameters:

• z (jnp.ndarray) – Redshift.

• h (float) – Dimensionless Hubble constant.

• omega_m (float) – Total matter density.

• w0 (float) – Dark energy equation-of-state today. Default \(-1\).

• wa (float) – CPL time-variation. Default \(0\).

Returns:

• D_L (jnp.ndarray) – Luminosity distance [Mpc].

• Accuracy

• ——–

• Exact identity D_L = χ(1+z) verified to < 1e-6 relative error pointwise

• for z ∈ [0.01, 3] (N=100, 2026-04-23). Consistent with Etherington (1933)

• reciprocity relation D_L = (1+z)² D_A.

• Timing

• ——

• ~ 70 µs / call (JIT-compiled, N=100 redshifts, CPU x86-64, 2026-04-23).

—

Linear Power Spectrum

(hod_mod.core.power_spectrum)

The dimensionless power spectrum is

\[\Delta^2(k, z) \equiv \frac{k^3 P(k, z)}{2\pi^2}.\]

For \(\Delta^2 \ll 1\) the density field is in the linear regime.

CAMB backend (default)

The full Boltzmann code CAMB (Lewis, Challinor & Lasenby 2000) [Lewis2002] is invoked via LinearPowerSpectrum.pk_linear(k, z, theta). A single CAMB run takes ~30 s; for MCMC use the CachedPkLinear wrapper that interpolates on a pre-computed \(k\)-grid keyed on \((\Omega_m, h, \ln 10^{10}A_s, z)\).

Eisenstein-Hu 1998 fitting formula (fast)

Eisenstein & Hu 1998 [EisensteinHu1998] provide an accurate analytical approximation to the transfer function:

\[P_{\rm EH}(k) \propto k^{n_s} T^2(k)\]

where the transfer function \(T(k)\) captures baryonic acoustic oscillations through a fitting formula involving the baryon-to-matter ratio, the matter-radiation equality scale, and the Silk damping scale. Implemented as eisenstein_hu_pk(k, theta) in JAX; differentiable with respect to all cosmological parameters.

Growth factor

Linear growth is encoded as \(P(k,z) = D^2(z) P(k,0)\) with the growth factor

\[D(z) = \frac{5\Omega_m}{2} H(z) / H_0 \int_z^\infty \frac{(1+z')}{[H(z')/H_0]^3}\,dz'.\]

Linear matter power spectrum via CAMB (Lewis, Challinor & Lasenby 2000).

class hod_mod.core.power_spectrum.CsstLinearPowerSpectrum

Bases: object

Linear P(k, z) via the CSST CEmulator (Chen+2025, v2.0).

The emulator is initialised once on instantiation. Cosmology is set via set_cosmos before each call to pk_linear. k is in h/Mpc and output P(k) is in (Mpc/h)^3.

Parameter ranges (will raise ValueError if exceeded):

• Omega_b ∈ [0.04, 0.06]

• Omega_m ∈ [0.24, 0.40]

• H0 ∈ [60, 80] (inferred as h * 100)

• n_s ∈ [0.92, 1.00]

• A_s ∈ [1.7e-9, 2.5e-9]

• w0 ∈ [−1.3, −0.7]

• wa ∈ [−0.5, 0.5]

• m_nu ∈ [0, 0.3] eV (default 0.06 eV when not in theta)

pk_linear(k: Array, z: float, theta: dict) → Array

Linear P(k) [(Mpc/h)^3] at redshift z via the CSST emulator.

Parameters:

• k (array_like, h/Mpc — interpolated onto emulator k-grid)

• z (float — must lie in [0, 3])

• theta (dict — hod_mod cosmological parameter dict)

class hod_mod.core.power_spectrum.LinearPowerSpectrum

Bases: object

Linear P(k, z) computed with CAMB.

Parameters:

(none — no pre-trained weights required)

static default_cosmology() → dict

Planck 2018 TT,TE,EE+lowE best-fit values (flat ΛCDM).

pk_linear(k: Array, z: float, theta: dict) → Array

Linear P(k) [(Mpc/h)^3] at redshift z (total matter).

Supports CPL dark energy via w0 and wa keys in theta (defaults to ΛCDM if absent). CAMB uses the PPF dark energy model (Hu & Sawicki 2007) which remains accurate for \(w < -1\).

Parameters:

• k (array_like, h/Mpc)

• z (float)

• theta (dict — keys: h, Omega_b, Omega_cdm, n_s, ln10^{10}A_s,) – w0 (default -1), wa (default 0)

pk_linear_b(k: Array, z: float, theta: dict) → Array

Baryon auto-power spectrum P_b(k) [(Mpc/h)^3] at redshift z.

\[P_b(k) = T_b^2(k)\,P_{\rm prim}(k)\]

Baryonic BAO wiggles are more pronounced here than in the total spectrum.

Parameters:

• k (array_like, h/Mpc)

• z (float)

• theta (dict — same keys as pk_linear())

pk_linear_cdm(k: Array, z: float, theta: dict) → Array

CDM auto-power spectrum P_CDM(k) [(Mpc/h)^3] at redshift z.

\[P_{\rm CDM}(k) = T_{\rm cdm}^2(k)\,P_{\rm prim}(k)\]

where \(T_{\rm cdm}\) is the CAMB CDM transfer function. To recover total matter: \(f_b P_b + f_c P_{\rm CDM} \approx P_{\rm tot}\) (exact in linear theory when cross-spectrum terms dominate, valid within ~5%).

Parameters:

• k (array_like, h/Mpc)

• z (float)

• theta (dict — same keys as pk_linear())

hod_mod.core.power_spectrum.eisenstein_hu_pk(k: Array, theta: dict) → Array

Eisenstein & Hu (1998) transfer function with BAO wiggles.

Implements eqs. (2)–(7), (10)–(24) of EH98. The total transfer function is \(T(k) = f_b T_b(k) + f_c T_c(k)\) and the power spectrum is \(P(k) \propto k^{n_s} T(k)^2\), normalised to unity at \(k = 0.05\,h\,\mathrm{Mpc}^{-1}\).

Parameters:

• k (array_like, h/Mpc)

• theta (dict — keys: h, Omega_m, Omega_b, n_s; optional T_cmb (K, default 2.7255))

• Accuracy

• ——–

• Normalised to P(k=0.05 h/Mpc) = 1.0 by construction. Shape agrees with

• CAMB P(k) to < 10% rms for k ∈ [0.01, 0.3] h/Mpc (Planck 2018, z=0).

• Large-scale slope d log P / d log k ≈ n_s to < 0.15 for k < 0.003 h/Mpc

• (2026-04-23).

• Timing

• ——

• ~ 242 µs / call (JIT-compiled, N=200 wavenumbers, CPU x86-64, 2026-04-23).

hod_mod.core.power_spectrum.eisenstein_hu_pk_nowiggle(k: Array, theta: dict) → Array

Eisenstein & Hu (1998) no-wiggle (smooth) power spectrum.

Implements eqs. (26), (28)–(31) of EH98. Captures the baryon-induced shape suppression through an effective shape parameter \(\Gamma_{\rm eff}(k)\), without acoustic oscillations. Useful as a smooth reference spectrum.

\(P(k) \propto k^{n_s} T_0(q_{\rm eff})^2\), normalised to unity at \(k = 0.05\,h\,\mathrm{Mpc}^{-1}\).

Parameters:

• k (array_like, h/Mpc)

• theta (dict — keys: h, Omega_m, Omega_b, n_s; optional T_cmb (K, default 2.7255))

hod_mod.core.power_spectrum.eisenstein_hu_pk_phys(k: Array, theta: dict) → Array

Eisenstein & Hu (1998) matter power spectrum in physical (Mpc/h)³ units.

Same transfer function as eisenstein_hu_pk() but returns \(P(k)\) in physical \((\mathrm{Mpc}/h)^3\) units with the correct amplitude derived from the primordial curvature spectrum via the Poisson equation in conformal-Newtonian gauge:

\[P(k_h, z=0) = D^2(z=0)\,\frac{8\pi^2}{25} \frac{(c/H_{100})^4}{\Omega_m^2} A_s \left(\frac{h}{k_*}\right)^{n_s-1} k_h^{n_s}\,T^2(k_h)\]

where \(c/H_{100} = 2997.924\;\mathrm{Mpc}/h\), \(k_* = 0.05\;\mathrm{Mpc}^{-1}\) (CAMB pivot, physical), and \(A_s = e^{\ln10^{10}A_s}\times10^{-10}\).

\(D(z=0)\) is the linear growth factor normalised so that \(D \to a\) during matter domination (EdS limit). For \(\Omega_m < 1\) flat \(\Lambda\)CDM, \(D(z=0) < 1\) because dark energy suppresses growth after matter–\(\Lambda\) equality. This factor is not encoded in the EH98 transfer function shape and must be included in the amplitude. It is computed via the exact numerical integral

\[D(z=0) = \frac{5\Omega_m}{2} \int_0^1 \frac{\mathrm{d}a}{[a\,H(a)/H_0]^3}\]

with \(H(a) = H_0\sqrt{\Omega_m a^{-3} + 1 - \Omega_m}\).

Parameters:

• k (array_like, h/Mpc)

• theta (dict — keys: h, Omega_m, Omega_b, n_s, ln10^{10}A_s;) – optional Omega_cdm (defaults to Omega_m − Omega_b), optional T_cmb [K] (default 2.7255)

Returns:

• P(k) (jnp.ndarray, (Mpc/h)³)

• Accuracy

• ——–

• \(\sigma_8\) computed from this spectrum agrees with CAMB to

• \(< 1\%\) for Planck 2018 parameters. Residual discrepancies

• reflect the EH98 transfer-function shape error (not the amplitude

• formula).

hod_mod.core.power_spectrum.rho_critical_0() → float

Critical matter density at z=0 for H₀ = 100 km/s/Mpc, in h-units.

\[\rho_{\mathrm{crit},0} = \frac{3H_{100}^2}{8\pi G} \approx 2.775\times10^{11}\;(M_\odot/h)\,(\mathrm{Mpc}/h)^{-3}\]

In h-unit conventions the h² from \(H_0 = 100h\) km/s/Mpc cancels the \(h^{-3}\) from the comoving volume, so this quantity is independent of h. The mean matter density follows as \(\bar{\rho}_m = \Omega_m\,\rho_{\mathrm{crit},0}\).

Physical constants used: \(G = 6.67430\times10^{-11}\) m³ kg⁻¹ s⁻², 1 Mpc = 3.085677581×10²² m, 1 M⊙ = 1.989×10³⁰ kg.

Returns:

rho_crit0 (float, (Msun/h) / (Mpc/h)³)

—

Non-linear Power Spectrum

(hod_mod.core.nonlinear)

The non-linear power spectrum \(P_{\rm nl}(k,z)\) is computed by the Aletheia emulator (Contreras et al. 2023) [Aletheia2025], a neural-network emulator trained on a suite of N-body simulations spanning a wide cosmological parameter space including massive neutrinos and dynamical dark energy.

Non-linear matter power spectrum via the Aletheia or CSST CEmulator, or CAMB HALOFIT.

class hod_mod.core.nonlinear.CachedPkNonlinear(pk_nl_obj, n_k: int = 512)

Bases: object

Caching wrapper for any non-linear power spectrum backend.

Duck-types any object exposing pk_nonlinear(k, z, theta) -> array_like and caches the result on a fixed log-spaced k grid, keyed by (z, Omega_m, ln10As, h). Subsequent calls with the same cosmology are returned via cheap log-log interpolation.

Compatible backends:

• NonLinearPowerSpectrum with backend='aletheia' — Aletheia emulator (arXiv:2511.13826), fast and JAX-based.

• NonLinearPowerSpectrum with backend='csst' — CEmulator v2.0.

• HALOFITSpectrum — CAMB HALOFIT / HMcode variants.

This mirrors the CachedPkLinear pattern documented in CLAUDE.md; use it in MCMC hot loops to avoid repeated CAMB / emulator calls.

Parameters:

• pk_nl_obj (object) – Any instance with a pk_nonlinear(k, z, theta) method.

• n_k (int) – Number of points in the internal interpolation grid (k ∈ [10⁻⁴, 20] h/Mpc).

pk_nonlinear(k, z: float, theta: dict) → Array

Non-linear P_nl(k) [(Mpc/h)^3], log-log interpolated from cache.

Parameters:

• k (array_like [h/Mpc])

• z (float)

• theta (dict — hod_mod cosmological parameter dict)

Returns:

jnp.ndarray

class hod_mod.core.nonlinear.HALOFITSpectrum(halofit_version: str = 'mead2020')

Bases: object

Non-linear P(k, z) via CAMB’s built-in HALOFIT / HMcode variants.

Uses CAMB (already installed) with NonLinear = NonLinear_pk. The default variant is halofit_mead2020 (HMcode-2020, arXiv:2009.01858), which includes a baryonic feedback option.

Parameters:

halofit_version (str) – CAMB HALOFIT variant name. Common choices:

"mead2020" — HMcode-2020 (Mead+2021, arXiv:2009.01858) default "mead2020_feedback" — HMcode-2020 with baryonic feedback "takahashi" — Takahashi+2012 (arXiv:1208.2701) "mead" — HMcode-2015 (arXiv:1602.02154) "original" — Smith+2003 original HALOFIT (arXiv:astro-ph/0207664)

Notes

Each call runs a full CAMB evaluation (~1–5 s). Use the CachedPkLinear pattern from CLAUDE.md to cache results in hot loops. Gradients through this class are not available.

pk_nonlinear(k: ndarray, z: float, theta: dict) → Array

Non-linear P(k) [(h^{-1} Mpc)^3] from CAMB HALOFIT.

Parameters:

• k (array_like, [h/Mpc])

• z (float)

• theta (dict — hod_mod cosmological parameter dict)

Returns:

jnp.ndarray, shape (len(k),) – P_nl(k) in (h^{-1} Mpc)^3, interpolated from the CAMB output grid.

class hod_mod.core.nonlinear.NonLinearPowerSpectrum(backend: str = 'aletheia')

Bases: object

Non-linear P(k, z) via the Aletheia or CSST CEmulator.

Both backends are numpy-based; this class loads the requested emulator once and exposes a JAX-compatible interface by converting outputs to jnp arrays. Gradients through the emulator are not available — use finite differences.

Parameters:

backend ({“aletheia”, “csst”}) – "aletheia" — Sanchez 2025 (arXiv:2511.13826), valid k ∈ [0.006, 2] Mpc^-1. "csst" — Chen+2025 CEmulator v2.0, valid k ∈ [0.005, 10] h/Mpc, z ∈ [0, 3], nonlinear via HMcode-2020 boost ratio.

boost_factor(k: Array, z: float, theta: dict, pk_lin: Array) → Array

Non-linear boost B(k, z) = P_nl / P_lin.

pk_nonlinear(k: ndarray, z: float, theta: dict) → Array

Non-linear P(k) [(Mpc/h)^3] at redshift z.

Parameters:

• k (array_like, h/Mpc)

• z (float)

• theta (dict — hod_mod cosmological parameter dict)

Notes

For the aletheia backend the emulator is valid for k ∈ [0.006, 2.0] Mpc⁻¹. k values outside this range fall back to the Eisenstein-Hu linear spectrum multiplied by the nonlinear boost evaluated at the nearest valid boundary. This ensures the function is defined over any k grid used by the clustering module without raising a ValueError from the emulator.

pk_nonlinear_jax(k: Array, z: float, theta: dict) → Array

JAX-native P_nl(k) [(Mpc/h)^3] for the Aletheia backend.

Unlike pk_nonlinear(), this method keeps all array operations in JAX so that jax.grad / jax.jacobian can differentiate through the emulator with respect to both k and the HOD/cosmological parameters encoded inside theta.

Only available for backend='aletheia'. For CSST or HMcode use CachedPkNonlinear wrapping the numpy-based method.

Parameters:

• k (jnp.ndarray [h/Mpc])

• z (float)

• theta (dict — hod_mod cosmological parameter dict)

Returns:

jnp.ndarray [(Mpc/h)^3]

class hod_mod.core.nonlinear.WHMSpectrum(whm_version: str = 'brieden2023')

Bases: object

Non-linear P(k, z) via the Web-Halo Model (WHM, Brieden et al. 2025).

Requires the WHM-CAMB fork installed in place of (or alongside) standard CAMB:

git clone https://github.com/SamuelBrieden/WHM
pip install -e WHM/WHM-CAMB

The WHM is parameter-free and combines the halo model with perturbation theory, modelling haloes, filaments, and sheets via physically-motivated window functions. Achieves ~2 % accuracy up to k = 0.4 h/Mpc⁻¹ across a wide range of cosmologies (arXiv:2508.10902).

Parameters:

whm_version (str) – WHM variant. Default 'brieden2023'. Available choices:

'brieden2023' — baseline WHM (recommended) 'brieden2023_feedback' — WHM with baryonic feedback 'brieden2023_halo' — halo term only 'brieden2023_fila' — filament term only 'brieden2023_sheet' — sheet term only 'brieden2023_halosphere', 'brieden2023_filasphere', 'brieden2023_sheetsphere' — spherical-profile variants

Notes

Each call runs a full CAMB evaluation (~1–5 s). Wrap with CachedPkNonlinear in MCMC hot loops. Gradients are not available through this class.

pk_nonlinear(k: ndarray, z: float, theta: dict) → Array

Non-linear P(k) [(h^{-1} Mpc)^3] from the WHM via CAMB.

Parameters:

• k (array_like [h/Mpc])

• z (float)

• theta (dict — hod_mod cosmological parameter dict)

Returns:

jnp.ndarray, shape (len(k),) – P_nl(k) in (h^{-1} Mpc)^3, log-log interpolated from the CAMB grid.

—

Halo Mass Function

(hod_mod.core.halo_mass_function)

The comoving number density of halos per unit logarithmic mass is

\[\frac{dn}{d\ln M} = f(\sigma) \frac{\bar{\rho}_m}{M} \left|\frac{d\ln\sigma^{-1}}{d\ln M}\right|\]

where \(\bar{\rho}_m\) is the mean comoving matter density and \(\sigma^2(M, z)\) is the variance of the linear density field smoothed on the Lagrangian radius \(R = (3M/4\pi\bar{\rho}_m)^{1/3}\):

\[\sigma^2(M, z) = \frac{D^2(z)}{2\pi^2} \int_0^\infty P_{\rm lin}(k, 0)\, W^2(kR)\, k^2\, dk\]

with the top-hat window function \(W(x) = 3(\sin x - x\cos x)/x^3\).

Tinker+2008 multiplicity function (library default; the project’s fitting pipelines use "csst" as their baseline instead — see Emulator backends below)

Tinker et al. 2008 [Tinker2008] calibrated the multiplicity function \(f(\sigma)\) against N-body simulations for halos defined by a fixed overdensity \(\Delta = 200\) relative to the mean background:

\[f(\sigma) = A\left[1 + \left(\frac{\sigma}{b}\right)^{-a}\right] \exp\left(-\frac{c}{\sigma^2}\right)\]

with best-fit parameters \(A=0.186\), \(a=1.47\), \(b=2.57\), \(c=1.19\) at \(z=0\). Redshift evolution of the parameters is also given in Table 2 of that paper.

Halo bias

The linear halo bias relates the halo overdensity to the matter overdensity on large scales. Using the Tinker+2010 prescription (Tinker et al. 2010) [Tinker2010]:

\[b(M, z) = 1 - A_b \frac{\nu^{a_b}}{\nu^{a_b} + \delta_c^{a_b}} + B_b\,\nu^{b_b} + C_b\,\nu^{c_b}\]

where \(\nu = \delta_c / \sigma(M,z)\) is the peak height and \(\delta_c \approx 1.686\) is the linear collapse threshold.

The effective galaxy bias is obtained by weighting over the occupation-weighted HMF:

\[b_{\rm eff} = \frac{\int b(M)\,\langle N(M)\rangle\,\frac{dn}{dM}\,dM} {\int \langle N(M)\rangle\,\frac{dn}{dM}\,dM}\]

Analytic backends (all JAX-native, differentiable)

make_hmf(backend, pk_func=pk_lin.pk_linear) accepts any key in _FSIGMA_MODELS, including "tinker08" (the library’s dependency-free default when no backend is requested), "bocquet16", "yung25", and 14 others.

Emulator backends

Two optional emulator backends replace the analytic HMF with a Gaussian-Process or neural-net prediction trained on N-body suites. Both expose the same .dndm(), .sigma(), .bias(), .n_eff() interface and are selected via make_hmf(backend). Halo bias always falls back to Tinker+2010 (neither emulator provides it).

The project’s fitting pipelines (hod_mod/scripts/fitting/*.py and the corresponding configs/hod_fit_*.yml / configs/fitting/*_example.yml) use "csst" (CSSTEMU) as their baseline HMF backend rather than the library default "tinker08", since its wide calibration range (\(M \in [10^{10},10^{16}]\,M_\odot/h\), \(z \leq 3\)) covers the full mass range these pipelines integrate over. "aemulusnu" is not used as a pipeline baseline: it is only calibrated for \(M \geq 10^{13}\,M_\odot/h\), so HOD samples with non-negligible occupation below that mass (e.g. low stellar-mass-threshold samples) get silently extrapolated, dominating the integrated predictions with unreliable values (see AemulusNuHaloMassFunction). Literature-validation scripts (validate_*.py, configs/benchmarks/*.yml) intentionally stay on "tinker08" to match the HMF used by the papers they reproduce.

Backend key	Class	Paper	Calibration range	Extra dependency
"csst"	CsstHaloMassFunction	[ChenCSST2025] (SCPMA 2025)	\(M \in [10^{10},10^{16}]\,M_\odot/h\), \(z \leq 3\), \(\Omega_m,h,n_s,A_s,w_0,w_a,m_\nu,\Omega_b\)	pip install git+https://github.com/czymh/csstemu
"aemulusnu"	AemulusNuHaloMassFunction	[ShenAemulus2025] (JCAP 2025)	\(M \geq 10^{13}\,M_\odot/h\), \(z \leq 2\), \(w_0w_a\nu{\rm CDM}\)	pip install git+https://github.com/DelonShen/aemulusnu_hmf

Both packages have numpy ≥ 2.0 and scipy ≥ 1.11 incompatibilities that are patched automatically at import time by halo_mass_function.py (scipy.integrate.simps alias, np.trapz alias, GPR predict squeeze, Cosmology.get_Omegam scalar fix).

Usage:

from hod_mod.core.power_spectrum import LinearPowerSpectrum
from hod_mod.core.halo_mass_function import make_hmf

pk  = LinearPowerSpectrum()

# Library default: analytic, differentiable, no extra dependency
hmf = make_hmf("tinker08", pk_func=pk.pk_linear)

# CSST emulator (pk_func ignored; uses CSST PkLin internally) —
# this is the baseline backend used by the fitting pipelines
hmf = make_hmf("csst")

# Aemulus-ν emulator (requires pk_func for bias; M >= 1e13 Msun/h only)
hmf = make_hmf("aemulusnu", pk_func=pk.pk_linear)

# All expose the same interface:
dn  = hmf.dndm(m_grid, z, theta)   # [h^4 Mpc^{-3} (M_sun/h)^{-1}]
b   = hmf.bias(m_grid, z, theta)   # Tinker 2010 (all backends)

Halo mass functions: multiple models implemented in JAX, following Colossus conventions.

The halo mass function describes the comoving number density of halos per unit mass interval:

\[\frac{dn}{dM} = f(\sigma)\, \frac{\bar{\rho}_0}{M^2}\, \left|\frac{d\ln\sigma}{d\ln M}\right|\]

where \(\bar{\rho}_0\) is the mean comoving matter density at z=0, \(\sigma(M)\) is the RMS linear density fluctuation in a sphere of radius \(R(M)\), and \(f(\sigma)\) is the multiplicity function.

The variance \(\sigma^2(M)\) is computed from the linear power spectrum:

\[\sigma^2(M) = \frac{1}{2\pi^2} \int_0^\infty P(k)\, W^2(kR)\, k^2\, dk\]

with top-hat window \(W(x) = 3(\sin x - x\cos x)/x^3\) and \(R(M) = (3M / 4\pi\bar{\rho}_0)^{1/3}\).

Redshift evolution of \(\sigma\): \(\sigma(M,z) = \sigma(M,0) \times D(z)/D(0)\) where D(z) is the linear growth factor (flat ΛCDM approximation).

All masses in M_sun/h, distances in Mpc/h, ρ̄₀ in (M_sun/h)/(Mpc/h)^3.

class hod_mod.core.halo_mass_function.AemulusNuHaloMassFunction(pk_func=None, rho_mean: float = 86010827748.85518, Delta: float = 200.0)

Bases: object

Halo mass function from the Aemulus-ν emulator (Shen+2025, arXiv:2410.00913).

Wraps aemulusnu_hmf.emulator.dn_dM and exposes the same interface as HaloMassFunction.

Calibration range: M ≥ 10^13 M_sun/h, z ∈ [0, 2], wνCDM cosmologies (w0waCDM + neutrino mass). A UserWarning is issued if dndm is called with masses below 10^13 M_sun/h, but the mass array is not clipped — values below the calibrated floor are passed straight to the GP emulator and silently extrapolated, which can return wildly wrong (huge, near-zero, or non-monotonic) dn/dM values.

Caveat for HOD predictions: FullHaloModelPrediction, HODBase and HaloModelCrossSpectra all integrate over a fixed M = 10^10–10^16 M_sun/h grid regardless of which HMF backend is plugged in. For galaxy samples whose central occupation is non-negligible below 10^13 M_sun/h — e.g. low stellar-mass-threshold samples such as Comparat+2025 S1 (log10m_star_thresh ≈ 9.5) — the resulting wp/n_gal/delta_sigma/C_ell^{gX} predictions are dominated by this unreliable extrapolated region and can look dramatically different from analytic backends like tinker08. This backend is only recommended for cluster-scale / high-mass-threshold samples where the occupation is concentrated above 10^13 M_sun/h.

Bias falls back to Tinker 2010 using σ(M) from a linear power spectrum callable supplied at construction time (pk_func).

Parameters:

• pk_func (callable) – (k, z, theta) → P_lin(k) [(Mpc/h)^3]. Required for sigma() and bias().

• rho_mean (float) – Mean comoving matter density at z=0 [M_sun/h / (Mpc/h)^3].

• Delta (float) – Overdensity for the Tinker 2010 bias (default 200, mean density).

bias(m_h: Array, z: float, theta: dict) → Array

Tinker 2010 large-scale halo bias b(M, z) [dimensionless].

dndm(m_h: Array, z: float, theta: dict) → Array

Halo mass function dn/dM [h^4 Mpc^{-3} (M_sun/h)^{-1}].

Parameters:

• m_h (jnp.ndarray Halo masses [M_sun/h].)

• z (float Redshift (must be ≤ 2).)

• theta (dict Cosmological parameters.)

n_eff(m_min: float, m_max: float, z: float, theta: dict) → Array

Integrated number density n(m_min < M < m_max) [h^3 Mpc^{-3}].

sigma(m_h: Array, z: float, theta: dict) → Array

RMS linear density fluctuation σ(M, z) via the supplied pk_func.

class hod_mod.core.halo_mass_function.CsstHaloMassFunction(massdef: str = 'RockstarM200m', rho_mean: float = 86010827748.85518)

Bases: object

Halo mass function from the CSST CEmulator (Chen+2025, v2.0).

Wraps HMF_CEmulator.get_dndlnM and exposes the same interface as HaloMassFunction so it can be used interchangeably.

Mass definition: RockstarM200m (200× mean, Rockstar halo finder). Bias is computed from the Tinker 2010 formula applied to σ(M) derived from the CSST linear power spectrum emulator.

Parameters:

• massdef ({“RockstarM200m”, “FoFM200c”, “RockstarMvir”})

• rho_mean (float) – Mean comoving matter density at z=0 [M_sun/h / (Mpc/h)^3].

bias(m_h: Array, z: float, theta: dict) → Array

Tinker 2010 large-scale halo bias b(M, z).

dndm(m_h: Array, z: float, theta: dict) → Array

Halo mass function dn/dM [h^4 Mpc^{-3} (M_sun/h)^{-1}].

Converts from the emulator’s dn/dlnM output: dn/dM = (dn/dlnM) / M.

n_eff(m_min: float, m_max: float, z: float, theta: dict) → Array

Integrated number density n(m_min < M < m_max) [h^3 Mpc^{-3}].

sigma(m_h: Array, z: float, theta: dict) → Array

RMS linear density fluctuation σ(M, z) via the CSST linear P(k) (CAMB backend).

class hod_mod.core.halo_mass_function.HaloMassFunction(pk_func, rho_mean: float = 86010827748.85518, model: str = 'tinker08', Delta: float = 200.0, n_k: int = 512, **fsigma_kwargs)

Bases: object

Halo mass function dn/dM and halo bias, JAX-accelerated.

Computes σ(M) from the linear power spectrum at z=0 (normalized to sigma8 if provided), applies the linear growth factor for z-evolution, then evaluates the chosen multiplicity function.

Parameters:

• pk_func (callable) – (k, z, theta) → P_lin(k) [(Mpc/h)^3]. Used only at z=0 to build σ(M).

• rho_mean (float) – Mean comoving matter density at z=0 [M_sun/h / (Mpc/h)^3]. Default: rho_critical_0() × 0.3100 (Planck 2018 Ω_m).

• model (str) – Multiplicity function. Any key from _FSIGMA_MODELS.

• Delta (float) – Overdensity threshold w.r.t. mean density (used by tinker08 and bias).

• n_k (int) – Number of points in the k integration grid.

• **fsigma_kwargs – Extra keyword arguments forwarded to the multiplicity function via functools.partial at construction time. Examples:

  • hydro=True for bocquet16

  • delta_ratio=<float> for despali16

bias(m_h: Array, z: float, theta: dict) → Array

Tinker 2010 large-scale halo bias b(M, z) [dimensionless].

Parameters:

• m_h (jnp.ndarray Halo masses [M_sun/h].)

• z (float Redshift.)

• theta (dict Cosmological parameters.)

dndm(m_h: Array, z: float, theta: dict) → Array

Halo mass function dn/dM [h^4 Mpc^{-3} (M_sun/h)^{-1}].

\[\frac{dn}{dM} = f(\sigma)\, \frac{\bar{\rho}_0}{M^2}\, \left|\frac{d\ln\sigma}{d\ln M}\right|\]

The logarithmic derivative is computed via finite differences on ln M with step δ=0.01.

Parameters:

• m_h (jnp.ndarray Halo masses [M_sun/h].)

• z (float Redshift.)

• theta (dict Cosmological parameters.)

n_eff(m_min: float, m_max: float, z: float, theta: dict) → Array

Effective number density n(M > m_min) integrated to m_max [h^3 Mpc^{-3}].

Parameters:

• m_min, m_max (float Mass limits [M_sun/h].)

• z (float Redshift.)

• theta (dict Cosmological parameters.)

sigma(m_h: Array, z: float, theta: dict) → Array

RMS linear density fluctuation σ(M, z) [dimensionless].

Computed at z=0 from the linear power spectrum, optionally rescaled to match theta['sigma8'] (needed when pk_func returns shape-only spectra such as eisenstein_hu_pk), then multiplied by the growth factor D(z)/D(0):

\[\sigma(M, z) = \sigma(M, 0) \times D(z) / D(0)\]

Parameters:

• m_h (jnp.ndarray Halo masses [M_sun/h].)

• z (float Redshift.)

• theta (dict Cosmological parameters. Omega_m required.) – sigma8 optional — triggers amplitude normalisation.

hod_mod.core.halo_mass_function.delta_vir_flat_jax(z: float, omega_m) → Array

Virial overdensity w.r.t. critical density (Bryan & Norman 1998).

\(\Delta_{\rm vir}(z) = 18\pi^2 + 82 x - 39 x^2\) where \(x = \Omega_m(z) - 1\) (valid for flat \(\Lambda\mathrm{CDM}\)).

Parameters:

• z (float Redshift (Python float).)

• omega_m (float or jnp.ndarray Matter density parameter Ω_m at z=0.)

• Accuracy

• ——–

• Δ_vir(z=0, Ω_m=1) = 18π² ≈ 177.7 (EdS limit; x=0 by definition); verified

• to < 2% (2026-04-23). For Planck 2018 at z=0 (Δ_vir ≈ 358.)

• Timing

• ——

• ~ 3 µs / call (scalar input, CPU x86-64, 2026-04-23).

hod_mod.core.halo_mass_function.fsigma_angulo12(sigma: Array, z: float = 0.0) → Array

Angulo et al. 2012 multiplicity function.

\[f(\sigma) = 0.201 \left[(2.08/\sigma)^{1.7} + 1\right] \exp(-1.172/\sigma^2)\]

Parameters:

• sigma (jnp.ndarray σ(M, z).)

• z (float Redshift (unused).)

hod_mod.core.halo_mass_function.fsigma_bhattacharya11(sigma: Array, z: float = 0.0) → Array

Bhattacharya et al. 2011 multiplicity function (z-dependent).

\[f(\sigma) = A(z) \sqrt{\frac{2}{\pi}} \exp\!\left(-\frac{a(z)\nu^2}{2}\right) \left(1 + (a(z)\nu^2)^{-p}\right) (\nu\sqrt{a(z)})^q\]

with ν=δ_c/σ, A=0.333(1+z)^{-0.11}, a=0.788(1+z)^{-0.01}, p=0.807, q=1.795.

Parameters:

• sigma (jnp.ndarray σ(M, z).)

• z (float Redshift.)

hod_mod.core.halo_mass_function.fsigma_bocquet16(sigma: Array, z: float = 0.0, hydro: bool = False) → Array

Bocquet et al. 2016 multiplicity function calibrated for Δ=200m.

Power-law z-evolution of the Tinker-type parameters (Equations 6–8 and Table 2 of Bocquet+2016 MNRAS 456, 2361). Separate fits for DM-only and hydro simulations are selected via the hydro flag.

\[f(\sigma) = A(z)\left[(\sigma/b(z))^{-a(z)} + 1\right] \exp\!\left(-c(z)/\sigma^2\right)\]

Parameters:

• sigma (jnp.ndarray σ(M, z).)

• z (float Redshift.)

• hydro (bool If True use hydro-simulation fit; else DM-only (default).)

hod_mod.core.halo_mass_function.fsigma_comparat17(sigma: Array, z: float = 0.0) → Array

Comparat et al. 2017 multiplicity function (virial mass, z=0).

Bhattacharya-type fit to the MultiDark-Planck simulation at z=0. The parameters are updated from the published version. Calibrated for the virial SO mass definition (Comparat+2017 MNRAS 469, 4157).

\[f(\sigma) = A \sqrt{\frac{2}{\pi}} \exp\!\left(-\frac{a\nu^2}{2}\right) \left(1 + (a\nu^2)^{-p}\right) (\nu\sqrt{a})^q, \quad \nu = \delta_c / \sigma\]

with A=0.324, a=0.897, p=0.624, q=1.589, δ_c=1.686.

Parameters:

• sigma (jnp.ndarray σ(M, z).)

• z (float Redshift (unused — calibrated at z=0 only).)

hod_mod.core.halo_mass_function.fsigma_courtin11(sigma: Array, z: float = 0.0) → Array

Courtin et al. 2011 multiplicity function.

Sheth-Tormen-type fit with \(\delta_c = 1.673\):

\[f(\sigma) = A \sqrt{\frac{2a}{\pi}} \frac{\nu'}{\sigma} \exp\!\left(-\frac{\nu'^2}{2}\right) \left(1 + \nu'^{-2p}\right)\]

with A=0.348, a=0.695, p=0.1 (FoF).

Parameters:

• sigma (jnp.ndarray σ(M, z).)

• z (float Redshift (unused).)

hod_mod.core.halo_mass_function.fsigma_crocce10(sigma: Array, z: float = 0.0) → Array

Crocce et al. 2010 multiplicity function (z-dependent).

\[f(\sigma) = A(z) (\sigma^{-a(z)} + b(z)) \exp(-c(z)/\sigma^2)\]

with coefficients evolving as power laws in (1+z).

Parameters:

• sigma (jnp.ndarray σ(M, z).)

• z (float Redshift.)

hod_mod.core.halo_mass_function.fsigma_despali16(sigma: Array, z: float = 0.0, delta_ratio: float = 1.0) → Array

Despali et al. 2016 multiplicity function for arbitrary SO mass definition.

Parameterises the Sheth-Tormen form using a polynomial in \(x = \log_{10}(\Delta / \Delta_{\rm vir})\) (Equation 12 of Despali+2016 MNRAS 456, 2486):

\[A(x) = -0.1362 x + 0.3292, \quad a(x) = 0.4332 x^2 + 0.2263 x + 0.7665, \quad p(x) = -0.1151 x^2 + 0.2554 x + 0.2488\]

with \(f(\sigma) = 2A \sqrt{\nu'/(2\pi)} e^{-\nu'/2} (1 + \nu'^{-p})\), \(\nu' = a\delta_c^2/\sigma^2\).

Parameters:

• sigma (jnp.ndarray σ(M, z).)

• z (float Redshift (unused here — mass-def conversion handles z-dep.).)

• delta_ratio (float) – \(\Delta_{\rm target} / \Delta_{\rm vir}\) in critical density units. Pass 1.0 (default) for the virial mass definition. Use delta_vir_flat_jax() to compute \(\Delta_{\rm vir}\).

hod_mod.core.halo_mass_function.fsigma_jenkins01(sigma: Array, z: float = 0.0) → Array

Jenkins et al. 2001 multiplicity function.

\[f(\sigma) = 0.315 \exp\!\left(-|\ln\sigma^{-1} + 0.61|^{3.8}\right)\]

Parameters:

• sigma (jnp.ndarray σ(M, z).)

• z (float Redshift (unused).)

hod_mod.core.halo_mass_function.fsigma_press74(sigma: Array, z: float = 0.0) → Array

Press & Schechter 1974 multiplicity function.

\[f(\sigma) = \sqrt{\frac{2}{\pi}} \nu \exp\!\left(-\frac{\nu^2}{2}\right), \quad \nu = \delta_c / \sigma\]

Parameters:

• sigma (jnp.ndarray RMS linear density fluctuation σ(M, z).)

• z (float Redshift (unused, kept for uniform interface).)

• Accuracy

• ——–

• ∫ f(σ) d ln σ⁻¹ = 1 (cloud-in-cloud normalisation) verified to < 0.5%

• over σ ∈ [0.1, 5] (numerical integration, 2026-04-23).

• Timing

• ——

• ~ 85 µs / call (JIT-compiled, N=200 σ values, CPU x86-64, 2026-04-23).

hod_mod.core.halo_mass_function.fsigma_rodriguezpuebla16(sigma: Array, z: float = 0.0) → Array

Rodriguez-Puebla et al. 2016 multiplicity function (virial mass, z-dep.).

Polynomial z-evolution of Tinker-type parameters calibrated for the Planck cosmology and virial SO mass definition (Table 2 of Rodriguez-Puebla+2016 MNRAS 462, 893).

Parameters:

• sigma (jnp.ndarray σ(M, z).)

• z (float Redshift (calibrated 0 ≤ z ≤ 7).)

hod_mod.core.halo_mass_function.fsigma_seppi20(sigma: Array, z: float = 0.0) → Array

Seppi et al. 2020 multiplicity function marginalized over xoff and spin.

The full model is a 3D distribution over (σ, x_off, λ) (Equation 21 of Seppi+2020 A&A 643, A17). This function returns the 1-D marginal \(f(\sigma)\) obtained by integrating over \(\log_{10}(x_{\rm off})\) and \(\log_{10}(\lambda)\).

Calibrated for M > 4×10¹³ M☉/h and the virial SO mass definition.

Parameters:

• sigma (jnp.ndarray σ(M, z).)

• z (float Redshift.)

hod_mod.core.halo_mass_function.fsigma_sheth99(sigma: Array, z: float = 0.0) → Array

Sheth & Tormen 1999 multiplicity function.

\[f(\sigma) = A \sqrt{\frac{2a}{\pi}} \nu' \exp\!\left(-\frac{a\nu'^2}{2}\right) \left(1 + (a\nu'^2)^{-p}\right), \quad \nu' = \delta_c / \sigma\]

with A=0.3222, a=0.707, p=0.3.

Parameters:

• sigma (jnp.ndarray σ(M, z).)

• z (float Redshift (unused).)

• Accuracy

• ——–

• Reproduces Sheth & Tormen 1999 Fig. 2 to < 1% for σ ∈ [0.2, 3].

• Timing

• ——

• ~ 15 µs / call (JIT-compiled, N=200 σ values, CPU x86-64, 2026-04-23).

hod_mod.core.halo_mass_function.fsigma_tinker08(sigma: Array, z: float = 0.0, Delta: float = 200.0) → Array

Tinker et al. 2008 multiplicity function with z-evolution and Δ interpolation.

Equations 2–5 and Table 2 of Tinker+2008 (ApJ 688, 709):

\[f(\sigma) = A(z)\left[\left(\frac{\sigma}{b(z)}\right)^{-a(z)} + 1\right] \exp\!\left(-\frac{c}{\sigma^2}\right)\]

where the parameters at z=0 are interpolated from Table 2 as a function of overdensity Δ (w.r.t. mean), and evolve with redshift as:

\[A(z) = A_0 (1+z)^{-0.14}, \quad a(z) = a_0 (1+z)^{-0.06}, \quad b(z) = b_0 (1+z)^{-\alpha}, \quad \alpha = 10^{-(0.75/\log_{10}(\Delta/75))^{1.2}}\]

Parameters:

• sigma (jnp.ndarray σ(M, z).)

• z (float Redshift.)

• Delta (float Overdensity with respect to mean matter density (default 200).)

• Accuracy

• ——–

• f(σ=1, z=0, Δ=200) ≈ 0.283 (from Table 2 params (A₀=0.186, a₀=1.47, b₀=2.57,)

• c₀=1.19); agrees with Tinker+2008 to < 5% for σ ∈ [0.3, 2] (2026-04-23).

• Timing

• ——

• ~ 116 µs / call (JIT-compiled, N=200 σ values, CPU x86-64, 2026-04-23).

hod_mod.core.halo_mass_function.fsigma_warren06(sigma: Array, z: float = 0.0) → Array

Warren et al. 2006 multiplicity function.

\[f(\sigma) = A (\sigma^{-a} + b) \exp(-c/\sigma^2)\]

with A=0.7234, a=1.625, b=0.2538, c=1.1982.

Parameters:

• sigma (jnp.ndarray σ(M, z).)

• z (float Redshift (unused).)

• Accuracy

• ——–

• Reproduces Warren et al. 2006 Table 3 values to < 2% for σ ∈ [0.3, 2].

• Timing

• ——

• ~ 18 µs / call (JIT-compiled, N=200 σ values, CPU x86-64, 2026-04-23).

hod_mod.core.halo_mass_function.fsigma_watson13(sigma: Array, z: float = 0.0) → Array

Watson et al. 2013 multiplicity function (Friends-of-Friends).

Parameters for FoF linking length b=0.2: A=0.282, a=2.163, b=1.406, c=1.210, γ=1.082.

Parameters:

• sigma (jnp.ndarray σ(M, z).)

• z (float Redshift (unused — FoF calibration).)

hod_mod.core.halo_mass_function.fsigma_yung24(sigma: Array, z: float = 0.0) → Array

Yung et al. 2024 multiplicity function (virial mass, z-dep.).

Calibrated from the GUREFT simulations at high redshift. Uses the Tinker functional form with polynomial z-dependence (Table 2 of Yung+2024 MNRAS 530, 4868).

Parameters:

• sigma (jnp.ndarray σ(M, z).)

• z (float Redshift (calibrated 0 ≤ z ≤ 20).)

hod_mod.core.halo_mass_function.fsigma_yung25(sigma: Array, z: float = 0.0) → Array

Yung et al. 2025 multiplicity function (virial mass, high-z calibrated).

Calibrated from the GUREFT simulations. Unlike yung24, this fit is optimised for z > 6 and not recommended at low redshift (Yung+2025 MNRAS 543, 3802).

Parameters:

• sigma (jnp.ndarray σ(M, z).)

• z (float Redshift (calibrated 6 ≤ z ≤ 30).)

hod_mod.core.halo_mass_function.growth_factor(z, theta: dict)

D(z)/D(0) dispatcher: CPL ODE when theta carries beyond-ΛCDM keys.

The branch is on dictionary membership (static under tracing): forecast parameter dicts include w0/wa and get the differentiable ODE growth; production ΛCDM dicts keep the Carroll+1992 formula bit-identical.

hod_mod.core.halo_mass_function.make_hmf(backend: str = 'tinker08', pk_func=None, rho_mean: float = 86010827748.85518, Delta: float = 200.0, **fsigma_kwargs)

Return a HaloMassFunction (or emulator wrapper) for the requested backend.

All backends expose: .dndm(), .bias(), .sigma(), .n_eff().

Parameters:

• backend (str) – Analytic multiplicity model (any key in _FSIGMA_MODELS, e.g. tinker08, bocquet16, yung25) or an emulator backend:

  • "csst" — CSST CEmulator HMF (Chen+2025, SCPMA 2025).

  • "aemulusnu" — Aemulus-ν HMF (Shen+2025, JCAP 2025, arXiv:2410.00913). Valid for M ≥ 10^13 M_sun/h, z ≤ 2; masses below that are silently extrapolated (see AemulusNuHaloMassFunction), so this backend is not recommended for low stellar-mass-threshold HOD samples whose occupation extends well below 10^13 M_sun/h.

• pk_func (callable) – (k, z, theta) → P_lin(k). Required for analytic backends and for bias() / sigma() when using emulator backends. Ignored for "csst" σ (uses CSST PkLin emulator internally).

• rho_mean (float) – Mean comoving matter density at z=0 [M_sun/h / (Mpc/h)^3].

• Delta (float) – Overdensity threshold w.r.t. mean density (tinker08 and bias).

• **fsigma_kwargs – Forwarded to the multiplicity function for analytic backends.

Examples

>>> hmf = make_hmf("tinker08", pk_func=my_pk)
>>> hmf = make_hmf("bocquet16", pk_func=my_pk, hydro=True)
>>> hmf = make_hmf("csst")
>>> hmf = make_hmf("aemulusnu", pk_func=my_pk)

hod_mod.core.halo_mass_function.tinker08_fsigma(sigma: Array, z: float = 0.0, Delta: float = 200.0) → Array

Tinker et al. 2008 multiplicity function with z-evolution and Δ interpolation.

Equations 2–5 and Table 2 of Tinker+2008 (ApJ 688, 709):

\[f(\sigma) = A(z)\left[\left(\frac{\sigma}{b(z)}\right)^{-a(z)} + 1\right] \exp\!\left(-\frac{c}{\sigma^2}\right)\]

where the parameters at z=0 are interpolated from Table 2 as a function of overdensity Δ (w.r.t. mean), and evolve with redshift as:

\[A(z) = A_0 (1+z)^{-0.14}, \quad a(z) = a_0 (1+z)^{-0.06}, \quad b(z) = b_0 (1+z)^{-\alpha}, \quad \alpha = 10^{-(0.75/\log_{10}(\Delta/75))^{1.2}}\]

Parameters:

• sigma (jnp.ndarray σ(M, z).)

• z (float Redshift.)

• Delta (float Overdensity with respect to mean matter density (default 200).)

• Accuracy

• ——–

• f(σ=1, z=0, Δ=200) ≈ 0.283 (from Table 2 params (A₀=0.186, a₀=1.47, b₀=2.57,)

• c₀=1.19); agrees with Tinker+2008 to < 5% for σ ∈ [0.3, 2] (2026-04-23).

• Timing

• ——

• ~ 116 µs / call (JIT-compiled, N=200 σ values, CPU x86-64, 2026-04-23).

hod_mod.core.halo_mass_function.tinker10_bias(nu: Array, Delta: float = 200.0) → Array

Tinker et al. 2010 large-scale halo bias b(ν).

Equation 6 of Tinker+2010 (ApJ 724, 878):

\[b(\nu) = 1 - A \frac{\nu^a}{\nu^a + \delta_c^a} + B \nu^b + C \nu^c\]

Parameters are from Table 2 evaluated at \(\Delta=200\) (mean density): A, a depend on Δ; B=0.183, b=1.5, C, c depend on Δ.

Parameters:

• nu (jnp.ndarray Peak height ν = δ_c / σ(M, z).)

• Delta (float Overdensity w.r.t. mean density (default 200).)

• Accuracy

• ——–

• b(ν=1) ≈ 0.9–1.1 (near-unity bias at characteristic mass M_*); verified

• to < 15% vs Tinker+2010 Table 2 for Δ=200, ν ∈ [0.5, 3] (2026-04-23).

• Timing

• ——

• ~ 26 µs / call (JIT-compiled, N=200 ν values, CPU x86-64, 2026-04-23).

—

Halo Profiles

(hod_mod.core.halo_profiles)

NFW profile

Navarro, Frenk & White 1997 [NFW1997] showed that the radial density profiles of dark matter halos are well described by

\[\rho_{\rm NFW}(r) = \frac{\rho_s}{(r/r_s)(1 + r/r_s)^2}\]

where \(r_s = r_{200}/c\) is the scale radius and \(c\) is the concentration parameter. The characteristic density is

\[\rho_s = \frac{M_{200}}{4\pi r_s^3 \left[\ln(1+c) - c/(1+c)\right]}.\]

Einasto profile

Einasto 1965 [Einasto1965] proposed a power-law logarithmic slope profile:

\[\rho_{\rm Ein}(r) = \rho_{-2} \exp\left\{ -\frac{2}{\alpha_E}\left[\left(\frac{r}{r_{-2}}\right)^{\alpha_E} - 1\right] \right\}\]

where \(r_{-2}\) is the radius where the logarithmic slope equals \(-2\) and \(\alpha_E \approx 0.18\) for typical halos.

Concentration–mass relation

The concentration is obtained from colossus using the Diemer & Joyce 2019 relation (Diemer & Joyce 2019) [DiemerJoyce2019]:

\[c(M, z) = c_0 \left(\frac{M}{M_{\rm piv}}\right)^{-\kappa_c} \left(1 + z\right)^{-\mu_c}\]

Projected quantities

The surface mass density is the line-of-sight projection:

\[\Sigma(R) = 2 \int_0^\infty \rho\!\left(\sqrt{R^2 + \ell^2}\right) d\ell\]

The mean surface density within radius \(R\) is

\[\bar{\Sigma}(<R) = \frac{2}{R^2} \int_0^R R'\,\Sigma(R')\,dR'\]

and the excess surface density (the weak-lensing observable) is

\[\Delta\Sigma(R) = \bar{\Sigma}(<R) - \Sigma(R).\]

Fourier transform of the NFW profile (for the halo model):

\[u(k|M) = \frac{4\pi\rho_s r_s^3}{M} \left[\cos(kr_s)\left({\rm Ci}((1+c)kr_s) - {\rm Ci}(kr_s)\right) + \sin(kr_s)\left({\rm Si}((1+c)kr_s) - {\rm Si}(kr_s)\right) - \frac{\sin(ckr_s)}{(1+c)kr_s}\right]\]

where Ci and Si are the cosine and sine integrals.

NFW and Einasto halo profiles plus Fourier-space window functions.

Provides 3D density, projected surface density, lensing ΔΣ (all in JAX), the NFW normalized Fourier transform needed for the full halo model (Cooray & Sheth 2002), and the Einasto (1965) alternative profile (Asgari+2023 Eq. 47).

References

Bartelmann 1996; Wright & Brainerd 2000 — NFW projected Σ and ΔΣ Cooray & Sheth 2002, Phys.Rep. 372, 1 — NFW Fourier transform (Eq. 11) Einasto 1965; Asgari+2023 arXiv:2303.08752 Eq. 47 — Einasto profile

class hod_mod.core.halo_profiles.HaloProfile(cosmo_params: dict, cm_relation: str = 'diemer19', mdef: str = '200m')

Bases: object

Concentration–mass relation and NFW profile parameters.

Supports two backends:

• cm_relation='dutton14' — JAX-native Dutton & Macciò 2014 power-law (requires mdef='200c'). Fully differentiable w.r.t. halo mass.

• Any colossus key (e.g. 'diemer19') — wraps colossus; not autodiff-capable but supports all mass definitions and models.

Parameters:

• cosmo_params (dict) – Colossus-style cosmological parameters (ignored for cm_relation='dutton14').

• cm_relation (str) – 'dutton14' for the JAX-native backend, or any colossus model name.

• mdef (str) – Mass definition, e.g. '200m' or '200c'. Must be '200c' when cm_relation='dutton14'.

concentration(m_h: Array, z: float) → Array

Concentration parameter c(M, z) from the chosen c-M relation.

delta_sigma(R_proj: Array, m_h: Array, z: float, theta_cosmo: dict) → Array

ΔΣ(R) [M_sun h / Mpc^2] for a single halo of mass m_h.

rho_s_and_rs(m_h: Array, z: float, theta_cosmo: dict) → tuple[Array, Array]

Characteristic density ρ_s and scale radius r_s [Mpc/h] for NFW.

r_delta = (3 M / 4π delta rho_ref)^{1/3} with (delta, rho_ref) from the mass definition mdef set at construction time. c = r_delta / r_s.

Parameters:

• m_h (jnp.ndarray — halo mass [Msun/h])

• z (float — redshift)

• theta_cosmo (dict — cosmological parameters (needs Omega_m))

hod_mod.core.halo_profiles.concentration_dutton14_jax(m_h: Array, z: float) → Array

Concentration \(c_{200c}(M, z)\) from Dutton & Macciò 2014 (MNRAS 441, 3359).

\[ \begin{align}\begin{aligned}\log_{10}(c_{200c}) = a(z) + b(z)\, \log_{10}\!\left(\frac{M_{200c}}{10^{12}\,h^{-1}M_\odot}\right)\\a(z) = 0.520 + 0.385\,\exp(-0.617\,z^{1.21})\\b(z) = -0.101 + 0.026\,z\end{aligned}\end{align} \]

Valid for \(M_{200c} \in [10^{10}, 10^{15}]\,h^{-1}M_\odot\) and \(z \in [0, 5]\). Use with HaloProfile(mdef='200c', cm_relation='dutton14'). Fully differentiable w.r.t. m_h.

Parameters:

• m_h (jnp.ndarray — halo mass \(M_{200c}\) [M_sun/h])

• z (float — redshift (static; JIT-specialised per redshift value))

Returns:

c (jnp.ndarray — concentration \(c_{200c}\), same shape as m_h)

hod_mod.core.halo_profiles.einasto_rho(r: Array, rho_s: float, r_s: float, alpha: float = 0.18) → Array

Einasto (1965) density profile ρ(r) [M_sun h² / Mpc³].

\[\rho(r) = \rho_s \exp\!\left[-\frac{2}{\alpha} \left(\left(\frac{r}{r_s}\right)^\alpha - 1\right)\right]\]

(Asgari+2023 Eq. 47; Einasto 1965)

α ≈ 0.18 gives a profile close to NFW for cluster-mass halos (Klypin+2001, Merritt+2006). Smaller α → steeper inner cusp.

Parameters:

• r ([Mpc/h], shape (Nr,))

• rho_s (characteristic density [M_sun h² / Mpc³])

• r_s (scale radius [Mpc/h]; ρ(r_s) = ρ_s exp(0) = ρ_s)

• alpha (shape parameter (default 0.18))

Returns:

• rho ([M_sun h² / Mpc³], shape (Nr,))

• Accuracy

• ——–

• ρ(r_s) = ρ_s exactly (by construction; exp argument = 0 at r = r_s).

• Monotonically decreasing verified analytically; numerical normalisation

• ∫ 4πr² ρ dr (N=2000 log nodes) matches einasto_uk (k→0) to < 2%

• for c ∈ [5, 20] (2026-04-23).

• Timing

• ——

• ~ 21 µs / call (JIT-compiled, N=100 radii, CPU x86-64, 2026-04-23).

hod_mod.core.halo_profiles.einasto_uk(k_arr: ndarray, r_s_arr: ndarray, c_arr: ndarray, alpha: float = 0.18, n_r: int = 200) → Array

Einasto normalized Fourier transform û_m(k, M) via Gauss-Legendre quadrature.

\[\hat{u}_m(k|M) = \frac{ \int_0^{r_h} \rho_{\rm Ein}(r)\,j_0(kr)\,r^2\,\mathrm{d}r }{ \int_0^{r_h} \rho_{\rm Ein}(r)\,r^2\,\mathrm{d}r }\]

where \(r_h = c\,r_s\) is the truncation radius and

\[\rho_{\rm Ein}(r) = \rho_s\exp\!\left[ -\frac{2}{\alpha}\left(\left(\frac{r}{r_s}\right)^\alpha - 1\right) \right]\]

The ratio is independent of \(\rho_s\) and satisfies \(\hat{u}_m(k\to 0) = 1\). Integrals are evaluated by n_r-point Gauss-Legendre quadrature on \([0, c]\).

Parameters:

• k_arr (array_like, shape (Nk,), wavenumbers [h/Mpc])

• r_s_arr (array_like, shape (NM,), Einasto scale radii [Mpc/h])

• c_arr (array_like, shape (NM,), concentration c = r_h / r_s)

• alpha (float) – Einasto shape parameter (default 0.18, close to NFW for clusters).

• n_r (int) – Number of Gauss-Legendre quadrature nodes (default 200).

Returns:

• uk (jnp.ndarray, shape (Nk, NM), dimensionless, in (0, 1])

• Accuracy

• ——–

• k→0 limit û→1 verified to < 1% (n_r=200 nodes, α=0.18). Converges to

• < 0.1% relative error vs n_r=1000 benchmark for k ∈ [0.01, 100] h/Mpc

• (2026-04-23).

• Timing

• ——

• ~ 22 ms / call (not JIT-compiled, Nk=50 × NM=10, n_r=200, CPU x86-64,

• 2026-04-23).

hod_mod.core.halo_profiles.nfw_delta_sigma(R: Array, rho_s: float, r_s: float) → Array

NFW excess surface density ΔΣ(R) = Σ_bar(<R) − Σ(R) [M_sun h / Mpc^2].

This is the galaxy-galaxy lensing observable.

hod_mod.core.halo_profiles.nfw_mass(r: Array, rho_s: float, r_s: float) → Array

NFW enclosed mass M(<r) [M_sun/h].

hod_mod.core.halo_profiles.nfw_mean_sigma(R: Array, rho_s: float, r_s: float) → Array

Mean projected surface density Σ_bar(<R) inside radius R (analytic).

Σ_bar(<R) = (2/R²) ∫₀^R Σ(R’) R’ dR’ Uses Wright & Brainerd 2000 Eq. 13.

hod_mod.core.halo_profiles.nfw_rho(r: Array, rho_s: float, r_s: float) → Array

NFW 3D density profile ρ(r) [M_sun h^2 / Mpc^3].

ρ(r) = ρ_s / [(r/r_s)(1 + r/r_s)²]

hod_mod.core.halo_profiles.nfw_sigma(R: Array, rho_s: float, r_s: float) → Array

Projected NFW surface density Σ(R) [M_sun h / Mpc^2] (analytic).

Uses the Bartelmann 1996 / Wright & Brainerd 2000 closed form.

hod_mod.core.halo_profiles.nfw_uk(k_arr: ndarray, r_s_arr: ndarray, c_arr: ndarray) → Array

NFW normalized Fourier transform û_m(k, M) (Cooray & Sheth 2002, Eq. 11).

\[\hat{u}_m(k|M) = \frac{1}{M}\int_0^{r_h} \rho_{\rm NFW}(r)\,j_0(kr)\,4\pi r^2\,dr\]

The analytic result for a truncated NFW profile (truncation at r_h = c r_s):

\[\hat{u}_m = \frac{ \cos(K)[{\rm Ci}(K(1+c)) - {\rm Ci}(K)] + \sin(K)[{\rm Si}(K(1+c)) - {\rm Si}(K)] - \sin(cK) / [(1+c)K] }{\ln(1+c) - c/(1+c)},\quad K = k\,r_s\]

(derivation: IBP on ∫₀^c sin(Kx)/(1+x)² dx, substitute t = K(1+x))

û_m(k → 0) = 1 by l’Hôpital (verified analytically). Not JIT-compatible: uses scipy.special.sici.

Parameters:

• k_arr (array_like, shape (Nk,), wavenumbers [h/Mpc])

• r_s_arr (array_like, shape (NM,), NFW scale radii [Mpc/h])

• c_arr (array_like, shape (NM,), concentration c = r_h / r_s)

Returns:

• uk (jnp.ndarray, shape (Nk, NM), dimensionless, in (0, 1])

• Accuracy

• ——–

• k→0 limit û→1 verified to < 1% for K < 1e-6 (L’Hôpital guard applied).

• Shape agrees with direct numerical quadrature (200 nodes) to < 0.1% for

• k ∈ [0.01, 100] h/Mpc, c = 10, r_s = 0.3 Mpc/h (2026-04-23).

• Timing

• ——

• ~ 196 µs / call (not JIT-compiled, Nk=50 × NM=10, CPU x86-64, 2026-04-23).

hod_mod.core.halo_profiles.nfw_uk_jax(k_arr: Array, r_s_arr: Array, c_arr: Array) → Array

NFW normalized Fourier transform û_m(k, M), JAX-native (autodiff-compatible).

Same analytic formula as nfw_uk() (Cooray & Sheth 2002, Eq. 11) but replaces scipy.special.sici with a pure-JAX series/asymptotic implementation via _si_jax() / _ci_jax(). Fully JIT-compiled and differentiable w.r.t. r_s_arr and c_arr.

See nfw_uk() for the analytic formula and accuracy notes.

Parameters:

• k_arr (jnp.ndarray, shape (Nk,))

• r_s_arr (jnp.ndarray, shape (NM,))

• c_arr (jnp.ndarray, shape (NM,))

Returns:

• uk (jnp.ndarray, shape (Nk, NM), in (0, 1])

• Accuracy

• ——–

• Agrees with scipy-based nfw_uk to < 0.1% for K ∈ [10⁻⁴, 100] h/Mpc,

• c ∈ [3, 20], r_s ∈ [0.01, 5] Mpc/h (verified 2026-05-19).

hod_mod.core.halo_profiles.satellite_nfw_uk(k_arr: ndarray, r_s_arr: ndarray, c_arr: ndarray, r_vir_arr: ndarray, b_sat_conc: float = 1.0, f_cut: float = 0.0, gamma: float = 0.0, n_r: int = 100, n_k_coarse: int = 128) → Array

Satellite normalized FT combining three inner-profile extensions (GL quadrature).

The satellite number density profile:

\[n_{\rm sat}(r) \propto \left(\frac{r}{r_{\rm vir}}\right)^{\gamma} \left[1 - \exp\!\left(-\frac{r}{f_{\rm cut}\,r_{\rm vir}}\right)\right] \rho_{\rm NFW}(r;\,c_{\rm sat}), \quad 0 \le r \le r_{\rm vir}\]

with \(c_{\rm sat} = b_{\rm sat\_conc}\,c_{\rm DM}\).

—

Gas Profiles

(hod_mod.gas)

Two parametric halo gas profiles are provided for computing galaxy × tSZ (thermal Sunyaev-Zel’dovich Compton-\(y\)) and galaxy × soft X-ray cross-correlations within the halo model.

M₂₀₀ → M_500c conversion

Both the A10 pressure profile and the DPM density profile use overdensity \(\Delta = 500c\). The static halo model cache stores halo masses and radii at \(\Delta = 200m\). The helper function m200_to_m500c() performs the conversion analytically using the NFW enclosed-mass formula (Navarro, Frenk & White 1997):

\[M_{\rm NFW}(r | M_{200}, c_{200}) = 4\pi\rho_s r_s^3 \left[\ln\!\left(1 + \frac{r}{r_s}\right) - \frac{r/r_s}{1 + r/r_s}\right]\]

A bisection (scipy.optimize.brentq) finds \(r_{500c}\) such that \(M_{\rm NFW}(r_{500c}) = (4\pi/3)\,500\,\rho_{\rm crit}(z)\,r_{500c}^3\) and returns \((M_{500c},\,R_{500c})\).

Arnaud+2010 Pressure Profile (A10) — for tSZ

Reference: Arnaud et al. 2010, A&A 517, A92, Table 1.

The generalised NFW (gNFW) shape function:

\[p(x) = \frac{P_0}{\left(c_{500}\,x\right)^\gamma \left[1 + \left(c_{500}\,x\right)^\alpha\right]^{(\beta-\gamma)/\alpha}}\]

where \(x = r/R_{500c}\). Universal parameters (A10 Table 1): \(P_0 = 8.403\), \(c_{500} = 1.177\), \(\gamma = 0.3081\), \(\alpha = 1.0510\), \(\beta = 5.4905\), \(\alpha_p = 0.12\).

Physical electron pressure (A10, Eq. 11):

\[P_e(r|M_{500c},z) = 1.65\times10^{-3}\,h(z)^{8/3} \left[\frac{M_{500c}}{3\times10^{14}\,h_{70}^{-1}M_\odot}\right]^{2/3+\alpha_p} p\!\left(\frac{r}{R_{500c}}\right) \quad [h_{70}^2\,{\rm keV\,cm}^{-3}]\]

where \(h(z) = H(z)/H_0\).

The Fourier transform of the y-profile per halo:

\[\tilde{y}(k|M,z) = \frac{\sigma_T}{m_e c^2} \int_0^{r_{\rm max}} P_e(r|M,z)\,\frac{\sin(kr)}{kr}\,4\pi r^2\,dr \quad [({\rm Mpc}/h)^2]\]

where \(\sigma_T = 6.6524\times10^{-25}\,{\rm cm}^2\) and \(m_e c^2 = 511\,{\rm keV}\). The integral is computed via Gauss-Legendre quadrature with \(r_{\rm max} = 5\,R_{500c}\) and 200 nodes by default.

DPM Electron Density Profile — for soft X-ray

Reference: Oppenheimer et al. 2025, arXiv:2505.14782.

The gNFW-shaped density profile:

\[f(x|\boldsymbol{\alpha}) = x^{-\alpha_{\rm in}} \left(1 + x^{\alpha_{\rm tr}}\right)^{(\alpha_{\rm in}-\alpha_{\rm out})/\alpha_{\rm tr}}\]

where \(x = r/R_s\) with scale radius \(R_s = R_{200}/c_{\rm DPM}\), \(c_{\rm DPM} = 2.772\).

Electron density:

\[n_e(r|M_{200},z) = n_{e,03}\left(\frac{M_{200}}{10^{12}M_\odot}\right)^\beta E(z)^\gamma\,f\!\left(\frac{r}{R_s}\right)\]

where \(n_{e,03}\) is the normalisation at \(r = 0.3\,R_{200}\) for \(M_{200}=10^{12}\,M_\odot/h\) at \(z=0\). Three calibrated DPM variants are provided (model=1, 2, 3; see Oppenheimer et al. 2025, Table 2).

The X-ray emissivity Fourier transform per halo is

\[\tilde{\varepsilon}(k|M,z) = \int_0^{r_{\rm max}} n_e^2(r|M,z)\,\frac{\sin(kr)}{kr}\,4\pi r^2\,dr \quad [({\rm Mpc}/h)^3\,{\rm cm}^{-6}]\]

with \(r_{\rm max} = 3\,R_{200}\) and 200 GL nodes.

Hot-gas / ICM fields shared by the X-ray and thermal-SZ cross-correlations.

Pressure profiles (Arnaud+2010, DPM) feed the tSZ Compton-y signal; the gas density, cooling/emissivity (APEC) and metallicity profiles feed the soft X-ray emissivity. The eROSITA instrument response converts intrinsic emission to observed count rates. These objects are consumed by hod_mod.observables.cross_spectra.HaloModelCrossSpectra.

class hod_mod.gas.ApecCoolingTable(emin: float = 0.5, emax: float = 2.0, n_T: int = 60, T_min: float = 0.08, T_max: float = 20.0, n_Z: int = 15, Z_min: float = 0.05, Z_max: float = 3.0, apec_vers: str = None, nbins: int = 1000)

Bases: object

Band-integrated APEC cooling function Λ(T, Z) precomputed as a 2D table.

Uses soxs.ApecGenerator to compute the X-ray emission spectrum from the AtomDB APEC CIE plasma model for a grid of temperatures and metallicities, integrates each spectrum over the requested energy band, and stores the result as a 2D log-log interpolator for fast evaluation.

The output follows the APEC normalization convention:

\[\varepsilon(r) = n_e(r)\,n_H(r)\,\Lambda(T(r), Z(r))\]

with \(n_H \approx 0.83\,n_e\) (fully ionized solar-abundance plasma). The stored table gives \(\Lambda_{n_e^2}(T, Z) = 0.83\,\Lambda_{\rm APEC}\), so that \(\varepsilon = n_e^2\,\Lambda_{n_e^2}\).

Parameters:

• emin, emax (float) – Energy band edges [keV]. Default: 0.5–2.0 (soft X-ray).

• n_T (int) – Number of temperature grid points (log-spaced). Default: 60.

• T_min, T_max (float) – Temperature range [keV]. Default: 0.08–20.

• n_Z (int) – Number of metallicity grid points (log-spaced). Default: 15.

• Z_min, Z_max (float) – Metallicity range [Z_sun]. Default: 0.05–3.0.

• apec_vers (str) – APEC version string (default from soxs config, currently “3.1.3”).

• nbins (int) – Number of spectral bins used for the band integration. Default: 1000.

Notes

Requires soxs (pip install pyxsim soxs) and the APEC spectral tables, which are downloaded once via soxs.download_spectrum_tables("apec").

Initialisation takes a few seconds (60×15 = 900 APEC evaluations). Cache the instance across calls.

class hod_mod.gas.ErositaResponse(response_npz: str | None = None)

Bases: object

Combined TM1-7 eROSITA DR1 response → energy-conversion factors.

Parameters:

response_npz (str | None) – Path to the distilled response artifact (ARF + in-band RMF efficiency).

ecf_apec(kT: float, Z: float = 0.3, z: float = 0.0, nH: float = 0.03, apec=None, absorb: bool = True) → float

ecf_apec_table(z: float, nH: float = 0.03, Z: float = 0.3, kT_grid=None)

Return (kT_grid, ecf_grid) and a log-kT interpolator for the gas ECF.

ecf_from_rest(e_rest, s_rest, z: float, nH: float = 0.03, absorb: bool = True) → float

ECF for a rest-frame photon spectrum s_rest = dN/dE (any norm).

e_rest : keV ; nH : 1e22 cm^-2 ; returns cts/s per erg/s/cm^2.

ecf_powerlaw(photon_index: float = 1.9, z: float = 0.0, nH: float = 0.03, absorb: bool = True) → float

class hod_mod.gas.GasDensityDPM(model: int = 2, r_max_over_r200: float = 3.0, n_gl: int = 200, sigma_scatter: float = 0.0, concentration_model: str = 'diemer19')

Bases: object

DPM electron density profile for X-ray emissivity (Oppenheimer+2025).

Reference: Oppenheimer et al. 2025, arXiv:2505.14782, Table 1.

The profile uses a generalized NFW (gNFW) shape function:

\[f(x|\alpha) = x^{-\alpha_{\rm in}} \left(1 + x^{\alpha_{\rm tr}}\right)^{(\alpha_{\rm in} - \alpha_{\rm out})/\alpha_{\rm tr}}\]

where \(x = r/R_s\) and \(R_s = R_{200}/c_{\rm DPM}\) with \(c_{\rm DPM} = 2.772\) (Table 1 of arXiv:2505.14782).

The electron density is:

\[n_e(r, M_{200}, z) = n_{e0}\,f(x|\alpha^{n_e}) \,E(z)^{\gamma^{n_e}}\,M_{12}^{\beta^{n_e}}\]

where \(M_{12} = M_{200}/(10^{12}\,M_\odot)\), \(E(z) = H(z)/H_0\), and \(n_{e0}\) is normalised so that \(n_e(0.3 R_{200}, 10^{12}\,M_\odot, z=0) = n_{e,0.3}\).

Three calibrated models are provided (Table 1 of arXiv:2505.14782):

• Model 1 — self-similar (β=0)

• Model 2 — cluster-reduced slope (β=0.36)

• Model 3 — slope-changing outer profile

Parameters:

• model (int (1, 2, or 3), default 2)

• r_max_over_r200 (float (default 3.0))

• n_gl (int (default 200))

density_3d(r: ndarray, m200: float, r200: float, z: float, omega_m: float, c200c: float | None = None) → ndarray

Electron number density n_e(r|M₂₀₀, z) [cm⁻³].

Parameters:

• r (radii [Mpc/h])

• m200 (M₂₀₀ [Msun/h])

• r200 (R₂₀₀ [Mpc/h])

• z (redshift)

• omega_m (matter fraction Ω_m)

• c200c (concentration c₂₀₀c for this halo (pre-computed by the caller).) – If None, falls back to the fixed class constant _C_DPM.

density_uk(k_arr: ndarray, m200_arr: ndarray, r200_arr: ndarray, z: float, theta_cosmo: dict) → ndarray

FT of the electron density: ñ_e(k|M) = 4π ∫ n_e(r) j₀(kr) r² dr.

Output units: (Mpc/h)³ cm⁻³. Multiply by (Mpc_cm/h)³ to convert to dimensionless (but this is done at the power-spectrum level when needed).

Parameters:

• k_arr ((Nk,) [h/Mpc])

• m200_arr ((NM,) [Msun/h])

• r200_arr ((NM,) [Mpc/h])

• z (redshift)

• theta_cosmo (dict with ‘Omega_m’)

Returns:

uk ((Nk, NM) [(Mpc/h)³ cm⁻³])

emissivity_full_uk(k_arr: ndarray, m200_arr: ndarray, r200_arr: ndarray, z: float, theta_cosmo: dict, pressure_profile: PressureProfileDPM, metallicity_profile: MetallicityProfileDPM, cooling_fn: ApecCoolingTable) → ndarray

FT of n_e²(r) × Λ_APEC(T(r), Z(r)) — X-ray surface brightness emissivity.

Evaluates the full temperature- and metallicity-dependent APEC cooling function at each quadrature node:

\[\varepsilon(r) = n_e^2(r) \times \Lambda_{n_e^2}(T(r), Z(r))\]

where \(T(r) = P_{\rm DPM}(r) / n_e(r)\) [keV] (ideal gas law), \(Z(r)\) comes from MetallicityProfileDPM [Z_sun], and \(\Lambda_{n_e^2}\) is the band-integrated APEC emissivity from ApecCoolingTable (0.83 × Λ_{\rm APEC} converting \(n_e n_H \to n_e^2\)).

Parameters:

• pressure_profile (PressureProfileDPM)

• metallicity_profile (MetallicityProfileDPM)

• cooling_fn (ApecCoolingTable) – Precomputed APEC cooling table. Instantiate once and reuse.

Returns:

uk ((Nk, NM) [erg cm³ s⁻¹ × (Mpc/h)³ cm⁻⁶])

emissivity_full_uk_bands(k_arr: ndarray, m200_arr: ndarray, r200_arr: ndarray, z: float, theta_cosmo: dict, pressure_profile: PressureProfileDPM, metallicity_profile: MetallicityProfileDPM, cooling_fns: list) → ndarray

Multi-band version of emissivity_full_uk().

Computes the emissivity FT X̃_b(k|M) for a LIST of energy-band cooling tables in ONE batched spherical-Bessel FT. The bands share n_e(r,M), T(r,M), Z(r,M) and the j₀ geometry; only Λ_b differs, so this is ≈ the cost of a single emissivity_full_uk() (the per-band table eval is cheap). Used by the energy-band (temperature-resolved) joint fit.

Parameters:

cooling_fns (list of ApecCoolingTable) – One per band (e.g. 15 × ApecCoolingTable(emin, emax)).

Returns:

uk ((Nb, Nk, NM) [erg cm³ s⁻¹ × (Mpc/h)³ cm⁻⁶])

emissivity_uk(k_arr: ndarray, m200_arr: ndarray, r200_arr: ndarray, z: float, theta_cosmo: dict) → ndarray

FT of n_e²(r) without cooling function weighting.

Deprecated since version Use: emissivity_full_uk() with PressureProfileDPM, MetallicityProfileDPM, and ApecCoolingTable.

Output units: (Mpc/h)³ cm⁻⁶.

class hod_mod.gas.MetallicityProfileDPM

Bases: object

DPM gas metallicity profile (Oppenheimer+2025, arXiv:2505.14782, Eq. 4).

All three DPM models share the same metallicity profile (Table 1):

\[Z(r, M, z) = Z_0 \, f(r/R_s \mid \alpha^Z)\]

with \(\alpha_{\rm in}^Z = 0\), \(\alpha_{\rm tr}^Z = 0.5\), \(\alpha_{\rm out}^Z = 0.7\), \(\beta^Z = 0\), \(\gamma^Z = 0\) (no mass or redshift dependence). The normalisation is \(Z(0.3 R_{200}) = 0.3\,Z_\odot\).

The same _C_DPM = 2.772 scale radius convention is used.

This profile is used by GasDensityDPM.emissivity_full_uk() to evaluate the metallicity-dependent X-ray cooling function Λ(T, Z).

metallicity_3d(r: ndarray, r200: float) → ndarray

Gas metallicity Z(r) [Z_sun].

No mass or redshift dependence (β^Z = γ^Z = 0).

Parameters:

• r (radii [Mpc/h])

• r200 (R₂₀₀ [Mpc/h])

class hod_mod.gas.PressureProfileA10(r_max_over_r500c: float = 6.0, n_gl: int = 200)

Bases: object

Arnaud+2010 generalized NFW electron pressure profile for tSZ.

Reference: Arnaud, Pratt, Piffaretti et al. 2010, A&A 517, A92 (arXiv:0910.1234), Eq. 11 and Table 1.

The “universal pressure profile” is:

\[P_e(r|M_{500c}, z) = 1.65 \times 10^{-3}\,h_{70}^2\,E(z)^{8/3} \left[\frac{M_{500c}}{3 \times 10^{14}\,h_{70}^{-1}\,M_\odot} \right]^{2/3 + \alpha_p} p(r/R_{500c}) \quad [\text{keV cm}^{-3}]\]

with shape function:

\[p(x) = \frac{P_0}{(c_{500}\,x)^\gamma \left[1 + (c_{500}\,x)^\alpha\right]^{(\beta-\gamma)/\alpha}}\]

Universal parameters from Table 1 of arXiv:0910.1234: P₀=8.403, c₅₀₀=1.177, γ=0.3081, α=1.0510, β=5.4905, α_p=0.12.

Parameters:

• r_max_over_r500c (float) – Integration truncation radius as a multiple of R₅₀₀c (default 6).

• n_gl (int) – Gauss-Legendre quadrature nodes (default 200).

pressure_uk(k_arr: ndarray, m200_arr: ndarray, r200_arr: ndarray, c200_arr: ndarray, z: float, theta_cosmo: dict) → ndarray

Pressure-profile Fourier transform ỹ(k|M) in (Mpc/h)².

Defined as:

\[\tilde{y}(k|M,z) = \frac{\sigma_T}{m_e c^2} \frac{\mathrm{Mpc\_cm}}{h} \times 4\pi \int_0^{r_{\max}} P_e(r|M,z)\, \frac{\sin(kr)}{kr}\,r^2\,\mathrm{d}r\]

with r in Mpc/h and P_e in keV cm⁻³. The prefactor (σ_T/m_e c²)×(Mpc_cm/h) has units cm³/(keV·Mpc/h) so that

\[[\tilde{y}] = \frac{\mathrm{cm}^3}{\mathrm{keV}\cdot(\mathrm{Mpc}/h)} \times \frac{\mathrm{keV}}{\mathrm{cm}^3} \times (\mathrm{Mpc}/h)^3 = (\mathrm{Mpc}/h)^2\]

The 3D galaxy×y cross-power P_{gy}(k) then has units (Mpc/h)², and the projected Σ_y(r_p) = (1/π) ∫ P_{gy}(k) J₀(k r_p) k dk is dimensionless (Compton-y parameter).

Parameters:

• k_arr ((Nk,) [h/Mpc])

• m200_arr ((NM,) [Msun/h])

• r200_arr ((NM,) [Mpc/h])

• c200_arr ((NM,) concentration at the overdensity stored in the static cache)

• z (redshift)

• theta_cosmo (dict with keys ‘h’, ‘Omega_m’)

Returns:

uk ((Nk, NM) [(Mpc/h)²])

class hod_mod.gas.PressureProfileDPM(model: int = 2, r_max_over_r200: float = 3.0, n_gl: int = 200)

Bases: object

DPM electron pressure profile for tSZ (Oppenheimer+2025, arXiv:2505.14782).

Reference: Table 1 of arXiv:2505.14782 — 3 calibrated models for the generalized NFW pressure profile.

The profile uses the same gNFW shape as GasDensityDPM (Eq. 1), with the addition of a mass-dependent outer slope (Eq. 5):

\[\alpha_{\rm out}(M) = \alpha_{\rm out,12} + \alpha_{\rm out,var} \log_{10}(M_{200} / 10^{12}\,M_\odot/h)\]

The pressure profile is (Eq. 2):

\[P(r, M, z) = P_0 \, f(r/R_s \mid \alpha(M)) \, E(z)^{\gamma^P} \, M_{12}^{\beta^P}\]

normalised so that \(P(0.3 R_{200}, 10^{12}\,M_\odot/h, z=0) = P_{0.3}\).

The pressure_uk method uses the same unit convention as PressureProfileA10 and outputs in (Mpc/h)².

Parameters from Table 1 (DPM paper arXiv:2505.14782), converted to keV cm⁻³:

Param	Model 1	Model 2	Model 3
P_0.3	4.09e-4	1.15e-4	7.10e-5
α_in^P	0.3	0.3	−0.6
α_tr^P	1.3	1.3	0.2
α_out^P	4.1	4.1	2.0
β^P	2/3	0.85	0.92
γ^P	8/3	8/3	8/3

Note

The paper (arXiv:2505.14782 Table 1) lists P_0.3 as 409, 115, 71 in meV cm⁻³. The values stored here have been converted to keV cm⁻³ (factor 10⁻⁶) so that pressure_uk and _pressure_3d return physically correct units. Sanity check: T = P_0.3 / ne_0.3 gives 0.70, 2.36, 1.46 keV for models 1–3 at M=10¹² M☉/h, z=0 — consistent with observed group/cluster temperatures at those masses.

Parameters:

• model (int (1, 2, or 3), default 2)

• r_max_over_r200 (float (default 3.0))

• n_gl (int (default 200))

pressure_uk(k_arr: ndarray, m200_arr: ndarray, r200_arr: ndarray, z: float, theta_cosmo: dict) → ndarray

DPM pressure-profile Fourier transform ỹ(k|M) in (Mpc/h)².

Same interface and unit convention as PressureProfileA10.pressure_uk(). The tSZ Compton-y prefactor σ_T/(m_e c²) × (Mpc_cm/h) is applied assuming P_0.3 is in keV cm⁻³.

Parameters:

• k_arr ((Nk,) [h/Mpc])

• m200_arr ((NM,) [Msun/h])

• r200_arr ((NM,) [Mpc/h])

• z (redshift)

• theta_cosmo (dict with keys ‘h’, ‘Omega_m’)

Returns:

uk ((Nk, NM) [(Mpc/h)²])

hod_mod.gas.load_ecf_tables(sample: str)

Load the precomputed per-component ECF tables for a GALxEVT sample.

Returns (ecf_gas_interp, ecf_agn, ecf_fixed) where ecf_gas_interp is a callable T_keV -> ECF_gas(T) [cts/s per erg/s/cm²] (log-T interpolation, flat extrapolation), ecf_agn is the AGN power-law ECF, and ecf_fixed is the GALxEVT pipeline’s fixed conversion ARF_1keV/C. Built by scripts/galaxies/build_ecf_tables.py.

hod_mod.gas.m200_to_m500c(m200_arr: ndarray, c200_arr: ndarray, r200_arr: ndarray, rho_crit_z: float) → tuple[ndarray, ndarray]

Convert M₂₀₀ to M₅₀₀c and R₅₀₀c using NFW enclosed-mass bisection.

Given that the NFW halo profile has an enclosed mass:

\[M(<r) = M_{200} \frac{g(c_{200} r/r_{200})}{g(c_{200})}\]

we solve for \(r_{500c}\) such that:

\[M(<r_{500c}) = \frac{4\pi}{3} \times 500 \times \rho_{\rm crit}(z) \times r_{500c}^3\]

Parameters:

• m200_arr ((NM,) [Msun/h])

• c200_arr ((NM,) NFW concentration at Δ=200 (any ρ_ref))

• r200_arr ((NM,) [Mpc/h])

• rho_crit_z (float [(Msun/h)/(Mpc/h)³] comoving critical density at z)

Returns:

• m500c ((NM,) [Msun/h])

• r500c ((NM,) [Mpc/h])

hod_mod.gas.temperature_from_dpm(pressure_profile, density_profile, r: ndarray, m200: float, r200: float, z: float, theta_cosmo: dict) → ndarray

Gas temperature T(r, M, z) [keV] from DPM pressure and density profiles.

Uses the ideal gas law T = P_e / (n_e k_B). Convenience wrapper that evaluates both DPM profiles at the same radii.

Parameters:

• pressure_profile (PressureProfileDPM instance)

• density_profile (GasDensityDPM instance)

• r (radii [Mpc/h])

• m200 (M₂₀₀ [Msun/h])

• r200 (R₂₀₀ [Mpc/h])

• z (redshift)

• theta_cosmo (dict with key ‘Omega_m’)

Returns:

T ((Nr,) [keV])

hod_mod.gas.temperature_from_profiles(pressure: ndarray, n_electron: ndarray) → ndarray

Gas temperature T = P / (n_e k_B) [keV].

Used internally by GasDensityDPM.emissivity_full_uk() to build the temperature map from DPM pressure and density profiles (Section 3.1.1 of arXiv:2505.14782).

Parameters:

• pressure (P_e [keV cm⁻³] — from PressureProfileDPM)

• n_electron (n_e [cm⁻³] — from GasDensityDPM)

Returns:

T ([keV])

hod_mod.gas.xray_cooling_function(T_keV: ndarray, Z_solar: ndarray, alpha_T: float = 0.5, alpha_Z: float = 1.0, Lambda_0: float = 3e-23) → ndarray

Simplified power-law cooling function Λ(T, Z) [erg cm³ s⁻¹].

Deprecated since version Use: ApecCoolingTable instead, which evaluates the full APEC plasma code tables (soxs) over the specified energy band.

Models the 0.5–2 keV soft X-ray volume emissivity coefficient:

\[\Lambda(T, Z) = \Lambda_0 \left(\frac{T}{1\,\text{keV}}\right)^{\alpha_T} \left(\frac{Z}{0.3\,Z_\odot}\right)^{\alpha_Z}\]

—

Matter Power Spectrum (Halo Model)

(hod_mod.core.halo_model)

The 1-halo + 2-halo decomposition of the nonlinear matter power spectrum is (Cooray & Sheth 2002) [CooraySheth2002]:

\[P_{\rm mm}(k) = P^{1h}_{\rm mm}(k) + P^{2h}_{\rm mm}(k)\]

The 1-halo term (pairs within the same halo) dominates at \(k \gtrsim 1\,h\,{\rm Mpc}^{-1}\):

\[P^{1h}_{\rm mm}(k) = \int \frac{dn}{dM} \left(\frac{M}{\bar{\rho}_m}\right)^2 u^2(k|M)\,dM\]

The 2-halo term (pairs in different halos) dominates at large scales:

\[P^{2h}_{\rm mm}(k) = P_{\rm lin}(k) \left[\int \frac{dn}{dM} b(M)\,\frac{M}{\bar{\rho}_m}\, u(k|M)\,dM\right]^2\]

Full halo model matter power spectrum P_mm(k) = P^{1h}_mm + P^{2h}_mm.

Implements the Asgari et al. (2023) halo model for the matter auto-power spectrum. The 1-halo term captures intra-halo (shot-noise) clustering; the 2-halo term recovers linear clustering on large scales.

\[ \begin{align}\begin{aligned}P^{1h}_{mm}(k) = \frac{1}{\bar{\rho}_m^2} \int M^2\, \hat{u}_m^2(k,M)\, n(M)\, dM\\P^{2h}_{mm}(k) = P_{\rm lin}(k)\, \left[\frac{1}{\bar{\rho}_m} \int M\, \hat{u}_m(k,M)\, b(M)\, n(M)\, dM\right]^2\end{aligned}\end{align} \]

References

Asgari et al. 2023, arXiv:2303.08752 — halo model review (Eqs. 34–35) Cooray & Sheth 2002, Phys.Rep. 372, 1 — NFW window function (Eq. 11)

class hod_mod.core.halo_model.HaloModelPowerSpectrum(hmf, halo_profile, pk_lin, m_min: float = 10000000000.0, m_max: float = 1e+16, n_m: int = 100)

Bases: object

Matter power spectrum P_mm(k) from the halo model.

Combines a 1-halo and 2-halo term using NFW profile window functions and the chosen halo mass function and bias.

\[ \begin{align}\begin{aligned}P^{1h}_{mm}(k) = \frac{1}{\bar{\rho}_m^2} \int M^2\, \hat{u}_m^2(k,M)\, n(M)\, dM\\P^{2h}_{mm}(k) = P_{\rm lin}(k)\, \left[\frac{1}{\bar{\rho}_m} \int M\, \hat{u}_m(k,M)\, b(M)\, n(M)\, dM\right]^2\end{aligned}\end{align} \]

where \(\hat{u}_m(k,M)\) is the NFW normalized Fourier transform (Cooray & Sheth 2002 Eq. 11, implemented in nfw_uk), \(n(M)\) is the halo mass function, \(b(M)\) is the linear halo bias, and \(\bar{\rho}_m = \Omega_m \rho_{\rm crit,0}\).

On large scales (k → 0): \(\hat{u}_m → 1\) so the 2-halo integral → 1 (by the mass-weighted bias normalization), recovering \(P^{2h}_{mm} → P_{\rm lin}\) as expected.

Parameters:

• hmf (HaloMassFunction) – Provides dndm(m, z, theta) and bias(m, z, theta).

• halo_profile (HaloProfile) – Provides rho_s_and_rs and concentration (colossus c–M relation).

• pk_lin (LinearPowerSpectrum) – Provides pk_linear(k, z, theta) for the 2-halo term.

• m_min, m_max (float [M_sun/h]) – Mass integration limits.

• n_m (int) – Number of log-spaced mass bins.

pk_1h_mm(k_arr: ndarray, z: float, theta: dict) → Array

1-halo matter power spectrum (Asgari+2023 Eq. 34).

\[P^{1h}_{mm}(k) = \frac{1}{\bar{\rho}_m^2} \int M^2\, \hat{u}_m^2(k,M)\, \frac{dn}{dM}\, dM\]

Parameters:

k_arr ([h/Mpc], shape (Nk,))

Returns:

p1h ([(Mpc/h)³], shape (Nk,))

pk_2h_mm(k_arr: ndarray, z: float, theta: dict) → Array

2-halo matter power spectrum (Asgari+2023 Eq. 35).

\[P^{2h}_{mm}(k) = P_{\rm lin}(k)\, \left[\frac{1}{\bar{\rho}_m} \int M\, \hat{u}_m(k,M)\, b(M)\, \frac{dn}{dM}\, dM\right]^2\]

Parameters:

k_arr ([h/Mpc], shape (Nk,))

Returns:

p2h ([(Mpc/h)³], shape (Nk,))

pk_mm(k_arr: ndarray, z: float, theta: dict) → Array

Total matter power spectrum P_mm = P^{1h}_mm + P^{2h}_mm.

Parameters:

k_arr ([h/Mpc], shape (Nk,))

Returns:

pk ([(Mpc/h)³], shape (Nk,))

—

Key references

[PressSchechter1974], [ShethTormen1999], [Jenkins2001], [Nishimichi2019], [SeljakWarren2004], [WrightBrainerd2000], [BryanNorman1998].