---
title: 'PRyMordial: Precision BBN Calculations'
url: https://www.emergentmind.com/topics/prymordial
type: topic
---

# PRyMordial: Precision BBN Calculations

PRyMordial is a public numerical package for precision calculations of Big Bang Nucleosynthesis (BBN) observables within the Standard Model and in a broad class of beyond-the-Standard-Model scenarios. In its original formulation, it was designed to compute primordial light-element abundances together with \(N_{\rm eff}\), including non-instantaneous neutrino decoupling effects, while remaining fast enough for Monte Carlo and Bayesian parameter inference and flexible enough to accommodate modified thermal histories, weak rates, and nuclear networks [2307.07061]. Subsequent work has used it as both a high-precision Standard BBN engine and a modifiable framework for neutron-lifetime studies, baryon-density inference, varying-coupling analyses, nonstandard expansion histories, CPT-violating plasma physics, and other MeV-era deformations of early-universe microphysics [2210.12031].

## 1. Definition and scientific scope

PRyMordial occupies an intermediate position between precision Standard-BBN codes such as PArthENoPE and PRIMAT and more BSM-oriented but less precision-focused tools such as AlterBBN [2307.07061]. Its core remit is the first few minutes of cosmic evolution, over a temperature interval from \(\mathcal{O}(10)\,\mathrm{MeV}\) down to \(\mathcal{O}(\mathrm{keV})\), with explicit treatment of the radiation bath, weak \(n\leftrightarrow p\) conversion, and the thermonuclear network.

The package returns the principal BBN observables relevant for cosmological inference: the helium-4 mass fraction \(Y_P\), deuterium abundance \(D/H\), helium-3 abundance \({}^3\mathrm{He}/H\), lithium-7 abundance \({}^7\mathrm{Li}/H\), and \(N_{\rm eff}\) including incomplete neutrino decoupling [2307.07061]. It also computes related neutrino quantities such as the present-day relativistic neutrino abundance. This makes it directly usable for internal BBN consistency studies, BBN–CMB comparisons, and parameter inference in scenarios with extra radiation, lepton asymmetry, modified weak interactions, altered nuclear rates, or additional thermal sectors.

A central design choice is that PRyMordial computes the BBN thermal background from first principles rather than relying only on pretabulated Standard-Model backgrounds. This is the main reason it has been repeatedly adopted in later work as a general-purpose numerical laboratory for MeV-scale cosmology, rather than only as a fixed Standard-BBN abundance interpolator [2307.07061].

## 2. Physical and numerical formulation

PRyMordial is organized around three coupled components: the thermodynamics of the radiation bath, the \(n\leftrightarrow p\) weak rates, and the nuclear reaction network [2307.07061]. The baseline expansion is determined by the Friedmann equation
\[
H^2 \equiv \left(\frac{d\log a}{dt}\right)^2 = \frac{8\pi}{3M_{\rm Pl}^2}\rho_{\rm tot},
\]
with \(\rho_{\rm tot}\) built from photons, \(e^\pm\), neutrinos, and optionally new sectors. The code distinguishes the photon/plasma temperature \(T_\gamma\) from the neutrino temperature \(T_\nu\), and evolves them through coupled Boltzmann equations rather than assuming instantaneous decoupling.

The thermal sector is governed by
\[
(\rho_\gamma' + \rho_{e^\pm}')\frac{dT_\gamma}{dt} = -4H\rho_\gamma - 3H(\rho_{e^\pm}+p_{e^\pm}) + \delta C_{e^\pm},
\]
\[
\rho_{\nu,\rm tot}'\frac{dT_\nu}{dt} = -4H\rho_{\nu,\rm tot} + \delta C_\nu,
\]
together with total energy conservation,
\[
\frac{d\rho_{\rm tot}}{dt} + 3H(\rho_{\rm tot}+p_{\rm tot})=0.
\]
From the solved thermal history, PRyMordial computes
\[
N_{\rm eff} \equiv \frac{8}{7}\left(\frac{11}{4}\right)^{4/3} \left(\frac{\rho_{\rm rad}-\rho_\gamma}{\rho_\gamma}\right),
\]
recovering the Standard-Model prediction \(N_{\rm eff}^{\rm SM}=3.044\) [2307.07061].

For the weak sector, the package includes the six standard neutron–proton interconversion processes and implements them beyond the Born approximation, with QED radiative corrections, finite nucleon-mass effects, weak magnetism, and finite-temperature corrections following the PRIMAT treatment [2307.07061]. In the Born limit, the rates are written as
\[
\Gamma_{n\to p}^{\infty} = \widetilde G_F^2 \int_0^\infty dE_e\, E_e\sqrt{E_e^2-m_e^2}\,(E_\nu^-)^2 \left[ f_\nu(E_\nu^-)f_e(-E_e) + f_\nu(-E_\nu^-)f_e(E_e) \right],
\]
\[
\Gamma_{p\to n}^{\infty} = \widetilde G_F^2 \int_0^\infty dE_e\, E_e\sqrt{E_e^2-m_e^2}\,(E_\nu^+)^2 \left[ f_\nu(E_\nu^+)f_e(-E_e) + f_\nu(-E_\nu^+)f_e(E_e) \right],
\]
with \(E_\nu^\pm = E_e \pm \mathcal Q\), \(\mathcal Q=m_n-m_p\), and either neutron-lifetime or direct weak-parameter normalization [2307.07061].

The nuclear sector evolves abundance yields \(Y_i=n_i/n_B\) through a coupled reaction network. PRyMordial solves nuclei up to \(A=8\), \(Z=5\), with two practical network modes: a reduced 12-reaction network sufficient for accurate \(Y_P\) and D/H, and a 63-reaction network used when improved lithium predictions are needed [2307.07061]. The general abundance evolution equation is
\[
\frac{dY_i}{dt} = \sum_R \mathcal S_{i,R} \left[ \Gamma^{(R)}_{\dots\to i\dots} \prod_j \left( \frac{Y_j^{\mathcal S_{j,R}}}{\mathcal S_{j,R}!} \right) - \Gamma^{(R)}_{i\dots\to\dots} \prod_k \left( \frac{Y_k^{\mathcal S_{k,R}}}{\mathcal S_{k,R}!} \right) \right].
\]
Observable definitions follow standard BBN conventions,
\[
Y_P \equiv 4Y_{^4\mathrm{He}} \simeq \rho_{^4\mathrm{He}}/\rho_B, \qquad i/\mathrm H \equiv \frac{Y_i}{Y_p}, \quad i=\mathrm D,\ ^3\mathrm{He},\ ^7\mathrm{Li}.
\]
The package therefore implements the full chain from plasma thermodynamics to weak freeze-out to final abundance synthesis within a single stiff ODE framework [2307.07061].

## 3. Software architecture, interfaces, and execution model

PRyMordial is written in Python 3 and is publicly available on GitHub, with optional acceleration through Julia and the SciML ecosystem via `diffeqpy` and `PyJulia` [2307.07061]. Its code structure is explicitly modular. `PRyM_init.py` contains constants, baseline inputs, boolean flags, and rate tables; `PRyM_thermo.py` handles thermal quantities and collision terms; `PRyM_nTOp.py` and `PRyM_evalnTOp.py` govern the weak-rate layer; `PRyM_nuclear_net12.py` and `PRyM_nuclear_net63.py` implement the reduced and full nuclear networks; and `PRyM_main.py` provides the user-facing interface [2307.07061].

A minimal run is
```python
import PRyM.PRyM_main as PRyMmain
res = PRyMmain.PRyMclass().PRyMresults()
```
which returns an array containing \(N_{\rm eff}\), relic neutrino quantities, \(Y_P\), \(D/H\), \({}^3\mathrm{He}/H\), and \({}^7\mathrm{Li}/H\) [2307.07061]. The code also supports persistence of precomputed thermal backgrounds and weak-rate tables, which is important for repeated likelihood evaluations.

User-facing flags govern both physics and performance. Standard options include `smallnet_flag` for the 12-reaction network, `nacreii_flag` for NACRE II-style key rates, `tau_n_flag` for weak-rate normalization via neutron lifetime, and `julia_flag` for Julia-based ODE solving [2307.07061]. For BSM work, the code exposes switches such as `NP_thermo_flag`, `NP_nu_flag`, `NP_e_flag`, `NP_nTOp_flag`, and `NP_nuclear_flag`, together with explicit new-physics parameters such as `NP_delta_nTOp` and reaction-by-reaction `NP_delta_R` shifts [2307.07061].

Later analyses show that this modularity is operational rather than merely nominal. The sensitivity atlas study computed response coefficients for 14 fundamental or cosmological parameters and 63 thermonuclear rates using PRyMordial under two weak-rate normalization schemes and two nuclear-rate compilations, turning the code into a controlled response-function generator for BBN observables [2603.22414]. This suggests that PRyMordial’s architecture is well suited both to forward prediction and to local Jacobian-based uncertainty and sensitivity studies.

## 4. Uncertainty propagation, rate libraries, and sensitivity structure

A distinctive feature of PRyMordial is its explicit treatment of nuclear-rate uncertainties. Forward thermonuclear rates are modeled as log-normal distributed,
\[
\log \widetilde\Gamma^{(R)} = \log \langle \widetilde\Gamma^{(R)} \rangle + p^{(R)} \log \Delta \widetilde\Gamma^{(R)},
\]
with temperature-independent Gaussian nuisance parameters \(p^{(R)}\) [2307.07061]. This provides a natural interface for Monte Carlo propagation of nuclear systematics. The 2024 baryon-abundance update made this capability central: PRyMordial was used to marginalize explicitly over reaction-rate uncertainties with log-normal priors, yielding a conservative BBN baryon abundance
\[
\Omega_b h^2 = 0.02218 \pm 0.00055
\]
for PDG light-element abundances in \(\Lambda\)CDM, and
\[
\Omega_b h^2 = 0.02196 \pm 0.00063
\]
in \(\Lambda\)CDM+\(N_{\rm eff}\), with light-element-only constraints on extra radiation at the \(\sim 0.2\) level [2401.15054].

The dominant nuclear issue in precision deuterium prediction is the treatment of deuterium destruction, especially
\[
{}^2\mathrm{H}+p\to{}^3\mathrm{He}+\gamma,\qquad
{}^2\mathrm{H}+{}^2\mathrm{H}\to{}^3\mathrm{He}+n,\qquad
{}^2\mathrm{H}+{}^2\mathrm{H}\to{}^3\mathrm{H}+p.
\]
PRyMordial supports two practically important key-rate choices: a NACRE II-based mode and a PRIMAT-based mode [2307.07061]. Later studies showed that this choice is not a minor implementation detail. In neutron-lifetime scans, helium-4 remained a robust probe of \(\tau_n\), whereas deuterium conclusions changed substantially between NACRE II and PRIMAT choices [2210.12031]. In baryon-abundance inference, experimentally driven rates and ab-initio rates differed at roughly the \(1\sigma\) level in inferred \(\Omega_b h^2\), with NACRE II marginalization preferred as a deliberately conservative synthesis [2401.15054].

The 2026 sensitivity atlas sharpened this picture. With PRIMAT-like rates, D/H is dominated by \(\Omega_b h^2\) and the three deuterium-destruction reactions above; with NACRE-II, the \(ddn\) and \(ddp\) reactions dominate the D/H uncertainty budget even more strongly [2603.22414]. By contrast, \(Y_p\) is largely insensitive to nuclear rates and is controlled mainly by weak freeze-out inputs such as \(Q\), \(\tau_n\), \(g_A\), \(V_{ud}\), \(m_e\), and, when allowed to vary, \(\Delta N_{\rm eff}\) [2603.22414]. The same study also showed that if \(\Delta N_{\rm eff}\) is treated as free with current external uncertainty, it dominates the theoretical uncertainty in \(Y_p\), overwhelming most purely nuclear contributions [2603.22414].

A common misconception is that BBN theory errors are set only by observational abundance precision or by the baryon density. PRyMordial-based analyses instead indicate a more heterogeneous structure: helium-4 is weak-sector and expansion-rate dominated, deuterium is compilation- and rate-library-sensitive, and lithium-7 is controlled mainly by the \({}^7\mathrm{Be}\)-production and destruction network [2401.15054].

## 5. Extension mechanisms for beyond-standard BBN

PRyMordial was explicitly designed for nonstandard BBN, and later work demonstrates several distinct extension patterns. The simplest is background modification through the new-physics thermodynamic interface. In the Weylian-boundary study, the authors solved an external modified-gravity background, converted it into temperature-dependent functions
\[
\rho_{\rm NP}(T),\qquad p_{\rm NP}(T),\qquad \frac{d\rho_{\rm NP}}{dT},
\]
and passed these to the `PRyMclass` constructor via `rho_NP`, `p_NP`, and `drho_NP_dT`, leaving the internal nuclear network unchanged [2509.01162]. The same structural strategy was used in the noncommutative-spacetime study, where a deformed photon gas modified the radiation equation of state and therefore the BBN expansion rate, with PRyMordial acting as the downstream abundance solver [2510.10685]. The Early Dark Energy analysis likewise embedded model-dependent \(\rho_{DE}(T)\) or \(\rho_{DE}(a)\) into a PRyMordial-based pipeline through a wrapper code and nested-sampling inference [2605.26749].

A second pattern is direct modification of internal microphysics. The varying-weak-scale study altered the Higgs-vev-dependent weak and hadronic inputs, including
\[
\frac{\delta v}{v} \equiv \frac{v_{\rm BBN}-v_0}{v_0},\qquad
m_f \to m_f \left(1+\frac{\delta v}{v}\right),\qquad
G_F \propto v^{-2},
\]
together with a modified neutron–proton mass difference
\[
\frac{m_n-m_p}{\rm MeV}=2.493\left(1+\frac{\delta v}{v}\right)-1.2,
\]
and a \(v\)-dependent \(\Lambda_{\rm QCD}\) and deuteron binding energy, all of which were propagated through the full PRyMordial evolution rather than treated as post-processing shifts [2402.08626].

A third pattern is correlated parameter-response embedding. The unification-scenario extension introduced a perturbative module in which selected Standard-Model and cosmological inputs were updated according to
\[
\frac{\Delta Q}{Q}=f_Q(R,S)\frac{\Delta\alpha}{\alpha},\qquad
Q=Q_0\left(1+f_Q\frac{\Delta\alpha}{\alpha}\right),
\]
with explicit response coefficients for \(m_e\), \(m_p\), \(Q_N\), \(\tau_n\), \(G_F\), \(m_W\), \(m_Z\), and \(G_N\), under two distinct gravitational-sector assumptions [2604.04870]. That work modified selected variables in `PRyM_init.py` and used PRyMordial as a self-consistent propagator from correlated coupling variations to \(Y_p\) and D/H.

The most invasive example is the temperature-dependent CPT-violation analysis, which used a modified PRyMordial branch to treat unequal electron and positron masses controlled by
\[
b_0(T)=\alpha T^2.
\]
There the code solved a temperature-dependent electron chemical potential from charge neutrality,
\[
n_{e^-}(m_e,\mu,T)=n_{e^+}(m_{e^+},-\mu,T),
\]
and recomputed separate \(e^\pm\) thermodynamics, weak \(n\leftrightarrow p\) rates, neutrino decoupling collision terms, and QED plasma corrections [2601.06259]. This shows that PRyMordial can be extended beyond background-only deformations to genuinely asymmetric finite-temperature plasma physics.

Taken together, these extensions establish PRyMordial less as a single fixed BBN calculation and more as a platform whose primary abstraction layers are the thermal background, weak kernels, and nuclear network. A plausible implication is that its long-term value lies in this modular decomposition rather than in any single baseline Standard-BBN configuration.

## 6. Representative applications and scientific impact

PRyMordial has been used in both narrowly targeted Standard-BBN analyses and broader BSM studies. Representative examples are summarized below [2210.12031].

| Study | PRyMordial role | Main inference |
|---|---|---|
| Neutron lifetime anomaly [2210.12031] | Scan \(840\text{ s} \le \tau_n \le 1050\text{ s}\) at fixed \(N_{\rm eff}=3.044\) | \(Y_p\) is strongly sensitive to \(\tau_n\); D/H is only mildly sensitive and compilation-dependent |
| 2024 baryon update [2401.15054] | Explicit marginalization over nuclear-rate uncertainties | Conservative BBN baryon abundance \(\Omega_b h^2 = 0.02218 \pm 0.00055\) in \(\Lambda\)CDM |
| Weak-scale variation [2402.08626] | Modified weak and hadronic inputs with \(v\)-dependent masses and \(G_F\) | \(|\delta v/v|\sim 0.01\) bounds; EMPRESS-like low helium can be improved only at the cost of worsening D/H |
| Weylian boundary / noncommutativity / EDE [2509.01162; 2510.10685; 2605.26749] | New-physics background injection via \(\rho_{\rm NP}(T)\), \(p_{\rm NP}(T)\), \(d\rho_{\rm NP}/dT\) | BBN strongly limits nonstandard expansion histories and deformed radiation sectors |
| Sensitivity atlas [2603.22414] | Uniform response study over 14 parameters and 63 rates | Provides a model-independent ranking of dominant inputs for \(Y_p\), D/H, \({}^7\mathrm{Li}/H\), and \(N_{\rm eff}\) |
| Correlated varying couplings [2604.04870] | Perturbative response module for \(\Delta\alpha/\alpha\), \(R\), \(S\) | Constrains \(\Delta\alpha/\alpha\) at tens-of-ppm level during BBN |

Two broad scientific lessons recur across these applications. First, PRyMordial repeatedly identifies helium-4 as the cleaner direct probe of weak freeze-out and expansion-rate physics. In the neutron-lifetime study, \(Y_p\) favored the bottle value and slightly disfavored the beam value, whereas D/H could not robustly distinguish them because of reaction-rate systematics [2210.12031]. In the sensitivity atlas, \(Y_p\) remained dominated by weak-sector and \(N_{\rm eff}\)-type inputs rather than by nuclear rates [2603.22414].

Second, the code has sharpened rather than eliminated controversies internal to BBN. The most important is the rate-compilation dependence of deuterium. PRyMordial-based work does not remove the difference between empirically fitted and ab-initio deuterium-burning treatments; instead it quantifies how that difference propagates into \(\Omega_b h^2\), D/H consistency, and the apparent robustness of BBN–CMB agreement [2401.15054]. Likewise, PRyMordial extensions aimed at the lithium problem have not produced a full solution. Varying-coupling models allowed by \(Y_p\) and D/H were found to move lithium in the right direction only weakly, far short of the factor-of-three reduction required [2604.04870].

In this sense, PRyMordial’s significance is twofold. It is a precision BBN engine capable of reproducing state-of-the-art Standard-Model observables, and it is a modular inference tool for MeV-era cosmology whose outputs are sufficiently differential to expose where present theoretical limitations lie: weak normalization conventions, deuterium-burning systematics, and the incomplete closure of the lithium sector.

Source: https://www.emergentmind.com/topics/prymordial