First-principles modeling of three-body interactions in highly compressed solid helium

First-principles modeling of three-body interactions in highly compressed solid helium

Claudio Cazorla School of Materials Science and Engineering, UNSW Australia, Sydney NSW 2052, Australia
Integrated Materials Design Centre, UNSW Australia, Sydney NSW 2052, Australia
   Jordi Boronat Departament de Física i Enginyeria Nuclear, Universitat Politècnica de Catalunya, Campus Nord B4-B5, E-08034, Barcelona, Spain

We present a new set of three-body interaction models based on the Bruch-McGee (BM) potential that are suitable for the study of the energy, structural and elastic properties of solid He at high pressure. Our ab initio three-body potentials are obtained from the fit to total energies and atomic forces computed with the van der Waals density functional theory method due to Grimme, and represent an improvement with respect to previously reported three-body interaction models. In particular, we show that some of the introduced BM parametrizations reproduce closely the experimental equation of state and bulk modulus of solid helium up to a pressure of  GPa, when used in combination with standard pairwise interaction models in diffusion Monte Carlo simulations. Importantly, we find that recent predictions reporting a surprisingly small variation of the kinetic energy and Lindeman ratio on quantum crystals under increasing pressure are likely to be artifacts produced by the use of incomplete interaction models. Also, we show that the experimental variation of the shear modulus, , at  GPa can be quantitatively described with the new set of three-body BM potentials. At higher pressures, however, the agreement between our results and experiments deteriorates and thus we argue that higher order many-body terms in the expansion of the atomic interactions probably are necessary in order to better describe elasticity in very dense solid He.


I Introduction

The electronic structure of a single He atom is among the simplest in the periodic table of elements. Likewise, the atomic interactions in liquid and solid helium can be reproduced accurately with simple analytical functions that solely depend on the distance between particles taken in pairs. Examples of successful He–He interaction models include the Lennard-Jones and Aziz-type semiempirical potentials. kalos81 ; boronat94 ; aziz Yet, under conditions of large pressures and strain deformations the interparticle interactions become more complex due to the strong electronic repulsion experienced by neighboring atoms. Consequently, pairwise potentials, which work reasonably well under near-equilibrium conditions, turn out to be unreliable. This is, for instance, the case of the Aziz-II potential, aziz which at high pressure provides too repulsive atomic forces and a significant overestimation of the He molar volume and bulk modulus. cazorla08

A recently proposed straightforward way to correct for such modeling drawbacks consists in modifying the repulsive part of standard pairwise potentials by means of an exponential attenuation factor. moraldi12 This possibility has already been explored in highly compressed solid He azizbc and molecular hydrogen omiyinka13 with quantum Monte Carlo simulations, producing equations of state which are in very good agreement with experiments. Nevertheless, the use of modified pairwise potentials in very dense crystals poses a series of issues and open questions. For instance, a surprisingly small variation of the kinetic energy upon increasing pressure have been reported in works [azizbc, ] and [omiyinka13, ], and, owing to the lack of experimental data in the thermodynamic regime of interest, it remains to be demonstrated whether such predictions can be fully ascribed to genuine quantum nuclear effects or not. Also, pairwise potentials are in general not recommended for the study of elasticity in hcp crystals at high pressure since they inevitably lead to null values of the Cauchy discrepancy (defined as the difference between the two elastic constants and ), in contrast to what is observed in experiments. wallace72 ; pechenik08 ; freiman09 ; zha04

An alternative route to improve the description of quantum solids under extreme stress-strain conditions is to consider higher order terms, beyond pairwise additivity, in the approximation to the atomic interactions. In this context, several three-body interatomic models have already been proposed like, for instance, the Axilrod-Teller (AT), Bruch-McGee (BM), and Cohen-Murrel (CM) potentials. boronat94 ; bm73 ; cohen96 However, improvements resulting from the use of those three-body interaction models so far have been reported to be only marginal. For instance, three decades ago Loubeyre claimed, based on the outcomes of self-consistent phonon and classical Monte Carlo simulations, that the BM three-body interaction could bring into good agreement calculations and experiments performed on the equation of state of solid helium up to  GPa. loubeyre86 However, Chang et al. chang01 have shown more recently that when either the BM or CM three-body potentials are considered in quantum Monte Carlo simulations the resulting He molar volumes are significantly underestimated, already at few GPa. Similar discouraging results have been reported also by other authors who have employed analogous three-body interaction models. herrero06 ; tian06

In this article, we present new work done on the modeling of three-body interactions in highly compressed solid helium up to pressures of  GPa. We introduce a new set of BM potential parametrizations obtained from fits to ab initio energies and atomic forces calculated with the van der Waals corrected density functional theory method due to Grimme (DFT-D2). grimme06 We show that an overall improved description of the energy, elastic and structural properties of solid helium can be achieved with some of the introduced BM three-body interatomic potentials, when used in combination with pairwise potentials in quantum Monte Carlo simulations. Our work also brings new insight into the physics of quantum crystals at high pressure. For instance, we show that previously reported small variations of the kinetic energy, , and Lindeman ratio, , in solid helium under pressure azizbc are likely to be artifacts deriving from the use of incomplete atomic interaction models. Moreover, we quantify the role of quantum nuclear effects on the estimation of the shear modulus, , and conclude that they become secondary when pressure is raised. Finally, at  GPa we find that the agreement between our results and experiments starts to worsen. Therefore, we argue that higher order many-body terms in the expansion of the atomic interactions probably are necessary in order to describe elasticity in dense solid helium more accurately.

The organization of this article is as follows. In the next section, we outline the employed computational methods and provide the technical details in our calculations. In Sec. III, we explain the fitting strategy that we have followed to obtain the new set of three-body interaction models. Next, we present our results on the equation of state, kinetic energy, and structural and elastic properties in solid helium, together with some discussion. Finally, we summarize our main findings in Sec. V.

Ii Computational Methods

We used the density functional theory method including van der Waals corrections due to Grimme, grimme06 to compute the interactions and forces between helium atoms in the hexagonal close package (hcp) crystal structure, from equilibrium up to a pressure of  GPa. (Details of our ab initio DFT-D2 calculations can be found in elsewhere, azizbc hence we highlight here only the main technical features.) Subsequently, we found a series of three-body interaction models that, when used in combination with the pairwise Aziz-II potential aziz (hereafter denoted as ), reproduced very closely the obtained DFT-D2 results. The details of our fitting strategy are comprehensively explained in Sec. III. Next, we performed diffusion Monte Carlo (DMC) simulations in which the new three-body interaction models were used to calculate the energy, structural, and elastic properties of solid helium under pressure. In this section, we explain the specific implementation of the DFT-D2 and DMC methods in our work.

ii.1 Density functional theory

We chose the generalized gradient approximation to density functional theory proposed by Perdew, Burke, and Ernzerhof (GGA-PBE), pbe96 as is implemented in the VASP package. vasp Van der Waals interactions were taken into account by adding an attractive energy term to the exchange-correlation energy of the form (where indexes and label different particles, is a constant, and a damping factor is introduced at short distances to avoid divergences). grimme06 ; cazorla15 The projector-augmented-wave technique blochl94 ; kresse99 was employed to represent the core electrons since this approach has been shown to provide very accurate total energies and is computationally very efficient. cazorla07 ; taioli07 The electronic wave functions were represented in a plane-wave basis truncated at  eV, and for integrations within the first Brillouin zone (BZ) we employed dense -centered -point grids of . By using these parameters we obtained interaction energies that were converged to within  K per atom. Geometry relaxations were performed by using a conjugate-gradient algorithm that kept the volume of the unit cell fixed and permitted variations of its shape. The imposed tolerance on the atomic forces was  eVÅ. With such a DFT-D2 setup we calculated the total energy and shear modulus in solid He in the volume interval  Å/atom.

Additionally, we computed the vibrational phonon spectrum in solid He at eight different volumes by means of the “direct approach”. In the direct approach the force-constant matrix is directly calculated in real-space by considering the proportionality between atomic displacements and forces when the former are sufficiently small. cazorla13b ; shevlin12 ; cazorla08c In this case, large supercells have to be simulated in order to guarantee that the elements of the force-constant matrix have all fallen off to negligible values at their boundaries, a condition that follows from the use of periodic boundary conditions. alfe09 Once the force-constant matrix is obtained, we Fourier-transform it to obtain the phonon spectrum at any -point. The quantities with respect to which our DFT-D2 phonon calculations need to be converged are the size of the supercell and atomic displacements, and the numerical accuracy in the atomic forces. The following settings were found to fulfill our convergence requirement of correct zero-point energy corrections to within  K/atom: azizbc ; cazorla13b supercells (that is, repetitions of the hcp unit cell containing a total of atoms), and atomic displacements of  Å . Regarding the calculation of the atomic forces with VASP, we found that the density of -points had to be increased slightly with respect to the value used in the energy calculations (i.e., from to ) and that computation of the non-local parts of the pseudopotential contributions needed to be performed in reciprocal, rather than real, space.

ii.2 Diffusion Monte Carlo

In our DMC simulations, we used a guiding wave function, , that accounts simultaneously for the atomic periodicity and Bose-Einstein quantum symmetry in He crystals. This model wave function is expressed as cazorla09b


where indexes and run over particles and perfect lattice positions, respectively. In previous works we have shown that provides an excellent description of the ground-state properties of bulk hcp He and other similar quantum systems. cazorla09b ; cazorla08b ; cazorla10 ; gordillo11 The correlation factors in Eq. (1) were expressed in the McMillan,  , and Gaussian, , forms. Parameters and were optimized at each density point by using the variational Monte Carlo (VMC) method. For instance, at  Å  we obtained  Å  and  Å , and at  Å ,  Å  and  Å . We note that our choice of the guiding function was motivated by an interest in studying the possible effects of quantum atomic exchanges on the energetic and elastic properties of dense helium. However, we realised by direct comparison to the results obtained with non-symmetric wave function models in analogous DMC simulations, azizbc that such effects can be totally neglected in practice.

The technical parameters in our calculations were set to ensure convergence of the total energy per particle to less than  K. The value of the mean population of walkers was and the length of the imaginary time-step ()  K . We used simulation cells containing atoms. Numerical bias stemming from the finite size of the simulation box were minimised by following the variational correction approach explained in works [cazorla08, ] and [azizbc, ]. Statistics were accumulated over DMC steps performed after system equilibration, and the approximation used for the short-time Green’s function, , was accurate to second order in boronat94 ; chin90 The computational strategy that we followed to calculate the shear modulus was the same than in Refs. [cazorla12, ; cazorla13, ; cazorla12b, ].

Figure 1: (Top) Energy differences between the DFT-D2 method and potentials calculated on a reference set of configurations (see text). Details are magnified in the inset. (Bottom) Results of our fit obtained in the case of the atomic forces. stands for the difference in the atomic forces between the DFT-D2 method and many-body potentials, for the variance of the atomic forces computed with the DFT-D2 method, and for the average performed over particles and Cartesian components.

Figure 2: Phonon spectrum obtained with the DFT-D2 method (dashed lines) and the (BM-F) interaction model (solid lines), which was determined considering only the atomic forces in the corresponding fit (see text).

Iii Fitting strategy and three-body potential models

Our three-body potential matching algorithm ercolessi94 ; marco12 ; marco14 is based on a least square fit to the DFT-D2 reference data, that consists of total energies and atomic forces. The objective function to be minimized is given by


where is the number of reference configurations, the number of particles on each configuration, and and a weight assigned to the energy, , and force, , contributions to , respectively. With this definition of the objective function we ensure that despite different magnitudes are expressed in different units all them are normalized and contribute equally to . Subscripts “DFT” and “FF” refer to the DFT-D2 and classical potential results, respectively.

The set of reference configurations in our fit comprised the structures used in the calculation of the He vibrational phonon spectra in the interval  Å/atom by means of the “direct approach” (see Sec. II.1). cazorla13b ; shevlin12 ; cazorla08c Such atomic arrangements were generated by taking the relaxed hcp lattice supercells (, space group ) at different volumes and displacing one of the atoms sitting in an inequivalent Wyckoff position a distance of  Å  first along the direction (where represent the normalised Cartesian vectors), and then along (that is, we created two different atomic configurations at each volume). The reason for our choice was that we wanted to reproduce simultaneously the energy and elastic properties in highly compressed solid He. In fact, the atomic forces are defined as minus the first derivative of the total energy with respect to the atomic positions, whereas the elastic constants involve the second derivative of the total energy with respect to strain deformations. In spite of this apparent disconnection, atomic forces and elastic constants are indirectly related by the corresponding spectrum of vibrational phonon frequencies. Namely, on one side, phonons can be calculated from the variation of the atomic forces upon the displacement of atoms away from their equilibrium positions, and, on the other side, elastic constants can be estimated from the slope of specific acoustic branches in the vicinity of the point in reciprocal space (that is, in the limit). Therefore, even though we did not explicitly consider second derivatives in our definition of the objective function , we expected to achieve an acceptable description of elasticity in solid helium. We shall come back to this point later on this section.

The classical potential adopted in this study, denoted as “FF” in Eq. 2, is given by , where represents the pairwise Aziz-II interaction model aziz and the three-body Bruch-McGee (BM) potential given by bm73


where , and , , and , are the interior angles of the triangle formed by the atoms labelled , , and . is an attractive potential term representing triple dipole and three-body exchange interactions. Parameters , , and were varied during the minimization of the objective function (see Eq. 2). For this, we used a quadratic polynomial interpolation line-search with the directions found using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) formula. nocedal06 The gradient of the objective function was calculated analytically since otherwise numerical bias developed that impeded convergence. Actually, the typical size of the involved atomic forces is very small, of the order of  eV/Å , hence they needed to be calculated very precisely. The minimizations were stopped when all the gradients of the objective function in absolute value were smaller than . Typically, this was achieved within minimization loops when starting from a reasonable initial guess of the , , and parameters (e.g., the original values proposed by Bruch and McGee [bm73, ]).

 (K)  (K)  ()
(BM) [bm73, ]
Table 1: Bruch-McGee three-body potential parameters corresponding to the original model, (BM), and the new ones introduced in the present work. The (BM-E) set has been obtained by considering exclusively DFT-D2 energies on the fit [, ], the (BM-F) the atomic forces [, ], and (BM-EF) a combination of ab initio energies and atomic forces [, ] (see text).

Table I shows the values of the parameters obtained in our fits, in which we considered three different possibilities based on the choice of the relative energy and forces weights:   and , hereafter denoted as (BM-E),   and , (BM-F), and   and , (BM-EF). Our results differ appreciably from the original values proposed by Brunch and McGee [which hereafter are denoted as (BM)]. For instance, becomes negative when the energies are taken into account in the fit, and and systematically turn out to be larger.

In Figure 1, we illustrate the quality of our fits by plotting the energies and forces calculated on each reference configuration. For comparison purposes, we also enclose the results obtained with the original (BM) potential (i.e., with function ). For the sake of simplifying the notation, we only indicate the three-body part in the corresponding many-body potential. This convention will be adopted throughout the text if not stated otherwise. As is appreciated in the figure, (BM-E) reproduces the DFT-D2 energies more closely than any other model (as expected) whereas (BM) provides the worst description. The energies obtained with the (BM-EF) potential can be regarded also as fairly good. As for the atomic forces, (BM-F) produces the best results, as expected, and (BM), again, turns out to be the less reliable. In this latter case, the forces obtained with the (BM-EF) and, surprisingly also, (BM-E) potentials are not too distant from the reference DFT-D2 data.

In Fig 2, we plot the vibrational phonon spectra obtained with the DFT-D2 method and the (BM-F) potential in solid He at the highest analysed pressure (probably the most challenging case to be reproduced with a potential function, see Fig. 1). We note that the agreement between the two sets of data can be regarded as fairly good. The largest differences are found on the optical branches, which correspond to the highest vibrational frequency values. The DFT-D2 acoustic phonon modes in the vicinity of the point, however, are reasonably well reproduced by (BM-F). These outcomes demonstrate that, as we suggested above, by considering the atomic forces in the definition of in principle one can obtain a reasonable description of the elasticity in the reference system.

Figure 3: Zero-temperature equation of state calculated in helium with the DFT-D2 and DMC methods. In the DMC case, different pairwise and three-body interaction models have been employed. Experimental data from Ref. [loubeyre93, ] are shown for comparison. Inset: The high- region in the curves are magnified in order to appreciate better the differences.

Iv Results and Discussion

iv.1 Equation of state

We show the results of our calculations on the equation of state, , of solid helium in Fig. 3, together with experimental data from work [loubeyre93, ]. The DFT-D2 series was obtained with the ab initio methods explained in Sec. II.1, including quantum zero-point energy corrections. The other results were obtained in diffusion Monte Carlo (DMC) simulations using the indicated interaction potentials, as explained in Sec. II.2 and elsewhere [azizbc, ]. Labels “” and “(BC)” stand respectively for the pairwise potential due to Aziz aziz and a modified version of the former that we have recently introduced in work [azizbc, ]. The DMC (DFT-D2) calculations were performed at  () different volumes spanned in the interval  Å/atom. In each case, the resulting total energies were fitted to a third order Birch-Murnaghan equation of the form birch78 ; cazorla09


where is the value of the bulk modulus at the equilibrium volume , with , and all the derivatives are calculated at zero pressure. For reproducibility purposes, we enclose the , , and parameters obtained in all our fits in Table II.

 (Å)  (eV/Å)
Table 2: Parameters corresponding to the fit of our equation of state results to Birch-Murnaghan functions, see Eq. (4), as obtained with different computational approaches. In the DMC case, different pairwise and three-body potentials have been considered for the description of the interatomic forces.

Figure 4: Atomic kinetic energy calculated in He with the DFT-D2 and DMC methods and expressed as a function of pressure. In the DMC case, different pairwise and three-body interaction models have been considered for the description of the atomic interactions.

Very good agreement is obtained between our DFT-D2 results and experiments. This outcome justifies in part our choice of the DFT-D2 results as reference data in modeling of the many-body interactions. Likewise, the curves obtained with the (BC), (BM-E), and (BM-EF) potentials are also very close to the observations. We notice that the (BC) model introduced in Ref. [azizbc, ] was constructed to reproduce the equation of state calculated with the DFT-D2 method and that the good agreement displayed in Fig. 3 is not a new result. Contrarily, the , (BM), and (BM-F) potentials provide a poor description of the variation of the volume under pressure. In particular, we find that the (BM) potential systematically underestimates at pressures equal or larger than  GPa, in accordance with previous results reported by other authors. chang01 ; herrero06 Meanwhile, the and (BM-F) interaction models significantly overestimate the same quantity at pressures also close to or larger than  GPa. In this latter case, we notice a surprising resemblance between the two calculated curves.

The main conclusion emerging from this part of our study is that the new (BM-E) and (BM-EF) three-body potentials reproduce very accurately the equation of state of solid helium up to a pressure of  GPa (and possibly beyond). To the best of our knowledge, such a good agreement between theory and experiments has not been reported before for any known potential in solid He (see work [chang01, ]).

iv.2 Kinetic energy

Our kinetic energy, , results are shown in Fig. 4. In our DFT-D2 calculations, the kinetic energy was estimated within the quasiharmonic approximation through the expression


where are the vibrational phonon frequencies in the crystal calculated at wave vector and phonon branch , which depend on the volume, and the total number of wave vectors used for integration within the first Brillouin zone (see Sec. II.1 and works [azizbc, ; cazorla13b, ]). usually is referred to as the “zero-point energy” (ZPE) and in many computational studies turns out to be crucial for predicting accurate solid-solid phase transitions. shevlin12 ; cazorla08c ; cazorla09 Regarding our DMC calculations, we computed first the exact potential energy, , by means of the pure estimator technique barnett91 ; casulleras95 and subsequently obtained the exact kinetic energy by subtracting to the corresponding total energy. In all the cases, spline interpolations were applied to the calculated data points in order to obtain smooth -dependent energy curves (lines in Fig. 4).

As is appreciated in the figure, the DFT-D2 results differ enormously from the rest of series obtained with pairwise and three-body potentials in our DMC simulations. At the highest analysed pressure, for instance, the DFT-D2 kinetic energy is a factor of two larger than the obtained DMC value. Given the lack of experimental data in the thermodynamic regime of interest, we can not rigorously conclude which type of calculation is providing the most realistic description. Nevertheless, we think that the DFT-D2 results are overestimating severely because they have been obtained using the quasiharmonic approximation. In fact, it has been already demonstrated that the quasiharmonic approximation is not appropriate for studying crystals that behave much more classically than solid helium like, for instance, molecular hydrogen, geneste12 ; mcmahon12 ; morales13 ammonia, datchi06 ; pickard08 and some alkali metals. errea11 ; feng15 It is worth noticing here that although the quasiharmonic DFT-D2 approach can produce equations of state that are in very good agreement with experiments (as it has been shown in Sec. IV.1), the accompanying ZPE corrections have a lot of margin for error since at high these are always several orders of magnitude smaller than the energy of the perfect crystal lattice. We shall comment again on this point in the next paragraph.

Figure 5: (Top) Atomic density profile around the perfect lattice positions calculated with the DMC method considering different pairwise and three-body interaction models ( Å/atom). Solid lines correspond to Gaussian curves fitted to the results. The corresponding tails are magnified in the inset in order to better appreciate the differences. (Bottom) Lindeman ratio calculated in solid He with the DFT-D2 and DMC methods, expressed as a function of pressure.

It is interesting to analyse the differences found between the (full quantum) DMC results obtained with different pairwise and three-body potential models. The (BC) curve shows a plateau around  K at pressures equal and beyond  GPa. In a recent work, azizbc we identified such an infinitesimal variation in the kinetic energy with the presence of extreme quantum nuclear effects. However, calculations performed with the new set of three-body potentials introduced in this work bring new light into our previous interpretation of the (BC) results. As is observed in Fig. 4, the (BM-E), (BM-F), and (BM-EF) curves consistently display a small but steady increase in the kinetic energy under compression. At pressures below  GPa the pairwise and three-body interaction models roughly provide equivalent results however at  GPa the differences between them are as a large of  K, with the potentials providing always the largest values. Several conclusions can be drawn from these results. First, although attenuated pairwise potentials based on exponential prefactors moraldi12 can fairly reproduce experimental data, azizbc ; omiyinka13 they are likely to introduce unwanted bias on the calculation of the kinetic energy. And second, the large discrepancies observed between the DFT-D2 and results do not seem to be originated by the absence of four-, five- and so on many-body interactions in the DMC calculations. Actually, by comparing the energy curves obtained in the and cases one realizes that the effect of considering three-body interactions on is rather small [only in the (BM) case those effects are not negligible, although certainly minor]. Therefore, it is reasonable to expect similar trends when eventually one would add higher order many-body terms in the description of the atomic interactions. In regard to this last point, we notice that one of the main conclusions presented in work [azizbc, ], namely that the quasiharmonic DFT approach exceedingly overestimates in dense He, appears to be valid.

iv.3 Structural properties

An analysis of the atomic structure in solid He at high pressure will allow us to understand better the origins of the discrepancies found so far between the (BC) and potentials. Figure 5 shows the atomic density profiles, , and Lindeman ratio, , calculated using the DMC method and several atomic interaction models. The results (see top panel) were obtained at volume  Å/atom and subsequently were fitted to Gaussian functions (solid lines in the figure). As is observed there, the (BC) curve is noticeably broader than all the others, and its value at the origin is about  % of that calculated with the (BM) potential. Meanwhile, the (BM-E) and (BM-EF) profiles are practically indistinguishable and slightly higher near zero than the one obtained in the (BM) case. Clearly, the (BC) potential produces a much larger atomic delocalization than the rest of interaction models, which is consistent with the kinetic energy results explained in the previous section.

As for the Lindeman ratio (see bottom panel in Fig. 5), we have estimated the corresponding dependence on pressure for each analysed potential. In the DFT-D2 case, was computed within the quasiharmonic approximation using the formula , see Eq. (5) and works [cazorla08b, ; arms03, ]. The results obtained in the (BC) case are already known: a plateau around appears at pressures larger than  GPa. azizbc However, all the other interaction models, including and (BM), provide much smaller values of at similar conditions. Moreover, the computed Lindeman ratio curves get depleted when compression is raised [with the exception of (BM), in which saturates around at pressures larger than  GPa]. This latter trend is also observed in the DFT-D2 series, which systematically lies below the DMC predictions.

The results presented in this section show that the (BC) potential produces an unusually large delocalization of the atoms, which is at odds with the trends realised in the rest of cases. Such a huge particle dispersion effect is the responsible for the flat kinetic energy curve appearing in Fig. 4, which is likely to be an artifact deriving from the use of exponential attenuation factors at short distances.

Figure 6: Calculated bulk modulus in He using the DFT-D2 and DMC methods and expressed as a function of volume. In the DMC case, different pairwise and three-body interaction models have been considered. Experimental data from work [zha04, ] are shown for comparison. Inset: The high- region in the curves are magnified in order to appreciate better the differences.

Figure 7: Calculated shear modulus in He using the DFT-D2 method and several force fields considering the atoms immobile in the perfect lattice positions. Experimental data are from work [zha04, ].

iv.4 Elastic properties

In Figs. 6 and 7, we show the bulk and shear modulii, and respectively, calculated in solid helium under pressure. The bulk modulus was directly obtained from the Birch-Murnaghan fits explained in Sec. IV.1, and in the case spline interpolations were applied to the calculated data points in order to obtain smooth -dependent curves.

Concerning the analysis of our results, this is very much similar to the conclusions presented for the equation of state in Sec. IV.1. Essentially, the DFT-D2, (BC), (BM-E), and (BM-EF) curves are in good agreement with experiments whereas the , (BM-F), and (BM) curves are not. In this latter case, both and (BM-F) series are very similar and significantly overestimate the bulk modulus at small volumes. Likewise, the (BM) potential provides unrealistically small values of at large densities.

 [aziz, ]
(BC) [azizbc, ]
(BM) [bm73, ]
Table 3: Summary of the performance of the pairwise and three-body atomic interaction models analysed in this work in describing the energy, structural, and elastic properties of solid He at high pressure. Symbol () indicates correct (incorrect) description of the considered quantity, whereas means quantitatively correct up to a certain pressure. Question mark “” denotes a certain hesitation due to lack of experimental data in the high pressure regime of interest.

Let us now comment on the results shown in Fig. 7. All the values have been obtained considering the atoms fixed on their perfect lattice positions, that is, totally neglecting likely quantum nuclear effects (hence the employed subscript). This is done for the sake of comparison since it is technically difficult to account for quantun nuclear effects in the DFT-D2 calculations in an exact manner, that is, to go beyond the quasiharmonic approximation. Nevertheless, later on this section we will show that according to our DMC simulations quantum nuclear effects become secondary on at high pressure. As is observed in the figure, the DFT-D2 curve is in overall good agreement with the ambient temperature measurements performed by Zha and collaborators. zha04 Again, these findings justify our choice of the benchmark data for the modeling of three-body interactions. Regarding the performance of the original and new three-body BM potentials, we find that in general they reproduce quite satisfactorily the experimental data obtained at volumes larger than  Å/atom (i.e.,  GPa). This is especially true in the (BM-F) case where, as expected (see Sec. III), the calculated shear modulii follow closely those obtained with the DFT-D2 method. However, at volumes smaller than  Å/atom (i.e.,  GPa) we find that the differences between the BM curves (including the BM-F case), on one side, and the DFT-D2 results and experiments, on the other, become increasingly larger. We recall that the (BM-E) and (BM-EF) potentials provide a very good description of the equation of state and bulk modulus, whereas the (BM-F) potential does not. This appreciation let us to conclude that is very difficult to provide simultaneously a good account of the energy and elastic properties in solid helium by using an effective three-body approach. Higher order many-body contributions in the description of the atomic interactions probably are necessary in order to attain an overall correct description of solid helium at high pressure. As for the pairwise potentials, performs very similarly to the (BM-F) model, as we have also noted in the total energy (see Sec. IV.1) and bulk modulus cases. The (BC) model, however, remarkably fails in reproducing the variation of the shear modulus under pressure. Moreover, it predicts the occurrence of unrealistic mechanical instabilities (i.e., sinko02 ; grimvall12 at small volumes. Therefore, the use of the (BC) potential is strongly not recommended for the simulation of solid helium at high pressure.

In order to quantify the importance of quantum nuclear effects on the calculation of the shear modulus, we carried out additional quantum DMC calculations (see Sec. II.2 and works [cazorla12, ; cazorla13, ; cazorla12b, ] for details). To our surprise, we found that the quantum and classical shear modulii results are very similar. For instance, in the (BM-F) case the difference (where subscript “quantum” means calculated with the DMC method) amounts only to  GPa at  GPa. Similar results were obtained also in the rest of and cases. We note that the sign of the differences is always positive, thus the inclusion of quantum nuclear effects tends to lower the classical values, although in a small fraction (i.e.,  %). This last finding appears to be consistent with conclusions presented in a recent quantum Monte Carlo study by Borda et al.borda14 in which the ideal shear strength on the basal plane of hcp He was found to behave analogously than in classical solids.

V Conclusions

In Table III we summarise the performance of the analysed pairwise and three-body potentials in describing the energy, elastic and structural properties of solid He at high pressure. A number of tips can be drawn from our results. First of all, the use of pairwise potentials in general is not recommended. These either fail to reproduce the equation of state and bulk modulus, i.e., , or the kinetic energy, and structural and elastic features, i.e., (BC), in highly compressed quantum crystals. In this context, we urge to employ more versatile many-body interaction models. This is the case, for instance, of the new three-body BM potentials introduced in this work, which represent an improvement with respect to previously reported similar models. Overall, we recommend to consider the (BM-E) and (BM-EF) parametrizations in prospective simulation studies because they provide the most satisfactory general description of dense solid He. Indeed, those interaction models can be safely employed, for instance, in atomistic high- high- simulations (either classical or quantum), which are of relevance to planetary sciences. Nevertheless, we must note that it remains a challenge to attain a precise description of elasticity at high pressure by using effective three-body potentials, thus in this latter case consideration of higher order many-body terms appears to be necessary.

Importantly, we have shown that the addition of three-body forces corrects for the artificially large atomic delocalization found with modified pairwise potentials based on exponential attenuation factors. Nevertheless, given the lack of structural and kinetic energy measurements performed at high pressure, we have not been able to quantify the accuracy of our and DMC results obtained with the (BM-E) and (BM-EF) potential models. In this regard, advanced computational studies in which both the nuclear and electronic degrees of freedom in the crystal were to be treated at the quantum level are highly desirable.

This research was supported under the Australian Research Council’s Future Fellowship funding scheme (projects number RG134363 and RG151175), and MICINN-Spain (Grants No. MAT2010-18113, CSD2007-00041, and FIS2014-56257-C2-1-P).


  • (1) M. H. Kalos, M. A. Lee, and P. A. Whitlock, Phys. Rev. B 24, 115 (1981).
  • (2) J. Boronat and J. Casulleras, Phys. Rev. B 49, 8920 (1994).
  • (3) R. A. Aziz, F. R. W. McCourt, and C. C. K. Wong, Mol. Phys. 61, 1487 (1987).
  • (4) C. Cazorla and J. Boronat, J. Phys.: Condens. Matt. 20, 015223 (2008).
  • (5) M. Moraldi, J. Low Temp. Phys. 168, 275 (2012).
  • (6) C. Cazorla and J. Boronat, Phys. Rev. B 91, 024103 (2015).
  • (7) T. Omiyinka and M. Boninsegni, Phys. Rev. B 88, 024112 (2013).
  • (8) D. C. Wallace, Thermodynaimcs of Crystals (Wiley, New York, 1972).
  • (9) E. Pechenik, I. Kelson, and G. Makov, Phys. Rev. B 78, 134109 (2008).
  • (10) Y. A. Freiman, S. M. Tretyak, A. Grechnev, A. F. Goncharov, J. S. Tse, D. Errandonea, H.-K. Mao, and R. J. Hemley, Phys. Rev. B 80, 094112 (2009).
  • (11) C.-S. Zha, H.-K. Mao, and R. J. Hemley, Phys. Rev. B 70, 174107 (2004).
  • (12) L. W. Bruch and I. J. McGee, J. Chem. Phys. 59, 409 (1973).
  • (13) M. J. Cohen and J. N. Murrel, Chem. Phys. Lett. 260, 371 (1996).
  • (14) P. Loubeyre, Phys. Rev. Lett. 58, 1857 (1986).
  • (15) S.-Y. Chang and M. Boninsegni, J. Chem. Phys. 115, 2629 (2001).
  • (16) C. P. Herrero, J. Phys.:Condens. Matt. 18, 3469 (2006).
  • (17) C.-L. Tian, F.-S. Liu, F.-Q. Jing, and L.-C. Cai, J. Phys.:Condens. Matt. 18, 8103 (2006).
  • (18) S. Grimme, J. Comp. Chem. 27, 1787 (2006).
  • (19) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  • (20) G. Kresse and J. Fürthmuller, Phys. Rev. B 54, 11169 (1996).
  • (21) C. Cazorla, Coord. Chem. Rev. 300, 142 (2015).
  • (22) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  • (23) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  • (24) C. Cazorla, M. J. Gillan, S. Taioli, and D. Alfè, J. Chem. Phys. 126, 194502 (2007).
  • (25) S. Taioli, M. J. Gillan, C. Cazorla, and D. Alfè, Phys. Rev. B 75, 214103 (2007).
  • (26) C. Cazorla and J. iguez, Phys. Rev. B 88, 214430 (2013).
  • (27) S. A. Shevlin, C. Cazorla, and Z. X. Guo, J. Phys. Chem. C 116, 13488 (2012).
  • (28) C. Cazorla, D. Alfè, and M. J. Gillan, Phys. Rev. Lett. 101, 049601 (2008).
  • (29) D. Alfè, Comp. Phys. Commun. 180, 2622 (2009).
  • (30) C. Cazorla, G. Astrakharchick, J. Casulleras, and J. Boronat, New Journal of Phys. 11, 013047 (2009).
  • (31) C. Cazorla and J. Boronat, Phys. Rev. B 77, 024310 (2008).
  • (32) C. Cazorla, G. Astrakharchick, J. Casulleras, and J. Boronat, J. Phys.: Condens. Matter 22, 165402 (2010).
  • (33) M. C. Gordillo, C. Cazorla, and J. Boronat, Phys. Rev. B 83, 121406(R) (2011).
  • (34) S. A. Chin, Phys. Rev. A 42, 6991 (1990).
  • (35) C. Cazorla, Y. Lutsyshyn, and J. Boronat, Phys. Rev. B 85, 024101 (2012).
  • (36) C. Cazorla, Y. Lutsyshyn, and J. Boronat, Phys. Rev. B 87, 214522 (2013).
  • (37) R. Rota, Y. Lutsyshyn, C. Cazorla, and J. Boronat, J. Low Temp. Phys. 168, 150 (2012).
  • (38) F. Ercolessi and J. B. Adams, Europhys. Lett. 26, 583 (1994).
  • (39) J. Sala, E. Guàrdia, J. Martí, D. Spångberg, and M. Masia, J. Chem. Phys. 136, 054103 (2012).
  • (40) M. Masia, E. Guàrdia, and P. Nicolini, Int. J. Quantum. Chem. 114, 1036 (2014).
  • (41) J. Nocedal and S. J. Wright, Numerical Optimization (Springer, Berlin, 2006).
  • (42) P. Loubeyre, R. LeToullec, J. P. Pinceaux, H. K. Mao, J. Hu, and R. J. Hemley, Phys. Rev. Lett. 71, 2272 (1993).
  • (43) F. Birch, J. Geophys. Res. 83, 1257 (1978).
  • (44) C. Cazorla, D. Errandonea, and E. Sola, Phys. Rev. B 80, 064105 (2009).
  • (45) R. Barnett, P. Reynolds, and W. A. Lester Jr., J. Comput. Phys. 96, 258 (1991).
  • (46) J. Casulleras and J. Boronat, Phys. Rev. B 52, 3654 (1995).
  • (47) G. Geneste, M. Torrent, F. Bottin, and P. Loubeyre, Phys. Rev. Lett. 109, 155303 (2012).
  • (48) J. M. McMahon, M. A. Morales, C. Pierleoni, and D. M. Ceperley, Rev. Mod. Phys. 84, 1607 (2012).
  • (49) M. A. Morales, J. M. McMahon, C. Pierleoni, and D. M. Ceperley, Phys. Rev. B 87, 184107 (2013).
  • (50) F. Datchi, S. Ninet, M. Gauthier, A. M. Saitta, B. Canny, and F. Decremps, Phys. Rev. B 73, 174111 (2006).
  • (51) C. J. Pickard and R. J. Needs, Nature Materials 7, 775 (2008).
  • (52) I. Errea, B. Rousseau, and A. Bergara, Phys. Rev. Lett. 106, 165501 (2011).
  • (53) Y. Feng, J. Chen, D. Alfè, X-Z. Li, and E. Wang, J. Chem. Phys. 142, 064506 (2015).
  • (54) D. A. Arms, R. S. Shah, and R. O. Simmons, Phys. Rev. B 67, 094303 (2003).
  • (55) G. V. Sin’ko and N. A. Smirnov, J. Phys.: Condens. Matter 14, 6989 (2002).
  • (56) G. Grimvall et al., Rev. Mod. Phys. 84, 945 (2012).
  • (57) E. J. L. Borda, W. Cai, and M. de Koning, Phys. Rev. Lett. 112, 155303 (2014).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description