title: Thermonuclear reaction rates and primordial nucleosynthesis
authors: Christian Iliadis Alain Coc
text: 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Christian Iliadis
  iliadis@unc.edu
  
  0000-0003-2381-0412]Christian Iliadis
  Department of Physics & Astronomy, University of North Carolina at Chapel Hill, NC 27599-3255, USA
  Triangle Universities Nuclear Laboratory (TUNL), Durham, North Carolina 27708, USA
  
  
  CNRS/IN2P3, IJCLab, Université Paris-Saclay, Bâtiment, 104, F-91405 Orsay Campus, France
  
  
  
  Assuming the best numerical value for the cosmic baryonic density and the existence of three neutrino flavors, standard big bang nucleosynthesis is a parameter-free model. It is important to assess if the observed primordial abundances can be reproduced by simulations. Numerous studies have shown that the simulations overpredict the primordial ^7Li abundance by a factor of ≈ 3 compared to the observations. The discrepancy may be caused by unknown systematics in ^7Li observations, poorly understood depletion of lithium in stars, errors in thermonuclear rates that take part in the lithium and beryllium synthesis, or physics beyond the standard model. Here, we focus on the likelihood of a nuclear physics solution. The status of the key nuclear reaction rates is summarized. Big bang nucleosynthesis simulations are performed with the most recent reaction rates and the uncertainties of the predicted abundances are established using a Monte Carlo technique. Correlations between abundances and reaction rates are investigated based on the metric of mutual information. The rates of four reactions impact the primordial ^7Li abundance: ^3He(α,γ)^7Be, d(p,γ)^3He, ^7Be(d,p)2α, and ^7Be(n,p)^7Li. We employ a genetic algorithm to search for simultaneous rate changes in these four reactions that may account for all observed primordial abundances. When the search is performed for reaction rate ranges that are much wider than recently reported uncertainties, no acceptable solutions are found. Based on the currently available evidence, we conclude that it is highly unlikely for the cosmological lithium problem to have a nuclear physics solution.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § INTRODUCTION
   
  
  
  The standard model of the big bang rests on the three observational pillars of big bang nucleosynthesis <cit.>, cosmic expansion <cit.>, and cosmic microwave background radiation <cit.>. Primordial nucleosynthesis took place during the first 20 minutes after the big bang and probes the cosmological model at early time. Nucleosynthesis in the early universe proceeded at temperatures and densities below 1 GK and 10^-5 g/cm^3, respectively. It produced the nuclides ^1H (or H), ^2H (or D), ^3He, ^4He, and ^7Li, which represent the building blocks for the subsequent evolution of matter in the universe.
  
  Primordial abundances have been measured in the proper astrophysical site with high precision in recent years. The cosmological deuterium abundance is determined by observing damped Ly-α systems at high redshift, resulting in a primordial number abundance ratio relative to hydrogen of (D/H)_p^obs = (2.527±0.030)× 10^-5 <cit.>. The nuclide ^3He has not been observed outside of the Galaxy since its abundance is very small. Also, the galactic chemical evolution of ^3He is uncertain because it is both produced and destroyed inside stars. Therefore, only an upper limit has been established from observations, (^3He/H)_p^obs ≤ (1.1±0.2)× 10^-5 <cit.>. The primordial ^4He abundance is deduced from measurements of metal-poor extragalactic H II regions, resulting in a nucleon fraction of Y_p^obs = 0.2449±0.0040 <cit.>. The cosmological ^7Li abundance is estimated from observations of low-metallicity stars in the halo of the Galaxy, where the lithium abundance is nearly independent of metallicity. The result is (^7Li/H)_p^obs = (1.58±0.31)× 10^-10 <cit.>. The uncertainties for the measured primordial ^2H, ^4He, and ^7Li abundances correspond to values of 1.2%, 1.6%, and 19.6%, respectively. These values are uncorrelated since each of these nuclides is observed in different astronomical objects by independent measurement techniques. The observations will be most useful for constraining cosmology if predicted primordial abundance uncertainties can be reduced to a level similar to that of the experimental ones. 
  
  The precision in the cosmological parameters entering the model of big bang nucleosynthesis has also been improved significantly in recent years. In particular, the measurement of anisotropies in the cosmic microwave background radiation has determined the cosmic baryonic density, Ω_b · h^2 = 0.02233±0.00015 <cit.>, with a precision of less than 1%. If we assume, in addition, N_ν = 3 for the number of neutrino flavors, the standard model of big bang nucleosynthesis becomes a parameter-free theory. These results establish the expansion rate of the early universe and the prevailing thermodynamic conditions (i.e., temperature and density). All that is required to simulate the primordial abundances is to numerically solve the network of coupled differential equations describing the abundance evolution of the light nuclides. The numerical results depend sensitively on the magnitudes and associated uncertainties of the thermonuclear reaction and weak interaction rates (see, e.g., <cit.>). 
  
  Many studies <cit.> have found that the predicted primordial abundances, with one exception, conform broadly with the experimental values. This result is significant because the abundances span nine orders of magnitude. Consequently, big bang nucleosynthesis is a sensitive tool for probing the physics of the early universe as well as for constraining physics beyond the standard model <cit.>. The exception is the primordial ^7Li abundance for which the simulated value exceeds the observational result by a factor of about three. This is the cosmological lithium problem and it has not found a satisfactory solution yet. The discrepancy may be caused by unknown systematics in the ^7Li observations, poorly understood depletion of lithium in stars, errors in thermonuclear rates of reactions that take part in the lithium and beryllium synthesis, or physics beyond the standard model. Solutions involving exotic physics almost always lead to an overproduction of deuterium incompatible with observations <cit.>. Solutions involving stellar physics require a uniform in situ lithium depletion <cit.> over wide metallicity
  and efffective surface temperature ranges <cit.>. See <cit.> for reviews.
  
  
  The goal of this work is to investigate the nuclear physics aspects of this persistent problem. In Section <ref>, we summarize the most recent information regarding the thermonuclear rates of key reactions. Our reaction network, thermodynamic conditions, and initial abundances are discussed in Section <ref>. Results from a network Monte Carlo procedure are presented in Section <ref>. In Section <ref>, we study correlations between reaction rates and simulated abundances. A genetic algorithm is employed in Section <ref> to investigate whether simultaneous rate changes of key nuclear reactions can account for the observed primordial abundances. A concluding summary is given in Section <ref>. In the Appendix, we provide information on the rates of selected nuclear reactions. 
  
  
  
  
  § NUCLEAR PROCESSES: BAYESIAN REACTION RATES
   
  
  
  The twelve most relevant nuclear processes taking part in primordial nucleosynthesis are listed in column 1 of Table <ref>, together with their rate uncertainty at 1 GK (column 2), the reference for the rates adopted here (column 3), and comments on the data treatment (column 4). 
  
  
  lcll
  Main nuclear processes in primordial nucleosynthesis.a 
  2
  
   Process        Uncertainty (%)b      Most recent reference     Comments  n ↔ p                     0.09           <cit.>                  theory and neutron lifetime 
   
     p(n,γ)d                            0.4            <cit.>                    MCMC analysis of data  
   
     d(p,γ)^3He                       3.7            <cit.>                   Bayesian data analysis  
   
     d(d,p)t                                   1.1            <cit.>                     Bayesian data analysis  
   
     d(d,n)^3He                              1.1            <cit.>                     Bayesian data analysis  
   
     ^3He(n,p)t                              1.6            <cit.>            χ^2 data fitting  
   
     ^3He(d,p)α                       1.2            <cit.>              Bayesian data analysis  
   
     t(d,n)α                            0.3            <cit.>              Bayesian data analysis  
   
     ^3He(α,γ)^7Bec           2.4            <cit.>                   Bayesian data analysis  
   
     t(α,γ)^7Li                4.3            <cit.>            χ^2 data fitting 
   
     ^7Be(n,p)^7Li                         2.1            <cit.>               Bayesian data analysis  
   
     ^7Li(p,α)α                2.5            <cit.>            χ^2 data fitting  
   
  aThe symbols p, d, t, and α denote the nuclides ^1H, ^2H, ^3H, and ^4He, respectively.
  bRate uncertainty at 1 GK.
  cSee Section <ref>.
  
  
  
  The rates of the weak interactions that convert neutrons to protons, and vice versa, n ↔ p, have recently been discussed and computed by <cit.>, who took into account the effects of radiative corrections (including the effects of finite nucleon mass, finite temperature, weak magnetism, quantum electrodynamics, and incomplete neutrino decoupling). The weak rates scale inversely with the free-neutron lifetime[This value is based on the results of seven experiments compiled by the Particle Data Group <cit.> and three more recent measurements <cit.>. A slightly different value of τ = 879.5±0.8 s was used in <cit.>.], for which we adopt a value of τ = 879.7±0.8 s, and, thus, have an uncertainty of ≤ 0.1%. 
  
  
  The p(n,γ)d reaction rate has been calculated precisely using effective field theory <cit.>. The model parameters were constrained, using a Markov chain Monte Carlo procedure, by cross section data for neutron capture, deuteron photodissociation, and photon analyzing powers. The reported rate uncertainty is 0.4%.
  
  Ten reactions listed in Table <ref> have been measured directly in the laboratory at the energies of astrophysical interest. The rates for the ^3He(n,p)t, t(α,γ)^7Li, and ^7Li(p,α)α reactions are adopted from <cit.>. These were obtained using frequentist statistics and χ^2 optimization, with the implicit assumption of normal likelihoods and approximate treatments of systematic uncertainties. The reported rate uncertainties for these three reactions amount to 1.6%, 4.3%, and 2.5%, respectively. These  are small enough that they do not impact big bang nucleosynthesis predictions.
  
  
  The cross section data for the remaining seven reactions have been analyzed recently using Bayesian hierarchical models. For d(p,γ)^3He, d(d,n)^3He, d(d,p)t, and ^3He(α,γ)^7Be, microscopic theories of nuclear reactions provided the  cross section energy dependence, while the absolute cross section normalization was found from fitting the data within the Bayesian framework <cit.>. For ^3He(d,p)α, t(d,n)α, and ^7Be(n,p)^7Li, R-matrix theory was implemented into the Bayesian model to fit the data <cit.>. The Bayesian method provides statistically rigorous rate probability densities that can be used in Monte Carlo reaction network studies to derive meaningful simulated primordial abundance uncertainties. 
  
  
  
  
  <cit.> prefer to “use smooth polynomial fits to the data...by construction this ensures that the fits match the data as well as possible given the uncertainties.” In the present work, we instead prefer to use for our adopted rates (Table <ref>) well-established nuclear theory (microscopic models or R-matrix theory) as a foundation for the data analysis. We are not implying that one method is superior to the other, only that there is a range of reasonable options. We also prefer to adopt rates obtained using Bayesian techniques because they allow for a rigorous inclusion of statistical and systematic sources of uncertainties in the data analysis. For a recent example of such a study, see the analysis of ^7Be(n,p)^7Li data by <cit.>.
  
  A few other reaction rates that impact primordial nucleosynthesis, e.g., for d(n,γ)t, ^3He(t,d)α, and ^7Be(d,p)2α, have been adopted from the literature. They are listed in Table <ref> and will be discussed in the text below and in the Appendix. The rates of all other reactions in our network (see Section <ref>) have been adopted from statistical nuclear model calculations. See <cit.> for details. We assigned a factor of 10 uncertainty to the latter rates.
  
  
  lcll 
  Other nuclear reactions in primordial nucleosynthesis considered in the present work.a
  
  
  
  Process        2cUncertaintyb      Most recent reference 
  
                 Reportedc            Adoptedd              
  d(n,γ)te                   8.0%         factor 2            <cit.>   
   
     ^3He(t,d)αf              30%          factor 10           <cit.>           
   
    ^7Be(d,p)2αg             factor 3       factor 3            <cit.>   
  
  aThe symbols p, d, t, and α denote the nuclides ^1H, ^2H, ^3H, and ^4He, respectively.
  bRate uncertainty at 1 GK.
  cReported in original work.
  dAdopted here.
  eSee Section <ref>.
  fSee Section <ref>.
  gSee Section <ref>.
  
  
  
  
  
  § REACTION NETWORK, THERMODYNAMIC CONDITIONS, AND INITIAL ABUNDANCES
  
  
  
  We compute primordial nucleosynthesis using a reaction network consisting of 50 nuclides, ranging from p, n, ^4He, to ^22F. These are linked by 466 nuclear processes, such as particle captures, weak interactions, light-particle reactions, etc. Thermonuclear reaction rates are adopted  from STARLIB v6.8 (04/2020).[The STARLIB site has moved to <https://starlib.github.io/Rate-Library/>.] The library is in tabular format and lists reaction rates and rate factor uncertainties on a grid of temperatures between 1 MK and 10 GK <cit.>. Except for a few recently updated reaction rates (see text), the network is similar to the one used by <cit.>.
  
  
  Assuming that the baryonic density is known (Section <ref>), the thermodynamic conditions (temperature and density) in the early universe are determined by the expansion rate. The evolutions of temperature and density, adopted from <cit.>, are presented in Figure <ref>. We start the network calculations at a time of 1.0 s after the big bang, indicated by the vertical dashed lines in both panels. At that time, the temperature and density amount to T = 9.93 GK and ρ = 0.0537 g/cm^3, respectively. We end each calculation at a time of 10^5 s after the big bang. 
  
  
  
  
  The last ingredients needed for setting up the reaction network are the neutron and proton nucleon fractions at the beginning of the calculation. For times t ≤ 0.1 s, with temperatures of T ≥ 30 GK, all particles (photons, electrons, positrons, protons, neutrons, neutrinos, and antineutrinos) are in thermal equilibrium and the neutron-to-proton number abundance ratio is given by the Boltzmann distribution, N_n/N_p = exp(-Q/kT), where Q = 1.29333205(48) MeV is the energy equivalent of the neutron-proton mass difference <cit.> and k is the Boltzmann constant. As the universe expands and cools, neutrons and protons fall out of equilibrium because the weak rates (n + ν ↔ p + e^-, n ↔ p + e^- + ν, n + e^+ ↔ p + ν) 
  
  become comparable to, and then slower than, the expansion rate. Consequently, the neutron and proton abundances freeze out and their ratio continues to change as a result of free neutron decay. Figure <ref> displays the nucleon fractions of protons and neutrons as solid and long-dashed lines, respectively. The vertical dashed line indicates the time at which we start the network calculations, when the neutron and proton nucleon fractions are X_n = 0.23948 and X_p = 0.76052, respectively.
  
  
  
  
  Primordial nucleosynthesis starts about 100 s after the  big bang, when the temperature has declined to about 1 GK. This temperature is sufficiently low for the first nuclear reaction, p(n,γ)d, to become faster than its reverse photodisintegration, which until this point prevented the production of heavier nuclei. The nuclear transformations result mainly in the synthesis of ^4He, because it has the highest binding energy among all hydrogen and helium nuclides. It becomes the most abundant species (after ^1H), while all other nuclides are produced with much smaller abundances. The nucleosynthesis continues until about 1000 s after the big bang (shaded regions in Figure <ref>). For later times, the temperature and density in the early universe were too small for thermonuclear reactions to change primordial abundances. 
  
  
  
  
  § SIMULATED ABUNDANCES: MONTE CARLO STUDIES
   
  
  
  We will next investigate the answer to the following question: Given a set of recommended reaction rates and their uncertainties, and assuming that these were obtained by taking into account all known statistical and systematic effects in the measured data, what are the best estimates of the simulated primordial abundances? This question can best be answered using a Monte Carlo network technique <cit.>. The following steps are involved. First, the rates of all 466 reactions in the network are randomly sampled (see below). Second, the network is numerically solved using this set of sampled rates. Third, the first two steps are repeated n times to collect an ensemble of final abundance yields; n must be sufficiently large for statistical fluctuations to become much smaller than the widths of the obtained abundance distributions. Fourth, final primordial abundances are extracted from the accumulated abundance probability densities. 
  
  The nuclear interaction rates are sampled using lognormal rate probability densities <cit.>. For a given reaction, j, and temperature, T, the rate probability density is expressed as
  
  
      f[x(T)_j] = 1/σ√(2π)1/x(T)_j e^-[ln x(T)_j - μ(T)_j]^2/[2σ(T)_j^2]
  
  
  where the lognormal parameters μ and σ determine the location and the width of the distribution, respectively. For a lognormal probability density, rate samples, x_i, are drawn using <cit.> 
  
  
      x(T)_ij = x(T)_med,j [f.u.(T)]_j^p(T)_ij
  
  
  where x_med and f.u. are the median rate value and the rate factor uncertainty, respectively. Both of these are listed in columns 2 and 3 of STARLIB, respectively. The quantity p_ij is a normally distributed random variable, i.e., it is given by a Gaussian distribution with an expectation value of zero and a standard deviation of unity. Note that f.u.(T)^p(T), and not p(T), is the factor by which the sampled reaction rate is modified compared to its median value. We will refer to p(T) as the rate variation factor.
  
  We sample this rate variation factor for a given network run and nuclear reaction only once, i.e., p(T)_ij = p_ij has the same value at all temperatures. This assumption was found to reproduce the abundance uncertainties arising from more complex sampling schemes <cit.>. It implies that the rate samples still depend on temperature, since Equation (<ref>) takes the temperature dependence of the factor uncertainty, f.u.(T), into account. Also, the rates of corresponding forward and reverse reactions are not  sampled independently, because they are subject to the same value of the rate variation factor in a given network calculation.
  
  The Monte Carlo procedure has major advantages compared to varying rates one-by-one in sequential network runs. It is straightforward to estimate primordial abundances by adopting the 16th, 50th, and 84th percentiles of the resulting abundance probability densities. The impact of a given reaction rate uncertainty on the nucleosynthesis can be quantified by storing the values of p_ij for each sample reaction network run. A scatter plot for the final abundance of a given nuclide versus the sampled value of p_ij can then be inspected for correlations. Below, we will present results obtained by computing 10,000 Monte Carlo reaction network samples.
  
  Our simulated primordial abundances are listed in Table <ref> and displayed in Figure <ref>, together with those obtained using the code PRIMAT <cit.> and the most recently published measured values. 
  
  
  
  
  llll
  Observed and predicted primordial abundances.a
   
  
         Observedb      Presentc     PRIMAT18d   
  
  Y_p                                    0.2449(40)            e           0.24709(17) 
   
  (D/H)_p × 10^-5             2.527(30)             2.445(37)             2.459(36)   
   
  (^3He/H)_p × 10^-5        ≤1.1(2)          1.041(18)             1.074(26)   
   
  (^7Li/H)_p × 10^-10       1.58(31)              5.52(22)              5.62(25)    
   
  aFor ^4He, the nucleon fraction (Y_p) is given, while for the other species the number abundance ratio relative to ^1H is listed. Uncertainties are given in parenthesis, e.g., “0.2449(40)” stands for “0.2449±0.0040.”
  bFrom <cit.>.
  cRecommended values and uncertainties are obtained from the 16th, 50th, and 84th percentiles of the abundance probability densities (shaded blue in Figure <ref>).
  dFrom <cit.>.
  eNo improvement was obtained in the present work compared to the Y_p value of <cit.>, see text.
  
  
  
  
  The ^4He abundance depends sensitively on the rates of the weak interactions, n ↔ p, for which we adopted the theoretical results of <cit.> (see Table <ref>). The width of the abundance distribution, and, thus, the uncertainty of the predicted ^4He abundance (Table <ref>), depends mainly on the experimental uncertainty of the neutron lifetime since this value was used to calibrate the theoretical weak interaction rates. Since we adopted the weak interaction rates of <cit.>, we do not obtain an improved value for the ^4He nucleon fraction. Therefore, we do not list a present value of Y_p in Table <ref>.
  
  For ^2H (top panel in Figure <ref>), present and previous abundance predictions are in agreement. Interestingly, the measured value is larger by about 3% compared to the simulated ones. This “tension” had already been noted by <cit.>. The predicted deuterium abundance depends sensitively on the d(p,γ)^3He reaction rate for which we adopt the results of a recent Bayesian analysis (Table <ref>). New cross section data for this reaction are expected to be published soon by the LUNA collaboration[The recent experimental cross sections of <cit.> have not yet been used in d(p,γ)^3He rate evaluations. The uncertainties of their four data points are larger than previous results.]. 
  
  For ^3He (middle panel in Figure <ref>), the simulated abundances agree with the observed upper limit <cit.>. Our result is slightly lower than the simulation of <cit.>, which is explained by our adoption of the most recent ^3He(d,p)α rate <cit.> compared to their use of the older rate of <cit.>.
  
  Finally, the simulated present and previous ^7Li abundance values (bottom panel in Figure <ref>) are in agreement, despite the fact that both studies employ different ^7Be(n,p)^7Li reaction rates. The recent rate of <cit.>, used here, has a significantly larger uncertainty than the rate of <cit.> that was used by <cit.>. Regardless, the simulated ^7Li abundance exceeds the observed values by a factor of ≈ 3, as has been pointed out before <cit.>. We will consider correlations next.
  
  
  
  
  § CORRELATIONS BETWEEN ABUNDANCES AND REACTION RATES: MUTUAL INFORMATION
   
  
  
  We mentioned in Section <ref> that the network Monte Carlo procedure lends itself to studying correlations between simulated abundances and reaction rates. Since final abundances and rate variation factors are stored after each network run, the analysis of correlation scatter plots will reveal the impact of the variation of each rate on the abundance of every nuclide in the network. 
  The obvious challenge is to quantify which rates have the largest impact on a specific nuclide abundance. We need to adopt a useful metric to quantify the correlation, compute the metric for 50 isotopes × 466 reactions, and then sort the results according to impact. Pearson's (product-moment correlation coefficient) r, which was used by <cit.>, is a measure for the linear correlation between two random variables. <cit.> suggested to use Spearman's (rank-order correlation coefficient) r_s, which quantifies how well the relationship between two variables is described by a monotonic function. However, these metrics are not without problems in the present context, where correlations are frequently neither linear nor monotonic.
  
  Here, we will adopt the mutual information metric, which originates from information theory <cit.>. It quantifies how much information is communicated, on average, in one random variable about another when both are sampled simultaneously. For two random variables, Y and Z, with values of {y_1, y_2, y_3,...} and {z_1, z_2, z_3,...}, respectively, their mutual information is defined by
  
  
      MI = ∑_y ∑_z P(y,z) log[ P(y,z)/P(y)P(z)]
  
  
  where P(y) and P(z) are marginal distributions of y and z, respectively, and  P(y,z) is the joint probability density. An important theorem from information theory states that the mutual information between two variables is zero if, and only if, the two random variables are statistically independent. Unlike Pearson's r and Spearman's r_s coefficients, mutual information has no finite upper bound and, therefore, its absolute magnitude has no straightforward interpretation. However, here we are mainly interested in identifying the most important reactions for a given nuclide. For this purpose, the relative magnitude of the mutual information will suffice[<cit.> suggested to transform the mutual information by introducing the informational coefficient of correlation, defined as IC ≡ √(1-e^-2·MI). The coefficient IC lies between 0 and 1, equals zero when the two random variables are statistically independent, and equals unity when they are fully correlated. Again, there is no need in the present work to normalize the mutual information, MI.]. 
  
  
  
  
   §.§ ^4He abundance
  
  
  
  The weak interactions, n ↔ p, have the largest impact on the primordial ^4He abundance (Section <ref>). The correlation between Y_p and the rate variation factor for n → p is presented in Figure <ref>. The mutual information amounts to a value of MI = 1.1. None of the other reactions in the network are noticeably correlated with Y_p, i.e., all other mutual information values are MI ≤ 0.02. As already noted in Section <ref>, the weak interaction rate uncertainty is given by the present uncertainty in the neutron decay constant.
  
  Interestingly, if one would assume a factor of 10 uncertainty for the d(n,γ)t rate, the ^4He abundance distribution would display a pronounced tail on the right side of the peak. The published rate uncertainty <cit.> for this reaction is only 8%, although it is not clear to us how this value was derived. In the present work, we adopt a conservative uncertainty of a factor of 2. With this assumption, the d(n,γ)t reaction rate uncertainty has only a negligible impact on the primordial ^4He abundance. More information about the experiment of <cit.> is given in Appendix <ref>.
  
  
  
  
  
  
   §.§ ^2H abundance
  
  
  
  The primordial deuterium abundance is mainly impacted by the rate uncertainties of the d(p,γ)^3He, d(d,n)^3He, and d(d,p)t reactions. The correlations are depicted in Figure <ref>. The d(p,γ)^3He reaction is by far the most important among the three processes, as can be seen from its much larger mutual information value (MI = 0.59) compared to the d(d,n)^3He (MI = 0.09) and d(d,p)t (MI = 0.07) reactions. The rate uncertainties are listed in Table <ref> and have been derived using Bayesian techniques for all three reactions. 
  
  As expected, the d(p,γ)^3He rate is negatively correlated with the deuterium abundance: a smaller rate will result in a larger abundance of surviving deuterium nuclei. The present rate uncertainty for the d(p,γ)^3He reaction is 3.7%. Decreasing the recommended rate by 9% would shift the centroid of the simulated deuterium abundance distribution into the center of the observed range (top panel in Figure <ref>). Clearly, a new and improved measurement of the d(p,γ)^3He reaction is of significant interest (see Section <ref>). 
  
  Figure <ref> also displays the correlation for the ^3He(t,d)α reaction (MI = 0.04). The vast majority of network samples reveals an uncorrelated ^2H abundance, as can be seen from a distribution of points that is symmetric about a rate variation factor of p = 0. However, a few samples, for large values of p, do reveal a correlation. For this reaction, we adopted a rate uncertainty factor of f.u. = 10. An arbitrary increase of the recommended rate by a factor of 400 would also shift the peak of the simulated deuterium abundance distribution into the center of the observed range. The status of the rate of this reaction is discussed in Appendix <ref>. 
  
  
  
  
  
  
   §.§ ^3He abundance
  
  
  
  The primordial ^3He abundance is mainly influenced by the two reactions that produce and destroy ^3He, d(p,γ)^3He and ^3He(d,p)α, respectively. Correlations can be found in Figure <ref>. The mutual information values are MI = 0.58 and MI = 0.14, respectively, for these two reactions. The reaction rate uncertainties are listed in Table <ref> and have been derived using Bayesian techniques. None of the other reactions in the network are noticeably correlated with the ^3He abundance, i.e., all other mutual information values are MI ≤ 0.02.
  
  
  
  
  
  
   §.§ ^7Li abundance
  
  
  
  Of the primordial lithium synthesized in the early universe, about 95% originated from radioactive ^7Be during primordial nucleosynthesis, which later decayed to ^7Li when the expansion of the universe gave rise to temperatures and densities insufficient to sustain nuclear reactions. The final ^7Li abundance is mainly impacted by the rate uncertainties of the ^3He(α,γ)^7Be, d(p,γ)^3He, ^7Be(d,p)2α, and ^7Be(n,p)^7Li reactions. The correlations are provided in Figure <ref>. The mutual information values for these four reactions amount to MI = 0.20, 0.19, 0.11, and 0.07, respectively. The mutual information values of all other reactions in the network are MI ≤ 0.03.
  
  Recently estimated Bayesian rates are available for the ^3He(α,γ)^7Be, d(p,γ)^3He, and ^7Be(n,p)^7Li reactions (see Table <ref>). A more recent ^3He(α,γ)^7Be rate has been published by <cit.>, which is based on the calculation of <cit.>. However, these results are highly problematic and should not be used in big bang nucleosynthesis simulations, for reasons explained in Appendix <ref>. The status of the d(p,γ)^3He reaction rate has already been mentioned in connection with the primordial deuterium abundance (Section <ref>). The uncertainty of the most recently published ^7Be(n,p)^7Li rate <cit.> is 2.1% at big bang nucleosynthesis temperatures (Table <ref>).
  
  Figure <ref> indicates a noticeable tail towards smaller primordial ^7Li abundances for large values of the ^7Be(d,p)2α reaction rate variation factor, p. In our simulations, we adopted the recently published rate of <cit.>, which has a reported uncertainty of a factor of ≈3 at big bang nucleosynthesis temperatures. Increasing the recommended rate by a factor of ≳200 would lead to a significant overlap between the predicted and measured abundance values, and, thus, would solve the cosmological lithium problem. Considering the large uncertainties involved in the experiment of <cit.> (30% in the absolute cross section normalization and ±45 keV in the experimental energy calibration), it is important to remeasure this reaction in the future with improved techniques. More information regarding the present status of the ^7Be(d,p)2α reaction is provided in Appendix <ref>. 
  
  
  
  
  
  
  § UNKNOWN SYSTEMATIC BIAS IN REACTION RATES: GENETIC ALGORITHM
   
  
  
  In Section <ref>, we discussed how to estimate primordial abundances, given a set of thermonuclear rates for the key nuclear reactions (Table <ref>) that were derived from measured data, by taking into account all known statistical and systematic effects. In agreement with previous work, we find that the simulated primordial lithium abundance exceeds the measured value by a factor of ≈3 (Section <ref>). However, suppose that one or more nuclear reaction rates are subject to some unknown systematic bias, so that their true magnitudes differ significantly from the reported values. Therefore, in this section we ask a different question: By what factors would one need to multiply key reaction rates to reproduce simultaneously all observed primordial abundances? 
  
  The answer to this question is not simply obtained by independently changing the rates of individual nuclear reactions. We already discussed in Section <ref> that, for example, the d(p,γ)^3He reaction impacts the primordial abundances of ^2H, ^3He, and ^7Li (Figures <ref>, <ref>, and <ref>, respectively). Therefore, we must consider the simultaneous rate changes of several nuclear reactions.
  
  We will focus on the reactions d(p,γ)^3He, ^3He(α,γ)^7Be, ^7Be(d,p)2α, and ^7Be(n,p)^7Li, which correlate most strongly with the primordial ^7Li abundance (Figure <ref>). We require that simultaneous rate changes of these reactions reproduce the measured ^2H, ^3He, and ^7Li abundances (Table <ref>). Although <cit.> prefer to report an upper limit for the ^3He abundance inferred from their observations of distant metal-poor galactic H II regions, their upper limit value agrees with the simulated primordial ^3He abundance (Table <ref>). Therefore, for the discussion in this section, we will follow the suggestion of <cit.> and adopt the upper limit of <cit.> as an estimate of the primordial ^3He abundance.
  
  To investigate this optimization problem, we use a genetic algorithm, which provides a robust means of finding the global extremum in a multi-dimensional parameter space. Genetic algorithms are inspired by the mechanism of biological evolution through natural selection. See, for example, <cit.> for details. In brief, the idea is to start with an initial population of unbiased random sets of model parameter values. Each such set represents an individual in the overall population. Depending on their fitness, pairs of individuals are selected for the operations of crossover and mutation to produce offspring until the current generation of individuals has been replaced by a new generation of offsprings. The fitness of each offspring is assessed, and the process repeats itself by evolving through subsequent generations. As a result, the average population fitness increases over time, i.e., a better solution is produced with each generation. The process is terminated when a suitable solution (e.g., a parameter set of predefined fitness) is found or a given number of generations is reached. A major benefit of genetic algorithms is their ability to explore different parts of parameter space simultaneously <cit.>.  
  
  For the genetic algorithm-based optimizer, we used the general-purpose subroutine PIKAIA v1.2[The public domain source code and accompanying documentation is available at <https://www.hao.ucar.edu/modeling/pikaia/pikaia.php>.], which was first presented in <cit.>. It uses a stochastic rank-based selection scheme, a uniform one-point crossover (occurring with a probability of 0.85), and a uniform one-point adjustable mutation rate based on fitness. We set the initial, minimum, and maximum mutation probabilities to values of 0.005, 0.0005, and 0.25, respectively. For the reproduction plan, we chose the “steady-state-delete-worst” option, where the least fit individual of the parent population  is eliminated and replaced by the offspring.
  
  Fitness is the only point of contact between the genetic algorithm and the problem being solved. We will adopt a figure of merit of the form
  
  
      χ^2(u) = ∑_i=1^3 [y_obs,i - y_i(u)]^2/σ_obs,i^2
  
  
  where y_i and σ_i denote the mean values and uncertainties, respectively, of the ^2H, ^3He, and ^7Li measured abundances (Table <ref>), and y_i(u) are the corresponding predictions that depend on the vector of model parameters, u. In our context, the model parameters, u_1, u_2, u_3, and u_4, are the rate multiplication factors of the reactions d(p,γ)^3He, ^3He(α,γ)^7Be, ^7Be(d,p)2α, and ^7Be(n,p)^7Li, respectively. The fitness can then be defined as f = χ^-2. If the predictions for all three abundances, (D/H)_p, (^3He/H)_p, and (^7Li/H)_p, agree with their respective observations within uncertainties, then f ≳ 0.33 according to Equation (<ref>). 
  
  First, we restricted the search to the parameter regions 0.9 < u_1 ≤ 1.1, 0.9 < u_2 ≤ 1.1, 0.033 < u_3 ≤ 30, 0.5 < u_4 ≤ 2. These ranges significantly exceed the reported rate uncertainties listed in column 2 of Tables <ref> and <ref>. The result of the evolutionary run is shown in Figure <ref>. The top panel displays the value of f for the fittest individual in each generation. The system is seen to converge after ≈100 generations. After 1000 generations, the maximum fitness value found by the algorithm was only f = 0.11. In other words, changing the rates of these four reactions within the ranges specified above cannot reproduce the observational abundances. The bottom panel presents the corresponding evolution of the rate multiplication factors. While the factors for d(p,γ)^3He and ^3He(α,γ)^7Be are near unity (u_1 ≈ u_2 ≈ 1), those for ^7Be(d,p)2α and ^7Be(n,p)^7Li are at the upper values of their search ranges (u_3 ≈ 30, u_4 ≈ 2). Similar results are obtained with different random number seeds.
  
  Since the ^4He nucleon fraction depends mainly on the weak interaction rates, we have three data points (observations of ^2H, ^3He, and ^7Li) and four parameters (reaction rate multiplication factors). Because our system is underdetermined, we expect, in general, more than one possible combination of parameter values that could reproduce the observations, depending on the constraints (parameter ranges) applied. However, the important point is that the genetic algorithm does not find any acceptable solution at all if the rates of key reactions are varied simultaneously over ranges that are large compared to the reported uncertainties.
  
  Once the search ranges for ^7Be(d,p)2α and ^7Be(n,p)^7Li (u_3 and u_4, respectively) are widened even more, the genetic algorithm finds, within only a few generations, many different acceptable solutions (f ≳ 0.33), depending on the magnitude of the search ranges. For example, agreement between observed and simulated abundances is achieved if the ^7Be(d,p)2α and ^7Be(n,p)^7Li rates are simultaneously multiplied by factors of u_3 ≈ 0.2 and u_4 ≈ 3.4, respectively, or if the ^7Be(d,p)2α and ^7Be(n,p)^7Li rates alone are multiplied by factors of u_3 ≈ 410 and u_4 ≈ 4.0, respectively. Based on the current status of the nuclear physics data, we find such large rate changes to be unlikely. For more information about the ^7Be(d,p)2α reaction, see Appendix <ref>.
  
  
  
  
  
  
  § CONCLUDING SUMMARY
   
  
  This work focused on a discussion of thermonuclear reaction rates for big bang nucleosynthesis simulations. We started by summarizing the most recently published reaction rates. Among the twelve nuclear processes that are key to primordial element synthesis, the rates of seven reactions, d(p,γ)^3He, d(d,p)t, d(d,n)^3He, ^3He(d,p)α, t(d,n)α, ^3He(α,γ)^7Be, and ^7Be(n,p)^7Li, have been analyzed using Bayesian techniques (Table <ref>). Such methods provide results that are less biased when compared to previously applied techniques. Subsequently, we discussed other input to our reaction network, such as the evolution of temperature and density in the early universe, and initial abundances. We then presented simulated primordial abundances, (D/H)_p, (^3He/H)_p, and (^7Li/H)_p, including their uncertainties. These were obtained using a reaction network Monte Carlo method, i.e., by simultaneously sampling the rates of all reactions in our network according to their rate probability densities. In agreement with previous works, our simulated ^7Li abundance exceeds the measured primordial value by a factor of ≈3. Correlations between predicted abundances and reaction rates were analyzed using the metric of mutual information, which originates from information theory. The most important reactions for the synthesis of ^7Li are, in decreasing order, according to their mutual information values, d(p,γ)^3He, ^3He(α,γ)^7Be, ^7Be(d,p)2α, and ^7Be(n,p)^7Li. These conclusions are based on thermonuclear rates that were derived by taking into account all known statistical and systematic effects in the data analysis. Finally, for these four reactions, we account for unknown systematic bias in the data by using a genetic algorithm to search for rate multiplication factors that would reconcile all simulated and measured primordial abundances. No such solutions were found within reasonable search ranges, adding weight to previous findings that the solution of the cosmological lithium problem is unlikely to be found within the realm of nuclear physics.
  
  
  
  We would like to thank Rafael S. de Souza for his help with the genetic algorithm. We would also like to express our gratitude to Robert Janssens for providing valuable feedback. This work was supported in part by NASA under the Astrophysics Theory Program grant 14-ATP14-0007 and by U.S. DOE under contracts DE-FG02-97ER41041 (UNC) and DE-FG02-97ER41033 (TUNL).
  
  
  
  
  
  
  
  § COMMENTS ON SPECIFIC NUCLEAR REACTIONS
  
  
  
  
  
  
  
  
   §.§ d(n,γ)t
  
  
  
  The d(n,γ)t has been measured by <cit.> in the energy range important for primordial nucleosynthesis. Their three measured cross section data points, at neutron energies of 30.5 keV, 54.2 keV, and 521 keV, have uncertainties between 10% and 15%, including counting statistics, detector response, correction factors, and absolute normalizations. Their data, including previous results of thermal-neutron capture measurements <cit.>, were fit using a theoretical calculation based on the Faddeev approach. However, <cit.> provide insufficient information to assess how this input translates to a reaction rate uncertainty of only 8%. For the present simulations, we adopted a more conservative rate factor uncertainty of f.u. = 2.  
  
  
  
  
   §.§ ^3He(t,d)α
  
  
  
  The ^3He(t,d)α recommended rate used in the present work was adopted from <cit.>. This rate is presumably the same as the one published in <cit.>, who quote as reference <cit.>. The latter work measured the ^3He + t cross section at triton energies in the range from 150 to 970 keV. Their result is in reasonable agreement with <cit.>, who measured the cross section at a higher triton energy of 1.9 MeV, but disagrees with the earlier measurement of <cit.>, who obtained a cross section higher by a factor of ≈3. <cit.> assigned a rate uncertainty of ±30% to the ^3He(t,d)α reaction at all temperatures below T ≤ 10 GK. Although the data are sparse and exhibit significant inconsistencies, it is unlikely that this rate could change by a factor of ≈400 required to impact the primordial deuterium abundance significantly, as discussed in Section <ref>. In the present work, we adopted a conservative uncertainty of a factor 10 for this rate.
  
  
  
  
   §.§ The ^3He(α,γ)^7Be rate
  
  
  
  For a long time, measurements of the ^3He(α,γ)^7Be reaction have been plagued by systematic effects. Depending on the detection geometry, either the prompt γ rays, the ^7Be decay, or the recoiling ^7Be nuclei were detected. However, the cross section data from different measurements displayed significant differences. Consequently, the reaction rate based on the available data was rather uncertain <cit.>. The situation at the time called for new, careful measurements. The S-factors measured subsequently  <cit.>, which are displayed in Figure <ref> as black data points, are in overall agreement. Also, some groups have applied all three different experimental methods (prompt, decay, and recoil detection), yielding consistent results. Compared to the older measurements (see, for example, Figure 3 of <cit.>), the modern data are systematically higher and have a significantly smaller dispersion. 
  
  
  These modern ^3He(α,γ)^7Be data were used in a number of studies to derive reaction rates <cit.>. The solid black line in Figure <ref> displays the fit based on a hierarchical Bayesian model <cit.>. It resulted in a reaction rate with an uncertainty of 2.4% in the temperature range of big bang nucleosynthesis (Table <ref>). This rate was adopted in the present work. 
  
  Unfortunately, recent publications <cit.> have led to confusing claims that cannot be substantiated, as will be explained below. The fit of <cit.> did not take into account the single data point of <cit.>, who inferred the ^3He(α,γ)^7Be S-factor at solar energy from the solar neutrino flux. This indirect determination is not on par with the direct measurements mentioned above because it is difficult to assess the systematic uncertainties involved. This data point, shown in green in Figure <ref>, reveals large uncertainties and is located far away from the energy range important to primordial nucleosynthesis. However, motivated by this result, a new theoretical S-factor for the ^3He(α,γ)^7Be reaction has been proposed by <cit.>, which is shown as the dashed blue line. It can be seen that the fit of <cit.> (and those of <cit.>) and the theoretical prediction of <cit.> both agree with the (green) data point at solar energy. However, importantly, the theoretical S-factor of <cit.> clearly disagrees with all of the modern data at the energies most relevant for big bang nucleosynthesis. In conclusion, the S-factor of <cit.> should not be used in primordial nucleosynthesis calculations. It underpredicts the data by about 15%, implying a reduction in the simulated primordial lithium abundance by a similar amount.   
  
  
  
  The reason why we have discussed the results of <cit.> in some detail is because their S-factor was adopted in the big bang nucleosynthesis simulations of <cit.>. The latter authors reported a smaller predicted lithium abundance, (^7Li/H)_p = (4.447±0.067)×10^-10, which appears to be in lesser disagreement with the observed value (Table <ref>). However, as explained above, their result is erroneous because it rests on an inappropriate S-factor extrapolation from the solar region to the energy range of primordial nucleosynthesis. 
  
  
  
  
  
  
  
   §.§ The ^7Be(d,p)2α rate
  
  
  
  The ^7Be(d,p)2α reaction was recognized for some time as the most promising one for solving the cosmological lithium problem <cit.>. Its rate required an increase by about two orders of magnitude compared to the rate of <cit.> to significantly reduce the ^7Li yield; see Figure 4 in <cit.> and Table I in <cit.>. The <cit.> rate was adopted from <cit.>, who assumed a constant S-factor of 100 MeVb based on the experimental data of <cit.>. <cit.> writes: “Experimental measurements of the differential cross-section for this reaction were made by Kavanagh (1960) for deuteron energies from 700 to 1700 keV. Lacking complete angular distributions, these data can be approximately converted to total cross sections by multiplying by 4π and (2) by multiplying by a factor of 3 to take into account contributions from higher excited states in ^8Be.” Unobserved resonances that could significantly increase the cross section were proposed by <cit.>. Until very recently, this possibility seemed unlikely based on measurements of the average cross section <cit.> and investigations of the properties of candidate resonances <cit.>. Notice that the experiment by <cit.> was only sensitive to the ^7Be(d,p)^8Be(2α) reaction channel. 
  
  In a recent experiment, <cit.> found a new resonance in the previously undetected ^7Be(d,α)^5Li(pα) channel. Their experimental S-factor is much higher compared to the results of <cit.> and <cit.>. This conclusion was supported by preliminary results of <cit.>, showing an increase of the S-factor around E_cm ≈ 0.3 MeV. 
  
  Although the experiment by <cit.> provided valuable data at temperatures relevant to primordial nucleosynthesis, the magnitude of the reaction rates in the temperature range of primordial nucleosynthesis changed only modestly compared to earlier estimates. <cit.> pointed out that the recommended rate of <cit.> is close to the result of <cit.>, who had introduced an approximate rate multiplication factor of 3. Also, the rate uncertainty reported in <cit.> is close to the factor of 3 uncertainty adopted by <cit.>. Figure <ref> displays the rates of <cit.> as blue lines and those adopted by <cit.> as red lines. The shaded region indicates the relevant temperature range. The rate difference is much smaller than the two orders of magnitude required for a significant change in the predicted primordial lithium abundance; see Table <ref> and <cit.>. We also note that the claim in <cit.> of a ≈5% reduction in the primordial lithium abundance when the ^7Be(d,p)2α rate is set equal to zero is most likely erroneous. In our network calculations, the reduction amounts to only 0.75%.
  
  It is, nevertheless, of interest to remeasure the ^7Be(d,p)2α reaction and to reduce the uncertainties in resonance energy and cross section, as suggested by <cit.>.
  
  
  
  
  aasjournal
  
  