Uncertainty quantification for proton-proton fusion in chiral effective field theory

We compute the $S$-factor of the proton-proton ($pp$) fusion reaction using chiral effective field theory ($\chi$EFT) up to next-to-next-to-leading order (NNLO) and perform a rigorous uncertainty analysis of the results. We quantify the uncertainties due to (i) the computational method used to compute the $pp$ cross section in momentum space, (ii) the statistical uncertainties in the low-energy coupling constants of $\chi$EFT, (iii) the systematic uncertainty due to the $\chi$EFT cutoff, and (iv) systematic variations in the database used to calibrate the nucleon-nucleon interaction. We also examine the robustness of the polynomial extrapolation procedure, which is commonly used to extract the threshold $S$-factor and its energy-derivatives. By performing a statistical analysis of the polynomial fit of the energy-dependent $S$-factor at several different energy intervals, we eliminate a systematic uncertainty that can arise from the choice of the fit interval in our calculations. In addition, we explore the statistical correlations between the $S$-factor and few-nucleon observables such as the binding energies and point-proton radii of $^{2,3}$H and $^3$He as well as the $D$-state probability and quadrupole moment of $^2$H, and the $\beta$-decay of $^{3}$H. We find that, with the state-of-the-art optimization of the nuclear Hamiltonian, the statistical uncertainty in the threshold $S$-factor cannot be reduced beyond 0.7%.


I. INTRODUCTION
In main sequence stars such as the Sun, the conversion of hydrogen to helium proceeds predominantly through the pp chain, which is primarily triggered by the weak pp-fusion process [1,2], The accurate determination of the rate of this reaction is a critical ingredient for our understanding of many stellar processes. The reaction rate is conventionally parametrized in terms of the S-factor, which is related to the cross section by where E is the center of mass energy and η = m p /E α/2 is the Sommerfeld parameter. This cross section can only be measured in experiments performed at rather high energies in order to overcome the Coulomb barrier between the protons. Extrapolations of S(E) to relevant energy domains, E < 10 keV, which is where the Gamow peak of the Sun lies, yield extremely large uncertainties. Therefore, we are forced to rely on theoretical calculations to provide a precise prediction. This situation turns the accompanying uncertainty analysis into an absolute necessity. However, the quantification of theoretical uncertainties, which can have many different origins, is a very difficult task and a rigorous uncertainty analysis is still lacking. It is the purpose of this Letter to significantly advance the state-of-the-art of theoretical uncertainty quantification for the pp-fusion process.
Although calculations of S(E) can be performed at first order in the weak coupling constant, they still require a detailed knowledge of the nuclear interaction, which has to be treated nonperturbatively. Potential models therefore provided the first insights into this process [3][4][5][6]. Still, obtaining reliable uncertainty estimates has always proven to be hard with such models for the nucleon-nucleon interaction. Effective field theories (EFTs) are systematic low-energy expansions in a small parameter and are capable to provide an estimate of the inherent systematic error. An effort to provide uncertainty estimates was first carried out in [7,8] using the so-called hybrid approach in which current operators were obtained from χEFT and wave functions from phenomenological potentials. The first complete χEFT calculation of the S-factor was carried out by Marcucci et al. [9], who obtained S = (4.030 ± 0.006) × 10 −23 MeV fm 2 , with an uncertainty (see below) seven times smaller than that of the previously recommended value [2] S = (4.01 ± 0.04) × 10 −23 MeV fm 2 .
Pionless EFT calculations [10][11][12][13], which use only the nucleons as explicit degrees of freedom, have obtained consistent values, albeit with slightly larger theory errors. Since uncertainty analysis was not the main objective of Ref. [9], their error estimates were based on simple assumptions. The current operators and the wavefunctions were calculated to high precision -O(Q 3 ) and O(Q 4 ), respectively -by employing χEFT interactions widely used in the literature. The uncertainty was then estimated from the range of S-factor values obtained by using two different short-distance regulators in the twobody current and the potential. Thus, the error reported in Ref. [9] only reflects the resulting spread in the S-factor at a mixed order in the EFT expansion.
In this Letter we study the pp-fusion process in χEFT at next-to-next-to-leading order (NNLO). In particular, we perform an accompanying uncertainty quantification that builds on recent progress in mathematical optimization and statistical analysis of chiral nuclear forces and ab initio nuclear theory [14][15][16][17][18]. We quantify both statistical uncertainties, associated with the determination of the relevant low-energy constants (LECs), as well as several systematic ones, related to the computational method and to the chiral expansion. In this process we aim for consistency, e.g. by regulating the weak-current operator of the pp-fusion process in a manner that is consistent with the β-decay used in the fit of the nuclear potential.
This paper is organized as follows. In Sec. II, we discuss the weak current operator and the initial and final state wave functions. In Sec. III, we present an analysis of the uncertainties in the S-factor calculation. In Sec. IV, we present our final results along with a discussion.

II. FORMALISM
The cross section for the reaction in Eq. (1) can be written in the center of mass frame as where p e,ν are the positron and neutrino momenta, E e,ν their energies, m d is the deuteron mass, v rel is the pp relative velocity, and q is the momentum of the recoiling deuteron. The function F (Z, E e ) accounts for the distortion of the positron wave function due to the Coulomb field of the deuteron. Its classical expression, which can be found in Ref. [19], is increased by 1.62% due to radiative corrections [20]. The summation runs over the spin projections of all the initial and the final state particles. The initial state |i and the final state |f are direct products of leptonic and nuclear states. At nuclear energies, the weak interaction Hamiltonian can be written in terms of the leptonic and the nuclear weak current operators aŝ where G V is the vector coupling constant. The matrix element of the leptonic weak current operator j µ between the leptonic wave functions is l µ e −iq·x , where l µ satisfies l σ l τ * = 8 p ν σ p e τ + p ν τ p e σ − g στ p ν ρ p eρ − i ρσδτ p ν ρ p eδ , with the summation running over the spin projections of the leptons. Note that Eqs. (5) and (7) look different from the ones in Ref. [9,21] because we include the 2E e and 2E ν denominators in Eq. (5) to make the phase space volume element manifestly Lorentz invariant. The nuclear wave functions are calculated up to NNLO in χEFT. The nuclear weak current operators are consistently derived from the same effective Lagrangian up to the same order in Q, i.e. O(Q 3 ). In this approach, the strong interaction dynamics between the nucleons involved is based on the same theoretical grounds as their coupling to the leptonic current.

A. Nuclear weak current operators
The current operators for charge-changing weak interactions were derived in χEFT in Refs. [7,8,22,23].
The one-body (1B) operators that give non-vanishing matrix elements between an S-wave pp wave function and the S-and D-wave configurations of the deuteron include, up to O(Q 3 ), the Gamow-Teller operator, at leading order and the "weak-magnetism" operator which is O(Q). Here the index l runs over both nucleons, g A is the axial vector coupling constant, τ − l is the isospin lowering operator (τ x l − iτ y l )/2. Formally, there are additional operators [8,24] that enter according to the χEFT power counting scheme 1 . The matrix elements of those operators are, however, kinematically suppressed for the extremely small proton energies being considered here. The χEFT 1B currents operators used here have the same structure as the ones obtained phenomenologically in earlier studies (see, e.g. Ref. [25]).
The expression for the axial-vector two-body (2B) current, which is O(Q 3 ), reads Here k = (p 2 − p 2 − p 1 + p 1 )/2 and p = (p 1 + p 1 − p 2 − p 2 )/4, where p l (p l ) is the nucleon momentum in the initial (final) state, τ − × = (τ 1 × τ 2 ) x − i(τ 1 × τ 2 ) y , σ × = σ 1 × σ 2 and F π is the pion decay constant. The constants d 1 and d 2 that appear in the short-range 2B current are constrained by the Pauli principle such that only one linear combination, d 1 + 2d 2 = c D /(g A Λ EFT ) enters. The constants c 3 and c 4 , which accompany the one-pion exchange current, also appear in the N N and πN interactions, as well as along with c D in the N N N interaction. For the current, it is customary to define the counterterm to parametrize the strength of the meson-exchange current.

B. Nuclear wave functions
The NNLO momentum-space potentials that we employ are non-local. We obtain the S-state (L = 0) and D-state (L = 2) components of the deuteron wave function in coordinate space by diagonalizing in a harmonicoscillator basis. We find that it is necessary to correctly reproduce the deuteron wave function beyond 50 fm in order to achieve infrared convergence of the radial integrals.
This bound-state problem is easily solved equally well in either momentum space or coordinate space. Unfortunately, this dual approach is not as trivial for the relative pp-scattering wave function, ψ(r; E) as the presence of the long-ranged Coulomb potential implies that the momentum space representation of the scattering potential becomes singular for equal incoming and outgoing relative momenta. In order to facilitate a numerical solution we follow the prescription of Vincent and Phatak [26] and introduce a cutoff radius R c . For radii r < R c , the Coulomb potential is well-defined in momentum space. In this region we can therefore apply ordinary methods to find ψ(r < R c ; E). Furthermore, if we choose R c > 35 fm, such that the short-ranged nuclear interaction becomes negligible, we can smoothly match the solution at r = R c to the asymptotic Coulomb wave function, where δ 0 is the S-wave phase shift with respect to the nuclear and Coulomb potential, F 0 and G 0 are, respectively, the regular and the irregular Coulomb wave functions, and χ 0 (r; E) is the radial S-wave pp wave function, which, upon ignoring higher partial wave contributions, is related to ψ(r; E) by The pp wave function obtained is correct to NNLO in χEFT and first order in the electromagnetic coupling constant, α. We do not include higher-order electromagnetic contributions explicitly in this work. These electromagnetic interaction terms mainly lower the central value of the pp-fusion cross section [6,7,9], as discussed in Sec. IV, and leave the corresponding error due to uncertainties in the description of the strong nuclear force unchanged.

C. Radial matrix elements
The earliest calculations of the S-factor were performed using the Gamow-Teller operator only. In this approximation, the S-factor is proportional to Λ 2 (E), the square of the overlap between the pp wave function at energy E and the deuteron wave function, given by where C 2 0 = 2πη/(e 2πη − 1). However, the S-factor we calculate includes deuteron recoil effects and mesonexchange current contributions. The former modifies not only the phase space but also the matrix element in Eq. (5). The meson-exchange current, which turns out to be the dominant correction, is conventionally quantified as δ 2B , defined as the ratio of the matrix element of J − 2B to that of J − GT [8]. While the dependence of Λ(E) on the LECs of χEFT is solely through the N N interaction in the initial and the final state wave functions, δ 2B depends explicitly on many of the LECs through the current operator, and therefore deserves a careful scrutiny in spite of its relatively small contribution to the overall value of the S-factor.

III. CALCULATIONS
In this section we compute the quantities that are relevant for describing low-energy pp fusion. In addition, we will explore several sources of statistical and systematic uncertainties.

A. Axial-vector coupling constant
The S-factor depends predominantly on g A explicitly, by the relation S ∝ g A 2 . The most recent world average value for g A , calculated by the Particle Data Group (PDG), is 1.2701(25) [27]. The values obtained in some of the recent experiments [28,29] are higher than the world average. For our uncertainty analysis, we use the PDG recommended value which in turn is very close to the 1.2695(29) value used by Marcucci et al. in Ref. [9].

B. Chiral interactions and currents
In this work we employ a large set of 42 different NNLO sim interactions [18] to describe the nuclear force. This family of potentials are derived from χEFT up to NNLO. However, each of them was constructed using one of the seven different regulator cutoffs Λ = [450, 475, 500, 525, 550, 575, 600] MeV, and their LECs were constrained using one of six different truncations of the maximum scattering energy in the world database of N N scattering cross sections. It is also particularly useful that the covariance matrix of the LECs have been precisely determined for each NNLOsim interaction. Therefore, utilizing this large set of systematically optimized interactions, with known statistical properties, allows us to better uncover systematic uncertainties and explore the statistical uncertainties and correlations in the pp-fusion process. It should be pointed out that an equally good description of the fit data is attained with all NNLO sim interactions, and the numerical value of all constrained LECs are of natural order. Thus, all NNLO sim interactions are expected to perform very well in the few-nucleon sector of nuclear physics. In detail, the strengths of all relevant LECs were simultaneously optimized to reproduce the selected N N − scattering cross sections, πN −scattering cross sections, the binding energies and charge radii of 2,3 H and 3 He, the quadrupole moment of 2 H, as well as the comparative β-decay halflife of 3 H. The systematic uncertainty stemming from the excluded higher-order contributions in the chiral expansion was also accounted for in the determination of the LECs. A detailed expose of the optimization protocol is given in Ref. [18].
Furthermore, and similarly to the work in Ref. [8], we regulate the ultraviolet behavior of the matrix elements of all current operators using a Gaussian regulator on the form exp[−k 2 /2Λ 2 EFT ], where k is defined in Sec. II A. The N N -sector of the NNLO sim potentials is regulated in similar way; exp[−(p/Λ EFT ) 2n ], where n = 3 and p is the relative momentum of the two interacting nucleons. We use the same value for the cutoff Λ EFT in the currents as we use in the interactions.

C. Central values
The family of 42 optimized χEFT interactions at NNLO is employed to construct an averaged central value for each one of the computed quantities using the arithmetic mean of the separate calculations. The magnitudes of the statistical uncertainties in each calculation are nearly identical. Thus it is not necessary to explore weighted average schemes.
The energy dependence of the astrophysical S-factor, is usually parametrized using a polynomial expansion A second-or third-order polynomial is most common.
Higher-order polynomials turn out to be ill-conditioned extrapolants attributed with large statistical uncertainties. The normalized n-th order derivative S n (0)/S(0) often serve as input to e.g. neutrino flux computations using solar models. We fit theoretical S(E)-values to a third-order polynomial across an energy range E = 1−30 keV. From this we extract the S-factor at zero energy.

D. Uncertainty analysis
We start by considering systematic uncertainties associated with the methods used to compute the cross section and to extrapolate to zero energy.
A polynomial fit of S(E) is uninformed by the underlying physics of the pp fusion process. The polynomial fit and subsequent extrapolation to zero energy will depend on the limits of the fit-interval I E = [E min , E max ], for which it is not evident a priori what limits to use. The numerical precision of our computational code allows us to safely set E min = 1 keV. To find the optimal E max , we construct a penalty function χ 2 that describes how well a polynomial function f fit describes the correponding χEFT predictions f calc in I E , and minimize this with respect to E max . Note that we have also used the shorthand f i to denote f (E i ), and evaluate energies in I E at 1 keV intervals. Overall, we find that a cubic polynomial fits the χEFT predictions better than a quadratic polynomial.
For a cubic polynomial, and by varying E max between 4-100 keV, we find that E max = 30 keV is optimal. For this choice, the statistical uncertainties of S(0), S (0)/S(0), S (0)/S(0), due to the polynomial fit, attain a minimum. The extrapolated S(0)-values, as a function of E max are shown in Fig. 1. The central value for S(0) decreases and the corresponding statistical uncertainty of the polynomial fit increases when E max 30 keV. The same trend is observed for all 42 NNLO interactions. The statistical error of the fit is ten times larger at E max = 100 keV than it is at E max = 30 keV. Therefore, we choose to fit all computed S(E)-values to a third order polynomial across the energy interval I E = [1,30] keV. The resulting statistical uncertainty for S(0) that comes from the cubic polynomial fit is 0.0002 × 10 −23 MeV fm 2 , which corresponds to a relative uncertainty of ∼ 0.005%. Obviously, the systematic uncertainty due to a possibly ill-determined E max could be much larger, typically ∼ 1%. In order to reliably extract the derivatives of the S-factor with respect to the energy, a sufficiently large fit interval has to be chosen. However, if the fit interval is made larger than approximately 30 keV the necessity of including P -waves in the incoming channel becomes apparent in Fig. 1. The statistical uncertainties, due to the polynomial fit, for the first and second logarithmic derivatives of S(0) are 0.02% and 1.18%, respectively. The polynomial fit uncertainties for remaining quantities are negligible.
In addition, we varied the Vincent-Phatak (VP) matching radius R cut between 5-50 fm to extract uncertainties related to our numerical approach to pp scattering.
We find that R cut = 35 fm provides robust solutions, and we conclude that the uncertainty in S(0) due to the Vincent-Phatak procedure is 0.002 × 10 −23 MeV fm 2 , which corresponds to a relative uncertainty of ∼ 0.05%.
The three remaining sources of uncertainties that we explore are related to the χEFT description of the ppfusion process. First, there exists a statistical uncertainty due to the non-zero variances of the LECs in the NNLO sim interactions. Second, the actual choice of input N N data is not uniquely defined. Indeed, it is not clear what maximum kinetic energy T Lab one should consider in the selection of scattering data. This ambiguity gives rise to a systematic uncertainty. Third, variations in the regulator cutoffs, Λ EFT , of the χEFT description of the interaction and the currents will induce another systematic uncertainty.
In the process of propagating statistical uncertainties we employ well-founded methods for error propagation (see Ref. [18] for details). We also calculate the corresponding statistical covariance matrix of the observables and quantities that we compute. To find the covariance Cov(A, B) between two observables, A and B, due to the statistical covariances Cov(α) of the LECsα, we employ the linear approximation where J A is the Jacobian vector of partial derivatives J A,i = ∂A ∂αi , and similarly for J B . We compute ppfusion in the S-wave channel only, thus the relevant LECs areα = (c 1 , c 3 , c 4 ,C pp 1 S0 ,C3 S1 , C1 S0 , C3 S1 , C3 S1− 3 D1 , c D ). Further, we define the uncertainty in any observable A, due to the uncertainties in the LECs, as We use the covariance matrices Cov(α) from Ref. [18], and obtain the necessary derivatives ∂A ∂αi via a straightforward univariate spline fit to 10 function evaluations of the observable A; changing the LECα i in the neighborhood of its optimal value. We benchmarked this spline approximation using the known derivative values of the deuteron binding energy [18], and found that it was accurate to at least ∼ 0.001%.
From this procedure we find that the error in S(0) due to the statistical uncertainties in the LECs is 0.009 × 10 −23 MeV fm 2 , which implies a relative uncertainty of ∼ 0.2%. This result is very stable for each of the simultaneously optimized NNLO interactions. Propagating also the statistical errors of g A increases this statistical uncertainty to 0.019 × 10 −23 MeV fm 2 . The statistical uncertainty of G V has a negligible impact on the statistical uncertainty of S(0). Also, the statistical uncertainties of all LECs have negligible impact on the logarithmic derivatives of the zero-energy S-factor. The relative uncertainties are; 0.05% in S (0)/S(0) and slightly more , 0.2%, in S (0)/S(0). Furthermore, the derivatives of S(0) are very insensitive to changes in the cutoffs Λ EFT and T max Lab . The relative uncertainty of the squared overlap of the deuteron and proton-proton wavefunctions, Λ 2 , is ∼ 0.05%, and this error is dominated by the LEC variations. The Λ 2 is only somewhat sensitive to variations in Λ EFT and T max Lab . For all 42 NNLO interaction that we employed, we only observed Λ 2 in the range 7.064−7.101. Not surprisingly, the corresponding uncertainty and variation in δ 2B , which is more sensitive to short distance physics, is much larger. The typical relative uncertainty due to variations in the LECs is ∼ 5%. Furthermore, δ 2B varies between 0.30% and 0.52% for all NNLO interactions we explored.
As previously mentioned, systematic uncertainties in the chiral expansion is probed using the family of 42 different simultaneously optimized NNLO potentials. The range of NNLO predictions for S(0), including the total statistical uncertainty, is shown in Fig. 2. For a given cutoff Λ EFT , the width of the green band indicates the magnitude of all considered uncertainties. The total error budget is dominated by the statistical uncertainties in the sub-leading LECs of the chiral expansion and the axial-vector coupling constant g A . We operate with currents and interactions that have been optimized simultaneously and at the same chiral order, which ensures consistent renormalization.

E. Correlation analysis
In addition to the diagonal variances, we also compute the statistical correlations between all relevant pp-fusion quantities and observables. This study includes masses, radii, and half-lifes of A = 2, 3, 4 nuclei. Correlations can possibly reveal more information, but this exercise also serves as a sanity check of the entire uncertainty analysis. We should recover correlations expected from physical arguments. We employ the Jacobian and covariance matrices of A = 2, 3, 4 observables with respect to the NNLO LECs published in Ref. [18] and contract those with the spline Jacobians extracted in this work. A graphical representation of the relevant correlations is shown in Fig. 3. This particular correlation matrix is based on the NNLO interaction with Λ EFT = 500 MeV and T max Lab = 290 MeV. The same pattern emerges with any of the 42 different interactions employed in this work. As expected from the Q-value dependence of the phase space volume, the Sfactor strongly anticorrelates with the deuteron ground state energy. It is noteworthy that the squared radial overlap Λ 2 of the deuteron and relative-proton wave functions does not correlate significantly with S(0). This indicates that the dependence of the S-factor on binding energy indeed occurs predominantly through the phase space. We also observe that an increase in the deuteron radius would increase the radial overlap with the protonproton wave function. The quadrupole moment of the deuteron and its D-state probability anti-correlate with Λ 2 . Here, it is important to point out that our squared radial overlap only contains the 1B piece of the current operator. Thus it only measures the overlap between S-wave components. A smaller D-state probability implies a larger S-state probability. Consequently, the anticorrelation between Λ 2 and Q( 2 H)/D( 2 H) mostly traces the same underlying S-wave component of the deuteron wave function. Finally, we observe a strong correlation between the strength of the 2B current and the reduced axial-vector current of the triton β-decay. In fact, the LEC c D plays a dominant role for both currents. In conclusion, we quantify all expected correlations and confirm that they emerge in our statistical analysis.

IV. RESULTS AND DISCUSSION
We have calculated the pp-fusion S-factor using χEFT and carried out a state-of-the-art uncertainty analysis by employing a family of mathematically optimized chiral potentials at NNLO with consistently renormalized currents. We focused on the threshold S-factor and have therefore only considered initial S-wave pp scattering. To O(α), we obtain a threshold S-factor where we combined, for simplicity, all uncertainties by adding them in quadrature, and then taking the min/max values of the green band in Fig 2. This error represents all uncertainties originating from χEFT, the computational method, and the statistical extrapolation to obtain the threshold value. The effects of higher order electromagnetic contributions that are proportional to α 2 remains to be accounted for. These corrections lower the threshold S-factor by about a percent [6,7,9]. From the energy dependence of these corrections, calculated in Ref. [6], we estimate a 0.84% reduction in S(0).
For comparison, the uncertainty presented here is four times larger than the estimate reported in the pioneering χEFT calculation in Ref [9]. The comparison of the central values, however, is not so straightforward since their calculation includes additional terms in the current operator involving additional LECs, namely g 4S and g 4V , and the relativistic correction to the axial one-body current. We estimate the contribution of these missing terms to be of the order of 0.1% . However, our value given in Eq. (20) seems to be slightly higher than the result obtained by Marcucci et al. even when these corrections are considered. This issue could be related to the infrared convergence of the matrix elements between a bound state and a scattering state wave function, as mentioned briefly in Sec. II B. Additional details will be communicated through a separate publication [30].
Furthermore, our work displays that great care has to be taken in order to obtain reliable uncertainty estimates in EFT calculations. While it had previously been understood that a cutoff variation is one necessary part in the quantification of uncertainties, our analysis shows that potentials that use the same regulator but are optimized to a different energy range of scattering data can lead to different results for electroweak matrix elements. This is even more remarkable as the various potentials give the same low-energy observables. We also find that the choice of the fit interval is an important consideration when using polynomial extrapolation to find the threshold S-factor. The appropriate interval has to be chosen by comparing the statistical errors of the fit for several different fit intervals.
In this work, we have only considered S-wave initial state interaction. In the calculation of Ref. [9], capture in the P -wave channel contributed about 1% to S(E) at low E, which is roughly the same size as the errors stemming from our χEFT description of the interactions and from higher order electromagnetic processes. More generally, the effect of higher partial waves on a low-energy scattering process has the generic supression given by δ l ∝ k 2l+1 , which suggests that the P -wave contribution to the fusion cross section formally enters at O(Q 4 ), which is beyond the order we work to in this paper. It would be interesting to combine our type of analysis with the results obtained using pionless EFT. There, the S-factor is expressed as a low-energy expansion that contains only effective range parameters and a few lowenergy constants. Such a parameterization should be particularly useful in extracting reliable uncertainties for the energy dependence of the S-factor. Future work also involves to consider muon-capture by the deuteron. This two-nucleon process contains information on thed R -term that is strongly tied to the three-nucleon force. This study could provide additional important information to constrain the nuclear many-body Hamiltonian.