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. [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 pionexchange contributions which, per construction, maintains the analytic structure of the long-range 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 threshold 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 highprecision 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-toleading 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 related work on electroweak currents by the JLab-Pisa group and Ref. [29] for a review article. 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 few-body 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 the 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 [39,40,38] in the three-nucleon sector. Last but nor least, we explore the role of selected N 4 LO short-range 3NF terms in Nd scattering.
Our paper is organized as follows. In section 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 section 3, while the role of selected N 4 LO 3NF operators is discussed in section 4. The main results of our study are summarized in section 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-/manynucleon sectors see Refs. [42,43,44,45] and [6,39,46,47,40,38,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. [49,41], a general Bayesian approach to calculate 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 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,49,41] 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,49,41], 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) 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 i∈A pr(c i |c) , (6) where the set A is defined as A = {n ∈ N 0 | n ≤ k ∧ n = 1 ∧ n = m} and 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 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] with c < = andc > = 1/ , which makes no assumption about either the maximum and 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 noninformative, the prior set C is generally expected to yield conservative estimates for truncation errors. 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 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.
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 as 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 X (0) as X (0) = αX (0) with α → 0 being a dimensionless parameter, one finds   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 partic-ular 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 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 a 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 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 customary to express the LECs D and E in terms of the corresponding dimensionless coefficients via where, following [56,38], we use Λ χ = 700 MeV. For the LECs c i , we employ the central values from the Roy-Steiner-equation analysis of Ref. [12] at the corresponding chiral order, namely c 1 = −0.74 GeV −1 , c 3 = −3.61 GeV −1 and c 4 = 2.44 GeV −1 . 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 the 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 )/Λ 2 ), where p 12 = (p 1 − p 2 )/2, p 12 = (p 1 −p 2 )/2. For the contact interaction proportional to the LEC E, we apply a nonlocal Gaussian regulator in momentum space where k 3 = 2(p 3 − (p 1 + p 2 )/2)/3 and k 3 = 2(p 3 − (p 1 + p 2 )/2)/3 are the corresponding Jacobi momenta.
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 H binding energy and the nucleon-deuteron differential cross section minimum at E N lab = 70 MeV. Specifically, we fit to the experimental data in the range of θ CM = 107 . . . In Figs. 4-7 we show the results for the differential cross section and selected polarization observables in elastic nucleon-deuteron scattering at NLO and N 2 LO at energies of E lab = 10 MeV, E lab = 70 MeV and E lab = 135 MeV 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 section 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 nucleon-nucleus 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   [15] (CD Bonn NN potential in combination with the Tucson-Melbourne 3NF [57]). Open circles are neutron-deuteron from Ref. [58] and proton-deuteron data from Ref. [59,60,61], corrected for the Coulomb effects, see Ref. [56] for details.
where m is the nucleon mass and " " refers to the nonrelativistic approximation. Identifying the scale p ≡ |p| in Eq.
that ensures that p coincides with p CM in the NN system. Here, s is the usual Mandelstam variable. One can thus with s = m 2 (A + 1) 2 + 2AmE 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) to define the expansion parameter parameter Q in Eq.  Fig. 5, the results at N 2 LO can be regarded as parameter-free predictions. It is reassuring to see that the calculated observables are in a reasonably good agreement with the experimental data, which in most cases lie within the 95% DoB intervals. One should, however, keep in mind that the estimated truncation errors depend on the Bayesian model, assumed prior sets and the values of parameters M eff π and Λ b . While model dependence of uncertainty estimates is expected to decrease at high chiral orders, it may still be significant at N 2 LO.
While we only show the uncertainty bands for the cutoff value of Λ = 500 MeV, the results for different values of Λ are similar. To illustrate this point we plot in 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-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 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 finetuned 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  Fig. 6. 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 Ayy, Axz and Axx in elastic nucleondeuteron 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. 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. Fig. 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   Predictions for polarization transfer coefficients K y xx , K y y , K y yy , K y xz , K y xx − K y yy and the induced polarization P y 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. [63]. For remaining notation see Fig. 4.  [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 three-nucleon forces. The authors of Ref. [68] have succeeded to fit the coefficients of the short-range op-erators 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 iT 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 Ei according to For the coefficients c Ei , we consider in this study the fixed values of c Ei = ±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 short-range 3NF terms can be quantified by comparing the results for c Ei = ±2 with those for c Ei = 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 E1 , c E7 . To allow for a meaningful interpretation of the obtained results, we follow in each case exactly the same fitting procedure as explained in section 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 E1 and c E7 . 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 E7 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 can be estimated based on naive dimensional analysis via | V 3N | ∼ Q 3 | V 2N | ∼ 0.3 3 ×40 MeV ∼ 1 MeV. For c E7 = ±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 described below.
We are now in the 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 Table 1. Expectation values in the triton state (calculated using the N 4 LO + NN force alone) of the various short-range terms in the 3NF (in MeV) for the cutoff Λ = 450 MeV.  Fig. 9, which indicates a significantly larger shifts in the LECs c D , c E induced by changing δc E7 = ±2 as compared with the ones induced by δc E1 = ±2. Another interesting observation is that both the c E1and c E7 -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 section 3, and indicate that the apparent A y -puzzle 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 E1 also significantly affects A y . A comparison of effects due to the  Fig. 10. Results for the differential cross section, nucleon and deuteron analyzing powers A n y and A d y as well as deuteron tensor analyzing powers Ayy, Axz and Axx in elastic nucleondeuteron scattering at laboratory energy of E N lab = 135 MeV based on the SMS NN potentials of Ref. [7] at N 4 LO + in combination with the 3NF at N 2 LO using Λ = 450 MeV. Blue light-(dark-) shaded bands show the expected truncation uncertainty for a complete N 3 LO calculation and are obtained by multiplying the N 2 LO truncation error corresponding to 95% (68%) DoB intervals for the modelC 650 0.5−10 with the expansion parameter Q 0.45. Short-dashed-dotted and long-dasheddotted purple lines show the impact of the N 4 LO central shortrange 3NF ∝ cE 1 with cE 1 = −2 and cE 1 = 2, respectively. Similarly, short-dashed and long-dashed blue lines show the impact of the N 4 LO spin-orbit short-range 3NF ∝ cE 7 with cE 7 = −2 and cE 7 = 2, respectively. For remaining notation see Fig. 6.
c Ei -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 configuration 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 Ei -terms has a negligible effect on the cross section in this breakup configuration, and the ob-  Fig. 11. Same as Fig. 12 but for polarization transfer coefficients K y xx , K y y , K y yy , K y xz , K y xx − K y yy and the induced polarization P y in elastic nucleon-deuteron scattering at laboratory energy of E N lab = 135 MeV. For remaining notation see Figs. 6 and 10.
tained results essentially coinside with the one presented in Ref. [73]. The observed discrepancy between the theoretical calculations and the neutron-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 [39,40,38] 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.

Acknowledgments
This study has been performed within Low Energy Nu- 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 noninformative prior C withc < = andc > = 1/ , the expression for the posterior, after taking the limit → 0 for c 2 k = 0 simplifies to B Partial wave decomposition of the N4LO 3NF contact terms For the E 1 -term we choose the Faddeev component V 3N |pkα = 32π 2 E 1 δ s s δ l 0 δ l0 δ sj δ sj δ T T δ M T M T × δ t t (k 2 + k 2 ) δ λ 0 δ λ0 δ I 1 2 δ I 1 2 − 2 3 k k δ λ 1 δ λ1 δ I I , with p, p , k and k denoting the corresponding Jacobi momenta.
For the E 7 -term we choose the Faddeev component