Atomic nuclei from quantum Monte Carlo calculations with chiral EFT interactions

The study of nuclear properties as they emerge from the individual interactions among protons and neutrons is a fascinating long-standing problem, subject of both theoretically and experimentally research activities. From a theoretical point of view, a truly ab-initio description of nuclei is still very challenging at present. The underlying theory of strong interactions, Quantum Chromodynamics (QCD), that describes how quarks and gluons interact to form nucleons and nuclei, in the low-energy regime is non-perturbative in its coupling constant. Despite remarkable progresses [1, 2], realistic computations of many-body nuclear systems in terms of the fundamental degrees of freedom of QCD – quarks and gluons – are still extremely challenging.

A more feasible approach to the problem consists in assuming that at the energy regime relevant to the description of atomic nuclei, quarks, and gluons are confined within hadrons. The latter are the active degrees of freedom at soft scales, and they interact among themselves through non-relativistic effective potentials that are consistent with the symmetries of QCD. The solution of nuclear many-body problems requires two main ingredients: an Hamiltonian that accurately models the interactions among the nucleons, and reliable numerical many-body methods to solve the corresponding Schrödinger equation.

Microscopic nuclear Hamiltonians, capable of reproducing nucleon-nucleon scattering data and the properties of few-body systems, have been successfully used to describe light nuclei. For example, the highly-realistic Argonne $v_{18}$ two-body potential [3] combined with the phenomenological Illinois-7 three-body force have been employed to predict several properties of nuclei up to $A=12$ with great accuracy [4]. Several calculations of energies, rms radii, transitions, and densities turn out to be in excellent agreement with experimental data. The main limitation of these phenomenological Hamiltonians is that it is not clear how they can be systematically improved, and how to quantify theoretical, i.e., systematic, uncertainties related to the specific interaction model. Another approach that became very popular in the last two decades consist in deriving nuclear interactions within the framework of chiral Effective Field Theory ($\chi$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 $\chi$EFT Hamiltonians were typically written in momentum space. Over the past few years, the situation has drastically changed with the development of local $\chi$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.

The microscopic model of nuclear theory assumes that nuclear systems can be described as point-like nucleons, whose dynamics is characterized by a non-relativistic Hamiltonian

where $T_{i}$ is the one-body kinetic energy operator, $v_{ij}$ is the nucleon-nucleon ($N\!N$) 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 $\sim 100\,\rm keV$ in \isotope[4]He [14] or pure neutron matter [15], and therefore are negligible compared to $N\!N$ and $3N$ components. Hence, current formulations of the microscopic model do not typically include them (see, for example, reference [4]).

The $N\!N$ interaction term in the nuclear Hamiltonian is the most studied of all, with thousands of experimental data points at laboratory energies ($E_{\rm lab}$) from essentially zero to hundreds of MeV. It consists of a long-range component, for inter-nucleon separation $r\gtrsim 2$ fm, due to one-pion exchange (OPE) [16], and intermediate- and short-range components, for, respectively, $1\,{\rm fm}\lesssim r\lesssim 2\,\rm fm$ and $r\lesssim 1\,\rm fm$, derived, up to the mid 1990’s, almost exclusively from meson-exchange phenomenology [17, 3, 18]. These models fit the large amount of empirical information about $N\!N$ scattering data contained in the Nijmegen database [19], available at the time, with a $\chi^{2}/{\rm datum}\simeq 1$ for $E_{\rm lab}$ up to pion-production threshold. Two well-known and still widely used examples in this class of $N\!N$ interactions are the CD-Bonn [18] and the Argonne $v_{18}$ (AV18) [3] potentials.

The AV18 interaction is a local, configuration-space $N\!N$ potential that has been extensively and successfully used in a number of QMC calculations. It is expressed as a sum of electromagnetic and OPE terms and phenomenological intermediate- and short-range parts:

The electromagnetic term $v^{\gamma}_{ij}$ has one- and two-photon exchange Coulomb interaction, vacuum polarization, Darwin-Foldy, and magnetic moment terms, with appropriate form factors that keep terms finite at $r=0$ (see reference [3] for more details). The OPE part includes the charge-dependent (CD) terms due to the difference in neutral ($m_{\pi_{0}}$) and charged pion ($m_{\pi_{\pm}}$) masses, and in coordinate-space it reads

where ${\bm{\sigma}}$ adn ${\bm{\tau}}$ are the Pauli matrices that operate over the spin and isospin of particles, and $S_{ij}=3\,{\bm{\sigma}}_{i}\cdot\hat{\textbf{r}}_{ij}\,\,{\bm{\sigma}}_{j}\cdot\hat{\textbf{r}}_{ij}-{\bm{\sigma}}_{i}\cdot{\bm{\sigma}}_{j}$ and $T_{ij}=3\,\tau_{iz}\tau_{jz}-{\bm{\tau}}_{i}\cdot{\bm{\tau}}_{j}$ are the tensor and isotensor operators, respectively. The functions, $v^{\pi}_{\sigma\tau}(r)$, $v^{\pi}_{t\tau}(r)$, $v^{\pi,\rm}_{\sigma T}(r)$, and $v^{\pi}_{tT}(r)$ are defined as

where $Y_{\alpha}(r)$ and $T_{\alpha}(r)$ are the Yukawa and tensor functions given by

with $x_{\alpha}=m_{\pi_{\alpha}}r$, and $g_{A}=1.267$, $f_{\pi}=92.4\,\rm MeV$ being the axial-vector coupling constant of the nucleon and the pion decay constant, respectively.

The intermediate-range region, $v^{I}_{ij}$, is parametrized in terms of two-pion exchange (TPE), based on, but not consistently derived from, a field-theory analysis of box diagrams with intermediate nucleons and $\Delta$ isobars [20]. The short-range region, $v^{S}_{ij}$, is instead represented by spin-isospin and momentum-dependent operators multiplied by Woods-Saxon radial functions [3].

The AV18 model can be written as an overall sum of eighteen operators $\left(\mathcal{N}=18\right)$

where the first eight are given by

with the spin-orbit contribution expressed in terms of the relative angular momentum $\mathbf{L}=\frac{1}{2i}(\mathbf{r}_{i}-\mathbf{r}_{j})\times(\bm{\nabla}_{i}-\bm{\nabla}_{j})$ and the total spin $\mathbf{S}=\frac{1}{2}({\bm{\sigma}}_{i}+{\bm{\sigma}}_{j})$ of the pair. There are six additional charge-independent operators corresponding to $p=9-14$ that are quadratic in $L$

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

The AV18 model has a total of 42 independent parameters. A simplex routine [21] was used to make an initial fit to the phase shifts of the Nijmegen partial-wave analysis (PWA) [19], followed by a final fit direct to the database, which contains 1787 $pp$ and 2514 $np$ observables for $E_{\rm lab}$ up to $350\,\rm MeV$. The $nn$ scattering length and deuteron binding energy were also fit. The final $\chi^{2}/{\rm datum}=1.1$ [3]. While the fit was made up to $350\,\rm MeV$, the phase shifts are qualitatively good up to much larger energies, $E\leq 600\,\rm 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, $\mathcal{O}^{p=1-8}_{ij}=\big{[}1,{\bm{\sigma}}_{i}\cdot{\bm{\sigma}}_{j},S_{ij},\mathbf{L}\cdot\mathbf{S}\big{]}\otimes\big{[}1,{\bm{\tau}}_{i}\cdot{\bm{\tau}}_{j}\big{]}$, of the full $N\!N$ potential, constructed to preserve the potential in all $S$ and $P$ waves as well as the ${}^{3}D_{1}$ and its coupling to the ${}^{3}S_{1}$, while over-binding the deuteron by $18\,\rm 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 $N\!N$ interactions did not provide enough binding for the three-body nuclei, \isotope[3]H and \isotope[3]He [24]. In the late 1990s and early 2000s this realization was also extended to the spectra (ground and low-lying excited states) of light $p$-shell nuclei, for instance, in calculations based on QMC methods [25] and in no-core shell-model studies [26]. Consequently, the microscopic model with only $N\!N$ interactions fit to scattering data, without the inclusion of a $3N$ interaction, is no longer considered realistic.

In addition to $N\!N$ 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 \isotope[3]H and the saturation-point of symmetric nuclear matter, while the IL7 potential involves five parameters constrained on the low-lying spectra of nuclei in the mass range $A=3-10$.

Despite their success in predicting a wide range of nuclear properties [4], the phenomenological potentials suffer from several drawbacks. For example, the resulting AV18+IL7 Hamiltonian leads to predictions of $\approx 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 $N\!N$ 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 $\chi$EFT [30, 31, 32]. $\chi$EFT is a low-energy effective theory of QCD and provides the most general scheme accommodating all possible interactions among nucleons and pions ($\Delta$-less $\chi$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 $\Delta$ isobar ($\Delta$-full $\chi$EFT), because the $\Delta$-nucleon mass splitting is only $300\,\rm MeV\sim 2m_{\pi}$.

By its own nature, the $\chi$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 $\Lambda_{b}\sim\rm GeV$, $Q=p/\Lambda_{b}$, where the breakdown scale denotes momenta at which the short distance structure becomes important and cannot be neglected and absorbed into contact interactions anymore (see references [33, 34, 35, 36] for recent review articles). It is important mentioning that alternative power-counting schemes have been also suggested [37, 38, 39, 40, 41, 42] but not fully explored.

This expansion introduces an order by order scheme, defined by the power $\nu$ of the expansion scale $Q$ associated with each interaction terms: leading order (LO) for $\nu=0$, next-to-leading order (NLO) for $\nu=2$, next-to-next-to-leading order (N²LO) for $\nu=3$ and so on. Similarly as for nuclear interactions, such a scheme can also be developed for electroweak currents. Therefore, $\chi$EFT provides a rigorous scheme to systematically construct many-body forces and consistent electroweak currents, and tools to estimate their uncertainties [43, 44, 45, 46, 47, 48].

Chiral nuclear forces are comprised of both pion-exchange contributions and contact terms. Pion-exchange contributions represent the long-range part of nuclear interactions and some of the pion-nucleon ($\pi N$) couplings entering sub-leading diagrams can consistently be obtained from low-energy $\pi N$ scattering data [49, 50, 51, 52, 53, 54, 55, 56]. On the other end, contact terms encode the short-range physics, and their strength is specified by unknown low-energy constants (LECs), obtained by fitting experimental data. Similarly to the phenomenological interactions, the LECs entering the $N\!N$ component are obtained by fitting $N\!N$ scattering data up to $300\,\rm 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 $N\!N$ 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 $N\!N$ and $3N$ forces to low-energy $N\!N$ data, deuteron binding energy, and binding energies and charge radii of hydrogen, helium, carbon, and oxygen isotopes.

Within $\chi$EFT, many studies have been carried out dealing with the construction and optimization of $N\!N$ and $3N$ interactions [58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 34, 49, 42, 69, 70, 71, 72, 73, 44, 74, 75, 76, 77, 57, 78, 79, 80, 81] and accompanying isospin-symmetry-breaking corrections [82, 83, 84]. These interactions are typically formulated in momentum space, and include cutoff functions to regularize their behavior at large momenta. This makes them strongly non-local when Fourier-transformed in configuration space, and therefore unsuitable for use with QMC methods. In this context, an interaction is local if it depends solely on the momentum transfer $\textbf{q}={\textbf{p}}-{\textbf{p}}^{\prime}$ (p and ${\textbf{p}}^{\prime}$ 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 ${\textbf{k}}=({\textbf{p}}^{\prime}+{\textbf{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 $\textbf{p}^{\prime}$, and (ii) contact interactions that explicitly depend on k.

In recent years, local configuration-space chiral $N\!N$ interactions have been derived by two groups. On the one side, the authors of references [12, 85] constructed $N\!N$ local chiral potentials within $\Delta$-less $\chi$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 $\chi^{2}$ fit to $N\!N$ phase shifts from the Nijmegen PWA up to $150\,\rm MeV$ lab energy. The fitting procedure is carried out for different values of the cutoff $R_{0}$ in the range of $R_{0}=1.0-1.2\,\rm 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 $N\!N$ 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 \isotope[4]He as well as the $P$-wave $n-\alpha$ elastic scattering phase shifts. In addition, they explore different parametrization for the $3N$ force, accordingly to Fierz identities [87, 88, 89]. In the present work, we are referring to a set of these local chiral interactions, specifically the $(D2,E\tau)$ model with $R_{0}=1.0\,\rm fm$ of reference [89], as GT+E$\tau$-1.0.

On the other side, the authors of references [10, 90] developed a different set of $N\!N$ local chiral interactions by i) including diagrams with the virtual excitation of $\Delta$-isobars in the TPE contributions up to N²LO ($\Delta$-full $\chi$EFT), ii) retaining contact terms up to N³LO. The LECs entering the $N\!N$ contact interactions in this model are constrained to reproduce $N\!N$ scattering data from the most recent and up-to-date database collected by the Granada group [91, 92, 93]. The contact terms are implemented via a Gaussian representation of the three-dimensional delta function with $R_{S}$ as the Gaussian parameter [90, 10, 94]. The pion-range operators are regularized at high-value of momentum transfer via a special radial function characterized by the cutoff $R_{L}$ [90, 10, 94]. There are two classes of these potentials. Class I (II) are fit to data up to $125\,\rm MeV$ ($200\,\rm MeV$). For each class, two combinations of short- and long-range regulators have been used, namely $(R_{S},R_{L})=(0.8,1.2)\,\rm fm$ (models NV2-Ia and NV2-IIa) and $(R_{S},R_{L})=(0.7,1.0)\,\rm fm$ (models NV2-Ib and NV2-IIb). Class I (II) fits about 2700 (3700) data points with a $\chi^{2}/{\rm datum}\lesssim 1.1$ ($\lesssim 1.4$) [90, 10]. In conjunction with these models, two distinct sets of $\Delta$-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 $\beta$-decay [94]. The resulting Hamiltonians were labelled as NV2+3-Ia/b and NV2+3-IIa/b (or Ia/b and IIa/b for short) in the first case, and as NV2+3-Ia^∗/b^∗ and NV2+3-IIa^∗/b^∗ (or Ia^∗/b^∗ and IIa^∗/b^∗) in the second.

The interactions between external electroweak probes – electrons and neutrinos – and interacting nuclear systems is described by a set of effective nuclear currents and charge operators. Analogously to the nuclear interactions, electroweak currents can also be expressed as an expansion in many-body operators that act on nucleonic degrees of freedom. Electroweak currents have been developed in both meson-exchange and $\chi$EFT approaches. We refrain to discuss them in this work, redirecting the interested reader to dedicated reviews [4, 96, 97, 98] and references therein.

The $\chi$EFT Hamiltonians and the consistent electroweak currents discussed in the previous Section are the main input of sophisticated many-body methods aimed at solving with controlled approximations the nuclear many-body Schrödinger equation

This is a highly non-trivial problem, mainly because of the non-perturbative nature and the strong spin-isospin dependence of realistic nuclear forces. In this work, we will focus on QMC techniques, namely the variational Monte Carlo (VMC), the Green’s function Monte Carlo (GFMC), and the auxiliary-field diffusion Monte Carlo (AFDMC) methods.