Linear perturbations

solve_perturbations integrates the full scalar Boltzmann hierarchy for one Fourier mode: photons (temperature and polarization), massless neutrinos, baryons, cold dark matter, momentum-resolved massive neutrinos, and — when the cosmology contains them — decaying dark matter with its dark-radiation product, and dark-energy fluid perturbations. The integration is in the conformal Newtonian gauge with the metric potential $\phi$ evolved from the energy constraint (numerically stable outside the horizon, where the algebraic alternative cancels catastrophically) and $\psi$ obtained from the anisotropic stress.

Two design choices differ from the standard codes. There is no tight-coupling approximation: the pre-recombination era, where the Thomson rate exceeds every other scale in the problem, is handled by a stiff implicit integrator instead of an analytic expansion, so there is no switching error to tune. And the hierarchy truncations were set by convergence against CAMB rather than tradition — the defaults $\ell_{\max,\gamma} = 25$, $\ell_{\max,\nu} = 32$ bring $P(k)$ to 0.04% between $k = 0.1$ and $2\,\mathrm{Mpc}^{-1}$, where the textbook $\ell_{\max} = 12$ leaves a 2% error from power reflected off the top of the tower.

The derivations of the evolution equations, including the curved-space generalizations, are collected in Perturbation equations.

Initial conditions

Each mode starts when it is both far outside the horizon and deep in radiation domination — a per-mode start time, since a fixed early start would be wasted effort at large scales and wrong at small ones. The adiabatic growing mode is normalized to unit primordial curvature $\mathcal{R} = 1$; the system being linear, the primordial spectrum multiplies afterwards.

The four regular isocurvature modes are implemented alongside: cold-dark-matter (:cdi), baryon (:bi), neutrino-density (:nid) and neutrino-velocity (:niv) isocurvature, each normalized to unit initial entropy. The super-horizon solutions are the Bucher–Moodley–Turok expansions written in the synchronous gauge and rotated to Newtonian through the Einstein constraints. Pure-mode CMB spectra reproduce the characteristic isocurvature acoustic-peak positions of CLASS (first CDI peak at $\ell \approx 300$ where the adiabatic peak sits at 220), and correlated mixtures are handled by PrimordialMatrix (see Primordial spectra).

Dark-energy perturbations

Any non-$\Lambda$ dark energy is perturbed automatically. The fluid equations with a rest-frame sound speed $\hat c_s^2$ (default 1) are evolved in the variables $(\delta, Q \equiv (1+w)\theta)$, which are algebraically identical to the textbook pair away from $w = -1$ but remain regular there — necessary for quintessence, which starts frozen at exactly $w = -1$. For a canonical scalar field these fluid equations with $\hat c_s^2 = 1$ are the exact field perturbations, no Klein–Gordon block required.

Parameterized fluids that genuinely cross $w = -1$ (the region current baryon-acoustic-oscillation data mildly favour) leave the single-fluid description ill-defined at the crossing; there the PPF construction of Fang, Hu and Lewis takes over automatically (de_scheme = :auto), evolving one regular variable $\Gamma$ matched to scalar-field behaviour on both sides. The crossing response of the CMB agrees with CLASS's PPF implementation.

Gauges

The solution is stored in the Newtonian gauge, and gauge_shift, δ_species, θ_species translate to synchronous, comoving, uniform-density and spatially-flat gauges; the gauge-invariant curvature perturbations curvature_ℛ, curvature_ζ and the Bardeen potentials are provided directly. $\mathcal{R}$ conservation outside the horizon is part of the test suite.

Reference

Cosmic.MassiveNuGrid — Type
MassiveNuGrid

Gauss-Legendre nodes in comoving momentum q, with the Fermi-Dirac weight folded in.

Massive neutrinos are the one species that cannot be reduced to a fluid. A fluid has a single velocity at each point. A free-streaming gas of massive particles has particles of different momenta moving at different speeds through the same point, and that spread in velocity is precisely what erases small-scale structure. So the perturbation has to be tracked as a function of momentum – a whole hierarchy per momentum bin – which is why this is expensive and why it is the last thing to get implemented in any Boltzmann code.

source
Cosmic.initial_conditions — Method
initial_conditions(c, bg, k, L, x0; ic = :adiabatic)

Adiabatic growing mode, deep in radiation domination and far outside the horizon, normalised to a primordial curvature perturbation ℛ = 1. Since the system is linear, the primordial spectrum is applied afterwards rather than here.

Derived in docs/src/perturbation_equations.md. With fν = ρν/(ργ+ρν):

ψ = 1/(3/2 + 2f_ν/5),  φ = (1 + 2f_ν/5) ψ
δ_γ = δ_ν = -2ψ,       δ_c = δ_b = -(3/2)ψ
θ (all species) = (1/2) k² η ψ
σ_ν = (1/15) k² η² ψ
source
Cosmic.solve_perturbations — Method
solve_perturbations(c, bg, k; lmax_γ = 12, lmax_ν = 12)

Integrate the Boltzmann hierarchy for one mode k (1/Mpc).

Defaults for the hierarchy truncations are set by convergence against CAMB, not by tradition: ℓmaxγ = 12 (a common textbook choice) leaves a 2% error in P(k) at k = 1/Mpc, because power reflected off the top of the truncated hierarchy leaks back into the low multipoles that source the metric. ℓmaxγ = 25, ℓmaxν = 32 brings P(k) to 0.04% of CAMB across 0.1 < k < 2.

The start time is chosen per mode so that the mode is comfortably outside the horizon (kη ≈ 1e-2) and deep in radiation domination – both are what make the adiabatic initial conditions valid. A fixed early start for every k would waste effort on large scales and be outright wrong for small ones.

Before recombination the system is stiff: κ̇ exceeds k and ℋ by orders of magnitude, which is exactly what "tight coupling" means. Boltzmann codes usually swap in an analytic tight-coupling expansion at that point for speed; we instead pay for an implicit solver, which is slower but has no extra approximation to get wrong.

Pass mg = stable_basis_solve(c, spec) to evolve linear Horndeski gravity (αT = 0) alongside the matter hierarchy: the scalar's velocity potential vX gets its own second-order equation of motion (Bellini & Sawicki 1404.3713, reduced), and the modified Einstein constraints feed the altered metric potentials φ, ψ back into every species. Omit it (the default) and the GR code path is untouched, bit for bit. Massive neutrinos and decaying CDM are supported: their density, pressure, momentum and (time-derivative of) anisotropic-stress moments — plus, for decaying CDM, the dark radiation's decay-sourced background pressure derivative — enter the scalar EOM. The B&S 3.17 Hamiltonian constraint, which the system does not evolve, stays satisfied along the solution to the same level as the massless GR limit under the same model (≈1e-5 relative for a physically-motivated spec; the residual is set by the model's activation transient, not by the massive-ν or dcdm sources). The MG branch is flat space only and still errors on a dark-energy fluid.

source
Cosmic.Φ — Method
Φ(p, a)

Spatial curvature perturbation. Constant outside the horizon in both radiation and matter domination, and famously drops by a factor 9/10 across the transition between them.

source
Cosmic.Ψ — Method
Ψ(p, a)

Newtonian potential. Equal to Φ up to the anisotropic stress, which is carried by free-streaming neutrinos and, after decoupling, photons.

source
Cosmic.δ_matter — Method
δ_matter(p, a)

Total matter density contrast in the conformal Newtonian gauge, ρ-weighted over CDM and baryons. Baryons are ~16% of the matter and do not trace CDM inside the sound horizon, so this is not δ_cdm.

For the matter power spectrum you almost certainly want δ_matter_comoving instead – see the warning there.

source
Cosmic.δ_matter_comoving — Method
δ_matter_comoving(p, a)

Total matter density contrast in the comoving (synchronous) gauge,

δ_sync = δ_Newtonian + 3ℋ θ_c / k²

This – not the Newtonian-gauge δ – is what "the matter power spectrum" means by convention, and what CAMB and CLASS report.

The distinction is not academic. A density contrast is gauge dependent, and the two definitions agree only well inside the horizon, where the correction term 3ℋθ/k² is negligible. As k approaches ℋ they diverge badly: at k = 5e-4 /Mpc the Newtonian-gauge P(k) is 70% larger than the comoving one. Using the wrong one looks fine at k ≳ 0.01 and is quietly, badly wrong on the largest scales – the ones surveys are increasingly trying to measure.

The synchronous gauge is defined by CDM being at rest, which is why the correction is built from θ_c specifically.

source
Cosmic.δ_matter_total_comoving — Method
δ_matter_total_comoving(p, a)

Total-matter density contrast in the comoving (synchronous) gauge, summing cold dark matter, baryons and massive neutrinos, weighted by their physical densities at a. For a massless-ν cosmology this is identical to δ_matter_comoving; the two diverge only when Σm_ν > 0, where the free-streaming neutrinos suppress the total-matter spectrum by a few percent below the cold (cdm+baryon) one on small scales.

The synchronous correction 3ℋθc/k² is applied to every component with the proper (1+w) weight, so it reduces exactly to the cold expression when there are no massive neutrinos. Use [`δmatter_comoving`](@ref) (cold) for the halo mass function and σ(R); this total version is what the HMcode two-halo term wants.

source
Cosmic.bardeen_Φ — Method
bardeen_Φ(p, a), bardeen_Ψ(p, a)

The Bardeen potentials. These are gauge-invariant by construction, and they happen to coincide with the metric potentials in the conformal Newtonian gauge – which is the deeper reason that gauge is so convenient: its variables are already the gauge-invariant ones, so nothing is hiding in the coordinates.

source
Cosmic.curvature_ζ — Method
curvature_ζ(p, a)

The uniform-density curvature perturbation,

ζ = φ - δρ_tot / (3(ρ̄ + P̄))

Also gauge invariant, and equal to ℛ on super-horizon scales (they differ only by terms of order (k/ℋ)²). ζ is the more natural variable when the equation of state is changing, ℛ when it is not; both are conserved outside the horizon.

source
Cosmic.curvature_ℛ — Method
curvature_ℛ(p, a)

The comoving curvature perturbation,

ℛ = φ + ℋ θ_tot / k²

Gauge invariant, and conserved on super-horizon scales for adiabatic perturbations – which is what lets inflation predict an initial condition at all, since ℛ set during inflation survives untouched through reheating, through equality, right up to horizon re-entry.

Cosmic normalises the initial conditions to ℛ = 1 per mode, so on large scales this should return 1 at every time. That is not a tautology: ℛ is reconstructed from the evolved φ and θ_tot, so it only comes out to 1 if the hierarchy and the Einstein equations are consistent with each other.

source
Cosmic.gauge_shift — Method
gauge_shift(p, a, gauge)

The time shift T that takes the conformal-Newtonian solution into gauge. Returns T in Mpc (the same units as conformal time).

source
Cosmic.δ_species — Function
δ_species(p, a, which, gauge)

Density contrast of one species in the requested gauge. which is one of :cdm, :baryon, :photon, :neutrino, :matter.

source
Cosmic.θ_species — Function
θ_species(p, a, which, gauge)

Velocity divergence of one species in the requested gauge.

source