# Atomic Nuclei From Quantum Monte Carlo Calculations With Chiral EFT Interactions

^{1}Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM, United States^{2}Facility for Rare Isotope Beams, Michigan State University, East Lansing, MI, United States^{3}Physics Division, Argonne National Laboratory, Argonne, IL, United States^{4}Trento Institute of Fundamental Physics and Applications (TIFPA), National Institute for Nuclear Physics (INFN), Trento, Italy^{5}Physics Department, Washington University, St. Louis, MO, United States

Quantum Monte Carlo methods are powerful numerical tools to accurately solve the Schrödinger 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.

## 1. 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 potentials that are consistent with the symmetries of QCD. The solution of nuclear many-body problems requires two main ingredients: an Hamiltonian that accurately models the interactions among the nucleons, and reliable numerical many-body methods to solve the corresponding Schrödinger equation.

Microscopic nuclear Hamiltonians, capable of reproducing nucleon-nucleon scattering data and the properties of few-body systems, have been successfully used to describe light nuclei. For example, the highly-realistic Argonne *v*_{18} two-body potential [3] combined with the phenomenological Illinois-7 three-body force have been employed to predict several properties of nuclei up to *A* = 12 with great accuracy [4]. Several calculations of energies, rms radii, transitions, and densities turn out to be in excellent agreement with experimental data. The main limitation of these phenomenological Hamiltonians is that it is not clear how they can be systematically improved, and how to quantify theoretical, i.e., systematic, uncertainties related to the specific interaction model. Another approach that became very popular in the last two decades consist in deriving nuclear interactions within the framework of chiral Effective Field Theory (χEFT). The advantage of this approach is that it provides the necessary tools to systematically improve the interaction models, to estimate uncertainties related to the truncation of the chiral expansion, and to consistently derive electroweak currents.

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.

## 2. 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 (3*N*) 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 3*N* 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 [3, 17, 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}^{\gamma}$ 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{\sigma}_{i}\xb7{\widehat{\text{r}}}_{ij}{\sigma}_{j}\xb7{\widehat{\text{r}}}_{ij}-{\sigma}_{i}\xb7{\sigma}_{j}$ and *T*_{ij} = 3τ_{iz}τ_{jz} − **τ**_{i} · **τ**_{j} are the tensor and isotensor operators, respectively. The functions, ${v}_{\sigma \tau}^{\pi}(r)$, ${v}_{t\tau}^{\pi}(r)$, ${v}_{\sigma T}^{\pi ,}(r)$, and ${v}_{tT}^{\pi}(r)$ are defined as

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}_{ij}^{I}$, 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}_{ij}^{S}$, 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 $({N}=18)$

where the first eight are given by

with the spin-orbit contribution expressed in terms of the relative angular momentum $\text{L}=\frac{1}{2i}({\text{r}}_{i}-{\text{r}}_{j})\times ({\nabla}_{i}-{\nabla}_{j})$ and the total spin $\text{S}=\frac{1}{2}({\sigma}_{i}+{\sigma}_{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 1,787 *pp* and 2,514 *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}^{\prime}$ (AV8′) contains a charge-independent eight-operator projection, ${{O}}_{ij}^{p=1-8}=[1,{\sigma}_{i}\xb7{\sigma}_{j},{S}_{ij},\text{L}\xb7\text{S}]\otimes [1,{\tau}_{i}\xb7{\tau}_{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 3*N* interaction, is no longer considered realistic.

In addition to *NN* forces, sophisticated phenomenological 3*N* 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-pion-exchange contribution, and a 3*N* contact [4]. More specifically, two families of 3*N* 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–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 ~ 2*m*_{π}.

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–36] for recent review articles). It is important mentioning that alternative power-counting schemes have been also suggested [37–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–48].

Chiral nuclear forces are comprised of both pion-exchange contributions and contact terms. Pion-exchange 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–54, 54–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 3*N* terms are fixed by reproducing properties of light-nuclei. This optimization procedure involves separate fit of the *NN* and 3*N* 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 3*N* 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 3*N* interactions [34, 42, 44, 49, 57–81] and accompanying isospin-symmetry-breaking corrections [82–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 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**.

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 3*N* force is straightforward. In their models, the unknown 3*N* 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 3*N*, force accordingly to Fierz identities [87–89]. In the present work, we are referring to a set of these local chiral interactions, specifically the (*D*2, *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–93]. The contact terms are implemented via a Gaussian representation of the three-dimensional delta function with *R*_{S} as the Gaussian parameter [10, 90, 94]. The pion-range operators are regularized at high-value of momentum transfer via a special radial function characterized by the cutoff *R*_{L} [10, 90, 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 2,700 (3,700) data points with a χ^{2}/datum ≲ 1.1 (≲ 1.4) [10, 90]. In conjunction with these models, two distinct sets of Δ-full 3*N* interactions have also been constructed up to N^{2}LO. In the first, the 3*N* unknown LECs were determined by simultaneously reproducing the experimental trinucleon ground-state energies and neutron-deuteron (*nd*) doublet scattering length for each of the NN 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 labeled 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–98] and references therein.

## 3. 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 spin-isospin 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.

### 3.1. Variational Monte Carlo

The variational Monte Carlo method is routinely used to obtain approximate solutions to the many-body 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 $\left(\begin{array}{c}A\\ 4\end{array}\right)$ partitions of the *A* particles into four *s*-shell and (*A* − 4) *p*-shell states. As suggested by standard shell-model studies, the independent-particle 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}_{1}^{c}$ and ${q}_{2}^{c}$ are variational parameters. In addition 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 energy of the interaction $\stackrel{-}{v}-\lambda $, with

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}_{1}^{p}$ and ${q}_{2}^{p}$ being variational parameters [25]. The three-body correlation operator *F*_{ijk} turns out to be particularly relevant for when 3*N* 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 3*N* force. It has been shown that the vast majority of the 3*N* correlations can be recovered by omitting the commutator term ${\u03f5}_{C}{V}_{ijk}^{C}$, 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 a 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 3*A*-dimensional integration is facilitated by introducing a probability distribution, *W*_{pq}(**R**), such that

In standard VMC calculations, one usually takes ${W}_{pq}(\text{R})=|\text{Re}({\text{\Psi}}_{T,p}^{\u2020}(\text{R}){\text{\Psi}}_{T,q}(\text{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}=\left(\begin{array}{c}A\\ Z\end{array}\right)$. 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 ${\sigma}_{i}\xb7{\sigma}_{j}=2{{P}}_{ij}^{\sigma}-1$, where $2{{P}}_{ij}^{\sigma}$ 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 in 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}|{\text{\Psi}}_{T,p}\rangle $. A simple scalar product of this quantity with 〈Ψ_{T, q}|, provides the numerator of the local estimate ${\text{\Psi}}_{T,q}^{\u2020}(\text{R}){O}{\text{\Psi}}_{T,p}(\text{R})/{W}_{pq}(\text{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, 6*A* + 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 3*A*(*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*.

### 3.2. 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 Equation (29) on the left by ${\text{\Psi}}_{I}^{\u2020}({\text{R}}_{i+1})$ yields

The quantity within squared brackets is the importance-sampled propagator ${G}_{\delta \tau}^{I}({\text{R}}_{i+1},{\text{R}}_{i})$. Note that a set of walkers can be sampled from ${\text{\Psi}}_{I}^{\u2020}({\text{R}}_{i})\text{\Psi}(\tau +\delta \tau ,{\text{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 ${\text{\Psi}}_{I}^{\u2020}({\text{R}}_{i})\text{\Psi}(\tau +\delta \tau ,{\text{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}^{*}\simeq {\delta}_{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 ${\lambda}^{2}=4\frac{{\hslash}^{2}}{2m}\delta \tau $. 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}_{jk}^{0}({\text{r}}_{jk,i},{\text{r}}_{jk,i+1})$ is the two-body free- particle propagator [108]. At variance with the propagator of Equation (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 spin-isospin space. To deal with it, first a scalar approximation to the importance sampled Green's function, denoted as ${\stackrel{~}{G}}_{\delta \tau}^{I}({\text{R}}_{i+1},{\text{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 ${\stackrel{~}{\text{\Psi}}}_{I}(\text{r})=\sqrt{{\text{\Psi}}_{J}^{\u2020}(\text{R}){\text{\Psi}}_{J}(\text{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}_{\delta \tau}^{0}({\text{R}}_{i+1},{\text{R}}_{i})$ by sampling a 3*A*-dimensional vector from a gaussian distribution to shift the spatial coordinates. To remove the linear terms coming from the exponential of Equation (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}_{\pm}/({\sum}_{\pm}{w}_{\pm})$ and the average weight ${\sum}_{\pm}{w}_{\pm}/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 to ${\stackrel{~}{G}}_{\delta \tau}^{I}({\text{R}}_{i+1},{\text{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 ${\text{\Psi}}_{I}^{\u2020}({\text{R}}_{i})\text{\Psi}(\tau ,{\text{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. 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 ${\text{\Psi}}_{T,p}^{\u2020}({\text{R}}_{i})\text{\Psi}(\tau ,{\text{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)

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–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.

### 3.3. 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 [89, 111–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

The state vector is fully specified by a set of 4*A* complex coefficients. As opposed to the many-body spin-isospin 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}\left(\begin{array}{c}A\\ Z\end{array}\right)$ amplitudes of the many-body spin-isospin basis. For this reason, the trial wave function typically used in AFDMC calculations [89, 115] 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 [89, 111, 112] have been investigated utilizing the local χEFT interactions of references [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/isospin-independent *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 ${\psi}_{i\alpha}^{n}$ 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 fields ${\stackrel{-}{x}}_{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 3*A*-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 heat-bath sampling among the four normalized weights *w*_{i}/*W*, with $W={\sum}_{i=1}^{4}{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 [89, 120]

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 [89, 111] and infinite nucleonic matter [116, 121]. 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.

**Figure 1**. Unconstrained evolution in ^{6}Li for the GT+Eτ-1.0 local chiral Hamiltonian. Blue dots with error bars bands are AFDMC energies. The red curve is a single-exponential function fit to the AFDMC results. The green band represents the final energy 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].

## 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) [89, 111]. 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

For the local chiral interaction GT+Eτ-1.0, results are presented at N^{2}LO (*i* = 3) considering *Q* = *m*_{π}/Λ_{b}, with *m*_{π} ≈ 140 MeV and Λ_{b} = 600 MeV [89, 111].

**Figure 2**. Ground-state energies in *A* ≤ 16 nuclei. For each nucleus, experimental results [122] are shown in green at the center. GFMC (AFDMC) results for the NV2+3-Ia [11] (GT+Eτ-1.0 [89]) potential are shown in red (blue) to the left (right) of the experimental values. For the NV2+3-Ia (GT+Eτ-1.0) potential, the colored bands include statistical (statistical plus systematic) uncertainties.

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

**Figure 3**. Energy ratio between the calculated binding energies and the experimental data. The color scheme is the same as Figure 2.

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 $\langle {R}_{p}^{2}\rangle =0.770(9){\mathrm{\text{fm}}}^{2}$ is the proton radius [127], $\langle {R}_{n}^{2}\rangle =-0.116(2){\mathrm{\text{fm}}}^{2}$ is the neutron radius [127], $(3{\hslash}^{2})/(4{M}_{p}^{2}{c}^{2})\approx 0.033{\mathrm{\text{fm}}}^{2}$ is the Darwin-Foldy correction [128], and $\langle {r}_{\text{so}}^{2}\rangle $ is a spin-orbit correction due to the anomalous magnetic moment in halo nuclei [129]. The point-nucleon radius *r*_{pt} is calculated as

where **r**_{i} is the intrinsic coordinate of Equation (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 $\langle {r}_{\text{ch}}^{2}\rangle $ 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.

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}O. As indicated by the cluster VMC analysis of reference [100], the 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.

**Figure 5**. Proton density in ^{12}C. 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).

**Figure 6**. Same as Figure 5 but for ^{16}O. 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].

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 ${\stackrel{~}{\rho}}_{N}(q)$ is the Fourier transform of the single-nucleon density defined in Equation (77), and ${Q}_{\text{el}}^{2}={\mathbf{\text{q}}}^{2}-{\omega}_{\text{el}}^{2}$ is the four-momentum squared, with ${\omega}_{\text{el}}=\sqrt{{q}^{2}+{m}_{A}^{2}}-{m}_{A}$ the energy transfer corresponding to the electron scattering elastic peak, *m*_{A} being the mass of the target nucleus. ${G}_{E}^{N}({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 in ^{16}O from cluster VMC calculations [100]. Green stars are the available experimental results [131, 132, 136–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 [130, 135, 139], such operators give a non-negligible contribution only for *q* > 2 fm^{−1}, as they basically include relativistic corrections.

**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].

**Figure 8**. Same as Figure 7 but for ^{12}C. Experimental data are taken from reference [132]. Adapted from reference [89].

**Figure 9**. Same as Figure 7 but for ^{16}O. Experimental data are from Sick, based on references [136–138]. Adapted from reference [89].

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) *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.

## 5. 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, nucleon-nucleon 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 [116, 121], 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].

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.

## Funding

The work of SG was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under contract No. DE-AC52-06NA25396, by the NUCLEI SciDAC program, by the LDRD program at LANL, and by the DOE Early Career research Program. The work of DL was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under Contract No. DE-SC0013617, and by the NUCLEI SciDAC program. The work of AL was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under contract No. DE-AC02-06CH11357. The work of MP was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under the FRIB Theory Alliance award DE-SC0013617. Computational resources have been provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under Contract No. 89233218CNA000001, and by the National Energy Research Scientific Computing Center (NERSC), which is supported by the U.S. Department of Energy, Office of Science, under contract No. DE-AC02-05CH11231.

## Conflict of Interest

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

## References

1. Winter F, Detmold W, Gambhir AS, Orginos K, Savage MJ, Shanahan PE, et al. First lattice QCD study of the gluonic structure of light nuclei. *Phys Rev D*. (2017) **96**:094512. doi: 10.1103/PhysRevD.96.094512

2. Chang E, Davoudi Z, Detmold W, Gambhir AS, Orginos K, Savage MJ, et al. Scalar, axial, and tensor interactions of light nuclei from lattice QCD. *Phys Rev Lett*. (2018) **120**:152002. doi: 10.1103/PhysRevLett.120.152002

3. Wiringa RB, Stoks VGJ, Schiavilla R. An Accurate nucleon-nucleon potential with charge independence breaking. *Phys Rev C*. (1995) **51**:38–51. doi: 10.1103/PhysRevC.51.38

4. Carlson J, Gandolfi S, Pederiva F, Pieper SC, Schiavilla R, Schmidt KE, et al. Quantum Monte Carlo methods for nuclear physics. *Rev Mod Phys*. (2015) **87**:1067–118. doi: 10.1103/RevModPhys.87.1067

5. Hagen G, Papenbrock T, Hjorth-Jensen M, Dean DJ. Coupled-cluster computations of atomic nuclei. *Rep Progr Phys*. (2014) **77**:096302. doi: 10.1088/0034-4885/77/9/096302

6. Hagen G, Papenbrock T, Ekström A, Wendt KA, Baardsen G, Gandolfi S, et al. Coupled-cluster calculations of nucleonic matter. *Phys Rev C*. (2014) **89**:014319. doi: 10.1103/PhysRevC.89.014319

7. Barrett BR, Navrátil P, Vary JP. *Ab initio* no core shell model. *Prog Part Nucl Phys*. (2013) **69**:131–81. doi: 10.1016/j.ppnp.2012.10.003

8. Bogner S, Furnstahl R, Schwenk A. From low-momentum interactions to nuclear structure. *Prog Part Nucl Phys*. (2010) **65**:94–147. doi: 10.1016/j.ppnp.2010.03.001

9. Dickhoff WH, Barbieri C. Self-consistent Green's function method for nuclei and nuclear matter. *Prog Part Nucl Phys*. (2004) **52**:377–496. doi: 10.1016/j.ppnp.2004.02.038

10. Piarulli M, Girlanda L, Schiavilla R, Navarro Pérez R, Amaro JE, Ruiz Arriola E. Minimally nonlocal nucleon-nucleon potentials with chiral two-pion exchange including Δ resonances. *Phys Rev C*. (2015) **91**:024003. doi: 10.1103/PhysRevC.91.024003

11. Piarulli M, Baroni A, Girlanda L, Kievsky A, Lovato A, Lusk E, et al. Light-nuclei spectra from chiral dynamics. *Phys Rev Lett*. (2018) **120**:052503. doi: 10.1103/PhysRevLett.120.052503

12. Gezerlis A, Tews I, Epelbaum E, Gandolfi S, Hebeler K, Nogga A, et al. Quantum Monte Carlo calculations with chiral effective field theory interactions. *Phys Rev Lett*. (2013) **111**:032501. doi: 10.1103/PhysRevLett.111.032501

13. Tews I, Gandolfi S, Gezerlis A, Schwenk A. Quantum Monte Carlo calculations of neutron matter with chiral three-body forces. *Phys Rev C*. (2016) **93**:024305. doi: 10.1103/PhysRevC.93.024305

14. Rozpedzik D, Golak J, Skibinski R, Witala H, Glockle W, Epelbaum E, et al. A first estimation of chiral four-nucleon force effects in He-4. *Acta Phys Polon B*. (2006) **37**:2889–904.

15. Krüger T, Tews I, Hebeler K, Schwenk A. Neutron matter from chiral effective field theory interactions. *Phys Rev C*. (2013) **88**:025802. doi: 10.1103/PhysRevC.88.025802

16. Yukawa H. On the interaction of elementary particles I. *Proc Phys Math Soc Jap*. (1935) **17**:48–57.

17. Stoks VGJ, Klomp RAM, Terheggen CPF, de Swart JJ. Construction of high quality N N potential models. *Phys Rev C*. (1994) **49**:2950–62. doi: 10.1103/PhysRevC.49.2950

18. Machleidt R. The High precision, charge dependent Bonn nucleon-nucleon potential (CD-Bonn). *Phys Rev C*. (2001) **63**:024001. doi: 10.1103/PhysRevC.63.024001

19. Stoks VGJ, Klomp RAM, Rentmeester MCM, de Swart JJ. Partial wave analaysis of all nucleon-nucleon scattering data below 350-MeV. *Phys Rev C*. (1993) **48**:792–815. doi: 10.1103/PhysRevC.48.792

20. Smith RA, Pandharipande VR. Nucleon-nucleon potentials including pi n delta coupling effects. *Nucl Phys A*. (1976) **256**:327–48. doi: 10.1016/0375-9474(76)90376-6

21. Nelder JA, Mead R. A simplex method for function minimization. *Comput J*. (1965) **7**:308–13. doi: 10.1093/comjnl/7.4.308

22. Gandolfi S, Carlson J, Reddy S, Steiner AW, Wiringa RB. The equation of state of neutron matter, symmetry energy, and neutron star structure. *Eur Phys J A*. (2014) **50**:10. doi: 10.1140/epja/i2014-14010-5

23. Wiringa RB, Pieper SC. Evolution of nuclear spectra with nuclear forces. *Phys Rev Lett*. (2002) **89**:182501. doi: 10.1103/PhysRevLett.89.182501

24. Friar JL, Gibson BF, Payne GL. Recent progress in understanding trinucleon properties. *Ann Rev Nucl Part Sci*. (1984) **34**:403–33. doi: 10.1146/annurev.ns.34.120184.002155

25. Pudliner BS, Pandharipande VR, Carlson J, Pieper SC, Wiringa RB. Quantum Monte Carlo calculations of nuclei with A < = 7. *Phys Rev C*. (1997) **56**:1720–50. doi: 10.1103/PhysRevC.56.1720

26. Navratil P, Vary JP, Barrett BR. Large basis ab initio no-core shell model and its application to C-12. *Phys Rev C*. (2000) **62**:054311. doi: 10.1103/PhysRevC.62.054311

27. Carlson J, Pandharipande VR, Wiringa RB. Three-nucleon interaction in 3-body, 4-body, and infinite-body systems. *Nucl Phys A*. (1983) **401**:59–85. doi: 10.1016/0375-9474(83)90336-6

28. Pieper SC, Pandharipande VR, Wiringa RB, Carlson J. Realistic models of pion exchange three nucleon interactions. *Phys Rev C*. (2001) **64**:014001. doi: 10.1103/PhysRevC.64.014001

29. Maris P, Vary JP, Gandolfi S, Carlson J, Pieper SC. Properties of trapped neutrons interacting with realistic nuclear Hamiltonians. *Phys Rev C*. (2013) **87**:054318. doi: 10.1103/PhysRevC.87.054318

30. Weinberg S. Nuclear forces from chiral Lagrangians. *Phys Lett B*. (1990) **251**:288–92. doi: 10.1016/0370-2693(90)90938-3

31. Weinberg S. Effective chiral Lagrangians for nucleon–pion interactions and nuclear forces. *Nucl Phys B*. (1991) **363**:3–18. doi: 10.1016/0550-3213(91)90231-L

32. Weinberg S. Three body interactions among nucleons and pions. *Phys Lett B*. (1992) **295**:114–21. doi: 10.1016/0370-2693(92)90099-P

33. Epelbaum E, Hammer HW, Meissner UG. Modern theory of nuclear forces. *Rev Mod Phys*. (2009) **81**:1773–825. doi: 10.1103/RevModPhys.81.1773

34. Machleidt R, Entem DR. Chiral effective field theory and nuclear forces. *Phys Rept*. (2011) **503**:1–75. doi: 10.1016/j.physrep.2011.02.001

35. Machleidt R, Sammarruca F. Chiral EFT based nuclear forces: achievements and challenges. *Phys Scripta*. (2016) **91**:083007. doi: 10.1088/0031-8949/91/8/083007

36. Machleidt R. Historical perspective and future prospects for nuclear interactions. *Int J Mod Phys E*. (2017) **26**:1730005. doi: 10.1142/S0218301317300053

37. Kaplan DB, Savage MJ, Wise MB. A new expansion for nucleon-nucleon interactions. *Phys Lett B*. (1998) **424**:390–6. doi: 10.1016/S0370-2693(98)00210-X

38. Kaplan DB, Savage MJ, Wise MB. Two nucleon systems from effective field theory. *Nucl Phys B*. (1998) **534**:329–55. doi: 10.1016/S0550-3213(98)00440-4

39. Nogga A, Timmermans RGE, van Kolck U. Renormalization of one-pion exchange and power counting. *Phys Rev C*. (2005) **72**:054006. doi: 10.1103/PhysRevC.72.054006

40. Pavon Valderrama M, Ruiz Arriola E. Renormalization of NN interaction with chiral two pion exchange potential. central phases and the deuteron. *Phys Rev C*. (2006) **74**:054001. doi: 10.1103/PhysRevC.74.054001

41. Long B, Yang CJ. Renormalizing chiral nuclear forces: triplet channels. *Phys Rev C*. (2012) **85**:034002. doi: 10.1103/PhysRevC.85.034002

42. van Kolck U. Few nucleon forces from chiral Lagrangians. *Phys Rev C*. (1994) **49**:2932–41. doi: 10.1103/PhysRevC.49.2932

43. Furnstahl RJ, Phillips DR, Wesolowski S. A recipe for EFT uncertainty quantification in nuclear physics. *J. Phys. G Nucl. Part. Phys*. (2015) **42**:034028. doi: 10.1088/0954-3899/42/3/034028

44. Epelbaum E, Krebs H, Meissner UG. Improved chiral nucleon-nucleon potential up to next-to-next-to-next-to-leading order. *Eur Phys J A*. (2015) **51**:53. doi: 10.1140/epja/i2015-15053-8

45. Furnstahl RJ, Klco N, Phillips DR, Wesolowski S. Quantifying truncation errors in effective field theory. *Phys Rev C*. (2015) **92**:024005. doi: 10.1103/PhysRevC.92.024005

46. Wesolowski S, Klco N, Furnstahl RJ, Phillips DR, Thapaliya A. Bayesian parameter estimation for effective field theories. *J. Phys. G Nucl. Part. Phys*. (2016) **43**:074001. doi: 10.1088/0954-3899/43/7/074001

47. Melendez JA, Wesolowski S, Furnstahl RJ. Bayesian truncation errors in chiral effective field theory: nucleon-nucleon observables. *Phys Rev C*. (2017) **96**:024003. doi: 10.1103/PhysRevC.96.024003

48. Wesolowski S, Furnstahl RJ, Melendez JA, Phillips DR. Exploring Bayesian parameter estimation for chiral effective field theory using nucleon-nucleon phase shifts. *J Phys G Nucl Part Phys*. (2018) **46**:045102. doi: 10.1088/1361-6471/aaf5fc

49. Krebs H, Epelbaum E, Meissner UG. Nuclear forces with delta-excitations up to next-to-next-to-leading order. I. Peripheral nucleon-nucleon waves. *Eur Phys J A*. (2007) **32**:127–37. doi: 10.1140/epja/i2007-10372-y

50. Ditsche C, Hoferichter M, Kubis B, Meissner UG. Roy-Steiner equations for pion-nucleon scattering. *J. High Energ. Phys*. (2012) **6**:043. doi: 10.1007/JHEP06(2012)043

51. Hoferichter M, Ruiz de Elvira J, Kubis B, Meissner UG. High-precision determination of the pion-nucleon σ term from Roy-Steiner equations. *Phys Rev Lett*. (2015) **115**:092301. doi: 10.1103/PhysRevLett.115.092301

52. Hoferichter M, Ruiz de Elvira J, Kubis B, Meissner UG. Roy-Steiner-equation analysis of pion-nucleon scattering. *Phys Rept*. (2016) **625**:1–88. doi: 10.1016/j.physrep.2016.02.002

53. Hoferichter M, Ruiz de Elvira J, Kubis B, Meissner UG. Matching pion-nucleon Roy-Steiner equations to chiral perturbation theory. *Phys Rev Lett*. (2015) **115**:192301. doi: 10.1103/PhysRevLett.115.192301

54. Siemens D, Bernard V, Epelbaum E, Gasparyan A, Krebs H, Meissner UG. Elastic pion-nucleon scattering in chiral perturbation theory: a fresh look. *Phys Rev C*. (2016) **94**:014620. doi: 10.1103/PhysRevC.94.014620

55. Siemens D, Ruiz de Elvira J, Epelbaum E, Hoferichter M, Krebs H, Kubis B, et al. Reconciling threshold and subthreshold expansions for pion-nucleon scattering. *Phys Lett B*. (2017) **770**:27–34. doi: 10.1016/j.physletb.2017.04.039

56. Yao DL, Siemens D, Bernard V, Epelbaum E, Gasparyan AM, Gegelia J, et al. Pion-nucleon scattering in covariant baryon chiral perturbation theory with explicit Delta resonances. *JHEP*. (2016) **5**:038. doi: 10.1007/JHEP05(2016)038

57. Ekström A, Jansen GR, Wendt KA, Hagen G, Papenbrock T, Carlsson BD, et al. Accurate nuclear radii and binding energies from a chiral interaction. *Phys Rev C*. (2015) **91**:051301. doi: 10.1103/PhysRevC.91.051301

58. Ordonez C, Ray L, van Kolck U. The two nucleon potential from chiral Lagrangians. *Phys Rev C*. (1996) **53**:2086–105. doi: 10.1103/PhysRevC.53.2086

59. Kaiser N, Brockmann R, Weise W. Peripheral nucleon-nucleon phase shifts and chiral symmetry. *Nucl Phys A*. (1997) **625**:758–88. doi: 10.1016/S0375-9474(97)00586-1

60. Kaiser N, Gerstendorfer S, Weise W. Peripheral NN scattering: role of delta excitation, correlated two pion and vector meson exchange. *Nucl Phys A*. (1998) **637**:395–420. doi: 10.1016/S0375-9474(98)00234-6

61. Epelbaum E, Gloeckle W, Meissner UG. Nuclear forces from chiral Lagrangians using the method of unitary transformation. 1. Formalism. *Nucl Phys A*. (1998) **637**:107–34. doi: 10.1016/S0375-9474(98)00220-6

62. Epelbaum E, Gloeckle W, Meissner UG. Nuclear forces from chiral Lagrangians using the method of unitary transformation. 2. The two nucleon system. *Nucl Phys A*. (2000) **671**:295–331. doi: 10.1016/S0375-9474(99)00821-0

63. Kaiser N. Chiral 3 pi exchange N N potentials: Results for representation invariant classes of diagrams. *Phys Rev C*. (2000) **61**:014003. doi: 10.1103/PhysRevC.61.014003

64. Kaiser N. Chiral three pi exchange N N potentials: results for diagrams proportional to *g*(*A*)^{4} and *g*(*A*)^{6}. *Phys Rev C*. (2000) **62**:024001. doi: 10.1103/PhysRevC.62.024001

65. Kaiser N. Chiral 3 pi exchange N N potentials: results for dominant next-to-leading order contributions. *Phys Rev C*. (2001) **63**:044010. doi: 10.1103/PhysRevC.63.044010

66. Kaiser N. Chiral 2 pi exchange N N potentials: two loop contributions. *Phys Rev C*. (2001) **64**:057001. doi: 10.1103/PhysRevC.64.057001

67. Kaiser N. Chiral 2 pi exchange NN potentials: relativistic 1/*M*^{2} corrections. *Phys Rev C*. (2002) **65**:017001. doi: 10.1103/PhysRevC.65.017001

68. Entem DR, Machleidt R. Accurate charge dependent nucleon nucleon potential at fourth order of chiral perturbation theory. *Phys Rev C*. (2003) **68**:041001. doi: 10.1103/PhysRevC.68.041001

69. Epelbaum E, Nogga A, Gloeckle W, Kamada H, Meissner UG, Witala H. Three nucleon forces from chiral effective field theory. *Phys Rev C*. (2002) **66**:064001. doi: 10.1103/PhysRevC.66.064001

70. Navratil P. Local three-nucleon interaction from chiral effective field theory. *Few Body Syst*. (2007) **41**:117–40. doi: 10.1007/s00601-007-0193-3

71. Bernard V, Epelbaum E, Krebs H, Meissner UG. Subleading contributions to the chiral three-nucleon force II: short-range terms and relativistic corrections. *Phys Rev C*. (2011) **84**:054001. doi: 10.1103/PhysRevC.84.054001

72. Girlanda L, Kievsky A, Viviani M. Subleading contributions to the three-nucleon contact interaction. *Phys Rev C*. (2011) **84**:014001. doi: 10.1103/PhysRevC.84.014001

73. Krebs H, Gasparyan A, Epelbaum E. Chiral three-nucleon force at N^{4}LO I: longest-range contributions. *Phys Rev C*. (2012) **85**:054006. doi: 10.1103/PhysRevC.85.054006

74. Epelbaum E, Krebs H, Meissner UG. Precision nucleon-nucleon potential at fifth order in the chiral expansion. *Phys Rev Lett*. (2015) **115**:122301. doi: 10.1103/PhysRevLett.115.122301

75. Entem DR, Kaiser N, Machleidt R, Nosyk Y. Peripheral nucleon-nucleon scattering at fifth order of chiral perturbation theory. *Phys Rev C*. (2015) **91**:014002. doi: 10.1103/PhysRevC.91.014002

76. Entem DR, Kaiser N, Machleidt R, Nosyk Y. Dominant contributions to the nucleon-nucleon interaction at sixth order of chiral perturbation theory. *Phys Rev C*. (2015) **92**:064001. doi: 10.1103/PhysRevC.92.064001

77. Ekström A, Baardsen G, Forssén C, Hagen G, Hjorth-Jensen M, Jansen GR, et al. Optimized chiral nucleon-nucleon interaction at next-to-next-to-leading order. *Phys Rev Lett*. (2013) **110**:192502. doi: 10.1103/PhysRevLett.110.192502

78. Ekström A, Hagen G, Morris TD, Papenbrock T, Schwartz PD. Δ isobars and nuclear saturation. *Phys Rev C*. (2018) **97**:024332. doi: 10.1103/PhysRevC.97.024332

79. Reinert P, Krebs H, Epelbaum E. Semilocal momentum-space regularized chiral two-nucleon potentials up to fifth order. *Eur Phys J A*. (2018) **54**:86. doi: 10.1140/epja/i2018-12516-4

80. Binder S, et al. Few-nucleon and many-nucleon systems with semilocal coordinate-space regularized chiral nucleon-nucleon forces. *Phys Rev C*. (2018) **98**:014002. doi: 10.1103/PhysRevC.98.014002

81. Krebs H, Gasparyan AM, Epelbaum E. Three-nucleon force in chiral EFT with explicit Δ(1232) degrees of freedom: longest-range contributions at fourth order. *Phys Rev C*. (2018) **98**:014003. doi: 10.1103/PhysRevC.98.014003

82. Friar JL, van Kolck U. Charge independence breaking in the two pion exchange nucleon-nucleon force. *Phys Rev C*. (1999) **60**:034006. doi: 10.1103/PhysRevC.60.034006

83. Friar JL, van Kolck U, Rentmeester MCM, Timmermans RGE. The nucleon-mass difference in chiral perturbation theory and nuclear forces. *Phys Rev C*. (2004) **70**:044001. doi: 10.1103/PhysRevC.70.044001

84. Friar JL, Payne GL, van Kolck U. Charge-symmetry-breaking three-nucleon forces. *Phys Rev C*. (2005) **71**:024003. doi: 10.1103/PhysRevC.71.024003

85. Gezerlis A, Tews I, Epelbaum E, Freunek M, Gandolfi S, Hebeler K, et al. Local chiral effective field theory interactions and quantum Monte Carlo applications. *Phys Rev C*. (2014) **90**:054323. doi: 10.1103/PhysRevC.90.054323

86. Fierz M. Zur Fermischen Theorie des β-Zerfalls. *Z Phys*. (1937) **104**:553–65. doi: 10.1007/BF01330070

87. Lynn JE, Tews I, Carlson J, Gandolfi S, Gezerlis A, Schmidt KE, et al. Chiral three-nucleon interactions in light nuclei, neutron-α scattering, and neutron matter. *Phys Rev Lett*. (2016) **116**:062501. doi: 10.1103/PhysRevLett.116.062501

88. Lynn JE, Tews I, Carlson J, Gandolfi S, Gezerlis A, Schmidt KE, et al. Quantum Monte Carlo calculations of light nuclei with local chiral two- and three-nucleon interactions. *Phys Rev C*. (2017) **96**:054007. doi: 10.1103/PhysRevC.96.054007

89. Lonardoni D, Gandolfi S, Lynn JE, Petrie C, Carlson J, Schmidt KE, et al. Auxiliary field diffusion Monte Carlo calculations of light and medium-mass nuclei with local chiral interactions. *Phys Rev C*. (2018) **97**:044318. doi: 10.1103/PhysRevC.97.044318

90. Piarulli M, Girlanda L, Schiavilla R, Kievsky A, Lovato A, Marcucci LE, et al. Local chiral potentials with Δ-intermediate states and the structure of light nuclei. *Phys Rev C*. (2016) **94**:054007. doi: 10.1103/PhysRevC.94.054007

91. Navarro Pérez R, Amaro JE, Ruiz Arriola E. Coarse-grained potential analysis of neutron-proton and proton-proton scattering below the pion production threshold. *Phys Rev C*. (2013) **88**:064002. doi: 10.1103/PhysRevC.88.064002

92. Navarro Pérez R, Amaro JE, Ruiz Arriola E. Coarse grained NN potential with chiral two pion exchange. *Phys Rev C*. (2014) **89**:024004. doi: 10.1103/PhysRevC.89.024004

93. Navarro Pérez R, Amaro JE, Ruiz Arriola E. Statistical error analysis for phenomenological nucleon-nucleon potentials. *Phys Rev C*. (2014) **89**:064006. doi: 10.1103/PhysRevC.89.064006

94. Baroni A, et al. Local chiral interactions, the tritium Gamow-Teller matrix element, and the three-nucleon contact term. *Phys Rev C*. (2018) **98**:044003. doi: 10.1103/PhysRevC.98.044003

95. Baroni A, Girlanda L, Kievsky A, Marcucci LE, Schiavilla R, Viviani M. Tritium β-decay in chiral effective field theory. *Phys Rev C*. (2016) **94**:024003. doi: 10.1103/PhysRevC.94.024003

96. Bacca S, Pastore S. Electromagnetic reactions on light nuclei. *J. Phys. G Nucl. Part. Phys*. (2014) **41**:123002. doi: 10.1088/0954-3899/41/12/123002

97. Marcucci LE, Gross F, Pena MT, Piarulli M, Schiavilla R, Sick I, et al. Electromagnetic structure of few-nucleon ground states. *J. Phys. G Nucl. Part. Phys*. (2016) **43**:023002. doi: 10.1088/0954-3899/43/2/023002

98. Lynn JE, Tews I, Gandolfi S, Lovato A. Quantum Monte Carlo methods in nuclear physics: recent advances. *Annu Rev*. (2019) **69**:279–305. doi: 10.1146/annurev-nucl-101918-023600

99. Pieper SC, Varga K, Wiringa RB. Quantum Monte Carlo calculations of A = 9, A = 10 nuclei. *Phys Rev C*. (2002) **66**:044310. doi: 10.1103/PhysRevC.66.044310

100. Lonardoni D, Lovato A, Pieper SC, Wiringa RB. Variational calculation of the ground state of closed-shell nuclei up to *A* = 40. *Phys Rev C*. (2017) **96**:024326. doi: 10.1103/PhysRevC.96.024326

101. Wiringa RB, Pieper SC, Carlson J, Pandharipande VR. Quantum Monte Carlo calculations of A = 8 nuclei. *Phys Rev C*. (2000) **62**:014001. doi: 10.1103/PhysRevC.62.014001

102. Lomnitz-Adler J, Pandharipande VR, Smith RA. Monte Carlo calculations of triton and 4 He nuclei with the Reid potential. *Nucl Phys A*. (1981) **361**:399–411. doi: 10.1016/0375-9474(81)90642-4

103. Pudliner BS. *Green's Function Monte Carlo Calculations of Few-Nucleon Systems*. Champaign, IL: University of Illinois at Urbana-Champaign (1996).

104. Schiavilla R, Pandharipande VR, Wiringa RB. Momentum distributions in a A = 3 and 4 nuclei. *Nucl Phys A*. (1986) **449**:219–42. doi: 10.1016/0375-9474(86)90003-5

105. Carlsson BD, Ekström A, Forssen C, Stromberg DF, Jansen GR, Lilja O, et al. Uncertainty analysis and order-by-order optimization of chiral nuclear interactions. *Phys Rev X*. (2016) **6**:011019. doi: 10.1103/PhysRevX.6.011019

106. Foulkes WMC, Mitas L, Needs RJ, Rajagopal G. Quantum Monte Carlo simulations of solids. *Rev Mod Phys*. (2001) **73**:33–83. doi: 10.1103/RevModPhys.73.33

107. Carlson J. Alpha particle structure. *Phys Rev C*. (1988) **38**:1879–85. doi: 10.1103/PhysRevC.38.1879

108. Schmidt KE, Lee MA. High-accuracy Trotter-formula method for path integrals. *Phys Rev E*. (1995) **51**:5495–8. doi: 10.1103/PhysRevE.51.5495

109. Zhang S, Carlson J, Gubernatis JE. Constrained Path Quantum Monte Carlo method for fermion ground states. *Phys Rev Lett*. (1995) **74**:3652–5. doi: 10.1103/PhysRevLett.74.3652

110. Schmidt KE, Fantoni S. A quantum Monte Carlo method for nucleon systems. *Phys Lett B*. (1999) **446**:99–103. doi: 10.1016/S0370-2693(98)01522-6

111. Lonardoni D, Carlson J, Gandolfi S, Lynn JE, Schmidt KE, Schwenk A, et al. Properties of nuclei up to *A* = 16 using local chiral interactions. *Phys Rev Lett*. (2018) **120**:122502. doi: 10.1103/PhysRevLett.120.122502

112. Lonardoni D, Gandolfi S, Wang XB, Carlson J. Single- and two-nucleon momentum distributions for local chiral interactions. *Phys Rev C*. (2018) **98**:014322. doi: 10.1103/PhysRevC.98.014322

113. Lynn JE, Lonardoni D, Carlson J, Chen JW, Detmold W, Gandolfi S, et al. *Ab initio* short-range-correlation scaling factors from light to medium-mass nuclei. *J Phys G Nucl Part Phys*. (2020) **47**:045109. doi: 10.1088/1361-6471/ab6af7

114. Tews I, Carlson J, Gandolfi S, Reddy S. Constraining the speed of sound inside neutron stars with chiral effective field theory interactions and observations. *Astrophys J*. (2018) **860**:149. doi: 10.3847/1538-4357/aac267

115. Gandolfi S, Lovato A, Carlson J, Schmidt KE. From the lightest nuclei to the equation of state of asymmetric nuclear matter with realistic nuclear interactions. *Phys Rev C*. (2014) **90**:061306. doi: 10.1103/PhysRevC.90.061306

116. Lonardoni D, Tews I, Gandolfi S, Carlson J. Nuclear and neutron-star matter from local chiral interactions. *arXiv*. (2019) 1912.09411 nucl-th.

117. Sarsa A, Fantoni S, Schmidt KE, Pederiva F. Neutron matter at zero temperature with auxiliary field diffusion Monte Carlo. *Phys Rev C*. (2003) **68**:024308. doi: 10.1103/PhysRevC.68.024308

118. Zhang S, Krakauer H. Quantum Monte Carlo method using phase-free random walks with Slater determinants. *Phys Rev Lett*. (2003) **90**:136401. doi: 10.1103/PhysRevLett.90.136401

119. Zhang S, Carlson J, Gubernatis JE. A Constrained path Monte Carlo method for fermion ground states. *Phys Rev B*. (1997) **55**:7464. doi: 10.1103/PhysRevB.55.7464

120. Pederiva F, Sarsa A, Schmidt KE, Fantoni S. Auxiliary field diffusion Monte Carlo calculation of ground state properties of neutron drops. *Nucl Phys A*. (2004) **742**:255–68. doi: 10.1016/j.nuclphysa.2004.06.030

121. Piarulli M, Bombaci I, Logoteta D, Lovato A, Wiringa RB. Benchmark calculations of pure neutron matter with realistic nucleon-nucleon interactions. *arXiv*. (2019) 1908.04426 nucl-th.

122. Wang M, Audi G, Kondev FG, Huang WJ, Naimi S, Xu X. The AME2016 atomic mass evaluation (II). Tables, graphs and references. *Chinese Phys C*. (2017) **41**:030003. doi: 10.1088/1674-1137/41/3/030003

123. Fricke G, Heilig K. Nuclear charge radii.1 In: Schopper H, editor. *Elementary Particles, Nuclei and Atoms.* Springer (2004).

124. Sick I. Precise radii of light nuclei from electron scattering. *Lect Notes Phys*. (2008) **745**:57–77. doi: 10.1007/978-3-540-75479-4_4

125. Mueller P, Sulai IA, Villari ACC, Alcántara-Núñez JA, Alves-Condé R, Bailey K, et al. Nuclear charge radius of ^{8}He. *Phys Rev Lett*. (2007) **99**:252501. doi: 10.1103/PhysRevLett.99.252501

126. Nörtershäuser W, Neff T, Sánchez R, Sick I. Charge radii and ground state structure of lithium isotopes: experiment and theory reexamined. *Phys Rev C*. (2011) **84**:024307. doi: 10.1103/PhysRevC.84.024307

127. Beringer J, Arguin JF, Barnett RM, Copic K, Dahl O, Groom DE, et al. Review of particle physics. *Phys Rev D*. (2012) **86**:010001. doi: 10.1103/PhysRevD.86.010001

128. Friar JL, Martorell J, Sprung DWL. Nuclear sizes and the isotope shift. *Phys Rev A*. (1997) **56**:4579–86. doi: 10.1103/PhysRevA.56.4579

129. Ong A, Berengut JC, Flambaum VV. Effect of spin-orbit nuclear charge density corrections due to the anomalous magnetic moment on halonuclei. *Phys Rev C*. (2010) **82**:014320. doi: 10.1103/PhysRevC.82.014320

130. Lovato A, Gandolfi S, Butler R, Carlson J, Lusk E, Pieper SC, et al. Charge form factor and sum rules of electromagnetic response functions in ^{1}2C. *Phys Rev Lett*. (2013) **111**:092501. doi: 10.1103/PhysRevLett.111.092501

131. Li GC, Sick I, Whitney RR, Yearian MR. High-energy electron scattering from ^{6}Li. *Nucl Phys A*. (1971) **162**:583–92. doi: 10.1016/0375-9474(71)90257-0

132. de Vries H, de Jager CW, de Vries C. Nuclear charge-density-distribution parameters from elastic electron scattering. *At Data Nucl Data Tables*. (1987) **36**:495. doi: 10.1016/0092-640X(87)90013-1

133. McVoy KW, Van Hove L. Inelastic electron-nucleus scattering and nucleon-nucleon correlations. *Phys Rev*. (1962) **125**:1034–43. doi: 10.1103/PhysRev.125.1034

134. Kelly JJ. Simple parametrization of nucleon form factors. *Phys Rev C*. (2004) **70**:068202. doi: 10.1103/PhysRevC.70.068202

135. Wiringa RB, Schiavilla R. Microscopic calculation of ^{6}Li elastic and transition form factors. *Phys Rev Lett*. (1998) **81**:4317–20. doi: 10.1103/PhysRevLett.81.4317

136. Sick I, McCarthy JS. Elastic electron scattering from ^{1}2C and ^{1}6O. *Nucl Phys A*. (1970) **150**:631–54. doi: 10.1016/0375-9474(70)90423-9

137. Schütz W. Elastische Elektronenstreuung an14N,15N,16O und18O bei kleiner Impulsübertragung. *Z Phys A*. (1975) **273**:69–75. doi: 10.1007/BF01435760

139. Mihaila B, Heisenberg JH. Microscopic calculation of the inclusive electron scattering structure function in ^{1}6O. *Phys Rev Lett*. (2000) **84**:1403–6. doi: 10.1103/PhysRevLett.84.1403

Keywords: quantum Monte Carlo methods, variational Monte Carlo, Green's function Monte Carlo, auxiliary field diffusion Monte Carlo, chiral effective field theory, nuclear Hamiltonians, nuclear structure

Citation: Gandolfi S, Lonardoni D, Lovato A and Piarulli M (2020) Atomic Nuclei From Quantum Monte Carlo Calculations With Chiral EFT Interactions. *Front. Phys.* 8:117. doi: 10.3389/fphy.2020.00117

Received: 06 January 2020; Accepted: 25 March 2020;

Published: 24 April 2020.

Edited by:

Carlo Barbieri, University of Surrey, United KingdomReviewed by:

Andrea Celentano, Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Genova, ItalyMarco La Cognata, Laboratori Nazionali del Sud (INFN), Catania, Italy

Copyright © 2020 Gandolfi, Lonardoni, Lovato and Piarulli. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Stefano Gandolfi, stefano@lanl.gov