Atomic Nuclei From Quantum Monte Carlo Calculations With Chiral EFT Interactions

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₁₈ 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

$\begin{array}{ll}H={\displaystyle \sum _i}{T}_i+{\displaystyle \sum _{i<j}}{v}_{ij}+{\displaystyle \sum _{i<j<k}}{V}_{ijk}+\cdots ,& (1)\end{array}$

where T_i is the one-body kinetic energy operator, v_ij is the nucleon-nucleon (NN) interaction between particles i and j, V_ijk is the three-nucleon (3N) interaction between particles i, j, and k, and the ellipsis indicate interactions involving more than three particles. There are indications that four-nucleon interactions may contribute at the level of only ~ 100 keV in ⁴He [14] or pure neutron matter [15], and therefore are negligible compared to NN and 3N components. Hence, current formulations of the microscopic model do not typically include them (see, for example, reference [4]).

The NN interaction term in the nuclear Hamiltonian is the most studied of all, with thousands of experimental data points at laboratory energies (E_lab) from essentially zero to hundreds of MeV. It consists of a long-range component, for inter-nucleon separation r ≳ 2 fm, due to one-pion exchange (OPE) [16], and intermediate- and short-range components, for, respectively, 1 fm ≲ r ≲ 2 fm and r ≲ 1 fm, derived, up to the mid 1990's, almost exclusively from meson-exchange phenomenology [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 χ²/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₁₈ (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:

$\begin{array}{ll}{v}_{ij}=v_{ij}^{\gamma }+v_{ij}^{\pi }+v_{ij}^I+v_{ij}^S.& (2)\end{array}$

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

$\begin{array}{ll}{v}_{ij}^{\pi }=\left[v_{\sigma \tau }^{\pi }(r){\sigma }_i·{\sigma }_j+v_{t\tau }^{\pi }(r)S_{ij}\right]{\tau }_i·{\tau }_j+\left[v_{\sigma T}^{\pi }(r){\sigma }_i·{\sigma }_j+v_{tT}^{\pi }(r)S_{ij}\right]T_{ij},& (3)\end{array}$

where σ adn τ are the Pauli matrices that operate over the spin and isospin of particles, and $S_{ij}=3{\sigma }_i·{\widehat{\text{r}}}_{ij}{\sigma }_j·{\widehat{\text{r}}}_{ij}-{\sigma }_i·{\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

$\begin{array}{ll}{v}_{\sigma \tau }^{\pi }(r)=\frac{Y_0(r)+2Y_{+}(r)}{3},& v_{t\tau }^{\pi }(r)=\frac{T_0(r)+2T_{+}(r)}{3},\\ v_{\sigma T}^{\pi }(r)=\frac{Y_0(q)-Y_{+}(r)}{3},& v_{tT}^{\pi }(r)=\frac{T_0(r)-T_{+}(r)}{3},& (4)\end{array}$

where Y_α(r) and T_α(r) are the Yukawa and tensor functions given by

$\begin{array}{ll}{Y}_{\alpha }(r)=\frac{g_A^2}{12\pi }\frac{m_{{\pi }_{\alpha }}^3}{{(2f_{\pi })}^2}\frac{e^{-x_{\alpha }}}{x_{\alpha }},T_{\alpha }(r)=Y_{\alpha }(r)\left(1+\frac{3}{x_{\alpha }}+\frac{3}{x_{\alpha }^2}\right),& (5)\end{array}$

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)$

$\begin{array}{ll}{v}_{ij}={\displaystyle \sum _{p=1}^{{N}}}{v}^p(r_{ij}){{O}}_{ij}^p,& (6)\end{array}$

where the first eight are given by

$\begin{array}{ll}{{O}}_{ij}^{p=1-8}=[1,{\sigma }_i·{\sigma }_j,S_{ij},\text{L}·\text{S}]\otimes [1,{\tau }_i·{\tau }_j],& (7)\end{array}$

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

$\begin{array}{ll}{{O}}_{ij}^{p=9-14}=[{\text{L}}^2,{\text{L}}^2{\sigma }_i·{\sigma }_j,{(\text{L}·\text{S})}^2]\otimes [1,{\tau }_i·{\tau }_j],& (8)\end{array}$

while the p = 15 − 18 are charge-independence breaking terms

$\begin{array}{ll}{{O}}_{ij}^{p=15-18}=[T_{ij},T_{ij}{\sigma }_i·{\sigma }_j,T_{ij}{S}_{ij},{\tau }_{i,z}+{\tau }_{j,z}].& (9)\end{array}$

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 χ²/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·{\sigma }_j,S_{ij},\text{L}·\text{S}]\otimes [1,{\tau }_i·{\tau }_j]$, of the full NN potential, constructed to preserve the potential in all S and P waves as well as the ³D₁ and its coupling to the ³S₁, 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, ³H and ³He [24]. In the late 1990s and early 2000s this realization was also extended to the spectra (ground and low-lying excited states) of light p-shell nuclei, for instance, in calculations based on QMC methods [25] and in no-core shell-model studies [26]. Consequently, the microscopic model with only NN interactions fit to scattering data, without the inclusion of a 3N interaction, is no longer considered realistic.

In addition to NN forces, sophisticated phenomenological 3N interactions have been then developed. They are generally expressed as a sum of a TPE P-wave term, a TPE S-wave contribution, a three-pion-exchange contribution, and a 3N contact [4]. More specifically, two families of 3N interactions were obtained in combination with the AV18 potential: the Urbana IX (UIX) [27] and Illinois 7 (IL7) [28] models. The UIX potential contains two parameters fit to reproduce the ground-state energy of ³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 ~ 2m_π.

By its own nature, the χEFT formulation has an expansion in powers of pion momenta as its organizing principle. Most chiral interactions employed in recent nuclear structure and reaction calculations are based on Weinberg power counting. Within Weinberg power counting, the interactions are expanded in powers of the typical momentum p over the breakdown scale Λ_b ~ GeV, Q = p/Λ_b, where the breakdown scale denotes momenta at which the short distance structure becomes important and cannot be neglected and absorbed into contact interactions anymore (see references [33–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²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 3N terms are fixed by reproducing properties of light-nuclei. This optimization procedure involves separate fit of the NN and 3N terms. Recently, a different strategy has been introduced by Ekström et al. [57]. This new approach is based on a simultaneous fit of the NN and 3N forces to low-energy NN data, deuteron binding energy, and binding energies and charge radii of hydrogen, helium, carbon, and oxygen isotopes.

Within χEFT, many studies have been carried out dealing with the construction and optimization of NN and 3N interactions [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²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 χ² 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₀ in the range of R₀ = 1.0 − 1.2 fm. The motivations why the authors of references [12, 85] truncated the chiral expansion of these local potentials at N²LO is because at this order it is (i) possible to have a fully local representation of the NN chiral interactions and (ii) the inclusion of consistent 3N force is straightforward. In their models, the unknown 3N LECs are obtained by reproducing the binding energy of ⁴He as well as the P-wave n − α elastic scattering phase shifts. In addition, they explore different parametrization for the 3N, force accordingly to Fierz identities [87–89]. In the present work, we are referring to a set of these local chiral interactions, specifically the (D2, Eτ) model with R₀ = 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²LO (Δ-full χEFT), ii) retaining contact terms up to N³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 χ²/datum ≲ 1.1 (≲ 1.4) [10, 90]. In conjunction with these models, two distinct sets of Δ-full 3N interactions have also been constructed up to N²LO. In the first, the 3N unknown LECs were determined by simultaneously reproducing the experimental trinucleon ground-state energies and neutron-deuteron (nd) doublet scattering length for each of the 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

$\begin{array}{ll}H|{\text{Ψ}}_n\rangle =E_n|{\text{Ψ}}_n\rangle .& (10)\end{array}$

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

$\begin{array}{ll}\begin{array}{l}\frac{\langle {\Psi }_T|H|{\Psi }_T\rangle }{\langle {\Psi }_T|{\Psi }_T\rangle }=E_T\ge E_0\hfill \end{array}& (11)\end{array}$

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

$\begin{array}{ll}|{\text{Ψ}}_T\rangle =(1-{\displaystyle \sum _{i<j<k}}{F}_{ijk})\left({S}{\displaystyle \prod _{i<j}}{F}_{ij}\right)|{\Phi }_J\rangle ,& (12)\end{array}$

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]:

$\begin{array}{l}|{\Phi }_J\rangle ={A}[{\displaystyle \prod _{i<j<k}{f}_{ijk}^c}{\displaystyle \prod _{i<j\le 4}{f}_{ss}}(r_{ij}){\displaystyle \prod _{k\le 4<l\le A}{f}_{sp}}(r_{kl})\\ \times {\displaystyle \sum _{LS[n]}(}{\beta }_{LS[n]}{\displaystyle \prod _{4<l<m\le A}{f}_{pp}^{[n]}}(r_{lm})|{\Phi }_A{(LS[n]J{J}_z{T}_z)}_{1234:5\dots A}\rangle )].& (13)\end{array}$

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_zT_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_zT_z)_1234:5…A〉 is chosen to be independent of the center of mass as it is expressed in terms of the intrinsic coordinates

$\begin{array}{ll}{\text{r}}_i\to {\text{r}}_i-{\text{R}}_{\text{CM}},{\text{R}}_{\text{CM}}=\frac{1}{A}{\displaystyle \sum _{i=1}^A}{\text{r}}_i.& (14)\end{array}$

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

$\begin{array}{ll}{f}_{ijk}^c=1-q_1^c({\text{r}}_{ij}·{\text{r}}_{ik})({\text{r}}_{ij}·{\text{r}}_{jk})({\text{r}}_{ik}·{\text{r}}_{jk})e^{-q_2^c(r_{ij}+r_{ik}+r_{jk})},& (15)\end{array}$

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)

$\begin{array}{ll}{F}_{ij}=(1+U_{ij})=(1+{\displaystyle \sum _{p=2}^6}{u}^p(r_{ij}){{O}}_{ij}^p).& (16)\end{array}$

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

$\begin{array}{ll}\stackrel{-}{v}-\lambda ={\displaystyle \sum _{p=1}^{18}}\left({\alpha }_p{v}^p(r_{ij}){{O}}_{ij}^p-{\lambda }_p(r_{ij})\right).& (17)\end{array}$

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]

$\begin{array}{ll}{u}^p(r_{ij})\to \left[{\displaystyle \prod _{k\ne i\ne j}}{f}_{ijk}^p({\text{r}}_{ij},{\text{r}}_{ik})\right]u^p(r_{ij}),& (18)\end{array}$

where the three-body operator-dependent correlation induced by the NN interaction is usually expressed as

$\begin{array}{ll}{f}_{ijk}^p({\text{r}}_{ij},{\text{r}}_{ik})=1-q_1^p(1-{\widehat{\text{r}}}_{ik}·{\widehat{\text{r}}}_{jk})e^{-q_2^p(r_{ij}+r_{ik}+r_{jk})},& (19)\end{array}$

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 3N interactions are present in the nuclear Hamiltonian. In this case, its form is suggested by perturbation theory

$\begin{array}{ll}{F}_{ijk}={\displaystyle \sum _q}{ϵ}_q{V}_{ijk}^q(y_q{r}_{ij},y_q{r}_{ik},y_q{r}_{jk}),& (20)\end{array}$

where y_q is a scaling parameter, and ϵ_q a small, constant. The superscript q indicates the various terms of the 3N force. It has been shown that the vast majority of the 3N correlations can be recovered by omitting the commutator term ${ϵ}_C{V}_{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

$\begin{array}{ll}\langle {O}\rangle =\frac{\int d\text{R}{\text{Ψ}}_T^{†}(\text{R}){O}{\text{Ψ}}_T(\text{R})}{\int d\text{R}{\text{Ψ}}_T^{†}(\text{R}){\text{Ψ}}_T(\text{R})}.& (21)\end{array}$

A deterministic integration of the above integral is computationally prohibitive, therefore Metropolis Monte Carlo techniques are employed to stochastically evaluate it. The order of operators in the symmetrized product of Equation (12), denoted by p and q for the left and right hand side wave functions, respectively, is also sampled. The 3A-dimensional integration is facilitated by introducing a probability distribution, W_pq(R), such that

$\begin{array}{ll}\langle {O}\rangle =\frac{{\displaystyle {\sum }_{p,q}}\int d\text{R}\frac{{\text{Ψ}}_{T,p}^{†}(\text{R}){O}{\text{Ψ}}_{T,q}(\text{R})}{W_{pq}(\text{R})}{W}_{pq}(\text{R})}{{\displaystyle {\sum }_{p,q}}\int d\text{R}\frac{{\text{Ψ}}_{T,p}^{†}(\text{R}){\text{Ψ}}_{T,q}(\text{R})}{W_{pq}(\text{R})}{W}_{pq}(\text{R})}.& (22)\end{array}$

In standard VMC calculations, one usually takes $W_{pq}(\text{R})=|\text{Re}({\text{Ψ}}_{T,p}^{†}(\text{R}){\text{Ψ}}_{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

$\begin{array}{ll}|{\text{Ψ}}_T\rangle ={\displaystyle \sum _{i_s\le n_s,i_t\le n_t}}a(i_s,i_t;\text{R})|{\chi }_{i_s}{\chi }_{i_t}\rangle .& (23)\end{array}$

The $n_s=2^A$ many-body spin states can be written as

$\begin{array}{l}|{\chi }_1\rangle =|{\downarrow }_1,{\downarrow }_2,\dots ,{\downarrow }_A\rangle \\ |{\chi }_2\rangle =|{\uparrow }_1,{\downarrow }_2,\dots ,{\downarrow }_A\rangle \\ |{\chi }_3\rangle =|{\downarrow }_1,{\uparrow }_2,\dots ,{\downarrow }_A\rangle \\ \dots \\ |{\chi }_{n_s}\rangle =|{\uparrow }_1,{\uparrow }_2,\dots ,{\uparrow }_A\rangle & (24)\end{array}$

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_zT_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 σ₁ · σ₂ on a three-body spin state (for simplicity we neglect the isospin components). Noting that ${\sigma }_i·{\sigma }_j=2{{P}}_{ij}^{\sigma }-1$, where $2{{P}}_{ij}^{\sigma }$ exchanges the spin of particles i and j, we obtain:

$\begin{array}{ll}{\sigma }_1·{\sigma }_2\left(\begin{array}{c}{a}_{\uparrow \uparrow \uparrow }\\ a_{\uparrow \uparrow \downarrow }\\ a_{\uparrow \downarrow \uparrow }\\ a_{\uparrow \downarrow \downarrow }\\ a_{\downarrow \uparrow \uparrow }\\ a_{\downarrow \uparrow \downarrow }\\ a_{\downarrow \downarrow \uparrow }\\ a_{\downarrow \downarrow \downarrow }\end{array}\right)=\left(\begin{array}{c}{a}_{\uparrow \uparrow \uparrow }\\ a_{\uparrow \uparrow \downarrow }\\ 2a_{\downarrow \uparrow \uparrow }-a_{\uparrow \downarrow \uparrow }\\ 2a_{\downarrow \uparrow \downarrow }-a_{\uparrow \downarrow \downarrow }\\ 2a_{\uparrow \downarrow \uparrow }-a_{\downarrow \uparrow \uparrow }\\ 2a_{\uparrow \downarrow \downarrow }-a_{\downarrow \uparrow \downarrow }\\ a_{\downarrow \downarrow \uparrow }\\ a_{\downarrow \downarrow \downarrow }\end{array}\right).& (25)\end{array}$

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 σ₁ · σ₂ can then be expressed as

$\begin{array}{l}{\sigma }_1·{\sigma }_2{\displaystyle \sum _{i_s,i_t}}a(i_s,i_t;\text{R})|{\chi }_{i_s}{\chi }_{i_t}\rangle \\ ={\displaystyle \sum _{i_s,i_t}}[2a(i_s,i_t;\text{R})-a(n_{\text{exch}}(k_{ij},i_s),i_t;\text{R})]|{\chi }_{i_s}{\chi }_{i_t}\rangle .& (26)\end{array}$

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{Ψ}}_{T,p}\rangle$. A simple scalar product of this quantity with 〈Ψ_{T, q}|, provides the numerator of the local estimate ${\text{Ψ}}_{T,q}^{†}(\text{R}){O}{\text{Ψ}}_{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, 6A + 1 evaluations of Ψ_T(R) are needed. Finally, using the trick described in reference [104], we can evaluate the action of the angular momentum dependent terms in the potential evaluating Ψ_T(R) an additional 3A(A − 1)/2 times.

Not only does the size of the wave vector grows exponentially with the number of nucleons, but so does the number of evaluations necessary to calculate the energy, limiting the applicability of the VMC method to A ≤ 12 nuclei. Sampling the spin-isospin state and evaluating the trial wave function's amplitude for that sampled state still requires a number of operations exponential in the particle number, bringing little savings in terms of computing time. Extending VMC calculations to larger nuclear systems requires devising trial wave functions that can capture most of physics of the system while requiring computational time that scales polynomially with A.

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

$\begin{array}{ll}|{\text{Ψ}}_0\rangle \propto {\displaystyle \underset{\tau \to \infty }{lim}}{e}^{-(H-E_T)\tau }|{\text{Ψ}}_T\rangle .& (27)\end{array}$

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

$\begin{array}{l}\text{Ψ}(\tau ,{\text{R}}_N)={\displaystyle \int }{\displaystyle \prod _{i=0}^{N-1}}\text{d}{\text{R}}_i\langle {\text{R}}_N\left|e^{-(H-E_T)\delta \tau }\right|{\text{R}}_{N-1}\rangle \cdots \\ \langle {\text{R}}_1\left|e^{-(H-E_T)\delta \tau }\right|{\text{R}}_0\rangle \langle {\text{R}}_0\rangle {\text{Ψ}}_T.& (28)\end{array}$

The GFMC wave function at imaginary time τ+δτ can be written in an integral form

$\begin{array}{ll}\text{Ψ}(\tau +\delta \tau ,{\text{R}}_{i+1})={\displaystyle \int }d{\text{R}}_i{G}_{\delta \tau }({\text{R}}_{i+1},{\text{R}}_i)\text{Ψ}(\tau ,{\text{R}}_i),& (29)\end{array}$

where we defined the short-time propagator, or Green's function,

$\begin{array}{ll}{G}_{\delta \tau }({\text{R}}_{i+1},{\text{R}}_i)=\langle {\text{R}}_{i+1}|e^{-H\delta \tau }|{\text{R}}_i\rangle .& (30)\end{array}$

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{Ψ}}_I^{†}({\text{R}}_{i+1})$ yields

$\begin{array}{l}{\text{Ψ}}_I^{†}({\text{R}}_{i+1})\text{Ψ}(\tau +\delta \tau ,{\text{R}}_{i+1})\\ ={\displaystyle \int }d{\text{R}}_i\left[{\text{Ψ}}_I^{†}({\text{R}}_{i+1})G_{\delta \tau }({\text{R}}_{i+1},{\text{R}}_i)\frac{1}{{\text{Ψ}}_I^{†}({\text{R}}_i)}\right]{\text{Ψ}}_I^{†}({\text{R}}_i)\text{Ψ}(\tau ,{\text{R}}_i).& (31)\end{array}$

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{Ψ}}_I^{†}({\text{R}}_i)\text{Ψ}(\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{Ψ}}_I^{†}({\text{R}}_i)\text{Ψ}(\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

$\begin{array}{ll}{\text{Ψ}}_I({\text{R}}_i)={\displaystyle \sum _n}{c}_n{\text{Ψ}}_n(\text{R}).& (32)\end{array}$

The Green's function can also be expressed in terms of the same eigenstates:

$\begin{array}{ll}{G}_{\delta \tau }({\text{R}}_{i+1},{\text{R}}_i)={\displaystyle \sum _n}{\text{Ψ}}_n({\text{R}}_{i+1})e^{-(E_n-E_T)\delta \tau }{\text{Ψ}}_n^{†}(\text{R}).& (33)\end{array}$

Inserting the last two relations into Equation (29) and integrating over R_i+1, we get

$\begin{array}{l}{\displaystyle \sum _n}{c}_n^{*}{\displaystyle \int }d{\text{R}}_{i+1}{\text{Ψ}}_n^{†}({\text{R}}_{i+1})\text{Ψ}(\tau +\delta \tau ,{\text{R}}_{i+1})\\ ={\displaystyle \sum _n}{c}_n^{*}{\displaystyle \int }d{\text{R}}_i{\text{Ψ}}_n^{†}({\text{R}}_i)e^{-(E_n-E_T)\delta \tau }\text{Ψ}(\tau ,{\text{R}}_i).& (34)\end{array}$

If the importance-sampling function closely resembles the ground-state wave function, then $c_n^{*}\simeq {\delta }_{n0}$ and E_T ≃ E₀, implying

$\begin{array}{ll}{\displaystyle \int }d{\text{R}}_{i+1}{\text{Ψ}}_0^{†}({\text{R}}_{i+1})\text{Ψ}(\tau +\delta \tau ,{\text{R}}_{i+1})\simeq {\displaystyle \int }d{\text{R}}_i{\text{Ψ}}_0^{†}({\text{R}}_i)\text{Ψ}(\tau ,{\text{R}}_i).& (35)\end{array}$

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

$\begin{array}{ll}{G}_{\delta \tau }({\text{R}}_{i+1},{\text{R}}_i)=e^{-V({\text{R}}_{i+1})\delta \tau /2}\langle {\text{R}}_{i+1}|e^{-T\delta \tau }|{\text{R}}_i\rangle e^{-V({\text{R}}_i)\delta \tau /2}+o(\delta {\tau }^3).& (36)\end{array}$

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

$\begin{array}{ll}\langle {\text{R}}_{i+1}|e^{-T\delta \tau }|{\text{R}}_i\rangle =G_{\delta \tau }^0({\text{R}}_{i+1},{\text{R}}_i)={\left[\frac{1}{{\lambda }^3{\pi }^{3/2}}\right]}^A{e}^{-{({\text{R}}_{i+1}-{\text{R}}_i)}^2/{\lambda }^2},& (37)\end{array}$

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

$\begin{array}{ll}{G}_{\delta \tau }({\text{R}}_{i+1},{\text{R}}_i)=\left({S}{\displaystyle \prod _{j<k}}\frac{g_{jk}({\text{r}}_{jk,i},{\text{r}}_{jk,i+1})}{g_{jk}^0({\text{r}}_{jk,i},{\text{r}}_{jk,i+1})}\right)G_{\delta \tau }^0({\text{R}}_{i+1},{\text{R}}_i),& (38)\end{array}$

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 ${\tilde{G}}_{\delta \tau }^I({\text{R}}_{i+1},{\text{R}}_i)$, is introduced. Recalling the form of the importance sampled Green's function

$\begin{array}{ll}{G}_{\delta \tau }^I({\text{R}}_{i+1},{\text{R}}_i)=\frac{{\text{Ψ}}_I^{†}({\text{R}}_{i+1})}{{\text{Ψ}}_I^{†}({\text{R}}_i)}{G}_{\delta \tau }({\text{R}}_{i+1},{\text{R}}_i),& (39)\end{array}$

constructing its scalar counterpart requires defining a scalar approximation for the importance-sampling function, which can be taken to be ${\widetilde{\text{Ψ}}}_I(\text{r})=\sqrt{{\text{Ψ}}_J^{†}(\text{R}){\text{Ψ}}_J(\text{R})}$. As for the potential, we can use the average of the central parts in the ¹S₀ and ³S₁ channels, thus

$\begin{array}{ll}\begin{array}{l}{\tilde{G}}_{\delta \tau }^I(R_{i+1},R_i)=\frac{{\tilde{\Psi }}_I(R_{i+1})}{{\tilde{\Psi }}_I(R_i)}{e}^{-[V_{10}(R_{i+1})+V_{01}(R_{i+1})]\delta \tau /4}{G}_{\delta \tau }^0\\ (R_{i+1},R_i)e^{-[V_{10}(R_i)+V_{01}(R_i]\delta \tau /4}.\end{array}& (40)\end{array}$

At each time-step, the walkers are propagated with $G_{\delta \tau }^0({\text{R}}_{i+1},{\text{R}}_i)$ by sampling a 3A-dimensional vector from a gaussian distribution to shift the spatial coordinates. To remove the linear terms coming from the exponential of Equation (37), we use two mirror points R_i+1 = R_i ± δR and we consider the corresponding two weights

$\begin{array}{ll}{w}_{\pm }=\frac{{\tilde{\Psi }}_I(R_i\pm \delta R)}{{\tilde{\Psi }}_I(R{R}_i)}{e}^{-[V_{10}(R_i\pm \delta R)+V_{01}(R_i\pm \delta R)+V_{10}(R_i)+V_{01}(R_i]\delta \tau /4}{e}^{E_T\delta \tau }.& (41)\end{array}$

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

$\begin{array}{ll}\text{Ψ}(\tau +\delta \tau ,{\text{R}}_{i+1})={\displaystyle \int }d{\text{R}}_i\left[\frac{G_{\delta \tau }({\text{R}}_{i+1},{\text{R}}_i)}{{\tilde{G}}_{\delta \tau }^I({\text{R}}_{i+1},{\text{R}}_i)}\right]{\tilde{G}}_{\delta \tau }^I({\text{R}}_{i+1},{\text{R}}_i)\text{Ψ}(\tau ,{\text{R}}_i).& (42)\end{array}$

Since the new positions are sampled according to ${\tilde{G}}_{\delta \tau }^I({\text{R}}_{i+1},{\text{R}}_i)$, we can conveniently define

$\begin{array}{ll}\text{Ψ}(\tau +\delta \tau ,{\text{R}}_{i+1})=\frac{G_{\delta \tau }({\text{R}}_{i+1},{\text{R}}_i)}{{\tilde{G}}_{\delta \tau }^I({\text{R}}_{i+1},{\text{R}}_i)}\text{Ψ}(\tau ,{\text{R}}_i).& (43)\end{array}$

The imaginary-time evolution of the walker density is given by

$\begin{array}{l}{\text{Ψ}}_I^{†}({\text{R}}_{i+1})\text{Ψ}(\tau +\delta \tau ,{\text{R}}_{i+1})=\\ {\displaystyle \int }d{\text{R}}_i\left[\frac{{\text{Ψ}}_I^{†}({\text{R}}_{i+1})G_{\delta \tau }({\text{R}}_{i+1},{\text{R}}_i)\text{Ψ}(\tau ,{\text{R}}_i)}{{\text{Ψ}}_I^{†}({\text{R}}_i){\tilde{G}}_{\delta \tau }^I({\text{R}}_{i+1},{\text{R}}_i)\text{Ψ}(\tau ,{\text{R}}_i)}\right]{\tilde{G}}_{\delta \tau }^I\\ ({\text{R}}_{i+1},{\text{R}}_i){\text{Ψ}}_I^{†}({\text{R}}_i)\text{Ψ}(\tau ,{\text{R}}_i).& (44)\end{array}$

Iterations of Equation (44) amount to multiple matrix multiplications

$\begin{array}{l}\text{Ψ}(\tau ,{\text{R}}_N)=\left[\frac{G_{\delta \tau }({\text{R}}_N,{\text{R}}_{N-1})}{{\tilde{G}}_{\delta \tau }^I({\text{R}}_N,{\text{R}}_{N-1})}\right]\left[\frac{G_{\delta \tau }({\text{R}}_{N-1},{\text{R}}_{N-2})}{{\tilde{G}}_{\delta \tau }^I({\text{R}}_{N-1},{\text{R}}_{N-2})}\right]\\ \cdots \left[\frac{G_{\delta \tau }({\text{R}}_1,{\text{R}}_0)}{{\tilde{G}}_{\delta \tau }^I({\text{R}}_1,{\text{R}}_0)}\right]{\text{Ψ}}_T({\text{R}}_0),& (45)\end{array}$

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{Ψ}}_I^{†}({\text{R}}_i)\text{Ψ}(\tau ,{\text{R}}_i)$, we sample the walkers from the positive-definite density distribution

$\begin{array}{l}I({\text{R}}_i)=\left|{\displaystyle \sum _{i_s,i_t}}\langle {\text{Ψ}}_I({\text{R}}_i)|{\chi }_{i_s}{\chi }_{i_t}\rangle \langle {\chi }_{i_s}{\chi }_{i_t}|\text{Ψ}(\tau ,{\text{R}}_i)\rangle \right|\\ +ϵ{\displaystyle \sum _{i_s,i_t}}\left|\langle {\text{Ψ}}_I({\text{R}}_i)|{\chi }_{i_s}{\chi }_{i_t}\rangle \langle {\chi }_{i_s}{\chi }_{i_t}|\text{Ψ}(\tau ,{\text{R}}_i)\rangle \right|.& (46)\end{array}$

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

$\begin{array}{ll}\langle {O}(\tau )\rangle =\frac{\langle {\text{Ψ}}_T|{O}|\text{Ψ}(\tau )\rangle }{\langle {\text{Ψ}}_T|\text{Ψ}(\tau )\rangle }=\frac{{\displaystyle {\sum }_{R_i}}\langle {\text{Ψ}}_T({\text{R}}_i)|O|\text{Ψ}(\tau ,{\text{R}}_i)\rangle /I({\text{R}}_i)}{{\displaystyle {\sum }_{R_i}}\langle {\text{Ψ}}_T({\text{R}}_i)|\text{Ψ}(\tau ,{\text{R}}_i)\rangle /I({\text{R}}_i)}.& (47)\end{array}$

For all other observables, we compute the mixed estimates

$\begin{array}{ll}\langle {O}(\tau )\rangle \simeq 2\frac{\langle {\text{Ψ}}_T|{O}|\text{Ψ}(\tau )\rangle }{\langle {\text{Ψ}}_T|\text{Ψ}(\tau )\rangle }-\frac{\langle {\text{Ψ}}_T|{O}|{\text{Ψ}}_T\rangle }{\langle {\text{Ψ}}_T|{\text{Ψ}}_T\rangle },& (48)\end{array}$

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

$\begin{array}{ll}{\text{Ψ}}_0^{†}({\text{R}}_i)\text{Ψ}(\tau ,{\text{R}}_i)=0,& (49)\end{array}$

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 Ψ₀:

$\begin{array}{ll}\langle {\text{Ψ}}_T|{\text{Ψ}}_d\rangle =0.& (50)\end{array}$

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{Ψ}}_{T,p}^{†}({\text{R}}_i)\text{Ψ}(\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

$\begin{array}{ll}{O}_{T,p}=\Re \left[{\text{Ψ}}_{T,p}^{†}({\text{R}}_i)\text{Ψ}(\tau ,{\text{R}}_i)\right].& (51)\end{array}$

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)

$\begin{array}{l}P\left[{\text{Ψ}}_{T,p}^{†}({\text{R}}_i),\text{Ψ}(\tau ,{\text{R}}_i)\right]=\{\begin{array}{cc}0& O/I>{\alpha }_c\\ \frac{{\alpha }_C-O/I}{{\alpha }_c-{\beta }_c}& {\alpha }_c>O/I>{\beta }_c\\ 1& O/I<{\beta }_c\end{array}\end{array}$

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

$\begin{array}{ll}|{\chi }_{i_s}{\chi }_{i_t}\rangle \to |S\rangle \equiv |s_1\rangle \otimes |s_2\rangle \otimes \cdots \otimes |s_A\rangle ,& (52)\end{array}$

where

$\begin{array}{ll}|s_i\rangle =a_{i,\uparrow p}|\uparrow p\rangle +a_{i,\downarrow p}|\downarrow p\rangle +a_{i,\uparrow n}|\uparrow n\rangle +a_{i,\downarrow n}|\downarrow n\rangle .& (53)\end{array}$

The state vector is fully specified by a set of 4A 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 σ₁ · σ₂ on a three-body spin state

$\begin{array}{ll}\begin{array}{c}{\sigma }_1·{\sigma }_2[\left(a_{1,\uparrow }|\uparrow \rangle +a_{1,\downarrow }|\downarrow \rangle \right)\otimes \left(a_{2,\uparrow }|\uparrow \rangle +a_{2,\downarrow }|\downarrow \rangle \right)\otimes \hfill \\ \left(a_{3,\uparrow }|\uparrow \rangle +a_{3,\downarrow }|\downarrow \rangle \right)]=2[\left(a_{2,\uparrow }|\uparrow \rangle +a_{2,\downarrow }|\downarrow \rangle \right)\otimes \hfill \\ \left(a_{1,\uparrow }|\uparrow \rangle +a_{1,\downarrow }|\downarrow \rangle \right)\otimes \left(a_{3,\uparrow }|\uparrow \rangle +a_{3,\downarrow }|\downarrow \rangle \right)]\hfill \\ -[\left(a_{1,\uparrow }|\uparrow \rangle +a_{1,\downarrow }|\downarrow \rangle \right)\otimes \left(a_{2,\uparrow }|\uparrow \rangle +a_{2,\downarrow }|\downarrow \rangle \right)\otimes \hfill \\ \left(a_{3,\uparrow }|\uparrow \rangle +a_{3,\downarrow }|\downarrow \rangle \right)].\hfill \end{array}\hfill & (54)\end{array}$

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

$\begin{array}{ll}|{\Psi }_T\rangle =\left(1-{\displaystyle \sum _{i<j}{F}_{ij}}-{\displaystyle \sum _{i<j<k}{F}_{ijk}}\right)|{\Phi }_J\rangle ,& (55)\end{array}$

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),

$\begin{array}{ll}\begin{array}{l}|{\Phi }_J\rangle ={\displaystyle \prod _{i<j}{f}_{ij}^c}{\displaystyle \prod _{i<j<k}{f}_{ijk}^c}|{\Phi }_A(J^{\pi },J_z,T_z)\rangle ,\hfill \end{array}& (56)\end{array}$

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,

$\begin{array}{ll}\langle X|\Phi (J^{\pi },J_z,T_z)\rangle ={\displaystyle \sum _n}{c}_n{\left[{\displaystyle \sum _{J{J}_z}}{C}_{J{J}_z}{A}[{\varphi }_{{\alpha }_1}(x_1)\dots {\varphi }_{{\alpha }_A}(x_A)]\right]}_{J{J}_z}.& (57)\end{array}$

In the above equation we have introduced X = {x₁, …, x_A}, where the generalized coordinate x_i≡{r_i, s_i} represents both the position R = r₁, …, r_A and the spin-isospin coordinates S = s₁, …, 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

$\begin{array}{ll}{\varphi }_{\alpha }(x_i)=R_{nl}(r_i)Y_{l{l}_z}({\widehat{r}}_i){\chi }_{s{s}_z}(\sigma ){\chi }_{t{t}_z}(\tau ),& (58)\end{array}$

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)

$\begin{array}{ll}\text{Ψ}(\tau +\delta \tau ,X_{i+1})={\displaystyle \sum _{S_i}}{\displaystyle \int }d{\text{R}}_i{G}_{\delta \tau }(X_{i+1},X_i)\text{Ψ}(\tau ,X_i).& (59)\end{array}$

Using the Trotter decomposition of Equation (36), the short-time Green's function factorizes as

$\begin{array}{l}{G}_{\delta \tau }(X_{i+1},X_i)=G_{\delta \tau }^0({\text{R}}_{i+1},{\text{R}}_i)\langle S_{i+1}|e^{-(V({\text{R}}_{i+1})/2+V({\text{R}}_i)/2-E_T)\delta \tau }|S_i\rangle \\ +o(\delta {\tau }^3).& (60)\end{array}$

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

$\begin{array}{ll}{e}^{-\lambda {{O}}^2\delta \tau /2}=\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }}dx{e}^{-x^2/2}{e}^{x\sqrt{-\lambda \delta \tau }{O}},& (61)\end{array}$

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

$\begin{array}{ll}{V}_{SD}=\frac{1}{2}{\displaystyle \sum _{i\alpha j\beta }}{A}_{i\alpha ,j\beta }{\sigma }_i^{\alpha }{\sigma }_j^{\beta }=\frac{1}{2}{\displaystyle \sum _{n=1}^{3A}}{{O}}_n^2{\lambda }_n,& (62)\end{array}$

where the operators ${{O}}_n$ are defined as

$\begin{array}{ll}{{O}}_n={\displaystyle \sum _{i,\alpha }}{\sigma }_i^{\alpha }{\psi }_{i\alpha }^n.& (63)\end{array}$

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

$\begin{array}{ll}{e}^{-V(R_i)\delta \tau /2}|S_i\rangle ={\displaystyle \prod _n\frac{1}{\sqrt{2\pi }}}{\displaystyle \int d}{x}_n{e}^{-x_n^2/2}{e}^{x_n\sqrt{-\lambda \delta \tau }{O}_n}|S_i\rangle ,& (64)\end{array}$

and the imaginary-time propagation is performed by sampling the auxiliary fields ${\stackrel{-}{x}}_n$ from the Gaussian probability distribution

$\begin{array}{ll}|S_{i+1}\rangle ={\displaystyle \prod _n{e}^{{\overline{x}}_n\sqrt{-\lambda \delta \tau }{O}_n}}|S_i\rangle .& (65)\end{array}$

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

$\begin{array}{ll}{G}_{\delta \tau }^I(X_{i+1},X_i)=G_{\delta \tau }(X_{i+1},X_i)\frac{{\text{Ψ}}_I(X_{i+1})}{{\text{Ψ}}_I(X_i)},& (66)\end{array}$

and we typically take Ψ_I(X) = Ψ_T(X). At each time step, each walker is propagated sampling a 3A-dimensional vector to shift the spatial coordinates and a set of auxiliary fields ${X}$ from Gaussian distributions. To remove the linear terms coming from the exponential of both Equations (37) and (64), in analogy to the GFMC method, we consider four weights, corresponding to separately flipping the sign of the spatial moves and spin-isospin rotations

$\begin{array}{ll}{w}_i=\frac{{\text{Ψ}}_I(\pm {\text{R}}_{i+1},S_{i+1}(\pm {X}))}{{\text{Ψ}}_I({\text{R}}_i,S_i)}.& (67)\end{array}$

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

$\begin{array}{ll}W\to W{e}^{-\left[V_{SI}({\text{R}}_i)/2+V_{SI}({\text{R}}_{i+1})/2-E_T\right]\delta \tau },& (68)\end{array}$

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

$\begin{array}{ll}\langle {O}(\tau )\rangle =\frac{\langle {\text{Ψ}}_T|{O}|\text{Ψ}(\tau )\rangle }{\langle {\text{Ψ}}_T|\text{Ψ}(\tau )\rangle }=\frac{{\displaystyle {\sum }_{X_i}}\langle {\text{Ψ}}_T(X_i)|{O}|\text{Ψ}(\tau ,X_i)\rangle /{\text{Ψ}}_I(X_i)}{{\displaystyle {\sum }_{X_i}}\langle {\text{Ψ}}_T(X_i)|\text{Ψ}(\tau ,X_i)\rangle /{\text{Ψ}}_I(X_i)},& (69)\end{array}$

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

$\begin{array}{ll}\frac{{\text{Ψ}}_I(X_{i+1})}{{\text{Ψ}}_I(X_i)}\to \text{Re}\left\{\frac{{\text{Ψ}}_I(X_{i+1})}{{\text{Ψ}}_I(X_i)}\right\},& (70)\end{array}$

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]

$\begin{array}{ll}{\text{Ψ}}_I(X)=\left|\text{Re}\left\{{\text{Ψ}}_T(X)\right\}\right|+\alpha \left|\text{Im}\left\{{\text{Ψ}}_T(X)\right\}\right|,& (71)\end{array}$

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

Figure 1. Unconstrained evolution in ⁶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 ¹⁶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⁽ⁱ⁾ at order i = 0, 2, 3, …, the theoretical uncertainty δX⁽ⁱ⁾ is estimated according to the prescription of Epelbaum et al. [74]:

$\begin{array}{l}\delta X^{(0)}=Q^2\left|X^{(0)}\right|,\\ \delta X^{(i)}={\displaystyle \underset{2\le j\le i}{max}}\left(Q^{i+1}\left|X^{(0)}\right|,Q^{i+1-j}|\Delta X^{(j)}|\right)\text{ for }i\ge 2,\\ \delta X^{(i)}\ge max(|X^{(j\ge i)}-X^{(k\ge i)}|),& (72)\end{array}$

where

$\begin{array}{l}\Delta X^{(2)}\equiv X^{(2)}-X^{(0)},\\ \Delta X^{(i)}\equiv X^{(i)}-X^{(i-1)}\text{ for }i\ge 3.& (73)\end{array}$

For the local chiral interaction GT+Eτ-1.0, results are presented at N²LO (i = 3) considering Q = m_π/Λ_b, with m_π ≈ 140 MeV and Λ_b = 600 MeV [89, 111].

FIGURE 2

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 (¹⁰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 (¹⁰He, ¹¹B) and some other sightly overbound (⁹Be, ¹⁰B, ¹²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

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 ⁸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 ¹²C and ¹⁶O, for the GT+Eτ-1.0 potential is higher than the expected result. However, the binding energy of ¹⁶O is still compatible with the experimental value within the fully uncertainty estimate. As discussed in reference [89], the discrepancy found for ¹²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

$\begin{array}{ll}\langle r_{\text{ch}}^2\rangle =\langle r_{\text{pt}}^2\rangle +\langle R_p^2\rangle +\frac{A-Z}{Z}\langle R_n^2\rangle +\frac{3{\hslash }^2}{4M_p^2{c}^2}+\langle r_{\text{so}}^2\rangle ,& (74)\end{array}$

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)/(4M_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

$\begin{array}{ll}\langle r_N^2\rangle =\frac{1}{{N}}\langle \text{Ψ}|{\displaystyle \sum _i}{{P}}_{N_i}|{\text{r}}_i{|}^2|\text{Ψ}\rangle ,& (75)\end{array}$

FIGURE 4

Figure 4. Same as Figure 2 but for charge radii. Experimental data are from references [123–126].

where r_i is the intrinsic coordinate of Equation (14), ${N}$ is the number of protons or neutrons, and

$\begin{array}{ll}{{P}}_{N_i}=\frac{1\pm {\tau }_{z_i}}{2}& (76)\end{array}$

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 ⁸Be and ¹⁰Be. The charge radius of ⁹Li is underpredicted, whereas that of ¹²C is overestimated. The GT+Eτ-1.0 potential works remarkably well in predicting the charge radius of ¹²C and ¹⁶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 ⁶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 ⁶Li.

In QMC methods, single-nucleon densities are calculated as

$\begin{array}{lll}{\rho }_N(r)& =\frac{1}{4\pi r^2}\langle \text{Ψ}\left|{\displaystyle \sum _i}{{P}}_{N_i}\delta (r-|{\text{r}}_i|)\right|\text{Ψ}\rangle ,& (77)\end{array}$

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 ¹²C and ¹⁶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 ¹²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 ¹⁶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 ¹⁶O AFDMC density for the GT+Eτ-1.0 potential is instead in better agreement with the experimental curve.

FIGURE 5

Figure 5. Proton density in ¹²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

Figure 6. Same as Figure 5 but for ¹⁶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:

$\begin{array}{ll}{F}_L(q)=\frac{1}{Z}\frac{G_E^p(Q_{\text{el}}^2){\tilde{\rho }}_p(q)+G_E^n(Q_{\text{el}}^2){\tilde{\rho }}_n(q)}{\sqrt{1+Q_{\text{el}}^2/(4m_N^2)}},& (78)\end{array}$

where ${\tilde{\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 ⁶Li, ¹²C, and ¹⁶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 ⁶Li from VMC calculations [135], AV18+IL7 in ¹²C from GFMC calculations [130], and AV18+UIX in ¹⁶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⁻¹, as they basically include relativistic corrections.

FIGURE 7

Figure 7. Longitudinal elastic form factor in ⁶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

Figure 8. Same as Figure 7 but for ¹²C. Experimental data are taken from reference [132]. Adapted from reference [89].

FIGURE 9

Figure 9. Same as Figure 7 but for ¹⁶O. Experimental data are from Sick, based on references [136–138]. Adapted from reference [89].

In ⁶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⁻¹, 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⁻¹. 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 ¹²C and ¹⁶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 ¹²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²(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 ¹²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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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.

Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

64. Kaiser N. Chiral three pi exchange N N potentials: results for diagrams proportional to g(A)⁴ and g(A)⁶. Phys Rev C. (2000) 62:024001. doi: 10.1103/PhysRevC.62.024001

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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.

Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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 ¹2C. Phys Rev Lett. (2013) 111:092501. doi: 10.1103/PhysRevLett.111.092501

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

138. Sick I. (1975).

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

PubMed Abstract | CrossRef Full Text | Google Scholar

140. Abbott BP, Abbott R, Abbott TD, Acernese F, Ackley K, Adams C, et al. GW170817: Measurements of neutron star radii and equation of state. Phys Rev Lett. (2018) 121:161101. doi: 10.1103/PhysRevLett.121.161101

PubMed Abstract | CrossRef Full Text | Google Scholar