Gaussian process error modeling for chiral effective-field-theory calculations of $np\leftrightarrow d\gamma$ at low energies

We calculate the energy-dependent cross section of the $np\leftrightarrow d\gamma$ process in chiral effective field theory and apply state-of-the-art tools for quantification of theory uncertainty. We focus on the low-energy regime, where the magnetic dipole and the electric dipole transitions cross over, including the range relevant for big-bang nucleosynthesis. Working with the leading one- and two-body electromagnetic currents, we study the order-by-order convergence of this observable in the chiral expansion of the nuclear potential. We find that the Gaussian process error model describes the observed convergence very well, allowing us to present Bayesian credible intervals for the truncation error with correlations between the cross sections at different energies taken into account. We obtain a 1$\sigma$ estimate of about 0.2\% for the uncertainty from the truncation of the nuclear potential. This is an important step towards calculations with statistically interpretable uncertainties for astrophysical reactions involving light nuclei.

We calculate the energy-dependent cross section of the np ↔ dγ process in chiral effective field theory and apply state-of-the-art tools for quantification of theory uncertainty. We focus on the lowenergy regime, where the magnetic dipole and the electric dipole transitions cross over, including the range relevant for big-bang nucleosynthesis. Working with the leading one-and two-body electromagnetic currents, we study the order-by-order convergence of this observable in the chiral expansion of the nuclear potential. We find that the Gaussian process error model describes the observed convergence very well, allowing us to present Bayesian credible intervals for the truncation error with correlations between the cross sections at different energies taken into account. We obtain a 1σ estimate of about 0.2% for the uncertainty from the truncation of the nuclear potential. This is an important step towards calculations with statistically interpretable uncertainties for astrophysical reactions involving light nuclei.
Over the past couple of decades, effective field theories (EFTs) have emerged as powerful tools to study nuclear structure and reactions within systematic frameworks where a well-justified quantification of theory uncertainties is feasible [1]. An EFT is a model-independent description of interacting particles at some low-momentum scale Q in terms of the minimal degrees of freedom relevant below a breakdown momentum scale, Λ. Physical observables are calculated as systematic expansions in Q/Λ, with undetermined parameters called low-energy constants (LECs). The EFT can predict an observable at a given order once all the LECs that appear up that order have been fixed, e.g., by fitting to experimental data for other observables. The truncation of the EFT expansion at a finite order results in a theory error that can be estimated. Recent adoption of Bayesian methods, which have found increasing prominence in nuclear theory , have enabled a statistically meaningful estimation of truncation errors which can be easily combined with other uncertainties that arise from experimental data, calculational methods, and parameter fitting (see, e.g., Ref. [20]). Such a complete and consistent accounting of theory uncertainties is of paramount importance not only when the theoretical results have to be compared to precise experiments, but also when theoretical predictions need to be used as proxies for experimental data that are imprecise or unavailable. The latter situation is common for nuclear reactions in astrophysics and cosmology, for example in those relevant to big-bang nucleosynthesis (BBN).
There is, in general, a good agreement between simulations of nuclear reaction networks and astronomical observations for abundances of light nuclei that were produced during BBN. This agreement stands not only as a historic cornerstone of big-bang cosmology, but also as a crucial test for future extensions to the Standard Model. The ability of BBN to constrain key quantities such as cosmic baryon density rests heavily on the pri-mordial Deuterium to Hydrogen abundance ratio, whose uncertainty is currently dominated by the uncertainties of the Deuterium-burning processes: 2 H(p, γ) 3 He, 2 H(d, n) 3 He and 2 H(d, p) 3 H. Improvements in precision on their rates, exemplified by the recent 3% measurement of the 2 H(p, γ) 3 He cross section by the LUNA collaboration [26], will eventually call for an update of the primordial Deuterium production reaction, np → dγ. This reaction is the first step in the BBN network and marks the end of the so-called Deuterium bottleneck. BBN simulations [27][28][29][30] rely on theoretical predictions for this important cross section because available experimental data are too sparse and not sufficiently precise in the relevant energy regime. It is therefore imperative to study it using different theoretical approaches and to rigorously quantify uncertainties. Furthermore, ever increasing precision of the constraints on the primordial Deuterium abundance from cosmic microwave background analysis [31] and observational astronomy of low-metallicity stars [32] need to be independently matched by nuclear physics determinations.
This reaction has traditionally been studied using phenomenological models for nuclear interactions and currents [33][34][35]. Pionless EFT, the low-energy EFT of quantum chromodynamics (QCD) that uses nucleons as the only dynamical degrees of freedom and has a breakdown scale Λ at the pion mass m π , was applied to np → dγ at BBN energies in Refs. [36][37][38]. By fitting an LEC in the E1 amplitude to photodissociation data and an LEC in the M 1 amplitude to p(n, γ)d measurement [39] at energy E = 1.2625 × 10 −8 MeV in the np center-of-mass frame, Ref. [37] quoted a sub-percentage precision for the cross section at BBN energies. The pionless EFT results of Refs. [37,38], along with the accompanying uncertainties, are adopted by modern BBN simulations [28,30].
This cross section has also been calculated at E = 1.2625 × 10 −8 MeV in Ref. [40] using lattice QCD, in Ref. [41] using chiral EFT (χEFT) and in Ref. [42] using a hybrid approach that utilized χEFT currents and phenomenological potentials. Based on the scale hierarchy Q ∼ m π < Λ, χEFT provides a description of hadronic and nuclear phenomena in terms of nucleons and pions as effective degrees of freedom, and is therefore well suited for investigating this cross section at energies beyond the thermal capture, including the cosmological [20 keV-200 keV] range [43] and higher.
In this work, we present the first χEFT results for the energy-dependent np ↔ dγ cross section for a range of energies from threshold to MeV scale, that encompasses the BBN regime. We extend the Bayesian procedure for estimating the χEFT truncation error, which was earlier applied to various other observablesnucleon-mass expansion [2,5], nucleon-nucleon elastic scattering [4,7,8,13,19], nucleon-deuteron scattering [16,18], nuclear-matter equation of state [44,45], pion-photoproduction on the nucleon [21,25], and properties of light nuclei [20,22]-to an electromagnetic reaction cross section. We adopt the error model developed by Melendez et al. [13] that used machine learning with Gaussian processes (GPs) calibrated by physics-based hyperparameters to determine the convergence of EFT predictions that may be correlated across independent variables (e.g., energy). It allows us to quantify, for the first time, the error in the np → dγ cross section, σ np (E), and in that of the reverse process, σ γd (E), from truncation of the χEFT potential as a Bayesian degree-of-belief band.
We begin by presenting χEFT results for the np ↔ dγ cross section and related observables in Sec. I. We then briefly introduce the Gaussian process error model in Sec. II. The fitted model, diagnostic checks and predicted uncertainties are presented in Sec. III. We close with a brief summary and outlook in Sec. IV.

I. THE np ↔ dγ CROSS SECTION IN χEFT
The detailed-balance principle relates the np → dγ cross section σ np to the deuteron photodissociation dγ → np cross section σ γd : where m d is the deuteron mass, m is the isospin-averaged nucleon mass, and δm is the neutron-proton mass difference. The Mandelstam variable s can be conveniently expressed in terms of the np relative energy E, or the neutron energy E n in the rest frame of the proton, where m p is the proton mass. The photodissociation cross section can be expressed as a function of the photon energy ν (in the rest frame of the deuteron) through in terms of the transverse response function of the deuteron, defined as Here, j λ are the spherical components of the electromagnetic current operator j at four-momentum transfer (ν, q). The deuteron ground state is denoted by |ψ M , where M is the projection of the total angular momentum, whereas |φ k,S S z ,T denotes the pn scattering state with the relative momentum, total isospin, total spin and spin projection given by k, T , S and S z respectively. E ± are the energies of the final-state nucleons in the rest frame of the deuteron, i.e., E ± = (q/2 ± k) 2 + m 2 . The deuteron and the np scattering-state wave functions are obtained by solving their Lippmann-Schwinger equations using the momentum-space regularized semilocal potentials of Ref. [46] at various orders in the χEFT expansion. The electromagnetic currents were first derived in χEFT in Refs. [47][48][49]. We use multipole expansions [50,51] of these operators [62] that contribute at orders (Q/Λ) −3,−2,1,0 , namely the one-body currents consisting of the convection and spin-magnetization terms, and the leading one-pion exchange currents. Although we retain the quadrupole and octupole operators in the multipole expansions, their contributions are negligible compared to the dominant M 1 and E1 operators at low energies. To compare the relative strengths of the M 1 and E1 transitions, we calculate the analyzing power, defined as Σ(θ) ≡ (N − N ⊥ )/(N + N ⊥ ), where N (N ⊥ ) is the number of outgoing neutrons parallel (perpendicular) to the photon-polarization plane in a photodissociation experiment. It is related to the M 1 and E1 contributions to the photodissociation cross section through [38] where θ is the angle of the neutrons with respect to the photon beam axis. We now present the χEFT predictions for observables related to np ↔ dγ starting with the p(n, γ)d cross section for thermal neutrons in Table I. Here (and throughout this work), we denote cross sections calculated with an nth-order (N n LO) potential of Ref. [46] (and electromagnetic currents fixed to one-body plus two-body onepion exchange) as y n . The contribution y 1 is zero because there are no corrections to the χEFT potential at this order. At this threshold energy, our results undershoot the experimental value of Ref. [39] by a few percent. It was shown in Ref. [41], that percentage-level agreement with the experiment of Ref. [39] required current operators at the order (Q/Λ) 1 not considered in this work; these included leading two-pion exchange, sub-leading corrections to one-pion exchange and contact operators. Since the goal of this paper is to apply machine learning tools for the first time to uncertainty quantification in nuclear electroweak reactions in χEFT and including these operators significantly complicates the calculations, we work with fixed currents and focus on analyzing convergence in the χEFT potential only, leaving the inclusion of higher order currents to future work. We discuss the implications of our findings for such a study further below. In Fig. 1, we compare our χEFT predictions at different orders for various observables related to np ↔ dγ with the pionless EFT results of Ref. [37] and available experimental data at energies relevant for BBN and beyond. The χEFT predictions agree with the experimental data for all orders. Since Ref. [37] uses the experimental cross section given in Table I as input, their M 1 contribution is larger than our χEFT predictions at low energies.

II. GAUSSIAN PROCESS MODEL FOR CORRELATED EFT ERRORS
Following the formalism introduced by Melendez et al. [13], we consider an observable y as a function of the kinematic variable p, which, in our case is the np relative momentum. The order-k EFT prediction is written as and its EFT truncation error as Here y ref (p) is a dimensionful quantity that sets the overall scale such that the a priori unknown dimensionless coefficients c n (p) are smooth naturally-sized (i.e., of order 1) curves, provided that the EFT is converging in the expected manner. The Bayesian model for EFT error quantification seeks to build upon this prior knowledge with the data, the values of c n≤k (p) evaluated from order-by-order EFT calculations up through order k, to refine our expectations for c n>k (p), and thereby obtain an estimate for δy k (p). Refs. [3-5, 7, 8, 13] developed a pointwise error model that is applicable when (i) y is a number and not a function of one or more kinematic variables; (ii) y(p) is analyzed at only one value of p or at multiple values of p that are sufficiently far apart such that the y(p) values can be safely assumed to be uncorrelated. Ref. [13] extended this framework to study curvewise convergence by employing GPs. Below we first introduce the GP error model and then present application to np ↔ dγ. The GP error model-The basic idea of the error model is to use GPs to build a stochastic representation (an emulator ) to serve as a surrogate for the sequence of deterministic calculations (the simulator ) that yields order-by-order EFT predictions up through order k. The statistical properties of the emulator are then exploited to yield estimates for the EFT truncation errors at various orders, subsuming the contributions of even those terms that have never been calculated. Specifically, we start with the assumption that c n (p) are independent draws from an underlying GP, i.e., they follow a multivariate normal distribution for every finite set of p. The GP is completely specified by the mean µ and covariance functionc 2 r(p, p ; ), The correlation function r(p, p ; ) is commonly chosen to have a squared-exponential form, The correlation length , the mean µ and the marginal variancec 2 are the hyperparameters of the GP, which are learned from the training data set that comprises orderby-order EFT calculations up through order k. Remarkably, it follows from Eq. (8) that the truncation error defined by Eq. (7) has the distribution and With point estimates for and Q(p), the normal-inverseχ 2 prior serves as the conjugate prior for the hyperparameters (µ,c 2 ) of the Gaussian processes above, i.e., their Bayesian posteriors can be analytically derived and have the same functional forms as the priors (see Ref. [13]). The assumptions made above are known to impose certain limitations [57]; however their validity can be assessed from several diagnostic metrics on the validation data set (see below). The photon analyzing power for photodissociation versus its energy. Experimental data are from Refs. [52] (triangles), [53] (circle), [54] (crosses), [55] (square) and [56] (diamonds). Beam-energy resolution errors are not shown. Purple hexagons are calculated using pionless EFT results of Ref. [37].
Application to σ np and σ γd -We now apply the tools described above to the np ↔ dγ cross section. We begin by partitioning the data, i.e., the order-by-order χEFT results depicted in Fig. 1, into training and validation sets. This data set consists of 32 values of p in 1 MeV increment starting form the smallest that corresponds to thermal-neutron p(n, γ)d. These values span a range of approximately 0 to 1 MeV (2 MeV) in E (E n ), and encompass the 4 MeV < p < 14 MeV interval most relevant for BBN. We use every fifth point for training and the rest for validation. We will see further below that our maximum a posteriori (MAP) value of vindicates this grid and training-validation splitting. Similar results are obtained for other reasonable choices. As in Ref. [13], we take Λ = 600 MeV and Q(p) = (p 8 + m 8 π )/(p 7 + m 7 π ) , which give an expansion parameter of approximately 0.23 with a very weak dependence on p. We take y ref (p) = y 0 (p). The coefficient c 5 (p) turns out to be rather unnaturally sized, particularly at smaller values of p. We therefore exclude it from the statistical analysis. This leaves c 2,3,4 (p) for GP modeling. Finally, we need a value for the nugget σ 2 n , the variance of the Gaussian white noise added to the data to stabilize matrix inversions during fitting. We find that σ 2 n > 10 −8 is required to avoid singularities and that σ 2 n > 10 −3 can introduce noises similar in size to our precision goal of 0.1-0.2% on the predicted cross sections. We therefore pick σ 2 n = 10 −5 as this value produces the best performing model on the validation set under the diagnostic criteria discussed below.

III. MODEL CALIBRATION, VALIDATION AND PREDICTION
We now present the results of the GP modeling, obtained using the package gsum [13]. Fig. 2(a) shows the coefficients c n (p) for the observable σ np , along with their GP emulators. The MAP value of is found to be MAP =10.4 MeV. Our choices of the grid and trainingvalidation splitting places the training points at 5 MeV intervals. This results in 3 training data points within one correlation length, which is optimal for GP modeling. Furthermore, our data set spans approximately 3 MAP which, as a rule of thumb, is the range beyond which the data points become uncorrelated. Fig. 2(a) provides a beautiful visual indication that the GP model has done an excellent job of emulating the actual χEFT calculations. However, detailed diagnostic checks are needed to quantitatively assess the adequacy of the model. To this end, we use the diagnostic metrics proposed by Ref. [57], and originally implemented in EFT error model by Ref. [13]. Fig. 2(b) shows the squared Mahalanobis distances, defined as where we have used the notation f for the vector of validation data points, m for the vector of means of the emulator at these points and K for its covariance matrix. D 2 MD is a generalization of the sum of squared residuals to the case of data points that are correlated across the independent variable. Values much larger than its reference distribution, a χ 2 , indicate conflict between the emulator and the simulator. Values much smaller than the reference, as we see in the case of c 2 to some extent, indicate that, given statistical fluctuations, there is an unusually close match between the emulator and the simulator.
We now look at a more informative metric, the pivoted Cholesky diagnostic D PC . For a "correct" emulator, it returns draws from a standard Gaussian at every index, which represents the validation points arranged such that the first element is the one with the largest predictive variance, the second element is the one with the largest variance conditioned on the first element, and so on. A group of unusually large or small D PC values across all indices indicates a misestimated variance whereas a group of unusually large or small D PC values in the latter part of the sequence indicates an inappropriate correlation structure. Overall, Fig. 2(c) shows that the points are distributed as expected, e.g., 4 out of the 72 points lie outside the −2 < D PC < 2 range. We also notice that there is a slight indication that the variance on c 2 (c 3 ) might have been somewhat overestimated (underestimated) by observing the spread of the corresponding data points.
The credible interval diagnostic involves constructing uncertainty bands at each order and checking whether it actually encompasses the correction that enters at the next order. The claimed (1 − α)100% credible intervals are then plotted against the percentage of validation data points found within the interval-emulators that output credible intervals containing too few data points compared to the reference distribution are overconfident and those that contain too many are underconfident. For uncorrelated data points, the reference distribution is a binomial. For correlated data points, the reference distribution is numerically estimated by sampling a large number of emulators from the underlying process. Fig. 2(d) shows that the model is performing as expected and that it is important to account for correlations while assigning truncation errors. Now that we have demonstrated that the coefficients c 2,3,4 can be appropriately described by a GP, we can now use Eq. (10) to compute the truncation errors. We list the np → dγ cross section values at order k = 4 along with their 1σ truncation errors at several energies in Table II. We note that these errors are different from, and vary much more smoothly with E, than naive estimates [58] based on multiplying the largest order-to-order shift with the appropriate power of the expansion parameter in a pointwise manner at each value of E. In Fig. 3(a), we plot σ np v n , for which we showed order-by-order results earlier in Fig. 1(a), versus E n . This quantity is proportional to the reaction rate in BBN. The bands represent 2σ truncation errors, i.e., 95% Bayesian credible intervals for σ np v n . In Fig. 3(b), we show these bands for σ γd and compare them with photodissociation data shown earlier in Fig. 1(b). Reassuringly, the experimental data as well as the highest-order theoretical calculation, y 5 , are compatible with the truncation error estimates. The uncertainties quoted in Table II and Fig. 3, which amount to about 0.2%, only include truncation error in the χEFT potential. The full theory uncertainty also comprises statistical error from fitting the LECs to experimental data. The framework we have adopted allows one to also incorporate fitting uncertainties on the LECs that appear up to the EFT order we have calculated [13]. Another missing source of uncertainty, which is important in the M 1-dominated regime, is the truncation of the current. Inclusion of the M 1 operator at order (Q/Λ) 1 is crucial for obtaining agreement with the experimental value for threshold capture given in Table I [41,59], although it introduces several new LECs and thus poses significant challenges for rigorous uncertainty analysis. The fitting strategy that uses minimal assumptions about the short-distance behavior of the current operator, among several explored by Ref. [41], is the one that constrains the LECs d V 1,2 simultaneously to σ np (E = 1.2625 × 10 −08 MeV) and the isovector combination of the A = 3 magnetic moments. However, it was found that this yields unnatural values for d V 1,2 . This fine-tuning can be mitigated by including the theory uncertainty we calculated in this paper, as well as the experimental error, in the fit. Such a strategy for performing parameter estimation with χEFT truncation error included as a guard against overfitting was recently successfully pursued by Ref. [20] in the context of constraining 3N interactions from properties of A = 3, 4 nuclei. A calculation of np ↔ dγ along these lines is a subject for future work.

IV. SUMMARY AND OUTLOOK
We performed the first χEFT calculations of the energy-dependent np ↔ dγ cross section at low energies, including the range relevant to BBN, and the first Bayesian analysis of χEFT truncation error to a nuclear reaction cross section. Working with fixed one-and twobody electromagnetic current operators, we studied the convergence of this observable in the EFT expansion of the potential. By harnessing recent progress in Bayesian analysis of EFT uncertainty, we were able to provide statistical estimates, amounting to 0.2%, for the theory uncertainty that stems from truncation of the χEFT potential.
At the χEFT order up to which we work, our calculations are pure predictions as no new LECs enter that need to be fixed. Inclusion of the next order in the current operator adds new parameters that are not well constrained at present and will most likely require fitting to electromagnetic observables in A = 2, 3 systems. To make predictions for this cross section with subpercentage-level precision even at threshold, we will therefore need further investigation into the nature of the current operator at short distances and calculation of the electromagnetic observables for A = 3 nuclei with NN and 3N interactions truncated consistently. Finally, uncertainty analysis in theoretical calculations of Deuterium burning processes will be an important development because these strongly affect the primordial Deuterium abundance and there is currently some discrepancy between theory and experiments, most notably for 2 H(p, γ) 3 He [60].