Ab initio informed evaluation of the radiative capture of protons on $^7$Be

The radiative capture of protons by $^7$Be, which is the source of $^8$B that $\beta$-decays emitting the majority of solar neutrinos measured on earth, has not yet been measured at astrophysically relevant energies. The recommended value for its zero-energy S-factor, $S_{17}$(0) = 20.8$\pm$(0.7)exp$\pm$(1.4)theory eV$\cdot$b, relies on theoretical extrapolations from higher-energy measurements, a process that leads to significant uncertainty. We performed a set of first-principle (or, ab initio) calculations of the $^7$Be($p$, $\gamma$)$^8$B reaction to provide an independent prediction of the low-energy S-factor with quantified uncertainties. We demonstrate underlying features in the predicted S-factor allowing the combination of theoretical calculations and measurements to produce an evaluated S-factor of $S_{17}$(0) = 19.8$\pm$0.3 eV$\cdot$b. We expect the calculations and uncertainty quantification process described here to set a new standard for the evaluation of light-ion astrophysical reactions.

Astrophysical reactions powering low-mass stars such as our sun have been at the center of theoretical and experimental attention ever since nuclear reactions were proposed as a mechanism for nucleosynthesis and energy generation in stellar interiors [1,2].As a result, solar fusion reactions are amongst the most precisely measured and thoroughly evaluated nuclear reactions; see for example Refs.[3,4] and references therein.Occurring at the tail end of the proton-proton chain, the radiative capture of a proton by a 7 Be nucleus to produce an 8 B nucleus (or 7 Be(p, γ) 8 B reaction) is key in determining the solar neutrino flux measured in terrestrial observatories [5,6].Given its importance, it has been measured multiple times over the years with various techniques [7][8][9][10][11][12][13][14][15].However, due to Coulomb repulsion between the proton and the 7 Be nucleus, a direct measurement at the astrophysically relevant energies is still missing, and theory calculations [16,17] are used to extrapolate from higher energy experimental data.As a result of this extrapolation process, the uncertainty in the currently recommended [4] value of the zero-energy S-factor (S 17 (0)), S 17 (0) = 20.8 ± 0.7(expt) ± 1.4(theory) eV•barn, is dominated by theoretical contributions.
First-principle (or ab initio) theoretical approaches provide an independent prediction of nuclear reaction observables, with the interaction between nucleons being their sole input.Consequently, the bulk of the theoretical uncertainty of ab initio calculations will come from the nuclear interaction employed.In this Letter we present first-principle calculations of the 7 Be+p system, including the 7 Be(p, γ) 8 B reaction, using nucleon-nucleon (NN) and three-nucleon (3N) interactions derived from chiral effective field theory (χEFT) [18,19], with the goal of extracting universal features of the system, and removing (in part) the uncertainty that stems from the choice of a specific interaction and its parametrization.
The no-core shell model with continuum (NCSMC), first introduced in [20,21], is a first-principle technique that has been successful in delivering predictive calculations of nuclear properties of light nuclei by combining bound and dynamic descriptions of an A-nucleon system (see Ref. [22] for an in depth review of results).In the NCSMC, the A-body Schrödinger equation for a total angular momentum J and parity π is solved for bound and scattering boundary conditions by means of an ansatz that takes the form The states |AλJ π are obtained from the no-core shell model (NCSM), a many-body harmonic oscillator expansion method [23], and represent the λ-th bound-like solution to the A-body Schrödinger equation.The channel basis states A ν |Φ J π νr correspond to totally antisymmetric binary-cluster states where the interacting nuclei (in this case 7 Be and p) are a distance r apart.The collective index ν corresponds to all asymptotic quantum numbers (internal states, spins, and parities of the fragments, relative angular momentum , and spin s).The unknown discrete parameters c λ and continuous amplitudes γ J π ν (r) are then determined via the microscopic Rmatrix method [24].In the NCSMC, the A nucleons are treated as point-like, 'active', degrees of freedom and the same NN+3N interaction determines both the intrinsic wave functions (belonging to the projectile, target and aggregate system) obtained in the NCSM, as well as the reaction dynamics.
Nuclear interactions derived from χEFT come with a systematic organization in powers of a small expansion parameter Q < 1.The truncation of the expansion at some order gives rise to errors that can be quantified [19,[25][26][27][28][29][30], by combining consistent calculations performed at different truncation orders.In this work we use the NN interactions at next-to-next-to-next-to-leading (4th) order of the chiral expansion defined in Ref. [31], denoted N 3 LO * , and those at 3rd, 4th and 5th order of Ref. [32], denoted respectively N 2 LO, N 3 LO, and N 4 LO.These NN interactions are supplemented by the 3N interaction of Ref. [33] with both local (3N loc ) [34,35] and local plus non-local (3N lnl ) regulators [30,36].Finally, a 3N interaction employing a local plus non-local regulator but with an added sub-leading contact term enhancing the strength of the spin-orbit interaction (E 7 term), as described in [37], was also employed (3N * lnl ).A list of all NN+3N combinations that are employed in this work can be found in Table I.The interactions described in Ref. [32] are a consistent set and we expect truncation errors extracted from their differences to be similar to errors for other parametrizations.
In all calculations, the interaction was softened using a similarity renormalization group (SRG) transformation [38], including induced forces up to the three-body level.Four-and higher-body induced terms are small at the λ SRG =2.0 fm −1 resolution scale used in present calculations [39].The NCSM calculations of the 7 Be and 8 B nuclei were carried out allowing for up to nine quanta of excitation (N max = 9) for natural parity states and N max = 10 for the negative parity states of 8 B. The relative motion between the two fragments was also computed within the same N max = 9(10) space for positive (negative) parity channels.The harmonic oscillator parameter Ω was chosen at 20 MeV, which is close to the value that minimizes the ground-state energies of investigated nuclei [40].The reaction channel basis states are constructed by taking into account the first five states of 7 Be (J π = 3/2 − , 1/2 − , 7/2 − , 5/2 − , 5/2 − ).Earlier work [41] has demonstrated that this choice is sufficient to reach convergence in the channel basis expansion.The discrete part of the NCSMC ansatz is spanned by the ten lowest energy states of 8 B for each parity value, with spin values ranging from J = 0 to J = 4.
The positive parity eigenphase shifts for 7 Be+p scattering obtained using the N 4 LO+3N * lnl chiral interaction up to 10 MeV in center of mass energy (E c.m. ) show the well-established 1 + and 3 + 8 B resonances as well as several other, yet unobserved, broader resonances (Fig. 1).The NCSMC S-wave phase shifts manifest scattering length signs consistent with those determined in recent measurements (negative for 5 S 2 , positive for 3 S 1 ) [42], unlike previous calculations [40] that disregarded the discrete portion of the ansatz of Eq. (1).Using this parametrization of the NN+3N interaction, the 2 + state in 8 B is slightly unbound.The very narrow nearthreshold 2 + resonance is not visible in the figure.We find it is difficult to produce a bound 8 B ground state, with the exception of the N 3 LO+3N loc interaction, which reproduces the ground state energy at the 30 keV level.Owing to the parity difference between the ground state of 8 B (2 + ) and that of 7 Be (3/2 − ), the low-energy cross section of the 7 Be(p, γ) 8 B radiative capture proceeds via an E1 transition [43].While the bare E1 operator has a one-body form, a consistent SRG evolution to the same scale as the interaction, will induce many-body contributions.These contributions were found to be of a shortrange nature [44], thus we expect them to be small in a calculation involving a loosely-bound system such as 8 B. Operators of the M1/E2 types contribute at higher energies and are treated using a closure approximation of the 8 B NCSM states and 7 Be+p channels [22].
The 7 Be(p, γ) 8 B S-factor is sensitive to the spatial extent of the 8 B wave function [45], that is in turn determined by the ground state binding energy of 136 (1) keV [46].While it is impossible to exactly reproduce this binding energy in all NN+3N interaction models used in this study, especially when considering all sources of uncertainty, it is possible to adjust the NCSMC ansatz so that the overall binding energy reproduces the experimental value, as done for example in [47].This phenomenological correction, dubbed NCSMC-pheno, is performed by shifting the NCSM eigenenergies of 7 Be so that the excitation energies (and therefore thresholds) match the experimental ones exactly.This ensures that decaying states have the correct phase space available, corresponding to their energy.Furthermore, the 8 B NCSM eigenenergies in the 2 + , 1 + , 3 + channels are also modified bringing the bound and unbound NCSMC states in the experimentally observed positions.The resulting scattering lengths, asymptotic normalization coefficients (ANCs), and S 17 (0) are shown in Table I.We consider this approach to be an ab initio guided evaluation process where experimental data are fed into the theoretical prediction to correct small deficiencies of the nuclear interaction, resulting in greater predictive capability.
The astrophysical S-factors obtained within this NCSMC-pheno approach starting from the N  mentally observed sharp peak at about 0.6 MeV and, to a lesser extent, the broader structure around 2.2 MeV that are caused by magnetic dipole (M1) and electric quadrupole (E2) transitions of, respectively, the 1 + and 3 + resonances to the 2 + ground state of 8 B (Fig. 2).On the one hand, the calculation using the N 4 LO+3N * lnl interaction matches well the direct measurement from Junghans et al. [48] starting at the 1 + resonance and up until the 3 + bump, but slightly underestimates them at low energies.On the other hand, the N 3 LO * +3N lnl results follow the Junghans et al. data from threshold and up to the 1 + peak, but overestimate them somewhat at higher energies.More in general, none of the calculations we carried out is able to describe the S-factor experimental data over the entire range of energies.However, The multiple calculations of the 7 Be + p system allow for a more systematic look at its inherent properties without focusing on a specific interaction.Indeed, the use of various chiral order truncations and differing regulators gives us a window into the universal properties of the system as described by χEFT.As previously pointed out in Refs.[49,50], there is a linear relationship between S 17 (0) and the sum of the squares of the ANCs (C 2 p 1/2 + C 2 p 3/2 ).We observe such a relationship (Fig. 3a), with all NCSMC calculations lying along a line with slope 38.53 ± 1.45 eV•b•fm.The error bars in the points of Figs.3a,b, propagated through the linear regression, correspond to an estimate of the chiral truncation uncertainty for each interaction.We note that the error bars shown for the theoretical calculations in Fig. 3a should not be treated as uncorrelated errors; the relation between the ANCs and S 17 (0) is inherent in the equations being solved, not the specific interaction.The tight uncertainty band on the fitted line implies that if ANCs were accurately measured, the resulting theoretical uncertainty would be orders of magnitude smaller than the currently recommended value.Nevertheless, we can still look for correlations with al- ready available experimental data, that could lead to an accurate determination of S 17 (0).
We also find a nearly linear correlation between S 17 (0) and the S-factor at energies E c.m. < 0.5 MeV.As an example, in Fig 3b we show the correlation between the S-factor at zero and at 250 keV; the fit is obtained with a cubic polynomial incorporating deviations that are bound to exist when comparing different energies.We verified that the degree of the polynomial used does not significantly impact our results.While the 1σ uncertainty of the fit is slightly larger than in the linear fit of the ANCs, it is still an order of magnitude smaller than the recommended theoretical uncertainty in the region of the evaluation of Ref. [4].
We combine the tight relationship we have identified between the S-factor at zero and higher energies with experimental measurements for energies E c.m. ≤ 0.5 MeV to provide an evaluated prediction for S 17 (0) that includes both theoretical and experimental uncertainties.Specifically, we i) carry out a regression of the predicted S 17 (E c.m. ) versus S 17 (0) relationship for each E c.m. value where an experimental measurement exists (similar to Fig. 3b); and ii) infer the S 17 (0) (ordinate) corresponding to the measured S 17 (E c.m. ) (abscissa) from the fitted polynomials (Fig. 3c).The calculations at N 4 LO are essential in reducing the regression uncertainty, resulting in a small theoretical error bar for the inferred S 17 (0).We find the experimental contribution dominates the uncertainty budget, with the theoretical part being about an order of magnitude smaller.A regression to a horizontal line using all inferred data points yields S 17 (0) = 21.1 ± 0.2 eV•b (Fig. 3c), which agrees with the currently recommended value within 1σ.Extending this

FIG. 3. (a)
S-factor at zero energy versus the sum of the squares of ANCs for all interactions used in this work.Theoretical uncertainties are estimated using Q = 0.4.The linear fit (magenta line), along with its 1σ uncertainty band (dashed lines), suggests a reduced theoretical uncertainty compared to the current evaluation (grey shaded area, dark grey band for experimental uncertainty only).(b) Correlation between the S-factor at 0 keV and at 250 keV.The fitted line is a cubic polynomial (see text for details).The 1σ uncertainty of the fit is significantly smaller at the region of interest compared to the edges.(c) Resulting values for S17(0) using the regression method described in the text for each data point with Ec.m. ≤ 0.5 MeV.Error bars correspond to 1σ (theoretical plus experimental) uncertainty, while the cap points denote the experimental uncertainty only.The black line corresponds to a best fit value for S17(0) with a 3σ uncertainty band.The red line uses a modified fit measure (see text for details).(d) Same as (c) but using also data with 0.8 ≤ Ec.m. ≤ 1.5 MeV.
protocol to include energies with 0.8 MeV ≤ E c.m. ≤ 1.5 MeV brings the evaluation down to 20.4 ± 0.1 eV•b (Fig. 3d).However, we find a discrepancy between the S 17 (0) values predicted with the lower-and higher-energy data from Junghans et al. [14] (green diamonds).This discrepancy can be further illustrated by using only the higher-energy data for the fit, in which case we obtain S 17 (0)=20.2±0.1 eV•b.
The discrepancy in the values of S 17 (0) extracted from the three regressions described above can be traced back to our treatment of all measurements in an experiment as uncorrelated.This treatment leads to fits dominated by a single experiment with multiple measurements and small error bars, as the one of Ref. [14].This no correlation assumption is not necessarily true when looking at, for example, data points at nearby energies.Indeed, computing the covariance matrix for the theoretical calculations at various energies, we find large correlation coefficients (>0.9) across the full energy range in question.While, there are a few options (from guessing a covariance form, to using theoretical covariances), directly addressing this issue goes beyond the scope of the present work.Instead, we assign a factor normalizing the contribution of each experiment based on the number of data points.Such a normalization is equivalent to a assuming a fully correlated covariance matrix, where only one 'average' data point is used from each experiment.Following this process to fit to the E c.m. ≤ 0.5 MeV data yields S 17 (0) = 19.8 ± 0.2 eV•b, which again agrees with the recommended value within uncertainties, but disagrees with the value obtained from the fit that did not consider covariances at the 3σ level (Fig. 3c).The fit to the 0.8 ≤ E c.m. ≤ 1.5 MeV data gives a value of 19.8±0.2 eV•b, while fitting over both energy ranges yields 19.8±0.1 eV•b (Fig. 3d).Based on the robustness of the fit with modified experimental weightswe arrive at a suggested value for the 7 Be(p, γ) 8 B S-factor at zero energy of 19.8±0.3 eV•b, with the 1σ uncertainty conservatively increased to account for the different data fitting procedures.
In conclusion, we have conducted an extensive ab initio analysis of the 7 Be+p system, providing theoretical predictions for phase shifts and the radiative capture cross section.To this end, we performed NCSMC calculations with a variety of χEFT interactions, all of which provide a reasonable description of the 7 Be and 8 B nuclei, and further refined them to reproduce observed binding energies and resonance positions.Finally, we introduced a new protocol that combines predictive ab initio calculations (with quantified uncertainties) and experimental data to evaluate observables in regions where experimental measurements are not feasible, and demonstrated the necessity for high-order (larger or equal to fifth order) calculations with reduced chiral truncation uncertainties.Alongside such calculations, high-precision experiments both at low and higher energies can help in further pinning down the value of S 17 (0).We believe that this work sets a new standard for the evaluation of thermonuclear reaction rates that will lead to more robust predictive capabilities for astrophysical models.It also highlights the necessity of incorporating estimates of experimental and theoretical covariances in the evaluation process.

FIG. 2 .
FIG.2.Astrophysical S-factor of the 7 Be(p, γ) 8 B radiative capture obtained from the NCSMC-pheno approach with the N 3 LO * +3N lnl (blue line) and the N 4 LO+3N * lnl (red line) interactions compared to experimental data.The inset shows the low-energy part and shows the evaluation of S17(0) with its uncertainty[4] (orange box).
Computing support for this work came from the LLNL institutional Computing Grand Challenge program and from an INCITE Award on the Summit supercomputer of the Oak Ridge Leadership Computing Facility (OLCF) at ORNL.Prepared in part by LLNL under Contract DE-AC52-07NA27344.This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under Work Proposal No. SCW0498, LLNL LDRD Project No. 20-LW-046, under the FRIB Theory Alliance award no.DE-SC0013617, and by the NSERC Grants SAPIN-2016-00033 and SAPPJ-2019-00039.TRIUMF receives federal funding via a contribution agreement with the National Research Council of Canada.

TABLE I .
[42]es for asymptotic normalization coefficients (ANCs) in fm −1/2 , scattering lengths in fm, and zero-energy S-factor in eV•b obtained from the set of interactions used in this work, after applying a phenomenological correction (see text) to the 8 B bound-state as well as the 1 + and 3 + resonance energies.The experimentally extracted values from Ref.[42]are also included for comparison.