Bayesian estimation of the low-energy constants up to fourth order in the nucleon-nucleon sector of chiral effective field theory

We use Bayesian methods and Hamiltonian Monte Carlo (HMC) sampling to infer the posterior probability density function (PDF) for the low-energy constants (LECs) up to next-to-next-to-next- to-leading order (N3LO) in a chiral effective field theory ($\chi$EFT) description of the nucleon-nucleon interaction. In a first step, we condition the inference on neutron-proton and proton-proton scattering data and account for uncorrelated $\chi$EFT truncation errors. We demonstrate how to successfully sample the 31-dimensional space of LECs at N3LO using a revised HMC inference protocol. In a second step we extend the analysis by means of importance sampling and an empirical determination of the neutron-neutron scattering length to infer the posterior PDF for the leading charge-dependent contact LEC in the $^{1}S_0$ neutron-neutron interaction channel. While doing so we account for the $\chi$EFT truncation error via a conjugate prior. We use the resulting posterior PDF to sample the posterior predictive distributions for the effective range parameters in the $^{1}S_0$ wave as well as the strengths of charge-symmetry breaking and charge-independence breaking. We conclude that empirical point-estimate results of isospin breaking in the $^{1}S_0$ channel are consistent with the PDFs obtained in our Bayesian analysis and that, when accounting for $\chi$EFT truncation errors, one must go to next-to-next-to-leading order to confidently detect isospin breaking effects.


I. INTRODUCTION
In chiral effective field theory (χEFT) [1][2][3][4] the nuclear interaction is parametrized in terms of low-energy constants (LECs) that capture unresolved physics and must be determined from data.The number of LECs grows with the order of the chiral expansion, and the parameter estimation becomes a challenging inference problem that is directly connected with the precision of the theory.Moreover, in χEFT the long-range pion-nucleon (πN ) interaction appears as a subprocess of the nuclear interaction.We can therefore constrain the long-range part of the nuclear interaction rather well using measured πN scattering data [5], albeit less so in the delta-full sector of χEFT [6].In a Bayesian context, this knowledge can be straightforwardly accounted for as a prior when learning more about the nuclear interaction from new data, as demonstrated in, e.g., Refs.[7,8].In fact, Bayesian inference methods allow us to account for any prior beliefs about χEFT, most importantly its truncation error [9].This highlights some of the central synergies of χEFT and Bayesian methods.In this paper, we focus on two challenges tied to a Bayesian analysis of the nucleon-nucleon (N N ) interaction: (i) how to reliably sample the highdimensional Bayesian posterior probability density functions (PDFs) of the LECs up to next-to-next-to-next-toleading order (N3LO) in χEFT, and (ii) how to extend these posterior PDFs, here inferred from neutron-proton (np) and proton-proton (pp) scattering data, to account for the uncertainty in the isospin breaking (IB) and leading neutron-neutron (nn) short-range LEC acting in the 1 S 0 partial wave.Throughout the paper we will often * isak.svensson@chalmers.seuse the short-hand labels posterior and prior to indicate posterior and prior PDFs, respectively.
Sampling a high-dimensional PDF poses a significant challenge in any Bayesian analysis.In a previous paper [8] we explicitly demonstrated the efficiency and manageable dimensional scaling of the Hamiltonian Monte Carlo (HMC) [10] algorithm applied to χEFT up to nextto-next-to-leading order (NNLO).HMC exploits the geometry of the parameter space and Hamiltonian dynamics to draw virtually independent Markov chain Monte Carlo (MCMC) samples concentrated to the bulk of the probability mass.In this paper, we revise our HMC protocol to sample the LEC posterior at the next chiral order, i.e., N3LO, where we encounter a significantly more challenging inferential problem in 31 dimensions, one per LEC.
Furthermore, the fundamental effect of IB leads to slight differences in the strong interaction between neutrons, protons, and between protons and neutrons.It originates from differences in the masses and electromagnetic charges of the up-and down-quarks [11] and is expected to be weaker in the irreducible interaction between three nucleons [12,13].Although weak, the IB of the strong nuclear interaction plays an important role in ab initio, and mean-field, analyses of infinite nuclear matter and finite nuclei, in particular towards the driplines with pronounced proton-to-neutron ratios (see, e.g., Refs.[14][15][16][17][18][19][20][21][22][23]).It is therefore important to quantify the uncertainties pertaining to the LECs governing the strengths of the IB effects.Unfortunately, inferences conditioned solely on the world database of N N scattering data [24,25], which does not contain nn cross sections, leaves the posterior of LECs acting (only) in the nn isospin channel unconstrained.To handle this, a point estimate of the leading charge-dependent LEC in the nn channel is typically obtained using the empirical value for the corresponding scattering length, primarily in the 1 S 0 partial wave.This latter quantity parametrizes the total nn cross section in the limit of zero scattering energy and its value is estimated from data on hadronic reactions that involve two neutrons in the initial and/or final state [26].A future prospect is to employ results from realistic lattice quantum chromodynamics to directly infer the values of the LECs [27].
In χEFT, isospin is an exact symmetry at leading order (LO) [28].At next-to-leading order (NLO), one introduces realistic charged-to-neutral pion mass splittings in the one-pion exchange potential (OPEP).The ≈ 3% mass splitting of the rather light pions induces IB in the S waves and beyond.We also have charge-dependent and nonderivative LECs in the 1 S 0 partial wave.Higherorder IB effects can be accounted for systematically by introducing, e.g., charge-dependent πN LECs, (charged) pion-photon interactions, and charge-dependent LECs in partial waves with nonzero angular momentum.However, many of those IB effects are estimated to be negligible compared to the OPEP mass splitting and nonderivative S-wave LECs [2][3][4].Indeed, the leading IB πN LECs were recently [29] inferred using N N scattering data and fifth-order χEFT, i.e., N4LO, and found to exhibit no significant charge dependence.In fact, higherorder IB effects are often neglected in quantitative chiral interactions [30][31][32][33][34].
In this paper, we employ N N interactions from χEFT up to N3LO in Weinberg power counting as defined in Ref. [3], and with the IB effects due to pion mass-splitting in the OPEP and charge-dependent leading S-wave contacts.We use nonlocal regulators in relative momenta p according to f (p) = exp −p 2n /Λ 2n with a cutoff Λ = 450 MeV and n = 3. Two-pion exchanges are spectral-function regulated [35,36] with a cutoff of 700 MeV.At N3LO, to complete the subleading two-pion exchange, we include all two-loop diagrams with some of them evaluated using numerical integration.To remedy the on-shell redundancy [32,37] in the S-wave contact potential at N3LO we set the contact fourth-order contact LECs D1 S0 , D3 S1 , and D3 S1− 3 D1 (in the notation of Ref. [3]) to zero.For describing pp and np low-energy scattering data we append the standard electromagnetic interactions up to second order in the fine-structure constant as outlined in, e.g., Ref. [38], at all chiral orders.
Within this χEFT framework, we perform a Bayesian study of the N N interaction, up to N3LO, conditioned on N N scattering data and subsequently extend the LEC posterior using an empirical value for the 1 S 0 nn scattering length a nn exp .This allows us to quantify the uncertainties of the IB effects due to the short-range LECs in χEFT.In the process of doing so, we also give an example of the flexibility of the Bayesian framework to straightforwardly expand existing results by introducing new parameters and conditioning on new data.We also test the robustness of a commonly employed model [9,37] for estimating truncation errors in χEFT.
This paper is organized as follows.In Sec.II, we out-line the statistical model upon which we base all inferences and demonstrate how to draw samples from the posterior PDFs and analyze the consistency of our inferences.In Sec.III, we draw samples from the posterior predictive distributions (PPDs) for scattering lengths and effective ranges in the 1 S 0 partial wave.We summarize our findings in Sec.IV.

II. STATISTICAL METHOD
In this section we explain our method for inferring the joint posterior PDFs for all the LECs in the N N sector of χEFT up to N3LO.After the specification of the prior and likelihood, our inference is performed in two stages.In a first step, we apply the HMC algorithm from Ref. [8] to quantify the posterior PDFs for the N N LECs in the np and pp sectors conditioned on np and pp scattering data.The posteriors presented in this paper account for an uncorrelated χEFT error model conditioned on the order-by-order convergence pattern up to N3LO.In a second step, we marginalize-in the nn nonderivative LEC C nn 1S0 and condition the inference on a single datum: the nn scattering length in the 1 S 0 partial wave, a nn exp = −18.9± 0.4 fm [26,39].We note that this is one of the currently accepted values for this scattering length and that there are conflicting experimental values for which the experimental uncertainties are not fully understood.However, we do not account for this additional level of uncertainty.See, e.g., Ref. [40] for a summary of the present status on this topic and the proposal of a novel method to measure the nn scattering length.In this paper we only operate with empirical scattering lengths and effective ranges for which electromagnetic (EM) effects have been removed by the originators of those values.Similarly, our theoretical predictions for effective range parameters do not include any EM effects either.We do however include EM effects when we compute scattering cross sections during the HMC sampling.

A. Finding an expression for the joint posterior
Using the notation pr(A|B) for the conditional probability that proposition A is true given B, the posterior PDF of interest pr(⃗ α|D, I) can, according to Bayes' theorem, be evaluated as pr(⃗ α|D, I) ∝ pr(D|⃗ α, I) × pr(⃗ α|I), i.e., as a product of the data likelihood pr(D|⃗ α, I) times the prior PDF pr(⃗ α|I).Here, ⃗ α denotes the vector of LECs to be inferred, D is employed data, and I encompasses all other assumptions and given information.The given information includes, for example, the chiral order, the regulator cutoff, masses, etc., and we make assumptions about the truncation error, data selection, and so on.Throughout this paper, we omit the overall normal-ization factor pr(D|I) as it does not play a central role in parameter estimation.
In a first step we use the np and pp scattering data D np,pp contained in the Granada database [24,25].In doing so we obtained posteriors for the πN and contact N N LECs, excluding the C nn 1S0 LEC.To simplify the notation we often denote this latter, explicitly charge-dependent, nn LEC as α nn while all other LECs are collectively de where we have assumed conditional independence in the final equality.Thus far, no strong assumptions have been made regarding the relation between α nn and ⃗ α np,pp , and the analysis is quite general.
We use HMC to sample pr(⃗ α np,pp |D np,pp , I), and the strategy that we employ for this is explained in Ref. [8].In brief, we employ order-by-order differences up to N3LO, omitting the LO results, to estimate the variance of the (uncorrelated) χEFT truncation error for describing scattering data.We place a multivariate Gaussian prior on the subleading πN LECs at NNLO and N3LO using the results from a Roy-Steiner analysis of πN scattering amplitudes [6].Furthermore, we place a rather weak prior on the N N contact LECs at all orders using an uncorrelated Gaussian PDF with zero mean and standard deviation of 5 × 10 4 GeV −(k+2) for the LECs belonging to the k = 0, 2, 4 (LO, NLO, N3LO) contact Lagrangian as defined in Weinberg power counting.The full prior factorizes into independent N N and πN PDFs since we assume no correlation between these sectors.In Sec.II B 1 we present further details about the sampling and how we revised the sampling protocol to reach N3LO.In the remainder of this section we assume that pr(⃗ α np,pp |D np,pp , I) is known to us.
Next, we extend this posterior by incorporating α nn .Using Bayes' theorem, we rewrite the first factor in the final row of (2) according to pr(α nn |⃗ α np,pp , a nn exp , I) ∝ pr(a nn exp |α nn , ⃗ α np,pp , I) where we again used that the prior for α nn is conditionally independent of ⃗ α np,pp One could certainly argue, using previous knowledge of IB, that we are in the right to place a narrow prior for α nn based on the marginal PDFs for C np 1S0 and C pp 1S0 .
However, to avoid building in such expectations on IB we selected a rather weak, and normally distributed, prior according to of width ᾱnn = 5 ×10 4 GeV −2 .We develop the likelihood in ( 4) by relating a theoretically computed value a nn th of the nn scattering length to the experimental result via a stochastic model This relation introduces the experimental error δa nn exp and the χEFT error δa nn th , which we assume is dominated by the truncation of the χEFT series.Equation ( 6) implies that other errors-such as numerical errors-are negligible.The assumption that δa nn exp and δa nn th are independent and normally distributed random variables leads to a Gaussian likelihood for the scattering length, where σ 2 exp and σ 2 th denote the variances of the experimental and theoretical errors, respectively.We use σ exp = 0.4 fm [26,40,41] as the experimental error.
To model the truncation error, we use the procedure from Refs.[9,37,42].We therefore assume that we can express the order-by-order predictions for a nn th as the sum where k is the chiral order, i.e., k = 0 is LO and k = 2, 3, 4, . . .corresponds to NLO, NNLO, N3LO, a ref is a dimensionful reference value for the scattering length, c i are dimensionless χEFT expansion coefficients, and the χEFT expansion parameter is assumed to be given by Q = m π /Λ b , which is reasonable for a quantity defined in the zero-momentum limit and analyzed in a pionful theory.To simplify notation, we omit explicit reference to the k dependence of a nn th .We employ a breakdown scale Λ b = 600 MeV in accordance with the first analysis we performed in [8].Assuming that the expansion coefficients c i (between and within chiral orders) are independent and normally distributed yields the following form for the truncation error: where c2 is the variance that characterizes the magnitude of the χEFT expansion coefficients for the scattering length and effective range.One can show [9,37] that For the analysis of the scattering length, we place a (conjugate) inverse-gamma (IG) prior c2 , leading to a marginal posterior for c2 which also follows an IG distribution but with updated values for the parameters [43].Marginalization of the variance leads to a Student's t distribution [44] for the truncation error (see Appendix A).Thus, we have that the IG prior for c2 , with initial parameters α 0 and β 0 , is updated by a set of observed expansion coefficients ⃗ c according to with where n c is the length of the vector ⃗ c.In practice, we are of course rather limited in the amount of data that we have available to infer α nn and learn about the corresponding χEFT error for a nn th .To that end we here exploit the order-by-order convergence pattern of a pp th and a np th to learn about the χEFT truncation error in a nn th .As will be further discussed in Sec.II B, we consider the NNLO-N3LO shift of a pp th to be an outlier and therefore omit this from the order-by-order data used to learn about the magnitude of the EFT truncation error.In detail, we have n c = 5 expansion coefficients for informing the marginal posterior for c2 ; these are calculated (see, e.g., Ref. [8]) from the LO-NLO, NLO-NNLOshifts of a np th and a pp th , and the NNLO-N3LO shift of a np th , all evaluated at the maximum a posteriori (MAP) locations for the PDF pr(⃗ α np,pp |D np,pp , I); see Sec.II B for details.However, with our priors, we find that the inference does not change dramatically if we include the NNLO-N3LO shift of a pp th .This is described further in Sec.III and Appendix D. We choose (α 0 , β 0 ) = (1.0,12.0), which yields the prior and resulting marginal posterior shown in Fig. 1.Our prior is predicated on our previous experience of these expansion coefficients, albeit for N N scattering data: we see it as unlikely that c2 < 4, but otherwise acknowledge our ignorance of the size of the truncation error.The mode of the prior, given by β 0 /(α 0 + 1), is located at c2 = 6.Exposure to data ⃗ c shifts the bulk of the PDF towards larger truncation errors, with the mode of the marginal posterior falling at c2 = 11.5.
To account for the marginal PDF on the χEFT truncation error, the posterior in (4) is expanded, viz., where we now explicitly indicate that the posterior is conditional on ⃗ c.In slightly more detail, the joint posterior of ⃗ α and c2 , conditioned on scattering data, scattering lengths, and order-by-order information, is the object of interest that we end up evaluating numerically.Concluding this section we list possible extensions to our analysis that we leave for future work: (1) incorporating a finite correlation length between the χEFT expansion coefficients as function of energy, (2) allowing for LEC variability, c i = c i (α nn ), (3) modifying the χEFT expansion parameter Q to a slightly greater value m π /500 MeV [7,8] or, better yet, (4) account for the uncertainty in Q by using an accompanying prior that ensures smooth matching to external (soft) momenta > m π in χEFT.

B. Evaluating posteriors
In this section we expound on our sampling of the posterior for ⃗ α np,pp using HMC and how we combine this posterior with α nn to produce a joint posterior (15) for all LECs at a given order.For clarity, we now make explicit that the inference of ⃗ α np,pp is conditional on a fixed variance c2 np,pp of the χEFT expansion coefficients.A detailed account of how to efficiently sample pr ⃗ α np,pp |D np,pp , c2 np,pp , I using HMC was given in Ref. [8].Here, we will mainly remark on new developments and results.Our procedure for extracting the joint posterior can be considered a two-step process: We follow the steps laid out in Ref. [8], with two modifications.First, we consistently use all available chiral orders, i.e., up to N3LO, to learn about the variance of the truncation error.Second, we employ a new method for rapidly tuning the parameters of the HMC algorithm such that we achieve sufficiently high sampling efficiency to perform reliable sampling at N3LO.This yields an N3LO posterior that passes all imposed MCMC convergence tests.
Let us first focus on how we learn about the truncation error.In our previous work [8] we limited ourselves to only use information about the χEFT convergence pattern up to the order at which we were sampling to estimate the variance c2 of the truncation error, and we included the zeroth-order coefficients ⃗ c 0 in the estimation.We then concluded that this procedure typically leads to an underestimation of the size of the truncation error.Here, we follow our own advice and infer the variance c2 from all orders available to us, and exclude the rather uninformative and biasing zeroth-order coefficients.We use the same grid of observables, laboratory energies, and scattering angles, as in our previous work and arrive at cnp,pp = 4.1.This value is somewhat greater than the nominally expected natural scale of c ≈ 1, yet nothing too alarming and certainly in line with what we have observed before.As discussed in Sec.II A, this is also the basis for us shifting the prior for the χEFT truncation error for the effective range towards slightly greater values.
The importance of tuning the HMC parameters, along with various strategies for achieving efficient sampling, is detailed in Ref. [8].In particular, we stress the importance of a suitable so-called mass matrix: a matrix that accounts for differences in scale between parameters in the sampled PDF.Previously we extracted a performant mass matrix by executing a short preliminary sampling and inverting the covariance matrix of the samples.Here, we instead optimize the LECs with respect to the posterior, i.e., we solve for and estimate the parameter covariance matrix at that optimum.We use the BFGS optimization algorithm [45][46][47][48] which is a so-called quasi-Newton method that relies on first-order gradients to update approximations to the Hessian employed in the Newton algorithm.We employ automatic differentiation (AD) [38,49,50] to obtain the first-order gradients.In the end, we find that the approximate Hessian at the optimum, as found by the BFGS algorithm, makes for a suitable mass matrix in our application.In our case we have access to higher-order gradients via AD and could, in principle, obtain ⃗ α ⋆ using any expedient optimization method and compute an exact-tomachine-precision Hessian.However, we find that BFGS is sufficient for tuning the HMC sampler.Bypassing the need for the ≈ 20 000 function evaluations required during a preliminary HMC sampling is certainly a great improvement, as the BFGS algorithm typically terminates after at most a few hundred function evaluations.Better yet, the optimization-based method is in our experience more reliable and ultimately yields slightly higher performance in the subsequent sampling of the posterior.
We gather N ≈ 10 000 samples in each of the ten chains we run at each chiral order.We find integrated autocorrelation times τ (τ NLO = 0.71, τ NNLO = 0.93, and τ N3LO = 3.33) that reach a plateau within very few samples, thus indicating converged HMC chains.The τ value should not be too great since it is inversely proportional to the effective sample size that measures the number of independent samples in the HMC chain.We note that the NLO and NNLO autocorrelation times are similar to our previous results [8], achieved with the more cumbersome method of preliminary samplings.We should reiterate the concerns that always surround the topic of convergence in MCMC.In real-world applications it is not possible for us to declare a chain "converged".This is because we can never explore the entire parameter space in finite time, and we can miss non-negligible (or even dominating) probability regions.All we can do is to probe the MCMC chains for signs of non-convergence.This issue becomes progressively more pressing as the dimensionality and domain of the parameter space grows, and is especially concerning at N3LO in our case.With that said, we have searched the parameter domains to the best of our abilities by initializing the BFGS optimization, and the HMC algorithm, at multiple locations and we did not find any signs of multimodality and non-convergence.
In Table I we predict the scattering lengths and effective ranges for the np and pp isospin channels in the 1 S 0 partial wave at the MAP points ( 16) of the LEC posteriors at different orders.Overall, the predictions are reasonably close to the experimental, or rather the empirical, values [3].This builds confidence in our model inference.These point estimates give a first indication of our model's performance in the low-energy limit.Overall, we observe convergence towards empirical results as we move to higher chiral orders.A notable exception is the N3LO prediction of the scattering length in the pp isospin channel.It differs from the empirical result by 1.9 standard deviations, compared to just 1.2 standard deviations for the NNLO prediction.Furthermore, the difference is in the opposite direction.This is possibly caused by a lack of informative data during the inference of the LECs.

Step 2: Evaluating the joint posterior
Equipped with an MCMC chain to represent pr ⃗ α np,pp |D np,pp , c2 np,pp , I , we proceed to evaluate the posterior pr ⃗ α, c2 |D, ⃗ c, I in Eq. (15).Due to the expected weakness of IB effects in our χEFT we can be fairly certain that the posterior probability mass of the latter PDF will not be very far from the one of the former.We take advantage of this expectation and proceed by using the principles of sampling/importance sampling [51] where we first define a sampling distribution with π(α nn ) a simple, bounded uniform distribution.Since the likelihood for the single nn scattering length poses no significant computational challenge, for each chiral order we first explore the α nn space to identify an interval of values encompassing the bulk of the probability mass.These intervals are (in units of 10 i ) from the sampling distribution (17).In practice, we go through the entire MCMC chain of ⃗ α np,pp sequentially since we know [8] that these are independent and random samples from pr ⃗ α np,pp |D np,pp , c2 np,pp , I .For each ⃗ α np,pp sample we draw a sample of c2 from its marginal posterior pr(c 2 |⃗ c, I), see Eq. ( 12), and a sample of α nn from π(α nn ).

Evaluate the ratio ω
Due to the form of our sampling distribution (17) and the factorization of the joint posterior (15), this ratio is simply the product of the Gaussian likelihood for the scattering length (7) and the prior for α nn (5).
Once we have the lists of samples (⃗ α i , c2 i ) we compute the normalized (importance) weights q i = ω i / j ω j .This provides us with a weighted chain of (⃗ α, c2 ) values distributed as the target posterior pr ⃗ α, c2 |D, ⃗ c, I (15).The method described here is simple to apply but may not work well in all cases, e.g., if the marginal posterior of α nn is relatively unknown and/or defined for a highdimensional parameter domain.We have validated our results obtained using importance sampling by approximating pr ⃗ α np,pp |D np,pp , c2 np,pp , I with a multivariate normal distribution and directly sampling the full distribution (15) using HMC.Multivariate normal approximations to all posteriors are provided in the Supplemental Material [52].
We report bivariate and univariate marginals of the joint LEC posteriors in Appendix B. In Table II we make point-estimate predictions for selected few-nucleon ground-state observables using the MAP point of the joint posteriors.This is a first check of the model inference and the results are reasonable.We note that the N3LO predictions deviate more from experiment than NNLO.However, we do not incorporate N N N interactions in this analysis and for that reason we also refrain from estimating χEFT truncation errors for the results presented in Table II.On the other hand, we perform a full analysis of the N N effective range expansion in Sec.III.
Apart from the inclusion of C nn 1S0 and c, the NLO and NNLO posteriors are overall similar to their counterparts in Ref. [8].The differences that do occur arise prior to the inclusion of C nn 1S0 as a result of different characterizations of the χEFT truncation error.The correlation patterns are nearly identical, and the majority of the marginal posteriors overlap at the 68% level.However, some parameters in higher partial waves have somewhat shifted values.For example, in Ref. [8] we reported C 3P 2 = −0.162(1)× 10 4 GeV −4 at NLO, while Fig. 5 reveals that C 3P 2 = −0.175(2)× 10 4 GeV −4 .
The posteriors presented here are confined to rather small volumes in parameter space and, as expected, the marginal posteriors for the C 1S0 LEC acting in the different isospin channels are similar to each other.The TABLE II.MAP predictions of few-nucleon ground-state energies E (in MeV) and point-proton radii R (in fm) for nuclei with mass number A = 2, 3, 4 obtained using the Jacobi no-core shell-model [53] in a harmonic oscillator basis with frequency 22 MeV/ℏ in model spaces with 251, 41, and 21 oscillator shells, respectively.Note that these predictions are not including a three-nucleon (N N N ) interaction.For 2 H we also list the D-state probability (PD) in % and the electric quadrupole moment (Q) in units of eb.The experimental and empirical values are the same as in Ref. [38].shows stronger correlations with other LECs than at NLO and NNLO.The fourthorder contact LECs, D, are sensitive to high momentum observables and therefore less constrained by the data.We find that several of the D LECs return the prior if we do not condition the posterior on data in the T lab ∈ [40.5, 290] MeV region.Furthermore, the posteriors indicate that some D LECs are notably unnatural; for instance, we have D 1S0 = −29.1(3)× 10 4 GeV −6 , which is of particular interest in this case as it acts in the same partial wave as-and is strongly correlated with-our explicitly isospin breaking LECs.The inference of this LEC, then, is influential on our predictions of 1 S 0 scattering lengths and may explain the stark difference between the N3LO predictions for a np and a pp in Table I.There are precedents for large values of D 1S0 ; see, e.g., Ref. [54].Unlike earlier works, we also find rather unnaturally sized MAP values for D 3S1 = −32.2(3)× 10 4 GeV −6 and D 3P 1 = −20.4(1)× 10 4 GeV −6 , where the uncertainties denote 68% credible intervals.At N3LO, the various C and C parameters are akin to their values inferred at NNLO, with the exception of C 3P 1 which is both several times larger and has the opposite sign, i.e., C 3P 1 = −0.956(5)× 10 4 GeV −4 at NNLO and C 3P 1 = 6.13(1) × 10 4 GeV −4 at N3LO.
We have also investigated how the posteriors are affected if we do not exclude the unexpectedly large NNLO-N3LO shift in a pp th from the estimation of the truncation error (see Sec. II A).We find that the width of the marginal C nn 1S0 posterior roughly doubles at NLO, widens by a small but noticeable amount at NNLO, and is virtually unaffected at N3LO.From this we draw the conclusion that the theoretical uncertainty dominates at NLO, the theoretical and experimental uncertainties are roughly equal at NNLO, and the experimental uncertainty dominates at N3LO.We find nonoverlapping priors and posteriors for the πN LECs at both NNLO and N3LO.This tension was previously seen at NNLO in Ref. [8].It turns out that if the posteriors are conditioned on only a low-energy data set, T lab ∈ [0, 40.5]MeV, we return our priors from a Roy-Steiner analysis.This is summarized for the NNLO case in Fig. 2.Although there are significant deviations between the πN prior and posterior when we condition on high-energy data, the discrepancies are rather small on a naturalness scale.As such, the statistical model we have set up to relate low-energy data and χEFT at NNLO appears to preserve the long-range physics rather well.Moreover, the credible intervals of the πN posterior cannot be straightforwardly compared with the variance of the prior determined in a Roy-Steiner analysis [5], for which it is difficult to estimate the truncation uncer-tainty.At N3LO we similarly find that conditioning on the energy-truncated data set returns the prior for the πN LECs.The prior and posterior at N3LO (conditioned on the full data set) are shown side by side in Appendix B. It is peculiar that the marginal and univariate posterior for the πN LEC c 4 sits on top of the prior.

III. POSTERIOR PREDICTIVE DISTRIBUTIONS OF 1 S0 ISOSPIN BREAKING
Equipped with samples from the joint posterior PDFs of the N N LECs up to N3LO we proceed to sampling the posterior predictive distributions (PPDs) for the S-wave scattering length a and effective range r defined by the effective range expansion (ERE) Here, p is the relative momentum and δ(p) is the scattering phase shift.The PPD is the PDF of unobserved values for a and r conditioned on np and pp scattering data as well as the empirical nn scattering length.We quantify the PPD in all three isospin channels t z =pp,np,nn of the 1 S 0 partial wave.The full PPDs are given by the finite set of samples and with (⃗ α, c2 ) ∼ pr ⃗ α, c2 |D, ⃗ c, I .Furthermore, the conjugacy of the IG prior for the χEFT truncation error variance enables a pointwise and closed-form evaluation of the associated χEFT errors δa th and δr th ; see Ref. [43] and Appendix A. As a result of our sampling procedure (see Sec. II B) the sets in Eqs.(19) and (20) are composed of approximately 10 5 weighted samples.In the following we present results based on an inference of c where we excluded the NNLO-N3LO shift of the a pp th (see Table I).It turns out that this does not impact any predictions beyond NLO; see Appendix D for comparison.
The results at NLO, NNLO, and N3LO are shown in Fig. 3. Clearly, our model is not fine-tuned to reproduce any of the ERE parameters in Fig. 3 except a nn , which therefore serves as another model check, but all PPDs agree with the empirical results within their uncertainties and the different distributions at different orders overlap reasonably with each other.The exception here is the N3LO prediction of a pp , that we also omitted from the estimation of the χEFT truncation error.In contrast, for a np , both NNLO and N3LO agree perfectly with the rather precise empirical value.For the empirical values of the pp ERE parameters, the error bands emerge predominantly from the model dependence in the analysis of scattering data and removal of electromagnetic effects.The PPDs for the effective ranges show a great deal of congruity.This implies that the N N contact C1 S0 that enters at NLO with a quadratic momentum dependence is sensibly inferred.This contact is missing at LO and the predictions of effective ranges at that order are thus rather poor, as seen in Table I.Overall, one must go beyond NNLO to achieve predictions that are more precise than the corresponding empirical uncertainties.
To quantify the strength of IB effects, we convert the PPDs for the ERE parameters in the 1 S 0 partial wave to standard measures of the CIB and charge symmetry breaking (CSB), where the latter amounts to a π rotation around the y axis in isospace, i.e., The PPDs for the CIB and CSB effects, including χEFT errors, are summarized in Fig. 4. We first note that it is only at NNLO and beyond that we can detect, with confidence, the overall magnitudes of IB effects in the 1 S 0 partial wave.Only CIB in the scattering length can be said to be different from zero with any confidence at NLO.We also find that our PPDs for CSB and CIB agree, within uncertainties, with existing empirical data and a range of point estimates using well-known chiral potentials at NNLO and N3LO.These point estimates are all rather close to each other and fall within the empirical uncertainties at all orders.The CSB and CIB in the effective range is somewhat underestimated and overestimated, respectively, compared to the empirical values.The outlier N3LO result for a pp of course propagates to the results for CIB and CSB and induces a comparatively large isospin breaking at this order.Yet, the results are in line with other chiral potentials.

IV. SUMMARY AND OUTLOOK
In this paper, we sampled high-dimensional posteriors up to N3LO in χEFT and studied the effects of IB in the N N sector.We used Bayesian inference, conditioned on np and pp scattering data as well as the empirical value a nn exp for the 1 S 0 scattering length in the nn channel, to infer posterior PDFs for the LECs at LO, NLO, NNLO, and N3LO.We split the inference in two steps.First we employed HMC and conditioned the LEC posteriors on np and pp scattering data and accounted for uncorrelated χEFT truncation errors.A new approach to tuning the mass matrix, based on posterior optimization, enabled us to extract a converged 31-dimensional posterior also at N3LO.We find such advancements pivotal for enabling robust MCMC sampling of the N3LO posterior.In the next step of the inference, we included a nn exp and employed importance sampling to marginalize  4. CSB and CIB in the 1 S0 scattering lengths and effective ranges; see (21).Colors and shadings as in Fig. 3. Empirical results (ann from Ref. [26]; all others from Ref. [3]) are shown in black with gray error bands.Calculations using three widely used potentials are also shown: NNLOsat [31], Idaho-N3LO [54], and ∆-NNLOGO(394) [55].Note that the three latter potentials agree almost exactly regarding ∆rCSB,CIB and are thus difficult to visually distinguish.Also note that the histograms have been scaled such that they have the same peak height.
in the nn contact LEC.For the variance of the truncation error in the ERE parameters in the 1 S 0 partial wave we employed a conjugate IG prior.In the end, we find that the resulting LEC posteriors at NLO and NNLO match our previous [8] posteriors in the np and pp isospin sectors and that the LEC C nn 1S0 is roughly twice as broad but exhibits a correlation pattern similar to C np 1S0 and C pp 1S0 .Next, we sampled the PPDs for the ERE parameters and found that our results are consistent with existing point estimates using well-known chiral potentials, with the exception of a pp at N3LO.This outlier might be traced back to insufficient information in the N N scattering data set we conditioned the first part of the inference on.When accounting for χEFT truncation errors, we find that one must go to NNLO to confidently detect IB effects.
One way to improve a Bayesian analysis of the strong interaction in the low-energy region of relevance to the ERE might be to mix pionless EFT [4,56] and χEFT.Indeed, χEFT harbors an intrinsic uncertainty due to the low-energy scale set by m π .Since the exact form of the χEFT expansion parameter Q for scattering amplitudes is obscured by the non-perturbative resummation in the Lippmann-Schwinger equation its uncertainty pr(Q|I) should be accounted for.It could then be interesting to apply Bayesian mixture models to combine pionless EFT and χEFT predictions.
When we condition the inference on low-energy N N scattering data with T lab ∈ [0, 40.5]MeV we return the Roy-Steiner priors for the πN LECs at NNLO and N3LO.This is in accordance with χEFT being a low-energy theory with long-ranged physics governed by the πN interaction.When conditioning the inference on all N N data up to the pion-production threshold at 290 MeV, the marginal πN LEC posteriors are significantly shifted with respect to the priors.Still, the discrepancies are small on a naturalness scale.
Ab initio studies [55,57] indicate that the inclusion of the ∆ isobar yields more realistic predictions for bulk properties of nuclei and nuclear matter.However, a Roy-Steiner prior for the πN LECs with explicit ∆ isobars is less precise [6].As such, it will also become more important to employ additional nuclear structure and reaction data in the inference and employ error models that account for correlations between the χEFT expansion coefficients.Although it rapidly becomes computationally challenging to evaluate likelihoods encompassing nuclear data for increasing mass numbers, the development of fast and accurate emulators [58][59][60][61] appears to provide sufficient leverage.In particular, the present work serves as an example how to sequentially incorporate multiple sets of low-energy nuclear structure and reaction data in a Bayesian framework to perform inference and quantify the theoretical precision.uncertainties as a dark (light) gray area.The empirical results are without electromagnetic effects and gathered from Ref. [3], except for ann which we take from Ref. [26].Note that the histograms have been scaled such that they have the same peak height.
FIG. 1.The prior pr(c 2 |I) (dashed line) for c2 and the resulting marginal posterior pr(c 2 |⃗ c, I) (solid line).The latter is obtained from allowing the order-by-order predictions of a np th and a pp th to flow through the prior via the observed expansion coefficients ⃗ c.The individual (squared) elements of ⃗ c are shown as vertical lines, with black (orange) representing np (pp) elements.The center dot in the y-axis label pr(c 2 |•) acts as a placeholder for "I" and "⃗ c, I".

FIG. 2 .
FIG. 2. Prior and posterior PDFs for the πN LECs c1, c3, c4, in units of GeV −1 .The posteriors were obtained using HMC.The inner (outer) ellipses enclose approximately 39% (86%) of the probability mass.Black lines indicate the prior as provided by the Roy-Steiner analysis [6].Blue lines indicate the posterior conditioned on low-energy scattering data only, i.e., with T lab up to 40.5 MeV.Red lines indicate the posterior conditioned on scattering data up to T lab = 290 MeV.

TABLE I .
[3] predictions of np and pp scattering lengths (a) and effective ranges (r) in units of fm.Empirical results are from Ref.[3].
[26]3.PPDs of scattering lengths and effective ranges at NLO (blue), NNLO (purple), and N3LO (red).Empirical results are shown as black lines, with corresponding 1σ (2σ) uncertainties as a dark (light) gray area.The empirical results are without electromagnetic effects and gathered from Ref.[3], except for ann which we take from Ref.[26].Note that the histograms have been scaled such that they have the same peak height. FIG.