Perturbations: Structure Formation¶

The Perturbations module computes structure formation in the universe.

Building on the Background module, Perturbations calculates how galaxies, clusters, and cosmic web structures form and evolve.

Overview¶

This module computes the growth and evolution of density fluctuations over cosmic time, including:

Matter power spectra at various redshifts and scales
Growth rates of structures at different epochs
Key cosmological parameters like σ₈

This module bridges smooth background cosmology and the observed clustered universe.

The Perturbations Protocol¶

Protocol Definition: cloelib.cosmology.cosmology.Perturbations

The Perturbations protocol defines what every perturbation calculator must provide. This module requires a Background instance as it depends on cosmological distances and densities.

Units convention

All wavenumbers within cloelib are always in 1/Mpc (not h/Mpc). Power spectra are therefore in Mpc³. Any class implementing the Perturbations protocol must adopt this convention. If the underlying code uses h/Mpc internally, the implementation wrapper is responsible for converting inputs and outputs.

Required Property¶

background: Reference to the associated Background object

This is key. Perturbations always needs a Background to compute distances, densities, etc.

Required Methods¶

`matter_power_spectrum(zs, ks)`¶

Compute the matter power spectrum for total matter Pmm(k, z).

Inputs:

zs: Redshifts (can be array)
ks: Wavenumbers in 1/Mpc (can be array)

Returns: Power spectrum in Mpc³

Note: Can be linear or non-linear depending on implementation.

`matter_power_spectrum_cb(zs, ks)`¶

Compute the CDM+baryons matter power spectrum Pcb(k, z).

Inputs:

zs: Redshifts (can be array)
ks: Wavenumbers in 1/Mpc (can be array)

Returns: Power spectrum in Mpc³

Note: Can be linear or non-linear depending on implementation.

cb x matter cross spectrum: `cb_cross_prescription`¶

Pairing a PositionsTracer(use_Pcb=True) with a tracer of the total matter (lensing, magnification, CMB lensing) needs the cb x matter cross spectrum, which a backend does not provide directly. Every backend implementing matter_power_spectrum_cb declares how it is built through the class attribute cb_cross_prescription, one of the constants in cloelib.cosmology.cosmology:

Value	Cross spectrum	Exact when
`CB_CROSS_CORRELATED` (`"correlated"`)	\(\sqrt{P_{cb} P_{mm}}\)	the backend applies the same non-linear boost to cb and total matter, \(P_{cb}^{NL}/P_{cb}^{L} = P_{mm}^{NL}/P_{mm}^{L}\) (always true if linear)
`CB_CROSS_LINEAR_NU` (`"linear_neutrinos"`)	\(f_{cb}\,(P_{cb}^{NL} - P_{cb}^{L}) + \sqrt{P_{cb}^{L} P_{mm}^{L}}\)	massive neutrinos, and the cb x neutrino spectrum, stay linear, \(P_{mm}^{NL} - P_{mm}^{L} = f_{cb}^2\,(P_{cb}^{NL} - P_{cb}^{L})\); needs `linearperturbations`

with \(f_{cb} = \Omega_{cb,0} / \Omega_{m,0}\). A backend declares the prescription whose condition its own non-linear \(P_{mm}\) and \(P_{cb}\) satisfy. How much the two prescriptions differ depends on the backend: at \(m_\nu = 0.15\) eV, about \(3 \times 10^{-4}\) in the cross spectrum for CAMB, at most \(2 \times 10^{-3}\) for the emulators and about \(1\%\) for CLASS with HMcode2020.

Prescription	Backends
`CB_CROSS_CORRELATED`	all linear backends; CAMB, HMcode2020Emu and CosmoPowerJAX (`NonLinear`, `NonLinearHalofit`) non-linear
`CB_CROSS_LINEAR_NU`	EE2, BACCOemu, JAX; CLASS, hi_class, MGCLASS and mochi_class non-linear
that of the baseline	MG-CosmoPowerJAX non-linear and emantis

Emulators (HMcode2020Emu, CosmoPowerJAX): \(P_{mm}^{NL}\) and \(P_{cb}^{NL}\) are emulated separately, and at \(m_\nu \lesssim 0.15\) eV the two conditions cannot be told apart within the emulation accuracy (\(\sim 10^{-3}\)). At higher masses the HMcode-based emulators follow the uniform boost of CAMB, so they are correlated. For the CosmoPowerJAX halofit sector the check is inconclusive; it shares the label of the HMcode one.
CLASS family (CLASS, hi_class, MGCLASS, mochi_class): the non-linear model (halofit/HMcode) is run on \(P_{cb}\) separately from \(P_{mm}\), so neither condition holds exactly. CB_CROSS_LINEAR_NU is used since it only relies on the backend's own non-linear \(P_{cb}\) and on linear spectra. Pass the matching *LinearPerturbations as linearperturbations to get the cross spectrum.
Boosts of a baseline (MG-CosmoPowerJAX, emantis): the same boost rescales the ΛCDM baseline \(P_{mm}\) and \(P_{cb}\), so cb_cross_prescription is a property returning the baseline's. emantis is calibrated without massive neutrinos, and its matter_power_spectrum_cb raises a ValueError when \(m_\nu > 0\).

A new backend providing matter_power_spectrum_cb must declare it, as a class attribute or a property (and be added to tests/test_cb_cross_prescription.py): computing the cross spectrum for a backend without a valid cb_cross_prescription raises a ValueError. The attribute is not part of the Perturbations protocol, since only the cross spectrum needs it.

`growth_factor(zs, ks)`¶

Calculate the linear growth factor D(z, k).

Normalized such that D(z=0) ≈ 1 in matter-dominated era.

`growth_rate(zs=None, ks=None)`¶

Calculate the linear growth rate f(z) = d ln D / d ln a. Every implementation follows the same convention:

zs=None returns the growth rate on the redshift grid the instance was built on; otherwise it is evaluated at zs.
ks=None returns an array of shape (nz,); otherwise an array of shape (nz, nk). A scale-independent growth rate is broadcast along k.

Scale-dependent implementations evaluate f(z, k) at the given ks. Without ks, hi_class and mochi_class use k = 1 Mpc⁻¹, while BACCOemu, which has no such default, raises a ValueError.

MGCLASS computes the growth rate by finite differences on its redshift grid, so it needs to be built on at least 3 redshifts to provide it.

`sigma8_0()`¶

Compute σ₈ at redshift z=0.

The RMS matter fluctuation in 8 Mpc/h spheres—a key cosmological parameter.

Constructors¶

All implementations are built with the same keyword arguments, so they can be swapped without changing the calling code:

lin = LinearPerturbations(background=bg, redshifts=zs)
nonlin = NonLinearPerturbations(background=bg, linearperturbations=lin, redshifts=zs)

This is described by the LinearPerturbationsFactory and NonLinearPerturbationsFactory protocols in cloelib.cosmology.cosmology, which can be used to type a class passed around as a parameter. Their type parameters are the background and, for the nonlinear one, the linear perturbations the class is built from. For example, NonLinearPerturbationsFactory[CAMBBackground, HMemuLinearPerturbations] is satisfied by HMemuNonLinearPerturbations, while NonLinearPerturbationsFactory[CAMBBackground, CAMBLinearPerturbations] is not, since HMcode2020Emu needs the cached Pk/Pk_cb grid of its linear perturbations (WithLinearSpectrumGrid).

Options specific to a backend (e.g. nonlinear_model) are passed as additional keyword arguments and are not part of these protocols. In particular, the HMcode2020 baryonic feedback parameter log10TAGN is accepted by the HMcode-based backends (CAMB, CLASS, mochi_class, HMcode2020Emu and CosmoPower-JAX). CAMB, CLASS and mochi_class select their HMcode2020 feedback model when log10TAGN is given without a nonlinear model, and raise a ValueError if the chosen model ignores it.

Backends that compute the nonlinear spectrum themselves (CAMB, CosmoPower-JAX) accept but do not use linearperturbations; CLASS, hi_class, MGCLASS and mochi_class only use it for the cb x matter cross spectrum (see cb_cross_prescription). EmantisFofrNonLinearPerturbations additionally requires the ΛCDM nonlinear perturbations and fR0.

Existing Implementations¶

CAMB Perturbations¶

CAMB provides separate linear and non-linear perturbation implementations: CAMBLinearPerturbations and CAMBNonLinearPerturbations.

Interfaces with CAMB for perturbation calculations.

Location: cloelib/cosmology/camb_cosmology.py

When to use: accurate, production-ready, but sometimes slow

Features:

Linear and non-linear power spectra
Multiple non-linear models (Halofit, HMCode, ...)
Neutrino effects
Modified gravity

Example:

from cloelib.cosmology.camb_cosmology import (
    CAMBBackground,
    CAMBLinearPerturbations,
    CAMBNonLinearPerturbations,
)
import numpy as np

# First, create background
bg = CAMBBackground(
    H0=67.5,
    Omega_b0=0.0492,
    Omega_cdm0=0.2650,
    As=2.1e-9,
    ns=0.965,
    alpha_s=0.0,   # running of the spectral index
    # ... other parameters
)

# Then create linear and non-linear perturbations
redshifts = np.array([0.0, 0.5, 1.0])
linear_pert = CAMBLinearPerturbations(background=bg, redshifts=redshifts)
pert = CAMBNonLinearPerturbations(
    background=bg,
    linearperturbations=linear_pert,
    redshifts=redshifts,
)

# Compute power spectrum
z = np.array([0.0, 0.5, 1.0])
k = np.logspace(-3, 1, 100)  # k in 1/Mpc
P_k_z = pert.matter_power_spectrum(z, k)

# Get σ₈
sigma8 = pert.sigma8_0()
print(f"σ₈ = {sigma8:.4f}")

CLASSPerturbations¶

Interfaces with CLASS.

Location: cloelib/cosmology/class_cosmology.py

When to use: CLASS-specific features, comparison studies

Example:

from cloelib.cosmology.class_cosmology import (
    CLASSBackground,
    CLASSLinearPerturbations,
    CLASSNonLinearPerturbations,
)

bg = CLASSBackground(...)
lin = CLASSLinearPerturbations(background=bg, redshifts=zs)

# HMcode2020 with baryonic feedback, as with CAMB and HMcode2020Emu
nonlin = CLASSNonLinearPerturbations(
    background=bg, linearperturbations=lin, redshifts=zs, log10TAGN=7.8
)

HMCode2020EmuPerturbations¶

Fast emulator for non-linear power spectra using HMCode2020Emu.

Location: cloelib/cosmology/HMcode2020Emu_cosmology.py

When to use: Fast non-linear predictions, MCMC sampling

Features:

Lightning-fast emulator
Accurate non-linear P(k)
Limited parameter range

EE2Perturbations¶

Simulation-based emulator for non-linear power spectra using euclidemu2.

Location: cloelib/cosmology/EE2_cosmology.py

When to use: Fast non-linear predictions based on simulations, MCMC sampling

Features:

Lightning-fast (emulator.)
Accurate DM-only non-linear boost for P_mm
More limited parameter range

BACCOemuPerturbations¶

Accurate and fast emulators of the linear, non-linear, and baryonic power spectra using BACCOemu.

Location: cloelib/cosmology/baccoemu_cosmology.py

When to use: Predictions for the matter clustering, in the linear and non linear regime. Predictions of baryonic effects.

Features:

Fast predictions of linear power spectra, growth factors, growth rates, and amplitude of fluctuations;
Accurate total and cold non-linear matter power spectrum emulated from high-resolution simulations;
Inclusion of baryonic effects through baryonification;
Large cosmological parameter range;
Neural network evaluation with JAX;

FlamingoBaryonResponseEmulatorPerturbations¶

CAMB nonlinear power spectrum with a baryonic suppression correction from the FLAMINGO hydrodynamical simulation suite.

Location: cloelib/cosmology/FlamingoBaryonResponseEmulator_cosmology.py

When to use: Accurate baryonic feedback modelled on FLAMINGO hydrodynamical simulations; can vary AGN jet fraction and ICM gas/stellar parameters.

Features:

CAMB-based nonlinear baseline
FLAMINGO emulator response \(B(k,z) = P_\text{hydro}/P_\text{DMO}\) applied at evaluation time
Three baryonic parameters: fgas_sigma, Mstar_sigma, jet_fraction
Valid for \(z \leq 3\) and \(k \leq 10^{1.5}\,h/\text{Mpc}\); outside these ranges the response is automatically clamped to 1
Standalone FlamingoBaryonBoostMixin composable with any nonlinear perturbations class

Example:

from cloelib.cosmology.FlamingoBaryonResponseEmulator_cosmology import (
    CAMBNonLinearFLAMINGOBaryonicPerturbations,
)

pert = CAMBNonLinearFLAMINGOBaryonicPerturbations(
    background=bg,
    linearperturbations=lin_pert,
    redshifts=zs,
    baryon_kwargs=dict(
        fgas_sigma=0.0,    # 0σ from calibrated gas fraction
        Mstar_sigma=0.0,   # 0σ from stellar mass function
        jet_fraction=0.0,  # pure thermal AGN feedback
    ),
)

k = np.logspace(-2, 1, 100)  # k in 1/Mpc
Pk = pert.matter_power_spectrum(np.array([0.5]), k)  # shape (n_k,)

ReACTEmu / BoostedPerturbations¶

Modified-gravity nonlinear boost module based on ReACT and the MGEmu emulator package.

Location: cloelib/cosmology/ReACTEmu_cosmology.py

When to use: Fast nonlinear modified-gravity corrections on top of an existing LCDM perturbation pipeline.

Features:

Supports fr, dgp, gamma, mu, wCDM, w0waCDM, and ide model tags
Builds a boost interpolator B(z, k) and applies it to a baseline nonlinear matter spectrum
Includes configurable high-redshift behavior via high_z_policy and z_decay
Exposes boosted matter_power_spectrum, growth_factor, and sigma8_0

Installation:

pip install ".[react]"

This extra installs MGEmu together with the TensorFlow support it expects. In environments where you want to drive the baseline spectra with CAMB as well, pip install ".[react,camb]" is a sensible default.

Example:

import numpy as np

from cloelib.cosmology.camb_cosmology import (
    CAMBBackground,
    CAMBLinearPerturbations,
    CAMBNonLinearPerturbations,
)
from cloelib.cosmology.ReACTEmu_cosmology import (
    MGemuNonlinearBoost,
    BoostedPerturbations,
)

zs = np.array([0.0, 0.5, 1.0])
k = np.logspace(-3, 1, 100)

background = CAMBBackground(...)
linear_pert = CAMBLinearPerturbations(background, zs)
nonlinear_pert = CAMBNonLinearPerturbations(
    background, linearperturbations=linear_pert, redshifts=zs
)

boost = MGemuNonlinearBoost(
    background,
    linear_pert,
    zs,
    gravity_model="fr",
    mgpars={"fr0": 1e-5},
)
mg_pert = BoostedPerturbations(linear_pert, nonlinear_pert, boost.MGboost_interp)
P_mg = mg_pert.matter_power_spectrum(zs, k)

Usage notes:

MGemuNonlinearBoost expects an existing background plus a linear LCDM perturbation object. BoostedPerturbations then combines that boost with a baseline nonlinear spectrum.
Internal wavenumbers follow the standard cloelib convention in 1/Mpc, even when external emulator data are defined in h/Mpc.
For models other than wCDM, w0waCDM, ide, and mu, the current implementation assumes a LCDM background with w0 = -1 and wa = 0.

Emulator ranges:

fr: Omega_m 0.24-0.35, Omega_b 0.04-0.06, H0 63-75, ns 0.9-1.01, As 1.7e-9-2.5e-9, Omega_nu 1e-11-0.00317, fR0 1e-10-1e-4, z 0-2.0
dgp: Omega_m 0.2899-0.3392, Omega_b 0.04044-0.05686, H0 62.9-73.1, ns 0.9432-0.9862, As 1.5e-9-2.7e-9, Omega_nu 1e-11-0.00317, omegarc 1e-3-100, z 0-2.4
gamma: Omega_m 0.2899-0.3392, Omega_b 0.04044-0.05686, H0 63.8-73.1, ns 0.9432-0.9862, As 1.9511e-9-2.2669e-9, Omega_nu 1e-11-0.00317, gamma 0-1, q1 -5 to 5, z 0-2.4
mu: Omega_m 0.24-0.4, H0 60-84, As 1.7e-9-2.5e-9, w0 -1.5 to -0.5, wa -0.5 to 0.5, mu0 -0.999 to 3, c1 -0.3333 to 1, lam 0-2, q1/q2/q3 -2 to 2, z 0-2.5
wCDM, w0waCDM, ide (shared ds backend): Omega_m 0.22-0.37, Omega_b 0.03-0.08, H0 63-84, ns 0.8-1.1, As 1.7e-9-2.5e-9, Omega_nu 1e-11-0.00317, w0 -1.3 to -0.5, wa -2 to 0.5, xi 0-150, z 0-2.5

Stay within these training ranges when sampling. Outside them, the implementation clips emulator inputs and applies its configured high-redshift policy, which is convenient for robustness but should not be treated as a new calibration region.

TabulatedBoost / TabulatedBoostedPerturbations¶

Lightweight wrapper for applying a tabulated beyond-LCDM nonlinear boost to an existing nonlinear matter power spectrum.

Location: cloelib/cosmology/TabulatedBoost_cosmology.py

When to use: You already have boost data from simulations or another external pipeline and want to interpolate it onto the cloelib perturbation grid.

Features:

Reads tabulated boost files with columns [k, B(z_1), B(z_2), ...]
Converts input k values from h/Mpc to the internal 1/Mpc convention
Reuses cloelib's wavenumber extrapolation machinery
Supports power_law, freeze, one, and taper high-redshift policies
Exposes a TabulatedBoostedPerturbations wrapper with boosted matter_power_spectrum, growth_factor, and sigma8_0

Input format:

Column 1: k in h/Mpc
Remaining columns: boost values B(k, z_i) at each tabulated redshift
z_cols must list the redshifts corresponding to those boost columns, in the same order as the file

Example:

import numpy as np

from cloelib.cosmology.camb_cosmology import (
    CAMBBackground,
    CAMBLinearPerturbations,
)
from cloelib.cosmology.HMcode2020Emu_cosmology import HMemuNonLinearPerturbations
from cloelib.cosmology.TabulatedBoost_cosmology import (
    TabulatedNonlinearBoost,
    TabulatedBoostedPerturbations,
)

zs = np.linspace(0.0, 4.0, 100)
z_cols = [
    3.017980,
    2.479559,
    2.161320,
    2.013288,
    1.609499,
    1.259818,
    1.000000,
    0.823800,
    0.677100,
    0.552300,
    0.444200,
    0.349200,
    0.264800,
    0.188900,
    0.120200,
    0.057540,
    0.000031,
    0.0,
]

background = CAMBBackground(...)
linear_perturbations = CAMBLinearPerturbations(background, zs)
nonlinear_perturbations = HMemuNonLinearPerturbations(
    background,
    linear_perturbations,
    zs,
)

tabulated_boost = TabulatedNonlinearBoost(
    background,
    linear_perturbations,
    zs,
    "m11_p01_boosts.txt",
    z_cols,
    high_z_policy="freeze",
)
boosted_perturbations = TabulatedBoostedPerturbations(
    linear_perturbations,
    nonlinear_perturbations,
    tabulated_boost.MGboost_interp,
)

CosmoPowerJAXPerturbations¶

Fast JAX-based emulator for linear and nonlinear power spectra using cosmopower-jax.

Location: cloelib/cosmology/cosmopower_jax_cosmology.py

When to use: Fast predictions with full JAX compatibility, gradient-based inference, MCMC sampling

Features:

Full JAX compatibility — automatic differentiation and JIT compilation
Linear and nonlinear P(k) and P_cb(k)
Two nonlinear prescriptions: HMcode2020 (with baryonic feedback) and halofit (dark-matter-only)
σ₈(z), fσ₈(z), growth factors D(z,k) and D_cb(z,k), growth rate f(z)
Baseline (ΛCDM/wCDM/w0waCDM) and extended (curvature, running spectral index) cosmologies
Emulator files downloaded automatically on first use from the Euclid-DR1-matter-emulators GitHub repository; the k-mode grid is read directly from each emulator

Available classes¶

Every Perturbations class exposes the same inner classes, each providing both the total-matter spectrum (matter_power_spectrum) and the CDM+baryon spectrum (matter_power_spectrum_cb, built at construction from the paired cb emulator), with the cached Pk/Pk_cb grids and growth_factor/growth_factor_cb:

Linear — linear P(k) and P_cb(k)
NonLinear — nonlinear P(k) and P_cb(k) using HMcode2020, with baryonic feedback controlled by log10TAGN
NonLinearHalofit — nonlinear P(k) and P_cb(k) using halofit (Takahashi 2012); dark-matter-only, so log10TAGN is accepted for interface compatibility but ignored.

Class	Cosmology	N_mnu	Nonlinear recipes
`CosmoPowerJAXLCDMPerturbations`	ΛCDM	0,1,2,3	HMcode2020, halofit
`CosmoPowerJAXwCDMPerturbations`	wCDM	0,1,2,3	HMcode2020, halofit
`CosmoPowerJAXw0waCDMPerturbations`	w0waCDM	0,1,2,3	HMcode2020, halofit
`CosmoPowerJAXLCDMCurvaturePerturbations`	ΛCDM + curvature (Ω_k)	0,1,3	HMcode2020
`CosmoPowerJAXw0waCurvaturePerturbations`	w0waCDM + curvature	0,1,3	HMcode2020
`CosmoPowerJAXLCDMRunningIndexPerturbations`	ΛCDM + running (α_s)	0,1,3	HMcode2020
`CosmoPowerJAXw0waRunningIndexPerturbations`	w0waCDM + running	0,1,3	HMcode2020

The neutrino configuration is selected automatically from background.N_mnu: 0 (massless), 1 (one massive), 2 (two degenerate), 3 (three degenerate). Baseline cosmologies cover 0–3; the extended (curvature / running) cosmologies cover 0, 1 and 3 (an unsupported value raises ValueError).

Selecting the nonlinear prescription¶

The recipe is chosen by which nonlinear class you use — the constructor signature is identical, so it is a drop-in swap:

# HMcode2020 (default) — baryonic feedback via log10TAGN
nl = CosmoPowerJAXLCDMPerturbations.NonLinear(
    background=bg, linearperturbations=lin, redshifts=zs, log10TAGN=7.8
)

# halofit — dark-matter-only (log10TAGN ignored)
nl = CosmoPowerJAXLCDMPerturbations.NonLinearHalofit(
    background=bg, linearperturbations=lin, redshifts=zs
)

Parameter ranges¶

Inputs are validated against a single global guardrail (CP_EMULATOR_BOUNDS); a parameter outside its range raises ValueError. The redshift z is not bounds-checked (it is clipped to the training range and extrapolated; z=0 is always added to the emulator grid, so σ₈(z=0) and the D(z=0)=1 normalisation are exact even when the requested redshifts start above zero). These global bounds are deliberately loose — each emulator's actual training box is tighter and varies by family.

Parameter	Enforced range	Applies to
ombh2	[0.001, 0.1]	all
omch2	[0.05, 0.9]	all
H0	[20, 100]	all
ns	[0.6, 1.3]	all
lnAs	[1.61, 5.0]	all
mnu	[0, 1] eV	massive-neutrino runs
w0	[-3, 1]	w0waCDM
wa	[-3, 3]	w0waCDM
w	[-3, 0]	wCDM
Omega_k0	[-0.3, 0.3]	curvature
alpha_s	[-0.3, 0.3]	running
log10TAGN	[7.3, 8.5]	HMcode nonlinear

Example¶

import numpy as np
from cloelib.cosmology.camb_cosmology import CAMBBackground
from cloelib.cosmology.cosmopower_jax_cosmology import CosmoPowerJAXLCDMPerturbations

zs = np.array([0.0, 0.5, 1.0, 2.0])
ks = np.logspace(-4, 1, 200)

bg = CAMBBackground(
    H0=67.0, Omega_b0=0.049, Omega_cdm0=0.270,
    As=2.1e-9, ns=0.96, mnu=0.06, N_mnu=1,
    w0=-1.0, wa=0.0, Omega_k0=0.0, gamma_MG=0.0,
)

# Linear P(k)
lin = CosmoPowerJAXLCDMPerturbations.Linear(background=bg, redshifts=zs)
Pk = lin.matter_power_spectrum(0.0, ks)   # shape (1, len(ks))

# Nonlinear P(k) with baryonic feedback (HMcode2020)
nl = CosmoPowerJAXLCDMPerturbations.NonLinear(
    background=bg, linearperturbations=lin, redshifts=zs, log10TAGN=7.6
)
Pk_nl = nl.matter_power_spectrum(0.0, ks)

# sigma8 and fsigma8 as a function of redshift
print(f"sigma8(z=0) = {lin.sigma8[0]:.4f}")
print(f"fsigma8(z=0) = {lin.fsigma8[0]:.4f}")

JAXPerturbations¶

Pure JAX implementation for automatic differentiation.

Location: cloelib/cosmology/jax_cosmology.py

When to use: Computing gradients, Fisher forecasts, HMC sampling

MG CosmoPower-JAX boost perturbations¶

Applies a CosmoPower-JAX modified-gravity boost \(B(k,z) = P*\mathrm{MG}/P*{\Lambda\mathrm{CDM}}\) as a multiplicative operator on top of an external LCDM baseline perturbation object, so that the GR limit (\(\mu = \eta = 1\)) reproduces the baseline spectrum up to the boost emulators' own residual at GR (\(|B-1| \lesssim 3\times10^{-4}\) for \(k \le 1\,h/\mathrm{Mpc}\), up to \(2.5\times10^{-3}\) at \(k \sim 10\,h/\mathrm{Mpc}\) for the nonlinear multi-bin emulator); no regridding error is introduced.

Location: cloelib/cosmology/mg_cosmopower_jax_cosmology.py

When to use: Phenomenological \(\mu(z), \eta(z)\) modified gravity where the nonlinear MG power spectrum is obtained by boosting a fast LCDM emulator (e.g. CosmoPowerJAXLCDMPerturbations) rather than solving MG directly.

Two modes, selected by MGParams.bin_index:

Single-bin (bin_index is an int): \(\mu, \eta\) vary in one redshift bin (per-bin emulators mg-boost-{linear,nonlinear}-bin{0..4}.npz);
Multi-bin (bin_index=None): \(\mu\) (and \(\eta\), linear only) vary in all five bins simultaneously (joint emulators mg-boost-{linear,nonlinear}-multibin.npz).

Features:

Emulator-based MG boost applied to any external LCDM Perturbations solver;
Provides the modified lensing parameter \(\Sigma(z) = \mu(1+\eta)/2\) (a step function over the redshift bins), applied to every lensing kernel: cosmic shear, magnification bias and CMB lensing (photo.shear, photo.positions, cmb);
Scale-independent growth rate \(f(z)\) and \(\sigma_8\) from an internally-built linear MG \(P(k)\); the boost is applied on the baseline k-grid, so the GR limit is recovered up to the emulators' own residual (no regridding error);
Mutable MGParams holder for injecting sampled \(\mu, \eta\) each likelihood call;
Implemented in module-level MGLinearPerturbations / MGNonLinearPerturbations (statically checked against the Perturbations protocol); mg_perturbations returns subclasses with the MG parameters and LCDM baseline bound;
Emulators are downloaded on first use from the Euclid-DR1-matter-emulators GitHub repository and cached locally, mirroring CosmoPowerJAXPerturbations (no local model directory needed);

Example:

import numpy as np
from cloelib.cosmology.cosmopower_jax_cosmology import CosmoPowerJAXLCDMPerturbations as LCDM
from cloelib.cosmology.mg_cosmopower_jax_cosmology import MGParams, mg_perturbations

# single-bin (bin 4); use MGParams(mu=..., eta=...) with bin_index=None for multi-bin
mg = MGParams(mu=1.0, eta=1.0, bin_index=4)
Lin, NonLin = mg_perturbations(
    mg, baseline_linear=LCDM.Linear, baseline_nonlinear=LCDM.NonLinear
)
# in the sampling wrapper, before each loglike:  mg.mu, mg.eta = ...

Emulator ranges (enforced; predictions outside raise ValueError):

Parameter	Single-bin	Multi-bin
\(\Omega_m\)	0.25–0.35	0.25–0.40
\(\Omega_b\)	0.040–0.055	0.040–0.055
\(h\)	0.65–0.73	0.65–0.75
\(n_s\)	0.95–1.00	0.80–1.20
\(\ln(10^{10}A_s)\)	2.996–3.091	2.944–3.219
\(\mu\) (per bin)	0.9–1.1	0.9–1.1
\(\eta\) (per bin)	0.9–1.1	0.9–1.1
\(z\)	0.01–3	0–3

The \(z\) upper edge is enforced (and handled by the redshift clamp); the lower edge is not, so \(z = 0\) remains available for the growth/\(\sigma_8\) normalisation.

MGCLASSPerturbations¶

Interfaces with MGCLASS.

Location: cloelib/cosmology/mgclass_cosmology.py

When to use: MGCLASS-specific features, comparison studies

Example:

from cloelib.cosmology.mgclass_cosmology import MGCLASSBackground, MGCLASSPerturbations

bg = MGCLASSBackground(...)
pert = MGCLASSPerturbations(
    background=bg,
    # other parameters
)

hi_classPerturbations¶

Interfaces with hi_class.

Location: cloelib/cosmology/hi_class_cosmology.py

When to use: hi_class-specific features, comparison studies

Example:

from cloelib.cosmology.hi_class_cosmology import hi_classBackground, hi_classPerturbations

bg = hi_classBackground(...)
pert = hi_classPerturbations(
    background=bg,
    # other parameters
)

EmantisPerturbations¶

Accurate and fast emulator of the nonlinear matter power spectrum in modified gravity using e-MANTIS.

Location: cloelib/cosmology/emantis_cosmology.py

When to use: Predictions for the nonlinear matter clustering in f(R) gravity.

Features:

Fast predictions of the nonlinear matter power spectrum in f(R) gravity;
Accurate emulation of the nonlinear modified gravity boost based on N-body simulations;
Limited to the Hu & Sawicki model (n=1) with fR0 as free parameter;

Baryonic Boost Mixin Framework¶

cloelib provides a composable BaryonBoostMixin Protocol that decouples the baryonic suppression from the underlying nonlinear solver. This makes it easy to combine any nonlinear backend with any baryonic model.

Protocol location: cloelib.cosmology.cosmology.BaryonBoostMixin

`BaryonBoostMixin` Protocol¶

Any class implementing this protocol must expose:

def baryonic_suppression(self, zs, ks, k_hunit=False) -> np.ndarray:
    """Return B(k,z) = P_hydro / P_DMO, shape (n_z, n_k)."""

Three ready-made mixins are provided:

Mixin	Backend	Baryonic parameters
`FlamingoBaryonBoostMixin`	FLAMINGO emulator	`fgas_sigma`, `Mstar_sigma`, `jet_fraction`
`BACCOemuBaryonBoostMixin`	BACCOemu	`M_c`, `eta`, `beta`, `M1_z0_cen`, `theta_out`, `theta_inn`, `M_inn`
`HMcode2020BaryonBoostMixin`	HMcode2020emu	`log10TAGN`

Each mixin __init__ must be called after the base nonlinear __init__ (it reads self.background and emulator internals set by the base).

`with_baryon_boost` Factory¶

cloelib.cosmology.cosmology.with_baryon_boost(NonLinearClass, BaryonMixinClass) creates a combined class with zero boilerplate:

from cloelib.cosmology.cosmology import with_baryon_boost
from cloelib.cosmology.baccoemu_cosmology import BACCOemuNonLinearPerturbations
from cloelib.cosmology.FlamingoBaryonResponseEmulator_cosmology import FlamingoBaryonBoostMixin

BACCOemuFLAMINGO = with_baryon_boost(BACCOemuNonLinearPerturbations, FlamingoBaryonBoostMixin)
pert = BACCOemuFLAMINGO(
    background, linear_pert, redshifts,
    nonlinear_model_name="Arico2023",
    baryon_kwargs=dict(fgas_sigma=0.0, Mstar_sigma=0.0, jet_fraction=0.0),
)

The resulting object is fully Perturbations-compatible and can be passed directly to any cloelib tracer or likelihood.

Total vs. cold matter boost

The same suppression \(B(k,z)\) is applied to matter_power_spectrum and to matter_power_spectrum_cb. The backends do not all define the response the same way — FLAMINGO's is \((b + \mathrm{cdm} + \nu)/(\mathrm{cdm} + \nu)\) while BACCOemu's is the cold \((b + \mathrm{cdm})/(\mathrm{cdm} + \nu)\) — and one could refine this by adding the neutrino contribution at linear order. The difference is ~\(10^{-5}\) even for \(m_\nu = 0.4\) eV, and FLAMINGO simulations show total and cold boosts to be equivalent at that level, so cloelib treats them as the same.

Return shape

All matter_power_spectrum methods (including baryonic ones) return shape (n_z, n_k), also for a single redshift, which gives (1, n_k). baryonic_suppression returns the same shape, so the boosted spectra keep it.

Adding Your Own Perturbations Implementation¶

To add a new Perturbations implementation, follow these steps.

Step 1: Create Your Class¶

Create cloelib/cosmology/my_solver_cosmology.py:

from cloelib.cosmology.cosmology import Perturbations, Background
import numpy as np

class MySolverPerturbations:
    """Interface to MySolver for perturbation calculations."""

    def __init__(
        self,
        background: Background,
        linearperturbations: Optional[object] = None,
        redshifts: np.ndarray = None,
        nonlinear_model: str = "halofit",
        kmax: float = 10.0,
        zmax: float = 5.0,
        **kwargs
    ):
        """
        Initialize perturbation calculator.

        Args:
            background: Background object (any implementation)
            linearperturbations: Linear perturbations object. Accepted for interface
                compatibility with cloelike, which always passes this as the second
                positional argument when constructing NonLinPerturbations. Unused by
                codes that compute nonlinear corrections internally (e.g. CAMB, CLASS).
            redshifts: Array of redshifts for the calculations.
            nonlinear_model: Which non-linear model to use
            kmax: Maximum wavenumber in 1/Mpc
            zmax: Maximum redshift
        """
        self._background = background
        self.nonlinear_model = nonlinear_model

        # Initialize your solver
        self._solver = MySolver(
            H0=background.H0,
            Omega_b=background.Omega_b0,
            # ... pass cosmological parameters from background
        )

        # Set up k and z grids
        self._solver.set_k_grid(kmin=1e-4, kmax=kmax, nk=200)
        self._solver.set_z_grid(zmin=0.0, zmax=zmax, nz=50)

        # Compute power spectrum
        self._solver.compute_power_spectrum(nonlinear=nonlinear_model)

    @property
    def background(self) -> Background:
        """Return the background object."""
        return self._background

    def matter_power_spectrum(self, zs: np.ndarray, ks: np.ndarray) -> np.ndarray:
        """
        Compute the total matter power spectrum P_mm(k, z).

        Args:
            zs: Redshifts (1D array)
            ks: Wavenumbers in 1/Mpc (1D array)

        Returns:
            P(k, z) in Mpc³, shape (len(zs), len(ks))
        """
        # Ensure inputs are arrays
        zs = np.atleast_1d(zs)
        ks = np.atleast_1d(ks)

        # Call solver and interpolate
        P_k_z = self._solver.get_power_spectrum(zs, ks)

        return P_k_z

    def matter_power_spectrum_cb(self, zs: np.ndarray, ks: np.ndarray) -> np.ndarray:
        """
        Compute the CDM+baryons matter power spectrum P_cb(k, z).

        Args:
            zs: Redshifts (1D array)
            ks: Wavenumbers in 1/Mpc (1D array)

        Returns:
            P_cb(k, z) in Mpc³, shape (len(zs), len(ks))
        """
        zs = np.atleast_1d(zs)
        ks = np.atleast_1d(ks)

        # If the solver distinguishes cb from total matter, use the dedicated method;
        # otherwise fall back to the total matter power spectrum
        P_cb_k_z = self._solver.get_cb_power_spectrum(zs, ks)

        return P_cb_k_z

    def growth_factor(self, zs: np.ndarray, ks: np.ndarray) -> np.ndarray:
        """Calculate linear growth factor D(z, k)."""
        zs = np.atleast_1d(zs)
        ks = np.atleast_1d(ks)

        # Most codes compute scale-independent D(z)
        # If scale-dependent, use ks; otherwise ignore
        D_z = self._solver.get_growth_factor(zs)

        # Return shape (len(zs), len(ks)) even if scale-independent
        return D_z[:, np.newaxis] * np.ones((len(zs), len(ks)))

    def growth_rate(self, zs: np.ndarray | None = None, ks: np.ndarray | None = None) -> np.ndarray:
        """Calculate growth rate f(z) = d ln D / d ln a."""
        if zs is None:
            raise ValueError("zs must be provided")

        zs = np.atleast_1d(zs)
        f_z = self._solver.get_growth_rate(zs)

        # If ks provided, expand to 2D array
        if ks is not None:
            ks = np.atleast_1d(ks)
            f_z = f_z[:, np.newaxis] * np.ones((len(zs), len(ks)))

        return f_z

    def sigma8_0(self) -> float:
        """Compute σ₈ at z=0."""
        return self._solver.get_sigma8()

Step 2: Connect to Background¶

The key pattern: your Perturbations wraps around a Background object.

# Users do this:
bg = CAMBBackground(...)  # Or any Background implementation
pert = MySolverPerturbations(background=bg)

# Now pert can access:
H_z = pert.background.hubble_parameter(z)
chi = pert.background.comoving_distance(z)

This means implementations can be mixed and matched — for example, using a CAMB background with a custom perturbations class.

Step 3: Handle Array Shapes¶

Pay attention to array shapes—it's easy to get confused! cloelib will always return arrays!

def matter_power_spectrum(self, zs, ks):
    """
    Output shape: (len(zs), len(ks))

    Example:
        zs = [0.0, 0.5, 1.0]  # 3 redshifts
        ks = [0.1, 0.2, ...]  # 100 wavenumbers
        P_k_z has shape (3, 100)
    """
    pass

Step 4: Add Tests¶

Create tests/test_my_solver_cosmology.py:

import pytest
import numpy as np
from cloelib.cosmology.camb_cosmology import CAMBBackground
from cloelib.cosmology.my_solver_cosmology import MySolverPerturbations

@pytest.fixture
def background():
    """Create a test background."""
    return CAMBBackground(
        H0=67.5,
        Omega_b0=0.0492,
        Omega_cdm0=0.2650,
        As=2.1e-9,
        ns=0.965,
        # ... minimal parameters
    )

def test_initialization(background):
    """Test Perturbations initialization."""
    pert = MySolverPerturbations(background=background)
    assert pert.background is background

def test_matter_power_spectrum(background):
    """Test P(k, z) calculation."""
    pert = MySolverPerturbations(background=background)

    z = np.array([0.0, 0.5, 1.0])
    k = np.logspace(-2, 1, 50)
    P_k_z = pert.matter_power_spectrum(z, k)

    # Check shape
    assert P_k_z.shape == (len(z), len(k))

    # Check positivity
    assert np.all(P_k_z > 0)

    # Check power decreases at high z (structure less developed)
    assert np.all(P_k_z[0, :] > P_k_z[-1, :])

def test_sigma8(background):
    """Test σ₈ calculation."""
    pert = MySolverPerturbations(background=background)
    sigma8 = pert.sigma8_0()

    # Reasonable range for σ₈
    assert 0.5 < sigma8 < 1.2

def test_growth_factor(background):
    """Test growth factor calculation."""
    pert = MySolverPerturbations(background=background)

    z = np.array([0.0, 1.0, 2.0])
    k = np.array([0.1, 1.0])
    D = pert.growth_factor(z, k)

    # Check shape
    assert D.shape == (len(z), len(k))

    # Check D decreases with z (earlier = less grown)
    assert np.all(D[0, :] > D[1, :])
    assert np.all(D[1, :] > D[2, :])

Also add the new classes to the shared interface tests:

tests/test_perturbations_interface.py: add a builder to BACKENDS, so the class is checked at runtime against the common growth_rate convention.
tests/typing/perturbations_conformance.py: add the class, so ty checks that it implements Perturbations and that its constructor matches LinearPerturbationsFactory/NonLinearPerturbationsFactory.

Step 5: Update Documentation¶

Add to docs/api.md
Update supported codes table in README.md
Add usage examples within the playrgound repository

Interface Patterns¶

Working with External Codes¶

Most Perturbations implementations wrap external Boltzmann codes:

class MyPerturbations:
    def __init__(self, background, ...):
        # Extract parameters from background
        params = {
            'h': background.h,
            'Omega_b': background.Omega_b0,
            'Omega_c': background.Omega_cdm0,
            # ... more parameters
        }

        # Initialize external code
        self._external = ExternalCode(**params)
        self._external.compute()

        # Store for later use
        self._z_grid = ...
        self._k_grid = ...
        self._P_k_z_grid = ...

Tips & Tricks¶

Non-linear vs. Linear Power Spectra¶

Make it clear what you are computing linear or non-linear power spectra at the level of the protocol. The protocols should always return this method and their corresponding CDM and baryon only matter power spectrum _cb.

def matter_power_spectrum(self, zs, ks):
    """
    Compute matter power spectrum.
    """
    pass

Parameter Validation¶

Check inputs early:

def matter_power_spectrum(self, zs, ks):
    if np.any(ks <= 0):
        raise ValueError("Wavenumbers must be positive")
    if np.any(zs < 0):
        raise ValueError("Redshifts must be non-negative")
    if np.any(ks > self.kmax):
        raise ValueError(f"k exceeds kmax={self.kmax}")

Next Steps¶

Now that you have a grounding in Perturbations, explore:

Observables – Connect structure to survey measurements
Summary Statistics – Compute final data products
Background – Review the foundation
API Reference – Full technical details

Perturbations: Structure Formation¶

Overview¶

The Perturbations Protocol¶

Required Property¶

Required Methods¶

matter_power_spectrum(zs, ks)¶

matter_power_spectrum_cb(zs, ks)¶

cb x matter cross spectrum: cb_cross_prescription¶

growth_factor(zs, ks)¶

growth_rate(zs=None, ks=None)¶

sigma8_0()¶

Constructors¶

Existing Implementations¶

CAMB Perturbations¶

CLASSPerturbations¶

HMCode2020EmuPerturbations¶

EE2Perturbations¶

BACCOemuPerturbations¶

FlamingoBaryonResponseEmulatorPerturbations¶

ReACTEmu / BoostedPerturbations¶

TabulatedBoost / TabulatedBoostedPerturbations¶

CosmoPowerJAXPerturbations¶

Available classes¶

Selecting the nonlinear prescription¶

Parameter ranges¶

Example¶

JAXPerturbations¶

MG CosmoPower-JAX boost perturbations¶

MGCLASSPerturbations¶

hi_classPerturbations¶

EmantisPerturbations¶

Baryonic Boost Mixin Framework¶

BaryonBoostMixin Protocol¶

with_baryon_boost Factory¶

Adding Your Own Perturbations Implementation¶

Step 1: Create Your Class¶

Step 2: Connect to Background¶

Step 3: Handle Array Shapes¶

Step 4: Add Tests¶

Step 5: Update Documentation¶

Interface Patterns¶

Working with External Codes¶

Tips & Tricks¶

Non-linear vs. Linear Power Spectra¶

Parameter Validation¶

Next Steps¶

`matter_power_spectrum(zs, ks)`¶

`matter_power_spectrum_cb(zs, ks)`¶

cb x matter cross spectrum: `cb_cross_prescription`¶

`growth_factor(zs, ks)`¶

`growth_rate(zs=None, ks=None)`¶

`sigma8_0()`¶

`BaryonBoostMixin` Protocol¶

`with_baryon_boost` Factory¶