1 Introduction

1.1 Lithium Problem

Big Bang Nucleosynthesis (BBN) is one of nucleosynthesis categories, which currently faces a great challenge called the lithium problem. The lithium problem is the difference between the initial lithium abundance obtained from the observations of population II stars in the galaxy compared to its theoretical value in the Big Bang nucleosynthesis calculations (Fig. 1). Establishing an agreement between these two values is considered as one of the proofs of the correctness of the Big Bang Theory (Iliadis et al. 2016). This discrepancy is such that the theoretical value is three times the estimated experimental value (Iliadis and Coc 2020).

Fig. 1
figure 1

Schematic comparison of the observed abundances (Black Box) of (a). Deuterium (Left panel), (b). 3He (Middle Panel), and (c). 7Li (Right Panel) with the simulated abundances of BBN (Red diagonal-hatched box). A different scheme of lithium problem can be seen in (Fields, & Sarkar 2006)

Full size image

In Figure 1, the black box at each panel shows statistical uncertainty. The green vertical-hatched band shows the amount of baryonic matter in the universe measured by Cosmic Microwave Background Radiation (CMB) as the first pillar of Big Bang theory, and the red diagonal-hatched band shows the predictions or calculations of Big Bang theory for the abundance of elements. According to Figure 1, matching of the narrow bands of the confidence limit (CL), the range of changes of the Big Bang nucleosynthesis predictions and the observed abundance of light nuclei means the accuracy of the abundance calculations. For 4He and D, the mentioned match can be observed, but for the abundance of 7Li, there is no matching area, and vertical bands of either BBN or CMB fall outside the observed area. Among the explanations for the lithium problem are the lack of exact understanding of lithium processes in stellar environments, the systematic errors in 7Li observations, errors related to calculated thermonuclear rates that take part in the lithium and beryllium synthesis, or the need of a physics beyond the standard model (Iliadis and Coc 2020).

Solving the lithium problem may be possible through improving the reaction rate calculations of big bang nucleosynthesis. Hence, using the so-called Bayesian inference to perform the calculations, this study will focus on the study of d(p,γ)3He and 3He(d,p)4He reactions in this field, which affect the abundance of lithium.

1.2 Importance of the Investigated Reactions

Among twelve Big Bang nucleosynthesis reactions (Fig. 2) that eventually end up with the 7Li nucleus and somehow affect the lithium problem, the d(p,γ)3He and 3He(d,p)4He reactions have been selected for this study.

Fig. 2
figure 2

The network of twelve Big Bang Nucleosynthesis reactions between the rectangular points according to Cook et al. (2013). The small round points are not part of the network

Full size image

The 3He(d,p)4He reaction affects the abundance of primary deuterium and has a noticeable effect on the primary abundance of 3He and 7Li. Therefore, despite the indirect connection between this reaction and the lithium problem, the improvement of its reaction rate could have a meaningful impact on the abundances of the Big Bang nucleosynthesis. The d(p,γ)3He reaction as the second step in the pp hydrogen burning chain in stars, with its very high rate compared to its rival reactions, plays an important role in the early stages of star evolution and also in Big Bang nucleosynthesis by effecting the abundance of D, 3He and 7Li at astrophysical energies of about 100 keV.

1.3 Bayes’ Principal

In the theoretical context, among the most recent methods (such as methods on the basis of χ2 minimization) to improve the accuracy of the theoretical extrapolations is the Bayes' principle to infer an estimated value of a quantity in experimentally un-achievable areas using the existing experimental values of achievable areas. This principle follows as (Brewer 2020):

$$ P\left( {\theta |y} \right) = P\left( {y|\theta } \right)\frac{P\left( y \right)}{{P\left( \theta \right)}} $$
(1)

in which \(P\left( y \right)\) is a certain coefficient called prior that is our early belief about a quantity, \(P\left( \theta \right)\) is the probability density function (PDF) for the studied parameter \(\theta\) (evidence or normalization coefficients), \(P\left( {y{|}\theta } \right)\) is the probability of the quantity y provided that θ has occurred (likelihood), and \(P\left( {\theta {|}y} \right)\) is the updated result that is our final belief about the parameter (posterior). The likelihood, in nuclear astrophysics can be defined as a product of the PDF for the astrophysical S-factor of the model and the statistical uncertainty (or σstat standard deviation) (Esmaeili et al. 2022).

According to the selection of the Lognormal probability density function in this study to calculate the desired quantities, its relation is stated in Eq. 2, which includes the two components of “quantity location” (μ) and “shape parameter” (σ) (defined, respectively in Eqs. 3 and 4):

$$ f\left( x \right) = \frac{1}{{\sigma \sqrt {2\pi } }}\frac{1}{x}e^{{\frac{{ - \left( {\ln x - \mu } \right)^{2} }}{{\left( {2\sigma } \right)^{2} }}}} ,\,\,\, x > 0 $$
(2)
$$ \mu = \ln \left( {E\left[ x \right]} \right) - \frac{1}{2}{\text{ln}}\left( {1 + \frac{V\left[ x \right]}{{E\left[ x \right]^{2} }}} \right) $$
(3)
$$ \sigma = \sqrt {{\text{ln}}\left( {1 + \frac{V\left[ x \right]}{{E\left[ x \right]^{2} }}} \right)} $$
(4)

where eμ is the median value of the distribution, V[x] and E[x] are the variance of the natural logarithm distribution and the expected mean value for the variable x, respectively, and the factor uncertainty “f.u.” of the studied quantity is eσ (Hilbe et al. 2017). An example for a lognormal probability distribution having the quantity location (μ) and the shape parameter (σ) can be seen in Fig. 3.

Fig. 3
figure 3

A random lognormal probability distribution with two components: location of quantity (μ) and shape parameter (σ) according to Thomopoulos and Johnson (2003)

Full size image

The lognormal distribution follows (Hilbe, et al. 2017):

$$ \pi \left( f \right) = \frac{1}{{\ln \left( {f.u.} \right)\sqrt {2\pi f} }}e^{{ - \frac{{\left[ {\ln f} \right]^{2} }}{{2\left[ {\ln \left( {f.u.} \right)} \right]^{2} }}}} $$
(5)

where f is the normalization factor (normal probability distribution of systematics) and f.u. is the factor uncertainty that was explained earlier. The above Eq. 5 is written and recognized in R-programming as follows:

$$ f\sim {\text{LogNormal}}\left( {0,\left[ {\ln \left( {f.u.} \right)} \right]^{2} } \right) $$
(6)

The symbol “ ~ ” means that the desired quantity is sampled from the stated function. When one writes Eq. 6 in R, the program itself recognizes the mathematical form of it as Eq. 5, and implements it using the imported numbers of “zero” and “f.u.” in the Eq. 5 for further calculations. Therefore, there would be no need to be aware of the exact Eq. 5.

1.4 Uncertainties

All experiments, including the reaction rate, face with statistical (σstat) and systematic (ξsys) uncertainties. In general, the result of a reaction rate measurement or astrophysical S-factor can be expressed as \(S_{{ \pm {\text{sys}}}}^{{ \pm {\text{stat}}}}\), in which the statistical uncertainty (or standard deviation) and systematic uncertainty appear as indices in front of the reported quantity.

2 Literature Review

As a brief literature review on this field, Table 1 summarizes some of most important pioneer studies performed using Bayesian approach in this field up to year 2022.

Table 1 Summary of nuclear astrophysics studies using Bayesian method for reaction rates
Full size table

All of the above studies in this field have used relatively identical relations and methods, and only minor parts of their studies differ from each other, such as the number of experimental data points as inputs of the model, as well as the type of their MCMC algorithm. Of course, in all these studies, compared to other studies using other calculation methods, better and more accurate results (having less uncertainty in simulations) have been obtained.

3 Relations and Method

In this study, to calculate the exact value of a quantity, the existing experimental data for that quantity are used and then by attributing a specific probability distribution to that quantity, the behavior of that quantity (probability distribution of the quantity) is updated using Bayesian inference. The Bayes method is performed with the help of Monte Carlo sampling for numerous experimental input quantities (for example, resonance energy and intensity, partial widths and reduced widths). The process output, which is a probability density for the studied quantity, is used to statistically extract meaningful estimates. In this study, similar to the mentioned studies of Table 1, the used MCMC algorithm is on the basis of the Metropolis–Hastings algorithm.

The models of this study for both reactions include two sets of parameters: physical parameters (resonance energy, reduced widths and channel radii) and experimental parameters (data uncertainty and normalization coefficients). Like De Souza et al. (2020), the Normal densities of the priors of the physical parameters are truncated at zero. For interacting nuclei, the radius of the channel is obtained by:

$$ {\text{a}}_{c} = {\text{r}}_{0} \left( {A_{1}^{\frac{1}{3}} + A_{2}^{\frac{1}{3}} } \right) $$
(7)

where A1 and A2 are the mass numbers of two interacting nuclei and r0 is a constant number between 1 and 2 fm (De Souza et al. 2019a). The Equations necessary to fit the data for astrophysical S-factor including S(0), S'(0) and S"(0) are considered as below in energies less than 1.1 MeV (Iliadis et al. 2016):

$$ S\left( E \right) = S_{{{\text{bare}}}} \left( E \right).\exp \left( {\pi \eta \frac{{U_{e} }}{E}} \right) $$
(8)
$$ S_{{{\text{bare}}}} \left( E \right) = S\left( 0 \right) + S^{\prime}\left( 0 \right)E + \frac{1}{2}S^{\prime\prime}\left( 0 \right)E^{2} $$
(9)
$$ 2\pi \eta = 0.98951013Z_{0} Z_{1} \sqrt {\frac{{M_{0} M_{1} }}{{M_{0} + M_{1} }}\frac{1}{E}} $$
(10)

where Sbare is the S-factor of a nucleus without electron screening effect and Ue is the energy-independent electron screening potential (Iliadis et al. 2016). In the end, by numerical integration of the reaction rate equation [Eq. (11)] (Iliadis 2015) for the value of resulted S-factor from previous steps, the reaction rate is calculated for a specified temperature range.

$$ N_{A} \sigma \upsilon = \left( {\frac{8}{{\pi m_{01} }}} \right)^{\frac{1}{2}} \frac{{N_{A} }}{{\left( {kT} \right)^{\frac{3}{2}} }}\mathop \smallint \limits_{0}^{\infty } S\left( E \right).\exp \left( { - 2\pi \eta } \right)\exp \left( { - \frac{E}{kT}} \right){\text{d}}E $$
(11)

where m01 is the reduced mass of the interacting particles \(\left( {\frac{{m_{0} m_{1} }}{{m_{0} + m_{1} }}} \right)\), and kT is in terms of MeV, which corresponds to Giga Kelvin temperatures.

The cross section of the investigated reaction is obtained by fitting Eq. 9 for which we have (De Souza et al. 2019a):

$$ S_{{{\text{bare}}}} \left( E \right) = E.\exp \left( {2\pi \eta } \right).\sigma_{ab} \left( E \right) $$
(12)
$$ \sigma_{ab} \left( E \right) = \frac{\pi }{{k_{d}^{2} }}\frac{2J + 1}{{\left( {2J_{1} + 1} \right)\left( {2J_{2} + 1} \right)}}\left| {S_{ab} } \right|^{2} $$
(13)

where kd is the deuteron wave number, η is the Sommerfeld parameter, and Sab is the scattering matrix element with the following form for input nucleus of a and output b (De Souza et al. 2019a):

$$ \left| {S_{ab} } \right|^{2} = \frac{{{\Gamma }_{a} {\Gamma }_{b} }}{{\left( {E_{0} + {\Delta } - E_{r} } \right)^{2} + \left( {{\Gamma }/2} \right)^{2} }} $$
(14)

where Γa and Γb are the partial widths of the reaction input channel (a) and the output channel (b), Γ is the total width, Δ is the level shift, E0 is the energy eigenvalue of the level, and Er is the value of the resonance energy of the studied reaction, which are the main parameters of the matrix R and are obtained from the following equations (De Souza et al. 2019a):

$$ {\Gamma } = \mathop \sum \limits_{c} {\Gamma }_{c} , {\Gamma }_{c} = 2\gamma_{c}^{2} {\text{P}}_{c} $$
(15)
$$ {\Delta } = \mathop \sum \limits_{c} {\Delta }_{c} = {\Delta }_{1} + {\Delta }_{2} , {\Delta }_{c} = - \gamma_{c}^{2} {\text{P}}_{c} \left( {{\text{S}}_{c} - {\text{B}}_{c} } \right) $$
(16)

where γ2c is the reduced width and Bc is the boundary condition parameter for interacting nuclei. The energy-dependent quantities Sc and Pc are, respectively, the shift and the penetration factors for channel c (each of the two allowed reaction channels) that can be obtained using the Coulomb wave functions Fℓ and Gℓ by:

$$ {\text{P}}_{c} = \frac{{{\text{k}}_{d} {\text{a}}_{c} }}{{F_{\ell }^{2} + G_{\ell }^{2} }}, {\text{S}}_{c} = \frac{{{\text{k}}_{d} {\text{a}}_{c} \left( {F_{\ell }^{{}} F_{\ell }{\prime} + G_{\ell }^{{}} G_{\ell }{\prime} } \right)}}{{F_{\ell }^{2} + G_{\ell }^{2} }} $$
(17)

where ℓ is the orbital angular momentum of the reaction channel and ac is the radius of the reaction channel (De Souza et al. 2019a). The likelihood for this study, considering both statistical and systematic uncertainty, is defined by nested equations as follow:

$$ S_{i} ^{\prime}\sim {\text{Normal}} \left( {S_{i} \left( \theta \right),\sigma_{{{\text{extr}}}}^{2} } \right) $$
(18)
$$ S_{{\text{exp;i}}} \sim {\text{Normal}} \left( {\xi_{i} S_{i} ^{\prime},\sigma_{{{\text{stat}};i}}^{2} } \right) $$
(19)

Normal means a normal distribution function that is called in the R program environment. The symbol “ ~ ” means that the desired quantity is sampled from the stated function. The indices “extr” stand for the total unknown sources of the uncertainties including statistical uncertainty, systematic uncertainty or etc. The normalization factor for data set j has the following form:

$$ \xi_{j} \sim {\text{lognormal}} \left[ {0,\ln \left( {f.u._{j} } \right)^{2} } \right] $$
(20)

4 Results and Discussion

There are a few experimental data for S-factor in some limited areas of energy, which due to the nuclear coulomb barrier, cannot be obtained in other areas using the same experiments. Since each data set of a reaction is limited to a special energy interval, all the available datasets for a reaction are simultaneously used to find the curve of the studied quantity and also to estimate the value of S-factor in the missed areas using statistical fitting methods, such as extrapolations that are based on Bayesian principal. It must be mentioned that the mentioned data-points are all of the available valid data that are based on experiments.

However, to fit the astrophysical S-factor (equivalent to the cross section) of the d(p,γ)3He reaction in terms of energy, the numerical data of Ma et al. (1997), Schmid et al. (1997), Bystritsky et al. (2008) and Casella (2002) were used as the only data that have both statistical and systematic uncertainties. After the final fitting, the value of \(S_{bare - 0} = 0.23_{ \pm 0.089}^{ \pm 0.058} eV.b\) obtained, which is in very good agreement with the results of Iliadis et al. (2016) (\(S_{0} = 0.2156_{ \pm 0.036}^{ \pm 0.038} eV.b\)), with a slightly higher the amount of uncertainty in this study, in comparison, for example, with \(S_{0} = 0.214_{ - 0.16}^{ + 0.17} eV.b\) (Adelberger et al. 2011), which obviously represents higher uncertainties, while the results of the Bayesian method yields lower uncertainties, hence a better result. The used data to fit the astrophysical S-factor can be seen in Fig. 4.

Fig. 4
figure 4

The astrophysical S-factor (equivalent to the cross-section) for d(p,γ)3He

Full size image

As can be seen, the behavior of the S- is dominated factor by a broad resonance whose peak is around 100 keV energy, and as the energy decreases, its value drops to about 0.2 eV b. According to above figure, it can be seen that at an energy of about 30 keV, the S-factor data of Casella (2002) shows a much intense drop than the data of Schmid et al. (1997) in the same energy region, which is caused by the difference in their measurement conditions as well as the difference in the model considered by them, both of which can be found in the original works of them.

Using the resulting S-factor and placing it in the reaction rate formula (Eq. 11), the reaction rate values can be obtained as the diagram of Fig. 5 indicates.

Fig. 5
figure 5

Reaction rate changes at GK for d(p,γ)3He at temperatures less than 10 GK. The declined red point value of the reaction rate at about 9 GK indicates the non-linearity of the reaction rate at higher temperatures

Full size image

As can be seen in Fig. 5, with the increase in temperature, due to the increase in the kinetic energy of the interacting particles, the studied reaction particles will have a higher velocity and thus a higher rate and the graph grows. But at lower temperatures and energies, due to the energy drop of the interacting particles, which have less ability to penetrate the Coulomb repulsion barrier, the reaction rate decreases. These changes show a nonlinear behavior considering more complex models with more parameters. The different value of the reaction rate at the temperature of about 9 GK indicates the non-linear behavior of the reaction rate at higher temperatures, which depends on the conditions of the star environment and the abundance of the isotopes.

However, according to the numerical values obtained for the studied reaction {increase in the cross section and as a result, the increase of the reaction rate of d(p,γ)3He}, and by referring to Fig. 2, it can be seen that this increase will require more deuteron consumption, therefore, the amount of the conversion of the nuclei into lithium will reduce to a negligible extent, which confirms the reduction of the discrepancy between the theoretical and experimental values of lithium abundances, hence leading to an improvement in lithium problem. But due to the slight decrease of this value, it cannot be claimed that the existing problem has been completely resolved, but only an insignificant reduction (in the range of only a few percent compared to the difference of three times between the theoretical and experimental values) will be resulted, and the lithium problem will still remain, the reason for which should be found in other cases that were explained in introduction.

4.1 Results of the 3He(d,p)4He Reaction

To fit the astrophysical S-factor of the 3He(d,p)4He in terms of energy in the low-energy regions, numerical data of Zhichang (1997), Krauss et al. (1987), Möller and Besenbacher (1980), Geist et al. (1999), Costantini et al. (2000) and Aliotta et al. (2001) were used. After the final fitting, the value of \(S_{{{\text{bare}} - 0}} = 5.82_{ \pm 0.098}^{ \pm 0.098} \,\,{\text{MeV}}\,{\text{b}}\) was obtained, with a very small difference in comparison to the results of De Souza et al. (2019b) (\(S_{0} = 5.729_{ \pm 0.088}^{ \pm 0.097} \,\,{\text{MeV}}\,{\text{b}}\)), which shows a very good agreement, but the difference that the uncertainty value in this study is slightly higher. However, the result can be compared, for instance, with \(S_{0} = 5.9_{ - 0.3}^{ + 0.3} \,\,{\text{MeV}}\,{\text{b}}\) (Descouvemont et al. 2004), which obviously has a higher uncertainty. The behavior of S-factor by the data used to fit the astrophysical S-factor can be seen in Fig. 6.

Fig. 6
figure 6

The astrophysical S-factor (equivalent to the cross-section) for 3He(d,p)4He

Full size image

According to the Fig. 6, the behavior of the S-factor is dominated by a broad resonance at an energy of about 200 keV (the resonance peak at the energy of about 0.2 MeV) and as the energy decreases, its value drops to about 6 MeV.b (Sbare-factor without electron screening effect). But the remarkable thing about this reaction is that the effect of electron screening in the Bayesian model of this reaction causes the value of S-factor to increase significantly at energies outside the experimental range. The effect of electron screening resulting from this research is in full agreement with the results reported in several studies, including the study of Aliotta et al. (2001) and De Souza et al. (2019b) [the reported values for electron screening effects are about 120 eV to 220 eV for 3He-d and d-3He reactions, respectively), and that is why this reaction is considered as one of the most important reactions for investigating the electron screening effects in nuclear astrophysics.

It should be noted that the data of Zhichang (1997) and Möller & Besenbacher (1980) at very low energies, face a more severe drop than the data of other authors in the same energy region, which is caused by the difference in their measurement conditions, as well as the difference in the model considered by them. In other words, for example, the effect of electron screening can be the most important factor in this significant difference.

In the end, by using the resulting S-factor and placing it in the reaction rate formula (Eq. 11), the reaction rates can be obtained in terms of temperature as the diagram of Fig. 7.

Fig. 7
figure 7

Variations of the reaction rate at GK for 3He(d,p)4He

Full size image

As can be seen in Fig. 7, with the increase in the temperature, the kinetic energy of the interacting particles increases, and because of the high density of deuteron in the environment, the rate of the reaction as a result increases and the graph grows. But at lower temperatures and energies, due to the energy drop of the interacting particles, which have less ability to penetrate the Coulomb repulsion barrier, the reaction rate decreases. Of course, the value close to zero shown in the Fig. 7 is actually due to the insignificant rate of the reaction, and in reality it does not simply mean that the desired reaction does not occur. The rate changes of 3He(d,p)4He can be considered in a relatively exponential form, which has a non-linear behavior when considering more complex models with more components. Also, according to different amounts of 3He in the environment, a final limit can be considered for this reaction, or the reaction rate diagram can have different final behaviors (decrease, increase, or approaching a specific final value).

However, again, according to the resulted increase in the cross section and the reaction rate of 3He(d,p)4He, and by referring to Fig. 2, one can see that this increase will require more 3He consumption, therefore, the amount of the conversion of the interacting nuclei into lithium will decrease, which confirms the reduction of the discrepancy between the theoretical and experimental values of lithium abundances, hence leading to an improvement in lithium problem. But due to the negligible resulted decrease, it cannot be claimed that the mentioned discrepancy has been completely resolved. In contrary, only an insignificant reduction (in the range of only a few percent compared to the difference of three times between the theoretical and experimental values) will be resulted, and the lithium problem still remains.

5 Conclusion

Considering the numerical values obtained for the rate of the two studied reactions, which have a slight difference (only a few percent) compared to the values obtained from previous non-Bayesian calculations by other researchers, despite the innovativeness and efficiency of Bayes method, it cannot be claimed that the lithium problem has been completely resolved, but only an insignificant reduction (compared to the three times difference between the theoretical and experimental values) will occur in this difference, and the lithium problem will still remain, the reason for which must be sought in other cases, including the fact that these calculations have no direct impact on the abundance of lithium, and as a result, the reason for the mentioned discrepancy between the theoretical and experimental values of the initial abundance of lithium should be sought in other areas, including a new correct knowledge of the lithium destruction mechanism in stars or the need for new physical theory beyond the standard model.

However, despite the validated accuracy of the current calculations in comparison with results of the all references, the need to improve or complete these results by other models is still felt, and there is a need to include other factors in the defined models in order to try to fully solve the lithium problem.