Atomic nuclei from quantum Monte Carlo calculations with chiral EFT interactions

Quantum Monte Carlo methods are powerful numerical tools to accurately solve the Schr\"odinger equation for nuclear systems, a necessary step to describe the structure and reactions of nuclei and nucleonic matter starting from realistic interactions and currents. These ab-initio methods have been used to accurately compute properties of light nuclei -- including their spectra, moments, and transitions -- and the equation of state of neutron and nuclear matter. In this work we review selected results obtained by combining quantum Monte Carlo methods and recent Hamiltonians constructed within chiral effective field theory.


INTRODUCTION
The study of nuclear properties as they emerge from the individual interactions among protons and neutrons is a fascinating long-standing problem, subject of both theoretically and experimentally research activities. From a theoretical point of view, a truly ab-initio description of nuclei is still very challenging at present. The underlying theory of strong interactions, Quantum Chromodynamics (QCD), that describes how quarks and gluons interact to form nucleons and nuclei, in the low-energy regime is non-perturbative in its coupling constant. Despite remarkable progresses [1,2], realistic computations of many-body nuclear systems in terms of the fundamental degrees of freedom of QCD -quarks and gluons -are still extremely challenging.
A more feasible approach to the problem consists in assuming that at the energy regime relevant to the description of atomic nuclei, quarks, and gluons are confined within hadrons. The latter are the active degrees of freedom at soft scales, and they interact among themselves through non-relativistic effective Several many-body methods have been developed to numerical solve the many-body Schrödinger equation. Most of them rely on basis expansions, for example the coupled cluster method [5,6], the no core shell model [7], the similarity renormalization group [8], and the self consistent Green's function [9]. Each of these methods has distinct advantages, and many are able to treat a wider variety of nuclear interaction models. These many-body techniques are very effective and achieve a good convergence only when relatively soft potentials are used.
Quantum Monte Carlo (QMC) methods are ideally suited to study strongly correlated many-body systems, and have no difficulties in treating "stiff" nuclear interactions, but are limited to nearly local nuclear potentials. For this reason, until fairly recently, the applicability of QMC methods was limited to phenomenological interactions, as χEFT Hamiltonians were typically written in momentum space. Over the past few years, the situation has drastically changed with the development of local χEFT potentials, both with [10,11] and without explicit delta degrees of freedom [12,13], that have provided a way to combine an EFT-based description of nuclear dynamics with precise QMC techniques. In this work we will review selected results of nuclei obtained using QMC methods and chiral Hamiltonians.

NUCLEAR INTERACTIONS
The microscopic model of nuclear theory assumes that nuclear systems can be described as point-like nucleons, whose dynamics is characterized by a non-relativistic Hamiltonian where T i is the one-body kinetic energy operator, v ij is the nucleon-nucleon (NN ) interaction between particles i and j, V ijk is the three-nucleon (3N ) interaction between particles i, j, and k, and the ellipsis indicate interactions involving more than three particles. There are indications that four-nucleon interactions may contribute at the level of only ∼ 100 keV in 4 He [14] or pure neutron matter [15], and therefore are negligible compared to NN and 3N components. Hence, current formulations of the microscopic model do not typically include them (see, for example, reference [4]).
The NN interaction term in the nuclear Hamiltonian is the most studied of all, with thousands of experimental data points at laboratory energies (E lab ) from essentially zero to hundreds of MeV. It consists of a long-range component, for inter-nucleon separation r 2 fm, due to one-pion exchange (OPE) [16], and intermediate-and short-range components, for, respectively, 1 fm r 2 fm and r 1 fm, derived, up to the mid 1990's, almost exclusively from meson-exchange phenomenology [17,3,18]. These models fit the large amount of empirical information about NN scattering data contained in the Nijmegen database [19], available at the time, with a χ 2 /datum 1 for E lab up to pion-production threshold. Two well-known and still widely used examples in this class of NN interactions are the CD-Bonn [18] and the Argonne v 18 (AV18) [3] potentials.
The AV18 interaction is a local, configuration-space NN potential that has been extensively and successfully used in a number of QMC calculations. It is expressed as a sum of electromagnetic and OPE terms and phenomenological intermediate-and short-range parts: The electromagnetic term v γ ij has one-and two-photon exchange Coulomb interaction, vacuum polarization, Darwin-Foldy, and magnetic moment terms, with appropriate form factors that keep terms finite at r = 0 (see reference [3] for more details). The OPE part includes the charge-dependent (CD) terms due to the difference in neutral (m π 0 ) and charged pion (m π ± ) masses, and in coordinate-space it reads where σ adn τ are the Pauli matrices that operate over the spin and isospin of particles, and S ij = 3 σ i ·r ij σ j ·r ij − σ i · σ j and T ij = 3 τ iz τ jz − τ i · τ j are the tensor and isotensor operators, respectively. The functions, v π στ (r), v π tτ (r), v π, σT (r), and v π tT (r) are defined as v π στ (r) = Y 0 (r) + 2 Y + (r) 3 , v π tτ (r) = T 0 (r) + 2 T + (r) 3 , v π σT (r) = where Y α (r) and T α (r) are the Yukawa and tensor functions given by with x α = m πα r, and g A = 1.267, f π = 92.4 MeV being the axial-vector coupling constant of the nucleon and the pion decay constant, respectively.
The intermediate-range region, v I ij , is parametrized in terms of two-pion exchange (TPE), based on, but not consistently derived from, a field-theory analysis of box diagrams with intermediate nucleons and ∆ isobars [20]. The short-range region, v S ij , is instead represented by spin-isospin and momentum-dependent operators multiplied by Woods-Saxon radial functions [3].
The AV18 model can be written as an overall sum of eighteen operators where the first eight are given by with the spin-orbit contribution expressed in terms of the relative angular momentum L = 1 2i (r i − r j ) × (∇ i − ∇ j ) and the total spin S = 1 2 (σ i + σ j ) of the pair. There are six additional charge-independent operators corresponding to p = 9 − 14 that are quadratic in L while the p = 15 − 18 are charge-independence breaking terms The AV18 model has a total of 42 independent parameters. A simplex routine [21] was used to make an initial fit to the phase shifts of the Nijmegen partial-wave analysis (PWA) [19], followed by a final fit direct to the database, which contains 1787 pp and 2514 np observables for E lab up to 350 MeV. The nn scattering length and deuteron binding energy were also fit. The final χ 2 /datum = 1.1 [3]. While the fit was made up to 350 MeV, the phase shifts are qualitatively good up to much larger energies, E ≤ 600 MeV [22].
Simplified versions of these interactions, including only a subset of the operators in Equation (7), are available. For instance, the Argonne v 8 (AV8 ) contains a charge-independent eight-operator projection, O p=1−8 ij = 1, σ i · σ j , S ij , L · S ⊗ 1, τ i · τ j , of the full NN potential, constructed to preserve the potential in all S and P waves as well as the 3 D 1 and its coupling to the 3 S 1 , while over-binding the deuteron by 18 keV due to the omission of electromagnetic terms [23]. The main missing features of these simplified interactions is the lack of terms describing charge and isospin symmetry breaking, as well as a slightly poorer description of nucleon-nucleon scattering data in higher partial waves. However, these contributions are very small, as outlined in reference [23].
Already in the 1980s, accurate three-body calculations showed that contemporary NN interactions did not provide enough binding for the three-body nuclei, 3 H and 3 He [24]. In the late 1990s and early 2000s this realization was also extended to the spectra (ground and low-lying excited states) of light p-shell nuclei, for instance, in calculations based on QMC methods [25] and in no-core shell-model studies [26]. Consequently, the microscopic model with only NN interactions fit to scattering data, without the inclusion of a 3N interaction, is no longer considered realistic.
In addition to NN forces, sophisticated phenomenological 3N interactions have been then developed. They are generally expressed as a sum of a TPE P -wave term, a TPE S-wave contribution, a three-pionexchange contribution, and a 3N contact [4]. More specifically, two families of 3N interactions were obtained in combination with the AV18 potential: the Urbana IX (UIX) [27] and Illinois 7 (IL7) [28] models. The UIX potential contains two parameters fit to reproduce the ground-state energy of 3 H and the saturation-point of symmetric nuclear matter, while the IL7 potential involves five parameters constrained on the low-lying spectra of nuclei in the mass range A = 3 − 10.
Despite their success in predicting a wide range of nuclear properties [4], the phenomenological potentials suffer from several drawbacks. For example, the resulting AV18+IL7 Hamiltonian leads to predictions of ≈ 100 ground-and excited-state energies up to A = 12 in good agreement with the corresponding empirical values. However, when used to compute the neutron-star equation of state, such Hamiltonian does not provide sufficient repulsion to guarantee the stability of the observed stars against gravitational collapse [29]. On the other end, the AV18+UIX model, while providing a reasonable description of s-shell nuclei and nuclear matter properties, it somewhat underbinds light p-shell nuclei.
Thus, in the context of the phenomenological nuclear interactions, we do not have a Hamiltonian that can explain the properties of all nuclear systems, from NN scattering to dense nuclear and neutron matter. Furthermore, this phenomenological approach does not provide a rigorous scheme to consistently derive two-and many-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.
These shortcomings were addressed when a new phase in the evolution of microscopic models began in the early 1990's with the emergence of χEFT [30,31,32]. χEFT is a low-energy effective theory of QCD and provides the most general scheme accommodating all possible interactions among nucleons and pions (∆-less χEFT) compatible with the relevant symmetries and symmetry breakings -in particular chiral symmetry -of low-energy QCD. In some modern approaches, the choice of degrees of freedom also includes the ∆ isobar (∆-full χEFT), because the ∆-nucleon mass splitting is only 300 MeV ∼ 2m π .
By its own nature, the χEFT formulation has an expansion in powers of pion momenta as its organizing principle. Most chiral interactions employed in recent nuclear structure and reaction calculations are based on Weinberg power counting. Within Weinberg power counting, the interactions are expanded in powers of the typical momentum p over the breakdown scale Λ b ∼ GeV, 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 references [33,34,35,36] for recent review articles). It is important mentioning that alternative power-counting schemes have been also suggested [37,38,39,40,41,42] but not fully explored.
This expansion introduces an order by order scheme, defined by the power ν of the expansion scale Q associated with each interaction terms: 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, χEFT provides a rigorous scheme to systematically construct many-body forces and consistent electroweak currents, and tools to estimate their uncertainties [43,44,45,46,47,48]. Chiral nuclear forces are comprised of both pion-exchange contributions and contact terms. Pionexchange contributions represent the long-range part of nuclear interactions and some of the pion-nucleon (πN ) couplings entering sub-leading diagrams can consistently be obtained from low-energy πN scattering data [49,50,51,52,53,54,55,56]. On the other end, contact terms encode the short-range physics, and their strength is specified by unknown low-energy constants (LECs), obtained by fitting experimental data. Similarly to the phenomenological interactions, the LECs entering the NN component are obtained by fitting NN scattering data up to 300 MeV lab energies, while the LECs involved in the 3N terms are fixed by reproducing properties of light-nuclei. This optimization procedure involves separate fit of the NN and 3N terms. Recently, a different strategy has been introduced by Ekström et al. [57]. This new approach is based on a simultaneous fit of the NN and 3N forces to low-energy NN data, deuteron binding energy, and binding energies and charge radii of hydrogen, helium, carbon, and oxygen isotopes.
Within χEFT, many studies have been carried out dealing with the construction and optimization of NN and 3N interactions [58,59,60,61,62,63,64,65,66,67,68,34,49,42,69,70,71,72,73,44,74,75,76,77,57,78,79,80,81] and accompanying isospin-symmetry-breaking corrections [82,83,84]. These interactions are typically formulated in momentum space, and include cutoff functions to regularize their behavior at large momenta. This makes them strongly non-local when Fourier-transformed in configuration space, and therefore unsuitable for use with QMC methods. In this context, an interaction is local if it depends solely on the momentum transfer q = p − p (p and p are the initial and final relative momenta of the two nucleons), which, upon Fourier transform, leads to dependencies solely on r. However, interactions in momentum-space can also depend on the momentum scale k = (p + p)/2, whose Fourier transform introduces derivatives in coordinate space. These k dependencies, and thus nonlocalities, 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.
In recent years, local configuration-space chiral NN interactions have been derived by two groups. On the one side, the authors of references [12,85] constructed NN local chiral potentials within ∆-less χEFT by including one-and two-pion exchange contributions and contact terms up to N 2 LO in the chiral expansion. The contact terms are regularized in coordinate space by a cutoff function depending only on the relative distance between the two nucleons, and use Fierz identities [86] to remove completely the dependence on the relative momentum of the two nucleons, by selecting appropriate combinations of contact operators. Their strength is characterized by 11 LECs, fixed by performing order by order χ 2 fit to NN phase shifts from the Nijmegen PWA up to 150 MeV lab energy. The fitting procedure is carried out for different values of the cutoff R 0 in the range of R 0 = 1.0 − 1.2 fm. The motivations why the authors of references [12,85] truncated the chiral expansion of these local potentials at N 2 LO is because at this order it is i) possible to have a fully local representation of the NN chiral interactions and ii) the inclusion of consistent 3N force is straightforward. In their models, the unknown 3N LECs are obtained by reproducing the binding energy of 4 He as well as the P -wave n − α elastic scattering phase shifts. In addition, they explore different parametrization for the 3N force, accordingly to Fierz identities [87,88,89]. In the present work, we are referring to a set of these local chiral interactions, specifically the (D2, Eτ ) model with R 0 = 1.0 fm of reference [89], as GT+Eτ -1.0.
On the other side, the authors of references [10,90] developed a different set of NN local chiral interactions by i) including diagrams with the virtual excitation of ∆-isobars in the TPE contributions up to N 2 LO (∆-full χEFT), ii) retaining contact terms up to N 3 LO. The LECs entering the NN contact interactions in this model are constrained to reproduce NN scattering data from the most recent and up-to-date database collected by the Granada group [91,92,93]. The contact terms are implemented via a Gaussian representation of the three-dimensional delta function with R S as the Gaussian parameter [90,10,94]. The pion-range operators are regularized at high-value of momentum transfer via a special radial function characterized by the cutoff R L [90,10,94]. There are two classes of these potentials. Class I (II) are fit to data up to 125 MeV (200 MeV). For each class, two combinations of short-and long-range regulators have been used, namely (R S , R L ) = (0.8, 1.2) fm (models NV2-Ia and NV2-IIa) and (R S , R L ) = (0.7, 1.0) fm (models NV2-Ib and NV2-IIb). Class I (II) fits about 2700 (3700) data points with a χ 2 /datum 1.1 ( 1.4) [90,10]. In conjunction with these models, two distinct sets of ∆-full 3N interactions have also been constructed up to N 2 LO. In the first, the 3N unknown LECs were determined by simultaneously reproducing the experimental trinucleon ground-state energies and neutron-deuteron (nd) doublet scattering length for each of the N N models considered, namely NV2-Ia/b and NV2-IIa/b [11,95]. In the second set, these LECs were constrained by fitting, in addition to the trinucleon energies, the empirical value of the Gamow-Teller matrix element in tritium β-decay [94]. The resulting Hamiltonians were labelled as NV2+3-Ia/b and NV2+3-IIa/b (or Ia/b and IIa/b for short) in the first case, and as NV2+3-Ia * /b * and NV2+3-IIa * /b * (or Ia * /b * and IIa * /b * ) in the second.
The interactions between external electroweak probes -electrons and neutrinos -and interacting nuclear systems is described by a set of effective nuclear currents and charge operators. Analogously to the nuclear interactions, electroweak currents can also be expressed as an expansion in many-body operators that act on nucleonic degrees of freedom. Electroweak currents have been developed in both meson-exchange and χEFT approaches. We refrain to discuss them in this work, redirecting the interested reader to dedicated reviews [4,96,97,98] and references therein.

QUANTUM MONTE CARLO METHODS
The χEFT Hamiltonians and the consistent electroweak currents discussed in the previous Section are the main input of sophisticated many-body methods aimed at solving with controlled approximations the nuclear many-body Schrödinger equation This is a highly non-trivial problem, mainly because of the non-perturbative nature and the strong spinisospin dependence of realistic nuclear forces. In this work, we will focus on QMC techniques, namely the variational Monte Carlo (VMC), the Green's function Monte Carlo (GFMC), and the auxiliary-field diffusion Monte Carlo (AFDMC) methods.

Variational Monte Carlo
The variational Monte Carlo method is routinely used to obtain approximate solutions to the manybody Schrödinger equation for a wide range of strongly interacting nuclear systems, including few-body nuclei, light closed-shell nuclei, and nuclear and neutron matter [4]. The VMC algorithm relies on the Rayleigh-Ritz variational principle to find the optimal set of variational parameters defining the trial wave function Ψ T . As far as the nuclear many-body problem is concerned, it is customary to assume that the trial state factorizes into long-and short-range components where F ij and F ijk are two-and three-body correlations, respectively. The symbol S indicates a symmetrized product over nucleon pairs since, in general, the F ij do not commute. VMC calculations explicitly account for the underlying strong alpha-cluster structure of light nuclei. For instance, the totally antisymmetric Jastrow wave function of p-shell nuclei is constructed from a sum over independent-particle terms, Φ A , each having four nucleons in an α-like core and the remaining (A − 4) nucleons in p-shell orbitals [99]: The operator A stands for an antisymmetric sum over all possible A 4 partitions of the A particles into four s-shell and (A − 4) p-shell states. As suggested by standard shell-model studies, the independentparticle wave function |Φ A (LS[n]JJ z T z ) 1234:5...A with the desired JM value of a given nuclear state is obtained using LS coupling, which is most efficient for nuclei with up to A = 12. The symbol [n] is the Young pattern that indicates the spatial symmetry of the angular momentum coupling of the p-shell nucleons [25]. Note that |Φ A (LS[n]JJ z T z ) 1234:5...A is chosen to be independent of the center of mass as it is expressed in terms of the intrinsic coordinates The pair correlation for particles within the s-shell, f ss , arises from the structure of the α particle. The f sp is similar to the f ss at short range, but with a long-range tail that goes to unity at large distances, allowing the wave function to develop a cluster structure. Finally, f pp is set to give the appropriate cluster structure outside the α core. The three-body central correlations, induced by the two-body potential has the following operator independent form where q c 1 and q c 2 are variational parameters. In addition the the scalar correlations of Equation (13), VMC trial wave functions include spin-dependent nuclear correlations, whose operator structure reflects the one of the NN potential of Equation (6) More sophisticated trial wave functions can be constructed by explicitly accounting for spin-orbit correlations, as, for instance, in the cluster variational Monte Carlo calculations of reference [100]. However, the computational cost of these additional terms is significant, while the gain in the variational energy is relatively small [101]. The radial functions u p (r ij ) are generated by minimizing the two-body cluster The variational parameters α p simulate the quenching of spin-isospin interactions between particles i and j due to interactions of these particles with others in the system. The Lagrange multipliers λ p (r ij ) account for short-range screening effects, and are fixed at large distances by the asymptotic behavior of the correlation functions, which is encoded by an additional set of variational parameters. The quality of the trial wave function is improved by reducing the strength of the spin-and isospin-dependent correlation functions u p (r ij ) when a particle k comes close to the pair ij [102] where the three-body operator-dependent correlation induced by the NN interaction is usually expressed as with q p 1 and q p 2 being variational parameters [25]. The three-body correlation operator F ijk turns out to be particularly relevant for when 3N interactions are present in the nuclear Hamiltonian. In this case, its form is suggested by perturbation theory where y q is a scaling parameter, and q a small constant. The superscript q indicates the various terms of the 3N force. It has been shown that the vast majority of the 3N correlations can be recovered by omitting the commutator term C V C ijk , provided that the strength of the anticommutator term A is opportunely adjusted. This allows to save a significant amount of computing time, since anticommutators involving pairs ij and jk can be expressed as generalized tensor operators involving the spins of nucleons i and k only. Hence, the computing time scales as the number of pairs rather than the number of triplets [25].
The expectation values of the form of Equation (11) contain multi-dimensional integrals over all particle positions A deterministic integration of the above integral is computationally prohibitive, therefore Metropolis Monte Carlo techniques are employed to stochastically evaluate it. The order of operators in the symmetrized product of Equation (12), denoted by p and q for the left and right hand side wave functions, respectively, is also sampled. The 3A-dimensional integration is facilitated by introducing a probability distribution, W pq (R), such that

Frontiers
In standard VMC calculations, one usually takes W pq (R) = | Re(Ψ † T,p (R)Ψ T,q (R))|, even though simpler choices might be used to reduce the computational cost. The Metropolis algorithm is used to stochastically sample the probability distribution W pq (R) and obtain a collection of uncorrelated or independent configurations.
Since the nuclear interaction is spin-isospin dependent, the trial state is a sum of complex amplitudes for each spin-isospin state of the system The n s = 2 A many-body spin states can be written as and the isospin ones can be recovered by replacing ↓ with n and ↑ with p. Note that, because of charge conservation, the number of isospin states reduces to n t = A Z . To construct the trial state, one starts from the mean-field component |Φ A (LS[n]JJ z T z ) 1234:5...A . For fixed spatial coordinates R, the spin-isospin independent correlations needed to retrieve |Φ J are simple multiplicative factors, common to all spin amplitudes. The symmetrized product of pair correlation operators is evaluated by successive operations for each pair, sampling their ordering as alluded to earlier. As an example, consider the application of the operator σ 1 · σ 2 on a three-body spin state (for simplicity we neglect the isospin components). Noting that σ i · σ j = 2 P σ ij − 1, where 2 P σ ij exchanges the spin of particles i and j, we obtain: Hence, the many-body spin-isospin basis is closed under the action of the operators contained in the nuclear Hamiltonian.
Most of the computing time is spent on spin-isospin operations like the one just described. They amount to an iterative sequence of large sparse complex matrix multiplications that are performed on-the-fly using explicitly coded subroutines, which mainly rely on three useful matrices. The first matrix m(i, i s ) gives the z-component of the spin of particle i associated to the many-body spin-state i s . A second useful matrix is n exch (k ij , i s ), that provides the number of the many-body spin state obtained by exchanging the spins of particles i and j, belonging to the pair labeled k ij in the state i s . The matrix n flip (i, i s ) yields the number of the spin state obtained by flipping the spin of particle i in the spin state. The action of the operator σ 1 · σ 2 can then be expressed as By utilizing this representation, we only need to evaluate 2 A operations for each pair, instead of the 2 A × 2 A operations that are required using a simple matrix representation in spin space. The tensor operator is slightly more complicated to evaluate and requires both matrices m(i, i s ) and n flip (i, i s ) [103]. Analogous matrices are employed to perform operations on the isospin space, as the two representations are practically identical.
The expectation values of Equation (21) are evaluated by having the operators act entirely on the right hand side of the trial wave function. The matrix machinery used to apply the spin-dependent correlation operators is also used to evaluate O|Ψ T,p . A simple scalar product of this quantity with Ψ T,q |, provides the numerator of the local estimate Ψ † T,q (R)OΨ T,p (R)/W pq (R) and W pq (R) is computed in a similar fashion. The first and second derivatives of the wave function are numerically computed by means of the two-and three-point stencil, respectively. Hence, to determine the kinetic energy, 6A + 1 evaluations of Ψ T (R) are needed. Finally, using the trick described in reference [104], we can evaluate the action of the angular momentum dependent terms in the potential evaluating Ψ T (R) an additional 3A(A − 1)/2 times.
Not only does the size of the wave vector grows exponentially with the number of nucleons, but so does the number of evaluations necessary to calculate the energy, limiting the applicability of the VMC method to A ≤ 12 nuclei. Sampling the spin-isospin state and evaluating the trial wave function's amplitude for that sampled state still requires a number of operations exponential in the particle number, bringing little savings in terms of computing time. Extending VMC calculations to larger nuclear systems requires devising trial wave functions that can capture most of physics of the system while requiring computational time that scales polynomially with A.

Green's function Monte Carlo
Green's function Monte Carlo overcomes the limitations intrinsic to the variational ansatz by using an imaginary-time projection technique to enhance the true ground-state component of a starting trial wave function In the above equation, τ is the imaginary time, and E T is a parameter used to control the normalization. In addition to ground-states properties, excited states can be computed within GFMC. The imaginary-time diffusion yields the lowest-energy eigenstate with the same quantum numbers as |Ψ T . Thus, to obtain an excited state with distinct quantum numbers from the ground state, one only needs to construct a trial wave function with the appropriate quantum numbers. If the excited-state quantum numbers coincide with those of the ground state, more care is needed, but precise results for such states can still be obtained [105].
Except for some specific cases, the direct computation of the propagator e −Hτ for arbitrary values of τ is typically not possible. For small imaginary times δτ = τ /N with N large, the calculation is tractable, and the full propagation to large imaginary times τ can be recovered through the following path integral The GFMC wave function at imaginary time τ + δτ can be written in an integral form where we defined the short-time propagator, or Green's function, Monte Carlo techniques are used to sample the paths by simultaneously evolving a set of configurations -dubbed walkers -in imaginary time, until the distribution converges to the ground-state wave function [106]. To avoid the large statistical errors arising from configurations that diffuse into regions where they make very little contribution to the ground-state wave function, the diffusion process is guided by introducing an importance-sampling function Ψ I (R), which has the same quantum numbers as the ground-state. The importance function is typically taken to coincide with the variational wave function, but different choices are possible. Multiplying Eq. (29) on the left by Ψ † The quantity within squared brackets is the importance-sampled propagator G I δτ (R i+1 , R i ). Note that a set of walkers can be sampled from Ψ † I (R i )Ψ(τ + δτ, R i ) only if this density is positive definite. In this case, the latter can be interpreted as a probability density distribution and its integral determines the size of the population, i.e., the number of walkers. In Fermion systems, however, the positiveness of Ψ † I (R i )Ψ(τ + δτ, R i ) is only granted for exact importance-sampling functions. In general, the nodal surface of the ground state can be different from that of Ψ I . We will return to this point later on. The importance function can be expanded in terms of eigenstates of the Hamiltonian as The Green's function can also be expressed in terms of the same eigenstates: Inserting the last two relations into Equation (29) and integrating over R i+1 , we get If the importance-sampling function closely resembles the ground-state wave function, then c * n δ n0 and E T E 0 , implying Therefore, having accurate importance function reduces the fluctuations in the population size from one time step to the next, thereby reducing the statistical errors in the calculation.
A common approximation for the short-time propagator is based upon the Trotter-Suzuki expansion Here, T is the kinetic energy giving rise to the free-particle propagator that, for non-relativistic systems, can be expressed as a simple Gaussian in configuration space with λ 2 = 4 2 2m δτ . The exponentials of the two-body potentials can be approximated to first order by turning the sums over pairs in the exponent into a symmetrized product of exponentials of the individual pair potentials. The first six terms of the potential can be easily exponentiated, while momentum dependent terms cannot be treated this way. A simple way to include them consists in expanding the exponential of the momentum dependent terms to first order in δτ and use integration by parts to let the derivatives act on the free-particle Green's function. This approach can only be successfully applied to the terms in the potential that are linear in momentum, such as L · S and (L · S) τ i · τ j [107]. On the other hand, contributions to the potential that are quadratic in the momentum cannot be evaluated to first order in this manner. For this reason we use approximations to the full NN potentials, such as the AV8 interaction, that only contain the first eight operators. The difference between AV18 and AV8 is treated in perturbation theory.
More sophisticated alternatives of reducing the time-step error exist and are routinely used in GFMC calculations. The most common one consists in building the Green's function operator as a product of exact two-body propagators where g jk (r jk, i , r jk, i+1 ) is the exact two-body propagator and g 0 jk (r jk, i , r jk, i+1 ) is the two-body freeparticle propagator [108]. At variance with the propagator of Eq. (36), terms quadratic in the angular momentum can in principle be accounted for into the exact pair propagator. However, the inclusion of these terms requires the sampled distribution to have the same locality structure to keep statistical errors under control. Thus, simplified AV8 potentials are also used in the pair propagator, even though in this case no approximations in treating L · S and (L · S) τ i · τ j terms are necessary.
Since the matrix V is the spin/isospin-dependent interaction, the propagator is in turn a matrix in spinisospin space. To deal with it, first a scalar approximation to the importance sampled Green's function, denoted asG I δτ (R i+1 , R i ), is introduced. Recalling the form of the importance sampled Green's function constructing its scalar counterpart requires defining a scalar approximation for the importance-sampling function, which can be taken to beΨ I (R) = Ψ † J (R)Ψ J (R). As for the potential, we can use the average of the central parts in the 1 S 0 and 3 S 1 channels, thus At each time-step, the walkers are propagated with G 0 δτ (R i+1 , R i ) by sampling a 3A-dimensional vector from a gaussian distribution to shift the spatial coordinates. To remove the linear terms coming from the exponential of Eq. (37), we use two mirror points R i+1 = R i ± δR and we consider the corresponding two weights One of the two walkers is kept in the propagation according to a heat-bath sampling among the two normalized weights w ± /( ± w ± ) and the average weight ± w ± /2 is associated to the propagated configuration.
In terms of the scalar Green's function, the propagation of Equation (29) reads Since the new positions are sampled according toG I δτ (R i+1 , R i ), we can conveniently define The imaginary-time evolution of the walker density is given by Iterations of Equation (44) amount to multiple matrix multiplications that are performed using the same matrices used to construct |Ψ T . It has to be stressed that Ψ(τ, R N ) is not the ground-state wave function. It rather represents a spin-isospin set of amplitudes that, when taken in product with the Hermitian conjugate of the importance function, gives an overlap for each component of the wave function. Are the changes in these overlaps that drive the distribution of walkers toward that of the true ground state.
To avoid sign fluctuations in Ψ † I (R i )Ψ(τ, R i ), we sample the walkers from the positive-definite density distribution The first term simply measures the magnitude of the overlap of the wave functions, while the second, with a small coefficient 0.01, ensures a positive definite importance function to allow diffusion across nodal surfaces. This choice for the sampling distribution is monitored by checking how much this estimate of the population size deviates from the actual number of configurations. As alluded to earlier, similarly to standard Fermion diffusion Monte Carlo algorithms, the GFMC method suffers from the Fermion sign problem. Since the configurations are distributed according to I(R i ) defined in Equation (46), the expectation values of observables that commute with the Hamiltonian are estimated as For all other observables, we compute the mixed estimates where the first and the second term correspond to the DMC and VMC expectation value, respectively.
As in standard Fermion diffusion Monte Carlo algorithms, the GFMC method suffers from the Fermion sign problem that arises from stochastically evaluating the matrix elements in Equation (47). The imaginary-time propagator is a local operator, while antisymmetry is a global property of the system. As a consequence, |Ψ(τ ) can contain bosonic components that have much lower energy than the Fermionic ones and are exponentially amplified during the propagation. When the dot product with the antisymmetric Ψ T is taken, the desired Fermionic component is projected out in the expectation values, but the variance -and hence the statistical error -grows exponentially with τ . Because the number of pairs that can be exchanged grows with A, the sign problem also grows exponentially with the number of nucleons. Already for A = 8, the statistical errors grow so fast that convergence cannot be achieved.
To control the sign problem, we adopt the so-called "contrained-path" method [101], originally developed to study condensed matter systems [109]. This method is based on discarding those configurations that in future generations will contribute only noise to expectations values. If we knew the exact ground state, we could discard any walker for which where a sum over spin-isospin states is implied. The sum of these discarded configurations can be written as a state |Ψ d , which has zero overlap with the ground state. Disregarding |Ψ d is justified because it only contains excited-states components and should decay away as τ → ∞. However, in general, the exact ground state is not known, and the constraint is approximately imposed using Ψ T in place of Ψ 0 : The GFMC wave function evolves smoothly in imaginary time and changes can be made arbitrarily small by reducing δτ . Hence, if the wave function is purely scalar, any configuration which yields a negative overlap must first pass through a point at which Ψ T -and hence the overlap -is zero. Discarding these configurations is then sufficient to stabilize the simulation and produce "fixed-node" variational solutions, to the many-Fermion problem. However, the GFMC trial wave function is a vector in spin-isospin space, and there are no coordinates for which all the spin-isospin amplitudes vanish. In addition, the overlap Ψ † T, p (R i )Ψ(τ, R i ) is complex and depends on the particular sampled order p. As a consequence, it does not evolve smoothly and can pass through zero. The constraint of Equation (50) cannot be satisfied for individual configurations, but rather only on average for the sum of discarded configurations. To circumvent these difficulties, we define the overlap We can then introduce a probability for discarding a configuration in terms of the ratio O T, p /I T, p where I T, p corresponds to choosing the ordering p in Ψ I as defined in Equation (46) P The constants α c and β c are adjusted such that the average of the overlap O T, p /I T, P is zero within statistical errors.
In a few cases the constrained propagation converges to the wrong energy (either above or below the correct energy). Therefore, a small number, n u = 10 to 80, of unconstrained steps are made before evaluating expectation values. These few unconstrained steps appear to be sufficient to remove the bias introduced by the constraint but do not greatly increase the statistical error.

Auxiliary field diffusion Monte Carlo
Over the last two decades, the auxiliary field diffusion Monte Carlo method [110] has become a mainstay for studying atomic nuclei [111,89,112,113] and infinite neutron matter [13,87,114]. The AFDMC overcomes the exponential scaling with the number of nucleons of the GFMC by using a spin-isospin basis given by the outer product of single-nucleon spinors where |s i = a i,↑p | ↑ p + a i,↓p | ↓ p + a i,↑n | ↑ n + a i,↓n | ↓ n .
The state vector is fully specified by a set of 4A complex coefficients. As opposed to the many-body spinisospin basis defined in Equation (23), the single-particle one is not closed under the action of two-body operators. To see this, lets apply again the operator σ 1 · σ 2 on a three-body spin state In general, the action of all pairwise spin/isospin operators needed to construct the trial state defined in Equation (12) generates all the 2 A A Z amplitudes of the many-body spin-isospin basis. For this reason, the trial wave function typically used in AFDMC calculations [115,89] is simpler than the one of the GFMC and takes the form where F ij and F ijk are defined in Equations (16) and (20), respectively. Since it contains a linearized version of spin/isospin-dependent two-body correlations, this wave function is significantly cheaper to evaluate than the one used in GFMC, as it scales polynomially with the number of nucleons rather than exponentially. However, because only pairs of nucleons are correlated at a time, the cluster property is violated. Nevertheless, the use of these linearized spin-dependent correlations has enabled a number of remarkably accurate AFDMC calculations, in which properties of atomic nuclei up to A = 16 [111,89,112] have been investigated utilizing the local χEFT interactions of Refs. [12,87]. Very recently, the AFDMC trial wave function has been improved by including quadratic pair correlations [89,116].
The Jastrow component of |Ψ T is also simpler than the one of Equation (13), where the two-body scalar correlation are obtained consistently with the u p (r ij ) minimizing the two-body cluster energy. The three-body scalar correlation is the one defined in Equation (15). The mean-field component is modeled by a sum of Slater determinants, In the above equation we have introduced X = {x 1 , . . . , x A }, where the generalized coordinate x i ≡ {r i , s i } represents both the position R = r 1 , . . . , r A and the spin-isospin coordinates S = s 1 , . . . , s A of the A nucleons. The determinants are coupled with Clebsch-Gordan coefficients C JJz in order to reproduce the total angular momentum, total isospin, and parity. The single-particle orbitals are given by where R nl (r) is the radial function, Y llz is the spherical harmonic, and χ ssz (σ) and χ ttz (τ ) are the complex spinors describing the spin and isospin of the single-particle state.
The AFDMC imaginary-time propagation can be broken up in small time steps similarly to what is done in Equation (28) for the GFMC method. This time however, the generalized coordinate X is used instead of R and the spin-isospin degrees of freedom are also sampled. The AFDMC wave function at imaginary time τ + δτ can be written in an integral form analogous to the one of Equation (29) Using the Trotter decomposition of Equation (36), the short-time Green's function factorizes as Quadratic spin-isospin operators contained in the nuclear potential can connect a single spin-isospin state |S i to all possible |S i+1 states. In order to preserve the single-particle representation, the short-time propagator is linearized utilizing the Hubbard-Stratonovich transformation where x are the auxiliary fields and the operators O are obtained as follows. The first six terms defining the NN potential can be conveniently separated in a spin/isospin-dependent V SD and spin/isospinindependent V SI contributions. To see this in more details, lets consider purely neutron systems, where τ i · τ j = 1, since the extension to isospin-dependent terms is trivial [89]. In this case, V SD can be cast in the form where the operators O n are defined as In the above equations λ n and ψ n iα are the eigenvalues and eigenvectors of the matrix A. Hence, applying the exponential of the spin-dependent terms of the NN interaction amounts to rotating the spin-isospin states of nucleons and the imaginary-time propagation is performed by sampling the auxiliary fieldsx n from the Gaussian probability distribution The spin-orbit term of the NN potential -p = 7 in Equation (6) -is implemented in the propagator as described in reference [117], and appropriate counter terms are included to remove the spurious contributions of order δτ . Presently, the isospin-dependent spin-orbit term of the NN potential, corresponding to p = 8 in Equation (6), cannot be properly treated within the AFDMC algorithm, as its counter term contains cubic spin-isospin operators, preventing the straightforward use of the Hubbard-Stratonovich transformation.
Importance sampling techniques are also routinely implemented in the AFDMC method -in both the spatial coordinates and spin-isospin configurations -to drastically improve the efficiency of the algorithm. To this aim, the propagator of Equation (60) is modified as and we typically take Ψ I (X) = Ψ T (X). At each time step, each walker is propagated sampling a 3A-dimensional vector to shift the spatial coordinates and a set of auxiliary fields X from Gaussian distributions. To remove the linear terms coming from the exponential of both Equations (37) and (64), in analogy to the GFMC method, we consider four weights, corresponding to separately flipping the sign of the spatial moves and spin-isospin rotations In the same spirit as the GFMC algorithm, only one of the four configurations is kept according to a heatbath sampling among the four normalized weights w i /W , with W = 4 i=1 w i /4 being the cumulative weight. The latter is then rescaled by and associated to this new configuration for branching and computing observables. This "plus and minus" procedure, first implemented in the AFDMC method in reference [115] significantly reduces the dependence of the results on δτ .
Expectation values are estimated during the imaginary-time propagation in a similar fashion as for the GFMC To alleviate the sign problem, as in reference [118], we implement an algorithm similar to the constrained-path approximation [119], but applicable to complex wave functions and propagators. The weights w i of Equation (67) are evaluated with and they are set to zero if the ratio is negative. Unlike the fixed-node approximation, which is applicable for scalar potentials and for cases in which a real wave function can be used, the solution obtained from the constrained propagation is not a rigorous upper-bound to the true ground-state energy [101]. To remove the bias associated with this procedure, the configurations obtained from a constrained propagation are further evolved using the following positive-definite importance sampling function [120,89] Ψ I (X) = Re{Ψ T (X)} + α Im{Ψ T (X)} , where we typically take 0.1 < α < 0.5. Along this unconstrained propagation, the expectation value of the energy is estimated according to Equation (69). The asymptotic value is found by fitting the imaginary-time behavior of the unconstrained energy with a single-exponential function, as in reference [25]. Unconstrained propagations have been performed in the latest AFDMC studies of atomic nuclei [111,89] and infinite nucleonic matter [121,116]. An example of unconstrained propagation in 6 Li for the GT+Eτ -1.0 local chiral Hamiltonian is reported in Figure 1, where the blue dots with error bars are the AFDMC unconstrained energies, the red curve is the exponential fit, and the green band represents the final result including the uncertainty coming from the fitting procedure.
In summary, the VMC method is used to find the best possible guess for the wave function for a given nucleus, i.e., it is used to optimize the wave function variational parameters. VMC energies are usually above the ones coming from GFMC and AFDMC calculations, while other observables such as radii and density distributions are in closer agreement. The variationally optimized wave function is then used as input for the (statistically) exact GFMC and AFDMC algorithms. The difference between these two methods relies in their accuracy and limitations. The GFMC method is very accurate in predicting several observables with very small statistical error bars, but its applicability is limited up to 12 nucleons. The AFDMC method can tackle larger systems, but its precision is somewhat reduced and it is currently limited to somewhat simplified interactions [4].

NUCLEAR STRUCTURE RESULTS
GFMC and AFDMC are complimentary methods that have been extensively used in the past to accurately calculate ground-state properties of light nuclei (A 16). In the following we will present results obtained using the GFMC method for ∆-full χEFT potentials, and using the AFDMC method for ∆-less χEFT interactions. In Figure 2 we show the binding energies of nuclei up to 16 O as calculated with GFMC for the NV2+3-Ia potential (red, left) [11], and with AFDMC for the GT+Eτ -1.0 interaction (blue, right) [111,89]. The central green bars are the experimental data. GFMC results only carry Monte Carlo statistical uncertainties, while for AFDMC results, theoretical uncertainties coming from the truncation of the chiral expansion are also included. For an observable X (i) at order i = 0, 2, 3, . . ., the theoretical uncertainty δX (i) is estimated according to the prescription of Epelbaum et al. [74]: where The NV2+3-Ia interaction provides an overall good description of the ground-state energy of light nuclei, including neutron-rich systems with isospin asymmetry as large as 0.6 ( 10 He). This can be appreciated even more by looking at Figure 3, where the ratio between QMC results and experimental data is shown. Above A = 8, the NV2+3-Ia description of binding energies looks slightly less accurate, with some nuclei slightly underbound ( 10 He, 11 B) and some other sightly overbound ( 9 Be, 10 B, 12 C). However, the difference with the experimental values is always less than 0.2 MeV/A, discrepancy that we expect to be fully covered by the uncertainty coming from the truncation of the chiral expansion (i.e., theoretical uncertainty from the interaction model), currently not available for the NV2+3-Ia potential.
The binding energy of very light nuclei is also well reproduced by the GT+Eτ -1.0 interaction, with 8 He slightly underbound (0.37 MeV/A difference compared to the experimental value), but compatible with observations within the estimated statistical plus systematic uncertainties, see Figure 3. Differently from GFMC calculations, AFDMC results for 8 ≤ A ≤ 11 open-shell nuclei are currently not available. The ground-state energy of heavier closed-shell systems, such as 12 C and 16 O, for the GT+Eτ -1.0 potential is higher than the expected result. However, the binding energy of 16 O is still compatible with the experimental value within the fully uncertainty estimate. As discussed in reference [89], the discrepancy found for 12 C is due to the somewhat too simplistic A = 12 AFDMC wave function, that only includes couplings in the p-shells, rather than a deficiency of the interaction itself. It has to be noted that AFDMC results for the GT+Eτ -1.0 interaction carry larger overall uncertainties compared to GFMC results for the NV2+3-Ia potential. This is because the full uncertainty evaluation includes both statistical and theoretical errors. Both QMC methods imply statistical uncertainties of the order of few percent. For the ∆-less potential, the theoretical errors coming from the truncation of the chiral expansion dominate compared to the statistical errors. Considering the next order in the chiral expansion should reduce theoretical uncertainties, and work is currently being done in developing such potentials. Figure 4 shows the charge radii of A ≤ 16 nuclei for the NV2+3-Ia and GT+Eτ -1.0 potentials, with respect to the available experimental data. The expectation value of the charge radius is derived from the point-proton radius r pt using the relation where R 2 p = 0.770(9) fm 2 is the proton radius [127], R 2 n = −0.116(2) fm 2 is the neutron radius [127], (3 2 )/(4M 2 p c 2 ) ≈ 0.033 fm 2 is the Darwin-Foldy correction [128], and r 2 so is a spin-orbit correction due to the anomalous magnetic moment in halo nuclei [129]. The point-nucleon radius r pt is calculated as  Figure 3. Energy ratio between the calculated binding energies and the experimental data. The color scheme is the same as Figure 2.
where r i is the intrinsic coordinate of Eq. (14), N is the number of protons or neutrons, and is the projector operator onto protons (+) or neutrons (−). The charge radius is a mixed expectation value, and it requires the calculation of both VMC and DMC point-proton radii, according to Equation (48). Even though mixed expectation values typically depend on the quality of the employed trial wave functions, for the highly-accurate wave functions employed in the GFMC and AFDMC methods, the extrapolation of the mixed estimate r 2 ch is always small.
Both chiral interactions nicely reproduce the charge radius of helium isotopes. The NV2+3-Ia potential also reproduces the radius of lithium, beryllium, and boron isotopes, with new predictions for 8 Be and 10 Be. The charge radius of 9 Li is underpredicted, whereas that of 12 C is overestimated. The GT+Eτ -1.0 potential works remarkably well in predicting the charge radius of 12 C and 16 O, even though theoretical uncertainties, that dominate over the statistical one, are large. As discussed in the previous paragraphs, going to the next order in the chiral expansion will reduce such theoretical uncertainties. For the GT+Eτ -1.0 interaction, the charge radius of 6 Li turns out to be smaller compared to the experimental value. Once again, this is not a feature of the employed interaction, rather a deficiency of the AFDMC wave function. In fact, differently from GFMC, the current AFDMC wave function does not include dedicated α-deuteron-like correlations, necessary to capture the structural properties of 6 Li.   Figure 2 but for charge radii. Experimental data are from references [123,124,125,126].
In QMC methods, single-nucleon densities are calculated as where P N i is the projector operator of Equation (76) and ρ N integrates to the number of nucleons. In Figures 5, 6 we show the QMC proton density in 12 C and 16 O for the available phenomenological (black) and chiral EFT (blue) potentials. Error bars correspond to statistical uncertainties only. The green bands are the experimental single-nucleon densities, obtained from the "sum-of-Gaussians" parametrization of the charge densities given in reference [132] by unfolding the nucleon form factors and subtracting the small contribution of the neutrons. As can be seen, both phenomenological and chiral EFT interactions provide a good description of the proton density in 12 C. The small discrepancy with the experimental curve at short distance is due to two-body meson exchange currents, not included in the proton density presented here. As shown in reference [130], such currents have little effect on the single-nucleon density for A ≥ 12, slightly reducing its value at small r. The phenomenological AV18+UIX potential underestimates the proton density a short distance in 16  Black triangles are GFMC results for the AV18+IL7 potential [130]. Blue dots are AFDMC results for the GT+Eτ -1.0 interaction [89]. The green band corresponds to the experimental results, unfolded from electron scattering data (see text for details).
three-body potential UIX introduces repulsion in the system, pushing nucleons far away from the nucleus center of mass, and thus resulting in larger radius and smaller central density. The 16 O AFDMC density for the GT+Eτ -1.0 potential is instead in better agreement with the experimental curve.
As opposed to the charge radius, densities are not observables themselves. However, the single-nucleon density can be related to the longitudinal elastic (charge) form factor, physical quantity experimentally accessible via electron-nucleon scattering processes. In fact, the charge form factor can be expressed as the ground-state expectation value of the one-body charge operator [133], which, ignoring small spin-orbit contributions in the one-body current, results in the following expression: whereρ N (q) is the Fourier transform of the single-nucleon density defined in Eq. (77), and Q 2 el = q 2 −ω 2 el is the four-momentum squared, with ω el = q 2 + m 2 A − m A the energy transfer corresponding to the electron scattering elastic peak, m A being the mass of the target nucleus. G N E (Q 2 ) are the nucleon electric form factors, for which we adopt Kelly's parametrization [134].
In Figures 7-9 we show the charge form factor in 6 Li, 12 C, and 16 O. Lines with bands correspond to chiral interactions, solid red for NV2+3-Ia from GFMC calculations and dotted blue for GT+Eτ -1.0 from AFDMC calculations. The black triangles are the results for the phenomenological potentials: AV18+UIX in 6 Li from VMC calculations [135], AV18+IL7 in 12 C from GFMC calculations [130], and AV18+UIX Black triangles are cluster VMC results for the AV18+UIX potential [100]. Blue dots are AFDMC results for the GT+Eτ -1.0 interaction [89].
in 16 O from cluster VMC calculations [100]. Green stars are the available experimental results [131,132,136,137,138]. Note that for all QMC calculations of the charge form factor only one-body charge operators are considered, i.e., no two-body electromagnetic currents are included. However, as shown in references [135,130,139], such operators give a non-negligible contribution only for q > 2 fm −1 , as they basically include relativistic corrections.
In 6 Li all interactions provide a consistent description of the charge form factor, with NV2+3-Ia and AV18+UIX compatible with the experimental results up to q ≈ 2 fm −1 , where two-body currents start playing a role. In the same range, the GT+Eτ -1.0 results are slightly higher, as already indicated by the too small charge radius (see Figure 4). Interestingly, only the phenomenological potential is capable of reproducing the kink in the experimental data, while chiral interactions predict a smooth charge form factor also above q ≈ 3 fm −1 . The inclusion of two-body currents could improve the description of the charge form factor at high momentum. However, this is a momentum range roughly corresponding to the characteristic cut-off of chiral potentials, hence their description of observables in such regime is not supposed to hold. Similar conclusions can be drawn for the charge form factor in 12 C and 16 O, where chiral forces produce results compatible with the experimental data, in particular for the position of the first diffraction peak. This is slightly underestimated for 12 C with the NV2+3-Ia potential, but we expect it to be consistent with the experimental results once the uncertainties coming from the truncation of the chiral expansion are taken into consideration.
Note that the "zero" in the form factor is due to the presence of a term like sin 2 (qR), where R is related to the nucleus charge radius. The zero is obtained when qR = π. Therefore, a smaller (larger) q value for the zero compared to the experimental data suggests a larger (smaller) R value, i.e., a larger (smaller)  Figure 7. Longitudinal elastic form factor in 6 Li for different nuclear potentials. For the NV2+3-Ia (solid red line) and AV18+UIX (black triangles) potentials, errors correspond to statistical Monte Carlo uncertainties. The blue band for the GT+Eτ -1.0 potential also includes the uncertainties coming from the truncation of the chiral expansion. Green stars are the experimental values [131]. Adapted from reference [89].
r ch value. This is indeed verified by QMC calculations. For instance, in Figure 8, the NV2+3-Ia potential predicts a smaller q value for the zero of the charge form factor in 12 C, hence a larger value for the charge radius, as confirmed by Figure 4.

CONCLUSIONS
In this work we have reviewed recent advancements in the development of realistic nuclear interactions and of ab-initio many-body methods for nuclear physics. In particular, we have discussed the recent integration of nearly-local interactions derived within chiral effective field theory, both with and without the inclusion of ∆ degrees of freedom, in quantum Monte Carlo methods, namely variational Monte Carlo, Green's function Monte Carlo, and auxiliary field diffusion Monte Carlo. Such a successful combination lead to accurate and realistic calculations of ground-and excited-state properties of nuclei, that include but is not limited to spectra, charge radii, and longitudinal elastic form factors. Even though the chiral interactions discussed in this work have been constructed using few-body observables only, nucleonnucleon scattering data and properties of nuclei up to A = 5, they provide a remarkable description of the physics of nuclei up to, at least, A = 16, with excellent agreement with experimental data.
The same techniques and nuclear potentials reviewed here have also been used to calculate the equation of state of infinite nuclear and neutron matter [121,116], and to infer properties of neutron stars, with results compatible with astrophysical observations including constraints extracted from gravitational waves of the neutron-star merger GW170817 by the LIGO-Virgo detection [140].  Figure 7 but for 12 C. Experimental data are taken from reference [132]. Adapted from reference [89].
Future efforts will be dedicated to i) further improve the employed local chiral interactions, by extending to higher order in the chiral expansion, ii) calculate electroweak properties in nuclear systems, by consistently deriving electroweak currents, and iii) extend the calculations to heavier nuclei, by improving the AFDMC variational wave functions and the scaling of the algorithm.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.