Galaxy × Thermal Sunyaev–Zel’dovich (tSZ)

The thermal SZ effect measures the line-of-sight integral of the electron pressure through the hot intracluster and circumgalactic medium. Cross-correlating a galaxy sample with a Compton-\(y\) map probes the pressure–halo connection of the galaxies’ host haloes.

This pipeline reuses the same halo-model engine as the galaxy clustering and galaxy × X-ray pipelines: HaloModelCrossSpectra wraps an existing FullHaloModelPrediction (reusing its cached halo mass function, bias, linear power spectrum and dark-matter profile FT) and adds an electron-pressure field (PressureProfileA10, Arnaud+2010, or the DPM variant).

The model

The galaxy × Compton-\(y\) cross-power spectrum has the usual 1-halo + 2-halo decomposition,

\[P_{g,y}(k, z) = P_{g,y}^{1h}(k,z) + P_{g,y}^{2h}(k,z),\]

where the galaxy leg is the occupation \(\langle N_g \rangle(M)\) weighted by the halo number density, and the \(y\) leg is the Fourier transform of the electron pressure profile \(\tilde{y}(k|M, z)\) (dimensionless Compton-\(y\)). The projected and angular observables follow by:

  • projected_gy() — the stacked Compton-\(y\) profile \(\Sigma_y(r_p)\) via Abel projection;

  • angular_cl_gy() — the angular cross-spectrum \(C_\ell^{g,y}\) via the Limber approximation.

Worked example

import numpy as np
from hod_mod.core.power_spectrum import LinearPowerSpectrum
from hod_mod.core.halo_mass_function import make_hmf
from hod_mod.core.halo_profiles import HaloProfile
from hod_mod.connection.hod import MoreHODModel
from hod_mod.observables.clustering import FullHaloModelPrediction
from hod_mod.observables.cross_spectra import HaloModelCrossSpectra
from hod_mod.gas import PressureProfileA10

theta = {"h": 0.6774, "Omega_m": 0.3089, "Omega_b": 0.0486,
         "n_s": 0.9667, "sigma8": 0.8159}

pk_lin = LinearPowerSpectrum()
hmf    = make_hmf("tinker08", pk_func=pk_lin.pk_linear)
hod    = MoreHODModel(hmf, hmf.bias)
hp     = HaloProfile()

fhmp   = FullHaloModelPrediction(pk_lin, hod, hp)
pp     = PressureProfileA10(r_max_over_r500c=5.0, n_gl=150)
cross  = HaloModelCrossSpectra(fhmp, pressure_profile=pp)

rp      = np.logspace(-1, 1.3, 20)           # Mpc/h
z       = 0.3
params  = hod.default_params()
sigma_y = cross.projected_gy(rp, z, theta, params)

# angular C_ell^{g,y} integrates over the galaxy redshift distribution n(z)
ell     = np.logspace(2, 4, 30)
z_arr   = np.linspace(0.2, 0.5, 16)
nz_g    = np.exp(-0.5 * ((z_arr - 0.3) / 0.05) ** 2)
cl_gy   = cross.angular_cl_gy(ell, z_arr, nz_g, theta, params)

The validation figures (A10 pressure profile, \(P_{g,y}(k)\) decomposition, \(\Sigma_y(r_p)\), \(C_\ell^{g,y}\)) are produced by:

hod-mod validate sz-xray

References

  • Arnaud et al. 2010, arXiv:0910.1234 — A10 generalized-NFW pressure profile.

  • Amodeo et al. 2021, arXiv:2009.05557 — ACT × BOSS CMASS/LOWZ stacked tSZ.

  • Pandey et al. 2025, arXiv:2506.07432 — DES Y3 shear × ACT DR6 tSZ.

API

class hod_mod.observables.cross_spectra.HaloModelCrossSpectra(fhmp, pressure_profile: PressureProfileA10 | PressureProfileDPM | None = None, density_profile: GasDensityDPM | None = None, metallicity_profile=None, cooling_function=None, agn_model=None, ecf_gas_table=None, ecf_agn=None)[source]

Bases: object

Halo model galaxy × gas cross-power spectra and projected observables.

Wraps a FullHaloModelPrediction and reuses its static cache (HMF, bias, DM profile FT, linear P(k)), adding gas-profile Fourier transforms for the y-field (tSZ) and X-ray emissivity field.

Parameters:

  • fhmp (FullHaloModelPrediction) – Already-instantiated prediction object whose static cache is reused.

  • pressure_profile (PressureProfileA10, optional) – Electron pressure profile for tSZ. If None, tSZ methods raise.

  • density_profile (GasDensityDPM, optional) – Electron density profile for X-ray. If None, X-ray methods raise.

angular_cl_XX(ell_arr: ndarray, z_arr: ndarray, nz_X: ndarray, theta_cosmo: dict, psf_fwhm_arcsec: float | None = None, psf_king_theta_c_arcsec: float | None = None, psf_king_alpha: float = 1.5, return_components: bool = False, n_workers: int = 1, beta_gas: float | None = None, beta_pressure: float | None = None) → ndarray | dict[source]

Angular auto-power spectrum C_ℓ^{X,X} of the total X-ray emission.

Includes the 1-halo and 2-halo gas–AGN cross-terms that vanish in angular_cl_gX() (which is exact and linear in X). See _pk_tables_XX() for the underlying P_{X,X}(k) decomposition.

\[C_\ell^{X,X} = \int_0^{\chi_{\max}} \frac{\mathrm{d}\chi}{\chi^2} W_X(\chi)^2\,P_{X,X}\!\left(k=\frac{\ell+\tfrac{1}{2}}{\chi}, z(\chi)\right)\]
Parameters:

  • ell_arr ((Nell,) angular multipoles)

  • z_arr ((Nz,) redshift array bracketing the X-ray source distribution)

  • nz_X ((Nz,) X-ray window function (e.g. emissivity-weighted dV/dz,) – or a matched source dN/dz). Normalized internally like nz_g in angular_cl_gX(), but appears squared in the Limber integral since both legs of the correlation are the X-ray field.

  • return_components (bool) – If True return {"total", "gas_gas", "cross", "agn_agn"}.

Returns:

cl_XX ((Nell,) or dict when return_components=True)

angular_cl_gX(ell_arr: ndarray, z_arr: ndarray, nz_g: ndarray, theta_cosmo: dict, hod_params: dict, psf_fwhm_arcsec: float | None = None, psf_king_theta_c_arcsec: float | None = None, psf_king_alpha: float = 1.5, return_components: bool = False, n_workers: int = 1, beta_gas: float | None = None, beta_pressure: float | None = None, agn_kwargs: dict | None = None, x_uk_override: list | ndarray | None = None) → ndarray | dict[source]

Angular cross-power spectrum C_ℓ^{g,X} via the Limber approximation.

Identical structure to angular_cl_gy() but for the X-ray emissivity field (DPM Model). The returned spectrum has units of the emissivity power spectrum [(Mpc/h)³ cm⁻⁶] divided by χ² [(Mpc/h)²], giving [(Mpc/h) cm⁻⁶].

\[C_\ell^{g,X} = \int_0^{\chi_{\max}} \frac{\mathrm{d}\chi}{\chi^2} W_g(\chi)\,P_{g,X}\!\left(k=\frac{\ell+\tfrac{1}{2}}{\chi}, z(\chi)\right)\]
Parameters:

  • ell_arr ((Nell,) angular multipoles)

  • z_arr ((Nz,) redshift array bracketing the galaxy distribution)

  • nz_g ((Nz,) dN/dz of the galaxy sample (normalized internally))

  • theta_cosmo (cosmological parameters)

  • hod_params (HOD parameters)

  • psf_fwhm_arcsec (float | None) – eROSITA PSF FWHM [arcsec]. If given, multiply C_ℓ by the Gaussian PSF window B_ℓ = exp(−ℓ²σ²/2) (single-field convolution). Use 30.0 for the eROSITA soft X-ray PSF.

  • psf_king_theta_c_arcsec (float | None) – King core radius [arcsec] for the analytic PSF window. If given, multiply C_ℓ by B_ℓ = exp(−ℓ θ_c) (α=3/2) or the general Bessel-K form. Fitted to eROSITA TM CalDB on-axis: 8.64”. Cannot be used together with psf_fwhm_arcsec.

  • psf_king_alpha (float) – King slope for the analytic PSF window. Default 1.5.

  • return_components (bool) – If True return a dict {"total", "gas", "agn"} instead of the total array.

  • n_workers (int) – Number of threads for parallel evaluation of _pk_tables_gX at each redshift. -1 (default) uses os.cpu_count(). Set to 1 to disable parallelism. The z-points are independent, so thread-based parallelism is safe because JAX releases the GIL during computation.

Returns:

cl_gX ((Nell,) [(Mpc/h) cm⁻⁶] or dict when return_components=True)

angular_cl_gy(ell_arr: ndarray, z_arr: ndarray, nz_g: ndarray, theta_cosmo: dict, hod_params: dict, psf_fwhm_arcsec: float | None = None, psf_king_theta_c_arcsec: float | None = None, psf_king_alpha: float = 1.5) → ndarray[source]

Angular cross-power spectrum C_ℓ^{g,y} via the Limber approximation.

Under Limber (Limber 1953; LoVerde & Afshordi 2008):

\[C_\ell^{g,y} = \int_0^{\chi_{\max}} \frac{\mathrm{d}\chi}{\chi^2} W_g(\chi)\,P_{g,y}\!\left(k=\frac{\ell+\tfrac{1}{2}}{\chi}, z(\chi)\right)\]

where \(W_g(\chi) = \mathrm{d}N_g/\mathrm{d}\chi\) (normalized). The y-field window is unity (already a LOS integral).

Parameters:

  • ell_arr ((Nell,) angular multipoles)

  • z_arr ((Nz,) redshift array for n(z) [must bracket the galaxy distribution])

  • nz_g ((Nz,) dN/dz of the galaxy sample (will be normalized internally))

  • theta_cosmo (cosmological parameters)

  • hod_params (HOD parameters)

  • psf_fwhm_arcsec (float | None) – If given, multiply C_ℓ by the Gaussian PSF window B_ℓ (single field).

  • psf_king_theta_c_arcsec (float | None) – If given, multiply C_ℓ by the analytic King-profile PSF window B_ℓ = exp(−ℓ θ_c) for α=3/2, or the general Bessel-K form. Cannot be used together with psf_fwhm_arcsec.

  • psf_king_alpha (float) – King slope for the analytic PSF window. Default 1.5.

Returns:

cl_gy ((Nell,) [(Mpc/h)²] (dimensionless after h-unit cancellation with χ²))

emissivity_xuk_bands_per_z(z_arr, theta_cosmo, hod_params, cooling_fns)[source]

Multi-band version of emissivity_xuk_per_z() for the energy bands.

Computes the per-z, per-band raw emissivity FT X̃_b(k|M)/Λ_ref,b in one batched FT per z (GasDensityDPM.emissivity_full_uk_bands). Each band is divided by its OWN reference cooling Λ_ref,b = cooling_fns[b](1 keV, 0.3 Z⊙) so band b feeds angular_cl_gX(x_uk_override=...) exactly like the broad-band path (which divides by the broad-band Λ_ref).

Returns a list over z; each entry is a list over bands of (Nk, NM) arrays (so out[iz][b] is a valid x_uk_override for band b at z-index iz). self._dp/_pp/_mp must be set (the cooling tables are passed explicitly, not taken from self._cooling_fn).

emissivity_xuk_per_z(z_arr, theta_cosmo, hod_params)[source]

Precompute the per-z raw emissivity FT X̃(k|M)/Λ_ref for the emulator.

Returns a list (one entry per z) of (Nk, NM) arrays — exactly what angular_cl_gX() consumes via x_uk_override. This isolates the expensive full-APEC FT (emissivity_full_uk, ~1.3 s/z) so a joint fit can cache it on a (p2, r_max) grid and skip it at evaluation time. The cached value is at this instance’s self._dp._ne_03; the caller scales by (n_e,0.3 / n_e,0.3_ref)² before passing it back.

Requires the full-APEC profiles (_dp/_pp/_mp/_cooling_fn).

projected_gX(rp_arr: ndarray, z: float, theta_cosmo: dict, hod_params: dict) → ndarray[source]

Projected galaxy × X-ray emissivity w_{g,X}(r_p).

Same Abel projection as projected_gy() but for P_{g,X}(k). Units: (Mpc/h)³ cm⁻⁶ × (h/Mpc)² = (Mpc/h) cm⁻⁶. Multiply by the effective cooling function Λ_eff [erg cm³ s⁻¹] to compare to surface-brightness data.

Parameters:

  • rp_arr ((NR,) [Mpc/h])

  • z, theta_cosmo, hod_params (as for _pk_tables_gX)

Returns:

wgX ((NR,) [(Mpc/h) cm⁻⁶])

projected_gy(rp_arr: ndarray, z: float, theta_cosmo: dict, hod_params: dict) → ndarray[source]

Projected galaxy × y signal Σ_y(r_p) [dimensionless Compton-y].

Computes the Abel projection of P_{g,y}(k):

\[\Sigma_y(r_p) = \frac{1}{2\pi^2} \int_0^\infty k\,P_{g,y}(k)\,J_0(k r_p)\,dk\]
Parameters:

  • rp_arr ((NR,) projected separations [Mpc/h])

  • z, theta_cosmo, hod_params (as for _pk_tables_gy)

Returns:

sigma_y ((NR,) [dimensionless])