PRyMordial: the first three minutes, within and beyond the standard model

In this work we present PRyMordial: A package dedicated to efficient computations of observables in the Early Universe with the focus on the cosmological era of Big Bang Nucleosynthesis (BBN). The code offers fast and precise evaluation of BBN light-element abundances together with the effective number of relativistic degrees of freedom, including non-instantaneous decoupling effects. PRyMordial is suitable for state-of-the-art analyses in the Standard Model as well as for general investigations into New Physics active during BBN. After reviewing the physics implemented in PRyMordial, we provide a short guide on how to use the code for applications in the Standard Model and beyond. The package is written in Python, but more advanced users can optionally take advantage of the open-source community for Julia. PRyMordial is publicly available on GitHub.


Introduction
The snapshot of the Universe approximately three minutes after the Big Bang [1] can be regarded as one of the most remarkable predictions of the Standard Model (SM) of Particle Physics in conjunction with the (so-called) concordance model of Cosmology, ΛCDM.
While a theory for the origin of chemical elements based on an epoch of high-energy densities and pressures was already formulated by Alpher, Bethe, and Gamow more than seventy years ago [2], the discovery of the quasi-black body spectrum of the Cosmic Microwave Background (CMB) [3,4] paved the road for the modern formulation of the theory of Big Bang Nucleosynthesis (BBN) [5].Indeed, thanks to the CMB, we know today that the SM particle species were in a thermal state during an epoch dominated by radiation.Extrapolating this cosmological picture back in time when the Universe was not yet transparent to light, within the standard lore of Cosmology and of Particle Physics we can accurately predict [6][7][8][9][10][11]: 1) The evolution of the number of relativistic degrees of freedom until recombination, N eff ; 2) The cosmological abundance of light nuclides synthesized from protons and neutrons, as a function of the number density of baryons relative to photons, η B ≡ n B /n γ .
Regarding 1), given the current knowledge of neutrino oscillations [12], N eff is predicted in the SM via solving a set of integro-differential equations for the neutrino density matrix at finite temperature [13], yielding N SM eff = 3.044 with an error estimated to be below the level of per mil [14][15][16].Concerning 2), a detailed analysis of CMB anisotropies in temperature and polarization currently constrains η B with 1% accuracy or better [17], anchoring the primordial asymmetry between baryons and anti-baryons to be O(10 −10 ) [18].Assuming no large asymmetry in the lepton sector as well, see e.g.[19], standard BBN turns into an extremely predictive theory, often dubbed "parameter free".
On the observational side, multi-wavelength astronomical campaigns have been able to provide rich spectroscopic information about emission and absorption lines of gas clouds in metal-poor extragalactic environments, see e.g.[20][21][22][23], bringing us today to a percent-level determination of the abundance of primordial deuterium and helium-4.Given the predictions of the standard theory and the precision of those measurements, together with the strong constraints on the thermal history provided by the CMB [24,25], the study of the Early Universe around the BBN epoch offers unique insight on New Physics (NP) [26][27][28][29][30][31][32][33][34][35].Looking at the exciting prospects of next-gen CMB experiments [36][37][38], and at the expected future sensitivity in the field of observational astronomy [39,40], it is therefore very timely to have tools at our disposal that allow for numerically efficient, yet precise computations that test the SM in the Early Universe, and that are flexible enough to broadly explore NP scenarios.
A few packages have already been developed to accurately investigate the BBN era.A publicly available version of the historical code of Ref. [41] (whose most up-to-date version is currently adopted by the PDG [42]) is described in [43].At the same time, publicly released codes dedicated to stateof-the-art BBN analyses are also available; in particular: • PArthENoPE [44][45][46] is a code originally written in FORTRAN 77 that in its latest re-incarnation also enjoys a graphical user interface; it offers a very efficient evaluation of BBN light-element abundances based on fitting formulae worked out for both weak rates and nuclear cross sections.
• PRIMAT [47,48] is an user-friendly Mathematica package containing all the inputs and ingredients for an ab-initio computation of neutron freeze-out and of weak rates; moreover, it has tabulated the largest nuclear network at hand in order to track the abundance of heavy nuclides as well.
Both codes include a few built-in options to account for the study of some specific NP scenarios.
AlterBBN [49,50] is a C++ open-source software developed for broad investigation of Physics Beyond the SM (BSM) in the BBN era.However, while allowing for fast numerical evaluations, AlterBBN does not implement the same level of detail and accuracy in its computation of light primordial abundances compared to PArthENoPE or PRIMAT.In fact, these two packages may currently represent the best tools to perform precision cosmological analyses [24,51].While powerful and flexible, these public codes nevertheless suffer from a few limitations and/or missing features.A precision tool for Cosmology, able to handle BSM Particle Physics should: -Allow for the evaluation of the physics of the thermal bath in a fast but precise way, following, e.g., the approach highlighted in [34,52,53], and implemented in the standalone code NUDEC BSM; -Interconnect a first-principle computation of the thermal background with an ab-initio precise calculation of the neutron-to-proton (n ↔ p) conversion, as the one implemented in PRIMAT [47]; -Render easily accessible exploration of the impact of the input parameters characterizing the BBN era and the uncertainties in the set of thermonuclear rates on the basis of more modeldependent/more data-driven approaches available in literature, see [35,54,55]; -Adopt a user-friendly, modern programming language compatible with numerical efficiency of the computations, while smoothly interfacing with standard libraries for statistically advanced analyses like Monte Carlo (MC) ones, see e.g.[56,57].
In this work, we introduce PRyMordial: A new public tool for the community of Particle Physics and Cosmology that precisely aims at filling in the above gaps for precision studies on the physics of the Early Universe both within and beyond the SM.The package is written and runs entirely with Python 3.Moreover, for the most advanced users, the resolution of the set of stiff differential equations for the BBN nuclear-reaction network can be further optimized with the optional switch to some routines of the SciML kit [58], the open-source software for scientific machine learning in Julia.
Our article is organized as follows: In section 2 we present all the key ingredients of the physics implemented in PRyMordial; In section 3 we discuss in detail how PRyMordial is structured and we provide several examples on the usage of the code; In section 4 we comment on future directions for further development of PRyMordial along with possible interesting applications.We finally collect in Appendix A a set of instructions for the installation of the package and its dependencies.

Physics in PRyMordial
In this section we present the key equations present in PRyMordial, which stand out as a reference for the physics implemented within the code as well as representing a guideline regarding its use (see section 3).We organize the presentation in three distinct topics: the thermodynamics of the plasma; the weak rates for n ↔ p conversion; and the set of thermonuclear rates for the key reactions responsible of the non-zero primordial abundance of deuterium, helium-3 and -4, and lithium-7.

Thermodynamics with Non-instantaneous Decoupling
The description of the thermal background during the BBN era in ΛCDM follows from an isotropic, homogeneous Universe modelled by the Einstein field equation: where H is the Hubble rate of space-time expansion, a the scale factor of the Friedmann-Lemaître-Robertson-Walker metric, ρ tot the total energy density present in the Universe, and Within an axiomatic characterization of the Early Universe provided by local thermodynamic equilibrium [59,60], SM species are described according to the spin-statistics theorem and the temperature T γ of the thermal bath (provided chemical potentials µ can be neglected, i.e., µ/T γ ≪ 1).Standard BBN takes place during radiation domination, and thus features contributions to ρ tot largely from relativistic species, i.e. ρ tot ≃ ρ rad ∝ T 4 γ .This observation dramatically simplifies the investigation of BBN, allowing one to decouple the study of the thermal background from the nucleon dynamics.Indeed, after the QCD crossover takes place [61] protons and neutrons are already non-relativistic, i.e. they are highly Boltzmann-suppressed well before the MeV scale temperatures characteristic of the BBN era.
Hence, for temperatures T γ ≲ O(10) MeV, one can accurately describe ρ tot in the SM as a sum of just three contributions: where x e ≡ m e /T γ and we distinguish the temperature of the electron-positron-photon system, T γ , from that of neutrinos, T ν . 1 Indeed, while the initial condition T ν = T γ must hold at early times for the two systems to be in thermal (more precisely, in chemical and kinetic) equilibrium, around the MeV scale neutrinos are expected to freeze out from the thermal bath as weakly-interacting relativistic species [63].Neglecting tiny departures from a Fermi-Dirac distribution in ν phase space, one can study the evolution of the two systems according to the momentum-integrated Boltzmann equations: with ′ ≡ d/dT , p the pressure density (equal to ρ/3 for a relativistic species), δC the (momentumintegrated) collision term, and where we have conveniently traded energy densities for temperatures in light of Eq. (2.2).Due to energy-momentum conservation, the sum over all δCs must vanish, so that one recovers the continuity equation for the total energy density of the Universe: In the SM, where Eq. (2.3) holds, such a constraint implies: δC ν = −δC e ± .The collision term δC ν has been evaluated in [63] under Maxwell-Boltzmann approximation, nicely refined in [52,53] taking into account relativistic corrections as well as finite mass effects from m e ̸ = 0, and more recently recomputed independently in [34].Including finite temperature QED corrections to the electromagnetic plasma [64], one can solve the system of coupled differential equations in Eq. (2.3), to find T γ (t), T ν (t), and, as a byproduct, T ν (T γ ). 2 Such a treatment naturally includes non-instantaneous decoupling effects, and allows one to perform a numerically fast, but accurate prediction of the effective number of relativistic degrees of freedom from first principles, yielding (in the SM) at T γ ≪ MeV: while also opening up novel explorations of BSM physics in the Early Universe [33,34,53]. 3 Based on these results, one can also easily evaluate the relic density of neutrinos (neglecting phase space spectral distortions).From the CMB we know the photon temperature today is T γ,0 = 0.2348 meV; plugging this value into the solution of Eq. (2.3) yields the temperature T ν,0 = 0.1682 meV, corresponding to the cosmological abundance of SM neutrinos: ) 1 While Te = Tγ follows from e ± being tightly coupled to photons via fast QED processes, the approximation underlying Tν , namely Tν e ≃ Tν µ ≃ Tν τ , can be motivated by the active flavor mixing of ν oscillations at Tγ of few MeV [62]. 2 In the current version of PRyMordial we adopt the computation of δCν as well as the next-to-leading (NLO) QED corrections to the electromagnetic pressure of the plasma directly from the numerical results tabulated in NUDEC BSM [53]. 3 Eq.(2.3) can be easily generalized to include new sectors.This contrasts with typical existing BBN codes which compute the thermodynamic background by interpolating the tabulated result of the (numerically intensive) integrodifferential Boltzmann equation, solved for the neutrino phase-space density in the SM.
which reproduces the relic neutrino abundance computed, e.g., in Ref. [65] to the per mil level.
In order to obtain T γ (t) and T ν (t) from Eq. (2.3), we have made use both of Eq. (2.1) together with Eq. (2.2).At this point, to complete the study of the thermodynamic background, we must extract the scale factor a as a function of time t and temperature T γ .This can be accomplished by applying (again) the notion of local thermodynamic equilibrium, which allows one to introduce the entropy density for each species i as: s i = (ρ i + p i − µ i n i )/T i , where n i is the number density of the species with associated chemical potential µ i .
For negligible chemical potentials, the total entropy density of the Universe s tot per comoving volume must be conserved as a consequence of energy-momentum conservation, Eq. (2.4).Then, during radiation domination s tot roughly scales as T 3 γ , underlying the approximate relation a ∝ 1/T γ .Nevertheless, even under the assumption of µ i /T i ≪ 1, the entropy of each species is generally not separately conserved due to heat exchanges related to the interactions with other species.The Boltzmann equation for s i generally follows (see, e.g., the discussion in Refs.[47,66]): where the first collision term (divided by the temperature) is the one appearing in the Boltzmann equation for the density ρ i , while the second collision term has been rewritten using the Boltzmann equation for the number density n i . 4In the SM, in the limit5 µ e /T γ ≪ 1, we use Eq.(2.7) for the electromagnetic bath to pin down the relation between a and T γ ; with spl ≡ (s γ + s e ± )/T 3 γ , we get: Knowing all the thermodynamical quantities as a function of T γ in the integrand above, Eq. (2.8) allows one to extract a(T γ ) up to the scale-factor value of today, a 0 , customarily defined as 1.Note that for T γ ≲ m e one has s ′ pl = 0, and taking the limit N ν → 0, the expected scaling set by d(s γ a 3 )/dt = 0 is easily recovered.The solution in Eq. (2.8) precisely tracks the relation between the scale factor and T γ in the case of non-instantaneous decoupling of neutrinos.While in the SM these effects are tiny (since N ν /3 ≪ spl ), they could become non-negligible in a BSM scenario.
It is worth noting that given T γ (t) from the solution of Eq. (2.3) and a(T γ ) from Eq. (2.8), one obtains a(t) as a byproduct, which allows to assess the evolution of the number density of baryons in t or T γ during the BBN era, since by definition: n B ∝ 1/a 3 .

Neutron Freeze Out beyond the Born Approximation
Shortly after hadrons form, neutrons and protons are non-relativistic species that do not contribute appreciably to the total energy budget stored in the thermal bath.Nevertheless, their abundance is eventually responsible for the tiny fraction of light primordial elements relative to hydrogen which are observable today in pristine astrophysical environments.
According to local thermodynamic equilibrium, the relative number density of nucleons is initially given by the Maxwell-Boltzmann distribution: where Q = m n − m p , µ Q = µ n − µ p , m n,p and µ n,p are the mass and chemical potential of neutrons and protons.For clarity, we have used T ν = T γ (valid for temperatures well above MeV) in the Q term, but retain T ν explicitly in the µ Q term.Assuming µ n ≃ µ p (e.g. a negligible contribution from lepton chemical potentials), Eq. (2.9) implies that at equilibrium n n ≃ n p .Indeed, fast electroweak processes efficiently convert n ↔ p: and govern the Boltzmann equations for the nucleon yields Y n,p ≡ n n,p / n B = n n,p / (n n + n p ): ) whose equilibrium solutions: , are in agreement with Eq. (2.9).These reactions guarantee chemical equilibrium among the involved species, implying µ Q ≃ −µ ν .Eq. (2.9) thus demonstrates that a primordial non-zero lepton asymmetry in the neutrino sector [67,68] can impact the initial conditions for BBN by altering the neutron-to-proton ratio, with notable cosmological consequences [35,69].At temperatures close to neutrino decoupling, n ↔ p conversion falls out of equilibrium, freezing out the neutron-to-proton ratio to ∼ 1/6 (in the SM), up to finite neutron lifetime effects [59,60].The weak rates for neutron freeze out require the evaluation of an involved multi-dimensional phase-space integral: e.g. for n e + → p ν (and similarly for the others) [70]: where dΠ i and P i are the Lorentz-invariant phase-space element and 4-momentum of the particle i, f i is the relativistic thermal distribution of the species i in the rest frame of the thermal bath, and M is the full matrix element of the process summed over initial and final spins.The latter can be computed from the weak effective theory for β decay [71]: where G F is the Fermi constant [42], V ud corresponds to the Cabibbo angle [72], g A and κ are the axialcurrent and weak-magnetism constant of the nucleon of mass m N [73], and The computation of |M| 2 can be found in detail in Appendix B of Ref. [47] (see also [70,74]).While expressions like Eq. (2.11) can be reduced to a five-dimensional integral in phase space by exploiting the symmetries of the problem, a dramatic simplification is obtained in the limit of infinite nucleon-mass at fixed Q [70,74].This is the so-called Born approximation, in which the kinetic energy of the 'infinitely' heavy neutrons and protons may be neglected, leading to the simplification: In that limit the n ↔ p rates read: where The outcome of Eq. (2.13) are rates that generally depend on both background temperatures and chemical potentials (i.e.T γ , T ν and µ ν ).For T ν = T γ (and negligible chemical potentials) detailed balance follows as: The dimensionful factor G F depends on V ud , g A , and G F , whose value is precisely determined by the muon lifetime.However, this factor is often more conveniently extracted from neutron decay in the vacuum, since in the SM: where F n incorporates a phase-space statistical factor for the neutron decay at zero temperature [75] plus electroweak radiative corrections [76].For a precise calculation of F n , see the very recent reassessment in Ref. [77] and references therein.This approach allows one to trade the combination V 2 ud (1 + 3g 2 A ) for the measured τ n . 6Using Eq. (2.14), in PRyMordial one can choose to adopt either a normalization of the weak rates based on the determination of the neutron lifetime, or one involving the knowledge of the modified Fermi constant G F .
In the SM the Born approximation predicts a neutron freeze-out temperature of slightly below 1 MeV.At smaller temperatures, the neutron-to-proton ratio is still affected by β decay until the Universe cools down sufficiently enough to preclude photo-dissociation of deuterium: for a binding energy B D = 2.2 MeV, this happens at temperatures around B D / log(1/η B ) ∼ 0.1 MeV [59,60].At that point, virtually all of the neutrons experience two-body nuclear reactions, ultimately resulting in their binding in helium-4, the most stable light element.As a result, the uncertainty on the Born-level theory prediction for helium-4 is only a few % (see Table 5 in [47]).
PRyMordial implements all of these corrections, following the treatment in PRIMAT (see Appendix B of [47]), where particular care was taken to attempt to combine several existing state-of-the-art recipes for electroweak rates beyond the Born approximation.
It is worth noticing that in the context of the SM, the corrections to the Born rates due to the incomplete neutrino decoupling are only marginal [87,88].Nevertheless, NP could dramatically alter T ν (T γ ), a(T γ ) and a(t), and the departure from the standard value for the weak rates can impact the final BBN abundances in a non-negligible way [31,33].As a result, the approach undertaken in subsection 2.1 is particularly useful not only for the study of neutrino decoupling, but also for a careful assessment of the neutron-to-proton ratio in BSM scenarios.

BBN Thermonuclear Reactions
Local thermodynamical equilibrium implies that at temperatures above neutron decoupling, a nuclear species i of atomic number Z i , mass number A i , spin s i , and binding energy B i follows a Boltzmann distribution with internal degrees of freedom: In terms of the yield Y i ≡ n i /n B , this equilibrium distribution reads: where we made use of the relation: ).This expression holds for the nucleons (A N = 1, B N = 0) themselves, and is consistent with Eq. (2.9).Importantly, it offers another handle on the estimate for the start of nucleosynthesis as the time in which the relative abundance of neutrons after freeze out becomes comparable to deuterium as dictated by Eq. (2.15), and pointing again to a temperature of about 0.1 MeV.
Starting from the initial conditions, abundances are determined by a network of Boltzmann equations that generalize Eq. (2.10) (see, e.g., Refs.[89,90]) to include the relevant nuclei: where the sum R is performed over all reactions involving the nuclear species i; S i,R is the stoichiometric coefficient for the species i in the nuclear reaction R; and the products j and k run over all of the initial and final states of the reaction with nuclear rate Γ (R) . Given the range of energies characterizing the BBN era, the nuclear reaction rates of interest can be measured in the laboratory, and are often tabulated as Γ i...l→j...m ≡ N S i ...S l −1 A ⟨σ i...l→j...m v⟩ [91], where N A is Avogadro's number (typically expressed in units of mol −1 ), and the velocity averaged cross section is obtained by weighting the appropriate cross section by the Maxwell-Boltzmann velocity distribution for the non-relativistic species (see e.g.Ref. [92] for a detailed description).By definition, for a given number-density rate ⟨σ (2.17) A priori, Eq. (2.16) includes the rates of both forward and reverse reactions in the evolution of the abundance of the nuclear species i. Nevertheless detailed balance implies since local thermodynamical equilibrium ensures that the forward and reverse reactions should balance.Thus, it is easy to evaluate the reverse reaction rates given the forward ones.It is customary to parameterize the relationship as: where the constants α R , β R , and γ R for a given process R from e.g. the up-to-date nuclear database of Ref. [93] via Eq.(2.15).
PRyMordial solves the general system of equations Eq. (2.16) following the strategy of Ref. [47] which conveniently breaks nucleosynthesis into three steps: We analyze n ↔ p conversion by solving Eq. (2.10) from an initial temperature of O( 10) MeV (and initial conditions from Eq. (2.9)) down to standard neutron freeze out, around MeV; 2) We use the values of Y n,p obtained from 1) together with Eq. (2.15) and evolve with a network comprised of the 18 key thermonuclear rates for the abundance of n, p together with all of the nuclides up to A = 8 and Z = 57 down to the temperature where deuterium photo-dissociation becomes inefficient, around 0.1 MeV; 3) We further evolve the network with the full set of thermonuclear processes and with initial conditions given by the nuclide yields obtained in step 2), evolving the abundances of the aforementioned nuclides down to O(keV) (i.e., well below e ± annihilation), when BBN is over.
The output of Step 3) is the abundances of the light-element originating from BBN.To compare with data, it is customary to quote helium-4 in terms of the primordial mass fraction:8 The other primordial elements under the lamppost of astrophysical observations are deuterium, helium-3 and lithium-7 (see, e.g., [11] for a recent report on the status of these measurements), which are usually quoted in terms of the relative number densities with respect to hydrogen: 3 He, 7 Li . (2.21) Notice that the final yield of primordial helium-3 receives a contribution from unstable species such as tritium; likewise, the final amount of lithium-7 includes the decay of beryllium-7.
The literature contains several publicly accessible compilations of the thermonuclear rates relevant for BBN.It is important to note that there are several different parameterizations of these rates adopted in BBN studies, and they differ not only with respect to the theoretical approach, but also with respect to the measured nuclear reaction data included in fitting them.To highlight a few of the more important approaches: • The NACRE II database [95] collects an extended evaluation of reaction rates of charged-particle induced reactions on target nuclides with mass number A < 16, adopting the so-called potential model [91] to describe nuclear cross sections in the energy range of interest.
• PArthENoPE implements semi-analytic expressions resulting from polynomial fits to nuclear data including theory modeling of screening and thermal effects [92,101]; data-oriented analyses relevant for BBN rates can be also found in the work of Refs.[102,103].
If one limits the scope to precise predictions of the helium-4 and deuterium abundances, the relevant portion of the nuclear network simplifies considerably, contracting to O(10) processes [104].Thus, PRyMordial offers the option of restricting the BBN analysis to a small network of 12 key reactions [105], implemented according to two different sets of thermonuclear rates: the first is largely based on the NACRE II compilation, whereas the second is based on the tabulated rates in PRIMAT.These two sets differ marginally in their predictions for helium-4, but lead to relevant differences in the prediction for deuterium, as discussed at length in Ref. [54], after the important measurement carried out by the LUNA collaboration [106]. 9For the most precise prediction of lithium-7, PRyMordial offers the possibility to solve a nuclear network including the 51 additional reactions listed in Appendix B, by adopting part of the network in Ref. [100] included in the PRIMAT database.
PRyMordial handles uncertainties on the tabulated thermonuclear rates Γ (R) by providing (for each forward10 nuclear reaction) a set of median values, ⟨ Γ (R) ⟩ together with an uncertainty factor ∆ Γ (R) , corresponding to a sample of temperatures.Following the method outlined in Refs.[107,108], to perform a MC analysis with PRyMordial one should treat the provided thermonuclear rates as lognormal distributed, implying that for each nuclear process R a random realization of the thermonuclear rate will be: log where p (R) is a temperature-independent coefficient following a normal distribution [109].Hence, in order to properly take into account the uncertainties of the thermonuclear rates in a MC analysis of BBN, one should independently vary the nuisance parameters p (R) for all the reactions R included in the study, see, e.g., the work carried out in Ref. [35] and the MC examples presented in section 3.

How to Use PRyMordial
In this section we provide some example code that demonstrates the use of PRyMordial.We start by detailing the modules of the code including their inputs and key parameters.We show how to implement a state-of-the-art analysis of the BBN era within the SM.Finally, we provide a concise description on how to use the code for the study of NP, and discuss how to implement and analyze generic BSM scenarios.
3.1 Structure of the Code and Hello, World!PRyMordial is a numerical tool dedicated to efficiently and accurately evaluate in the SM and beyond all the key observables related to the BBN era, discussed in section 2, namely: • The number of effective relativistic degrees of freedom, N eff , Eq. (2.5) ;

PRyM_init.py:
Initializes all the parameters and options for the study of the BBN era PRyM_main.py:Takes information from all other modules and computes all the key observables at the end of BBN  • The cosmic neutrino abundance today, Ω ν h 2 , Eq. (2.6) ; • The helium-4 mass fraction (both for BBN and CMB), Y P , Eq. (2.20) ; • The relative number density of deuterium, helium-3 and lithium-7, Eq. (2.21) .
In contrast to other BBN codes available, PRyMordial begins by computing the thermal background from first principles.As a byproduct of the determination of N eff and Ω ν h 2 , the relationship between time, scale factor and temperature of relativistic species is determined precisely, including effects from non-instantaneous decoupling within and beyond the Standard Model.
Next, PRyMordial evaluates the weak rates for neutron freeze out via a state-of-the-art implementation that includes nucleon finite-mass effects, one-loop QED corrections and finite-temperature effects.While the latter are typically negligible within current observational precision and can be conveniently stored between runs, the remainder are generally recomputed for each iteration of a generic BBN analysis.
Finally, PRyMordial solves a network of nuclide reactions for their yields within three different physical regimes: i) a high-temperature era in which one can restrict the study to nucleons with an initial temperature of O(10) MeV and a final temperature close to neutrino decoupling; ii) a mid-temperature era from O(1) MeV down to O(0.1) MeV, during which photo-dissociation of nuclear bound states is relevant; iii) and a low temperature era starting at O(0.1) MeV during which PRyMordial follows all of the nuclear species of interest, which ends at a temperature well below e ± heating of the thermal bath, i.e. down to O(1) keV.Local thermal equilibrium sets the initial nuclide abundances and detailed balance determines all of the reverse reactions included in the chosen set of nuclear reactions.These three regimes are matched such that the solution for each one provides the initial conditions for the successive period.
PRyMordial is a Python package with optional dependencies which allow more advanced users to speed up execution by exploiting the Julia programming language.The recommended libraries and general requirements are tabulated in Appendix A. As highlighted in Figure 1, PRyMordial is organized in five primary modules: • PRyM init.py is an initialization module where physical constants and Boolean flags for usercontrolled options are defined; in particular, three main blocks for input parameters are found: ⋆ Fundamental constants, masses (in natural units), initialized according to the PDG [42]; 11⋆ Additional parameters needed for the evaluation of the n ↔ p rates beyond the Born level; ⋆ Cosmological inputs including the CMB temperature and the abundance of baryonic matter [24].
Boolean flags allow the user to switch on/off the following options: • verbose flag: Allows the user to run the code with all of the internal messages enabled; • numba flag: If True, speeds up some numerical integrations, if the Numba library is installed; • numdiff flag: If True, performs numerical derivatives using Numdifftools library; • aTid flag: Controls the inclusion of incomplete-decoupling effects in the determination of the scale factor as a function of time and temperature; • compute bckg flag: If True, recomputes the thermodynamical background as presented in subsection 2.1 (via save bckg flag the outcome can be stored in a file for future runs); • NP thermo flag: If True, includes the contribution(s) of new (interacting) species to the dynamics of the thermal bath (by default, one must also provide a NP temperature); • NP nu flag: If True, includes new species thermalized with the neutrino bath; • NP e flag: If True, includes new species thermalized with the plasma; • compute nTOp flag: If True, recomputes weak rates beyond Born as discussed in subsection 2.2 (via save nTOp flag the outcome can be stored in a file for future runs); • nTOpBorn flag: If True, adopts the Born approximation for the neutron freeze out; • compute nTOp thermal flag: If True, recomputes thermal corrections to n ↔ p rates via Vegas (since this is numerically intensive, we recommend save nTOp thermal flag = True); • tau n flag: If True, uses the neutron lifetime to normalize the weak rates, see subsection 2.2; • NP nTOp flag: If True, includes NP affecting n ↔ p rates in units of the Born rates; • smallnet flag: If True, restricts the nuclear network to the set of 12 key nuclear processes collected in Table 1 of Appendix B; • nacreii flag: If True, the key nuclear rates adopted in PRyMordial will be mostly based on NACRE II compilation rather than those of PRIMAT, see subsection 2.3; • NP nuclear flag: If True, shifts the nuclear rates due to NP in units of the standard ones; • julia flag: If True, solves all of the systems of ordinary differential equations using routines in the SciML kit [58] developed for the Julia programming language; the optional dependencies described in Appendix A are then required.
This module also loads the tabulated nuclear rates (as well as the coefficients of Eq. (2.19)).
• PRyM thermo.py is the module where all of the thermodynamical quantities for the species contributing to the expansion of the Universe during radiation domination are defined, together with all the collision terms that enter in Eq. ( 2.3) and Eq.(2.7).
• PRyM nTOp.py is the module which imports the weak rates for n ↔ p conversion described in subsection 2.2, either relying on the additional module PRyM evalnTOp.py-where the actual computation of the rates is performed from scratch -or by loading pre-stored rates from a file.
• PRyM nuclear net12.pyand PRyM nuclear net63.pyare the modules which set up the systems of ordinary differential equations -see Eq. (2.16) -involving the nuclear rates loaded by PRyM init.py.The Boolean flag smallnet flag controls whether PRyMordial sets up and solves the smaller network of 12 key reactions or the full set of 63 nuclear processes.
• PRyM main.py is the main module, which calls the other modules to solve for the thermodynamical background, compute N eff and the cosmic neutrino abundance, and solve for the nuclide yields.It contains the Python class PRyMclass(), designated to return all the cosmological observables implemented in the package.
• PRyM jl sys.py is an optional module which allows the user to solve all of the systems of ordinary differential equations in PRyM main.py by taking advantage of the numerically efficient routines that are part of the SciML kit [58] developed in Julia.In some cases, this significantly speed up the execution time of the code (to a degree depending on both the adopted precision of the computation as well as the specific choice of differential-equation solver).
After downloading PRyMordial, the code can be used immediately.To run a Hello, World!style example, the user would enter the package folder, start an interactive Python session, and type: # Hello, World! of PRyMordial import PRyM.PRyM_main as PRyMmain res = PRyMmain.PRyMclass().PRyMresults() which executes a BBN computation and fills the array res with the values of: Located in the same folder are: • a folder PRyM in which all of the modules described above reside; • a folder PRyMrates in which all the essential thermal, weak and nuclear rates are present, and where new evaluations of them can be stored; • a script named runPryM julia.pythat provides a simple example for the user as to how to use the package, with execution-time benchmarking in both standard and Julia modes.
In the following subsections we present more sophisticated examples illustrating some of PRyMordial's capabilities.

SM examples: the PDG Plot and Monte Carlo Analysis
In an interactive session in Python, any default value in PRyM init.py can be changed using the syntax: import PRyM.PRyM_init as PRyMini # New assignment x for parameter X PRyMini.X = x This includes the Boolean flags listed in the previous subsection.Hence -to perform a run with: i) the computation of the thermal background from scratch, including non-instantaneous decoupling effects; ii) the ab-initio evaluation of the weak rates for neutron freeze out; and iii) the inclusion of key nuclear processes based on the tabulated rates of the NACRE II compilation -one should type: Li/H Figure 2: Primordial abundances of helium-4, deuterium, helium-3, and lithium-7 as predicted by PRyMordial within the SM, as a function of the cosmic baryon density.Central predictions are shown without theory uncertainties (i.e. using the nominal nuclear rates for the largest set implemented in the package with the NACRE II compilation for the key processes) and at the central values of all of the inputs.Measurements of light-element abundances (orange) as well as the CMB constraint on the baryon-to-photon ratio (cyan) follow from Figure 24.1 of the PDG [42].
# No need to recompute n <--> p rates as well PRyMini.compute_nTOp_flag= False PRyMini.save_nTOp_flag= False # Compute PRyMordial observables: now faster!import PRyM.PRyM_main as PRyMmain res = PRyMmain.PRyMclass().PRyMresults()While it may be necessary in general to recompute the thermal background and/or the rates for neutron freeze out, there are cases for which storing the outcome of these computations can be computationally advantageous.An example is the classic PDG review BBN plot of the primordial abundances as a function of the baryon-to-photon ratio η B [42].Once thermal background and weak rates have been stored, the behaviour of the abundances in the PDG Figure 24.1 can be reproduced with PRyMordial: # PDG plot npoints = 50 import numpy as np etabvec = np.logspace(-10,-9,npoints)# Initialization of array of observables YP_vec, DoH_vec, He3oH_vec, Li7oH_vec = np.zeros((4,npoints))for i in range(npoints): # Update value of baryon-to-photon ratio and store new obs The outcome of this code is illustrated in Figure 2, which adopts the largest nuclear network for the most accurate prediction of the relative abundance of lithium-7.It is worth noting that the BBN prediction for deuterium matches observations of quasar absorption systems, and is also in line with the cosmological abundance of baryons independently determined from the CMB (without a BBN prior).As pointed out in Ref. [54] and further scrutinized in Ref. [35], this test of concordance would fail if the PRIMAT rates were to be adopted, i.e. nacreii flag = False.
To perform a Monte Carlo analysis of the SM predictions taking into account uncertainties (similar to the one presented in Ref. [35]): # SM MC run num_it = 10000 # number of iterations import numpy as np Yp_vec, YDoH_vec, YHe3oH_vec, YLi7oH_vec = np.zeros((4,num_it))# Import PRyM modules import PRyM.PRyM_init as PRyMini import PRyM.PRyM_main as PRyMmain # Baryon eta from Planck 18 (no BBN prior) mean_eta0b = PRyMini.eta0bstd_Omegabh2 = 2*1.e-4std_eta0b = PRyMini.Omegabh2_to_eta0b*std_Omegabh2 # Neutron lifetime from PDG 2023 mean_tau_n = PRyMini.tau_nstd_tau_n = 0.5 Figure 3: 1D probability distributions (and 2D joint 68% and 95% probability regions) for the light primordial abundances predicted in the SM with PRyMordial.Predictions are obtained using a Gaussian prior for the neutron lifetime τ n = 878.4±0.5 s (comprising the eight best measurements from ultra-cold neutron experiments combined in Ref. [42]), and the cosmic baryon density, Ω B h 2 = 0.02230 ± 0.00020 (from Table 5 of Ref. [24] for the analysis with an uninformative Y P prior).The large network of nuclear reactions has been used, implying an additional 63 nuisance parameters varied with a log-normal distribution.Two different sets of key nuclear rates have been considered on the basis of the Boolean flag nacreii flag, and the statistics of the marginalized distributions for each case is presented.

# Compute primordial abundances at each iteration def ComputeAbundances(i):
# for i in range(num_it)) Yp_vec,YDoH_vec,YHe3oH_vec,YLi7oH_vec = np.array(FinalAbundances).transpose() The output maps out the probability distributions, shown in Figure 3, where the light elements at the end of the BBN era are predicted within the SM via a MC analysis that involves: i) a cosmological prior on the cosmic baryon abundance; ii) a particle-physics measurement prior on the neutron lifetime; and iii) a dedicated treatment of the uncertainties in the rates of the nuclear processes.Figure 3 displays the "deuterium anomaly" present for the PRIMAT compilation of the key nuclear rates, and further shows that it is completely washed out when one employs the NACRE II database 12 .Figure 3 suggests that the "primordial lithium problem" stands out as statistically significant, regardless of the approach undertaken for the nuclear network.However, the up-to-date analysis of the lithium problem in Ref. [94] points out that the predicted primordial abundance of lithium-7 could be depleted via stellar (and cosmic-ray) nucleosynthesis.Given this argument, the observational inference of Figure 2 and Figure 3, in which the observations lie below the theoretical prediction for primordial lithium-7, are consistent with a resolution for this long-standing puzzle.

NP examples: New Interacting Sectors and BBN
PRyMordial allows the user to perform state-of-the-art analyses for Physics beyond the SM in the Early Universe.A few options already built-in to the current release include: • additional relativistic degrees of freedom contributing to the expansion rate of the Universe in the form of a shift of N eff , see Eq. (2.5); • a non-zero chemical potential for neutrinos, influencing both the cosmological expansion rate as well as the equilibrium distributions in the weak processes for neutron-to-proton conversion; • Boolean flags specific to the study of new species interacting with the plasma and/or neutrino bath, as well as flags implementing a new entire sector with temperature T NP ̸ = T γ,ν ; • a Boolean flag and a dedicated parameter encoding NP effects as a phenomenological modification of n ↔ p conversion rates (in units of the Born rates); • a set of parameters that allow one to similarly investigate NP effects in the nuclear processes as a simple shift in terms of the median rate of each process.
The first two have been extensively investigated in Ref [35], and thus we focus here on the others.
The following is code demonstrating how to implement an electrophilic species in thermal equilibrium with the SM bath during BBN:  Figure 4: Investigation of the cosmological impact at the end of the BBN era from a new relativistic species X with degrees of freedom corresponding to a real / complex scalar (light / dark-blue lines), a real massive vector (magenta), or a Majorana / Dirac fermion (red / green); X is assumed to be in thermal equilibrium with either the electron-positron-photon plasma (left panels) or with the SM neutrino thermal bath (right panels).The orange bands represent the observational constraints at the 2σ level from Refs.[24,42].Predictions with PRyMordial are obtained at nominal inputs and rates.
In Figure 4 we present the results for NP scenarios of this type, reproducing the qualitative features already well-discussed, e.g., in Ref. [31].In particular, we observe three primary NP effects: i) a change in the cosmological expansion rate, affecting the time-temperature relation; ii) an impact on the evolution of the neutrino-to-photon temperature ratio, relevant for both neutrino and neutron decoupling; and iii) additional entropy released in the plasma, altering the number of baryons per a given baryon-to-photon ratio.Note that in Figure 4 we use the set of nuclear reactions from PRIMAT (nacreii flag = False) and as a result a neutrinophilic species around ∼ 10 MeV in mass appears to be favored by current observations of primordial D/H while remaining compatible with the other cosmological NP probes based on helium-4 and N eff .
In contrast to the previous scripts, this code calls PRyMclass() with three functions (of temperature) as arguments: the contribution to the energy density, its derivative, and the pressure of the new species added to the bath.More generally, one can include a new interacting sector with its own temperature T NP and non-trivial collision term δC NP along the lines of the recent work in Ref. [34].In PRyMordial one may study such "dark sectors" consistently by generalizing the set of equations in Eq. ( 2.3) to follow T NP together with T γ,ν , and solving for the entropy density involved in Eq (2.8) taking into account the effect of the NP.To do this, one switches on the Boolean flag NP thermo flag and codes all of the relevant contributions to the energy density, its derivative (which can optionally be evaluated numerically via Numdifftools), pressure and collision term for the NP sector, and passes them to PRyMresults.
One can also study NP resulting in changes to the weak rates for neutron freeze out and/or any of the implemented thermonuclear rates.To modify the weak rates, one sets the Boolean flag NP nTOp flag = True and change the parameter NP delta nTOp from its default of zero.Analogously, for the nuclear rates one switches on the flag NP nuclear flag and modifies the value of NP delta R with R being the reaction of interest.
As an example, we consider NP which results in a small change to the n ↔ p conversion rates.We perform a Bayesian fit to Y P and D/H (as quoted by the PDG [42]) and allowing τ n , Ω B h 2 , and the other key nuclear rates to vary within their uncertainties (in line with the SM MC analysis of the previous subsection):  3) and flags smallnet flag and nacreii flag are both switched on.Helium-4, deuterium measurements correspond to the recommended values from the PDG [42].
Figure 5 shows the resulting 2D joint (68% and 95%) probability regions for NP delta nTOp correlated with the measurements of primordial helium-4 and deuterium.To perform the statistical analysis, we adopt the emcee package [56].For the sake of computational efficiency, we restrict the analysis to the network of 12 key reactions (with nacreii flag = True), as is sufficient given the focus on helium-4 and deuterium.Figure 5 indicates that BBN is consistent with NP in the n ↔ p conversion rates at the level of at most a few percent relative to the standard Born rates.The tight correlation with Y P illustrates the importance of neutron freeze out in determining the primordial helium-4 abundance.

Outlook
In this work we have presented PRyMordial: A new tool to explore the physics of BBN in great detail, with an unprecedented eye toward applications for physics beyond the SM.The package also allows for fast, user-friendly precision analyses of the BBN era within the SM of Particle Physics, reaching the same level of accuracy as the state-of-the-art codes already publicly available.
In section 2 we have provided in some detail a review of the BBN era, highlighting the physics implemented in the code.The main novelties in PRyMordial are that it is: • A package entirely written in Python, easy to install, run and modify, efficient in the evaluation of the key quantities for the study of BBN; moreover, an optional dependence on Julia allows the user to make the code run even faster; • A computation of the thermal background based on the Boltzmann equations governing the evolution of the relativistic species present at that time.This allows for an accurate prediction of N eff from first principles and opens up new avenues for the study of BSM Physics; • A fast and accurate evaluation of the weak rates including QED, nucleon-finite mass and thermal corrections for a prediction of the neutron-to-proton ratio that confronts the precision of current and next-generation measurements; • A BBN code that easily allows exploration of uncertainties and changes in all of the input parameters and most importantly, includes by default different treatments for the nuclear rates in order to give to the user a better handle on the overall theoretical systematics.
In section 3 we describe the structure of the code and provide several examples of its usage within the Standard Model and for a few interesting scenarios of NP.
There are many directions that can be pursued in the future to make PRyMordial an even more compelling and flexible tool for the community.One important aspect we plan to expand upon is the characterization of the thermal background.At the moment, only a single common temperature for neutrinos is considered and no evolution equation for primordial chemical potentials is given by default.All of these can be easily implemented along the lines of Ref. [53].
Also relevant for precision studies would be an approach to efficiently include effects from phasespace spectral distortions of relativistic species.In this regard, we plan to further enrich the physics in PRyMordial with a dedicated framework for neutrino decoupling that includes effects from oscillations at non-zero lepton chemical potentials, see Ref. [110].
It would be a very interesting (though formidable) task to further improve the current next-toleading order computation of neutron freeze out in the Early Universe, filling in the gaps of some of the approximations undertaken in the literature (see Appendix B of [111] as well as the improvements brought by the recent effective-field-theory study at zero temperature of Ref. [77]).We would also eventually like to include higher-order QED corrections such as the ones available in Refs [64] and [112], as well as the NLO QED corrections to e + e − ↔ ν ν matrix elements recently inspected in Ref. [113].
Finally, in the future we would like to enlarge the nuclear network beyond the 63 nuclear reactions currently implemented, which encode all of the processes involving nuclides up to boron-8 in atomic and mass number (needed for an accurate prediction of lithium-7 in the Standard Model).
With the public release of PRyMordial we hope to provide to the community an important new tool to address fundamental questions about the Early Universe, whose study remains central to further progress in our understanding of Nature.In the wise words of a giant of our time [1]: "[Human beings] are not content to comfort themselves with tales of gods and giants, or to confine their thoughts to the daily affairs of life; they also build telescopes and satellites and accelerators, and sit at their desks for endless hours working out the meaning of the data they gather.The effort to understand the universe is one of the very few things that lifts human life a little above the level of farce, and gives it some of the grace of tragedy."import julia julia.install()import diffeqpy diffeqpy.install()At this point the user will be able to exploit the SciML routines developed in Julia to solve the nuclearreaction network in PRyMordial, speeding up the execution of time by a factor of two or more, and with the possibility of cherry-picking from a large collection of differential-equation solvers built-in in the package, see the documentation at https://docs.sciml.ai/DiffEqDocs/stable/solvers/odesolve.
To use the SciML routines, the user must set the flag PRyM init.flagjulia = True.In some systems, the very first call of PRyM main.PRyMresults() might need to be in Python and therefore requires initially PRyM init.flagjulia = False.Also, the first call in Julia will inevitably be slow, since it will compile PRyM jl sys.py.As a concise example of the dedicated script runPRyM julia.pycoming with the present release, here below is how things should work in the Julia mode: import PRyM.PRyM_init as PRyMini import PRyM.PRyM_main as PRyMmain # Initialization call in Python: PRyMini.julia_flag= False res = PRyMmain.PRyMclass().PRyMresults() # First call in Julia will be slow: PRyMini.julia_flag= True res = PRyMmain.PRyMclass().PRyMresults() # From here on, any call will be fast!

B Nuclear Processes in PRyMordial
In this appendix we collect the 12 key reactions necessary to accurately predict helium-4 and deuterium, see Table 1, as well as the 51 additional reactions comprising the full set recommended for a more robust prediction of lithium-7, Table 2.For the general aspects of the evaluation of the nuclear rates in the Early Universe as well as the theoretical and statistical details behind the compilation of the nuclear rates present in PRyMordial, we refer the interested reader to Ref. [92] and Refs.[107,108].The key nuclear reactions adopted in PRyMordial, with corresponding references.The red (blue) column refers to the option nacreii flag = True (False), see subsection 2.3 for further details.Notice that the compilation of the blue column is present also in the code PRIMAT [47].

Figure 1 :
Figure 1: PRyMordial in a nutshell: Schematic of the modules making it up and their inter-relations.
[56]traint on a relative change of the weak n ↔ p conversion rates from NP, based on a Bayesian fit performed with PRyMordial with the use of the emcee[56]package.Gaussian priors on the neutron lifetime and the cosmic baryon abundance are assumed (as for Figure SMFigure 5: