Posterior predictive distributions of neutron-deuteron cross sections

We quantify the posterior predictive distributions (PPDs) of elastic neutron-deuteron ($nd$) scattering cross sections using nucleon-nucleon ($NN$) interactions from chiral effective field theory ($\chi$EFT) up to and including next-to-next-to-next-to-leading order (N$^3$LO). These PPDs quantify the spread in $nd$ predictions due to the variability of the low-energy constants (LECs) inferred from $NN$ scattering data. We use the wave-packet continuum discretization method to solve the Alt-Grassberger-Sandhas form of the Faddeev equations for elastic scattering. We draw 100 samples from the PPDs of $nd$ cross sections up to 67 MeV in scattering energy, i.e., in the energy region where the effects of three-nucleon forces are expected to be small. We find that the uncertainty about $NN$ LECs inferred from $NN$ scattering data, when assuming uncorrelated errors, does not translate to significant uncertainty in the low-energy $nd$ continuum. Based on our estimates, the uncertainty of $nd$ predictions are dominated by the $\chi$EFT truncation error, at least below N$^3$LO. At this order, the 90% credible interval of the PPD and the truncation error are comparable, although both are very small on an absolute scale.


I. INTRODUCTION
Chiral effective field theory (χEFT) [1][2][3] promises a systematically improvable description of the nuclear interaction grounded in the symmetries of low-energy quantum chromodynamics. Two-nucleon (N N ) and three-nucleon (N N N ) interactions from χEFT are used extensively in modern ab initio predictions of atomic nuclei and nuclear matter, see, e.g., Refs. [4][5][6] for recent overviews. To make quantitative predictions of the properties of nuclear systems, the numerical values of the low-energy constants (LECs) that govern the strengths of the pion-nucleon (πN ) and nucleon-contact couplings must first be inferred from low-energy data. For this, the Bayesian approach to statistics [7] provides a natural framework since it yields a (posterior) probability density function (PDF) that quantifies our uncertainty about the values of the LECs. Propagating this uncertainty when making theoretical predictions amounts to averaging the distribution of predictive samples over the LEC posterior PDF. The result of this is called a posterior predictive distribution (PPD). This type of distribution sits at the center of the scientific process whereby we try to predict future data based on previous data and theory.
There are existing efforts to quantify Bayesian PPDs for various nuclear observables, e.g., N N scattering cross sections [8] and scattering lengths [9], few-nucleon [10] and many-nucleon [11][12][13] energies, radii, and decays, as well as nuclear mass models [14,15]. These probability distributions quantify our degree-of-belief, and facilitate a meaningful comparison with experimental data. For example, the PPD finds use in model checking [7], such as posterior predictive checks. There one simulates data, using a fitted model, and compares to observed data. The simulated data corresponds to draws from the PPD and it should look roughly like the observed data if the model did indeed contain all relevant physics and there has been a sufficient amount of calibration data.
In this work, we sample the PPDs of selected neutron-deuteron (nd) scattering cross sections arising from the variability of the LEC posterior when conditioned on N N scattering data. We use χEFT descriptions of the N N interaction at all orders up to next-to-next-to-next-toleading order (N 3 LO) in Weinberg power counting. To the best of our knowledge there exists only frequentist statistical analyses encompassing a subset of nucleondeuteron (N d) scattering cross sections and scattering lengths [16][17][18], for which various estimates for dispersion have been quantified. Our analysis is rooted in Bayesian methodology and therefore provides probability densities for the predicted observables of interest. As such, the results of this work facilitates a quantitative measure of the predictive power in the low-energy N N N continuum using χEFT interactions carefully calibrated using N N scattering data. This work is part of an ongoing effort towards a full Bayesian analysis of χEFT conditioned also on experimental data in the N d continuum [19].
To sample the PPDs of elastic nd-scattering cross sections, we repeatedly solve the Alt-Grassberger-Sandhas [20] (AGS) form of the Faddeev equations using the wave-packet continuum discretization (WPCD) method [19,21]. This method is parallelizable with respect to the scattering energy, denoted with E Lab , in the laboratory frame of reference. Therefore, it is particularly suitable for sampling PPDs across a range of E Lab values. Still, the collection of samples from the PPDs is limited by the number of times we can solve the AGS equation. For this reason, we currently neglect N N N forces (3NFs) and focus our analysis on cross sections and polarization observables with E Lab ≤ 67 MeV, for which N N -only models typically perform well [16,[22][23][24]. The low-energy vector analyzing power, A y , is a possible exception to this statement and we therefore place a special focus on the analysis of this polarization observable.
In addition to the inherent uncertainty of inferred LEC values, there are also other sources of theoretical uncertainty. The model discrepancy due to the omission of higher chiral orders is an obvious one. Neglecting this un-arXiv:2209.06501v1 [nucl-th] 14 Sep 2022 certainty can lead to biased and over-confident inferences and predictions [25]. Fortunately, χEFT is designed to be an order-by-order improvable description of the nuclear interaction, and as such the theory itself provides valuable information about the magnitude of the truncation error [26]. Indeed, there exists several efforts to quantify the truncation error in effective field theory predictions of nuclear systems, see, e.g., Refs. [8,10,[27][28][29][30][31]. Although our focus is to quantify the PPDs of nd scattering observables due to variability in the N N LECs, we will also contrast our findings with estimates of the truncation error.
In Sec. II we define the general structure of the PPDs we sample in this work. In Sec. III we present the essential elements of the WPCD method we use to produce elastic nd cross sections. In Sec. IV we present the strategy for sampling the PPDs, with particular focus on achieving computational speedup, and the results of the sampling. We also compare the degree-of-belief intervals of the PPDs with some of the other components of the total error budget; the χEFT truncation error in particular. We end with a summary and outlook in Sec. V.

II. SETTING UP THE POSTERIOR PREDICTIVE DISTRIBUTION
The PPD is a PDF pr(y|D, M, I) for a quantity y as predicted by a model M . This distribution quantifies the uncertainty about y given previous data D and any other assumptions or information I. Here, we focus on the uncertainty of the numerical values of the LECs, denoted α, present in the underlying χEFT N N interaction. As such, we must marginalize over the LECs by evaluating (1) In the second line we introduced a short-hand y( α) for a deterministic model prediction given numerical values for α from some parameter domain Ω. We also used that y is conditionally independent of D. The proportionality indicates that we are only interested in the width and shape of the PPD, and not the overall normalization constant.
We will refer to the χEFT description of the N N interaction at a chiral order ν as a "model", and denote this as M ν . The chiral orders are defined according to Weinberg power counting with ν = 0, 2, 3, 4, and as is common, we refer to them to as leading order (LO), next-to-leading order (NLO), next-to-next-to-leading order (N 2 LO), and N 3 LO, respectively. The values of α depend on the chiral order ν, but to simplify notation we do not index α by ν.
The PPD is a probabilistic generalization of the familiar point-estimate value y = y( α ), obtained by evaluating the model M ν at some preferred parameter value α , such as a local parameter-optimum. We will in some cases resort to evaluating the PPD at the maximum a posteriori (MAP) value of the LEC posterior α ≡ argmax α pr( α|D, M ν , I). ( Note that the PPD does not necessarily attain its maximum for α . Indeed, the evaluation of y( α), through the AGS equation, is neither linear nor monotonic. Evaluating the integral in Eq. (1) requires knowledge about the PDF, pr( α|D, M ν , I). We utilize the available LEC posteriors up to and including N 3 LO published in Ref. [9]. These posteriors were sampled using Hamiltonian Monte Carlo (HMC) while accounting for uncorrelated χEFT truncation errors, and were conditioned on the Granada database [32,33] of N N scattering cross sections for scattering energies E Lab ≤ 290 MeV. The leading neutron-neutron (nn) isospin-breaking LEC was inferred using an empirical value for the nn scattering length in the 1 S 0 partial-wave channel. We note that other methods accounting for correlated χEFT truncation errors exist, see, e.g., Ref. [34], which may change the inferred, and rather narrow, distributions of LEC values we use here.
The HMC algorithm is particularly well-suited for sampling high-dimensional PDFs and yields virtually uncorrelated draws from pr( α|D, M ν ). A detailed analysis [8] suggests that the HMC chains we employ in this work to represent the LEC posteriors are sufficiently converged at all orders, unimodal, and rather concentrated in parameter space. As such, we have in-depth knowledge of the location of the posterior mass, which helps tremendously when evaluating the integral in Eq. (1).

III. WAVE-PACKET CONTINUUM DISCRETIZATION
In this section we summarize the WPCD method [21] for solving the AGS equation in momentum space. Our results are based on the implementation presented in Ref. [19] 1 . The AGS equation for nd scattering, without 3NFs, can be written aŝ whereÛ is the transition matrix between asymptotic scattering states, is the resolvent of the free N N N Hamiltonianĥ 0 , E is the total energy, 1 The implementation, named "Tic-tac", is available under a GNU open-source license (GPLv3) on https://github.com/seanbsm/Tic-tac t 1 denotes the scattering T -matrix for the pair-system (23) as written in standard odd-man-out notation, and P ≡ 2P 123 whereP 123 is the permutation matrix acting on partially-antisymmetric N N N states 2 . The large dimensionality of the N N N Hilbert space makes it challenging to apply matrix-inversion type methods to solve Eq. (3). Instead, one usually resorts to expanding the AGS equation in a Neumann series that is subsequently re-summed using a Padé approximant [35] to handle the divergence originating from the integral-kernelĜ 0v1 with Weinberg eigenvalues [36] outside the unit circle.
It is well understood how to obtain converged solutions for U in a standard plane-wave basis, see, e.g., Ref. [37]. In this basis,Û is obtained for a specific value of the onshell energy E, and the resolventĜ 0 and N N T -matrix t 1 depend explicitly on E. This dependency inflicts several complications such as moving singularities in the resolvent operator, and a requirement for antisymmetrizing N N T -matrices at many energies when evaluating the AGS integral kernel, which is typically handled using splines [38].
In this work, we use the WPCD method [21] for solving the AGS equation. This is one of many boundstate approaches [39] for describing scattering processes. In WPCD, we discretize the continuum using a wavepacket basis. Doing so simplifies the numerical analysis of the AGS equation. First, one can derive a closed-form expression of the channel-resolvent, treating the associated singularities analytically. Second, the P -matrix has no need for splining. Third, it factorizes the on-shell energy dependence out of the matrix multiplications associated with the terms of the Neumann series expansion, providing significant speedup of the most time-consuming parts of the numerical solution.
As a downside, the WPCD method entails large matrix dimensionalities compared with the plane-wave representation. However, scattering amplitudes can be calculated at multiple scattering energies with minor extra computational cost per energy. This makes WPCD particularly suitable for sampling Bayesian PPDs across ranges of energies. In fact, we find that calculating scattering amplitudes at multiple scattering energies only doubles the computational cost compared to computing the amplitude at a single energy [19] We define a wave packet |x as a finite integral of continuum states |p , e.g., plane-wave states, within a mo- where f (p) is a weighting function and N is the normalization constant. An A-body wave packet can be straightforwardly defined using wave-packet discretization for each Jacobi coordinate. A N N N wave packet is given by |X ≡ |x ⊗ |x , where |x corresponds to the pair-system p-momentum and |x corresponds to the spectator q-momentum. The eigenstates of the N N Hamiltonianĥ 1 in a (planewave) wave-packet basis can be used to approximate "scattering" N N N wave packets rather well. In this basis, it is also possible to evaluate the channel-resolvent whereÛ now depends on E only viaĜ 1 . This is the starting point for solving the AGS equation in the WPCD method. Here, as in Ref. [19], we use an equal number of wave packets, N WP , to discretize the p and q continua, yielding matrices in Eq. (5) that scale in size as O(N 4 WP ). We find that the runtime of the code follows this quartic scaling with N WP quite closely. Note, however, that the calculations at N 3 LO are ∼ 10% more costly since the Padé resummation of the Neumann series typically requires more terms to converge.

IV. EVALUATING POSTERIOR PREDICTIVE DISTRIBUTIONS
We sample the PPD of a scattering observable by evaluating Eq. (1) numerically. This is done by computing the nd scattering observable of interest for a finite set of LEC values drawn from the posterior PDF, pr( α|D, M ν ). In practice, we use the Markov chains obtained in Ref. [9].
For every sample that we draw from the PPD we must solve the AGS equation. Fortunately, with the WPCD method we get access to all scattering cross sections at all angles and energies without any significant computational overhead. Also, since the permutation operatorP does not depend on the LECs, we only have to compute this once and re-use it throughout the sampling process. However, we have to setup the Neumann series for every new sample, and this is the most time-consuming part.
In all calculations done here, we use a spin-angular basis of N N N partially-antisymmetric partial-waves with total angular momentum J ≤ 17/2, using both parities, and using N N total angular momentum J ≤ 3. We also explicitly account for the charge dependence of the strong N N interaction in the 1 S 0 channel. This state space provides sufficiently converged U -matrix elements for E Lab ≤ 100 MeV when using the chiral potentials defined in Ref. [8,9]. Note that our study is limited to E Lab ≤ 67 MeV due to the omission of 3NFs. It has been shown that, at low scattering energies, the scattering amplitudes are likely dominated by N N forces [16,[22][23][24].
We discuss our general strategy to quantify the PPD in Sec. IV A, present results for the PPDs of the differ-ential nd cross section in Sec. IV B, relate this to estimates of the χEFT truncation errors in Sec. IV C, and discuss spin-polarization observables, focusing on A y (n), in Sec. IV D.
A. Trading wave-packets for computational speedup In the limit N WP → ∞, the WPCD results converge towards the results from an exact calculation, e.g., a continuous plane-wave solution [37] of Eq. (5). However, the computational cost increases quartically with N WP , and larger values for N WP will significantly increase the PPD sampling cost. Balancing cost and accuracy, we found it sufficient to draw N = 100 samples from each PPD that we study, since we are quantifying univariate distributions. Also, we noticed that the shapes and widths of the PPDs studied here did not change visibly when varying N WP , and as such we could limit ourselves to N WP ≤ 75 and extrapolate to larger values. This will be discussed in the next section.
At present, using N WP = 75, it takes roughly 12 nodehours (384 core-hours 3 ) to compute all necessary scattering amplitudes at ∼50 scattering energies below 100 MeV for a single configuration of values for the LECs at a specific chiral order ν. This translates to roughly 150k core-hours to compute all scattering amplitudes for 100 different LEC values at four chiral orders. The same calculation with N WP = 150 would be 16 times more expensive and cost roughly 2.5M core-hours. To monitor the reduced method accuracy at N WP = 75, we repeat the PPD sampling, with copies of the same LEC samples, at every chiral order with N WP = 30 and 50. We also use a restricted set of 10 posterior samples with N WP = 100. In addition, we evaluate the PPD at the MAP value α of the LEC PDF using N WP = 30, 50, 75, 100, and 150. The N WP = 75, 150 MAP predictions will be used for extrapolation in the next section.
We had little cost-related reason to restrict calculations to E Lab ≤ 67 MeV. Instead, we computed the on-shell U -matrices at all wave-packet N N Hamiltonian eigenenergies below E Lab = 100 MeV, which was roughly two thirds of the wave-packet basis size. Between these energies we perform linear interpolation of the U -matrix elements to virtually any E Lab < 100 MeV. Consequently, we obtained 100 samples from the PPD of any elastic scattering cross section at every order up to, and including, N 3 LO. Of course, with the neglect of 3NFs, we consider our predictions above E Lab = 67 MeV to be incomplete and have therefore been omitted from the present study. Nonetheless, they allowed us to check on the width and shape of PPDs all the way to E Lab = 100 MeV. Although the HMC-chains of LEC posterior samples are virtually uncorrelated, this does not imply that ensuing samples from the nd cross section PPD are equally uncorrelated. Unfortunately, a chain of 100 samples is typically too short to quantify, e.g., an integrated autocorrelation time or reliably determine the autocorrelation itself. Nevertheless, an inspection of the trace plots of the PPD samples, as shown in Fig. 1, does not indicate any hints of strong correlation between samples. In the event of observing strongly correlated samples, the information content of the PPD chain, as measured by its effective sample size, will drop inversely to the integrated autocorrelation time and we would have to increase the number of samples accordingly [8] B. The differential cross section The convergence of the differential elastic nd cross section at E Lab = 12 MeV with respect to N WP is shown in Fig. 2. Clearly, with N WP ≈ 100, the results begin to stabilize, at least for subleading orders. The somewhat reduced convergence rate for the LO results might be caused by the rather coarse wave-packet representation of the N N potential for low relative momenta [40]. To remedy this one should either re-distribute the discretization boundaries to improve the coverage of the lower momentum region, or simply increase N WP if possible. Since we detect a sufficient convergence at subleading orders, we opt for keeping the discretization mesh the same throughout all calculations and at all chiral orders.
Next, we study the convergence of the PPD with respect to N WP . In Fig. 3  75 are very similar in terms of shape and width. In fact, for all observables that we study in this work 4 , the width and shape of the PPD remains approximately constant as we vary N WP , and the main effect is a shift of the entire distribution. Therefore, we shift the mean of the samples obtained with N WP = 75 using the difference between the MAP predictions obtained with N WP = 75 and N WP = 150. This makes a comparison with experimental data more meaningful.
We did not detect a robust exponential or power-law convergence pattern with respect to N WP and leave further analysis of the N WP -convergence and the WPCD method uncertainty to future work. As such, there might be additional corrections to the PPDs when using N WP > 150 that we do not account for. However, assuming that the widths and shapes of the PPDs remain unchanged, our main conclusions in this work will not be affected.
After shifting the differential cross section obtained with N WP = 75 to N WP = 150 we obtain the result shown in Fig. 4. At all chiral orders and energies we study, the PPD is rather narrow. At LO, the PPD width is comparable to the experimental uncertainty, while at subleading orders the experimental uncertainty is typically greater than the width of the PPD.
To quantify the width of the PPDs, we compute the 90% highest posterior density interval (HPDI), normalize it to the mean of the PPD, and average over θ c.m. . This way, we find that the average HPDI for the differential cross section at E Lab = 12 MeV is, 5.7%, 2.3%, 0.7%, and 0.5% at LO, NLO, N 2 LO, and N 3 LO, respectively. The decreasing values reflects the increasingly narrow LEC posterior densities obtained at higher chiral orders [8,9]. Moving to higher scattering energies we find that the PPDs remain very narrow still. Apart from LO, the average HPDI values are comparable to frequentist estimates of dispersion quantified in Ref. [17], where a similar increase in uncertainty was noted at higher scattering energies.
Recently it was shown that N 2 LO potentials with 3NFs yield an excellent description of differential cross section data [24]. It was suggested in Ref. [42] that 3NFs are necessary to reproduce the differential cross section minimum in the vicinity of E Lab = 65 MeV. Here, however, we see similar reproduction of data at NLO and N 2 LO without 3NFs. Going to N 3 LO, the reproduction of experimental data deteriorates. As shown in Ref. [9], the 3 H and 3 He ground state energies and radii at N 3 LO are also markedly worse compared to N 2 LO. This trend is a testament to the importance of inferring LECs in the N N -and N N N -sectors of χEFT simultaneously [43].
We conclude, based on the inference of N N LECs made in Ref. [9], that the discrepancies between experimental low-energy nd cross section data and theoretical predictions are not due to the uncertainties stemming from the LEC variability. Given the very narrow PDFs  for the LECs, an opposite finding would be a testament to a tremendous fine tuning of scattering observables in the N N N continuum relative to the N N continuum.

C. The EFT truncation error
The truncation of the χEFT expansion used to describe the nuclear interaction leads to a model discrepancy referred to as an EFT truncation error. Following Ref. [26], we assume that the theoretical prediction at chiral order ν for some observable y can be written as where x denotes the kinematic variables E Lab and θ c.m. and y ref is a reference value for the observable in question. This expression renders the expansion coefficients c k dimensionless quantities, which we also expect to be of natural size, i.e., c k ∼ O(1). We assume a χEFT expansion parameter of the form and set the χEFT breakdown scale to Λ b = 600 MeV as in Ref. [9] from where we also obtain the LEC posteriors. We set the c.m. momentum, q, according to the kinetic energy, E Lab , of the incoming nucleon. The χEFT truncation error, δy ν , is the expected magnitude of the sum of contributions from terms beyond the order ν. Under the assumption of having independent and normally distributed expansion coefficients, c k , it is shown in, e.g., Ref. [34], that δy ν is also normally distributed and given by wherec 2 denotes the variance of the expansion coefficients. Thus, knowingc 2 enables us to quantify the (variance of the) χEFT truncation error. For this purpose, we follow the procedure of, e.g., Ref. [8] and employ the root-mean-square (RMS) value of order-by-order differences to estimatec 2 . The order-by-order differences are computed from the mean values of the PPDs at each order ν, thus averaging over a possible LEC dependence. We wish to compare the magnitude of the χEFT truncation error with the theoretical error in y stemming from the uncertainty about the numerical values of the LECs. Let us take the differential cross section at E Lab = 12 MeV as an example and inspect it closer. Limiting ourselves to this low value of E Lab , the effect of 3NFs are expected to be small [42]. Therefore, we retain the expansion in Eq. (6) and use Eq. (8) to quantify the χEFT truncation error. We set y ref to the LO prediction. At this scattering energy, we also have Q = m π /Λ b ≈ 0.23. An RMS estimate from the expansion coefficients at θ c.m. = 30, 90 and 150 degrees (omitting LO results due to their role in the definition of y ref ) givesc = 14.8. This is a fairly unnatural value which arises from an oscillating convergence when including higher chiral orders. 11 12 13 14 dσ/dΩ (mb) To set the absolute scale, we included the experimental measurement (gray) from Ref. [44] The PPDs due to the LEC variabilities and the χEFT truncation errors are compared in Fig. 5. Clearly, the truncation error (solid lines) is typically much greater than the error due to the uncertain values of the N N LECs (histograms) up to an including N 2 LO. Therefore, we find it unnecessary to account for a possible LEC variability in the expansion coefficients in Eq. (6). At N 3 LO, the two errors are becoming comparable. However, at this order, both of the errors are tiny, 4%, on an absolute scale. In fact, they are both smaller than typical experimental errors, indicated as the gray area in Fig. 5. In addition to the RMS estimate ofc 2 we also show the truncation error (dashed line) based on a naturalness assumption where we setc 2 = 1. In this limit, the two errors become comparable for this observable already at NLO.
At higher energies, we see in Fig, 4 that the predictions at N 3 LO deviates from the ones at NLO and N 2 LO. When analyzing the truncation errors at E Lab = 36 MeV, we obtainc = 65.1, which signals the presence of an unnaturally large contribution in the χEFT expansion. We find that omitting the shift between N 2 LO and N 3 LO has a significant impact and yields a more reasonable value ofc = 15.1. Doing the same at E Lab = 12 MeV yields c = 11.3, i.e., a relatively small change from when including the shift. The truncation error is expected to increase with the on-shell energy, and thus it should become greater than the LEC uncertainty for E Lab > 12 MeV, but we leave a more detailed study for future work.

D. Spin-polarization observables
There are many different possibilities to form observables related to spin-polarization in the initial and/or fi-nal states of the N d reactants [45]. The fine details of the angular dependence of these observables can depend sensitively on the spin structure of the N N and N N N interactions. A well-known example is the low-energy vector analyzing power A y . This observable depends sensitively on the 3 P partial waves of the N N -interaction [46,47]. There are indications that it also depends sensitively on parts of the subleading 3NF [24]. It has turned out to be very challenging to reproduce the experimental data for this observable at laboratory scattering energies E Lab 30 MeV [37,48].
Given the possibly fine-tuned nature of A y , it is particularly interesting to quantify the PPD due to the variability in the N N LECs of χEFT. In Fig. 6, we show the PPDs for A y at NLO, N 2 LO, and N 3 LO as well as the average 90% credibility intervals. At E Lab = 10 MeV we do not reproduce the experimental data at any chiral order. We note that the N 3 LO calculation appears to improve the description of the data at the polarization maximum. However, the low-angle description is markedly worse compared to the result at N 2 LO. For E Lab ≈ 35 − 67 MeV it appears sufficient to use N Nonly forces at N 2 LO to describe presently available data.
It is clear that the variability due to the LECs inferred from N N data does not give rise to any significant uncertainty nor does it explain discrepancies between theory and data. We refrain from quantifying the χEFT truncation error for this observable since our calculation omits 3NFs, which may very well play a significant role in explaining the low-energy A y values. Nevertheless, a crude estimate to account for the χEFT truncation error with missing 3NFs can be obtained by pulling out factors of Q in Eq. (8), starting at N 2 LO [50]. We found that this procedure induced rather large χEFT uncertainties that covered the experimental data at all orders.
As for the remaining spin-polarization observables, their N N PPDs exhibit similar patterns and widths as presented above for the differential cross section and A y (n), i.e., the vastly dominating source of uncertainty is the χEFT truncation error, at least below N 3 LO.

V. SUMMARY AND OUTLOOK
We sampled the PPDs for the nd differential cross section dσ/dΩ at E Lab = 12, 36, and 65 MeV scattering energy, and neutron analyzing power A y (n) at E Lab = 10, 35, and 67 MeV. The underlying samples from the LEC posterior were obtained from a previous analysis of N N data [9]. The HMC algorithm used in that analysis yields virtually uncorrelated samples which we find most likely persists for the elastic nd observables. The main conclusion from this work is that the uncertainty about N N LECs, when conditioned on N N scattering data and uncorrelated estimates of the χEFT truncation errors, does not entail significant uncertainties in the low-energy nd continuum. Although we only show results for selected observables, we find them to be representative of all elas- tic nd scattering observables, at least for E Lab 67 MeV. When compared with estimates of the χEFT truncation error, we find that the uncertainty stemming from the numerical values of the N N LECs are negligible, at least up to (and including) N 2 LO in Weinberg power counting. At N 3 LO, the width of the PPD and the credible interval of the truncation error are starting to become comparable. However, these uncertainties are very small and, in fact, are comparable to typical experimental errors.
In this work we have not quantified the errors due to having a finite number of wave-packets in the WPCD method. Instead, we extrapolated all results to N WP = 150 and relied on the fact that the widths and shapes of all studied PPDs remain the same when using fewer wave-packets, i.e., N WP = 50 and 75. Future work should be dedicated to understanding the scaling of the WPCD method-error with respect to the discretization of the continuum.
Throughout our analysis, the PPDs were conditioned on N N scattering data. For the predicted differential cross section, we find reasonable agreement with experimental N d scattering data. The same observation was made for many polarization observables, not shown explicitly in this paper. However, less accuracy is observed in the low-energy A y analyzing power. A natural next step would therefore be to simultaneously infer the N N and N N N LECs from N N plus N d scattering data. This would shed more light on the necessity of including 3NFs to explain this data.
The inference of LECs in χEFT is not restricted to use only scattering observables. In fact, any low-energy nuclear data can be utilized (and will be relevant given that it has a high information content). On the other hand, the abundant sets of experimentally measured N N [32,33], πN [51], and N d [23] scattering cross sections provide data where theoretical predictions do not rely on many-body interactions beyond 3NFs. In addition, a scattering cross section can be tied to a welldefined (external) momentum, providing a clear interpretation of the soft scale entering the χEFT expansion parameter Q and the associated truncation error. This identification of a soft scale is more ambiguous in bound states of nuclear many-body systems.
A Bayesian analysis of LECs in χEFT conditioned on N d data requires efficient solutions to the AGS equations. Indeed, traversing larger domains of the multidimensional LEC parameter-spaces would require orders of magnitude more samples than what we employed in this work. Fortunately, recent advances in model reduction methods [52], utilizing singular value decomposition [53] and eigenvector continuation [54][55][56] methods, show great promise in delivering accurate and fast solutions to the Faddeev equations. Some of these methods appear compatible with our existing implementation for solving the AGS equations with the WPCD method.