Introduction

The hot Big Bang theory, first proposed in 1946 by Gamow [1], is now the most widely accepted cosmological model of the universe, where it has expanded from a very high-density state dominated by radiation. This theory has been confirmed by the observations of the cosmic microwave background [2-4], the expansion of the universe, and good global agreement between the predictions and observations of the primordial abundances of the lightest elements in nature: hydrogen, helium, and lithium. According to the Big Bang theory, the universe began with a fireball approximately 13.8 billion years ago. Following inflation and cooling, primordial big-bang nucleosynthesis (BBN) began when the universe was approximately 3 min old (when the temperature was reduced to approximately 1 GK, i.e., particle energy E≈kT≈0.1 MeV) and ended less than half an hour later when the nuclear reactions were quenched by the low temperature and density conditions in the expanding universe. Only the lightest nuclides were synthesized in appreciable quantities through BBN: approximately 75% ¹H and 25% ⁴He, with small amounts of ²H, ³He, and ⁷Li. These relics provide us with a unique window into the early universe. More comprehensive reviews on BBN can be found in the literature [5-7].

In general, the primordial abundances of ²H (referred to as D) and ⁴He inferred from observational data agree with predictions, except for the lithium problem [8, 9]. Deuterium, a fragile isotope, was destroyed after the BBN. Its most primitive abundance was determined from observing cosmological clouds at a high redshift on the line of sight of distant quasars. Recently, the precision of deuterium observations in cosmological clouds has dramatically improved, reaching an accuracy of 1.19% for primordial deuterium abundance (D/H), that is, D/H = (2.527±0.030)×10^-5[10]. Therefore, the nuclear cross section data relevant to the deuterium involved in the BBN network needs to be known with similar precision to constrain further the cosmological parameter, that is, the cosmic baryon density. In this study, we focused on BBN deuterium abundance.

In BBN, deuterium synthesis is the first to occur, and the accumulation of primitive deuterium significantly affects the rate of subsequent reactions. Therefore, accurate determination of the abundance of deuterium in BBN calculations is important. Eextensive studies relevant to BBN deuterium abundance have been conducted [11-13] but only a few mention the impact of the uncertainties of the relevant reaction rates on D/H uncertainty. In this study, we used the Monte Carlo simulation code PRIMAT [14] to study the D/H uncertainty within BBN reaction networks and demonstrated the reaction rate sensitivity for some major reactions. In addition, the uncertainties of both the nuclear physics input and cosmological parameters were proposed to satisfy the accuracy of the observed deuterium abundance.

BBN model & Monte Carlo method

For a standard BBN model, the evolution of the nucleosynthetic abundance can be obtained by solving the following system of differential equations [15]: $\frac{\text{d}{Y}_i}{\text{d}t}={\displaystyle \sum _{j\mathrm{,}k\mathrm{,}l}{N}_i}\left({\Gamma }_{kl\to ij}\frac{Y_l^{N_l}{Y}_k^{N_k}}{N_l!N_k!}-{\Gamma }_{ij\to kl}\frac{Y_i^{N_i}{Y}_j^{N_j}}{N_i!N_j!}\right),$ (1) where Ni represents the mass number of the corresponding nuclide, Yi is the abundance of the corresponding nuclide, and Γ denotes the reaction rate, which is usually obtained using the following formula [16, 9, 17]: $\langle \sigma v\rangle ={\left[\frac{8}{\pi \mu }\right]}^{1/2}{[kT]}^{-3/2}{\displaystyle \int }\sigma (E)E{e}^{-E/kT}\text{ d}E.$ (2) At the same time, the standard Big Bang model is also a single-parameter model, with only the parameter η used, that is, the ratio of the number of baryons to the number of photons; different η values correspond to different nuclide abundance evolution curves or frozen abundances. In the study of BBN, we usually obtain the corresponding abundance value by entering the observed η into the model as a parameter and then comparing it with the observed abundance. To perform such calculations, input parameters such as nuclear reaction rates, element abundances, and baryon density are required. However, because reaction network models incorporate an increasing number of nuclides, conventional numerical calculations have become increasingly computationally expensive when computing uncertainties. Furthermore, according to the research conducted by Longland et al. [18], the traditional method of error propagation fails to incorporate the statistical significance of errors. In this study, we adopted the Monte Carlo sampling technique to calculate the abundance of the targeted nuclide. The abundance uncertainty can be directly determined by generating a data distribution through the sampling process.

The Monte Carlo method is a statistical method based on the large number theorem. Converting the solutions of mathematical problems into random samples can significantly reduce the difficulty in solving complex models. Because of the sampling of specific physical quantities following a certain distribution (such as Poisson, Gaussian, and log-normal distributions), it also preserves the statistical significance of the physical quantities and their associated uncertainties and applies them realistically in the calculation results.

In this study, we used the PRIMAT code [14] for BBN calculations. PRIMAT is mainly divided into the following parts: first, it needs to determine some cosmological parameters, such as the parameter a(t) of the Friedmanan-Lemaitre (FL) spacetime, which in the code is obtained from the thermodynamics of the plasma, and its change with temperature is obtained by numerical inversion; second, it calculates the effects of some weak interactions and stores the relationship between them and the temperature in a hard disk; finally, a reaction network model is established to calculate the nuclear reaction processes of ≥10 GK, 10–1.25 GK, and ≤1.25 GK according to the temperature change. During the entire process, the code uses random sampling of the distribution of relevant parameters to obtain the distribution function of abundance and then determines the uncertainty of the abundances.

In particular, the main focus of BBN is nuclear physics input quantities, specifically the reaction rates of the relevant reactions. The PRIMAT code assumes that the reaction rate follows a lognormal distribution [14]: $f(x)=\frac{1}{\sigma \sqrt{2\pi }}\frac{1}{x}{e}^{-\frac{{\left(\ln x-\mu \right)}^2}{2{\sigma }^2}},$ (3) where x denotes the reaction rate, that is, $x=N_{\text{A}}\langle \sigma v\rangle$, μ and σ are parameters of the lognormal distributions. For this distribution, the corresponding low, median, and high rates are expressed as follows: $\begin{array}{ccc}{x}_{\text{low}}=e^{\mu -\sigma }\mathrm{,}& x_{\text{med}}=e^{\mu }\mathrm{,}& x_{\text{high}}=e^{\mu +\sigma },\end{array}$ (4) where $x_{\text{low}}$, $x_{\text{med}}$ and $x_{\text{high}}$ represent the reaction rates with probabilities of 16%, 50%, and 84%, respectively. Further information is available in Refs. [18].

Model parameters & Simulation

The main parameters in the PRIMAT code that can be adjusted are the reaction rates and the associated uncertainties of the 11 reactions of primary importance [19] involved in the BBN network to evaluate their impact on deuterium abundance. The reaction rates were adopted from the default rates used in the PRIMAT code, and the data sources are listed in Table 1. In addition, the default values set in the PRIMAT code were adopted, that is, the cosmological parameter was Ω_bh² = (0.02225±0.00016) [4] for the baryon-to-photon ratio, and the neutron lifetime was τ_n = (879.5±0.8) s [20]. Notably, in the original PRIMAT calculations [14], a calculated D/H uncertainty of 1.46% was obtained using the default reaction rates and their uncertainties, but without considering the uncertainties of the cosmological parameter Ω_bh² at that time. In our calculations, a similar D/H uncertainty of 1.47% was obtained under the same conditions, confirming the correctness of our calculations.

First, we studied the effect of uncertainties in the 11 reaction rates, Ω_bh² and τ_n on the calculated D/H uncertainties. The results are listed in Table 2. This shows that the uncertainty in the neutron lifetime $\text{Δ}{\tau }_{\text{n}}$ negligibly impacts the D/H uncertainty, whereas the uncertainties in the reaction rates and Ω_bh² have a much larger impact. Conversely, $\text{Δ}{\tau }_{\text{n}}$ had a significant impact on the ⁴He abundance. The impact of the neutron lifetime on Big Bang nucleosynthesis was studied in Ref. [26] and is not discussed here.

In the following, we investigate the sensitivity of the reaction-rate uncertainties to the primordial deuterium abundance (D/H) uncertainties by considering the uncertainties of Ω_bh² and τ_n simultaneously. The default reaction rate uncertainties adopted in PRIMAT for the four most important reactions are presented in Table 3.

To examine the sensitivity, we multiplied the default uncertainty of just one reaction rate by arbitrary factors of M.F. = 0.1, 0.8, 0.9, 1, 2, 4 and 6, respectively, and set the remaining 10 reactions with their default rate uncertainties. We found that the d(p,γ)³He, d(d,n)³He, d(d,p)t and p(n,γ)d reactions are the four major reactions that have the largest influence on the D/H uncertainty, compared to the remaining reactions. The calculated results are listed in Table 4 and shown in Fig. 1.

Fig. 1Full-screen View

Calculated D/H uncertainties v.s. multiplying factor of the reaction-rate uncertainty for the four most important reactions. The connecting lines are used just to guide the eyes. Table 4 can be referenced for details. The inset shows the enlarged plot for the multiplicative factors of M.F. = 0.1, 0.8 and 0.9. Here, the “Original” line (corresponding to a D/H uncertainty of 1.83%) indicates the results calculated with all the default rates and the associated uncertainties adopted in the PRIMAT code [14]

For multiplying factors of less than one (i.e., reducing the current rate uncertainty), the decreasing trend for the four reactions of interest is shown in the inset of Fig. 1. This demonstrates that the calculated D/H uncertainties can reach 1.41%, 1.76%, 1.77%, and 1.82% when the reaction rate uncertainties of d(p,γ)³He, d(d,n)³He, d(d,p)t and p(n,γ)d are reduced by a factor of 10 (i.e., $M\mathrm{.}F\mathrm{.}$ = 0.1, nearly no rate uncertainty). In general, the current default $d{(p\mathrm{,}\gamma )}^3$ He rate uncertainty (i.e., 3.7%) dominates the calculated D/H uncertainty because the default rate uncertainties adopted for the latter three reactions are already quite small (see Table 3). Thus, reducing their uncertainties further would not significantly affect the calculated D/H uncertainty. Therefore, the main source of D/H uncertainty originates from the d(p,γ)³He reaction rate. This conclusion is consistent with the findings of the recent LUNA studies [29]. Here, the “Original” line (i.e., M.F. = 1 case) shown indicate the results calculated with the default rates and their uncertainties adopted in the original PRIMAT code [14], i.e., corresponding to a D/H uncertainty value of 1.83%.

In addition, we found that changing the uncertainty of the reaction rate also affected the central value of deuterium abundance. We believe that this is because the sampling of the reaction rates in the code was performed according to the log-normal distribution. This distribution is a typical skewed distribution, where the larger the parameter σ, the more the most likely position is to be skewed to the left. We found that the change in the central value of the deuterium abundance caused by increasing uncertainty was equivalent to the change in the central value caused by the shift in the most likely position associated with uncertainty.

Access to the observational precision

In recent astronomical observations [10], the recommended value of the deuterium abundance is D/H = (2.527±0.030)×10^-5 with an accuracy of approximately 1.19%, which is better than the above PRIMAT calculated accuracy of 1.83%, that is, (2.459±0.045)×10^-5. Therefore, we focus here on the four reactions that have the largest influence on the D/H uncertainty and use a binary search to confirm whether reducing their rate uncertainties can improve the BBN calculated accuracy to the level of the observation. However, we found that the accuracy of the current deuterium abundance observations could not be achieved by reducing only the uncertainties of the reaction rates, and the uncertainty of Ω_bh² should also be reduced. Note that the d(p,γ)³He reaction rate is adopted from Ref. [22] for all the abovementioned calculations.

Most recently, Mossa et al. [29, 30] measured the d(p,γ)³He cross section in the $E_{\text{c}\text{.m}\text{.}}$ = 33–263 keV energy region using the LUNA 400 kV accelerator to an unprecedented precision of better than 3% by exploiting the million-fold reduction in cosmic-ray muons at Gran Sasso. Their new astrophysical S factor remarkably improved the evaluation of the present-day baryon density, Ω_bh², using the standard BBN model alone. In this study, we utilized the updated reaction rate of d(p,γ)³He from LUNA [29] to calculate the abundances and uncertainties of these primordial light nuclides. The results are listed in Table 5, where the abundance of ⁴He is expressed using its mass fraction $Y_{\text{p}}$, whereas the abundances of the other nuclides are expressed as the ratio of their number density to that of ¹H. We find that when utilizing the reaction rate of d(p,γ)³He reported by LUNA, the calculated D/H uncertainty can be reduced from 1.83% to 1.56%, but it does not reach observational precision. We found that this goal can be achieved by reducing the current uncertainty in Ω_bh² by approximately 50%.

Conclusion & outlook

In this study, we investigated the primordial deuterium abundance and its uncertainty using the BBN code PRIMAT. We found that the predicted deuterium abundance uncertainties were dominated by the reaction rate uncertainties in the four most important reactions of d(p,γ)³He, d(d,n)³He, d(d,p)t, and p(n,γ)d, as well as that in Ω_bh². Although the current BBN calculation can reach a deuterium abundance precision level of 1.56% with the recent precise LUNA d(p,γ)³He rate, there is still a gap of 0.4 % in the observational precision. We found that this gap cannot be largely reduced by only reducing the uncertainties of the reaction rates and the uncertainty of Ω_bh² should also be reduced. If the uncertainty of Ω_bh² adopted from the Planck 2015 results [4] is reduced by approximately 50%, the calculated D/H uncertainty can reach the observational level.

In nuclear physics, the reaction rate uncertainties for the remaining three reactions must be significantly reduced, except for the d(p,γ)³He reaction. For instance, relevant reactions can be measured directly at the China Jinping Underground Laboratory (CJPL) [31-33], which is the deepest operational underground laboratory for particle and nuclear physics experiments worldwide. With such a unique superlow background environment[34], several successful experimental campaigns [35-42] were conducted at the Jinping Underground Nuclear Astrophysics Experimental Facility (JUNA) [43, 44].

Figure 2 shows the observational and BBN model calculated primordial deuterium abundances and their uncertainties. The observed deuterium abundance was recommended as D/H = (2.527±0.030)×10^-5 [10, 14], whereas the currently calculated value was (2.471±0.039)×10^-5 and there was an approximately 1.5σ deviation in the central abundance value. In addition, the calculated values were less accurate than the observed values. Therefore, From a nuclear physics perspective, these important reactions still need to be measured precisely. Furthermore, the BBN calculated precision for other primordial nuclides can also be effectively reduced using more precise reaction rates, with which the cosmic parameters can be strictly constrained.

Fig. 2Full-screen View

(Color online) Observational and model calculated primordial deuterium abundances and the associated uncertainties. The data (from left to right) are adopted from Refs. [7, 10-12, 14, 45-53]. The current result shown is calculated by using the recent LUNA d(p,γ)³He rate, and by considering the uncertainties in Ω_bh² and τ_n

For instance, we collected all the available data for the important production reaction of p(n,γ)d, as shown in Fig. 3. This shows that there are still few experimental data in the 0.1–1.0 MeV energy region of the BBN interest. Furthermore, Hara’s et al. photoneutron data [58] still had relatively large error bars of (7–14)%. However, the PRIMAT code adopts only a very small (0.45–0.73)% reaction rate uncertainty for this p(n,γ)d reaction based on theoretical calculations [21]. Such a small uncertainty requires further validation using more precise experimental data in the BBN energy region. We are now planning to conduct photodisintegration measurements at the Shanghai Laser Electron Gamma Source (SLEGS) facility [63, 64], which uses the back-Compton scattering of electrons from the 3.5 GeV electron storage ring of the Shanghai Synchrotron Radiation Facility (SSRF) and the CO₂ laser to generate the γ beam in the energy range of 0.4~20 MeV [65]. The γ flux is expected to reach 10⁵~10⁷/s, depending on both the γ-ray energy and the collimator size [66, 67].

Fig. 3Full-screen View

(Color online) The reaction rate data for p(n, γ)d. The black solid line, blue dotted line, red solid line, green dash-dotted line, and purple dashed line are from FCZI [54], Hale in 1991 [55] and in 2003 [56], JENDL [57], and Hara in 2003 [58], respectively. The black square [58], the red circle [59], the blue triangle [60], the green inverted triangle [61] and the purple diamond [62] represent the available experimental data. The shaded area indicates the goal region with the SLEGS facility

In the proposed experiment, a 4π flat-efficiency ³He neutron detector array was intended for use as a photoneutron detector with a detection efficiency of approximately 42% [68] in the energy region of interest. An LaBr₃ scintillator detector was used for γ flux monitoring and waveform acquisition. We aim to measure the cross section in the energy region of 0.01–10 MeV (as shown in the shaded region in Fig. 3) with precision greater than 5% to further reduce the deuterium abundance uncertainty associated with the p(n, γ)d rate. As stated in Ref. [21], “the results of the R-matrix theory [69] at E = 0.1 and 1 MeV differ significantly from the other theoretical estimations by ≈4.6%. Therefore, it would be important to experimentally measure the np-capture cross sections at these energies to resolve this significant discrepancy”. Thus, precise data can also help constrain nuclear reaction models.

Finally, it is observed that different BBN codes generated quite different abundance values, maybe mainly owing to the different nuclear-physics inputs used. Therefore, the consistencies of different BBN codes can be verified using the same inputs. In addition, the Monte Carlo-type PRIMAT code always yields much smaller uncertainties than other codes; hence, the Monte Carlo method should be implemented in other codes to perform a cross-check. Moreover, it is important not only to check how the uncertainties in the reaction rates affect the D/H uncertainty but also how a shift in the central rate values may influence the D/H values. All these issues, which are beyond the scope of this study, need to be studied in the future.