High-Precision Nuclear Forces From Chiral EFT: State-of-the-Art, Challenges, and Outlook

1. 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¹ 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.

2. 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⁴LO, require input from the following effective Lagrangians (with each line containing the contributions with a fixed number of the nucleon fields)

where M_π and F_π are the pion mass and decay constant², 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}}_{\pi N}$ is given in Bernard et al. [40], and Fettes et al. [41], ${{L}}_{\text{NN}}^{(0)}$ was introduced in Weinberg [1, 2], ${{L}}_{\text{NN}}^{(2)}$ can be found in Epelbaum et al. [42], and Girlanda et al. [43], the minimal form of ${{L}}_{\text{NN}}^{(4)}$ is given in Reinert et al. [17], ${{L}}_{\pi NN}^{(1)}$ and ${{L}}_{\text{NNN}}^{(0)}$ are discussed in Epelbaum et al. [44] while ${{L}}_{\text{NNN}}^{(2)}$ was constructed in Girlanda et al. [45]. Notice further that the chiral symmetry breaking terms $\propto M_{\pi }^2$ are not shown explicitly in ${{L}}_{\text{NN}}$ and ${{L}}_{\text{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}}_{\pi 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}}_{\pi NN}$, which have not been worked out yet and would be needed to finalize the derivation of the 3NF at N⁴LO.

The long-range parts of the nuclear forces emerge from pion exchange diagrams and can be derived from ${{L}}_{\pi }$ and ${{L}}_{\pi 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⁴LO, namely c_{1, …, 4} from ${{L}}_{\pi N}^{(2)}$, d₁ + d₂, d_{3, 5, 18}, and d₁₄ − d₁₅ from ${{L}}_{\pi N}^{(3)}$ and e_{14, …, 18} from ${{L}}_{\pi N}^{(4)}$. Here, we made use of the fact that the contributions from the LECs l₃, 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 short-range part of the nuclear forces, and the corresponding LECs have to be determined from NN scattering and three- or more-nucleon observables.

In the single-nucleon sector, the effective Lagrangian ${{L}}_{\pi }+{{L}}_{\pi N}$ can be used to systematically compute the scattering amplitude in perturbation theory by applying the chiral expansion, a simultaneous expansion in particles' external three-momenta $p\equiv \left|\overrightarrow{p}\right|$ 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³. 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 ¹S₀ 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~{\Lambda }_{\text{b}}^2/M_{\pi }\gg {\Lambda }_{\text{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 A-body problem

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 scheme-dependent 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 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 three-nucleon potentials guarantees the validity of matching for the scattering amplitude. This can indeed be verified explicitly using the expressions for the 3NFs $V_{3\text{N}}^{2\pi ,1/m}$ from Equations (4.9) to (4.11) of Bernard et al. [55] and $V_{3\text{N}}^{2\pi -1\pi }$ from Equations (2.16) to (2.20) of Bernard et al. [54] and employing DR to evaluate loop integrals⁴.

Clearly, DR is impractical for a numerical solution of the A-body 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 non-renormalizable 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–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 many-body 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}(Λ), C_i(Λ), D_i(Λ), D(Λ), E(Λ), E_i(Λ), …, 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_i^{\text{r}}(\mu )+R_i(\mu )$. 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⁵, which leaves the implicit renormalization outlined above as the only available option. Notice that contrary to the renormalized LECs $l_i^{\text{r}}(\mu )$, $d_i^{\text{r}}(\mu )$, …, 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 residual Λ-dependence 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]⁶. 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].

3. 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₀ 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 (12) 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 (14) 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_{\text{eff}}^{\text{UT}}$ from the effective chiral Lagrangian in Equation (1) one needs to solve the non-linear decoupling Equation (13) 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 $\nu =-2+\sum _i{V}_i{\kappa }_i$, where V_i denotes the number of vertices of type i and κ_i is the inverse mass dimension of the corresponding coupling constant, ${\kappa }_i=d_i+\frac{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⁷. 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 (15) 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⁰. 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² (NLO) involve the leading two-pion exchange potential (TPEP) [31, 96, 97] and 7 short range interactions ∝C_i. At order Q³ (N²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⁴ (N³LO), the NN potential receives the contributions from the leading three-pion exchange [99–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⁵ (N⁴LO), the NN potential receives corrections to the three-pion 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⁴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 short-range operators. Work is still in progress to derive the remaining 3NF and 4NF at N⁴LO. We further emphasize that all calculations mentioned above are carried out using DR or equivalent schemes.

The effective potential $V_{\text{eff}}^{\text{UT}}$ 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 on-shell physics. One example is an energy-independent potential defined by

The proof that $V_{\text{eff}}^{\text{EI}}$ of Equation (17) 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_{\text{eff}}^{\text{EI}}$ is manifestly non-hermitian. 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–111]). This is carried out in perturbation theory by counting the nucleon mass via m ~ Λ_b

where the superscripts indicate the chiral order Qⁿ. 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_{2\pi }^{(3)}(\nu =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⁸, it can be transformed to the hermitian potential $v_{2\pi }^{(3)}(\nu =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.

4. Regularization

4.1. 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 momentum-space matrix elements of the contact potential are multiplied by a simple nonlocal Gaussian regulator

Here and in what follows, $p\equiv \left|\overrightarrow{p}\right|$ and $p^{\prime }\equiv \left|{\overrightarrow{p}}^{\prime }\right|$. 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–115])⁹.

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 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=\left|\overrightarrow{l}\right|$ and $\overrightarrow{l}$ denoting the three-momentum of the exchanged pion. The introduction of the Gaussian form factor in the pion propagators leads to properly regularized long-range potentials that are finite at short distances in coordinate space. In order to have a clean separation of the long-range 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 coordinate-space 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\equiv \left|\overrightarrow{q}\right|\equiv |{\overrightarrow{p}}^{\prime }-\overrightarrow{p}|$ and ${\overrightarrow{\sigma }}_i$ (τ_i) are the Pauli spin (isospin) matrices of the i-th nucleon. Here, the static pion propagator has been regularized according to Equation (23) and a likewise-regularized LO contact interaction has been added to the OPEP. Its coefficient C(M_π),

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(\overrightarrow{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 $\sqrt{M_{\pi }^2+{\lambda }^2}$

where ${\overrightarrow{l}}_1$ and ${\overrightarrow{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 (23) 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 $\rho (\mu )=\Im (V_{2\pi }(q)){|}_{q=0^{+}-i\mu }$ which are readily available up to N⁴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 (27) whose number depends on the chiral order of the contribution at hand. Introducing the pion propagator regulators in Equation (26), the regularized generic spectral integral of Equation (27) is replaced by

see Reinert et al. [17] for more details. The resulting potential $V_{2\text{N},\Lambda }^{2\pi }(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 (24) or Equation (28), 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 short-range terms polynomial in q². 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 (29) convergent and according to our convention, we have additionally fixed the subtraction coefficients $C_{C,1}^2(\mu )$ and $C_{C,2}^2(\mu )$ by the requirement that ${W_{C,\Lambda }^{(2)}(r)|}_{r=0}={\frac{d^2}{d{r}^2}{W}_{C,\Lambda }^{(2)}(r)|}_{r=0}=0$. [The first derivative of $W_{C,\Lambda }^{(2)}(r)$ vanishes at the origin regardless of the subtraction coefficients]. Figure 3 shows the ratio of the regularized and unregularized expressions in Equation (29) 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_{C,1}^2(\mu )=C_{C,2}^2(\mu )=0$ does not vanish at the origin¹⁰.

4.2. 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³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³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]

with ${\overrightarrow{q}}_i={\overrightarrow{p}}_i^{\prime }-{\overrightarrow{p}}_i$, ${\overrightarrow{k}}_i=1/2({\overrightarrow{p}}_i^{\prime }-{\overrightarrow{p}}_i)$, and ${\overrightarrow{p}}_i$, (${\overrightarrow{p}}_i^{\prime }$) the initial (final) momenta of the i-th nucleon. The complete expression for the relativistic corrections to the 3NF at N³LO can be found in Bernard et al. [55]. We now consider the first iteration of $V_{3\text{N}}^{2\pi ,1/m}$ with the static OPEP

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_{3\text{N}}^{2\pi -1\pi }$ 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 three-nucleon 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 $\propto {\overrightarrow{q}}_3·{\overrightarrow{\sigma }}_3{\overrightarrow{q}}_3·{\overrightarrow{\sigma }}_1$ is canceled by the D counter term in the second 3NF diagram at N²LO in Figure 2. To cancel the linearly divergent contribution $\propto {\overrightarrow{q}}_2·{\overrightarrow{\sigma }}_2$ one would, however, need to introduce a vertex in ${{L}}_{\pi NN}^{(1)}$ 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_{\pi }^2$, cannot appear in ${{L}}_{\pi NN}^{(1)}$. 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) on-shell 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_{3\text{N}}^{2\pi -1\pi }$ with a cutoff regularization for the iterative contributions in Equation (34), see Krebs [77] for another example with the NN axial vector current operator at N³LO. Indeed, recalculating the loop integral in $V_{3\text{N}}^{2\pi -1\pi }$ 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_{3\text{N}}^{2\pi -1\pi }$ with all pion propagators being regularized according to Equation (23). This would indeed solve the problem with the cancelation of linear divergencies at N³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³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.

5. 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⁽ⁱ⁾, $\Delta 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 $\delta X^{(k)}=\sum _{n>k}\Delta 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_{\pi }^{\text{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_{\pi }^{\text{eff}}=M_{\pi }$ 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ⁱ, 0 ≤ i ≤ k, are then estimated via

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–126] as well as few- and many-nucleon [8, 22, 127–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 model ${\bar{C}}_{0.5-10}^{650}$ 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 (37) in terms of dimensionless expansion coefficients c_i by introducing a (generally dimensionfull) scale X_ref via

where¹¹

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_{\text{ref}}=\left|X^{(0)}\right|$ 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 $\text{pr}(c_i|\bar{c})$ with a hyperparameter $\bar{c}$. 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 ${\Delta }_k\equiv \sum _{n=k+1}^{\infty }{c}_n{Q}^n\simeq \sum _{n=k+1}^{k+h}{c}_n{Q}^n$ to take a value Δ_k = Δ, given the knowledge of {c_{i ≤ k}}, is given by Melendez et al. [132]

where the set A is defined as A = {n ∈ ℕ₀|n ≤ k ∧ n ≠ 1 ∧ n ≠ m} and

The model ${\bar{C}}_{0.5-10}^{650}$ utilizes the Gaussian prior of “set C” from Melendez et al. [132],

for which the integrals in Equation (43) can be performed analytically [132], and uses the values of h = 10, ${\bar{c}}_{<}=0.5$ and ${\bar{c}}_{>}=10$. Following Epelbaum [134], the scales that control the expansion parameter are set to $M_{\pi }^{\text{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.

6. Selected Results

6.1. 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_{{\pi }^{\pm }}=139.57$ MeV and $M_{{\pi }^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 magnetic-moment 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¹². 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⁴LO⁺ version of the potential. As the proper description of these data requires a precise reproduction of F-waves, the N⁴LO⁺ potential adds the four leading F-wave contact interactions

formally appearing at N⁵LO and entering the ³F₂, ¹F₃, ³F₃, and ³F₄ partial waves, to the N⁴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⁴LO⁺¹³. 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⁴LO and N⁴LO⁺, while we reduced the energy to E_max = 25, 100, 125, and 200 MeV at the orders LO, NLO, N²LO, and N³LO, respectively. Notice that balancing the tradeoff between these two competing features can be handled using Bayesian methods (see e.g., [143]). From N³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 χ²/datum values at all considered orders for the cutoff Λ = 450 MeV¹⁴. 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 adjustable parameters at each order which also includes isospin-breaking LECs contributing to the ¹S₀ partial wave. It should be noted that no new contact interactions are added when going from NLO to N²LO and that the observed improvement of the χ²/datum values is entirely due to the N²LO contributions to the parameter-free TPEP. A similar situation occurs when going from N³LO to N⁴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³LO, a satisfactory description of the neutron-proton 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⁴LO potential improves on this, especially at intermediate and higher energies, it does not achieve a χ²/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⁴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 χ²/datum value of 27.88.

The situation is much improved once we switch to the N⁴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 χ²/datum ~ 1 description of the complete scattering database. Throughout the orders LO − N⁴LO the χ²/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⁴LO⁺, the reproduction of proton-proton data becomes very accurate while the slightly larger χ²/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⁴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⁴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⁴LO⁺ potential of Entem et al. [23] at the intermediate cutoff Λ = 500 MeV.

In general, there is good agreement between the shown N⁴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 ³P₀ phase shift where the N⁴LO⁺ prediction for the proton-proton 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 neutron-proton 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⁴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 ³P₀, ³P₂, and ³D₂ 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 ¹P₁ and ³P₁ 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 ¹G₄, ³H₄, and ϵ₄ phase shifts and mixing angle shown in Figure 5 since they do not involve any adjustable short-range parameters at N⁴LO⁺ but are solely determined by the long-range pion-exchange potential. Here, the local regulator of Equation (23) 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⁴LO⁺ potential as discussed above. The data are well described within the given truncation error for all considered orders, but the N⁴LO⁺ clearly allows for a proper quantitative description. Likewise, the reproduction of the spin observables in Figure 6 is excellent already at N³LO with a good convergence pattern. Notice however, that the error bands at lower orders for D (A) at the minimum (maximum) around ${\Theta }_{\text{CM}}=150^{\bullet }$ do not overlap with the ones for N^≥3LO 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³LO while lower-order corrections are small, such that the overall scale in Equation (42) 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²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]

where the LECs ${\tilde{C}}_i$, C_i, D_i, and E_i start to contribute at order Q⁰, Q², Q⁴, and Q⁶, 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 partial-wave 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 (47), we consequently should obtain values of unit magnitude. Figure 7 shows the absolute values of the LECs at N⁴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⁴ 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⁶ LECs E_i included at N⁴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 χ²/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 from ${\tilde{C}}_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 neutron-proton 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 ¹S₀ and ³S₁ − ³D₁ channels. Nevertheless, all LECs are well-constrained as can already be seen by looking at the errors in Figure 7. We can further look at the largest eigenvalue of the covariance matrix of the natural LECs as a measure of how well-determined the parameters are. Throughout the considered range of the cutoff Λ = 400 − 550 MeV, the largest eigenvalue of the covariance matrix does not exceed 0.1 and is ~ 0.08 for Λ = 450 MeV.

From the point of view of data fitting, another check concerns the statistical assumptions underlying a χ² fit. One usually assumes that the residuals $r_i=(O_i^{\text{exp}}-O_i^{\text{th}})/\Delta O_i$ follow a normal distribution ${N}(0,1)$ with zero mean and unit standard deviation. Here $O_i^{\text{exp}}$ and ΔO_i are the experimental value and its error of an observable and $O_i^{\text{th}}$ 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 χ² per degree of freedom. For the N⁴LO⁺ fit with Λ = 450 MeV we get χ² = 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 χ²/ν = 1.026 with ν = N_dat − N_par. If the residuals are indeed normal-distributed then χ²/ν should follow the χ²-distribution and yields ${\chi }^2/\nu =1\pm \sqrt{2/\nu }=1\pm 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 “tail-sensitive 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⁴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 quantile-quantile 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 cutoff-artifacts 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⁴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.

• Statistical uncertainties of πN LECs: In addition to the central values, the authors of Hoferichter et al. [50] also give the covariance matrix as determined from πN scattering data. Propagation of these uncertainties to NN observables is, however, less straightforward, because the values of the NN contact interactions depend on the values of πN LECs. We thus resort to some Monte Carlo sampling of the multivariate Gaussian probability distribution of the πN LECs given by their central values and their covariance matrix. For each of the sampled sets of LECs, we refit the NN contact LECs before calculating any observables. The uncertainty of a given observable can then be estimated in a standard way from the variance of the results calculated with different πN LEC 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⁴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⁴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-of-belief. 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_{\pi }^{\text{eff}}/{\Lambda }_{\text{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⁴LO⁺ and for Λ = 450 MeV the values of $A_s=0.8847_{(-3)}^{(+3)}(5)(0)(1)$ and $\eta =0.02553_{(-9)}^{(+11)}(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_{{\pi }^{\pm }}$ and $M_{{\pi }^0}$ and charge dependence of the short-range potential in the ¹S₀ partial wave. These are the dominant and well-understood isospin-breaking effects necessary to arrive at e.g. a correct description of the charge-dependence of the ¹S₀ scattering length. For the calculation of scattering observables in the two-nucleon system, the isospin-breaking due to long-range 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 isospin-breaking 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 determination). Last but not least, charge-dependence in the short-range potential entering P-waves should also be taken into account starting from N⁴LO [163].

6.2. Three-Nucleon Scattering

As discussed in the previous subsection, the N⁴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 high-precision 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⁴LO⁺ ones of [17]) tend to underestimate the ³H and ³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 (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–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²LO in the chiral expansion and are, therefore, suppressed by Q³ 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). For ³H and ⁴He, one can use the typical expectation values of the NN potential energy of ${\langle V_{\text{NN}}\rangle }_{3_{\text{ H}}}~50$ MeV and ${\langle V_{\text{NN}}\rangle }_{4_{\text{ He}}}~100$ MeV [83], along with the estimation of the expansion parameter $Q~M_{\pi }^{\text{eff}}/{\Lambda }_{\text{b}}$ with $M_{\pi }^{\text{eff}}=200$ MeV and Λ_b = 650 MeV, in order to estimate the expected 3NF contributions to the binding energy to be ${\langle V_{\text{3N}}\rangle }_{3_{\text{ H}}}~Q^3{\langle V_{\text{NN}}\rangle }_{3_{\text{ H}}}~1.5$ MeV and ${\langle V_{\text{3N}}\rangle }_{4_{\text{ He}}}~Q^3{\langle V_{\text{NN}}\rangle }_{4_{\text{ He}}}~3$ MeV. These qualitative estimations agree well with the actual underprediction of the ³H and ⁴He by the NN interactions alone which, using the AV18 [176], CD Bonn [146], N²LO [113], and Idaho N³LO [112] potentials as representative examples, amounts to 0.5…0.9 MeV and 2.1…4.1 MeV, respectively. 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% 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 (38) is dominated by the momentum scale $p=\sqrt{2/3m{E}_{\text{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-of-envelope 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²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 ${\overrightarrow{q}}_i={\overrightarrow{p}}_i\prime -{\overrightarrow{p}}_i$ with ${\overrightarrow{p}}_i\prime$ and ${\overrightarrow{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_{\pi }^2{\Lambda }_{\chi })$ and $E=c_E/(F_{\pi }^4{\Lambda }_{\chi })$ [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²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 ³H, ⁴He, and heavier nuclei, nucleon-deuteron doublet scattering length ²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 ³H binding energy of B_3H 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 ²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 ³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²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 ³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²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²LO have obviously rather low precision at this energy. In fact, the N²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 dashed-dotted 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⁴LO. At the energy of E_N = 135 MeV, the uncertainty bands at N⁴LO are expected to become 4-5 times more narrow as compared with the N²LO ones shown in Figure 11.

So where do we stand in terms of efforts to include 3NF corrections beyond N²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³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³LO and N⁴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⁴LO, is the determination of the LECs appearing in the 3NF at this order. While the N³LO contributions do not involve free parameters, the short-range part of the 3NF at N⁴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⁴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²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⁴LO can be tuned to reproduce the ³H binding energy, nd scattering lengths, cross section and polarization observables of pd scattering at 2 MeV center-of-mass 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₁ and E₇, respectively,

where ${\overrightarrow{K}}_i=({\overrightarrow{p}}_i\prime -{\overrightarrow{p}}_i)/2$. Parameterizing the dimension-full LECs E₁, E₇ in terms of the corresponding dimensionless parameters via $E_i=c_{E_i}/(F_{\pi }^4{\Lambda }_{\chi }^3)$ with Λ_χ = 700 MeV, we studied the impact of these N⁴LO terms on selected Nd scattering observables for the fixed values of the LECs of c_Ei = ±2. Based on the experience in the NN sector and with the N²LO 3NF, we expect the actual values of these LECs to lie well within this range. The expectation values of various contributions to the 3NF in the triton state indicate that the employed values c_E7 = ±2 may already overestimate the expected natural range of this LEC.

In order to compute the contributions of the c_Ei-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²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_Ei. The calculations have been performed using the N⁴LO⁺ SMS NN potential from Reinert et al. [17] in combination with the SMS N²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³LO, obtained by rescaling the N²LO Bayesian truncation uncertainty with the expansion parameter Q¹⁵, and visualize the expected impact of N⁴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²LO (see also [30, 45] for a similar conclusion). While the A_y is well-known to be particularly sensitive to spin-orbit types of 3NFs [191] such as the one proportional to c_E7, our results also show an unexpectedly strong sensitivity to the subleading central interaction of the c_E1-type.

At higher energies, the effects of the considered N⁴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_Ei-terms on Nd scattering observables is, in general, consistent with the estimated N³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.

6.3. 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²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 ⁴⁸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–195] and employing Similarity Renormalization Group (SRG) transformed interactions [196–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–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²LO and underbinding at N³LO and N⁴LO. Notice that the strongly repulsive nature of the N³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³LO. No large gap between the N²LO and N³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⁴LO⁺.

In Epelbaum et al. [8], a complete N²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²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²LO but appear to be significantly larger at NLO. We also calculated in that paper the excitation energies for selected states of A = 6 − 12 nuclei and the point-proton radius of ⁴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 ¹⁶O, which appears to be overbound. Notice that the deviation between the predicted and experimental values of the ¹⁶O binding energy is comparable to the 95% DoB Bayesian error at N²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.

7. 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³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, which—from 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⁴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⁴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²LO [84]. The LECs c_D and c_E, which enter the 3NF at this order, have been determined from the ³H binding energy and the very precise pd cross section data at $E_{\text{lab}}^N=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⁴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⁴LO.

• The novel semilocal nuclear forces, derived in the finite-cutoff 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⁴ 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 ³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³ (i.e., N²LO), this improvement results solely from taking into account the parameter-free subleading TPEP contributions. Notice that certain alternative power counting schemes suggest that some of the contact interactions that appear at order Q⁴ in the Weinberg scheme are enhanced and should yield contributions to observables larger than the order-Q³ TPEP (see e.g., Table 1 of [81]). The clear evidence of the chiral TPEP at orders Q³ and Q⁵ 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³LO and manifests itself in their dependence on arbitrary real phases ${\bar{\beta }}_8$, ${\bar{\beta }}_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 ¹S₀ and ³S₁-³D₁ channels proportional to the LECs $D_{1S0}^{\text{off}}$, $D_{3S1}^{\text{off}}$, and $D_{ϵ1}^{\text{off}}$ [17, 204, 205]. The SMS potentials of Reinert et al. [17] make use of the standard choice for ${\bar{\beta }}_{8,9}$ namely ${\bar{\beta }}_8=-{\bar{\beta }}_9=1/4$, which minimizes the amount of 1/m²-corrections to the OPEP, and employ $D_{1S0}^{\text{off}}=D_{3S1}^{\text{off}}=D_{ϵ1}^{\text{off}}=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 ${\bar{\beta }}_{8,9}$ and $D_i^{\text{off}}$, therefore, probes the impact of neglected higher-order terms and should lie within the estimated truncation errors. We have redone the fits at N⁴LO⁺ for Λ = 450 MeV using alternative choices of $D_i^{\text{off}}$ [17] and also developed a version of the potential with ${\bar{\beta }}_8={\bar{\beta }}_9=1/2$ [206]. The letter choice is motivated by the vanishing isoscalar exchange charge density operator at N³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²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³LO as detailed in section 4.2. This issue has not been paid attention to in the recent calculations involving the 3NFs [208–211] and exchange electroweak currents [110, 212, 213] at N³LO. As explained in section 4.2, using ad hoc regularization approaches at N³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⁴LO. While the N³LO contributions to the 3NF and 4NF do not involve unknown parameters, the N⁴LO corrections to the 3NF involve 10 LEC accompanying purely short-range operators [45] and one or more LECs entering the one-pion-exchange-contact 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_{\text{str}}=1.9731_{-0.0018}^{+0.0013}$ fm, by pushing the chiral expansion of the electromagnetic exchange charge density beyond N³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³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³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.

Funding

This work was supported by BMBF (Grant No. 05P18PCFP1) and by DFG through funds provided to the Sino-German CRC 110 110 Symmetries and the Emergence of Structure in QCD (Grant No. TRR110).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We are grateful to Ashot Gasparyan, Jambul Gegelia, and Ulf-G. Meißner for a careful reading of the manuscript and useful comments and to all members of the LENPIC Collaboration for sharing their insights into the topics addressed in this review article. We also thank Pieter Maris for providing us with a source file of Figure 8 from Epelbaum et al. [8].

Footnotes

References