The extended model: implementation of the missing physics

The implemented counterpart of What the model does not yet contain.

Between 2026-07-03 and 2026-07-04 (branch feature/missing-physics) the propositions of What the model does not yet contain were implemented in four waves, growing the differentiable parameter vector from the tier-2 61 entries to 90 (waves 1–3), then to 102 (the Tier-3 forecast: multi-wavelength maps, band LFs, z < 2, M* > 10⁹ SED calibrations), 106 (wave 4, galaxy morphology) and 111 (Tier-4 forecast: the morphology observables, the morphology observables) — every addition fiducial-preserving, so the Tier-2 forecast: nothing fixed (90 parameters) predictions are reproduced exactly (or to a stated tolerance) at the fiducial, and every mechanism carries an exact invariant that is enforced by tests/test_missing_physics.py. This page documents the architecture, the physics of each sector with its equations and parameters, the survey noise models, the invariants, usage, and the approximations deliberately made.

Architecture

Three design rules govern every extension:

1. Append-only, fiducial-preserving parameters. New entries land at the end of PARAM_NAMES (MISSING_PHYSICS = PARAM_NAMES[61:], 29 names) with fiducials equal to the previously hard-coded behaviour: promoted constants keep their values, switches sit at their “off” value (eta_w_norm = 0, sum_mnu = 0 eV, dlx_quenched = 0), and physical relations take their literature centres. Tier-1/tier-2 scripts pin the whole extension by default (fix=params.TIER2_EXTENSION).

2. Static dispatch, no traced branches. Where behaviour must change with the parameter set rather than parameter values, the branch is on dictionary membership or a constructor flag — static under JAX tracing. The growth factor switches to the CPL ODE because the forecast cosmology dict has a w0 key; production ΛCDM dicts keep the Carroll formula bit-identical. Sample definitions (sfq, ssfr_cut, cm_relation, pk_correction) are constructor arguments.

3. Distill-to-table for non-traceable dependencies. Anything that cannot be traced (CAMB, soxs/APEC) is evaluated once in numpy, stored as an npz shipped with the package, and consumed through differentiable jnp interpolation — the pattern established by apec_bands and extended by pk_camb_ratio.

The parameter block

The 29 MISSING_PHYSICS parameters (fiducial → prior width)

Name	Fiducial	Prior σ	Sector	Enters through / constrained by
eps_sn	0.1	0.3	baryon	SN channel of the energy closure; low-mass SMF, group L_X–M
w0, wa	−1, 0	1.0, 2.0	cosmology	CPL growth ODE + all Limber distances + E(z); BAO/FS, 3×2pt
sum_mnu [eV]	0	0.25	cosmology	ν suppression of the P(k) shape (σ8-anchored); cluster counts + CMB lensing
log10_Mq_cen, mu_q_cen	11.78, 0.41	1.0, 1.0	sfq	ZM16 central quenching Weibull; f_Q(M*, z), split w_p/ΔΣ
log10_Mq_sat, mu_q_sat	12.19, 0.24	1.0, 1.0	sfq	ZM16 satellite quenching; group catalogues
dlx_quenched [dex]	0	1.0	sfq	L_X–M offset of quenched samples; eROSITA SF/Q CGM stacks
agn_xi_rx, agn_xi_rm	0.60, 0.78	0.3, 0.3	agn	fundamental-plane slopes; radio LF + FP priors
agn_b_r, agn_sig_r	7.33, 0.88	2.0, 0.3	agn	FP zero point and scatter
log10_M0_hi, log10_Mmin_hi, alpha_hi	10.63, 12.3, 0.24	1, 1, 0.3	coldgas	M_HI(M_h); HIMF, HI HOD, 21 cm crosses
eta_w_norm, alpha_w	0, 1.0	0.5, 0.7	baryon	wind mass loading in η(M); kSZ, group baryon census
ssfr_ms_norm, ssfr_ms_slope, ssfr_ms_zs	−9.9, −0.25, 0	0.5, 0.3, 2.0	sfq	the star-forming main sequence + evolution; MS measurements
dhi_quenched [dex]	0	1.0	coldgas	HI deficit of quenched centrals; xGASS conditional gas fractions
sigma_ms, dssfr_q [dex]	0.30, −1.6	0.15, 0.5	sfq	the double-lognormal p(log sSFR | M*); sSFR distributions
loii_norm	40.85	0.3	sfq	log L_[OII]/SFR calibration [Kennicutt1998]; [OII] LFs
f_loud0, beta_loud, b_jet	0.01, 2.0, 41.0	0.05, 1.5, 2.0	agn	radio-loud jet population; bright-end radio LF
agn_bc_ir	−0.52	0.3	agn	L_IR(6 μm)/L_bol; IR AGN LF and counts

Wave 1 — cosmology, concentration, quenching, fundamental plane, HI

Beyond-ΛCDM growth and neutrinos

The linear growth factor is obtained from the growth ODE in \(x = \ln a\),

\[D'' + \Big[2 + \frac{{\rm d}\ln E}{{\rm d}\ln a}\Big] D' - \tfrac{3}{2}\,\Omega_m(a)\,D = 0, \qquad E^2(a) = \Omega_m a^{-3} + (1-\Omega_m)\, a^{-3(1+w_0+w_a)}\,e^{-3 w_a (1-a)},\]

integrated by a fixed-grid RK4 through jax.lax.scan (160 steps from \(\ln a = -7\); jit/vmap/jacfwd-safe, no extra dependency) with the matter-domination growing mode \(D \propto a\) as the initial condition (_growth_factor_cpl_jax()). The dispatcher growth_factor() selects it whenever the cosmology dict carries a w0/wa key; at \((w_0, w_a) = (-1, 0)\) it reproduces the Carroll et al. (1992) fitting formula to a few \(\times 10^{-4}\) and the exact ΛCDM quadrature \(D \propto E \int {\rm d}a\,(aE)^{-3}\) to <0.2% (both tested). The CPL parameters also propagate through geometry: every Limber comoving distance and every \(E(z)\) factor in the forecast (halo definitions, L_X/kT/P500 self-similar scalings, the energy closure) uses the CPL expressions of hod_mod.core.distances / ForwardModel._e2z.

Massive neutrinos suppress the shape of the EH98 spectrum (_nu_suppression()):

\[\frac{P_\nu(k)}{P_0(k)} = 1 - 8 f_\nu \cdot \tfrac{1}{2}\Big[1 + \tanh\frac{\ln (k/k_{\rm nr})}{1.5}\Big], \quad f_\nu = \frac{\Sigma m_\nu / 93.14\,{\rm eV}}{\Omega_m h^2},\]

with \(k_{\rm nr} \simeq 0.018\sqrt{\Omega_m m_{\nu,i}/1\,{\rm eV}}\) h/Mpc (three degenerate species) and a small mass floor inside the square root keeping the jacfwd column finite at the massless fiducial. Because the suppression multiplies the shape used everywhere (the HMF path, the σ8 anchor and the 2-halo term), σ8 keeps its measured-amplitude meaning and only the shape responds — exactly massless (bit-identical) at \(\Sigma m_\nu = 0\). The wave-2 CAMB-ratio derivative row (below) corrects the residual of this first-order form.

Cosmology-dependent concentration

ForwardModel(cm_relation="diemer15") replaces the Planck13-frozen Dutton14 fit with the Diemer & Kravtsov / Diemer & Joyce model [DiemerKravtsov2015], [DiemerJoyce2019],

\[\begin{split}c_{\rm 200c} = \frac{c_{\min}}{2} \left[\left(\frac{\nu}{\nu_{\min}}\right)^{-\alpha} + \left(\frac{\nu}{\nu_{\min}}\right)^{\beta}\right], \qquad \begin{aligned} c_{\min} &= \phi_0 + \phi_1\, n_{\rm eff},\\ \nu_{\min} &= \eta_0 + \eta_1\, n_{\rm eff}, \end{aligned}\end{split}\]

with the DJ19 median parameters (φ₀, φ₁, η₀, η₁, α, β) = (6.58, 1.27, 7.28, 1.56, 1.08, 1.77), peak height \(\nu = 1.686/\sigma(M, z)\) from the JAX σ(M) machinery, and \(n_{\rm eff} = {\rm d}\ln P/{\rm d}\ln k\) evaluated at \(k_R = \kappa\, 2\pi/R(M)\) with the DJ19 calibration κ = 0.42 (ForwardModel._neff_eh98, differentiable). Both inputs carry cosmology, so c(M) finally responds to σ8/n_s/h — and, through the growth factor, to w0/wa/Σm_ν.

Important

Implementing this exposed a factor-~2 bug in the pre-existing hod_mod.core.concentration.c_diemer15(): the second power law was applied as \(c_{\min} x^{-\alpha}(1+x^{\beta})\) (an effective \(\beta-\alpha\) exponent, no ½) instead of the DK15 Eq. 9 form above, and the n_eff helper used κ = 1. Both are fixed (core and forecast); values now match the COLOSSUS [Colossus2018] anchors — \(c_{\rm 200c}(10^{12}, z{=}0) = 8.8\) (~9.5), \(c_{\rm 200c}(10^{15}) = 4.7\) (~5.0) — enforced by test_dk15_concentration_is_cosmology_dependent.

The SF/quiescent split

ForwardModel(sfq="sf"|"q") weights the occupations with the Zu & Mandelbaum halo-quenching Weibull fractions [ZuMandelbaum2016] (f_red_cen_zu16() / f_red_sat_zu16, the in-repo MAP fiducials):

\[f_{\rm Q}(M_h) = 1 - \exp\!\left[-\big(M_h/M_{\rm q}\big)^{\mu_{\rm q}}\right],\]

separately for centrals and satellites. By construction the star-forming and quenched samples sum exactly to the unsplit sample — the regression invariant, tested to \(10^{-12}\) on occupations, abundances and the SFR density. dlx_quenched shifts the L_X–M relation of quenched samples only (_lx_kt_of), turning the eROSITA CGM star-forming-vs-quiescent measurements [Zhang2025CGM] into a one-parameter test through the split-cell \(C_\ell^{gX}\). In the tier-2 assembly, Tier2Forecast(split_sfq=True) doubles the cell blocks while the shell X-ray observables (XLF, band \(C_\ell^{XX}\)) keep a dedicated unsplit model, so the quenched offset cannot leak into the total-gas X-ray auto.

The fundamental plane of black-hole activity

The radio luminosity attaches to the Powell chain’s per-halo \((M_{\rm BH}, L_X)\) through [MerloniHeinzDiMatteo2003] (priors from [Gultekin2019]):

\[\log L_R = \xi_{RX} \log L_X + \xi_{RM} \log M_{\rm BH} + b_R .\]

Because the \(M_{\rm BH}\) scatter \(\sigma_{lm}\) is the same deviate in \(L_X\) and \(M_{\rm BH}\), it enters the radio kernel with coefficient \((\xi_{RX}+\xi_{RM})\), and the FP scatter adds in quadrature:

\[\sigma_R^{\rm eff} = \sqrt{(\xi_{RX}+\xi_{RM})^2 \sigma_{lm}^2 + \sigma_R^2}.\]

The rlf observable (5 GHz νL_ν luminosity function) reuses the XLF kernel machinery; the exact identity — at \((\xi_{RX}, \xi_{RM}, b_R, \sigma_R) = (1, 0, 0, 0)\) and jets off, the rlf equals the hard-band XLF on the same abscissas — is tested to \(10^{-10}\).

The HI sector

The Villaescusa-Navarro et al. halo model [VillaescusaNavarro2018HI]:

\[M_{\rm HI}(M_h) = M_0 \left(\frac{M_h}{M_{\min}}\right)^{\alpha} \exp\!\left[-\left(\frac{M_{\min}}{M_h}\right)^{0.35}\right],\]

feeding (i) the himf HI mass function through a 0.35 dex lognormal conditional (the SMF/XLF kernel pattern), and (ii) the 21 cm × galaxy cross cl_gHI — the \(C_\ell^{gX}\) machinery with the NFW-distributed HI mass as the tracer field, so \(C_\ell^{g{\rm HI}} \propto M_0\) exactly (tested). dhi_quenched (wave 2) dims the HI around quenched centrals, following the NeutralUniverseMachine phenomenology [GuoNUM2023].

Supernova coupling

eps_sn promotes the formerly fixed \(\varepsilon_{\rm SN} = 0.1\) of the energy-closure SN channel into the vector. Note for testing: at the default fiducial the log10_M_pivot slot (reinterpreted as \(\log_{10}\varepsilon_{\rm AGN}\) in closure mode) saturates the energy budget’s min(), so exercising the SN derivative requires a physical coupling (the gate uses \(\log_{10}\varepsilon_{\rm AGN} = -2\)).

Wave 2 — CAMB-quality P(k), winds, main sequence, radio/HI surveys

The linearized CAMB ratio

A Fisher forecast needs the spectrum and its first derivatives to be accurate at the fiducial — nothing more. That reduces “CAMB-quality P(k)” from an emulator problem to eleven CAMB evaluations (hod_mod.forecast.pk_camb_ratio):

\[\ln R(k;\theta) \simeq \ln R_0(k) + \sum_i \frac{\partial \ln R}{\partial \theta_i}(k)\,(\theta_i - \theta_i^{\rm fid}), \qquad R = \frac{P^{\rm shape}_{\rm CAMB}}{P^{\rm shape}_{\rm EH98}},\]

with both shapes pivot-normalised at k = 0.05 h/Mpc and derivative rows for (h, Ω_b, Ω_m, n_s, Σm_ν). The fiducial ratio confirms the expected EH98 error: \(\max|\ln R_0| = 3.8\%\) around the BAO scale. The Σm_ν row is computed against the EH98 shape including the wave-1 tanh suppression, so it corrects that form’s residual to CAMB accuracy. w0/wa need no rows (they do not alter the z = 0 transfer shape). The table ships as hod_mod/data/pk_ratio/camb_eh98_ratio.npz; evaluation is two jnp.interp calls and an exp, applied inside pk_shape so the HMF, σ8 anchor and 2-halo term stay mutually consistent. Opt-in: ForwardModel(pk_correction="camb_linear"); rebuild with python -m hod_mod.forecast.pk_camb_ratio (needs CAMB). Beyond the tabulated k range the log-ratio is edge-clamped (deep power-law tail).

Wind mass loading

The Muratov-style [Muratov2015] SN wind,

\[\eta_w(M_h) = \eta_{w,0} \left(\frac{V_c(M_h)}{200\,{\rm km/s}}\right)^{-\alpha_w}, \qquad \eta_{\rm eff}(M) = \frac{\eta_{\rm sigmoid}(M)}{1 + \eta_w(M)},\]

couples into the existing gas-concentration slot: stronger winds puff out the low-mass hot gas, which the ΔΣ baryon split, the X-ray emissivity and the tSZ pressure all see. \(\alpha_w\) interpolates the momentum-driven (1) to energy-driven (2) scalings [SomervilleDave2015]. At the fiducial \(\eta_{w,0} = 0\) the tier-2 sigmoid is reproduced bit-identically (tested), and the \(\eta_{w,0}\) derivative is live; \(\alpha_w\) is flat exactly at the fiducial (the prior bounds it) — a documented property, not a bug.

Surveys: radio, HI, and the local-HIMF lesson

RadioSurvey (LoTSS-like: \(f_{\rm sky} = 0.13\), a νL_ν(5 GHz)-equivalent detection threshold of \(3\times10^{-22}\) erg/s/cm² driving the L_lim(z) completeness of the radio LF) and HISurvey (ALFALFA-like: \(f_{\rm sky} = 0.16\), the standard 21 cm mass–flux relation \(M_{\rm lim} = 2.36\times10^5 d_L^2 S_{\rm int}\), plus an effective noise-to-signal recipe for the 21 cm intensity-mapping cross) enter Tier2Forecast(include_radio=True, include_hi=True).

One design lesson worth recording: attaching the HIMF to the Δz = 0.1 shells fails physically — at \(z_{\rm hi} = 0.3\) the flux limit gives \(M_{\rm lim} \approx 2.5\times10^{11} M_\odot\), flagging every bin. Blind HIMFs are local measurements, so the assembly builds one dedicated hi_local block (z ≤ 0.06), matching how ALFALFA [Jones2018ALFALFA] and MIGHTEE-HI [Ponomareva2023] actually operate.

Wave 3 — continuous sSFR, [OII], jets, infrared

The continuous sSFR distribution and selection

The conditional sSFR distribution is the double lognormal

\[p(\log {\rm sSFR}\,|\,M_*) = f_{\rm Q}\,\mathcal{N}\big(\mu_{\rm MS} + \Delta_{\rm Q},\, \sigma_{\rm Q}\big) + (1-f_{\rm Q})\,\mathcal{N}\big(\mu_{\rm MS},\, \sigma_{\rm MS}\big),\]

with the Speagle-like main sequence [Speagle2014] \(\mu_{\rm MS} = {\rm ssfr\_ms\_norm} + {\rm ssfr\_ms\_slope}\, (\log M_* - 10.5)\) (evolution through the standard _zs mechanism), free MS scatter sigma_ms and quenched offset dssfr_q (\(\sigma_{\rm Q} = 0.5\) dex fixed). ForwardModel(ssfr_cut=...) selects sSFR-thresholded samples (ELG-like) by the per-component Gaussian survival \(S = \tfrac12\,{\rm erfc}[({\rm cut}-\mu)/\sqrt{2}\sigma]\) applied inside the occupation (_sfq_weights). The selection composes with the SF/Q split: SF-cut + Q-cut ≡ mixture-cut exactly, and a cut at \(-\infty\) recovers the unsplit sample exactly (both tested).

Two observables consume it:

• sfrd — the cell’s SFR density \(\rho_{\rm SFR} = \sum_{\rm pop} n_{\rm pop}\,\langle{\rm SFR}\rangle_{\rm pop}\) with the lognormal means \(\langle{\rm SFR}\rangle = 10^{\mu + \log M_*} e^{(\sigma\ln 10)^2/2}\); \(\rho_{\rm SFR} \propto 10^{\rm ssfr\_ms\_norm}\) exactly and SF+Q partition it exactly (both tested). Data: the cosmic star-formation history [MadauDickinson2014].

• oiilf — the z-resolved [OII] luminosity function through the Kennicutt-like calibration [Kennicutt1998] \(\log L_{\rm [OII]} = {\rm loii\_norm} + \log{\rm SFR}\) on the ZM15 SHMR + main sequence, with kernel width \(\sqrt{\sigma_{\rm MS}^2 + \sigma_{\rm [OII]}^2 + (1+{\rm slope})^2\sigma_{M_*}^2}\) and the star-forming fraction as weight (centrals-only v1). Data: the [OII] LFs of [Comparat2015OII]. A designed degeneracy is made explicit and testable: the loii_norm and ssfr_ms_norm Jacobian columns are identical (both shift the same abscissa) — only SFR-side data (sfrd, ssfr) break it.

Radio-loud jets

The fundamental plane describes the radiatively efficient population. The jetted HERG/LERG population [BestHeckman2012] is a second rlf component from all central black holes — deliberately not tied to the ERDF-active fraction:

\[f_{\rm loud}(M_{\rm BH}) = \min\!\big[f_{\rm loud,0}\, 10^{\beta_{\rm loud}(\log M_{\rm BH} - 8)},\, 1\big], \qquad \log L_R^{\rm jet} \sim \mathcal{N}\big(b_{\rm jet} + (\log M_{\rm BH} - 8),\, \sigma_{\rm jet}\big).\]

Consequently the \(\Phi_R \propto 10^{f_{\rm ERDF}}\) amplitude identity of the FP-only rlf breaks by design once jets are on — tested in both directions (exact with \(f_{\rm loud,0} = 0\); a deficit with jets on).

The infrared AGN luminosity function

ilf maps the Powell chain’s bolometric output to 6 μm, \(L_{\rm IR} = 10^{\rm agn\_bc\_ir}\,L_{\rm bol}\) ([Hopkins2007]-style correction), with no obscuration suppression: the IR LF is obscuration-robust by construction (the agn_fabs Jacobian column is exactly zero — tested), so together with the f_abs-dimmed soft-X-ray XLF it is the cross-band consistency check of the obscured fraction, the WISE-selected halo statistics [Donoso2014], [Petter2023] probe from the clustering side. Noise: IRSurvey (WISE/SPHEREx-like, νL_ν(6 μm) completeness limit).

Wave 4 — galaxy morphology

The last roadmap topic with no code at all. Vector 102 → 106 (WAVE4_MORPHOLOGY, sector morphology); the TIER3_EXTENSION slice is frozen at [90:102].

The conditional early-type fraction

hod_mod.connection.morphology mirrors the ZM16 halo-quenching pattern: a Weibull early-type fraction of centrals,

\[f_{\rm early,c}(M_h) = 1 - \exp\!\left[-\left(M_h/M_{\rm morph} \right)^{\beta_{\rm morph}}\right],\]

with fiducial \(\log_{10} M_{\rm morph} = 12.5\), \(\beta_{\rm morph} = 0.8\), and a satellite boost toward early types \(f_{\rm early,s} = f_{\rm early,c} + f_{\rm morph,sat} (1 - f_{\rm early,c})\) (environmental transformation; ∈ [0, 1] by construction). ForwardModel(morph="early"|"late") weights the occupations exactly like the SF/Q split — EARLY + LATE ≡ unsplit, and the four-way SF/Q × early/late partition sums exactly to the unsplit sample (both tested to 1e-12; morphology and quenching selections are treated as independent at fixed M_h, a documented simplification).

The f_early observable

The roadmap’s cheap route: one Euclid-VIS-like datum \(f_{\rm early}(M_*, z)\) per (z, M*) cell — the occupation-weighted mean early-type fraction of the cell’s sample (Tier2Forecast(include_morph=True) / --include-morph; default ON in Tier3Forecast). Noise is binomial counting over the cell (negligible for wide cells) plus a morphological-calibration floor (SpectroSurvey.fmorph_err = 0.02 absolute). Data: Euclid VIS morphologies, COSMOS-Web/HSC at higher z, group-catalogue morphology–halo-mass trends [Yang2007groups], [Tinker2021groups].

The black-hole–bulge coupling

mbh_bt_slope inserts the bulge proxy into the Powell chain (\(M_{\rm BH} \propto (B/T \cdot M_*)\)-like coevolution [Yang2019BHbulge]), with \(B/T\) proxied by the mean early-type fraction of the halo:

\[\langle \log M_{\rm BH} \rangle \mathrel{+}= {\rm mbh\_bt\_slope} \cdot \log_{10}\!\big[f_{\rm early,c}(M_h) + 10^{-4}\big] .\]

At the mbh_bt_slope = 0 fiducial the chain is exactly unchanged (zero morphology response of the XLF — tested), while the mbh_bt_slope Jacobian column is live at the fiducial and, off-fiducial, routes log10_M_morph/beta_morph into the XLF, radio and IR LFs — morphology becomes testable through the AGN sector, as the roadmap proposed. (The assembly-bias hook — ordering B/T by formation time through the cosmology-dependent c(M) — remains a documented refinement.)

Observables and noise: the full menu

Opt-in observables added by the three waves

Name	Physical content	Block / flag	Noise	Completeness cut
rlf	5 GHz AGN radio LF (FP + jets)	shell / --include-radio	Poisson over the radio footprint	νL_ν ≥ L_lim(z_hi), LoTSS-like
ilf	6 μm AGN IR LF	shell / --include-ir	Poisson over the IR footprint	νL_ν ≥ L_lim(z_hi), WISE/SPHEREx-like
oiilf	[OII] emission-line LF	shell / --include-ssfr	Poisson over the spectroscopic volume	L ≥ 4π d_L² F_[OII],lim (DESI-like)
himf	HI mass function	hi_local (z ≤ 0.06) / --include-hi	Poisson, ALFALFA footprint	M_HI ≥ 2.36×10⁵ d_L² S_int
cl_gHI	21 cm × galaxy cross	cell / --include-hi	Knox with the effective IM recipe	—
ssfr	mean MS log sSFR of the cell	non-Q cell / --include-ssfr	0.05 dex absolute	—
f_early	early-type fraction of the cell	cell / --include-morph	binomial + 0.02 calibration floor	—
sfrd	SFR density of the cell	non-Q cell / --include-ssfr	12% relative (MD14-like)	—

All abundance-type observables (rlf, ilf, oiilf, himf, ssfr, sfrd) bypass the r_p/ℓ scale cuts; cl_gHI follows the angular cut of its cell. Rows failing a completeness limit receive σ = ∞ (zero Fisher weight) and are reported by the driver.

Exact invariants (the test contract)

Every mechanism ships with at least one exact identity in tests/test_missing_physics.py — the extension’s regression contract:

Invariant	Test
Growth ODE ≡ Carroll (ΛCDM) to <0.5%; ≡ exact quadrature to <0.2%	test_growth_ode_lcdm_matches_carroll_and_exact
\(\partial D/\partial w_0 \ne 0\), correct sign	test_growth_responds_to_dark_energy
ν suppression ≡ 1 at Σm_ν = 0 (bit-identical); small scales only	test_nu_suppression_shape
_c_dk15 ≡ core c_diemer15; COLOSSUS anchors; ∂c/∂σ8 ≠ 0 (and = 0 for Dutton14)	test_dk15_*
SF + Q ≡ unsplit (occupations, n_gal, ρ_SFR) to 1e-12	test_sfq_split_sums_to_unsplit, test_sfrd_observable
dlx_quenched/dhi_quenched touch ONLY quenched samples	test_d*_quenched_only_touches_*
rlf ≡ hard-band XLF at (ξ_RX, ξ_RM, b_R, σ_R) = (1,0,0,0), jets off	test_rlf_identity_with_xlf
Φ ∝ 10^f_ERDF exactly for FP-rlf and ilf; broken (deficit) by jets	test_rlf_amplitude_identity..., test_ilf_observable
\(C_\ell^{g{\rm HI}} \propto M_0\) exactly (∂ln C/∂log₁₀M₀ = ln 10)	test_hi_sector
SF-cut + Q-cut ≡ mixture-cut; cut → −∞ ≡ no cut	test_ssfr_cut_selection
ρ_SFR ∝ 10^{sSFR_MS} exactly	test_sfrd_observable
∂Φ_[OII]/∂loii_norm ≡ ∂Φ_[OII]/∂ssfr_ms_norm (the designed degeneracy)	test_oiilf_observable
ilf has exactly zero agn_fabs response (obscuration-robust)	test_ilf_observable
wind η₀ = 0 ≡ tier-2 sigmoid bit-identically	test_wind_loading
CAMB ratio: bounded (<10%), differentiable, moves ∂/∂n_s	test_camb_ratio_correction
tier-2 assembly: split cells partition n_gal; shell/hi_local block structure; all new rows carry noise	test_tier2_split_sfq_assembly, test_tier2_wave*_observables

Usage

Forward model (all opt-in; every default reproduces tier-2):

from hod_mod.forecast.forward_jax import ForwardModel

m = ForwardModel(
    z_eff=0.35,
    cm_relation="diemer15",        # cosmology-dependent c(nu, n_eff)
    pk_correction="camb_linear",   # linearized CAMB/EH98 ratio
    sfq="sf",                      # or "q", or None
    ssfr_cut=-10.5,                # ELG-like sSFR threshold [log10 yr^-1]
)
out = m.predict(theta, ["wp", "rlf", "ilf", "oiilf", "himf",
                        "cl_gHI", "ssfr", "sfrd"])

Tier-2 forecast with everything enabled:

JAX_PLATFORMS=cpu python -m hod_mod.scripts.forecasts.run_tier2_forecast \
    --rmin 0.1 0.5 2.5 --n-bands 6 --split-sfq \
    --include-radio --include-hi --include-ssfr --include-ir

Rebuild the CAMB-ratio table after changing the fiducial (needs CAMB):

python -m hod_mod.forecast.pk_camb_ratio

Approximations and fixed constants

Deliberate simplifications, each a candidate for promotion later:

Constant / choice	Value	Meaning and caveat
\(\sigma_{\rm Q}\)	0.5 dex	width of the quenched sSFR lognormal (fixed)
\(\sigma_{\rm jet}\)	0.7 dex	jet-luminosity scatter at fixed M_BH; \(\xi_{\rm jet} = 1\) fixed
\(\sigma_{\rm [OII]}\)	0.2 dex	extra [OII]-calibration scatter beyond the MS width
\(\sigma_{M_{\rm HI}}\)	0.35 dex	HI-mass conditional scatter of the HIMF kernel
ν suppression	first order	the −8f_ν tanh form + the CAMB-ratio derivative row; scale-dependent growth of ν is neglected
CAMB ratio	linearized	valid near the fiducial (a Fisher application); rebuild for a new fiducial; edge-clamped beyond k ≈ 28 h/Mpc
\(\alpha_w\)	flat at fiducial	the wind slope decouples exactly at \(\eta_{w,0} = 0\) (prior-bounded)
[OII]/AGN kernels	centrals only	the Powell-chain convention; satellites are a documented refinement
21 cm IM noise	effective (rN, aN)	the calibrated-recipe precedent (tSZ); a physical T_sys model is the upgrade
N_H template, ZM16 forms, V₀ = 200 km/s, z_pivot = 0.3	fixed	inherited tier-2 / wave conventions

Results

The 90-parameter production forecast (six X-ray bands, SF/Q-split cells, radio + IR + [OII] + HI observables — 177 blocks, 22 539 data rows) extends the Tier-2 forecast: nothing fixed (90 parameters) decomposition. Headline numbers at \(r_\mathrm{min} = 0.1\,h^{-1}\) Mpc:

Parameter	All 90 free	Astrophysics pinned	Degradation
\(\Omega_\mathrm{m}\)	\(1.82\times10^{-4}\)	\(8.01\times10^{-5}\)	×2.3
\(\sigma_8\)	\(2.93\times10^{-4}\)	\(1.22\times10^{-4}\)	×2.4
\(w_0\)	\(3.17\times10^{-3}\)	\(1.46\times10^{-3}\)	×2.2
\(w_a\)	\(1.36\times10^{-2}\)	\(6.22\times10^{-3}\)	×2.2
\(\sum m_\nu\) [eV]	\(3.23\times10^{-5}\)	—	—

Despite 29 additional free parameters, the 90-parameter run beats the 61-parameter tier-2 baseline (\(\sigma(\Omega_\mathrm{m}) = 2.85\times 10^{-4}\), \(\sigma(\sigma_8) = 4.39\times10^{-4}\)) by ~35%: the new per-cell observables (sSFR, SFRD, \(C_\ell^{g\mathrm{HI}}\)) and the radio/[OII]/IR/HI shells add more information than the new freedom removes. The cumulative probe attribution is dominated by galaxy grid → +lensing → +X-ray/tSZ (\(3.1 \to 2.7 \to 1.9 \times 10^{-4}\) on \(\Omega_\mathrm{m}\)), with the wave observables contributing percent-level refinements thereafter. The three analytically flat directions (eps_sn at the saturated energy-closure fiducial, alpha_w at the wind-off fiducial, dssfr_q without an sSFR cut) return inf marginals exactly as designed — they are prior-bounded, not data-constrained. The dominant residual degeneracies are the intra-sector ERDF triangle (agn_sig_bh–agn_rho–agn_sig_mstar, \(|\rho| = 1.000\)) and the gas metallicity-evolution pair (z_gas_norm–z_gas_zs, \(\rho = -1.000\)); the worst-constrained Fisher directions are combinations of z_gas_norm/kt_slope/kt_zs and the jet duo f_loud0/beta_loud. Completeness pruning drops 49 rows (high-z [OII]/IR/radio/X-ray LF bins below the flux limits and the two lowest HIMF masses). SUMMARY, npz products, per-sector information gains and the nine diagnostic figures are written to $HOD_MOD_RESULTS/tier2_forecast/ (*_nb6).

API reference

Physical survey-noise models for the tier-2 forecast.

Unlike the tier-1 effective (rN, aN) recipe (run_stage4_forecast), every noise term here is derived from explicit survey specifications:

• ShearSurvey — Euclid+LSST: shape noise σ_e²/n̄ per tomographic bin.

• CMBLensingSurvey — S4-like flat N_L^{κκ}.

• AthenaAllSky — hypothetical all-sky X-ray survey pinned by the completeness argument: F_lim(0.5–2 keV) = 2e-16 erg/s/cm² is exactly the depth that makes an L_X > 1e42 erg/s AGN sample complete to z = 1. Photon (CXB) noise for the band C_ℓ, PSF beam, and the L_lim(z) completeness check.

• SpectroSurvey — DESI/4MOST-like: pair-count + cosmic-variance errors for the per-cell w_p and ΔΣ, Poisson errors for the abundances/XLF.

All functions are numpy and evaluated once at the fiducial — consistent with fisher.fisher_matrix treating the covariance as parameter-independent. Model-unit conversions (the forecast X-ray field carries arbitrary but self-consistent units) are handled by expressing every noise as a noise-to-signal ratio against fiducial model quantities.

class hod_mod.forecast.noise.AthenaAllSky(f_lim: float = 2e-16, psf_hew: float = 5.0, f_sky: float = 0.65, n_det_lim: float = 10.0, e_mean: float = 1.63e-09, i_cxb: float = 2.6e-08, gamma_cxb: float = 1.4, shell_cxb_frac: float = 0.05)

Bases: object

Hypothetical Athena all-sky survey (completeness-pinned depth).

f_lim = 2e-16 erg/s/cm² (0.5–2 keV) makes L_X > 1e42 complete to z = 1 (l_lim()); the implied all-sky depth t·A_eff = n_det·ē/F_lim ≈ 8e7 s·cm² is the stated optimistic premise. The 5” HEW PSF matters through source detection/confusion (the XLF), not the ℓ ≤ 3000 band C_ℓ (beam ≈ 1).

band_flux_fractions(bands)

CXB energy-flux fraction per band (Γ_cxb power law), Σ over 0.5–2 = 1.

beam(ell)

Gaussian beam b_ℓ for the PSF (HEW = 2.355 σ).

e_mean: float = 1.63e-09

property exposure_area

Effective survey depth t·A_eff [s·cm²] implied by the flux limit.

f_lim: float = 2e-16

f_sky: float = 0.65

gamma_cxb: float = 1.4

i_cxb: float = 2.6e-08

l_lim(z, h, Om)

Faintest detectable L_X(0.5–2) [erg/s] at redshift z (completeness).

n_det_lim: float = 10.0

noise_cl_model(ibar_model, bands=None)

White photon-noise power in MODEL units, N_X = Ī_model²/(n_γ·s_x²).

ibar_model is the fiducial mean intensity of the modeled shell field in the forecast’s own X-ray units; s_x = shell_cxb_frac anchors that field to its physical share of the CXB, whose photon statistics set the map noise (per band when bands is given).

photon_density(bands=None)

CXB photon surface density [sr⁻¹], total or per band.

psf_hew: float = 5.0

shell_cxb_frac: float = 0.05

class hod_mod.forecast.noise.BandLFSurvey(f_sky: float, nulnu_lim: float)

Bases: object

Generic broad-band luminosity-function survey: an f_sky footprint with a νL_ν-equivalent flux limit driving the L_lim(z) completeness (the Athena/radio/IR pattern, instantiated per band: UV, opt, NIR, AGN UV/opt).

f_sky: float

l_lim(z, h, Om)

Faintest detectable νL_ν [erg/s] at redshift z.

nulnu_lim: float

class hod_mod.forecast.noise.CMBLensingSurvey(f_sky: float = 0.4, n0: float = 7e-09)

Bases: object

S4-like CMB lensing: flat N_L^{κκ} (a good approximation at L < 3000).

f_sky: float = 0.4

n0: float = 7e-09

class hod_mod.forecast.noise.HISurvey(f_sky: float = 0.16, s_int_lim: float = 0.6, z_himf: float = 0.06, f_sky_im: float = 0.1, rn_im: float = 2.0, an_im: float = 0.5)

Bases: object

Blind HI survey + 21 cm intensity mapping (ALFALFA/MIGHTEE → SKA1 era).

The HIMF uses Poisson counts with an M_HI detection limit M_lim(z) = 2.36×10⁵ d_L² S_int (the standard 21 cm mass–flux relation); the 21 cm × galaxy cross uses the calibrated effective (rN, aN) recipe (the tSZ precedent) — a CHIME/MeerKLASS-era noise-to-signal.

an_im: float = 0.5

f_sky: float = 0.16

f_sky_im: float = 0.1

mhi_lim(z, h, Om)

Faintest detectable M_HI [Msun/h] at redshift z (2.36e5 d_L² S).

rn_im: float = 2.0

s_int_lim: float = 0.6

z_himf: float = 0.06

class hod_mod.forecast.noise.IRMapSurvey(f_sky: float = 0.65, bands: tuple = (3.4, 4.9, 12.0), rn: tuple = (0.3, 0.3, 0.4), an: float = 0.5)

Bases: object

WISE/SPHEREx-like infrared intensity maps in a few μm bands.

Same effective (rn, an) noise recipe as SKASurvey; zodiacal and stellar residuals dominate, hence the larger rn at 12 μm.

an: float = 0.5

bands: tuple = (3.4, 4.9, 12.0)

f_sky: float = 0.65

noise_cl(ell, cl_band, i)

Effective noise power for band i against its fiducial auto.

rn: tuple = (0.3, 0.3, 0.4)

class hod_mod.forecast.noise.IRSurvey(f_sky: float = 0.65, nulnu_lim: float = 1e-14)

Bases: object

WISE/SPHEREx-like all-sky infrared AGN survey.

nulnu_lim is the 6 μm νL_ν-equivalent detection threshold — the IR analogue of the Athena/radio flux limits, driving the L_lim(z) completeness of the AGN IR luminosity function.

f_sky: float = 0.65

l_lim(z, h, Om)

Faintest detectable 6 μm νL_ν [erg/s] at redshift z.

nulnu_lim: float = 1e-14

class hod_mod.forecast.noise.RadioSurvey(f_sky: float = 0.13, nulnu_lim: float = 3e-22)

Bases: object

LOFAR/LoTSS-like radio continuum survey with redshift counterparts.

nulnu_lim is the νL_ν detection threshold expressed at the fundamental plane’s 5 GHz reference (a 144 MHz flux limit of ~0.8 mJy scaled with a ν^{-0.7} synchrotron spectrum) — the radio analogue of the Athena F_lim, driving the L_lim(z) completeness of the radio luminosity function.

f_sky: float = 0.13

l_lim(z, h, Om)

Faintest detectable 5 GHz νL_ν [erg/s] at redshift z.

nulnu_lim: float = 3e-22

class hod_mod.forecast.noise.SKASurvey(f_sky: float = 0.5, bands: tuple = (0.95, 1.4, 3.0), rn: tuple = (0.2, 0.2, 0.3), an: float = 0.5)

Bases: object

SKA-like radio continuum intensity maps in a few GHz bands.

Map noise follows the calibrated effective recipe (the tSZ / 21 cm IM precedent): N_b(ℓ) = rn_b (ℓ/100)^an · C_b(ℓ) against the fiducial band auto — thermal noise, calibration residuals and bright-source masking are absorbed into (rn, an), the documented upgrade path being a physical T_sys/confusion model.

an: float = 0.5

bands: tuple = (0.95, 1.4, 3.0)

f_sky: float = 0.5

noise_cl(ell, cl_band, i)

Effective noise power for band i against its fiducial auto.

rn: tuple = (0.2, 0.2, 0.3)

class hod_mod.forecast.noise.ShearSurvey(n_eff: float = 30.0, sigma_e: float = 0.26, f_sky: float = 0.5)

Bases: object

Euclid+LSST combined shear: 30 sources/arcmin², f_sky = 0.5.

f_sky: float = 0.5

n_bin_sr(n_bins=1)

Source density per tomographic bin [sr⁻¹] (equal-number split).

n_eff: float = 30.0

noise_cl(n_bins=1)

Shape-noise power N_κ = σ_e²/n̄ per bin [sr].

sigma_e: float = 0.26

class hod_mod.forecast.noise.SpectroSurvey(f_sky: float = 0.5, f_cv0: float = 0.005, pi_max: float = 100.0, ssfr_err: float = 0.05, sfrd_rel: float = 0.12, foii_lim: float = 1e-16, fha_lim: float = 1e-16, fmorph_err: float = 0.02, size_err: float = 0.02, fmorph_agn_err: float = 0.05, mstar_lim0: float = None, mstar_lim_slope: float = 0.0)

Bases: object

DESI/4MOST-like spectroscopy over the shear footprint.

Tier-3 stellar-mass completeness: with mstar_lim0 set, a (z, M*) cell is complete iff its lower edge satisfies m_lo >= mstar_lim0 + mstar_lim_slope * z_hi (a magnitude-limited selection rises roughly linearly in log M* with z). None (default) keeps the tier-2 behaviour: complete everywhere.

complete_for(m_lo, z_hi)

Is a cell with lower M* edge m_lo complete to z_hi?

cv_rel(volume)

Relative cosmic-variance floor for a cell of volume (Mpc/h)³.

f_cv0: float = 0.005

f_sky: float = 0.5

fha_lim: float = 1e-16

fmorph_agn_err: float = 0.05

fmorph_err: float = 0.02

foii_lim: float = 1e-16

lha_lim(z, h, Om)

Faintest detectable L_Hα [erg/s] at redshift z.

loii_lim(z, h, Om)

Faintest detectable L_[OII] [erg/s] at redshift z.

mstar_lim0: float = None

mstar_lim_slope: float = 0.0

pi_max: float = 100.0

sfrd_rel: float = 0.12

size_err: float = 0.02

ssfr_err: float = 0.05

hod_mod.forecast.noise.athena_noise_cl_model(ath: AthenaAllSky, bands, z_eff, chi1, chi2, h)

White photon-noise power for a shell’s band C_ℓ^XX in MODEL units, (Nb,).

The forecast’s X-ray amplitudes are (L/10^45 erg/s) comoving emissivity densities, so the model→physical intensity conversion per unit comoving depth is conv = 10^45/(4π(1+z)⁴ (Mpc/h → cm)²). Anchoring the shell field with a thin-shell window of depth Δχ = χ₂−χ₁ gives the (signal-independent) white-noise level

\[N^{\rm model}_b = \frac{I_{\rm CXB,b}\,\bar e_b} {t A_{\rm eff}\; {\rm conv}^2\, \Delta\chi} .\]

Divide by ath.beam(ell)**2 for the beam-deconvolved noise (≈1 here).

hod_mod.forecast.noise.band_mean_photon_energy(bands, gamma_cxb=1.4)

Mean photon energy per band [erg] for an E^{-Γ} photon spectrum.

hod_mod.forecast.noise.chi_of(z, h, Om)

Comoving distance χ(z) [Mpc/h] (numpy scalar/array).

hod_mod.forecast.noise.delta_sigma_noise(rp, ds, z_l, ngal, volume, h, Om, zs_grid, nz_src, shear: ShearSurvey, spectro: SpectroSurvey, dz_buffer=0.1)

Absolute σ on ΔΣ(r_p) [h M⊙/pc²]: lensing shape noise (+ CV floor).

σ_shape = σ_e ⟨Σ_crit⁻¹⟩⁻¹ / √(n_src^eff · A_ann · N_lens), with the source density and ⟨Σ_crit⁻¹⟩ restricted to z_s > z_l + dz_buffer — high-z lens cells become noise-dominated automatically as the background empties.

hod_mod.forecast.noise.knox_auto(ell, cl, noise_cl, f_sky)

Absolute Gaussian σ on an auto spectrum: √(2/N_modes)·(C+N).

hod_mod.forecast.noise.knox_cross(ell, cl, cl_a, noise_a, cl_b, noise_b, f_sky)

Absolute Gaussian σ on a cross spectrum A×B.

hod_mod.forecast.noise.n2d_of(ngal, chi1, chi2)

Projected galaxy density [sr⁻¹] of a shell: n̄_g·(χ₂³−χ₁³)/3.

hod_mod.forecast.noise.n_modes(ell, f_sky)

Gaussian mode count per log-spaced ℓ bin (the stage4 convention).

hod_mod.forecast.noise.poisson_relerr(density, volume)

Relative Poisson error 1/√N for a count density × volume.

hod_mod.forecast.noise.shell_volume(z1, z2, h, Om, f_sky)

Comoving volume of the z-shell [(Mpc/h)³].

hod_mod.forecast.noise.wgp_noise(rp, wgp, z_l, ngal, volume, h, Om, shear: ShearSurvey, spectro: SpectroSurvey)

Absolute σ on w_g+(r_p) [Mpc/h]: shape noise per annulus + CV floor.

The delta_sigma_noise geometry WITHOUT the Σ_crit lensing weight — the intrinsic-alignment estimator correlates density tracers with the SHAPES of the same (or an overlapping) sample, so the noise per r_p bin is σ_e/√(n_shape·A_ann·N_dens) projected over 2π_max. An effective recipe (no IA–clustering cross-covariance), the tSZ/IM documentation precedent.

hod_mod.forecast.noise.wp_pair_sigma(rp, wp, ngal, volume, survey: SpectroSurvey)

Absolute σ on w_p(r_p) [Mpc/h]: pair-count Poisson + cosmic variance.

N_pair = ½ n̄² V · π(r₂²−r₁²) · 2π_max per r_p bin; σ_shot = 2π_max (1 + w_p/2π_max) / √N_pair.

hod_mod.forecast.noise.xlf_relerr(phi, volume, dloglx=0.5)

Relative Poisson σ per XLF point: N = Φ·V·Δlog L_X.

CAMB-quality linear P(k) for the forecast: a linearized EH98-ratio table.

The differentiable forecast uses the analytic EH98 shape, which deviates from CAMB by up to ~4% around the BAO scale. For a Fisher application the requirement is precise and modest: the spectrum and its first derivatives must be CAMB-accurate near the fiducial. That is exactly what a first-order expansion of the log shape-ratio delivers:

\[\ln R(k;\theta) \simeq \ln R_0(k) + \sum_i \frac{\partial \ln R}{\partial \theta_i}(k)\, (\theta_i - \theta_i^{\rm fid}), \qquad R = \frac{P^{\rm shape}_{\rm CAMB}}{P^{\rm shape}_{\rm EH98}}\]

with both shapes pivot-normalised at k = 0.05 h/Mpc (the σ8 anchoring makes the overall normalisation irrelevant). The EH98 side includes the forecast’s ν suppression, so the Σm_ν derivative row corrects the residual of the first-order tanh form against CAMB. Eleven CAMB evaluations (fiducial + central/forward differences over h, Ω_b, Ω_m, n_s, Σm_ν) are distilled once into hod_mod/data/pk_ratio/camb_eh98_ratio.npz (the apec_bands pattern); the JAX side is two jnp.interp calls and an exp.

Beyond the tabulated k range the log-ratio is edge-clamped (deep power-law tail where EH98 is accurate). w0/wa need no rows: they do not change the z = 0 transfer-function shape (only growth/geometry, handled elsewhere).

hod_mod.forecast.pk_camb_ratio.apply_ratio(pk, k, theta: dict, tab: dict)

Multiply an EH98 shape by the linearized CAMB ratio (differentiable).

hod_mod.forecast.pk_camb_ratio.build(npz_path: str = '/home/docs/checkouts/readthedocs.org/user_builds/hod-mod/checkouts/stable/hod_mod/data/pk_ratio/camb_eh98_ratio.npz', verbose: bool = True) → dict

Build the ratio table with CAMB (once; ~1 min) and cache it as npz.

hod_mod.forecast.pk_camb_ratio.load(npz_path: str = '/home/docs/checkouts/readthedocs.org/user_builds/hod-mod/checkouts/stable/hod_mod/data/pk_ratio/camb_eh98_ratio.npz') → dict

Load the table as jnp arrays; raise with build instructions if absent.

Tier-2 forecast assembly: the (z, M*) cell grid + AGN + global lensing blocks.

Unlike TomographicForecast (per-bin HOD copies), the tier-2 design uses ONE shared global parameter vector — redshift evolution is carried by the explicit *_zs slope parameters inside each ForwardModel (_theta_eff) — so the global vector is simply PARAM_NAMES (61 entries) and the Jacobian is the row-stack of independent per-block Jacobians:

• cell blocks — volume-limited (z, M*) samples: Δz = 0.1 shells × 0.2-dex M* bins, each predicting (wp, ds, cl_gX×bands, cl_gy, cl_gkCMB, n_gal) (the per-bin n_gal IS the stellar-mass-function datum, so smf is not a separate observable);

• shell blocks — per-Δz observables that do not depend on the M* split: the soft-band AGN XLF and the per-band X-ray auto cl_XX;

• global block — tomographic cosmic shear (cl_kk pairs), CMB lensing and their cross (cl_kCMB, cl_shear_kCMB);

• wp_agn blocks — projected clustering of complete soft-L_X-selected AGN samples in 0.5-dex bins at a few redshifts.

Each block gets its own jax.jacfwd (never one monolithic Jacobian) and a per-block npz cache, so re-runs are incremental. Noise is physical (hod_mod.forecast.noise): pair counts + cosmic variance for the projected statistics, shape noise for lensing, CXB photon noise + the completeness-pinned Athena spec for X-rays, Poisson counts for the XLF.

class hod_mod.forecast.tier2.Tier2Forecast(z_edges=None, mstar_edges=None, n_bands=6, n_shear_bins=5, agn_lx_bins=None, agn_z_centers=(0.1, 0.3, 0.5, 0.7, 0.9), shear=None, cmbl=None, athena=None, spectro=None, tsz=(0.25, 0.9, 0.3), split_sfq=False, include_radio=False, include_hi=False, include_ssfr=False, include_ir=False, include_morph=False, radio=None, hi=None, ir=None, **model_kw)

Bases: object

Shared-vector multi-block tier-2 forecast (see module docstring).

Parameters:

• z_edges, mstar_edges (array) – Cell grid: Δz = 0.1 shells over 0 < z < 1 × 0.2-dex bins over 10.0 ≤ log10 M* ≤ 11.6 by default (80 cells).

• n_bands (int) – X-ray energy bands over 0.5–2 keV: 1 (broad), 6 (default) or 15 (the validated production 100 eV grid); any explicit list of (emin, emax) pairs is also accepted.

• n_shear_bins (int) – Tomographic shear source bins (equal-number Smail split).

• agn_lx_bins, agn_z_centers – Soft-band L_X bins (0.5 dex, complete above 1e42) and the redshifts of the AGN clustering samples (Δz = 0.2 windows).

• shear, cmbl, athena, spectro (noise.py survey dataclasses)

• tsz ((rN, aN, f_sky)) – The calibrated stage-4 effective tSZ noise recipe (kept in v1).

• split_sfq (bool) – Split every (z, M*) cell into star-forming and quiescent samples (ZM16 Weibull quenching; missing-physics extension). Doubles the cell blocks; the shell X-ray observables stay unsplit (total gas).

• model_kw (forwarded to every ForwardModel (grids etc.).)

data_and_jacobian(fid, cache_dir=None, verbose=True)

Assemble (d0, J, meta) block by block; per-block npz caching.

meta is a dict of per-row arrays: block, kind, zeff, z_lo, z_hi, m_lo, m_hi, obs, x, sub (band / L_X-bin / shear-pair sub-index).

fiducial()

noise_sigma(fid, d0, meta, verbose=True)

Per-row absolute Gaussian σ from the physical survey noise models.

Completeness: XLF / wp_agn rows whose L_X bin dips below the Athena detection limit L_lim(z_hi) get σ = inf (zero weight) and are reported.

precompute_blocks(fid, cache_dir, jobs=1, verbose=True, max_tasks_per_child=4)

Fill the per-block npz cache with jobs worker processes.

Spawned workers rebuild this exact forecast from self._ctor_kwargs (matching the parent’s x64 mode) and write each block atomically (tempfile + os.replace), so a subsequent serial data_and_jacobian() assembles bit-identical results from the cache — the parallel == serial invariant. Returns the labels that were missing on entry.

max_tasks_per_child bounds worker memory: a worker’s JAX compilation cache grows with every distinct block shape it touches (several GB after a handful of blocks), and unbounded workers OOM the host on a full tier-3 run. Implemented as BATCHED POOLS — a fresh executor per chunk of jobs × max_tasks_per_child blocks — rather than the executor’s own max_tasks_per_child, whose worker-respawn path deadlocks on CPython 3.11 (observed: pool alive, all workers exited, no respawn). Pool teardown between batches frees the caches identically, at one forecast rebuild (~seconds) per worker per batch.

prior(add_planck=False, fix=('log10DC',))

Regularizing prior; the retired log10DC is pinned by default (agn_emission=”powell” removes it from the emissivity entirely).

scale_cut_mask(meta, rmin)

Per-block mask: r_p > rmin (projected), ℓ < χ(z_eff)/rmin (angular).

Tier-3 forecast assembly: multi-wavelength maps, band LFs, z < 2 × M* > 10⁹.

Extends Tier2Forecast (whose block/cache/noise machinery is reused unchanged through its extension hooks) with

• a coarse exploratory grid: Δz = 0.2 shells over 0 < z < 2 × 0.2-dex M* bins over 9.0 ≤ log10 M* ≤ 11.6 (130 cells, ×2 with the SF/Q split), each cell on the extended mass grid (log10m_min = 8.5);

• radio and IR intensity maps (SKA-like GHz bands; WISE/SPHEREx-like μm bands): per-cell galaxy crosses cl_gR/cl_gI, per-shell autos cl_RR/cl_II and AGN crosses cl_aR/cl_aI/cl_ag;

• galaxy band LFs (uvlf, optlf, nirlf, the half Hα LF) and AGN UV/optical LFs (qlf_uv, qlf_opt) per shell;

• a per-shell wide-M* SFRD(z) block (the Madau–Dickinson measurement);

• the four extras: tSZ auto cl_yy, 21 cm auto cl_HIHI, X-ray cluster counts ncl (per shell) and AGN lensing ds_agn (per AGN z-block).

Two-tier galaxy completeness: the wide spectroscopic survey carries a stellar-mass limit log10 M*_lim(z) = mstar_lim0 + mstar_lim_slope·z; cells below it fall back to a small deep field (spectro_deep), and cells complete in neither tier are not built (recorded in skipped_cells).

class hod_mod.forecast.tier3.Tier3Forecast(z_edges=None, mstar_edges=None, agn_z_centers=(0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 1.7, 1.9), spectro=None, spectro_deep=None, ska=None, irmap=None, lf_uv=None, lf_opt=None, lf_nir=None, lf_quv=None, lf_qopt=None, include_maps=True, include_bandlfs=True, include_extras=True, cell_log10m_min=8.5, cell_n_m=256, **kw)

Bases: Tier2Forecast

Tier-3 multi-wavelength forecast (see module docstring).

Beyond Tier2Forecast (whose wave flags default ON here):

Parameters:

• spectro, spectro_deep (noise.SpectroSurvey) – Wide tier (f_sky = 0.5, M*-complete to 10^{9+z}) and deep tier (f_sky = 0.004, flat at 10^9).

• ska, irmap (noise.SKASurvey, noise.IRMapSurvey) – Radio/IR intensity-map surveys (bands + effective noise recipes).

• lf_uv, lf_opt, lf_nir, lf_quv, lf_qopt (noise.BandLFSurvey) – Band-LF footprints and νL_ν flux limits (GALEX/Rubin/WISE-like and quasar-survey-like defaults).

• include_maps, include_bandlfs, include_extras (bool) – Toggle the three tier-3 observable families.

• cell_log10m_min, cell_n_m (float, int) – Mass grid of the cell/sfrd models (the 8.5 floor resolves M* = 10^9 occupations; n_m = 256 converges to <1e-3).

noise_sigma(fid, d0, meta, verbose=True)

Tier-2 noise for the inherited rows, then the tier-3 families.

Band-integrated APEC tables + AGN spectral templates for the tier-2 X-ray layer.

The differentiable forecast cannot trace soxs/AtomDB, so this module distills hod_mod.gas.cooling.ApecCoolingTable (the machinery behind the production fit_xray_joint_bands energy-band fit) into static log10 Λ_b(log10 T, log10 Z) grids, cached as one npz per band set under hod_mod/data/apec_bands/. The forecast then only ever does JAX bilinear interpolation on those arrays — soxs is needed once per new band configuration.

Also provides the Morrison & McCammon (1983) ISM photoelectric cross-section and per-band transmission templates for the obscured-AGN fraction (agn_fabs).

hod_mod.forecast.apec_bands.band_tables(bands, n_T=48, T_min=0.05, T_max=60.0, n_Z=13, Z_min=0.05, Z_max=3.0)

log10 Λ(log10 T, log10 Z) tables for bands + the 0.5–2 broad band.

Returns a dict of numpy arrays: lt (n_T,), lz (n_Z,), tables (n_bands+1, n_T, n_Z) with the broad band LAST, and edges. Cached to an npz keyed on the full configuration; building needs soxs once.

hod_mod.forecast.apec_bands.band_transmission(bands, nh=1e+22, gamma=1.8, n_e=256)

Energy-flux-weighted transmission of a Γ power law through N_H, per band.

A static template for the obscured-AGN fraction: the intra-band Γ dependence is second order, so the template is evaluated at the fiducial Γ and the free agn_gamma only moves the band fractions, not the transmissions.

hod_mod.forecast.apec_bands.mm83_sigma(e_kev)

MM83 ISM photoelectric cross-section σ(E) [cm² per H atom], 0.03–10 keV.