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

TL;DR

This paper reviews the progress in high-order chiral effective field theory calculations for three-nucleon scattering, demonstrating good agreement with data up to third order and discussing future challenges for fifth order accuracy.

Contribution

It provides selected results for nucleon-deuteron scattering using regularized chiral potentials and explores the impact of higher-order three-nucleon forces on observables.

Findings

Good description of experimental data at third chiral order

Higher-order effects are significant for polarization observables

Pushing to fifth order is necessary for high-precision modeling

Abstract

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.

Tables

Table 1. 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 Λ 450 \Lambda=450  MeV.

$(c_{E_1},c_{E_7})$	$⟨V_D⟩$	$⟨V_E⟩$	$⟨V_{E_1}⟩$	$⟨V_{E_7}⟩$
$(0,\mathrm{\hspace{0.17em}0})$	$0.46$	$-0.53$	$0$	$0$
$(2,\mathrm{\hspace{0.17em}0})$	$-0.52$	$-0.06$	$0.69$	$0$
$(-2,\mathrm{\hspace{0.17em}0})$	$1.16$	$-0.30$	$-0.69$	$0$
$(0,\mathrm{\hspace{0.17em}2})$	$-1.37$	$1.43$	$0$	$0$
$(0,-2)$	$1.42$	$-1.32$	$0$	$0$

Equations

$$X$$

$$Q=\max \left(\frac{p}{{\Lambda }_b},\frac{M_{\pi }^{\mathrm{eff}}}{{\Lambda }_b}\right),$$

$$X_{\mathrm{ref}}=X^{(0)}$$

$$X_{\mathrm{ref}}=\max \left(|X^{(0)}|,\frac{|\Delta X^{(2)}|}{Q^2}\right)$$

$$X_{\mathrm{ref}}=\max \left(|X^{(0)}|,\frac{|\Delta X^{(2)}|}{Q^2},\frac{|\Delta X^{(3)}|}{Q^3}\right)$$

$${\mathrm{pr}}_h(\Delta |\{c_{i\le k}\})=\frac{{\int }_0^{\infty }d\bar{c}{\mathrm{pr}}_h(\Delta |\bar{c})\mathrm{pr}(\bar{c}){\prod }_{i\in A}\mathrm{pr}(c_i|\bar{c})}{{\int }_0^{\infty }d\bar{c}{\prod }_{i\in A}\mathrm{pr}(c_i|\bar{c})},$$

$${\rm pr}_{h}(\Delta|\bar{c})\equiv\left[\prod_{i=k+1}^{k+h}\int_{-\infty}^{\infty}dc_{i}{\rm pr}(c_{i}|\bar{c})\right]\,\delta\bigg{(}\Delta-\sum_{j=k+1}^{k+h}c_{j}Q^{j}\bigg{)}\,.$$

$$\mathrm{pr}(c_i|\bar{c})=\frac{1}{\sqrt{2\pi }\bar{c}}{e}^{-c_i^2/(2{\bar{c}}^2)}$$

$$\mathrm{pr}(\bar{c})=\frac{1}{\ln ({\bar{c}}_{>}/{\bar{c}}_{<})}\frac{1}{\bar{c}}\theta (\bar{c}-{\bar{c}}_{<})\theta ({\bar{c}}_{>}-\bar{c}).$$

$$V^{3N}$$

$$D=\frac{c_D}{F_{\pi }^2{\Lambda }_{\chi }},E=\frac{c_E}{F_{\pi }^4{\Lambda }_{\chi }},$$

$$\frac{1}{{\vec{q}}_i^2+M_{\pi }^2}\to \frac{1}{{\vec{q}}_i^2+M_{\pi }^2}{e}^{-\frac{{\vec{q}}_i^2+M_{\pi }^2}{{\Lambda }^2}}.$$

$$V_{\mathrm{cont}}^{3N}\to V_{\mathrm{cont}}^{3N}{e}^{-\frac{4{\vec{p}}_{12}^2+3{\vec{k}}_3^2}{4{\Lambda }^2}}{e}^{-\frac{4{{\vec{p}}_{12}^{\prime }}^2+3{{\vec{k}}_3^{\prime }}^2}{4{\Lambda }^2}},$$

$${\vec{p}}_{\mathrm{CM}}^2=\frac{E_{\mathrm{lab}}^N{A}^2{m}^2(E_{\mathrm{lab}}^N+2m)}{m^2(A+1{)}^2+2Am{E}_{\mathrm{lab}}^N}\simeq \frac{2A^2}{(A+1{)}^2}m{E}_{\mathrm{lab}}^N,$$

$$\sqrt{s}-(A+1)m=:2\sqrt{{\vec{p}}^2+m^2}-2m,$$

$${\vec{p}}^2=\frac{s-2(A-1)m\sqrt{s}+(A+1)(A-3)m^2}{4}$$

$${\vec{p}}^2=\frac{A}{A+1}m{E}_{\mathrm{lab}}^N.$$

$$\delta c_1$$

$$\delta c_3=-\delta c_4$$

$$V_{3N}$$

$$E_i=\frac{c_{E_i}}{F_{\pi }^4{\Lambda }_{\chi }^3}.$$

$$\displaystyle{\rm pr}_{h}^{C}(\Delta|\{c_{i\leq k}\})=\frac{1}{\sqrt{\pi\bar{q}^{2}\mbox{\boldmath$c$}_{k}^{2}}}\left(\frac{\mbox{\boldmath$c$}_{k}^{2}}{\mbox{\boldmath$c$}_{k}^{2}+\Delta^{2}/\bar{q}^{2}}\right)^{k/2}$$

$$\displaystyle\times\frac{{\rm\Gamma}\left[\frac{k}{2},\,\frac{1}{2\bar{c}_{>}^{2}}\left(\mbox{\boldmath$c$}_{k}^{2}+\frac{\Delta^{2}}{\bar{q}^{2}}\right)\right]-{\rm\Gamma}\left[\frac{k}{2},\,\frac{1}{2\bar{c}_{<}^{2}}\left(\mbox{\boldmath$c$}_{k}^{2}+\frac{\Delta^{2}}{\bar{q}^{2}}\right)\right]}{{\rm\Gamma}\left[\frac{k-1}{2},\,\frac{\mbox{\boldmath$c$}_{k}^{2}}{2\bar{c}_{>}^{2}}\right]-{\rm\Gamma}\left[\frac{k-1}{2},\,\frac{\mbox{\boldmath$c$}_{k}^{2}}{2\bar{c}_{<}^{2}}\right]}\,,$$

$$\mathrm{\Gamma }(s,x)={\int }_x^{\infty }dt{t}^{s-1}{e}^{-t}.$$

$${\rm pr}_{h}^{C_{\epsilon}}(\Delta|\{c_{i\leq k}\})=\frac{1}{\sqrt{\pi\bar{q}^{2}\mbox{\boldmath$c$}_{k}^{2}}}\frac{{\rm\Gamma}\left(\frac{k}{2}\right)}{{\rm\Gamma}\left(\frac{k-1}{2}\right)}\left(\frac{\mbox{\boldmath$c$}_{k}^{2}}{\mbox{\boldmath$c$}_{k}^{2}+\Delta^{2}/\bar{q}^{2}}\right)^{k/2}\,.$$

$$V_{3N}=E_1\sum _{i\ne j\ne k}{\vec{q}}_i^2,$$

$$⟨p^{\prime }{k}^{\prime }{\alpha }^{\prime }|V_{3N}^{(1)}|pk\alpha ⟩=32{\pi }^2{E}_1{\delta }_{s^{\prime }s}{\delta }_{l^{\prime }0}{\delta }_{l0}{\delta }_{s{j}^{\prime }}{\delta }_{sj}{\delta }_{T^{\prime }T}{\delta }_{M_T^{\prime }{M}_T}$$

$$\times {\delta }_{t^{\prime }t}\left[(k^2+k^{\prime 2}){\delta }_{{\lambda }^{\prime }0}{\delta }_{\lambda 0}{\delta }_{I^{\prime }\frac{1}{2}}{\delta }_{I\frac{1}{2}}-\frac{2}{3}k{k}^{\prime }{\delta }_{{\lambda }^{\prime }1}{\delta }_{\lambda 1}{\delta }_{I^{\prime }I}\right],$$

$$V_{3N}=i{E}_7\sum _{i\ne j\ne k}{\vec{q}}_i\times ({\vec{K}}_i-{\vec{K}}_j)\cdot ({\vec{\sigma }}_i+\overrightarrow{{\sigma }_j}),$$

$$V_{3N}^{(1)}$$

$$⟨p^{\prime }{k}^{\prime }{\alpha }^{\prime }|V_{3N}^{(1)}|pk\alpha ⟩=8{\pi }^2{E}_7{\delta }_{T^{\prime }T}{\delta }_{M_T^{\prime }{M}_T}{\delta }_{t^{\prime }t}$$

$$\displaystyle\times\bigg{[}-\sqrt{2}\left(1-(-1)^{s^{\prime}+s}\right)(-1)^{J+\tfrac{1}{2}}$$

Full text

1 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. Epelbaum:2008ga ; Epelbaum:2012vx ; Machleidt:2011zz for review articles. The fifth-order (N4LO) contributions to the nucleon-nucleon (NN) force have been worked out in Ref. Entem:2014msa and a new generation of accurate and precise chiral EFT NN potentials up through N4LO has been developed in Refs. Epelbaum:2014efa ; Epelbaum:2014sza ; Reinert:2017usi ; Entem:2017gor , see also Refs. Ekstrom:2013kea ; Gezerlis:2013ipa ; Piarulli:2014bda for related developments. The N4LO interactions of Refs. Reinert:2017usi ; Entem:2017gor utilize the values of the pion-nucleon ($\pi$N) low-energy constants (LECs) determined from matching chiral perturbation theory to the Roy-Steiner-equation analysis of $\pi$N scattering at the subthreshold point Hoferichter:2015tha , but differ substantially by the regularization procedure. The potentials developed by our group in Ref. Reinert:2017usi build upon our earlier studies Epelbaum:2014efa ; Epelbaum:2014sza and employ a local regulator for pion-exchange contributions which, per construction, maintains the analytic structure of the long-range interaction. Differently to the nonlocally regularized potentials of Ref. Entem:2017gor , the interactions constructed in Ref. Reinert:2017usi do not produce long-range artifacts at any finite order of expansion in inverse powers of the momentum-space cutoff $\Lambda$. The resulting semilocal111The term “semilocal” refers to a local regularization approach for long-range interactions in combination with a nonlocal cutoff for the short-range part of the nuclear force. momentum-space regularized (SMS) potentials of Ref. Reinert:2017usi are currently the most precise chiral EFT interactions and provide, at the highest order N4LO+222The N4LO+ potentials of Ref. Reinert:2017usi include four sixth-order (i.e. N5LO) short-range operators contributing to F-waves, which are needed to describe certain very precisely measured proton-proton scattering observables at intermediate and high energies. The same operators are included in the N4LO potentials of Ref. Entem:2017gor ., a nearly perfect and $\Lambda$-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 Perez:2013jpa . 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 Wiringa:1994wb , CD Bonn Machleidt:2000ge and Nijm I, II Stoks:1994wp 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 (N2LO) in the chiral expansion. What makes high-accuracy calculations of such systems/reactions beyond N2LO 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 (N2LO) and fourth (N3LO) 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 Bernard:2007sp ; Bernard:2011zr , 4NF Epelbaum:2006eu ; Epelbaum:2007us and electroweak charge and current operators Kolling:2009iq ; Kolling:2011mt ; Krebs:2016rqz ; Krebs:2019aka completely up through N3LO using dimensional regularization (DR) to compute pion loop contributions, see also Refs. Pastore:2008ui ; Pastore:2009is ; Piarulli:2012bn ; Baroni:2015uza for a related work on electroweak currents by the JLab-Pisa group and Ref. Hermann_review for a review article. Furthermore, selected N4LO contributions to the 3NF have been worked out in Refs. Krebs:2012yv ; Krebs:2013kha ; Girlanda:2011fh , and the longest-range part of the 3NF was also analyzed in the framework involving $\Delta(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. Golak:2009ri ; Hebeler:2015wxa . 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 N3LO) and beyond CD18-Hermann ; EE-NTSE . 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. Slavnov:1971aw in the NN sector. Work along these lines is in progress.

In this paper we update our recent study Epelbaum:2018ogq based on the semilocal coordinate-space regularized (SCS) NN forces of Refs. Epelbaum:2014efa ; Epelbaum:2014sza and analyze nucleon-deuteron (Nd) scattering observables using the new SMS chiral NN potentials of Ref. Reinert:2017usi in combination with the SMS 3NF at N2LO. We also refine our previous estimations of truncation uncertainties by employing a Bayesian approach instead of the algorithm proposed in Ref. Epelbaum:2014efa and used in our earlier studies Binder:2015mbz ; Binder:2018pgl ; Epelbaum:2018ogq in the three-nucleon sector. Last but nor least, we explore the role of selected N4LO 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. Melendez:2017phj . Our results for Nd scattering observables at N2LO are presented in section 3, while the role of selected N4LO 3NF operators is discussed in section 4. The main results of our study are summarized in section 5.

2 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. Epelbaum:2014efa , 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. Siemens:2016hdi ; Yao:2016vbz ; Siemens:2017opr ; Blin:2018pmj and Epelbaum:2014sza ; Binder:2015mbz ; Hu:2016nkw ; Skibinski:2016dve ; Binder:2018pgl ; Epelbaum:2018ogq ; NevoDinur:2018hdo , 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. Furnstahl:2015rha ; Melendez:2017phj , 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. Epelbaum:2014efa , the EKM error estimations for the neutron-proton (np) total cross section at selected energies were found in the Bayesian approach of Ref. Furnstahl:2015rha 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. Furnstahl:2015rha . 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,\ldots$, (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

[TABLE]

with $\Delta X^{(2)}:=X^{(2)}-X^{(0)}$ and $\Delta 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_{\rm 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 Epelbaum:2014efa ; Epelbaum:2014sza ; Furnstahl:2015rha ; Melendez:2017phj

[TABLE]

where $\Lambda_{b}$ is the breakdown scale of the chiral expansion. The quantity $M_{\pi}^{\rm eff}$ serves to model the expansion of NN observables around the chiral limit. In Refs. Epelbaum:2014efa ; Epelbaum:2014sza ; Furnstahl:2015rha ; Melendez:2017phj , this scale was set to the physical pion mass $M_{\pi}$. However, as pointed out in Ref. CD18-EE , the error plots in Ref. Epelbaum:2014efa indicate that the transition between the two expansion regimes in the NN sector actually appears at a scale $M_{\pi}^{\rm eff}$ higher than $M_{\pi}$. On the other hand, both Bayesian model parameters $\Lambda_{b}$ and $M_{\pi}^{\rm eff}$ can be determined/tuned empirically by calculating the success rates for a given set of observables and/or energies. In particular, Ref. Furnstahl:2015rha confirmed the EKM estimation $\Lambda_{b}\sim 600$ MeV based on the results for the total np cross section and using the SCS potentials of Ref. Epelbaum:2014efa , but it also found somewhat larger values of $\Lambda_{b}$ to be statistically consistent, see also a related discussion in Melendez:2019izc . A similar empirical analysis was performed in Ref. CD18-EE for both $\Lambda_{b}$ and $M_{\pi}^{\rm eff}$ yielding the values of $\Lambda_{b}\sim 650\ldots 700$ MeV and $M_{\pi}^{\rm eff}\sim 200$ MeV.

Suppose the results for the observable $X(p)$ are available up through the order $X^{(k)}$, $k\geq 2$. The goal is then to estimate the truncation error $\delta X^{(k)}\equiv\sum_{i>k}\Delta X^{(i)}$ resulting from neglecting the unknown higher-order contributions, i.e. to compute the posterior probability distribution function (pdf) for $\delta X^{(k)}$ given the explicit knowledge of $\{X^{(0)},\,X^{(2)},\,\ldots,X^{(k)}\}$. The Bayesian model of Ref. Melendez:2017phj uses the leading-order (LO) result for $X(p)$ to set the overall scale

[TABLE]

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. Melendez:2019izc . In such circumstances it is advantageous to set the overall scale from the next-to-leading order (NLO) contribution $\Delta X^{(2)}$ via $X_{\rm ref}=\Delta X^{(2)}/Q^{2}$ in order to avoid underestimating $X_{\rm ref}$. To have an approach applicable when both $X^{(0)}$ and/or $\Delta X^{(2)}$ are accidentally small, we set the scale $X_{\rm ref}$ via

[TABLE]

for $k=2$ and

[TABLE]

for $k\geq 3$. Let $c_{m}=1$, $m\in\{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 ${\rm pr}(c_{i}|\bar{c})$ with a hyperparameter $\bar{c}$ and performing marginalization over $h$ chiral orders $k+1,\ldots,k+h$, which are assumed to dominate the truncation error, the probability distribution for the dimensionless residual $\Delta_{k}\equiv\sum_{n=k+1}^{\infty}c_{n}Q^{n}\simeq\sum_{n=k+1}^{k+h}c_{n}Q^{n}$ to take a value $\Delta_{k}=\Delta$, given the knowledge of $\{c_{i\leq k}\}$, is given by Melendez:2017phj

[TABLE]

where the set $A$ is defined as $A=\{n\in\mathbb{N}_{0}\,|\,n\leq k\,\land\,n\neq 1\,\land\,n\neq m\}$ and

[TABLE]

Here and in what follows, we employ the Gaussian prior of “set C” from Ref. Melendez:2017phj , namely

[TABLE]

and assume a log-uniform probability distribution Schindler:2008fh

[TABLE]

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 ${\rm pr}_{h}(\Delta|\{c_{i\leq k}\})$ from Ref. Melendez:2017phj . The posterior pdf is an even function of $\Delta$, and for any given DoB interval, the corresponding value of the residual $\Delta_{k}$ and the resulting truncation error $\delta X^{(k)}=X_{\rm ref}\Delta_{k}$ can be readily obtained by numerically integrating ${\rm pr}_{h}(\Delta|\{c_{i\leq k}\})$ over $\Delta$.

As a first application, we employ the Bayesian model of Ref. Melendez:2017phj and set the overall scale $X_{\rm 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_{\rm lab}=143$ MeV for the non-informative prior $C_{\epsilon}$ Furnstahl:2015rha with $\bar{c}_{<}=\epsilon$ and $\bar{c}_{>}=1/\epsilon$, which makes no assumption about either the maximum and minimum size of $\bar{c}$ by taking the limit $\epsilon\to 0$, see Eq. (23). Here, the breakdown scale was set to $\Lambda_{b}=700$ MeV, and the resulting Bayesian model is referred to as $C_{\epsilon}^{700}$.

Being noninformative, the prior set $C_{\epsilon}$ is generally expected to yield conservative estimates for truncation errors. However, setting $\bar{c}_{<}\to 0$ yields a $\delta$-function-like posterior ${\rm pr}_{h}^{C_{\epsilon}}(\Delta|\{c_{i\leq k}\})$ for \mbox{\boldmathc}_{k}^{2}\to 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 $\Delta^{(2)},\;\ldots,\;\Delta^{(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\sigma/d\Omega$ around $\theta_{\rm CM}\sim 85^{\circ}$, for the polarization transfer $D_{t}$ around $\theta_{\rm CM}\sim 140^{\circ}$, for the spin correlation coefficients $C_{kp}$ at $\theta_{\rm CM}\sim 65^{\circ}$ and for the coefficient $C_{kk}$ at $\theta_{\rm CM}\sim 25^{\circ}$, and $\theta_{\rm CM}\sim 115^{\circ}$. In all these cases, the corresponding functions $\Delta X^{(2)}(\theta_{\rm CM})$ change their sign. To circumvent the problem with the underestimation of the truncation error in such kinematical regions, the authors of Ref. Melendez:2017phj suggested to use a more informative (but not too restrictive) prior set $C_{0.25-10}$ corresponding to $\bar{c}_{<}=0.25$ and $\bar{c}_{>}=10$. Here and in what follows, we make the choice $\bar{c}_{<}=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_{0.5-10}^{650}$ indeed yields reasonable estimates of the truncation errors for $d\sigma/d\Omega$. However, the more informative prior with $\bar{c}_{<}\neq 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 $\theta_{\rm CM}\sim 100^{\circ}$ and for $C_{kk}$ at $\theta_{\rm CM}\sim 10^{\circ}$, $\theta_{\rm CM}\sim 75^{\circ}$ and $\theta_{\rm CM}\sim 180^{\circ}$, 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)}=\alpha\tilde{X}^{(0)}$ with $\alpha\to 0$ being a dimensionless parameter, one finds $\Delta_{k}\sim\alpha^{-1}$ for the prior set $C_{\epsilon}$ while $\Delta_{k}\sim\alpha^{0}$ leading to $\delta X^{(k)}=\alpha\tilde{X}^{(0)}\Delta_{k}\sim\alpha$ for the prior set $C_{\bar{c}_{<}-\bar{c}_{>}}$ with $\bar{c}_{>}<\infty$. The problem can be easily fixed by replacing Eq. (3) with Eq. (4) as shown in the second row of Figs. 1 and 2. Here and in what follows, the resulting Bayesian model is referred to as $\tilde{C}$. 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\geq 3$, Eq. (4) with Eq. (5) and refer to the resulting Bayesian model as $\bar{C}$. 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_{\pi}^{\rm eff}=M_{\pi}$ and $\Lambda_{b}=600$ MeV).

As pointed out in Ref. Furnstahl:2015rha , 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_{\rm lab}=50$ MeV with the models $C$ and $\bar{C}$ as compared to the original EKM approach are mainly due to the larger value of $M_{\pi}^{\rm 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\ldots 75\%$ ($100\%$), which are statistically consistent with the DoB intervals of $68\%$ ($95\%$).333For the EKM approach, the success rate equals $100\%$ per construction. We have verified that this conclusion also holds true for a larger set of energies considered in Ref. CD18-EE . 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. Melendez:2019izc , which is based on Gaussian processes and encodes correlation structure of coefficients $c_{i}(\theta_{\rm CM})$.

3 Nucleon-deuteron scattering at N2LO using SMS nuclear potentials

We now turn to the main topic of our study and consider nucleon-deuteron scattering in chiral EFT. The N2LO three-nucleon force is given by vanKolck:1994yi ; Epelbaum:2002vt

[TABLE]

where the subscripts refer to the nucleon labels, $\vec{q}_{i}=\vec{p}_{i}\,^{\prime}-\vec{p}_{i}$ with $\vec{p}_{i}\,^{\prime}$ and $\vec{p}_{i}$ being the final and initial momenta of the nucleon $i$ and $\vec{\sigma}_{i}$ (\mbox{\boldmath\tau}_{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_{\pi}$ 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

[TABLE]

where, following Epelbaum:2002vt ; Epelbaum:2018ogq , we use $\Lambda_{\chi}=700$ MeV. For the LECs $c_{i}$, we employ the central values from the Roy-Steiner-equation analysis of Ref. Hoferichter:2015tha 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. Reinert:2017usi at N2LO. Differently to Ref. Epelbaum:2018ogq , we employ the same semilocal momentum-space regulator for the 3NF as in the NN potentials of Ref. Reinert:2017usi by replacing the pion propagators via

[TABLE]

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(-(\vec{p}_{12}^{2}+{\vec{p}_{12}^{\prime}}^{2})/\Lambda^{2})$, where $\vec{p}_{12}=(\vec{p}_{1}-\vec{p}_{2})/2$, ${\vec{p}}_{12}^{\,\prime}=(\vec{p}_{1}^{\,\prime}-\vec{p}_{2}^{\,\prime})/2$. For the contact interaction proportional to the LEC $E$, we apply a nonlocal Gaussian regulator in momentum space

[TABLE]

where $\vec{k}_{3}=2(\vec{p}_{3}-(\vec{p}_{1}+\vec{p}_{2})/2)/3$ and $\vec{k}_{3}^{\,\prime}=2(\vec{p}_{3}^{\,\prime}-(\vec{p}_{1}^{\,\prime}+\vec{p}_{2}^{\,\prime})/2)/3$ are the corresponding Jacobi momenta.

It is important to emphasize that the SMS NN potentials of Ref. Reinert:2017usi 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. Epelbaum:2018ogq , we determine the LECs $c_{D}$ and $c_{E}$ from the 3H binding energy and the nucleon-deuteron differential cross section minimum at $E_{\rm lab}^{N}=70$ MeV. Specifically, we fit to the experimental data in the range of $\theta_{\rm CM}=107\ldots 141^{\circ}$. The resulting $c_{D}$- and $c_{E}$-values are strongly dependent on the cutoff $\Lambda$ and change from $c_{D}=8.9$ ($c_{E}=1.15$) for $\Lambda=400$ MeV to $c_{D}=-5.4$ ($c_{E}=-0.25$) for $\Lambda=550$ MeV. Such a strong cutoff dependence has to be expected given that these LECs refer to bare quantities. The calculated observables, on the other hand, show only a weak residual cutoff dependence consistent with the estimated truncation uncertainty.

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 N2LO at energies of $E_{\rm lab}=10$ MeV, $E_{\rm lab}=70$ MeV and $E_{\rm lab}=135$ MeV for the cutoff $\Lambda=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 model $\bar{C}_{0.5-10}^{650}$ 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_{\rm lab}^{N}$ 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_{\rm lab}^{N}$ via

[TABLE]

where $m$ is the nucleon mass and “$\simeq$” refers to the nonrelativistic approximation. Identifying the scale $p\equiv|\vec{p}|$ in Eq. (2) with $p_{\rm CM}\equiv|\vec{p}_{\rm CM}|$ results in $A$-dependent values of the expansion parameter $Q$ corresponding to the same excess energy. For example, the pion production threshold in the NN (Nd) system with $E_{\rm lab}^{N}\sim 290$ MeV ($E_{\rm lab}^{N}\sim 215$ MeV) corresponds to $p_{\rm CM}\sim 370$ MeV ($p_{\rm CM}\sim 425$ MeV), leading to the expansion parameter of $Q=0.57$ ($Q=0.65$). Alternatively, one can define the momentum scale $p$ in terms of the Lorentz-invariant excess energy $\sqrt{s}-\sqrt{s_{0}}=\sqrt{s}-(A+1)m$ available in the $A+1$-nucleon system and define the momentum scale $p$ via the relation

[TABLE]

that ensures that $p$ coincides with $p_{\rm CM}$ in the NN system. Here, $s$ is the usual Mandelstam variable. One can thus express the scale $p$ in terms of $E_{\rm lab}^{N}$ via

[TABLE]

with $s=m^{2}(A+1)^{2}+2AmE_{\rm lab}^{N}$. In the nonrelativistic approximation, this relation simplifies to

[TABLE]

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. (2). The breakdown scale $\Lambda_{b}=650$ MeV then corresponds to the excess energy of $\sim 400$ MeV independently of the number of nucleons $A$ in the target nucleus. Notice that the employed model leads to less conservative error estimates for $A>1$ as compared with the assignment of $p=p_{\rm CM}$ in Eq. (2).

We now turn to the results for Nd scattering observables at $E_{\rm lab}^{N}=10\ldots 135$ MeV shown in Figs. 4-7. Except for the differential cross section at $E_{\rm lab}^{N}=70$ MeV shown in the upper left panel of Fig. 5, the results at N2LO 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_{\pi}^{\rm eff}$ and $\Lambda_{b}$. While model dependence of uncertainty estimates is expected to decrease at high chiral orders, it may still be significant at N2LO.

While we only show the uncertainty bands for the cutoff value of $\Lambda=500$ MeV, the results for different values of $\Lambda$ are similar. To illustrate this point we plot in Fig. 5 for the intermediate energy of $E_{\rm lab}^{N}=70$ MeV the N2LO results for all available cutoff values in the range of $\Lambda=400\ldots 550$ MeV. Notice that the residual cutoff dependence of the considered observables at NLO is similar to the one at N2LO and is, in most cases, comparable with the N2LO $68\%$ DoB intervals. These results demonstrate that our calculations for different values of $\Lambda$ are consistent with each other within errors. Notice further that the softest cutoff of $\Lambda=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 Machleidt:2000ge with and without the Tuscon-Melbourne (TM) 3NF Coon:2001pv . 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_{\rm lab}^{N}=70$ and $135$ MeV, our N2LO 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 Gloeckle:1995jg and a similar structure of the TM and the leading chiral 3NF at N2LO Epelbaum:2007sq ; Meissner:2008zza largely driven by the intermediate $\Delta$(1232) excitation.

Similarly to our findings in Ref. Epelbaum:2018ogq based on the SCS interactions of Refs. Epelbaum:2014efa ; Epelbaum:2014sza , the nucleon and deuteron vector analyzing powers are also properly described (within errors) at the lowest considered energy of $E_{\rm lab}^{N}=10$ MeV showing no evidence of the so-called $A_{y}$-puzzle Gloeckle:1995jg 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. Epelbaum:2018ogq 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. Binder:2015mbz ; Binder:2018pgl , we have calculated Nd scattering observables based on the SCS NN chiral potentials of Refs. Epelbaum:2014efa ; Epelbaum:2014sza and estimated the truncation errors using the EKM approach. While these calculations are incomplete starting from N2LO due to the missing 3NF, they have demonstrated that the expected accuracy of chiral EFT at high orders such as N4LO 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. Reinert:2017usi , including the 3NF at N2LO and replacing the EKM approach to estimating truncation errors by the Bayesian model $\bar{C}_{0.5-10}^{650}$. Specifically, the incomplete N3LO and N4LO results shown in this figure are based on the N3LO and N4LO+ NN potentials accompanied with the N2LO 3NF with the LECs $c_{D}$ and $c_{E}$ being readjusted to the 3H binding energy and the differential cross section at $E_{\rm lab}^{N}=70$ MeV in exactly the same way as done at N2LO. In the 3NF, we have used the values of the LECs $c_{i}$ from Ref. Hoferichter:2015tha consistent with the NN interactions at the corresponding chiral order, namely $c_{1}=-1.07$ GeV*-1*, $c_{3}=-5.32$ GeV*-1* and $c_{4}=3.56$ GeV*-1* at N3LO and $c_{1}=-1.10$ GeV*-1*, $c_{3}=-5.54$ GeV*-1* and $c_{4}=4.17$ GeV*-1* at N4LO, subject to the additional shifts of

[TABLE]

generated by the pion loop contributions to the 3NF at N3LO Bernard:2007sp . Since we do not have complete results beyond N2LO, the error bands in Fig. 8 are obtained by just rescaling the corresponding $68\%$ and $95\%$ N2LO DoB intervals. The incomplete N3LO and N4LO+ results may, of course, still be regarded as complete N2LO predictions. At $E_{\rm lab}^{N}=200$ MeV, the N3LO uncertainty bands are still quite sizable indicating that the N4LO contributions to the 3NF could play a significant role. Thus, fully in line with the findings of Ref. Binder:2015mbz ; Binder:2018pgl , 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 N4LO. 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 N4LO (or even N4LO+) Reinert:2017usi .

4 Subleading short-range 3NF: An exploratory study

To include the 3NF contributions beyond N2LO one needs to regularize the corresponding pion loop contributions consistently with the NN interactions of Reinert:2017usi in a chirally symmetric manner, see Refs. CD18-Hermann ; EE-NTSE for discussion. Such consistently regularized pion-exchange contributions to the 3NF are not yet available beyond N2LO. In addition to long- and intermediate-range interactions generated by pion-exchange diagrams, the chiral 3NF involves ten purely short-range operators at N4LO, which have been worked out in Ref. Girlanda:2011fh . In the exploratory study of Ref. Girlanda:2018xrw , the effects of these subleading short-range terms are investigated in proton-deuteron scattering below $E_{\rm lab}^{p}=3$ MeV within the hybrid approach based on phenomenological two- and three-nucleon forces. The authors of Ref. Girlanda:2018xrw have succeeded to fit the coefficients of the short-range operators to obtain a good description of experimental data at $E_{\rm lab}^{p}=3$ MeV. However, except for fine-tuned observables like the neutron-deuteron doublet scattering length, which is well known to be correlated with the 3H 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 Gloeckle:1995jg . 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 KalantarNayestanaki:2011wz .

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

[TABLE]

where $\vec{K}_{i}=(\vec{p}_{i}^{\prime}+\vec{p}_{i})/2$ while $E_{1}$ and $E_{7}$ denote the corresponding LECs. We apply the same nonlocal Gaussian regulator as employed in the N2LO short-range part of the 3NF and defined in Eq. (13), and restrict ourselves to the cutoff $\Lambda=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

[TABLE]

For the coefficients $c_{E_{i}}$, we consider in this study the fixed values of $c_{E_{i}}=\pm 2$. Based on the variation of the LEC $c_{E}$, which, depending on the cutoff $\Lambda=400\ldots 550$ MeV, takes the values in the ranges of $c_{E}=-0.3\ldots 1.2$, $c_{E}=-0.3\ldots 0.8$, $c_{E}=-0.6\ldots 0.7$ and $c_{E}=-0.5\ldots 0.7$ when using the the NN potential at orders N2LO, N3LO, N4LO and N4LO+, respectively, we expect that the actual values of $c_{E_{i}}$ should lie well within the range spanned by $c_{E_{i}}=\pm 2$. The adopted naturalness estimates are, however, subject to convention-dependent ambiguities444For example, our notation for the spin-orbit term $\propto E_{7}$ in the 3NF differs from the one adopted in Ref. Girlanda:2018xrw by a factor of $1/2$. 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. Reinert:2017usi for the NN potentials. This, however, would require the inclusion of a complete set of independent contact operators in the 3NF at N4LO, 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_{E_{i}}=\pm 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 section 3 by adjusting $c_{D}$ and $c_{E}$ to the 3H binding energy and the Nd cross section minimum at $E_{\rm lab}^{N}=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 $\Lambda=450$ MeV, the expectation value of the two-pion exchange 3NF contributions amounts to $\langle V_{2\pi}\rangle=-0.19$ MeV.555We 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. Reinert:2017usi by explicitly subtracting the short-range pieces of the two-pion exchange, the expectation value changes to $\langle V_{2\pi}\rangle=-0.61$ MeV, showing that the N2LO 3NF is actually dominated by the long-range pieces. Notice that the spin-orbit 3NF term $\propto 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 3H binding energy. The apparent contradiction with the findings of Ref. Girlanda:2018xrw 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 $|\langle V_{3N}\rangle|\sim Q^{3}|\langle V_{2N}\rangle|\sim 0.3^{3}\times 40\mbox{ MeV}\sim 1$ MeV. For $c_{E_{7}}=\pm 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. Reinert:2017usi at N4LO+ together with the 3NF at N2LO.

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_{\rm lab}^{N}=135$ MeV. The light- and dark-shaded blue bands show the $95\%$ and $68\%$ DoB intervals for the truncation error at N3LO. These error bands do not properly reflect the uncertainty of our calculation, which is only complete through N2LO, but show the expected size of N4LO corrections estimated within the employed Bayesian approach. Next, we repeat the calculations by switching on the N4LO short-range terms in the 3NF. The resulting predictions for $c_{E_{1}}=\pm 2$ and $c_{E_{7}}=0$ ($c_{E_{1}}=0$ and $c_{E_{7}}=\pm 2$) are shown in Figs. 10, 11 by the purple dashed-dotted (blue dashed) lines. As expected from the estimated truncation uncertainty at N3LO, the considered N4LO 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}}=\pm 2$ corrections compares well with the width of the N3LO 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}}=\pm 2$. On the other hand, the contributions of the $c_{E_{7}}$-term lie in most cases outside of the N3LO truncation bands, which suggests that the employed values of $c_{E_{7}}=\pm 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 $\delta c_{E_{7}}=\pm 2$ as compared with the ones induced by $\delta c_{E_{1}}=\pm 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. Melendez:2019izc for a similar conclusion.

At the lowest considered energy of $E_{\rm lab}^{N}=10$ MeV, the effects of the considered N4LO 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_{y}^{N}$ and $A_{y}^{d}$. 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 N4LO level by the corresponding short-range contributions to the 3NF Girlanda:2016neb , see also Ref. Margaryan:2015rzg 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 Kievsky:1999nw , we found, quite surprisingly, that the isoscalar central short-range 3NF term $\propto 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_{\rm lab}^{N}=10$ MeV leads to the same conclusions as reached at the higher energy of $E_{\rm lab}^{N}=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_{\rm lab}^{N}=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. Witala:2019ffj for recent results based on the SCS nuclear potentials of Refs. Epelbaum:2014efa ; Epelbaum:2014sza ; Epelbaum:2018ogq . 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. Witala:2019ffj . 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 N3LO and remaining N4LO contributions to the 3NF.

5 Summary and conclusions

In this paper we analyzed selected Nd scattering observables at N2LO in chiral EFT based on the SMS interactions of Ref. Reinert:2017usi . The main results of our study can be summarized as follows.

• •

Following the approach of Ref. Melendez:2017phj , 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. Reinert:2017usi accompanied with the N2LO 3NF regularized in the same way, we have determined the LECs $c_{D}$ and $c_{E}$ entering the 3NF from the 3H binding energy and the Nd cross section data of Ref. Sekiguchi:2002sf at $E_{\rm lab}^{N}=70$ MeV. The resulting N2LO results for elastic Nd scattering observables agree within errors with our earlier N2LO calculations based on the SCS interactions Epelbaum:2018ogq and with experimental data.

• •

The truncation errors for various Nd scattering observables estimated in Binder:2015mbz ; Binder:2018pgl ; Epelbaum:2018ogq using the approach of Ref. Epelbaum:2014efa are found to be consistent with $68\%$ DoB intervals for the employed Bayesian model $\bar{C}_{0.5-10}^{650}$. In particular, we confirm our earlier findings in Refs. Binder:2015mbz ; Binder:2018pgl , obtained using the SCS NN forces of Refs. Epelbaum:2014efa ; Epelbaum:2014sza , 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 N4LO. This conclusion is in line with the convergence pattern of the chiral expansion in the NN sector Reinert:2017usi and is further substantiated by an exploratory study of the short-range 3NF contributions $\propto 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 N2LO are indeed potentially capable of resolving the discrepancies between theory and experiment observed in Nd scattering at intermediate and higher energies KalantarNayestanaki:2011wz .

As a next step, this study should be extended to N3LO, 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 Nuclear Physics International Collaboration (LENPIC) project and was supported by BMBF (contract No. 05P18PCFP1), by DFG and NSFC (TRR 110), by DFG (project No. 279384907 - CRC 1245), by the Polish National Science Centre under Grants No. 2016/22/M/ST2/00173 and 2016/ 21/D/ST2/01120, by the Chinese Academy of Sciences (CAS) Presidents International Fellowship Initiative (PIFI) (Grant No. 2018DM0034) and by VolkswagenStiftung (Grant No. 93562). Numerical calculations were performed on the supercomputer cluster of the JSC, Jülich.

Appendix A Analytic expressions for the posterior pdf

For the Gaussian prior pdf in Eq. (6), the posterior pdf takes the form

[TABLE]

where $\bar{q}^{2}\equiv\sum_{i=k+1}^{k+h}Q^{2i}$, \mbox{\boldmathc}_{k}^{2}\equiv\sum_{i\in A}c_{i}^{2} and the incomplete gamma function is defined as

[TABLE]

For the noninformative prior $C_{\epsilon}$ with $\bar{c}_{<}=\epsilon$ and $\bar{c}_{>}=1/\epsilon$, the expression for the posterior, after taking the limit $\epsilon\to 0$ for \mbox{\boldmathc}_{k}^{2}\neq 0 simplifies to

[TABLE]

Appendix B Partial wave decomposition of the N$4$LO 3NF contact terms

For the $E_{1}$-term

[TABLE]

we choose the Faddeev component $V_{3N}^{(1)}=2\,E_{1}\,\vec{q}_{1}^{2}$ invariant with respect to the interchange of nucleons $2$ and $3$ and obtain

[TABLE]

with $\vec{p}$, $\vec{p}^{\prime}$, $\vec{k}$ and $\vec{k}^{\prime}$ denoting the corresponding Jacobi momenta.

For the $E_{7}$-term

[TABLE]

we choose the Faddeev component

[TABLE]

and obtain

[TABLE]

with $\hat{X}\equiv 2X+1$. For more details on notation see Ref. Gloeckle:1995jg .