# Progress in Monte Carlo Calculations of Fermi Systems: Normal Liquid 3He

## Abstract

The application of the diffusion Monte Carlo method to a strongly interacting Fermi system as normal liquid <sup>3</sup>He is explored. We show that the fixed-node method together with the released-node technique and a systematic method to analytically improve the nodal surface constitute an efficient strategy to improve the calculation up to a desired accuracy. This methodology shows unambiguously that backflow correlations, when properly optimized, are enough to generate an equation of state of liquid <sup>3</sup>He in excellent agreement with experimental data from equilibrium up to freezing.

Liquid <sup>3</sup>He has been for many years a benchmark in the field of quantum-many body physics. The Fermi statistics of its atoms, combined with the strong correlations induced by the hard core of their interatomic potential, has turned it into a paradigm of strongly-correlated Fermi systems. At zero temperature, an approximate microscopic description has been achieved by means of variational methods, both for the Fermi liquid <sup>3</sup>He  as well as for its bosonic counterpart liquid <sup>4</sup>He . From a Monte Carlo viewpoint, the quantum many-body problem can be tackled in a more ambitious way with the aid of the Green’s function Monte Carlo (GFMC) and the diffusion Monte Carlo (DMC) methods . The variational wave function can be used as an input for the Monte Carlo method which, for boson systems, is able to solve the Schrödinger equation of the $`N`$-body system providing exact results. In Fermi systems as liquid <sup>3</sup>He, the exactness of the method is lost due the involved sign problem that makes a straightforward interpretation of the wave function not possible.

The cancellation methods developed up to now to solve this intricate problem have proved their efficiency in model problems or with very few particles but become unreliable for real many-body systems . In the meantime, the approximate fixed-node (FN)  method has become a standard tool. In the FN-DMC method, the antisymmetry is introduced in the trial wave function used for importance sampling imposing its nodal surface as a boundary condition. This approach provides upper bounds to the exact eigenvalues, the quality of which is related to the accuracy of the nodal surface of the trial wave function. The main drawback of the FN-DMC method is the lack of control over the influence of the imposed nodal surface on the results obtained, not to say the impossibility of properly correcting for such effect. In the present Letter, we come back to this problem using the FN-DMC as a main approach but crucially combined with two auxiliary methods: the released-node (RN) estimation technique  and an analytical method able to enhance the quality of any given nodal surface. The combination of the above methods provides information on the bias due to the imposed nodal surface and a procedure which can evaluate and bring down, in principle to any arbitrary required precision, this influence. This complete program has been applied to the study of liquid <sup>3</sup>He bringing, as we will show, the systematic error under control and down to levels below the current statistical errors. This has allowed a very accurate microscopic calculation of the equation of state of liquid <sup>3</sup>He, including a prediction for the negative-pressure region and the spinodal density.

DMC  and GFMC  calculations have provided up to now the best upper bounds to the ground-state energy of liquid <sup>3</sup>He . These calculations have unambiguously shown the relevance of the Feynman-Cohen-type backflow correlations in the improvement of the energy. However, the energy gain appears too small to recover the experimental data and, what is more conclusive, the density dependence of the pressure which stresses the curvature of the equation of state shows clear differences with the experiment. A conclusion that naturally emerged from those results was that the nodal surface, originated by backflow correlations, is not accurate enough and probably new state-dependent correlations ought to be considered. Our present results prove that backflow correlations, when properly optimized, effectively do provide very accurate nodal surfaces.

In the DMC method, the Schrödinger equation written in imaginary time is translated into a diffusion-like differential equation which can be stochastically solved in an iterative procedure. Specific information on the implementation of the DMC method is given in Ref. . As far as the FN framework  is concerned, the choice of the trial wave function $`\psi (\text{R})`$ used for importance sampling is a key point. The simplest model is the Jastrow-Slater wave function

$$\psi =\psi _\mathrm{J}D_{}D_{},$$
(1)   
with $`\psi _\mathrm{J}=_{i<j}f(r_{ij})`$ a Jastrow wave function and $`D_{}`$ ($`D_{}`$) a Slater determinant of the spin-up (spin-down) atoms with single-particle orbitals $`\phi _{\alpha _i}(𝐫_j)=\mathrm{exp}(i𝐤_{\alpha _i}𝐫_j)`$.

In this variational description (1), the dynamical correlations induced by the interatomic potential are well modelled by the Jastrow factor, and the statistical correlations, implied by the antisymmetry, are introduced with a Slater determinant of plane waves which is the exact wave function of the free Fermi sea. The two factors account well for the dynamical correlations and the Fermi statistics when these effects are independently considered but their product is only a relatively poor approximation for a strongly correlated Fermi liquid. It is well known, from previous variational and GFMC/DMC calculations , that a significant improvement on the Jastrow-Slater model is achieved by introducing backflow correlations in $`\phi _{\alpha _i}(\text{r}_j)`$, a name which is taken from the Feynman-Cohen famous work on the microscopic description of the phonon-roton spectrum in liquid <sup>4</sup>He .

At this point, it becomes essential to set up a method for analytically enhancing a given model. Such a procedure is already contained in the imaginary-time Schrödinger equation. Let us consider a time-dependent wave function $`\varphi (𝐑,t)`$, with $`\varphi (𝐑,t=0)=\psi (𝐑)`$ the initial guess for the trial wave function, satisfying

$$\frac{\varphi (𝐑,t)}{t}=H\varphi (𝐑,t).$$
(2)   
A natural choice for a more accurate trial wave function is obtained solving Eq. (2) at first order in $`t`$. Near the nodes, which is the relevant region to our purposes, one readily captures the main correction to the original $`\psi _\mathrm{A}(𝐑)D_{}D_{}`$ in the form $`\varphi (𝐑,t)=\psi _\mathrm{J}(𝐑)\psi _\mathrm{A}(\stackrel{~}{𝐑})`$ with $`\stackrel{~}{𝐑}=𝐑(t)`$ and $`𝐑/t=D𝐅_\mathrm{J}(𝐑)`$, $`𝐅_\mathrm{J}(𝐑)`$ being the drift force coming from the Jastrow wave function $`\psi _\mathrm{J}(𝐑)`$ and $`D=\mathrm{}^2/(2m)`$. In this form, the new nodal surface is described by the original antisymmetric wave function but with arguments that are shifted due to the effect of dynamical correlations. It is worth noting that this approach generates the Feynman-Cohen backflow in $`\phi _{\alpha _i}(𝐫_j)`$ ($`\psi _\mathrm{A}^{\mathrm{BF}}(𝐑)`$) as a first order correction to the plane-wave orbitals. Recursively, entering with $`\psi _\mathrm{A}^{\mathrm{BF}}(𝐑)`$ the next order correction can be analytically obtained (see Eq. 5).

Finally, once a specific model for the nodal surface has been chosen it is necessary to establish a method to test its quality. This can be accomplished by means of the released-node technique . In the RN approach a superposition of a small boson component in the wave function is allowed, with the primary effect of resetting the nodal surface to the exact position. This is technically accomplished by introducing a positively-defined guiding wave function $`\psi _\mathrm{g}(\text{R})`$ so that the walkers are not confined into a region of definite sign of $`\psi (\text{R})`$. The basic requirements on choosing $`\psi _\mathrm{g}(\text{R})`$ are twofold: proximity to $`|\psi (\text{R})|`$ away from the nodal surface and being positive-defined at the nodes. The choice we have made is

$$\psi _\mathrm{g}(𝐑)=(\psi (𝐑)^2+a^2)^{1/2},$$
(3)   
$`a`$ being a parameter which controls the crossing frequency. Other different choices can be considered but the specific details of $`\psi _\mathrm{g}(\text{R})`$ are not relevant since the RN energy is calculated projecting out its antisymmetric component. The use of $`\psi _\mathrm{g}(\text{R})`$ does not introduce any systematic bias in the RN energies, which approach the exact eigenvalue when $`t_\mathrm{r}\mathrm{}`$ , $`t_\mathrm{r}`$ being the maximum allowed lifetime after the first crossing. However, the variance of the energy grows exponentially with $`t_\mathrm{r}`$ due to the boson component, and thus in general the asymptotic value cannot be obtained. In contrast, what is straightforwardly available is the slope of the energy versus $`t_\mathrm{r}`$ at $`t_\mathrm{r}0`$, which provides a direct measure of the quality of the input nodal surface (the true antisymmetric ground-state wave function would generate a zero slope), and constitutes a means of comparing different trial wave functions. In particular, it provides feedback information on whether the next analytical correction to $`\psi _\mathrm{A}^{\mathrm{BF}}(𝐑)`$ is necessary.

We have applied all the above methodology to the study of normal liquid <sup>3</sup>He at zero temperature. The results reported have been obtained with $`N=66`$ particles, but we have made size checks using also $`N=54`$ and $`N=114`$. In Fermi systems, the kinetic energy includes statistical contributions that show an oscillating behavior with $`N`$. We have observed that this behavior follows very closely that of a discretized Fermi-gas energy, a fact that could be expected since such a term appears explicitly in the local kinetic energy. It is worth noticing that the case $`N=66`$ is specially well suited for MC calculations as the correction amounts only 0.015 K.

As in previous calculations, we use a short-ranged backflow in the form $`\phi _{\alpha _i}(𝐫_j)=\mathrm{exp}(i𝐤_{\alpha _i}\stackrel{~}{𝐫}_j^{\mathrm{BF}})`$, with

$$\stackrel{~}{𝐫}_j^{\mathrm{BF}}=𝐫_j+\lambda _\mathrm{B}\underset{kj}{}\eta (r_{jk})𝐫_{jk},$$
(4)   
and $`\eta (r)=\mathrm{exp}(((rr_\mathrm{B})/\omega _\mathrm{B})^2)`$. The two-body correlation factor has been chosen of McMillan type, $`f(r)=\mathrm{exp}(0.5(b/r)^5)`$, and the pairwise HFD-B(HE) Aziz potential , which has proved high accuracy in liquid <sup>4</sup>He calculations , has modelled the atomic interactions. At the experimental equilibrium density $`\rho _0^{\mathrm{expt}}=0.273\sigma ^3`$ ($`\sigma =2.556`$ Å), we have started the calculation with $`b=1.15\sigma `$ and the backflow parameters optimized in Ref.  ($`\lambda _\mathrm{B}=0.14`$, $`r_\mathrm{B}=0.74\sigma `$, $`\omega _\mathrm{B}=0.54\sigma `$). With this initial set of parameters, the results obtained are clearly biased by the trial wave function. Even though the RN approach  corrects numerically the shortcomings of $`\psi (𝐑)`$ and, in some applications, allows for an exact estimation of the eigenvalue, this is not the case for liquid <sup>3</sup>He. In Fig. 1, the RN energies as a function of the released times are shown for the cases $`\lambda _\mathrm{B}=0`$ and $`\lambda _\mathrm{B}=0.14`$. As one can see, at small imaginary times the RN method reveals the presence of corrections to the FN energies but a common asymptotic regime is far beyond the scope of the available MC data.

The next step then was looking for the next order correction to backflow correlations, as well as for a possibly better set of backflow parameters. We have found that the ones we were using correspond to a local minimum of the FN energy, and that a narrower but deeper minimum exists with $`\lambda _\mathrm{B}=0.35`$ and $`r_\mathrm{B}`$ and $`\omega _\mathrm{B}`$ unchanged. The resulting energy versus released time is also plotted in Fig. 1. The relation of initial slopes, $`1:0.27:0.016`$ for $`\lambda _\mathrm{B}=0,0.14,0.35`$, provides information on the accuracy of $`\psi (𝐑)`$. In the optimal case, $`\lambda _\mathrm{B}=0.35`$, the slope is practically inexistent and the energy correction would be $`0.01`$ K if the asymptotic regime could be reached. In order to get additional evidence on the size of this correction, and as a closing checkmark of the reliability of our results, we have included corrections to the backflow trial wave function using the analytical method previously described. It can be shown that these new terms incorporate explicit three-body correlations in $`\phi _{\alpha _i}(\stackrel{~}{𝐫}_j)`$ of the form

$$\stackrel{~}{𝐫}_j^{\mathrm{BFT}}=\stackrel{~}{𝐫}_j^{\mathrm{BF}}+\lambda _{\mathrm{BT}}\underset{kj}{}\eta (r_{jk})(_j_k),$$
(5)   
with $`_i=_{li}\eta (r_{il})𝐫_{il}`$. We have carried out a FN-DMC calculation with this new trial wave function at $`\rho _0^{\mathrm{expt}}`$ and the result for the energy correction has been found $`<0.01`$ K. Both this analytical check and the numerical findings provided by the RN method point out the excellent description that backflow correlations make of the nodal surface in liquid <sup>3</sup>He.

The FN energies with $`\lambda _\mathrm{B}=0`$ (no backflow), $`\lambda _\mathrm{B}=0.14`$, and $`\lambda _\mathrm{B}=0.35`$ are reported in Table I, together with the corresponding kinetic energy obtained as the difference between the total energy and a pure estimation  of the potential energy. The comparison with the experimental energy , also contained in the table, shows the successive improvement of the FN-DMC result until an excellent agreement with $`\lambda _\mathrm{B}=0.35`$. Concerning the kinetic energy, a sizeable difference between theory and experiment  survives, a fact that has been generally attributed to long-range wings in the high-$`q`$ inelastic response that are difficult to incorporate effectively in the experimental analysis.

The FN-DMC calculation has been extended to a wide range of densities ranging from the spinodal point up to a maximum value $`\rho =0.403\sigma ^3`$, located near to the experimental freezing density $`\rho _\mathrm{f}^{\mathrm{expt}}=0.394\sigma ^3`$. We have used $`N=114`$ only at the highest density and below that $`N=66`$ have proved to be accurate enough. Among the three variational parameters entering in the backflow wave function (4) only $`\lambda _\mathrm{B}`$ shows a density dependence which is nearly linear in the range studied ($`\lambda _\mathrm{B}=0.42`$ at $`\rho =0.403\sigma ^3`$). The results are displayed in Fig. 2 in comparison with the experimental data of Ref. . The solid line in the same figure is a third-degree polynomial fit to our data with $`\chi ^2/\nu =1.2`$. According to this fit, the equilibrium density is $`\rho _0=0.274(1)\sigma ^3`$ and the energy at this density $`(E/N)_0=2.464(7)`$ K, in close agreement with experimental data.

The quality of the equation of state is even more stressed by looking at its derivatives. In Fig. 3, the behavior of the pressure and the sound velocity with the density is shown in comparison with experimental data from Refs. . The theoretical prediction for both quantities, derived form the polynomial fit to $`E/N(\rho )`$ (Fig. 1), shows again an excellent agreement with the experimental data from equilibrium up to freezing. The sound velocity that at $`\rho _0^{\mathrm{expt}}`$ is $`c=182.2(6)`$ m$`/`$sec, in close agreement with the experimental value $`c^{\mathrm{expt}}=182.9`$ m$`/`$sec , goes down to zero at the spinodal point. The location of this point has been previously obtained both from extrapolation of experimental data at positive pressures  and from density-functional theories . The present microscopic calculation allows for an accurate calculation, free from extrapolation uncertainties, that locates the spinodal point at a density $`\rho _\mathrm{s}=0.202(2)\sigma ^3`$ corresponding to a negative pressure $`P_\mathrm{s}=3.09(20)`$ atm, much closer to the equilibrium than in liquid <sup>4</sup>He where $`P_\mathrm{s}=9.30(15)`$ atm .

In conclusion, we have analyzed the possibilities of the diffusion Monte Carlo method in the study of Fermi systems. The FN-DMC method, combined with the RN technique to evaluate the effect of the nodal surface model used in the trial function, and with a systematic method to analytically improve the nodal surface, constitutes a closed loop able to improve the quality of the antisymmetric wave function up to a required precision. The fact that this general approach has allowed to deal with a strongly interacting Fermi liquid suggests that it could be also useful for tackling other Fermi systems. We have applied this methodology to the study of liquid <sup>3</sup>He. The effect of corrections beyond the backflow terms has been evaluated and found to be less than 0.01 K. The equation of state so obtained presents an accuracy comparable to the result obtained in bosonic liquid <sup>4</sup>He. The precision proved by the method allows for posterior studies as the characterization of spin-polarized liquid <sup>3</sup>He. Preliminary calculations for the fully-polarized phase at $`\rho _0^{\mathrm{expt}}`$ indicate a less binded system with an energy $`E/N=2.22(4)`$ K.

This research has been partially supported by DGES (Spain) Grant N<sup>0</sup> PB96-0170-C03-02. We also acknowledge the supercomputer facilities provided by the CEPBA.