py-SP(k) - A hydrodynamical simulation-based model for the impact of baryon physics on the non-linear matter power spectrum

                      _____ ____   ____   _
    ____  __  __     / ___// __ \_/_/ /__| |
   / __ \/ / / /_____\__ \/ /_/ / // //_// /
  / /_/ / /_/ /_____/__/ / ____/ // ,<  / /
 / .___/\__, /     /____/_/   / //_/|_|/_/
/_/    /____/                 |_|    /_/

py-SP(k) (Salcido et al. 2023) is a python package aimed at predicting the suppression of the total matter power spectrum due to baryonic physics as a function of the baryon fraction of haloes and redshift.

Requirements

Runtime dependencies:

• numpy
• pydantic (required; used for input validation)
• scipy

Optional dependencies:

• astropy (required only for the cosmology-based relations; see Methods 2 and 3)

Example (notebook) dependencies:

• pandas, matplotlib, cycler
• emcee (MCMC sketch)
• corner (posterior corner plot)

Installation

Using pip:

pip install pyspk

If you need a user install:

pip install --user pyspk

If you want the cosmology-based relations (Methods 2 and 3):

pip install "pyspk[cosmology]"

Using uv (recommended for working from source / running examples):

uv sync

To run the example notebook dependencies:

uv sync --extra examples

Usage

The main entrypoints are:

• pyspk.sup_model: compute suppression $P_\mathrm{hydro}(k)/P_\mathrm{DM}(k)$
• pyspk.get_limits: fitting limits for $f_b$ as a function of mass/redshift
• pyspk.optimal_mass: optimal halo mass as a function of scale/redshift

All user-facing inputs are validated via Pydantic models (clear errors are raised if required parameters are missing or inconsistent).

py-SP(k) is not restrictive to a particular shape of the baryon fraction – halo mass relation. To provide flexibility, there are 4 supported ways to specify the required $f_b$ - $M_\mathrm{halo}$ relation. A Jupyter notebook with more detailed examples is available at examples/pySPk_Examples.ipynb.

Quickstart

import pyspk as spk

k, sup = spk.sup_model(SO=200, z=0.125, fb_a=0.4, fb_pow=0.3, fb_pivot=10**13.5)

Method 1: Using a power-law fit to the $f_b$ - $M_\mathrm{halo}$ relation

py-SP(k) can be provided with power-law fitted parameters to the $f_b$ - $M_\mathrm{halo}$ relation using the functional form:

$$f_b/(\Omega_b/\Omega_m)=a\left(\frac{M_{SO}}{M_{\mathrm{pivot}}}\right)^{b},$$

where $M_{SO}$ could be either $M_{200c}$ or $M_{500c}$ in $\mathrm{M}_ \odot$, $a$ is the normalisation of the $f_b$ - $M_\mathrm{halo}$ relation at $M_\mathrm{pivot}$, and $b$ is the power-law slope. The power-law can be normalised at any pivot point in units of $\mathrm{M}_ {\odot}$. If a pivot point is not given, spk.sup_model() uses a default pivot point of $M_{\mathrm{pivot}} = 1 \mathrm{M}_ \odot$. $a$, $b$ and $M_\mathrm{pivot}$ can be specified at each redshift independently.

Next, we show a simple example using power-law fit parameters:

import pyspk as spk

z = 0.125
fb_a = 0.4
fb_pow = 0.3
fb_pivot = 10**13.5

k, sup = spk.sup_model(SO=200, z=z, fb_a=fb_a, fb_pow=fb_pow, fb_pivot=fb_pivot)

Method 2: Redshift-dependent power-law fit to the $f_b$ - $M_\mathrm{halo}$ relation

For the mass range that can be relatively well probed in current X-ray and Sunyaev-Zel'dovich effect observations (roughly $10^{13} \leq M_{500c} [\mathrm{M}_ \odot] \leq 10^{15}$), the total baryon fraction of haloes can be roughly approximated by a power-law with constant slope (e.g. Mulroy et al. 2019; Akino et al. 2022). Akino et al. (2022) determined the of the baryon budget for X-ray-selected galaxy groups and clusters using weak-lensing mass measurements. They provide a parametric redshift-dependent power-law fit to the gas mass - halo mass and stellar mass - halo mass relations, finding very little redshift evolution.

We implemented a modified version of the functional form presented in Akino et al. (2022), to fit the total $f_b$ - $M_\mathrm{halo}$ relation as follows:

$$f_b/(\Omega_b/\Omega_m)= \left(\frac{e^\alpha}{100}\right) \left(\frac{M_{500c}}{10^{14} \mathrm{M}_ \odot}\right)^{\beta - 1} \left(\frac{E(z)}{E(0.3)}\right)^{\gamma},$$

where $\alpha$ sets the power-law normalisation, $\beta$ sets power-law slope, $\gamma$ provides the redshift dependence and $E(z)$ is the usual dimensionless Hubble parameter. For simplicity, we use the cosmology implementation of astropy to specify the cosmological parameters in py-SP(k).

Note that this power-law has a normalisation that is redshift dependent, while the the slope is constant in redshift. While this provides a less flexible approach compared with Methods 1 (simple power-law) and Method 4 (binned data), we find that this parametrisation provides a reasonable agreement with our simulations up to redshift $z=1$, which is the redshift range proved by Akino et al. (2022). For higher redshifts, we find that simulations require a mass-dependent slope, especially at the lower mass range.

In the following example we use the redshift-dependent power-law fit parameters with a flat LambdaCDM cosmology. Note that any astropy cosmology could be used instead.

import pyspk as spk
from astropy.cosmology import FlatLambdaCDM

H0 = 70
Omega_m = 0.2793

cosmo = FlatLambdaCDM(H0=H0, Om0=Omega_m)

alpha = 4.189
beta = 1.273
gamma = 0.298
z = 0.5

k, sup = spk.sup_model(SO=500, z=z, alpha=alpha, beta=beta, gamma=gamma, cosmo=cosmo)

Method 3: Redshift-dependent double power-law fit to the $f_b$ - $M_\mathrm{halo}$ relation

We implemented a redshift-dependent double power-law fit to the total $f_b$ - $M_\mathrm{halo}$ relation as follows:

$$f_b/(\Omega_b/\Omega_m)= \frac{1}{2} \epsilon \left[\left(\frac{M_{500c}}{M_{\mathrm{pivot}}}\right)^{\alpha} + \left(\frac{M_{500c}}{M_{\mathrm{pivot}}}\right)^{\beta}\right] \left(\frac{E(z)}{E(0.3)}\right)^{\gamma},$$

where $\epsilon$ sets the normalisation at the pivot point, $M_{\mathrm{pivot}}$, in units of $\mathrm{M}_ {\odot}$, $\alpha$ and $\beta$ are the power-law slopes at low and high mass respectively, $\gamma$ provides the redshift dependence and $E(z)$ is the usual dimensionless Hubble parameter. For simplicity, we use the cosmology implementation of astropy to specify the cosmological parameters in py-SP(k).

We find that this redshift-dependent double power-law form provides a good fit to the whole range of baryon fractions in the ANTILLES simulations.

In the following example we use the redshift-dependent double power-law fit parameters with a flat LambdaCDM cosmology. Note that any astropy cosmology could be used instead.

import pyspk as spk
from astropy.cosmology import FlatLambdaCDM

H0 = 68.0
Omega_m = 0.306

cosmo = FlatLambdaCDM(H0=H0, Om0=Omega_m)

epsilon = 0.410
alpha = -0.385
beta = 0.451
gamma = 0.125
m_pivot = 1e13
z = 0.2

k, sup = spk.sup_model(
    SO=500,
    z=z,
    epsilon=epsilon,
    alpha=alpha,
    beta=beta,
    gamma=gamma,
    m_pivot=m_pivot,
    cosmo=cosmo,
)

Method 4: Binned data for the $f_b$ - $M_\mathrm{halo}$ relation

The final, and most flexible method is to provide py-SP(k) with the baryon fraction binned in bins of halo mass. This could be, for example, obtained from observational constraints, measured directly form simulations, or sampled from a predefined distribution or functional form. For an example using data obtained from the BAHAMAS simulations (McCarthy et al. 2017), please refer to the examples provided.

Structured (Pydantic) API

If you prefer an explicit, validated input object (e.g., for configuration-driven workflows), the Pydantic models are available in pyspk.api:

relation.kind supports the same four relation kinds used throughout the docs: "power_law", "cosmo_power_law", "double_power_law", and "binned".

from pyspk.api import SupModelRequest

req = SupModelRequest(
    SO=200,
    z=0.125,
    relation={"kind": "power_law", "fb_a": 0.4, "fb_pow": 0.3, "fb_pivot": 10**13.5},
    k_min=0.1,
    k_max=8,
    n=100,
)

You can also inspect relation-specific requirements via help(pyspk.sup_model).

MCMC / fast inner loops (errors=False)

Relation kinds (relation_kind)

The fast evaluator API requires an explicit relation_kind to avoid per-call validation. Supported values are:

• "power_law" (Method 1): requires fb_a, fb_pow (optional fb_pivot).
• "cosmo_power_law" (Method 2; Akino et al. 2022): requires alpha, beta, gamma and either cosmo or efunc.
• "double_power_law" (Method 3): requires epsilon, alpha, beta, gamma, m_pivot and either cosmo or efunc.
• "binned" (Method 4): requires M_halo, fb (optional extrapolate).

If you are calling py-SP(k) in a tight loop (e.g., an MCMC likelihood) and you typically use errors=False, you can build a fast evaluator. This avoids per-call Pydantic validation and caches the k-grid and fitting-limit interpolators.

import pyspk

evaluator = pyspk.build_sup_model_evaluator(
    SO=500,
    relation_kind="cosmo_power_law",
    k_max=8,
    n=100,
)

# In your MCMC loop:
k, sup = evaluator(z=z, alpha=alpha, beta=beta, gamma=gamma, cosmo=cosmo)

For cosmology-based relations, you may also pass efunc directly (a callable returning $E(z)$) instead of an astropy cosmology object.

Optional: emcee sketch

If you use emcee (or a similar sampler), the key idea is to build the evaluator once and call it inside log_prob. Install with pip install emcee (or uv sync --extra examples when working from source).

import numpy as np
import pyspk as spk

# Optional (for Method 2/3):
from astropy.cosmology import FlatLambdaCDM

cosmo = FlatLambdaCDM(H0=70, Om0=0.2793)
evaluator = spk.build_sup_model_evaluator(
    SO=500,
    relation_kind="cosmo_power_law",
    k_max=8,
    n=100,
)

k_data = np.logspace(-1, np.log10(8.0), 60)
sup_data = np.ones_like(k_data)  # replace with your measured/target suppression
sigma = 0.05

def log_prob(theta: np.ndarray) -> float:
    alpha, beta, gamma = theta
    z = 0.5
    k, sup = evaluator(z=z, alpha=alpha, beta=beta, gamma=gamma, cosmo=cosmo)
    sup = np.interp(k_data, k, sup)
    if not np.all(np.isfinite(sup)):
        return -np.inf
    return -0.5 * np.sum(((sup - sup_data) / sigma) ** 2)

Priors (Method 2: Akino-like Single Power-Law)

These priors apply to the Method 2 redshift-dependent single power-law (Akino-like) functional form only (relation_kind="cosmo_power_law"; parameters alpha, beta, gamma).

While py-SP(k) was calibrated using a wide range of sub-grid feedback parameters, some applications may require a more limited range of baryon fractions that encompass current observational constraints. For such applications, we used the gas mass - halo mass and stellar mass - halo mass constraints from the fits in Table 5 in Akino et al. (2022), and find the subset of simulations from our 400 models that agree to within $\pm 2$ or $3 \times \sigma$ of the inferred baryon budget at redshift $z=0.1$. We note that for our simulations, we include all stellar and gas particles within a spherical overdensity radius. Hence, in order to make reasonable comparisons with the fits in Akino et al. (2022), we included an additional 15% contribution to the total stellar masses from the contribution of blue galaxies, and 30% additional stellar mass to the brightest cluster galaxies (BCGs) to account for the diffuse intracluster light (ICL, see Akino et al. 2022).

We utilised the simulations satisfying these restrictions to determine the redshift-dependent power-law parameters for the $f_b$ - $M_\mathrm{halo}$ relation up to redshift $z=1$ (Method 2), and then utilised these parameters to infer suitable priors. We limited the fitting range to $6 \times 10^{12} \leq M_{500c} [\mathrm{M}_ \odot] \leq 10^{14}$.

Priors inferred from simulations that fall within $\pm 2 \times \sigma$ of the inferred baryon budget:

Parameter	Description	Prior
$\alpha$	Normaliasation	$\mathcal{N}$(4.16, 0.07)
$\beta$	Slope	$\mathcal{N}$(1.20, 0.05)
$\gamma$	Redshift evolution	$\mathcal{N}$(0.39, 0.09)

where $\mathcal{N}(\mu,\sigma)$ is a Gaussian distribution with mean $\mu$ and standard deviation $\sigma$.

Priors inferred from simulations that fall within $\pm 3 \times \sigma$ of the inferred baryon budget:

Parameter	Description	Prior
$\alpha$	Normaliasation	$\mathcal{N}$(4.18, 0.12)
$\beta$	Slope	$\mathcal{N}$(1.26, 0.08)
$\gamma$	Redshift evolution	$\mathcal{N}$(0.42, 0.10)

where $\mathcal{N}(\mu,\sigma)$ is a Gaussian distribution with mean $\mu$ and standard deviation $\sigma$.

Acknowledging the code

Please cite py-SP(k) using:

```bibtex
@ARTICLE{SPK_Salcido_2023,
    author = {Salcido, Jaime and McCarthy, Ian G and Kwan, Juliana and Upadhye, Amol and Font, Andreea S},
    title = "{SP(k) – a hydrodynamical simulation-based model for the impact of baryon physics on the non-linear matter power spectrum}",
    journal = {Monthly Notices of the Royal Astronomical Society},
    volume = {523},
    number = {2},
    pages = {2247-2262},
    year = {2023},
    month = {05},
    issn = {0035-8711},
    doi = {10.1093/mnras/stad1474},
    url = {https://doi.org/10.1093/mnras/stad1474},
    eprint = {https://academic.oup.com/mnras/article-pdf/523/2/2247/50512773/stad1474.pdf},
}
```

For any questions and enquires please contact me via email at jaime.salcido@gmail.com