Quantum Monte Carlo methods for nuclear physics
Abstract
Quantum Monte Carlo methods have proved very valuable to study the structure and reactions of light nuclei and nucleonic matter starting from realistic nuclear interactions and currents. These abinitio calculations reproduce many lowlying states, moments and transitions in light nuclei, and simultaneously predict many properties of light nuclei and neutron matter over a rather wide range of energy and momenta. We review the nuclear interactions and currents, and describe the continuum Quantum Monte Carlo methods used in nuclear physics. These methods are similar to those used in condensed matter and electronic structure but naturally include spinisospin, tensor, spinorbit, and threebody interactions. We present a variety of results including the lowlying spectra of light nuclei, nuclear form factors, and transition matrix elements. We also describe lowenergy scattering techniques, studies of the electroweak response of nuclei relevant in electron and neutrino scattering, and the properties of dense nucleonic matter as found in neutron stars. A coherent picture of nuclear structure and dynamics emerges based upon rather simple but realistic interactions and currents.
Contents:
 I INTRODUCTION
 II HAMILTONIAN
 III Quantum Monte Carlo methods
 IV Light nuclei

V Electroweak currents
 V.1 Conventional nuclear electroweak currents
 V.2 Electromagnetic currents in chiral effective field theory
 V.3 Elastic and inelastic form factors
 V.4 Second state of C: Hoyle state
 V.5 Magnetic moments and electroweak transitions
 V.6 Electroweak Response of Light Nuclei
 V.7 Sum rules of electroweak response functions
 V.8 Euclidean response functions
 VI The equation of state of neutron matter
 VII Conclusions and Future Directions
I Introduction
Nuclei are fascinating few and manybody quantum systems, ranging in size from the lightest nuclei formed in the big bang to the structure of neutron stars with 10 km radii. Understanding their structure and dynamics starting from realistic interactions among nucleons has been a longstanding goal of nuclear physics. The nuclear quantum manybody problem contains many features present in other areas such as condensed matter physics, including pairing and superfluidity and shell structure, but also others that are less common including a very strong coupling of spin and spatial degrees of freedom, clustering phenomena, and strong spinorbit splittings. The challenge is to describe diverse physical phenomenon within a single coherent picture.
This understanding is clearly important to describe nuclear properties and reactions, including reactions that synthesized the elements and the structure of neutronrich nuclei. An accurate picture of interactions and currents at the nucleonic level is critical to extend this understanding to the properties of dense nucleonic matter as occurs in neutron stars, and to use nuclei as probes of fundamental physics through, for example, beta decay, neutrinoless doublebeta decay, and neutrinonucleus scattering.
Over the last three decades it has become possible using Quantum Monte Carlo (QMC) methods to reliably compute the properties of light nuclei and neutron matter starting from realistic nuclear interactions. While many of the most basic properties of nuclei can be obtained from comparatively simple meanfield models, it has been a challenge to relate the two and threenucleon interactions inferred from experiments to the structure and reactions of nuclei. This challenge arises because the scale of the nuclear interactions obtained by examining nucleonnucleon phase shifts is of order (50100) MeV or more, significantly larger than a typical nuclear binding energy of 8 MeV per nucleon.
In addition, the nucleonnucleon interaction is much more complex than the Coulomb force used in molecular and atomic physics, the van der Waals potential between atoms used, for example, in studies of liquid helium systems, or the contact interaction that dominates dilute coldatom physics. The primary force carrier at large nucleon separations is the pion, which couples strongly to both the spin and isospin of the nucleons with a strong tensor component. In addition there are significant spinorbit forces. As a consequence, there is strong coupling between the spin and isospin and spatial degrees of freedom.
These features lead to complex nuclear phenomena. The interactions are predominantly attractive at low momenta, resulting in large pairing gaps in nuclei and associated superfluidity in matter. In light nuclei, there is further clustering of neutrons and protons into alphaparticle like configurations that are very evident in the lowlying excitations of some nuclei. At moderate nucleon separations, the tensor character of the neutronproton interaction produces significant highmomentum components in the nuclear wave function that impact the electroweak response observed in electron and neutrino scattering. The nuclear correlations also significantly quench the singleparticle description of nucleon knockout and transfer reactions. A major challenge has been to include both the shortrange highmomentum phenomena and the longrange superfluid and clustering properties of nuclei and matter in a consistent framework.
QMC methods based upon Feynman path integrals formulated in the continuum have proved to be very valuable in attacking these problems. The sampling of configuration space in variational (VMC) and Green’s function (GFMC) Monte Carlo simulations gives access to many of the important properties of light nuclei including spectra, form factors, transitions, lowenergy scattering and response. Auxiliary Field Diffusion Monte Carlo (AFDMC) uses Monte Carlo to also sample the spinisospin degrees of freedom, enabling studies of, for example, neutron matter that is so critical to determining the structure of neutron stars. In this review we concentrate on continuum Monte Carlo methods. Lattice QMC methods have also recently been employed to study both neutron matter Muller et al. (2000); Lee and Schäfer (2006); Seki and van Kolck (2006); Abe and Seki (2009); Wlazłowski et al. (2014); Roggero et al. (2014) and certain nuclei Lee (2009); Epelbaum et al. (2012). Other Monte Carlo methods combined with the use of effective interactions and/or space models like the shell model have been also developed to study properties of larger systems; see for example Koonin et al. (1997); BonettMatiz et al. (2013); Otsuka et al. (2001); Abe et al. (2012); Bonnard and Juillet (2013).
Other manybody methods, many of which have direct analogues in other fields of physics, have also played important roles in the study of nuclei. These include the coupled cluster method Hagen et al. (2014, 2014), the no core shell model Barrett et al. (2013), the similarity renormalization group Bogner et al. (2010), and the Self Consistent Green’s Function Dickhoff and Barbieri (2004). Each of these methods has distinct advantages, and many are able to treat a wider variety of nuclear interaction models. Quantum Monte Carlo methods, in contrast, are more able to deal with a wider range of momentum and energy and to treat diverse phenomenon including superfluidity and clustering.
Progress has been enabled by simultaneous advances in the input nuclear interactions and currents, the QMC methods, increasingly powerful computer facilities, and the applied mathematics and computer science required to run efficiently these calculations on the largest available machines Lusk et al. (2010). Each of these factors have been very important. QMC methods have been able to make use of some of the most powerful computers available, through extended efforts of physicists and computer scientists to scale the algorithms successfully. The codes have become much more efficient and also more accurate through algorithmic developments. The introduction of Auxiliary Field methods paved the way to scale these results to much larger nuclear systems than would otherwise have been possible. Equally important, advances in algorithms have allowed to expand the physics scope of our investigations. Initial applications were to nuclear ground states, including energies and elastic form factors. Later advances opened the way to study lowenergy nuclear reactions, the electroweak response of nuclei and infinite matter.
Combined, QMC and other computational methods in nuclear physics have allowed us, for the first time, to directly connect the underlying microscopic nuclear interactions and currents with the structure and reactions of nuclei. Nuclear wave functions that contain the manynucleon correlations induced by these interactions are essential for accurate predictions of many experiments. QMC applications in nuclear physics span a wide range of topics, including lowenergy nuclear spectra and transitions, lowenergy reactions of astrophysical interest, tests of fundamental symmetries, electron and neutrinonucleus scattering, and the properties of dense matter as found in neutron stars. In this review we briefly present the interactions and currents and the Monte Carlo methods, and then review results that have been obtained to date across these different diverse and important areas of nuclear physics.
Ii Hamiltonian
Over a substantial range of energy and momenta the structure and reactions of nuclei and nucleonic matter can be studied with a nonrelativistic Hamiltonian with nucleons as the only active degrees of freedom. Typical nuclear binding energies are of order 10 MeV per nucleon and Fermi momenta are around 1.35 fm. Even allowing for substantial correlations beyond the mean field, the nucleons are essentially nonrelativistic. There is a wealth of nucleonnucleon () scattering data available that severely constrains possible interaction models. Nuclear interactions have been obtained that provide accurate fits to these data, both in phenomenological models and in chiral effective field theory. This is not sufficient to reproduce nuclear binding, however, as internal excitations of the nucleon do have some impact. The lowest nucleon excitation is the resonance at 290 MeV. Rather than treat these excitations as dynamical degrees of freedom, however, it is more typical to include them and other effects as threenucleon () interactions.
Therefore, in leading order approximation, one can integrate out nucleon excitations and other degrees of freedom resulting in a Hamiltonian of the form
(1) 
where is the kinetic energy and is an effective interaction, which, in principle, includes nucleon potentials, with :
(2) 
The interaction term is the most studied of all, with thousands of experimental data points at laboratory energies from essentially zero to hundreds of MeV. Attempts are now being made to understand this interaction directly through lattice QCD, though much more development will be required before it can be used directly in studies of nuclei Beane et al. (2013); Ishii et al. (2007). Traditionally the scattering data has been fit with phenomenological interactions that require a rather complicated spinisospin structure because of the way the nucleon couples to the pion, other heavier mesons, and nucleon resonances. More recently, advances have been made using chiral effective field theory, which employs chiral symmetry and a set of lowenergy constants to fit the scattering data. This has led to an understanding of why chargeindependent terms are larger than isospinbreaking ones, why interactions are a small fraction () of interactions, and has provided a direct link between interactions and currents.
In what follows we will focus on potentials developed in coordinate space, which are particularly convenient for QMC calculations. Many phenomenological models are primarily local interactions (although often specified differently in each partial wave) and local interactions can be obtained within chiral effective theory, which is an expansion in the nucleon’s momentum. The interaction is predominantly local because of the nature of onepion exchange, but at higher orders derivative (momentum) operators must be introduced. Local interactions are simpler to treat in continuum QMC methods because the propagator is essentially positive definite, a property that is not always true in nonlocal interactions. The Monte Carlo sampling for such positive definite propagators is much easier, reducing statistical errors in the simulation.
A number of very accurate potentials constructed in the 1990s reproduce the longrange onepionexchange part of the interaction and fit the large amount of empirical information about scattering data contained in the Nijmegen database Stoks et al. (1993b) with a for lab energies up to 350 MeV. These include the potentials of the Nijmegen group Stoks et al. (1994), the Argonne potentials Wiringa et al. (1995); Wiringa and Pieper (2002) and the CDBonn potentials Machleidt et al. (1996); Machleidt (2001). Of those potentials derived more recently by using chiral effective field theory, the most commonly used is that of Entem and Machleidt (2002). The most practical choice for QMC calculations is the Argonne potential Wiringa et al. (1995), which is given in an rspace operator (nonpartial wave) format and has a very weak dependence on nonlocal terms. The latter are small and hence are tractable in QMC calculations. Another less sophisticated interaction that, apart from chargesymmetry breaking effects, reproduces the gross features of Argonne is the Argonne . These are the potentials adopted in most of the QMC calculations.
However all of these interactions, when used alone, underestimate the triton binding energy, indicating that at least forces are necessary to reproduce the physics of H and He. A number of semiphenomenological potentials, such as the Urbana Carlson et al. (1983); Pudliner et al. (1996) series, were developed to fit three and fourbody nuclear ground states. The more recent Illinois Pieper et al. (2001); Pieper (2008a) potentials reproduce the ground state and lowenergy excitations of light shell nuclei (). More sophisticated models may be required to treat nucleonic matter at and above saturation density . Particularly in isospinsymmetric nuclear matter, the manybody techniques for realistic interactions also need to be improved. Effective field theory techniques and QMC methods may help to provide answers to these questions.
ii.1 The nucleonnucleon interaction
Among the realistic interactions, the Argonne (AV18) potential Wiringa et al. (1995) is a finite, local, configurationspace potential that is defined in all partial waves. AV18 has explicit chargeindependence breaking (CIB) terms, so it should be used with a kinetic energy operator that keeps track of the protonneutron mass difference by a split into chargeindependent (CI) and chargesymmetry breaking (CSB) pieces:
(3)  
where and are the proton and neutron mass, and is the operator that selects the third component of the isospin. AV18 is expressed as a sum of electromagnetic and onepionexchange (OPE) terms and phenomenological intermediate and shortrange parts, which can be written as an overall operator sum
(4) 
The electromagnetic term has one and twophotonexchange Coulomb interaction, vacuum polarization, DarwinFoldy, and magnetic moment terms, with appropriate form factors that keep terms finite at =0. The OPE part includes the chargedependent (CD) terms due to the difference in neutral and charged pion masses:
(5) 
where the coupling constant is , are the Pauli matrices that operate over the isospin of particles, and is the isotensor operator. The radial functions are
(6)  
(7)  
(8) 
where or , , the scaling mass , are Pauli matrices that operate over the spin of nucleons, and is the tensor operator. The and are the normal Yukawa and tensor functions with a shortrange cutoff with fm.
The intermediate and shortrange stronginteraction terms have eighteen operators and are given the functional forms
(9)  
(10) 
where is constructed with the average pion mass, , and is a WoodsSaxon potential with radius fm and diffuseness fm. Thus the former has twopionexchange (TPE) range, while the shortrange part remains finite and is constrained to have zero slope at the origin, except for tensor terms which vanish at the origin. The first fourteen operators are CI terms:
(11) 
where is the relative angular momentum of the pair , and is the total spin. The remaining operators include three CD and one CSB terms:
(12)  
(13) 
The maximum value of the central (=1) potential is GeV.
The AV18 model has a total of 42 independent parameters , , and . A simplex routine Nelder and Mead (1965) was used to make an initial fit to the phase shifts of the Nijmegen PWA93 analysis Stoks et al. (1993a), followed by a final fit direct to the data base, which contains 1787 and 2514 observables for MeV. The scattering length and deuteron binding energy were also fit. The final = 1.1 Wiringa et al. (1995). While the fit was made up to 350 MeV, the phase shifts are qualitatively good up to much larger energies MeV Gandolfi et al. (2014).
The CD and CSB terms are small, but there is clear evidence for their presence. The CD terms are constrained by the longrange OPE form and the differences between and scattering in the channel. The CSB term is shortranged and constrained by the difference in and scattering lengths, and is necessary to obtain the correct He–H mass difference.
Direct GFMC and AFDMC calculations with the full AV18 potential are not practical because the spinisospindependent terms which involve the square of the orbital momentum operator have very large statistical errors. However, these terms in AV18 are fairly weak and can be treated as a firstorder perturbation. Using a wave function of good isospin also significantly reduces the cost of calculations in GFMC. Hence it is useful to define a simpler isoscalar AV8 potential with only the first eight (central, spin, isospin, tensor and spinorbit) operators of Eq. (II.1); details are given in Pudliner et al. (1997); Wiringa and Pieper (2002). The AV8 is not a simple truncation of AV18, but a reprojection that preserves the isoscalar average of the strong interaction in all and partial waves as well as the deuteron. It has been used in benchmark calculations of He by seven different manybody methods, including GFMC Kamada et al. (2001).
It has proved useful to define even simpler reprojections of AV8, particularly an AV6 potential without spinorbit terms that is adjusted to preserve deuteron binding. The AV6 has the same CI OPE potential as AV8 and preserves deuteron binding and wave and partial wave phase shifts, but partial waves are no longer properly differentiated. Details are given in Wiringa and Pieper (2002), where the evolution of nuclear spectra with increasing realism of the potentials was investigated.
ii.2 Threebody forces
The Urbana series of potentials Carlson et al. (1983) is written as a sum of twopionexchange wave and remaining shorterrange phenomenological terms,
(14) 
The structure of the twopion wave exchange term with an intermediate excitation (Fig. 1a) was originally written down by Fujita and Miyazawa (1957); it can be expressed simply as
(15) 
where is constructed with the average pion mass and is a sum over the three cyclic exchanges of nucleons . For the Urbana models , as in the original FujitaMiyazawa model, while other potentials like the TucsonMelbourne Coon et al. (1979) and Brazil Coelho et al. (1983) models, have a ratio slightly larger than . The shorterrange phenomenological term is given by
(16) 
For the Urbana IX (UIX) model Pudliner et al. (1995), the two parameters and were determined by fitting the binding energy of H and the density of nuclear matter in conjunction with AV18.
While the combined AV18+UIX Hamiltonian reproduces the binding energies of shell nuclei, it somewhat underbinds light shell nuclei. A particular problem is that the twoparameter Urbana form is not flexible enough to fit both He and Be at the same time. A new class of potentials, called the Illinois models, has been developed to address this problem Pieper et al. (2001). These potentials contain the Urbana terms and two additional terms, resulting in a total of four strength parameters that can be adjusted to fit the data. The general form of the Illinois models is
(17) 
One term, , is due to wave scattering as illustrated in Fig. 1b and is parametrized with a strength . It has been included in a number of potentials like the TucsonMelbourne and Brazil models. The Illinois models use the form recommended in the latest Texas model Friar et al. (1999), where chiral symmetry is used to constrain the structure of the interaction. However, in practice, this term is much smaller than the contribution and behaves similarly in light nuclei, so it is difficult to establish its strength independently just from calculations of energy levels.
A more important addition is a simplified form for threepion rings containing one or two s (Fig. 1c,d). As discussed by Pieper et al. (2001), these diagrams result in a large number of terms, the most important of which are used to construct the Illinois models:
(18) 
Here the and are operators that are symmetric or antisymmetric under any exchange of the three nucleons, and the subscript or indicates that the operators act on, respectively, spin or isospin degrees of freedom.
The is a projector onto isospin triples:
(19) 
To the extent isospin is conserved, there are no such triples in the shell nuclei, and so this term does not affect them. It is also zero for scattering. However, the term is attractive in all the shell nuclei studied. The has the same structure as the isospin part of the anticommutator part of , but the term is repulsive in all nuclei studied so far. In shell nuclei, the magnitude of the term is smaller than that of the term, so the net effect of the is slight repulsion in shell nuclei and larger attraction in shell nuclei. The reader is referred to the appendix of Pieper et al. (2001) for the complete structure of .
The first series of five Illinois models (IL15) explored different combinations of the parameters , , , and , and also variation of the OPE cutoff function . One drawback of these models is that they appear to provide too much attraction in dense neutron matter calculations Sarsa et al. (2003). To help alleviate this problem, the latest version Illinois7 (IL7) introduced an additional repulsive term with the isospin projector:
(20) 
After fixing at the Texas value, and taking from AV18, the four parameters , , , and were searched to obtain a best fit, in conjunction with AV18, for energies of about 20 nuclear ground and lowlying excited states in nuclei Pieper (2008a).
ii.3 Nuclear Hamiltonians from chiral effective field theory
Chiral effective field theory (EFT) has witnessed much progress during the two decades since the pioneering papers by Weinberg (1990, 1991, 1992). In EFT, the symmetries of quantum chromodynamics (QCD), in particular its approximate chiral symmetry, are employed to systematically constrain classes of Lagrangians describing, at low energies, the interactions of baryons (in particular, nucleons and isobars) with pions as well as the interactions of these hadrons with electroweak fields. Each class is characterized by a given power of the pion mass and/or momentum, the latter generically denoted by , and can therefore be thought of as a term in a series expansion in powers of , where GeV specifies the chiralsymmetry breaking scale. Each class also involves a certain number of unknown coefficients, called lowenergy constants (LEC’s), which are determined by fits to experimental data. See, for example, the review papers Bedaque and van Kolck (2002) and Epelbaum et al. (2009), and references therein. Thus EFT provides a direct connection between QCD and its symmetries and the strong and electroweak interactions in nuclei. From this perspective, it can be justifiably argued to have put lowenergy nuclear physics on a more fundamental basis. Just as importantly, it yields a practical calculational scheme, which can, at least in principle, be improved systematically.
Within the nuclear EFT approach, a variety of studies have been carried out in the stronginteraction sector dealing with the derivation of and potentials Ordonez et al. (1996); Epelbaum et al. (1998); Entem and Machleidt (2003); Machleidt and Entem (2011); Navratil (2007); Epelbaum et al. (2002); van Kolck (1994); Bernard et al. (2011); Girlanda et al. (2011) and accompanying isospinsymmetrybreaking corrections Friar and van Kolck (1999); Epelbaum and Meissner (1999); Friar et al. (2004, 2005). In the electroweak sector additional studies have been made dealing with the derivation of parityviolating potentials induced by hadronic weak interactions Haxton and Holstein (2013); Zhu et al. (2005); Girlanda (2008); Viviani et al. (2014) and the construction of nuclear electroweak currents Park et al. (1993); Pastore et al. (2009, 2011); Piarulli et al. (2013); Kölling et al. (2009, 2011).
Recently chiral nuclear interactions have been developed that are local up to nexttonexttoleading order (NLO) Gezerlis et al. (2013). These interactions employ a different regularization scheme from previous chiral interactions, with a cutoff in the relative momentum . They are therefore fairly simple to treat with standard QMC techniques to calculate properties of nuclei and neutron matter Lynn et al. (2014); Gezerlis et al. (2013).
As explained in Gezerlis et al. (2014), up to NLO, the momentumdependent contact interactions can be completely removed by choosing proper local operators. For example, at LO there are several operators that are equivalent for contact interactions: , , , and . Similarly, interactions at NLO and NLO can be constructed by adding extra operators that include the , , and . The shortrange regulators are also chosen to be local, i.e., . In this way, by fitting the lowenergy constants, the chiral potentials are completely local up to NLO. At the next order NLO nonlocal operators start to appear, but their contributions are expected to be very small Piarulli et al. (2015).
Iii Quantum Monte Carlo methods
There is a large variety of Quantum Monte Carlo algorithms, and it would be out of the scope of this review to cover all of them. We will limit ourselves to describing a specific subset of QMC algorithms that has been consistently applied to the many nucleon problem, namely algorithms that are based on a coordinate representation of the Hamiltonian, and that are based on recursive sampling of a probability density or of a propagator. This set of methods includes the standard Variational Monte Carlo (VMC), Green’s Function Monte Carlo (GFMC) and Diffusion Monte Carlo methods.
These methods have been successfully applied to a broad class of problems. The major fields of application of this set of algorithms are quantum chemistry and materials science B.J. Hammond (1994); Nightingale and Umrigar (1999); Foulkes et al. (2001), where QMC is a natural competitor of methods such as Coupled Cluster theory and standard Configuration Interaction methods that are very accurate for problems where the uncorrelated or HartreeFock state provides already a good description of the manybody ground state. In these fields several software packages have been developed with the aim of making the use of QMC methods more and more widespread across the community. Other applications in condensed matter theory concern the physics of condensed helium systems, both He and He Schmidt and Ceperley (1992); Ceperley (1995). Several QMC calculations have been extensively performed to investigate properties of both bosonic and fermionic ultracold gases; see for example Carlson et al. (2003b); Giorgini et al. (2008).
Because of the strong correlations induced by nuclear Hamiltonians, QMC methods have proved to be very valuable in understanding properties of nuclei and nucleonic matter. Variational Monte Carlo methods were introduced for use with nuclear interactions in the early 1980s LomnitzAdler et al. (1981). VMC requires an accurate understanding of the structure of the system to be explored. Typically, a specific class of trial wave functions is considered, and using Monte Carlo quadrature to evaluate the multidimensional integrals, the energy with respect to changes in a set of variational parameters is minimized.
GFMC was introduced in nuclear physics for spin isospindependent Hamiltonians in the late 1980s Carlson (1987, 1988). It involves the projection of the ground state from an initial trial state with an evolution in imaginary time in terms of a path integral, using Monte Carlo techniques to sample the paths. GFMC works best when an accurate trial wave function is available, often developed through initial VMC calculations. This method is very accurate for light nuclei, but becomes increasingly more difficult moving toward larger systems. The growth in computing time is exponential in the number of particles because of the number of spin and isospin states. The largest nuclear GFMC calculations to date are for the C nucleus Lovato et al. (2013, 2014, 2015), and for systems of 16 neutrons Gandolfi et al. (2011); Maris et al. (2013) (540,672 and 65,536 spinisospin states, respectively).
Auxiliary Field Diffusion Monte Carlo (AFDMC) was introduced in 1999 Schmidt and Fantoni (1999). In this algorithm the spin and isospindependence is treated using auxiliary fields. These fields are sampled using Monte Carlo techniques, and the coordinatespace diffusion in GFMC is extended to include a diffusion in the spin and isospin states of the individual nucleons as well. This algorithm is much more efficient at treating large systems. It has been very successful in studying homogeneous and inhomogeneous neutron matter, and recently has been shown to be very promising for calculating properties of heavier nuclei, nuclear matter Gandolfi et al. (2014), and systems including hyperons Lonardoni et al. (2013, 2014, 2015). It does require the use of simpler trial wave functions, though, and is not yet quite as flexible in the complexity of nuclear Hamiltonians that can be employed. Extending the range of interactions that can be treated with AFDMC is an active area of research.
iii.1 Variational Monte Carlo
In VMC, one assumes a form for the trial wave function and optimizes variational parameters, typically by minimizing the energy and/or the variance of the energy with respect to variations in the parameters. The energy of the variational wave function
(21) 
is greater than or equal to the groundstate energy with the same quantum numbers as . Monte Carlo methods can be used to calculate and to minimize the energy with respect to changes in the variational parameters.
For nuclear physics, the trial wave function has the generic form:
(22) 
With this form, a factorization of the wave function into longrange lowmomentum components and shortrange highmomentum components is assumed. The shortrange behavior of the wave functions is controlled by the correlation operator , and the quantum numbers of the system and the longrange behavior by . In nuclei the separation between the shortdistance correlations and the lowmomentum structure of the wave function is less clear than in some systems. For example, alpha particle clusters can be very important in light nuclei, and their structure is of the order of the interparticle spacing. Also the pairing gap can be a nontrivial fraction of the Fermi energy, and hence the coherence length may be smaller than the system. Nevertheless this general form has proved to be extremely useful in both light nuclei and nuclear matter.
Shortrange structure: F
The correlation operator is dominated by Jastrowlike correlations between pairs and triplets of particles:
(23) 
where is the symmetrization operator, is a twobody and is a threebody correlation. The twobody correlation operator can include a strong dependence upon spin and isospin, and is typically taken as:
(24) 
where
(25) 
and the are functions of the distance between particles and . The pair functions are usually obtained as the solution of Schrödingerlike equations in the relative distance between two particles:
(26) 
The pair functions are obtained by solving this equation in different spin and isospin channels, for example , , and can then be recast into operator form. For =1 channels the tensor force enters and this equation becomes two coupled equations for the components with and .
The are functions designed to encode the variational nature of the calculation, mimicking the effect of other particles on the pair in the manybody system. Additional variational choices can be incorporated into boundary conditions on the . For example, in nuclear and neutron matter the pair functions are typically shortranged functions and the boundary condition that and at some distances , which may be different in different channels, is enforced. Usually it is advantageous for the tensor correlation to be finite out to longer distances because of the onepionexchange interaction. The distances are variational parameters, and the equations for the pair correlations are eigenvalue equations; the eigenvalues are contained in the . See Pandharipande and Wiringa (1979) for complete details.
For the lightest shell nuclei (A= 3 and 4), on the other hand, the asymptotic properties of the wave function are encoded in the pair correlation operators . To this end the are determined by requiring the product of pair correlations to have the correct asymptotic behavior as particle is separated from the system. These boundary conditions are described in Schiavilla et al. (1986) and Wiringa (1991).
It has been found advantageous to reduce the strength of the spin and isospindependent pair correlation functions when other particles are nearby, with the simple form above altered to
(27) 
where the central (spinisospin independent) quenching factor is typically 1, while for other operators it is parametrized so as to reduce the pair correlation when another particle is near the pair Pudliner et al. (1997).
The becomes particularly important when the Hamiltonian includes a force. A good correlation form is:
(28) 
with , a scaling parameter, and a (small negative) strength parameter. The superscript denotes various pieces of the force like and , so Eq. (28) brings in all the spinisospin dependence induced by that piece of the potential. In practice the in Eq. (23) is usually replaced with a sum which is significantly faster and results in almost as good a variational energy. For three and fourbody nuclei and nuclear matter, pair spinorbit correlations have also been included in Eq. (23), but they are expensive to compute and not used in the work reviewed here.
The typical number of variational parameters for sshell nuclear wave functions is about two dozen for a twobody potential like AV18, as shown in Wiringa (1991) and Pudliner et al. (1997). Another four to six parameters are added if a threebody potential is included in the Hamiltonian. One can also add a few additional parameters to break charge independence, e.g., to generate components in the trinucleon wave functions, but these are generally used only for studies of isospin violation. For shell nuclei, the alphaparticle pair and triplet correlations are varied only minimally, and most optimization is done with the longrange correlations discussed below.
The variational parameters have generally been optimized by hand. Variational wave functions with significantly larger numbers of parameters and more sophisticated optimization have since been developed Usmani et al. (2009, 2012), but are not in general use. However, they have provided useful insight for improving the simpler parameter sets. The calculation of light nuclei is now sufficiently fast that automated optimization programs might be profitably employed in the future.
LongRange Structure:
The quantum numbers and longrange structure of the wave function are generally controlled by the term in Eq. (22). For nuclear and neutron matter this has often been taken to be an uncorrelated Fermi gas wave function. Recently, the crucial role of superfluidity has been recognized, particularly in lowdensity neutron matter. In such cases the trial wave function includes a of BCS form. For the wave pairing relevant to lowdensity neutron matter, this can be written:
(29) 
where the finite particle number projection of the BCS state has been taken, with the individual pair functions, and the unprimed and primed indices refer to spinup and spindown particles respectively. These pair states are functions of the distance between the two nucleons in the pair. The operator is an antisymmetrization operator Carlson et al. (2003a); Gezerlis and Carlson (2008). For a more general pairing, a Pfaffian wave function is needed (see for example Gandolfi et al. (2008a, 2009a) and references therein).
For light nuclei, the simplest can be written as the sum of a few Slater determinants, essentially those arising from a very modest shellmodel treatment of the nucleus. The singleparticle orbitals in such calculations are written in relative coordinates so as to avoid introducing any spurious centerofmass (CM) motion. An explicit antisymmetrization of the wave function summing over particles in wave, wave, etc., orbitals is required to compute .
Improved wave functions can be obtained by considering the significant cluster structures present in light nuclei. For example the ground state of Be has a very large overlap with two wellseparated alpha particles. Alphacluster structures are important in many light nuclei, for example states in helium and carbon. To this end, it is useful to use a “Jastrow” wave function which includes spinisospin independent two and threebody correlations and the clusterstructure for the :
(30) 
This wave function must be explicitly antisymmetrized as it is written in a particular cluster structure, with particles being in an alphaparticle cluster, summed over the possible partitions. The spinisospin independent twobody correlations , , and are different for pairs of particles where both are in the shell, both in the shell, or one in each. The comes from the structure of an alpha particle, the is constructed to go to unity at large distances. The is set to give the appropriate cluster structure outside the particle core, for example it is similar to a deuteron for Li and to a triton for Li; see Pudliner et al. (1997) for more details.
Except for closedshell nuclei, the complete trial wave function is constructed by taking a linear set of states of the form in Eq. (30) with the same total angular momentum and parity. Typically these correspond to the lowest shellmodel states of the system. QMC methods are then used to compute the Hamiltonian and normalization matrix elements in this basis. These coefficients are often similar in magnitude to those produced by a very small shellmodel calculation of the same nucleus. In light nuclei coupling is most efficient; examples of the diagonalization may be found in Pudliner et al. (1997); Wiringa et al. (2000); Pieper et al. (2002) and compared to traditional shell model studies such as Kumar (1974). The VMC calculations give very good descriptions of inclusive observables including momentum distributions, but the energies and other observables can then be improved, using the results of the VMC diagonalization to initiate the GFMC calculations.
Computational Implementation
The spatial integrals in Eq. (21) are evaluated using Metropolis Monte Carlo techniques Metropolis et al. (1953). A weight function is first defined to sample points in 3Adimensional coordinate space. The simplest choice is , where the brackets indicate a sum over all the spin isospin parts of the wave function. For spinisospin independent interactions the particle wave function is a function of the coordinates of the system only, and the weight function is the square of the wave function. The Metropolis method allows one to sample points in largedimensional spaces with probability proportional to any positive function W through a suitable combination of proposed (usually local) moves and an acceptance or rejection of the proposed move based upon the ratio of the function W at the original or proposed points. Iterating these steps produces a set of points in 3A dimensional space with probability proportional to W(R).
For spinisospin dependent interactions, the wave function is a sum of complex amplitudes for each spinisospin state of the system:
(31) 
and the spin states are:
(32) 
and similarly for the isospin states with and instead of and . The isospin states can be reduced by using charge conservation to states and, by assuming the nucleus has good isospin , further reduced to
(33) 
components. The weight function in this case is the sum of the squares of the individual amplitudes: .
Given a set of coordinates , to calculate the wave function one must first populate the various amplitudes in the trial state by calculating the Slater determinant, BCS state, or Jastrow wave function . Spinisospin independent operators acting on are simple multiplicative constants for each amplitude . Pair correlation operators then operate on the ; these are sparse matrix multiplications for each pair. The sparse matrices are easily computed onthefly using explicitly coded subroutines Pieper (2008b). The product over pair correlations is built up by successive operations for each pair. For example, the effect of the operator on the wave function of threeparticles can be written as follows (The notation means the amplitude for nucleon 1 being spin up and nucleons 2 and 3 being spin down; the isospin components have been omitted for simplicity):
(34) 
Metropolis Monte Carlo is used to sample points in the 3dimensional space by accepting and rejecting trial moves of the particles. Enforcing detailed balance ensures that the asymptotic distribution of such points will be distributed according to the weight . The energy can then be computed as the average over the points in the random walk:
(35) 
where the angled brackets imply the sum over spin and isospin states for each set of spatial coordinates . The matrix elements of the Hamiltonian are evaluated using the same techniques as those used for the pair correlation operators.
The computational time for the VMC method scales exponentially with the particle number. At first glance, this may seem to be because of the explicit sums over exponentially large number of spinisospin amplitudes calculated from the trial wave function. If that were the only reason, it would be trivial to sample the spinisospin state and evaluate the trial wave function’s amplitude for that sampled spinisospin state. This sampling can in fact be done but the fundamental problem remains that good trial wave functions constructed as described in Eqs. (22–24), require exponential in the particle number operations to evaluate either a single spinisospin amplitude or all of them. Evaluating a single amplitude provides negligible savings, so the computational time is reduced by explicitly summing over the amplitudes, which removes any variance that would occur from sampling. If trial wave functions could be constructed which capture the important physics, while requiring computational time that scales polynomially with particle number for a single spinisospin amplitude, VMC calculations would be straightforward for all nuclei.
In reality one does not usually compute the full wave function with all orders of pair operators implied by the symmetrization operator in the definition of the wave function. One can sample the orders of the pairs independently for the left and right (bra and ket) wave functions of Eq. (35), and define a slightly more complicated positive definite form for the weight function in terms of the two sets of amplitudes and for the order of pair operators and in the left and righthand wave functions. From several thousand to several tens of thousands of points are sufficient for a typical evaluation of the energy, and statistical errors are obtained using standard techniques.
To search for optimal variational parameters embedded in , it is very useful to first generate a Monte Carlo walk with configurations and weights for a given parameter set. Then one can change one or more parameters and reuse the same set of configurations to evaluate the change in the energy. The correlated energy difference will have a much smaller statistical error than differencing two large energies obtained from independent random walks. In this manner, a chain of small incremental improvements can be developed that leads to a lower variational energy. When the norm of the improved wave function starts to differ significantly from the original walk, a new reference walk can be made and the search continued from that set.
One way to overcome the exponential growth in computational requirements and access larger nuclei is to use a cluster expansion. Cluster expansions in terms of the operator correlations in the variational wave function were developed more than two decades ago and used in the first QMC calculations of O Pieper et al. (1992). In these calculations a full dimensional integral was done for the Jastrow part of the wave function while up to fournucleon linkedclusters were used for the operator terms. Earlier versions of the Argonne and Urbana interactions were used. Given the tremendous increase in computer power since then, this method might profitably be reconsidered for calculations of much bigger nuclei.
iii.2 Green’s function Monte Carlo
GFMC methods are used to project out the ground state with a particular set of quantum numbers. GFMC methods were invented in the 1960s Kalos (1962) and have been applied to many different problems in condensed matter, chemistry, and related fields. They are closely related to finitetemperature algorithms which calculate the density matrix Ceperley (1995), but they use trial wave functions on the boundaries of the paths to project out the quantum numbers of specific states.
GFMC typically starts from a trial wave function and projects:
(36) 
where is a parameter used to control the normalization. For stronglyinteracting systems one cannot compute directly, however one can compute the hightemperature or shorttime propagator, and insert complete sets of states between each shorttime propagator,
(37) 
and then use Monte Carlo techniques to sample the paths in the propagation. The method is accurate for small values of the time step , and the accuracy can be determined by simulations using several different values of the time step and extrapolating to zero. In the GFMC method, Monte Carlo is used to sample the coordinates ; Eq. (37) also has an implied sum over spin and isospin states at each step of the walk which is calculated explicitly.
ImaginaryTime Propagator
In the simplest approximation the propagator:
(38)  
where is the nonrelativistic kinetic energy:
(39) 
with , yielding a Gaussian diffusion for the paths. The matrix is the spin and isospindependent interaction:
(40) 
where indicates a symmetrization over orders of pairs. Each pair interaction can be simply evaluated as the exponent of a small spinisospin matrix. This treatment is adequate for static spindependent interactions.
In practice one needs to include momentumdependent spinorbit interactions as well as interactions. It is more efficient to calculate the propagator explicitly, storing the radial and spinisospin dependence on a grid for each initial and final state. This is done by calculating the propagator independently in each partial wave and then summing them to create the full propagator. This was first done in studies of liquid Helium Ceperley (1995); Schmidt and Lee (1995) and then adapted to the nuclear physics case Pudliner et al. (1997). This has the advantage of summing all interactions for each pair explicitly, allowing for larger time steps in the pathintegral simulation. The propagator is defined:
(41) 
where and are the initial and final relative coordinates, is the Hamiltonian including relative kinetic energy and the interaction, and and are initial and final spin and isospin states, respectively. The pair propagator is calculated for the AV8 Hamiltonian, denoted as . At present higher order terms in the momenta are treated perturbatively. Though the pair propagator can be calculated for these interactions, the Monte Carlo sampling can lead to large variance Lynn and Schmidt (2012).
The pair propagators are then combined to produce the full propagation matrix for the system. The interaction is included symmetrically, and the full propagation matrix for each step can then be written as:
(42) 
The spinorbit interaction in the product of propagators with the full v interaction yields spurious interactions resulting from quadratic terms in the difference from different pairs. One can correct for this but in practice the effect is not significant. Using the calculated propagators allows for a factor of 510 larger time steps than the simple approximation in Eq.(III.2.1) Pudliner et al. (1997).
Implementation
Once the propagator for each step is specified, an algorithm must be chosen to sample over all possible paths. A branching random walk algorithm very similar to that used in standard diffusion Monte Carlo (DMC) Foulkes et al. (2001) is used. This random walk does not sample the entire path at once; it uses Markov Chain Monte Carlo to perform each step given the present coordinates and amplitudes in the propagated wave function. One difference with standard DMC is that the importance sampled Green’s function is explicitly sampled rather than using a small timestep extrapolation for the wave functions.
A positive definite “weight” is first defined as a function of the trial function and the propagated wave function . Typically the form used is
(43) 
where is a small parameter. Sampling of the paths and branching for the importance function are then implemented with the scalar function . Given the present position , several different possible final states are sampled from the free propagator . For each sample of the corresponding configuration is included in the sample. The weight function is then calculated for each of the possible new points , and the final point is chosen according to the relative weights and scaled with the ratio of the average to the actual . Branching is performed with the ratio of weight functions after and before the step, or typically after several steps. The weights of different paths used to calculate observables will eventually diverge, yielding the entire contribution from only a few paths that dominate. This is commonly avoided by using the branching technique, in which the configurations are redistributed by killing or making copies of each one according to
(44) 
where is the weight of the ith configuration obtained by multiplying the weight of Eq. (43) by ( is the spin/isospin independent part of the potential), is a random number with uniform distribution between 0 and 1, and in the above equation means the truncated integer number of the argument. Different random number seeds are given to new copies generated from the same walker. This procedure guarantees that the configurations with small weight, contributing by generating only noise to the observables, are dropped. The full procedure is described in Pudliner et al. (1997).
After every typically 20 to 40 steps, the energy as a function of imaginary time is calculated as:
(45) 
where the sum over indicates the sum over samples of the wave function. The brackets in the numerator and denominator of the last expression indicate sums over spins and isospins for each sample. The initially decrease rapidly from the VMC () energy but then stabilizes and just fluctuates within the statistical errors; examples of this are shown in Fig. 2, discussed below, and also in Sec. V.4. These stable values are averaged to get the converged GFMC results.
In principle, the GFMC propagation should converge to the lowestenergy state of given quantum numbers . The nuclei considered here may have a few particlestable and multiple particleunstable excited states of the same quantum numbers. In practice, GFMC propagation can obtain good energy estimates for many of these additional states.
First, a set of orthogonal VMC trial functions are generated that are diagonalized in the small singleparticle shell basis of differing and spatial symmetry combinations that can make a given . These are pseudobound wave functions that fall off exponentially at long range, with matter radii not much larger than the ground state. Then independent GFMC propagations are carried out starting from each of these trial functions. An example is shown in Fig. 2 for the four shell states in Li, all of which are particleunstable Pieper et al. (2004). The GFMC propagations stay nearly orthogonal to fairly large MeV, as shown by the solid symbols in the figure. The overlaps between different states can be evaluated, and an explicit reorthogonalization made, shown by the open symbols. The states remain wellseparated in energy.
The first state in Fig. 2 is physically wide ( keV) because it has the spatial symmetry of alpha plus triton and is several MeV above the threshold for breakup into separated clusters. Consequently, a GFMC propagation is expected to eventually drop to that threshold energy, and the figure shows, after a rapid initial drop from 26 to 32 MeV by MeV, a slowly decreasing energy as increases, reaching 33 MeV at MeV. In cases like this, the energy is quoted at the small value of where the rapid initial improvement over the variational starting point has saturated. The second state in Fig. 2 is physically narrow ( keV) because it has a spatial symmetry like Li+ and is only 20 keV above that breakup threshold. The GFMC propagation shows the same rapid initial drop in energy, and then no appreciable further decline, allowing us to identify a clear energy for this state. The third and fourth states are not experimentally identified, but from the GFMC propagation behavior we would expect the third state to be physically narrow, and the fourth to be fairly broad. An alternative approach to calculate systems in the continuum by imposing specific boundary conditions is presented in sec. IV.6.
In general the GFMC method suffers from the fermion sign problem, in that the numerator and denominator of Eq. (45) tend to have an increasing ratio of error to signal for a finite sample size and large imaginary times . Other than for a few special cases such as purely attractive interactions, Hubbard models at halffilling, or lattice QCD at zero chemical potential, QMC methods typically all have this difficulty. This is basically because when is not exact it contains contamination from the Bosonic groundstate that will be unavoidably sampled. For scalar potentials, or in any case where a real wave function can be used, the sign problem is avoided by using the fixednode approximation, and the problem is solved in a restricted (Bosonic) subspace, where the trial wave function always maintains the same sign. In this case the problem would be exactly solved if the nodes of the true groundstate were known. Because this is not the case, the solution obtained is a rigorous upperbound to the true groundstate energy Moskowitz et al. (1982). For spinisospin dependent Hamiltonians a complex wave function must be used, and the general fixednode approximation does not apply. Instead the sign problem is circumvented by using a ‘constrained path’ algorithm, essentially limiting the original propagation to regions where the propagated and trial wave functions have a positive overlap. This approximation, like the fixednode algorithm for spinindependent interaction, involves discarding configurations that have zero overlap with the trial wave function. As such, they are exact for the case when the trial wave function is exact and are therefore variational. However, unlike the fixednode case, the constrained path method does not provide upper bounds Wiringa et al. (2000).
To address the possible bias introduced by the constraint, all the configurations (including those that would be discarded) for a previous number of steps are used when evaluating energies and other expectation values. is chosen to be as large a number of time steps as feasible with reasonable statistical error (again typically 20 to 40 steps). Tests using different trial functions and very long runs indicate that energies in shell nuclei are accurate to around one per cent using these methods. This has been tested in detail in Wiringa et al. (2000), where the use of different wave functions is discussed.
Expectation values other than the energy are typically calculated from “mixed” estimates; for diagonal matrix elements this is:
(46) 
The above equation can be verified by assuming that the true ground state is well represented by the variational wave function and a small perturbation, i.e., , and is a small parameter. Since the variational wave functions are typically very good the extrapolation is quite small. This can be further tested by using different trial wave functions to extract the same observable, or using the HellmanFeynman theorem. For the case of simple static operators, improved methods are available that propagate both before and after the insertion of the operator Liu et al. (1974), i.e. directly calculating operators with on both sides. However these techniques might be very difficult to apply for nonlocal operators.
Because a Hamiltonian commutes with itself, the total energy of the Hamiltonian used to construct the propagator [Eq. (42)] is not extrapolated; thus this total energy is not the sum of its extrapolated pieces, rather the sum differs by the amount the energy was improved. As was noted above, the full AV18 potential cannot be used in the propagator; rather an containing the AV8 approximation to AV18 is used. In practice AV8 gives slightly more binding than AV18 so the the repulsive part of the potential is increased to make . The difference must be extrapolated by Eq. (46). The best check of the systematic error introduced by this procedure is given by comparing GFMC calculations of H and He energies with results of more accurate fewnucleon methods; this suggests that the error is less than 0.5% Pudliner et al. (1997).
In the case of offdiagonal matrix elements, e.g., in transition matrix elements between initial and final wave functions, Eq. (46) generalizes to:
(47)  
Technical details can be found in Pervin et al. (2007).
Recently, the capability to make correlated GFMC propagations has been added Lovato et al. (2015). In these calculations, the values of for every time step, the corresponding weights , and other quantities are saved during an initial propagation. Subsequent propagations for different initial or different nuclei (such as isobaric analogs) then follow the original propagation and correlated differences of expectation values can be computed with much smaller statistical errors than for the individual values.
iii.3 Auxiliary Field Diffusion Monte Carlo
The GFMC method works very well for calculating the low lying states of nuclei up to C. Its major limitation is that the computational costs scale exponentially with the number of particles, because of the full summations of the spinisospin states. An alternative approach is to use a basis given by the outer product of nucleon position states, and the outer product of single nucleon spinisospin spinor states. An element of this overcomplete basis is given by specifying the Cartesian coordinates for the nucleons, and specifying four complex amplitudes for each nucleon to be in a spinisospin state. A basis state is then defined
(48) 
The trial functions must be antisymmetric under interchange. The only such functions with polynomial scaling are Slater determinants or Pfaffians (BCS pairing functions), for example,
(49) 
or linear combinations of them. Operating on these with the product of correlation operators, Eq. (23), again gives a state with exponential scaling with nucleon number. In most of the AFDMC calculations, these wave functions include a stateindependent, or central, Jastrow correlation:
(50) 
Calculations of the Slater determinants or Pfaffians scale like when using standard dense matrix methods, while the central Jastrow requires operations if its range is the same order as the system size. These trial functions capture only the physics of the gross shell structure of the nuclear problem and the stateindependent part of the twobody interaction. Devising trial functions that are both computationally efficient to calculate and that capture the statedependent two and threebody correlations that are important would greatly improve both the statistical and systematic errors of QMC methods for nuclear problems.
The trial wave functions above can be used for variational calculations. However, the results are poor since the functions miss the physics of the important tensor interactions. More recently the improved form