The matter power spectrum

MatterPowerSpectrum

Tabulated transfer function δm(k) at one redshift, ready to be turned into P(k). Built by [`matterpower_spectrum`](@ref).

bao_k_grid(rec, kmin, kmax; per_decade = 10, per_decade_bao = 70,
           bao_center = 3.0, bao_width = 4.0)

Wavenumber grid with the sampling density boosted across the BAO region.

A uniform grid in ln k is the obvious choice and it is wrong, quietly. The baryon acoustic oscillations in P(k) are periodic in linear k with period Δk = 2π/rs ≈ 0.04 /Mpc, so a log grid that comfortably resolves the turnover at keq gives only two or three points per wiggle by k ~ 0.2, and the interpolated spectrum comes out visibly kinked – the wiggles get aliased, not smoothed.

CLASS solves this by varying the density (perturbations.c, and the defaults in precisions.h): 10 points per decade in general, rising to 70 per decade inside a window centred on 3·krec with logarithmic width 4, where krec = 2π/rs(z*). That window is exactly where the oscillations live, so the cost is paid only where it buys something.

dln k = ln10 / [ n_pk + (n_bao - n_pk)·(1 - tanh( ((ln k - ln(3 k_rec))/ln 4)^4 )) ]

dimensionless_power(P, k)

Δ²(k) = k³P(k)/(2π²) = δm(k)² Δ²ℛ(k). The variance of the matter field per logarithmic interval in k, which is the form σ₈ is an integral of.

growth_factor(c, rec, z; k = 1e-3)

Linear growth factor D(z), normalised to D(0) = 1.

Measured from an actual mode rather than from the usual scale-independent growth ODE. That matters: with massive neutrinos (or any component that free-streams) growth is scale dependent, and the textbook ODE quietly assumes it is not. The default k = 1e-3 /Mpc is large-scale enough to recover the scale-independent answer in plain ΛCDM.

growth_rate(c, rec, z; k = 1e-3)

f = dln D / dln a, the quantity redshift-space distortions measure. In ΛCDM it is close to Ω_m(a)^0.55.

matter_power_spectrum(c, rec; z = 0, kmin = 1e-4, kmax = 10, nk = nothing, total = false)

Solve the hierarchy on a grid of wavenumbers and tabulate the matter transfer function at redshift z. k is in 1/Mpc.

By default this returns the cold (cdm+baryon) spectrum, the convention CAMB and CLASS report and the input the halo mass function and σ(R) want. Pass total = true for the total-matter spectrum that also includes massive-neutrino density perturbations (identical to cold for massless ν); HMcode's two-halo term uses this one.

With nk = nothing (the default) the grid comes from bao_k_grid, which refines itself across the acoustic oscillations. Passing an integer nk forces a plain log-spaced grid of that size instead – simpler, faster, and liable to alias the BAO wiggles, so it is for quick looks rather than for results.

One Boltzmann solve per mode: this is the expensive step, and it is embarrassingly parallel, so the modes are threaded.

power(P, k)

Matter power spectrum P(k) in Mpc³:

P(k) = δ_m(k)² · A_s (k/k_pivot)^(n_s-1) · 2π²/k³

primordial_power(c, k)

Dimensionless primordial curvature power spectrum Δ²ℛ(k). By default the power law Δ²ℛ(k) = As (k/kpivot)^(ns − 1), but if the cosmology carries a primordial override (any callable `k ↦ Δ²ℛ(k)`, e.g. an inflaton solution or a tabulated broken-scale-invariance spectrum) that is used instead. The override lets features, steps, running, or a computed inflationary spectrum enter without touching the transfer-function machinery, since the spectrum is applied linearly on top of the ℛ = 1 transfer functions.

top_hat_window(x)

Fourier transform of a real-space spherical top hat, W(x) = 3(sin x - x cos x)/x³, with W(0) = 1.

transfer(P, k)

Matter transfer function δ_m(k) at the spectrum's redshift, ℛ = 1 normalised.

σ8(P)

σ₈: the RMS matter fluctuation in spheres of radius 8 h⁻¹ Mpc. Note the radius is in h⁻¹ Mpc while the code works in Mpc, so R = 8/h.

σ_R(P, R)

RMS linear matter fluctuation in spheres of comoving radius R (Mpc):

σ_R² = ∫ dln k  Δ²(k)  W(kR)²

HalofitSpectrum

Non-linear P(k) from Takahashi et al. (2012). Callable via power: power(hf, k) returns the non-linear P(k) in Mpc³. Also carries the linear spectrum it was built from and the non-linear-scale diagnostics.

Dimensionless non-linear power Δ²(k) = k³P(k)/2π².

halofit_power(P)

Build the Takahashi 2012 + Bird 2012 non-linear spectrum on top of a linear MatterPowerSpectrum P. The linear spectrum must reach well past the non-linear scale (kmax ≳ 10/Mpc) so the Gaussian variance integral converges.

nonlinear_scale(P)

Return (k_σ, n_eff, C) — the non-linear wavenumber (1/Mpc), the effective slope, and the curvature of the linear spectrum at the non-linear scale. k_σ is found from σGauss(1/kσ) = 1 by bisection in log R, exactly as CLASS does. Throws if the linear spectrum does not reach far enough into the non-linear regime (σ never reaches 1 over the available k-range).

power(hf, k)

Non-linear matter power spectrum P(k) in Mpc³ at wavenumber k (1/Mpc). Below the non-linear scale it returns the linear spectrum; above it, the Takahashi halo (one-halo) plus quasi-linear (two-halo) sum.

hmcode_power(P; P_total = P)

Build the HMcode-2020 non-linear spectrum on top of a linear cold (cdm+baryon) MatterPowerSpectrum P (which must reach kmax ≳ 10/Mpc). Evaluate with power: power(hm, k) returns P_nl(k) in Mpc³.

The entire halo model — σ(R), the mass function, the one- and two-halo terms — runs on the cold spectrum P, exactly as HMcode-2020 prescribes (neutrinos do not cluster into haloes). For a massive-ν cosmology the total-matter nonlinear spectrum is then reconstructed from the cold halo-model result by the linear total/cold ratio,

P_mm(k) = P_cb^NL(k) · [P_mm^lin(k) / P_cb^lin(k)],

because the neutrinos stay linear and carry that scale-dependent suppression (→ fcb² at high k) through the nonlinear evolution unchanged. Pass `Ptotalfrommatterpowerspectrum(...; total = true)to activate it; for massless νP_total = Pmakes the ratio 1 and the result is the pure cold spectrum. Validated against CAMBmead2020` (Σmν = 0.15 eV) to <0.1% across 0.05 < k[h] < 8.

power(hm, k)

HMcode-2020 non-linear matter power spectrum P(k) in Mpc³, k in 1/Mpc.

eft_power(P; cs2 = 0.0, ir_resum = true, rd = 100/h Mpc)

Build the 1-loop EFT spectrum on a linear MatterPowerSpectrum. cs2 is the effective sound-speed counterterm in Mpc² (a free parameter — fit it to data; 0 gives pure 1-loop SPT). ir_resum applies BAO IR resummation. Evaluate with power.

p13(P, k; qmax_factor = 100)

The P₁₃ one-loop integral (Mpc³) at wavenumber k (1/Mpc), built from the linear MatterPowerSpectrum P. Negative, as it should be at high k.

p22(P, k)

The P₂₂ one-loop integral (Mpc³) at k (1/Mpc). Positive.

power(eft, k)

1-loop EFT matter power spectrum P(k) in Mpc³ at k (1/Mpc). With IR resummation the BAO wiggle in both the linear and loop pieces is damped by e^{-k²Σ²}; the counterterm -2 cs² k² Plin is always included.