Local nucleon-nucleon and three-nucleon interactions within chiral effective field theory

To obtain an understanding of the structure and reactions of nuclear systems from first principles has been a long-standing goal of nuclear physics. In this respect, few- and many-body systems provide a unique laboratory for studying nuclear interactions. During the past decades, the development of accurate representations of the nuclear force has undergone substantial progress. Particular emphasis has been devoted to chiral effective field theory (EFT), a low-energy effective representation of quantum chromodynamics (QCD). Within chiral EFT, many studies have been carried out dealing with the construction of both the nucleon-nucleon ($N\!N$) and three-nucleon ($3N$) interactions. The aim of the present article is to give a detailed overview of the chiral interaction models that are local in configuration space, and show recent results for nuclear systems obtained by employing these local chiral forces.


INTRODUCTION
The last few decades have marked the emergence of the basic model of nuclear theory in which nuclear systems -particularly atomic nuclei and infinite nucleonic matter -can be described as a collection of point-like particles, the nucleons, interacting with each other in terms of two-and many-body effective interactions, and with external electroweak probes via effective current operators. This approach, in conjunction with a computational method of choice to solve the many-body Schrödinger equation, can then be used to study the structure and dynamics of nuclear systems in a fully microscopic way, which is commonly referred to as ab-initio calculations. Examples of such calculations are based on the no-core shell model (NCSM) [1,2], the coupled cluster (CC) [3,4] or hyperspherical harmonics (HH) [5] expansions, similarity renormalization group (SRG) approaches [6,7], self-consistent Green's function techniques [8,9], quantum Monte Carlo (QMC) methods [10], and nuclear lattice effective field theory (NLEFT) [11]. Although significant progress has been made in recent years, these ab-initio techniques remain challenging and their domain of applicability is, at present, limited to provide quantitative description of light and medium-mass nuclei [1,4,7,9,8,10,12] and their reactions [13,14,15,16]. A special but related challenge is the development of microscopic models that include continuum couplings which are mandatory to describe, for instance, weakly bound nuclear systems [17,18].
One might argue that nucleons are not the fundamental building blocks of the nuclear systems at hand, and that one should instead start from Quantum Chromodynamics (QCD). QCD provides the theoretical framework to describe strong interactions which governs the dynamics and properties of quarks and gluons. However, while strong interactions are weak and perturbative at high energies, i.e., short distances (asymptotic freedom), quarks are strongly interacting at low energies or long distances, of relevance for nuclear physics, which makes a nonperturbative treatment necessary. In addition, at these energies quarks are confined into colorless objects called hadrons (baryons, consisting of three quarks, e.g., the nucleon, and mesons consisting of a quark and an anti-quark, e.g. the pion). Hence, while QCD is responsible for the complex inter-nucleon forces in nuclear systems, which can be thought of as residual interactions among quarks, a description in terms of nucleon degrees of freedom is particularly valid at sufficiently low energies.
How the interactions among nucleons emerge from the fundamental theory, QCD, has kept nuclear theorists occupied for many decades. Since QCD is nonperturbative at low energies of interest in nuclear systems, one may try to solve QCD with brute computing power on a discretized Euclidean space-time lattice (known as lattice QCD) However, in spite of many advances [19,20,21,22], lattice QCD calculations are still limited to small nucleon numbers and/or large pion masses, and hence, at the present time, can only be used to address a limited set of representative key-issues.
As a consequence, most theoretical studies of nuclear systems have to resort to using the basic model of nuclear theory, i.e., assuming pointlike nucleons to be the relevant degrees of freedom instead of quarks. In this review, we will briefly introduce this basic model and discuss the current state-of-the-art for nuclear interactions, chiral effective field theory (EFT). We will then focus on a particular subclass of chiral EFT interactions, local chiral EFT interactions, intended for the use in QMC methods.
The review is structured as follows. In Section 2, we discuss the general features of nuclear interactions starting with the phenomenological ones and moving to those obtained in chiral EFT. In Section 3, we provide many details about the theoretical derivation of local interactions in both delta-full and delta-less chiral EFT, i.e., when explicitly including the delta resonance or not. In Section 4, we briefly discuss finite cutoff and regulator artifacts that can appear in calculations with local interactions. Finally, in Section 5, we report selected results for light and medium-mass nuclei and the equation of state of pure neutron matter using QMC methods.

NUCLEAR HAMILTONIANS
The basic model of nuclear theory assumes that a nuclear system can be described by a non-relativistic Hamiltonian that contains interactions among nucleons, i.e., protons and neutrons. The individual nucleons mostly interact via two-body (NN ) interactions. However, nucleons can also interact via three-body (3N ) and higher many-body interactions. The way these many-body interactions appear is twofold. First, nucleons are compound particles and, hence, treating them as point-like particles induces effective manybody interactions even if only two-quark interactions were to be considered. This is similar to describing tides on Earth, where the three-body system given by Earth, Moon, and Sun is relevant, even though gravity is only a two-body force. Second, since quarks themselves can have multi-quark interactions, this immediately leads to the appearance of 'true' 3N forces among nucleons, where, for example, single quarks in each of the 3 nucleons interact with each other.
The resulting Hamiltonian can then be written as a sum of the non-relativistic one-body kinetic energy (T i ), NN interactions between particle i and j (V ij ), 3N interactions between particle i, j, and k (V ijk ), and additional many-body interactions, and provides a good approximation for interacting nucleons in a given nuclear system: There are indications that four-body interactions may contribute at the level of only ∼ 100 keV in 4 He [23] or pure neutron matter [24], and therefore are negligible compared to NN and 3N interactions. Hence, current formulations of the basic model do not typically include them (see, for example, Ref. [10]).
In order to derive two-and three-body nuclear forces, one has to take into account some general considerations, specify the theoretical framework in which such interactions are formulated, and the experimental inputs necessary to determine possible unknown parameters of the theory.

General considerations for nuclear interactions
To accurately describe nuclear systems that are governed by QCD, nuclear interactions need to obey all the relevant symmetries of QCD. Hence, nuclear potentials need to have the following properties (we will focus on NN forces here, but the statements remain true for all parts of the interaction): • V is hermitian, because the Hamiltonian is hermitian, • V is symmetric under the permutation of identical particles, i.e., V ij = V ji , • V is translationally and rotationally invariant, • V is Lorentz invariant (for nonrelativistic interactions this reduces to Galilean invariance), • V is invariant under parity transformations and time reversal, • V has to conserve baryon and lepton number, • V has to be approximately isospin symmetric and charge independent, • and V has to include the properties of spontaneously and explicitly broken chiral symmetry. Chiral symmetry is a symmetry of the QCD Lagrangian with massless quarks under independent rotations of left-and right-handed quarks. Considering only u and d quarks, this symmetry can be written as SU(2) L ×SU(2) R . This expression contains two symmetries: the first (vector) one represents isospin symmetry, i.e., symmetry under the exchange of u and d quarks, and the second (axial) one is the socalled chiral symmetry. These two symmetries imply degenerate fermions under isospin and spin-parity transformations. While isospin symmetry is approximately fulfilled in nature, i.e., the neutron and proton have similar masses, nucleons with spin 1/2 + and 1/2 − have very different masses (940 MeV vs 1535 MeV). This signals that chiral symmetry is broken in nature.
In fact, chiral symmetry is broken twofold. First, it is broken spontaneously, leading to the formation of Goldstone bosons, that can be identified with the pions. Second, chiral symmetry is also explicitly broken by the finite quark masses, which leads to the pion being pseudo-Goldstone bosons with finite but small mass. In contrast, isospin symmetry remains a good symmetry, because the ratio (m d − m u )/Λ QCD is very small, where m u 2.4 MeV and m d 4.8 MeV.
These symmetries only allow certain operator structures for nuclear interactions. Galilean invariance, for instance, implies that nuclear interactions depend only on relative momenta between two nucleons, p = p i − p j , while symmetry under parity transformations implies that nuclear interactions cannot be linear in p, and charge independence requires that the nuclear interactions can be written as and so on. In addition, the spin dependencies are included through operators like 1, σ i · σ j , spin-orbit interactions given by L · S with L = r × p, where r = r i − r j , or tensor interactions with the tensor operator S ij (r) = σ i ·r σ j ·r − σ i · σ j . As a consequence, interactions typically have a spin-isospin operator structure given by where the individual operators carry momentum-dependent functions consistent with all required symmetries.

Phenomenological interactions
Historically, N N interactions were derived using phenomenological insight. They were characterized by a long-range component characterizing the interaction for inter-nucleon separations r 1/m π , due to one-pion exchange (OPE) [25], and intermediate-and short-range components describing the interactions at 1 fm r 2 fm and r 1 fm, respectively. The intermediate-and short-range components were included to simulate intermediate-range attraction as well as short-range repulsion.
Up to the mid-1990's, nuclear interactions were based almost exclusively on meson-exchange phenomenology. Interactions of the mid-1990's [26,27,28] were constrained by fitting nucleon-nucleon (NN ) elastic scattering data up to laboratory energies of 350 MeV, with χ 2 /datum 1 relative to the database available at the time [29]. Two well-known and still widely used examples in this class are the Argonne v 18 (AV18) [27] and CD-Bonn [28] interactions. These are so-called phenomenological interactions.
Already in the 1980's, accurate three-body calculations showed that contemporary NN interactions alone did not provide sufficient binding to reproduce experimental numbers for nuclei with nucleon number A = 3, 3 H and 3 He [30]. This realization was later on extended to the spectra (ground and low-lying excited states) of light p-shell nuclei, for instance, in calculations based on quantum Monte Carlo (QMC) methods [31] and in no-core shell-model (NCSM) studies [32]. Consequently, the basic model with only NN interactions fit to scattering data, without the inclusion of a three-nucleon (3N ) interaction, was found to be unsatisfactory. However, because of the composite nature of the nucleon and, in particular, the dominant role of the ∆ resonance, a spin-3/2, isospin-3/2 nucleon resonance, in pion-nucleon scattering, many-body interactions arise quite naturally in meson-exchange phenomenology.
For example, the Illinois 3N interaction [33] consists of a dominant two-pion exchange (TPE) -the Fujita-Miyazawa interaction [34] -and smaller multi-pion exchange components resulting from the excitation of intermediate ∆'s. The most recent version, Illinois-7 (IL7) [35], also contains phenomenological isospindependent central terms. The parameters characterizing this 3N potential have been determined by fitting the low-lying spectra of nuclei in the mass range A = 3-10. The resulting AV18+IL7 Hamiltonian, generally utilized with QMC methods, then leads to predictions of 100 ground-and excited-state energies up to A = 12, including the 12 C ground-and Hoyle-state energies, in good agreement with the corresponding experimental values [10]. However, when used to compute the neutron-star equation of state, these interactions do not provide sufficient repulsion to guarantee the stability of the observed neutron stars with masses larger than two solar masses against gravitational collapse [36]. Thus, in the context of the phenomenological nuclear interactions, we do not have a Hamiltonian that can predict the properties of all nuclear systems, from NN scattering to dense nuclear and neutron matter.
Furthermore, high-precision phenomenological potentials suffer from several limitations, most notably the missing connection with the low-energy QCD and, hence, the absence of a guiding principle for the construction of interactions. As a consequence, phenomenological interactions do not provide rigorous schemes to consistently derive two-and three-body forces and compatible electroweak currents. In addition, there is no clear way to properly assess the theoretical uncertainty associated with the nuclear potentials and currents.

Chiral effective field theory
These drawbacks were addressed when a new phase in the evolution of the basic model began in the early 1990's with the emergence of chiral effective field theory (EFT) [37,38,39].
Chiral EFT is a low-energy effective theory of QCD based on the choice of baryons as effective degrees of freedom: in chiral EFT one chooses pions and nucleons. At typical momenta in nuclei, p ∼ m π ∼ O(100 MeV), this choice is accurate, because shorter-range structures, e.g., the quark substructure, or heavier meson exchanges, e.g., exchanges of the ρ-meson, are not resolved, and can be absorbed in short-range nucleon contact interactions. This separation of scales between typical momenta p and scales of the same order, i.e., the pion mass m π ∼ 140 MeV, and larger scales, e.g., the mass of the ρ, m ρ ∼ 770 MeV, can then be used to systematically derive an effective and most general scheme accommodating all possible interactions among the relevant degrees of freedom consistent with the symmetries of QCD. In some modern approaches, the choice of degrees of freedom also includes the ∆ isobar (delta-full chiral EFT), because the ∆-nucleon mass splitting is only 300 MeV∼ 2m π .
The starting point in chiral EFT is the most general Lagrangian in terms of the chosen degrees of freedom, which contains all allowed interaction mechanisms in accordance with the considerations in Sec. 2.1. As a consequence, this Lagrangian contains an infinite number of terms and needs to be truncated using a given power-counting scheme. Most chiral interactions used in nuclear structure calculations are based on Weinberg power counting, which itself is based on naive dimensional analysis of interaction contributions. Within Weinberg power counting, the interactions are expanded in powers of the typical momentum p over the breakdown scale Λ b , Q = p/Λ b , where the breakdown scale denotes momenta at which the short distance structure becomes important and cannot be neglected and absorbed into contact interactions anymore (see Refs. [40,41,42,43] for recent review articles). It is worthwhile mentioning that alternative power-counting schemes have been also suggested as in Refs. [44,45,46,47,48,49].
This expansion defines an order by order scheme, defined by the power ν of the expansion scale Q in each interaction contribution: leading order (LO) for ν = 0, next-to-leading order (NLO) for ν = 2, next-to-next-to-leading order (N 2 LO) for ν = 3 and so on. Similarly as for nuclear interactions, such a scheme can also be developed for electroweak currents. Therefore, chiral EFT provides a rigorous scheme to systematically construct many-body forces and consistent electroweak currents, and tools to estimate their uncertainties [50,51,52,53,54,55]. From this perspective, it can be justifiably argued that chiral EFT has put the basic model on a more fundamental basis, by providing a link between QCD with all its symmetries, and the strong and electroweak interactions in nuclei. Figure 1 shows the state of the art of chiral contributions to the NN and 3N interactions in the delta-less and delta-full chiral EFT. Higher many-body forces, such as four-nucleon (4N ) or five-nucleon (5N ) interactions, can naturally also be derived within this framework [42], but they will not be discussed here. Nuclear forces in chiral EFT are separated into pion-exchange contributions and contact terms. Pion-exchange contributions represent the long-and intermediate-range parts of nuclear interactions and contain all chiral physics. Contact terms, on the other hand, encode the unresolved short-range physics and their strength is specified by unknown low-energy constants (LECs), that need to be adjusted to experimental data.
At LO, besides the already mentioned OPE potential, there are two NN contact terms with no momentum dependence that contribute only to the S-wave. They are identified by the four-nucleon-leg diagram with a momentum-independent vertex denoted by a small dot in the first row of Fig. 1. The interaction at LO is a very simple approximation, but already takes into account some of the important features of the NN force. For instance, the OPE generates the tensor component of the nuclear force known to be crucial to properly describe the two-nucleon bound state (deuteron).
The leading NN two-pion-exchange (TPE) contributions appear at NLO. Diagrams involving virtual excitations of the ∆-isobars [56,57,58,59] also appear at NLO in the delta-full chiral-EFT approach. Most importantly, seven new momentum-dependent contact terms can be constructed at this order, which are denoted by the four-nucleon-leg graph with a solid square in the second row of Fig. 1. These additional contact terms are important to correctly describe NN scattering in the Sand P -waves. More details about these contributions are presented in the next sections. Another important contribution at NLO is the leading 3N force, which can be described by the well-known Fujita-Miyazawa diagram [34], which involves intermediate excitation of the ∆-isobars between three nucleons. While this contribution has to be considered in the ∆-full approach, it can be shown that the net contribution of 3N forces vanishes in the delta-less chiral EFT [39,49] at this order. At next order, N 2 LO, the sub-leading NN TPE diagrams contain vertices (large solid dots) proportional to the so-called c i coefficients. The values of these parameters can be obtained by pion-nucleon (πN ) [60,61,62,63,64,65,66,67] or NN scattering data [41]. In the delta-less chiral EFT, these coefficients mimic the effect of the ∆-isobar (or some other meson resonances) through a mechanism known as resonance saturation. Hence, they are enhanced in magnitude and found to be "unnaturally" large. The explicit inclusion of the ∆ isobar in the delta-full theory reduces the strength of the c i 's and promotes the corresponding contributions to a lower order (see gray arrows in Fig. 1). As a consequence, the convergence of the expansion in the delta-full theory improves considerably at these orders. In the delta-full approach, additional sub-leading TPE contributions appear that have also been worked out at this order [60].
In addition to the NN sector, additional 3N diagrams appear at N 2 LO in both approaches. They involve a 3N TPE, a OPE-contact interaction, and a true 3N contact diagram. The 3N TPE potential also involves the c i parameters already present in the TPE NN force. As in the case of the NN force, these contributions absorb the presence of the ∆-isobar in the delta-less approach, while some of their strength is promoted to lower order in the form of the already discussed Fujita-Miyazawa diagram in the delta-full approach. The OPE-contact and 3N contact diagrams include two purely three-body LECs that have to be adjusted to A ≥ 3 data. Finally, the are no additional diagrams due to ∆ contributions to the 3N force at N2LO [68].
At higher orders, the number of contributions to the NN force dramatically increases. In Fig. 1 only a few representative diagrams are displayed. For instance, at N 3 LO more TPE contributions occur -in both delta-less and delta-full chiral EFT -involving leading two-loop and relativistic corrections. In addition, leading three-pion (3π) exchange contributions arise at this order but they are found to be negligible. The main feature at N 3 LO is the presence of additional contact interactions represented by the four-nucleon-leg with a solid diamond. Since these interactions are ∼ p 4 , p 4 , they have a relevant impact up to the D waves. Their full operator structure will be discussed in the next section. Additional complicated 3N diagrams appear at N 3 LO, as well as the first contributions to four-nucleon forces (4N ). We will not discuss these diagrams here and refer the reader to Refs. [69,70,71,72]. For additional contributions at N 4 LO and N 5 LO, we refer the interested reader to Refs. [73,74,75,76].
An important aspect of nuclear interactions (and currents) in the basic model is that they suffer from ultraviolet (UV) divergences which need to be removed by a proper regularization and renormalization procedure. There are two sources of UV divergences that require regularization: first, UV divergences appear in loop corrections, and second when solving the Schrödinger or Lippmann-Schwinger equations or when calculating matrix elements involving nuclear currents. Loop divergences can be treated via dimensional regularization (DR) or spectral-function regularization (SFR), where the latter is implemented through the inclusion of a finite cutoff in the spectral functions. To cure divergences when solving the Schrödinger or Lippmann-Schwinger equations, the nuclear potential is multiplied by regulator functions that remove large-momentum contributions above a chosen cutoff scale. The regularization of the potential (and current) operators is followed by a renormalization procedure, i.e., dependencies on the regularization scheme and cutoff are reabsorbed, order by order, by the LECs entering the potential (and currents).
Nucleon-nucleon scattering has been extensively studied in chiral EFT in the past two decades following the pioneering work by Weinberg [37,38,39] and Ordonez et al. [58]. In particular, NN potentials at N 3 LO in the chiral expansion are available since the early 2000's [77,78] and have served as a basis for numerous ab initio calculations of nuclear structure and reactions. Recently, accurate and precise chiral EFT potentials up to fifth order in the chiral expansion, i.e. N 4 LO, have been developed [73,74,75,76], and provide an extremely accurate description of NN data bases up to laboratory energies of 300 MeV with a χ 2 per datum close to one. These databases have been provided by the Nijmegen group [29,26], the VPI/GWU group [79], and more recently the Granada group [80,81,82]. In the standard optimization procedure, the NN potentials are first constrained through fits to neutron-proton (np) and proton-proton (pp) phase shifts, and then refined by minimizing the total χ 2 obtained from a direct comparison with the NN scattering data. However, new optimization schemes are being explored in Refs. [83,84]. For instance, the optimization strategy of the N2LO sat interaction of Ref. [84] is based on a simultaneous fit of the twoand three-nucleon forces to low-energy NN data, the deuteron binding energy, and binding energies and charge radii of hydrogen, helium, carbon, and oxygen isotopes using consistent NN and 3N interactions at N 2 LO. However, despite the good description of properties of 16 O and 40 Ca, the NN component of this interaction shows deficiencies in reproducing the pp and np scattering data even at very low energy.
Three-nucleon forces and their impact on nuclear structure and reactions has become an important frontier in nuclear physics, see Refs. [85,86] for review articles. As shown in Fig. 1, chiral contributions to the 3N interaction have been derived up to N4LO in the chiral expansion [69,70,87,88,89]. However, few-and many-nucleon calculations are, with very few exceptions, still limited to chiral 3N forces at N 2 LO. At this order, as we have mentioned above, 3N forces are characterized by the presence of two unknown LECs that have to be determined. The two LECs -namely c D in the OPE-contact and c E in the 3N contact interaction -have been constrained either by fitting exclusively strong-interaction observables [90,91,92,93] or by relying on a combination of strong-and weak-interaction observables [94,95,96]. This last approach is made possible by the relation between c D in the OPE-contact interaction and the LEC in the NN contact axial current [94,95], established in chiral EFT [97]. This connection allows one to use nuclear properties governed by either strong or weak interactions to constrain simultaneously the 3N interaction and NN axial current.
As chiral EFT is a low-momentum expansion of nuclear interactions, many of the chiral interactions available in the literature are formulated in momentum space and have the feature of being strongly nonlocal in coordinate space. This makes them not well-suited for certain numerical algorithms, for example QMC methods. In this context, an interaction is local if it depends solely on the momentum transfer q = p − p , which Fourier transforms to dependencies on r. However, interactions in momentum-space can also depend on the momentum scale k = (p + p)/2, which Fourier transform to derivatives in coordinate space. These k dependencies, and thus non-localities, come about because of (i) the specific functional choice made to regularize the momentum space potentials in terms of the two momentum scales p and p , and (ii) contact interactions that explicitly depend on k. QMC methods, for example variational (VMC) and Green's Function Monte Carlo (GFMC) [10,98] techniques, provide reliable solutions of the many-body Schrödinger equation -presently for up to A = 12 nucleons -with full account of the complexity of the many-body, spin-and isospin-dependent correlations induced by nuclear interactions. The sampling of configuration space in VMC and GFMC simulations gives access to many important properties of light nuclei such as spectra, form factors, transitions, low-energy scattering, and response functions. Auxiliary Field Diffusion Monte Carlo (AFDMC) [10,98] uses Monte Carlo techniques to additionally sample the spin-isospin degrees of freedom, enabling studies of, for example, nuclei up to A = 16 [99,100] and neutron matter [90,91,101,102,103] that is so critical to determining the structure of neutron stars. QMC simulations have surely proved to be very valuable in attacking many nuclear-structure problems over the last three decades but require local chiral interactions as input. Therefore, there is a need to develop local chiral interactions for the use in QMC methods in order combine these accurate many-body methods with systematic nuclear interactions and to test to what extent the chiral EFT framework impacts our knowledge of few-and many-body systems.

Local two-nucleon interactions
A major thrust of our work is based on the theoretical derivation, optimization, and implementation of chiral interactions suitable for QMC methods. In recent years, local configuration-space chiral NN interactions have been derived by two groups [104,105,106,107]. In this section, we will introduce these two families of interactions, that are either derived in the delta-less [104,105] or delta-full [106,107] approach. We begin by introducing general features of both approaches and then describe the specifics. We will be stating general considerations in momentum-space, where q dependencies indicate local parts of interactions and k dependencies indicate non-localities, and then switch to coordinate-space where interactions are local if they only depend on the relative distance r = r i − r j . Fourier transformations connect interactions in momentum-and coordinate-space, with q and r being associated variables, while k leads to appearances of gradient terms.
As discussed before, nuclear interactions can generally be separated into different interaction channels depending on their operator structure. Obviously, chiral interactions can also be separated into long-range physics, mediated by pion-exchange interactions, and short-range physics, which is described by a set of operators consistent with all symmetries and accompanied by LECs adjusted to reproduce experimental data: Each of these components can then be expanded in chiral order ν as discussed before: At LO, ν = 0, both delta-less and delta-full chiral EFT have the same operator structure. At this order, only the leading contact interactions as well as the one-pion exchange (OPE) interaction contribute, see Fig. 1. Generally, pion-exchange interaction can be written as with central, spin, tensor, spin-orbit and quadratic spin-orbit components, respectively. In the local chiral interactions discussed in this review, the spin-orbit and quadratic spin-orbit terms are not included as they Frontiers are of higher order. The one-pion exchange interaction is given in momentum space as where g A , f π = 92.4 MeV, and m π denote the axial-vector coupling constant of the nucleon, the pion decay constant, and the pion mass, respectively. As a consequence, the OPE contributes to the W T channel.
Including isospin-symmetry breaking effects induced by the mass difference between charged and neutral pions, the OPE interaction can be rewritten as Hence, when including isospin-symmetry breaking, the OPE adds to the W S and W T parts of Eq. (7).
with Y α (q) and T α (q) given by Here, m πα denotes the neutral (m π 0 ) and charged (m π ± ) pion masses. When Fourier-transformed, the coordinate-space OPE is given by (11) where the individual functions can be obtained from Eq. (9) with q → r and with the functions Y α (r) and T α (r) given by Here, x α = m πα r. Note that Eq. (11) only holds in the case r > 0. In addition, upon Fourier transformation a δ-function appears, which has been dropped from Eq. (11), because it can be reabsorbed in the short-range contact terms at LO, which we will discuss next.
The LO contact interactions are momentum-independent and can be described by the most general operator set allowed by all symmetries: As these terms describe the interactions of nucleons, i.e., fermions, these interactions are used between antisymmetrized wave functions. One can define the antisymmetrized interaction by applying the antisymmetrizer, given by One then finds It follows immediately that only two out of these four couplings are linearly independent, describing the two possible S-wave scattering channels. The two commonly chosen LO contact operators are but in principle any different two of the four contact interactions can be chosen and lead to the same physical description for fermionic systems. This is analogous to Fierz ambiguities and in the following we will call this freedom to choose operators Fierz rearrangement freedom.
Additionally, there are isospin breaking corrections to the LO contact interactions that have to be taken into account. These are due to different masses of u and d quarks, and account for differences in neutron-neutron (nn), np, and pp S-wave scattering lengths: At higher orders, the description of the potential changes depending on the choice of delta-less or delta-full approach. In the following, we will describe both approaches as pursued by individual research groups.

Without Delta isobars
At NLO in chiral EFT, additional momentum-dependent contact interactions as well as TPE interactions appear. For the TPE, we give the expressions within the spectral-function representation (SFR) as detailed in Ref. [108], with spectral functions ρ i and η i : Here,Λ is the SFR cutoff. Similar expressions are valid for W C , W S , and W T in terms of η C , η S , and η T . The TPE spectral functions at NLO are given by [109] ρ (2) For the NLO contact interactions, the most general set of operators is given by Using the same arguments as for the LO contact interactions, only 7 out of these 14 operators are linearly independent. To construct local interactions, one typically chooses the 6 local operators (proportional to γ 1 -γ 4 , γ 11 , and γ 12 ) as well as the spin-orbit operator (proportional to γ 9 ): In coordinate space, this translates to Piarulli et al.
The relativistic 1/m N corrections, with m N being the nucleon mass, have been omitted here since, in the counting employed here, they would appear at N 3 LO, provided the nucleon mass is counted according to Q/m N ∼ Q 2 /Λ 2 b as suggested in Ref. [38].
The delta-less chiral EFT approach has been used to construct local interactions up to N 2 LO. At next higher order, N 3 LO, contact interactions cannot be written down in a purely local fashion, as only 8 out of 30 possible operators are local. A possible way forward is the definition of "maximally local" N 3 LO potentials, which has been pursued in the delta-full approach and will be discussed in the next section.
Finally, it is necessary to specify a regulator scheme. For the delta-less local interactions of Refs. [104,105], the following long-and short-range regulators are used: The long-range regulator multiplies each function Y (r), while the short-range regulator replaces all δ functions. The regulator functions depend on the cutoff scale R 0 , that determines how long-and short-range physics are separated. For a smaller cutoff R 0 (i.e., for a larger momentum-space cutoff), the interactions is probed at shorter distances, and typically shows stronger short-range repulsion. We show the delta-less local chiral interactions in the 1 S 0 channel in Fig. 2 for different values of the cutoff. Introducing a local regulator function leads to the appearance of regulator artifacts that brake Fierz-rearrangement freedom. We will address this topic in detail in Section 4.1.

With Delta isobars
In the delta-full local chiral interactions, coordinate-space expressions for the TPE terms at NLO and N 2 LO are obtained by using the spectral function representation [109,108] but with dimensional regularization (DR) [59]. This implies taking the cutoffΛ in Eqs.
where three values for the radius R L are considered: R L = (0.8, 1.0, 1.2) fm with the diffuseness a L fixed at a L = R L /2 in each case.
Another difference between the delta-less and delta-full coordinate-space interactions lies in the operator structure of their short-range components. In the delta-full potentials, selected contact terms at N 3 LO are also retained in addition to the LO and NLO contributions given in Eqs. (16) and (26). The contact potential at order N3LO, V N3LO cont (q, k), which involves four gradients acting on the nucleon fields, is expressed in terms of 15 independent operators [41] after considering the Fierz rearrangement freedom. Its standard parametrization, adopted in momentum-space potentials, is given by However terms proportional to k 2 and k 4 in those expressions, upon Fourier transformation, would lead to gradient operators in coordinate-space (p −→ −i∇ is the relative momentum operator), making the NN potential strongly non-local.
The number of non-localities can be reduced by reconsidering the Fierz rearrangement freedom. However some of these non-local terms still persist at N 3 LO leading to the definition of "minimally nonlocal" contact interactions: In coordinate space, this reads as where O l=1,...,11 referred to as c, τ , σ, στ , t, tτ for the first six operators, and b, bτ , bb, q, qσ for the remaining five operators. The four additional terms, denoted as p, pσ, pt, and ptτ , in the anti-commutator of Eq. (36) are p 2 -dependent. For the definition of the radial functions v l S (r) as well as those multiplying the p 2 -terms, we refer the reader to Appendix A.
A comment is now in order. The strict adherence to power counting would require the inclusion of additional one-loop as well as two-loop TPE and three-pion exchange contributions at N 3 LO. For the time being, these contributions have been neglected, since part of their strength is promoted at lower orders due to the inclusion of the ∆ resonance, and some of the remaining diagrams are also known to be small (see, for example, Ref. [41]). Furthermore it is the D i LEC's at N 3 LO that are critical for a good reproduction of phase shifts in lower partial waves, particularly D-waves, and a good fit to the NN database. However, the consistency between the long-and short-range part at higher orders in the delta-full chiral EFT is work in progress.
The local versions of these "minimally nonlocal" NN potentials have been defined by dropping the terms proportional to p 2 in the anti-commutator when the optimization procedure for estimating the LECs is carried out [107]. In Ref. [107] we observed that the inclusion of the p 2 -dependent terms would have improved the fits to the database in the laboratory energy range up to 200 MeV only marginally. However, apart from the small improvement that the p 2 -dependent terms would bring to the total χ 2 in the fit to the NN scattering data, the effect of these terms on nuclear observables has not been studied.
Lastly, the delta-full local interactions contain additional isospin breaking terms at NLO. They are parametrized by the following operators O l=12,...,16 referred to as τ z, T , σT , tT , bT . The radial functions multiplying these operators are also reported in Appendix A.
The short-range part of these potentials involve the local regulator given in Eq. (32) with n = 2, where we consider, in combination with R L = (0.8, 1.0, 1.2) fm, R S = (0.6, 0.7, 0.8) fm, corresponding to typical momentum-space cutoffs Λ S = 2/R S ranging from about 660 MeV down to 500 MeV. Hereafter, we will denote the potential with cutoffs (R L , R S ) = (1.2, 0.8) fm as model a, that with (1.0, 0.7) fm as model b, and that with (0.8, 0.6) fm as model c.
There are 26 LECs in the definition of the delta-full local interactions. Of these, 20 LECs describe the charge-independent part of the interaction: 2 at LO (Q 0 ), 7 at NLO (Q 2 ), and 11 at N 3 LO (Q 4 ). The remaining 6 LECs describe its charge-dependent part: 2 at LO (one each from CIB and CSB), and 4 at NLO from CIB. The optimization procedure to fix these 26 LECs uses pp and np scattering data (including normalizations), as assembled in the Granada database [80], the NN scattering length, and the deuteron binding energy. The minimization of the χ 2 objective function with respect to the LECs is carried out with the Practical Optimization Using no Derivatives (for Squares) routine, POUNDerS [110]. For each of three different sets of cutoff radii (R S , R L ), two classes of local interactions have been developed, which only differ in the range of laboratory energy over which the fits were carried out, either 0-125 MeV in class I or 0-200 MeV in class II. The χ 2 /datum achieved by the fits in class I (II) was 1.1( 1.4) for a total of about 2700 (3700) data points. In the literature, we are referring to these NN interactions generically as the Norfolk potentials (NV2s), and designate those in class I as NV2-Ia, NV2-Ib, and NV2-Ic, and those in class II as NV2-IIa, NV2-IIb, and NV2-IIc.
The NV2 interactions were found to provide insufficient attraction in calculations of the ground-state energies of nuclei with A = 3-6 [107]. To remedy this and similar shortcomings, 3N interactions at N 2 LO have to be included in both approaches. This will be described in the next section.

Local three-nucleon interactions
Three-nucleon forces are very important ingredients for the correct description of physical systems. They naturally appear within chiral EFT and are consistent with the NN sector. The exact description of the 3N interactions depends on the choice of delta-less vs. delta-full approach. In the following, we review 3N forces in both approaches.

Without Delta isobars
In the delta-less chiral EFT approach, the leading 3N contributions appear at N 2 LO in the power counting. They an be separated into three topologies: (i) a long-range TPE interaction named V C depending on the pion-nucleon LECs c 1 , c 3 , and c 4 , that already appear in the NN sector, (ii) a one-pion-exchange-contact interaction V D dependent on a new LEC c D , and (iii) a 3N contact interaction V E dependent in a new LEC c E . The LECs c D and c E solely describe 3N physics and need to be adjusted to properties of A ≥ 3 systems. In momentum space, these interactions are defined as where we sum over all permutations of the particles i, j, and k, where the first pion carries a momentum q i from nucleon i to j, while the second pion carries q k from j to k, and where F αβ ijk is given by As one can easily see, all of these interactions are local, as long as local regulator functions are applied. To obtain expressions in coordinate space, these interactions have to be Fourier transformed. For the part of V C proportional to c 1 , we find This results in where we have used and Y (r) = exp(−m π · r) m π r , U (r) = 1 + 1 m π r .

Frontiers
For the other parts of V C we find where and For the one-pion-exchange-contact part V D we find and for the three-nucleon-contact interaction V E we find To regularize these 3N topologies, we choose consistent regulators with the NN sector, i.e., we replace δ functions by f short (r) and multiply Yukawa functions with f long (r). The cutoff scale for 3N interactions does not necessarily have to be the same as for the NN sector, and we call it R 3N in the following.
To adjust the appearing 3N couplings to experimental data, one should select few-body observables that are uncorrelated. In the delta-less approach, these observables have been chosen to be the 4 He binding energy and n-α scattering P wave phase shifts; see Fig. 3, where we show parameter curves for the 3N LECs for different 3N cutoffs R 3N , chosen similar to R 0 , and for different parametrizations that we will discuss in the next section. Stars in the parameter curves mark fits that also describe neutron-alpha scattering, shown in the right panel. For more details, see Ref. [91].

With Delta isobars
In the delta-full chiral EFT approach, the structure of the 3N force at N 2 LO is similar to the 3N force in the delta-less approach. We still have the three topologies V C , V D and V E at N 2 LO but, in addition, the well-known Fujita-Miyazawa interaction [34] (V ∆ ), which in the delta-less approach is absorbed by V C , appears already at NLO in the power counting. In momentum space, it reads as where S, S † and T , T † are the transition spin and isospin operators: The operator S (T ) converts a spin (isospin) 1/2 into a spin (isospin) 3/2 particle.
The configuration-space expression follows from where the following definitions have been introduced: and the dimensionless functions Y (r) and T (r) defined before. The term [ · · · ] in Eq. (56) can be written as and the transition-spin and transition-isospin operators can be eliminated using the identities to obtain where the function X ij was defined in the previous section. In the definitions above, the δ(r)-function terms have been dropped.
In analogy to the 3N delta-less chiral EFT, we regularize the 3N contributions in the delta-full chiral EFT by replacing the δ functions with f ∆ short (r) and multiplying the Yukawa functions with f ∆ long (r). Note that the implementation of V C and V D in the delta-full chiral EFT does not retain the terms proportional to σ i · σ j in the definition of X ik , in Eq. (50), and in Eq. (52). They can be reabsorbed in the redefinition of the short-range contact terms.
In the delta-full chiral EFT, two different sets for the values of c D and c E were obtained, leading to two different parametrization of the 3N interaction [93,96]. In the first, these LECs were determined by simultaneously reproducing the experimental trinucleon ground-state energies and neutron-deuteron (nd) doublet scattering length, as shown in the left panel of Fig. 4. In the second set, these c D and c E were constrained by fitting, in addition to the trinucleon energies, the empirical value of the Gamow-Teller matrix element in tritium β decay [96], see right panel of Fig. 4. Because of the much reduced correlation between binding energies and the GT matrix element, the second fit procedure leads to a more robust determination of c D and c E then attained in the first one. Note that these observables have been calculated with hyperspherical-harmonics (HH) expansion methods [5] as described in Refs. [93,96].

FINITE CUTOFF AND REGULATOR ARTIFACTS
The derivations of local interactions in the last sections did not include any of the local regulator functions that necessarily have to be applied to the interactions to make them suitable for the use in nuclear manybody methods. Generally, when introducing a regulator function, terms beyond the order at which one is working are affected. Hence, the use of such regulator functions with finite values for the cutoff leads to the appearance of regulator artifacts, that might influence calculations of many-body observables. In this section, we will address the different regulator artifacts that can appear in calculations with local interactions.

Violation of Fierz-rearrangement freedom
The first regulator artifact for local interactions affects short-range operators. In previous sections we had shown how only half of the operators at each order are linearly independent due to their insertion between antisymmetric fermionic states, see, e.g., Eq. (15). However, this argument changes when a regulator function is applied. The discussion in this section will follow Ref. [111].
In general, a regulator function can depend on two momentum scales, f R (q, k). Local regulators, on the other hand, only depend on q, f R,loc (q). The derivation of Eq. (15) remains valid if the regulator function commutes with the antisymmetrizer and, hence, reduces to a simple prefactor in Eq. (15), i.e., when We can immediately see, that a purely local regulator can never fulfill this condition while typical nonlocal regulators of the form [84,83,76,77,78] do. As a consequence, and the Fierz-rearrangement freedom is violated. For general regulator functions as defined in the previous sections, this leads to where V f corr (p · p ) captures all the regulator artifacts that are of higher-order in the EFT. It depends on the functional form of the regulator and the cutoff value. One can also see, that the corrections can be angle-dependent, which leads to a mixing of different partial waves. As a consequence, when applying these regulators to a three-neutron system, for example, pure contact interactions, that otherwise would vanish due to the Pauli principle, start to contribute. This mixing of partial waves complicates the fitting procedure, increases theoretical uncertainties, and makes calculated observables dependent on the operator structure that was chosen.
In Fig. 5 we show results for the 4 He ground-state energy for different LO operator choices. As one can see, the ground-state energies can vary by approximately 10 MeV at LO, depending on the operator choice. However, when going to higher order and including subleading contact operators, regulator artifacts get partially absorbed and corrected. Then, only higher-order artifacts remain, which improves the situation considerably, as can be seen for the NLO results. In this case, the spread originating from different choices of LO operators reduces to approximately 4 MeV.
A similar effect appears in the 3N sector, where the V E contact interaction suffers from a similar violation of the Fierz freedom when local regulators are applied. While 3N forces are typically fit to symmetric systems where this dependence can then be approximately accounted for, in triples with S = 3/2 or T = 3/2 (where typically no 3N contact force can contribute due to the Pauli principle) regulator artifacts appear, and lead to a finite contribution from 3N contact interactions that depend on the operator choice. We show this behavior in Fig. 5 in the right panel in the case of pure neutron matter, where all triples have T=3/2. The three different bands explore 3 choices for the 3N contact operators. At nuclear saturation density, we find that the regulator artifacts introduce a spread of approximately 5 MeV. Unfortunately, higher-order correction terms only appear at N 4 LO and, to date, are not systematically included in any calculation.
Finally, we mention that the finite cutoff also introduces an ambiguity in the V D term, that depends on the choice of the initial spin-isospin structure when Fourier transforming: Both expressions are identical for true δ functions (infinite cutoff) but differ when a finite cutoff is applied.

Weaker pion exchanges
A second regulator artifact for local regulators affects the pion exchanges. In Ref. [90] it was shown that locally regulated pion exchanges lead to less 3N repulsion than nonlocally regulated pion exchanges. At the Hartree-Fock level, for a typical cutoff of 2.5fm −1 , when applying nonlocal regulators ≈ 97% of the infinite cutoff result is recovered, while local regulators only recover ≈ 60%. To reproduce the momentum-space results, the cutoff has to be considerably increased.
Local regulators for pion exchanges have been investigated in detail in Ref. [112] in both the NN and 3N sector. The fact that the contribution due to pion exchanges is weaker for local than for nonlocal regulator functions is easy to understand in the Hartree-Fock approximation. At the Hartree-Fock level, there are both a direct and an exchange term. The momentum transfer q = p − p vanishes in the direct term because p = p , but it is q = 2p in the exchange term because p = −p . A typical local regulator of the form exp − q Λ n , thus, evaluates to 1 in the direct term, but to exp − p Λ/2 n in the exchange term. Therefore, compared to nonlocal regulators for which both terms are identical, exp − p Λ n , local regulators have a very different behavior. In particular, local regulators have an effectively lower cutoff in the exchange channel. In the Hartree-Fock approximation, where the direct term vanishes for spindependent interactions like pion exchanges, only the exchange term contributes and, hence, is weaker for local than for nonlocal regulators.
While the situation is more complicated when abandoning the Hartree-Fock approximation, this reasoning qualitatively remains valid and locally regulated pion exchanges are weaker than nonlocally regulated pion exchanges.

SELECTED RESULTS
In this section, we will briefly show the successes of Quantum Monte Carlo calculations with local chiral interactions for light atomic nuclei and infinite matter.   [99] with permission from the American Physical Society.

Light nuclei
Local chiral interactions, both in the delta-less and delta-full approach, have been used to successfully describe properties of light nuclei using QMC methods. In Fig. 6, we show GFMC results for groundand excited states for nuclei up to 12 C within the delta-full approach compared to experimental data. In addition, the results obtained with chiral EFT are compared to results with phenomenological interactions. The results clearly show that chiral interactions describe spectra of light nuclei with great success and are compatible to the accuracy of phenomenological interaction in these systems. In addition, we also show ground-state energies obtained in the AFDMC method for nuclei up to 16 O for delta-less chiral interactions. Results are given at LO, NLO, and N 2 LO for two different 3N parameterizations to explore regulator artifacts. Again, chiral interactions agree well with experimental results, which are shown as green points.
In addition to energies, local chiral interactions describe charge radii well. In Fig. 7, we present orderby-order AFDMC results for the charge radii of nuclei up to 16 O, compared to experiment. Again, the description is accurate. In addition, as mentioned before, delta-less chiral interactions have been adjusted to reproduce neutron-α scattering phase shifts, see Fig. 3. While NN interactions alone cannot reproduce the P wave splitting in this system (NLO calculations in Fig. 3), chiral Hamiltonians at N 2 LO, including 3N interactions, reproduce the neutron-α P wave scattering phase shifts accurately.

Infinite Matter
In addition to properties of atomic nuclei, local chiral interactions have been used to study infinite matter, and in particular, pure neutron matter. In the right panel of Fig. 5, we have already shown results for the energy per particle of pure neutron matter. Results are shown for 3 Hamiltonians at N 2 LO, that explore the uncertainty due to regulator artifacts and the truncation of the chiral series. While uncertainties in pure neutron matter are enhanced due to the local regulator artifacts discussed before, indicated by the differences between the three bands, the resulting neutron-matter equation of state (EOS) is consistent with other ab initio determinations within uncertainties.
These calculations have been successfully used to study the EOS of neutron stars, and it has been found that the resulting equations of state are consistent with astrophysical observations of pulsar masses. The EOS have also been used to study gravitational waves from neutron-star mergers [102,113,114].

CONCLUSION AND OUTLOOK
The quest to understand properties of nuclear systems in terms of forces acting between the nucleons has been considered one of the most challenging efforts of nuclear theory. During the past quarter century, particular emphasis has been devoted to the systematic framework provided by chiral EFT. This approach allows for a consistent description of the two-and many-body interactions and ensuing many-body currents, and a quantification of the theoretical uncertainty due to the truncation error in the chiral expansion.
In this review, we have presented a comprehensive description of the two families of local chiral interactions that have been developed for the use in QMC methods: one within the delta-less and one within the delta-full approach. We provided many details about the theoretical derivation and optimization of these nuclear models addressing their similarities and differences. For completeness, we also presented selected QMC results for light nuclei and neutron matter. These results show that the combination of local chiral EFT interactions with powerful QMC many-body methods can accurately describe ground-and excited-state energies, radii of nuclei up to 16 O, and n-α scattering, as well as the equation of state of neutron matter.
These local chiral interactions have also been used to calculate the distribution of nucleons in a nucleus in both momentum and coordinate space which are related to experimental observations [99,100,115,116], in benchmark calculations of the energy per particle of pure neutron matter as a function of the baryon density [103] and in studies of neutrinoless double-beta decays [117].
In future, local chiral interactions will continue to serve as input for precise QMC methods to systematically study, for example, electroweak reactions, along the lines of Refs. [118,119,120,121,122,123,124] and infinite matter also at finite proton fractions.
Improvements to the interactions that reduce uncertainties due to the scheme and scale dependence of the interactions, e.g., the inclusion of higher orders in the chiral expansion in both the N N and 3N sectors, will provide exciting prospects and permit precision studies of many nuclear systems.
where each of these functions is multiplied by the cutoff f ∆ long (r) defined in Section 3.1.2. Below, we give the coordinate-space representation of the TPE radial functions at NLO and N2LO that contribute to V C,π , W C , V S , W S , V T , and W T .
At NLO, the terms without ∆'s read as where g A , f π , and m π denote the axial-vector coupling constant of the nucleon, the pion decay constant, and the pion mass, respectively, x = m π r and K n are modified Bessel functions of the second kind. The NLO terms with a single ∆ intermediate state are given by v 2π,NLO c (r; 1∆) = − 1 6π 2 r 5 y v 2π,NLO στ (r; 1∆) = 1 54π 2 r 5 y v 2π,NLO tτ (r; 1∆) = − 1 54π 2 r 5 y where h A is the N -to-∆ axial coupling constant, y = ∆r (∆ is the ∆-nucleon mass difference) and the parametric integral over µ is carried out numerically. The NLO terms with 2 ∆ intermediate states are

Contact terms
The contact expressions introduced in Section 3.1.2 are defined by the functions v l S (r). These functions are given by