Quantum Monte Carlo Methods in Nuclear Physics: Recent Advances

TL;DR

Recent advances in quantum Monte Carlo methods combined with effective field theories have significantly improved understanding of nuclear systems, neutron matter, and electroweak processes, enabling studies of heavier nuclei and neutron stars.

Contribution

This review highlights recent progress in QMC techniques with chiral EFT interactions, expanding applicability to complex nuclear systems and astrophysical phenomena.

Findings

Enhanced QMC methods enable studies of medium-mass nuclei.

Development of local chiral EFT interactions improves accuracy.

Insights into neutron matter and neutron star properties.

Abstract

In recent years, the combination of precise quantum Monte Carlo (QMC) methods with realistic nuclear interactions and consistent electroweak currents, in particular those constructed within effective field theories (EFTs), has lead to new insights in light and medium-mass nuclei, neutron matter, and electroweak reactions. This compelling new body of work has been made possible both by advances in QMC methods for nuclear physics, which push the bounds of applicability to heavier nuclei and to asymmetric nuclear matter and by the development of local chiral EFT interactions up to next-to-next-to-leading order and minimally nonlocal interactions including $\Delta$ degrees of freedom. In this review, we discuss these recent developments and give an overview of the exciting results for nuclei, neutron matter and neutron stars, and electroweak reactions.

Figures

Equations

$$H=T+\sum _{i<j}{V}_{ij}^{\text{NN}}+\sum _{i<j<k}{V}_{ijk}^{\text{3N}}+\cdots .$$

$$V_{ij}^{\text{NN},\pi }=\frac{f_{\pi }^2}{4\pi }\frac{m_{\pi }}{3}\left[Y(m_{\pi }r){\sigma }_i\text{dotproduct}{\sigma }_j+T(m_{\pi }r)Y(m_{\pi }r)S_{ij}\right]{\tau }_i\text{dotproduct}{\tau }_j,$$

$$O_{ij}^{1-8}$$

$$O_{ij}^{9-14}$$

$$O_{ij}^{15-18}$$

$$V=\sum _{\nu =0}^{\infty }{V}^{\nu }(C_i^{\nu }){\left(\frac{p}{{\Lambda }_b}\right)}^{\nu },$$

$$V_{\text{cont}}^{\nu =0}=C_{\text{mathds}1}\text{mathds}1+C_{\sigma }{\sigma }_1\text{dotproduct}{\sigma }_2+C_{\tau }{\tau }_1\text{dotproduct}{\tau }_2+C_{\sigma \tau }{\sigma }_1\text{dotproduct}{\sigma }_2{\tau }_1\text{dotproduct}{\tau }_2,$$

$$V_{\pi }^{\nu =0}(\text{p},{\text{p}}^{\prime })=-{\left(\frac{g_A}{2f_{\pi }}\right)}^2\frac{{\sigma }_1\text{dotproduct}\mathbf{q}{\sigma }_2\text{dotproduct}\mathbf{q}}{{\mathbf{q}}^2+m_{\pi }^2}{\tau }_1\text{dotproduct}{\tau }_2,$$

$$f_{\text{short}}=\alpha \exp \left(-{\left(\frac{r}{R_0}\right)}^n\right),f_{\text{long}}={\left(1-\exp \left(-{\left(\frac{r}{R_0}\right)}^{n_1}\right)\right)}^{n_2},$$

$$J_{\mathrm{NC}}^{\mu }=-2{\sin }^2{\theta }_W{J}_{\gamma ,S}^{\mu }+(1-2{\sin }^2{\theta }_W)J_{\gamma ,z}^{\mu }+J_{5,z}^{\mu },$$

$$J^{\mu }=\sum _i{j}^{\mu }(i)+\sum _{i<j}{j}^{\mu }(ij)+\cdots$$

$$j_{\gamma ,S}^0(i)$$

$${\mathbf{j}}_{\gamma ,S}^{\perp }(i)$$

$$j_{\gamma ,S}^{\parallel }(i)$$

$$j_{5,z}^0(i)=-\frac{G_A(Q^2)}{4m}{\tau }_{i,z}{\sigma }_i\text{dotproduct}(\mathbf{q}+{\mathbf{p}}_i),{\mathbf{j}}_{5,z}(i)=-\frac{G_A(Q^2)}{2}{\tau }_{i,z},$$

$$j_{5,\text{PS}}^{\mu }(i)=\frac{G_A(Q^2)}{m_{\pi }^2+Q^2}{\tau }_{i,\pm }{q}^{\mu }{\sigma }_i\text{dotproduct}\mathbf{q}.$$

$$\begin{array}{rl}{E}_V& =\frac{\text{expectationvalue}H{\Psi }_T}{\text{innerproduct}{\Psi }_T{\Psi }_T}=\frac{{\sum }_{\sigma \tau }\int \text{differential}\mathbf{R}{\Psi }_T^{\ast }(\mathbf{R},\sigma ,\tau )H{\Psi }_T(\mathbf{R},\sigma ,\tau )}{{\sum }_{\sigma \tau }\int \text{differential}\mathbf{R}{\Psi }_T^{\ast }(\mathbf{R},\sigma ,\tau ){\Psi }_T(\mathbf{R},\sigma ,\tau )},\end{array}$$

$$E_V=\frac{{\sum }_{\sigma \tau }\int \text{differential}\mathbf{R}P(\mathbf{R},\sigma ,\tau )H{\Psi }_T(\mathbf{R},\sigma ,\tau )/{\Psi }_T(\mathbf{R},\sigma ,\tau )}{{\sum }_{\sigma \tau }\int \text{differential}\mathbf{R}P(\mathbf{R},\sigma ,\tau )},$$

$$\hat{F}=\left(\prod _{i<j}{f}_c(r_{ij})\right)\left[\mathcal{S}\prod _{i<j}\left(1+F_{ij}\right)\right],$$

$$F_{ij}=f_{\tau }(r_{ij}){\tau }_i\text{dotproduct}{\tau }_j+f_{\sigma }(r_{ij}){\sigma }_i\text{dotproduct}{\sigma }_j+f_{\sigma \tau }(r_{ij}){\sigma }_i\text{dotproduct}{\sigma }_j{\tau }_i\text{dotproduct}{\tau }_j+f_t(r_{ij})S_{ij}+f_{t\tau }(r_{ij})S_{ij}{\tau }_i\text{dotproduct}{\tau }_j.$$

$$\displaystyle\ket{\Phi}=\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),\hskip 60.00009pt{\bf\it\sigma}_{1}\dotproduct{\bf\it\sigma}_{2}\ket{\Phi}=\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)\,.$$

$$|{\Psi }_0⟩\propto \underset{\tau \to \infty }{\lim }\exp [-(H-E_0)\tau ]|{\Psi }_T⟩,$$

$$\text{innerproduct}{\mathbf{R}}_N{\Psi }_T=\int \prod _{i=0}^{N-1}\text{differential}{\mathbf{R}}_i\text{matrixelement}{\mathbf{R}}_N\exp [-H\delta \tau ]{\mathbf{R}}_{N-1}\cdots \text{matrixelement}{\mathbf{R}}_1\exp [-H\delta \tau ]{\mathbf{R}}_0\text{innerproduct}{\mathbf{R}}_0{\Psi }_T,$$

$$\begin{array}{rl}{G}_{\delta \tau }(\mathbf{R},{\mathbf{R}}^{\prime })& =\text{matrixelement}{\mathbf{R}}^{\prime }\exp (-V\delta \tau /2)\exp (-T\delta \tau )\exp (-V\delta \tau /2)\mathbf{R}+\mathcal{O}(\delta {\tau }^3),\end{array}$$

$$\text{matrixelement}\mathbf{R}\exp (-V\delta \tau )\mathbf{R}\approx \mathcal{S}\prod _{i<j}\exp [-V_{ij}\delta \tau ].$$

$$|s_i⟩=a_i|p\uparrow ⟩+b_i|p\downarrow ⟩+c_i|n\uparrow ⟩+d_i|n\downarrow ⟩,s_i\equiv \{a_i,b_i,c_i,d_i\}$$

$$|S⟩=|s_1⟩\otimes \cdots \otimes |s_A⟩.$$

$${\sigma }_i^{\alpha }|S⟩=|s_1⟩\otimes \cdots \otimes {\sigma }_i^{\alpha }|s_i⟩\otimes \cdots \otimes |s_A⟩=|s_1⟩\otimes \cdots \otimes |s_i^{\prime }⟩\otimes \cdots \otimes |s_A⟩=|S^{\prime }⟩$$

$$\exp \left(-\frac{1}{2}\lambda {\hat{O}}^2\right)=\frac{1}{\sqrt{2\pi }}\int \text{differential}x\exp \left(-\frac{x^2}{2}+\sqrt{-\lambda }x\hat{O}\right),$$

$$\text{innerproduct}SR{\Psi }_V=⟨SR|\left[\prod _{i<j}{f}_c(r_{ij})\right]\left[1+\sum _{i<j}{U}_{ij}\right]|\Phi ⟩.$$

$${\rho }_{1,N}(A,r)\equiv \frac{1}{4\pi r^2}\text{expectationvalue}\sum _{i=1}^A\frac{1\pm {\tau }_{z,i}}{2}\delta (r-|{\mathbf{r}}_i-{\mathbf{R}}_{\text{cm}}|){\Psi }_0,$$

$${\rho }_{2,\mathcal{O}}(A,r)\equiv \frac{1}{4\pi r^2}\text{expectationvalue}\sum _{i<j}^A{\mathcal{O}}_{ij}\delta (r-|{\mathbf{r}}_{ij}|){\Psi }_0,$$

Full text

††journal: Annu. Rev. Nucl. Part. Sci.

1. 1 INTRODUCTION

2. 2 NUCLEAR INTERACTIONS

3. 3 ELECTROWEAK CURRENTS

4. 4 QUANTUM MONTE CARLO METHODS

5. 4.1 Variational Monte Carlo

6. 4.2 Green’s Function Monte Carlo

7. 4.3 Auxiliary Field Diffusion Monte Carlo

8. 5 RESULTS IN LIGHT AND MEDIUM-MASS NUCLEI

9. 5.1 Energies, Radii, and Scattering

10. 5.2 Distributions and Short-Range Correlations

11. 6 THE EQUATION OF STATE OF NEUTRON MATTER AND NEUTRON STARS

12. 7 RESULTS IN ELECTROWEAK REACTIONS

13. 8 SUMMARY AND OUTLOOK

1 INTRODUCTION

Over the past decade, significant progress has been made in the theoretical description of strongly-interacting nuclear systems. This progress is reflected in an increasingly accurate prediction of nuclear-structure observables for heavier nuclei, including, e.g. radii, masses, and neutron-separation energies. In particular, compelling progress has been made in ab initio nuclear structure, where the many-body Schrödinger equation is solved with controlled approximations and protons and neutrons are assumed to be the relevant degrees of freedom. This progress includes, for example, the discovery of new shell closures in neutron-rich nuclei, new studies of doubly magic nuclei, and the role of short-range correlations in weak transitions [1, 2, 3, 4, 5, 6]. In addition, the description of nuclear and neutron matter, and therefore the symmetry energy, has also become more accurate, and now includes reliable uncertainty estimates that are important for the extrapolation to the density regime encountered in neutron stars [7, 8]. These advances have been steered by the developments of systematic nuclear interactions and reliable many-body methods.

Systematic Hamiltonians from nuclear effective field theories (EFTs) [9, 10, 11] have played a key role in obtaining reliable results for nuclear systems. These EFT Hamiltonians are rooted in the symmetries of the fundamental theory of strong interactions, quantum chromodynamics (QCD), but describe the dynamics of nuclear system in terms of nucleonic degrees of freedom.{marginnote}[]\entryEFTeffective field theory \entryQCDquantum chromodynamics They are based on a power-counting scheme that allows for the derivation of nuclear interactions and consistent electroweak currents in a systematically improvable fashion. Furthermore, nuclear EFTs provide a recipe to estimate theoretical uncertainties, a key ingredient for a meaningful comparison with experimental data. In addition, nuclear EFTs naturally predict many-body forces, which are necessary for the correct description of nuclear systems. Pionless EFT and chiral (pion-full) EFT, both used in QMC calculations, will be briefly discussed in this review.

Thanks to the increasing availability of computing resources and the development of new algorithms, nuclear ab initio methods have extended their reach to medium-heavy nuclei. Among these many-body methods, quantum Monte Carlo (QMC) techniques are known for their accuracy in describing properties of light nuclei up to \isotope[12]C, see e.g. References [12, 13]. Among these QMC approaches are the Green’s Function Monte Carlo (GFMC) and the auxiliary field diffusion Monte Carlo (AFDMC) methods. They solve the Schrödinger equation by exploiting an imaginary-time evolution to enhance the ground-state component from a starting trial wave function. While some approximations are made during this evolution, mainly to cope with the fermion sign problem, the final results can be considered “stochastically exact”, as expectation values are estimated on finite Monte Carlo samples.{marginnote}[]\entryQMCquantum Monte Carlo\entryGFMCGreen’s function Monte Carlo\entryAFDMCauxiliary field diffusion Monte CarloThe GFMC method has been used to successfully predict the spectra and electroweak processes of nuclei with $A\leq 12$, where $A$ is the number of nucleons, with a percent-level accuracy. Because it sums over all spin/isospin states, the GFMC scales exponentially with $A$, which presently prevents its applicability to $A>12$ nuclei. The AFDMC method, on the other hand, uses Hubbard-Stratonovich transformations to sample the spin/isospin degrees of freedom and achieve a polynomial scaling in $A$. This has enabled the computation of systems with larger $A$, such as \isotope[16]O and neutron matter, at the cost of using somewhat simplified wave functions.

In this review, we present state-of-the-art QMC results for nuclei up to \isotope[16]O and neutron matter using nuclear EFT interactions, and electroweak processes using realistic phenomenological potentials, and discuss future directions. The review is structured as follows. In Section 2 we discuss nuclear interactions: Starting with phenomenological ones for context and moving to those based on chiral EFT. In Section 3 we briefly discuss electroweak currents. In Section 4 we introduce the GFMC and AFDMC methods for nuclear physics. In Sections 5, 6 and 7 we present recent results for light and medium-mass nuclei, neutron matter and neutron stars, as well as for electroweak reactions.

2 NUCLEAR INTERACTIONS

The fundamental degrees of freedom for nuclear systems are quarks and gluons, whose dynamics are primarily governed by the QCD Lagrangian. However, a description of, e.g., atomic nuclei, in terms of these degrees of freedom requires nonperturbative lattice techniques. Because of their tremendous computing cost, they are currently not practical for $A\gtrsim 2$ at a physical value of the pion mass [14].

Instead, at the energy regime relevant for the description of nuclear systems, the effective degrees of freedom are point-like nucleons, whose dynamics are dictated by the nonrelativistic Hamiltonian

[TABLE]

In the above equation, $T$ denotes the kinetic energy, $V_{ij}^{\text{NN}}$ is the two-nucleon (NN) interaction between nucleons $i$ and $j$, $V_{ijk}^{\text{3N}}$ is the three-nucleon (3N) interaction between nucleons $i$, $j$, and $k$, and the ellipsis indicate interactions involving more than three particles. As suggested by nuclear matter studies [15], these four-nucleon (and beyond) interactions are small compared to the current level of precision and can safely be omitted.{marginnote}[]\entryNNtwo-nucleon\entry3Nthree-nucleon

Traditionally, nuclear interactions have been constructed relying on meson-exchange models, e.g. in the CD-Bonn potential [16], or with the goal of reproducing scattering data with very high accuracy, e.g., nuclear interactions of the Argonne type [17]. Phenomenological Argonne NN interactions have been extensively and successfully used in a number of GFMC and AFDMC calculations. They describe the NN interaction by explicitly including the long-range one-pion exchange (OPE) interaction and a set of intermediate- and short-range terms that model the more complicated multi-pion exchanges and short-range dynamics. The OPE is given by{marginnote}[]\entryOPEone-pion exchange

[TABLE]

where $f_{\pi}$ is the $\pi\text{N}$ coupling constant, $m_{\pi}$ is the average pion mass, $S_{ij}=3{\bf\it\sigma}_{i}\dotproduct\mathbf{\hat{r}}\,{\bf\it\sigma}_{j}\dotproduct\mathbf{\hat{r}}-{\bf\it\sigma}_{i}\dotproduct{\bf\it\sigma}_{j}$ is the tensor operator in coordinate space with the Pauli matrices ${\bf\it\sigma}$, $Y(x)=\exp(-x)/x$ is the Yukawa function and $T(x)=(1+3/x+3/x^{2})$. The short-range divergent behavior of both $T$ and $Y$ is regulated by multiplying them by $f_{R}(x)=1-\exp(-cx^{2})$ where $c$ is a cutoff parameter, typically taken to be $c=2.1$ $\text{fm}^{-2}$. The intermediate- and short-range parts are modeled by a set of spin-/isospin- and momentum-dependent operators multiplied by $T^{2}$ (to approximate the two-pion-exchange) and Woods-Saxon-like radial functions respectively. The latest version in this class of potentials, denoted as Argonne $v_{18}$ (AV18) [17], is expressed in terms of 18 operators:

[TABLE]

In the above equations, $\mathbf{L}$ is the relative angular momentum of the pair, $\mathbf{S}$ is the total spin, and $T_{ij}=3\tau_{i}^{z}\tau_{j}^{z}-{\bf\it\tau}_{i}\dotproduct{\bf\it\tau}_{j}$ is the isotensor operator. All parameters of AV18 have been fit to the Nijmegen NN scattering database with $\chi^{2}/{\rm datum}\simeq 1$. Simplified versions of these interactions, comprising only a subset of the operators reported in Equation 5 are available. For instance, the Argonne $v_{8}^{\prime}$ ($\text{AV8}^{\prime}$) interaction, widely used in neutron-matter studies, only contains the first 8 operators, and other even simpler interactions have been explored [18]. In addition to the NN forces, phenomenological 3N interactions have been developed. They are generally expressed as a sum of a two-pion-exchange $P$-wave term (Fujita-Miyazawa), a two-pion-exchange $S$-wave contribution, a three-pion-exchange contribution, and a 3N contact. More specifically, the Urbana IX (UIX) [19] interaction contains only the first and last terms, while the Illinois 7 (IL7) [20] potential contains all four contributions. The UIX 3N interaction is fit to reproduce the ground-state energies of \isotope[3]H and \isotope[4]He and the saturation-point of symmetric nuclear matter, while the IL7 interaction was fit to the low-lying spectra of nuclei in the mass range $A=\text{3--10}$. Phenomenological interactions have been successfully used in a multitude of QMC calculations of nuclear systems, and results have been reviewed in e.g., Reference [12].

However, these interactions suffer from important shortcomings. Since they are constructed in an empirical way without a clear guiding principle, it is not possible to assess theoretical uncertainties associated with modeling nuclear dynamics. Also, it is not clear how to improve these interactions, especially in the 3N sector. For example, though the AV18+IL7 Hamiltonian leads to a description of more than 100 ground- and excited-state energies up to $A=12$ in good agreement with experimental data, it fails to provide sufficient repulsion in pure neutron matter [21]. On the other hand, the AV18+UIX model, while providing a reasonable description of nuclear matter properties, does not satisfactorily reproduce the spectrum of light nuclei. In addition, the derivation of consistent electroweak currents is not straightforward.

A solution to the previously discussed limitations of phenomenological interactions has been presented with the advent of nuclear EFTs [9, 10, 11]. Nuclear EFTs exploit the symmetries of QCD and enable a systematic approach to nuclear forces based on a low-momentum expansion. Within the nuclear EFT approach, one starts from hadronic degrees of freedom relevant for the system at hand. Additional degrees of freedom, e.g. heavier mesons or even the nucleon substructure, relevant only at higher energy scales beyond those treated within the EFT, are integrated out. This so-called “separation of scales” determines the breakdown scale, $\Lambda_{b}$, of the theory, and can be used to construct a systematic EFT: One writes down the most general Lagrangian consistent with all symmetries of QCD, and uses a power-counting scheme to arrange the terms according to their importance, typically in powers of $p/\Lambda_{b}$, where $p$ is a typical momentum scale in the nuclear system. The resulting scheme is valid only when $p\ll\Lambda_{b}$.

The most general EFT Lagrangian contains an infinite series of interaction terms,

[TABLE]

where $V^{\nu}(C_{i}^{\nu})$ is the contribution at order $\nu$ which depends on low-energy couplings (LECs) $C_{i}^{\nu}$. The LECs encode the unresolved physics that is integrated out and are determined by fitting experimental data.{marginnote}[]\entryLECslow-energy couplingsIn a converging EFT, the LECs are natural, i.e. of order $1$ and, hence, the higher-order contributions to $V$ decrease in magnitude. This permits the truncation of the series expansion at a certain order $\nu$. By going to higher orders, one can work to a desired accuracy at the cost of computing more diagrams. (In this review we denote leading order by LO, next-to-leading order by NLO, and $\text{next-to-}\cdots\text{-leading order}$ by NxLO, with $x$ the number of orders beyond LO). {marginnote}[]\entryLOleading order\entryNLOnext-to-leading order\entryNxLO$\overbrace{\text{next-to-}\cdots}^{x\ \text{times}}$ -leading order This systematic expansion can be used to estimate meaningful theoretical uncertainties. Another advantage of nuclear EFTs is that the procedure described above leads to the natural appearance of many-body forces.

At very low momentum scales, $p\ll m_{\pi}$, pions can be integrated out and nuclear interactions reduce to contact interactions with different numbers of derivatives. QMC calculations with pionless EFT interaction have been used to analyze lattice QCD calculations with great success [22, 23], e.g., a pionless-EFT Hamiltonian was used in AFDMC calculations aimed at extending lattice QCD predictions to \isotope[16]O [22]. Leading-order results indicate that for $m_{\pi}=805$ MeV and $m_{\pi}=510$ MeV, \isotope[16]O is not stable against breakup into four \isotope[4]He nuclei. We refer the reader to References [9, 10] for more details on pionless EFT.

While pionless EFT has been used successfully in low-energy nuclear physics, see, e.g., Reference [24], typical momenta in nuclear many-body systems are of the order of $m_{\pi}$, and therefore larger than its breakdown scale. Chiral EFT is based on the observation that pions naturally emerge as pseudo Goldstone bosons associated with the spontaneous breaking of the approximate chiral symmetry of QCD. Within this “pionfull” chiral EFT [10, 11], nuclear interactions are comprised of both contact terms, written in a general operator basis, and one- and multi-pion-exchange interactions. Modern chiral EFT interactions are based on Weinberg power counting [25, 26], but alternative power-counting schemes have been suggested; see, e.g., References [27, 28, 29, 30, 31].

At LO, the contact interactions are given by the momentum-independent contributions

[TABLE]

and the pion-exchange interactions are given by the well-known OPE interaction,

[TABLE]

where $\mathbf{p}$ and $\mathbf{p}^{\prime}$ are the relative nucleon momenta before and after the interaction, and $\mathbf{q}=\mathbf{p}-\mathbf{p}^{\prime}$ is the momentum transfer. At higher orders, more complicated interaction pieces contribute, i.e., momentum-dependent contacts $\sim p^{\nu},p^{\prime\nu}$, tensor contacts, and multiple-pion exchanges, and we refer the reader to References [10, 11] for more details. For example, at N3LO the NN operator basis includes a set of operators similar to the phenomenological ones of Equation 5, but the operators in the contact and pion sector appear in a systematic fashion. A chief advantage of the chiral EFT formulation is that 3N interactions are consistent with the NN potential, i.e., the same vertices in both sectors have the same LECs. The leading 3N forces are given by a two-pion exchange interaction ($V_{C}$), an OPE–contact interaction ($V_{D}$) and a 3N contact contribution ($V_{E}$), with only two unknown LECs to be determined [32].

Chiral EFT Hamiltonians have been extensively employed in recent years in a host of nuclear many-body methods. At the same time, new strategies to improve chiral interactions and reduce the theoretical uncertainties have been proposed. Interactions up to fifth order in the chiral expansion have recently been developed [33, 34]. These interactions reproduce the pp and np scattering data from the Granada-2013 database with a $\chi^{2}/\text{datum}\sim 1$, matching the precision of phenomenological potentials. Furthermore, potentials with explicit $\Delta$ degrees of freedom have been constructed [35, 36] which in principle enable a detailed comparison of the order-by-order convergence in both the $\Delta$-less and $\Delta$-full theories. New optimization schemes are being explored [37, 38], which may improve the multidimensional fits of the LECs in the chiral Hamiltonians. Such schemes may play an important role when working at N3LO or beyond, where 24 or more LECs have to be simultaneously determined. Also, new forms of uncertainty estimates using Bayesian statistical tools are being explored [39, 40]. These tools are necessary to enable a meaningful comparison of theoretical calculations with experimental data and to study the convergence of the chiral expansion in a systematic way. Lastly, new regularization schemes are being explored. As for phenomenological interactions, when employing nuclear EFT Hamiltonians in many-body calculations, the high-momentum components of the interactions have to be regularized to prevent divergences. This is achieved by introducing a regulator function $f_{R}$ that is $\mathcal{O}(1)$ at low momenta and $\mathcal{O}(0)$ at high momenta, e.g. $f_{R}(p)=\exp\left(-\left(p/\Lambda\right)^{n}\right)$, where $p$ is the regulated momentum scale, $\Lambda$ is the cutoff scale that determines which contributions are discarded, and $n$ is an integer. Typically, chiral interactions were regulated nonlocally, but recently, local and semilocal regularization schemes have been proposed, see e.g., References [41, 33] for more details.

In this review, we focus on the application of chiral interactions within QMC methods. Chiral interactions are often constructed in momentum space and contain various momentum dependencies in terms of the average incoming and outgoing nucleon momenta $\mathbf{p}$ and $\mathbf{p}^{\prime}$. However, the GFMC and AFDMC methods are best suited to use local interactions as input to solve the many-body Schrödinger equation, i.e., the interactions should depend on the momentum transfer $\mathbf{q}$. This is true for pion-exchange interactions up to N2LO but not generally true for contact interactions. However, it is possible to construct fully local chiral interactions up to N2LO by choosing a local contact operator basis using Fierz ambiguities and local regulators. Local chiral interactions have recently been developed for use in the GFMC and AFDMC methods within $\Delta$-less [41, 42] and $\Delta$-full [35] chiral EFT.

The local regulators for short- and long-range interactions typically are of the form

[TABLE]

where the exponents $n$, ${n_{1}}$, and ${n_{2}}$ describe the sharpness of the regulator and $\alpha$ is a normalization factor. To cut off the divergences at the origin also of multi-pion-exchange interactions, the exponents ${n_{1}}$ and ${n_{2}}$ need to be chosen sufficiently large. Note that the minimally nonlocal interactions of References [35, 43] use a Woods-Saxon–like functional form for the long-range regulator, which however, in practice differs only nominally from Equation 9. More details on the construction of local chiral interactions can be found in the above-mentioned references.

Nevertheless, local chiral interactions suffer from regulator artifacts that do not appear for typical nonlocal chiral interactions. When constructing local contact potentials, Fierz ambiguities among different contact operators at a given order are exploited to eliminate nonlocal terms. However, local regulators break this Fierz rearrangement freedom and short-range regulator artifacts arise [44, 45]. It can be shown that the latter are of the same order as higher-order terms in the chiral expansion and, therefore, are cured by explicitly including higher-order contact interactions [46, 45]. While this effect is not very dramatic in the NN sector at, e.g. N2LO, it is quite sizable in the 3N sector. These regulator artifacts lead to an ambiguity in the shorter-range $V_{D}$ and $V_{E}$ topologies [47], which increases the theoretical uncertainty at typical cutoff scales. While local regulators preserve the analytic structure of the partial-wave amplitude near threshold for long-range pion-exchange interactions in NN scattering [48], they also lead to larger regulator artifacts for the 3N two-pion exchange interaction at typical cutoff scales; see References [49, 46]. Hence, local chiral interactions together with QMC methods lead to exciting insights, but one has to carefully consider regulator artifacts that appear and their influence on the results.

3 ELECTROWEAK CURRENTS

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. Those associated with neutral-current transitions can be written as [50]

[TABLE]

where $\theta_{W}$ is the Weinberg angle ($\sin^{2}\theta_{W}=0.23122$ [51]), $J^{\mu}_{\gamma,S}$ and $J^{\mu}_{\gamma,z}$ are the isoscalar and isovector pieces of the electromagnetic current $J^{\mu}_{\text{EM}}=J^{\mu}_{\gamma,S}+J^{\mu}_{\gamma,z}$, and $J^{\mu}_{5,z}$ denotes the isovector term of the axial current.

Analogously to the nuclear interaction, electroweak currents can also be expressed as an expansion in many-body operators that act on nucleonic degrees of freedom

[TABLE]

The one-body charge and current operators have the standard expressions [52] obtained from the nonrelativistic reduction of the covariant single-nucleon current, and include terms proportional up to $1/m^{2}$, $m$ being the nucleon mass. The transverse ($\perp$) and longitudinal ($\parallel$) components to the momentum transfer $\mathbf{q}$ of the isoscalar term read

[TABLE]

where $\mathbf{p}_{i}$ is the momentum of the $i$th nucleon and current conservation has been used to relate $j^{\parallel}_{\gamma,S}(i)$ to $j^{0}_{\gamma,S}(i)$. The corresponding isovector components of $j^{\mu}_{\gamma,z}(i)$ are obtained by $G_{E,M}^{S}(Q^{2})\longrightarrow G_{E,M}^{V}(Q^{2})\,\tau_{i,z}$, with $G_{E}^{S/V}(Q^{2})$ and $G_{M}^{S/V}(Q^{2})$ being the isoscalar/isovector combinations of the proton and neutron electric ($E$) and magnetic ($M$) form factors.

Omitting for brevity terms proportional to $1/m^{2}$, the isovector components of the axial weak neutral current $j^{\mu 5}_{z}$ are given by

[TABLE]

where $G_{A}(Q^{2})$ is the axial form factor of the nucleon, which is usually parametrized by a dipole $G_{A}(Q^{2})=g_{A}/(1+Q^{2}/M_{A}^{2})$. The nucleon axial-vector coupling constant is taken to be $g_{A}=1.2723$ [51] and the axial mass $M_{A}=1.03$ GeV [53], as obtained from an analysis of pion electroproduction data [54] and measurements of the reaction $\nu_{\mu}+p\to n+\mu$ [55]. Uncertainties in the $Q^{2}$ dependence of the axial form factor have a significant impact upon neutrino-nucleus cross-section predictions. In particular, the dipole parametrization has been the subject of intense debate: An alternative “z-expansion” analyses [56] has been proposed and dedicated lattice-QCD calculations of $G_{A}(Q^{2})$ have been carried out [57].

The charge-changing weak current is written as the sum of polar- and axial-vector components $J^{\mu}_{\text{CC}}=J^{\mu}_{\gamma\pm}+J^{\mu}_{5,\pm}$, whose one-body contributions can be obtained from $j^{\mu}_{\gamma,z}(i)$ and $j^{\mu}_{5,z}(i)$ by replacing $\tau_{i,z}/2\longrightarrow\tau_{i,\pm}=(\tau_{i,x}\pm\mathrm{i}\tau_{i,y})/2$. In addition, one has to retain the induced pseudoscalar contribution [58, 59]

[TABLE]

The gauge invariance of the theory imposes that the electromagnetic charge and current operators must satisfy the continuity equation $\mathbf{q}\dotproduct\mathbf{j}_{\text{EM}}=\commutator{H}{\rho_{\text{EM}}}$ where $\rho_{\text{EM}}\equiv J^{0}_{\text{EM}}$, hence providing an explicit connection between the nuclear interactions and the longitudinal component of the current operators. For instance, the isospin and momentum dependence of the NN interactions leads to nonvanishing commutators with the one-body charge operator and hence to the emergence of two-body terms in the current operator. In QMC calculations, both the phenomenological “Standard Nuclear Physics Approach” (SNPA) and chiral EFT have been exploited to derive many-body current operators.

The SNPA isoscalar and isovector components of the nuclear electromagnetic current $J^{\mu}_{\gamma,S}$ and $J^{\mu}_{\gamma,z}$, whose explicit expressions can be found in Reference [50], lead to a satisfactory description of static properties (charge radii, quadrupole moments, and M1 transition widths), charge and magnetic form factors of nuclei with $A\leq 12$ [52, 60, 61, 62], and electromagnetic response functions [63, 64], which will be discussed in Section 7. They consist of “model-independent” and “model-dependent” terms [65]. The former are obtained from the NN interaction, and by construction satisfy current conservation. The leading operator is the isovector “$\pi$-like” current but important contributions also arise from $\rho$-like terms. The additional two-body currents arising from the momentum-dependence of the NN interaction have been numerically proven to be much smaller [61].

The transverse components of the two-body currents cannot be directly linked to the nuclear Hamiltonian. In the latest applications of the SNPA formalism [50, 62, 63, 64], they include the isoscalar $\rho\pi\gamma$ transition and the isovector current associated with the excitation of intermediate $\Delta$-isobar resonances. The $\rho\pi\gamma$ couplings are extracted from the widths of the radiative decay $\rho\to\pi\gamma$ [66] and the $Q^{2}$ dependence of the electromagnetic transition form factor is modeled assuming vector-meson dominance [67]. Among the model-dependent currents, those associated with the $\Delta$ isobar are the most important ones.

One of the chief advantages of the chiral EFT formulation is that electroweak currents are constructed in a consistent fashion with the nuclear interaction. Since the chiral Lagrangian is gauge invariant, nuclear electromagnetic currents automatically satisfy the continuity equation, order by order, with the corresponding chiral potentials [68]. In particular, an important advantage of chiral EFT over the SNPA is the explicit connection between the 3N interaction and the two-body axial current. For example, the LEC $c_{D}$ entering the 3N potential at N2LO is related to the LEC of the two-body contact axial current [69, 70]. Similarly to the EFT interactions discussed in Section 2, chiral EFT currents can also be systematically organized in powers of $(p/\Lambda_{b})^{\nu}$, where the generic low-momentum scale includes the momentum transferred by the external electroweak probe. Note that single-nucleon structure effects have to be accounted for by introducing appropriate form factors. Because of the shortcomings of chiral EFT nucleonic form factors for $Q^{2}\gtrsim 0.1$ GeV2 [71, 72], even in chiral EFT formulations, parametrized versions of the latter are usually employed.

Over the last decade, extensive work to construct two-body electromagnetic current operators has been carried out by the JLab-Pisa and by the Bochum-Bonn groups using standard time-ordered perturbation theory [73, 74, 75] and the method of unitary transformations [76, 77], respectively. The LO vector current $\mathbf{j}_{\gamma,S}$ and $\mathbf{j}_{\gamma,z}$, corresponding to $\nu=-2$, are the same as those obtained within the SNPA, and are reported in the second and third lines of Equation 12. At N2LO ($\nu=0$) one needs to account for relativistic corrections to the one-body currents, while at N3LO ($\nu=1$) there are pure two-pion exchange and short-range one-loop contributions. At this order, additional “minimal” and “nonminimal” contact diagrams, defined in terms of two new LECs, need to be accounted for. The former originate from the contact chiral EFT NN potential at NLO through the minimal substitution $\mathbf{p}\to\mathbf{p}-\mathrm{i}e\mathbf{A}$, where $e$ is the electric charge and $\mathbf{A}$ is the vector photon field. Consequently, these contributions, needed for the continuity equation to be satisfied, involve the same LECs as the contact NN term – and hence can be determined by fitting NN scattering data. The nonminimal contributions arise from the field strengths $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$, which transform covariantly under chiral symmetry [74], and their LECs need to be fixed against electromagnetic observables. It should be noted that the isoscalar and the isovector contributions of the minimal terms correspond to the model-dependent $\rho\pi\gamma$ and $\Delta$-excitation transverse currents of the SNPA.

Axial currents were also recently derived within chiral EFT up to one-loop in References [78, 79], including pion-pole contributions. The latter are crucial for the current to be conserved in the chiral limit and are suppressed in low-momentum transfer processes. In the axial current, OPE contributions enter at N3LO $(\nu=0)$ and involve the LECs $c_{3}$, $c_{4}$, and $c_{6}$, which also enter the NN interaction. The situation is different for the axial charge, as pion-range contributions already enter at NLO $(\nu=-1)$. One-loop corrections to the axial current appear at N4LO $(\nu=1)$ and were first utilized in the calculation of the tritium Gamow-Teller matrix element [80]. On the other hand, at N3LO $(\nu=1)$ the calculation of the one-loop contributions has been carried out in Reference [81], aimed at studying the inclusive neutrino scattering off the deuteron at low energies. Electroweak currents which explicitly include the $\Delta$ excitation, consistently with the nuclear interactions discussed in Section 2, were derived in Reference [82] and applied to the calculation of the the Gamow-Teller matrix element contributing to tritium $\beta$ decay.

4 QUANTUM MONTE CARLO METHODS

Quantum Monte Carlo methods provide powerful tools to solve for the ground state of strongly interacting many-body systems. These methods have been used for problems in quantum chemistry and materials with very high accuracy, see e.g. References [83, 84, 85, 86]. Several different QMC implementations exist, both for bosons and fermions; here we will limit our description to the particular methods recently used to calculate properties of nuclear systems.

4.1 Variational Monte Carlo

The Variational Monte Carlo (VMC) method is used to calculate observables (e.g. the energy) of a many-body system once a suitable guess for its wave function $\Psi_{T}$ (the “trial” wave function) is provided.{marginnote}[]\entryVMCVariational Monte Carlo The variational energy $E_{V}$ of an $A$-nucleon system is given by

[TABLE]

where $\mathbf{R}=\{\mathbf{r}_{1}\ldots\mathbf{r}_{N}\}$, $\sigma=\{\sigma_{1}\dots\sigma_{N}\}$, and $\tau=\{\tau_{1}\dots\tau_{N}\}$ include all particles’ positions $\mathbf{r}_{i}$, spins $\sigma_{i}$, and isospins $\tau_{i}$, and $H$ is the nuclear Hamiltonian. The energy $E_{V}$ provides an upper bound to the ground-state energy $E_{0}$ and is equal to $E_{0}$ only if $\Psi_{T}$ coincides with the true ground-state wave function of the system, $\ket{\Psi_{T}}=\ket{\Psi_{0}}$. The calculation of $E_{V}$ requires the numerical evaluation of a multidimensional integral, but the high dimensionality limits standard numerical integration techniques to very small systems.

Monte Carlo integration is a natural solution to this limitation. Equation 15 can be rewritten as

[TABLE]

where the function $P(\mathbf{R},\sigma,\tau)$ is a probability distribution, and one natural choice is $P(\mathbf{R},\sigma,\tau)=\Psi_{T}^{\dagger}(\mathbf{R},\sigma,\tau)\Psi_{T}(\mathbf{R},\sigma,\tau)$. In the VMC method $P$ is used to sample a set of $M$ configurations in $\{\mathbf{R},\sigma,\tau\}$ space that are used to solve the integral above. A common way to generate such configurations is provided by the Metropolis algorithm, but many others are available. See e.g., Reference [86].

For strongly interacting systems, a common ansatz for variational wave function is $\ket{\Psi_{T}}=\hat{F}\ket{\Phi}$. The correlation operator $\hat{F}$, modeling the short-range correlations induced by the Hamiltonian, can generically be written as

[TABLE]

where we have omitted three-body correlations just for simplicity. In the above equation, $f_{c}(r)$ is a spin/isospin-independent correlation, and

[TABLE]

Evaluating the symmetrization operator $\mathcal{S}$ would require a factorial number of operations. In practice, the order of pairs for the left and right wave functions is instead sampled for each configuration. The radial correlations $f_{i}(r)$ include variational parameters that are chosen in order to minimize $E_{V}$.

The long-range antisymmetric part $\ket{\Phi}$ is typically a Slater determinant of single-particle orbitals, appropriate for the nuclear system of interest. For homogeneous matter, the orbitals can be plane waves or also include pairing correlations [87]. For nuclei, the single-particle orbitals are generally states written in the $ls$ or $jj$ basis that are properly combined to give the desired total angular momentum $J$ and isospin $T$ of the nucleus [88]. Within the GFMC method, $\ket{\Phi}$ consists of a set of amplitudes, each representing a particular spin/isospin configuration of the many-body state. For example, the spin amplitudes for 3 neutrons and the amplitudes after a spin/spin operator has acted on the state are represented by:

[TABLE]

The isospin is treated in a similar fashion, except that in this case the number of elements is smaller due to charge and/or total isospin conservation. For this reason, in GFMC calculations the number of many-body spin/isospin states (in the charge basis) is equal to $2^{A}\binom{A}{Z}$. Note that the coefficients associated with a given many-body spin/isospin state are such that the wave-function is fully anti-symmetric.

4.2 Green’s Function Monte Carlo

In the GFMC method, the ground state of the system is obtained with an imaginary-time projection

[TABLE]

where $\tau$ is the imaginary time, and $E_{0}$ is a parameter used to control the normalization (that we set to [math] in the following). The direct computation of the propagator $\exp[-H\tau]$ for arbitrary $\tau$ is typically not possible, but for small imaginary times $\delta\tau=\tau/N$ with $N$ large, the calculation is tractable, and the full propagation to large imaginary times $\tau$ can be obtained through the path integral

[TABLE]

where Monte Carlo techniques are used to sample the paths ${\mathbf{R}_{i}}$. In practice, a set of configurations, typically called walkers, are simultaneously evolved in imaginary time, and then used to calculate observables once convergence is reached. In the GFMC method, each walker contains the nucleon positions and a complex amplitude for each spin/isospin state of the nucleus, implying an unfavorable exponential scaling with the number of nucleons.

The most common and easiest approximation for the short-time propagator $G_{\delta\tau}(\mathbf{R},\mathbf{R}^{\prime})\equiv\matrixelement{\mathbf{R}^{\prime}}{\exp[-H\delta\tau]}{\mathbf{R}}$ is obtained by using the Trotter-Suzuki expansion:

[TABLE]

but more sophisticated and accurate ways to reduce the time-step error above are available [89]. Here, $T$ is the nonrelativistic kinetic energy giving rise to the free-particle propagator $\matrixelement{\mathbf{R}^{\prime}}{\exp[-T\delta\tau]}{\mathbf{R}}\propto\exp[-(\mathbf{R}-\mathbf{R}^{\prime})^{2}/\lambda^{2}]$, with $\lambda^{2}=4\frac{\hbar^{2}}{2m}\delta\tau$, yielding a Gaussian diffusion for the particles. The matrix $V$ is the spin/isospin-dependent interaction:

[TABLE]

Each pairwise interaction can be simply evaluated by exponentiating a small spin/isospin matrix. This treatment is adequate for static spin/isospin-dependent NN interactions, as they are diagonal in coordinate space. In practice one also needs to include momentum-dependent spin-orbit (LS) NN interactions as well as 3N interactions. For these and other details, see Reference [12].{marginnote}[]\entryLSspin-orbitThe method is exact in the limit of $\delta\tau=0$, and in practice small values of the time step are used to extrapolate to $\delta\tau=0$.

In addition to ground states, excited states have been accessed in GFMC calculations. The diffusion $\lim_{\tau\to\infty}\mathrm{e}^{-H\tau}\ket{\Psi_{T}}\to\ket{\Psi_{0}}$ drives $\ket{\Psi_{T}}$ to the lowest-energy eigenstate with the same quantum numbers as $\ket{\Psi_{T}}$. Thus, to obtain an excited state with distinct quantum numbers from the ground state, one need only construct a trial wave function with the appropriate quantum numbers. If the excited-state quantum numbers coincide with the ground state, more care is needed, but results for such states can still be obtained [90].

4.3 Auxiliary Field Diffusion Monte Carlo

The basic idea of the AFDMC method [91] is to achieve a better scaling with $A$ than GFMC by sampling the spin/isospin states rather than explicitly considering all of them. Let us define the single-nucleon spinor as

[TABLE]

where $a_{i}$, $b_{i}$, $c_{i}$ and $d_{i}$ are complex numbers, and $\{\ket{p\uparrow},\ket{p\downarrow},\ket{n\uparrow},\ket{n\downarrow}\}$ is the proton-up, proton-down, neutron-up and neutron-down basis. The spin/isospin states employed in AFDMC are products of single-particle states

[TABLE]

whose dimensionality scales linearly with the number of particles. On the other hand, computing the Slater determinant of the mean-field part of the wave-function, $\innerproduct{S}{\Phi}$, scales polynomially with the number of nucleons, as it requires $A^{3}$ operations. The main issue associated with the single-particle spin states of Equation 41 is that they are not closed with respect to the application of a quadratic spin (or isospin) operators. For instance, it can be easily shown that ${\bf\it\sigma}_{i}\dotproduct{\bf\it\sigma}_{j}\ket{S}\neq\ket{S^{\prime}}$. For a linear operator instead, one finds

[TABLE]

Realistic nuclear Hamiltonians contain quadratic spin/isospin operators. Thus, the imaginary-time propagation in AFDMC is carried out using the Hubbard-Stratonovich transformation, suitable to linearize such quadratic spin/isospin dependence:

[TABLE]

where the $x$ are called auxiliary fields. For example, for a spin-dependent interaction, $v(r_{ij}){\bf\it\sigma}_{i}\dotproduct{\bf\it\sigma}_{j}\delta\tau$, it is possible to define new operators such that $\sigma_{i}^{\alpha}\sigma_{j}^{\alpha}\to\hat{O}^{2}$ and $\lambda\to 2v(r_{ij})\delta\tau$ [12]. The strategy of the AFDMC method is to propagate particle positions in the continuum as commonly done in the GFMC method, but also sample the spin states of the nucleons in the continuum using the auxiliary fields. For more details see References [92, 12, 93].

The AFDMC trial (variational) wave function must be antisymmetric under the exchange of pairs. Omitting again three-body correlations, it is usually written as

[TABLE]

The long range part is given by $\innerproduct{SR}{\Phi}={\cal A}\{\phi_{\alpha_{1}}(r_{1},s_{1})\dots\phi_{\alpha_{N}}(r_{N},s_{N})\}$ with $\phi_{\alpha_{i}}(\mathbf{r}_{j},s_{j})=\innerproduct{\mathbf{r}_{j},s_{j}}{\phi_{\alpha_{i}}}=\innerproduct{\mathbf{r}_{j}}{f_{n_{i}}(r)}\innerproduct{s_{j}}{\xi_{i}}$ being single-particle orbitals, which are either constructed from a Wood-Saxon Hamiltonian or opportunely rescaled from mean-field calculations. While a single Slater determinant suffices for closed-shell nuclei, a sum of them is required to described open-shell configurations. For homogeneous matter, the orbitals can be either plane waves or pairing correlations.

5 RESULTS IN LIGHT AND MEDIUM-MASS NUCLEI

The combination of ab initio QMC methods with interactions derived from chiral EFT is an exciting development and yields interesting new insights into nuclear structure. In this section, we present a selection of recent results highlighting the rapidly growing reach of this novel combination.

5.1 Energies, Radii, and Scattering

One of the most interesting findings that emerge from QMC simulations with chiral EFT interactions is the ability to simultaneously describe three different types of nuclear systems: The binding energies and charge radii of light nuclei, the LS splitting in the $P$-wave phase shifts of elastic $n$–$\alpha$ scattering, and the equation of state (EOS) of pure neutron matter.{marginnote}[]\entryEOSequation(s) of state Historically, these three systems could not be described by a single combination of phenomenological potentials. While the combination $\text{AV}18+\text{UIX}$ gives a good description of light nuclei and neutron matter, it cannot reproduce the LS splitting in the $J^{\pi}=3/2^{-}$ and $1/2^{-}$ partial-wave elastic $n$–$\alpha$ scattering phase shifts. The conclusion reached was that certain topologies in the 3N interaction (namely three-pion exchange ring diagrams) were necessary to account for this LS splitting [94]. On the other hand, the combination $\text{AV}18+\text{IL}7$, which includes such three-pion exchange ring diagrams, gives an excellent description of light nuclei and the LS splitting in the $P$ waves of elastic $n$–$\alpha$ scattering but produces a too soft EOS for neutron matter inconsistent with physical expectations  [21]. The three-pion-exchange ring diagrams used in the IL7 3N potential appear also in chiral EFT; however, not until N3LO, N4LO, or N5LO (depending on the number of intermediate $\Delta$ states and $\Delta$-full or $\Delta$-less EFT). Hence, the question emerges how well local chiral N2LO interactions perform in the three benchmark systems discussed above. We present results from References [47, 95, 93] to answer this question.

The binding energies and radii of light nuclei up to \isotope[16]O are presented in Figure 1 from AFDMC calculations reported in Reference [95, 93]. The two undetermined LECs appearing in the N2LO 3N interaction are fit to the \isotope[4]He binding energy and the LS splitting in the $P$-wave phase shifts of elastic $n$–$\alpha$ scattering. While the agreement for $A=4$ and 5 is by construction, the very good agreement between the AFDMC simulations and experiment persists up to \isotope[6]Li in the energies, after which it deteriorates somewhat. For the radii, the agreement is very good up to \isotope[16]O with an exception for \isotope[6]Li, which is also found in simulations using phenomenological potentials [12].

The LS splitting between different states is an important feature in light nuclei. For local chiral EFT interactions at N2LO, in Reference [47] the two 3N LECs were fit to the two $P$-wave phase shifts extracted from an $R$-matrix analysis of the data. While at NLO, the (NN-only) interaction generates too little splitting between the two partial waves (the nonresonant $1/2^{-}$ partial wave is well reproduced, but not the resonant $3/2^{-}$ wave), at N2LO with the addition of the 3N interactions, agreement with both partial waves can be reproduced well even for different parameterizations of the 3N interaction. In Reference [96], ground and some excited states of light nuclei up to \isotope[12]C have been calculated in GFMC using interactions derived from $\Delta$-full chiral EFT [43]: See Figure 2. Overall, the agreement with experiment is very good, with an RMS deviation from experiment of $<1$ MeV.

5.2 Distributions and Short-Range Correlations

In addition to energies and radii, QMC methods can provide detailed information on the distribution of nucleons in a nucleus in both coordinate and momentum space. These distributions are connected to experimental results in several ways. For example, the one-body point-proton and -neutron densities, defined as

[TABLE]

with $+$ for the proton ($N=p$) density and $-$ for the neutron ($N=n$) density, are related via Fourier transform to the longitudinal electric form factor $F_{L}(Q)$; see Figure 3. Overall, the comparisons of these electric charge form factors with experiment is very good, with the first diffraction minima well reproduced.

Two-body coordinate-space distributions,

[TABLE]

can also be related to experimentally observable quantities. One of the most interesting results to arise from the novel combination of EFT with QMC methods is the relation of the so-called two-body short-range correlation (SRC) scaling factors $a_{2}$ to the short-distance behavior of the ratio of the two-body central correlations.{marginnote}[]\entrySRCshort-range correlationThat is,

[TABLE]

with $A$ representing a nucleus with $A$ nucleons and $d$ representing the deuteron ($A=2$); See Reference [98] for more details. These SRC scaling factors are extracted from quasielastic scattering from nuclei at intermediate Bjorken $x$ values. In Figure 4 we compare experimental values for $a_{2}(A/d$) in light nuclei with values extracted from two-body distributions from GFMC calculations using local chiral EFT interactions at N2LO and phenomenological potentials as well as VMC calculations using phenomenological potentials. This novel idea not only sheds new light on two-body SRCs in nuclei, but may also help clarify the nature of the so-far elusive 3N SRCs and on the isospin dependence of the EMC effect through the EMC-SRC linear relationship; see [98] for more details.

So far, we have discussed distributions calculated in coordinate space (and their Fourier transforms). However, it is also possible to calculate distributions directly in momentum space, such as the two-nucleon momentum distribution (the probability of finding a pair of nucleons in a nucleus with relative momentum $\mathbf{q}$ and total center-of-mass momentum $\mathbf{Q}$)

[TABLE]

where $\mathcal{P}_{\text{NN}}(ij)=(1/4)(1\pm\tau_{z,i})(1\pm\tau_{z,j})$ is an isospin projector: See Reference [99] for more details. Such distributions with $\mathbf{Q}=0$ (so-called “back-to-back” pairs) are valuable to compare with exclusive electron scattering experiments ${}^{A}Z(\text{e},\text{e}^{\prime}\text{pp})/^{A}Z(\text{e},\text{e}^{\prime}\text{np})$, where the dramatic dominance of np pairs over pp pairs has been observed [100, 101, 102]; see Figure 5. The agreement between the experimentally extracted ratios and those calculated from two-nucleon momentum distributions using local chiral EFT interactions at N2LO is very good. In addition, the agreement between the phenomenological results using $\text{AV}18+\text{UX}$ (solid black line in Figure 5) is notable.

6 THE EQUATION OF STATE OF NEUTRON MATTER AND NEUTRON STARS

The EOS of PNM is closely related to the nuclear symmetry energy that can be studied at higher densities only in heavy-ion collisions, and to the EOS of neutron stars (NS), where proton fractions in the core are typically of the order of only few %.{marginnote}[]\entryPNMpure neutron matter\entryNSneutron starTherefore, PNM provides a natural bridge between astrophysical observations of NS and terrestrial nuclear experiments. Also, the EOS of PNM is an interesting model system to study nuclear interactions because interactions in PNM are simpler than in systems containing also protons. For example, only the $T=3/2$ isospin channel contributes to PNM (where $T$ is the total isospin) while the presence of protons also permits contributions from the $T=1/2$ channel. Nevertheless, the $T=3/2$ isospin channel is only weakly accessible by studying properties of nuclei.

While a complete calculation of nuclear matter with arbitrary proton fractions up to $x=0.5$ is still not possible with QMC methods, the AFDMC method has been widely used to calculate the EOS of PNM for many different nuclear interactions in the past years. In practice, in QMC methods the infinite system is simulated by a fixed number of neutrons in a periodic box at a given baryon density. In particular, simulations using 66 neutrons (33 spin up and 33 spin down) give results very close to the thermodynamic limit [103, 104].

In Figure 6 we present results for PNM using the AFDMC method with local chiral interactions up to N2LO. The three different bands correspond to using different short-range operator structures for the 3N contact interaction $V_{E}$ at N2LO as described in Reference [47] and discussed in Section 2; the differences are due to finite-cutoff effects and vanish in the limit of large (momentum-space) cutoffs. Each band depicts a truncation uncertainty estimate based on the order-by-order results at LO, NLO (both also shown in the figure), and N2LO. The results are compared to calculations for the phenomenological $\text{AV8}^{\prime}$ NN interactions and when additionally including the UIX 3N forces. Note, that the blue (lower) band produces an EOS that is very soft and leads to negative pressure at $\approx 1.5n_{0}$, which is unphysical. The other two bands, instead, lead to an EOS that is compatible with calculations using phenomenological Hamiltonians, but provide uncertainty estimates.

To describe NS, PNM calculations have to be extended to both $\beta$ equilibrium as well as to higher densities. While nuclear Hamiltonians have been used in QMC calculations at all densities encountered in NSs [106], it is not clear if a description in terms of nucleonic degrees of freedom remains valid at high densities. Therefore, a more conservative approach is to use results based on realistic Hamiltonians at small densities, where uncertainties are under control, and extrapolate to higher densities using general extrapolation schemes, e.g. polytropic extensions [107, 108] or speed-of-sound extensions [8, 109]. With the EOS specified, the structure of an idealized spherically symmetric NS can be calculated by integrating the Tolman-Oppenheimer-Volkoff equations.

Such a general extrapolation has been performed for the PNM EOS from AFDMC calculations in References [8, 105]. In the left panel of Figure 7 we show the resulting mass-radius uncertainty when the QMC input is used up to $2n_{0}$, for smooth EOS models based on an expansion of nuclear matter around saturation density [minimal model (MM), black band] and for EOS models that are based on a general extension of the speed of sound and also allow phase transitions [speed-of-sound model (CSM), red band]. As one can see, the radius of a typical neutron star has an uncertainty of $\approx 4$ km. While future measurements of neutron star radii from NICER [110] or eXTP [111] might strongly constrain the EOS, currently the strongest constraints come from NS mass measurements. Observations in the past eight years [112, 113] have found two NSs with masses near $2M_{\odot}$. These two observations provide some of the strongest constraints on the nature of the EOS above the nuclear saturation density and have been taken into account in Figure 7.

Finally, we address the recently observed NS merger, GW170817 [114, 115]. The gravitational-wave signal from this event can be used to constrain the tidal polarizabilities $\Lambda_{i}$ of the two NSs in the binary system. The tidal polarizability measures to what extent an NS deforms under an external gravitational field, and depends on the compactness of the NS. Using the previously defined EOS, one can compute the tidal polarizabilities of the two NSs in GW170817 and compare the results to the observation by the LIGO-Virgo collaboration [105]. We show the results in the right panel of Figure 7 and compare the two extrapolations to the observation (blue). One finds that nuclear physics calculations up to $2n_{0}$ are more constraining for tidal polarizabilities than the observation of GW170817 with its uncertainty, and that GW170817 does not improve our understanding of the EOS.

7 RESULTS IN ELECTROWEAK REACTIONS

The description of neutrino interactions with nuclei provides an essential input for current and planned accelerator-neutrino experiments [116, 117, 118, 119, 120]. Since the energy of neutrino beams, produced as secondary decay products, is not monochromatic, neutrino-oscillation experiments are sensitive to a variety of reaction mechanisms, whose contributions depend on the energy and momentum transfer. The low energy-transfer regime is dominated by coherent scattering, excitations of low-lying nuclear states, and collective modes. At energies of the order of hundreds of MeVs, the leading mechanism is quasielastic scattering, in which the probe interacts primarily with individual bound nucleons. Corrections to this leading mechanism arise from processes in which the lepton couples to pairs of interacting nucleons. At higher energies, neutrinos can also excite a struck nucleon to a baryon resonance state that quickly decays into pions. In this regime, a description solely based on nucleonic degrees of freedom starts to fail, as deep-inelastic scattering (DIS) processes need to be accounted for.{marginnote}[]\entryDISdeep inelastic scattering Achieving a unified description of all these reaction mechanisms is a formidable nuclear theory challenge. Valuable information can be inferred from the analysis of the large wealth of available electron-nucleus scattering data [121]. The experiments recently carried out at JLab [122] are particularly relevant in this regard, as they probe phenomena occurring at small internuclear distances, such as the reduction in the DIS cross-section ratios for heavier nuclei relative to deuterium (EMC effect), and its connection to NN SRCs in nuclei [98, 123]; see also Section 5.2.

The interactions of an external electroweak probe with a nucleus are described by the response functions, encoding the strong-interaction dynamics of the nucleons, and their coupling to these external fields. The response functions – two for the electromagnetic processes, and five for the neutral or charge-changing weak processes – can be schematically written as

[TABLE]

where ${\bf q}$ and $\omega$ are the momentum and energy transfers injected by the external field into the nucleus, $|0\rangle$ and $|f\rangle$ represent respectively its initial ground state of energy $E_{0}$ and final states of energy $E_{f}$ (possibly in the continuum), and $J_{\alpha}$ denotes the appropriate components of the nuclear electroweak current operator [50, 124].

Even at intermediate values of the momentum transfer ($q\lesssim 0.5$ GeV) and for energy-transfer corresponding to the quasielastic region, the calculation of the response functions involves severe difficulties, as it requires summation over the entire excitation spectrum of the nucleus and the inclusion of one- and many-body terms in the electroweak currents.

Integral properties of the responses can be studied by means of their sum rules

[TABLE]

where $C_{\alpha\beta}$ are $q$-dependent normalization factors. Fixing the $\omega$-dependence of the current operators at the quasielastic peak, $\omega_{\text{qe}}=\sqrt{q^{2}+m^{2}}-m$, the sum rules can be expressed as ground-state expectation values $S_{\alpha\beta}(q)=\expectationvalue{J_{\alpha}^{\dagger}(\omega_{\text{qe}},\mathbf{q})J_{\beta}(\omega_{\text{qe}},\mathbf{q})}{0}$.

GFMC calculations of the sum rules of the electromagnetic and neutral-current response functions of \isotope[12]C have been reported in References [62, 125]. Processes involving two-body currents substantially increase ($\simeq 30$%) the one-body sum rules even down to small momentum transfers. At low momentum transfers terms involving two-body currents only dominate, primarily with the same pair contributing in both $J_{\alpha}^{\dagger}$ and $J_{\beta}$. At higher momentum transfers the interference between one- and two-nucleon currents plays a more important role, as noted in Reference [126]. Consistently with the analysis carried out by experiments at Brookhaven National Laboratory [127] and JLab [100, 102] on exclusive measurements of back-to-back pairs in \isotope[12]C, the contribution of np pairs is the most important one. This has to be ascribed to the tensor component of the nuclear interaction, which plays a larger role in np pairs where it can act in relative $S$-waves, while it acts only in relative $P$-waves (and higher partial waves) in nn and pp pairs [128, 129, 130, 131].

The sum rules calculations are not capable of identifying the energy-transfer dependence of the calculated excess strength induced by two-body currents. A direct GFMC calculation of $R_{\alpha\beta}(q,\omega)$ is impractical, because it would require evaluating each individual transition amplitude $\ket{0}\longrightarrow\ket{f}$ induced by the current operators. To circumvent this difficulty, the use of integral transform techniques has proved to be quite helpful [132]. Such an approach is based on the Laplace transform of the response functions, i.e. the Euclidean response [133, 134], defined as

[TABLE]

The lower integration limit $\omega_{\rm el}=q^{2}/2M_{A}$, $M_{A}$ being the mass of the target nucleus, is the elastic scattering threshold – corresponding to the $\ket{f}=\ket{0}$ term in the sum of Equation 49 – whose contribution is excluded. Using the same procedure as in the sum-rule calculations to fix the $\omega$ dependence of the current operators, the Euclidean responses can be expressed as ground-state expectation values,

[TABLE]

The Euclidean response functions reported in References [63, 64, 135] are computed with the variational wave function, $\ket{0}=\ket{\Psi_{V}}$. This is justified by the fact that the sum rules computed with $\ket{\Psi_{V}}$ for \isotope[12]C are very close (within less than $5\%$) to those computed with the GFMC wave function [62, 125]. The calculation of the matrix element above proceeds in two steps [136]. First, an unconstrained imaginary-time propagation of the VMC state $\ket{\Psi_{V}}$ is performed and saved. Next, the states $J_{\beta}(\omega_{\text{qe}},\mathbf{q})\ket{\Psi_{V}}$ are evolved in imaginary time following the path previously saved. During this latter imaginary-time evolution, scalar products of $\exp\left[-\left(H-E_{0}\right)\tau_{i}\right]J_{\beta}(\omega_{\text{qe}},\mathbf{q})\ket{\Psi_{V}}$ with $J_{\alpha}(\omega_{\text{qe}},\mathbf{q})\ket{\Psi_{V}}$ are evaluated on a grid of $\tau_{i}$ values, and from these scalar products estimates for $E_{\alpha\beta}(q,\tau_{i})$ and for the associated statistical error are obtained [133, 134].

Retrieving the energy dependence of the response functions requires a numerical inversion of the Laplace transform of Equation 51, a notoriously ill-posed problem. The GFMC calculations of the electroweak response functions carried out over the last few years exploit maximum entropy techniques [137, 138] to perform the analytic continuation of the Euclidean response function. More specifically, the so called “historic maximum entropy” technique, employed in Reference [63], has been augmented to better propagate the statistical errors associated with $E_{\alpha\beta}(q,\tau)$. By exploiting GFMC and maximum entropy techniques, the authors of Reference [64] have demonstrated that accurate calculations of the response, based on a realistic correlated nuclear wave function and containing one- and two-body currents, can reproduce the \isotope[12]C electromagnetic response functions in the quasielastic region. In the top two panels of Figure 8, the GFMC response functions of \isotope[12]C at $q=570$ MeV in which only one-body or both one- and two-body terms are included in the electromagnetic current operators – denoted by (red) dashed and (black) solid lines and labeled GFMC-J1b and GFMC-J1b+2b, respectively – are compared to the experimental world data analysis of Reference [139]. The red and gray shaded areas show the uncertainty of the inversion procedure, ultimately associated with the statistical error of the corresponding Euclidean responses. While the contributions from two-body charge operators tend to slightly reduce the longitudinal response in the threshold region, those from two-body currents generate a large excess of strength in the transverse channel, significantly improving the agreement with experimental data. The absence of explicit pion production mechanisms restricts the applicability of the GFMC method to the quasielastic region of the transverse response. Within this picture, the so-called quenching of the longitudinal response near the quasielastic peak emerges as a result of initial- state correlations and final-state interactions, as opposed to the in-medium modification of the nucleon form factors advocated in Reference [140].

The $\nu$ and $\overline{\nu}$ differential cross sections and the $\nu/\overline{\nu}$ ratios for a fixed value of the three-momentum transfer ($q\,$= 570 MeV/c) as function of the energy transfer for a number of scattering angles are displayed in the bottom two panels of Figure 8. In terms of the response functions, they are given by

[TABLE]

where $-$ ($+$) refers to $\nu$ ($\overline{\nu}$), $k^{\prime}$ and $E^{\prime}$ are the momentum and energy of the outgoing neutrino, $q$ and $\omega$ are the momentum and energy transfers with $Q^{2}=q^{2}-\omega^{2}$ being the four-momentum transfer, $\theta$ is the outgoing neutrino scattering angle relative to the direction of the incident neutrino beam, and $G_{F}=1.1803\times 10^{-5}$ GeV*-2* [141]. GFMC results for the response functions and cross sections relevant to neutrino scattering off 12C induced by neutral-current transitions are reported in Reference [135]. These calculations are based on the same dynamical model employed in the electromagnetic case: the nucleons interact with each other via AV18 and IL7 nuclear potentials and with electroweak fields via the phenomenological currents presented is Section 3. Because of the cancellation between the $R_{xx}$ and $R_{xy}$ response functions in Equation 53, the $\overline{\nu}$ cross section decreases rapidly relative to the $\nu$ one as the scattering angle changes from the forward to the backward hemisphere. For analogous reasons, two-body currents, which mostly impact the $R_{xx}$ and $R_{xy}$ responses, are larger for the $\nu$ than for the $\overline{\nu}$ cross section, becoming almost negligible in the latter case for backward-angle cross section. Their contributions significantly increase the magnitude of the cross sections over the entire quasielastic region, and particularly the ratio of neutrino to anti-neutrino cross sections. The analysis of the five response functions entering Equation 53 reveals that this enhancement is mostly due to constructive interference between the one- and two-body current matrix elements, and is consistent with that expected on the basis of sum-rule analyses discussed earlier. It has to be noted that, at variance to the electromagnetic case, two-body terms in the weak neutral charge also produce excess strength in $R_{00}$ and $R_{0z}$ beyond the quasielastic peak.

A major limitation of the GFMC response functions comes from the nonrelativistic nature of the calculation. In Reference [142] the applicability of GFMC has been extended in the quasielastic region to intermediate momentum transfers by performing the calculations in a reference frame that minimizes the momenta of the struck nucleon. Additional relativistic effects in the kinematics are accounted for employing the two-fragment model, which relies on the assumption that the quasielastic reaction is dominated by the break-up of the nucleus into a knocked-out nucleon and a remaining $(A-1)$ system. This assumption enables one to connect, in a relativistically correct way, the energy transfer to the excitation energy of the nucleus, entering the energy-conserving delta function of the scattering process. It has to be noted that the two-fragment model has been adopted only for determining the kinematic input of a calculation where the full nuclear dynamics of the system is taken into account. This is achieved by interpolating the GFMC response function at energy and momentum transfer that fulfill a relativistic energy-conserving delta function. Despite the two-fragment model does not contain tunable parameters, it improves the agreement between experimental data and GFMC calculations for the inclusive electron-4He cross sections, especially for relatively large values of the incoming lepton energy – see Figure 7 of Reference [142].

8 SUMMARY AND OUTLOOK

In this review we have presented recent advances in QMC methods for nuclear physics. Most of these advances have been made possible both by new developments of many-body methods themselves and by the implementation of systematic interactions from chiral EFT. We have presented results showing that this fruitful combination can accurately describe energies of ground and excited states, radii, momentum distributions of nuclei up to \isotope[16]O, and $n$-$\alpha$ scattering. At the same time, QMC methods with chiral interactions give a reliable description of neutron matter for astrophysical applications, such as neutron stars and neutron-star mergers. We have also reviewed exciting results on electroweak reactions, which are important for the understanding of interactions between neutrinos and matter.

Looking forward to the future, there is still interesting and important work to be done. A remaining milestone is the combination of accurate QMC methods with consistent chiral EFT interactions and electroweak currents and the study of heavier systems with this combined approach. To achieve these goals, both QMC methods and local chiral EFT interactions need to be improved.

For the methods, one question that must be addressed is how to build a wave function for AFDMC calculations of nuclei, which is sophisticated enough to capture important correlations for larger nuclear systems such as \isotope[40]Ca, while still maintaining favorable scaling in $A$. In addition, the implementation of different boundary conditions might allow one to access the thermodynamic limit in nuclear matter for smaller particle numbers, which might permit systematic computations of asymmetric matter. On the interaction side, extending minimally nonlocal chiral interactions to higher orders in the EFT power counting, including consistent many-body forces, will be an important step forward. For example, complete chiral interactions at N3LO will help to reduce the systematic uncertainties and allow to determine how well the chiral expansion is converging in the local chiral EFT approach. Within this context, it would be desirable to make a systematic comparison, order by order, of the $\Delta$-less and $\Delta$-full local chiral interactions to better understand the effect of $\Delta$ degrees of freedom on the order-by-order convergence. The effect of regulator artifacts and a wider range of cutoffs also needs to be studied in few- and many-body systems. Finally, explicitly including pion fields in QMC methods [143] might offer new insights on the chiral expansion.

DISCLOSURE STATEMENT

The authors are not aware of any affiliations, memberships, funding, or financial holdings that might be perceived as affecting the objectivity of this review.

ACKNOWLEDGMENTS

We wish to give our special thanks to the late Steven C. Pieper, who contributed so much to the progress of nuclear QMC methods. We also thank J. Carlson, D. Lonardoni, R. Schiavilla, A. Schwenk, and R.B. Wiringa. The work of J.E.L. was supported by the BMBF under contract No 05P15RDFN1. The work of I.T. and S.G. was supported by the U.S. DOE under contract DE-AC52-06NA25396, and by the LANL LDRD program. S.G. was also supported by the NUCLEI SciDAC program and by the DOE Early Career Research Program. The work of A.L. was supported by the U.S. DOE under contract DE-AC02-06CH11357. This research used resources 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. Computational resources have been provided 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. Calculations for this research were also conducted on the Lichtenberg high performance computer of the TU Darmstadt.