text:
# Quantum Monte Carlo Calculations of Light Nuclei

## 1 INTRODUCTION

A major goal of nuclear physics is to understand the stability, structure, and reactions of nuclei as a consequence of the interactions between individual nucleons. This is an extremely challenging many-body problem, exacerbated by the fact that we do not know in detail what the interactions are. Quantum chromodynamics, the fundamental theory for strong interactions, is so difficult to solve in the nonperturbative regime of low-energy nuclear physics that no one is currently able to quantitatively describe the interactions between two nucleons from first principles. However, we do have a tremendous amount of experimental data on nucleon-nucleon ($`NN`$) scattering and we can construct accurate representations of $`NN`$ interactions by means of two-body potentials. These potentials are complicated, depending on the relative positions, spins, isospins, and momenta of the nucleons, which makes the problem of finding accurate quantum-mechanical solutions of many-nucleon bound and scattering states a demanding task. Nevertheless, a variety of precise calculations have been carried out for the three- and four-nucleon systems which reveal an additional complication: two-nucleon forces alone are inadequate to quantitatively explain either the bound-state properties or a variety of $`Nd`$ and $`Nt`$ scattering data. This is not surprising, because meson-exchange theory and the composite nature of nucleons indicate that many-nucleon forces, at least three-nucleon forces, are important. However, these are even more difficult to construct from first principles. In like manner, meson-exchange introduces two-body charge and current operators, and calculations of electroweak transitions show that these are also necessary to achieve agreement with data.

Despite these difficulties, tremendous progress has been made in the past decade, both in the characterization of nuclear forces and currents, and in the development of accurate many-body techniques for evaluating them. A multienergy partial-wave analyses of elastic $`NN`$ scattering data below $`T_{lab}=350`$ MeV was produced by the Nijmegen group in 1993 (?), which demonstrated that over 4300 data points could be fit with a $`\chi ^2`$ per datum $``$ 1. This was achieved by careful attention to the known long-range part of the $`NN`$ interaction and a rigorous winnowing of inconsistent data. Their work inspired the construction of a number of new $`NN`$ potential models which gave comparably good fits to the database.

Concurrently, several essentially exact methods have been developed for the study of few-nucleon ($`A`$ = 3,4) systems with such realistic $`NN`$ potentials. The first accurate, converged calculations of three-nucleon ($`NNN`$) bound states were obtained by configuration-space Faddeev methods in 1985 (?), while the first accurate 4He bound states were obtained by the Green’s function Monte Carlo method in 1987 (?). Other methods of comparable accuracy have been developed since then for both $`A`$ = 3,4 bound states and for $`NNN`$ scattering. A comprehensive review of many aspects of the few-nucleon problem can be found in the article by Carlson and Schiavilla (?). The result of all the progress in the last decade is that we have obtained an excellent understanding of few-nucleon systems, including both structure and reactions, although a few unanswered problems remain.

Until recently, nuclei larger than $`A=4`$ have been described within the framework of the nuclear shell model or mean-field pictures such as Skyrme models. The classic description of p-shell ($`5A16`$) nuclei was developed by Cohen and Kurath in the 1960s (?). In that work, a large number of experimental levels were fit in terms of a relatively small number of two-body matrix elements, but no direct connection was made to the underlying $`NN`$ force. Beginning with the work of Kuo and Brown (?), approximate methods for relating the shell-model matrix elements to the underlying forces were developed and applied to even larger nuclei. However, in most calculations, only the valence nucleons, i.e., the nucleons in the last unfilled shell, are treated as active, while the core is inert. Recently the first ”no-core” shell model calculations, in which all the nucleons are active, have been made in the p-shell with effective interactions that have been derived by a G-matrix procedure from a realistic underlying $`NN`$ force (?). These calculations agree with the other exact methods for the $`A=3,4`$ nuclei, but fully-converged calculations for larger systems have not yet been obtained.

The present article focuses on recent developments in quantum Monte Carlo methods that are making the light p-shell nuclei accessible at a level of accuracy close to that obtained for the s-shell ($`A4`$) nuclei. The quantum Monte Carlo methods include variational Monte Carlo (VMC) and Green’s function Monte Carlo (GFMC) methods. The VMC is an approximate method (?, ?, ?, ?) that is used as a starting point for the more accurate GFMC calculations. Because of its relative simplicity, the VMC has also been used to make first studies of various reactions, such as the nuclear response to electron (?, ?) or pion scattering (?), and for electroweak capture reactions of astrophysical interest (?, ?). The GFMC is exact in the sense that binding energies can be calculated with an accuracy of better than 2%; it has been used to compute ground states and low-lying excitations of all the $`A8`$ nuclei (?, ?, ?, ?). These studies of light p-shell nuclei have shown that nuclear structure can be predicted from the elementary nuclear forces. As presently implemented, both the VMC and GFMC methods require significant additional computer resources for every nucleon added, so they probably will be limited to systems of $`A12`$ nuclei. A new method, whose accuracy lies between that of the VMC and GFMC methods, is auxiliary-field diffusion Monte Carlo (AFDMC); it has been implemented for large pure neutron systems (?) and with further development should be capable of handling nuclei above $`A=12`$.

Most of the calculations discussed here have been made using Hamiltonians containing the Argonne $`v_{18}`$ $`NN`$ potential (?) alone or with one of the Urbana (?) or Illinois (?) series of $`NNN`$ potentials. Argonne $`v_{18}`$ (AV18) is representative of the modern $`NN`$ potentials that give accurate fits to scattering data, while the Urbana and Illinois models are modern $`NNN`$ potentials based on meson-exchange interactions and fit to the binding energies of light nuclei. These models are described in some detail in Section 2 to show the complexity of modern nuclear forces, and the consequent challenge for many-body theory. The VMC method is presented in Section 3; it starts with the construction of a variational wave function of specified angular momentum, parity and isospin, $`\mathrm{\Psi }_V(J^\pi ;T)`$, using products of two- and three-body correlation operators acting on a fully antisymmetrized set of one-body basis states coupled to specific quantum numbers. Metropolis Monte Carlo integration (?) is used to evaluate $`\mathrm{\Psi }_V|H|\mathrm{\Psi }_V`$, giving an upper bound to the energy of the state. The GFMC method, described in Section 4, is a stochastic method that systematically improves on $`\mathrm{\Psi }_V`$ by projecting out excited state contamination using the Euclidean propagation $`\mathrm{\Psi }(\tau )=\mathrm{exp}[(HE_0)\tau ]\mathrm{\Psi }_V`$. In the limit $`\tau \mathrm{}`$, this leads to the exact $`|H|`$.

Energy results for the ground states and many low-lying excited states of the $`A8`$ nuclei are presented in Section 5. With the new Hamiltonians, the excitation spectra are reproduced very well, including the relative stability of different nuclei and the splittings between excited states. Energy differences between isobaric analog states and isospin-mixing matrix elements have also been computed to study the charge-independence-breaking of the strong and electromagnetic forces. A variety of density and momentum distributions, electromagnetic moments and form factors, transition densities, spectroscopic factors, and intrinsic deformation are presented in Section 6. One of the chief benefits of this approach is that absolute normalizations can be computed, and effective charges so commonly needed in traditional shell-model calculations of electroweak transitions are not required.

While the quantum Monte Carlo methods have been used primarily for stationary states, it is also possible to study scattering states, as discussed in Section 7. Some first calculations of astrophysical radiative capture reactions in the p-shell are presented in Section 8. Both of these applications are in their early stages, but promise to be major fields of endeavor for future work. In Section 9 we discuss neutron drops, systems of interacting neutrons confined in an artificial external well, that can serve as reference points for the development of Skyrme interactions for use in very large neutron-rich nuclei. Finally, in Section 10, we consider the outlook for future developments.

## 2 HAMILTONIAN AND CURRENTS

The Hamiltonian used in this review use includes nonrelativistic kinetic energy, the AV18 $`NN`$ potential (?) and either the Urbana IX (?) or one of the new Illinois series of $`NN`$ potential (?):

$$H=\underset{i}{}K_i+\underset{i<j}{}v_{ij}+\underset{i<j<k}{}V_{ijk}.$$
(1)
The kinetic energy operator has charge-independent (CI) and charge-symmetry-breaking (CSB) components, the latter due to the difference between proton and neutron masses,

$$K_i=K_i^{CI}+K_i^{CSB}\frac{\mathrm{}^2}{4}(\frac{1}{m_p}+\frac{1}{m_n})_i^2\frac{\mathrm{}^2}{4}(\frac{1}{m_p}\frac{1}{m_n})\tau _{zi}_i^2,$$
(2)
where $`𝝉_i`$ is the isospin of nucleon $`i`$.

The AV18 potential is written as a sum of electromagnetic and one-pion-exchange terms and a shorter-range phenomenological part,

$$v_{ij}=v_{ij}^\gamma +v_{ij}^\pi +v_{ij}^R.$$
(3)
The electromagnetic terms include one- and two-photon-exchange Coulomb interaction, vacuum polarization, Darwin-Foldy, and magnetic moment terms, all with appropriate proton and neutron form factors. The one-pion-exchange part of the potential includes the charge-dependent (CD) terms due to the difference in neutral and charged pion masses. It can be written as

$$v_{ij}^\pi =f^2\left(\frac{m}{m_s}\right)^2\frac{1}{3}mc^2\left[X_{ij}𝝉_i𝝉_j+\stackrel{~}{X}_{ij}T_{ij}\right],$$
(4)
where $`T_{ij}=3\tau _{zi}\tau _{zj}𝝉_i𝝉_j`$ is the isotensor operator and

$`X_{ij}=\frac{1}{3}\left(X_{ij}^0+2X_{ij}^\pm \right),`$ (5)
$`\stackrel{~}{X}_{ij}=\frac{1}{3}\left(X_{ij}^0X_{ij}^\pm \right),`$ (6)
$`X_{ij}^m=\left[Y(mr_{ij})𝝈_i𝝈_j+T(mr_{ij})S_{ij}\right].`$ (7)
Here $`Y(x)`$ and $`T(x)`$ are the normal Yukawa and tensor functions:

$`Y(x)`$ $`=`$ $`{\displaystyle \frac{e^x}{x}}\xi _Y(r),`$ (8)
$`T(x)`$ $`=`$ $`\left({\displaystyle \frac{3}{x^2}}+{\displaystyle \frac{3}{x}}+1\right)Y(x)\xi _T(r).`$ (9)
with short-range cutoff functions $`\xi _Y(r)`$ and $`\xi _T(r)`$, and the $`X^{\pm ,0}`$ are calculated with $`m=m_{\pi ^\pm }`$ and $`m_{\pi ^0}`$. The coupling constant in Eq.(4) is $`f^2=0.075`$ and the scaling mass $`m_s=m_{\pi ^\pm }`$.

The remaining terms are of intermediate (two-pion-exchange) and short range, with some 40 adjustable parameters. The one-pion-exchange and the remaining phenomenological part of the potential can be written as a sum of eighteen operators,

$$v_{ij}^\pi +v_{ij}^R=\underset{p=1,18}{}v_p(r_{ij})O_{ij}^p.$$
(10)
The first fourteen are CI operators,

$`O_{ij}^{p=1,14}=[1,𝝈_i𝝈_j,S_{ij},𝐋𝐒,𝐋^2,𝐋^2(𝝈_i𝝈_j),(𝐋𝐒)^2][1,𝝉_i𝝉_j],`$ (11)
while the last four,

$$O_{ij}^{p=15,18}=[1,𝝈_i𝝈_j,S_{ij}]T_{ij},\mathrm{and}(\tau _{zi}+\tau _{zj}),$$
(12)
are CD and CSB terms. The potential was fit directly to the Nijmegen $`NN`$ scattering data base (?) containing 1787 $`pp`$ and 2514 $`np`$ data in the range 0–350 MeV, to the $`nn`$ scattering length, and to the deuteron binding energy, with a $`\chi ^2`$ per datum of 1.09.

While the CD and CSB terms are small, there is a clear need for their presence. The Nijmegen group studied many $`NN`$ potentials from the 1980s and before, and found that potentials fit to $`np`$ data in $`T=1`$ states did not fit $`pp`$ data well even after allowing for standard Coulomb effects, and vice versa (?). All five modern $`NN`$ potentials, Reid93, Nijm I and II (?), CD Bonn (?), and AV18, which fit $`pp`$ and $`np`$ scattering data with a $`\chi ^2/`$point near 1, have CD components. The CSB term is required to fit the difference between the $`pp`$ and $`nn`$ scattering lengths, and is consistent with the mass difference between 3H and 3He, as shown below. The identification of the proper CI, CD, and CSB $`NN`$ force components is important in setting the proper baseline for the introduction of $`NNN`$ forces, which are required to fit the 3H and 3He binding energies as a primary constraint.

The Urbana series of three-nucleon potentials is written as a sum of two-pion-exchange P-wave and remaining shorter-range phenomenological terms,

$$V_{ijk}=V_{ijk}^{2\pi ,P}+V_{ijk}^R.$$
(13)
The structure of the two-pion P-wave exchange term with an intermediate $`\mathrm{\Delta }`$ excitation (Fig. 1a) was originally written down by Fujita and Miyazawa (?); it can be expressed simply as

$$V_{ijk}^{2\pi ,P}=\underset{cyc}{}\left(A_{2\pi }^P\{X_{ij}^\pi ,X_{jk}^\pi \}\{𝝉_i𝝉_j,𝝉_j𝝉_k\}+C_{2\pi }^P[X_{ij}^\pi ,X_{jk}^\pi ][𝝉_i𝝉_j,𝝉_j𝝉_k]\right),$$
(14)
where $`X_{ij}^\pi `$ is constructed with an average pion mass, $`m_\pi =\frac{1}{3}m_{\pi ^0}+\frac{2}{3}m_{\pi ^\pm }`$, and $`_{cyc}`$ is a sum over the three cyclic exchanges of nucleons $`i,j,k`$. For the Urbana models $`C_{2\pi }^P=\frac{1}{4}A_{2\pi }^P`$, as in the original Fujita-Miyazawa model (?), while other potentials like the Tucson-Melbourne (?) and Brazil (?) models, have a ratio slightly larger than $`\frac{1}{4}`$. The shorter-range phenomenological term is given by

$$V_{ijk}^R=\underset{cyc}{}A_RT^2(m_\pi r_{ij})T^2(m_\pi r_{jk}).$$
(15)
For the Urbana IX (UIX) model, the parameters were determined by fitting the binding energy of 3H and the density of nuclear matter in conjunction with AV18.

As shown below, while the combined AV18/UIX Hamiltonian reproduces the binding energies of s-shell nuclei, it does not do so well for light p-shell nuclei. Recently a new class of $`NNN`$ potentials, called the Illinois models, has been developed to address this problem (?). These potentials contain the Urbana terms and two additional terms, resulting in a total of four coupling constants that can be adjusted to fit the data.

One term, $`V_{ijk}^{2\pi ,S}`$, is due to $`\pi N`$ S-wave scattering as illustrated in Fig. 1b. It has been included in a number of $`NNN`$ potentials like the Tucson-Melbourne (?) and Brazil (?) models. The Illinois models use the form recommended in the latest Texas model (?), where chiral symmetry is used to constrain the structure of the interaction. However, in practice, this term is much smaller than the $`V_{ijk}^{2\pi ,P}`$ 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 three-pion rings containing one or two deltas (Fig. 1c,d). As discussed in Ref. (?), these diagrams result in a large number of terms, the most important of which are used to construct the Illinois models:

$$V_{ijk}^{3\pi ,\mathrm{\Delta }R}=A_{3\pi }^{\mathrm{\Delta }R}\left[\frac{50}{3}S_\tau ^IS_\sigma ^I+\frac{26}{3}A_\tau ^IA_\sigma ^I\right].$$
(16)
Here the $`S_x^I`$ and $`A_x^I`$ are operators that are symmetric or antisymmetric under any exchange of the three nucleons, and the subscript $`x=\sigma `$ or $`\tau `$ indicates that the operators act on, respectively, spin and space or just isospin degrees of freedom.

The $`S_\tau ^I`$ is a projector onto isospin-$`\frac{3}{2}`$ triples:

$$S_\tau ^I=2+\frac{2}{3}\left(𝝉_i𝝉_j+𝝉_j𝝉_k+𝝉_k𝝉_i\right)=4P_{T=3/2}.$$
(17)
To the extent isospin is conserved, there are no such triples in the s-shell nuclei, and so this term does not affect them. It is also zero for $`Nd`$ scattering. However, the $`S_\tau ^IS_\sigma ^I`$ term is attractive in all the p-shell nuclei studied. The $`A_\tau ^I`$ has the same structure as the isospin part of anticommutator part of $`V^{2\pi ,P}`$, but the $`A_\tau ^IA_\sigma ^I`$ term is repulsive in all nuclei studied so far. In p-shell nuclei, the magnitude of the $`A_\tau ^IA_\sigma ^I`$ term is smaller than that of the $`S_\tau ^IS_\sigma ^I`$ term, so the net effect of the $`V_{ijk}^{3\pi ,\mathrm{\Delta }R}`$ is slight repulsion in s-shell nuclei and larger attraction in p-shell nuclei. The reader is referred to Ref. (?), and its appendix, for the complete structure of $`V_{ijk}^{3\pi ,\mathrm{\Delta }R}`$.

Relativistic corrections to the Hamiltonian of Eq.(1) have been studied in three- and four-body nuclei (?, ?, ?). We only briefly outline the results here. If a square-root kinetic energy,

$$K_i^{rel}=\sqrt{p_i^2+m_i^2}m_i,$$
(18)
where $`p_i`$ is the momentum of nucleon $`i`$, is introduced into the Hamiltonian, then the $`NN`$ potential must be readjusted to refit the $`NN`$ data. This can be done for AV18 with only small adjustments of the parameters, to produce a $`v^{rel}`$. VMC calculations show that the resulting $`K^{rel}+v^{rel}K+v`$ for 3H and 4He; thus the square-root kinetic energy may be neglected in light nuclei.

A second relativistic correction, the so-called “boost correction,” arises from the fact that the two-nucleon interaction $`v_{ij}`$ depends both on the relative momentum $`𝐩=(𝐩_i𝐩_j)/2`$ and the total momentum $`𝐏=𝐩_i+𝐩_j`$ of the interacting nucleons,

$$v_{ij}(𝐩,𝐏)=\stackrel{~}{v}_{ij}(𝐩)+\delta v_{ij}(𝐩,𝐏),$$
(19)
where $`\delta v_{ij}(𝐩,𝐏=0)=0`$. Fitting the two-body potential to two-nucleon scattering data determines only $`\stackrel{~}{v}_{ij}`$. The $`\delta v_{ij}`$ can be determined from $`\stackrel{~}{v}_{ij}`$; a suitable approximation (?, ?, ?) is

$$\delta v_{ij}(𝐏)=\frac{P^2}{8m^2}\stackrel{~}{v}_{ij}+\frac{1}{8m^2}[𝐏𝐫𝐏\mathbf{},\stackrel{~}{v}_{ij}],$$
(20)
where only the first six terms (the static terms) of $`\stackrel{~}{v}_{ij}`$ are retained and the gradient operators do not act on $`\stackrel{~}{v}_{ij}`$. The $`\delta v_{ij}`$ is meant to be used only in first-order perturbation in a non-relativistic wave function. Its expectation values in light nuclei are nearly proportional to $`V_{ijk}^R`$ (?), and thus one can consider that part of $`A_R`$, Eq.(15), comes from $`\delta v_{ij}`$. This proportionality does not, however, extend to nuclear matter or pure neutron systems. There one must use a reduced $`A_R^{}`$ that does not contain the part ascribable to $`\delta v_{ij}`$, and explicitly evaluate $`\delta v_{ij}`$; doing so results in a softer equation of state for dense matter.

## 3 VARIATIONAL MONTE CARLO

The variational method can be used to obtain approximate solutions to the many-body Schrödinger equation, $`H\mathrm{\Psi }=E\mathrm{\Psi }`$, for a wide range of nuclear systems, including few-body nuclei, light closed-shell nuclei, nuclear matter, and neutron stars (?). A suitably parameterized wave function, $`\mathrm{\Psi }_V`$, is used to calculate an upper bound to the exact ground-state energy,

$$E_V=\frac{\mathrm{\Psi }_V|H|\mathrm{\Psi }_V}{\mathrm{\Psi }_V|\mathrm{\Psi }_V}E_0.$$
(21)
The parameters in $`\mathrm{\Psi }_V`$ are varied to minimize $`E_V`$, and the lowest value is taken as the approximate solution. Upper bounds to excited states are also obtainable, if they have different quantum numbers from the ground state, or from small-basis diagonalizations if they have the same quantum numbers. The corresponding $`\mathrm{\Psi }_V`$ can be used to calculate other properties, such as particle density or electromagnetic moments, or to start a Green’s function Monte Carlo calculation. In this section we describe the ansatz for $`\mathrm{\Psi }_V`$ for light nuclei and briefly review how the expectation value is evaluated and the parameters of $`\mathrm{\Psi }_V`$ are fixed.

The best variational wave function has the form (?)

$$|\mathrm{\Psi }_V=\left[1+\underset{i<j<k}{}(U_{ijk}+U_{ijk}^{TNI})+\underset{i<j}{}U_{ij}^{LS}\right]|\mathrm{\Psi }_P,$$
(22)
where the pair wave function, $`\mathrm{\Psi }_P`$, is given by

$$|\mathrm{\Psi }_P=\left[𝒮\underset{i<j}{}(1+U_{ij})\right]|\mathrm{\Psi }_J.$$
(23)
The $`U_{ij}`$, $`U_{ij}^{LS}`$, $`U_{ijk}`$, and $`U_{ijk}^{TNI}`$ are noncommuting two- and three-nucleon correlation operators, and $`𝒮`$ is a symmetrization operator. The form of the antisymmetric Jastrow wave function, $`\mathrm{\Psi }_J`$, depends on the nuclear state under investigation. For the s-shell nuclei the simple form

$$|\mathrm{\Psi }_J=\left[\underset{i<j<k}{}f_{ijk}^c(𝐫_{ik},𝐫_{jk})\underset{i<j}{}f_c(r_{ij})\right]|\mathrm{\Phi }_A(JMTT_3)$$
(24)
is used. Here $`f_c(r_{ij})`$ and $`f_{ijk}^c`$ are central (spin-isospin independent) two- and three-body correlation functions and $`\mathrm{\Phi }_A`$ is an antisymmetrized spin-isospin state, e.g.,

$`|\mathrm{\Phi }_3(\frac{1}{2}\frac{1}{2}\frac{1}{2}\frac{1}{2})=𝒜|ppn,`$ (25)
$`|\mathrm{\Phi }_4(0000)=𝒜|ppnn,`$ (26)
with $`𝒜`$ the antisymmetrization operator.

The two-body correlation operators $`U_{ij}`$ and $`U_{ij}^{LS}`$ are sums of spin, isospin, tensor, and spin-orbit terms:

$`U_{ij}`$ $`=`$ $`{\displaystyle \underset{p=2,6}{}}\left[{\displaystyle \underset{ki,j}{}}f_{ijk}^p(𝐫_{ik},𝐫_{jk})\right]u_p(r_{ij})O_{ij}^p,`$ (27)
$`U_{ij}^{LS}`$ $`=`$ $`{\displaystyle \underset{p=7,8}{}}\left[{\displaystyle \underset{ki,j}{}}f_{ijk}^p(𝐫_{ik},𝐫_{jk})\right]u_p(r_{ij})O_{ij}^p,`$ (28)
where the $`O_{ij}^p`$ were introduced in Eq.(11). The spin-orbit correlations are only summed in Eq.(22) because of the extra computational expense of evaluating powers of $`𝐋𝐒`$ that would occur if it was inserted in the symmetrized product of Eq.(23).

The central $`f_c(r)`$ and noncentral $`u_p(r)`$ pair correlation functions reflect the influence of the two-body potential at short distances, while satisfying asymptotic boundary conditions of cluster separability. Reasonable functions are generated by minimizing the two-body cluster energy of a somewhat modified $`NN`$ interaction; this results in a set of eight coupled, Schrödinger-like, differential equations corresponding to linear combinations of the first eight operators i of $`v_{ij}`$, with a number of embedded variational parameters (?). The $`f_c(r)`$ is small at short distances, to reduce the contribution of the repulsive core of $`v_{ij}`$, and peaks at an intermediate distance corresponding to the maximum attraction of $`v_{ij}`$, as illustrated in Fig. 2 for several nuclei. For the s-shell nuclei, $`f_c(r)`$ falls off at larger distances to keep the system confined. For example, in 4He, $`[f_c(r)]^3\mathrm{exp}(\kappa r)/r`$ at large r, as shown by the dashed line in Fig. 2, where $`\kappa `$ corresponds to an $`\alpha t+p`$ separation energy. The noncentral $`u_p(r)`$ are all relatively small; the most important is the long-range tensor-isospin part $`u_{t\tau }(r)`$, also shown in Fig. 2 for several nuclei, which is mainly induced by the one-pion-exchange part of $`v_{ij}`$.

The $`f_{ijk}^c`$, $`f_{ijk}^p`$, and $`U_{ijk}`$ are three-nucleon correlations also induced by $`v_{ij}`$ (?). The $`U_{ijk}^{TNI}`$ are three-body correlations induced by the three-nucleon interaction, which have the form suggested by perturbation theory:

$$U_{ijk}^{TNI}=\underset{x}{}ϵ_xV_{ijk}(\stackrel{~}{r}_{ij},\stackrel{~}{r}_{jk},\stackrel{~}{r}_{ki}),$$
(29)
with $`\stackrel{~}{r}=yr`$, $`y`$ a scaling parameter, and $`ϵ_x`$ a (small negative) strength parameter. A somewhat simpler wave function to evaluate than $`\mathrm{\Psi }_V`$ is given by:

$`|\mathrm{\Psi }_T=\left[1+{\displaystyle \underset{i<j<k}{}}\stackrel{~}{U}_{ijk}^{TNI}\right]|\mathrm{\Psi }_P,`$ (30)
where $`\stackrel{~}{U}_{ijk}^{TNI}`$ is a slightly truncated TNI correlation. The $`\mathrm{\Psi }_T`$ gives about $`12`$ MeV less binding than the full $`\mathrm{\Psi }_V`$, but costs less than half as much to construct, making it a more efficient starting point for the GFMC calculation (?).

The Jastrow wave function for the light p-shell nuclei is more complicated, as a number of nucleons must be placed in the unfilled p-shell. The $`LS`$ coupling scheme is used to obtain the desired $`JM`$ value of a given state, as suggested by the shell-model studies of light p-shell nuclei (?). Different possible $`LS`$ combinations lead to multiple components in the Jastrow wave function. The possibility that the central correlations $`f_c(r_{ij})`$ could depend upon the shells (s or p) occupied by the particles and on the $`LS`$ coupling is also allowed for:

$`|\mathrm{\Psi }_J`$ $`=`$ $`𝒜\{{\displaystyle \underset{i<j<k}{}}f_{ijk}^c{\displaystyle \underset{i<j4}{}}f_{ss}(r_{ij}){\displaystyle \underset{k4<lA}{}}f_{sp}(r_{kl})`$
$`{\displaystyle \underset{LS[n]}{}}(\beta _{LS[n]}{\displaystyle \underset{4<l<mA}{}}f_{pp}^{LS[n]}(r_{lm})|\mathrm{\Phi }_A(LS[n]JMTT_3)_{1234:56\mathrm{}A})\}.`$
The operator $`𝒜`$ indicates an antisymmetric sum over all possible partitions of the $`A`$ particles into 4 s-shell and $`(A4)`$ p-shell ones. The central correlation $`f_{ss}(r)`$ is the $`f_c(r)`$ from the 4He wave function. The $`f_{sp}(r)`$, shown in Fig. 2 for Li nuclei, is similar to the $`f_c(r)`$ at short range, but with a long-range tail going to unity; this helps the wave function factorize to a cluster structure like $`\alpha +d`$ in 6Li or $`\alpha +t`$ in 7Li at large cluster separations. The $`f_{pp}^{LS[n]}(r)`$ is similar to the deuteron (triton) $`f_c(r)`$ in the case of 6Li (7Li).

The $`LS`$ components of the single-particle wave function are given by:

$`|\mathrm{\Phi }_A(LS[n]JMTT_3)_{1234:56\mathrm{}A}=|\mathrm{\Phi }_\alpha (0000)_{1234}{\displaystyle \underset{4<lA}{}}\varphi _p^{LS[n]}(R_{\alpha l})`$ (32)
$`\{[{\displaystyle \underset{4<lA}{}}Y_{1m_l}(\mathrm{\Omega }_{\alpha l})]_{LM_L[n]}\times [{\displaystyle \underset{4<lA}{}}\chi _l(\frac{1}{2}m_s)]_{SM_S}\}_{JM}\times [{\displaystyle \underset{4<lA}{}}\nu _l(\frac{1}{2}t_3)]_{TT_3}.`$
The $`\varphi _p^{LS[n]}(R_{\alpha k})`$ are p-wave solutions of a particle in an effective $`\alpha +N`$ potential and are functions of the distance between the center of mass of the $`\alpha `$ core and nucleon $`k`$; they may be different for different $`LS[n]`$ components. A Woods-Saxon potential well,

$$V_{\alpha N}(r)=V_p^{LS}[1+exp(\frac{rR_p}{a_p})]^1,$$
(33)
where $`V_p^{LS}`$, $`R_p`$, and $`a_p`$ are variational parameters, and a Coulomb term if appropriate is used.

For $`A7`$ nuclei, with three or more p-shell particles there are multiple ways of coupling the orbital angular momentum, $`L`$, spin, $`S`$, and isospin, $`T`$, to obtain the total $`(J^\pi ;T)`$ of the system. The spatial permutation symmetry, denoted by the Young pattern $`[n]`$, is used to enumerate the different possible terms, as discussed in Appendix 1C of Ref. (?). The different possible contributions to $`A`$ = 6–8 nuclei are given in Table 1, with the corresponding spin states of highest spatial symmetry for each nucleus.

After other parameters in the trial function have been optimized, a series of calculations are made in which the $`\beta _{LS[n]}`$ of Eq.(3) may be different in the left- and right-hand-side wave functions to obtain the diagonal and off-diagonal matrix elements of the Hamiltonian and the corresponding normalizations and overlaps. The resulting N$`\times `$N matrices are diagonalized to find the $`\beta _{LS[n]}`$ eigenvectors, using generalized eigenvalue routines because the correlated $`\mathrm{\Psi }_V`$ are not orthogonal. This allows us to project out not only the ground states, but excited states of the same $`(J^\pi ;T)`$ quantum numbers. For example, the ground states of 6Li, 7Li, and 8Li are obtained from 3$`\times `$3, 5$`\times `$5, and 7$`\times `$7 diagonalizations, respectively.

The energy expectation value of Eq.(21) is evaluated using Monte Carlo integration. A detailed technical description of the methods used can be found in Refs. (?, ?). Monte Carlo sampling is done both in configuration space and in the order of operators in the symmetrized product of Eq.(23) by following a Metropolis random walk. The expectation value for an operator $`O`$ is computed with the expression

$$O=\frac{_{p,q}𝑑𝐑\left[\mathrm{\Psi }_p^{}(𝐑)O\mathrm{\Psi }_q(𝐑)/W_{pq}(𝐑)\right]W_{pq}(𝐑)}{_{p,q}𝑑𝐑\left[\mathrm{\Psi }_p^{}(𝐑)\mathrm{\Psi }_q(𝐑)/W_{pq}(𝐑)\right]W_{pq}(𝐑)},$$
(34)
where samples are drawn from a probability distribution, $`W_{pq}(𝐑)`$. The subscripts $`p`$ and $`q`$ specify the order of operators in the left- and right-hand-side wave functions, while the integration runs over the particle coordinates $`𝐑=(𝐫_1,𝐫_2,\mathrm{},𝐫_A)`$. The probability distribution is constructed from the $`\mathrm{\Psi }_P`$ of Eq.(23):

$$W_{pq}(𝐑)=|\mathrm{Re}(\mathrm{\Psi }_{P,p}^{}(𝐑)\mathrm{\Psi }_{P,q}(𝐑))|.$$
(35)
This is much less expensive to compute than using the full wave function of Eq.(22) with its spin-orbit and operator-dependent three-body correlations, but it typically has a norm within 1–2% of the full wave function.

Expectation values have a statistical error which can be estimated by the standard deviation $`\sigma `$:

$$\sigma =\left[\frac{O^2O^2}{N1}\right]^{1/2},$$
(36)
where $`N`$ is the number of statistically independent samples. Block averaging schemes can also be used to estimate the autocorrelation times and determine the statistical error.

The wave function $`\mathrm{\Psi }`$ can be represented by an array of $`2^A\times (_Z^A)`$ complex numbers,

$$\mathrm{\Psi }(𝐑)=\underset{\alpha }{}\psi _\alpha (𝐑)|\alpha ,$$
(37)
where the $`\psi _\alpha (𝐑)`$ are the complex coefficients of each state $`|\alpha `$ with specific third components of spin and isospin. This gives vectors with 96, 1280, 4480, and 14336 complex numbers for 4He, 6Li, 7Li, and 8Li, respectively. The spin, isospin, and tensor operators $`O_{ij}^{p=2,6}`$ contained in the two-body correlation operator $`U_{ij}`$, and in the Hamiltonian are sparse matrices in this basis. For forces that are largely charge-independent, as is the case here, this charge-conserving basis can be replaced with an isospin-conserving basis that has $`N(A,T)=2^A\times I(A,T)`$ components, where

$$I(A,T)=\frac{2T+1}{\frac{1}{2}A+T+1}\left(\begin{array}{c}A\\ \frac{1}{2}A+T\end{array}\right).$$
(38)
This reduces the number of vector elements to 32, 320, 1792, and 7168 for the cases given above — a significant savings. In practice, isospin operators are more expensive to evaluate in this basis, but the overall savings in computation is still large. Furthermore, if for even-$`A`$ nuclei the $`M=0`$ state is computed, only half the spin components need to be evaluated; the other ones can be found by time-reversal invariance.

Expectation values of the kinetic energy and spin-orbit potential require the computation of first derivatives and diagonal second derivatives of the wave function. These are obtained by evaluating the wave function at $`6A`$ slightly shifted positions of the coordinates $`𝐑`$ and taking finite differences, as discussed in Ref. (?). Potential terms quadratic in L require mixed second derivatives, which can be obtained by additional wave function evaluations and finite differences.

The rapid growth in the size of $`\mathrm{\Psi }`$ as $`A`$ increases is the chief limitation, both in computer time and memory requirements, on extending the current method to larger nuclei. The cost of an energy calculation for a given configuration also increases as the number of pairs, $`P=\frac{1}{2}A(A1)`$, because of the number of pair operations required to construct $`\mathrm{\Psi }`$, and roughly as the number of particles, $`A`$, due to the number of wave function evaluations required to compute the kinetic energy by finite differences. These factors are shown in Table 2 for a number of nuclei, with a final overall product showing the cost relative to that for an 8Be calculation. Some initial $`A=9,10`$ calculations have been made, but are at the limit of current computer resources, so while 12C should be reached in a few years, it may be the practical limit for this approach. Calculations of nuclei like 16O or 40Ca, will require other methods such as auxiliary-field diffusion Monte Carlo, which samples the spin and isospin components of the wave function by introducing auxiliary fields, as the plain Monte Carlo samples the spatial part. Also shown in the table is the cost for several pure neutron problems, such as an eight-body drop (with external confining potential) or 14 or 38 neutrons in a box (with periodic boundary conditions) as neutron matter simulations. The 8n drop has been calculated (?, ?) by VMC and GFMC, and initial box simulations have been done for 14 neutrons by GFMC (?) and up to 54 neutrons by AFDMC (?).

A major problem arises in minimizing the variational energy for p-shell nuclei using the above wave functions: there is no variational minimum that gives reasonable rms radii. For example, the variational energy for 6Li is slightly more bound than for 4He, but is not more bound than for separated 4He and 2H nuclei, so the wave function is not stable against breakup into $`\alpha +d`$ subclusters. Consequently, the energy can be lowered toward the sum of 4He and 2H energies by making the wave function more and more diffuse. Such a diffuse wave function would not be useful for computing other nuclear properties, or as a starting point for the GFMC calculation, so the search for variational parameters is constrained by requiring the resulting point proton rms radius, $`r_p`$, to be close to the experimental values for 6Li and 7Li ground states. For other $`A`$ = 6–8 ground states, and all the excited states, the trial functions contain minimal changes to the 6Li and 7Li wave functions, with the added requirement that excited states should not have smaller radii than the ground states. The final step is always the diagonalization of the Hamiltonian in the $`\beta _{LS[n]}`$ mixing parameters.

## 4 GREEN’S FUNCTION MONTE CARLO

The GFMC method (?, ?) projects out the exact lowest-energy state, $`\mathrm{\Psi }_0`$, for a given set of quantum numbers, using $`\mathrm{\Psi }_0=lim_\tau \mathrm{}\mathrm{exp}[(HE_0)\tau ]\mathrm{\Psi }_T`$, where $`\mathrm{\Psi }_T`$ is an initial trial function. If the maximum $`\tau `$ actually used is large enough, the eigenvalue $`E_0`$ is calculated exactly while other expectation values are generally calculated neglecting terms of order $`|\mathrm{\Psi }_0\mathrm{\Psi }_T|^2`$ and higher (?). In contrast, the error in the variational energy, $`E_V`$, is of order $`|\mathrm{\Psi }_0\mathrm{\Psi }_T|^2`$, and other expectation values calculated with $`\mathrm{\Psi }_T`$ have errors of order $`|\mathrm{\Psi }_0\mathrm{\Psi }_T|`$. In the following we present a brief overview of modern nuclear GFMC methods; much more detail may be found in Refs. (?, ?).

We start with the $`\mathrm{\Psi }_T`$ of Eq.(30) and define the propagated wave function $`\mathrm{\Psi }(\tau )`$

$`\mathrm{\Psi }(\tau )=e^{(HE_0)\tau }\mathrm{\Psi }_T=\left[e^{(HE_0)\mathrm{}\tau }\right]^n\mathrm{\Psi }_T,`$ (39)
where we have introduced a small time step, $`\tau =n\mathrm{}\tau `$; obviously $`\mathrm{\Psi }(\tau =0)=\mathrm{\Psi }_T`$ and $`\mathrm{\Psi }(\tau \mathrm{})=\mathrm{\Psi }_0`$. The $`\mathrm{\Psi }(\tau )`$ is represented by a vector function of $`𝐑`$ using Eq.(37), and the Green’s function, $`G_{\alpha \beta }(𝐑,𝐑^{})`$ is a matrix function of $`𝐑`$ and $`𝐑^{}`$ in spin-isospin space, defined as

$$G_{\alpha \beta }(𝐑,𝐑^{})=𝐑,\alpha |e^{(HE_0)\mathrm{}\tau }|𝐑^{},\beta .$$
(40)
It is calculated with leading errors of order $`(\mathrm{}\tau )^3`$. Omitting the spin-isospin indices $`\alpha `$, $`\beta `$ for brevity, $`\mathrm{\Psi }(𝐑_n,\tau )`$ is given by

$$\mathrm{\Psi }(𝐑_n,\tau )=G(𝐑_n,𝐑_{n1})\mathrm{}G(𝐑_1,𝐑_0)\mathrm{\Psi }_T(𝐑_0)𝑑𝐑_{n1}\mathrm{}𝑑𝐑_1𝑑𝐑_0,$$
(41)
with the integral being evaluated stochastically.

The short-time propagator should allow as large a time step $`\mathrm{}\tau `$ as possible, because the total computational time for propagation is proportional to $`1/\mathrm{}\tau `$. Earlier calculations (?, ?, ?) used the propagator obtained from the Feynman formulae in which the kinetic and potential energy terms of $`H`$ are separately exponentiated. The main error in this approximation to $`G_{\alpha ,\beta }`$ comes from terms in $`e^{H\mathrm{}\tau }`$ having multiple $`v_{ij}`$, like $`v_{ij}Kv_{ij}(\mathrm{}\tau )^3`$, where $`K`$ is the kinetic energy, which can become large when particles $`i`$ and $`j`$ are very close due to the large repulsive core in $`v_{ij}`$. This requires a rather small $`\mathrm{}\tau 0.1`$ GeV-1.

It has been found in studies of bulk helium atoms (?) that including the exact two-body propagator allows much larger time steps. This short-time propagator is

$`G_{\alpha \beta }(𝐑,𝐑^{})`$ $`=`$ $`G_0(𝐑,𝐑^{})\alpha |\left[𝒮{\displaystyle \underset{i<j}{}}{\displaystyle \frac{g_{ij}(𝐫_{ij},𝐫_{ij}^{})}{g_{0,ij}(𝐫_{ij},𝐫_{ij}^{})}}\right]|\beta ,`$ (42)
where

$`G_0(𝐑,𝐑^{})`$ $`=`$ $`𝐑|e^{K\mathrm{}\tau }|𝐑^{}=\left[\sqrt{{\displaystyle \frac{m}{2\pi \mathrm{}^2\mathrm{}\tau }}}\right]^{3A}\mathrm{exp}\left[{\displaystyle \frac{(𝐑𝐑^{})^2}{2\mathrm{}^2\mathrm{}\tau /m}}\right],`$ (43)
$`g_{ij}`$ is the exact two-body propagator,

$`g_{ij}(𝐫_{ij},𝐫_{ij}^{})=𝐫_{ij}|e^{H_{ij}\mathrm{}\tau }|𝐫_{ij}^{},`$ (44)
and $`g_{0,ij}`$ is the free two-body propagator. All terms containing any number of the same $`v_{ij}`$ and $`K`$ are treated exactly in this propagator, as we have included the imaginary-time equivalent of the full two-body scattering amplitude. It still has errors of order $`(\mathrm{}\tau )^3`$, however they are from commutators of terms like $`v_{ij}Kv_{ik}(\mathrm{}\tau )^3`$ which become large only when both pairs $`ij`$ and $`ik`$ are close. Because this is a rare occurrence, a five times larger time step, $`\mathrm{}\tau 0.5`$ GeV-1, can be used (?).

Including the three-body forces and the $`E_0`$ in Eq.(40), the complete propagator is given by

$`G_{\alpha \beta }(𝐑,𝐑^{})`$ $`=`$ $`e^{E_o\mathrm{}\tau }G_0(𝐑,𝐑^{})\mathrm{exp}[{\displaystyle (V_{ijk}^R(𝐑)+V_{ijk}^R(𝐑^{}))\frac{\mathrm{}\tau }{2}}]`$ (45)
$`\alpha |I_3(𝐑)|\gamma \gamma |\left[𝒮{\displaystyle \underset{i<j}{}}{\displaystyle \frac{g_{ij}(𝐫_{ij},𝐫_{ij}^{})}{g_{0,ij}(𝐫_{ij},𝐫_{ij}^{})}}\right]|\delta \delta |I_3(𝐑^{})|\beta ,`$
where

$`I_3(𝐑)=\left[1{\displaystyle \frac{\mathrm{}\tau }{2}}{\displaystyle V_{ijk}^\pi (𝐑)}\right].`$ (46)
and $`V_{ijk}^\pi =V_{ijk}^{2\pi }+V_{ijk}^{3\pi }`$ represents, in general, all non-central $`V_{ijk}`$ terms. With the exponential of $`V_{ijk}^\pi `$ expanded to first order in $`\mathrm{}\tau `$, there are additional error terms of the form $`V_{ijk}^\pi V_{i^{}j^{}k^{}}^\pi (\mathrm{}\tau )^2`$. However, they have negligible effect because $`V_{ijk}^\pi `$ has a magnitude of only a few MeV.

Quantities of interest are evaluated in terms of a “mixed” expectation value between $`\mathrm{\Psi }_T`$ and $`\mathrm{\Psi }(\tau )`$:

$`O_{\mathrm{Mixed}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Psi }_T|O|\mathrm{\Psi }(\tau )}{\mathrm{\Psi }_T|\mathrm{\Psi }(\tau )}}`$ (47)
$`=`$ $`{\displaystyle \frac{𝑑𝐏_n\mathrm{\Psi }_T^{}(𝐑_n)OG(𝐑_n,𝐑_{n1})\mathrm{}G(𝐑_1,𝐑_0)\mathrm{\Psi }_T(𝐑_0)}{𝑑𝐏_n\mathrm{\Psi }_T^{}(𝐑_n)G(𝐑_n,𝐑_{n1})\mathrm{}G(𝐑_1,𝐑_0)\mathrm{\Psi }_T(𝐑_0)}},`$
where $`𝐏_n=𝐑_0,𝐑_1,\mathrm{},𝐑_n`$ denotes the ‘path’, and $`d𝐏_n=d𝐑_0d𝐑_1\mathrm{}d𝐑_n`$ with the integral over the paths being carried out stochastically. The desired expectation values would, of course, have $`\mathrm{\Psi }(\tau )`$ on both sides; by writing $`\mathrm{\Psi }(\tau )=\mathrm{\Psi }_T+\delta \mathrm{\Psi }(\tau )`$ and neglecting terms of order $`[\delta \mathrm{\Psi }(\tau )]^2`$, we obtain the approximate expression

$`O(\tau )={\displaystyle \frac{\mathrm{\Psi }(\tau )|O|\mathrm{\Psi }(\tau )}{\mathrm{\Psi }(\tau )|\mathrm{\Psi }(\tau )}}O(\tau )_{\mathrm{Mixed}}+[O(\tau )_{\mathrm{Mixed}}O_T],`$ (48)
where $`O_T`$ is the variational expectation value. More accurate evaluations of $`O(\tau )`$ are possible (?), essentially by measuring the observable at the mid-point of the path. However, such estimates require a propagation twice as long as the mixed estimate and require separate propagations for every $`O`$ to be evaluated. The nuclear calculations published to date use the approximation of Eq.(48).

A special case is the expectation value of the Hamiltonian. The $`H(\tau )_{\mathrm{Mixed}}`$ can be re-expressed as (?)

$`H(\tau )_{\mathrm{Mixed}}={\displaystyle \frac{\mathrm{\Psi }_T|e^{(HE_0)\tau /2}He^{(HE_0)\tau /2}|\mathrm{\Psi }_T}{\mathrm{\Psi }_T|e^{(HE_0)\tau /2}e^{(HE_0)\tau /2}|\mathrm{\Psi }_T}}E_0,`$ (49)
since the propagator $`\mathrm{exp}[(HE_0)\tau ]`$ commutes with the Hamiltonian. Thus $`H(\tau )_{\mathrm{Mixed}}`$ approaches $`E_0`$ in the limit $`\tau \mathrm{}`$, and the expectation value obeys the variational principle for all $`\tau `$.

The AV18 interaction contains terms that are quadratic in orbital angular momentum, $`L`$. These terms are, in essence, state- and position-dependent modifications of the mass of the nucleons. If they are included in the calculation of $`g_{ij}`$, then the ratio $`g_{ij}/g_{0,ij}`$ in Eq.(43) becomes unbounded for large $`|𝐫_{ij}|`$ or $`|𝐫_{ij}^{}|`$ and the Monte Carlo statistical error will also be unbounded. Hence the GFMC propagator is constructed with a simpler isoscalar interaction, $`H^{}`$, with a $`v_{ij}^{}`$ that has only eight operator terms, $`[1,𝝈_i𝝈_j,S_{ij},𝐋𝐒][1,𝝉_i𝝉_j]`$, chosen such that it equals the CI part of the full interaction in all S- and P-waves and in the deuteron. The $`v_{ij}^{}`$ is a little more attractive than $`v_{ij}`$, so a $`V_{ijk}^{}`$ adjusted so that $`H^{}H`$ is also used. This ensures the GFMC algorithm will not propagate to excessively large densities due to overbinding. Consequently, the upper bound property applies to $`H^{}(\tau )`$, and $`HH^{}`$ must be evaluated perturbatively.

Another complication that arises in the GFMC algorithm is the “fermion sign problem”. This arises from the stochastic evaluation of the matrix elements in Eq.(47). The $`G(𝐑_i,𝐑_{i1})`$ is a local operator and can mix in the boson solution. This has a (much) lower energy than the fermion solution and thus is exponentially amplified in subsequent propagations. In the final integration with the antisymmetric $`\mathrm{\Psi }_T`$, the desired fermionic part is projected out, but in the presence of large statistical errors that grow exponentially with $`\tau `$. Because the number of pairs that can be exchanged grows with $`A`$, the sign problem also grows exponentially with increasing $`A`$. For $`A8`$, the errors grow so fast that convergence in $`\tau `$ cannot be achieved.

For simple scalar wave functions, the fermion sign problem can be controlled by not allowing the propagation to move across a node of the wave function. Such “fixed-node” GFMC provides an approximate solution which is the best possible variational wave function with the same nodal structure as $`\mathrm{\Psi }_T`$. However, a more complicated solution is necessary for the spin- and isospin-dependent wave functions of nuclei. In the last few years a suitable “constrained path” approximation has been developed and extensively tested, first for condensed matter systems (?) and more recently for nuclei (?). The basic idea of the constrained-path method is to discard those configurations that, in future generations, will contribute only noise to expectation values. If the exact ground state $`|\mathrm{\Psi }_0`$ were known, any configuration for which

$$\mathrm{\Psi }(𝐏_n)^{}\mathrm{\Psi }_0(𝐑_n)=0,$$
(50)
where a sum over spin-isospin states is implied, could be discarded. Here $`𝐏_n`$ designates the complete path, $`(𝐑_0,\mathrm{},𝐑_n)`$ that has led to $`𝐑_n`$. The sum of these discarded configurations can be written as a state $`|\mathrm{\Psi }_d`$, which obviously has zero overlap with the ground state. The $`\mathrm{\Psi }_d`$ contains only excited states and should decay away as $`\tau \mathrm{}`$, thus discarding it is justified. Of course the exact $`\mathrm{\Psi }_0`$ is not known, and so configurations are discarded with a probability such that the average overlap with the trial wave function,

$$\mathrm{\Psi }_d|\mathrm{\Psi }_T=0.$$
(51)
Many tests of this procedure have been made (?) and it usually gives results that are consistent with unconstrained propagation, within statistical errors. However a few cases in which the constrained propagation converges to the wrong energy (either above or below the correct energy) have been found. Therefore a small number, $`n_u=10`$ to 20, unconstrained steps are made before evaluating expectation values. These few unconstrained steps, out of typically 400 total steps, appear to be enough to damp out errors introduced by the constraint, but do not greatly increase the statistical error.

Figure 3 shows the progress with increasing $`\tau `$ of typical constrained GFMC calculations, in this case for various states of 8Be. The values shown at $`\tau =0`$ are the VMC values using $`\mathrm{\Psi }_T`$; VMC values with the best known $`\mathrm{\Psi }_V`$ are about 1.5 MeV below these. The GFMC very rapidly makes a large improvement on these energies; by $`\tau =0.01`$ MeV-1, the $`\mathrm{\Psi }_T`$ energies have been reduced by 6.5 to 8 MeV. This rapid improvement corresponds to the removal of small admixtures of states with excitation energies $`1`$ GeV from $`\mathrm{\Psi }_T`$. Averages over typically the last nine $`\tau `$ values are used as the GFMC energy. The standard deviation, computed using block averaging, of all of the individual energies for these $`\tau `$ values is used to compute the corresponding statistical error. The solid lines show these averages; the corresponding dashed lines show the statistical errors. The g.s., $`1^+`$, and $`3^+`$ states of 8Be have widths less than 250 keV and their $`E(\tau )`$ appear to be converged and constant over the averaging regime. The $`2^+`$ state has a width of 1.5 MeV and the $`4^+`$ state’s width is 3.5 MeV; the $`E(\tau )`$ for the $`2^+`$ state might be converged, but that of the $`4^+`$ state is clearly steadily decreasing. Reliable estimates of the energies of broad resonances requires the use of scattering-state boundary conditions (see Sec. 7) which have not yet been implemented for more than five nucleons.

Tests with different $`\mathrm{\Psi }_T`$ as starting points for the GFMC calculation are described in Ref. (?). The most crucial aspect in choosing $`\mathrm{\Psi }_T`$ is that it have the correct mix of spatial symmetries, i.e., the $`\beta _{LS[n]}`$ admixtures. This is because the limited propagation time is not sufficient for the GFMC algorithm to filter out low-lying excitations with the same quantum numbers. This issue is unimportant in 4He, where the first 0+ excited state is near 20 MeV, but in 6Li, the first 1+ excited state is at only 5.65 MeV; other $`p`$-shell nuclei have similar low-lying excitations. Otherwise, the GFMC algorithm is able to rapidly correct for some very poor $`\mathrm{\Psi }_T`$ whose variational energies are actually positive. However, a good $`\mathrm{\Psi }_T`$ helps to keep the error bars small at larger $`\tau `$. The most efficient balance between speed of construction of $`\mathrm{\Psi }_T`$ and smallest number of samples needed to achieve a given error bar seems to be given by Eq.(30). A number of other tests of the GFMC algorithm and its ability to determine nuclear radii are described in Refs. (?, ?).

The calculations described here are computationally intense, and would not have been possible without the advent of parallel supercomputers. The initial studies (?) of 6Li required $`2,000`$ node hours on an IBM SP1 to propagate 10,000 configurations to $`\tau `$ = 0.06 MeV-1. Today, with improvements in the algorithm and after much effort to optimize the computer codes, the same calculation requires about 15 node hours on a third generation IBM SP. However, a present calculation for 8Li, requiring 10,000 configurations to get a reasonable error bar, takes 1,000 node hours, so forefront computer resources remain essential for this program.

## 5 ENERGY RESULTS

There are a number of accurate many-body methods for evaluating the binding energy of the s-shell nuclei, 3H, 3He, and 4He, using realistic nuclear forces. These include Faddeev in configuration space (?) (FadR), Faddeev and Faddeev-Yakubovsky in momentum space (?) (FadQ), and hyperspherical harmonics (?) (HH), correlated hyperspherical harmonics (?) (CHH) and pair-correlated hyperspherical harmonics (?) (PHH). Some results of these methods for the AV18 and AV18/UIX Hamiltonians are shown in Table 3. We observe that the VMC upper bounds are generally 1.5–2.0% above the GFMC energies, which in turn are 0.25–1.0% above the Faddeev results. The disagreement between GFMC and the other methods may be attributable to the fact that GFMC is propagated with an $`H^{}`$ and a small piece of the Hamiltonian is computed in perturbation. At this fine level of comparison, one needs to worry about how features such as $`T=\frac{3}{2}`$ admixtures and $`np`$ mass differences in the trinucleon ground state are treated in each calculation.

Table 3 also shows the necessity of including (at least) a three-nucleon potential in the Hamiltonian in order to reproduce the A=3,4 experimental binding energies. This is true for all the other modern $`NN`$ potentials, as shown in Refs. (?, ?). The local potentials AV18, Reid93, and Nijm II, give very similar results, while the slightly nonlocal Nijm I gives 1–3% more binding, and the more nonlocal CD-Bonn gives 5–8% more binding, or about halfway between the local potential values and experiment. At present, no two- plus three-nucleon potential combination gives an exact fit to both trinucleon and 4He energies, but several combinations, like AV18/UIX, come quite close. It also appears that it may be easier to fit 3He and 4He simultaneously, rather than 3H and 4He, perhaps because of the better data constraints on the $`pp`$ interaction compared to the $`nn`$ interaction (?). In principle, there should also be four-body forces, but their contribution must be quite small, and it is impossible to unambiguously identify a need for such terms until a more thorough survey of possible three-nucleon potentials is made and the discrepancies between the various many-body methods are resolved.

Figure 4 compares VMC (using the full $`\mathrm{\Psi }_V`$ of Eq.(22)) and GFMC calculations for the AV18/UIX Hamiltonian, and also shows experimental energies. All GFMC energies in this review are from Ref. (?); VMC energies are from Refs. (?, ?). All particle-stable or narrow-width ($`\mathrm{\Gamma }<150`$ keV) states for $`4A8`$ are shown, except isobaric analogs. We see that the VMC calculations get progressively worse in the p-shell compared to the s-shell, being on the order of 10–15% above the final GFMC results, for nuclei with $`NZ`$. In absolute terms, $`\mathrm{\Psi }_V`$ misses roughly an extra 1.5 MeV of binding for each p-shell nucleon that is added. The VMC results also fail to reproduce important qualitative features of the GFMC calculations; for example 6Li is stable against breakup into $`\alpha `$+d for this Hamiltonian but the VMC calculation shows it unbound by $``$2 MeV (the black dashed lines show the indicated thresholds computed using the sub-cluster energies appropriate to each calculation or experiment). Clearly, there is some significant feature of light p-shell nuclei that is not yet included in the trial wave functions.

Table 4 shows the GFMC energy values for several Hamiltonians. The table and Fig. 4 show that the AV18/UIX Hamiltonian, which gives excellent energies, compared to experiment, for the s-shell nuclei, significantly underbinds the p-shell nuclei. The Li isotopes are stable against breakup into subclusters with AV18/UIX, but progressively more underbound as $`A`$ increases. There is also an isospin problem, in that the He isotopes are even further off from the experimental values, and in fact do not show stability against breakup into subclusters.

The new Illinois three-nucleon potentials were constructed to solve the p-shell binding problems. Figure 5 and Table 4 compare GFMC calculations for AV18 with no $`V_{ijk}`$ and AV18/IL2 to the same experimental energies shown in Fig. 4. The AV18/IL2 Hamiltonian does a good job of reproducing these energies; the rms deviation from experiment for these levels is only 360 keV, while it is 2.3 MeV for AV18/UIX. The AV18 values with no $`V_{ijk}`$ show the large contribution that three-nucleon potentials make to these binding energies; for $`A=8`$ the IL2 increases the binding energy by more than 10 MeV.

Table 5 shows some splittings for states that in a simple shell-model picture have the same $`L`$ and $`S`$ but different total $`J`$. In all cases the lowest state for each $`J`$ is used. Some of the states, particularly the 5He states, have large experimental widths; calculations such as those described in Sec. 7 would be more reliable. The splittings are computed for the AV18, AV18/UIX, and AV18/IL2 Hamiltonians and are compared to experimental values. The AV18 with no $`V_{ijk}`$ significantly underpredicts all the splittings except the $`\frac{1}{2}^{}\frac{3}{2}^{}`$ doublet in 7Li; however in this case the relatively large statistical error for the small splitting makes any conclusion difficult (the computed splittings are the difference of two independent GFMC calculations whose statistical errors must be added in quadrature). The AV18/UIX Hamiltonian substantially improves the 5He splitting, and a VMC study showed that an earlier Urbana $`V_{ijk}`$ significantly increased the computed splitting in 15N, resulting in much better agreement with experiment (?). However the other splittings in the table are still underpredicted with AV18/UIX. The AV18/IL2 results in good predictions of all the splittings except for the 1$`{}_{}{}^{+}2^+`$ doublet in 8Li which is overpredicted.

A detailed breakdown of the GFMC ground-state energies for AV18/IL2 into kinetic, two- and three-nucleon interactions is given in Table 6. Because of the extrapolation of the mixed expectation values, Eq.(48), these components are not as accurate as the total energy, and do not add up to the full amount; indeed the sum of the individual kinetic and potential energies differs from $`H`$ by the same amount that the GFMC has improved the $`\mathrm{\Psi }_T`$ energy. Nevertheless, they give a good idea about the relative size of the different terms involved. There is a big cancellation between the kinetic and two-body potential terms. Consequently, while $`V_{ijk}`$ is less than 8% of $`v_{ij}`$ in magnitude, its expectation value is up to 50% of the net binding (however the net effect of the IL2 $`V_{ijk}`$, defined as the difference of the binding energies computed without and with $`V_{ijk}`$, is at most 30%, as can be seen in Table 4). This difference between expectation value and net effect is due to the large change that $`V_{ijk}`$ induces in $`K+v_{ij}`$.

Among the subcomponents of $`v_{ij}`$, the one-pion exchange dominates, being 70–80% of the total $`v_{ij}`$ for $`A3`$ nuclei. Similarly, the two-pion exchange is the dominant component of $`V_{ijk}`$. The three-pion exchange term is small and repulsive in s-shell nuclei, but attractive in the p-shell nuclei, which have $`T=\frac{3}{2}`$ triples. It is this term which results in the large improvement of the Illinois models over the Urbana models for p-shell nuclei; however its contributions are always less than 15% of the two-pion $`V_{ijk}`$. Finally, we note that the electromagnetic $`v_{ij}^\gamma `$ is dominated by the Coulomb interaction between protons, but about 17% (8%) of its total contribution comes from the magnetic moment and other terms in He (Li) isotopes.

Energy differences among isobaric analog states are probes of the charge-indepen-dence-breaking parts of the Hamiltonian. The understanding of these “Nolen-Schiffer energies” has been a theoretical problem for 30 years (?). The energies for a given isomultiplet of states can be expanded as

$$E_{A,T}(T_z)=\underset{n2T}{}a_n(A,T)Q_n(T,T_z),$$
(52)
where $`Q_0=1`$, $`Q_1=T_z`$, and $`Q_2=\frac{1}{2}(3T_z^2T^2)`$ are isoscalar, isovector, and isotensor terms (?). The isovector and isotensor coefficients $`a_n(A,T)`$, and various contributions to them, are given in Table 7. These were calculated as expectation values in the $`T_z=T`$ GFMC wave functions for AV18/IL2. The individual terms are: 1) the effect of the neutron-proton mass difference on the kinetic energy, $`K^{CSB}`$; 2) the proton-proton Coulomb potential, including the AV18 form factor, $`v_{C1}(pp)`$; 3) all other electromagnetic terms such as vacuum polarization and magnetic-moment terms, $`v^{\gamma ,R}`$; and 4) the strong-interaction contributions, $`v^{CSB}`$ which contributes to the isotensor coefficients and $`v^{CD}`$ which contributes to the isovector coefficients.

The 3H–3He mass difference is 757 keV with the AV18/IL2 Hamiltonian, in excellent agreement with the experimental value of 764 keV. The bulk of the difference is the $`pp`$ Coulomb energy, but 108 keV comes from the other terms. Because of its very long range, the $`v_{C1}(pp)`$ dominates the $`a_1(A,T)`$ of heavier nuclei even more and small errors in the rms radius of the nucleus can result in changes in the $`v_{C1}(pp)`$ contribution that can totally mask the effects of the other, much smaller, terms.

Two (2+;0+1) states occur very close together in the spectrum of 8Be at 16.6 and 16.9 MeV excitation; these isospin-mixed states come from blending the (2+;1) isobaric analog of the 8Li ground state with the second (2+;0) excited state. There are also fairly close (1+;0,1) and (3+;0,1) pairs at slightly higher energies in the 8Be spectrum. The isospin-mixing matrix elements that connect these pairs of states,

$$E_{01}(J)=\mathrm{\Psi }(J^+;0)|H|\mathrm{\Psi }(J^+;1),$$
(53)
have been computed using VMC wave functions for the AV18/UIX Hamiltonian. Results for the $`E_{01}(J)`$ are given in Table 8. The experimental values are determined from the observed decay widths and energies (?). The dominant contribution, from the Coulomb potential, typically accounts for less than half of the matrix element. We see that the remaining part of the electromagnetic interaction and the strong CSB interaction can provide a significant boost, although the experimental mixing is still underpredicted by $``$20%. It appears that these mixing elements are a more sensitive test of the small CSB components of the Hamiltonian than are the isobaric analog energy differences.

## 6 DENSITY AND MOMENTUM DISTRIBUTIONS

The one- and two-nucleon distributions of light $`p`$-shell nuclei are interesting in a variety of experimental settings. For example, the Borromean 6He and 8He nuclei are popular candidates for study as ‘halo’ nuclei whose last neutrons are weakly bound. In addition, the polarization densities of 6Li and 7Li are important because of possible applications in polarized targets.

Point proton rms radii and quadrupole moments are shown in Table 9, as computed by GFMC for the AV18/IL2 model, along with experimental values. The experimental charge radii have been converted to point proton radii by removing the proton and neutron $`r^2`$ of 0.743 and –0.116 fm2, respectively. The radii are in good agreement with experiment, thanks at least in part to the fact that the AV18/IL2 model reproduces the binding energies of these nuclei very well. The quadrupole moments have been calculated in impulse approximation; two-body charge contributions are expected to provide only a few percent correction. The agreement with experiment is again fairly good, with the exception of the 6Li quadrupole moment, which involves a delicate cancellation between the contributions from the deuteron quadrupole moment and the D-wave part of the $`\alpha `$+$`d`$ relative wave function. In general, quadrupole moments are difficult to calculate accurately with quantum Monte Carlo methods because they are dominated by the long-range parts of the wave functions, which contribute very little to the total energy that VMC and GFMC are both designed to optimize.

Calculated and experimental isoscalar and isovector magnetic moments are shown in Table 10, as defined with the convention used in Eq.(52); thus the isovector values for the $`T=\frac{1}{2}`$ cases are $``$2 times those often quoted in the literature. (For the $`A=8`$, $`T=1`$ nuclei the experimental isoscalar and isovector moments are obtained from the sum and difference of the values for B and Li, since the magnetic moment of the $`T=1`$ $`J^\pi =2^+`$ state in Be is not measured.) The moments were calculated as expectation values in the GFMC wave functions with $`T_z=T`$ in impulse approximation. These impulse isoscalar moments are quite close to experiment and it is expected that two-body current contributions are small (?). For isovector moments, however, there can be significant corrections from meson-exchange contributions. We note that the corrections needed for the $`A`$=7,8 isovector moments are the same sign and magnitude as for $`A`$=3, so an exchange-current model that fixes the s-shell (?) is likely to work for the light p-shell nuclei also.

In Fig. 6, we present proton and matter (proton$`+`$neutron) densities for the stable helium isotopes, as calculated with GFMC for AV18/IL2. Nucleon densities are calculated as simple $`\delta `$-function expectation values, with possible spin and/or isospin projectors: for example, the proton density is given by

$$\rho _p(r)=\frac{1}{4\pi r^2}\mathrm{\Psi }|\underset{i}{}\frac{1+\tau _{zi}}{2}\delta (r|𝐫_i𝐑_{cm}|)|\mathrm{\Psi }.$$
(54)
The blue symbols and curves show results for 6He, red ones for 8He; open symbols give the proton distributions which are also interpreted as the alpha “core” density, full symbols are the total matter densities. It can be seen that, as more neutrons are added, the tails of the matter distributions broaden considerably because of the relatively weak binding of the p-shell neutrons. In addition, the central neutron and proton densities decrease rather dramatically. This effect does not necessarily require any changes to the $`\alpha `$ core, but can be understood at least partially from the fact that the $`\alpha `$ no longer sits at the center of mass of the entire system. The motion relative to the center of mass spreads out the mass distribution relative to that of 4He.

The figure also shows two attempts to extract the densities of 6,8He from scattering of beams of these short-lived nuclei from proton targets (?). The analysis of Alkhazov, et al. (dashed curves) used Glauber theory and assumed proton (core) and matter density distributions which were varied to fit the measured cross sections. The later work of Al-Khalili and Tostevin (?) (dot-dashed curves) improved on this by doing the Glauber theory using model wave functions which contained correlations between the valence neutrons. The GFMC calculations definitely prefer the latter analysis.

The proton-proton distribution function is defined by

$$\rho _{pp}(r)=\frac{1}{4\pi r^2}\mathrm{\Psi }|\underset{i<j}{}\frac{1+\tau _{zi}}{2}\frac{1+\tau _{zj}}{2}\delta (r|𝐫_i𝐫_j|)|\mathrm{\Psi }.$$
(55)
These distributions are directly related to the Coulomb sum measured in inclusive longitudinal electron scattering; such measurements in 3He have been used to put constraints on the $`\rho _{pp}(r_{ij})`$, and realistic calculations agree with the experimental results(?). The behavior of $`\rho _{pp}(r)`$ at short distances is largely determined by the repulsive core of the $`NN`$ potential and is nearly independent of the nucleus, but at larger distances it is determined by the size of the nucleus.

It is interesting to compare $`\rho _{pp}`$ for 4He to that of 6He and 8He. These distribution functions are shown in Fig. 7, again calculated with GFMC for the AV18/IL2 model. These nuclei each have just one $`pp`$ pair which presumably is in the “alpha core” of 6,8He. Unlike the one-body densities, these distributions are not sensitive to center of mass effects, and thus if the alpha core of 6,8He is not distorted by the surrounding neutrons, all three $`\rho _{pp}`$ distributions in the figure should be the same. We see that the $`pp`$ distribution spreads out slightly with neutron number in the helium isotopes, with an increase of the pair rms radius of approximately 4% in going from 4He to 6He, and 8% to 8He. While this could be interpreted as a swelling of the alpha core, it might also be due to the charge-exchange ($`\tau _i\tau _j`$) correlations which can transfer charge from the core to the valence nucleons. Since these correlations are rather long-ranged, they can have a significant effect on the $`pp`$ distribution. VMC calculations of 4He with wave functions modified to give $`\rho _{pp}`$ distributions close to those of 6,8He suggest that the alpha cores of 6,8He are excited by $`80`$ and $`350`$ keV, respectively.

In general, VMC calculations give one-body densities very similar to the GFMC results, although the two-body densities may differ by up to 10% in the peak, with the GFMC having a somewhat sharper structure. Because the VMC wave functions are simpler and less expensive to construct they have been used in a number of applications where the single-particle structure is dominant. All these calculations were carried out with wave functions for the AV18/UIX Hamiltonian, but we expect that GFMC calculations with the improved AV18/IL2 model will not qualitatively alter the results. (Remember that the VMC wave functions are constrained to reproduce available experimental rms radii.)

VMC calculations of elastic and inelastic electromagnetic form factors for 6Li are shown in Fig. 8. These have been made in impulse approximation (IA) and with meson-exchange contributions (MEC) to the charge and density operators (?). The comparison with data for the elastic longitudinal form factor, $`F_L(q^2)`$, is excellent, as is the E2 transition to the 3+ first excited state. The MEC corrections are small, but stand out noticeably in the first minimum where they significantly improve the fit to data. The elastic transverse form factor, $`F_T(q^2)`$, is good up to the first zero, but is a little too large in the region of the second maximum. The M1 transition to the 0+ isobaric analog state is also reproduced reasonably well. This kind of quantitative agreement with data has not been achieved in the past with either shell model or $`\alpha `$+$`d`$ cluster wave functions. In fact, shell model calculations normally require the introduction of effective charges, typically adding a charge of $`0.5e`$ to both the proton and neutron, to obtain reasonable transition rates (?). No effective charges are used in the VMC calculations.

While electron scattering from nuclei is primarily sensitive to the proton distributions in nuclei, $`\pi ^{}`$ scattering is most sensitive to the neutron distributions. Pion inelastic scattering at medium energy ($`80E_{lab}300`$ MeV) is dominated by strong absorption due to excitation of the $`\mathrm{\Delta }`$ resonance, and is well described in distorted-wave impulse approximation (DWIA). An analysis of meson factory data for p-shell nuclei based on the Cohen-Kurath shell model (?) showed that reasonable agreement with data could be achieved if the quadrupole excitation operator was enhanced by a factor $``$ 1.75–2.5 (?). Figure 9 shows a calculation of the differential cross sections for 7Li($`\pi ,\pi ^{}`$) scattering to the first three excited states using VMC transition densities as input, with no enhancement factors (?). Solid red lines show the full results, which are in good agreement with the data, while the dashed blue lines give the contribution coming from protons only. Because of the strong isospin dependence of the $`\pi N`$ scattering $`t`$-matrix, this reaction can be a stringent test of the neutron densities predicted by the quantum Monte Carlo calculations.

The VMC wave functions for the AV18/UIX model have also been used to calculate single-nucleon momentum distributions in many nuclei, and a variety of cluster-cluster overlap wave functions, such as $`dp|t`$, $`dd|\alpha `$, and $`\alpha d|^6\mathrm{Li}`$ (?, ?). Recently the overlaps $`^6\mathrm{He}(J^\pi )+p(\mathrm{}_j)|^7\mathrm{Li}`$ for all possible p-shell states in 6He were studied (?). The spectroscopic factors obtained from these overlaps are 0.41 to the ground state of 6He and 0.19 to the 2+ first excited state. These factors are significantly smaller than the predictions of the Cohen-Kurath shell model (?), which gives values of 0.59 and 0.40, respectively. The CK shell model requires that the possible 6He$`+p`$ states sum to unity within the p-shell, whereas in the VMC calculation, the correlations in the wave function push significant strength to higher momenta that cannot be represented as a 6He state plus p-wave proton. The VMC overlaps were used as input to a Coulomb DWIA analysis of recent 7Li$`(e,e^{}p)^6`$He data taken at NIKHEF (?) and found to give an excellent fit to the data, as shown in Fig. 10.

The VMC wave functions are based on one-body parts that have a shell-model structure, namely four nucleons in an $`\alpha `$ core coupled to $`(A4)`$ one-body $`\mathrm{}=1`$ wave functions. However the low-lying states of 8Be exhibit a rotational spectrum and are believed to be well approximated as two $`\alpha `$’s rotating around their common center of mass. It is possible to recover this picture from the VMC wave functions by a modified Monte Carlo density calculation (?).

The standard Monte Carlo method for computing one-body densities, $`\rho (𝐫)`$, is to make a random walk that samples $`|\mathrm{\Psi }(𝐫_1,𝐫_2,\mathrm{},𝐫_A)|^2`$ and to bin $`𝐫_1,𝐫_2,\mathrm{},𝐫_A`$ for each configuration in the walk. The density is then proportional to the number of samples in each bin. In the case of a $`J`$ = 0 nucleus, this “laboratory” density will necessarily be spherically symmetric.

The intrinsic density in body-fixed coordinates can be approximated by computing the moment of inertia matrix, $``$, of the $`A`$ positions for each configuration. The eigenvalues and eigenvectors of $``$, are found and a rotation to those principal axes is made. The resulting $`𝐫_1^{},𝐫_2^{},\mathrm{},𝐫_A^{}`$ is then binned. The eigenvector with the largest eigenvalue is chosen as the $`𝐳^{}`$ axis. This procedure will not produce a spherically symmetric distribution, even if there is no underlying deformed structure, because almost every random configuration will have principal axes of different lengths and the rotation will always orient the longest principal axis in the $`𝐳^{}`$ direction. However no artificial structure is introduced (?).

When the above procedure is applied to the three lowest 8Be states, a dramatic intrinsic structure is revealed, as shown in Fig. 11 for the ground state. The figure shows contours of constant density plotted in cylindrical coordinates. The left side of the figure shows the standard, lab frame, density calculation. For the $`J`$ = 0 ground state, this is spherically symmetric as shown. The right side of the figure shows the intrinsic density. It is clear that the intrinsic density has two peaks, with the neck between them having only one-third the peak density; we regard these as two $`\alpha `$’s. This assignment is strengthened by making the same construction for the $`J=2^+`$ and $`4^+`$ states; although the laboratory densities for these states (in $`M=J`$ states) are quite different, the intrinsic densities are, within statistical errors, the same as the $`J=0`$ intrinsic density (?).

If the $`0^+`$, $`2^+`$, and $`4^+`$ states are generated by rotations of a common deformed structure, then their electromagnetic moments and transition strengths should all be related to the intrinsic moments which can be computed by integrating over the projected body-fixed densities. This is explored in Ref. (?) and a generally consistent picture is shown.

## 7 LOW-ENERGY NUCLEON-NUCLEUS SCATTERING

The calculations presented so far have treated resonant nuclear states as if they were particle stable; that is the VMC trial wave functions decay exponentially at large distances and the GFMC propagation does not impose a scattering-state boundary condition. As was shown in Fig. 3, this GFMC propagation converges for unstable states that are narrow, but for wide states, such as 8Be(4+), the $`E(\tau )`$ keeps decreasing for increasing $`\tau `$. Such states should be computed using scattering-wave boundary conditions.

The 5He($`\frac{1}{2}^{}`$, $`\frac{3}{2}^{}`$) states have been studied in VMC (?) and GFMC (?) using scattering-wave boundary conditions. The VMC calculations use a trial function that is a straight-forward generalization of the $`\mathrm{\Psi }_V`$ given in Sec. 3. In it the $`\mathrm{\Psi }_J`$ contains one p-wave single-particle function that goes to zero at a specified (large) radius, $`R_n`$. The variational energy is minimized subject to this boundary condition. The phase shift at the resulting scattering energy, $`E_s=E(^5\text{He}(R_n))E(^4\text{He})`$ is then obtained from

$$\mathrm{tan}(\delta _l)=\frac{j_l(kR_n)}{n_l(kR_n)},$$
(56)
where $`k=(2\mu E_s)^{1/2}/\mathrm{}`$, $`\mu `$ is the $`\alpha `$+$`n`$ reduced mass, and $`j_l`$ and $`n_l`$ are the spherical Bessel functions (in this case $`l=1`$). It is possible to generalize this formulation by specifying a logarithmic-derivative at $`R_n`$, rather than having the wave function go to zero there. This method requires the computation of the energy difference, $`E_s`$. If this is done with two different VMC calculations, one for 5He with the specified boundary condition, and one for 4He, then the statistical errors of these two calculations must be added in quadrature. Instead one can directly evaluate the energy difference in a VMC random walk that is controlled by the 5He wave function. In practice this gives $`E_s`$ with a smaller statistical error than either of the separate statistical errors (?). The GFMC calculations start with the scattering-state trial function and use a modified propagator that preserves the boundary condition (?). Both methods result in a phase shift and energy corresponding to a given boundary condition. By repeating the calculations for different boundary conditions, the phase shift as a function of energy can be mapped out.

The VMC calculations of Ref. (?) used older two- and three-nucleon potentials than have been used for the other calculations reported in this review. They reproduced the qualitative features of the experimental $`\alpha `$+$`n`$ $`\frac{1}{2}^{}`$ and $`\frac{3}{2}^{}`$ phase shifts, but, in particular, gave a spin-orbit splitting about 1 MeV too small. The GFMC calculations (?) used a different, but still somewhat old Hamiltonian and obtained values that are closer to experiment. These calculations need to be repeated with the modern Hamiltonians described in this review, and for other systems such as $`\alpha `$+$`d`$, 6He+$`n`$, and $`\alpha `$+$`\alpha `$.

## 8 ASTROPHYSICAL ELECTROWEAK REACTIONS

Most of the key nuclear reactions in primordial nucleosynthesis (?) and solar neutrino production (?) involve only the $`A8`$ nuclei. Many of these are electroweak capture reactions that are difficult or impossible to measure in the laboratory. With the high-precision few-body methods, such as PHH and Faddeev, or new effective field theory treatments, the reactions to s-shell final states have been calculated with unprecedented accuracy in the last few years. These include the 1H$`(p,e^+\nu _e)^2`$H (?) and 3He$`(p,e^+\nu _e)^4`$He (?) weak capture reactions, and the 1H$`(n,\gamma )^2`$H (?), 2H$`(n,\gamma )^3`$H and 2H$`(p,\gamma )^3`$He (?) radiative capture reactions.

Recently, the VMC method has been applied to some of the radiative captures to p-shell final states, including 2H($`\alpha ,\gamma )^6`$Li, 3H($`\alpha ,\gamma )^7`$Li, and 3He($`\alpha ,\gamma )^7`$Be (?, ?). As discussed in the previous section, many-nucleon scattering states can be constructed within the quantum Monte Carlo framework, but it has not yet been done for $`A6`$ nuclei. Thus for these first studies, the appropriate electromagnetic matrix elements (primarily E1 and E2) are evaluated between a VMC six- or seven-nucleon final state, and a correlated cluster-cluster scattering state which is not variationally improved.

Since much of the capture takes place at long range, an important ingredient in these calculations is an asymptotically-correct description of the six- or seven-body final state. This feature has to be implemented in the initial variational wave function because, as remarked before, the VMC and GFMC methods find their best wave functions by optimizing the energy, which is not sensitive to the long-range behavior. The asymptotic behavior can be imposed by requiring the one-body wave function in Eq.(3) to satisfy the condition:

$$[\varphi _p^{LS[n]}(r\mathrm{})]^nW_{km}(2\gamma r)/r,$$
(57)
where $`W_{km}(2\gamma r)`$ is the Whittaker function for bound-state wave functions in a Coulomb potential and $`n`$ is the number of p-shell nucleons. Here $`\gamma ^2=2\mu _{4n}B_{4n}/\mathrm{}^2`$, with $`\mu _{4n}`$ and $`B_{4n}`$ the appropriate two-cluster effective mass and binding energy. For 6Li, $`B_{42}=1.47`$ MeV; for 7Li or 7Be, capture can be to the $`\frac{3}{2}^{}`$ ground or $`\frac{1}{2}^{}`$ excited state with corresponding values for $`B_{43}`$.

The initial-state wave functions are taken as elastic-scattering states of the form

$$|\psi _{\alpha \tau };LSJM=𝒜\left\{\varphi _{\alpha \tau }^{JL}(r_{\alpha \tau })Y_{LM_L}(\widehat{𝐫}_{\alpha \tau })\underset{ij}{}G_{ij}|\psi _\alpha \psi _\tau ^{m_S}\right\}_{LSJM},$$
(58)
where curly braces indicate angular momentum coupling, $`𝒜`$ antisymmetrizes between clusters, $`\psi _\alpha `$ is the 4He ground state, and $`\psi _\tau ^{m_S}`$ is the deuteron or trinucleon ground state in spin orientation $`m_S`$. The $`G_{ij}`$ are identity operators if the nucleons $`i`$ and $`j`$ are in the same cluster, else, they are a set of pair correlation operators, including both central and spin-isospin dependent terms, which introduce distortions in each cluster, under the influence of individual nucleons from the other cluster. They are similar to the correlations discussed in Sec. 3, except that they revert to the identity operator at pair separations beyond about 2 fm. The correlations $`\varphi _{\alpha \tau }^{JL}`$ are generated from optical potentials that describe the experimental phase shifts in cluster-cluster scattering; see Refs. (?, ?) for details. Because the $`G_{ij}`$ go to unity, the $`\psi _{\alpha \tau }`$ has the same phase shifts as those generated by the optical potential.

Evaluation of electromagnetic matrix elements between the initial scattering and final bound states can be split into two parts by noting that all of the energy dependence is contained in the relative wave function $`\varphi _{\alpha \tau }^{JL}`$ and the transition operators. Thus, using a technique first applied in a VMC calculation of $`d`$+$`d`$ radiative capture (?), the matrix element for a given scattering partial wave can be written as

$`T_{\mathrm{}}^{LSJ_iJ_f}(q)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑xx^2\varphi _{\alpha \tau }^{J_iL}(x)`$
$`\times `$ $`\psi _A^{J_fm_f}|T_\mathrm{}\lambda (q)𝒜\{\delta (xr_{\alpha \tau })Y_L^{M_L}(\widehat{𝐫}_{\alpha \tau })){\displaystyle \underset{ij}{}}G_{ij}|\psi _\alpha \psi _\tau ^{m_S}\}_{LSJ_iM_i},`$
where $`T_\mathrm{}\lambda `$ denotes any of the standard $`E_{\mathrm{}}`$ or $`M_{\mathrm{}}`$ operators. The integration over all coordinates except $`x`$ can be calculated just once for each partial wave by Monte Carlo sampling, and the result can then be used to compute the full integral for as many energies as desired by recomputing $`\varphi _{\alpha \tau }^{J_iL}(x)`$ only.

The result for the 2H($`\alpha ,\gamma )^6`$Li capture $`S`$-factor is shown in Fig. 12, where it is compared to a collection of direct and indirect data. The $`S`$-factor falls rapidly at low-energy because a pseudo-orthogonality between the initial and final states suppresses normal S-wave capture. Whenever there is such a suppression, it is important to consider higher-order terms that might contribute, including relativistic corrections and two-body charge and current operators (?). This is possible with the VMC calculation because of the fully correlated $`A`$-body wave functions that are used. In this case, a relativistic center-of-energy correction leads to an E1 contribution that dominates at low energy, while the resonance region and above is primarily E2. Results for 3H($`\alpha ,\gamma )^7`$Li, and 3He($`\alpha ,\gamma )^7`$Be can be found in Ref. (?)

This kind of VMC calculation is only a first investigation into the p-shell astrophysical reactions, and should be improved in the future by using the more exact GFMC wave functions for both bound and scattering states, and the Hamiltonians with improved $`NNN`$ potentials.

## 9 NEUTRON DROPS

Neutron-rich nuclei are interesting both because of their importance in various astrophysical contexts, such as the r-process or neutron-star crusts, and the current interest in experiments with radioactive beams. These nuclei are often studied within the framework of mean-field models using Skyrme or other potential models. The parameters of these models are fit to experimentally known binding energies, that is for situations with $`NZ`$. In particular the isospin dependence of the spin-orbit component of such potentials is considered to be not strongly constrained. Neutron drops offer the possibility of theoretical guidance for the isotopic dependence of such parameters.

Neutron drops are systems of interacting neutrons confined in an artificial external well. Eight neutrons form a closed shell and single-particle spin-orbit splittings can be studied in drops of seven or nine neutrons. Pairing energies can be studied if six-neutron drops are also computed. Calculations of systems of seven and eight neutrons interacting with AV18/UIX were used as a basis for comparing Skyrme models of neutron-rich systems with microscopic calculations based on realistic interactions (?). The external one-body well used is a Woods-Saxon:

$$V_1(r)=\underset{i}{}\frac{V_0}{1+\mathrm{exp}[(r_ir_0)/a_0]};$$
(60)
the parameters are $`V_0=20`$ MeV, $`r_0=3.0`$ fm, and $`a_0=0.65`$ fm. Neither the external well nor the total internal potential ($`v_{ij}+V_{ijk}`$) are individually attractive enough to produce bound states of seven or eight neutrons; however the combination does produce binding.

Tables 4-6 show results for the neutron drops with this external well. The $`T=\frac{3}{2}`$ nature of the $`S_\tau ^IS_\sigma ^I`$ term of $`V^{3\pi }`$ results in large contributions in the neutron drops. As a result the seven-neutron drops computed with AV18/IL2 have double the spin-orbit splitting predicted by AV18/UIX. Thus the spin-orbit splitting in neutron-rich systems depends strongly on the Hamiltonian used.

## 10 CONCLUSIONS AND OUTLOOK

Quantum Monte Carlo methods can now be used to obtain accurate (within 2%) energies for light p-shell nuclei up to $`A`$ = 8. Initial calculations of $`A`$=9,10 nuclei have been made by the authors, and studies of $`A`$=11,12 nuclei should be feasible in the next few years by variational and Green’s function Monte Carlo methods. This progress in the nuclear many-body problem is due both to the rapid growth of computational power and the continuing evolution of algorithms. The new auxiliary-field diffusion Monte Carlo method, which samples spins and isospins by auxiliary fields, and space by standard diffusion Monte Carlo, may be the key to doing even larger nuclear systems (?).

Studies of the p-shell nuclei allow us to test nuclear forces in new ways not accessible in s-shell nuclei, particularly the odd partial waves of $`NN`$ scattering and the $`T=\frac{3}{2}`$ triples for $`NNN`$ forces. They also give us many more cases in which to examine charge-dependent and charge-symmetry-breaking interactions. With these calculations of nuclear spectra, we see for the first time that nuclear structure, including both single-particle and clustering aspects, really can be explained directly in terms of bare nuclear forces that fit $`NN`$ data. A crucial ingredient for quantitative agreement is the addition of realistic $`NNN`$ forces, including at least two terms beyond the standard long-range two-pion-exchange potential (?).

Many aspects of nuclear structure and reactions are described quantitatively in these studies. The energy differences between isobaric analog states are explained well once the complete electromagnetic and strong charge-independence-breaking forces deduced from $`NN`$ scattering are included. Charge radii and quadrupole moments agree with experiment, and we expect magnetic moments will also once the important two-body exchange currents are included, as they have been in s-shell nuclei (?). Elastic and transition form factors measured in electron scattering, and transition densities that are tested in pion scattering, are accurately predicted without the introduction of effective charges. Spectroscopic factors are naturally quenched by the correlations in $`A`$-body wave functions and we see that those correlations build up intrinsic cluster structure, as in the case of 8Be.

Initial calculations of nucleon-nucleus scattering states and electroweak capture reactions are encouraging. These two related problems will be a major activity for the quantum Monte Carlo studies over the next several years, with studies of scattering states in nuclei like 7He and 9He, and of the cluster-cluster scattering states that enter into the astrophysical capture reactions. Additional reactions like 7Be$`(p,\gamma )^8`$B and even 8Be$`(\alpha ,\gamma )^{12}`$C should become feasible. Still other projects will include the study of weak decays and, going beyond the p-shell, looking at the unnatural-parity intruder states which start to dip into the low-energy spectrum with the $`A=9`$ nuclei, and can be particle stable by $`A=10`$.

Refining the nuclear Hamiltonian will remain a major aspect of future work. Extensive new $`NN`$ scattering data taken since the Nijmegen partial-wave analysis has increased the database to more than 6000 points, and a new CD-Bonn 2000 potential has been constructed that again achieves a $`\chi ^2/`$datum $`1`$ (?). Also, intensive searches are being made to look for $`NNN`$ force signatures in $`N+d`$ scattering experiments (?, ?), which may well indicate a need for additional terms beyond those contained in the Illinois models. Our ability to test new force models, as they become available, in light nuclei by accurate quantum Monte Carlo methods will continue to be an important tool for nuclear physics.

## Acknowledgments

The authors thank J. Carlson, A. Kievsky, V. R. Pandharipande, R. Schiavilla, J. P. Schiffer and K. E. Schmidt for useful communications and discussions. This work was supported by the U.S. Department of Energy, Nuclear Physics Division, under contract No. W-31-109-ENG-38.