An Exact Monte Carlo Method for Continuum Fermion Systems M. H. Kalos1 and Francesco Pederiva2 1Lawrence Livermore National Laboratory Livermore, CA 94550 USA 2Dipartimento di Fisica, Universita` di Trento, I-38050 Povo-Trento, Italy (February 1, 2008) 0 0 0 2 n Abstract a J 3 1 1 We offer a new proposal for the Monte Carlo treatment of many-fermion v 3 9 systems in continuous space. It is based upon Diffusion Monte Carlo with 1 1 significant modifications: correlated pairs of random walkers that carry oppo- 0 0 site signs; different functions “guide” walkers of different signs; the Gaussians 0 / t a used for members of a pair are correlated; walkers can cancel so as to con- m - serve their expected future contributions. We report results for free-fermion d n systems and a fermion fluid with 14 3He atoms, where it proves stable and o c : correct. Its computational complexity grows withparticle number,butslowly v i X enough to make interesting physics within reach of contemporary computers. r a Typeset using REVTEX 1 Monte Carlo methods have provided powerful numerical tools for quantum many-body physics [1,2]. They include methods such as Green’s function Monte Carlo (GFMC) [3], Diffusion Monte Carlo (DMC) [4], or Path Integral Monte Carlo (PIMC) [5] that are ca- pable of giving, at least for moderate size bosonic systems, answers with no uncontrolled approximations. Such accurate treatment of fermionic systems has been made vastly more difficult by a “sign problem.” Progress in the application of Quantum Monte Carlo methods to condensed matter physics, to electonic structure, and to nuclear structure physics has been impeded for years by the lack of exact and efficient methods for dealing with fermions. This paper offers a new proposal for solving many-fermion systems by an extension of DMC. In the systems we have studied, the signal-to-noise ratio of the Monte Carlo esti- mates are constant at long imaginary times, by contrast to the behavior of ordinary DMC where they decay exponentially [2]. Except for the use of a short-time Green’s function, no approximations– physical, mathematical, or numerical– are made. The effect of a finite interval of imaginary time is easily controlled, or may be eliminated entirely. ItisnosurprisethatMonteCarlomethodscansolvethemany-bodySchrodingerequation ~ in imaginary time for bosonic systems. Let R denote all coordinates of an N-body system, and V(R~) be the potential at R~. h¯2 ∂ψ(R~,τ) 2 ~ ~ [− ∇ +V(R)]ψ(R,τ)+h¯ = 0 (1) 2m ∂τ This equation also describes the diffusion of an object (a “random walker”) in a 3N di- ~ mension space in which the potential V(R) serves as a generalized absorption rate. Because the potential in physical problems can be unbounded from above and below, a direct simula- tion of that diffusion, although straightforward, will be inefficient. Some form of importance sampling transformation has been found to be highly useful. In the standard DMC [3,4], this is a technical device for accelerating the convergence; in our new method it becomes an essential feature. DMC uses an “importance” or “guiding” function ψ (R~) and a trial eigenvalue E to G T constructarandomwalk. Asimpleversionisasfollows: Usingafixedstepinimaginarytime, 2 δτ, a walker at R~ is (a) moved to R~+δτ∇~ lnψ (R~); (b) then each coordinate is incremented G by an element of a random vector U~, a Gaussian with mean zero and variance δτ; finally, ˆ ~ ~ (c) each walker is turned into M walkers with < M >= exp δτ E −Hψ (R)/ψ (R) , T G G n h io where Hˆ is the Hamiltonian. The resulting random walk has expected density f(R~,τ) = ψ (R~) a exp[(E −E )τ]φ (R~) (2) G k T k k k X where φ (R~) are eigenfunctions of Hˆ with eigenvalues E , and a are expansion coefficients. k k k The limit of f(R~,τ) for large τ is dominated by the eigenfunction φ having the lowest 0 eigenvalue E . 0 We alter the structure of DMC in the following ways: (1) In order to represent an anti- symmetricwavefunctionthatisbothpositiveandnegative, weintroducewalkers, {R~+,R~−}, m m that respectively add or subtract their contributions to statistical expectations. The com- putation now involves ensembles of pairs of walkers carrying opposite signs. (2)Two distinct functions, ψ±(R~) are used to guide the ± walkers. (3) The Gaussians U~± for the paired G walkers are not independent; rather U~− is obtained by reflecting U~+ in the perpendicular bisector of the vector R~+ − R~−. Finally (4) the overlapping distributions that determine the next values of R~± are added algebraically so as to allow positive and negative walkers to cancel, while preserving the correct expected values. In a Monte Carlo calculation of this kind, we “project” quantities of interest by calculat- ing integrals weighted with some trial function, say ψ (R~). In DMC the energy eigenvalue, T E , can be determined from: 0 ˆ ~ Hψ (R ) T m E0 = HˆψT(R~)φ0(R~)dR~ = Xm ψG(R~m) (3) R ψT(R~)φ0(R~)dR~ ψT(R~m) ~ R m ψG(Rm) X replacing integrals by sums over positions of the random walk. In our modified dynamics, Eq.(3) now takes the form 3 Hˆψ (R~+) Hˆψ (R~−) [ T m − T m ] ψ+(R~+) ψ−(R~−) E0 = Xm G m G m . (4) ψ (R~+) ψ (R~−) [ T m − T m ] ψ+(R~+) ψ−(R~−) Xm G m G m If{R~+}and{R~−}followthesamedynamics using thesame ψ , thentheexpected values m m G of numerator and denominator of Eq. (4) decay exponentially at large τ. Some correlation among walkers is essential. This observation is reinforced by noting that the Pauli principle, which demands that fermion wave functions be antisymmetric, is a global condition, and cannot be satisfied by independent walkers. Put another way, it will be necessary to have dynamics that distinguish between walkers that carry different signs. These motivations underlie aspects (2) and (3) of our method outlined above. A second aspect of the difficulty in treating fermion systems is that the density that one obtains naturally from a random walk is the symmetric ground state. In order for Eq. (4) to have an asymptotically bounded signal-to-noise ratio, walkers of opposite signs must be able efficiently to cancel. This underlies the need for modification (4) given above. The need for some degree of cancellation has been a theme of previous research starting with the work of Arnow et al. [7]. The need for distinct dynamics for positive and negative walkers was stressed in [8]. That these two aims could be accomplished by appropriate correlation among walkers was pointed out by Liu, Zhang, and Kalos [9]. The use of distinct guiding functions is new and serves as the connection among the different algorithmic ideas that enables the treatment of general potentials. Stable results can be obtained using correlated pairs only. Correct results are ensured when the dynamics have the property that the random walk for either member of the pair is the same as that of a single free walker, except when they cancel, a condition satisfied here. The expectations of Eqs.(3) and (4) are linear in the walker densities and are unchanged by correlations. Furthermore, we have devised a method of canceling opposite walkers that also preserves these expectations. ~ ~ Let ϕ (R) be a trial function for the fermionic state. Let ϕ (R) be some approximation A S 4 to the symmetric ground state wave function of the same Hamiltonian. Define: ψ±(R~) = ϕ2(R~)+c2ϕ2(R~)±cϕ (R~) (5) G S A A q The following properties of these two functions are significant: (a) they are positive; (b) when c is small, they are dominated by ϕ , so that opposite walkers behave similarly; (c) S + ψ transforms under an odd permutation P as follows: G ψ+(PR~) = ψ−(R~). (6) G G As mentioned above we modify simple DMC in several ways. The “drift” is applied in the usual way to walkers assumed to be at R~±, using the two guiding functions: 0 R~+ = R~+ +δτ∇~ lnψ+(R~+) 0 G (7) R~− = R~− +δτ∇~ lnψ−(R~−). 0 G Diffusion of the walkers, however, is carried out in a correlated way: let U~+ be a vector of 3N Gaussian random variables each of mean zero and variance δτ. New trial positions R~± n are now given by R~+ = R~+ +U~+;R~− = R~− +U~−, (8) n n where the random vector U− is obtained by reflection in the perpendicular bisector of the vector R~+ −R~− as described in (4) above. Walker densities can be subtracted by computing their chances of arrival at a common point, but because they have different guiding functions, they do not exactly cancel. The analysis of “forward walking” [1,6,10] allows one to determine the expected future contribu- tion of a walker to any projected quantity. Thus, we can compute the change in expectations when a pair meets. For this change to be zero [11], a positive walker at R~+ must survive to n the next time step with probability P+(R~+;R~+,R~−) = n (9) B−(R~+|R~−)G(R~+ −R~−)ψ+(R~+) max 0,1− n n G n " B+(R~+|R~+)G(R~+ −R~+)ψ−(R~+)# n n G n 5 where exp[−(R~′ −R~)2/(2δτ)] G(R~′ −R~) = (10) (2πδτ)3N/2 is the Gaussian density used in DMC. The branching factors, B+(R~|R~+) and B−(R~|R~−) are Hψ±(R~) B±(R~) = exp δτ E − G (11)   T ψ±(R~)   G    An analogous expression is used for negative walkers.  An isolated walker may appear as a result of different branching factors at {R~+} and m {R~−}; if, with probability one half, one generates a walker of opposite sign by interchanging m the coordinates of two like-spin particles, then a pair is reconstituted that preserves future expectations. To determine the energy, we use the estimator of Eq. (4). A sharp indication of the stability of the calculation is the behavior of its denominator ψ (R~+) ψ (R~−) D = [ T m − T m ] (12) ψ+(R~+) ψ−(R~−) G m G m Inanaivecalculation, D decaystozeroinanimaginarytimeoforderτ = 1/(E −E )where c A S E and E are the fermion and boson energies. A stable method will show D asymptotically A S constant. Although a system of free fermions in a periodic box is analytically trivial, it presents an exigent test of this method. For this system, the lowest symmetric state is constant, and the exact fermionic wave function is a determinant of plane waves. We use ρ = 0.5 and set ψ±(R~) = 1+c2ϕ2(R~)±cϕ (R~) (13) G A A q where ϕ is a Slater determinant of one body orbitals χ~k of the following form: A ~ri χ~k = exp i~k · ~r +λ η(r )~r (14) ~ri   i B ij ij j6=i X    The parameter λ controls the departure of the nodal structure of this function from the B exact shape. The fact that these functions are modulated only a little from a constant by ϕ means that the polarization of the population of plus and minus walkers is small. A 6 In table I we report the results obtained for periodic systems of 7, 19, and 27 free fermions. The results agree with the analytic eigenvalues within the Monte Carlo estimates of the standarderror. It hasbeen conjectured that the computational complexity of Fermion Monte Carlo calculations will grow as N!, where N is the number of particles in the system. Since (27!/7!) = 2.16 ×1024, a calculation with 27 or even 19 bodies would be impossible were that conjecture to be true. We have also applied this algorithm to a system of 14 3He atoms in a periodic box at equilibrium density, ρ = 0.0216˚A−3. Energies are expressed in Kelvins, and lengths in ˚A. With interatomic potentials that have a hard core, we may use the same function ϕ as A for free fermions, but also need a Jastrow product. With ϕ = ϕ (R~) = exp[−(b/r )5], (15) S S ij i<j Y the guiding functions now have the form: ψ±(R~) = ϕ (R~) 1+c2ϕ2(R~)±cϕ (R~) (16) G S A A (cid:20)q (cid:21) In Fig. 1 we plot the cumulative denominator as a function of imaginary time for a typical run. As can be seen, the fundamental stability of the method is well demonstrated. Fig. 2 shows the decay of the same denominator, when the method is made unstable by setting c = 0. Table II exhibits the eigenvalues of various runs with our method applied to the periodic system with 14 3He atoms. They are all consistent, and yield a weighted average of - 2.2558(39). The run marked (b) is a continuation of the run labeled (a) separated by a long run with a longer time step. As a whole, including such continuations, the longest aggregate sequence comprises a total imaginary time of 1830 K−1. Using a total system energy difference of 20 K (as we have measured), that corresponds to 3.6 × 104 fermion decay times. An alternative measure of the length of the run, suggested by David Ceperley [12], is the ratio of the rms diffusion length of a particle to the mean spacing between particles. For this sequence of runs, that ratio is 19. Thus the observation of stable values of the sums in Eq. (4) is significant. 7 Space limitations preclude a complete description of the other checks that we have made that the results for 3He are correct: they include a fixed node calculation of exactly the same model problem, which yielded an eigenvalue of -2.08(1) K. A transient estimate (cf. Fig. 3), relaxing from the fixed node, is consistent with our result (shown as the dashed line.) Analysis of the results in Fig. 2 leads to a fermion-boson energy difference of 1.434(35) K per particle. This agrees well with a direct calculation of the energy of a 14-body mass-3 boson system that gave -3.68(1) K. By construction, the method proposed here introduces no approximations other than that of the short imaginary-time Green’s function. In other words, if the results are stable, then they are correct. Although we have have not yet proved the stability of the method (i.e. that the long-term average of the denominator of the eigenvalue quotient is not zero), we believe that we have convincingly demonstrated the stability. Perhaps the most impor- tant conclusion that we may draw is that the “sign problem” of Fermion Monte Carlo for continuous systems is not intractable; the search for elegant computational methods in this and related applications is justified. One of the authors (FP) has been supported, in part, by the National Science Founda- tion under grant ASC 9626329. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract No. W-7405-Eng.-48. We are happy to acknowledge helpful conversations with B.J. Alder, J. Carlson, D.M. Ceperley, G.V. Chester, R.Q. Hood, S.B. Libby, K.E. Schmidt, C.J. Umrigar, and S. Zhang. We thank particularly D.E. Post for his help and strong support. 8 REFERENCES [1] D.M. Ceperley and M.H. Kalos, “Quantum Many-body Problems” in Monte Carlo Methods in Statistical Physics, ed. by K. Binder, Springer-Verlag (1979). [2] K.E. Schmidt and M.H. Kalos, “Few- and Many-fermion Problems” in Monte Carlo Methods in Statistical Physics, vol. 2, ed. by K. Binder, Springer-Verlag (1984). [3] M.H. Kalos, D. Levesque, and L. Verlet, Rev. A 9, 2178 (1974). [4] P.J. Reynolds, D.M. Ceperley, B.J. Alder, and W.A. Lester Jr., J. Chem. Phys., 77, 5593 (1982). [5] D.M. Ceperley, Rev. Mod. Phys. 78, 279, 1995. [6] M.H. Kalos, Phys. Rev. A 2, 250 (1970). [7] D.M. Arnow, M.H. Kalos, M.A. Lee, K.E. Schmidt, J. Chem. Phys. 77, 5562 (1982). [8] S. Zhang and M.H. Kalos, Phys. Rev. Letts 67, 3074 (1991). [9] Z. Liu, S. Zhang, M.H. Kalos, Phys. Rev. E 50, 3220 (1994). [10] R.N. Barnett, P.J. Reynolds, and W.A. Lester, Jr., J. Comp. Phys. 96, 258 (1991). [11] An exposition of the basic method– without the current results– has been published in the proceedings of the NATO-ASI conference “Quantum Monte Carlo Methods in Physics and Chemistry”, Cornell University, M.P. Nightingale and C.J. Umrigar, eds., Kluwer 1999. [12] D.M. Ceperley, private communication, September, 1999. 9 TABLES N E E ex 7 2.91285(49) 2.912712 19 2.760(25) 2.757454 27 2.796(30) 2.763316 TABLE I. Energies and errors for a periodic system of N free fermions. The analytic result is E . ex b c λ E B 1.15 0.025 0 -2.251(08)a 1.145 0.025 0 -2.258(17) 1.135 0.025 0 -2.257(10) 1.15 0.016 0 -2.246(20) 1.15 0.010 0 -2.250(19) 1.15 0.025 0 -2.2559(84)b 1.15 0.025 -0.05 -2.249(12) 1.15 0.025 0.05 -2.268(10) TABLE II. Energies and errors for a periodic system of 14 3He atoms 10