Inference of the low-energy constants in $\Delta$-full chiral effective field theory including a correlated truncation error

We sample the posterior probability distributions of the low-energy constants (LECs) in $\Delta$-full chiral effective field theory ($\chi$EFT) up to third order. We use eigenvector continuation for fast and accurate emulation of the likelihood and Hamiltonian Monte Carlo to draw effectively independent samples from the posteriors. Our Bayesian inference is conditioned on the Granada database of neutron-proton ($np$) cross sections and polarizations. We use priors grounded in $\chi$EFT assumptions and a Roy-Steiner analysis of pion-nucleon scattering data. We model correlated EFT truncation errors using a two-feature Gaussian process, and find correlation lengths for $np$ scattering energies and angles in the ranges 45--83 MeV and 24--39 degrees, respectively. These correlations yield a non-diagonal covariance matrix and reduce the number of independent scattering data with a factor of 8 and 4 at the second and third chiral orders, respectively. The relatively small difference between the second and third order predictions in $\Delta$-full $\chi$EFT suppresses the marginal variance of the truncation error and the effects of its correlation structure. Our results are particularly important for analyzing the predictive capabilities in \textit{ab initio} nuclear theory.


I. INTRODUCTION
A chiral effective field theory (χEFT) description of the nuclear interaction [1][2][3][4][5] is endowed with a power counting (PC) to organize the order-by-order contributions of the strong-interaction dynamics to nuclear observables.This kind of organization, regardless of the particulars of the adopted PC [6], is a hallmark of EFT [7] and ab initio [8] approaches to nuclear theory since it promises a handle on the theoretical uncertainty coming from truncating the EFT expansion [9].Accounting for the truncation error is key to mitigating overfitting of the low-energy constants (LECs) as well as assessing the importance of discrepancies.Indeed, the modeling of EFT truncation errors can play a significant role in the robustness of LEC inferences and ensuing nuclear predictions [10][11][12][13][14][15][16].Melendez et al. [17] have proposed a Bayesian model for the truncation error that accounts for finite correlations across independent variables, e.g., the scattering energy and angle for nucleon-nucleon (N N ) scattering cross sections and polarizations.To date, LEC inference in χEFT typically accounts for truncation errors in the fully correlated or uncorrelated limits [18][19][20][21][22], and the robustness with respect to the correlation structure is not well known.
In this paper, we quantify a correlated truncation error and analyze its effects on a Bayesian estimation of the LEC posteriors for a ∆-full χEFT description of the neutron-proton (np) interaction up to next-to-nextto-leading order (NNLO) [23,24].This extends our previous work on Bayesian LEC estimation in ∆-less χEFT where we employed an uncorrelated truncation error [14,21,22].We also use eigenvector continuation (EC) [25], i.e., a reduced basis method [26], to efficiently and accurately emulate the scattering amplitudes entering the np scattering-data likelihood.Following our previous publications [21,22], we employ Hamiltonian Monte Carlo (HMC) [27] to draw effectively independent samples from the LEC posteriors.
The ∆(1232) resonance plays an important role in nuclear physics since it represents a rather low excitation energy and couples strongly to the pion-nucleon (πN ) system.This was recognized already in early χEFT descriptions of the N N interaction [28][29][30], and several modern χEFT interactions incorporate the ∆ as well [31][32][33][34].In ∆-full χEFT there are four subleading πN LECs up to NNLO, usually denoted c 1 , c 2 , c 3 , c 4 .They govern the strength of subleading 2π-exchange diagrams of the N N interaction and the leading threenucleon (N N N ) 2π exchange with an intermediate ∆ excitation, the so-called Fujita-Miyazawa force [35,36].A Roy-Steiner analysis of πN scattering amplitudes by Hoferichter et al. [37] has enabled a determination of the πN LECs.Unfortunately, the relatively unknown value of the axial πN ∆ coupling h A propagates to approximately five times greater uncertainties for c 2 , c 3 , c 4 (compared to a determination in ∆-less χEFT) when matching ∆-full χEFT to the subthreshold amplitudes in the Roy-Steiner formalism [38].A well-founded truncation error would therefore pave the way for learning more about the strength of subleading 2π exchange, and leading N N N forces, from N N scattering data.This paper is organized as follows: In Sec.II we present our statistical model for linking experiment and χEFT.In Sec.III we discuss our priors and likelihood, and in Sec.IV we introduce the two-feature Gaussian-process model of the correlated EFT truncation error.In Sec.V we discuss the training of EC np scattering emulators.In Sec.VI we present the results from HMC sampling of the LEC posteriors.A summary and an outlook are given in Sec.VII.

II. LINKING EXPERIMENT AND THEORY
Following our previous papers [14,21,22], we relate an experimental measurement y exp of some scattering observable with a theoretical prediction y (k) th , up to chiral order k, using an additive model to account for the respective uncertainties, δy exp and δy  ( The theory prediction y th depends deterministically on the vector of LECs ⃗ α and the independent variable ⃗ x = (T lab , θ), where T lab denotes the kinetic energy of the incoming nucleon in the laboratory frame and θ denotes the scattering angle in the center-of-mass frame.For the total scattering cross section we have ⃗ x = (T lab ) as this observable is integrated over all θ.We have suppressed the explicit ⃗ α-dependence of δy th as we will assume a fixed parameter value in the model of the theory error.
The composition of the LEC vector ⃗ α depends on the chiral order.For the np potentials in this paper, the LEC vector up to each order is given by where LO is leading order (k = 0), NLO is next-toleading-order (k = 2), and NNLO is next-to-next-toleading order (k = 3). 1 We employ units and a notation linked to the momentum partial-wave basis; see Refs.[4,23] for details.The potential is nonlocally regulated with a regulator cutoff Λ = 450 MeV as in Ref. [24].
In the following we will only refer to the generic vector ⃗ α, while the specific chiral order, if important, should be obvious from the context.The power counting of ∆-full χEFT allows us to express a prediction y (k) th (⃗ α; ⃗ x) up to chiral order k as a sum of order-by-order contributions where y ref is some characteristic scale to be decided, and c (i) are dimensionless expansion coefficients.The dimensionless expansion parameter Q = f (m π , δ, p)/Λ χ is a 1 Due to symmetries [3][4][5] the order k = 1 vanishes.
ratio composed of a soft scale given by some function of the pion mass (m π ), the ∆N mass splitting δ, the (external) relative N N momentum p, and a hard scale Λ χ for which we assume a point-estimate value of 600 MeV [18,39].We treat δ ∼ m π as a small scale [40], although it does not vanish in the chiral limit.We resum the potential to all orders which obfuscates the form of f , and in line with [17] we therefore assume the following functional form which facilitates a smooth transition across the soft scale m π .We find that the exact form of this function does not impact our inference significantly, and reverting to a canonical estimate f = max(m π , p) does not change any of our results.The upshot of having an order-byorder description of y th , as in Eq. ( 5), is a handle on the theoretical uncertainty via Clearly, we have not explicitly extracted EFT coefficients c (i>k) for any ⃗ x and as a consequence we are uncertain about their values.However, if we assume naturalness to hold, it is reasonable to have c (i) ∼ O(1) for all expansion coefficients, including the unseen ones.Furthermore, contributions can be either positive or negative.We will therefore consider the coefficients as drawn from some underlying probability distribution bounded by a finite variance c2 expected to be of order 1.In addition, there is most likely information about y th (⃗ x) at nearby values of the independent variable.As such, there is likely nonzero covariance between expansion coefficients c (i) (⃗ x m ) and c (i) (⃗ x n ) at different values ⃗ x m and ⃗ x n .That correlation structure is of primary interest in this paper.

III. SETTING UP THE LEC PARAMETER INFERENCE
We seek a posterior probability density function (PDF) pr(⃗ α|D, I) for the LECs ⃗ α in ∆-full χEFT up to NNLO conditioned on selected np scattering data D from the Granada database [41,42] and other information I as specified below.In the following, we use bold symbols to denote quantities that depend on a range of values for the independent variable ⃗ x, e.g., x = (⃗ x 1 , ⃗ x 2 , . ..) and D = (y exp (⃗ x 1 ), y exp (⃗ x 2 ), . ..),where we let y denote observable type (differential cross section, vector polarization, etc.).
We ignore the evidence term in this paper as it provides overall normalization and does not impact the shape of the posterior.

A. Prior
We assume independent priors for the LECs of the N N contact and πN potentials, i.e., pr where ⃗ α πN = (c 1 , c 2 , c 3 , c 4 )2 denotes the subleading πN LECs.For these, we place a normal and correlated prior based on the (nondiagonal) covariance matrix Σ πN obtained in the Roy-Steiner analysis of πN scattering data by Siemens et al. [38], Along the lines of naturalness, we place a normal and uncorrelated prior on the LECs ⃗ α N N of the contact potential, where Σ N N is a diagonal covariance matrix with standard deviations 5 × 10 4 GeV −(k+2) for the contact LECs first appearing at orders k = 0, 2. This is the same prior as in our previous papers [21,22].

B. Likelihood
We consider only strong interaction effects in elastic reaction channels, and neglect all electromagnetic interactions in our calculations.Thus, we condition our LEC inferences on low-energy np scattering data below the pion-production threshold T lab = 290 MeV, and omit np scattering data with T lab < 30 MeV as the nonzero magnetic moments of neutrons and protons can distort the low-energy np scattering amplitude [43,44].We partition the experimental data set as D = {D y } Ny y=1 , where N y = 18 is the number of unique types of scattering cross sections and polarizations in the considered energy range.In Fig. 1 we summarize in tabular form the data of measured np scattering observables that we include, as well as some of the results that will be discussed below.
We assume that experimental and theoretical errors are independent of each other such that the total covariance matrix can be written as The (diagonal) covariance matrix Σ exp is provided by the experimenters and we employ normalization factors from the Granada database.We model the covariance of the EFT errors independently for each type of np scattering observable y, and make specific assumptions and choices I y per observable type.The (block-diagonal) covariance matrix Σ th is given by where the (non-diagonal) covariance matrix Σ th,y at chiral order k contains elements that we model using a Gaussian process; see Sec.IV.
In accordance with the principle of maximum entropy we employ a normally distributed data likelihood, which factorizes as pr(D|⃗ α, I) = Ny y=1 pr(D y |⃗ α, I y ), (17) where (18) Here, r y (⃗ α) = [r y,1 (⃗ α), r y,2 (⃗ α), . . ., r y,N d,y (⃗ α)] is a (column) vector of N d,y residuals, each one given by the difference between experimental and order-k theoretical prediction values for the independent variable ⃗ x, i.e., All cross sections are computed from amplitudes obtained numerically via the nonrelativistic Lippmann-Schwinger equation [12,21].This approach is virtually exact and we do not account for any method uncertainties of numerical origin in this paper.Note that in Eq. (19) we assume that all errors taken into account have mean zero, which will be the case in this paper.[41,42], number of data N d,y and measurement energies NT lab ,y per observable type y, effective dimension n eff given theory error correlations, maximum a posteriori (MAP) values for GP correlation lengths across scattering energy ℓT lab and angle ℓ θ , and marginal variance c2 .The spin-dependent integrated cross sections σT and σL have polarizations parallel and anti-parallel to each other and transverse (T ) or longitudinal (L) to the beam direction.The σT observable was omitted in our analysis, see the text for details.For spin polarizations we list the polarization tensors, X srbt , using the notation of [45] where subscripts denote the polarization of the (s)catter, (r)ecoil, (b)eam, and (t)arget particles in the direction of the incoming projectile ( k), normal to the scattering plane (n), ŝ = n × k, and 0 is an unpolarized state.The illustrations show initial and final spin-polarizations in the laboratory frame using gray arrows for reactants, blue arrows for the polarization direction in the scattering plane, and circled dots for the outward-pointing polarization vector.The projectile impinges from the left, the target (not indicated) is located in the center, and the projectile (recoil) scatters upwards (downwards) an angle indicated by the arc (gray).In this table, and in the text, we leave the energy dependence of all observables implicit.

IV. CORRELATING THE EFT TRUNCATION ERROR IN N N SCATTERING
Following the suggestion by Melendez et al. [17], we model the EFT truncation error using Gaussian processes (GPs) [46,47].This results in a covariance matrix Σ th from Eq. ( 15) with off-diagonal elements.In this section, we limit our discussion to the GP modeling of a single covariance submatrix Σ th,y for a specific observable type y, e.g., the differential cross section σ(θ), and therefore omit the observable index y in the following.
We develop a GP model to handle two features, i.e., GP : R 2 → R 1 , for correlating the distribution of expansion coefficients at different values of the scattering energy and angle.For total cross sections we only operate with the scattering energy, as all angles are integrated over.We expect unseen coefficients c (i) (x) ≡ {c (i) (⃗ x j )} j to be distributed as where m is assumed to be a constant mean function and ℓ, c2 denote the hyperparameters for the GP correlation length(s) and marginal variance, respectively.The calibration of these parameters will be discussed in Sec.IV A.
The marginal standard deviation c is a scale factor that quantifies the average distance of the GP away from its mean, and the correlation length quantifies the typical extent of the GP fluctuations.Elastic np scattering observables typically exhibit a smooth behavior as a function of energy and angle.Therefore, we expect smooth variations of the expansion coefficients.Assuming stationarity, we employ a squared-exponential kernel to parametrize the correlation structure, where the length scale(s) are given by in an obvious notation.
We are mainly interested in the quantified hyperparameters ℓ and c2 that characterize the finite correlation length and variance of the EFT error.A GP is a linear-Gaussian model for which there are closed-form expressions [46,47] for the distribution of predicted data conditioned on calibration data.Once the posterior distribution, or point estimates, of the hyperparameters are determined it is straightforward to quantitatively use the GP to predict EFT expansion coefficients exhibiting a variance and correlation consistent with the calibration data.Furthermore, one can show [17] that the GP for the EFT truncation error pertaining to a scattering observable y th is given by where the respective mean (M ) and kernel (K) functions are given by and This lays the foundation for how we model the covariance matrix in Eq. ( 15) and sample the LEC posterior pr(⃗ α|D, I) via the likelihood and prior defined in Secs.III A and III B. As will be discussed shortly, we set m = 0 in this paper.

A. Optimizing the GP hyperparameters
We set up individual GPs for the block-diagonal covariance matrices Σ th,y in Eq. ( 15).The hyperparameters for each observable type y are optimized separately and we find maximum a posteriori (MAP) values for them using gsum [17].For each GP we use a set of training data c ⋆ consisting of N train expansion coefficients obtained as for a training range x ⋆ of values for the independent variable.We pick N train ≈ 25-50 points uniformly distributed across approximately 4-6 energies T lab ∈ [30,290] MeV and 5-10 angles θ ∈ [0, 180] degrees.The training grid varies somewhat depending on the observable type.The order-by-order differences in Eq. ( 26) are computed from theoretical predictions based on a MAP estimate for the LECs, ⃗ α = ⃗ α ⋆ , obtained in line with Ref. [21] in the uncorrelated limit using c = 1.
The LO predictions y th in ∆-full χEFT are identical to the ∆-less predictions and as such they are rather poor, as expected in Weinberg PC [12,48].Thus, the c (2) expansion coefficients incur rather unnatural values.One could argue that a Bayesian approach, with a prior for the expansion parameters, should be able to properly handle such outliers.However, a (likely) deficient LO will conflict with our initial assumptions in Eq. (7).Instead of using the truncation error to absorb deficiencies rooted in the PC, we decided to neglect LO predictions in this paper.Thus, all training data consist of c (3) expansion coefficients only.Despite this limited amount of data, we expect the available coefficients to carry valuable information about the correlated truncation error and we proceed with a Bayesian GP model.For total and differential cross sections we set y ref to the predictions of the Idaho-N3LO(500) potential [49], which is almost the same as using y ref = y exp .We encounter a zero crossing in y ref for σ T .This leads to a discontinuity in the c (3) coefficients as a function of T lab ; see Fig. 2. Such small-scale structures are virtually impossible to faithfully represent using an infinitely smooth squared-exponential kernel.As there are only three experimental data points for this observable type, we decided to exclude them from D rather than choosing an exceptional y ref for this one case.For polarization data we encounter an excess of small-scale structures when setting y ref to the predictions of Idaho-N3LO(500).Using y ref = 0.5 removes them entirely and we are spared further kernel design.We choose 0.5 as a representative value for these observables.Other choices could be y ref = 1.0 [18] or an average of experimental values [21].
All training data are pruned for outliers using a threeinterquartiles threshold [21].In the end, all c (3) (x) coefficients are of natural order, and a vast majority of them pass the outlier filter.Operating with a GP allows us to incorporate known symmetry constraints on the polarization observables.Indeed, some polarizations are identically equal to zero at θ = 0 and/or θ = 180 degrees.We impose such boundary constraints by adding zerovalued expansion coefficients for the endpoint(s) of the angular direction to the training data.In a future application we will formally incorporate continuous boundary constraints in the GP [50].
We employ a bounded and uniform prior for the length scales ℓ of the form pr(ℓ θ |I) = U(ϵ, 180) degrees (27)

and pr(ℓ
where ϵ = 10 −5 is introduced to avoid numerical issues in the optimization of the kernel (the exact value of ϵ is unimportant in this case since the posterior values are never close to the edge of the prior).We place a conjugate inverse-χ 2 prior on the variance c2 according to and a Dirac delta prior pr(m) = δ(m) on the mean m.This encodes our expectation of natural c2 values, although the heavy tail in our prior allows for some unnaturalness, and expansion coefficients symmetrically distributed with mean m ≡ 0. Our chosen prior for c2 is shown in orange in Fig. 3.We do not allow the mean to vary owing to the limited amount of informationone chiral order-we have available to learn from.This strict prior on m could be substituted for, e.g., a Gaussian prior if more information is available.If so, Eq. ( 19) must be updated to include the systematic contribution of a truncation error with non-zero mean.
The conjugacy of the χ −2 (c 2 ) prior implies that the posterior pr(c 2 |m = 0, ℓ, c ⋆ , x, I) is also an inverse-χ 2 PDF with updated parameters [17].The prior for the length scale ℓ is not conjugate, and to find the posterior MAP values ( ℓ, c2 ) of we use the numerical optimization routine L-BFGS-B [51] with multiple restarts.We include a so-called nugget, set to 10 −10 , to avoid numerical issues.The resulting MAP estimates of the hyperparameters are listed in the table in Fig. 1.The delta prior on m obviously yields a corresponding MAP value m = 0.The GP predictions for expansion coefficients c (3) for three of the most abundant observable types are shown in Fig. 4. The amount of training data is sufficient to make the posterior likelihood-dominated in most cases.We explore a number of priors, shown in Fig. 3, and find that the inference is rather robust as long as we allow for very small values c2 ≳ 0.2 2 .Indeed, a prior that is too restrictive in this regard can bias the result for total cross sections, where the training data are one dimensional.

B. Validating the GP hyperparameters
We validate the GP hyperparameters for each observable type y using a set of N val complementary validation data c.The validation data are generated by the same method as the training data, but with energies and angles shifted so that the validation values x for the independent variable do not overlap with ones used during training.In addition to visual inspections, which are certainly useful, we also employ a set of diagnostics inspired by Refs.[17,52].These diagnostics are the Mahalanobis distance, pivoted Cholesky errors, and the credible interval diagnostic.These, and a few more, are thoroughly discussed in Ref. [17].
The (squared) Mahalanobis distance, D 2 MD , is a multivariate analog to the sum of squared residuals.Here, it is computed as where m is the GP mean at x-equivalent to zero with our choice of prior-and This distance is commonly used to quantify the deviation of a prediction compared to data in a correlated setting.
Here, we use it to diagnose whether a set of validation data could reasonably have been drawn from a GP with covariance according to Eq. ( 32) and mean m = m = 0. Either too large or small values of D 2 MD , compared to the reference χ 2 distribution with N val degrees of freedom, point to a possible conflict between the GP and the validation data.The Mahalanobis distance is a scalar measure and does not provide detailed insight about the possible tension with respect to the validation data.
Furthermore, it is instructive to study the GP and the validation data point by point as a function of the independent variable.Such comparisons are not straightforward to interpret given the mutual correlation of the validation data.Therefore, we decorrelate and rescale the covariance matrix to independent unit variances using a Cholesky decomposition; Σ GP = GG T where G is a triangular standard deviation matrix.From this we define the Cholesky errors To order the vector D G in a meaningful way we pivot the decomposition in decreasing conditional variances.One should not detect any pattern when plotting the pivoted D G versus the index of the validation data.To reveal further information one can introduce a ratio scale also on the abscissa by plotting the pivoted Cholesky errors versus the conditional variances used for pivoting [53].
All but one of the GPs readily pass our diagnostics, with D 2 MD landing within 95% credible intervals of the respective reference distributions.We see little to no structure in the pivoted Cholesky decompositions, and empirical coverages that roughly match the corresponding credible intervals.In addition, visual inspections indicate that the inferred hyperparameters are plausible.Naturally, incorporating expansion coefficients from other orders would provide a stronger test of our GP model.We only reject σ T , and therefore also the three experimental data points.Neither D 2 MD nor the D G exhibit significantly poor performance for this observable given our validation data.However, upon visual inspection of the c (3) curve for this observable, shown in Fig. 2, we find that it is discontinuous for T lab ≈ 205 MeV because the reference value crosses zero at this point, leading to extremely large coefficients near the crossing point and an abrupt change in sign.Clearly, our GP model, based on a squared-exponential kernel, is ill suited to handling this discontinuity, and the inferred length scale is strongly dependent on the chosen training data and will approach zero as we increase the number of training data.If it were not for the existence of an experimental data point in the vicinity of the problematic scattering energy we would have considered excluding only the problematic region.This example clearly shows that one should not blindly trust the diagnostics as they are conditioned on the chosen training and validation data.

C. The structure of the correlated truncation error
Introducing a correlated truncation error reduces the number of independent directions in the data space.Indeed, two data residing within one correlation length of the scattering energy or angle carry joint information.To quantify the impact of this we compute the effective dimension n eff , per observable type y, using a measure [54,55] , where γ i denotes a normalized eigenvalue γ i = λ i /tr(C y ) of the respective correlation matrix.Here, C y = S −1 y Σ y S −1 y and S y = diag(Σ y ) with Σ y = Σ exp,y + Σ th,y .We could equally well compute n eff for Σ y directly.However, operating with the correlation matrix leads to the values n eff = 1 and n eff = N d,y in the limits of having a full off-diagonal correlation and zero off-diagonal correlation, respectively.Consequently, we interpret n eff as the effective number of data, and note that the logarithm of n eff corresponds to the discrete Shannon entropy of the spectrum of the normalized eigenvalues γ i .
The resulting n eff values at NNLO are summarized in the table in Fig. 1.They show that the correlations reduce the number of effective data from 2087 to 841, which is still plentiful.The main reason for the relatively weak impact of the correlated truncation error is the dominance of the experimental variances along the diagonal of the covariance matrix Σ.Indeed, for the correlation matrices of the GP kernel for EFT truncation errors alone we find n eff ≈ 5 − 10.In Fig. 5 we show correlation matrices for the total cross section C σtot , with and without the experimental covariance matrix accounted for.The diagonal dominance of the experimental errors is clearly visible.This dominance weakens for higher energies (higher indices in the figure) as the truncation error increases with the scattering energy T lab .At NLO, where the truncation error is greater by one chiral order Q, the corresponding n eff values are typically a factor 2 smaller as the correlations of the truncation error become relatively more important.The total number of effective data is 498.In Fig. 6(a) we compare the 68% credible interval ∆y th of the EFT truncation error at NLO and NNLO with the experimental errors ∆y exp for all data.
The EFT truncation error, and its correlations, is expected to be more important for inferences conditioned on very precise pp data.To estimate this effect, we inspected the pp A y (θ) data set more closely.With T lab ∈ [30,290] MeV there are 496 data points, which is almost the same as in the np sector.Reducing the experimental variances of the np A y (θ) data, to mimic the average level of pp variances, we observe n eff = 154 and 91 at NNLO and NLO, respectively.For other observable types, like σ(θ), the average pp variance is greater than the average np variance, and the opposite effect is likely observed.Note that the distribution of np and pp measurement energies and angles differ and we do not account for this in our estimate.

V. EMULATING np SCATTERING CROSS SECTIONS
We use an EC-based method [25] to construct accurate and efficient emulators for np scattering observables.Operating with emulators, instead of the exact simulators, helps reduce the computational cost of evaluating the likelihood.As an added bonus, once the emulators are trained they are straightforward to distribute in the scientific community.
We use Newton's functional formulation of on-shell np T -matrix elements for setting up the emulators [56].Technically, this leads to one emulator per np partialwave and unique scattering energy in the database D (table in Fig. 1).Truncating the partial-wave expansion at a maximum np angular-momentum quantum number J = 30 leads to 182 partial waves.With N T lab = 177 energies in D we end up with 32 214 T -matrix emulators per chiral order.The training values for the LECs are drawn according to a latin hypercube design in a sufficiently wide interval [−4, +4] (see Sec. III A for units).To simplify the setup, we employ the same training protocol for all emulators and find it sufficient to use 7 ( 8) training points at NLO (NNLO).As there are at most 3 (7) relevant LECs acting per partial wave, this approach leads to very accurate emulation of all relevant scattering observables.We estimate that emulator errors ∆y emu are at least 10, and typically 1000, times smaller than experimental errors.This is quantified using the difference between emulated and exact results for all np observables in D at ten different sets of LEC values following randomized latin hypercube designs within the training intervals specified above.In Fig. 6(b) we show the ratio ∆y emu /∆y exp at NNLO for a random set of LEC values.This is virtually identical to the NLO result.We conclude that the emulator errors are sufficiently small to neglect them in our inferences.In the following we refer to the collection of all amplitude-emulators, per chiral order, as the emulator for scattering observables y.

VI. SAMPLING POSTERIOR DISTRIBUTIONS
We sample the LEC posteriors pr(⃗ α|D, I) at NLO and NNLO using HMC, an advanced Markov chain Monte Carlo (MCMC) algorithm suitable for high-dimensional sampling problems.We have detailed the use and performance of HMC in a previous paper [21] and further improved its performance in Ref. [22].Here, we employ the same sampling strategy as in the latter paper.HMC needs tuning in order to perform well, and in particular a good so-called mass matrix is needed.The mass matrix should be a decent approximation to the inverse of the parameter covariance matrix.We extract such an approximation by maximizing the posterior through the standard BFGS algorithm [57][58][59][60].
HMC requires gradients of the (logarithm of the) posterior with respect to the sampled parameters.The underlying linearity of the emulators furnishes easy access to derivatives of T -matrix amplitudes with respect to the LECs.However, we employ the automatic differentiation (AD) and just-in-time (JIT) compilation tools for Python in the Google JAX [61] library to compute the required gradients and boost the execution speed of our Python code; it takes on the order of one second to evaluate the entire data likelihood, including derivatives.We opted for this approach due to its simplicity and speed.It also enables straightforward computation of derivatives of the posterior with respect to parameters other than the LECs, such as GP hyperparameters.A threaded Ccode implementation would likely be more efficient, but not pivotal for the present applications.
We diagnose the convergence of the MCMC chains using a criterion based on the integrated autocorrelation time τ int [62].This is a measure of the autocorrelation between MCMC samples, and it tends to be underestimated for short MCMC chains.In line with Ref. [21] we therefore declare convergence when (a) the estimation of τ int has stabilized, and (b) when N ≥ 50τ int , where N is the length of the MCMC chain.All our chains readily pass this test.Like other MCMC convergence tests, the τ int criterion may falsely declare convergence, e.g., if the posterior is multimodal; we therefore search for local optima using BFGS optimization initialized at random positions.We have not detected signs of multimodality and are confident that our chains are converged.

A. LEC posteriors
In a first step we compare the LEC posteriors obtained with correlated and uncorrelated EFT truncation errors.In the uncorrelated limit we determined c2 from the rootmean-square value of the NLO-NNLO expansion coefficients for the σ tot , σ(θ), A y (θ) observables, i.e., omitting the LO contribution and discarding outliers.We employed the same grids of scattering energies and angles as in the training of the correlated GP model for these observables, see Sec.IV A. This leads to a fixed c2 = 0.42 2 .
We find that the introduction of a finite correlation length in the EFT truncation error does not affect the marginal LEC posteriors much; see Figs. 9, 10, 11, and 12 in Appendix A. However, there are some differences and similarities worth commenting on.Most noticeably, the marginal posteriors for all LECs are approximately twice as wide when including a correlated truncation error.This is consistent with the observed reduction in the effective number of data.Other than that, the posterior correlation structure remains the same and the respective locations of the modes are largely the same, except for a significant shift in the value of the C 3S1−3D1 contact and the C 1P 1 contact LEC becoming (almost) consistent with zero in the correlated (uncorrelated) limit at NLO (NNLO).At NNLO the MAP values for the subleading πN LECs are substantially shifted with respect to the mean of their prior values in Eq. ( 10).This result does not change when neglecting a correlated truncation error.In the case of a correlated truncation error we obtain a MAP value for the πN LECs, ⃗ α πN = [−0.72(2),−0.96(5), −0.01(6), 0.69(5)] GeV −1 , (34) with 68% credible intervals indicated by error bars in parentheses.The squared Mahalanobis distance of this point with respect to the mean and covariance of the prior is D 2 MD ( ⃗ α πN , ⃗ µ πN ; Σ πN ) = 9.65, which is just far enough for the posterior MAP value to be outside 95% of the prior probability mass.This can be cast in terms of a p value (0.047) and traditional significance testing, which would lead us to reject the correctness of the Roy-Steiner prior (null hypothesis) on the 5% significance level.However, we are inclined to place doubt on our model for the truncation error and in particular its variance.As a side remark, the D 2 MD values for the MAP values of the πN LEC posteriors for ∆-less NNLO [21] and for the next order N3LO [22], inferred using an uncorrelated truncation error, are even greater and again point to significant tension between the respective LECs inferred using N N data and πN data.
We suspect that the truncation error is underestimated since the NLO-NNLO shift in the ∆-full theory is on the smaller side, which is not too surprising as the inclusion of the ∆-isobar pulls down higher-order contributions in the ∆-less theory.Also, we cannot rule out that the contributions at ∆-full N3LO are substantial owing to the introduction of a rich fourth-order contact potential.To shed some light on the possible underestimation of c2 we sampled the joint posterior pr(⃗ α, c2 |D, I) in the uncorrelated limit of the truncation error.For this, we employed the same LEC prior as defined in Sec.III A, and assumed the same c2 for all observable types.For the latter parameter we employed an inverse-χ 2 prior with ν = 23.75 and τ 2 = 0.19, as obtained via conjugacy of an inverse−χ 2 (hyper)prior with ν 0 = τ 2 0 = 1 updated using training data from expansion coefficients c (3) for σ(θ), A y (θ), σ tot on the same grid of θ and T lab values as used when inferring the GP model parameters for the correlated truncation error.This prior is sharply peaked at c2 ≈ 0.45 2 .The marginal LEC values from the pr(⃗ α, c2 |D, I) posterior are not noticeably different from before.However, the posterior modes for c2 are 0.96(5) and 3.5(2) at NLO, and NNLO, respectively, with 68% credible intervals indicated in parentheses.Clearly, conditioning on D has a significant impact and increases the EFT truncation error.

B. Posterior predictive distributions
In this section we quantify the posterior predictive distributions (PPDs) for selected np scattering observables at NLO (k = 2) and NNLO (k = 3).A PPD is a distribution of values for a, possibly unseen, observable y ⋆ conditioned on experimental data [63].Specifically, we will sample from pr(y ⋆ |D, ⃗ x ⋆ , I) with y ⋆ = y th (⃗ x ⋆ ) + δy th (⃗ x ⋆ ) at some value(s) ⃗ x ⋆ for the independent variable.Note that we include our estimate of the truncation error in all PPDs.To achieve this, we sample y th,y ) where the covariance matrix is informed by a GP with MAP values ℓ θ , ℓ T lab , c2 according to the values listed in the table in Fig. 1.
As a model check we first inspect the PPD for some of the seen data, i.e., the data already in D. In Fig. 7 we show 500 draws from the PPDs for the vector polarization A y (θ) at NLO with T lab = 40 MeV and T lab = 140 MeV.The model calculations were performed using EC emulators trained as above.The observed correlation length, ℓ θ = 28 degrees, is reasonable, and we hit most of the data within the theory uncertainty.Due to the symmetries of the strong interaction we must have A y (θ = 0) = A y (θ = 180) = 0.As we incorporated this type of constraint during the training the GP, the predictive variance indeed goes to zero for θ = 0, 180 degrees.Given the rather long correlation lengths, the variance constraints at the angular endpoints appears to propagate to the interior to further suppress the truncation error.We have not studied the empirical coverage of the truncation error as we expect it to extrapolate rather poorly to unseen data at higher energies.
In Fig. 8 we quantify PPDs for unseen data and com- Refs.[64,65] and [66], respectively.pare predictions based on correlated and uncorrelated truncation errors.We have drawn 500 samples from the PPDs for the differential cross section σ(θ) at T lab = 319 MeV, the vector polarization A y (θ) at T lab = 350 MeV, and the total cross section σ tot for T lab = 30 − 350 MeV, where the latter is at NLO and the former two at NNLO using correlated as well as uncorrelated truncation errors.
The smoothing effect of the correlated truncation error is clearly visible in panels (a) and (c).For σ(θ), in panels (a) and (b), the predictions in the low-θ region fare significantly worse when including correlations.We speculate that this is connected to the reduced weight of high-energy np data residing within one correlation length of the truncation error.In panels (c) and (d) we compare A y (θ) predictions at NNLO.Besides overall similarities, the predictive variance of A y (θ) is finite at the edges even when we include a correlated truncation error for which we have imposed the relevant symmetry constraints.However, at T lab = 350 MeV we are approximately one correlation length away from the training point where the constraint was imposed, and its effect has deteriorated.Finally, the PPD for the total cross section in panel (e) underestimates the experimental data for T lab > 100 MeV, and the correlated truncation error appears to be somewhat too small as well.

VII. SUMMARY AND OUTLOOK
We quantified a correlated truncation error for ∆-full χEFT predictions of np scattering observables at NLO and NNLO.The correlation structure was modeled using a GP with a squared-exponential kernel, and the resulting MAP values for the correlation lengths in the scattering energy and angle directions are in the ranges of 45-83 MeV and 24-39 degrees, respectively, for the 17 different observable types in the set of np scattering data that we consider.These are significant correlation lengths that, in principle, could reduce the effective number of independent np data in the likelihood by two orders of magnitude and therefore strongly impact the LEC inference.However, other than doubled widths of the univariate marginal posteriors for all LECs we find that the introduction of a correlated EFT truncation error does not change the structure of the LEC posteriors by much.This is explained well by the relatively small marginal variances of the truncation errors that we quantify in this paper.Indeed, we found effective dimensions of the data likelihood that were approximately 1/8 and 1/4 of the length of the experimental np database, at NLO and NNLO respectively.
The marginal variance of the uncorrelated truncation error increases up to four times when jointly sampling it with the LEC values.Thus, future inferences should attempt to marginalize predictions over the hyperparameters of the GP model and if possible also the breakdown scale, or the expansion parameter Q.For this, the HMC algorithm with AD should provide the necessary leverage to enable sampling of a posterior with, at least, doubled dimensionality.There is most likely useful information to learn about the truncation error from theoretical predictions at N3LO in ∆-full χEFT.Adding more chiral orders in the analysis of the truncation error might reveal challenging structures in the order-by-order data that call for more sophisticated GP-kernel design than we employ in this paper.
It will be interesting in future studies to explore the performance of the inferred interactions in predictions of many-nucleon observables.The general ab initio modeling capabilities of ∆-full chiral interactions at NNLO were studied recently in Refs.[15,16,69].These papers targeted, in particular, predictions for heavy nuclei up to 208 Pb and infinite nuclear matter.A common strategy was to first use history matching to identify interaction parametrizations that give acceptable model predictions for a suite of low-energy nuclear observables.Energies and radii of light nuclei, as well as selected np phase shifts, were included sequentially in an implausibility measure that led to the removal of large regions of the LEC parameter space.Although a lower regulator cutoff Λ = 394 MeV was used in the history-matching analyses, it is still interesting to compare the identified regions with the parameter posterior PDF inferred in this paper.First, we find that the posterior mode from this paper projected on c 2,3,4 is situated just outside the nonimplausible range.That is not surprising since the Roy-Steiner prior was employed without updates throughout history matching.Furthermore, the mode is located at the upper end of the non-implausible range in the C 1S0 direction and at the lower edges of the C 1P 1 and C 3P 2 ranges.The significance and origin of this tension should be explored in a joint analysis of N N , few-, and manynucleon observables.
The scattering emulators and HMC sampling presented here can be combined with history matching and reduced-basis methods for many-body observables to enable a joint analysis of chiral interactions.Such studies, however, will rely on an improved understanding of correlated EFT truncation errors and soft momentum scales in finite nuclei.
FIG. 2. (a) Calculated expansion coefficients for σT .(b) The NLO, NNLO, and reference values from which the coefficients are calculated.Note the zero crossing for the reference value at T lab ≈ 205 MeV.
FIG.3.The inverse-χ 2 prior pr(c 2 |I) (orange) for the GP hyperparameter c2 .Shown in black are alternative priors that we use to show the robustness of the inference.

FIG. 5 .
FIG. 5. Correlation matrices based on the NNLO covariance matrix Σσ tot of the 119 σtot data.In panels (a) and (b) we exclude and include, respectively, the experimental variances Σexp,σ tot .The data are sorted by energy in increasing order.

FIG. 6 .
FIG. 6.The EFT truncation and EC emulator errors compared to experimental errors for all considered np scattering data.(a) The standard deviations of the EFT truncation errors at NLO and NNLO.(b) The NNLO EC emulator errors (∆yemu).Colors as in Fig. 2. In both panels we divide by the the experimental 1σ errors (∆yexp) and the horizontal dashed lines indicate a ratio of 1.

FIG. 8 .
FIG.8.500 draws from the PPDs for unseen data for: (a) σ(θ) at NNLO and T lab = 319 MeV with a correlated error, and experimental data from Ref.[67].(b) Same as (a) but with an uncorrelated error.(c) Ay(θ) at NNLO and T lab = 350 MeV with a correlated error, and experimental data from Ref.[68].(d) Same as (c) but with a correlated error.(e) σtot at NLO and T lab = 30-350 MeV (see Ref.[42] for references to the experimental data).

FIG. 9 .FIG. 10 .FIG. 11 .FIG. 12 .
FIG. 9. Posterior PDF of the values for the LECs at NLO inferred using an EFT truncation error in the uncorrelated limit.The inner (outer) gray contour line encloses 39% (86%) of the probability mass.The dot-dashed vertical lines indicate a 68% credibility interval in the univariate marginals.See main text for units and conventions.