Bayesian Analysis of χ EFT at Leading Order in a Modified Weinberg Power Counting Approach

,


I. INTRODUCTION
The foundational principles of chiral effective field theory (χEFT) [1][2][3] promise a systematically improvable description of the strong force between nucleons that is consistent with quantum chromodynamics (QCD).However, establishing a power counting of chiral nuclear interactions that fulfills renormalization requirements presents several challenges-see, e.g., van Kolck [4] for an overview.Recently, Yang et al. [5] analyzed nuclear ground-state energies using renormalization-group (RG) invariant formulations of χEFT up to (perturbative) next-to-leading order (NLO) corrections.Apparently, the essential mechanism for nuclear-binding tends to fail already at leading order (LO) for selected A > 4 nuclei when using available RG-invariant power counting schemes.In fact, similar results had already been found using χEFT [6] based on the canonical Weinberg power counting (WPC) [7][8][9][10], and pionless EFT [11][12][13].Yang et al. [5] put forward three possible reasons for these shortcomings at LO in χEFT: i) One (or more) scales critical to the physical description of finite nuclei might not be correctly captured by the contact terms at LO [14]; ii) The LO nucleon-nucleon (N N ) interaction should possibly be complemented with other interaction terms such as sub-leading pion exchange [15] and manynucleon interactions [16][17][18][19]; iii) The description of the nuclear interaction might be finely tuned and therefore require careful calibration of the low energy constants (LECs).Indeed, Yang et al. [5] renormalized the relevant LECs by demanding exact reproduction of selected phase shifts at a single scattering energy.This procedure resulted in point estimates of the LECs without the possibility to analyze possible fine-tuning effects.* oliver.thim@chalmers.se In this work we tackle overfitting and expose possible fine tuning using Bayesian methods.Specifically, we infer a posterior probability density function (pdf) for the values of the LECs conditioned on neutron-proton (np) scattering data.The advantages of a Bayesian approach are several.First, we obtain a probabilistic measure of our uncertainty about the values of the LECs and subsequent predictions, something that is not obtained when doing maximum likelihood estimation [6].Second, with the Bayesian approach we can utilize the expected systematicity of χEFT as prior knowledge about the truncation error [20].There exist several Bayesian studies of χEFT interactions [21][22][23][24][25] and quantified posterior predictive distributions (ppds) of nuclear properties [26][27][28].However, so far all such studies are grounded in χEFT based on WPC.
In EFTs that are perturbative at all orders, power counting usually follows the momentum scaling of the different Feynman diagrams.This is known as naive dimensional analysis (NDA) [29].However, in χEFT we must account for bound multi-nucleon states and therefore face non-perturbative physics with the consequence of infrared enhancement of purely nucleonic intermediate states.To deal with this, Weinberg [2,3] suggested to apply NDA to the potential which is then iterated to all orders by solving, e.g., the Lippmann-Schwinger equation.This prescription assumes that the iteration to infinite order does not introduce additional divergences with the need for higher-order counterterms-an assumption that did not hold upon closer inspection.Indeed, taking the momentum cutoff Λ of the regulator very large, i.e., far beyond the anticipated breakdown scale of χEFT, Nogga et al. [30] found that WPC does not yield RG-invariant amplitudes in attractive spin-triplet partial-waves.By now there exist several proposals on how to modify WPC, see, e.g., Refs.[5,18,[31][32][33][34][35][36][37].We refer to all such proposals as modified Weinberg power counting (MWPC).One can argue that WPC provides a consistent EFT framework as long as Λ is kept in the vicinity of the breakdown scale Λ χ [38][39][40][41].In that scheme, all orders are typically summed up non-perturbatively and, starting at third order, it can provide realistic predictions for selected nuclei [42] and nuclear matter [43].However, this achievement does not imply that WPC provides the foundation for an EFT of QCD.
This paper is organized as follows: In Sec.II the LO potential and non-relativistic np scattering are addressed along with relevant background regarding power counting and the renormalization of singular potentials.In Sec.III we discuss our Bayesian approach and in Sec.IV we review the numerical sampling of the posterior pdf for the LECs.Finally, in Sec.V ppds are presented and analyzed for some np-scattering observables and the deuteron ground-state energy.A conclusion and outlook is presented in Sec.VI.

II. THEORY
To keep this work self-contained we first discuss how counterterms (and associated LECs) are introduced to renormalize the singular nature of the one-pion exchange (OPE) potential at LO in χEFT.We also study limitcycle-like behavior in some detail since they play an important role during the Bayesian inference.

A. Effective field theory expansion
The EFT approach promises an order-by-order improvable description of a nuclear observable y residing in the low-energy domain below the relevant breakdown scale.In general, up to some finite order n we can expect an expansion of the form [10,44] and we specialize to χEFT assuming a breakdown scale Λ χ = 600 MeV as per previous analyses [21,45] and denote the momentum cutoff by Λ.We also assume a soft scale Q = max(p, m π ) with p denoting the external momentum and m π denoting the pion mass.The unspecified function D n+1 depends on ratios of the relevant low energy scales and, importantly, absorbs the residual cutoff dependence.In this work we set the reference scale y ref to the corresponding experimental value y exp of y.Pulling out y ref leads to dimensionless expansion coefficients b k , and if the theory is renormalized order-by-order, these should not exhibit any cutoff dependence although they do depend on ratios of relevant scales.Also, along the lines of naturalness [46], we expect the b k -values to be of order unity.
For a perturbative EFT-where NDA applied to the Feynman diagrams carries over to the amplitudes and hence to observables-the power counting in Eq. ( 1) follows straightforwardly from the Lagrangian.The relation between the Lagrangian and observables in a nonperturbative theory like χEFT is less direct.The nonperturbative effects in combination with the need to treat the problem numerically pose challenges to finding a consistent power counting.When the amplitudes cannot be obtained perturbatively, the LO T -matrix must scale as [44,47,48].The summation in Eq. ( 1) can still start at k = 0, since the Q −1 dependence can be absorbed into y ref .Terms with k = 1, 2, . . .then correspond to higher-order corrections with respect to LO.We employ MWPC where amplitudes of the LO potential are computed non-perturbatively and sub-leading corrections should be accounted for using perturbation theory [4,33,[49][50][51].

B. Two-nucleon potential and scattering amplitudes
The momentum-space and isospin-symmetric LO potential considered in this work is adopted from Ref. [5], and using our conventions it reads The first term is the OPE potential where τ i and σ i , for i = 1, 2, are the isospin and spin operators for the respective nucleon, p (p ′ ) are the ingoing (outgoing) relative momenta with normalization ⟨p ′ |p⟩ = δ 3 (p ′ − p), q = p ′ − p is the momentum transfer, and we have p ≡ |p| and p ′ ≡ |p ′ |.The contact LECs: C1 S0 , C3 S1 , C3 P0 , C3 P2 , carry implicit projection operators such that they act in the indicated partial-waves, (2s+1) l j , where s, l, j denote the quantum numbers of the N N spin, orbital angular momentum, and total angular momentum, respectively.We use the PDG [52] values for the axial coupling constant g A = 1.275, the pion decay constant f π = 92.1 MeV, and averaged pion mass m π = 138.039MeV.Finally, for the partial-wave decomposed [53] potential, and operators, we use the notation As discussed in Sec.II C and Ref. [5], this LO potential is restricted to low partial waves: Consequently, it is set to zero for all partial waves with l > 1 that have no coupling to l ≤ 1.
To render the integrals of the Lippmann-Schwinger equation finite, the relative momenta of in-and out-going nucleons are regulated using the function, which was also used by Yang et al. [5] 1 , where We straightforwardly account for the small effects of relativistic kinematics using the minimal relativity prescription [9,54].Combined with the momentum regulators f Λ , the LO potential in Eq. ( 2) is thus modified according to where mp+mn is the nucleon mass with proton and neutron masses m p = 938.272MeV and m n = 939.918MeV, respectively [52].
We condition all inferences on np scattering data and must therefore solve the corresponding Lippmann-Schwinger equation for the T -matrix We solve this numerically in momentum space using a standard matrix-inversion method by first converting to a real equation for the reaction matrix [55].Furthermore, we sum the partial-wave amplitudes to construct the spin scattering matrix, M s m ′ s ms (p, θ c.m. ), see Appendix A. Here, m s , (m ′ s ) is the in-(out-)going total spin projection and θ c.m. is the center-of-mass scattering angle.All scattering observables can be computed from this spin scattering matrix [56].The types of observables that we condition our inferences on are listed in Table I and discussed in Sec.III A.

C. Singular potentials and limit-cycle-like behavior
An attractive potential is singular near the origin if it behaves as −λ/r n with n > 2 (or n = 2 with sufficiently large λ > 0) [57,58].Two particles interacting only via a singular potential will collapse towards the origin with increasing velocity.Akin to the infinities in quantum field theory, the singularity at the origin is cured via renormalization.The OPE potential is singular in attractive spin-triplet partial waves, e.g., 3 P 0 , 3 S 1 − 3 D 1 , 3 D 2 and 3 P 2 − 3 F 2 [30].From an EFT perspective, the singular character of the OPE becomes physically meaningful only with the addition of counterterms that parameterize, and cure, these short-range pathologies [59,60].
To avoid the introduction of infinitely many counterterms, the OPE potential in Eq. ( 2) is truncated to act only in partial waves with orbital angular momentum l < l c .It is not obvious, however, where the limit l c should be set.Higher-order terms can be included using the distorted-wave Born approximation and there is evidence that this does not necessitate the introduction of additional counterterms, see Refs.[37,41] for contrasting views.According to previous studies, OPE contributions to the N N scattering amplitude in partial waves with orbital angular momentum l > 1 can be treated perturbatively up to at least p ∼ 300 MeV [4,35,61].We therefore truncate the OPE potential such that it is nonzero only for channels that contain a partial wave component with l ≤ l c = 1.
LECs associated with counterterms that renormalize a singular potential with n > 2 usually exhibit a limitcycle-like behavior [30,59,60], i.e., provided that there is only one counterterm per partial wave, the corresponding LEC will exhibit periodic discontinuities as a function of the regulator value.This was extensively investigated for one-pion exchange in Ref. [30].We can reproduce those results exactly, but we also observe a slight shift in the location of the limit-cycles-like behavior when including the minimal relativity correction in Eq. ( 5).In Fig. 1 we show how the LECs in the 1 S 0 , 3 S 1 , 3 P 0 and 3 P 2 partial waves run with Λ when renormalized to reproduce the Nijmegen partial-wave phase shifts [62] for the laboratory kinetic energy T lab = 50 MeV of the projectile nucleon.There is no limit-cycle-like behavior observed in the 3 P 2 wave for the cutoff region studied here.
When this renormalization procedure is used, spurious and deeply bound states appear in the singular and attractive partial-waves as the cutoff Λ is increased.This is not a problem if the states are deeply bound and thus outside the applicable domain of the EFT.In practice, one can project these states out of the spectrum in a phase-equivalent fashion [30].In Fig. 2 we show how the spurious states in 3 P 0 and 3 S 1 − 3 D 1 appear at the threshold value of Λ corresponding to where the limit-cycle-like behavior is observed in Fig. 1.
In this work we infer the LEC posteriors at fixed values of the cutoff Λ.There exist threshold LEC values for which the total potential becomes sufficiently attractive such that a spurious bound state appears.The (Λdependent) LEC values for which this happens are tabulated in Appendix B for the respective channels.Note that a spurious state is not deeply bound in the immediate vicinity of these threshold values and we stress that the behavior depicted in Fig. 2 is obtained when requiring the exact reproduction of a specific phase shift.
As the bound state pole moves near the threshold when the corresponding LEC is varied, the phase shift will change dramatically by Levinson's theorem [63].Rapidly varying phase shifts lead to particular challenges in the inference process since scattering observables that are included in the likelihood, and that depend sensitively on partial-wave amplitudes, might then constrain the Bound-state energies in selected channels with the LECs renormalized to reproduce the phase shift from the Nijmegen partial-wave analysis [62] at T lab = 50 MeV.There are two bound states in the 3 S1− 3 D1 channel, where the upper one (orange dots) is the physical deuteron state while the lower one (blue dots) is spurious.The latter, deeply-bound state appears at Λ ≥ 1140 MeV.In the 3 P0 channel there is a spurious, deeply-bound state for Λ ≥ 770 MeV.
LEC values to very narrow regions of parameter space.In Fig. 3 we exemplify this by showing phase shifts at T lab = 50 MeV as a function of the LEC in the partial waves 3 P 0 and 3 P 2 for different cutoff values.The red, horizontal line indicates the empirical value of the phase shift according to the Nijmegen analysis [62].Fig. 3 also shows that for certain cutoffs we can obtain a wide range of LEC values that reproduce the empirical phase shift with reasonable precision.
Figures 1 and 3 provide complementary views on limitcycles-like behavior.For example, in Fig. 1 we see that C3 P0 changes from a large positive value to a large negative one as we approach Λ ≈ 780 MeV from below.The same information is contained in Fig. 3 where the LO 3 P 0 phase shift intersects with the Nijmegen result for a positive (negative) value of C3 P0 for Λ = 700 (1000) MeV.Enforcing the exact reproduction of a phase shift there- fore implies that the spurious bound state will be deeply bound exactly when the limit-cycle-like behavior appears as shown in Fig. 1.For the 3 P 2 wave, the phase shift only intersects with the Nijmegen result on one side of the discontinuity, at least for the low values of Λ shown here.Consequently, there is no spurious bound state, or limit-cycle-like behavior appearing.

III. BAYESIAN PARAMETER ESTIMATION
Bayes' theorem expresses the posterior pdf, pr (θ|D, I), for the relevant model parameters, θ, conditioned on experimental data, D, and other pertinent information, I, in terms of quantities that we can evaluate pr (θ|D, I) = pr (D|θ, I) × pr (θ|I) pr (D|I) .
Here, pr (θ|I) is the prior, pr (D|θ, I) is the likelihood, and pr (D|I) is the model evidence.The latter is independent of θ and can be omitted for parameter estimation purposes leaving us with the simpler expression pr (θ|D, I) ∝ pr (D|θ, I) × pr (θ|I) .
The set of model parameters in this study is collected in a vector where α denotes the LECs of χEFT at LO and c governs the magnitude of the EFT truncation error, as described below in Eq. ( 17).The elements of the vector α at LO are in units of 10 4 GeV −2 and 10 4 GeV −4 for the S-and P -wave LECs, respectively.

A. Likelihood
When relating theory y th and experiment y exp for some observable y we must account for both the theoretical and experimental uncertainties.This can be achieved with a statistical model in which the uncertainties are modeled by random variables δy th and δy exp , respectively, and where we suppressed all parameter dependencies for simplicity.One of the great benefits of working with EFTs is the expected systematicity of the truncation error.Assuming that QCD accounts for all relevant physics, and that our EFT is sound, i.e., it systematically approximates QCD in the low-energy regime, then we can express the value of some nuclear observable y as in the notation of Eq. ( 1).Thus, the LO theoretical prediction can be expressed as where the second term contains the residual cutoff dependence through D 1 .For the present case we assume that b 1 = 0 by time-reversal and parity invariance [8][9][10].
With this, Eq. ( 12) can be written as The coefficients c k are related to b k and the power series expansion of D 1 .Since b 1 = 0, it implies that c 1 = 0 if the first term in the power-series expansion of D 1 vanishes, which we assume it does by the same symmetry argument as for b 1 .The coefficients c k , like b k , are also assumed to be of natural size, i.e., of order one.Finally, we note that c k might have a residual cutoff dependence inherited from D 1 .
We can now express the LO truncation error as Quantitatively, this is unknown to us since we have not computed anything beyond LO.However, we can evaluate the probability distribution of the random variable δy th if we assume a distribution for the expansion coefficients c k guided by domain knowledge included in I.In this work we follow Ref. [20] and assume that where N (µ, σ 2 ) denotes a normal distribution with mean µ and variance σ 2 .Assuming independent and identically distributed c k coefficients the distribution for δy th at each kinematical point is given by [64] pr (δy th |I) = N 0, σ 2 th , Hence, the variance of the truncation error is governed by c2 , which is unknown a priori.However, operating with an EFT we expect c to be of order one.We will quantify this in Sec.III B.
The experimental error for each datum is assumed to follow a normal distribution with variance σ 2 exp , i.e., pr (δy Assuming that the EFT truncation error and the experimental errors are independent the likelihood for a single observation reads by Eq. ( 11), where we now explicitly indicate the relevant parameter dependencies.Following Refs.[24,45] we further assume that EFT errors of scattering observables at different energies and angles are independent.The experimental data set, D = {y exp,i } N i=1 , on which we condition the MCMC inference is listed in Table I and contains N = 1043 scattering observables.This is not all the data in the Granada database [65,66].We attempted an inference using all np observables with T lab < 100 MeV, but a model check revealed that most calibration data were poorly reproduced.For example, low-energy total cross sections deviated significantly from experimental values, which in turn yielded significantly overbound deuteron states for most cutoffs.This is probably due to a misspecified EFT truncation error, leading to unnaturally large values for c and thus overestimated EFT errors for low-energy cross sections (that are expected to be reproduced relatively well at LO).There are several suggestions [22] for how to construct more sophisticated EFT error models that would accommodate the incorporation of more data in the likelihood.We postpone such developments until we have incorporated sub-leading orders in MWPC.

B. Priors
We treat α and c as random variables and must therefore place priors on them.For the LECs we employ independent and normally distributed priors with a standard deviation of four times their naturalness estimates: the S-and P -wave LECs, respectively [9,10,39].Similarly, we expect c to be of natural size and always positive.As shown in Fig. 4, we employ an inverse gamma distribution with parameters α = 3 and β = 4.2.This distribution has a mode around one, but also a heavy tail that allows for a significant variation.

IV. SAMPLING THE POSTERIOR AT DIFFERENT CUTOFF VALUES
The parameter inference proceeds in two steps.First, we employ a Bayes linear approach known as history matching [28,[67][68][69] to identify relevant domains for the LECs.Second, we employ an affine invariant MCMC algorithm called emcee [70] to numerically draw samples from the posterior pdf.The history-matching domains identified in the first step are used to initialize the MCMC sampling.Furthermore, to avoid detrimental influence of the limit-cycle-like behavior during sampling we select momentum cutoffs based on the analysis in Sec.II C. TABLE I. Total and differential np scattering observables used to condition the Bayesian inference in this work.The total number of data points is 1043.This is a subset of the data in the Granada database [65,66], where the normalizations determined from the Granada analysis are included.The uncertainties in the normalizations are not taken into account.The definitions of the various scattering observables can be found in Ref. [56].

Observable
T History matching is an iterative scheme that can be employed to identify and exclude regions of parameter space that produce model outputs inconsistent with observational data, taking relevant uncertainties into account.See, e.g., Ref. [28] for details of the algorithm and an application in nuclear physics.At the end, history matching returns the region(s) of parameter space that could not be ruled out, the so-called non-implausible domain.
In this work, however, we employ history matching as a precursor to a full Bayesian analysis.A similar use of history matching can be found in Ref. [71].Our aim is then to identify regions of parameter space where we expect to find large contributions to the posterior probability mass.We will be less concerned with the observational property of the data that we utilize in this step (we will employ scattering phase shifts), or with the rigour of our uncertainty estimates, as the sole purpose is to pick promising starting points for the MCMC algorithm.For this reason we refer to excluded samples as inefficient starting points, whereas the final domain is at least non-inefficient.These designations replace the standard implausible and non-implausible labels that are found in the history-matching literature.
We employ an "inefficiency" measure I M (α) that gauges the performance of the theoretical predictions for a selected subset of data given values for the LECs The index, i, enumerates the set of experimental observations that is included in the history matching.The measure I M (α) is defined as the maximum over the observations {y exp,i } and is therefore governed by the datum that is reproduced the worst.The experimental and theoretical errors are incorporated via the variances σ 2 exp,i and σ 2 th,i , respectively.A second measure, I 2M (α), is sometimes used to construct an additional constraint.It is analogously defined as the second maximum over the observations and obviously fulfills I 2M (α) ≤ I M (α).
We perform two waves of history matching.In the first wave, the scattering phase shifts listed in Table II are considered as observational data.Since all LECs act in distinct partial waves, each LEC could be constrained independently within one-dimensional sub-waves.We estimate conservative phase-shift errors σ th from Eq. ( 17) using fixed c ≈ 6.This choice implies a 30% error for T lab < 40.6 MeV (i.e., p < m π ) and 80% for phase shifts with T lab ≈ 100 MeV.For each sub-wave we use 10 4 different LEC values in a space-filling Latin hypercube design [72] across a rather wide interval informed by the phase-shift analysis in Sec.II C and corresponding to the ranges shown on the y-axes in Fig. 5.Samples for which the inefficiency measures are larger than some cutoff val-ues are deemed as poor candidates for initializing the MCMC algorithm.In this analysis we use a sequence of cutoffs, I M > 3.0 combined with I 2M > 2.5, which is similar to Ref. [69] (therein based on Pukelsheim's three sigma rule [73]).In the end, the relevant parameter volume is reduced by a large fraction.The ratio of the noninefficient over initial volume was found to be as small as 10 −5 , although this ratio depends rather strongly on the value of the cutoff.
In the second wave, a set of 13 np scattering observables is used to construct the inefficiency measures.In this wave all four LECs are active simultaneously.Here we employ 10 5 samples, using the same space-filling design, in the non-inefficient domain resulting from the first wave.The larger number of samples, and the much reduced volume helps to provide sufficient resolution for detecting narrow domains.The set of observables and corresponding model uncertainties, estimated using the same prescription as in the first wave, can be found in Table III.
The resulting, relevant ranges for each of the LECs are shown in Fig. 5.The domains are classified, and colorcoded, according to if the LEC gives rise to a spurious (unphysical) bound state in the respective channel.The rapid variation of the phase shift by 180 • where a bound state appears-discussed and illustrated in connection to Fig. 3-gives rise to disjoint domains.Therefore, we can anticipate multimodal structures of the LEC posterior pdf.
Depending on the cutoff, we find spurious states in all channels.In the 1 S 0 channel there is always a region with an unphysical bound state, but it becomes more narrow for increasing values of the cutoff.The physical deuteron state is always present in 3 S 1 − 3 D 1 , and we observe additional (spurious) states in disjoint domains for several cutoffs.For Λ ≥ 3000 MeV there are two disjoint domains that contain at least one spurious bound state each.
While there are no spurious bound states in 3 P 0 for Λ ≤ 500 MeV, a second domain with lower values of the LEC contains a spurious bound state for Λ = 550, 600, 700 MeV.For larger cutoffs, Λ ≥ 2000 MeV, TABLE II.np scattering phase shifts in degrees, used as observational data within the first wave of history matching.Values for T lab are given in MeV.The phase shift at T lab = 40 MeV was obtained through spline interpolation.All other phase shifts were taken from Refs.[65,66].The standard deviation for the theoretical error in Eq. ( 21) is presented in the last column as a relative error.We can understand the Λ-dependence of the P -wave domains by studying Fig. 3 in light of the errors going into the inefficiency measure (21).In doing so we see that for 3 P 0 and Λ = 550 − 700 MeV the theoretical and empirical phase shifts match across wide regions of LEC values that qualitatively agree with the non-inefficient domains observed in Fig. 5.A similar argument holds for 3 P 2 and greater values of Λ.
Unfortunately, history-matching fails to reduce the initial domain of C3 P0 and C3 P2 for some momentum cutoffs and thus provides limited information in these cases.In the following we perform a Bayesian analysis.The LEC posteriors presented in Section IV B are conditioned on observational data (see Table I) and a more reliable error model.They provide more informative and credible inference compared to the domains shown in Fig. 5.

B. Markov chain Monte Carlo sampling
We sample the posterior pdf in Eq. ( 8) using the affine-invariant ensamble sampler from the python package emcee [70] conditioned on the experimental data in Table I.As noted in Sec.IV A, there are disconnected non-inefficient regions for some LECs at several values of the cutoff.These regions might correspond to multimodal structures of the posteriors.In infinite time, the MCMC walkers will explore all parameter space, but in finite time they are likely trapped in local modes.To handle such convergence challenges we initiate sampling in several of the regions identified using history matching.Ideally, all combinations of such regions should be explored in order to locate the dominant posterior mode(s) with some certainty.
We initialize 50 MCMC walkers using a random subset  of points from the non-inefficient domains for each cutoff and let them take 5000 steps each after a burn-in period of 1000 steps.This is typically sufficient for obtaining a good representation of each posterior mode of our fivedimensional posterior pdf.The initialization of walkers in the P -waves are chosen in the region of the domain without bound states, shown in orange in Fig. 5, except for the higher cutoffs in 3 P 0 .The same principle is applied in the 3 S 1 partial wave, where we initialize in the region containing the least amount of spurious states.The posterior turns out to be relatively flat in the direction of P -wave LECs.Thus, the MCMC sampling is less sensitive to where the walkers are initialized in the rather large start domains.Moving forward, we focus on the possible multimodalities originating from the 1 S 0 channel where we always found a parameter region giving rise to a shallow and unphysical bound state.The initialization of MCMC walkers in the 1 S 0 channel is chosen either in the unbound (orange) or bound (blue) region.Proceeding like this for selected cutoff values between 400 MeV and 4000 MeV we indeed find multimodal structures in the posterior pdf, as summarized in Table VI in Appendix C.
To analyze the probability mass gathered in each mode we calculate the model evidence Z K for all domains K enclosing a mode.We do this using the Laplace approximation, which is accurate if the probability distribution can be locally approximated with a multivariate Gaussian, i.e., where Σ −1 is the Hessian matrix of pr (θ|D, I) evaluated at the mode θ * .Each mode listed in Table VI in Appendix C is indeed found to be well approximated with a Gaussian.The numerical evaluation of Eq. ( 22) might be sensitive to the convergence of the MCMC chains.However, this sensitivity mainly impacts the evaluation of the determinant of the inverse Hessian which turns out to be a less significant contribution to the evidence compared to the density at the MAP point [74].The ln(Z K ) values for all investigated modes are summarized in Table VI of Appendix C. Fortunately, a well-defined "best" mode with highest evidence can be identified at all cutoffs except for Λ = 700 and Λ = 4000 MeV for which there are two P -wave modes with very similar evidences.All the significant modes identified in the evidence analysis are retained, and none of them contain an unphysical bound state in the 1 S 0 channel.
To ensure that we obtain a high-quality representation of the relevant posteriors we perform the final sampling using 50 walkers taking 5000 steps each initialized in the vicinity of the largest ln(Z K ) mode for each cutoff.This lays the foundation for accurate inferences of the ppds presented in Sec.V. We use an autocorrelation analysis [24] and visual inspection of the traces to confirm MCMC convergence.At the cutoffs Λ = 700 MeV and Λ = 4000 MeV the two P -wave modes with similar evidences remain.The resulting values for the model evidences and MAP locations at several cutoffs are shown in Table IV.The model evidences are of similar magnitude for all cutoffs.In fact, the biggest difference, on a logarithmic scale, is about 40, i.e., the maximal ratio of evidences is e 40 .The same difference can, for example, be obtained by shifting all predictions by about 5%, which is quite small compared to the EFT truncation error.
In Figs.6a and 6b, we show the marginal posterior pdfs for the cutoff values Λ = 450 MeV and Λ = 700 MeV.The posterior pdfs for the remaining cutoff values are shown in Appendix C. The posteriors are consistent with our naturalness assumptions for both the LECs and the scale of the EFT truncation error as no tensions with our priors are observed.
The running of the inferred LECs can now be monitored.The cutoff dependence of the marginal posteriors for the LECs is shown in Fig. 7.We conclude that C1 S0 (Λ) and C3 S1 (Λ) have a similar behavior as the running of those couplings determined from the phase shift fits in Sec.II C and Fig. 1.This finding can likely be attributed to the fact that these LECs are quite well constrained by the employed low-energy np scattering cross sections.However, the P -wave LECs C3 P0 and C3 P2 are not that well constrained by the np data used in the inference and run differently with Λ in the Bayesian analysis.
The inference of c is quite interesting in its own right.In Fig. 7 we show its marginal posterior and see that it is rather insensitive to cutoff variation, in particular TABLE IV.The model evidences and maximum a posteriori (MAP) estimates (θ * ) of the LECs for all cutoffs Λ (in MeV).A posterior bimodality with comparable probability mass per mode exists for both Λ = 700 and Λ = 4000 MeV.

V. POSTERIOR PREDICTIVE DISTRIBUTIONS
In this section we present ppds for scattering phase shifts as well as selected np observables.The ppd for an observable y is a marginalization of the model predictions over the posterior of the relevant parameters pr (y|D, I) = dθ pr (y|θ, D, I) pr (θ|D, I) .
When sampling the model prediction, pr (y|θ, D, I), we can choose to include the EFT truncation error, in which case it takes the form pr (y|θ, or we can exclude it, which corresponds to just propagating the parametric LEC uncertainty, by using a delta distribution pr (y|θ, D, I) = δ y − y The integral in Eq. ( 23) can be straightforwardly evaluated using the samples from the MCMC chains obtained in Sec.IV.In the following, ppds including (excluding) the EFT error are colored blue (purple).

A. Phase shifts
The Bayesian inference in this work differs from the point estimates obtained from a standard phase shift optimization.In Fig. 8 we compare the ppds from the present analysis with the optimized phase from Yang et al. [5] for Λ = 450 MeV and the Nijmegen partialwave analysis [62].The EFT error is not included in the ppd.One finds that the parametric uncertainty of the phase shifts, stemming from the LEC posterior, is rather small and the 95% credible intervals shown in the figure are mostly visible in the 3 P 2 partial wave.We find good agreement with the results by Yang et al. [5] except for the P -wave phase shifts which are more repulsive in the Bayesian analysis.Apparently, the P -wave LECs receive significant contributions when conditioning on scattering data.At Λ = 450 MeV, The P -wave phase shifts 1 P 1 and 3 P 1 , which are not part of the contact potential, are more attractive than the ones from the Nijmegen analysis, see Fig. 9.Note that the LO potential only includes low partial waves, for which the OPE is not perturbative.To prevent overfitting to higher order effects, which are better described by high partial waves, either the LO potential or the model for the LO EFT error likely needs Note that the EFT error is not included.Our results are compared to the phase shifts of Yang et al. [5] and the Nijmegen partialwave analysis [62].The labels Y40 and Y200 indicate the two types of fits that were done by Yang et al. [5] in the 3 P2− 3 F2 channel; renormalizing the corresponding LECs to reproduce the 3 P2 phase shifts at T lab = 40 MeV and T lab = 200 MeV, respectively.
to be revisited.We also foresee that the LO LECs will receive corrections at higher orders and this could significantly change the scattering amplitudes and the pole structure in a given channel.In Fig. 10 we compare experimental data for several np scattering observables and the corresponding ppds, including the EFT truncation error, at LO and using the cutoff Λ = 450 MeV.Total cross sections, σ, and differential cross sections, dσ/dΩ, are reproduced quite well, although the ppd error bands are rather wide.Using Our results are compared to the results from the Nijmegen partial-wave analysis [62].
Eq. ( 17) and c ≈ 2.8 from Fig. 6a we estimate that the EFT truncation error for p ≲ m π is about 15%.The ppd for the spin-polarization observable P b , which was not included in the MCMC sampling, reproduces experimental data rather poorly.This result is particularly striking in the ppd at T lab = 25.0MeV.The situation improves at higher energies, which is somewhat surprising for a low-energy EFT.As such, conditioning on low-energy P b data could improve the model.However, expanding the dataset in Table I requires an improved model for the EFT truncation error [22] and likely access to order-byorder calculations beyond LO to guide the inference of c.
In Fig. 11 we show the ppds for σ and P b as a function of Λ, including the EFT error, for a handful of values and θ c.m. .Except for a of P b for Λ < 700, no significant cutoff variation is observed.It can also be seen that the experimental value is well reproduced for σ and more poorly for P b .The LO EFT error for σ appears to be the conservative side and might shrink as we learn it better.
The ppd, without the EFT error, for the deuteron ground-state energy is shown in Fig. 12 and reveals an approximate 0.5 MeV underbinding, but no significant cutoff dependence.For some of the higher cutoffs additional spurious and deeply bound states appear at E ≈ −1 GeV.Employing the same error model as for the scattering observables one can attempt to make an uncertainty estimate for the deuteron ground-state Using the binding momentum p gs = m N (−E gs ) ≈ 50 MeV < m π as a proxy for the relevant soft scale we arrive at an estimated EFT error of around 15% of the predicted value.Including this error makes the lower end of the ppds touch the experimental value.
The ppds for the energies of the unphysical bound states in the channels 3 P 0 and 3 P 2 − 3 F 2 are shown in Fig. 13, except for a very deep (E ≈ −1 TeV) state in 3 P 0 for Λ = 4000 MeV.Some of the energies are around 25 MeV, which is well within the validity of the EFT.Note also that the ppds in Fig. 13 show a strong dependence on the cutoff, which is generally the case for spurious states, which was also seen in Sec.II C. In the 3 P 2 − 3 F 2 channel for the cutoff Λ = 700 the deeply bound state at E ≈ −1 GeV is produced by the very small tail of negative C3 P2 -values shown in the marginal pdf in Fig. 7.We see spurious states in all singular channels for Λ ≥ 2000 MeV.

VI. CONCLUSIONS AND OUTLOOK
In this work we have presented a detailed Bayesian study of a LO χEFT potential in MWPC.We have performed robust inferences of the LECs and the scale c of the EFT truncation error using history matching and MCMC sampling for several values of the momentum cutoff in the range from 400 to 4000 MeV.Below, we summarize the main takeaways.
• Multimodal LEC posteriors are induced via the existence of bound states, virtual states and resonances that move around near threshold for varying values of the LECs and cutoff Λ.This makes inference challenging.In particular, we could always find modes where the 1 S 0 state was bound at an energy of approximately 1 MeV.• We computed the model evidence for each mode in the Laplace approximation.At each cutoff we inia final round of MCMC sampling around the mode(s) with the largest evidence (see Table IV).Modes harboring a bound 1 S 0 state were always found to be insignificant.However, pairs of modes in the 3 P 2 − 3 F 2 channel for Λ = 700 MeV and in the 3 P 0 channel for Λ = 4000 MeV had similar evidences and could not be discerned.We expect that the poor constraints on the LECs might be improved with a better error model and with the inclusion of more types of scattering observables in the calibration data.
• Our inference shows that c ≈ 2.8 across all cutoffs investigated.Since the main posterior modes have similar evidences (see Table IV) this implies that the description of the low-energy scattering data is essentially equivalent for the different momentum cutoffs in the range 400 to 4000 MeV.
• Conditioning the LEC inference on scattering data leads to excessive repulsion in the 3 P 0 and 3 P 2 partial waves compared to the Nijmegen analysis [62] and to the phase-shift optimization performed by Yang et al. [5].Despite the poor accuracy of P -wave phase shifts (see Fig. 8), we find that scattering observables are reproduced surprisingly well (see Figs. 10 to 12) for a LO potential.We conclude that it might be problematic to use phase shifts for calibration of LO potentials in MWPC.These potentials, acting in the first few partial waves only, need to compensate for excluded higher-order contributions to reproduce observables.Consequently, the inference becomes crucially dependent on the specification of the EFT error model, of which we have limited information.This highlights the need for further studies at higher orders.
• We conclude that the ppds for various bound-state and scattering observables in the np sector exhibit no significant cutoff dependence.
It remains to be seen what happens as sub-leading orders are included perturbatively and when the model for the EFT truncation error is refined.The inferred LO potentials, with broad distributions for the P -wave LECs, will likely lead to significant predictive uncertainties in nuclei.Therefore, the ppds for ground-state energies and radii should also be quantified to assess the quality and physical relevance of existing RG-invariant χEFT formulations.for the S-and P -waves respectively.The median and the 68% equal-tailed credible interval are indicated for the univariate marginal pdfs.

4 GeV − 4 ]FIG. 1 .
FIG.1.LECs renormalized to reproduce the phase shift from the Nijmegen partial-wave analysis[62] at T lab = 50 MeV as a function of the cutoff Λ.The limit-cycle-like behavior in C3 S 1 is seen at Λ = 1150 MeV and in C3 P 0 at Λ = 780 MeV.

FIG. 3 .
FIG.3.The phase shifts δ3 P 0 and δ3 P 2 as a function of the respective LEC, for four different cutoff values Λ = 450, 550, 700, 1000 MeV at T lab = 50 MeV.The horizontal (red solid) line is the empirical value of the phase shift[62].

2 FIG. 5 .
FIG.5.The ranges of relevant starting points for each LEC, resulting from two waves of history matching for several cutoff values.The units used for the LECs are defined below Eq.(10).Non-inefficient ranges are separated and color-coded according to if the channel contains spurious bound states.The range of each LEC axis indicates the region that was searched in the first wave of history matching.

FIG. 6 .
FIG.6.Posterior pdfs for the parameters θ = (α, c) using cutoffs Λ = 450 MeV and Λ = 700 MeV.The units of the LECs are 10 4 GeV −2 and 10 4 GeV −4 for the S-and P -waves respectively.The median and the 68% equal-tailed credible interval are indicated for the univariate marginal pdfs.

FIG. 7 .
FIG.7.The marginal posterior pdfs for the parameters θ for cutoffs, Λ, in the range 400 to 4000 MeV.The units of the LECs are 10 4 GeV −2 and 10 4 GeV −4 the S-and P -waves, respectively.

FIG. 8 .
FIG.8.Predicted phase shifts (in degrees) with Λ = 450 MeV represented by the median of the ppd (solid purple line) and the 95% equal-tailed credible interval (purple band).Note that the EFT error is not included.Our results are compared to the phase shifts of Yang et al.[5] and the Nijmegen partialwave analysis[62].The labels Y40 and Y200 indicate the two types of fits that were done by Yang et al.[5] in the 3 P2− 3 F2 channel; renormalizing the corresponding LECs to reproduce the 3 P2 phase shifts at T lab = 40 MeV and T lab = 200 MeV, respectively.

FIG. 9 .
FIG.9.Predicted phase shifts (in degrees) in the P -waves not part of the contact potential using the cutoff Λ = 450 MeV.Our results are compared to the results from the Nijmegen partial-wave analysis[62].

FIG. 10 .
FIG.10.Median (solid blue line) of the ppds for selected np scattering observables using MWPC at LO with a cutoff Λ = 450 MeV.The EFT error is included.Shaded dark and light blue bands represent 68% and 95% equal-tailed credible intervals, respectively.Experimental data (exp) from Refs.[65,66].

10 ( 5 ( 4000 FIG. 11 . 3 S 1 − 3 D 1 FIG. 12 .
FIG.11.The ppds for the scattering observables σ and P b at two different lab energies as a function of the cutoff, Λ, and including the EFT error.The red horizontal line and band represent the experimental value and error, respectively[65,66].

2 − 3 F 2 Λ
FIG.13.The ppds for the energy of selected bound states in the 3 P0 (upper row) and 3 P2− 3 F2 (lower row) channels as a function of the cutoff, Λ, and excluding the EFT error.Note the different scales on the y-axes.The red horizontal line indicates the binding threshold.

FIG. 15 .FIG. 16 .
FIG.15.Posterior pdfs for the parameters θ = (α, c) at different cutoffs.The units of the LECs are 10 4 GeV −2 and 10 4 GeV −4 for the S-and P -waves respectively.The median and the 68% equal-tailed credible interval are indicated for the univariate marginal pdfs..
there is only one narrow domain and it gives rise to a spurious bound state.In 3 P 2 we observe larger domains and spurious bound states for Λ ≥ 2000 MeV.

TABLE III .
np scattering observables used within the second wave of history matching.Note that σ, without any subscripts, is denoting the total np cross section, while σ th /yexp is the (relative) model error for each observable.
The np phase shifts for the MAP estimate of the LECs for Λ = 450 MeV are tabulated in Table VII in Appendix D.
B. Observables in the np sector