Elastic constants of hcp 4 He: Path-integral Monte Carlo results versus experiment
Luis Aldemar Peña Ardila, Silvio A. Vitiello, Maurice de Koning
EElastic constants of hcp He: Path-integral Monte Carlo results versus experiment
L. A. Pe˜na Ardila, Silvio A. Vitiello and Maurice de Koning ∗ Instituto de F´ısica Gleb Wataghin, Universidade Estadual de Campinas,UNICAMP, 13083-859, Campinas, S˜ao Paulo, Brazil
The elastic constants of hcp He are computed using the path-integral Monte Carlo (PIMC)method. The stiffness coefficients are obtained by imposing different distortions to a periodic cellcontaining 180 atoms, followed by measurement of the elements of the corresponding stress tensor.For this purpose an appropriate path-integral expression for the stress tensor observable is derivedand implemented into the PIMC++ package. In addition to allowing the determination of the elasticstiffness constants, this development also opens the way to an explicit atomistic determination ofthe Peierls stress for dislocation motion using the PIMC technique. A comparison of the resultsto available experimental data shows an overall good agreement of the density dependence of theelastic constants, with the single exception of C . Additional calculations for the bcc phase, on theother hand, show good agreement for all elastic constants. PACS numbers: 67.80.B, 02.70.Ss, 62.20.D
I. INTRODUCTION
Despite the intensive research efforts developed overthe past six years, the remarkable results of the torsional-oscillator (TO) experiments on solid He by Kim andChan [1, 2] still elude a consistent explanation. Althoughit has been suggested that the observed nonclassical rota-tional inertia (NCRI) is a manifestation of superfluidityin the solid phase, [3–5] several questions remain unan-swered. One of them is the apparent correlation betweenthe NCRI data and the observation of elastic stiffeningupon cooling, both of which display similar temperatureand He-concentration dependences [6, 7]. The observedstiffening effect has been interpreted in terms of dislo-cation pinning due to He impurities [6, 8] and this hascontributed to ideas that the NCRI may not have an ex-clusive non-superfluid origin but rather that it might alsobe a manifestation of mechanical behaviour [8–12].A major difficulty is the fact that the insight obtainedfrom recent experiments [6, 8, 13] relies on the indirectinterpretation of data, a process that is inevitably basedon assumptions. The dislocation-pinning interpretation,for instance, conjectures the presence of dislocation net-works that are pinned by He impurities binding to thedislocation cores [6, 8]. Despite the elevated degree ofsophistication of recent experiments [6, 8, 13], however,it has not yet been possible to verify these assumptionsexplicitly, hampering a conclusive understanding of theobserved phenomenology and its relation to the NCRIdata obtained in TO experiments.In this context, theoretical modelling can serve as auseful complementary approach. In principle, it allows asystematic study of a hierarchy of well-controlled struc-tures not accessible to experiment, starting at the defect-free crystal and progressing through a series of configura-tions containing different defect geometries. Analysis ofthe results obtained in these different situations may thenassist in the interpretation of existing experimental data or even lead to new predictions. However, a first stepin this effort is to gauge the results of these modellingapproaches for quantities that can be directly comparedto experimental data. In the context of the recent ob-servations of the mechanical behaviour of solid He, itsintrinsic elastic properties are an evident target for suchan assessment.The purpose of the present paper is to compare theresults of state-of-the-art path-integral Monte Carlo [14,15] (PIMC) calculations based on the Aziz pair poten-tial [16] to the available experimental data for the elasticstiffness constants from the 1970s. The PIMC methodhas shown to give excellent agreement with experimentfor several properties of the liquid and solid states, in-cluding the energy, pair correlation functions, structurefactors, and superfluid density. Triggered by TO exper-iments of Kim and Chan, the PIMC methodology hasalso been applied extensively to the study of He super-solidity, focusing on a variety of properties, including thepossibility of superflow induced by lattice defects such asvacancies [17, 18] and dislocations [19, 20].The first path-integral approach to the computation ofelastic constants was based on a direct measurement ofthe elastic constants in terms of the second derivativesof the partition function with respect to strain compo-nents [21]. Here we adopt a different approach based onthe development of an expression for the stress-tensor ob-servable in the path integral formalism. Not only doesthis observable allow the determination of the elastic con-stants of the defect-free crystal, it is also a key observablein the characterization of plastic deformation in terms oflattice defect properties. In this light, this developmentis of interest in its own right, allowing for instance, anexplicit atomistic determination of the Peierls stress fordislocation motion [22]. Here we utilize it to computethe complete set of stress-strain relations by measuringthe elastic stress response to small homogeneous defor-mations, giving the five independent elastic stiffness con- a r X i v : . [ c ond - m a t . o t h e r] S e p stants of the hcp He crystal. In addition to a direct com-parison of elastic constants as a function of density, weanalyze the degree to which the Cauchy relation, whichmeasures to what extent non central forces and zero-pointeffects are important, is satisfied.The remainder of the paper has been organized asfollows. In Sec. II we derive an expression for thestress-tensor observable within the path-integral formal-ism based on the pair-action approximation and whichhas been implemented in the PIMC++ code [18]. SectionIII describes the employed computational setup and sum-marizes the parameters used in the simulations. Next, wedescribe and discuss the obtained results in Sec. IV andsummarize in Sec. V.
II. PATH-INTEGRAL EXPRESSION FOR THESTRESS-TENSOR OBSERVABLE.
Following the usual approach of Parrinello and Rah-man [23–25], we describe homogeneous distortions of aperiodic simulation cell in terms of the box matrix h,whose columns are the three periodic repeat vectors a, b,and c of the computational cell [25], h = a x b x c x a y b y c y a z b z c z (1)In terms of the h matrix, an absolute position r i withinin the cell is written as r i = h s i, (2)where s i is a relative coordinate vector whose com-ponents have values ranging between 0 and 1. In thisdescription, homogeneous deformations are described interms of variations of the c matrix at fixed relative coor-dinates s i .To obtain an expression for the stress tensor in thepath integral formalism we start with its thermodynamicdefinition, describing its components σ ij in terms ofderivatives of the appropriate thermodynamic potential.Specifically, we have [25] σ ij = − h ) (cid:88) k =1 h jk (cid:18) ∂F∂ h ik (cid:19) N,T , (3)where F = F ( N, h , T ) (4)is the Helmholtz free energy of a system containing N particles that are confined to a volume described by thematrix h and in thermal equilibrium with a heat bath attemperature T . As usual, the microscopic description for the stress tensor is obtained through the connection be-tween the thermodynamic potential and the correspond-ing partition function. Here, we have F ( N, h , T ) = − β ln Z ( N, h , T ) (5)where Z ( N, h , T ) is the canonical partition functionand β = ( K B T ) − . In this manner, the components ofthe stress tensor P become P ij = − h ) βZ (cid:88) k =1 h jk (cid:18) ∂Z∂ h ik (cid:19) N,T , (6)In the path-integral formalism, the canonical partitionfunction for a system of distinguishable particles can bewritten as [14] Z ( N, h , T ) = (cid:90) · · · (cid:90) dR · · · dR M − × exp (cid:2) − S path ( R , ··· , R M − , R ) (cid:3) (7)where R k = { r ,k , · · · , r N,k } represents the set of po-sition vectors of the N particles in the k th time sliceof theM-bead closed path R = { R , ··· , R M − , R } , and S path ( R , ··· , R M − , R ) is the path action. Proper sym-metrization for the case of indistinguishable particles isstraightforward [14]. Describing the position vectors interms of the h matrix, the expression becomes Z ( N, h , T ) = (cid:90) · · · (cid:90) [det( h )] NM dS · · · dS M − × exp (cid:2) − S path ( S , ··· , S M − , S ; h ) (cid:3) , (8)where S k = { s ,k , · · · , s N,k } represents the set of scaledposition vectors of the N particles in the k − th time slice.The derivatives of the partition function with respect tothe elements of the h matrix are then given by (cid:18) ∂Z∂ h ij (cid:19) N,T = [det( h )] NM (cid:90) · · · (cid:90) dS · · · dS M − × exp (cid:2) − S path (cid:3) (cid:20) N M ( H ) − ji − ∂S path ∂ h ij (cid:21) . (9)Substitution into Eq. (6) then gives P ij = 1det( h ) β (cid:34) N M δ ij − (cid:88) k =1 h jk (cid:28) ∂S path ∂ h ik (cid:29)(cid:35) (10)where the angular brackets indicate averaging overclosed paths. Given that the paths consist of M links, theabove expression can also be written in terms of averagesover link actions, namely, P ij = − M det( h ) β (cid:34) N δ ij − (cid:88) k =1 h jk (cid:28) ∂S Link ∂ h ik (cid:29)(cid:35) (11)The hydrostatic pressure P is then given by P = 13 Tr P = − M det( h ) β (cid:20) N − h ik (cid:28) ∂S Link ∂ h ik (cid:29)(cid:21) . (12)Next we determine the derivative of the link actionwith respect to the elements of the h matrix. Here, weare specifically interested in the pair approximation forthe action, S link = S ( R, R (cid:48) ; τ ) = S kin ( R, R (cid:48) ; τ ) + S pot ( R, R (cid:48) ; τ )= N (cid:88) n =1 K ( r n , r (cid:48) n ; τ ) + (cid:88) n 12 ( | r | + | r (cid:48) | ) = 12 ( | h s | + | h s (cid:48) | ) , (20) z ≡ | r | − | r (cid:48) | = | h s | − | h s (cid:48) | , (21) t ≡ | r − r (cid:48) | = | h ( r − r (cid:48) ) | , (22)the potential action then takes the form S pot = N (cid:88) n All calculations are based on an orthorhombic com-putational cell containing 180 He atoms arranged on ahcp lattice characterized by the ideal ratio c/a = 1 . . 3% of the experimental estimates for allconsidered densities [27]. The cell is constructed in sucha manner that the x , y , and z axes are parallel to thecrystallographic [1 2 1 0], [1 0 1 0], and [0 0 0 1] di-rections, respectively. Standard periodic boundary con-ditions are applied throughout. The simulations havebeen carried out using the P IM C ++ package [18], whichis a C + + implementation of the P IM C algorithms de-scribed in Ref. [14]. The used pair action was obtainedfrom a standard matrix squaring procedure [14, 15] us-ing the Aziz HFD − B3 − FCI1 pair potential [16] andan interaction cutoff of 8˚ A . All PIMC simulations em-ploy a time step τ = β/m = 1 / K − [28], with β theinverse temperature and M the number of beads in thering polymers. The determination of the elastic stiffnessconstants C ijkl is based on the definition [29]. σ ij = (cid:88) k,l C ijkl (cid:15) kl (29)where (cid:15) kl is the kl component of the strain tensor and σ ij is the ij component of the corresponding stress tensor.The stiffness constants are then determined by imposingdifferent kinds of strains (cid:15) kl and measuring the inducedstress responses and σ ij . Specifically, each simulation iscarried out using a cell in which only one of the six in-dependent strain components (measured with respect tothe undeformed cell reference) is different from zero. Inthis manner, there is only one term in the summation ofEq. (29). Exploring the linearity of Eq. (29), we deter-mine the stress response as a function of the magnitudeof a given strain component. The stiffness constants arethen given by the slopes of the σ ij ( (cid:15) kl ) graphs, C ijkl = dσ ij d(cid:15) kl . (30) IV. RESULTS AND DISCUSSION Figure [1] shows typical results for the stress compo-nents as a function of the imposed strain at a temperature FIG. 1: (color online). Stress response as a function of im-posed deformation. Magnitude of the error bars correspondsto size of symbols. Full lines depict linear fits to the PIMCresults. (a) Shear stress σ xz in response to shear strain (cid:15) xz ,scale of the left-hand side. (b) Tensile stress σ xx a functionof tensile strain (cid:15) zz , scale of the right-hand side. of 1 K and a molar volume of 20 cm . Figure [1-a] showsthe shear stress response σ xz as a function of the imposedshear strain (cid:15) xz , where the indications x and z refer tothe x and z directions of the computational cell. In allcases, we verified that the small distortions did not causeany disruptions of the crystal structure, leading only tosmall homogeneous deformations of the hexagonal struc-ture within the elastic limit. The behaviour is distinctlylinear up to absolute strain values of 2%, as attested bythe linear fit shown by the full line. As expected, theline passes through the origin, with no shear stress be-ing present at zero shear deformation. The slope of theline then determines the shear elastic constant C xzxz or,in Voigt notation [29], C , usually known as the shearmodulus µ . Figure [1-b] shows a similar response curve,plotting the tensile stress σ xx in response to the tensilestrain (cid:15) zz . In this case the stress response does not passthrough the origin given that the reference cell is in astate of hydrostatic compression. Once again the be-haviour is manifestly linear, with the slope giving thevalue of the elastic constant C xxzz or, in Voigt notation, C .Using stress-strain curves of the type shown in Fig. [1]obtained from the PIMC simulations, we determine the 66 elastic-constant matrix [29] of hcp solid He as a func-tion of the molar volume at a temperature of 1 K . Asexpected, the results are found to obey the symmetry re-lations for hexagonal crystals [29], leaving only five inde-pendent stiffness constants. These are shown in Fig. [2],which also includes sets of experimental data producedin the 1970s [30, 35]. A comparison between the theoret-ical and experimental data, however, must be conductedwith some care. Whereas the PIMC calculations give val-ues for the isothermal stiffness constants, the ultrasonicexperimental data typically probe the adiabatic elasticconstants. Accordingly, a direct comparison between thetwo data sets is meaningful only if the difference betweenthese two kinds of elastic constants is small. In its esti-mation, we applied the expressions for the difference be-tween the isothermal and adiabatic elastic constants [36]and used experimental thermodynamic data [37] for theisochoric heat capacity C ν and the isochoric pressure co-efficient ( ∂P/∂V ) ν . For a molar volume of 19.135 cm and a temperature of 1 K , we find the difference to beof the order of 101 bar, which is very small compared tothe absolute values of the elastic constants and their er-ror bars, thus justifying a direct comparison of the PIMCdata with experiment.Fig. [2-a] compares our PIMC results to experimen-tal data for the four elastic constants C , C , C , and C . For further comparison we have also included resultsfrom recent variational Monte Carlo (VMC) for T = 0 K based on the shadow wave function formalism [22]. Theagreement between PIMC and experiment is excellent for C , and C ,, with values essentially within each oth-ers error bars across the considered density range. Theagreement is also good for C , and C , with the PIMCvalues ∼ 10 - 20% below the experimental data. For thesefour constants, the PIMC results also show slightly betteragreement with experiment compared to the VMC data,which systematically overestimate all four constants by ∼ 20 - 30%.The situation is markedly different for the remainingindependent elastic constant, C , however. As shown bythe diamond symbols in Fig. [2-b] the PIMC data over-estimate experiment by 50 - 100%, with exception of thedata point measured by Greywall in 1971, although thelatter is characterized by an error bar of 300% [32]. Tofurther verify our result we performed additional com-putations at a molar volume of 20 cm . First, we car-ried out the PIMC calculations using a reduced timestep of τ = 1 / K − . Furthermore, we also determined C in an indirect way, computing the bulk modulus B = V ( ∂P/∂V ) T by a finite-difference derivative of thehydrostatic pressure with respect to volume changes andusing the relation C = (3 B − C ) / 2. As can be seen inFig. [2-b], neither led to significantly different results for C , lending further support to the internal consistencyof our calculations. Interestingly, the VMC result for C is actually in good agreement with experiment, differingby ∼ C + C = C + C , (31)which should be satisfied in case the c/a ratio is in-dependent of the density. This independence has been FIG. 2: (color online). (Color online) Elastic constants andsome of their specific combinations for hcp solid He as a func-tion of the molar volume at a temperature of 1K. (a) Filledsymbols depict PIMC results. Half-filled symbols representzero-temperature VMC results of Ref. [22]. Open symbolsrepresent corresponding experimental data.When error barsare not shown they are smaller than the symbol size. C (upward triangles), C (circles), C (squares), C (righttriangles). (b) C PIMC (filled diamonds, τ = 1 / − ,open upward triangle, τ = 1 / − , open downward trian-gle, result obtained from bulk modulus), C experimental(open diamonds), C VMC (half-filled diamond), C + C PIMC (o), C + C PIMC (pluses), C + C experimental(squares), C + C experimental (crosses). verified experimentally for a wide range of pressure val-ues [30]. This is consistent with the fact that the experi-mental values of the left- and right-hand sides of Eq. (32)are essentially equal, as shown by the open squares andcrosses, respectively, in Fig. [2-b] .Despite the overesti-mate for stiffness constant C , the PIMC results are infact consistent with this characteristic of solid He in thehcp phase as shown by the open circles and pluses in thesame plot.Another concerns the validity of the so-called Cauchyrelation [38, 39], C − P = C + P , where P is the ap-plied hydrostatic pressure. These are relations betweenelastic constants that are expected to hold when the (a)interatomic interactions can be described by purely two-body central forces and (b) when vibrational effects, zero-point or thermal in origin, are negligible. It is useful toquantify deviations from this relation through the pa-rameter [35–39] δ = C + P − ( C − P )( C − P ) . (32)For a classical crystal at T = 0K and characterized bypairwise central interaction forces one has δ = 0. Devia-tions from this value are then a measure of the magnitude FIG. 3: (color online). Cauchy deviation values for the hcp(squares) and bcc (circles) phases. Experimental results (opensymbols), PIMC results (full symbols). of vibrational effects and the importance of many-bodyinteraction forces. In the case of quantum crystals at lowtemperatures, the deviation thus is a probe for the roleof zero-point motion and/or many body interactions. Itseems plausible, however, that the latter are rather smallfor the condensed phases He [40].Fig. [3] shows δ as a function of the density for the ex-perimental data and our PIMC results for the hcp phase.The results show a distinct discrepancy. While the cal-culations give relatively small deviation values, rangingbetween 0.05 and 0.14, the experimental deviations are0.8-0.9. Again, this discrepancy can be traced back tothe difference in the value of stiffness constant C .To further assess these issues we also computed theelastic constants for the bcc phase, which exists at similardensities and temperatures as those considered for thehcp form. In this case, we find very good agreementbetween PIMC and experiment for all elastic constants.Specifically, at a density of 21 cm and a temperatureof 1.5 K , the PIMC calculations give C = 371 ± C = 330 ± 25, and C = 217 ± C =349 ± C = 301 ± 10, and C = 215 ± δ = ( C − C + 2 P ) / ( C − P ) in Fig. [3].The interpretation of the above results is challenging.In the present situation, our results agree very well withexperiment for four of the five independent elastic moduliof the hcp phase, but show a deviation for the remainingconstant, C . If this were to be interpreted as a flaw in our calculations, then discrepancies would also be ex-pected for other phases. This does not seem to be thecase. In addition to previous results for the liquid phasebased on very similar interaction models [14], our calcu-lations are in very good agreement with experiment forall elastic constants of the bcc phase. A further elementin this discussion is that, although conducted for zerotemperature, recent VMC and diffusion Monte Carlo re-sults [22, 41] do not seem to show this deviation for C in hcp He.To unveil this puzzle it would be extremely useful torenew the experimental data for the elastic constants ofhcp He. The data sets available today are more than30 years old and, with present experimental capabili-ties, it should be possible to obtain the elastic constantswith significantly greater accuracy. Not only would suchdata be useful toward gauging current theoretical mod-eling approaches, they would also be particularly useful,for instance, in understanding the elastic properties ofpolycrystalline He samples [11]. In addition, such infor-mation would also be of value in discriminating betweenchanges in elasticity and supersolidity, respectively, wheninterpreting frequency changes in TO experiments [11]. V. CONCLUSIONS In summary, in this paper we have reported values ofthe elastic stiffness constants of solid He in the hcp phasedetermined using the PIMC approach based on the Azizpair potential model. To this end we have developed anexpression for the stress observable in the path-integralformalism, allowing the direct measurement of the inter-nal stress state of a system in PIMC simulations. Thisdevelopment is of interest in its own right, allowing for in-stance, an explicit atomistic determination of the Peierlsstress for dislocation motion [22]. Here, we use it to com-pute the elastic stiffness constants by measuring the lin-ear stress response to imposed small strain conditions.Four of the five computed elastic stiffness constants asa function of density show good agreement with exper-iment. The stiffness coefficient C , which is ∼ C and C . The samecalculations for the bcc phase, on the other hand, showgood agreement between experiment and PIMC for allelastic constants. ACKNOWLEDGMENTS We gratefully acknowledge support from the Brazilianagencies CNPq, Fapesp, and Capes. We thank Prof.J. Beamish for helpful discussions. The calculationswere performed at CCJDR-IFGW-UNICAMP andCENAPAD-SP. ∗ Corresponding author: dekoning@ifi.unicamp.br [1] E. Kim and M. H. W. Chan, Nature (London) , 225(2004).[2] E. Kim and M. H. W. Chan, Science , 1543 (1970).[4] G. V. Chester, Phys. Rev. A , 256 (1970).[5] S. Balibar, Contemp. Phys. , 31 (2007).[6] J. Day and J. Beamish, Nature (London) , 853(2007).[7] M. H. W. Chan, Science , 1207 (2008).[8] X. Rojas, A. Haziot, V. Bapst, S. Balibar, and H. J.Maris, Phys. Rev. Lett. , 145302 (2010).[9] J. D. Reppy, Phys. Rev. Lett. , 255301 (2010).[10] J. Beamish, Physics , 51 (2010).[11] H. Maris and S. Balibar, J. Low Temp. Phys. , 5(2010).[12] H. Maris and S. Balibar, J. Low Temp. Phys. , 12(2011).[13] J. Day, O. Syshchenko, and J. Beamish, Phys. Rev. B , 214524 (2009).[14] D. M. Ceperley, Rev. Mod. Phys. , 279 (1995).[15] W. Krauth, Statistical Mechanics: Algorithms and Com-putations (Oxford University Press, Oxford, 2006).[16] R. A. Aziz, A. R. Janzen, and M. R. Moldover, Phys.Rev. Lett. , 1586 (1995).[17] M. Boninsegni, A. B. Kuklov, L. Pollet, N. V. Prokofev,B. V. Svistunov, and M. Troyer, Phys. Rev. Lett. 97,080401 (2006).[18] B. K. Clark and D. M. Ceperley, Comput. Phys. Com-mun. , 82 (2008).[19] M. Boninsegni, A. B. Kuklov, L. Pollet, N. V. Prokofev,B.V. Svistunov, and M. Troyer, Phys. Rev. Lett. ,035301 (2007).[20] S. G. S¨oyler, A. B. Kuklov, L. Pollet, N. V. Prokofev, andB. V. Svistunov, Phys .Rev. Lett. , 175301 (2009).[21] P. Sch¨offel and M. H. M¨user, Phys. Rev. B , 224108(2001).[22] R. Pessoa, S. A. Vitiello, and M. de Koning, Phys. Rev.Lett. , 085301 (2010). [23] M. Parrinello and A. Rahman, J. Appl. Phys. , 7182(1981).[24] M. Parrinello and A. Rahman, J. Chem. Phys. , 2662(1982).[25] M. E. Tuckerman, Statistical Mechanics: Theory andMolecular Simulation. Oxford University Press, Oxford,(2010).[26] K. P. Essler, Advancements in the Path Integral MonteCarloMethod for Many-body Quantum Systems at Fi-nite Temperature,Ph.D. thesis, University of Illinois atUrbana-Champaign (2006).[27] Y. A. Freiman, S. M. Tretyak, A. Grechnev, A. F. Gon-charov, J. S.Tse, D. Errandonea, H.-k. Mao, and R. J.Hemley, Phys. Rev. B , 094112 (2009).[28] B. Bernu and D. Ceperley, Quantum Simulations of Com-plex Many-Body Systems: From Theory to Algorithms,NIC Series,Vol. 10, edited by A. M. J. Grotendorst andD. Marx (John von Neumann Institute for Computing,J´’ulich, 2002), p. 51.[29] J. F. Nye, Physical Properties of Crystals: Their Rep-resentation by Tensors and Matrices (Oxford UniversityPress, Oxford, 1985).[30] J. P. Franck and R. Wanner, Phys. Rev. Lett. , 345(1970).[31] R. H. Crepeau, O. Heybey, D. M. Lee, and S. A.Strauss,Phys. Rev.A , 1162 (1971).[32] D. S. Greywall, Phys. Rev. A , 2106 (1971).[33] R. Wanner and J. P. Franck, Phys. Rev. Lett. , 365(1970).[34] D. S. Greywall, Phys. Rev. B 16, 5127 (1977).[35] J. Beamish, Handbook of Elastic Properties of Solids,Liquids,and Gases, Vol. 2, edited by M. Levy and L. Furr(Academic Press,San Diego, 2000).[36] D. C. Wallace, Thermodynamics of Crystals (Dover Pub-lications,New York, 1998), p. 26.[37] G. Ahlers, Phys. Rev. A , 1505 (1970).[38] M. Born and K. Huang, Dynamical Theory of CrystalLattices (Clarendon, Oxford, 1954).[39] C.-S. Zha, H.-k. Mao, and R. J. Hemley, Phys. Rev. ,174107(2004).[40] S. Ujevic and S. A. Vitiello, J. Chem. Phys. , 8482(2003).[41] C. Cazorla, Y. Lutsyshyn, and J. Boronat, Phys. Rev. B85