Recompute a Primordial Spectrum and Test the Solver

Predict before you run

Write down what you expect. For the $R^2$ plateau with $N_*=55$, the lowest-order formulas give $n_s\approx0.9636$ and $r\approx0.004$. For $\tfrac12m^2\phi^2$ they give $n_s\approx0.9640$ and $r\approx0.144$.

Also predict the size of the step’s effect. A step placed $59$ e-folds before the end should disturb modes about four e-folds larger than the pivot, near $k\approx10^{-3}\,\mathrm{Mpc}^{-1}$, and leave the pivot alone.

Reading: Planck Collaboration, Planck 2018 results. X. Constraints on inflation, arXiv:1807.06211v2 · Original background and mode calculations with an independent Python reference

Get the files and fix the environment

Download the ZIP below and extract it into a fresh directory. It contains the reference solver, the benchmark script, a comparison script, five input models and their expected outputs. There is no browser Engine in it.

Install the pinned versions in a virtual environment. The expected files were made with Python 3.9.6, NumPy 2.0.2 and SciPy 1.13.1, and every output records the versions that produced it.

python3 -m venv venv
. venv/bin/activate
pip install -r requirements.txt

Reading: Original background and mode calculations with an independent Python reference · CLASS explanatory.ini, lesgourg/class_public commit c5d93d01990116162b408b789b62ec5669fc4ea4

Test the solver on an exact answer first

Power-law inflation, $V=V_0e^{-\lambda\phi}$, has constant $\epsilon=\lambda^2/2$, and its scalar and tensor spectra are known in closed form from Stewart and Lyth’s equations 31, 41 and 43. The benchmark compares the numerical spectra with that formula, not with the slow-roll estimate.

Our run covered $\lambda$ from $0.02$ to $0.3$. The worst power error was $2.0\times10^{-6}$ and the worst Wronskian drift $4.0\times10^{-7}$. The leading slow-roll formula was off by $1.1\times10^{-4}$ at $\epsilon=2\times10^{-4}$ and by $2.6\%$ at $\epsilon=0.045$. That is why the closed form, not slow roll, is the test.

The script exits with an error if any check fails. Move on once it passes on your machine.

python verify_benchmarks.py my-benchmarks.json

Reading: Stewart & Lyth, A more accurate analytic calculation of the spectrum of cosmological perturbations produced during inflation (1993), gr-qc/9302019v1 · Original background and mode calculations with an independent Python reference

Recompute a fixed model

Each input names the potential, the amplitude $V_0$ in units of $M_{\mathrm{pl}}^4$, the starting field, $N_*$, the step and a grid of $27$ wavenumbers in $\mathrm{Mpc}^{-1}$. The step uses the tanh form of Adams, Cresswell and Easther’s equation 1, with its position and width given in field units.

The output holds the background, with total e-folds, pivot field and $\epsilon$ there, then one row per wavenumber. Each row gives the scalar and tensor power, the slow-roll estimate, the exit time and the field at exit. Each preset takes about one second.

On the same software versions our fresh run matched the expected file exactly. A different SciPy may move the last digits; the comparison script fails above $10^{-4}$, which is far larger than the solver error and far smaller than any physical effect discussed here.

python recompute_inflation_spectrum.py inputs/step-ripples.json my-step.json
python compare_with_expected.py my-step.json expected/step-ripples.json

Reading: Adams, Cresswell & Easther, Inflationary perturbations from a potential with a step (2001), astro-ph/0102236v2 · Original background and mode calculations with an independent Python reference

Measure the tilt and tensor ratio yourself

Make a three-point input around the pivot and difference the logarithm of the power. The snippet below does it for the plateau.

Our results: plateau $n_s=0.96490$, $r=0.00355$; with the default step, the same to four decimal places; quadratic $n_s=0.96338$, $r=0.14403$; shorter reheating, $N_*=46$, $n_s=0.95826$, $r=0.00498$.

The amplitude needs care. The plateau’s $V_0$ was chosen so that $A_s=2.0989\times10^{-9}$, matching Planck’s $\ln(10^{10}A_s)=3.044$. With the same $V_0$ the $N_*=46$ model gives $1.4880\times10^{-9}$. Since $\mathcal{P}_{\mathcal{R}}\propto V_0$ at fixed $N_*$, rescaling $V_0$ by $1.4105$ restores it without changing the tilt or $r$.

import json, math, subprocess
model = json.load(open('inputs/smooth-plateau.json'))
model['kGridMpc'] = [0.049, 0.05, 0.051]
json.dump(model, open('tilt-input.json', 'w'))
subprocess.run(['python', 'recompute_inflation_spectrum.py', 'tilt-input.json', 'tilt-output.json'], check=True)
low, mid, high = json.load(open('tilt-output.json'))['spectrum']
print('n_s =', 1 + math.log(high['scalar'] / low['scalar']) / math.log(0.051 / 0.049))
print('r   =', mid['tensor'] / mid['scalar'])

Reading: Original background and mode calculations with an independent Python reference · Planck Collaboration, Planck 2018 results. X. Constraints on inflation, arXiv:1807.06211v2

Compare with published constraints

Planck 2018 X gives $n_s=0.9649\pm0.0042$ and $\ln(10^{10}A_s)=3.044\pm0.014$; PDG 2025 gives $r<0.036$ at $95\%$ confidence. The plateau sits on the centre of the tilt band. The quadratic model fits the tilt but its $r$ is four times the bound. The $N_*=46$ plateau is $1.6\sigma$ low in tilt.

These are comparisons of primordial parameters, which were themselves fitted assuming a power-law spectrum. A model with a step must instead go through a Boltzmann code and the likelihood. CLASS can accept your own table through its external_Pk option, which is the natural next step; this course stops at the primordial spectrum.

Reading: Planck Collaboration, Planck 2018 results. X. Constraints on inflation, arXiv:1807.06211v2 · Ellis, Vennin & Wands, Review 23: Inflation, Review of Particle Physics (PDG 2025, revised August 2025) · CLASS explanatory.ini, lesgourg/class_public commit c5d93d01990116162b408b789b62ec5669fc4ea4

Failure cases we actually hit

Absolute tolerance. An early version used a fixed $\mathrm{atol}=10^{-30}$. For $N_*=46$ the starting $|\mathcal{R}|$ ranges from about $2\times10^{-29}$ to $10^{-24}$ across the grid in these units, so the tolerance was loose. Over the $27$-point grid the power came out up to $0.55\%$ wrong and the Wronskian drifted by up to $0.88\%$, while the integrator reported success. The plateau at $N_*=55$ stayed below $10^{-4}$. The fix scales the tolerance to each component’s starting size.

Step position. With the step crossed $58$ e-folds before the end, its ringing reached the pivot and moved the reported tilt. We moved the default step to $59$ e-folds and now check the pivot tilt against the smooth model.

Pivot choice for $r$. The default step leaves $r$ at $0.05\,\mathrm{Mpc}^{-1}$ unchanged, at $0.00355$, but raises it at $0.002\,\mathrm{Mpc}^{-1}$ from $0.00318$ to $0.00441$. The Lab therefore reports $r$ at $0.05\,\mathrm{Mpc}^{-1}$.

Reading: Original background and mode calculations with an independent Python reference · Planck Collaboration, Planck 2018 results. X. Constraints on inflation, arXiv:1807.06211v2

What this reproduction does not establish

It is one canonical field in Einstein gravity, with a Bunch–Davies vacuum and no dynamics after $\epsilon=1$. It produces primordial spectra only; no transfer function or likelihood is applied, so no fit to data is claimed.

Agreement between the Lab’s integrator and this reference, within $2.5\times10^{-4}$ on all five presets, shows the two codes solve the same equations.

Reading: Original background and mode calculations with an independent Python reference · Ellis, Vennin & Wands, Review 23: Inflation, Review of Particle Physics (PDG 2025, revised August 2025)