Nucleon-Nucleon Scattering Up to N5LO in Chiral Effective Field Theory

During the past few decades a large effort has been made toward describing the NN interaction in the framework of chiral Effective Field Theory (EFT). The main idea is to exploit the symmetries of QCD to obtain an effective theory for low energy nuclear systems. In 2003, the first accurate charge-dependent NN potential in this scheme was developed and it has been applied to many ab-initio calculations, opening the possibility to study nuclear systems in a systematic and accurate way. It was shown that the fourth order (N3LO) was necessary and sufficient to describe the NN scattering data with a χ2/d.o.f on the order of so-called high precision potentials. However the systematics of chiral EFT also allow to relate two- and many-body interactions in a well-defined way. Since many-body forces make their first appearance at higher order, they are substantially smaller than their two-body counterparts, but may never-the-less be crucial for some processes. Thus, there are observables where they can have a big impact and, for example, there are indications that they solve the long standing Ay puzzle of N-d scattering. The last few years, have also seen substantial progress toward higher orders of chiral EFT which was motivated by the fact that only three-body forces of rather high order may solve some outstanding issues in microscopic nuclear structure and reactions. In this chapter we will review the latest contributions of the authors to development of chiral EFT based potentials up to N4LO as well as first calculations conducted for NN scattering at N5LO.


INTRODUCTION
The modern view of the NN interaction is given in the framework of Chiral Effective Field Theory (χEFT). The concept of an Effective Field Theory (EFT) is not a new one. The main idea is to identify the relevant degrees of freedom and symmetries for a certain system at a certain scale, and use this to find a Quantum Field Theory that is able to describe the system. However the traditional renormalization condition used to build theories like QCD is not required and a renormalization order by order is used instead. Nowadays, this approach is widely applied in different areas of physics.
In the case of strong interactions, we know that the fundamental theory is given by Quantum Chromodynamics (QCD). However for nuclear systems, the relevant degrees of freedom are not quarks and gluons, but nucleons and pions. Applying the EFT concept to nuclear systems allows to build theories for nucleons and pions that are consistent with the symmetries of the underlying theory. In the case of QCD, a very important property for low energy dynamics is that the original approximate chiral symmetry is broken spontaneously. This effect makes the pion come into play as the pseudo-Goldstone boson of the theory, which naturally explains the low mass of the pion as compared to other scales in nuclear systems. Chiral Perturbation Theory (ChPT) uses these ideas to determine observables making a perturbative expansion in the pion mass or some low energy external momenta. The Goldstone-boson character of the pion allows for this perturbative expansion, having always derivative couplings. ChPT was first applied to ππ systems [1] and πN systems [2] with quite some success. Chiral EFT is essentially based on ChPT, however in the case of the NN interaction this perturbative expansion is inadequate and non-perturbative resummations are needed. The complicate structure of the amplitudes makes it difficult to resum these contributions using the techniques of Unitarized ChPT that are applied in two-meson systems [3]. However first attempts to use similar techniques using the so called N/D method have been made [4].
The use of χEFT for the two-nucleon system was introduced by Weinberg in two seminal papers [5,6]. Weinberg realized that reducible diagrams violate the chiral expansion and, therefore, proposed to determine the potential using the rules of ChPT and then insert it into a Schrödinger-like equation to conduct the non-perturbative resummation.
Soon after, the first nuclear potentials were obtained by Ordoñez and van Kolck [7][8][9]. These position-space potentials were developed up to next-to-next-to-leading order (N 2 LO) and regularized by a cutoff function. Momentum-space potentials up to N 2 LO using dimensional regularization were derived by the Bochum group [10,11]. The simple and transparent momentumspace expressions obtained in this type of derivation [12] made chiral potentials more popular. However it was not until 2003 that χEFT reached high precision when the first chiral potential at N 3 LO was developed by Entem and Machleidt [13,14] that was able to describe the NN scattering data with a χ 2 /d.o.f similar to what the high-precision potentials of the 90's had achieved [15][16][17][18].
Since then, many applications of N 3 LO NN potentials together with chiral three-nucleon forces (3NFs) have been reported. These investigations include few-nucleon reactions [19][20][21][22], structure of light-and medium-mass nuclei [23][24][25][26][27] and infinite matter [28][29][30][31][32][33]. Although satisfactory predictions have been obtained in many cases, persistent problems continue to pose serious challenges, as the overbinding in medium mass nuclei [25] or the descriptions of charge and matter radii [34]. There is also the well-known A y puzzle of nucleon-deuteron scattering [35]. In this case recent calculations including contact 3NFs at N 4 LO have been shown to be able to solve the puzzle [36]. This suggests that one may have to proceed to the next higher order, namely, N 4 LO, for the two-nucleon force.
Thus, during the past few years, chiral potentials up to N 4 LO have been developed by the Idaho-Salamanca group [37] as well as the Bochum group [38].
In the whole chapter we will be referring to the so called -less EFT, where degrees of freedom have been integrated out. There are recent advances in the -full theory [39,40]. We refer the interested reader to contributions on this topic in the present monograph.
The chapter is organized as follows. In section 2 we review the most important aspects of χEFT for the two-nucleon system. In section 3 we apply the perturbative amplitude obtained to study peripheral NN scattering up to N 5 LO. In section 4 we review NN potentials up to N 4 LO. We conclude with a summary in section 5.

Power Counting
In order to build an EFT for the two nucleon system, the Lagrangians for the involved degrees of freedom have to be constructed. However, there is an infinite number of terms in the Lagrangian compatible with the allowed symmetries. For this reason, it is necessary to order all terms by what we call power counting. Following power counting, the terms in the Lagrangian are arranged by order. Moreover, the diagrams representing an amplitude calculated from the Lagrangian are also of a well defined order. Since higher orders include loop diagrams that diverge, the power counting also needs to be such that all the infinities generated at a certain order can be reabsorbed into redefinitions of the coupling constants of the Lagrangian at the same order. With these ideas in mind Weinberg, proposed the so called Weinberg power counting which is based on naive dimensional analysis.
Following naive dimensional analysis, a nucleon propagator counts as Q −1 , where Q stands for a low momentum or pion mass, a pion propagator as Q −2 , each derivative or pion mass insertion counts as Q and each four momentum integration as Q 4 . The power of a diagram is then given by the simple formula [5,6,14] where A is the number of nucleons involved, C the number of connected pieces, L the number of loops, and the sum runs over all vertexes i with i the index of the vertex given by with d i the number of derivatives or pion mass insertions (chiral dimension) and n i the number of nucleon legs. In this way the contribution of a diagram goes as (Q/ b ) ν with b the breakdown scale.
In the heavy-baryon formalism, an expansion in terms of Q/M N is performed, with M N denoting the nucleon mass. It is used for low energy nucleon systems and we will count these contributions as Q/M N ∼ (Q/ b ) 2 for reasons explained in Weinberg [5,6].
An important property of chiral symmetry is that the index of the vertexes is always zero or positive i ≥ 0. This fact implies that for a fixed number of nucleons with A ≥ 2 and considering diagrams with one connected piece, the power of a diagram is always bounded from below. This fact is crucial for the convergence of the chiral expansion.
A very important aspect of the EFT is that it relates two-body forces with many-body forces. We know that two-body forces are the main contribution to nuclear forces, however, many-body forces should exist. If we consider lowest order diagrams with L = 0 and i = 0, for an m-body force in an A-nucleon system, the number of separately connected pieces is C = A − m + 1, and so the power of the diagram is given by ν = 2m − 4. This means that two-body forces (m = 2) appear at ν = 0, three-body forces (m = 3) at ν = 2, four-body (m = 4) at ν = 4 and so on. So the power counting explains in a simple way the hierarchy of nuclear forces. In Figure 1 we summarize this hierarchy up to N 5 LO or sixth order of the chiral expansion.

The Lagrangian
We will limit ourselves to the -less version of χEFT, and so the relevant degrees of freedom are pions and nucleons. The effective Lagrangian, subdivided in terms of the number of nucleon legs, is given by where L π π stands for the Lagrangian that deals with pion dynamics, L π N the interaction between pions and a nucleon, and L NN contains four nucleon legs and no pion fields. The ellipsis stands for terms that involve two nucleons plus pions and three or more nucleons with or without pions, not relevant for the two nucleon sector. All the pieces in the Lagrangian are then organized in terms of the chiral dimension (number of derivatives/pion mass insertions) of increasing order L π π = L (2) π π + L (4) π π + . . . , where the superscript refers to the chiral dimension and the ellipsis refers to terms of higher dimensions. We use the heavybaryon formulation of the Lagrangians, the explicit expressions of which can be found in Machleidt and Entem [14] and Krebs et al. [41]. Notice that only in the NN case the chiral dimension is the same as the index i .

The Scattering Amplitude
Having the Lagrangian, we can now calculate the NN scattering amplitude. The NN amplitude has contributions from irreducible as well as reducible diagrams. The reducible diagrams are those that we can separate into two diagrams by cutting only nucleon lines. In covariant perturbation theory the separation is well defined, however when we apply a three-dimensional reduction of the Bethe-Salpeter equation it depends on the way this reduction is performed. See Machleidt and Entem [14] for a discussion on this point. We will come back to this when we define the potential. The amplitude for diagrams involving pions is organized in terms of the number of pions exchanged by the two nucleons Then each piece is divided in terms of the power counting described previously as where the superscript denotes the order ν. Besides these diagrams, contributions coming from Lagrangian L NN are also present. These contributions are contact-like contributions and take into account the unknown short-distance dynamics. They are again organized using the power counting where the superscript is the order ν. Due to symmetry requirements these contributions come only in even powers. Then the order by order contributions are given by where LO stands for leading order, NLO next-to-leading order, etc. For the presentation of amplitudes we will use the following decomposition where p ′ and p denote the final and initial nucleon momenta in the center-of-mass system (CMS), respectively. Moreover, q = p ′ − p is the momentum transfer, k = ( p ′ + p)/2 the average momentum, and S = ( σ 1 + σ 2 )/2 the total spin, with σ 1,2 and τ 1,2 the spin and isospin operators, of nucleon 1 and 2, respectively. For on-shell scattering, V α and W α (α = C, S, LS, T, σ L) can be expressed as functions of q = | q | and p = | p ′ | = | p |, only.
FIGURE 1 | Hierarchy of nuclear forces up to N 5 LO or sixth order of the chiral expansion. Only some representative diagrams are included. Small dots, large solid dots, solid squares, triangles, diamonds, and stars denote vertexes of index i = 0, 1, 2, 3, 4, and 6, respectively. Reprinted figure with permission from Entem et al. [37], copyright (2017) by the American Physical Society.

Pion-Exchange Contributions
We now specify the contributions coming from pion exchanges which provide the long-range interactions. Contributions at LO, NLO, and NNLO are diagrammatically given by the graphs in Figure 2.

Leading Order
The leading order (LO) is just the charge-independent one-pionexchange (OPE). The expression is given by where g A , f π , and m π denoted the axial-vector coupling constant, pion-decay constant, and the pion mass, respectively. There are corrections at higher orders that renormalize the coupling constant. They are taken into account by using g A /f π = g π N /M N , with g π N the πNN coupling constant. Numerical values are given in Table 1. Note that, on-shell, there are no relativistic corrections. Charge dependence is taken into account using with I the isospin of the two-nucleon system and m π 0 denotes the mass of the neutral pion and m π ± the one of the charged pion. The charge dependence is an NLO effect [14], but we include it already at leading order to make comparison with phase-shifts more meaningful.

Next-to-Leading Order
The NLO contributions appear at order ν = 2. Symmetry requirements make the contributions at ν = 1 vanish. In the past, the expressions for these diagrams as obtained in dimensional regularization were used [14]. Here, we apply the socalled spectral-function regularization (SFR) [42]. The potentials are obtained using dispersion relations from the imaginary part of the amplitude in the left-hand cut. However a cut-off˜ is used in the dispersion relation to constrain the potentials to the low-energy region where χEFT is applicable. The contribution is given by with which agrees with the dimensional regularization expressions [14] when replacing L(˜ ; q) by L(q). In fact,

Next-to-Next-to-Leading Order
Here the diagrams that contribute include a vertex with i = 1 which is represented by a large solid dot in Figure 2. The NNLO contribution is with As in the case of the NLO contribution, dimensional regularization is recovered when using Notice that, here, we demote the relativistic corrections of the NLO diagrams to N 3 LO, while in Machleidt and Entem [14] they were counted NNLO.

N 3 LO Contributions
At this order the first 3π exchange contributions appear. However it was shown in Kaiser [43,44] that they give negligible There are three types of contributions given by the three classes represented in Figure 3. The first one is the football diagram (a). The contribution is [45], The second class (b) corresponds to the 2π-exchange two-loop diagrams.
Here as well as for the N 4 LO expressions (see below), we state contributions in terms of their spectral functions, from which the momentum-space amplitudes V α (q) and W α (q) are obtained via the subtracted dispersion integrals: and similarly for W C,S,T . The thresholds are given by n = 2 for two-pion exchange and n = 3 for three-pion exchange. For˜ → ∞ the above dispersion integrals yield the finite parts of loop-functions as in dimensional regularization, while for finite˜ >> nm π we employ the method known as spectralfunction regularization (SFR). The purpose of the finite scale˜ is to constrain the imaginary parts to the low-momentum region where chiral effective field theory is applicable. The spectral functions for class (b) are given by [45,46] Im V C = 3g 4 A (2m 2 π − µ 2 ) πµ(4f π ) 6 (m 2 π − 2µ 2 ) Here and below all imaginary parts are evaluated at iµ, because that is where they are needed for the calculation of the SFR integrals.
Finally the relativistic corrections of the NLO diagrams corresponding to class (c) are given by [14]

N 4 LO Contributions
The 2π-exchange contributions at N 4 LO have three different classes of diagrams shown in Figure 4. The contributions of  class (a) and (b) are given in terms of spectral functions and Equation (34). The spectral functions for class (a) are obtained by integrating the product of the leading one-loop πN amplitude and the subleading chiral ππNN vertex proportional to c i over the Lorentz-invariant 2π-phase space. The result for the nonvanishing amplitudes is given by [46] with the dimensionless variable u = µ/m π > 2 and the logarithmic function Class (b) is obtained in the same way but multiplying the oneloop πN amplitude proportional to c i (see [41] for details) and the leading-order chiral πN amplitude. The result is [46] where the only two independent LEC'sē 14 The 3π-exchange contributions at order N 4 LO are shown in Figure 5. The spectral functions have been calculated first in Kaiser [47] where the classification scheme applied in Figure 5 was introduced. Class XI vanishes while class X and part of class XIV give negligible contributions. Thus, we include in our calculations only class XII and XIII, and the V S contribution of class XIV. In Kaiser [47], the spectral functions were presented in terms of integrals over the invariant mass of a pion pair. These integrals have been solved analytically in Entem et al. [46], and the spectral functions are given by Im W (XII) Im W (XIII)

Going Beyond N 4 LO
The next order is N 5 LO or sixth order. At this order, no complete calculation exists; however, the presumed dominant contributions have been evaluated in Entem et al. [48].
As before, we will state contributions in terms of their spectral functions, from which the momentum-space amplitudes V α (q) and W α (q) are obtained via subtracted dispersion integrals which, for N 5 LO read: and similarly for W C,S,T . The thresholds are given by n = 2 for two-pion exchange and n = 3 for three-pion exchange. The 2π-exchange at N 5 LO is given by the diagrams of Figure 6. There are three different classes. Class (a) is obtained from the subleading one loop πN amplitude folded with the subleading ππNN vertex proportional to c i . The results for the non-vanishing spectral functions are The classification scheme of Kaiser [47] is applied. The same notation as in Figure 2 is use.
Reprinted figure with permission from Entem et al. [46], copyright (2015) by the American Physical Society.
with the dimensionless variable u = µ/m π > 2 and the logarithmic function B(u) defined in Equation (47). We give the result in terms of the independent pion LEC'sē 14 andē 18 . Class (b) is obtained from the leading one-loop πN amplitude folded by itself. The result is Class (c) is obtained from the leading two-loop πN amplitude with the tree-level πN amplitude. The two-loop πN amplitude has not been evaluated and we omit this class of diagrams. The next contribution is the 1/M 2 N correction to the leading one-loop chiral 2π-exchange diagrams. They were given in Kaiser [49] and are shown in Figure 7. The explicit expressions are The next contribution is given by 3π-exchange contributions. There are several classes of diagrams as shown in Figure 8. The class (a) diagrams are proportional to c 2 i . We use the same notation as in Kaiser [47] and Entem et al. [46].
The contributions to ImW S and Im(µ 2 W T − W S ) are split into three pieces according to their dependence on the isoscalar/isovector low-energy constants c 1,3 and c 4 : The next contribution is given by class (b). Each diagram includes the one-loop πN amplitude. Not all the contributions could be treated; only those contributions that are independent of the pion-nucleon CMS energy in the loop or linearly dependent could be included. The contributions are in general small. The omitted contributions are typically an order of magnitude smaller. Class Xb: Class XIb: Class XIIb: setting w = 1 + u 2 − 2uω 1 . Class XIIIb: setting again w = 1 + u 2 − 2uω 1 . Class XIVb: where the auxiliary function G(w) is defined as Finally 4π-exchange diagrams occur for the first time at N 5 LO. These diagrams are three loop diagrams with only leading vertices. As mentioned before, three-pion exchanges with just leading order vertices turned out to be negligible. For that reason, we expect the leading four-pion exchanges to be even smaller, and we leave them out.

NN Contact Terms
Contact terms are given by the NN piece of the Lagrangian Equation (6). They start at order ν = 0 with non-derivatives terms given by [5] They contribute to S waves, only. The next order is ν = 2 (NNLO), which introduces seven new contact terms, given by [11] V (2) ct ( p ′ , p) = C 1 q 2 + C 2 k 2 + C 3 q 2 + C 4 k 2 σ 1 · σ 2 The next order is ν = 4 (N 3 LO) which has 15 contributions given by We note that, on shell, there are only 12 independent operators. The redundancy on-shell has been shown to generate large correlations. Reinert et al. [38] and Wesolowski et al. [50] claim that removal of the three (on-shell) redundant operators improves the fit. The partial wave decomposition of all these terms can be found in Machleidt and Entem [14]. Contact contributions are polynomials in external momenta and they only give contributions to partial waves with L ≤ ν/2.

PERIPHERAL NN SCATTERING
Peripheral NN scattering is of special interest since it is less sensitive to the short distance dynamics. A way to study it is to consider partial waves with high angular momentum, since the centrifugal barrier prevents sensitivity to short distance forces.
In the framework of EFT, the short distance physics is mimicked by the contact terms. In momentum space, they are given by polynomial terms in external momenta. This has the property that they don't give contributions to all partial waves, but only to angular momenta L ≤ ν 2 . This means that, for example at N 5 LO, there are only contributions up to F-waves.
One important aspect of peripheral waves is that the interaction is weaker and perturbative calculations can be performed, so avoiding all the problems posed by singular interactions in the Lippmann-Schwinger equation. For these reasons, it can be viewed as a clean probe of chiral dynamics in the NN sector.
The calculation is conducted by using the K matrix perturbatively as with V π ( p ′ , p) the χEFT amplitude where the iteration of OPE has been subtracted, and V 2π ,it ( p ′ , p) representing the once iterated OPE given by where P denotes the principal value integral and E p ′′ = M 2 N + p ′′ 2 . There is no unique way to subtract the iterative part of OPE. The prescription given by Equation (112) is slightly different from the one used in Kaiser et al. [12]. The difference between them is reabsorbed in a redefinition of the irreducible part. See Appendix C of Machleidt and Entem [14] for more details. Now the order by order calculation is conducted as follows. At LO only OPE is included in V π and no iteration is included. At NLO V π up to order ν = 2 is included and V 2π ,it is included. Higher orders (NNLO, N 3 LO, etc) include V π up to this order and the once iterated OPE. N 3 LO and higher orders should also include the twice iterated OPE contribution. However the difference between the once iterated OPE and the infinitely iterated OPE is very small and can not be identified on the scale of the figures. For this reason, we omit iterations of OPE beyond what is contained in V 2π ,it .

Fifth-Order (N 4 LO) Results
The contributions at NNLO [12] and N 3 LO [51] are in general too attractive, especially when the c i LEC's obtained from πN scattering are used.
We analyze now the contributions at N 4 LO. In Figure 9 we show results for selected F and G waves. Curve (1) gives the results for the N 3 LO calculation. Curve (2) adds the relativistic corrections (proportional to c i /M N ) of the NNLO terms. In curve (3), the 2π-exchange two-loop contributions of class (a) (Figure 4 and section 2.4.5) are added. Curve (4) adds the twoloop contribution of class (b). Finally curve (5) adds 3π-exchange contributions giving the final result at N 4 LO. In all calculations a SFR cutoff˜ = 1.5 GeV is used.
One can see that 3π-exchange contributions are significantly smaller than 2π-exchanges which can be interpreted as a convergence in regard to the number of pions exchanged. The 3π contribution is the sum of individual contributions that can be sizable but they add up to a small final result.
The c i /M N and two-loop contributions are mainly repulsive which helps to overcome the excess of attraction at N 3 LO. An exception is the 1 F 3 partial wave where the two-loop contribution of class (b) gives attraction, resulting in too much attraction for the whole N 4 LO contribution at higher energies.
For F and G waves (except 1 F 3 ) the final N 4 LO result is in very good agreement with the empirical phase-shifts. An interesting case is the 3 G 5 that is a problem at N 3 LO [51]; however, the final result at N 4 LO is in almost perfect agreement with the phase-shift analysis.
Here we have used˜ = 1.5 GeV. It is interesting to note that other potentials constructed from dispersion relations like the Stony Brook [52] and the Paris [53] potentials cut the dispersion integral at µ 2 = 50m 2 π which is equivalent to a SFR cut-off of ∼ 1 GeV. In Figures 10, 11 we show the impact of the SFR cutoff on the results at different orders. In general the variations for N 3 LO are large and always too attractive while at N 4 LO variations are smaller and close to the data. We also include lower orders to compare the relative size of the order-by-order contributions. One would expect a convergence pattern going from NNLO to N 3 LO and further to N 4 LO; however, this is not the case as seen in Figures 10, 11.
Concerning the LECs used, note that in the calculations of this subsection, the "KH" set of LECs shown in Table 2 was applied, while in the calculations of the next subsection the "GW" set is employed.

Going Beyond Fifth Order
As mentioned before there is no complete calculation at sixth order (N 5 LO). However a study of peripheral NN scattering with the expected dominant contributions was performed in Entem et al. [48]. We present here the results at this order.
For N 5 LO we consider G and higher waves, since they are not affected by contact terms at this order. In Figure 12, we show how individual groups of diagrams contribute to two G waves. Curve (1) represents the N 4 LO result. Curve (2) adds the N 5 LO 2π-exchange contributions of class (a) and curve (3) adds also class (b) (Figure 6 and Section 2.4.6). 3π-exchange (Figure 8) of class (a) are included in curve (4) and class (b) is contained in curve (5). The final result at N 5 LO is given by curve (6) which includes the 1/M 2 N corrections. In all cases a SFR cutoff˜ = 800 MeV is used.
The two-loop 2π-exchange class (a) (Figure 6) generates a strong repulsive central force, while the spin-spin and tensor forces provided by this class are negligible. The fact that this class produces a relatively large contribution is not unexpected, since it is proportional to c 2 i . The 2π-exchange contribution class (b) creates a moderately repulsive central force and a noticeable tensor force, as the impact on 3 G 5 demonstrates.  The c i ,d i , andē i are the LECs of the second, third, and fourth order π N Lagrangian given in Krebs et al. [41] and are in units of GeV −1 , GeV −2 , and GeV −3 , respectively. GW refers to the LECs obtained fitting to the George Washington University partial wave analysis from Arndt et al. [54], while KH refers to the Karlsruhe-Helsinki analysis from Koch [55].
The 3π-exchange class (a) (Figure 8) is negligible in 1 G 4 , but noticeable in 3 G 5 and, therefore, it should not be neglected. This contribution is proportional to c 2 i , which suggests a nonnegligible size but it is typically smaller than the corresponding 2π-exchange contribution class (a). The 3π-exchange class (b) contribution turns out to be negligible [see the difference between curve (4) and (5) in Figure 12]. This may not be unexpected since it is a three-loop contribution with only leading-order vertexes. Finally the relativistic 1/M 2 N corrections to the leading 2π-exchange have a small but non-negligible impact, particularly in 3 G 5 .
The predictions for G and H waves are shown in Figure 13, with shaded bands corresponding to a variation of the SFR cut-of over the range 700-900 MeV. The N 5 LO contribution shows a moderately repulsive effect, reducing further the excess attraction at N 3 LO. The N 5 LO result is, in general, substantially smaller than the N 4 LO one, indicating a signature of convergence. At N 5 LO, there is excellent agreement with the data.
Concerning the values for the LECs, let us note again that, in this subsection, the "GW" set of LECs shown in Table 2 was used, while in the calculations of the previous subsection the "KH" set was applied. Figure 13 includes only the three highest orders. However, a comparison between all orders is also of interest. Therefore, we show in Figure 14 the contributions to phase shifts through all six chiral orders from LO to N 5 LO. Note that the difference between the LO prediction (one-pion-exchange) and the data (filled and open circles) is to be provided by two-and three-pion exchanges, i.e., the intermediate-range part of the nuclear force. How well that is accomplished is a crucial test for any theory of nuclear forces. NLO produces only a small contribution, but NNLO (denoted by N2LO in the figure) creates substantial intermediate-range attraction (most clearly seen in 1 G 4 , 3 G 5 , and 3 H 6 ). In fact, NNLO is the largest contribution among all orders. This is due to the one-loop 2π-exchange (2PE) triangle diagram which involves one ππNN-contact vertex proportional to c i . This vertex represents correlated 2PE as well as intermediate (1232)isobar excitation. It is well-known from the traditional meson theory of nuclear forces that these two features are crucial for a realistic and quantitative 2PE model. Consequently, the oneloop 2π-exchange at NNLO is attractive and assumes a realistic size describing the intermediate-range attraction of the nuclear force about right. At N 3 LO, more one-loop 2PE is added by the bubble diagram with two c i -vertices, a contribution that seemingly is overestimating the attraction. This attractive surplus is then compensated by the prevailingly repulsive two-loop 2πand 3π-exchanges that occur at N 4 LO and N 5 LO.
In this context, it is worth to note that also in conventional meson theory the one-loop models for the 2PE contribution always show some excess of attraction. In conventional meson theory, the surplus attraction is reduced by heavy-meson exchange (ρ-and ω-exchange) which, however, has no place in chiral effective field theory (as a finite-range contribution). Instead, in the latter approach, two-loop 2π-and 3π-exchanges provide the corrective action.

NN POTENTIALS UP TO N 4 LO
The starting point of all ab-initio calculations of nuclear systems is the NN potential. For that reason, it is necessary to define a potential.
We define the NN potential as the sum of the irreducible NN diagrams discussed in previous sections, which are calculated perturbatively. However, in reality, the NN system is characterized by the presence of a shallow bound state (the deuteron) and large (S-wave) scattering lengths that cannot be obtained perturbatively. Therefore, the potential has to be applied in a scattering equation to obtain the NN amplitude. Since our approach is in principal covariant (with relativity taken into account perturbatively), a proper equation would be the Bethe-Salpeter equation. However, it is more convenient, to use one of the three-dimensional reductions of that equation. We use the Blankenbeclar-Sugar (BbS) version of the equation [56] which reads where V is the potential and E p ′′ = M 2 N + p ′′2 . Since this is a relativistic equation, it includes relativistic kinematical corrections to all orders.
If we now definê The c i ,d i , andē i are the LECs of the second, third, and fourth order π N Lagrangian Krebs et al. [41] and are in units of GeV −1 , GeV −2 , and GeV −3 , respectively. The uncertainties in the last digits are given in parentheses after the values.
the BbS equation becomeŝ which is the Lippmann-Schwinger equation andV can be used like a non-relativistic potential. All the technical details to solve the Lippmann-Schwinger equation, including the case where the Coulomb interaction is included, can be found in Machleidt [18].
The amplitude V and the potentialV are built order-byorder following the Equations (12-16) with two exceptions. We add to V N 3 LO the 1/M N corrections of the NNLO 2π-exchange proportional to c i . This c i /M N correction is formally an N 4 LO contribution, however, in Entem et al. [46] it was shown that the football diagram proportional to c 2 i at N 3 LO was unrealistically attractive, while the c i /M N correction is large and repulsive. Therefore, it makes sense to group these diagrams together to arrive at a more realistic intermediate-range attraction at N 3 LO.
The other exception is to include, at N 4 LO, the four F-wave contacts that formally appear at N 5 LO, cf. Equation (17). This ensures an optimal fit of the NN data for the potential of the highest order to be constructed.

Regularization
The potentialV obtained previously is in most cases singular. Singular potentials are those that diverges in momentum space when the momentum goes to infinity, being more singular than 1/r 2 in coordinate space. For this reason they cannot be included in a Lippmann-Schwinger equation without further manipulation. The practical way to solve this problem is to cut the potential at a certain scale by multiplying with a regulator function f (p ′ , p)V where the function f (p ′ , p) can be taken to be This regularization allows to obtain finite results, however renormalization requires to have regularization independent results. The implicit assumption in Weinberg's proposal [5,6] was that the same contact interactions that renormalize loop diagrams would also renormalize the iterative loops of the (infinite) resummation in the Lippmann-Schwinger equation. This is not necessarily true and has given rise to a comprehensive discussion about non-perturbative renormalization. This is one of the key issues where the EFT community is divided, mainly, in two different points of view, one with the cut-off scale below the hard-scale of the EFT, and the other with a value above (let's say, infinity). This topic has been discussed by many authors [4,[57][58][59][60][61][62][63][64][65][66][67][68][69][70][71][72][73][74][75][76], and we refer the interested reader to contributions about this topic in the monograph. However, using cutoffs in the order of 450 − 550 MeV (first point of view) has been shown to give mild regularization dependence and to be phenomenologically successful at N 3 LO [77], although renormalization is not so clear. The parameter n is usually chosen in such a way that the corrections induced by the regulator are of an order that is higher than the given order. We choose n = 2 for 3PE and 2PE and n = 4 for OPE (except in LO and NLO, where we use n = 2 for OPE). For contacts of order ν, we choose 2n > ν.

Charge Dependence
In order to fit the np and pp databases, charge dependence has to be included. All orders include the charge dependence due to pion mass splitting in the one-pion exchange as was already discussed. Charge dependence is most important in the 1 S 0 partial wave at low energies, particularly in the scattering lengths. The charge dependence from OPE cannot explain it all. The remainder is accounted for by treating the 1 S 0 LO contact term parameterC1 S 0 ≡ 4π(C s −3C T ) in a charge-dependent way.
So, we distinguish betweenC . For pp at any order, the relativistic Coulomb interaction is included [78,79]. Finally at N 3 LO and N 4 LO, we take into account irreducible πγ exchange [80], which affects only the np potential. Also, the charge-dependent effects from n-p mass splitting are taking into account by using the correct values for the nucleon masses.
For a detailed discussion of possible sources for charge dependence of the NN interaction, see Machleidt and Entem [14].
The pion exchange contribution, V π , is fixed by the πN LECs for which we use the values from the very accurate analysis by Hoferichter et al. [86], Table 3. However, the short-range part given by V ct has to be determined from NN scattering. This was done by fitting the NN potentials to the NN database. The database includes all NN data below 350 MeV laboratory energy published in refereed physics journals between January 1955 and December 2016 that are not discarded when applying the Nijmegen rejection criteria [79]. There are alternative criteria   [87] which have been applied, e.g., in the Granada database [88], however we continue to use the Nijmegen criteria to be consistent with the pre-2000 part of our database.
The database finally consists of 3072 pp scattering data and 3569 np data. The 2013 Granada NN database [88] consists of 2996 pp and 3717 np data. The larger number of pp data in our base is mainly due to the inclusion of 140 pp data from The EDDA Collaboration [89] which are left out in the Granada base. On the other hand, the Granada base contains 148 more np data, which is a consequence of the modified rejection criteria applied by the Granada group, which allows for the survival of a few more np data.
In the fitting procedure, only data below 290 MeV were taken into account. One starts with the pp potential, since the pp data are more accurate than the np data. First, a fit to the pp phaseshifts is made, and then a rough minimization of the χ 2 is performed by using the Nijmegen error matrix [90]. In the end, the potential is fitted directly to the scattering data. For this the SAID software package [91] that includes all electromagnetic contributions necessary for the calculation of NN observables at low energy is used.
Then the I = 1 np potential is fixed by starting from the pp potential and applying charge dependence. For the 1 S 0 part of the np potential, theC (np) 1 S 0 LEC is adjusted to the np scattering length. The I = 0 part is then fitted in a similar way as the I = 1 part. After the I = 0 fit, some small variations of the I = 1 parameters were allowed to obtain a minimal over-all χ 2 .
The nn potential is obtained from the pp one by leaving out Coulomb, replacing the proton mass by the neutron mass, and fitting theC (nn) 1 S 0 LEC to the 1 S 0 nn scattering length.
The above procedure is basically the same as used in the construction of the so called high-precision potentials of the 1990s [15,16,18], which all have χ 2 /datum ≈ 1. This differs from the procedure applied in the recent construction of the NNLO sat potential [83] where NN data up to 35 MeV and the ground-state energies and radii of nuclei up to 16 O are taken into account to fix simultaneously the two-and three-nucleon forces. Our procedure also differs from the construction of some recent chiral NN potentials by the Bochum group [81,82], where only phase-shifts are fitted. However, in their most recent potential constructions, the Bochum group [38] does apply a procedure where the fitted potentials are directly confronted with the NN data.

Results for NN Scattering
The χ 2 /datum for the reproduction of the NN data is given in Table 4. For the close to 5000 pp plus np data below 290 MeV (pion-production threshold), the χ 2 /datum is 51.4 at NLO and 6.3 at NNLO, which is of special relevance since the number of NN contact terms is the same for both orders. The improvement comes entirely from a better description of the 2PE at NNLO. At N 3 LO, the χ 2 /datum further improves to 1.63. It, finally, reaches 1.15 at N 4 LO, in acordance with high precision potentials, showing a great convergence pattern.
np phase shifts are displayed in Figure 15, which reflect the same features as the χ 2 , namely, an excellent convergence when going from NNLO to N 3 LO and, finally, to N 4 LO. From Entem et al. [37]. However, at LO and NLO there are large discrepancies between the predictions and the empirical phase shifts as to be expected from the corresponding χ 2 values. This fact renders applications of the LO and NLO nuclear forces useless for any realistic calculation (but they could be used to demonstrate truncation errors).
It is important to be aware of the regulator dependence of the NN phase shifts and scattering observables. For this reason, potentials with cutoffs = 450, 500, and 550 MeV were constructed. We show in Figure 16 the phase shifts at NNLO (green curves, left panel) and N 4 LO (purple curves, right panel) for potentials with varying cutoffs. As   [18] for references; the empirical value for rstr is from Jentschura et al. [92]. For some of the notation, see Table 5, where also empirical information on the deuteron and triton can be found.
expected, the cutoff dependence diminishes with increasing order, being very small at N 4 LO. The cutoff window we selected is motivated by the fact that for values ≤ 450 MeV cutoff artifacts start to appear above 200 MeV as seen in the 1 D 2 and 3 D 2 partial waves. The upper limit is given by the fact that the breakdown scale occurs around b ∼ 600 MeV [82].

Deuteron and Triton
The deuteron binding energy is fitted at all orders to the empirical value of 2.224575 MeV using the nonderivative contact term in the 3 S 1 partial wave. Different observables of the deuteron and triton are given at all orders in Table 5. Notice that only the deuteron binding energy is fitted while all other observables are predictions. It is interesting to notice that already at NNLO all properties are close to the empirical values and vary little when going to higher orders, as one would expect, since they are low energy observables.
The triton binding energy is also given. A 34-channel charge dependent Faddeev calculation using only two-nucleon forces is used. The results show a smooth and steady convergence order by order toward a value around 8.1 MeV, giving some space to three-nucleon forces. The low deuteron D-state probabilities and the high triton binding energy predictions are due to the softness of the potentials.
In Table 6, we demonstrate, for order NNLO and N 4 LO, the cutoff dependence of the χ 2 /datum, the deuteron properties, and the triton binding energy. One observes a mild regulator dependence for most quantities. The exception is the deuteron D-state probability which, however, is not an observable. Linked to this (via the strength of tensor force) is the triton binding energy. This is due to the off-shell behavior of the two-nucleon force. This can be compensated by corresponding changes in the three-nucleon force.

SUMMARY
The past 25 years have seen great progress in our understanding of nuclear forces in terms of low-energy QCD. Key to this development was the realization that low-energy QCD is equivalent to an effective field theory which allows for a perturbative expansion that has become known as chiral perturbation theory. In this framework, two-and many-body forces emerge together and the empirical fact that nuclear manybody forces are substantially weaker than the two-nucleon force is explained naturally.
The main focus of this review, was on the two-nucleon force. We presented the order-by-order development from LO (∼ Q 0 ) to N 5 LO (∼ Q 6 ). Using low-energy constants (LECs) determined from πN scattering, our predictions for peripheral partial waves are parameter-free, except for the spectral function cutoff that regularizes the dispersion integrals which determine the NN amplitudes. This spectral-function regularization ensures that the calculated contributions are restricted to the longand intermediate range, where chiral effective field theory is applicable. Specifically, we have calculated perturbative NN scattering in peripheral partial-waves, which is dominated by one-, two-, and three-pion exchanges ruled by chiral symmetry. The order-by-order convergence is slow, but is ultimately achieved at N 5 LO, where predictions are in perfect agreement with empirical phase shifts.
Besides this, we have also discussed the construction of complete (i.e., including the lower partial waves) chiral NN potentials through all orders up to N 4 LO. The construction may be perceived as consistent, because the same power counting scheme as well as the same cutoff procedures are applied in all orders. The potential of the highest order (N 4 LO) reproduces the NN data below pion-production threshold with a χ 2 /datum of 1.15. This is among the highest precisions ever accomplished with any chiral NN potential to date. The NN potentials presented may serve as a solid basis for systematic ab initio calculations of nuclear structure and reactions that allow for a comprehensive error analysis. In particular, the order by order development of the potentials will make possible a reliable determination of the truncation error at each order.
In summary, this review presents the most comprehensive investigation of the implications of chiral symmetry for the NN system. The results provide the ultimate confirmation that chiral EFT is an adequate theory for nuclear forces.