Maximally local two-nucleon interactions at N$^3$LO in $\Delta$-less chiral effective field theory

We present new maximally local two-nucleon interactions derived in $\Delta$-less chiral effective field theory up to next-to-next-to-next-to-leading order (N$^3$LO), which include all contact and pion-exchange contributions to the nuclear Hamiltonian up to this order. Our interactions are fit to nucleon-nucleon phase shifts using a Bayesian statistical approach, and explore a wide cutoff range from $0.6-0.9$ fm ($\approx 660 - 440$ MeV). These interactions can be straightforwardly employed in quantum Monte Carlo methods, such as the auxiliary field diffusion Monte Carlo method. Together with local three-nucleon forces, calculations with these new interactions will provide improved benchmarks for the structure of atomic nuclei and serve as crucial input to analyses of astrophysical phenomena of neutron stars, such as binary neutron-star mergers.


I. INTRODUCTION
Astrophysical explosions involving neutron stars (NSs), like supernovae and NS mergers, are fascinating phenomena to study nuclear physics at the extremes.NSs and their mergers explore nuclear matter reaching the highest densities in the cosmos, making them ideal laboratories to elucidate strong nuclear interactions.These interactions manifest themselves in the form of the nuclear matter equation of state (EOS), which connects NSs with experiments of neutron-rich nuclei.Recently, multi-messenger analyses, combining state-of-theart nuclear theory with astrophysical observations, have provided a wealth of new information on the EOS [1][2][3][4][5][6][7][8].These analyses used constraints on the EOS of dense matter that were obtained from many-body calculations using interactions from chiral effective field theory (EFT) [9][10][11][12][13].Chiral EFT provides a systematic expansion of nuclear interactions and is connected to the fundamental theory of strong interactions, quantum chromodynamics [14][15][16].While there has been a lot of progress in recent years to improve chiral EFT constraints, it is crucial to reduce uncertainties of theoretical models for nuclear interactions to fully exploit the multitude of anticipated data from NS observations and nuclear experiments in the coming years.
One way to achieve this goal is to perform calculations at higher orders in the chiral EFT expansion.In this paper, we introduce a new family of local chiral EFT two-nucleon (N N ) interactions at next-to-next-to-nextto-leading order (N 3 LO) that can be employed in quantum Monte Carlo (QMC) computational methods.QMC methods are among the most precise many-body methods [17] but they require local interactions as input, see Refs. [18,19] for details.In the past, local interactions from chiral EFT have been developed up to next-to-nextto-leading order (N 2 LO) in the ∆-less approach [18,19] and including short-range pieces at N 3 LO in the ∆-full approach [20,21].Fully local chiral interactions have also recently been developed in the ∆-less approach where all available local operators up to N 3 LO were considered, even those that would be connected by antisymmetrization and the Fierz rearrangement freedom (FRF) [22].In particular, 4 leading-order (LO) operators, 8 local nextto-leading order (NLO) operators, and all 11 local N 3 LO operators were considered while all nonlocal operators and their associated physics were not included.
In this paper, we develop maximally local interactions for use in QMC methods while using FRF [23].By maximally local, we mean a complete set of operators in which the number of nonlocal operators is as low as possible, consistent with FRF.Local chiral EFT interactions up to N 2 LO, developed following our approach, have been used in QMC methods and provide a good description of nuclei with A ⪅ 20 [12,24,25] and dense matter [12,26,27].However, uncertainties in these calculations are still sizable and result from both the truncation of the chiral expansion and regulator artifacts.These uncertainties affect QMC predictions of nuclei and the nuclear equation of state, including the nuclear symmetry energy and related observables [7,[27][28][29].Forthcoming astrophysics constraints from gravitational-wave (GW) observatories will provide precision data on dense matter [30] and it is key that uncertainties in many-body calculations at low densities are reduced to enable best-possible analyses of these exciting new data.Here, we focus on developing novel maximally local N 3 LO N N interactions that can be used in QMC calculations of various nuclear systems.In a forthcoming paper, we will include the parameterfree N 3 LO three-nucleon (3N ) interactions [31,32], the charge-symmetry and charge-independence breaking corrections, and QMC calculations of many-body nuclear systems.Calculations at higher orders in the EFT ex-pansion are expected to reduce theoretical uncertainties by a factor of 2 − 3 [13,33].
The interactions developed here include a set of 21 contact operators, out of which 4 are nonlocal.All pionexchange interactions are local and fully included.The local interactions developed explore a wide range of cutoffs, R 0 = 0.6 − 0.9 fm (≈ 660 − 440 MeV), in order to reduce the impact of regulator artifacts.This is an important aspect of our calculation since previous studies have shown that regulator artifacts from local regulators are larger than for nonlocal regulators [12,23,34].Highcutoff interactions have smaller regulator artifacts and can easily be employed in QMC calculations, in contrast to most other many-body methods that require softer interactions for convergence.Our N 3 LO interactions, which we name N 3 LO LA -09 to N 3 LO LA -06, are fit to N N scattering phase shifts using Bayesian methods.This allows us to explicitly model EFT truncation uncertainties when performing the fits and we show how these uncertainties evolve with the cutoff R 0 .We also perform least-squares fits to the NN scattering phase shifts that do not incorporate EFT truncation uncertainties.The comparison between the two ways of fitting indicates the importance of modeling the EFT truncation uncertainties when chiral interactions are calibrated to data.We demonstrate that local high-cutoff interactions perform better than their softer counterparts, in the sense that they better reproduce N N scattering phase shifts and lead to smaller EFT truncation uncertainties.Finally, we show that although our interactions are not fit to the properties of the deuteron, our model predictions for these are in good agreement with experimental data, especially for our high-cutoff interactions.This paper is organized as follows.In Sec.II, we give the explicit form of the interactions that we use in this work.The couplings are fit to N N scattering phase shifts and the details of this fit are given in Sec.III, where we discuss both the Bayesian fit as well as the least-squares fit.We also show how the np phase shifts and their associated theoretical uncertainties change with increasing chiral order and varying the cutoff.In Sec.IV, we use our interactions to study the properties of the deuteron.Our main conclusions and summary are presented in Sec.V.

II. N 3 LO INTERACTIONS FROM CHIRAL EFT
In this section, we give the detailed expressions for the interactions along with the local regulators employed in this work.Since there are four nonlocal contact operators, we present the momentum-space expression for all the contacts that we use.The pion-exchange terms and the local regulators, on the other hand, are treated in coordinate space.

A. Chiral EFT and QMC methods
In atomic nuclei and nuclear matter below about twice the nuclear saturation density, chiral EFT is currently the main framework to describe nuclear interactions in a systematic order-by-order expansion [14][15][16].The chiral EFT framework provides consistent N N , 3N , and higher-body interactions, based on a low-momentum expansion of nuclear forces where the expansion parameter Q is the ratio of a typical momentum of the system under study with respect to the breakdown scale Λ b , which determines where chiral EFT becomes inapplicable.The effects of high momenta, that would be resolved above the breakdown scale, are absorbed into a set of low-energy couplings, the strengths of which are adjusted to reproduce experimental data.The advantages of chiral EFT over other approaches are that it (i) allows to quantify theoretical uncertainties [33,35,36] and (ii) provides consistent N N and many-body interactions, i.e., the same processes between different particles are described by the same low-energy constants (LECs) and operators.Order by order, predictions become more accurate and precise by a factor of 2 − 3 at the cost of more involved calculations.Chiral EFT is valid for relative nucleon momenta below ≈ 500 − 600 MeV, translating into densities below about twice the nuclear saturation density.
Solving the nuclear many-body problem is a challenging task that requires advanced computational tools.QMC methods are among the most precise nuclear manybody methods [17] and use stochastic techniques to extract ground-state properties of nuclear systems, providing quasi exact solutions with statistical uncertainties [24].This is the main benefit over other computational methods whose additional approximations can lead to systematic uncertainties.However, QMC methods require local interactions1 as input, i.e., interactions with no derivatives acting on the wave functions, to nonperturbatively solve the nuclear many-body problem.While chiral EFT is traditionally formulated in a nonlocal way and has been used to construct N N interactions to N 3 LO and beyond [37,38], local chiral EFT interactions have been introduced only in the past decade [18,20,22].Local interactions have so far been developed up to N 2 LO in chiral EFT [18,19] on the same footing as nonlocal interactions.In addition, maximally local interactions with selected N 3 LO contributions have been developed [20] but these typically do not include N 3 LO pion-exchange contributions.Finally, recent work saw the development of local N 3 LO interactions where all nonlocalities were neglected and replaced by local operators where possible, even when the latter are connected by FRF [22].Here, we develop complete maximally local N N interactions at N 3 LO in chiral EFT that are suited for QMC calculations, i.e., interactions that account for all short-and long-range contributions and all necessary local and nonlocal pieces up to that order, using FRF.While FRF is violated when local regulators are applied [22,23], this effect induces regulator artifacts that take the form of interaction pieces of higher order in the EFT expansion, and hence, decrease in size.We have shown that the most sizable regulator artifacts at LO can be absorbed by the regular interactions pieces at NLO [23].Similarly, artifacts at NLO and subleading artifacts at LO can be absorbed at N 3 LO.Hence, the remaining artifacts are of order Q 6 and are therefore expected to be small.Furthermore, we study interactions at large cutoffs where these regulator artifacts further decrease in size.For these reasons, regulator artifacts in the NN sector are small in this work.

B. Maximally local intactions at N 3 LO
Chiral EFT interactions are given in terms of a momentum expansion and can be decomposed into shortrange contact pieces and long-range pieces mediated by one and multiple pion exchanges, where ν is the chiral order, indicating the power Q ν .

Contact interactions
Up to N 2 LO, the contact interactions can be fully expressed using only local operators and the nonlocal spin-orbit interaction that can be treated by QMC methods [12,18,19,39].The leading-order (LO) momentumindependent contact interactions are given by where two out of four possible operators are chosen [19,23].The remaining two operators are linearly dependent due to the required antisymmetry of the wave function in nuclear systems.In other words, one can construct an antisymmetrized potential V A , where A is the exchange operator given as By antisymmetrizing the potential in this manner, it can be shown that any two out of four operators, describing both N N S-wave interaction channels, can be selected using FRF, see Refs.[18,19,23] for more details.
At NLO, the contact interaction is momentum dependent.For initial and final relative momenta p and p ′ , momentum transfer q = p ′ − p, and momentum transfer in the exchange channel k = (p ′ + p)/2, the NLO interaction is given by where, again, a subset of 7 independent out of 14 possible operators has been chosen.At NLO, the operators are selected such that the interaction is fully local (i.e., independent of k) except for the spin-orbit interaction.
Other choices are possible, too, which lead to nonlocal or partially local interactions [37,38,[40][41][42].There are no new contact operators that appear at N 2 LO.We now turn to N 3 LO where there are a total of 30 possible contact operators [41], Note, that terms ∼ (q • k) 2 can be expressed in terms of q 2 k 2 and (q × k) 2 , and hence, do not have to be included explicitly.Also, the last two operators are related to a quadratic spin-orbit operator and the angular momentum operator by Due to FRF, we can again choose a subset of 15 independent operators from this over-complete set.When selecting the subset, however, it is crucial to choose operators in such a way that upon antisymmetrization the complete operator set is recovered.Hence, for example, a nonlocal tensor operator cannot be replaced by a central local interaction piece as physics information at that order would be lost.One possible way to choose an independent set of operators would be to drop all terms in Eq. ( 6) that contain the product of isospin matrices τ 1 • τ 2 , as done in Ref. [41].Alternatively, in this work, we choose the following subset of operators: The We can further reduce the number of nonlocal contacts as it has been found that there are redundancies among the 15 contact operators at N 3 LO [38,43].By performing a unitary transformation (UT) on the Hamiltonian, we can decrease the number of independent nonlocal contacts from 7 to 4 [38].Following Ref. [38], we consider the unitary operator where T i are the three anti-Hermitian generators of the UT and γ i are the corresponding transformation angles.
Similarly to Ref. [38], we choose the following generators: Note that Eqs.(10), (11), and (12) constitute one choice of basis for the UT and other choices are possible due to FRF.For example, Ref. [38] replaces the operator τ 1 • τ 2 by σ 1 • σ 2 in Eq. ( 11) and, τ 1 • τ 2 with 1 1 in Eq. (12).Moreover, apart from the freedom in choosing the generator basis, there are only three two-nucleon contact operators that are anti-Hermitian and Galilean invariant at order Q 2 .Therefore, the unitary operator U considered in this work represents all possible UTs that can be generated at order Q 2 .We need to apply the UT only to the LO Hamiltonian since the UT, when applied to the higher-order interactions, will induce terms at order Q 6 and above which is beyond the desired accuracy for this work.The shift in the LO Hamiltonian is given by where the dots represent terms beyond the order of our calculation.It can be shown that [V cont , T i ] only induce shifts to contact operators of order Q 0 and Q 2 and therefore do not need to be considered explicitly [38].On the other hand, the commutator with the kinetic energy generates order Q 4 terms (see Ref. [38]): Using the identity and Eq. ( 10) of Ref. [38], we can express Eqs.( 13) and ( 14) as Note that some of the new operators, are not explicitly chosen in Eq. ( 8), but they are linearly dependent operators; see Eq. ( 6).Therefore, we do not need to consider them explicitly any further.
Having carried out the UT and including the shift to the Hamiltonian, the N3 LO contact interaction is now The variables γ 1 , γ 2 and γ 3 are completely arbitrary parameters of the UT and can be chosen to remove nonlocal operators.The parameter γ 1 can be chosen to be either −Λ 2 Therefore, we see that using the UT we can remove 3, leaving only 4 nonlocal operators.
In this paper, we set The nonlocal part of the N 3 LO contact interactions then consist of only 4 operators: where we have used the relative orbital angular momentum operator L = (q × k).We will refer to our maximally local N 3 LO interactions with this choice of 4 nonlocal operators by N 3 LO LA -09, N 3 LO LA -08, N 3 LO LA -07, and N 3 LO LA -06, where the number refers to the cutoff in coordinate space.The coordinate space expressions of the N 3 LO contact operators are given in Appendix A and the matrix elements of these operators in the partial wave basis is given in Appendix B. We have also considered other possible choices for the set of 4 nonlocal operators-see Appendix C-but found the set chosen in this work to be best suited for QMC methods because three operators can be directly mapped into the 18 operator channels of the phenomenological Argonne V18 (AV18) interaction [45] that has been used extensively in QMC simulations.
In summary, our N 3 LO LA interactions contain 21 contacts in total, out of which 4 are nonlocal.The corresponding 21 LECs are determined by fits to np phase shifts, see Sec.III for details.In future work, we will also include the isospin-breaking contributions to the EFT interaction, following the strategy of Refs.[18,19] to calculate neutron matter for the physical scattering length at all orders.

Pion-exchange interactions
The long-range and intermediate-range parts of chiral EFT interactions are mediated by pion exchanges.All pion-exchange interactions to N 3 LO are either fully local or accompanied by the spin-orbit operator, and thus, we directly give the coordinate space expressions here.
Without loss of generality, the pion-exchange contributions can be decomposed as [22,41] where is the tensor operator and S = (σ 1 + σ 2 )/2 is the total spin operator.At LO, only the one-pion exchange (OPE) contributes.It is given by [18,19,22]  where x = m π r, m π is the pion mass, g A is the axialvector coupling constant, and F π is the pion decay constant.For the constants m π , F π , g A , and the nucleon mass M N , we use the same values as in Ref. [19].Here, we use the charge-independence breaking form of the OPE as we have done before; see Ref. [19] for details.At NLO and beyond, V π receives contribution from two-pion exchange (TPE) diagrams.For these TPE pieces, we employ the expressions using spectral-function regularization (SFR) [18,22,41].Under this representation, the TPE potential is written in terms of spectral functions as and similarly for W C,S,T,LS .All TPE contributions to the Hamiltonian at NLO, N 2 LO, and N 3 LO can therefore be specified in terms of the spectral functions V C,S,T,LS (iµ).The explicit expressions for these spectral functions can be found in Appendix A of Ref. [22].We note that we add to V (4) π the 1/M N corrections to the N 2 LO TPE potential.These corrections nominally appear at fifth order in chiral EFT but, as is common  [18,19] whereas the couplings at N 3 LO are taken from Ref. [46].
practice [22,41], we include it at N 3 LO in order to arrive at an improved intermediate-range attraction.Also, note that unlike Ref. [22], we do not take the limit Λ → ∞ but leave the SFR cutoff finite.In this work, the SFR cutoff Λ is fixed at the value of 1 GeV, similarly to Refs.[18,41], and we have verified that changing it to Λ = 2 GeV does not significantly affect our results; see Appendix D.
The strength of the TPE potential is determined by the πN couplings which are constrained by chiral symmetry and can be determined by analyses of low-energy πNscattering.Here, we employ the values obtained from a Roy-Steiner (RS) analysis of πN scattering [46] at N 3 LO.At N 2 LO, we chose the same values as in Ref. [18,19] for consistency with our previous interactions.These values are very close to the values extracted from the RS analysis at N 2 LO except for c 4 .The values of the couplings employed in this work are given in Table I.R 0 = 0.9 fm FIG. 2. Phase shifts in the 1 S0 and 3 S1-3 D1 partial waves.In panels (a)-(d), the chiral EFT order is fixed to be N 3 LO and we show results for different cutoffs as indicated in the legend.In panels (e)-(h) [panels (i)-(l)], we show results for different chiral EFT orders at R0 = 0.6 fm [R0 = 0.9 fm].The bands correspond to the 95% CL.The solid lines represent the solutions that maximize the Bayesian posterior distributions while the dashed lines represent the least-squares fit results.For comparison, we show the Nijmegen partial-wave analysis (NPWA) [44] (black squares) and the phase shifts of the AV18 potential [45] (blue stars), which has been used in past QMC calculations.

C. Regulators
Chiral EFT interactions need to be regulated at short distances and/or high momenta: To define maximally local interactions, we choose regulators that are fully local.We choose Gaussian regulators in position space, that can be expressed as with n = n 1 = 2, n 2 = 6.We therefore, have one position-space cutoff R 0 that regulates both the shortand long-range pieces of the interaction.Upon Fourier transformation (FT) of the regulator functions, the position-space cutoff can be related to the momentumspace cutoff Λ c as R 0 = 2/Λ c .In this paper, we study four different cutoffs, R 0 = 0.9, 0.8, 0.7, and 0.6 fm, leading to four different interactions, N 3 LO LA -09 to N 3 LO LA -06.

III. ANALYSIS OF NN SCATTERING
In this work, the 21 operator LECs are determined by fits to np phase shifts, while in the future we will explore fits to scattering data.We perform the fits directly in momentum space using the formalism developed in Ref. [47].This allows us to fit interactions that include nonlocal pieces, which is crucial at N 3 LO.We take the phase-shift values from the Nijmegen partial-wave analysis (NPWA) [44], and incorporate EFT truncation uncertainties by performing Bayesian fits to these data.Bayes' theorem defines the posterior P as where Π is the prior distribution on the LECs, L is the likelihood function that incorporates information from the phase shifts, and the normalization constant, i.e., the evidence Z can, in principle, be used to perform model comparison, but we will not pursue this here.We take the prior Π to be uniform everywhere in parameter space, because the LECs can vary strongly in size as we vary the cutoff.We require our likelihood function to incorporate a model for the EFT truncation uncertainties.In order to do so, we assume that the EFT expansion holds for any scattering observable X(p), i.e., where the expansion parameter Q ≡ max(mπ,p) min(Λ b ,Λc) .Here, p is the relative momentum, Λ b is the breakdown scale which we take to be 600 MeV [35,38], and Λ c is the momentum-space cutoff.Note that in the Weinberg power counting scheme [41], c 1 = 0.For sufficiently large k, one might expect that the uncertainty stemming from the truncation of the EFT expansion at order k can be approximated by the contribution to the observable X at order k + 1, i.e., However, since c k+1 is unknown, Epelbaum, Krebs, and Meißner (EKM) [35] proposed approximating Eq. ( 35) as where j = max(2, k + 1).Note that this expression provides an estimate of the systematic truncation uncertainty but it does not have a strict probabilistic interpretation.However, it makes minimal assumptions regarding the nature of c n and the underlying probability distribution from which the c n are drawn.On the other hand, if one assumes that the expansion coefficients are drawn from a Gaussian distribution, the BUQEYE collaboration showed that the EFT truncation uncertainty, when summed to infinite order, also follows a Gaussian distribution whose variance is given as [43] where c2 is the root mean square of the expansion coefficients c n up to order k.Equation (37) holds in the limit where all correlations across different kinematic points are ignored, see Ref. [43] for the corresponding expression in the fully correlated limit.Assuming that experimental uncertainties are also Gaussian distributed, the BUQEYE approach leads to a Gaussian model for the likelihood function, where σ 2 i,BUQEYE = σ 2 i,exp + σ 2 i,theo,BUQEYE .The product over i indicates a product over all considered partial waves and kinematic variables (laboratory energies).The variable X denotes an observable which, in this case, is the phase shift for a given laboratory energy and partial wave.This likelihood function, along with priors that enforce naturalness of the LECs, was extensively used in Bayesian analyses of phase shifts in Ref. [43].
The BUQEYE prescription Eq. ( 38) holds under the assumption that the expansion coefficients c n are Gaussian distributed.The systematic uncertainties induced by truncating the EFT expansion does not necessarily need to follow such a distribution.Therefore, in this work, we turn to the EKM prescription (36) that makes fewer assumptions regarding the expansion coefficients c n , i.e., it assumes a uniform distribution of these parameters.As mentioned earlier the EKM expression (36) provides only an estimate of the EFT truncation uncertainty and does not have a robust probabilistic interpretation, but see Ref. [48] for a mapping of the EKM error bars to the BUQEYE statistical distributions.Therefore, as an ansatz, we choose ∆X k EKM to set the scale of our likelihood function, i.e., we take σ theo,EKM = α∆X k EKM , where α is a constant.We also add experimental uncertainties in quadrature to the theoretical truncation uncertainties, completing our ansatz for the variance of our likelihood function, We emphasize that Eq. ( 39) is merely a model for the variance of our likelihood function, i.e., we do not claim to have derived it and no assumptions have been made regarding the distribution followed by the expansion coefficients c n .We now turn to the form of the likelihood function.For this, we use the principle of maximum entropy which states that, among all real-valued functions with a specific variance σ 2 , the function that maximizes the differential entropy is a normal distribution with that variance.This fixes our likelihood function to be, where σ i is given by Eq. ( 39).As mentioned above, in this work we take the EFT truncation uncertainty to be σ i,theo,EKM = α∆X k EKM .In the limit σ 2 i,theo ≫ σ 2 i,exp (which we have verified for all i), the parameter α turns out to be only an overall constant in the log-likelihood function.In the case of uninformative uniform priors on the LECs, the choice of α, therefore, has negligible impact in our Bayesian analysis, eliminating the impact of differences of the EKM confidence interval at different orders.For simplicity, we take α = 1 and we have checked that changing it to α = 1/2 has no impact on our results.The experimental uncertainties σ i,exp are taken to be the uncertainties provided by the NPWA [44].
The likelihood employed in this work, Eq. ( 40), and that in Eq. ( 38) are similar, with the only difference being the contribution of the EFT truncation uncertainty to the variance of the Gaussian likelihood, i.e., σ theo,BUQEYE vs σ theo,EKM .In Appendix E, we compare these two quantities and show that these two choices are almost identical, leading to a similar Bayesian fit of our interactions.With this justification, we will henceforth focus on the posterior defined using Eq. ( 40) and leave a more detailed comparison with the BUQEYE approach for future work.Finally, we note that our approach to modeling theoretical truncation uncertainties does not incorporate the correlation of these uncertainties across different kinematic points [43,49].In Ref. [50], it has been shown that incorporating such correlations doubles the variance of the marginal posteriors of all LECs, but it leaves the structure of the LEC posteriors unchanged.Here, we thus limit ourselves to the uncorrelated case, and will perform a more sophisticated analysis that includes these correlations in the future.
For the fits at N 2 LO and N 3 LO, we set X LO = 0 because we have otherwise found very poor Bayesian fits at high laboratory energies, with many samples including resonances or spurious bound states.The reason for this is that, due to the large differences between the data and the LO predictions at high laboratory energies, the EKM uncertainty estimates are dominated by the poor LO predictions and are very large, removing the constraining power of any higher-energy data points and leading to spurious structures.Hence, we treat the LO contribution as an outlier in the expansion; see also Ref. [50] for a similar treatment of the LO predictions.
Evaluating the EKM uncertainty, Eq. ( 36), at a given order requires knowledge of the phase shifts at all lower orders.Therefore, as an initial step to performing Bayesian fits, we first determine the LECs via a leastsquares minimization of the objective function, where m is the number of experimental data points included in the fit.At LO, the least-squares optimizations are done by fitting to the 1 S 0 and 3 S 1 partial waves up to a laboratory energy of E max = 50 MeV.At NLO and N 2 LO, we fit the 1 S 0 , 3 S 1 , ϵ 1 , 1 P 1 , 3 P 0 , 3 P 1 , and 3 P 2 partial waves up to E max = 150 MeV.For these orders, we fit to phase shift values at energies specified in Ref. [19].At N 3 LO, we additionally include 3 D 1 , ϵ 2 , 1 D 2 , 3 D 2 , and 3 D 3 in the fit and we fit up to E max = 250 MeV, additionally including points at 200 and 250 MeV. 3 The results of these fits serve as an order-by-order estimate of the EFT convergence that is used to estimate the EKM uncertainty Eq. ( 36) which, in turn, is used as an input for our Bayesian fits.Note that the χ 2 function used in the least-squares fit does not incorporate theoretical EFT 3 For the least-squares fit at R 0 = 0.7 fm alone, we chose Emax = 350 MeV in order to remove a spurious resonance in the 1 P 1 channel at E ≈ 350 MeV.We verified that this change of Emax has a negligible effect on the phase shifts below 350 MeV in all other channels.
truncation uncertainties.Therefore, the least-squares fits also serve as complementary analyses to the Bayesian fits and can be used to estimate the importance of modeling EFT truncation uncertainties when chiral EFT interactions are calibrated to scattering phase shifts.In Fig. 1, the left panel shows the local component of our N 3 LO interactions in the 1 S 0 channel, with the LECs determined via the above mentioned least-squares fits.The interaction becomes increasingly hard for smaller R 0 because the LO 1 S 0 spectral LEC C1S0 increases rapidly with decreasing R 0 ; see right panel of Fig. 1.This growth is required to counter the attractive pionexchange potential that becomes increasingly singular for lower coordinate-space cutoffs [23].
The Bayesian analyses are carried out by a Markovchain Monte Carlo (MCMC) sampling of the posterior distribution defined in Eq. (33).The details of the Bayesian fits, i.e., the partial waves involved and E max at a given order, are the same as in the least-squares fit.The MCMC sampling results in a large number of samples drawn from the posterior distribution.We used the emcee Python package [51] to carry out the MCMC sampling and we used 5000 walkers.To check the convergence of the MCMC sampling, i.e., to ensure that sufficient iterations were performed, we imposed that the total number of sampler iterations was larger than 50 times the autocorrelation time for all the sampling parameters, i.e., for all the LECs.We also checked that our estimate of the autocorrelation time, as computed by the emcee package, remained stable (to within 1%) as a function of the sampler iteration.Once the MCMC sampling was terminated based on these convergence criteria, we discarded the first 2τ max samples, where τ max is the largest autocorrelation time across all sampling parameters, to make sure that we use only the "burnt-in" samples.We also thinned our sample chains by 0.5τ min to remove any spurious correlations among our samples.Finally, we found that a small percentage (< 5%) of the samples at N 3 LO contained spurious resonances.We discarded these samples and verified that this did not modify the posterior distribution significantly.In this manner, we obtained a posterior distribution over the entire parameter space, which at N 3 LO is a parameter space spanned by the 21 LECs.This posterior can be converted to posteriors over any two-body observable, such as np phase shifts.
In Fig. 2, we show our results for the np phase shifts in the 1 S 0 and 3 S 1 -3 D 1 partial waves for both the leastsquares fit (dashed lines) and the Bayesian fits (solid lines that represent the maximum of the posterior and the bands that correspond to the 95% confidence level (CL) given the definition of uncertainty discussed above).In the top row, we demonstrate the effect of varying the cutoff R 0 at N 3 LO, i.e., we show results for our interactions N 3 LO LA -09 to N 3 LO LA -06.In the middle (bottom) row, we vary the chiral EFT order with the cutoff fixed at R 0 = 0.6 fm (R 0 = 0.9 fm).Figs. 3 and 4 show similar results but for P and D waves respectively.In general, regardless of the cutoff, we find good convergence order-by-order, with our results approaching the NPWA data with reduced uncertainty bands when going to higher orders.By comparing results for different cutoffs, we see that the softer interactions deviate significantly from the data at higher laboratory energies in several partial waves, especially the D waves.Several comments are in order.First, note the large deviations of the LO phase shifts from the data in almost all channels.This effect enhances the contribution of the LO phase shifts to the EKM uncertainties, which is why it is important to set X LO = 0 in Eq. (36) for the fits at N 2 LO and N 3 LO.Second, it is important to keep in mind that the uncertainty bands shown here, reflecting the Bayesian posteriors obtained in the fit, describe the LEC uncertainties and are not necessarily the same as the bands obtained following the EKM prescription a posteriori, as was done in Ref. [35].This can be understood by considering an LO interaction where the EKM uncertainties widen with the laboratory energy due to the almost constant phase shift.However, parts of such an EKM uncertainty band would not be accessible in an actual fit.For example, describing phase shifts that decrease with the laboratory energy, as provided by the lower bound of an EKM uncertainty band, requires effective-range contributions which are only provided at NLO.Therefore, the results of a Bayesian fit will differ and show a less dramatic energy dependence.For this reason, we do not show the 95% CL bands for the LO results, and we plot only the least-squares and the maximum posterior results (dashed and solid lines, respectively).Finally, our interactions at lower orders are not fit to the D waves but the N 3 LO interactions are.As a consequence, the differences between N 2 LO and N 3 LO bands in the D waves are larger than the differences between NLO and N 2 LO.
The impact of the inclusion of the EFT truncation uncertainties that guide the Bayesian fits can be gauged by comparing the solid and dashed lines in all panels in Figs.2-4.Especially at lower cutoffs, the impact of the EFT truncation uncertainties is important and the Bayesian analysis yields a better reproduction of the data at lower energies as compared to the least-squares fits.At higher energies, the least-squares fits perform better generally.Note however that, since we do a global fit involving several partial waves, there are some channels in which the Bayesian fit gives better results at higher energies, e.g., for the mixing angle ϵ 1 for our N 3 LO LA interactions.
In general, we also find that the high-cutoff interactions give a better reproduction of the NPWA data than low-cutoff interactions, see the top rows in Figs. 2, 3 and 4 where the phase shifts at N 3 LO are compared for different cutoffs.The Bayesian fits clearly demonstrate the reduction of uncertainty when decreasing R 0 , which is equivalent to increasing the momentum-space cutoff Λ c .In order to stress this point, in Fig. 5 we show the objective function defined in Eq. ( 41) divided by the number of fit parameters (21 at N 3 LO) as a function of R 0 .We see the strong decrease in χ 2 as interactions become harder, before plateauing and eventually increasing at R 0 ≈ 0.55 fm.
The LECs that we obtain for our N 3 LO LA interactions are given in Table II for both fitting strategies: leastsquares fits and Bayesian fits (maximum posterior estimate).Similar results for the LECs at LO, NLO, and N 2 LO are given in Appendix F.

IV. DEUTERON
Finally, we study properties of the deuteron, which can be completely described by the N N interactions developed here.In this work, we do not fit the interactions to reproduce any properties of the deuteron, and instead choose experimental data on deuteron properties as a benchmark for our interactions.
In Fig. 6, we show the deuteron binding energy as a function of cutoff, for different orders.The bands represent the Bayesian posteriors at the 95% CL whereas the least-squares results are shown as dashed lines.We see that, at N 3 LO, the least-squares results which do not include any theoretical uncertainty estimates are compatible with the predictions obtained from the Bayesian fits.However, for NLO and N 2 LO, there is a clear difference between the two analyses.This, once again, illustrates the importance of incorporating theoretical uncertainties when EFTs are calibrated to experimental data.Also, note that the uncertainties decrease significantly with the chiral order as well as when decreasing R 0 , which is consistent with our results for the phase shifts.Our N 3 LO LA interactions reproduce the experimental deuteron binding energy within theoretical uncertainties for all considered cutoffs. .The results for different R0 correspond to the least-squares phase shifts analyses.The AV18 [45] analysis is shown as a reference.
In Fig. 7, we show results for the deuteron wave function.In the top row, we show the coordinate space wavefunction in the S wave, ψ L=0 (r), and D wave, ψ L=2 (r).Note that these wavefunctions are related to their components u(r) and w(r), which are sometimes reported in the literature [20,38,41], as ψ L=0 (r) = u(r)/r and ψ L=2 (r) = w(r)/r.In the bottom row, we show the momentum-space representation of the wavefunctions in the S wave, ψL=0 (p), and D wave, ψL=2 (p), showing the squared momentum-space wavefunctions in a logarithmic scale.As our fitting code is written in momentum space, we first solve for the momentum-space wavefunctions and then obtain the coordinate-space representations as, In coordinate space, our models do not contain any oscillations in the wave function at intermediate or large r, unlike some chiral interactions developed in the past [40,41].Interactions at different cutoffs are identical at large r but have different small-r behavior.In the D wave, the harder interactions, i.e., N 3 LO LA -06 and N 3 LO LA -07, lead to more pronounced peaks as compared to their softer counterparts N 3 LO LA -08 and N 3 LO LA -09.
We have used the deuteron wave functions at N 3 LO to calculate other deuteron observables, which we show in Table III.We give results for both our Bayesian fits (at the 95% CL) and least-squares fits.For the Bayesian fits, our results are in good agreement with the experimental data within theoretical uncertainties, except for the quadrupole moment which is to be expected since we do not take electromagnetic two-body currents into account.Note that the agreement between the median model predictions and the experiment is better for interactions with smaller coordinate space cutoffs.Furthermore, these harder interactions have smaller uncertainties, similar to what was observed in the phase shift analyses.The least-squares results are also in good agreement with the experimental data, and the agreement improves for harder interactions.We note that, since we fit our interactions to the phase shifts in the  II.LECs for our N 3 LOLA interactions for different cutoffs.We give results obtained from the Bayesian analyses (left) and least-squares fits (right).For the former, the quoted LECs correspond to the maximum of the posterior distribution.
it is not surprising that we get reasonable deuteron properties.However, this calculation still serves as a valuable model check.Finally, note that our interactions result in D-state probabilities P d ⪆ 6%, which is larger that what has been found for other local N 3 LO interactions [20,22].

V. SUMMARY AND OUTLOOK
We have constructed maximally local interactions at N 3 LO in ∆-less chiral EFT.Our interactions include a total of 21 contact operators at N 3 LO, out of which four are nonlocal.We have also included all pion-exchange contributions consistently up to N 3 LO.We have employed the method of Bayesian statistics in order to calibrate our maximally local interactions to np phase shifts, explicitly taking into account EFT truncation uncertainties.We have demonstrated the importance of incorporating these uncertainties in the fit.Finally, we have explored a range of cutoffs that is significantly larger than what is typically explored in the literature, showing explicitly that high-cutoff interactions lead to smaller EFT trun-cation errors and a better reproduction of experimental data.Furthermore, it is expected that these interactions significantly reduce uncertainties due to the violation of the Fierz rearrangement freedom, particularly from 3N interactions [12].We will demonstrate this explicitly in a forthcoming work.
Future QMC calculations using the interactions developed in this work will allow us to reduce theoretical uncertainties for nuclei and dense nuclear matter, which is important to improve nuclear-physics inputs to nuclear structure and reactions, and for analyses of multimessenger observations of neutron stars and related simulations.
where we have used the spectroscopic LECs that are related to the operator LECS by

4 b D 9
or Λ 4 b D 11 , removing either of the two corresponding contact operators.Similarly, γ 2 can be chosen as either −Λ 4 b D 10 or Λ 4 b D 12 .Finally, γ 3 can be set to −Λ 4 b D 14 .

FIG. 1 .
FIG. 1. (a):The local part of our N 3 LOLA interactions in the 1 S0 channel (i.e., 17 local contacts and all the pion-exchange pieces) for the least-squares fit as a function of particle separation r.The LECs are obtained via least-squares fits as described in Sec.III.We show the interactions for different cutoff values, i.e., N 3 LOLA-09 to N 3 LOLA-06.(b): Cutoff dependence of the LO spectral LEC in the 1 S0 partial wave at N 3 LO.The solid blue line shows the least-squares results whereas the band represents the 95% confidence level (CL) calculated using a Bayesian analysis; see Sec.III.

FIG. 5 .
FIG. 5. Objective function defined in Eq. (41) divided by the number of fit parameters N (21 at N 3 LO) as a function of cutoff for our N 3 LOLA interactions.The calculated values are written explicitly in the figure.

FIG. 6 .
FIG.6.Deuteron binding energy as a function of cutoff for different orders, as indicated in the legend.The bands correspond to the 95% CL obtained from the Bayesian analysis whereas the dashed lines represent the least-squares results.The dashed black line shows the experimental binding energy.

FIG. 7 .
FIG. 7. Deuteron wave function in coordinate space [panels (a) and (b)] and momentum space [panels (c) and (d)].The results for different R0 correspond to the least-squares phase shifts analyses.The AV18[45] analysis is shown as a reference.

4 FIG. 8 . 3 FIG. 9 .
FIG. 8. Differences between the phase shifts obtained using the full interaction and the predictions with all nonlocal parts set to zero.We show results for the N 3 LOLA-06 interaction as well as the interactions defined in Eqs.(C1) to (C4), i.e., Set 1 to Set 4.

FIG. 10 .
FIG.10.Comparison of the BUQEYE estimate of the EFT truncation uncertainties with the EKM estimate in different S-wave channels.This serves as input to our Bayesian analysis-see Eq. (40)-and similar error estimates will result in similar Bayesian fit posteriors.The blue distribution is obtained using Eq.(E1) and the solid black lines are the corresponding 68% CL.The EKM uncertainty, i.e., Eq. (36) with k = 4, is shown as dashed red lines.Different panels correspond to different partial waves and laboratory energies as written in the upper left of each panel.
operators D 1 to D 8 are local and D 9 to D 15 are nonlocal.Note that the terms D 5 and D 6 , while being nonlocal, can be treated by Monte Carlo methods since the momentum k appears only linearly.Hence we group them together with the local operators at N 3 LO, just as we have grouped C 5 along with the other 7 local operators at N 2 LO.

TABLE I
. πN couplings used in this work.The couplings at N 2 LO are taken from Refs.

TABLE III .
[56]erties of the deuteron for our N3LOLA interactions, for different cutoffs.For the Bayesian fit, the error bars are quoted at the 95% CL.Here, E d is the binding energy, Q d is the quadrupole moment, η d is the asymptotic D/S ratio, ⟨r 2 ⟩ d m is the root-mean-square matter radius, AS is the asymptotic S-wave normalization and P d is the D-state probability.The experimental value for E d is from Ref.[52], Q d is from Ref.[53], η d is from Ref.[54], ⟨r 2 ⟩ d m is from Ref.[55]and AS is from Ref.[56].