Recompute a small network from a frozen rate grid and propagate rate errors

What a frozen rate grid is

A BBN network spends most of its effort on quantities that depend only on temperature: thermonuclear reaction rates, the weak neutron-proton rates, and the energy density of the photon-electron plasma. None of them depends on the baryon density, the extra radiation or the neutron lifetime once the weak rates are written in units of $1/\tau_n$. The package computes them once on a grid of 173 temperatures, $\log_{10}T_9=-2.3,-2.275,\dots,2.0$, and saves them in frozen-rate-grid.json. The network then reads the grid, interpolates in $\log T_9$, and can rerun for any $\eta_{10}$, $\Delta N_{\mathrm{eff}}$ or $\tau_n$ in a fraction of a second. The Lab uses a byte-identical copy of the same grid.

The package is original code under the MIT licence. It contains no third-party code; the Serpico coefficients and the AME2020 masses are transcribed numbers, cited in the files that use them.

Download links: MANIFEST.json, input.json, expected.json, frozen-rate-grid.json, teaching-values.json and provenance.json. The Python files follow the same pattern, one URL per file name.

Units and conventions

Keep these straight before reading any number.

• Temperature: $T_9=T/10^9\,\mathrm K$, and $k_BT=0.08617\,T_9\,\mathrm{MeV}$. The network runs from $T_9=100$ ($8.6\,\mathrm{MeV}$) down to $T_9=0.01$ ($0.86\,\mathrm{keV}$), in steps uniform in $\ln T$.
• Reaction rates: $N_A\langle\sigma v\rangle$ in $\mathrm{cm^3\,mol^{-1}\,s^{-1}}$, stored as $\log_{10}$. Each reaction also stores $\delta f/f$, its fractional one-sigma error.
• Weak rates: $\tau_n\Gamma_{n\to p}$ and $\tau_n\Gamma_{p\to n}$, dimensionless.
• Abundances: $Y_i=n_i/n_b$ for n, p, d, t, ${}^3\mathrm{He}$, ${}^4\mathrm{He}$, ${}^7\mathrm{Li}$ and ${}^7\mathrm{Be}$. Baryon conservation means $\sum_iA_iY_i=1$, and every saved run records that sum.
• Outputs, following PRIMAT: $Y_p=4Y({}^4\mathrm{He})$; $\mathrm{D/H}=Y_d/Y_p$; ${}^3\mathrm{He/H}$ counts helium-3 plus tritium; ${}^7\mathrm{Li/H}$ counts lithium-7 plus beryllium-7, because tritium and beryllium-7 decay later.
• Inputs: $\eta_{10}$ is today's baryon-to-photon ratio times $10^{10}$; $\Delta N_{\mathrm{eff}}$ adds radiation counted in neutrino species on top of three; $\tau_n$ in seconds, default $878.4\,\mathrm s$, the value quoted in Section 24.2 of the 2025 Review of Particle Physics.

Step 1: download and check the files

Put the files in a new, empty folder. Replace the host name if you use a mirror.

mkdir bbn-frozen-rate-network && cd bbn-frozen-rate-network
base=https://sandboxphysics.com/stories/primordial-nucleosynthesis/materials
curl -fsSO "$base/MANIFEST.json"
for name in $(python3 -c "import json; print(' '.join(f['name'] for f in json.load(open('MANIFEST.json'))['files']))"); do
  curl -fsSO "$base/$name"
done
python3 -c "import hashlib, json; bad=[f['name'] for f in json.load(open('MANIFEST.json'))['files'] if hashlib.sha256(open(f['name'],'rb').read()).hexdigest()!=f['sha256']]; print('digests OK' if not bad else bad)"

The last line should print digests OK. The site's build stops if any served file no longer matches its size and digest in MANIFEST.json, so a mismatch means the copy changed after it left the site.

Step 2: replay the saved result

python3 -B replay.py actual-result.json

The run takes about two minutes and prints

wrote actual-result.json
PASS

replay.py rebuilds the grid from the Python modules, compares it with the saved grid, reruns every case in input.json, recomputes the benchmarks, sensitivities, error budget and Monte Carlo, and compares every number with expected.json at a relative tolerance of $10^{-7}$. The saved gridRebuildMaxDifference is exactly $0$: the saved grid is what the modules produce.

The package has been run only on CPython 3.9.6 on macOS. On another platform, the last bits of math.exp can differ, and a stiff solver can amplify that past $10^{-7}$. If you see FAIL with differences in the seventh or eighth significant digit, compare against the physical tolerances in the benchmarks below instead; differences in the third digit mean something is wrong.

input.json holds five cases: the standard universe ($\eta_{10}=6.12$), baryon-poor ($2$), baryon-rich ($15$), extra radiation ($\Delta N_{\mathrm{eff}}=1.5$) and missing radiation ($-1$). Each is run with 2000 steps and with the Lab's 1000 steps. For the standard case expected.json gives:

{
  "heliumMassFraction": 0.2423956078641409,
  "deuteriumToHydrogen": 2.5482022079050184e-05,
  "helium3ToHydrogen": 1.018552513743762e-05,
  "lithium7ToHydrogen": 4.5131264294261034e-10,
  "baryonNumberSum": 0.999999999999992
}

At the Lab's 1000 steps D/H is $2.5570\times10^{-5}$; at 4000 steps it is $2.5456\times10^{-5}$. The time-step error of the Lab is therefore about 0.4 percent on D/H and smaller on the other outputs.

Step 3: read the benchmarks one source at a time

The benchmarks block of expected.json compares the network with three published calculations. Each comparison tests something different; read them in this order.

Serpico et al., Tables 3 and 4

These rows use the same rate fits as the network, so they test the transcription and the network, not the physics. Serpico quotes results against $\omega_b$, with $\tau_n=885.7\,\mathrm s$ and $\eta_{10}=273.49\,\omega_b/(1-0.007\,Y_p)$; input.json records both. For $\omega_b=0.022$, $0.023$ and $0.024$ the network reproduces D/H within 0.4 percent, helium-3 within 0.8 percent and lithium-7 within 1.1 to 2.5 percent, below. A common misreading is to take this agreement as proof that the network is accurate. It shows that the fits were copied and integrated correctly; their physics is the 2004 physics.

PRIMAT, Table V, row Born

PRIMAT lists how each correction to the weak rates moves the abundances, starting from the Born approximation, the level the tutorial implements. At PRIMAT's inputs, $\eta_{10}=6.0913$ and $\tau_n=879.5\,\mathrm s$, the network gives $Y_p=0.24258$ against PRIMAT's $0.24262$, a difference of $-0.018$ percent. This tests the weak rates, the plasma thermodynamics and the clock. Only helium is compared. PRIMAT uses its own nuclear rates, so its D/H of $2.423\times10^{-5}$ and lithium of $5.635\times10^{-10}$ are not targets for a network built on the 2004 fits; the network gives $2.569\times10^{-5}$ and $4.47\times10^{-10}$. The same table shows the corrections the tutorial omits raise $Y_p$ to $0.24709$.

Fields et al., Table II

Fields lists logarithmic sensitivities at $\eta_{10}=6.129$ and $\tau_n=879.4\,\mathrm s$. The network reproduces all 60 entries for the inputs it shares with that table within $0.036$; the largest difference is the lithium sensitivity to $\eta$, $2.058$ against $2.094$. Sensitivities test the shape of the network's response rather than its absolute values, and they are what the error budget uses.

Two further checks compare the helium slope per extra neutrino species. The network gives $0.01341$; the PRIMAT-based table distributed with the CLASS code, sBBN_2025_primat.dat, gives $0.01345$, and the review quotes $0.013$.

Step 4: propagate the rate errors

propagate_rate_uncertainties.py treats each rate error as log-normal, fully correlated across temperature within one reaction and independent between reactions. For a standard normal number $z$ the rate becomes

The one-sigma budget sets $z=+1$ and $z=-1$ for one reaction at a time and records half the change in $\ln X$. The neutron lifetime is moved by $\pm0.5\,\mathrm s$. The quadrature total at $\eta_{10}=6.12$ is 0.69 percent for D/H, 0.75 percent for helium-3 and 1.76 percent for lithium-7. Helium-4 moves by 0.04 percent, almost all from the lifetime.

The Monte Carlo draws 200 sets of $z$ values with the fixed seed 20261112, so it reproduces exactly. Its relative standard deviations, 0.72 percent for D/H and 1.78 percent for lithium-7, agree with the quadrature totals. That agreement tests the linear approximation; with 200 samples the standard deviation itself carries a statistical error of about 5 percent.

To compute a budget at another baryon density, use the module directly:

from network_solver import FrozenGrid, load_grid
from propagate_rate_uncertainties import one_sigma_budget

grid = FrozenGrid(load_grid('frozen-rate-grid.json'))
budget = one_sigma_budget(grid, {'eta10': 5.0, 'neutron_lifetime': 878.4, 'steps': 2000})
print(budget['quadratureTotal'])

Step 5: regenerate the course numbers

python3 -B teaching_values.py my-teaching-values.json

This takes about fifteen seconds and should reproduce teaching-values.json: the freeze-out point at $T_9=7.71$, the clock from $T_9=10$ to $0.1$, helium half formed at 247 seconds, the equilibrium deuterium temperatures, the comparison runs and the baryon density inferred from the observed deuterium, $\eta_{10}=6.180$.

Failure cases worth reproducing

The failureCases block keeps the mistakes that this package had to avoid. Each one is a test you can run on another code.

• Raw fits outside their range. Evaluated on the full grid, every Serpico fit turns non-positive somewhere, between 6 and 65 of the 173 grid points depending on the reaction. The fits are only valid over their stated temperatures. rate_fit_windows.py continues each one below its window with the Coulomb-barrier form and holds it constant above. Moving every lower edge by a factor of $0.75$ or $1.33$ changes the outputs by less than $2\times10^{-4}$, so the result does not depend on that choice.
• A printed sign. Read literally, the printed $d(p,\gamma){}^3\mathrm{He}$ fit has a growing rather than a decaying barrier factor and gives $1452$ at $T_9=0.1$ instead of $5.59$, in the units above. The package uses the decaying form and records the check in serpico_rate_fits.py and provenance.json; a second sign, in the lower branch of the ${}^3\mathrm{He}(\alpha,\gamma){}^7\mathrm{Be}$ error fit, is treated the same way.
• A first-order time step. Replacing the step-doubled backward Euler scheme by plain backward Euler at the same 2000 steps shifts D/H by $+6.4$ percent and lithium-7 by $-8.5$ percent against a 4000-step reference. Agreement on helium alone would hide this.
• Too few steps. With 250 steps D/H is $2.651\times10^{-5}$, 4 percent high. Always rerun at double the steps before quoting a fourth digit.

What the package does not include

It stops at the Born approximation for the weak rates and decouples neutrinos instantly, so it underestimates $Y_p$ by about $0.005$ compared with full calculations. It uses the 2004 nuclear rates throughout, so its D/H at fixed $\eta$ is not the current best prediction. It has no QED corrections to the plasma and treats extra radiation as a constant $\Delta N_{\mathrm{eff}}$ present from the start. These limits matter when you compare with observations. Each check above compares like with like, so those checks still hold.

Checklist for the next network you meet

• Which rate compilation, which weak-rate corrections, and which neutrino decoupling treatment does it use?
• Does it conserve baryon number to rounding error, and does it record the check?
• Does it agree with its own source tables, and separately with an independent calculation at matched inputs?
• Has it been rerun at double the time steps?
• Are rate errors propagated per reaction with temperature dependence, and is the propagation checked by a seeded Monte Carlo?
• Are the absolute offsets of the model stated before any bound on new physics is quoted?