Zero-temperature equation of state of solid He at low and high pressures
We study the zero-temperature equation of state (EOS) of solid He in the hexagonal closed packet (hcp) phase over the GPa pressure range by means of the Diffusion Monte Carlo (DMC) method and the semi-empirical Aziz pair potential HFD-B(HE). In the low pressure regime (P GPa) we assess excellent agreement with experiments and we give an accurate description of the atomic kinetic energy, Lindemann ratio and Debye temperature over a wide range of molar volumes ( cm/mol). However, on moving to higher pressures our calculated curve presents an increasingly steeper slope which ultimately provides differences within 40 with respect to measurements. In order to account for many-body interactions arising in the crystal with compression which are not reproduced by our model, we perform additional electronic density-functional theory (DFT) calculations for correcting the computed DMC energies in a perturbative way. We explore both generalized gradient and local density approximations (GGA and LDA, respectively) for the electronic exchange-correlation potential. By proceeding in this manner, we show that discrepancies with respect to high pressure data are reduced to with few computational extra cost. Further comparison between our calculated EOSs and ab initio curves deduced for the perfect crystal and corrected for the zero-point motion of the atoms enforces the reliability of our approach.
The physics of helium at low temperatures is among the most intriguing and intensively studied topics in condensed matter science. Despite of being a rare gas element with one of the simplest possible electronic structures, helium constitutes a fundamental system which is challenging for the test and development of methods based on quantum theory. Because of its light atomic mass and weak interatomic interaction, helium is the only system that remains liquid under its own vapor pressure () at zero temperature. Below 2.17 K, liquid He features superfluidity and Bose-Einstein condensation, two striking and inherent quantum effects. With an external pressure of 25 bar, the fluid at crystallizes into the hexagonal closed packet structure (hcp), which remains the stable phase of solid helium at and high pressures, made the exception of an fcc loop along melting in between K. driessen86 (); loubeyre93 (); zha04 ()
Solid helium is manifestly a quantum crystal. In the regime of ultralow temperatures (few mK) this system possesses extraordinarily large atomic kinetic energy ( 24 K) and Lindemann ratio ( 0.26), and likewise anharmonic effects on it are of relevance for predicting and understanding its thermodynamic properties. glyde94 () Further testimony about the uniqueness of this solid is posed by the long-standing controversy sparked by recent experimental findings about whether perfect crystalline He may exhibit superfluid-like behavior and Bose-Einstein condensation (supersolid). kim04 (); rittner06 (); profkofev05 (); ceperley04 (); galli05 () From a technological side, solid helium also has some relevance since it is considered as the best quasi-hydrostatic medium hence modern technologies based on it have emerged and induced considerable progress in the field of high-pressure experiments. takemura01 (); bell81 (); downs96 (); takemura04 ()
In this paper, we study the zero-temperature equation of state of bulk solid He in the hcp phase over a wide pressure range ( GPa) with the Diffusion Monte Carlo (DMC) method and the HFD-B(HE) Aziz pair potential (hereafter referred to AzizII), aziz87 () and additionally with electronic density functional theory (DFT) to account for many-body interactions arising in the system with increasing pressure. In all the work we differentiate between two pressure regimes, namely low pressure ( GPa) and high pressure ( GPa). Quantum Monte Carlo (QMC) methods have proved among the most accurate and reliable tools for solving quantum many-body problems associated to condensed matter systems. hammond94 (); anderson99 (); guardiola98 (); ceperley95 () In particular, the DMC method is a zero-temperature approach which yields exact estimation (only subject to statistical uncertainty) of the ground-state energy and related properties of many-boson interacting systems. kalos66 (); kalos70 (); ceperley79 () During the last few decades this and other Monte Carlo techniques (mainly, the variational Monte Carlo -VMC- and Path Integral Monte Carlo -PIMC- methods) have been fruitfully applied to the study of noble gases and light elements and compounds like He, Ne, H, D, LiH and LiD in homogeneous and inhomogeneous phases and both in bulk and in reduced dimensionalities. sola06 (); boronat04 (); gordillo03 (); grau02 (); drummond06 (); cazorla05 (); cazorla04 () The great capability of the DMC method is related to the existence of accurate interatomic potentials, which are expressed in the form of many-body expansions, and are tuned to reproduce empirical and/or theoretical data. Interatomic potentials are of precious value because provide computational affordability by allowing one to model atoms as interacting points (thus avoiding direct treatment of the electronic degrees of freedom of the system), and also simplified understanding of the system under study. In the case of helium, the semi empirical pair-potential HFD-HE2 proposed by Aziz et al. aziz79 () more than twenty years ago has allowed for quite precise reproduction of the energetic and structural properties of the liquid and solid phases near equilibrium. kalos81 (); pandharipande83 () In this work, we use a newer version of this potential, namely the HFD-B(HE) one, aziz87 () which has demonstrated excellent performance in the description of the liquid boronat94 () but heretofore has not been tested in the crystal upon high pressure.
Anticipating some of the outcomes we are to present, excellent agreement between our results for EOS and experiments is assessed in the low pressure regime for volumes ranging from cm/mol to cm/mol; however, discrepancies start to develop at smaller volumes ( GPa). Within the low pressure regime, we provide exact estimation within some statistical error of the kinetic energy per atom, Lindemann ratio and Debye temperature of the system by means of the pure estimator technique. liu74 (); reynolds86 (); casulleras95 () In the high pressure regime, however, our equation of state systematic and increasingly overestimates the pressure. Discrepancies with respect to measurements amount to at GPa and to at GPa. Previous PIMC work on the EOS of solid He at ambient temperature ( K) and performed with akin model pair potentials arrived to similar disagreements. boninsegni01 (); herrero06 () With the aim of analyzing the possible causes of this large disagreement we first examine the influence of finite size effects in our results. Indeed, finite size effects become larger by increasing pressure because cut-off distances involved in the calculation of the atomic interactions within the system are continously reduced (generally these are chosen as half the length of the simulation box). Accordingly, the radial pair distribution function for crystals, , emerges progressively less smooth with compression. Therefore, customary corrections devised for dealing with finite size effects which are based on simple approximations for , might introduce appreciable deviations in the results (see Fig. 1). Because of these effects, we have regarded as essential to quote the energy tails accounting for the finite size corrections by means of two different approaches: (i) considering beyond a certain cut-off distance and then integrating the simplified analytic expressions for the tails, and (ii), computing the variational Monte Carlo energy of progressively enlarged systems and then estimating the energy of the corresponding infinite system by means of extrapolation to . Certainly, approach (ii) results computationally more demanding than (i) but also more accurate, and we find a pressure difference of GPa between both resulting EOSs at the smallest studied volume ( cm/mol). Nevertheless, this discrepancy by itself does not explain the large disagreement between our results and high-pressure data. In consequence, we turn our main concern to the characterization of the interatomic interactions.
It is well-known that the structural and electronic properties of the atomic and molecular systems may experience important arrangements by effect of pressure. cohen74 (); cohen84 (); christensen06 (); dewaele03 () Upon compression, overlappings between the electronic clouds within the system are promoted hence further correlations among the atoms (angular forces) emerge so as to lower their energies. In the case of solid He, it has been suggested and tested within the Self Consistent Phonon formalism that three-body exchange interactions become significant with increasing density. loubeyre87 () In Ref. boninsegni01, , Chang and Boninsegni included three-body effects into their high-pressure PIMC calculations performed with pair potentials, by computing the energy of several three-body interaction models over sets of configurations generated in their simulations (that is, perturbatively). In doing this, their agreement with experiments did not improve quantitatively, thus they suggested that higher order many-body contributions to the energy had to be considered. More recently, Herrero herrero06 () has adopted a similar but computationally more demanding approach to that of Chang and Boninsegni in which three-body interactions are explicitly included into the model Hamiltonian. For the case of a rescaled Bruch-McGee three-body interaction, Herrero’s work provides very notable agreement with experiments up to pressures of GPa and at room temperature.
In the present work, we propose a perturbative approach for quoting the many-body interactions happening within highly compressed solid He which are not accounted for by any atomic pair potential, and without increasing the computational cost significantly. Essentially, this consists in performing ab initio density functional calculations over sets of atomic configurations independently drawn from DMC simulations; subsequently, the energies previously computed with DMC are corrected according to the average difference between the ab initio interaction energy of the all-electron-ion system and the pair-potential energy. In this way, many-body interactions of order two and higher are included perturbatively into the EOS without requiring from the knowledge of any additional two-, three-, four-, and so on, body interatomic potential. We show that proceeding in this manner the agreement with high pressure experimental data is at the level with relatively few computational extra cost. Truly, the approach that here we present for helium can be extended to the study of other quantum crystals upon high pressure for which accurate pair potentials are available.
The remainder of this paper is as follows. Section II describes the methods that we have employed, the treatment of finite size effects and gives the technical details. In Section III, we present our results for the ground state of solid He at low and high pressure and yield further comparison with first-principles based calculations. Finally, in Section IV we analyze the pros and cons of the proposed perturbative scheme and conclude with the final remarks.
Ii Approach and Methods
ii.1 Diffusion Monte Carlo
DMC is a ground-state method which provides the exact energy within statistical errors of many-boson interacting systems of interest. hammond94 (); guardiola98 (); ceperley79 () This technique is based on a short-time approximation for the Green’s function corresponding to the imaginary time-dependent Schrdinger equation, which is solved up to a certain order of accuracy within an infinitesimal interval . Despite this method is algorithmically simpler than domain Green’s function Monte Carlo, ceperley79 (); kalos74 () it presents some bias coming from the factorization of the imaginary time propagator . Our implementation of DMC is quadratic, chin90 () hence the control of the time-step bias is efficiently controlled since the required extrapolation is nearly eliminated by choosing a sufficiently small time step. The Hamiltonian, , describing our system is
where is the mass of a He atom, the distance between atoms composing an , pair and the HFD-B(HE) Aziz interaction. aziz87 () The corresponding Schrdinger equation in imaginary time (),
with an arbitrary constant, can be formally solved by expanding the solution in the basis set of the energy eigenfunctions . It turns out that tends to the ground state wave function of the system for an infinite imaginary time as well as the expected value of the Hamiltonian tends to the ground-state value . The hermiticity of the Hamiltonian guarantees the equality
where is a convenient trial wave function which depends on the atomic coordinates of the system . Consequently, the ground-state energy of the system can be computed by calculating the integral
where (assuming it is normalized), and is the local energy defined as . The introduction of in is known as importance sampling and it certainly improves the way in which integral (4) is computed (for instance, by imposing when is smaller than the core distance of the interatomic interaction).
with and where and are variational parameters. This model is composed of two-body correlation functions accounting for the two-body correlations induced by , and one-body functions which localize each particle around a site of the equilibrium lattice of the crystal as given by the set of vectors . Initially, the parameters contained in are optimized by means of variational Monte Carlo at some molar volume near equilibrium; however, as we have explored the system over a wide range of volumes we have repeated this procedure at some other selected points along the EOS. For instance, the optimized value of the parameters and at the molar volume cm is and Å, respectively, while at cm results and Å . The parameters of the simulations, namely, the number of particles, , critical population of walkers, , and time step, , have been adjusted to eliminate any residual bias coming from them; their respective values are 180, 400 and 2.710 K. The parameter has been reduced progressively with increasing density in order to provide numerical stability.
ii.2 Finite size corrections
The description of an infinite system of interacting particles is obtained in practice through the simulation of a finite number of particles enclosed within a box. The difference between the scale of the real and simulated systems can be overcome by enlarging the size of the simulated system so much as possible and applying periodic boundary conditions to it. allen89 () Even so, several corrections to the energies quoted directly from the simulation must be done if correlations of longer range are present. Certainly, these corrections arise from the fact that the maximum distance involving correlations in the simulation coincides with the length-scale of the particle container. The expressions for the potential and kinetic energy corrections and , assuming a certain cut-off length R for the computation of the correlations (generally chosen as half the length of the simulation box), are :
where , and are the number, diffusion constant and density of particles, and , and the radial pair distribution function, pair potential and two-body correlation function entering the trial wave function, respectively. In the case of liquids, can be well-approximated to unity in equations (6) and (7), and consequently, and turn out to be analytically accessible (standard tail correction -STC-). Nevertheless, in the case of solids such approximation could result rather inaccurate owing that the pattern of the radial distribution function is still oscillating beyond the cut-off distance (see Fig. 1). In view of these facts and in order for the attained description of solid He to be as precise as possible, we have estimated and also by means of VMC (variational tail correction -VTC-) through the relation
where the superscripts in the energies refer to number of particles, is the number of particles used in the DMC simulations and . The limit , equivalent to in equations (6) and (7), is reached through successive enlargements of the simulation box at fixed density (up to 900 particles) and further linear extrapolation to infinite volume. Indeed, this procedure results computationally affordable within VMC but not within DMC. In Fig. 2, we shown the asymptotic agreement between standard and variational energy tail corrections for infinite solid He () within VMC.
ii.3 Ab initio calculations and perturbative approach
Density Functional Theory (DFT) is a first-principles quantum approach which has allowed for accurate and reliable knowledge of a great deal of materials and systems with exceptional computational affordability. A comprehensive description of DFT methods as applied to the modeling of condensed matter is given in recent books and reviews. martin04 (); kohanoff06 () In DFT, the ab initio free energy of an atomic system, given the positions and charges of its nuclei, is expressed as a functional of the electronic density, , as follows:
where is the electronic kinetic energy, and the atomic number and position of atom , respectively, and the electronic exchange-correlation energy (we have imposed and ). The other terms in Eq.(9) account for the Coulomb interactions between electrons, electrons and nuclei and nuclei. The Hohenberg-Kohn theorem states that the density which minimizes the functional corresponds to the true ground-state density of the system (thus ) and that this optimal solution is unique. It is demonstrated that DFT is an exact electronic ground-state method, whereas the electronic exchange-correlation functional is not known for most of the systems. Consequently, some approximations for it must be introduced in the calculations. The most widely used models for are the local density approximation (LDA) and the generalized gradient approximation (GGA), which have been parameterized by different groups. In this work, we use both Ceperley-Alder (CA) version of LDA ceperley80 () and Perdew-Wang (PW91) of GGA, perdew92 () since a priori one can not discern confidently which is going to result more reliable for the study. A completely independent issue from the choice of is the implementation of DFT that is used. This mainly concerns the way in which the electron orbitals are represented. Here, we use the projector augmented wave (PAW) framework developed by Blch blochl94 () and as implemented in the VASP program. kresse93 (); kresse94 (); kresse96 ()
The perturbative approach that we propose for correcting the DMC energies obtained with the pair potential , consists in averaging the quantity over sets of configurations drawn independently from the DMC simulations. According to this, the corrected energies result
The many-body correction includes two-, three- and so on many-body contributions to the total energy as can be seen by invoking a many-body expansion of the ab initio ground state energy
We note that the family of vectors here refer to the positions of the atoms (nuclei) and not to the sites of the perfect crystalline lattice. It turns out that all the many-body terms composing are evaluated for any arrangament of the atoms as generated according to the Hamiltonian in Eq.(1), and included into the total energy in a perturbative manner. Certainly, our many-body approach is not exact; firstly, it is noted that the full quantum Hamiltonian of the system expressed within the Born-Oppenheimer approximation (FQ-BO) might be written down as , where corresponds to the kinetic energy of the nuclei, and that Eq.(1) results a simplification of it, needless to be said, extremely accurate at low and moderate pressures. Nevertheless, using the DMC method for solving the ground-state of remains a future goal because of the numerous intricacies encountered in the treatment of the electronic degrees of freedom (i.e. choice of the trial wavefunction and sign problem) and large computational cost involved. Therefore, instead of abording the full quantum problem straightaway, we have opted for a simplified but affordable strategy: add and substract into , solve exactly the part of the Hamiltonian embodying most of the two-body interactions and account for the rest by means of first-order perturbation theory.
The ab initio calculations required for the computation of have been performed on supercells containing 180 particles and with Monkhorst-Pack -sampling of the Brillouin zone monkhorst76 () and cut-off energy 478.0 eV; these settings ensure energy convergence to better than 0.5 K per atom. On the other side, our criterion for the convergence of correction relies on the measure of its fluctuation, , over the same collection of DMC configurations than used for the average . Given a molar volume, we have requested to be less than 1 K per atom, that is approximately the total DMC energy obtained at large densities. In the process of drawing atomic configurations from the DMC runs, we have imposed the only constraint , where is the local energy of the considered configuration and the mean energy calculated over the population of walkers to which it corresponds. We have proceeded so for avoiding spurious configurations on the averages which otherwise are rejected within few steps in the DMC sampling. The number of atomic configurations required for the convergence of has proved smaller than initially expected in all the studied cases: about 15-25 were enough. This rapid convergence of the fluctuations reveals that despite two-body interactions by themselves are not sufficient to attain reliable description of very dense solid He they are still of leading relevance on it.
iii.1 Low pressure regime
The EOSs of solid He has been obtained by fitting a fourth-order polynomial to the DMC energies and subsequently performing the derivation respect to volume,
where is directly the equilibrium volume of the system and , constants. In Fig. 3, we compare the DMC energies at volumes close to equilibrium (), obtained with both STC and VTC, with experimental data. driessen86 () As one observes there, excellent agreement between measurements and VTC results is provided, however, differences with respect to STC results are quoted in an almost constant upwards shift of K at positive pressure. Despite these discrepancies will practically vanish when expliciting both VTC and STC EOSs because of the energy derivative involved (as it will be shown shortly), it is noted that for other magnitudes which explicitly depend on the internal energy, as for example the enthalpy or freezing and melting densities, VTC and STC lead to different results. In the same figure, we also display previous theoretical calculations performed with Green’s function Monte Carlo (GFMC) and the HFD-HE2 Aziz pair potential. kalos81 () The GFMC points perfectly coincide with our results obtained with VTC, however, they disagree with the STC ones in the same manner experimental points do. Since GFMC and DMC are exact ground-state methods, energy differences between both approaches should be due only to the model interaction. Assuming that the treatment of finite size effects adopted in Ref. kalos81, corresponds to the commonly used STC one, we may just conclude that both HFD-HE2 and HFD-B(HE) interatomic potentials are likely to produce equivalent curves at low pressures (likewise VTC and STC lead to practically identical EOSs) but not so total energies and other directly related properties.
Fig. 4 reports our results of the EOS of solid He at in the volume range cm/mol ( GPa). Curves obtained with VTC and STC are coincident as we anticipated in the previous paragraph. The parameters of the fit (12) also displayed in Fig. 4 (we provide the one obtained with VTC) are K, K and cm/mol (uncertainties are shown within parentheses). A glance at the plot reveals an excellent agreement between our results and experiments at low pressures, however, discrepancies become progressively larger as we move towards volumes smaller than cm/mol ( GPa). (For instance, at cm/mol our prediction of pressure overestimates the experimental value within .) It is worth noticing that the worsening of our curve roughly coincides with the interval in which the potential energy of the system becomes positive (see Table I). This fact indicates that the repulsive part of the HFD-B(HE) potential is probably too stiff. In the next subsection we will extensively deal with the shortcomings derived from the adopted model interaction, however now we continue with the description of other atomic magnitudes of interest that we have obtained at low pressures.
The zero-temperature atomic kinetic energy of solid He, , is an important (and challenging) quantity to measure and compute since it evidences the singular quantum nature of this crystal. It is well-known that the zero-point energy of solid helium is comparable in magnitude to its potential energy (cohesive energy), , and that the ratio between these two energies gives a qualitative idea about the relevance of anharmonic effects in the system (the larger is, the larger anharmonic contributions would result). glyde94 (); moleko85 () From a computational side, exact estimation of the expected ground state values of operators which do not commute with the Hamiltonian, as for instance the potential and kinetic energy operators, may be provided within the DMC scenario by means of the pure estimator technique. liu74 (); reynolds86 (); casulleras95 () In practice, this technique involves the introduction of additional weight factors into the customary DMC sampling which retain memory of the configurational replication processes occurring along the simulation. In order for our evaluation of the zero-temperature kinetic energy of solid He to be as reliable as possible, we first determine the exact potential energy of the system by means of the pure estimator method and then we subtract it to the total energy (we note that within DMC the estimator of the kinetic energy is slightly biased by the choice of the trial wavefunction). In Fig. 5 we display our results for and compare them to low temperature data provided by several authors: hilleke84 (); diallo04 (); adams07 (); celli98 () the overall agreement between them is excellent. In particular, we note the perfect agreement of our calculated value K at cm/mol with the very recent neutron scattering measurement of Diallo et al., diallo04 () K, performed at the same volume. In Table I, we enclose DMC results for the total, kinetic and potential energies of solid He including STC at some selected volumes within the interval - cm/mol. In Fig. 5, however, we have not refrained from including a further point at volume cm/mol for which we have shown that the description of the system attained with the model interaction seems to be not fully reliable. It should be mentioned that the treatment of finite size effects adopted in the calculations has little effect on the final values of since the largest contribution to the total energy tail correction stems from the interatomic interactions.
Another quantity of interest in the study of quantum solids is the atomic mean squared displacement, , which is directly measured in x-ray diffraction experiments. In connection to this magnitude, the Lindemann ratio is defined as , where is the distance between nearest neighbors in the perfect crystalline lattice. As we pointed out in the Introduction, the zero-temperature Lindemann ratio of solid He is uncommonly large (even if compared to other distinguished quantum solids like for example H which possesses ) as consequence of its light atomic mass and weak interatomic interaction. Using the pure estimator technique, we have studied the dependence of with volume over the range - cm/mol. We have depicted our results for in Fig. 6 and compared them to experimental data of different authors, and again the overall agreement between them is remarkable. Once is known, the Debye temperature of the system at , , is deduced straightforwardly through the relation . We have fitted our results for with the relation
where and which has been used previously to reproduce the density dependence of the phonon frequencies in solid H and He as well. driessen84 (); driessen86 () Our optimal coefficients for expression (13) plotted in Fig. 6 are: cm/mol, , , and . (Additional point at cm/mol in Fig. 6, has not been used in the fit.)
iii.2 High pressure regime
As we have already illustrated in Sec.III.1, the pair potential HFD-B(HE) performs excellently in the description of solid He up to volumes of cm/mol, however, it monotonically fails to reproduce its EOS as the density is increased beyond this point. In Fig. 7, we explicit the differences between measurements of Ref. driessen86, and loubeyre93, and our calculations performed with STC and VTC for finite size effects and over the pressure range GPa. The pressure difference between the EOSs obtained with VTC and STC at the highest studied density amounts to GPa, however this quantity results very small when compared to the discrepancy of both with experiments which is of the experimental value. This discrepancy is in overall agreement with the microscopic calculations of Refs. herrero06, and boninsegni01, .
Our results and other found in the literature, boninsegni01 (); herrero06 () pose the need for considering higher order many-body effects present on dense He instead of going in the search of improved pair potential models. As we pointed out in the Introduction, many authors have made efforts for elucidating the relevance of three-body and higher order (up to six-body) effects on the EOS of solid He with assorted degree of accuracy and success. boninsegni01 (); herrero06 (); tian06 () In this section, we present the curves calculated within our proposed scheme for correcting the DMC energies obtained with pair potentials as described in Sec. II.3. Just as we have explained there, all the many-body interactions not accounted for by are computed with ab initio DFT and included into the total energy in a perturbative way, without requiring from the knowledge of additional two-, three- and/or higher order many-body interaction models.
The fits to our results displayed in the figures of this and next subsections have been performed with the Vinet relation vinet87 ()
where , and are the equilibrium volume, equilibrium isothermal bulk modulus () and equilibrium , respectively. The experimental values of these parameters as provided by Ref. loubeyre93, are cm/mol, GPa and . We have also enclosed data points of Ref. driessen86, in the plots for additional comparison with our estimations. The improvement of our EOS when considering many-body effects computed with the proposed perturbative approach is substantial (see Fig. 8). For example, within LDA we obtain GPa at volume cm/mol which is quite close to the experimental value GPa and far below the non-corrected DMC result of GPa. The parameters of the fit corresponding to this case are cm/mol, GPa and , which as it is observed in Fig. 8 leads to constant underestimation of pressure within few Gpa respect to experimental data along the whole depicted range. On the contrary, the EOS obtained with GGA provides a notable description of the system near the equilibrium and few GPa above the experimental values at high pressures. Putting this into figures, cm/mol, GPa and , which in fact result closer to the experimental values of Ref. loubeyre93, than the LDA ones. It is worth mentioning that the observed tendency of GGA (LDA) for overestimating (underestimating) pressure in our results is a well-known outcome in the field of ab initio simulations.
Table II yields the values of the DMC energy and perturbative corrections for solid He at some selected volumes. Separately, we have shifted all the LDA and GGA corrections a same amount for providing null contributions at the largest enclosed volume so as to facilitate the comparison between them. Certainly, this can be done without any loss of generality since the zero of the LDA and GGA ab initio energies and that of the HFD-H(BE) interaction do not coincide and we are essentially interested in the pressure. Two main conclusions can be extracted from the values in Table II: (i) corrections performed with LDA decrease monotonically with compression, not so with GGA, and (ii) GGA corrections are smaller in absolute value than the LDA ones. Since the proposed approach for correcting the DMC energies obtained with pair potential is perturbative, the conclusion (ii) concedes more reliability to the results obtained with GGA than with LDA. Indeed, a conclusive answer about whether LDA figures are or not too large would be best provided by second order perturbative theory, however this is out from our scope. In the next subsection, we will comment again on this issue by supplying further comparison between results presented here and others obtained by means of ab initio procedures.
One known weakness of DFT calculations is that the usual approximations for may fail to capture the essence of the long-range forces present in the system. kohn98 (); basanta05 () In the case of rare gases, the van der Waals (vdW) energy, which physically accounts for Coulomb correlations between distant electrons, has notorious relevance in the cohesion of the system. With the aim of estimating the effect of this shortcoming in our corrections, we have added an effective two-body term accounting for the vdW interactions to the ab initio energy . This term is expressed as
where KÅ and for but for with Å (that is, as given by the HFD-B(HE) interaction), and it has been evaluated over the same sets of atomic configurations than used for the computation of . Following this receipt, the many-body correction devised for energies now can be redefined as . In Fig. 9, we plot the curves obtained with the correction , and Table III summarizes the parameters of all the fits that we have performed (for LDA and GGA corrections including and not including vdW contributions). On one hand, a glance at the figure reveals that considering vdW interactions as explained above has in general little effect on the results, just a slight and otherwise expected lowering of the curves within few GPa over the whole depicted range. On the other hand, the equilibrium properties of the system, as given by the parameters of the fits, change appreciably (see figures enclosed in Table III). This result seems to corroborate the accepted assumption that the effect of long-range interactions in the EOS of rare-gas solids becomes less important with increasing pressure. kim06 (); dewhurst02 (); iitaka02 (); kwon95 ()
In this subsection, we have not attempted to enclose any result for fcc He in the plots and/or tables since experiments indicate that hcp is the only stable phase of solid helium at high pressures ( GPa) and low temperatures, apart from a small fcc loop region around melting between and K. driessen86 (); loubeyre93 () Reassuringly, previous work based on first-principles calculations agrees to regard hcp as the most energetically favourable zero-temperature phase of He upon pressures up to GPa. nabi05 () In spite of this, we have carried out a series of calculations in highly compressed fcc He in order to check the predictability of our approach. Essentially, our results show no appreciable energy differences between the two phases within the numerical uncertainty. This outcome, however, appears to be not surprising since short-range interactions in helium are of leading importance, and the first and second shells of nearest neighbours in the fcc and hcp phases peak at practically indentical distances given a same density.
iii.3 Comparison with ab initio based calculations
Within the DFT formalism, the zero-temperature energy of a solid is usually written as a sum of two different contributions
where is the energy of the perfect crystal (atoms frozen on their sites) and accounts for the motion of the atoms and is expressed as a sum of harmonic and anharmonic terms. In practice, is obtained with standard DFT calculations and it involves affordable computations performed within one unit cell of the perfect crystal (apart from the summations involving periodic boundary conditions). On the other side, the estimation of requires from some knowledge of the phonon-related properties of the solid of interest. In the case of heavy-ion crystals, the quasiharmonic approximation in combination with finite displacement methods have allowed for an accurate description of the phonon frequency spectra. alfe01 (); vocadlo02 () The basic strategy underlying these methods consists in distorting the perfect crystal by displacing certain selected atoms slightly from their equilibrium positions and then evaluating the atomic forces arising on the system by means of the Hellman-Feynman theorem and DFT. This approach, however, fails to reproduce solid He since it provides negative (imaginary) phonon frequencies associated to its experimental stable phases at intermediate pressures. cazorla07 () Truly, the relevance of anharmonic effects in solid He makes the computation of its vibrational properties a tedious and complicate task which requires from approaches going beyond the harmonic and/or quasiharmonic approximations. It should be noted that within DMC this difficulty is circumvented since the phononic nature of the studied system is inherently cast into the method, hence further partition of the energy into static and vibrational parts is not required.
In order to contrast our results presented in Sec. III.2 with other obtained with ab initio based methods, we have computed the EOS of solid He through the relation (16) by evaluating with DFT and with the Mie-Gruneisen model hemley90 ()
where is the Debye temperature, the Gruneisen parameter which we approximate as (with ) and the gas constant. Indeed, we have used for the experimental relation provided by Driessen et al. driessen86 () since, as we have noted previously, the estimation of this or any other vibrational property of solid He at would result not trusty with customary ab initio strategies used in the study of normal (not quantum) solids. The resulting increases monotonically with compression and for instance it represents about of the total pressure of the system at volume cm/mol, thus it must not be neglected. In Fig. 9, we show the EOSs obtained with the already explained procedure using both LDA and GGA approximations for the exchange-correlation energy; also we include the curves quoted within DMC and corrected for the vdW energy and many-body interactions, and experimental data. As it is observed there, differences between both LDA and GGA perturbationally corrected curves and their ab initio counterparts are quite small. Moreover, these discrepancies are likely caused by the treatment of the long range interactions (vdW energy) described in Sec. III.2 and the approximation adopted for . This result is stimulating since it demonstrates that with the approach presented in Sec. II.3 one can obtain accurate EOSs for dense solid helium, or equivalently for any other light quantum solid, in excellent agreement with those which would be obtained by means of first-principles approaches but with the benefit of not requiring from the computation of the phonon dispersion curves of the crystal or experimental data.
Now we turn our attention to the concern posed over LDA in the previous subsection. As we noted there, a glance at Table II might lead to the conclusion that the LDA corrections are too large so as to be considered perturbative on top of the DMC energies (not so the GGA ones). Indeed, we do not dispose of a fair criterion for accepting or rejecting corrections in basis to their size, and in our opinion this is the most important shortcoming of our approach. Nevertheless, appealing to the good accordance between the LDA curve corrected for the atomic zero-point motion and the DMC one corrected with LDA, we may feel quite confident about the reliability of the latter.
Iv Discussion and Conclusions
In the Introduction we pointed out that Diffusion Monte Carlo is among the best suited methods for studying quantum solids. In the case of bosonic systems, this method provides the exact ground state energy and related properties without dependence on the choice of the trial wave function, which otherwise is related to the computational efficiency. Here, we have proved the excellent performance of DMC with the Aziz pair potential HFD-B(HE) in characterizing solid He at low pressures ( 0.65 GPa), by estimating its EOS, kinetic energy per particle and Debye temperature at different volumes, and comparing our results with experiments. Especial attention has been paid to the extend of finite size effects in our results. To this regard, we conclude that customary strategies devised to correct such effects based on the approximation of the pair radial distribution function to unity beyond certain cut-off distance, may result accurate enough for the derivation of the EOS but not so for the assessment of other magnitudes like the energy.
On the other side, solid helium under high compression (as most of the materials) undergoes important arrangements on electronic structure which lead to the appearance of angular correlations among the atoms. dewaele03 () This circumstance makes necessary to consider not only atomic pair interactions but also higher order many-body ones when investigating this crystal upon high pressure. Nevertheless, within DMC the minimal inclusion of three-body interactions on the model Hamiltonian has already the effect of drastically slow down the simulations. Furthermore, even in the supposed case computational cost was not a problem, first we should know much better than now the analytical form of these many-body interactions or alternatively to be able to devise them (which actually may result puzzling). According to this occurrence, ab initio methods emerge among the best candidates for quoting such contributions since they do not rely on potential models and, in general, are computationally affordable. However, fully ab initio analysis of crystals requires from the knowledge of the phonon-related properties, and for the case of solid He and other light quantum solids this is by no means a straightforward task.
In this work, we have presented an approach for the study of dense solid He at zero temperature which combines the versatility of the DMC method with the accuracy of ab initio calculations. On one hand, we naturally circumvent the calculation of the vibrational properties of helium thanks to the DMC strategy, and on the other, we account for the many-body interactions having place on the system by means of DFT. However, the way by which we enclose these many-body contributions to the DMC energy is not exact but perturbative and we do not dispose of rigorous tests for quoting the errors included on these corrections. Concerning results, we have yielded the EOSs corrected in this manner within LDA and GGA for the exchange-correlation energy and they have proved in fairly notable agreement with experiments over the pressure range GPa. Specifically, GGA provides a better description of the crystal near equilibrium than LDA. Further comparison of these curves with EOSs obtained trough DFT and corrected for the atomic zero-point motion by means of an approximate model, comes to support the reliability of our approach.
It should be mentioned that the zero-temperature scheme proposed in this work is also well suited for the study of other light quantum solids upon high pressure, like He, H, D and Li, for which accurate pair potentials are devised. Certainly, a further and promising improvement of the present framework would consist in going beyond the perturbative approach. This could be achieved by proper coupling of the DMC and DFT methods, as for instance, by considering the ab initio energy of the system within the branching weight of the DMC algorithm. Work on this direction is in progress.
Acknowledgements.We acknowledge financial support from DGI (Spain) Grant No. FIS2005-04181 and Generalitat de Catalunya Grant No. 2005GR-00779.
- (1) A. Driessen, E. van der Poll and I. F. Silvera, Phys. Rev. B 33, 3269 (1986).
- (2) P. Loubeyre, R. LeToullec, J. P. Pinceaux, H. K. Mao, J. Hu and R. J. Hemley, Phys. Rev. Lett. 71, 2272 (1993).
- (3) C-S. Zha, H-K. Mao and R. J. Hemley, Phys. Rev. B 70, 174107 (2004).
- (4) H. R. Glyde in Excitations in Liquid and Solid Helium (Clarendon Press, Oxford, 1994).
- (5) E. Kim and M. H. W. Chan, Science 305, 1941 (2004).
- (6) A. S. C. Rittner and J. D. Reppy, Phys. Rev. Lett. 97, 165301 (2006).
- (7) N. Profkof’ev and B. Svistunov, Phys. Rev. Lett. 94, 155302 (2005).
- (8) D. M. Ceperley and B. Bernu, Phys. Rev. Lett. 93, 155303 (2004).
- (9) D. Galli, M. Rossi and L. Reatto, Phys. Rev. B 71, 140506 (2005).
- (10) K. Takemura, J. Appl. Phys. 89, 662 (2001).
- (11) P. M. Bell and H. K. Mao, Year Book - Carnegie Inst. Washington 80, 404 (1981).
- (12) R. T. Downs, C. S. Zha, T. S. Duffy and L. W. Finger, Am. Mineral. 81, 51 (1996).
- (13) K. Takemura, K. Sato, H. Fujihisa and M. Onoda, Ferroelectrics 305, 103 (2004).
- (14) R. A. Aziz, F. R. W. McCourt and C. C. K. Wong, Mol. Phys. 61, 1487 (1987).
- (15) B. L. Hammond, W. A. Lester Jr. and P. J. Reynolds in Monte Carlo Methods in Ab initio Quantum Chemistry (World Scientific, 1994).
- (16) J. B. Anderson, Rev. Comp. Chem. 13, 133 (1999).
- (17) R. Guardiola in Microscopic Quantum Many-Body Theories and Their Applications ed. by J. Navarro and A. Polls (Springer, Berlin, 1998).
- (18) D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995).
- (19) M. H. Kalos, J. Comp. Phys. 1, 257 (1966).
- (20) M. H. Kalos, Phys. Rev. A 2, 250 (1970).
- (21) D. M. Ceperley and M. H. Kalos in Monte Carlo Methods in Statistical Physics (Springer, Berlin, 1979).
- (22) E. Sola, J. Casulleras and J. Boronat, Phys. Rev. B 73, 092515 (2006).
- (23) J. Boronat, C. Cazorla, D. Colognesi and M. Zoppi, Phys. Rev. B 69, 174302 (2004).
- (24) M. C. Gordillo, J. Boronat and J. Casulleras, Phys. Rev. B 68, 125421 (2003).
- (25) V. Grau, J. Boronat and J. Casulleras, Phys. Rev. Lett. 89, 045301 (2002).
- (26) N. D. Drummond and R. Needs, Phys. Rev. B 73, 024107 (2006).
- (27) C. Cazorla and J. Boronat, J. Low Temp. Phys. 139, 645 (2005).
- (28) C. Cazorla and J. Boronat, J. Low Temp. Phys. 134, 43 (2004).
- (29) R. A. Aziz, V. P. S. Nain, J. S. Carley, W. L. Taylor and G. T. Mc Conville, J. Chem. Phys. 70, 4330 (1979).
- (30) M. H. Kalos, M. A. Lee, P. A. Whitlock and G. V. Chester, Phys. Rev. B 24, 115 (1981).
- (31) V. R. Pandharipande, J. G. Zabolitsky, S. C. Pieper, R. B. Wiringa and V. Helmbrehct, Phys. Rev. Lett. 50, 1676 (1983).
- (32) J. Boronat and J. Casulleras, Phys. Rev. B 49, 8920 (1994).
- (33) K. S. Liu, M. H. Kalos and G. V. Chester, Phys. Rev. A 10, 303 (1974).
- (34) P. J. Reynolds, R. N. Barnett, B. L. Hammond and W. A. Lester Jr., J. Stat. Phys. 43, 1017 (1986).
- (35) J. Casulleras and J. Boronat, Phys. Rev. B 52, 3654 (1995).
- (36) S. Y. Chang and M. Boninsegni, J. Chem. Phys. 115, 2629 (2001).
- (37) C. P Herrero, J. Phys.: Condens. Matter 18, 3469 (2006).
- (38) S. G. Louie and M. L. Cohen, Phys. Rev. B 10, 3237 (1974).
- (39) K. J. Chang and M. L. Cohen, Phys. Rev. B 30, 004774 (1984).
- (40) N. E. Christensen and D. L. Novikov, Phys. Rev. B 73, 224508 (2006).
- (41) A. Dewaele, J. H. Eggert, P. Loubeyre and R. Le Toullec, Phys. Rev. B 67, 094112 (2003).
- (42) P. Loubeyre, Phys. Rev. Lett. 58, 1857 (1987).
- (43) M. H. Kalos, D. Levesque and L. Verlet, Phys. Rev. A 9, 2178 (1974).
- (44) S. A. Chin, Phys. Rev. A 42, 6991 (1990).
- (45) L. H. Nosanow, Phys. Rev. Lett. 13, 270 (1964).
- (46) J. P. Hansen and D. Levesque, Phys. Rev. 165, 293 (1968).
- (47) J. P. Hansen, Phys. Letters 30A, 214 (1969).
- (48) M. P. Allen and D. J. Tildesley in Computer Simulation of Liquids (Oxford University Press, 1989).
- (49) R. M. Martin in Electronic Structure (Cambridge University Press, 2004).
- (50) J. Kohanoff in Electronic Structure Calculations for Solids and Molecules: Theory and Computational Methods (Cambridge University Press, 2006).
- (51) D. M. Ceperley and B. I. Alder, Phys. Rev. Lett. 45, 566 (1980).
- (52) J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
- (53) P. E. Blchl, Phys. Rev. B 50, 17953 (1994).
- (54) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
- (55) G. Kresse and J. Hafner, Phys. Rev. B 49, 14251 (1994).
- (56) G. Kresse and J. Furthmller, Comput. Mat. Sci. 6, 15 (1996).
- (57) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
- (58) L. K. Moleko and H. R. Glyde, Phys. Rev. Lett. 54, 901 (1985).
- (59) R. O. Hilleke, P. Chaddah and R. O. Simmons, Phys. Rev. Lett. 52, 847 (1984).
- (60) S. O. Diallo, J. V. Pearce, R. T. Azuah and H. K. Glyde, Phys. Rev. Lett. 93, 075301 (2004).
- (61) M. A. Adams, J. Mayers, O. Kirichek and R. B. E. Down, Phys. Rev. Lett. 98, 085301 (2007).
- (62) M. Celli, M. Zoppi and J. Mayers, Phys. Rev. B 58, 242 (1998).
- (63) C. Stassis, D. Khatarnian and G. R. Kline, Solid State Comms. 25, 531 (1978).
- (64) C. T. Venkataraman and R. O. Simmons, Phys. Rev. B 68, 224303 (2003).
- (65) C. A. Burns and E. D. Isaacs, Phys. Rev. B 55, 5767 (1997).
- (66) A. Driessen and I. F. Silvera, J. Low Temp. Phys. 54, 361 (1984).
- (67) C. Tian, F. Liu, F. Jing and L. Cai, J. Phys.: Codens. Matt. 18, 8103 (2006).
- (68) P. Vinet, J. R. Smith, J. Ferrante and J. H. Rose, Phys. Rev. B 35, 1945 (1987).
- (69) W. Kohn, Y. Meir and D. E. Makarov, Phys. Rev. Lett. 80, 4153 (1998).
- (70) M. A. Basanta, Y. J. Dappe, J. Ortega and F. Flores, Europhys. Lett. 70, 355 (2005).
- (71) E. Kim, M. Nicol, H. Cynn and C. Yoo, Phys. Rev. Lett. 96, 035504 (2006).
- (72) J. K. Dewhurst, R. Ahuja, S. Li and B. Johansson, Phys. Rev. Lett. 88, 075504 (2002).
- (73) T. Iitaka and T. Ebisuzaki, Phys. Rev. B 65, 012103 (2002).
- (74) I. Kwon, L. A. Collins, J. D. Kress and N. Troullier, Phys. Rev. B 52, 16165 (1995).
- (75) Z. Nabi, L. Vitos, B. Johansson and R. Ahuja, Phys. Rev. B 72, 172102 (2005).
- (76) D. Alfè, G. D. Price and M. J. Gillan, Phys. Rev. B 64, 045123 (2001).
- (77) L. Vočadlo and D. Alfè, Phys. Rev. B 65, 214105 (2002).
- (78) As it has been checked by the authors.
- (79) R. J. Hemley, H. K. Mao, L. W. Finger, A. P. Jephcoat, R. M. Hazen and C. S. Zha, Phys. Rev. B 42, 6458 (1990).