Path Integral Monte Carlo methods for superfluids Introduction

Path Integral Monte Carlo methods for superfluids • Introduction to path integrals • Path integrals for superfluids • Application: – Helium liquid, droplets, supersolid, . . • Extensions to fermion systems and superfluids

Motivation for Path Integral MC • There are difficulties with VMC and DMC – Need to find good trial functions; this becomes increasing difficult as systems get more complex, especially if one doesn’t know the correct physics. – Mixed estimator problem for properties other than the energy. • Temperature is important: e. g. finite temperature phase transitions. • PIMC makes nice connection with DMC and with other theoretical approaches and leads to concepts such as Reptation MC, understanding of bose condensation, superfluidity, exchange … • Details given in : RMP 67, 279 (1995)

Imaginary Time Path Integrals

The thermal density matrix • Find exact many-body eigenstates of H. • Probability of occupying state is exp(- E ) • All equilibrium properties can be calculated in terms of thermal o-d density matrix • Convolution theorem relates high temperature to lower temperature.

Trotter’s theorem (1959) • We can use the effects of operators separately as long as we take small enough time steps. • n is number of time slices. • is the “time-step” • We now have to evaluate the density matrix for potential and kinetic matrices by themselves: • Do by FT’s • V is “diagonal” • Error at finite n is roughly: comes from communtator

Using this for the density matrix. • We sample the distribution: Where the “primitive” link action is: • Similar to a classical integrand where each particle turns into a “polymer. ” – K. E. is spring term holding polymer together. – P. E. is inter-polymer potential. • Trace implies R 1=Rm+1 closed or ring polymers

“Distinguishable” particles • Each atom is a ring polymer; an exact representation of a quantum wavepacket in imaginary time. • Trace picture of 2 D helium. The dots represent the “start” of the path. (but all points are equivalent) • The lower the real temperature, the longer the “string” and the more spread out the wavepacket.

Quantum statistics • For quantum many-body problems, not all states are allowed: allowed are totally symmetric or antisymmetric. Statistics are the origin of BEC, superfluidity, lambda transition. • Use permutation operator to project out the correct states: • Means the path closes on itself with a permutation. R 1=PRM+1 • Too many permutations to sum over; we must sample them. • PIMC task: sample path { R 1, R 2, …RM and P} with Metropolis Monte Carlo (MCMC) using “action”, S, to accept/reject.

Exchange picture • Average by sampling over all paths and over connections. • Trial moves involve reconnecting paths differently. • At the superfluid transition a “macroscopic” permutation appears. • This is reflection of bose condensation within PIMC.

3 boson example • Suppose the 2 particle action is exact. • Make Jastrow approximation for spatial dependance (Feynman form) • Spatial distribution gives an effective attraction (bose condensation). • For 3 particles we can calculate the “permanent, ” but larger systems require us to sample it. • Anyway permutations are more physical.

Main Numerical Issues of PIMC • How to choose the action. We don’t have to use the primitive form. Higher order forms cut down on the number of slices by a factor of 10. We can solve the 2 -body problem exactly. • How to sample the paths and the permutations. Single slice moves are too slow. We move several slices at once. Permutation moves are made by exchanging 2 or more endpoints. • How to calculate properties. There are often several ways of calculating properties such as the energy. If you use the simplest algorithm, your code will run 100 s or 1000 s of times slower than necessary. Calculations of 3000 He atoms can be done on a workstation-- if you are patient. Details see: RMP 67, 279 1995.

PIMC Sampling considerations • Metropolis Monte Carlo that moves a single variable is too slow and will not generate permutations. • We need to move many time slices together • Key concept of sampling is how to sample a “bridge”: construct a path starting at R 0 and ending at Rt. • How do we sample Rt/2? GUIDING RULE. Probability is: • Do an entire path by recursion from this formula. • Related method: fourier path sampling.

Bisection method 1. Select time slices ß 2. Select permutation from possible pairs, triplets, from: R’ 3. Sample midpoints 4. Bisect again, until lowest level R 5. Accept or reject entire move 0

Helium phase diagram • Because interaction is so weak helium does not crystallize at low temperatures. Quantum exchange effects are important • Both isotopes are quantum fluids and become superfluids below a critical temperature. • One of the goals of computer simulation is to understand these states, and see how they differ from classical liquids starting from nonrelativistic Hamiltonian:

Path Integral explanation of Boson superfluidity • Exchange can occur when thermal wavelength is greater than interparticle spacing • Localization in a solid or glass can prevent exchange. • Macroscopic exchange (long permutation cycles) is the underlying phenomena leading to: – Phase transition: bump in specific heat: entropy of long cycles – Superfluidity winding paths – Offdiagonal long range order--momentum condensation separation of cut ends – Absence of excitations (gaps) • Some systems exhibit some but not all of these features. • Helium is not the only superfluid. (2001 Nobel Prize for BEC)

T=2. 5 K T=1 K Normal “atomic” state “entangled” liquid

ENERGY Bose statistics have a small effect on the energy Below 1. 5 K 4 He is in the ground state. SPECIFIC HEAT • Characteristic shape when permutations become macroscopic • Finite size effects cause rounding above transition Kinetic term becomes smaller because Ncycke<N. Springs stretched more.

Superfluidity and PIMC 4 He rotating disks: Andronikashvili’s expt (1946) • We define superfluidity as a linear response to a velocity perturbation (the energy needed to rotate the system) “NCRI=nonclassical rotational inertia” • To evaluate with Path Integrals, we use the Hamiltonian in rotating frame: A = signed area of imaginary-time paths

Winding numbers in periodic boundary conditions • Distort annulus • The area becomes the winding (average center of mass velocity) • The superfluid density is now estimated as: • Exact linear response formula. (analogous to relation between ~<M 2> for Ising model. • Relates topological property of paths to dynamical response. Explains why superfluid is “protected. ” • Imaginary time dynamics is related to real time response. • How the paths are connected is more important than static correlations.

Superfluidity in pure Droplets • 64 atom droplet goes into the superfluid state in temperature range 1 K <T <2 K. NOT A PHASE TRANSITION! • But almost completely superfluid at 0. 4 K (according to response criteria. ) • Superfluidity of small droplets recently verified. Sindzingre et al 1990 Bulk experiment Droplet

Droplets and PIMC E. Draeger (LLNL) D. Ceperley(UIUC) • Provide precise microscopic probes for phenomenon such as superfluidity and vortices. • Provide a nearly ideal “spectroscopic matrix” for studying molecular species which may be unstable or weakly interacting in the gas phase. • PIMC can be used to simulate 4 He droplets of up to 1000 atoms, at finite temperatures containing impurities, calculating the density distributions, shape deformations and superfluid density. • Droplets are well-suited to take advantage of the strengths of PIMC: – Finite temperature (T=0. 38 K) – Bose statistics (no sign problem) – Finite size effects are interesting.

Experimental Setup for He droplets Toennies and Vilesov, Ann. Rev. Phys. Chem. 49, 1 (1998) • Adiabatic expansion cools helium to below the critical point, forming droplets. • Droplets then cool by evaporation to: T=0. 38 K, T=0. 15 K, (4 He) (3 He) • The droplets are sent through a scattering chamber to pick up impurities, and are detected either with a mass spectrometer with electron-impact ionizer or a bolometer. • Spectroscopy yields the rotational-vibrational spectrum for the impurity to accuracy of 0. 01/cm. Almost free rotation in superfluid helium but increase of MOI of rotating impurities.

Demonstration of droplet superfluidity Grebenev, Toennies, Vilesov: Science 279, 2083 (1998) • An OCS molecule in a 4 He droplet shows rotational bands corresponding to free rotation, with an increased moment of inertia (2. 7 times higher) 4 He are “coat” the impurity allowing it to freely rotate in the superfluid • They replaced boson 4 He with fermion 3 He. If Bose statistics are important, then rotational bands should disappear. they didn’t! • However, commercial 3 He has 4 He impurities, which would be more strongly attracted to an impurity. How much 4 He does it take to “coat” the impurity and get free rotation? • They found that it takes around 60 4 He atoms.

4 He Small impurities of 4 He in 3 He 4 He is more strongly attracted to impurity because of zero point effects, so it coats the impurity, insulating it from the 3 He.

Local Superfluid Density Estimator Although superfluid response is a non-local property, we can calculate the local contribution to the total response. Where A is the area. A 1 approximation: use only diagonal terms not positive definite could be noisy each bead contributes

Local Superfluid Reduction Average over region with cylindrical symmetry T=0. 38 K, N=500 (HCN)3 The local superfluid estimator shows a decrease in the superfluid response throughout the first solvation layer.

Bose condensation • BEC is the macroscopic occupation of a single quantum state (e. g. momentum distribution in the bulk liquid). • The one particle density matrix is defined in terms of open paths: • We cannot calculate n(r, s) on the diagonal. We need one open path, which can then exchange with others. • Condensate fraction is probability of the ends being widely separated versus localized. ODLRO (off-diagonal long range order) (The FT of a constant is a delta function. ) • The condensate fraction gives the linear response of the system to another superfluid.

How to calculate n(r) 1. Take diagonal paths and find probability of displacing one end. • advantage: – simultaneous with other averages, – all time slices and particle contribute. • disadvantage: unreliable for r>. 2. Do simulation off the diagonal and measure end-end distribution. Will get condensate when free end hooks onto a long exchange. • advantage: works for any r • Disadvantage: – Offdiagonal simulation not good for other properties – Normalization problem.

Comparison with experiment Single particle density matrix Neutron scattering cross section Condensate fraction

Surface of Liquid Helium 2 possible pictures of the surface • Dilute Bose gas model: Griffin and Stringari, PRL 76, 259 (1996). bulk: n 0 = 9% 4 He surface: n 0 = 100% • Ripplon model: Galli and Reatto, J. Phys. CM 12, 6009 (2000). Can smeared density profile be caused by ripplons alone? 4 He surface: n 0 = 50% Density profile does not distinguish

Condensate Fraction at the Surface of 4 He • PIMC supports dilute gas model. • Jastrow wavefunction • Shadow wavefunction

Dictionary of the Quantum-Classical Isomorphism Properties of a quantum system are mapped into properties of the fictitious polymer system Attention: some words have opposite meanings. Quantum Classical Bose condensation Delocalization of ends Boson statistics Joining of polymers Exchange frequency Free energy to link polymers Free energy Imaginary velocity Bond vector Kinetic energy Negative spring energy Momentum distribution FT of end-end distribution Particle Ring polymer Potential energy Iso-time potential Superfluid state Macroscopic polymer Temperature Polymer length

“Direct” Fermion Path Integrals • Path integrals map quantum mechanics into a system of crosslinking closed “polymers. ” 4 quantum paths R 0=PRM, P permutation, S(Ri, Ri+1) is “boltzmannon action” • Bosons are easy: simply sample P. • Fermions: sample the “action” and carry (-1)P as a weight. • Observable is even P - odd P. scales exponentially in N and T-1!

To get around the sign problem, why not just use the fixed -node method? What nodes? The ground state nodes are not necessarily the correct ones at T>>0. The nodes of the density matrix have an imaginary time dependence. High temperature Low temperature

Fixed-Node method with PIMC • Get rid of negative walks by canceling them with positive walks. We can do this if we know where the density matrix changes sign. Restrict walks to those that stay on the same side of the node. • Fixed-node identity. Gives exact solution if we know the places where the density matrix changes sign: the nodes. • Classical correspondence exists!! • Problem: fermion density matrix appears on both sides of the equation. We need nodes to find the density matrix. • But still useful approach. (In classical world we don’t know V(R). )

Proof of the fixed node method 1. The density matrix satisfies the Bloch equation with initial conditions. 2. 3. 4. 5. 6. One can use more general boundary conditions, not only initial conditions, because solution at the interior is uniquely determined by the exterior-just like the equivalent electrostatic problem. Suppose someone told us the surfaces where the density matrix vanishes (the nodes). Use them as boundary conditions. Putting an infinite repulsive potential at the barrier will enforce the boundary condition. Returning to PI’s, any walk trying to cross the nodes will be killed. This means that we just restrict path integrals to stay in one region. neg pos R 0 neg R PR 0

Ortho-para H 2 example In many-body systems it is hard to visualize statistics. • The simplest example of the effect of statistics is the H 2 molecule in electronic ground state. • Protons are fermions-must be antisymmetric. 1. Spins symmetric ( ). spatial wf antisymmetric (ortho) “fermions” 2. Spins antisymmetric ( - ). spatial wf symmetic (para) “bosons” 3. Non symmetrical case (HD) “boltzmannons” All 3 cases appear in nature! • Go to relative coordinates: r= r 1 -r 2 • Assume the bond length is fixed |r|=a. Paths are on surface of sphere of radius a. PIMC task is to integrate over such paths with given symmetries. For a single molecule there is no potential term, a “ring polymer” trapped on the surface of a sphere.

Paths on a sphere 1. “boltzmannons”Ring polymers on sphere O(r r) 2. “bosons” 2 types of paths allowed. + O(r -r) O(r r) 3. ”fermions” 2 types of paths allowed O(r r) - O(r -r) Low efficiency as

Restricted paths for ortho H 2 • Fix origin of path: the reference point. • Only allow points on path with a positive density matrix. paths staying in the northern hemisphere: r(t). r(0)>0 • Clearly negative paths are thrown out. • They have cancelled against positive paths which went south and then came back north to close. • The symmetrical rule in “t”: r(t). r(t’)>0 is incorrect. • Spherical symmetry is restored by averaging over the reference point: the north pole can be anywhere. • Can do many H 2 the same way. • Ortho H 2 is much more orientable than either HD or para H 2.

RPIMC with approximate nodes • In almost all cases, we do not know the “nodal” surfaces. • We must make an an ansatz. • This means we get a fermion density matrix (function with the right symmetry) which satisfies the Bloch equation at all points except at the node. • That is, it has all the exact “bosonic” correlation • There will be a derivative mismatch across the nodal surface unless nodes are correct. • In many cases, there is a free energy bound. (proved at high temperature and at zero temperature and when energy is always lower. ) • Maybe one can find the best nodes using the variational principle. (variational density matrix approach)

Exciton superfluidity • What system is the most appropriate to observe superfluidity of fermions? (strongest pairing) • Consider the simplest 1 -band model of particles and holes in a semiconductor. • Assume masses are isotropic and the same; only the charge is different. • At low temperature a particle and hole can bind together to form an exciton (like a hydrogen atom) which is a boson. • If the exciton density is high enough, they can bose condense. : T 3/2<2. 7 • Shumway-Ceperley (1999) observed this transition for excitons. • What do the paths look like? • Observed in Oct 2003 in atom traps! BCS region

Pairing Nodes • Free fermion nodes does not allow pairing because nodes of two species are independent. – Consider two pairs of fermions. – Possible exchanges are {I, Pa. Pb} and {Pa , Pb}. – The permutation Pa. Pb represents an exciton exchange but it is forbidden if nodes are independent since the path will cross “a” nodes or “b” nodes first. • Instead we used paired nodes: A{ g(a 1 -b 1)g(a 2 b 2)…. . } where g(r) is a pairing function (we used a Gaussian). • Nodes are time-independent winding number formula for superfluid response. • We can define 2 different responses. Let Wx be winding number of species x. – movement of walls: <(Wa+Wb)2> – magnetic field: <(Wa-Wb)2>

BEC of excitons Winding exchange (3, 6) [rs=6 T<Tc] Pair exchange (2, 2) Blue=electron lavendar=hole Superconductivity is (cooper pairing) of paths.

Theoretical • Restricted paths allow realistic calculations of many fermion systems. No sign problem. • Generalization of bosonic PIMC. • Unifies theory of bose and fermi systems: ring exchanges important for both. • Makes a nodal assumption which is only controlled for T>TF Practical • No Born-Oppenheimer approximation • Fully quantum protons • No density functional needed • No pseudopotential/k-space cutoff • No empirical potentials/chemical model • Paths get stuck at low temperature T<0. 1 TF unless you make other assumptions (e. g. ground state nodes. ) • CEIMC allows lower temperature PI simulations.