PRyMordial: Precision BBN Calculations

• PRyMordial is a public numerical package for precision Big Bang Nucleosynthesis that addresses both Standard Model and beyond-standard-model scenarios.
• It integrates first-principles radiation thermodynamics, detailed weak n↔p conversion, and a modular nuclear reaction network for robust cosmological inference.
• Its flexible design supports extensions to modified thermal histories, varying fundamental constants, and nonstandard plasma physics for comprehensive BBN studies.

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 (Burns et al., 2023). 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 (Chowdhury et al., 2022).

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 (Burns et al., 2023). 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 (Burns et al., 2023). 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 (Burns et al., 2023).

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 (Burns et al., 2023). The baseline expansion is determined by the Friedmann equation

$\mathcal{O}(10)\,\mathrm{MeV}$0

with $\mathcal{O}(10)\,\mathrm{MeV}$1 built from photons, $\mathcal{O}(10)\,\mathrm{MeV}$2, neutrinos, and optionally new sectors. The code distinguishes the photon/plasma temperature $\mathcal{O}(10)\,\mathrm{MeV}$3 from the neutrino temperature $\mathcal{O}(10)\,\mathrm{MeV}$4, and evolves them through coupled Boltzmann equations rather than assuming instantaneous decoupling.

The thermal sector is governed by

$\mathcal{O}(10)\,\mathrm{MeV}$5

$\mathcal{O}(10)\,\mathrm{MeV}$6

together with total energy conservation,

$\mathcal{O}(10)\,\mathrm{MeV}$7

From the solved thermal history, PRyMordial computes

$\mathcal{O}(10)\,\mathrm{MeV}$8

recovering the Standard-Model prediction $\mathcal{O}(10)\,\mathrm{MeV}$9 (Burns et al., 2023).

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 (Burns et al., 2023). In the Born limit, the rates are written as

$\mathcal{O}(\mathrm{keV})$0

$\mathcal{O}(\mathrm{keV})$1

with $\mathcal{O}(\mathrm{keV})$2, $\mathcal{O}(\mathrm{keV})$3, and either neutron-lifetime or direct weak-parameter normalization (Burns et al., 2023).

The nuclear sector evolves abundance yields $\mathcal{O}(\mathrm{keV})$4 through a coupled reaction network. PRyMordial solves nuclei up to $\mathcal{O}(\mathrm{keV})$5, $\mathcal{O}(\mathrm{keV})$6, with two practical network modes: a reduced 12-reaction network sufficient for accurate $\mathcal{O}(\mathrm{keV})$7 and D/H, and a 63-reaction network used when improved lithium predictions are needed (Burns et al., 2023). The general abundance evolution equation is

$\mathcal{O}(\mathrm{keV})$8

Observable definitions follow standard BBN conventions,

$\mathcal{O}(\mathrm{keV})$9

The package therefore implements the full chain from plasma thermodynamics to weak freeze-out to final abundance synthesis within a single stiff ODE framework (Burns et al., 2023).

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 (Burns et al., 2023). 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 (Burns et al., 2023).

A minimal run is $\mathcal{O}(10)\,\mathrm{MeV}$05 which returns an array containing $n\leftrightarrow p$0, relic neutrino quantities, $n\leftrightarrow p$1, $n\leftrightarrow p$2, $n\leftrightarrow p$3, and $n\leftrightarrow p$4 (Burns et al., 2023). 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 (Burns et al., 2023). 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 (Burns et al., 2023).

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 (Burns, 23 Mar 2026). 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,

$n\leftrightarrow p$5

with temperature-independent Gaussian nuisance parameters $n\leftrightarrow p$6 (Burns et al., 2023). 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

$n\leftrightarrow p$7

for PDG light-element abundances in $n\leftrightarrow p$8CDM, and

$n\leftrightarrow p$9

in $Y_P$0CDM+$Y_P$1, with light-element-only constraints on extra radiation at the $Y_P$2 level (Schöneberg, 2024).

The dominant nuclear issue in precision deuterium prediction is the treatment of deuterium destruction, especially

$Y_P$3

PRyMordial supports two practically important key-rate choices: a NACRE II-based mode and a PRIMAT-based mode (Burns et al., 2023). Later studies showed that this choice is not a minor implementation detail. In neutron-lifetime scans, helium-4 remained a robust probe of $Y_P$4, whereas deuterium conclusions changed substantially between NACRE II and PRIMAT choices (Chowdhury et al., 2022). In baryon-abundance inference, experimentally driven rates and ab-initio rates differed at roughly the $Y_P$5 level in inferred $Y_P$6, with NACRE II marginalization preferred as a deliberately conservative synthesis (Schöneberg, 2024).

The 2026 sensitivity atlas sharpened this picture. With PRIMAT-like rates, D/H is dominated by $Y_P$7 and the three deuterium-destruction reactions above; with NACRE-II, the $Y_P$8 and $Y_P$9 reactions dominate the D/H uncertainty budget even more strongly (Burns, 23 Mar 2026). By contrast, $D/H$0 is largely insensitive to nuclear rates and is controlled mainly by weak freeze-out inputs such as $D/H$1, $D/H$2, $D/H$3, $D/H$4, $D/H$5, and, when allowed to vary, $D/H$6 (Burns, 23 Mar 2026). The same study also showed that if $D/H$7 is treated as free with current external uncertainty, it dominates the theoretical uncertainty in $D/H$8, overwhelming most purely nuclear contributions (Burns, 23 Mar 2026).

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 $D/H$9-production and destruction network (Schöneberg, 2024).

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

${}^3\mathrm{He}/H$0

and passed these to the PRyMclass constructor via rho_NP, p_NP, and drho_NP_dT, leaving the internal nuclear network unchanged (Matei et al., 1 Sep 2025). 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 (Matei et al., 12 Oct 2025). The Early Dark Energy analysis likewise embedded model-dependent ${}^3\mathrm{He}/H$1 or ${}^3\mathrm{He}/H$2 into a PRyMordial-based pipeline through a wrapper code and nested-sampling inference (Matei et al., 26 May 2026).

A second pattern is direct modification of internal microphysics. The varying-weak-scale study altered the Higgs-vev-dependent weak and hadronic inputs, including

${}^3\mathrm{He}/H$3

together with a modified neutron–proton mass difference

${}^3\mathrm{He}/H$4

and a ${}^3\mathrm{He}/H$5-dependent ${}^3\mathrm{He}/H$6 and deuteron binding energy, all of which were propagated through the full PRyMordial evolution rather than treated as post-processing shifts (Burns et al., 2024).

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

${}^3\mathrm{He}/H$7

with explicit response coefficients for ${}^3\mathrm{He}/H$8, ${}^3\mathrm{He}/H$9, ${}^7\mathrm{Li}/H$0, ${}^7\mathrm{Li}/H$1, ${}^7\mathrm{Li}/H$2, ${}^7\mathrm{Li}/H$3, ${}^7\mathrm{Li}/H$4, and ${}^7\mathrm{Li}/H$5, under two distinct gravitational-sector assumptions (Dreyer et al., 6 Apr 2026). That work modified selected variables in PRyM_init.py and used PRyMordial as a self-consistent propagator from correlated coupling variations to ${}^7\mathrm{Li}/H$6 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

${}^7\mathrm{Li}/H$7

There the code solved a temperature-dependent electron chemical potential from charge neutrality,

${}^7\mathrm{Li}/H$8

and recomputed separate ${}^7\mathrm{Li}/H$9 thermodynamics, weak $N_{\rm eff}$0 rates, neutrino decoupling collision terms, and QED plasma corrections (Barenboim et al., 9 Jan 2026). 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 (Chowdhury et al., 2022).

Study	PRyMordial role	Main inference
Neutron lifetime anomaly (Chowdhury et al., 2022)	Scan $N_{\rm eff}$1 at fixed $N_{\rm eff}$2	$N_{\rm eff}$3 is strongly sensitive to $N_{\rm eff}$4; D/H is only mildly sensitive and compilation-dependent
2024 baryon update (Schöneberg, 2024)	Explicit marginalization over nuclear-rate uncertainties	Conservative BBN baryon abundance $N_{\rm eff}$5 in $N_{\rm eff}$6CDM
Weak-scale variation (Burns et al., 2024)	Modified weak and hadronic inputs with $N_{\rm eff}$7-dependent masses and $N_{\rm eff}$8	$N_{\rm eff}$9 bounds; EMPRESS-like low helium can be improved only at the cost of worsening D/H
Weylian boundary / noncommutativity / EDE (Matei et al., 1 Sep 2025, Matei et al., 12 Oct 2025, Matei et al., 26 May 2026)	New-physics background injection via $n\leftrightarrow p$0, $n\leftrightarrow p$1, $n\leftrightarrow p$2	BBN strongly limits nonstandard expansion histories and deformed radiation sectors
Sensitivity atlas (Burns, 23 Mar 2026)	Uniform response study over 14 parameters and 63 rates	Provides a model-independent ranking of dominant inputs for $n\leftrightarrow p$3, D/H, $n\leftrightarrow p$4, and $n\leftrightarrow p$5
Correlated varying couplings (Dreyer et al., 6 Apr 2026)	Perturbative response module for $n\leftrightarrow p$6, $n\leftrightarrow p$7, $n\leftrightarrow p$8	Constrains $n\leftrightarrow p$9 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, $\mathcal{O}(10)\,\mathrm{MeV}$00 favored the bottle value and slightly disfavored the beam value, whereas D/H could not robustly distinguish them because of reaction-rate systematics (Chowdhury et al., 2022). In the sensitivity atlas, $\mathcal{O}(10)\,\mathrm{MeV}$01 remained dominated by weak-sector and $\mathcal{O}(10)\,\mathrm{MeV}$02-type inputs rather than by nuclear rates (Burns, 23 Mar 2026).

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 $\mathcal{O}(10)\,\mathrm{MeV}$03, D/H consistency, and the apparent robustness of BBN–CMB agreement (Schöneberg, 2024). Likewise, PRyMordial extensions aimed at the lithium problem have not produced a full solution. Varying-coupling models allowed by $\mathcal{O}(10)\,\mathrm{MeV}$04 and D/H were found to move lithium in the right direction only weakly, far short of the factor-of-three reduction required (Dreyer et al., 6 Apr 2026).

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.