Path Integral Ground State study of 2D solid He
We have studied a two-dimensional triangular commensurate crystal of He with the exact K Path Integral Ground State (PIGS) Monte Carlo method. We have projected onto the true ground state both a Jastrow-Nosanow wave function, in which equilibrium positions are explicitly given and no Bose–Einstein (BEC) is present, and a translationally invariant shadow wave function, in which the solid phase emerges through a spontaneously broken symmetry process and it has BEC. We find a remarkable convergence to the same properties, both the diagonal ones as well as the off–diagonal one–body density matrix . This supplies a strong evidence that no variational bias are present in the PIGS method. We find no BEC in the commensurate 2D He crystal at K, shows an exponential decay in the large distance range. The structure found in is due to virtual vacancy–interstitial pairs and this shows up in the momentum distribution.
Quantum Monte Carlo methods have provided a very powerful tool in exploring the physics of strongly interacting many-body quantum systems. As far as the properties of Bose fluids and solids are concerned, Path Integral Monte Carlo (PIMC) methods have been proved to evaluate “exact” and “unbiased” expectation values on the thermal equilibrium state at finite temperature Ceperley (); “exact” means that the obtained results are the true expectation values within the statistical error, while “unbiased” means that the only required input is the interatomic potential. At zero temperature, different “exact” techniques are available, as the Green’s function Monte Carlo Kalos1 () (GFMC), the Diffusion Monte Carlo Kalos2 () (DMC), the Reptation Monte Carlo Baroni () or the Path Integral Ground State pigs () (PIGS) methods. All such methods, though “exact” in suitable limits, rely on variational models for the ground state wave function (wf) of the system. These variational models play a relevant role in the importance sampling employed by these methods; although, in principle, the final results should not be affected by the particular choice of the trial wf, in practice, the possibility that some bias could survive, especially in systems with complex broken symmetries like in the solid phase, has not been established. In this paper we will show the robustness of the results given by the PIGS method with respect to the choice of the variational wf by studying a model for He in two dimensions in the solid phase. In fact, we find the same physical properties of the system, within the statistical noise, starting from two radically different wfs. One, a Jastrow–Nosanow Nosanow () wf (JNWF), has built in the crystal lattice via Gaussian localization factors, it is not Bose–symmetric and has no Bose–Einstein condensation (BEC). The other, a shadow Viti () wf (SWF), is translationally invariant, with crystalline order arising as spontaneous broken symmetry, and it has BEC Masserini (); Galli (); Galli3 (). We have chosen to study this model in 2D for different reasons. The reduced dimensionality allows us to study correlations to much larger distances than in 3D, even more than 20 lattice parameters. Fluctuations are expected to be stronger in 2D so this is a more stringent test for convergence given the different symmetry properties of the starting variational wfs. Finally, the 2D system is a model which is relevant for adsorbed He atoms on a smooth substrate like graphite.
With the present study we address also the important question of the supersolid state of matter Chan (); Chan2 (); Reppy (); Shira (); Kubota (); Kojima (). PIMC computations give strong evidence that in a 3D commensurate solid He (number of atoms equal to the number of lattice sites) there is no superfluid response and no BEC Clark (); worm (). The PIMC computations are at finite temperature and the lowest is of order K which is above the experimental transition temperature of crystals of good quality Chan2 (). By computing the density matrix of crystalline He at K, we find that there is no BEC in a 2D crystal. Preliminary results indicate that this is true also in 3D.
Dealing with low temperature properties, He atoms are described as structureless zero–spin bosons, interacting through a realistic two–body potential, that we assume to be the HFDHE2 Aziz potential Aziz (). The aim of the PIGS method is to improve a variationally optimized trial wf by constructing a path in the Hilbert space of the system which connects the given wf to the true ground state of the system; during this “path”, the correct correlations among the particles arise through the “imaginary time evolution operator” , where is the Hamiltonian operator. Being a trial wf with non–zero overlap with the exact ground state , this is obtained as the limit of suitably normalized. This can be written analytically by discretizing the path in imaginary time and exploiting the factorization property . In this way, turns out to be expressed in term of convolution integrals which involve the “imaginary time propagator” for a , that can be small enough such that very accurate approximants are known Ceperley (); Suzuki (). This maps the quantum system into a classical system of open polymers pigs (). An appealing feature peculiar to the PIGS method is that, in , the variational ansatz acts only as a starting point, while the full path in imaginary time is governed by , which depends only on the Hamiltonian operator.
As a trial wf, we used both a JNWF and a SWF. The JNWF is written as the product of two–body correlations and of Gaussian one–body terms which localize the particles around the assumed lattice positions Nosanow (). In the SWF, beyond the explicit two–body factors, additional correlations are introduced via auxiliary (shadow) variables which are integrated out Viti (). Nowadays, SWF gives the best variational description of solid and liquid He Moroni (). In both the cases, the variational parameters have been chosen to minimize the expectation value of the Hamiltonian operator. In what follows, we will refer to PIGS when we deal with the projection of the JNWF, and to SPIGS spigs () (Shadow Path Integral Ground State) when we project the SWF.
Because of the Bose statistics obeyed by the atoms, when using as an approximation of the true ground state , one has, in principle, also to account for permutations in the propagator Ceperley (); pigs (). Permutation moves are necessary when the JNWF, which is not Bose–symmetric, is used. On the other hand permutation moves are not necessary whenever the trial wf in is already Bose–symmetric, as the SWF. However, also for SPIGS, adding permutation moves turns out to be useful in improving the efficiency and the ergodicity of the sampling, mainly in reaching the large–distance range of . In our algorithm we have introduced two different permutation samplings: cycles of inter–particles exchanges and swap moves. The first, which may involve an arbitrary number of particles, are described in detail in Ref. scambi, . In a PIGS or SPIGS calculation of , the efficiency can be further improved by introducing particular two–particles permutations cycles involving always one of the two positions and : the swap moves worm (). These moves improve very much the efficiency of the computation of and their acceptation rate is remarkably high: in the present 2D system, by using a stagingmolotrovo () method to sample the free (kinetic) part of the imaginary time propagator, we have found an acceptation rate for this swap move which is nearly 15% in the liquid phase and in the solid phase at densities close to the melting.
The 2D He system phase diagram is known from accurate finite temperature PIMC simulations Gordillo (); at zero temperature both DMC DMC () and GFMC GFMC () have been used to investigate its properties mainly for the liquid phase. We have performed SPIGS and PIGS simulations of a 2D He commensurate triangular crystal at Å, slightly above the melting density. In order to control the reliability of our results we have tested their dependence on both the “projection time” and the “time step” . For a fixed value of K, we have done calculations with , , and K and we have used the pair–product approximation Ceperley () for the imaginary time propagator. Reducing below K affects only marginally the results; for example, by using K, the obtained energy is only 1% lower then the one with K. So we have adopted in most of the computations K as a reasonable compromise between accuracy and computational effort. These tests provide also a robust check of the ergodicity of the sampling algorithm since a lower value of for given means a bigger number of small time projections (convolution integrals) so that one is dealing with polymers of increasing length. We have then increased till convergence in the results has been achieved.
Diagonal properties, like the energy, have been computed in a box which hosts particles with periodic boundary conditions. In Fig. 1 we give the energy per particle as a function of both for SPIGS and PIGS: at K the SPIGS result is already converged to the value K, while with the PIGS one reaches convergence at K, where the energy takes the value K. The potential energy values are respectively K and K. From Fig. 1, it is evident that the true ground state value is reached much more quickly (i.e. with a lower number of small-time projections) with SPIGS. The overlap between the trial wf and the true ground state plays a crucial role in the convergence: a large overlap ensures that is a good approximation even for not too large . It has been shown overlap () that where corresponds to the area of the region between and its convergence value . From our results, turns out to be 99.8% for the SWF and 97.9% for the JNWF. From the peaks in the static structure factor we have extracted the Debye-Waller factor; we have found that both SPIGS and PIGS converge to the same value (Fig. 1) and, as in the energy case, SPIGS shows a faster convergence. It is important to note that the energy convergence does not allow to deduce a priori the convergence of other physical properties: the convergence must be checked independently for each observable.
In order to study whether the 2D commensurate solid He has BEC, we have computed the one–body density matrix . A non–zero limit of for implies BEC and the Fourier transformation of gives the momentum distribution . We have computed in a 2D commensurate crystal at Å in a simulation box with particles (this allows us to explore distances up to about 39 Å). In our calculations we have sampled in two different ways. In the first one we force to lie on a nearest–neighbor direction, while, in the second and are allowed to explore the full plane; we have found perfectly consistent results. In Fig. 2 we show the convergence of computed along the nearest–neighbor direction with SPIGS and PIGS for increasing projection time. The given by SWF has a non–zero limit at large distance, i.e. there is BEC as in 3D Galli (). When projecting from SWF the large distance plateau decreases for increasing until, at K, it disappears up to the distance allowed by the simulation box. At convergence, displays an exponential decay, with a correlation length Å, with superimposed a small modulation reflecting the crystal symmetry. The given by the JNWF is essentially a Gaussian. Under projection, exchanges among the atoms greatly extended the range of . It is remarkable that at K the obtained starting from such different wfs coincide within the statistical uncertainty. Similar conclusions are reached also from the computation of in the whole plane.
In Fig. 3a we show obtained both with SPIGS ( range) and PIGS ( range) with the projection time K. Again, the results are indistinguishable within the statistical error. We have looked also for finite size effects by considering system at the same density but with different particles number (, 240 and 480) and we have found no appreciable differences. With PIGS we have studied up to a distance of about 60 Å, finding the same exponential decay.
We conclude that a 2D commensurate He crystal has no BEC, or, more precisely, the condensate fraction, if any, is below 10. We notice that exchange processes are rather significant in the system, so that the range of is significantly larger than the size of the unit cell. This manifests itself in the momentum distribution (Fig. 3b) as deviation of from the Gaussian corresponding to the kinetic energy. There is an excess of particles at low momenta up to Å whereas there is a deficit in the –space region around Å. The positions of the bumps of over the exponential decay suggest that the extension of beyond the unit cell can be interpreted as a signal of the appearence of vacancy-interstitial pairs (VIP): their position respect to the origin correspond to interstitial positions whereas the distance between the positions of two neighboring pumps is equal to the lattice parameter. In the case of the SWF these VIP are unbound allowing for a finite BEC Galli (), whereas these pairs are bounded in the exact . This might be due to the presence of some long range correlations in , possibly caused by the zero–point motion of transverse phonons, which are absent in the SWF. Since the ground state wf reflects the zero–point motion of any excitation in the system Reatto (), it is tempting to suggest that the evidence of VIP in is the manifestation of a new kind of excitation in the quantum solid beyond the phonon excitations.
Our results that the 2D commensurate crystal has no BEC at K and the similar results Clark (); worm () obtained in 3D at finite make plausible that also in 3D a commensurate crystal has no BEC at K and, in fact, we have some preliminary results for in 3D supporting this notion. Experimentally it is established that defects and He impurities have an important role on the non–classical rotational inertia (NCRI) of He crystal and dislocations are suspected to have a relevant role Chan2 (), at least in the case of good quality crystals. It is an open question what happens in a solid with less and less defects: will any NCRI effect go away? Our result suggests that it should go away unless some disorder is present even in the ground state as an intrinsic property of . A key question is therefore to establish whether zero–point defects are present in the ground state of a quantum solid. We believe that the presence of zero–point defects in solid He is still an open question Galli2 (). Since a SWF, which has BEC, is the exact ground state of a suitable Hamiltonian, the interesting question is for which class of interatomic potentials does the commensurate crystal have BEC and supersolidity?
The convergence in the results obtained starting from radically different wfs is a remarkable result because supplies strong evidence for the absence of any variational bias in the PIGS method. Moreover, even if not shown in this letter, we have also obtained convergence of diagonal properties, such as the static structure factor and the radial distribution function, by projecting a simple Jastrow wf which has just the minimal information on the short range behavior and displays no crystalline order at the considered density.
This work was supported by the INFM Parallel Computing Initiative, by the Supercomputing facilities of CILEA and by the Mathematics Department “F. Enriques” of the Università degli Studi di Milano.
- D.M. Ceperley Rev. Mod. Phys. 67, (1995).
- P.A. Whitlock, D.M. Ceperley, G.V. Chester and M.H. Kalos, Phys. Rev. B 19, 5598 (1979).
- M.H. Kalos, D. Levesque and L. Verlet, Phys. Rev. A 9, 2178 (1974).
- S. Baroni and S. Moroni, Phys. Rev. Lett. 82, 4745 (1999).
- A. Sarsa, K.E. Schmidt and W.R. Magro, J. Chem. Phys. 113, 1366 (2000).
- L.H. Nosanow, Phys. Rev. Lett. 13, 270 (1964).
- S.A. Vitiello, K. Runge, and M.H. Kalos, Phys. Rev. Lett. 60, 1970 (1988).
- L. Reatto and G.L. Masserini, Phys. Rev. B 38, 4516 (1988).
- D.E. Galli, M. Rossi and L. Reatto, Phys. Rev. B 71, 140506 (2005).
- D.E. Galli and L. Reatto, J. Low Temp. Phys. 124, 197 (2001).
- E. Kim and M.H.W. Chan, Science 305, 1941 (2004); Nature 427, 225 (2004); J. Low Temp. Phys. 138, 859 (2005); Phys. Rev. Lett. 97, 115302 (2006).
- X. Lin, A.C. Clark, and M.H.W. Chan, Nature 449, 1025 (2007).
- A.S. Rittner and J.D. Reppy, Phys. Rev. Lett. 97, 165301 (2006).
- M. Kondo, S. Takada, Y. Shibayama and K. Shirahama, J. Low Temp. Phys. 148, 695 (2007).
- A. Penzev, Y. Yasuta and M. Kubota, J. Low Temp. Phys. 148, 677 (2007).
- Y. Aoki, J. C. Graves, and H. Kojima, Phys. Rev. Lett. 99, 015301 (2007).
- M. Boninsegni, N.V. Prokof’ev and B.V. Svistunov, Phys. Rev. Lett. 96 105301 (2006); B.K. Clark and D.M. Ceperley, Phys. Rev.. Lett. 96, 105302 (2006).
- M. Boninsegni, N.V. Prokof’ev and B.V. Svistunov, Phys. Rev. Lett. 96, 070601 (2006); Phys. Rev. E 74, 036701 (2006).
- R.A. Aziz, V.P.S. Nain, J.S. Carley, W.L. Taylor and G.T. McConville, J. Chem. Phys. 70, 4330 (1979).
- M. Suzuki, Phys. Lett. A 201, 425 (1995).
- S. Moroni, D.E. Galli, S. Fantoni and L. Reatto, Phys. Rev. B 58, 909 (1998).
- D.E. Galli, L. Reatto, J. Low Temp. Phys. 124, 197 (2001).
- D.E. Galli and L. Reatto, Phys. Rev. Lett. 96, 165301 (2006).
- M. Boninsegni, J. Low Temp. Phys. 141, (2005).
- E.L. Pollock and D.M. Ceperley, Phys. Rev. B 30, 2555 (1984).
- M.C. Gordillo and D.M. Ceperley, Phys. Rev. B 58, 6447 (1998).
- S. Giorgini, J. Boronat and J. Casulleras, Phys. Rev. B 54, 6099 (1996).
- P.A. Whitlock, G.V. Chester and M.H. Kalos, Phys. Rev. B 38, 2418 (1988).
- C. Mora and X. Waintal, Phys. Rev. Lett. 99, 030403 (2007).
- L. Reatto and G.V. Chester, Phys. Rev. 155, 88 (1967).