High-precision nuclear forces from chiral EFT: State-of-the-art, challenges and outlook

We review a new generation of nuclear forces derived in chiral effective field theory using the recently proposed semilocal regularization method. We outline the conceptual foundations of nuclear chiral effective field theory, discuss all steps needed to compute nuclear observables starting from the effective chiral Lagrangian and consider selected applications in the two- and few-nucleon sectors. We highlight key challenges in developing high-precision tree-body forces, such as the need to maintain consistency between two- and many-body interactions and constraints placed by the chiral and gauge symmetries after regularization.


INTRODUCTION
Almost 30 years ago, Weinberg put forward his groundbreaking idea to apply chiral perturbation theory (ChPT), the low-energy effective field theory (EFT) of QCD, to the derivation of nuclear interactions [1,2]. This seminal work has revolutionized the whole field of nuclear physics by providing a solid theoretical basis and offering a systematically improvable approach to low-energy nuclear structure and reactions.
So where do we stand today in the implementation of the program initiated by Weinberg? Much has been learned about specific features of the nuclear interactions and currents and about the role of many-body forces from the point of view of the effective chiral Lagrangian, see Epelbaum [3], Epelbaum et al. [4], Machleidt and Entem [5], and Hammer et al. [6] for review articles covering different research areas, while some issues are still under debate [6,7]. Meanwhile, the interactions derived in chiral EFT, sometimes referred to as "chiral forces, " have largely replaced phenomenological potentials developed in the nineties of the last century. They are nowadays commonly used in ab initio nuclear structure calculations, see Epelbaum et al. [8], Piarulli et al. [9], Lonardoni et al. [10], Hagen et al. [11], Gebrerufael et al. [12], and Cipollone et al. [13] for recent examples using a variety of continuum ab initio methods and Epelbaum et al. [14], Elhatisari et al. [15], and Lähde and Meißner [16] for selected highlights from nuclear lattice simulations. With the most recent chiral nucleon-nucleon (NN) potentials [17] providing a nearly perfect description of the mutually consistent neutron-proton (np) and proton-proton (pp) scattering data below pion production threshold from the Granada-2013 database [18], the two-nucleon sector is already in a very good shape. On the other hand, three-nucleon forces (3NF) are much less understood at the quantitative level [19] and constitute an important frontier in nuclear physics [20].
In this article we focus on the latest generation of chiral nuclear forces based on an improved regularization approach [17,21,22], which allows one to maintain the long-range part of the interaction as will be described in section 4.1. We review our recent work along these lines in the two-nucleon sector, describe the ongoing efforts by the Low-Energy Nuclear Physics International Collaboration (LENPIC) toward developing consistent 1 many-body forces and solving the structure and reactions of nuclei, and discuss selected applications. For selected recent studies along these lines by other groups, see Entem et al. [23], Gezerlis et al. [24], Piarulli et al. [25], Ekström et al. [26,27], Li et al. [28], Lynn et al. [29], and Girlanda et al. [30], and references therein.
Our paper is organized as follows. In section 2, we outline the foundations of the employed theoretical framework. Section 3 gives an overview of various methods to derive nuclear forces and currents from the effective chiral Lagrangian. It also summarizes the available results for nuclear potentials derived using dimensional regularization (DR). In section 4, we present the improved semilocal regularization approach, which is utilized in the most accurate and precise NN potentials of Reinert et al. [17]. We also discuss the challenges that need to be addresses to construct consistently regularized 3NFs and exchange current operators beyond tree level, which are not restricted to any particular type of cutoff regularization. Section 5 is devoted to uncertainty quantification in chiral EFT. Selected results for the NN system, three-nucleon scattering and light nuclei are presented in section 6. We conclude with a short summary and outlook in section 7.

THE FRAMEWORK IN A NUTSHELL
Throughout this work, we restrict ourselves to the two-flavor case of the light up-and down-quarks and employ the simplest version of the effective chiral Lagrangian with pions and nucleons as the only active degrees of freedom. Contributions of the (1232) isobar to the nuclear potentials are discussed in Ordonez et al. [31], Kaiser et al. [32], Krebs et al. [33], Epelbaum et al. [34,35], and Krebs et al. [36]. The effective Lagrangian involves all possible interactions between pions and nucleons compatible with the symmetries of QCD and is organized in powers of derivatives and quark (or equivalently pion) masses. Pions correspond to the (pseudo) Nambu-Goldstone bosons of the spontaneously broken axial generators and thus transform nonlinearly with respect to chiral SU(2) L ×SU(2) R transformations. The effective Lagrangian can be constructed in a straightforward way using covariantly transforming building blocks defined in terms of the pion fields [37,38]. All applications reviewed in this paper rely on a non-relativistic treatment of the nucleon fields and make use of the heavy-baryon formalism to eliminate the nucleon mass m from the leading-order Lagrangian. The individual terms in the effective Lagrangian are multiplied by the corresponding coupling constants, commonly referred to as low-energy constants (LECs), which are not fixed by the symmetry and typically need to be determined from experimental data. The most accurate currently available nuclear potentials at fifth order in the chiral expansion, i.e., at N 4 LO, require input from the following effective Lagrangians (with each line containing the contributions with a fixed number of the 1 The precise meaning of consistency of many-body forces is defined in sections 2 and 4.2. nucleon fields) L eff = L (2) π (M π , F π ) + L (4) π (l 1,...,7 ) + L (1) πN (g A ) + L (2) πN (m, c 1,..., 7 where M π and F π are the pion mass and decay constant 2 , g A is the nucleon axial-vector coupling while l i , c i , d i , e i , C i , D, D i , E, and E i are further LECs. The superscript n of L (n) denotes the number of derivatives and/or M π -insertions and is sometimes referred to as the chiral dimension. Notice that we only show new LECs that appear in the corresponding Lagrangians and suppress the dependence on the LECs appearing at lower orders. The pionic Lagrangian can be found in Gasser and Leutwyler [39], L π N is given in Bernard et al. [40], and Fettes et al. [41], L (0) NN was introduced in Weinberg [1,2], L (2) NN can be found in Epelbaum et al. [42], and Girlanda et al. [43], the minimal form of L (4) NN is given in Reinert et al. [17], L (1) π NN and L (0) NNN are discussed in Epelbaum et al. [44] while L (2) NNN was constructed in Girlanda et al. [45]. Notice further that the chiral symmetry breaking terms ∝ M 2 π are not shown explicitly in L NN and L NNN . For calculations at the physical value of the quark masses, their contributions are absorbed into the LECs listed in Equation (1). We have, furthermore, restricted ourselves in this equation to isospin-invariant terms for the Lagrangians involving two and three nucleons. The single-nucleon Lagrangian L π N does involve isospin-breaking contributions due to the quark mass difference and can be extended to include virtual photon effects [46,47]. The ellipses in the second-to-last line of Equation (1) refer to higher-order Lagrangians L π NN , which have not been worked out yet and would be needed to finalize the derivation of the 3NF at N 4 LO.
The long-range parts of the nuclear forces emerge from pion exchange diagrams and can be derived from L π and L π N . Fortunately, only a very restricted set of (linear combinations of) LECs from these Lagrangians contributes to the πN → πN and πN → ππN scattering amplitudes, which enter as subprocesses when deriving the long-range nuclear interactions up to N 4 LO, namely c 1,...,4 from L (2) π N , d 1 + d 2 , d 3,5,18 , and d 14 − d 15 from L (3) π N and e 14,...,18 from L (4) π N . Here, we made use of the fact that the contributions from the LECs l 3 , e 19,..., 22 , and e 35,..., 38 can be absorbed into the appropriate shifts of the LECs c i [48]. All these πN LECs can nowadays be reliably extracted by matching the πN scattering amplitude from the recent Roy-Steiner equation analysis [49] with ChPT at the subthreshold point [50], see also Siemens et al. [51] for an alternative strategy. Thus, the long-range nuclear interactions are completely determined by the spontaneously broken approximate chiral symmetry of QCD and experimental/empirical information on the πN system in a parameter-free way. The two-and three-nucleon interactions in the last two lines of Equation (1) parameterize the shortrange part of the nuclear forces, and the corresponding LECs have to be determined from NN scattering and three-or morenucleon observables.
In the single-nucleon sector, the effective Lagrangian L π + L π N can be used to systematically compute the scattering amplitude in perturbation theory by applying the chiral expansion, a simultaneous expansion in particles' external threemomenta p ≡ | p | and around the chiral limit M π → 0. The importance of every Feynman diagram is estimated by counting powers of the soft scales and applying the rules of naive dimensional analysis (NDA). The expansion parameter Q ∈ {p/ b , M π / b } is determined by the breakdown scale b , which may (optimistically) be expected to be of the order of the ρmeson mass 3 . At every order in the chiral expansion only a finite number of Feynman diagrams need to be evaluated. For more details on ChPT in the 1N sector see the review article [53].
Contrary to the 1N case, the NN S-wave scattering amplitude exhibits poles in the near-threshold region corresponding to the bound state (deuteron) and the virtual state in the 1 S 0 channel, which signal the breakdown of perturbation theory. In this context, it was pointed out by Weinberg that the contributions of multi-nucleon ladder diagrams are enhanced compared to the estimation based on the chiral power counting due to the appearance of pinch singularities (in the m → ∞ limit) [1,2]. Weinberg also argued that the nucleon mass needs to be counted as m ∼ 2 b /M π ≫ b in order to formally justify the need to perform a non-perturbative resummation of the ladder contributions. Given that the ladder diagrams are automatically resummed by solving the few-nucleon Schrödinger equation, Weinberg's chiral EFT approach to low-energy nuclear systems, perhaps not surprisingly, resembles the quantum mechanical where i is the Laplace operator acting on the nucleon i. The nuclear potentials V 2N , V 3N , . . . receive contributions from diagrams that cannot be reduced to ladder iterations and are calculable in a systematically improvable way within ChPT.
Among the many attractive features, the approach outlined above allows one to maintain consistency between nuclear forces and exchange current operators which are schemedependent quantities. To illustrate the meaning and importance of consistency consider the Feynman diagram on the left-hand side (l.h.s.) of the equality shown in Figure 1 as an example. The corresponding (on-shell) contribution to the scattering amplitude features both a reducible (i.e., of a ladder-type) and irreducible pieces as visualized in the figure. Reducible contributions to the amplitude are resummed up to an infinite order when solving the Faddeev equation corresponding to Equation (2). In doing so, the diagrams corresponding to its zeroth and first iterations shown in the figure must match 3 An upper bound for b is set by the scale 4πF π emerging from pion loops [52]. the result obtained from the Feynman diagram when taken on the energy shell. The iterative contribution from the first graph on the right-hand side of the depicted equality, however, involves NN and 3N potentials, whose off-shell behavior is scheme dependent. Also the 3NF corresponding to the last diagram is scheme dependent (even on the energy shell) [54], and only a consistent choice of the involved two-and threenucleon potentials guarantees the validity of matching for the scattering amplitude. This can indeed be verified explicitly using the expressions for the 3NFs V Clearly, DR is impractical for a numerical solution of the Abody problem and is usually replaced by cutoff regularization. Renormalization of the Schrödinger equation in the context of chiral EFT is a controversial and heavily debated topic, see Lepage [56], Pavon Valderrama and Ruiz Arriola [57], Nogga et al. [58], Birse [59], Epelbaum and Meißner, [60], Epelbaum and Gegelia [61], Long and Yang [62], Valderrama [63], Epelbaum et al. [7], and Hammer et al. [6] for a range of opinions. The essence of the problem is related to the nonrenormalizable nature of the Lippmann-Schwinger (LS) equation for NN potentials truncated at a finite order in the chiral expansion. Except for a few cases, such as the leading-order (LO) equation in pionless EFT and in chiral EFT in spin-singlet channels, ultraviolet (UV) divergences emerging from the loop expansion of the scattering amplitude cannot be absorbed into redefinitions of parameters appearing in the truncated potentials [7,64]. The problem can be avoided by treating the one-pion exchange (OPE) and higher-order contributions to the potential in perturbation theory using e.g., the systematic power counting scheme proposed by Kaplan et al. [65], but the resulting approach unfortunately fails to converge (at least) in certain spin-triplet channels [66,67] (see also [68] for a recent discussion). A renormalizable framework with the one-pion exchange potential (OPEP) treated non-perturbatively was proposed in Epelbaum and Gegelia [69] (see also [70]), based on a manifestly Lorentz invariant form of the effective Lagrangian. This approach requires a perturbative inclusion of higher-order contributions to the potential in order to maintain renormalizability (which may lead to convergence issues in some channels [71]) but has not been systematically explored beyond LO yet.
Throughout this work we employ a finite-cutoff version of nuclear chiral EFT in the formulation of Lepage [56], which is utilized in most of the applications available today. This is so far the only scheme, that has been advanced to high chiral orders and successfully applied to a broad range of few-and many-nucleon systems. Below, we briefly summarize the basic steps involved in the calculation of nuclear observables within this framework. In the following sections, all four steps outlined below will be discussed in detail.
i. Derivation of nuclear forces and current operators from the effective chiral Lagrangian. This can be achieved by separating out irreducible contributions to the A-nucleon scattering amplitude that cannot be generated by iterations of the dynamical equation using various methods outlined in section 3. The derivations are carried out in perturbation theory using the standard chiral power counting. In contrast to ChPT for the scattering amplitude, special efforts are needed to arrive at renormalized nuclear potentials. This requires that all UV divergences from irreducible loop diagrams are canceled by the corresponding counter terms. The renormalizability requirement imposes strong constraints on the unitary ambiguity of nuclear forces and currents [48,[72][73][74][75]. ii. Introduction of a regulator for external (off-shell) momenta of the nucleons in order to make the A-body Schrödinger equation well-behaved. Given the lack of counter terms needed to absorb all UV divergences from iterations of the dynamical equation with a truncated potential, the (momentum-space) cutoff must not be set to arbitrarily high values but should be kept of the order of the breakdown scale, ∼ b [7,56,61]. The accessible cutoff window is, in practice, further restricted by the need to avoid the appearance of spurious deeply bound states which provide a severe complication for applications beyond the NN system [76] and a preference for soft interactions in order to optimize convergence of ab initio many-body methods. Given the rather restricted available cutoff window, it is important to employ regulators that minimize the amount of finite-cutoff artifacts, see section 4 for discussion. While the regulator choice for V 2N still features a high degree of ambiguity, maintaining the relevant symmetries and consistency with regularized manybody forces and exchange currents beyond tree level represents a highly nontrivial task [77,78] (see section 4 for an example and discussion). iii. Renormalization of the few-nucleon amplitude by fixing the short-range multi-nucleon interactions from low-energy experimental data (see section 6 for details). This allows one to express the calculated scattering amplitude in terms of observable quantities instead of the bare LECs C S,T ( ), . ., and amounts to implicit renormalization of the amplitude. Notice that in the pion and 1N sectors of ChPT, renormalization is usually carried out explicitly by splitting the bare LECs l i , d i , e i , . . ., into the (finite) renormalized ones and counter terms, e.g., d i = d r i (µ)+R i (µ). Here, µ denotes the renormalization scale while R i are the corresponding counter terms, which diverge in the limit of a removed regulator (i.e., → ∞ in the cutoff regularization or the number of dimensions d → 4 in DR). Such a splitting is not unique as reflected by the scale µ, and the appropriate choice of renormalization conditions is essential to maintain the desired power counting, i.e., to ensure the appropriate scaling behavior of renormalized contributions to the amplitude leading to a systematic and self-consistent scheme (see e.g., [79,80]). In the few-nucleon sector, the non-perturbative resummation of pion-exchange potentials via Equation (2) can only be carried out numerically 5 , which leaves the implicit renormalization outlined above as the only available option. Notice that contrary to the renormalized LECs l r i (µ), d r i (µ), . . ., the bare LECs C S,T ( ), C i ( ), . . ., must be re-determined at every order in the expansion. iv. Estimation of the truncation uncertainty and a-posteriori consistency checks of the obtained results. These include, among others, testing the naturalness of the extracted LECs [17], making error plots for phase shifts as suggested in Lepage [56], and Grießhammer [81], verifying a reduced residualdependence of observables (within a specified cutoff range) upon including higher-order short-range interactions, see e.g., Figure 4 of Epelbaum [82], and confronting the contributions of many-body interactions and/or exchange currents with estimations based on the assumed power counting [83,84]. Our approach to error analysis is outlined in section 5, while selected consistency checks are discussed in section 6.
Before closing this section, several remarks are in order. First, we emphasize that the approach outlined above is applicable at the physical quark masses. Quark mass dependence of nuclear observables can be studied more efficiently in the renormalizable chiral EFT framework of Epelbaum and Gegelia [69,85], see also Baru et al. [86,87], and Lähde et al. [88] for an alternative method. Secondly, the validity (in the EFT sense) of the finite-cutoff EFT formulation outlined above has been demonstrated numerically by means of the error plots [56,89] and analytically [61] for toy-models with long-range interactions. It can also be easily verified in pionless EFT. For the case of exactly known non-singular long-range potentials, the employed approach reduces in the NN sector to the well-known modified effective range expansion [90]. The relation between the choice of renormalization conditions and power counting is discussed within pionless EFT in Epelbaum et al. [91] 6 . That paper provides an explicit example of the choice of subtraction scheme (i.e., renormalization conditions), which leads to a self-consistent EFT approach for two particles with both a natural and unnaturally large scattering length, while respecting the NDA scaling of renormalized LECs. Notice that in all applications reviewed in this article, few-nucleon short-range interactions are counted according to NDA. A number of authors advocate alternative approaches, in particular by inferring the importance of short-range operators from the requirement of -independence of the scattering amplitude at arbitrarily large values of as articulated in detail in Hammer et al. [6]. However, performing the loop expansion of the solution of the LS equation in spin-triplet channels for the resummed OPEP shows that the scattering amplitude is only partially renormalized in spite of the fact that it admits, in some cases, a finite → ∞ limit at a fixed energy [7]. The danger of choosing ≫ b in such partially renormalized non-perturbative expressions is demonstrated using an exactly solvable model in Epelbaum and Gegelia [61].

CHIRAL PERTURBATION THEORY FOR NUCLEAR POTENTIALS
One method to decouple pion-nucleon and purely nucleonic subspaces of the Fock space, thereby reducing a quantum field theoretic problem to a quantum mechanical one, is the unitary transformation (UT) technique. Let η and λ be the projection operators onto the purely nucleonic subspace of the Fock space and the rest, respectively. The time-independent Schrödinger equation can be written in the form where E denotes the eigenenergy of the πN system. The idea is to apply a UT to the Hamilton operator H in order to block diagonalize the matrix on the l.h.s. of Equation (3) leading to The decoupling requirement is given by To construct the UT U we first introduce a Møller operator [92], which is defined by with the requirement = η.
Here, | refers to few-nucleon scattering states below pion production threshold. See Lindgre [93] for a discussion of the properties of the operator . The Møller operator reproduces the original low-energy state out of projected state. By projecting Equation (6) onto the model space η one obtains the identity Using Equation (6), we can write the time-independent Schrödinger equation in the form where H 0 denotes a free Hamiltonian. On the other hand, projecting the original Schrödinger equation Equation (3) onto the model space and applying on the resulting equation the operator , we obtain Subtracting Equation (10) from Equation (9) leads to This way we obtain a non-linear equation for the Møller operator Defining the operator A via = : η+A with A = λAη, as follows from Equations (7) and (8), we rewrite Equation (11) in the form The UT U was parameterized by Okubo [94] in terms of the operator A via The resulting transformed Hamiltonian leads to the effective potential defined via Obviously, the Okubo transformation in Equation (13) is not the only possibility to obtain a block-diagonalized Hamiltonian. On top of the transformation U one can always apply e.g., a UT acting nontrivially on the η-space, thus leaving the Hamiltonian block-diagonal. This freedom has been exploited in a systematic manner to construct renormalizable/factorizable 3NFs and four-nucleon forces (4NFs) in chiral EFT in Epelbaum [72], Bernard et al. [54,55], and Krebs et al. [48,95].
To derive the potential V UT eff from the effective chiral Lagrangian in Equation (1) one needs to solve the non-linear decoupling Equation (12) for the operator A. This can be done perturbatively using NDA [3] to count powers of three-momenta and pion masses, denoted collectively by Q. For the sake of definiteness, we restrict ourselves in the following to nuclear potentials in the absence of external sources. The extension to the current operators is straightforward and discussed in details in Krebs et al. [75]. The irreducible contributions of any connected Feynman diagram scale as Q ν with ν = −2 + i V i κ i , where V i denotes the number of vertices of type i and κ i is the inverse mass dimension of the corresponding coupling constant, κ i = d i + 3 2 n i + p i − 4. Here, d i is the number of derivatives and/or M π -insertions, while n i and p i denote the number of nucleon and pion fields, respectively 7 . This particular form of the power counting allows one to formulate the chiral expansion in the form that is completely analogous to the expansion in powers of coupling constants. It is thus particularly well suited for algebraic approaches such as the method of UT. Once the operator A is available, one can perform the chiral expansion of Equation (14) to construct the effective potential order-by-order.
The chiral expansion of the nuclear forces is visualized in Figure 2. Below, we briefly discuss isospin symmetric contributions starting from the leading order (LO) Q 0 . The only contributions at this order emerge from the OPEP and two contact interactions ∝ C S,T [1,2]. The first corrections at order Q 2 (NLO) involve the leading two-pion exchange potential (TPEP) [31,96,97] and 7 short range interactions ∝ C i . At order Q 3 (N 2 LO), further corrections to the TPEP ∝ c i need to be taken into account [96]. At the same order one has the first non-vanishing contributions to the 3NF. They are given by the two-pion exchange diagram involving the LECs c i and two shorter-range tree-level diagrams involving the LECs D and E [44,98]. At order Q 4 (N 3 LO), the NN potential receives the contributions from the leading threepion exchange [99][100][101], further corrections to the TPEP [102,103] and 12 new short-range interactions ∝ D i [17]. At the same order, there are various one-loop corrections to the 3NF [54,55,104] and the first contributions to the 4NFs [72,73], which do not involve unknown parameters. Finally, at order Q 5 (N 4 LO), the NN potential receives corrections to the threepion exchange ∝ c i [101] and further contributions to the TPEP [105]. No additional unknown parameters appear in the isospin-conserving part of the NN force at this order. The 3NF also receives corrections at N 4 LO, most of which have already been worked out using DR [45,48,95]. Notice that the 3NF involves at this order a number of new shortrange operators. Work is still in progress to derive the remaining 3NF and 4NF at N 4 LO. We further emphasize that all calculations mentioned above are carried out using DR or equivalent schemes.
The effective potential V UT eff leads, by construction, to the same spectrum and on-shell scattering matrix as the original untransformed potential V [106,107]. There are, however, other possibilities to define the effective potential without changing onshell physics. One example is an energy-independent potential defined by 7 Alternatively (but equivalently), the chiral order ν of a connected, N-nucleon irreducible diagram with L loops can be expressed as The proof that V EI eff of Equation (16) leads to the same spectrum is trivial: where we used Equations (6) and (8) in the first and second lines, respectively. Note that the potential V EI eff is manifestly nonhermitian. However, due to its simplicity, it is widely used in the literature [92]. This example shows that there is a considerable freedom to define nuclear potentials. Nuclear forces and current operators constructed by the Bochum-Bonn group (see e.g., [17, 21, 22, 48, 54, 55, 73-75, 95, 97]), are obtained using the method of UT. The JLab-Pisa group utilizes a different approach by starting with the on-shell transfer matrix T and "inverting" it to obtain the effective potential (see e.g., [108][109][110][111]). This is carried out in perturbation theory by counting the nucleon where the superscripts indicate the chiral order Q n . The same counting scheme is used to organize the contributions to effective potential: The inversion of the LS equation is carried out iteratively to yield Obviously, the knowledge of the on-shell transfer matrix is insufficient to perform the inversion, and one needs to specify its off-shell extension. Notice that the potentials constructed in this way are not necessarily hermitian, and thus there is no guarantee that they are unitarily equivalent to the ones derived using the UT technique. It should, however, always be possible to find a similarity transformation that relates one potential to another. This is exemplified with the potential v (3) 2π (ν = 1) in Equation (20) of Pastore et al. [109], where ν is an arbitrary phase, which is manifestly non-hermitian. Using the similarity transformation in Equation (28) of that paper 8 , it can be transformed to the hermitian potential v (3) 2π (ν = 0), that is actually employed in the current version of the interactions developed by the JLab-Pisa group. With this choice, their potentials are unitarily equivalent to the ones of the Bochum-Bonn group.

Semilocal Momentum-Space Regularization of the NN Potential
In this review article we focus on the semilocal regularization approach of the chiral nuclear potentials carried out in momentum space [17]. For the purpose of regularization we will consider the two-nucleon interaction consisting of two distinct parts: the short-range contact interaction part and the long-range pion-exchange part. In this context, the term "semilocal" refers to the application of a nonlocal regulator for the former and a local regulator for the latter. In particular, the momentumspace matrix elements of the contact potential are multiplied by a simple nonlocal Gaussian regulator Here and in what follows, p ≡ | p | and p ′ ≡ | p ′ |. Such kinds of nonlocal regulators (albeit with different powers of p, p ′ , and ) have been and still are employed as the main method of regularization for the entire potential including the long-range interactions (see e.g., [23,27,[112][113][114][115]) 9 .
However, in Epelbaum et al. [21,22] it was shown that the amount of long-range cutoff artifacts can be significantly reduced by employing a local regulator for pion-exchange potentials. Notice that pion-exchange contributions, except for 9 Notice that the aforementioned potentials (except the one of [112]) additionally employ spectral function regularization (SFR) [116,117] of the TPEP in the form of a sharply cut-off spectral integral in order to suppress its remaining unphysical short-distance behavior. Notice, however, that the application of a nonlocal regulator exp(−(p 2n + p ′2n )/ 2n ) with suitably chosen n is sufficient to arrive at UV-finite iterations of the potential. some relativistic corrections, give rise to local potentials. We require the regulator to preserve the long-range part of the interaction, which is unambiguously determined in chiral EFT. More precisely, for ≫M π , the regulator is required not to affect the large-distance behavior of the n-pion exchange potential V nπ (r) ∼ exp(−nM π r) f (r), with f (r) being an irrational function, up to exponentially small corrections that vanish in the limit → ∞. Inspired by Rijken [118], this is achieved in our momentum-space approach by regularizing the static propagators of pions exchanged between different nucleons with a local Gaussian cutoff via with l = | l| and l denoting the three-momentum of the exchanged pion. The introduction of the Gaussian form factor in the pion propagators leads to properly regularized longrange potentials that are finite at short distances in coordinate space. In order to have a clean separation of the longrange pion-exchange potential from the short-range contact interactions, we made use of the available contact interactions to subtract out the remaining (finite) admixtures of short-range interactions [17]. The fixed coefficients of these subtractions are determined from the requirement that the corresponding coordinate-space potential and as many derivatives thereof as allowed by power counting vanish at the origin. This convention leads to a qualitatively similar regularization as the coordinatespace regulator previously employed in Epelbaum et al. [21,22]. Application of these ideas to the OPEP is straightforward and leads, in the limit of exact isospin symmetry, to where q ≡ | q | ≡ | p ′ − p | and σ i (τ i ) are the Pauli spin (isospin) matrices of the i-th nucleon. Here, the static pion propagator has been regularized according to Equation (22) and a likewiseregularized LO contact interaction has been added to the OPEP.
with erfc(z) denoting the complementary error function, is fixed by the requirement that the spin-spin part of the OPEP in coordinate space vanishes at the origin. For the regularization of the TPEP, we start with a generic three-dimensional loop integral I( q ) arising in the derivation of the TPEP using e.g., the method of unitary transformation as detailed in the previous section or comparable approaches like time-ordered perturbation theory or S-matrix-based methods [96]. As discussed in Rijken [118], the pion energy denominators in the corresponding 1-loop expressions can always be rewritten into an integral over a mass parameter λ involving a product of two static pion propagators with mass M 2 π + λ 2 where l 1 and l 2 denote the three-momenta of the exchanged pions and the ellipses refer to additional momentum-spin-isospin structures arising from the vertices of a particular diagram. With the pion propagators factorized in this a way, we can regularize them by applying the prescription specified in Equation (22) to each of them. Although the introduction of the regulator obviously affects the resulting expression for the TPEP, there is no need to rederive them explicitly. Indeed, the scalar functions accompanying the spin-isospin operators in the unregularized TPEP can be expressed using the dispersive representation with the spectral functions ρ(µ) = ℑ(V 2π (q))| q=0 + −iµ which are readily available up to N 4 LO. For the explicit expressions of the TPEP, additional subtractions of short-range terms have to be performed to arrive at a convergent spectral integral in Equation (26) whose number depends on the chiral order of the contribution at hand. Introducing the pion propagator regulators in Equation (25), the regularized generic spectral integral of Equation (26) is replaced by see Reinert et al. [17] for more details. The resulting potential V 2π 2N, (q) coinsides with the one obtained by explicitly evaluating the loop integral with regularized pion propagators up to a short-range function.
Expanding the exponentials in inverse powers of the cutoff in either Equation (23) or Equation (27), one observes that the regulator indeed does not affect the long-range part of the potential to any order, but generates an infinite series of shortrange terms polynomial in q 2 . Since an increasing number of contact interactions of this form with freely adjustable LECs become available with increasing chiral order, the perturbative restoration of cutoff-independence is also obvious in this scheme.
The expressions of the regularized and subtracted TPEP can be found in Reinert et al. [17]. Here we restrict ourselves to the example of the isospin-independent central part of the leading TPEP at NLO which is given by with the spectral function Two subtractions have been performed in order to render the unregularized spectral integral in Equation (28) convergent and according to our convention, we have additionally fixed the subtraction coefficients C 2 C,1 (µ) and C 2 C,2 (µ) by the requirement that derivative of W (2) C, (r) vanishes at the origin regardless of the subtraction coefficients]. Figure 3 shows the ratio of the regularized and unregularized expressions in Equation (28) in coordinate space. As one can see, the behavior of the regularized potential is smoother when fixing the subtraction coefficients by the convention explained above. Also note that the potential with C 2 C,1 (µ) = C 2 C,2 (µ) = 0 does not vanish at the origin 10 .

Regularization and Consistency of Nuclear Forces
Having defined the regularization scheme in the NN sector, we now turn to regularization of the 3NF. The expressions for the 3NFs described in section 3 have been worked out completely through N 3 LO using DR. They are off-shell consistent with the unregularized NN interactions reviewed in that section in the way explained in section 2. To arrive at regularized 3NFs, it is tempting to apply some kind of multiplicative regulators to the expressions of the 3NF derived using DR. Such a naive approach, however, leads to a violation of the chiral symmetry at N 3 LO and destroys the consistency between two-and three-nucleon forces after regularization.
To illustrate the problem consider the diagrams shown in Figure 1, which have already been discussed in section 2. The 3NF entering the first graph on the right-hand side (r.h.s.) is given by Bernard et al. [55] V 2π , 1/m 3N The complete expression for the relativistic corrections to the 3NF at N 3 LO can be found in Bernard et al. [55]. We now consider the first iteration of V as shown by the first diagram on the r.h.s. of Figure 1. By simply counting the powers of momenta in the loop integration one observes that the loop integral is linearly divergent, which leads to a finite result in DR. As already pointed out in section 2, adding the DR expression for the 3NF V 2π −1π

3N
from Equations (2.16) to (2.20) of Bernard et al. [54] yields (on-shell) the same result as obtained from calculating the Feynman diagram on the l.h.s. of Figure 1 as expected for consistent two-and threenucleon forces.
We now repeat this exercise using the semilocally regularized nuclear potentials in the calculation of the first diagram on the r.h.s. of Figure 2. This leads to where the ellipses refer to all permutations of the nucleon labels and terms finite in the → ∞-limit. The linear divergence ∝ q 3 · σ 3 q 3 · σ 1 is canceled by the D counter term in the second 3NF diagram at N 2 LO in Figure 2. To cancel the linearly divergent contribution ∝ q 2 · σ 2 one would, however, need to introduce a vertex in L (1) π NN corresponding to a derivative-less coupling of the pion to the NN systems. Such vertices violate the chiral symmetry and, being suppressed by powers of M 2 π , cannot appear in L (1) π NN . As a consequence, this linear divergence can not be absorbed into redefinition of the LECs, and the amplitude on the r.h.s. of Figure 1 can seemingly not be renormalized (i.e., made finite in the → ∞ limit). The r.h.s. of the shown equation, therefore, apparently cannot match the (renormalizable) onshell scattering amplitude from the Feynman diagram on the l.h.s.. The problem can be traced back to mixing the DR when calculating the 3NF V 2π −1π 3N with a cutoff regularization for the iterative contributions in Equation (33), see Krebs [77] for another example with the NN axial vector current operator at N 3 LO. Indeed, recalculating the loop integral in V 2π −1π 3N using the cutoff-regularized pion propagators leads to where the ellipses refer to the finite terms in the → ∞-limit. The problematic linear divergence cancels exactly and the agreement with the on-shell amplitude from the Feynman diagram is restored when both consistently regularized contributions on the r.h.s. of Figure 1 are added together.
One may worry whether the regularization issues discussed above could also be relevant for NN interactions. Fortunately, this is not the case since the momentum structure of the NN contact interactions is not restricted by the chiral symmetry. UV divergences emerging from iterations of the LS equation can, therefore, always be absorbed into redefinition of the bare LECs C S,T ( ), C i ( ), . . ..
In the considered example with the 3N amplitude, the consistently regularized 3NF could be obtained by simply recalculating V 2π −1π 3N with all pion propagators being regularized according to Equation (22). This would indeed solve the problem with the cancelation of linear divergencies at N 3 LO, but it would still lead to a violation of the chiral symmetry in diagrams involving three-and four-pion vertices, which depend on the parametrization of the pion field. For vertices involving up to four pion fields, this freedom is represented by a single real parameter α. In the effective chiral Lagrangian, all pion fields are collected in an SU(2) matrix U(π ), whose most general expression, expanded in powers of the pion fields, takes the form Clearly, the on-shell amplitude must be independent of the arbitrary parameter α. Evaluating the 3NF and 4NF with the regularized pion propagators, however, leads to α-dependent expressions (for finite values of ). This shows, perhaps not surprisingly, that the simplistic approach by regularizing all pion propagators as described above violates the chiral symmetry. A possible solution of this problem is provided by the symmetry preserving higher derivative regularization method introduced by Slavnov [119], see also Djukanovic et al. [120] and Long and Mei [121] for recent applications in chiral EFT. To summarize, we have shown that a naive regularization of the three-and more-nucleon forces by multiplying the expressions derived in DR with regulator functions leads to inconsistencies starting from N 3 LO, see Krebs [77] for the same conclusion for two-and more-nucleon charge and current operators. This problem is by no means restricted to semilocal cutoffs. To derive many-body forces and currents regularized consistently with the NN potentials of Reinert et al. [17], the expressions for the 3NF of Bernard et al. [54,55] and Krebs et al. [48], 4NF of Epelbaum [73], and exchange charge and current operators of Kölling et al. [74,122] and Krebs et al. [75,123] need to be recalculated using e.g., an appropriately chosen higher derivative regulator at the level of the effective Lagrangian.

TRUNCATION ERROR ANALYSIS
Estimating the uncertainty associated with truncations of the EFT expansion, which typically dominates the error budget (see section 6), is an important task -in particular since chiral EFT is being developed into a precision tool. In the past, truncation errors were often estimated in few-nucleon calculations from a residual cutoff dependence. This approach, however, suffers from several drawbacks and does not allow for a reliable estimation of truncation errors [113]. In Epelbaum et al. [21], we have formulated a simple algorithm to estimate the size of neglected higher-order terms based on the available information about the EFT expansion pattern for any given observable. To be specific, consider an arbitrary NN scattering observable X at the center of mass momentum p, which is calculated in chiral EFT up to the order Q k The corrections X (i) , X (i) = O X (0) Q i , are assumed to be known explicitly up to the order i = k. The goal is to estimate the size of neglected higher-order terms δX (k) = n>k X (n) . We, furthermore, assume that the expansion parameter Q is given by This simple ansatz is motivated by the expectation that at very low energies, the errors are dominated by the expansion around the chiral limit. The scale M eff π , which will be specified below, is related to the pion mass and controls the convergence rate of the expansion around the chiral limit. At higher energies one would, on the other hand, expect the expansion to be dominated by powers of momenta. This simple picture is in qualitative agreement with the error plots for NN phase shifts [21], which show clearly the two different regimes mentioned above, see Epelbaum [82] for a discussion. It is less clear how to estimate the relevant momentum scale in bound-state observables.
The algorithm proposed by Epelbaum, Krebs and Meißner (EKM) in Epelbaum et al. [21] employs M eff π = M π and b = 600 MeV based on the estimation from the error plots. It also assumes the truncation error δX (k) to be dominated by the first neglected term. The truncation errors at orders Q i , 0 ≤ i ≤ k, are then estimated via (38) subject to the additional constraint which ensures that the estimated errors cannot be smaller than the known actual higher-order contributions. Notice that this relation leads, per construction, to overlapping errors at different orders. In Binder et al. [83], the method was adjusted to make it applicable to incomplete calculations of few-body observables based on NN interactions only. The EKM approach was applied to a broad range of low-energy reactions in the single-baryon [51,[124][125][126] as well as few-and many-nucleon [8,22,[127][128][129] sectors. The robustness of this method and some alternative algorithms are discussed in Binder et al. [130]. The obvious drawback of the EKM approach is that the estimated uncertainties do not offer a statistical interpretation. In Furnstahl et al. [131] and Melendez et al. [132,133], a more general and statistically well-founded Bayesian approach was developed to calculate the probability distribution function (pdf) for truncation errors in chiral EFT. The EKM approach was then shown to correspond to one particular choice of prior probability distribution for the coefficients in the chiral expansion of X(p). In Furnstahl et al. [131], the EKM uncertainties for the np total cross section were found to be consistent with 68% degree-of-belief (DoB) intervals. The authors of that paper, furthermore, found using the semilocal coordinate-space regularized (SCS) potentials of Epelbaum et al. [21,22] the assumed value of the breakdown scale of b = 600 MeV to be statistically consistent for not too soft regulator values, see also Melendez et al. [132] for a related discussion. Recently, a slightly modified version of the Bayesian approach developed in Furnstahl et al. [131] and Melendez et al. [132] was applied by the LENPIC Collaboration to study NN and 3N scattering [84]. Below, we briefly outline the Bayesian modelC 650 0.5−10 proposed in that paper, which will be employed throughout section 6. For more details on the Bayesian approach the reader is referred to the original publications [131,132].
We begin with rewriting Equation (36) in terms of dimensionless expansion coefficients c i by introducing a (generally dimensionfull) scale X ref via where 11 (41) This choice of the reference scale was found in Epelbaum et al. [84] to be more robust for observables that depend on continuously varying parameters, as compared with the choice of X ref = |X (0) | adopted in Melendez et al. [132]. Alternatively, correlations between observables (and thus the coefficients c i ) evaluated at different values of continuously varying parameters can be taken into account using Gaussian processes [133]. Except for the coefficient c m = 1, m ∈ {0, 2, 3}, used to set the scale X ref , the expansion coefficients c i are assumed to be distributed according to some common pdf pr(c i |c) with a hyperparameterc. 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 Melendez et al. [132] pr where the set A is defined as A = {n ∈ N 0 | n ≤ k ∧ n = 1 ∧ n = m} and 11 No meaningful uncertainty estimation can be carried out within the Bayesian approach at LO.
The modelC 650 0.5−10 utilizes the Gaussian prior of "set C" from Melendez et al. [132], for which the integrals in Equation (42) can be performed analytically [132], and uses the values of h = 10,c < = 0.5 and c > = 10. Following Epelbaum [134], the scales that control the expansion parameter are set to M eff π = 200 MeV and b = 650 MeV. The sensitivity of the estimated uncertainties to the choice of prior pdf is discussed in Epelbaum et al. [84], Furnstahl et al. [131], and Melendez et al. [132]. One generally finds minor dependence on the prior pdf if a sufficient amount of information on the coefficients c i is available.

The Two-Nucleon System
We now turn to the calculation of phase shifts and observables in the two-nucleon system. While the derivation and regularization of the nuclear forces have been outlined in the previous sections, we also need to specify the numerical values of all relevant physical quantities and LECs to perform actual calculations. For pion-exchange contributions to the potential, all LECs can be extracted from processes involving at most one nucleon, making it parameter-free. In the TPEP, we use the values of the πN LECs as determined recently by matching the πN scattering amplitude from chiral perturbation theory to a Roy-Steiner equations analysis of πN scattering data at the subthreshold point [50].
We account for the isospin-breaking effects due to the different pion masses in the OPEP and employ the physical masses of the charged and neutral pions M π ± = 139.57 MeV and M π 0 = 134.98 MeV, while we use the isospin-averaged value of M π = 138.03 MeV in the TPEP. We adopt an effective value of g A = 1.29 for the nucleon axial coupling constant which is slightly larger than the current experimental average value of g A = 1.2723(23) [135] because it already accounts for the Goldberger-Treiman discrepancy (see [136] for a related discussion). The employed value of the pion decay constant is F π = 92.4 MeV.
In contrast to the parameter-free long-range potential, the short-range contact interactions in the two-nucleon force have to be determined from experimental NN data. In order to achieve a proper reproduction of pp and, to a lesser extent, np scattering data, it is crucial to also include electromagnetic interactions between the nucleons. Although these interactions are accompanied by powers of a numerically small coupling constant α ∼ 1/137, they are enhanced at low energies and/or forward angles due to the infrared singularity of the photon propagator or, equivalently, due to their long-range nature. Here, we follow the treatment of the Nijmegen group [137] and include the so-called improved Coulomb potential [138], the magneticmoment interaction [139] as well as the vacuum-polarization potential [140] in the calculation of proton-proton observables. The magnetic moment interaction is also taken into account in neutron-proton scattering. To the best of our knowledge, these effects have been included in every partial-wave analysis (PWA) of or fit of a high-quality potential model from NN data since the Nijmegen PWA of Stoks et al. [137], so that differences in their predictions stem from modeling the strong interaction, the experimental input and/or details of the fitting procedure itself.
For scattering data we use the Granada-2013 database [18] which consists of experimental data for NN elastic scattering up to E lab = 350 MeV from 1950 up to 2013 12 . The database contains the data that have been found to be mutually compatible by means of a 3σ rejection criterion in the corresponding phase shift analysis of Navarro Pérez et al. [18]. The presence of very precisely measured proton-proton data in the database, such as those of Cox et al. [142], motivated us to introduce the N 4 LO + version of the potential. As the proper description of these data requires a precise reproduction of F-waves, the N 4 LO + potential adds the four leading F-wave contact interactions formally appearing at N 5 LO and entering the 3 F 2 , 1 F 3 , 3 F 3 , and 3 F 4 partial waves, to the N 4 LO potential. The fits have been performed for all cutoffs = 400, 450, 500, and 550 MeV as well as for all orders from LO up to N 4 LO +13 . When determining the values of the contact LECs, one has to decide up to which energy E lab the experimental data should be taken into account. One is faced with the two competing features: on the one hand, the inclusion of as many data as possible is desirable from a data fitting point of view. On the other hand, the truncation errors for the chiral interactions become larger at high energies. We therefore chose the maximum energy E lab of data to be included to be E max = 260 MeV for N 4 LO and N 4 LO + , while we reduced the energy to E max = 25, 100, 125, and 200 MeV at the orders LO, NLO, N 2 LO, and N 3 LO, respectively. Notice that balancing the tradeoff between these two competing features can be handled using Bayesian methods (see e.g., [143]). From N 3 LO on, we also adjust the deuteron binding energy B d and the coherent neutron-proton scattering length b np to reproduce their experimental values of B d = 2.224575(9) MeV [144] and b np = −3.7405(9) fm [145]. Table 1 shows the reproduction of neutron-proton and proton-proton scattering data in terms of χ 2 /datum values at all considered orders for the cutoff = 450 MeV 14 . We employ a standard definition of the objective function in terms of a sum of squared residuals as detailed in Reinert et al. [17]. As expected, a clear order-by-order improvement in the description of the scattering data can be seen. Table 1 also gives the number of 12 Strictly speaking, our database differs from the one of Navarro Pérez [18] by the omission of the data set from Daub et al. [141] (see [17] for more details). 13 In our paper Reinert et al. [17], also the cutoff = 350 MeV was considered. Given the sizable finite-artifacts for this very soft cutoff choice, we do not consider this case in the following discussion. 14 We have corrected the last figures in the values for χ 2 /datum for np data in the E lab bins of 0-100 and 0-200 MeV at N 3 LO and N 4 LO + quoted in Table 3 of Reinert et al. [17]. The numbers in brackets after the order indicate the number of parameters entering the neutron-proton and proton-proton potentials.
adjustable parameters at each order which also includes isospinbreaking LECs contributing to the 1 S 0 partial wave. It should be noted that no new contact interactions are added when going from NLO to N 2 LO and that the observed improvement of the χ 2 /datum values is entirely due to the N 2 LO contributions to the parameter-free TPEP. A similar situation occurs when going from N 3 LO to N 4 LO, although here we also allow for additional isospin-breaking of the C 1S0 contact LEC splitting it into two independently adjustable parameters for the neutron-proton and proton-proton/neutron-neutron systems. These improvements demonstrate both the importance of the chiral TPEP in the nuclear force and the predictive power of chiral perturbation theory, which allows to use LECs extracted in one process for making parameter-free predictions for (parts of) another. Starting from N 3 LO, a satisfactory description of the neutronproton data in the energy range of E lab = 0 − 200 MeV and the proton-proton data for E lab = 0−100 MeV is achieved. Although the N 4 LO potential improves on this, especially at intermediate and higher energies, it does not achieve a χ 2 /datum ∼ 1 description of the proton-proton data for E lab ≥ 100 MeV. In the intermediate region, this value is significantly affected by the already mentioned high-precision data which requires an accurate description of F-waves. At N 4 LO the differential cross section data set of Cox et al. [142] at E lab = 144.1 MeV, although well described within the Bayesian truncation errors, yields a χ 2 /datum value of 27.88.
The situation is much improved once we switch to the N 4 LO + potential and short-range interactions in F-waves are added. The description of scattering data at higher energies is generally improved and also the high-precision proton-proton data at intermediate energies is accurately reproduced leading to a χ 2 /datum ∼ 1 description of the complete scattering database. Throughout the orders LO − N 4 LO the χ 2 /datum value for proton-proton scattering up to 200 or 300 MeV has been larger than the one for neutron-proton scattering. This is plausible as proton-proton data is in general more precise than neutron-proton data and because only isovector partial waves contribute to it and hence only roughly half of the total number of parameters. However, at N 4 LO + , the reproduction of proton-proton data becomes very accurate while the slightly larger χ 2 /datum values for the neutron-proton data as compared to proton-proton data reflect the larger statistical fluctuations among different data sets. This can be seen as an indication for reaching the threshold where the model accuracy approaches the precision of the data. In fact, the description of the scattering data at N 4 LO + and = 450 MeV is comparable to or exceeds that of the high-quality semi-phenomenological potentials such as CD-Bonn [146], Nijm I, II [147], and Reid93 [147]. Thanks to the parameter-free effects of the TPEP this is achieved with only 27 adjustable short-range parameters instead of the ∼ 40 − 50 parameters used in those potentials.
Indeed, due to the excellent description of the data, the obtained results at = 450 MeV qualify to be considered a partial-wave analysis. In Figures 4, 5, we show the obtained N 4 LO + neutron-proton and proton-proton phase shifts for = 450 MeV, respectively. We compare them to the 2013 Granada analysis [18] and in the case of neutron-proton scattering also to the corresponding 2008 analysis by Gross and Stadler [148]. Furthermore, we also show the predictions from the N 4 LO + potential of Entem et al. [23] at the intermediate cutoff = 500 MeV.
In general, there is good agreement between the shown N 4 LO + phase shifts and the results obtained by the considered phase shift analyses. This is especially true for the case of proton-proton phase shifts which are more strongly constrained by the precise experimental data. Some discrepancies among the different results remain e.g., around the maximum of the 3 P 0 phase shift where the N 4 LO + prediction for the protonproton phase is slightly larger than the ones of the Nijmegen and Granada PWAs, resulting in a ∼ 3σ deviation from the former at E lab = 50 MeV. On the other hand, our neutronproton phase shifts fall in between the results of the two PWAs. The study of isospin-breaking effects in P-waves beyond the ones included in the two PWAs and the current version of the semilocal momentum-space regularized (SMS) interaction of Reinert et al. [17] is expected to shed some light on this issue. We can also compare our results at N 4 LO + to the ones of Entem et al. [23]. Similar to the comparison with the PWAs, agreement with proton-proton phases is better than with neutron-proton ones. There are, however, notable differences in the 3 P 0 , 3 P 2 , and 3 D 2 waves starting at low or intermediate energies. At higher energies around E lab = 250 − 300 MeV, a change in curvature of the phase shift as a function of energy is visible e.g., in the 1 P 1 and 3 P 1 waves, which is presumably caused by the regulator employed in Entem et al. [23]. The effects of regulator artifacts can be observed particularly well in the 1 G 4 , 3 H 4 , and ǫ 4 phase shifts and mixing angle shown in Figure 5 since they do not involve any adjustable short-range parameters at N 4 LO + but are solely determined by the long-range pion-exchange potential. Here, the local regulator of Equation (22) leads to an undistorted reproduction of the peripheral phase shifts.
Selected proton-proton scattering observables and their estimated truncation error at various orders are shown in Figure 6 for E lab around ∼ 143 MeV. In particular, we show our predictions for the differential cross section at E lab = 144.1 MeV and compare them with two high-precision data sets, most notably the one of Cox et al. [142], which motivated the introduction of the N 4 LO + potential as discussed above. The data are well described within the given truncation error for all considered orders, but the N 4 LO + clearly allows for a proper quantitative description. Likewise, the reproduction of the spin observables in Figure 6 is excellent already at N 3 LO with a good convergence pattern. Notice however, that the error bands at lower orders for D (A) at the minimum (maximum) around CM = 150 • do not overlap with the ones for N ≥3 LO and are indeed underestimating the uncertainty. Here we find that the value of the observable in that particular angular region is notably shifted starting at N 3 LO while lower-order corrections are small, such that the overall scale in Equation (41) is still underestimated. Using a more sophisticated Bayesian approach of Melendez et al. [133] would likely allow for a more reliable estimation of the truncation errors at LO-N 2 LO in these particular cases.
There are various a posteriori checks that can be performed to test the self-consistency and quality of the fit. First, the values of the LECs have to be of natural size assuming the cutoff is kept below the hard scale. The expected sizes of the spectroscopic contact LECs can be estimated to be [21] |C where the LECsC i , C i , D i , and E i start to contribute at order Q 0 , Q 2 , Q 4 , and Q 6 , respectively. b is the breakdown scale of the chiral expansion discussed in section 5. Furthermore, the factor of 4π emerges from the angular integration of the partialwave decomposition and has been included in the definition of the spectroscopic LECs. If we now divide the contact LECs obtained in the fit by their expected sizes in Equation (46), we consequently should obtain values of unit magnitude. Figure 7 shows the absolute values of the LECs at N 4 LO + in these natural units for all considered values of the cutoff using b = 650 MeV. As can be seen, all LECs are indeed of natural size with D 1S0 and D 3S1 being among the largest in magnitude. This is especially true for the softest cutoff = 400 MeV, for which also most of the other-Q 4 LECs turn out to be slightly larger than at higher values of the cutoff. This indicates that at = 400 MeV and below, finite-cutoff artifacts start to increase, leading to a lower effective breakdown scale compared to the other considered cutoffs. Notice further that the values for the Q 6 LECs E i included at N 4 LO + turn out to be of a perfectly natural size. Therefore, even though we have emphasized their importance in describing some high-precision proton-proton data and achieving a χ 2 /datum ∼ 1 description of the database, their actual contributions agree with the expectations from naive dimensional analysis (i.e., Weinberg) power counting, and there is no need to promote them to a lower order.
In addition to the absolute of the central values, Figure 7 also shows the statistical uncertainties of the contact LECs as determined from the covariance matrix of the fit (expressed in their natural units). When going fromC i , C i , D i to E i the statistical relative errors tend to increase. This is in accordance with the decreasing importance of higher-order contributions as predicted by power counting. One also notices that errors are smaller for LECs entering isovector partial waves, because these parameters are mainly constrained by the more precise proton-proton data. Since we perform a combined fit of neutronproton and proton-proton data, the isovector partial waves are not only constrained by more precise data but also by more data in general compared to the isoscalar partial waves which have to be extracted from neutron-proton data alone. The covariance matrix also gives access to the correlations among the LECs. As to be expected, correlations mostly occur among LECs entering the same partial waves with the largest ones arising in the channels with the most parameters, namely in the 1 S 0 and 3 S 1 − 3 D 1 channels. Nevertheless, all LECs are well-constrained as can already be seen by looking at the errors in and O i are the experimental value and its error of an observable and O th i is its calculated "theoretical" value. If the assumptions on the normally-distributed residuals can be verified, this confirms that the data are described sufficiently well by the theoretical model. An easy and often employed check is the value of χ 2 per degree of freedom. For the N 4 LO + fit with = 450 MeV we get χ 2 = 4708.65 in the fitting range of E lab = 0 − 260 MeV with the number of data N dat = 4616 and the number of parameters N par = 27. Consequently, we obtain χ 2 /ν = 1.026 with ν = N dat − N par . If the residuals are indeed normal-distributed then χ 2 /ν should follow the χ 2distribution and yields χ 2 /ν = 1 ± √ 2/ν = 1 ± 0.021 as the 68% confidence interval.
We can go one step beyond this simple check and plot the quantiles of the empirical distribution of residuals r i that we obtain against the quantiles of the assumed normal distribution N (0, 1). If they are the same, they should lie on the diagonal line x = y. In order to statistically quantify deviations from the diagonal, confidence bands have been derived with one of the most recent and most sensitive being the ones of the "tailsensitive test" by Aldor-Noiman et al. [154]. This graphical test for normal-distributed residuals has been first applied to the analysis of nucleon-nucleon scattering by Navarro Pérez et al. [155]. Figure 8 shows a rotated quantile-quantile plot for the N 4 LO + residuals at = 450 MeV where the theoretical quantiles have been subtracted from the empirical ones on the y-axis, turning the diagonal line into a horizontal one. As evident from the figure, the empirical distribution of residuals lies within the 68% confidence region of the tail-sensitive test signaling that the residuals are indeed normal-distributed. The quantilequantile plot for the other values of the cutoff turn out to be overall similar, but perform slightly worse. For = 500 MeV and = 550 MeV the quantiles that are already close to the edge of the 68% confidence region in Figure 8 cross these limits but still stay well within the 95% confidence region. The increased cutoffartifacts at = 400 MeV manifest themselves in a stronger deviation from normality as the plotted quantiles also cross the 95% confidence limits at the spike at Q th = 2 in Figure 8.
We now turn to the extended error analysis for observable predictions. While the truncation of the chiral expansion is clearly the dominant source of uncertainty at higher energies, other sources of uncertainty can become relevant at N 4 LO + . In particular we account for the following sources of uncertainty: • Statistical uncertainties of NN LECs: As already mentioned, Figure 7 shows the statistical errors of the contact LECs as determined from the covariance matrix of the fit. The uncertainties of the parameters can then be propagated from the covariance matrix to the observable of interest. While it is always possible to do this via a Monte Carlo sampling of the corresponding multivariate Gaussian probability distribution, it is computationally much more convenient to do a Taylor expansion of the desired observable with respect to the LECs and evaluate the moments of the LECs analytically. While a linear expansion is commonly employed, it has been argued in Carlsson et al. [114], that some observables require a second order expansion for an accurate reproduction of their uncertainties. In the case of large second-order contributions, the error bars become asymmetric and we usually give both the upper and lower error to accommodate for this possibility.
Notice that in such a case, the probability density of the observable is not Gaussian and the quoted uncertainties do not necessarily correspond to a 68% degree-of-belief.   [142] and Jarvis et al. [150]. The data sets have been corrected for their estimated norms of 0.988 and 1.001, respectively. Analyzing power P at E lab = 142 MeV with experimental data taken from Taylor et al. [151]. The data have been floated and multiplied by an estimated norm of 0.942. Depolarization D, rotation parameter A, polarization transfer coefficient D t , and spin-correlation parameter C kp at E lab = 143 MeV with experimental data taken from Bird et al. [152] and Jarvis et al. [153]. The light-(dark-) shaded green, blue, and red bands depict the 68% (95%) DoB truncation errors at N 2 LO, N 3 LO, and N 4 LO + , respectively. Open circles show the predictions of the Nijmegen partial-wave analysis [137]. sets. Due to the need to refit the contact interactions for each sampled set of πN LECs and the computational overhead related to it, we have restricted the total number of such sets to 50. Although this is a quite low statistics for a Monte Carlo approach, it should give an idea of the order of magnitude of the uncertainty. It indeed turns out that the uncertainty related to the statistical error of the πN LECs is small compared to the other sources of uncertainty. However, the aforementioned approach does not probe the systematic errors in the determination of the πN LECs emerging from the truncation of the chiral expansion and thus does not represent the full uncertainty related to these LECs. • Uncertainty due to the choice of the maximum fit energy: The extracted values of the contact LECs also depend on details of the fitting protocol. In particular, we probe the impact of the choice for the maximum laboratory energy E max = 260 MeV up to which scattering data is included in the N 4 LO + fit. This is achieved by performing additional fits with E max = 220 MeV and E max = 300 MeV and determining the maximum deviation of the observables from the E max = 260 MeV predictions. Unlike the aforementioned uncertainties, the error estimated via this simple procedure does not reflect any particular degree-of-belief.
As an example, Figure 9 shows the neutron-proton total cross section and the corresponding uncertainties in the energy range E lab = 0 − 300 MeV. The plot on the left in Figure 9 shows the ratio of our predictions using the N 4 LO + potential at = 450 MeV and the result of the Nijmegen partial-wave analysis [137]. In the right panel, the different relative errors stemming from the various sources discussed above are shown. For the case of the statistical errors of the NN contact interactions, second order effects and resulting asymmetries in the error bands turn out to be small for the total cross section, and the plotted uncertainty corresponds to the average of upper and lower statistical errors. As expected, the dominant contribution to the uncertainty at higher energies (E lab > 100 MeV) arises from the truncation of the chiral expansion. At lower energies, however, other sources of uncertainty become relevant and indeed both the statistical errors of the NN contact LECs and the uncertainty due to the maximum fitting energy are larger than the truncation error in the range of E lab = 30 − 100 MeV. When quantitatively comparing the different errors, one has to keep in mind that the uncertainty due to the maximum fitting energy does not correspond to a particular degree-ofbelief. The uncertainty arising from the statistical errors of the πN LECs is found to be significantly smaller throughout the whole considered energy range and is negligible for the total cross section. Finally, we would like to comment on the origin of the existing kinks in the right-hand-side plot of Figure 9. In particular, the kink in the E max -error at around 200 MeV arises because of the maximum operation. Below 200 MeV, the error is dominated by the deviation of the E max = 220 MeV fit while it is given by the deviation of the E max = 300 MeV fit above 200 MeV. The second kink present in the truncation error, on the other hand, is caused by the transition of Q from M eff π / b to p/ b . Table 2 shows the deuteron properties as predicted by various high-quality potentials. Clearly, the error analysis can also be applied to the bound state properties of Table 2. However, the obtained uncertainties are only meaningful for a complete calculation of an unconstrained observable. This excludes the deuteron binding energy B d (as it is a fitted quantity), the quadrupole moment Q and deuteron radius r d (as meson exchange currents and relativistic corrections are not taken into account) as well as the D-state probability P D (which is not observable). On the other hand, we can perform the uncertainty quantification for the asymptotic S-state normalization A S and the asymptotic D/S-state ratio η for which we obtain at N 4 LO + and for = 450 MeV the values of A s = 0.8847 (+3) (−3) (5)(0)(1) and η = 0.02553 (+11) (−9) (4) (3)(8), respectively. Here the first, second, third, and fourth error refer to the NN statistical, truncation, πN statistical, and E max uncertainty, respectively. Notice that the quoted truncation errors estimated using the Bayesian model of section 5 are fairly similar to the ones given in Reinert et al. [17], which were obtained using the EKM method. On the other hand, the πN statistical uncertainties are much smaller than the corresponding errors quoted in Reinert et al. [17], where an attempt was made to also include systematic effects by using the values of these LECs determined in the physical region of πN scattering.
Finally, let us discuss the treatment of isospin-breaking effects in the two-nucleon interaction. Like all modern high-precision potentials, the SMS interactions include isospin-breaking in the OPEP due to the different physical pion masses M π ± and M π 0 and charge dependence of the short-range potential in the 1 S 0 partial wave. These are the dominant and wellunderstood isospin-breaking effects necessary to arrive at e.g. a correct description of the charge-dependence of the 1 S 0 scattering length. For the calculation of scattering observables in the two-nucleon system, the isospin-breaking due to longrange electromagnetic interactions is taken into account as discussed at the beginning of this section. This treatment of strong and electromagnetic isospin-breaking effects is identical to the Nijmegen PWA [137]. On the other hand, chiral EFT allows for a systematic inclusion of isospin-breaking effects beyond the ones previously considered. In fact, expressions for the leading isospin-breaking TPEP [161,162], the subleading isospin-breaking TPEP [163], and irreducible πγ exchange [164], which are (mostly) parameter-free in the two-nucleon system, have been available for some time. The long-standing question regarding the charge-dependence of the πNN coupling constant also re-emerges in a systematic treatment of isospinbreaking effects in the framework of chiral EFT. While the Nijmegen group did not find evidence for charge-dependence, the issue does not seem to be settled (see [165] for a recent  [156] and have been corrected for their estimated norm of 0.999. The plot on the right shows the relative uncertainties as discussed in the text. determination). Last but not least, charge-dependence in the short-range potential entering P-waves should also be taken into account starting from N 4 LO [163].

Three-Nucleon Scattering
As discussed in the previous subsection, the N 4 LO + SMS potentials of Reinert et al. [17] lead to excellent and in fact a nearly perfect description of np and pp scattering data below pion production threshold. Moreover, an order-by-order comparison of the results for various observables along with the Bayesian error analysis indicate a generally good convergence of the chiral expansion in the NN sector. On the other hand, a description of nucleon-deuteron elastic and breakup scattering data at a comparable level of accuracy is not available yet. Extensive calculations performed in the last decades using highprecision phenomenological NN potentials and 3NF models in the framework of the Faddeev equations [166] and using other ab initio methods [167] have revealed the following picture (see [19] and references therein): -Calculations based on high-precision NN potentials alone (including the N 4 LO + ones of [17]) tend to underestimate the 3 H and 3 He binding energy by ∼ 0.5 MeV and generally lead to similar predictions in Nd scattering observables. -At low energies, the resulting description of Nd data appears to be rather good apart from a few exceptions such as the underprediction of the nucleon analyzing power A y , known as the A y puzzle [168], and the discrepancy for the cross section for the symmetric space star deuteron breakup configuration [169]. 3NF effects in this energy range are found to be small in agreement with qualitative arguments based on the chiral power counting as explained below. -Starting from E lab ∼ 50 MeV, discrepancies between theory and experimental data set in and become large at E lab ∼ 200 MeV and above. Except for the cross section, the inclusion of the phenomenological 3NFs like the Tucson-Melbourne The binding energy has been calculated with the non-relativistic energy-momentum relation for the potentials of Entem et al. [23] and with the relativistic relation for the SMS potential of Reinert et al. [17] and the CD Bonn potential [146]. ⋆ The deuteron binding energy has been taken as input in the fit. † This value corresponds to the so-called deuteron structure radius, which is defined as a square root of the difference of the deuteron, proton and neutron mean square charge radii.
(TM99) [170] and Urbana-IX [171] models does not globally reduce the discrepancies between theory and data [19]. Relativistic effects have also been studied, see Witała et al. [172] and references therein, and found to be small at energies below the pion production threshold.
Assuming that the discrepancies between theory and experimental data in the 3N system are to be resolved by 3NFs, these findings demonstrate that the currently available phenomenological models do not provide an appropriate description of the 3NF. This should not come as a surprise given the enormously rich and complex spin-isospin-momentum structure of a most general 3NF [95,[173][174][175]. Here, chiral EFT offers a decisive advantage over more phenomenological approaches by predicting the long-range part of the 3NF in a model-independent way, establishing a clear importance hierarchy of short-range terms and providing a solid theoretical framework for maintaining consistency between two-and three-nucleon forces and ensuring scheme independence of the calculated observables. As already mentioned in section 3, three-body contributions to the nuclear Hamiltonian first appear at N 2 LO in the chiral expansion and are, therefore, suppressed by Q 3 relative to the dominant pairwise NN interaction. It is instructive to estimate the expected magnitude of 3NF effects for various observables solely on the basis of the chiral power counting (i.e., using NDA). The shallow nature of few-nucleon bound states indicates that there are large cancelations between the kinetic and potential energies. Because of this fine tuning, 3NF contributions to the binding energies are enhanced beyond the naive estimation of Q 3 ∼ 3% and actually reach 10 . . . 15%. On the other hand, there is generally no reason to expect a similar enhancement for Nd scattering observables at low energy except for some fine-tuned polarization observables such as A y . It is well known that tiny changes of the NN interaction in the triplet P-waves amount to large relative changes in the Nd A y [168]. On the other hand, starting from E N ∼ 60 MeV, the expansion parameter Q in Equation (37) is dominated by the momentum scale p = √ 2/3mE N [84]. At e.g. the energies of E N ∼ 100 MeV and E N ∼ 200 MeV, the expansion parameter becomes Q ∼ 0.40 and Q ∼ 0.55, and the relative contributions of the 3NF to a generic scattering observable are expected to increase to ∼ 6 and ∼ 16%, respectively. Clearly, these simplistic back-ofenvelope estimations only yield qualitative insights into the role of the 3NF. Nevertheless, they agree remarkably well with the observed trend of discrepancies between theoretical predictions based solely on the NN interactions and experimental data, which tend to increase with energy. For further examples and a more quantitative analysis along this line of Nd scattering, selected properties of light and medium-mass nuclei and the equation of state of nuclear matter (see [83,127,130,177]). We further emphasize that it is not entirely clear how to estimate the relevant momentum scale, that determines the expansion parameter in heavy nuclei, and how to quantify truncation errors for excited states (see [130] for an extended discussion).
As discussed in section 3 and visualized in Figure 2, the leading contributions to the 3NF at N 2 LO emerge from the two-pion exchange, one-pion-exchange-contact and purely contact tree-level diagrams, leading to the well-known expressions [44,98] where q i = p i ′ − p i with p i ′ and p i being the final and initial momenta of the nucleon i. The LECs D and E are usually expressed in terms of the corresponding dimensionless coefficients c D and c E via D = c D /(F 2 π χ ) and E = c E /(F 4 π χ ) [44]. In Epelbaum et al. [8] and [84], semilocal coordinate-and momentum-space regularized 3NF expressions in combination with the corresponding chiral NN potentials from Epelbaum et al. [21,22] and Reinert et al. [17], respectively, were employed by the LENPIC Collaboration to analyze Nd scattering observables at N 2 LO. The numerical implementation of the 3NF in the Faddeev equations is carried out in the partial wave basis. Partial wave decomposition (PWD) of a general 3NF can be carried out numerically using the machinery developed in Golak et al. [178] by performing five-dimensional angular integrations. Given the required number of partial waves and grid points for the four Jacobi momenta to reach converged results for Nd scattering observables, such a numerical PWD requires substantial computational resources. In Hebeler et al. [179], a more efficient approach was introduced, that exploits the local nature of the bulk of the 3NF.
To make predictions for few-nucleon observables, one first needs to determine the LECs c D and c E entering the 3NF. A broad range of few-and many-body observables including the binding energies and radii of 3 H, 4 He, and heavier nuclei, nucleon-deuteron doublet scattering length 2 a, n-α scattering, triton β-decay, and the saturation properties of nuclear matter have been proposed and employed in the past to determine these two LECs [9,27,29,44,180,181]. A reliable determination of c D , c E is complicated by the existence of strong correlations between some of the low-energy observables (see e.g., [182]), which originate from the large S-wave scattering lengths in the NN system. Furthermore, going beyond the 3N system may require, depending on the observable and the chiral order, the inclusion of 4NF and exchange current contributions. In Epelbaum et al. [8], we therefore, restricted ourselves to 3N observables in the determination of c D , c E . Specifically, we employed the 3 H binding energy of B3 H = 8.482 MeV to fix the LECs c E for a given value of c D . The remaining LEC c D was determined by considering a number of observables including 2 a = 0.645 ± 0.008 fm [145], nd total cross section data from [183] and precisely measured pd differential cross section in the minimum region at E N = 70 MeV [184], 108 MeV [185], and 135 MeV [184]. In the left panel of Figure 10, we show the extracted values of c D for the SCS interactions with the cutoff R = 0.9 fm. It is reassuring to see that the considered 3N observables lead to consistent values of c D . In addition, these results show that the strongest constraint on c D , given the experimental and the estimated truncation uncertainty, is imposed by the pd differential cross section data at E N = 70 MeV from Sekiguchi et al. [184] as visualized in the right panel of Figure 10. We also found no correlations between this observable and the 3 H binding energy. In particular, the resulting value of the LEC c D is largely determined by the differential cross section and almost insensitive to a variation of the triton binding energy.
In Epelbaum et al. [84], we have analyzed Nd scattering observables using the most recent SMS NN potentials from Reinert et al. [17] in combination with the N 2 LO 3NF regularized in the same way. Motivated by the experience with the SCS interactions [8], the LECs c D and c E were determined from the 3 H binding energy and the pd cross section minimum at E N = 70 MeV. In Figure 11, we show, as a representative example, our N 2 LO predictions for selected Nd scattering observables at E N = 135 MeV, along with the experimental data and calculations based on the CD Bonn NN potential with and without the TM99 3NF model. As an important internal consistency check of the calculations, we have verified that the predictions obtained using different cutoff values are consistent with each other (within errors) (see Figure 5 of [84]).
It is reassuring to see that the experimental data are generally well described by the theory. On the other hand, while accurate, our predictions at N 2 LO have obviously rather low precision at this energy. In fact, the N 2 LO truncation errors are comparable with or even larger than the observed deviations between experimental data and calculations based on phenomenological high-precision NN and 3NF models, see the dotted and dasheddotted lines in Figure 11. Based on the experience in the NN sector as discussed in section 6.1, it is conceivable that a high-precision description of Nd scattering data will require the chiral expansion of the 3NF to be pushed to (at least) N 4 LO. At the energy of E N = 135 MeV, the uncertainty bands at N 4 LO are expected to become 4-5 times more narrow as compared with the N 2 LO ones shown in Figure 11.
So where do we stand in terms of efforts to include 3NF corrections beyond N 2 LO? As explained in section 4.2, the main obstacle for the inclusion of higher order contributions to the 3NF is the lack of their consistently regularized expressions. Starting from N 3 LO, it is not sufficient anymore to naively regularize the available expressions for the 3NF from Bernard and Epelbaum [54,55] and Krebs et al. [48,95] derived using DR, since such an approach violates constraints imposed by the chiral symmetry. Rather, the N 3 LO and N 4 LO corrections to the 3NF need to be re-derived using the consistent finite-cutoff regularization approach. Work along these lines is in progress. Another challenge, that will have to be addressed at N 4 LO, is the determination of the LECs appearing in the 3NF at this order. While the N 3 LO contributions do not involve free parameters, the short-range part of the 3NF at N 4 LO depends on 10 unknown LECs [45], from which 9 contribute to the isospin-1/2 channel and thus can, in principle, be determined in Nd scattering. Furthermore, the yet-to-be-derived one-pion-exchange-contact contributions to the 3NF at N 4 LO will also involve unknown LECs. Given the still rather significant computational cost of solving the Faddeev equations in the 3N continuum, the complicated treatment of the Coulomb interaction [186] and the lack of partial wave analyses in the 3N sector, the determination of these LECs from 3N scattering data will certainly be a challenging task.
While a complete analysis of Nd scattering is currently not available beyond N 2 LO, it is instructive to explore the role of subleading short-range 3NF interactions. In Girlanda et al. [30], it was shown within a hybrid phenomenological approach that the 3N contact operators at N 4 LO can be tuned to reproduce the 3 H binding energy, nd scattering lengths, cross section and polarization observables of pd scattering at 2 MeV center-ofmass energy. The resulting models were shown to lead to a satisfactory description of pd polarization observables below the deuteron breakup. On the other hand, 3NF effects are expected to be much more visible at intermediate and higher energies. In Epelbaum et al. [84], we explored the impact of the short-range 3NF operators of the central and spin-orbit types proportional to the LECs E 1 and E 7 , respectively, In order to compute the contributions of the c E i -terms to 3N observables in a meaningful way, one needs to perform (implicit) renormalization as explained in section 2. This was achieved in Epelbaum et al. [84] by simultaneously adjusting the values of the N 2 LO LECs c D , c E to the triton binding energy and the cross section minimum at E N = 70 MeV for all considered values of the LECs c E i . The calculations have been performed using the N 4 LO + SMS NN potential from Reinert et al. [17] in combination with the SMS N 2 LO 3NF. In Figure 12, we show the resulting predictions at the lowest considered energy of E N = 10 MeV. The blue bands show the estimated truncation error at N 3 LO, obtained by rescaling the N 2 LO Bayesian truncation uncertainty with the expansion parameter Q 15 , and visualize the expected impact of N 4 LO terms. In agreement with the expectations, 3NF effects generally appear to be rather small at such low energies. This figure also provides a nice illustration of the fine tuned nature of the nucleon vector analyzing power A y , which shows a strong sensitivity to small changes in the Hamiltonian. What has been traditionally referred to as the A y -puzzle thus appears to be just a consequence of the fine-tuned nature of this observable, and the "puzzle" may be expected to be resolved by 3NF contributions beyond N 2 LO (see also [30,45] for a similar conclusion). While the A y is  [184], while the yellow band also takes into account the quoted systematic uncertainty of 2%. All results are obtained using the N 2 LO SCS NN potential from Epelbaum et al. [21] in combination with the N 2 LO SCS 3NF for the coordinate-space cutoff R = 0.9 fm. FIGURE 11 | Predictions for the differential cross section, nucleon and deuteron analyzing powers A n y and A d y , deuteron tensor analyzing powers A yy , A xz , A xx , 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 Nd scattering at laboratory energy of E N = 135 MeV at NLO (yellow bands) and N 2 LO (green bands). The light-(dark-) shaded bands indicate 95% (68%) DoB intervals using the Bayesian modelC 650 0.5−10 introduced in section 5. Open circles are proton-deuteron data from Sekiguchi et al. [184]. The dotted (dashed-dotted) lines show the results based on the CD Bonn NN potential [146] (CD Bonn NN potential in combination with the Tucson-Melbourne 3NF [170]). All results are obtained using the N 2 LO SMS NN potential from Reinert et al. [17] in combination with the N 2 LO SMS 3NF for the momentum-space cutoff = 500 MeV. FIGURE 12 | Results for the differential cross section, nucleon analyzing powers A n y as well as deuteron tensor analyzing powers A xx and A xx in elastic nucleon-deuteron scattering at laboratory energy of E N lab = 10 MeV based on the SMS NN potentials of Reinert et al. [17] at N 4 LO + in combination with the SMS 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 corresponding value of the expansion parameter Q. Short-dashed-dotted and long-dashed-dotted red lines show the impact of the N 4 LO central short-range 3NF ∝ c E1 with c E1 = −2 and c E1 = 2, respectively. Similarly, short-dashed and long-dashed blue lines show the impact of the N 4 LO spin-orbit short-range 3NF ∝ c E7 with c E7 = −2 and c E7 = 2, respectively. Open circles are neutron-deuteron data from Howell et al. [187] and proton-deuteron data from Sagara et al. [188], Rauprich et al. [189], and Sperisen et al. [190], corrected for the Coulomb effects (see [44] for details).
well-known to be particularly sensitive to spin-orbit types of 3NFs [191] such as the one proportional to c E 7 , our results also show an unexpectedly strong sensitivity to the subleading central interaction of the c E 1 -type.
At higher energies, the effects of the considered N 4 LO 3NF terms become more significant as visualized in Figure 13 for the case of selected spin-correlation parameters. More results for the cross section, vector and tensor analyzing powers and polarization transfer coefficients at E N = 135 MeV can be found in Epelbaum et al. [84]. It is comforting to see that the impact of the c E i -terms on Nd scattering observables is, in general, consistent with the estimated N 3 LO truncation errors. One should, however, keep in mind that the employed Bayesian approach may, under certain circumstances, become unreliable. This is, in particular, the case for observables that depend on a continuously varying parameter in the kinematical regions where the LO results and higher-order corrections change sign (see [84] for a detailed discussion). One such failure of the Bayesian model is shown in Figure 13 for the spin-correlation coefficient C x,z at E N = 200 MeV around θ = 120 • . In such problematic cases, the approach proposed in Melendez et al. [133] and based on Gaussian processes is expected to provide more reliable estimations of the truncation uncertainty.

Light Nuclei
While no results for light nuclei using SMS chiral interactions are available yet, we briefly review here some recent highlights obtained by the LENPIC Collaboration using the SCS NN potentials of Epelbaum et al. [21,22] with and without the corresponding 3NFs at N 2 LO. In Binder et al. [83,130], we have calculated the ground state energies and selected properties of light and medium-mass nuclei up to 48 Ca using the SCS NN interactions at various chiral orders. Specifically, A = 3, 4 nuclei were analyzed in the framework of the Faddeev-Yakubovsky equations while light p-shell nuclei were calculated using the No-Core Configuration Interaction (NCCI) method [193][194][195] and employing Similarity Renormalization Group (SRG) transformed interactions [196][197][198][199] to improve the convergence. The results for 16,24 O and 40,48 Ca were obtained within the coupled cluster and in-medium SRG group frameworks (see [12,[200][201][202][203] and references therein). A qualitatively similar convergence pattern was observed in all considered cases, namely a significant overbinding at LO, results close to the experimental values at NLO and N 2 LO and underbinding at N 3 LO and N 4 LO. Notice that the strongly repulsive nature of the N 3 LO contributions to the SCS NN interactions of Epelbaum et al. [21,22] was shown to be caused by the employed unnaturally large values of the redundant short-range operators [17]. The SMS interactions of Reinert et al. [17] utilize a soft choice for these contact terms, which leads to more perturbative interactions at and beyond N 3 LO. No large gap between the N 2 LO and N 3 LO results for the ground state energies is, therefore, expected for the new SMS NN interactions. The calculated charge radii of the considered medium-mass nuclei were found to show a systematic improvement with the chiral order, but remain underestimated using the NN interaction at the highest available order N 4 LO + .
In Epelbaum et al. [8], a complete N 2 LO analysis of p-shell nuclei was presented by the LENPIC Collaboration using the SCS NN and 3N interactions. In Figure 14, we show the NLO and N 2 LO results from that paper for nuclei up to A = 16. We emphasize that since the Hamiltonian has been completely determined in the NN and 3N system as described in sections 6.1 and 6.2, the ground-state energies shown in that figure are parameter-free predictions. In Figure 14, we have updated the corresponding figure from Epelbaum et al. [8] by replacing the truncation errors, that have been estimated in that paper using the EKM approach of Epelbaum et al. [21] and Binder et al. [130], with the Bayesian uncertainties calculated as described in section 5. The 68% DoB Bayesian truncation errors are similar to those quoted in Epelbaum et al. [8] at N 2 LO but appear to be significantly larger at NLO. We also calculated in that FIGURE 13 | Same as Figure 12 but for the deuteron-nucleon spin-correlation parameters C z,x , C y,y , C z,z , and C x,z for Nd elastic scattering at E N = 135 MeV (left) and E N = 200 MeV (right). Open circles are proton-deuteron data from von Przewoski et al. [192]. paper the excitation energies for selected states of A = 6 − 12 nuclei and the point-proton radius of 4 He. For almost all considered cases, adding the 3NF to the NN interaction was found to lead to a significant improvement in the description of experimental data. The predicted ground state energies of p-shell nuclei show a good agreement with the data except for 16 O, which appears to be overbound. Notice that the deviation between the predicted and experimental values of the 16 O binding energy is comparable to the 95% DoB Bayesian error at N 2 LO. It will be very interesting to repeat the calculations for the newest SMS interactions and to extend them to higher orders and heavier nuclei.

SUMMARY AND OUTLOOK
In this review article we have presented a snapshot of the current state-of-the-art in low-energy nuclear theory with a focus on the latest generation of semilocal nuclear potentials from chiral EFT. We now summarize some of the key conclusions of our paper.
• We have presented a concise and self-contained introduction to the conceptual foundations of chiral effective field theory in the few-nucleon sector and described in some detail all steps needed to compute low-energy observables from the effective chiral Lagrangians (including error analysis). Special emphasis was given to clarify the notion of consistency of nuclear forces and current operators in terms of a perturbative matching to the unambiguously defined on-shell scattering amplitude.
In particular, few-nucleon potentials from Epelbaum [73], Bernard et al. [54,55], Krebs et al. [48,95] and electroweak current operators from Kölling et al. [74,122] and Krebs et al. [75,123] at N 3 LO and beyond, derived using DR, are off-shell consistent with each other provided DR is also used to compute loop integrals arising from iterations of the dynamical equations. • We have reviewed the semilocal momentum-space regularized potentials of Reinert et al. [17], which are currently the most precise chiral EFT NN forces on the market. These are the only NN interactions derived in chiral EFT, whichfrom the statistical point of view-qualify to be regarded as PWA of NN data below pion production threshold (see section 6.1 for details). At the highest considered order N 4 LO + , these interactions describe the np and pp data from the self-consistent Granada-2013 database with a precision that is at least comparable to the one reached by modern phenomenological potentials with a much larger number of adjustable parameters. The significantly better description of the scattering data by the SMS N 4 LO + interactions of Reinert et al. [17] as compared to the nonlocal potentials of Entem et al. [23] at the same chiral order, and their much smaller residual cutoff dependence (see Figure 17 of Reinert et al. [17]), can presumably be traced back to the improved semilocal regulator, which maintains the long-range part of the interaction as described in section 4.1. We also addressed in detail the issue of uncertainty quantification in the NN sector.
In particular, we discussed statistical uncertainties of NN and πN LECs and their propagation to selected observables as well as uncertainty introduced by fixing the maximum fit energy in the determination of the NN LECs. We also estimated truncation errors at various chiral orders using the Bayesian model specified in section 5. • Beyond the NN sector, calculations based on the SMS interactions have so far been carried out up to N 2 LO [84]. The LECs c D and c E , which enter the 3NF at this order, have been determined from the 3 H binding energy and the very precise pd cross section data at E N lab = 70 MeV from Sekiguchi et al. [184]. Using the employed Bayesian model to estimate truncation uncertainties, the predicted ground state energies of p-shell nuclei up to A = 16 are generally in a good agreement with the data. Also the predicted Nd scattering observables including the vector analyzing power A y are consistent with experimental data within errors. We performed an additional test of the employed Bayesian model for truncation errors by exploring the impact of selected short-range 3NF terms at N 4 LO on observables in Nd elastic scattering. Our results suggest that a high-precision description of Nd scattering data will likely require the chiral expansion of the 3NF to be pushed to N 4 LO.
• The novel semilocal nuclear forces, derived in the finitecutoff formulation of chiral EFT with short-range interactions counted according to NDA (i.e., the Weinberg scheme), have already been successfully confronted with few-nucleon data and passed a number of a-posteriori consistency checks as briefly summarized below: -Using the minimal basis of the order-Q 4 NN contact interactions as detailed in section 6, the LECs determined from the np and pp scattering data come out of a natural size (see Figure 7). The same is true for the LECs c D and c E entering the leading 3NF, as can be seen e.g., from the corresponding expectation values in the 3 H state [84]. -The residual cutoff-dependence of NN phase shifts is strongly reduced at N 3,4 LO as compared to N 1,2 LO within the considered -range (see e.g., Figure 4 of [82]). -There is a clear systematic improvement in the description of np and pp data with increasing chiral orders (see Table 1). At order Q 3 (i.e., N 2 LO), this improvement results solely from taking into account the parameterfree subleading TPEP contributions. Notice that certain alternative power counting schemes suggest that some of the contact interactions that appear at order Q 4 in the Weinberg scheme are enhanced and should yield contributions to observables larger than the order-Q 3 TPEP (see e.g., Table 1 of [81]). The clear evidence of the chiral TPEP at orders Q 3 and Q 5 observed in Epelbaum et al. [21,22] and Reinert et al. [17] does, however, not support such alternative scenarios. -The resulting convergence pattern of the EFT expansion for selected NN observables was scrutinized using Bayesian statistical methods (see section 5 for details). For not too soft cutoffs, the assumed breakdown scale of the EFT expansion of b ∼ 600 MeV [21] was found to be statistically consistent [131] (see also [84,132] for a related discussion). -Scheme-dependence of nuclear potentials offers yet another way to perform nontrivial consistency checks of the theoretical framework by explicitly verifying (approximate) scheme-independence of observables. In the formulation we employ, scheme dependence of the nuclear forces first appears at N 3 LO and manifests itself in their dependence on arbitrary real phasesβ 8 ,β 9 , which parameterize the unitary ambiguity of the leading relativistic corrections [21,74], and the appearance of three off-shell short-range operators in the 1 S 0 and 3 S 1 -3 D 1 channels proportional to the LECs D off 1S0 , D off 3S1 , and D off ǫ1 [17,204,205]. The SMS potentials of Reinert et al. [17] make use of the standard choice forβ 8,9 namelyβ 8 = −β 9 = 1/4, which minimizes the amount of 1/m 2 -corrections to the OPEP, and employ D off 1S0 = D off 3S1 = D off ǫ1 = 0. Different choices of these parameters lead to different off-shell behaviors of the potential. They are related to each other by unitary transformations which, however, also induce an infinite tower of higher-order terms beyond the order one is working. The residual dependence of observables onβ 8,9 and D off i , therefore, probes the impact of neglected higher-order terms and should lie within the estimated truncation errors. We have redone the fits at N 4 LO + for = 450 MeV using alternative choices of D off i [17] and also developed a version of the potential with β 8 =β 9 = 1/2 [206]. The letter choice is motivated by the vanishing isoscalar exchange charge density operator at N 3 LO. In all considered cases, we found negligibly small changes in observables in spite of strong changes at the interaction level.
-Calculations of three-and more-nucleon observables based on solely NN interactions are incomplete beyond second order. They do, however, provide information about the magnitude of the missing 3NF contributions by assessing the spread in results at different orders Q ≥3 and via a comparison of such incomplete predictions with experimental data. In Binder et al. [83], such an analysis was performed for Nd scattering observables and selected properties of light nuclei using the SCS NN interactions of Epelbaum et al. [21,22]. The sizes of the 3NF contributions required to bring such incomplete results in agreement with experimental data were found to agree well with expectations based on Weinberg's power counting. Furthermore, recent calculations by the LENPIC Collaboration which include the leading 3NF [8,84] show that the resulting N 2 LO predictions for observables that have not been used in the determination of the LECs c D , c E are generally in a good agreement with the data (see section 6). No indications are found for enhanced contributions of the 3NF in general and of the c D -term in particular as suggested in Birse [207].
To summarize, major progress has been achieved in recent years toward developing chiral EFT into a precision tool for low-energy nuclear physics. In the NN sector, the latest SMS interactions at fifth chiral order have already reached the accuracy at or even below permille level for low-energy observables such as e.g., the deuteron asymptotic S-state normalization A S (see section 6.1 for details and further examples). The only essential missing step in the NN sector concerns the inclusion of isospin-breaking interactions up to fifth chiral order. Work along this line is in progress.
Pushing the precision frontier beyond the NN system opens exciting perspectives for low-energy nuclear theory and will allow one to confront chiral EFT with currently unsolved problems, such as a quantitative description of 3N scattering observables [19]. This, however, will require to address the two core challenges: (i) Derivation of consistent regularized three-and four-nucleon forces and exchange charge and current operators at and beyond N 3 LO as detailed in section 4.2. This issue has not been paid attention to in the recent calculations involving the 3NFs [208][209][210][211] and exchange electroweak currents [110,212,213] at N 3 LO. As explained in section 4.2, using ad hoc regularization approaches at N 3 LO and beyond generally leads to incorrect results for the scattering amplitude and other observables due to the appearance of uncontrolled short-range artifacts, which violate chiral symmetry and are not suppressed by inverse powers of . This puts the findings of these studies into question. (ii) Determination of the LECs in the 3NF at N 4 LO. While the N 3 LO contributions to the 3NF and 4NF do not involve unknown parameters, the N 4 LO corrections to the 3NF involve 10 LEC accompanying purely short-range operators [45] and one or more LECs entering the one-pion-exchangecontact topology, which has not been worked out yet. As discussed in section 6.2, the determination of these LECs from 3N data will require a computationally challenging analysis.
As a first example of a precision calculation not restricted to NN scattering, we have recently determined the deuteron structure radius with an accuracy below the permille level, r str = 1.9731 +0.0013 −0.0018 fm, by pushing the chiral expansion of the electromagnetic exchange charge density beyond N 3 LO and performing a thorough analysis of various types of uncertainty [206]. By combining the predicted value for r str with the very accurate atomic data from isotope shift measurements, it was, for the first time, possible to extract the neutron charge radius from experimental data on light nuclei. This study was facilitated by the absence of loop contributions in the isoscalar exchange charge density at N 3 LO [74,122], which allowed for a trivial construction of the corresponding consistently regularized expressions for the charge operator. Rederivation of the contributions to 3NFs, 4NFs and exchange currents at and beyond N 3 LO using a regulator, consistent with the one employed in Reinert et al. [17], would open the way for performing similar precision calculations for a broad class of low-energy few-nucleon reactions.

AUTHOR CONTRIBUTIONS
EE has lead the development of this review. All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.