text:
# Backflow Correlations for the Electron Gas and Metallic Hydrogen

## I Introduction

This paper concerns the form of the ground state wavefunction of metallic hydrogen at high enough density so that all the hydrogen molecules are dissociated and the electrons are delocalized. Neglecting possible quantum effects on the protonic motion, the many-body wavefunction can be regarded as the ground state of an electron gas under the influence of an external potential due to the actual positions of the protons. Quantum Monte Carlo (QMC) techniques are currently one of the most powerful methods to calculate accurately the properties of such a many-body quantum systemRMPQMC . However, since ground state QMC is based on trial wavefunctions, QMC typically demands compact and accurate descriptions of the ground state wavefunction. In this paper we review different approaches to obtain and improve trial wavefunctions, compare the qualities of the resulting many-body wavefunctions with previous QMC calculations for the electron gas and metallic crystal hydrogen, and present first results using these wavefunctions for liquid metallic hydrogen.

Most of the work within QMC has been done using a pair product (PP) (or Slater-Jastrow) wavefunction: a Slater determinant of single electron spin orbits times a product of pair electron (Jastrow) factors. Notwithstanding certain deficiencies such as a lack of direct spin coupling, this wavefunction has proven to be quite accurate, in particular within fixed-node Diffusion Monte CarloRMPQMC  (DMC). The first calculation on many-body hydrogendmc81  used an even simpler form of this wavefunction; the single electron orbits were taken to be free electron plane waves. We refer to this as the SJ-PW trial function. Later, Natoli natoli93 ; natoli95a  found that determinants using these orbitals are inaccurate by 0.05eV/atom within the fixed-node DMC calculations at the density corresponding to the transition between molecular and metallic hydrogen ($`r_s=1.31`$). Hence, more accurate orbitals, computed from either density functional (LDA) or Hartree-Fock (HF) calculations, are required. Because these orbitals are calculated assuming fixed ionic positions, inclusion of ionic motions, such as those from the zero point motion of the ions in the crystal, is difficult.

Recently, there have been new attemptsdewing01 ; dmc03  to calculate properties of disordered systems such as liquid hydrogen within QMC. In the Coupled Electron Ion Monte Carlo (CEIMC) methoddewing01  the protons are moved based on the results of a QMC calculation of the electronic energy. This approach requires accurate trial functions that can be obtained quickly as the ionic positions are changed; methods involving the solution of mean field equations such as LDA and HF, or even optimizing a parameterized trial function, can greatly slow down the overall performance of the CEIMC simulationdewing01 . Further, combining the orbitals obtained from LDA or HF with a pair correlation (Jastrow) factor to improve the accuracy is not straightforward; substantial modification of the orbitals might be necessary requiring a reoptimization of the orbitals and correlation factor, in principle, at each ionic positionfahy . This optimization step creates a bottleneck to coupling the QMC calculations with the ionic Monte Carlo.

One could consider obtaining the trial wavefunction from other variational approaches like Fermi-hypernetted (FHNC) chain or correlated basis functions (CBF) methodskrotscheck  which would not have the problems of optimization. However, in these approaches based on explicit integration, one is in general limited in the form of the trial function by the ease performing the integration, and these are typically much more time-consuming than LDA calculations.

One of the biggest advantages of the QMC approach is that one can use an arbitrary wavefunction without changing the algorithm in an essential way. Fast algorithms will result if one can find concise and accurate forms. In this paper, instead of using one-body orbitals from mean field theory or integral equations, we propose to use trial functions which depend explicitly and continuously on the ionic variables. Such wavefunctions do not have to be reoptimized for movements of the ions, are easy to implement, and accurate for disordered systems. Calculation of ionic forces is simplified since the derivative of the trial function with respect to ionic configurations is a straightforward application of the chain rule. These trial functions are a generalization of the backflow three-body wavefunctions used very successfully in highly correlated homogeneous quantum liquids: liquid 3He and the electron gas. There, backflow trial functions show much improvement over the pair product getting approximately 75% of the energy missing at the PP level and even more when done with the fixed-node method.

Backflow wavefunctions were developed by Feynman and Cohen feynman  for a single 3He impurity in liquid 4He when it was recognized that without backflow, the mass of the impurity was equal to the bare mass. Pandharipande and Itoh panda73  showed that the backflow arises from the momentum dependence of the correlation between the impurity and the liquid. The backflow wavefunction was then extended to bulk liquid 3He schmidt79 ; manousakis  using an integral equation method to evaluate expectation values. The first use of backflow in QMC was by Lee et al.lee81  and others schmidt81 ; panoff89  with calculations on liquid 3He. Moroni et al.moroni95  further optimized the trial function within liquid 3He. Kwon et al.kwon93 ; kwon98  used backflow functions for the electron gas in both 2D and 3D, obtaining significantly lower energies and improved excitation energies. Vitiello et al.vitiello97  discuss an equivalence of backflow and spin-dependent correlations, an aspect we will not further consider in this paper.

Using different approaches, we generalize the backflow three-body wavefunction to a two component system of electrons and protons and derive approximate expressions for the correlated trial function. We first present an argument based on the generalized Feynman-Kacs formula which shows that backflow is the next order improvement beyond the pair product (PP) wave function. Using perturbation theory, we then discuss general features of the backflow functions and obtain explicit expressions for the homogeneous electron gas and for the electron-proton plasma. A similar analysis using the Bohm-Pines method has been recently performed by Gaudoin et al.gaudoin , however, without going beyond the Slater-Jastrow wavefunction. Studying the problem of a single electron in the potential generated by a simple cubic lattice of protons, we show that the exact one electron wavefunction can be approximately rewritten by a backflow function. Finally, we optimize numerically simple functional forms for the backflow functions in the full many-body problem by variational Monte Carlo. We compare the quality of the wavefunctions stemming from these different approaches for the electron gas and for liquid and crystal hydrogen at the level of variational and diffusion Monte Carlo.

In the following we consider the non-relativistic Hamiltonian of $`N`$ protons and $`N`$ electrons:

$$\widehat{H}=\underset{i}{}\lambda _i_i^2+\underset{i<j}{}\frac{e_ie_j}{r_{ij}}$$
(1)
where $`\lambda _i=\mathrm{}^2/(2m_i)`$, $`i=1,\mathrm{},2N`$ and $`m_i`$ and $`e_i`$ are the electron or proton mass and charge. The Fermi wavevector is $`k_F`$. Numerical results are given in atomic units where $`\lambda _e=1/2`$ and $`\lambda _p=0`$ for classical protons, $`|e_i|=m_i=1`$. The electron density $`n=N/V`$ is quoted in terms $`r_s=a/a_0`$, where $`a=(4\pi n/3)^{1/3}`$ and $`a_0=\mathrm{}^2/m_ee^2`$ is the Bohr radius. Energies of the QMC calculations are given in Rydbergs per electron.

## II The Feynman-Kacs approach to improving the wavefunction

The Feynman-Kacs formula expresses the exact wavefunction in terms of average over Brownian paths. We now review how it can be generalized to random walks with “drift”.

We define the “importance-sampled” Green’s function as $`\widehat{G}_t=\psi \mathrm{exp}(t\widehat{H})\psi ^1`$ in operator notation where $`\psi `$ is an unsymmetrical trial function. $`\widehat{G}_t`$ acting on a function has the effect of enhancing the component of lower energy states. Then the lowest energy (exact) fermi wavefunction $`\varphi _F(R)`$ is given by:

$$\varphi _F(R)\widehat{𝒜}\psi (R)\underset{t\mathrm{}}{lim}𝑑R^{}R^{}|\widehat{G}_t|R$$
(2)
assuming only that the trial function has a non-zero overlap with $`\varphi _F`$ and that $`\varphi _F`$ is non-degenerate. Here $`\widehat{𝒜}`$ is a projection operator for fermion symmetry defined as

$$\widehat{𝒜}f(R)=\frac{1}{N!}\underset{P}{}(1)^Pf(PR)$$
(3)
and $`R=\{𝐫_1,𝐫_2,\mathrm{}\}`$ is a point in configuration space. The electron spin is treated by restricting the permutation in Eq. (3) to be exclusively within spin up or spin down electrons.

Following the derivation in Diffusion Monte Carloreynolds , the Green’s function can be split into diffusion, drift and branching processes. To show this, the master equation for the Green’s function is written:

$$\frac{d\widehat{G}}{dt}=\psi \widehat{H}\psi ^1\widehat{G}=[\underset{i}{}\lambda _i_i(_i+2_i\mathrm{ln}\psi )+E(R)]\widehat{G}.$$
(4)
The local energy $`E`$, defined as $`E(R)=\psi ^1\widehat{H}\psi `$, is the residual error of the trial function, and becomes a constant function as $`\psi `$ approaches an exact eigenvalue. Trotter’s formula applies to the above master equation, allowing us to split up the evolution into the first two terms describing a stochastic process, and the final term which is a branching or “weighting” process. Thus we have the generalized Feynman-Kacs formula :

$$\varphi _F(R)𝒜\psi (R)\mathrm{exp}\left[_0^{\mathrm{}}𝑑tE(R(t))\right]$$
(5)
where the brackets imply averaging over all drifting random walks $`R(t)`$ beginning at a point $`R`$. The above relation is exact for any real trial function. For trial functions having an imaginary component of $`\mathrm{ln}\psi `$, the formalism goes through, however, the Green’s function is no longer real and positive and therefore cannot be treated as a probability. Other methods are more appropriate. For the moment we will ignore this case.

To make further analytical progress, we take the average into the exponent. For any stochastic process, one can write the average of the exponent as the exponential of the cumulant expansion, the first two terms of which are:

$$\varphi _F(R)𝒜\psi (R)\mathrm{exp}(E+(1/2)\delta E^2\mathrm{}].$$
(6)
The double brackets are defined as $`E=_0^{\mathrm{}}𝑑tE(R(t))`$ with walks $`R(t)`$ generated from the drift and diffusion starting at a point $`R`$. We truncate the cumulant expansion after the first term. We then have an approximate method of improving the trial function.

$$\psi ^{(n+1)}=\psi ^{(n)}e^{E^{(n)}_n}$$
(7)
with the subscript indicating that the drift is given by $`\mathrm{ln}\psi ^{(n)}`$. If we split the log of the trial function into its real and imaginary parts $`\psi ^{(n)}=\mathrm{exp}(U^{(n)}+iS^{(n)})`$ with $`U`$ and $`S`$ real, we are led to the following equations for a single iteration:

$`U^{(n+1)}`$ $`=`$ $`U^{(n)}`$
$`+`$ $`V`$ $`+{\displaystyle \underset{i}{}}\lambda _i[_i^2U^{(n)}(_iU^{(n)})^2+(_iS^{(n)})^2_n`$
$$S^{(n+1)}=S^{(n)}+\underset{i}{}\lambda _i[_i^2S^{(n)}2_iU^{(n)}_iS^{(n)}]_n.$$
(9)
Here $`V(R)`$ is the total potential energy.

Specializing to the case of a fermi liquid, we take as an initial wavefunction $`U^{(0)}=0`$ and $`S^{(0)}=_i𝐤_i𝐫_i`$, i.e. singly occupied free particle states. (The usual spin functions are assumed but not explicitly written.) Note that this function is an unsymmetrical trial function, with a non-zero overlap with a fermion state as long as all the $`𝐤_i`$’s are distinct. When the wavefunction is antisymmetrized, one gets a determinant of plane waves. However, the antisymmetrization will be done only once, after the trial function has gone through several iterations of Eq. (7). This way simplify the procedure, since the local energy of the unsymmetric trial function is much simpler than that of an antisymmetric trial function. Note that in Eq. (5) both the antisymmetrization and the averaging are linear operators and so can be interchanged.

After the first iteration, the wavefunction will have the form:

$$U^{(1)}=V(R)_0U(R)$$
(10)
$$S^{(1)}=S^{(0)}.$$
(11)
In the above equation and the following discussion we drop, without mention, constant normalization terms. If $`V(R)=_{i<j}v(r_{ij})`$ is a pair potential, with a Fourier transform $`v_k`$, the averaging can be carried out analytically with a result that $`U^{(1)}`$ will also be a pair potential and will have a Fourier transform given by $`v_k/(\lambda k^2)`$ where $`\lambda =\lambda _i+\lambda _j`$. For a Coulomb potential the real-space correlation (Jastrow) function will then have the form: $`u_C(r)=e^2r/2\lambda `$.

Hence, the form of the first order wavefunction is of the SJ-PW or Slater-Jastrow form, with free particle orbitals. The pair term will be denoted by $`U`$ with $`U=_{i<j}u(r_{ij})`$. Typically the form of $`u`$ is derived from a variational principle, chosen so such that either the total energy or variance is minimized. This will, of course, give a lower energy than the cumulant form derived above. The above derivation does give the correct cusp condition (the limit of $`u`$ at large $`k`$ or small $`r`$). However, it does not give the long wavelength limit correctly because of the neglect of the higher cumulants. Gaskellgaskell  proposed an analytic form based on the Random Phase Approximation (RPA) without any parameters. It was founddmc78  for the homogeneous electron gas that the RPA form does, as well as, or better than simple assumed forms with parameters. Figure 1 shows a comparison of these correlation functions.

Note that the cumulant approximation will not exist if the Fourier transform of the potential does not exist. Two examples of such potentials are the hard sphere and Lennard-Jones interactions. However, for the short range part of a soft potential which does have a Fourier transform such as the Yukawa potential, the cumulant approximation works quite well (see remarks concerning the situation at finite temperature in Ref. rmpi ).

We now perform the next iteration of this procedure. To minimize the fluctuations in the local energy so that the cumulant approximation will be more accurate, we assume that first order wavefunction has been optimized but it still has a pair product form. Using Eq. (9), neglecting constants and combining pair terms together, we get in second order a function of the form:

$$U^{(2)}=\stackrel{~}{U}^{(1)}\underset{i}{}\lambda _i(_iU^{(1)})^2_1$$
(12)
and

$$S^{(2)}=\underset{i}{}𝐤_i(𝐫_i2\lambda _i_iU)_1.$$
(13)
Here $`\stackrel{~}{U}^{(1)}`$ includes additional pair terms.

At second order, we cannot perform the averaging analytically, since it involves drift under the influence of the first order wavefunction: $`U^{(1)}`$. We make the assumption that the averaging will not change the functional form of the quantity being averaged but only smoothes out the individual functions. That is, our ansatz for the iterated wavefunction is:

$$U^{(2)}=U(R)\underset{i}{}(_iW)^2$$
(14)
and

$$S^{(2)}=\underset{i}{}𝐤_i(𝐫_i_iY)$$
(15)
where $`U`$, $`W`$ and $`Y`$ are three different pair “potentials” to be optimized. In the following, we have adopted the convention that pair functions have the same sign as $`v_{ij}(r)`$, so that, for example, a repulsive $`v`$ leads to a repulsive $`w`$ and $`y`$.

The two new functions appearing at second order are the backflow function $`Y`$ and the three-body or polarization term $`W`$. The backflow potential is

$$Y=\underset{i<j}{}y_{ij}(r_{ij})$$
(16)
where $`y(r)`$ is a spherically symmetric function and the sum extends over all pairs of particles, including both electrons and protons. The backflow displacement is defined as the gradient of the backflow potential with respect to a particle coordinate:

$$\mathrm{\Delta }𝐫_i=_iY=\underset{ji}{}\eta (r_{ij})(𝐫_i𝐫_j)$$
(17)
where

$$\eta (r)=\frac{1}{r}\frac{dy(r)}{dr}$$
(18)
corresponds to the definition in previous work for homogeneous systems kwon93 ; kwon98 .

With this ansatz, the antisymmetized trial function is a determinant composed of “quasi-particle” coordinates:

$$\psi _F^{(2)}=det\left[\mathrm{exp}i𝐤_j(𝐫_i+\mathrm{\Delta }𝐫_i)\right]e^{U(R)+(W)^2}.$$
(19)
Recall that in the fixed-node or fixed-phase diffusion Monte Carlo method, one obtains the exact energy subject to the imposed constraint RMPQMC ; ortiz93 . The assumed node or phase limits the ultimate accuracy for fermion systems. Since the correction to the real part, the three-body term, is already symmetric, it is the backflow which is responsible for the change of node or phase of the trial function and is, in that sense, more important than the Jastrow and polarization part.

In the above derivation we have neglected any effects of a complex drift velocity. However, as already shown in Ortiz and Ceperleyortiz95 , a complex drift velocity does not affect the corrections to the wavefunction to the order we have considered; the Eqs (9) are valid to improve the wavefunction.

Now we consider the long range properties of the pair functions appearing. In periodic boundaries (or “super-cells”) we need to perform Ewald summations of the functions $`V,U,W,Y`$. This is most convenient in Fourier space. We define the Fourier transform of a radial function as:

$$\stackrel{~}{y}_k=𝑑𝐫^3e^{i\mathrm{𝐤𝐫}}y(r).$$
(20)
Using the Poisson sum formula, the “potential” of the $`i^{th}`$ particle in periodic boundary conditions is:

$$y_i=\frac{1}{V}\underset{𝐤,j}{}\stackrel{~}{y}_ke^{i𝐤(𝐫_i𝐫_j)}$$
(21)
where $`V`$ is the volume of the supercell. For example, to find the backflow displacement, Eq.(17), we simply take the gradient of the pair function:

$$\mathrm{\Delta }𝐫_i=\frac{1}{V}\underset{𝐤,j}{}i𝐤\stackrel{~}{y}_ke^{i𝐤(𝐫_i𝐫_j)}$$
(22)
where $`𝐤`$ ranges over the reciprocal lattice vectors of the supercell.

The three-body potential, $`W`$ is defined analogously in terms of a pair polarization $`w(r)`$. This function is related to that used in previous QMC workpanoff89 ; moroni95 ; kwon93 ; kwon98  by:

$$\sqrt{|\lambda _T|}\xi (r)=\frac{1}{r}\frac{dw(r)}{dr}.$$
(23)
The overall sign of $`w`$ is not important because only its square appears in the trial function, but the relative sign of the electron-electron to the electron-proton interaction is significant.

One of the simple ways of deriving conditions on the backflow function is to look at the action of the Hamiltonian on the wavefunction, the local energy, and to minimize the fluctuations of the local energy. Here we focus on the imaginary part of the local energy and consider a single electron with phase $`S=𝐪(𝐫Y)`$. Setting to zero the imaginary part of the local energy we obtain:

$$^2y(r)+2u(r)2u(r)y(r)=0.$$
(24)
Neglecting the last term, since it is higher order in the interaction, we obtain: $`^2y(r)=2u(r).`$ This has a solution in Fourier space: $`y_k=2u_k/k^2`$. (Because we want a solution which is smoother than $`u(r)`$ at $`r=0`$, we neglect a term proportional to $`r^1`$.) We get the same smoothing ($`k^2`$) that we observed at first order for the pair function. Shown in figure 2 are the $`u(r)`$ and $`\eta (r)`$ function coming from this approach. Note that this approach is based on a single electron description and therefore does not correctly describe the long wavelength (large r) behavior where the collective motion dominates.

To obtain a simple form for the three-body potential, we note that the averages used in the definition of $`Y`$ are similar to those for $`W`$, see Eq.(12-15). Hence an estimate of the polarization potential is

$$W=\sqrt{\lambda }k_FY$$
(25)
where we have approximated $`(_iU)^2(_iU)^2/\tau `$, averaged over a “typical” time $`\tau (\lambda k_F^2)^1`$. This relates the three-body contribution to the backflow potential.

The GFK approach is good for suggesting corrections, but there are serious problems in using it to find a good backflow function since the averaging is difficult to carry out, the linear cumulant approximation may be inadequate, and the long-time effects of the imaginary drift are being ignored. If one cannot analytically perform the averaging, one does not know what time to multiply the local energy by to get a wavefunction, nor the relative corrections at large versus small distances. We now discuss several other approaches which allow us to directly evaluate the Jastrow, threebody and backflow functions and give more insight into their form.

## III Perturbation Theory/Analytic methods

In this section we follow another approach to obtaining improved estimates of the many-body wavefunction. Many-body perturbation theory is a well studied approach to understanding the effects of weak correlation. Encouraged by the use of the RPA gaskell  which gave an excellent analytic two-body correlation function, we will extend this wavefunction by perturbative expressions for the Jastrow, backflow and threebody potentials for the electron gas and for metallic hydrogen. Rather than performing a systematic low or high density expansion to derive analytical expressions for the variational wavefunction of the electron gas or metallic hydrogen, we concentrate on improving this correlation factor. The collective coordinate formulation of Bohm and Pines BP ; gaudoin  allows us to use Slater-Jastrow wavefunctions as zeroth-order starting point. We obtain improved potentials for the homogeneous electron gas and metallic hydrogen, which compare very well with numerically optimized forms.

Even if perturbation theory assumes a weak coupling (or high density) expansion, we expect the derived properties to be qualitatively valid as long as the corresponding perturbation expansion remains regular, e. g. until there is a phase transition to an insulating phase.

### III.1 Single particle perturbation theory

Consider a single electron interacting with an arbitrary external potential $`v(𝐫)`$ with Fourier transform $`\stackrel{~}{v}(𝐤)`$. To avoid the problems arising from the long-range behavior of the Coulomb interaction, we restrict the analysis to a potential with a Fourier transform which is finite at the origin, $`|\stackrel{~}{v}(0)|<\mathrm{}`$, e.g. a screened Coulomb potential. We use the continuum notation in this section ($`\frac{1}{V}_k\frac{d^3k}{(2\pi )^3}`$). The solution of the Schrödinger equation $`\varphi _k(𝐫)`$ of a particle with wavevector $`𝐤`$,

$$\varphi _k(𝐫)=\frac{1}{(2\pi )^3}d^3𝐩c_k(𝐩)\mathrm{exp}[i𝐩𝐫]$$
(26)
can be written as

$$c_k(𝐩)=(2\pi )^3\delta (𝐤𝐩)+\frac{4\pi f(𝐤,𝐩)}{k^2p^2+i\delta }$$
(27)
where the off-shell scattering amplitudes, $`f(𝐤,𝐩)`$, are given by the integral equation

$$4\pi \lambda f(𝐤,𝐩)=\frac{d^3k^{}}{(2\pi )^3}\stackrel{~}{v}(𝐩𝐤^{})c_k(𝐤^{}).$$
(28)
Using the Born approximation we can write down the wavefunction to first order in $`\stackrel{~}{v}`$

$`\varphi _k(𝐫)`$ $``$ $`\varphi _k^{(0)}(𝐫)+\varphi _k^{(1)}(𝐫)`$ (29)
$`=`$ $`e^{i𝐤𝐫}\left(1{\displaystyle \frac{1}{\lambda }}{\displaystyle \frac{d^3p}{(2\pi )^3}e^{i𝐩𝐫}\frac{\stackrel{~}{v}(p)}{𝐩(𝐩+2𝐤)}}\right)`$
If we expand the solution around $`𝐤=0`$ and assume the change in the wavefunction is small, we can write it in the pair-product and backflow form, Eq.(19). We obtain for the pair potential

$$u(𝐫)\frac{d^3p}{(2\pi )^3}e^{i𝐩𝐫}\frac{\stackrel{~}{v}(p)}{\lambda p^2},$$
(30)
and for the backflow potential

$$y(𝐫)2\frac{d^3p}{(2\pi )^3}e^{i𝐩𝐫}\frac{\stackrel{~}{v}(p)}{\lambda p^4}.$$
(31)
Note that the small $`𝐩`$ part of the integral is usually cut off by the finite size of the box. In addition, it is $`\eta (r)`$, the derivative of $`y(r)`$ (see Eq. (18)) which enters in the trial function.

Although the first order approximation is only reliable in the case of a weak potential, it becomes correct in the high momentum region and hence, gives the correct cusp conditions. The derived form is identical to that obtained from the Feynman-Kacs formula in the previous section. For an arbitrary weak potential, we further get the long range behavior, $`u\stackrel{~}{v}(0)/r`$ and $`\eta 1/r^3`$ for $`r\mathrm{}`$, provided the potential has a finite range ($`\stackrel{~}{v}(q)\stackrel{~}{v}(0)q^2`$ for $`q0`$) and there is no other singularity in the integrand.

To find an approximate form for the three-body function $`W(𝐫)`$ we must go to higher order in the interaction, but only at $`𝐤=0`$. Using Eqs. (27-28), we can write down the second order corrections in $`\stackrel{~}{v}`$ to the wavefunction, $`\varphi _k^{(2)}`$ at $`k=0`$:

$$\varphi _{k=0}^{(2)}(𝐫)=\frac{1}{\lambda ^2}\frac{d^3q}{(2\pi )^3}\frac{d^3p}{(2\pi )^3}\frac{\stackrel{~}{v}(q)e^{i𝐪𝐫}}{(𝐪+𝐩)^2}\frac{\stackrel{~}{v}(p)e^{i𝐩𝐫}}{p^2}$$

$$\frac{d^3q}{(2\pi )^3}\frac{d^3p}{(2\pi )^3}\frac{\stackrel{~}{v}(q)e^{i𝐪𝐫}}{\lambda q^2}\left(1\frac{2𝐪𝐩}{q^2}\right)\frac{\stackrel{~}{v}(p)e^{i𝐩𝐫}}{\lambda p^2}.$$
(32)
This is almost in the form of the three-body correlation obtained with the FK approach: $`(w)^2`$. Note, however, that Eq. (33) is unsymmetrical in $`q`$ and $`p`$ so that in r-space it will be written as: $`(w_u)(w_y)`$ with $`w_u(𝐫)u(r)`$ and $`w_y(𝐫)y(𝐫)`$. Therefore the polarization term is not a square but a product of the gradients of two different functions. (In second order one will also find a contribution $`[u(r)]^2`$ to the pair term.)

The perturbative expressions (30) and (31) are based on the Born approximation for scattering between free states. However, an attractive potential as the electron-proton (effective) interaction might also lead to bound states. To include the effects of a possible bound state we can use the non-perturbative expression (28) for the scattering amplitudes: given an approximate expression for the bound state wavefunction of energy $`ϵ_k=\lambda k^2`$, we can calculate the scattering amplitudes and obtain corrections from the bound state to the pair and backflow potential in the same way as shown above for the scattering states within the Born approximation. In a similar way one should proceed to obtain approximations for the pair and backflow potentials for systems where the interatomic potentials cannot be treated within the Born-approximation, for example potentials dominated by a hard core.

Of course, in the case of a single electron in an external potential we can solve the Schrödinger equation by other means and obtain the “best” pair and backflow potentials from the exact (numerical) solution. This is done below for a perfect crystal using a band structure calculation. However, the simple perturbative approach above provides an easy way to get some intuition for the pair and backflow potential, and is already good enough to determine their asymptotic properties. These properties are expected to hold in the many-body case: the short-range properties are typically determined by two-body collisions and the influence of the remaining particles on the long-range properties are usually well described by an effective single particle potential. Many-body perturbation theory, which we discuss next, leads to similar expressions.

### III.2 Many-body perturbation theory

We now make an expansion of the exact $`N`$ particle wavefunction $`|\varphi `$ of the interacting system around the non-interacting (ground) state $`|\varphi _0`$; the ground state without both the electron-electron and electron-proton interaction. Let $`a_k`$ ($`a_k^{}`$) be the annihilation (creation) operator for an electron of wavevector $`k`$. Expanding in particle-hole excitations, we have:

$$|\varphi \left(1+\underset{q,k_1,k_2}{}\alpha _{k_1,k_2,q}a_{k_2q}^{}a_{k_2}a_{k_1+q}^{}a_{k_1}+\mathrm{}\right)|\varphi _0.$$
(33)
The problem is reduced to determining the coefficients $`\alpha _{k_1,k_2,q}`$. Just as in the single particle case, a further expansions of $`\alpha _{k_1,k_2,q}`$ around $`k_1=0`$ or $`k_2=0`$ together with an exponentiation brings the wavefunction into the desired functional form, thereby determining the pair and backflow potentials. To avoid over counting, we assume that the summation in Eq.(33) goes only over distinct states so that it is sufficient to antisymmetrize the wavefunction at the very end, once we have calculated the perturbative corrections. We have limited the expansion in Eq.(33) to the leading order corrections, particle-hole excitations; the generalization to include higher order excitations is straightforward, but not necessary to calculate the pair and backflow terms in the wavefunction.

In order to determine the coefficients, we write $`|\mathrm{\Phi }_{k1,k2,q}=a_{k_2q}^{}a_{k_2}a_{k_1+q}^{}a_{k_1}|\varphi _0`$ and multiply these states by a constant phase

$$\alpha _{k1,k2,q}=\mathrm{\Phi }_{k1,k2,q}|\varphi \frac{\varphi |\varphi _0\mathrm{\Phi }_{k_1,k_2,q}^{ee}|\varphi }{\varphi |\varphi _0\varphi _0|\varphi }$$
(34)
Note that with this phase factor, the right hand side of Eq. (34) is given by the expectation value of an operator over the true ground state and the coefficients can therefore be identified as a $`N`$-particle Green’s function Greens . By considering only particle-hole excitations in Eq.(33), the $`N`$-particle Green’s function reduces to a connected two-particle Green’s function and the lowest order modifications to the ideal gas ground state of the homogeneous electron gas are therefore related to the two particle Green’s function $`G^{(2)}`$ at equal times Greens  or equivalently

$$\alpha _{k1,k2,q}\frac{\varphi |a_{k_2q}^{}a_{k_2}a_{k_1+q}^{}a_{k_1}|\varphi }{\varphi |\varphi }$$
(35)
Summing up particle-hole bubble diagrams (corresponding to the RPA approximation) results in an effective interaction, $`\stackrel{~}{v}_{RPA}(p,\omega )`$,

$$\stackrel{~}{v}_{RPA}(p,\omega )=\frac{\stackrel{~}{v}(p)}{ϵ(k,\omega )},ϵ(k,\omega )=1\stackrel{~}{v}(p)D(p,\omega )$$
(36)
where $`D(p,\omega )`$ is the Lindhard function. Perturbation theory can now be arranged to be regular Gell-Mann . We note that Eq. (36) already contains the correct short- and long-range limit of the effective interaction.

Neglecting for the moment any contributions from plasmon excitations coming from the poles where $`ϵ(k_p,\omega _p(k_p))=0`$, we get

$`\alpha _{k_1,k_2,q}`$ $`=`$ $`(1n_{k_1q})n_{k_1}(1n_{k_2+q})n_{k_2}`$
$`\times `$ $`{\displaystyle \frac{\stackrel{~}{v}_{RPA}(q,\epsilon _{k_1}\epsilon _{k_1q})+\stackrel{~}{v}_{RPA}(q,\epsilon _{k_2}\epsilon _{k_2+q})}{2(\epsilon _{k_1}+\epsilon _{k_2}\epsilon _{k_1q}\epsilon _{k_2+q})}}`$
where $`n_k`$ are the occupation numbers of state $`k`$ in lowest order. Expanding around $`k_1=k_2=0`$, we get the Jastrow and the backflow potential. Including the plasmon excitations will give an important long-range contribution. However, in the simplest approximation, this contribution describes only the long wavelength limit correctly, and destroys the correct short distance behavior. We will circumvent this problem in the next section using the formalism of collective coordinates.

As already shown in the previous section, we expect a more general form for the three-body potential,

$$\varphi det\left[\mathrm{exp}i𝐤_j(𝐫_i+\mathrm{\Delta }𝐫_i)\right]e^{U(R)+W}$$
(38)
with W=

$$\underset{j}{}(_jW_u)(_jW_y)2\underset{i<j}{}[_jw_u(r_{ij})][_jw_y(r_{ij})],$$
(39)
where

$$w_u(r)u(r)w_y(r)y(r).$$
(40)
For the interactions of the electrons with static protons, we can use the static dielectric function $`ϵ(k,0)`$ to obtain the effective electron-proton interaction, and use directly the results of the single particle perturbation theory of the previous section with this screened potential.

The disadvantage of perturbation theory is that one gets correct behavior at long and short distances, but it does not provide an unique way to interpolate between these limits. In Table 1 we summarize the asymptotic properties of the pair and backflow potentials for the 3D electron gas.

### III.3 The Bohm-Pines collective coordinate approach

Instead of replacing the established form for the Jastrow part proposed by Gaskell gaskell  by the direct use of Eq. (III.2), we prefer to improve the RPA form of Gaskell by extending it using perturbative formulas. This is most easily done within the framework of the collective coordinate description of Bohm and Pines using additional field variables BP . In this approach, the original Hamiltonian of electrons interacting with each other and with static protons is extended by an additional boson field with generalized momentum variables $`\mathrm{\Pi }_k`$ coupling to the electron and proton density fluctuations

$`H`$ $`=`$ $`{\displaystyle \underset{i}{}}\lambda p_i^2+{\displaystyle \frac{1}{2V}}{\displaystyle \underset{k}{}}\stackrel{~}{v}_k\left(\rho _k^e\rho _k^eN\right)`$ (41)
$`{\displaystyle \frac{1}{V}}{\displaystyle \underset{k}{}}\stackrel{~}{v}_k\rho _k^e\rho _k^p`$
$`+{\displaystyle \frac{1}{V}}{\displaystyle \underset{k}{}}\left({\displaystyle \frac{\mathrm{\Pi }_k^{}\mathrm{\Pi }_k}{2}}+M_k\mathrm{\Pi }_k^{}\rho _k^e+P_k\mathrm{\Pi }_k^{}\rho _k^p\right)`$
where $`\rho _k^e`$ ($`\rho _k^p`$) is the Fourier transform of the electron (proton) density, $`\rho _k=_ie^{i𝐤𝐫_i}`$, and $`M_k`$ and $`P_k`$ are variational parameters. By imposing the extra conditions $`\mathrm{\Pi }_k\mathrm{\Psi }=0`$ on the wavefunction, the ground state wavefunction of the new extended Hamiltonian will be identical to the original one. For a detailed description of this approach we refer to the original literature BP ; we will only describe the main steps.

Carrying out the following canonical transformation

$$\varphi _{old}=\mathrm{exp}[iS/\mathrm{}]\varphi _{new},S=\frac{1}{V}\underset{k}{}\left(M_k\rho _k^e+P_k\rho _k^p\right)Q_k,$$
(42)
where $`Q_k`$ represents the field coordinate conjugate to $`\mathrm{\Pi }_k`$, we obtain an equivalent Hamiltonian

$`H`$ $`=`$ $`{\displaystyle \underset{i}{}}\lambda _ep_i^2+{\displaystyle \frac{1}{V}}{\displaystyle \underset{k}{}}\left({\displaystyle \frac{\mathrm{\Pi }_k^{}\mathrm{\Pi }_k}{2}}+\lambda _enk^2M_k^2Q_k^{}Q_k\right)`$ (43)
$`+H_{sr}^{ee}+H_{int}+H_{rw}+H_{sr}^{ep}`$
where

$`H_{sr}^{ee}`$ $`=`$ $`{\displaystyle \frac{1}{2V}}{\displaystyle \underset{k}{}}\left(\stackrel{~}{v}_kM_k^2\right)\left(\rho _k^e\rho _k^eN\right)`$ (44)
$`H_{int}`$ $`=`$ $`i{\displaystyle \frac{1}{V}}{\displaystyle \underset{k,j}{}}\left({\displaystyle \frac{𝐤𝐩_j}{m}}+{\displaystyle \frac{\mathrm{}k^2}{2m}}\right)M_kQ_ke^{i𝐤𝐫_j}`$ (45)
$`H_{rw}`$ $`=`$ $`{\displaystyle \frac{\lambda }{V^2}}{\displaystyle \underset{kk^{},j}{}}Q_k^{}Q_k^{}M_kM_k^{}𝐤𝐤^{}e^{i(𝐤𝐤^{})𝐫_j}`$ (46)
$`H_{sr}^{ep}`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{k}{}}\left(\stackrel{~}{v}_kM_kP_k\right)\rho _k^e\rho _k^p`$ (47)
Now the ground state of the additional field in the zeroth order Hamiltonian, Eq.(43), is simply given by harmonic oscillator ground states of frequencies $`V_k=(nk^2M_k^2/m)^{1/2}`$,

$$\varphi _{new}^0=det[\mathrm{exp}[i𝐤_i𝐫_j]]\mathrm{exp}\left[\frac{1}{V}\underset{k}{}\frac{\mathrm{\Pi }_k^{}\mathrm{\Pi }_k}{2\mathrm{}V_k}\right]$$
(48)
Transforming back and applying the subsidiary conditions replaces the field operator $`\mathrm{\Pi }_k`$ by $`M_k\rho _k^e+P_k\rho _k^p`$ and the zeroth order wavefunction is in the Slater-Jastrow form

$`\varphi _{old}^0`$ $`=`$ $`det[\mathrm{exp}[i𝐤_i𝐫_j]]`$ (49)
$`\mathrm{exp}\left[{\displaystyle \frac{1}{V}}{\displaystyle \underset{k}{}}{\displaystyle \frac{M_k^2\rho _k^e\rho _k^e+2M_kP_k\rho _k^e\rho _k^p}{2\mathrm{}V_k}}\right]`$
up to a constant factor. Instead of using $`M_k=(\stackrel{~}{v}_k)^{1/2}\theta (k_ck)`$ for the long wavelength part up to $`k_c`$, and optimizing the cut-off $`k_c`$, as done in the original work of Bohm and Pines, we can use $`\stackrel{~}{u}_k^{ee}`$ and $`\stackrel{~}{u}_k^{ep}`$ for the electron-electron and electron-proton Jastrow part taken in the RPA approximationdmc81  and relate these functions to $`M_k`$ and $`P_k`$. The resulting residual electron-electron and electron-proton interaction is screened, since $`M_k^2\stackrel{~}{v}_k`$ and $`M_kP_k\stackrel{~}{v}_k`$ in the long wavelength limit $`k0`$.

A second unitary transformation using

$$S=\frac{1}{V}\underset{k,j}{}M_k\frac{𝐤𝐩_j}{m\omega _p(k)(\mathrm{}\omega _p(0)+ϵ_k)}\mathrm{\Pi }_k^{}e^{i𝐤𝐫_j}$$
(50)
eliminates $`H_{int}`$ to first order. Here, $`\omega _p(k)`$ is the plasmon frequency at wavevector $`k`$. Note that this transformation brings the wavefunction into the backflow form. Furthermore, we treat the remaining terms of the Hamiltonian perturbatively as shown in the previous subsection.

The detailed functions we used for the electron gas and for metallic hydrogen are given in the appendix and the numerical tests are given in Section V.

## IV Comparison with the Band Structure Wavefunction

In this section, we consider another approach of generating backflow functions. As in the discussion of the single particle perturbation theory in the last section, we consider a perfect lattice of protons in which a single electron moves. It is straightforward to expand the wavefunction in plane waves and obtain a precise numerical solution of the one electron problem by diagonalization of the Hamiltonian matrix. We study to what extent we can recast the “band structure” wavefunction into a backflow form. The advantage of this approach is that we are evaluating the entire non-linear effect of a lattice of protons on the electron wavefunction, or orbital, which for a perfect lattice is a Bloch wave. However, effects of electron correlation or screening are absent for this model.

As was done in Eq. (26), the exact one-electron wave function is expanded in plane waves:

$$\varphi _𝐤(𝐫)=\underset{G}{}c_{k,G}e^{i(𝐆+𝐤)𝐫}$$
(51)
where $`𝐆`$ is a reciprocal vector of the lattice and $`𝐤`$ the crystal momentum. We then obtain numerical values for $`c_{k,G}`$ by conventional diagonalization of the Hamiltonian in this basis.

First, we study the wavefunction at $`𝐤=0`$ to determine the pair part of the wavefunction, $`U`$. Neglecting the three-body term we have:

$$\underset{i}{}u(|𝐫𝐳_i|)=\mathrm{ln}(\varphi _0(𝐫))$$
(52)
where $`𝐳_i`$ are the proton positions. Then by fourier transforming and assuming a Bravais lattice :

$$u_𝐆=_Vd^3re^{i\mathrm{𝐆𝐫}}\mathrm{ln}(\varphi _0(𝐫)).$$
(53)
This is shown in Fig. 1 and compared to the RPA form (solid line) and cumulant form. Note that we only obtain information about $`u_𝐤`$ at values of $`𝐤`$ on the reciprocal lattice. It is seen that except for the first few reciprocal lattice vectors, the pair wavefunction is determined by the cusp behavior. The non-cusp behavior is due to the neglect of higher order terms in the cumulant expansion. Some effects are picked up by the three-body term of the wavefunction. We note that even for the largest lattice vector, the values seem to follow a smooth curve, independent of the lattice directions. The $`𝐤=0`$ component, though important, will not affect the many-body nodal structure or the correlation effects near the fermi surface.

Now, let us use the same procedure to estimate the backflow function. First, we divide out the wavefunction at $`𝐤=0`$ to define the backflow functions :

$$𝐤Y_k(𝐫)=𝐤𝐫+i\mathrm{ln}[\varphi _𝐤/\varphi _o].$$
(54)
Assuming $`Y_k(𝐫)`$ is the sum of contributions of proton-electron terms on a Bravais lattice we get:

$$𝐤Y_k(𝐫)=\frac{1}{V}\underset{𝐆}{}i𝐆𝐤y_𝐆^ke^{i𝐆𝐫}.$$
(55)
Setting these two expressions equal and taking the Fourier transform we arrive at:

$$y_q^k=\frac{i}{𝐤𝐪}𝑑𝐫e^{i\mathrm{𝐪𝐫}}\mathrm{ln}[\varphi _𝐤(𝐫)/\varphi _o(𝐫)].$$
(56)
In general, the function $`y_q^k`$ depends on both $`𝐤`$ and $`𝐪`$. For small values of $`k`$, the ratio approaches a limit, independent of both the magnitude and direction of $`𝐤`$. As with the pair term, we can only determine $`y_𝐤`$ at reciprocal lattice vectors, $`𝐪`$ and for $`𝐤`$ in the first Brillouin zone. Shown in Fig. (3) is the ratio for several values of $`𝐤`$ evaluated for a simple cubic lattice plotted versus $`𝐪`$. The dispersion of the values from a smooth curve is a test of the extent to which the band structure orbital can be cast into the form of a backflow function. Note that only for the smallest values of $`𝐪`$ is the backflow function appreciable. At intermediate values of $`q`$ one does seem some effect of ”non-backflow” behavior, however it is not clear how important these effects are. At large $`q`$, we see the behavior $`y_q8\pi e^2/\lambda q^6`$ shown as the solid line as expected from the results of Section II and III.

In figure 4 is shown the error in the band energies with a backflow wavefunction (BF) and the results for having no backflow effects. For the comparison we used a BF function $`y_q=y_0\mathrm{exp}(bq)`$ fitted to the low $`q`$ behavior. Since $`y_q`$ drops off rapidly w.r.t $`q`$, it is primarily the effects at small $`k`$ that are important to describefootnoteBT . By definition the energies are identical at $`k=0`$ and the curvature around $`𝐤=0`$ is exactly put in by the backflow ansatz, at least assuming cubic symmetry. We see that the errors in the band energy go as $`k^4`$ instead of $`k^2`$ for the non-backflow trial function. However, near the band edge there are serious problems because our assumed form does not have mixing of the bands required by lattice periodicity. We expect such an effect to be much reduced for a disordered system since such degeneracies will not occur.

This achieves our goal or showing that the dominant band structure effects can be interpreted as backflow corrections, particularly at small $`𝐤`$. This implies that the changes in the nodal surfaces due to an external potential of protons are well approximated by backflow functions.

The backflow form is a much more succinct description of the single body wavefunction than the expansion in plane waves. In the introduction, we emphasized that this improves performance because we no longer have to perform the band structure calculation. However, there is also an improvement in speed of calculation of the orbitals using backflow. The expansion in plane waves can be quite slow, since the accuracy versus number of terms decreases quite slowly. In previous work on metallic hydrogen natoli93 , we divided the band structure orbital by an electron-proton Jastrow factor as an approximation to $`\varphi _0(r)`$, and then re-expanded in plane waves. The resulting expansion is much more quickly convergent in the number of plane waves since the cusp at $`r_{ep}=0`$ is in the Jastrow factor. It takes the sum of many plane waves to recover this non-analytic behavior at $`𝐫_i=𝐳_j`$. Backflow takes this even further by using the fact that near $`𝐤=0`$ the wavefunction can be expanded in pair terms with a higher-order cusp. These pair terms can be conveniently and rapidly computed, since much of the computational effort is to map each pair of particles (ee or ep) onto a grid value for a table look-up. The distances and grid values are then used for all of the pair terms: the potential, the Jastrow, the backflow and the polarization terms.

The problems concerning degeneracies of the unperturbed plane wave functions near the edge of the Brillouin zone are common to all analytical approaches considered up to now. Without a separate treatment of (nearly) degenerate zeroth order (plane wave) states, neither the cumulant method (Section II) nor perturbation theory (Section III) are able to produce the resulting energy splitting at the band edge. A degenerate case will have to be treated by including all of the degenerate states in the unperturbed basis.

## V Quantum Monte Carlo Testing of Trial Function Forms

There are two principal simulation methods used to calculate the ground state energies of quantum many-body systems: Variational Monte Carlo (VMC) and Diffusion Monte Carlo (DMC). In VMC, one samples the square of the wavefunction, and, in DMC, one uses a trial wavefunction and the imaginary-time evolution to project onto the ground state. VMC is potentially very powerful because one can use any wavefunction, as long as one can easily compute its values. One can add correlation directly to the wavefunction, leading to a very compact accurate wavefunction. The resulting integrals are similar to that of the classical partition function and therefore demand a simulation algorithm to evaluate. The disadvantage of the variational approach is that one needs to use the right functional space in order to get satisfactory properties. Though DMC is much less dependent on details of the trial wavefunction than VMC, however, lacking an exact fermion algorithm, the results still depend to some extent on the positions of the node (or phase) of the trial wavefunction.

The most straightforward, and rigorous approach to determine the trial function is to propose a definite analytic form, containing some parameters, $`a`$. One then uses VMC to evaluate the variational energy $`E_V(a)`$, an upper bound to the exact energy as a function of $`a`$. One can use various techniques to optimize the parameters to obtain the lowest energy, the lowest variance or some combination of the two. Variational optimizationkwon93 ; kwon98  has determined good backflow and three-body trial functions for the electron gas in both 2 and 3 D . The disadvantage of optimization is that beyond general trends, it is hard to extract analytic behavior because of the noisy behavior of the optimization method and the restriction to a limited functional form.

Here we compare several different trial wavefunctions on two systems: the 3D electron gas and metallic hydrogen. We employ three estimators of the quality of the wavefunction: the variational energy $`E_v=\psi \widehat{}\psi `$, the variational variance $`\sigma ^2=\psi \widehat{}^2\psi E_v^2`$ and the DMC (fixed-node) energy. The first two properties are sensitive to all aspects of a wavefunction; the variance is particularly sensitive to short-range structure since the energy fluctuations are larger. However, the DMC energy is determined only by the positions of the trial function node, not by the “bosonic” part of the trial function. The VMC/DMC calculations we performed were standard onesRMPQMC . All calculations are done with Periodic Boundary Conditions (PBC), equivalent to the $`\mathrm{\Gamma }`$ point for a band structure calculation in a cubic unit cell. Hence all trial functions were real. Though twist-averaged boundary conditions (TABC) lin2002  are useful in reducing size effects, tests showed the relative accuracy of various trial functions can be determined with PBC using real trial functions.

First, we discuss the results using backflow and three-body wavefunctions on the 3D electron gas as shown in Table II. The results using analytic trial functions give results comparable to the numerically optimized backflow results of Kwon et alkwon98 . We find that for $`r_s<20`$ the analytic wavefunction have a lower VMC energy than the numerically optimized wavefunction. This is mainly due to the inclusion of the long range part of the backflow potential. For all values of $`r_s`$ the analytic wavefunctions have a lower DMC energy, implying a more accurate nodal surface than obtained by numerical optimization. For $`r_s=20`$ the numerically optimized VMC energy is lower than that of the analytic wavefunction, indicating that at least the 3-body part of the wavefunction becomes inaccurate at strong correlations.

Now we consider the use of these same trial functions for a system composed of electrons and protons. To determine the properties using the optimization method, we used the RPA form for both the electron-electron (ee) and electron-proton (ep) $`u(r)`$. We used optimized Gaussians for both the backflow and polarization terms:

$$\eta (r)=\lambda \mathrm{exp}((rr_0)^2/w^2).$$
(57)
Even though the optimal functions may have a long range tail, as shown earlier, the additional energy gained is small and we neglect the long-range terms in setting up the parameterized trial functions. An additional Gaussian (with $`r_0=0`$ so as not to change the cusp value) was added to the pair term. We did not include ee backflow or polarization terms in the wavefunction. The resulting 10-parameter wavefunction was then optimized to minimize a linear combination of its energy and variance. Shown in Fig. 3 are the optimized backflow functions compared with the cumulant value, with the analytic form and with the band structure determination. The magnitude and shape are similar, though differences are apparent.

We compare the results with three other wavefunctions. The simplest is the SJ-PW functionsdmc81 , which do not contain backflow, three-body terms and the orbitals are simple plane waves. We also used optimized Slater-Jastrow functions with orbitals from a LDA calculationnatoli93 . Finally shown are various analytic backflow calculations: one contains only ep backflow (and 3body), the others have alse ee backflow (3body) included.

Shown in Table III are both VMC and DMC calculations of various wavefunctions for metallic bcc hydrogen at $`r_s=1.31`$, a density very close to the molecular-metallic transition. While the detailed results depend on the number of particles, in general we find that the SJ-PW function is in error within VMC by about 15mH/atom while the BF is in error by about 4 mH/atom and the LDA trial function by about 2 mH/atom. Within DMC the SJ-PW is in error by 6mH/atom and the BF is as accurate at the LDA trialfunction within the statistical error. This analysis of errors is done with the assumption that the LDA-DMC energy is exact. As another indication of the quality, the VMC wavefunction variance is roughly a factor of 3 smaller with the BF wavefunction than with the SJ-PW wavefunction.

We see that for $`N=16`$ the DMC backflow results are even lower than the LDA function. One reason for this could be that $`N=16`$ has a degenerate ground state for a single Slater determinant; many-body effects break the degeneracy. It may be that the current simulations, though similar to those of Natoli, broke the degeneracy in a more favorable way and thus have a lower energy. The $`N=54`$ system has a non-degenerate ground state at the mean field level, a closed shell, so the results may be more typical. Finally degeneracy effects are probably less important at $`N=128`$ since $`N`$ is larger.

We have tested the relative importance of including ee backflow in the case of metallic hydrogen. Using ep backflow-3body only, the analytical wavefunctions give considerably higher energies compared to the numerically optimized ones. Including ee backflow in the analytical forms, they become comparable. One should note that the analytical approaches derive ee and ep backflow at the same order of approximation; dropping one of them alone is not justified and might explain the importance of including ee and ep backflow in the analytical functions. The inclusion of 3-body terms does not noticeably affect the energies. This is similar to the results for the electron gas at comparable densitieskwon98 . Since the density is close to the transition from metallic to molecular hydrogen we tried to improve our wavefunction by considering the effects of a simple electron-proton bound state on Jastrow and backflow in our analytical formulas (see appendix) and found significantly lower energies within VMC.

Note that when ee backflow is included, it becomes necessary to move all electrons together, and for reasonable acceptance ratios one must choose an increasingly smaller time step as the system size increases. However, the more accurate nodal surface gives both a quantitative improvement in properties and a qualitative changes in some properties such as Fermi liquid parameterskwon94 .

We also used the CEIMCdewing01  method to generate a collection of proton positions appropriate to liquid metallic hydrogen at 5000K, far above the melting temperature of the lattice. Using these configurations we tested the accuracy of the same trial functions described above. See Table IV. The values marked BF3-O are obtained minimizing local energy and variance for 1000 different equilibrium configurations. We compared to the other ways of determining the backflow functions, either the analytic formulas (see the Appendix) or optimized on the lattice. While the optimized BF3 functions have a slightly lower energy in some cases, this does not compensate for the difficultly and reliability of performing the optimization. We find that the BF3 wavefunctions are about 20mH/atom lower in energy that from the SJ-PW at the VMC level, and have a lower variance. This comparison shows that disorder weakly affects the determined functions at least in this experiment. This supports our belief that the BF3 wavefunction is “transferable” to a variety of protonic configurations. In addition, we expect the backflow wavefunction to be more effective in the disordered system, since the energy degeneracies caused by crystal symmetry of a perfect lattice are not present. Comparisons using optimized LDA functions to support this hypothesis will be reported in a future publication.

## VI Conclusion

What we have shown in this paper is that ideas from perturbation theory can be used to generate an explicit trial wavefunction beyond the pair level. This gives us both an insight into the form of the many-body wavefunction and a more efficient quantum Monte Carlo simulation for disordered systems. This approach has also given intuition on the effect of an external potential on the wavefunction, even for a single electron. We have shown that one can approximate the band wavefunction (a 3d table of numbers for each Bloch wave), with three 1D functions ($`u,w,`$ and $`y`$) valid for all Bloch waves achieving reasonable accuracy. It should be recalled that for the electron-proton system, there will be these 3 functions for the ee interaction and 3 functions for the ep interaction. We have found analytical representations of these functions accurate throughout most of the phase diagram of the electron gas and promising for metallic hydrogen.

An important consideration in Monte Carlo is computational efficiency. For electron-electron backflow, the code runs slower due to having to move all the particles together. For electron-proton backflow that is not the case. You can still move electrons one at a time since all the changes in the Slater matrix are confined to a single column; each such matrix value is given by a term of the same form as a classical force, allowing it to be quickly computed once the ep distance have been computed. Expansions of single body orbitals in a plane wave basis can be quite time-consuming, especially when pseudopotentials are not used.

But the most important advantage of the backflow wavefunctions is that the form can be easily extended to put in effects of electron-electron correlation on the nodes. The outstanding problem in the simulation of quantum systems is the “fermion sign problem.” If the nodal surfaces are accurately approximated, then the “fixed-node” method will give accurate results. The present work, establishes new analytic properties of the backflow functions and thus leads to important progress in understanding nodal surfaces. In particular, the effect of long-range interactions resulting in perturbations to the nodal surfaces is important to establish. Strong short-range effects can be captured either by energy minimization or by the nodal release algorithm, which can solve for the exact wavefunction for relatively short projection times or for small numbers of fermions.ca80 . Fixing the relationship between the long-wavelength collective coordinates and the nodal surfaces could be crucial in obtaining accurate simulations for fermion systems.

In the above, we have discussed the use of backflow functions for simple metals, using plane waves as the reference state. It is straightforward to apply the approaches explored here to an insulating state. In that case, the reference state will be a determinant of Wannier functions, in the simplest case, Gaussians. The backflow ideas are applicable for suggesting improvements to the resulting Slater-Jastrow function. This will be considered in future work. A related problem is how to treat bound states in metallic liquid hydrogen in a more accurate way.

Backflow ideas are also useful at finite temperature. In that case we need to know how density matrices will evolve going from high temperature to low temperaturedmc96 . One knows how to put in backflow at high temperature. The challenge is to smoothly interpolate to zero temperature since it is clear that the backflow potential must be a smooth function of temperature. In the variational density matrix methodvdm  one uses a Hartree-Fock approach with a Gaussian basis to determined the evolution of the nodal surface of the many body density matrix. The various approaches we have described here in particular the Bohm-Pines method, will be useful in understanding the temperature dependence.

Another important problem is to generalize these methods to treat electrons with core states. The formalism should generate good trial functions in the valence region and can be used with either all-electron methods, or pseudopotentials in that region. We hope that with some modification the procedures we have discussed will be useful in the core region as well.

The hospitality of M. Mareschal and CECAM is gratefully acknowledged. This research was funded by NSF DMR01-04399, NASA, DOE (CSAR) the CNRS-University of Illinois exchange agreement, the CNRS, ENS-Lyon and the Dept. of Physics at the University of Illinois. We acknowledge useful discussions with K. Schmidt, S.-W. Zhang, R. M. Martin, R. Needs and G. Bachelet. Computational resources were provided by CECAM, NCSA and CINECA (Italy) through the INFN Parallel Computing Initiative.

## Appendix: Analytic expressions of the trial wavefunction

In this section we summarize the analytic two-body, backflow and polarization functions which describe the trial functions. We start from the pair-product (Slater-Jastrow) wavefunction based on the RPA approximation, using

$$2n\stackrel{~}{u}_q^{ee}=1+\left(1+\frac{2n\stackrel{~}{v}_q}{\epsilon _q}\right)^{1/2}$$
(58)
and

$$2n\stackrel{~}{u}_q^{ep}=\frac{2n\stackrel{~}{v}_q}{\epsilon _q\left(1+2n\stackrel{~}{v}_q/\epsilon _q\right)^{1/2}},$$
(59)
where $`\epsilon _q=\mathrm{}^2q^2/2m\lambda q^2`$. Here $`m`$ is the electron mass and $`n`$ is the electronic density. Using a trial function with ee and ep Jastrow factors corresponds to the following extended Hamiltonian, Eq.(43), with:

$`M_q^2=(\stackrel{~}{u}_q^{ee})^22n\epsilon _q`$ (60)
$`M_qP_q=\stackrel{~}{u}_q^{ep}\stackrel{~}{u}_q^{ee}2n\epsilon _q`$ (61)
Applying the unitary transformation (50) to the wavefunction, it generates the backflow potentials,

$`y_q^{ee,\mathrm{int}}`$ $`=`$ $`{\displaystyle \frac{2\lambda M_q^2}{\omega _p(q)(\omega _p(0)+\epsilon _q)}}`$ (62)
$`y_q^{ep,\mathrm{int}}`$ $`=`$ $`{\displaystyle \frac{2\lambda M_qP_q}{\omega _p(q)(\omega _p(0)+\epsilon _q)}}`$ (63)
where we used $`\omega _p^2(q)=8\pi \lambda e^2n+2.4k_F^2\lambda \epsilon _q+\epsilon _q^2`$ for the plasma frequencies and $`k_F`$ is the Fermi vector. The screened interaction between electrons, Eq.(44), and between electrons and protons, Eq.(47), can be treated by perturbation theory. Summing up the particle-hole (bubble) diagrams, using only the zeroth order plane waves, leads to coefficients $`\alpha _{k_1,k_2,q}`$ for the electron gas, as given by Eq.(III.2), but with an effective interaction and dielectric constant:

$$\stackrel{~}{v}_{\mathrm{eff}}^{ee}(q)=\stackrel{~}{v}_qM_q^2,ϵ_{\mathrm{eff}}(q,\omega )=1\stackrel{~}{v}_{\mathrm{eff}}^{ee}(q)D(q,\omega )$$
(64)
where and $`D(q,\omega )`$ is the real part of the Lindhard function. As the pair product form already accounts for plasmons, we do not consider any additional plasmon contributions to Eq. (III.2). Expanding Eq. (III.2) around $`k_i=0`$, we obtain:

$$y_q^{ee,sr}=\frac{[S(q)]^2}{2q^2}\frac{\stackrel{~}{v}_{\mathrm{eff}}^{ee}(q)}{\epsilon _qϵ_{\mathrm{eff}}(q,\epsilon _q)}$$
(65)
to obtain this formula we have further approximated the sum over occupied (unoccupied) states by $`[S(q)]^2/4`$ where $`S(q)`$ is the ideal gas structure factor,

$`S(q)=\{\begin{array}{cc}\frac{1}{2}\left[\frac{3q}{2k_F}\left(\frac{q}{2k_F}\right)^3\right]:\hfill & q<2k_F\hfill \\ 1:\hfill & q2k_F.\hfill \end{array}`$ (68)
The screened electron-proton interaction gives a similar term

$$y_q^{ep,sr}=\frac{2}{\epsilon _qq^2}\frac{\stackrel{~}{v}_{\mathrm{eff}}^{ep}(q)}{ϵ_{\mathrm{eff}}(q,0)}$$
(69)
with the screened electron proton interaction $`\stackrel{~}{v}_{eff}^{ep}(k)=(\stackrel{~}{v}_kM_kP_k)`$.

Adding these two contributions, the total backflow is:

$$y_q^{ee}=y_q^{ee,\mathrm{int}}+y_q^{ee,sr},$$
(70)
$$y_q^{ep}=y_q^{ep,\mathrm{int}}+y_q^{ep,sr}.$$
(71)
We also performed calculations with an additional ee Jastrow function $`\stackrel{~}{v}_{\mathrm{eff}}^{ep}(q)/(\epsilon _qϵ_{\mathrm{eff}}(q,0))`$ but this form did not lower the energy. Assuming this form disturbs the already correct limiting behavior of the Jastrow part $`u_q^{ee}`$ and $`u_q^{ep}`$ for $`q0`$ and $`q\mathrm{}`$, we took only the portion around the logarithmic singularity at $`2k_F`$, by using the following additional Jastrow factor:

$$\stackrel{~}{u}_q^{ee,\mathrm{add}}=\frac{[S(q)]^2}{4\epsilon _q}\stackrel{~}{v}_{\mathrm{eff}}^{ee}(q)\left\{\frac{1}{ϵ_{\mathrm{eff}}(q,\epsilon _q)}\frac{1}{ϵ_{\mathrm{eff}}(0,0)}\right\}$$
(72)
$$\stackrel{~}{u}_q^{ep,\mathrm{add}}=\frac{1}{\epsilon _q}\stackrel{~}{v}_{\mathrm{eff}}^{ep}(q)\left\{\frac{1}{ϵ_{\mathrm{eff}}(q,0)}\frac{1}{ϵ_{\mathrm{eff}}(0,0)}\right\}$$
(73)
We used $`\stackrel{~}{u}_q^{ep}+\stackrel{~}{u}_q^{ep,\mathrm{add}}`$ for the total electron-proton Jastrow potential, but only $`\stackrel{~}{u}_q^{ee}`$, since the additional term $`\stackrel{~}{u}_q^{ee,\mathrm{add}}`$ did not improve the variational energies of the electron gas. We used the unsymmetrical form of the polarization with different left and right components given by Eq. (39) and Eq. (40):

$$w_u^{ep}(r)=u^{ep,add}(r)w_y^{ep}(r)=y^{ep}(r).$$
(74)
Analogous forms were used for the electron-electron part.

For the case of metallic hydrogen we tried to take into account the effects of a possible bound state on the electron-proton pair and backflow potential. The single electron wavefunction $`\varphi _b`$ considering only one bound state can be approximately written by

$$\varphi _b\frac{A}{\sqrt{N}}\underset{i}{}\phi _b(|𝐫𝐫_i|)$$
(75)
where $`𝐫_i`$ is the position of the $`i^{th}`$ proton and the sum extends over all $`N`$ protons. As single particle orbital we will take the hydrogen ground state, $`\phi _b=(\pi a_b^3)^{1/2}\mathrm{exp}(r/a_b)`$, with energy $`e_b=1/2ma_b^2`$; $`A1`$ is a normalization taking into account the non-zero overlap between orbitals on different sites. Using Eq. (28) we obtain for the scattering amplitude

$$f(e_b,𝐩)=\underset{i}{}\frac{Ae^{i𝐩𝐫_i}}{\sqrt{N\pi a_0^3}}\frac{4\pi e^2}{p^2+(a_0^1+k_{TF})^2}$$
(76)
where we have taken a screened Coulomb interaction $`v(r)=e^2e^{k_{TF}r}/r`$ with the Thomas-Fermi wavevector $`k_{TF}^2=2k_Fe^2/\pi \lambda `$, and we have neglected overlap effects from different sites. From Eq. (27) we can finally derive the corrections to the pair potential,

$`u_q^{ep,b}`$ $`=`$ $`{\displaystyle \frac{A}{\sqrt{n\pi a_b^3}(1+(q/2k_{TF})^2)}}`$ (77)
$`\times {\displaystyle \frac{4\pi e^2}{[q^2+(a_b^1+k_{TF})^2][e_b+e_q]}}`$
and

$`y_q^{ep,b}`$ $`=`$ $`{\displaystyle \frac{8\pi e^2A}{\sqrt{n\pi a_b^3}(1+(q/2k_{TF})^2)}}`$ (78)
$`\times ({\displaystyle \frac{1}{[q^2+(a_b^1+k_{TF})^2]^2[e_b+e_q]}}`$
$`+{\displaystyle \frac{\lambda }{[q^2+(a_b^1+k_{TF})^2][e_b+e_q]^2}})`$
for the backflow potential. We have cut-off the short-range part of the corrections by multiplying with $`[1+(q/2k_{TF})^2]^1`$ in order not to destroy the cusp conditions. Since we expect a higher energy for the ground state of the screened Coulomb interaction than for the pure Coulomb potential, we used $`a_b2a_0`$ and $`A1`$ in the numerical calculations.

All potentials were split into a short range and long range partNatoli95  in such a way as optimize the accuracy for a given r-space and k-space cutoff. The short range function is evaluated in real space and the long range part is then calculated by summing over Fourier components. Figure 5 shows numerical values of $`\eta (r)`$ for the 3D electron gas. Comparing with the same figure of Kwon et al. where these functions were numerically optimized, we see that the short ranged functions are very similar for $`r_s<10`$ but different at larger $`r_s`$. Figure 6 shows the three-body contribution of $`i^{th}`$ wavefunction. It is a rapidly increasing function of $`r_s`$ and is somewhat narrower and more structured than the numerically optimized form.