Towards high-order calculations of three-nucleon scattering in chiral effective field theory

We discuss the current status of chiral effective field theory in the three-nucleon sector and present selected results for nucleon–deuteron scattering observables based on semilocal momentum-space-regularized chiral two-nucleon potentials together with consistently regularized three-nucleon forces up to third chiral order. Using a Bayesian model for estimating truncation errors, the obtained results are found to provide a good description of the experimental data. We confirm our earlier findings that a high-precision description of nucleon–deuteron scattering data below pion production threshold will require the theory to be pushed to fifth chiral order. This conclusion is substantiated by an exploratory study of selected short-range contributions to the three-nucleon force at that order, which, as expected, are found to have significant effects on polarization observables at intermediate and high energies. We also outline the challenges that will need to be addressed in order to push the chiral expansion of three-nucleon scattering observables to higher orders.


Introduction
The past few years have seen remarkable advances towards pushing the precision frontier of chiral effective field theory (EFT) in the two-nucleon sector; see Refs. [1][2][3] for review articles. The fifth-order (N 4 LO) contributions to the nucleon-nucleon (NN) force have been worked out in Ref.
a e-mail: evgeny.epelbaum@rub.de (corresponding author) [4] and a new generation of accurate and precise chiral EFT NN potentials up through N 4 LO has been developed in Refs. [5][6][7][8], see also Refs. [9][10][11] for related developments. The N 4 LO interactions of Refs. [7,8] utilize the values of the pion-nucleon (π N) low-energy constants (LECs) determined from matching chiral perturbation theory to the Roy-Steiner-equation analysis of π N scattering at the subthreshold point [12], but differ substantially by the regularization procedure. The potentials developed by our group in Ref. [7] build upon our earlier studies [5,6] and employ a local regulator for pion-exchange contributions which, per construction, maintains the analytic structure of the longrange interaction. Differently to the nonlocally regularized potentials of Ref. [8], the interactions constructed in Ref. [7] do not produce long-range artifacts at any finite-order of expansion in inverse powers of the momentum-space cutoff Λ. The resulting semilocal 1 momentum-space regularized (SMS) potentials of Ref. [7] are currently the most precise chiral EFT interactions and provide, at the highest order N 4 LO +2 , a nearly perfect and Λ-independent (within the considered cutoff range) description of neutron-proton and proton-proton scattering data below pion production thresh-old from the self-consistent 2013 Granada data base [13]. The achieved description of the scattering data is comparable to or even better than that based on the so-called high-precision semi-phenomenological potentials like the AV18 [14], CD Bonn [15] and Nijm I, II [16] models.
In spite of this exciting progress in the NN sector, applications of chiral EFT to three-and more-nucleon systems and in processes involving external sources are, with very few exceptions, still limited to the next-to-next-to-leading order (N 2 LO) in the chiral expansion. What makes high-accuracy calculations of such systems/reactions beyond N 2 LO so difficult? The main technical and conceptual difficulties are related to the treatment of many-body forces and exchange current operators. Three-(3NF) and four-nucleon forces (4NF) start to contribute at third (N 2 LO) and fourth (N 3 LO) orders of the chiral expansion, respectively, while the first contributions to the exchange electroweak currents appear already at second order (NLO) relative to the dominant single-nucleon terms. Over the past decade, we have worked out off-shell-consistent expressions for the 3NF [17,18], 4NF [19,20] and electroweak charge and current operators [21][22][23][24] completely up through N 3 LO using dimensional regularization (DR) to compute pion loop contributions; see also Refs. [25][26][27][28] for a similar work on electroweak currents by the JLab-Pisa group and Ref. [29] for a related discussion. Furthermore, selected N 4 LO contributions to the 3NF have been worked out in Refs. [30][31][32], and the longest-range part of the 3NF was also analyzed in the framework involving Δ(1232) degrees of freedom. On the technical side, the implementation of the 3NFs and exchange currents in fewbody calculations requires their partial-wave decomposition. This nontrivial task can nowadays be accomplished numerically for a general 3NF specified in momentum space as described in Refs. [33,34]. The main conceptual challenge that still needs to be addressed concerns a consistent regularization of many-body forces and exchange operators. While this issue is irrelevant in the NN sector, a naive regularization of many-body interactions and exchange currents based on the expressions derived using DR violates chiral symmetry and leads to inconsistent/wrong results at the one-loop level (i.e. at N 3 LO) and beyond [35,36]. A possible solution could be provided by using higher-derivative regularization instead of DR, which has to be chosen in a way compatible with the regularization scheme of Ref. [37] in the NN sector. Work along these lines is in progress.
In this paper we update our recent study [38] based on the semilocal coordinate-space regularized (SCS) NN forces of Refs. [5,6] and analyze nucleon-deuteron (Nd) scattering observables using the new SMS chiral NN potentials of Ref. [7] in combination with the SMS 3NF at N 2 LO. We also refine our previous estimations of truncation uncertainties by employing a Bayesian approach instead of the algorithm proposed in Ref. [5] and used in our earlier studies [38][39][40] in the three-nucleon sector. Last but not least, we explore the role of selected N 4 LO short-range 3NF terms in Nd scattering.
Our paper is organized as follows. In Sect. 2, we specify our Bayesian model for truncation errors by following the approach proposed in Ref. [41]. Our results for Nd scattering observables at N 2 LO are presented in Sect. 3, while the role of selected N 4 LO 3NF operators is discussed in Sect. 4. The main results of our study are summarized in Sect. 5.

Bayesian model for truncation errors
A reliable estimation of theoretical uncertainties is an essential ingredient of any systematic approach such as chiral EFT. Cutoff variation offers one possibility to quantify the impact of contributions beyond the truncation level. However, in the few-and many-nucleon sectors, the available cutoff range is often rather limited. Furthermore, cutoff variation does not allow one to probe the impact of neglected long-range interactions. In Ref. [5], a more reliable, universally applicable algorithm for estimating truncation errors using the available information on the chiral expansion for any observable of interest without relying on cutoff variation was proposed. Here and in what follows, this algorithm will be referred to as the EKM approach. For applications of the EKM method to a broad range of low-energy reactions in the single-baryon and few-/many-nucleon sectors, see Refs. [42][43][44][45] and [6,[38][39][40][46][47][48], respectively. Being very simple and easy to implement, the EKM approach does, however, not directly provide a statistical interpretation of the estimated uncertainties. In Refs. [41,49], a general Bayesian approach to calculating the posterior probability distribution for truncation errors in chiral EFT was developed. The EKM approach was then shown to essentially correspond to a particular choice of prior probability distribution for the coefficients in the chiral expansion of observables. Using the chiral NN potentials from Ref. [5], the EKM error estimations for the neutron-proton (np) total cross section at selected energies were found in the Bayesian approach of Ref. [49] to be consistent with 68% degree-of-belief (DoB) intervals.
Throughout this paper, we employ a slightly modified version of the Bayesian model from Ref. [49]. Specifically, consider a two-nucleon scattering observable X ( p) with p referring to the center-of-mass (CM) momentum. Calculating X ( p) using chiral EFT potentials at various orders Q i , i = 0, 2, 3, 4, . . ., (but for a fixed cutoff value) yields the corresponding predictions X (i) ( p), and the chiral expansion of X ( p) can be written in the form X = X (0) + ΔX (2) + ΔX (3) + ΔX (4) with ΔX (2) := X (2) − X (0) and ΔX (i) := X (i) − X (i−1) for i > 2. The second equality serves as a definition of the dimensionless expansion coefficients c i . The reference value X ref , which sets the overall scale, will be defined below.
Here and in what follows, the expansion parameter of chiral EFT is assumed to have the form [5,6,41,49] where Λ b is the breakdown scale of the chiral expansion. The quantity M eff π serves to model the expansion of NN observables around the chiral limit. In Refs. [5,6,41,49], this scale was set to the physical pion mass M π . However, as pointed out in Ref. [50], the error plots in Ref. [5] indicate that the transition between the two expansion regimes in the NN sector actually appears at a scale M eff π higher than M π . On the other hand, both Bayesian model parameters Λ b and M eff π can be determined/tuned empirically by calculating the success rates for a given set of observables and/or energies. In particular, Ref. [49] confirmed the EKM estimation Λ b ∼ 600 MeV based on the results for the total np cross section and using the SCS potentials of Ref. [5], but it also found somewhat larger values of Λ b to be statistically consistent; see also a related discussion in [51]. A similar empirical analysis was performed in Ref. [50] for both Λ b and M eff π yielding the values of Λ b ∼ 650 . . . 700 MeV and M eff π ∼ 200 MeV. Suppose the results for the observable X ( p) are available up through the order X (k) , k ≥ 2. The goal is then to estimate the truncation error δ X (k) ≡ i>k ΔX (i) resulting from neglecting the unknown higher-order contributions, i.e. to compute the posterior probability distribution function (pdf) for δ X (k) given the explicit knowledge of {X (0) , X (2) , . . . , X (k) }. The Bayesian model of Ref. [41] uses the leading-order (LO) result for X ( p) to set the overall scale in order to define the dimensionless expansion coefficients c i with c 0 = 1 in Eq. (1). As we will argue below, this approach may, for certain choices of the prior pdf, be too restrictive in the kinematical regions near the points where the LO contribution vanishes e.g. by changing the sign. Such situations are not uncommon if one looks at observables which depend on continuously varying parameters such as energy or scattering angle, see also a discussion in Ref. [51]. In such circumstances it is advantageous to set the overall scale from the next-to-leading order (NLO) contribution ΔX (2) via X ref = ΔX (2) /Q 2 in order to avoid underestimating X ref .
To have an approach applicable when both X (0) and/or ΔX (2) are accidentally small, we set the scale X ref via for k = 2 and for k ≥ 3. Let c m = 1, m ∈ {0, 2, 3}, be the coefficient used to define the overall scale. Assuming that the remaining coefficients c i are distributed according to some common pdf pr(c i |c) with a hyperparameterc and performing marginalization over h chiral orders k + 1, . . . , k + h, which are assumed to dominate the truncation error, the probability distribution for the dimensionless residual Δ k ≡ ∞ n=k+1 c n Q n k+h n=k+1 c n Q n to take a value Δ k = Δ, given the knowledge of {c i≤k }, is given by [41] pr h (Δ|{c i≤k }) = ∞ 0 dc pr h (Δ|c) pr(c) i∈A pr(c i |c) ∞ 0 dc pr(c) i∈A pr(c i |c) , (6) where the set A is defined as Here and in what follows, we employ the Gaussian prior of "set C" from Ref. [41], namely and assume a log-uniform probability distribution, [52] pr(c) = 1 ln(c > /c < ) The nice feature of this prior is that all integrations in Eq. (6) can be carried out analytically. For the sake of completeness, we give in Appendix A the corresponding expressions for pr h (Δ|{c i≤k }) from Ref. [41]. The posterior pdf is an even function of Δ, and for any given DoB interval, the corresponding value of the residual Δ k and the resulting truncation error δ X (k) = X ref Δ k can be readily obtained by numerically integrating pr h (Δ|{c i≤k }) over Δ.
As a first application, we employ the Bayesian model of Ref. [41] and set the overall scale X ref solely from the corresponding LO contribution; see Eq. (3). In the bottom row of Figs. 1 and 2, we show the 68% and 95% DoB intervals for selected np scattering observables at E lab = 143 MeV for the non-informative prior C [49] withc < = and c > = 1/ , which makes no assumption as regards either the maximum or the minimum size ofc by taking the limit → 0; see Eq. (23). Here, the breakdown scale was set to Λ b = 700 MeV, and the resulting Bayesian model is referred to as C 700 .
Being non-informative, the prior set C is generally expected to yield conservative estimates for truncation errors. Open circles refer to the results of the Nijmegen partial-wave analysis [53]. Data for the cross section are at E lab = 142.8 MeV and taken from [54]. The first, second, third and fourth rows correspond to the Bayesian models C 650 0.5−10 , C 650 0.5−10 ,C 650 0.5−10 and C 700 . All results shown are based on the SMS NN potentials using the cutoff of Λ = 450 MeV However, settingc < → 0 yields a δ-function-like posterior pr C h (Δ|{c i≤k }) for c 2 k → 0 as can be seen from Eq. (23), i.e. this model fails to provide an adequate estimation of the truncation error if the corrections Δ (2) , . . . , Δ (k) happen to be accidentally small. For the examples shown in Figs. 1 and 2, this is the case at NLO for the differential cross section dσ/dΩ around θ CM ∼ 85 • , for the polarization transfer D t around θ CM ∼ 140 • , for the spin-correlation coefficients C kp at θ CM ∼ 65 • and for the coefficient C kk at θ CM ∼ 25 • , and θ CM ∼ 115 • . In all these cases, the corresponding functions ΔX (2) (θ CM ) change their sign. To circumvent the problem with the underestimation of the truncation error in such kinematical regions, the authors of Ref. [41] suggested to use a more informative (but not too restrictive) prior set C 0.25−10 corresponding toc < = 0.25 andc > = 10. Here and in what follows, we make the choicec < = 0.5, which we found to be more efficient in resolving the above-mentioned issue while still sufficiently general. As shown in the upper row of Figs. 1 and 2, the prior set C 650 0.5−10 indeed yields reasonable estimates of the truncation errors for dσ/dΩ. However, the more informative prior withc < = 0 suffers from another issue: it yields a vanishingly small truncation error at all orders in the cases when X (0) happens to be accidentally small. This is the case for D t at θ CM ∼ 100 • and for C kk at θ CM ∼ 10 • , θ CM ∼ 75 • and θ CM ∼ 180 • ; see the plots in the upper row of Figs. 1 and 2. The most extreme situation is observed for the coefficient C kp , for which the LO contribution appears to be small for all scattering angles. The problem can be traced back to the misidentification of the overall scale by Eq. (3) in such accidental cases. Writing Here and in what follows, the resulting Bayesian model is referred to asC. However, while highly unlikely, it is still possible that both the LO contribution and the NLO correction are simultaneously accidentally small. For the considered observables, this happens for D t at backward and for C kk at forward angles. To prevent underestimating the truncation errors in such rare cases, we replace, for k ≥ 3, Eq. (4) with Eq. (5) and refer to the resulting Bayesian model asC. This model is found to yield robust results for NN scattering observables in all kinematical regions; see e.g. the third row in Figs. 1 and 2 and will be employed in nucleon-deuteron scattering calculations considered in the next sections.
It is important to emphasize that the considered examples have been selected to visualize the possible issues in certain accidental situations, and that the differences between the considered Bayesian models are fairly minor for most observables. This is exemplified in Fig. 3, where the truncation errors for the total cross section are shown for a variety of considered models including the original EKM approach (with M eff π = M π and Λ b = 600 MeV). As pointed out in Ref. [49], the dependence on a particular Bayesian model and/or assumed prior set decreases with an increasing order, i.e. with the increasing amount of information about the actual pattern of the chiral expansion. Notice that differently to the EKM approach, the considered Bayesian models exploit only the information up to the order, at which the truncation error is estimated. The more conservative error estimations at E lab = 50 MeV with the  Fig. 1 but for the spin-correlation parameters C kp and C kk models C andC as compared to the original EKM approach are mainly due to the larger value of M eff π . For the considered set of the total cross section calculations, counting the success rate for the next-higher order result to lie within the estimated uncertainty as shown in Fig. 3 yields the values of 62.5 . . . 75% (100%), which are statistically consistent with the DoB intervals of 68% (95%). 3 We have verified that this conclusion also holds true for the larger set of energies considered in Ref. [50]. We refrain from performing similar statistical tests for angular distributions due to their correlated nature. This could be done using the method proposed in Ref. [51], which is based on Gaussian processes and encodes correlation structure of coefficients c i (θ CM ). 3 For the EKM approach, the success rate equals 100% per construction.

Nucleon-deuteron scattering at N 2 LO using SMS nuclear potentials
We now turn to the main topic of our study and consider nucleon-deuteron scattering in chiral EFT. The N 2 LO threenucleon force is given by [55,56] where the subscripts refer to the nucleon labels, q i = p i −p i with p i and p i being the final and initial momenta of the nucleon i and σ i (τ i ) are the Pauli spin (isospin) matrices. Further, c i , D and E denote the corresponding low-energy constants (LECs), while g A and F π refer to the nucleon axial coupling and pion decay constant, respectively. It is custom-ary to express the LECs D and E in terms of the corresponding dimensionless coefficients via where, following [38,56], we use Λ χ = 700 MeV. For the LECs c i , we employ the central values from the Roy-Steinerequation analysis of Ref. [12] at the corresponding chiral order, namely The same values are used in the SMS NN potentials of Ref. [7] at N 2 LO. Differently to Ref. [38], we employ the same semilocal momentum-space regulator for the 3NF as in the NN potentials of Ref. [7] by replacing the pion propagators via For the D-term, the contact interaction between nucleons 1 and 2 is, in addition, regularized by multiplying the matrix elements with a nonlocal Gaussian regulator exp(−(p 2 12 + p 12 2 For the contact interaction proportional to the LEC E, we apply a nonlocal Gaussian regulator in momentum space, where It is important to emphasize that the SMS NN potentials of Ref. [7] employ additional (local) short-range subtractions to ensure that the coordinate-space expressions of the regularized pion-exchange contributions and derivatives thereof vanish at the origin. This convention ensures that regularized pion-exchange contributions contain only long-range pieces. On the other hand, using the regulator in Eq. (12), the resulting TPE contributions still contain admixtures of short-range terms of the D-and E-types.
Following Ref. [38], we determine the LECs c D and c E from the 3 [15] (CD Bonn NN potential in combination with the Tucson-Melbourne 3NF [57]). Open circles are neutron-deuteron from Ref. [58] and protondeuteron data from Ref. [59][60][61], corrected for the Coulomb effects; see Ref. [56] for details for the cutoff Λ = 500 MeV, along with the truncation errors corresponding to the DoB intervals of 68% and 95%.
Notice that the truncation errors are symmetric, and the actual results of our calculation lie in the middle of the corresponding error bands.
To estimate the truncation uncertainty we have used the Bayesian modelC 650 0.5−10 introduced in Sect. 2. For scattering reactions involving three and more nucleons, we, however, also need to specify the pertinent momentum p scale that sets the expansion parameter Q in Eq. (2). Consider nucleonnucleus scattering and let E N lab be nucleon beam energy in the laboratory system. Neglecting the neutron-proton mass difference and the binding energy of the target nucleus, which is assumed to consist of A nucleons, the CM momentum is related to E N lab via that ensures that p coincides with p CM in the NN system. Here, s is the usual Mandelstam variable. One can thus express the scale p in terms of E N lab via Predictions for the differential cross section, nucleon and deuteron analyzing powers A n y and A d y as well as deuteron tensor analyzing powers A yy , A xz and A xx in elastic nucleon-deuteron scattering at laboratory energy of E N lab = 135 MeV at NLO (yellow bands) and N 2 LO (green bands) based on the SMS NN potentials of Ref. [7] for Λ = 500 MeV. Open circles are proton-deuteron data from Ref. [62]. For remaining notation, see Fig. 4 with s = m 2 (A + 1) 2 + 2 Am E N lab . In the nonrelativistic approximation, this relation simplifies to The nonrelativistic approximation holds at a sub-percent level for the energy range considered in this study and we use the relation (17)  Notice that the residual cutoff dependence of the considered observables at NLO is similar to the one at N 2 LO and is, in most cases, comparable with the N 2 LO 68% DoB intervals. These results demonstrate that our calculations for different values of Λ are consistent with each other within errors. Notice further that the softest cutoff of Λ = 400 MeV shows the largest deviation from the bulk behavior and from the experimental data, which presumably points to the increasing amount of finite-regulator artifacts.
We also show in Figs. 4, 5, 6 and 7 the results based on the CD Bonn NN potential [15] with and without the Tuscon-Melbourne (TM) 3NF [57]. In particular, for the cases where the TM 3NF is known to provide sizable corrections such as e.g. for the differential cross section around its minimum and for the deuteron analyzing powers at the intermediate energies of E N lab = 70 and 135 MeV, our N 2 LO results agree well with the CD Bonn NN plus TM 3NF calculations, while the predictions based on the CD Bonn NN force alone are often outside the 65% and sometimes even 95% DoB intervals. This should not come as a surprise given the known weak dependence of the Nd elastic scattering observables on the off-shell behavior of the NN potentials [64] and a similar structure of the TM and the leading chiral 3NF at N 2 LO [65,66] largely driven by the intermediate Δ(1232) excitation.
Similarly to our findings in Ref. [38] based on the SCS interactions of Refs. [5,6], the nucleon and deuteron vector analyzing powers are also properly described (within errors) at the lowest considered energy of E N lab = 10 MeV showing no evidence of the so-called A y -puzzle [64] at this chiral order. The much larger truncation uncertainty for the vector analyzing powers at this low energy as compared with other observables indicates their strongly fine-tuned nature; see also Ref. [38] for a related discussion. We further emphasize that the approximate subtraction of the Coulomb effects from the proton-deuteron data at this energy may lead to sizable uncertainties.
In Refs. [39,40], we have calculated Nd scattering observables based on the SCS NN chiral potentials of Refs. [5,6] and estimated the truncation errors using the EKM approach. While these calculations are incomplete starting from N 2 LO due to the missing 3NF, they have demonstrated that the expected accuracy of chiral EFT at high orders such as N 4 LO should be substantially smaller than the observed discrepancies between state-of-the-art calculations and experimental data. Figure 8 shows an update of these finding by using the new SMS NN potentials of Ref. [7], including the 3NF at N 2 LO and replacing the EKM approach to estimating truncation errors by the Bayesian modelC   Open circles are proton-deuteron data from Ref. [67]. For remaining notation, see Fig. 4 generated by the pion loop contributions to the 3NF at N 3 LO [17]. Since we do not have complete results beyond N 2 LO, the error bands in Fig. 8 are obtained by just rescaling the corresponding 68% and 95% N 2 LO DoB intervals. The incomplete N 3 LO and N 4 LO + results may, of course, still be regarded as complete N 2 LO predictions. At E N lab = 200 MeV, the N 3 LO uncertainty bands are still quite sizable indicating that the N 4 LO contributions to the 3NF could play a significant role. Thus, fully in line with the findings of Ref. [39,40], our results suggest that the accurate description of Nd scattering data below pion production threshold will likely require the chiral expansion of the 3NF to be pushed to N 4 LO. Notice that the accurate and precise description of neutron-proton and proton-proton data below pion production threshold also required the chiral expansion of the NN force to be pushed to N 4 LO (or even N 4 LO + ) [7].

Subleading short-range 3NF: an exploratory study
To include the 3NF contributions beyond N 2 LO one needs to regularize the corresponding pion loop contributions consistently with the NN interactions of [7] in a chirally symmetric manner; see Refs. [35,36] for discussion. Such consistently regularized pion-exchange contributions to the 3NF are not yet available beyond N 2 LO. In addition to long-and intermediate-range interactions generated by pion-exchange diagrams, the chiral 3NF involves ten purely short-range operators at N 4 LO, which have been worked out in Ref. [32]. In the exploratory study of Ref. [68], the effects of these subleading short-range terms are investigated in proton-deuteron scattering below E p lab = 3 MeV within the hybrid approach based on phenomenological two-and threenucleon forces. The authors of Ref. [68] have succeeded to fit the coefficients of the short-range operators to obtain a good description of experimental data at E p lab = 3 MeV. However, except for fine-tuned observables like the neutron-deuteron doublet scattering length, which is well known to be correlated with the 3 H binding energy, and the analyzing powers A y and i T 11 ; see Fig. 4, the scattering observables at such low energies are dominated by the NN interaction, and the 3NFs are expected to play a minor role [64]. On the other hand, large discrepancies between theory and data are observed in Nd elastic scattering at intermediate and higher energies, where the 3NFs are expected to play a prominent role [69].
To explore the role of the subleading short-range 3NF contributions we choose two out of the ten terms, namely the isoscalar central and spin-orbit interactions where K i = (p i + p i )/2, while E 1 and E 7 denote the corresponding LECs. We apply the same nonlocal Gaussian regulator as employed in the N 2 LO short-range part of the 3NF and defined in Eq. (13), and restrict ourselves to the cutoff Λ = 450 MeV. In this study we do not attempt to determine the LECs E 1 and E 7 from data but explore effects of these 3NF terms for fixed values of these LECs. Specifically, in a complete analogy with Eq. (11), we express E i in terms of dimensionless coefficients c E i according to For the coefficients c E i , we consider in this study the fixed values of c E i = ±2. The adopted naturalness estimates are, however, subject to convention-dependent ambiguities 4 and should be taken with care. A more meaningful and reliable assessment of the natural size of the LECs can be carried out in the spectroscopic basis as done in Ref. [7] for the NN potentials. This, however, would require the inclusion of a complete set of independent contact operators in the 3NF at N 4 LO, which goes beyond the scope of our study.
For a given observable, the impact of the subleading shortrange 3NF terms can be quantified by comparing the results for c E i = ±2 with those for c E i = 0 after renormalization. Since we are only able to perform implicit renormalization by expressing the bare LECs in terms of low-energy observables, this requires a readjustment of the LECs c D and c E for every considered set of c E 1 , c E 7 . To allow for a meaningful interpretation of the obtained results, we follow in each case exactly the same fitting procedure as explained in Sect. 3 by adjusting c D and c E to the 3 H binding energy and the Nd cross section minimum at E N lab = 70 MeV. In Fig. 9 we show the resulting dependence of c D and c E on c E 1 and c E 7 .
In addition to looking at the absolute values of the various dimensionless LECs, it is also instructive to compare the corresponding expectation values in the triton state, which are listed in Table 1. For the considered SMS regulator and Λ = 450 MeV, the expectation value of the two-pion exchange 3NF contributions amounts to V 2π = −0.19 MeV. 5 Notice that the spin-orbit 3NF term ∝ c E 7 does not contribute to the S-wave partial waves in the triton state (for the employed angle-independent regulator), and is found to provide a negligible contribution to the 3 H binding energy. The apparent contradiction with the findings of Ref. [68] regarding this term is presumably caused by a different regulator employed in that paper. The expected natural contribution of the 3NF 4 For example, our notation for the spin-orbit term ∝ E 7 in the 3NF differs from the one adopted in Ref. [68] by a factor of 1/2. 5 We emphasize again that the regularized two-pion exchange contributions contain admixtures of the short-range c D -and c E -like terms. When adopting the convention of Ref. [7] by explicitly subtracting the shortrange pieces of the two-pion exchange, the expectation value changes to V 2π = −0.61 MeV, showing that the N 2 LO 3NF is actually dominated by the long-range pieces.
can be estimated based on naive dimensional analysis via For c E 7 = ±2, the individual terms in the 3NF already start exceeding their expected natural size, thus indicating that the considered values for this LEC likely overestimate its natural range. This conclusion is supported by Nd scattering results, as described below.
We are now in a position to discuss the effects of the subleading short-range contributions to the 3NF in selected Nd scattering observables. To that aim, we first perform calculations based on the NN SMS potential of Ref. [7] at N 4 LO + together with the 3NF at N 2 LO.
The resulting predictions for the Nd elastic scattering observables lie in the middle of the blue bands as shown in Figs. 10, 11 for the intermediate energy of E N lab = 135 MeV. The light-and dark-shaded blue bands show the 95% and 68% DoB intervals for the truncation error at N 3 LO. These error bands do not properly reflect the uncertainty of our calculation, which is only complete through N 2 LO, but show the expected size of N 4 LO corrections estimated within the employed Bayesian approach. Next, we repeat the calculations by switching on the N 4 LO short-range terms in the 3NF. The resulting predictions for c E 1 = ±2 and c E 7 = 0 (c E 1 = 0 and c E 7 = ±2) are shown in Figs. 10, 11 by the purple dashed-dotted (blue dashed) lines. As expected from the estimated truncation uncertainty at N 3 LO, the considered N 4 LO 3NF terms yield sizable contribution to the Nd scattering observables at this rather high energy, especially in the region of the cross section minimum and at backward angles. In fact, the magnitude of the c E 1 = ±2 corrections compares well with the width of the N 3 LO error bands, especially with the ones corresponding to 95% DoB intervals. This finding supports our expectation that the actual value of this LEC should be well within the considered range of c E 1 = ±2. On the other hand, the contributions of the c E 7 -term lie in most cases outside of the N 3 LO truncation bands, which suggests that the employed values of c E 7 = ±2 overestimate the natural size of this LEC. This conclusion is in line with the pattern shown in Fig. 9, which indicates a significantly larger shifts in the LECs c D , c E induced by changing δc E 7 = ±2 as compared with the ones induced by δc E 1 = ±2. Another interesting observation is that both the c E 1 and c E 7 contributions tend to lie well within the truncation bands at forward angles while outside at backward angles. This could indicate a shortcoming of the employed Bayesian model, which relies on Eq. (2) and does not explicitly account for higher momentum scales being probed in backward scattering as compared with forward scattering; see also Ref. [51] for a similar conclusion.
At the lowest considered energy of E N lab = 10 MeV, the effects of the considered N 4 LO 3NF terms turn out to be small; see the right panel of Fig. 12 for a representative example, except for the nucleon and deuteron vector analyzing powers A N y and A d y . These results provide yet another confirmation of the fine-tuned nature of these observables, see the discussion in Sect. 3, and indicate that the apparent A ypuzzle could be naturally resolved at the N 4 LO level by the corresponding short-range contributions to the 3NF [70]; see also Ref. [71] for a related discussion within pionless EFT. While the spin-orbit 3NF is well known to have a strong impact on the vector analyzing power [72], we found, quite surprisingly, that the isoscalar central short-range 3NF term ∝ c E 1 also significantly affects A y . A comparison of effects due to the c E i -terms with the estimated truncation error bands at E N lab = 10 MeV leads to the same conclusions as reached at the higher energy of E N lab = 135 MeV. Last but not least, we have also looked at the differential cross section in the so-called symmetric space star con-figuration of the Nd breakup reaction at E N lab = 13 MeV, which is known to represent another low-energy puzzle in the three-nucleon continuum. Contrary to A y , this observable is dominated by the S-wave components of the NN force and turns out to be highly insensitive to the 3NFs considered so far; see Ref. [73] for recent results based on the SCS nuclear potentials of Refs. [5,6,38]. We found that the inclusion of the c E i -terms has a negligible effect on the cross section in this breakup configuration, and the obtained results essentially coinside with the one presented in Ref. [73]. The observed discrepancy between the theoretical calculations and the neutron-deuteron and proton-deuteron data thus indeed appears to be puzzling. It will be interesting to see if this puzzle can be resolved by the inclusion of the N 3 LO and remaining N 4 LO contributions to the 3NF.

Summary and conclusions
In this paper we analyzed selected Nd scattering observables at N 2 LO in chiral EFT based on the SMS interactions of Ref. [7]. The main results of our study can be summarized as follows.
-Following the approach of Ref. [41], we have explored several pointwise Bayesian models for quantifying truncation uncertainties in chiral EFT and tuned them by calculating angular distributions of neutron-proton scattering. -Using the SMS NN forces of Ref. [7] accompanied with the N 2 LO 3NF regularized in the same way, we have determined the LECs c D and c E entering the 3NF from the 3 H binding energy and the Nd cross section data of Ref. [62] at E N lab = 70 MeV. The resulting N 2 LO results for elastic Nd scattering observables agree within errors with our earlier N 2 LO calculations based on the SCS interactions [38] and with experimental data. -The truncation errors for various Nd scattering observables estimated in [38][39][40] using the approach of Ref. [5] are found to be consistent with 68% DoB intervals for the employed Bayesian modelC 650 0.5−10 . In particular, we confirm our earlier findings in Refs. [39,40], obtained using the SCS NN forces of Refs. [5,6], that Nd scattering at intermediate and higher energies provides a "golden window" to study higher-order contributions to the chiral 3NF.
-Based on the estimated truncation errors, we argue that a high-precision description of neutron-deuteron scattering below pion production threshold will likely require pushing the chiral expansion of the 3NF at least to N 4 LO. This conclusion is in line with the convergence pattern of the chiral expansion in the NN sector [7] and is further substantiated by an exploratory study of the short-range 3NF contributions ∝ E 1,7 , which, as expected from our Bayesian analysis, are found to have significant effect on Nd polarization observables at intermediate energies.
Our results show that the contributions to the 3NF beyond N 2 LO are indeed potentially capable of resolving the discrepancies between theory and experiment observed in Nd scattering at intermediate and higher energies [69].
As a next step, this study should be extended to N 3 LO, which will require the inclusion of consistently regularized 3NF contributions. Work along these lines is in progress by the LENPIC Collaboration.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: All relevant data are already given in the figures and in the table.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Analytic expressions for the posterior pdf
whereq 2 ≡ k+h i=k+1 Q 2i , c 2 k ≡ i∈A c 2 i and the incomplete gamma function is defined as For the non-informative prior C withc < = andc > = 1/ , the expression for the posterior, after taking the limit → 0 for c 2 k = 0 simplifies to