# Quasiharmonic elastic constants corrected for deviatoric thermal stresses

###### Abstract

The quasiharmonic approximation (QHA), in its simplest form also called the statically constrained (SC) QHA, has been shown to be a straightforward method to compute thermoelastic properties of crystals. Recently we showed that for non-cubic solids SC-QHA calculations develop deviatoric thermal stresses at high temperatures. Relaxation of these stresses leads to a series of corrections to the free energy that may be taken to any desired order, up to self-consistency. Here we show how to correct the elastic constants obtained using the SC-QHA. We exemplify the procedure by correcting to first order the elastic constants of MgSiO-perovskite and MgSiO-post-perovskite, the major phases of the Earth’s lower mantle. We show that this first order correction is quite satisfactory for obtaining the aggregated elastic averages of these minerals and their velocities in the lower mantle. This type of correction is also shown to be applicable to experimental measurements of elastic constants in situations where deviatoric stresses can develop, such as in diamond anvil cells.

###### pacs:

62.20.Dc, 65.40.-b, 91.35.-x, 91.60.Gf## I Introduction

The quasiharmonic approximation (QHA) Wallace ; Anderson is a computationally efficient method for evaluating thermal properties of materials within the density functional theory (DFT) from low to temperatures above the Debye temperature. It provides high quality high pressure-high temperature materials properties karki99 ; RenataPRL ; Sha ; Menendez ; Tsuchiya ; RenataPNAS in a continuous pressure-temperature (PT) domain in which anharmonic effects are negligible.Carrier However, it has a not well recognized shortcoming: the non-hydrostatic nature of thermal stresses in non-isotropic structures. Broadly speaking, these calculations start by obtaining the static internal energy of fully relaxed DFT structures at various pressures. After computations of the vibrational density of states, the thermal energy contribution to the Helmholtz free energy is added. This latter contribution has anisotropic strain gradients and produces deviatoric stresses. This straightforward procedure should be referred to as the statically constrained (SC) QHA. It has been used to compute the elastic constant tensor of isotropic karki99 and non-isotropic minerals RenataPRL ; Tsuchiya at high PT as well, even though pressure conditions were not precisely hydrostatic in the latter calculations. In general, relaxation of deviatoric stresses, irrespective of their origin, is essential in both experimentsBassett ; Menendez and theory,Carrier for generating realistic and reproducible structural and elastic properties.

Here we show how to correct the elastic constant tensor obtained using the SC-QHA. We exemplify
the procedure by correcting to first order the elastic constants of MgSiO-perovskite (PV) and
MgSiO-post-perovskite (PPV), the major phases of the Earth’s lower mantle, for which
elasticity *data* are essential to interpret seismic information of
this region.Lay We
show that this first order correction is quite satisfactory for obtaining the aggregated
elastic averages of these minerals and their acoustic velocities in the PT range of the
lower mantle.

This article is organized as follow: we first discuss the equations used for numerically determining the elastic constant tensor within the SC-QHA. We then describe the procedure for correcting it to first order for deviatoric thermal stresses. We then evaluate these corrections to the previously reported elastic constant tensors of PVRenataPRL and PPV.Tsuchiya

## Ii Elasticity within and beyond the statically constrained (SC) QHA

The present procedure builds on a related procedure to correct structural parameters and equations of state of non-isotropic solids at high PTs.Carrier The method introduced in Ref. Carrier, can correct the SC crystal structure at to infinite order as long as the SC elastic constant tensor is simultaneously corrected. However, this is a very demanding computational procedure and, fortunately, unnecessary. A first order correction to the crystal structure using SC elastic constant, appears to be sufficient. This conclusion was reached after examining the crystal structure of one of the most studied materials at high PT: MgSiO-perovskite.ExpPV This type of experimental data is quite limited and results on other materials with similarly complex crystal structures would be helpful strengthen this conclusion.

According to the (SC) QHA the Helmholtz free energy is given by:

(1) |

where and are respectively Boltzmann’s and Planck’s constants. The first term, , is the volume dependent static energy obtained after full structural relaxation under isotropic pressure, and is the corresponding phonon spectrum. Both phonon spectrum and static energy are here determined using the DFT within the local density approximation (LDA),KohnSham but the methodology is general and applicable to any first principles method. Structural relaxations are performed using a variable cell shape (VCS) algorithmVCS and phonon spectra are computed using the PWscf codequantumespresso as described in Ref. SISSA, , based the linear response theory. The second term in Eq. (1) is the zero point energy, , such that the sum of the terms in the bracket is the energy at =0 K. The last term in Eq. (1) is the thermal excitation energy, (see Ref. Anderson, for details).

Pressure, , is obtained from using the standard thermodynamics relation

(2) |

This procedure implicitly assumes that remains isotropic at all temperatures, but this is only true for static calculations, where structures were optimized at target pressures. The two frequency dependent terms in Eq. (1), the zero-point energy and the thermal energy, contribute to but their strain gradients are intrinsically anisotropic. This effect was recently quantifiedCarrier by the computation of deviatoric thermal stresses, , defined as the difference between the stress tensor and the nominal pressure (diagonal) tensor. In Voigt’s notation:

(3) |

where is the Heaviside step function, equal to 0 for strictly negative and 1 otherwise. Deviatoric thermal stresses are caused by the vibrational (zero-point and thermal) energies and are shown to be important at high pressures and temperatures. The larger the temperature, the more visible these stresses are.

We have previously shown that these deviatoric stresses can be relaxed to first order if one knows the elastic constant tensor, , calculated within the (SC) QHA.Carrier The latter are obtained from the Gibbs free energy, ,

(4) |

by calculating the second derivative of with respect to the strains and :RenataPRL ; Nye

(5) |

The adiabatic elastic constants, which are the relevant ones for interpretation of seismic data, are then computed using appropriate thermodynamics relations.RenataPRL ; Musgrave Below all calculated elastic constants, bulk and shear moduli, and velocities are adiabatic.

Lattice parameters at high pressures and temperatures under hydrostatic conditions can then be corrected to first order by evaluating the strains, , involved in the relaxation of the deviatoric thermal stresses given in Eq. (3):

(6) |

The Cartesian components of the relaxed lattice vectors are then:

(7) |

where

are respectively the matrices of lattice vectors () and Cartesian
strains (keeping up with Voigt’s notation).
Notice that increase in symmetry or symmetry break (phase transformations) may be induced
by deviatoric thermal stresses in the presence of soft phonon, i.e., h and
do *not* necessarily have to the same space group.

In Ref. Carrier, we pointed that attainment of zero deviatoric thermal stresses
within the QHA should involve a *self-consistent* cycle with simultaneous recalculation of
the elastic constant tensor under hydrostatic condition followed by new structural relaxation,
and so on.
However, such procedure is extremely computationally intensive given the need to recompute
vibrational density of states on a PT grid every step of the cycle.
We show next how to obtain the elastic constant tensor corrected to first order with knowledge
of (6) only.

The components of the elastic constant tensor expanded in a Taylor series of strains (in Voigt’s notation) defined by Eq. (6) are:

Neglecting second and higher order terms one has:

(9) | |||||

In the last step above we assumed that pressure is unaffected by shear stresses, i.e., for 4, 5, and 6, thus reducing the index summation from 6 to 3. The stress derivatives of in Eq. (9) are determined using the definition of the pressure as the trace of the stress tensor, . Taking the derivative of the pressure as function of each stress leads to

(10) |

where . This correction requires only knowledge of the pressure derivatives of ’s which are known from the statically constrained QHA calculation, and the deviatoric thermal stresses given by Eq. (3). It gives to first order the elastic constants corrected for deviatoric stresses without having to explicitly calculate Gibbs free energy at the relaxed lattice parameters.

As a final remark, we point that Eq. (10) could also be used and tested
on experimental data as a mean for correcting any type of deviatoric stresses, as long as
the stress deviations remain small compared to the hydrostatic pressure (in a limit for the
Taylor expansion to be valid).
The correction only requires knowledge of (i) the three components ,
, and at the same time (ii) the pressure variation of the elastic constants at specified and
: .
Principal strain deviations, and , are measurable quantities,
for instance, using diffraction ring measurementsBassett
and their corresponding stresses are therefore also available from experiments.
Pressure variation of the elastic constantsDaniels are measurable quantitiesSingh
that require only few additional runs for estimating experimentally the pressure derivative
of at given ’s.
Eventual experimental setting that combines *simultaneously* measurements of (i) and (ii) above
can be used to measure the correction to the elastic constants due to
deviatoric stresses in DAC apparatus, after applying Eq. (10).

## Iii Elastic constants of PV and PPV

We present in this section new results on the deviatoric thermal stresses of PPV and the correction to the elastic constants obtained using the (SC) QHA.Tsuchiya Since deviatoric thermal stresses of PV were recently published,Carrier we also give here the corresponding correction to the elastic constants of PV.

The PT dependent elastic constant tensors of PV and PPV determined using the (SC) QHA have been reported respectively in Ref. RenataPRL, and Ref. Tsuchiya, . These are the major phases of the Earth’s lower mantle and their elastic properties are central information for the interpretation of the seismic properties of this inaccessible region in terms of temperature, composition, and mineralogy. PV and PPV are both orthorhombic crystals respectively with symmetry and . This difference of symmetry group implies in particular, as stated in Ref. Tsuchiya, , that “the , , and directions in the structure correspond to the , , and in the structure, respectively,” corresponding to a rotation of 45 of the — reciprocal lattices. Lattice deformations and deviatoric thermal stresses between PV and PPV are thus comparable only through this transformation. Figure 1(a) shows the deviatoric thermal stresses for PPV. Equivalent results for PV have recently been reported in Ref. Carrier, along with the analysis of its crystalline structure at high PT. The deviatoric stresses and in PPV have opposite sign but similar magnitudes to that of PV (see Ref. Carrier, ), except along the crystalline axes. As stated above, deviatoric thermal stresses for PV and PPV induce distinct deformations along lattices and . The deviatoric thermal stresses in the direction of PPV is considerably larger than the corresponding one in PV leading to larger corrections in PPV than in PV, as shown below. Figure 1(b) shows the percentage corrections to the lattice parameters of PPV, based on Eq. (6). Interestingly, Fig. 1 shows that zero-point energy (the black zero Kelvin line in that figure) also produces deviatoric stresses. With increasing temperature, these stresses are enhanced but their origin is the anisotropic nature of the phonon dispersions.

Figure 2 shows the resulting summation of the three deviatoric thermal stresses [of Fig. 1(a)] for PPV (and see Ref. Carrier, for PV’s deviatoric thermal stresses). It represents the first of the two ingredients necessary for the correction given by Eq. (10). Clearly, the correction for PPV is considerably larger than the one for PV. This is mostly due to that is larger in PPV than in PV (see above). The correction for PPV is always negative, which has the effect of decreasing its elastic constants, while for PV, the correction can be negative (mostly at low temperature) or positive (mostly at high temperature). In principle, there are no reasons for having deviations of systematic nature and they should vary depending on the crystalline structure. One observation that remains true for all crystalline structures, however, is that positive deviations in one direction are to be compensated by a negative deviation in another direction, as observed in both PV and PPV.

Figure 3 shows the pressure derivatives, , of all the elastic constants of PV and PPV, which is the second ingredient required for the correction according to Eq. (10). The figure shows the variations of with pressure for only two temperatures, 0 K and 3000 K, the latter being close to the temperature of the D layer in the lower mantle, where the PPV phase is important in the geophysical models.Lay

Figure 4 shows the corrected bulk and shear moduli, after applying Hill’sVRH (arithmetic) average to the elastic constants, at several temperatures. The corrections are largest at high pressure and high temperature in both PV and PPV. The nature of the correction is also structure-dependent. Notice that the general aspect of the correction to the bulk moduli in Fig. 4 is similar to displayed in Fig. 2, indicating that the dominant term in the correction of Eq. (10) is the deviatoric thermal stress, and to a lesser extent the pressure derivatives of the elastic constants. However, all corrections remain relatively small, meaning the (SC) QHA calculation does not suffer from significant deviatoric thermal stresses, although they can very well be corrected to any level of accuracy.

Table 1 summarizes the corrections to the (SC) QHA for the elastic constants at = 3000K for two pressures, = 100 GPa and = 120 GPa. Corrections are given in parenthesis. Bulk and shear moduli calculated using Voigt (uniform strain), Reuss (uniform stress), and Hill (arithmetic average between Voigt and Reuss) are shown.VRH The volume correction, , as shown in Fig. 1(b), is reported as density, . Velocities are then evaluated from Voigt-Reuss-Hill moduli, since it provides a realistic estimation of the true moduli. Notice that velocities are only slightly modified, because moduli are corrected along with the density; therefore, their ratio remains relatively unaltered.

## Iv Conclusions

In summary, we have introduced a scheme to correct high PT elastic constants obtained using the statically constrained quasiharmonic approximation for deviatoric thermal stresses that develop in calculations of anisotropic structures. This self-consistent scheme was used to compute to first order the elastic constants of the geophysically important MgSiO-perovskite and MgSiO-post-perovskite phases of the lower mantle. The corrections introduced by relaxation of these deviatoric stresses are quite small at relevant conditions of the lower mantle and previous (SC) QHA results remain essentially unchanged. However, this might not be the general case and the current scheme may be used to arbitrary order for computing high PT elastic constants to the desired level of accuracy.

Acknowledgements

This work was supported by NSF grants EAR-0230319, EAR-0635990, and ITR-0428774. We especially thank Shuxia Zhang from the Minnesota Supercomputing Institute for her assistance with optimizing the PWscf code performance on the BladeCenter Linux Cluster and on the SGI Altix XE 1300 Linux Cluster, and Yonggang Yu for helpful discussions relative to PWscf. PC acknowledges partial support from a MSI research scholarship and JFJ from Brazilian agency CNPq.

## References

- (1) D. C. Wallace, Thermodynamics of Crystals (Dover Publications, Mineola, 1972).
- (2) O. L. Anderson, Equations of State of Solids for Geophysics and Ceramic Science (Oxford University Press, New York, 1995).
- (3) B. B. Karki, R. M. Wentzcovitch, S. de Gironcoli, and S. Baroni, Science 286, 1705 (1999).
- (4) R. M. Wentzcovitch, B. B. Karki, M. Cococcioni, and S. de Gironcoli, Phys. Rev. Lett. 92, 018501 (2004).
- (5) X. Sha and R. E. Cohen, Phys. Rev. B 74, 064103 (2006).
- (6) T. Tsuchiya, J. Tsuchiya, K. Umemoto, and R. M. Wentzcovitch, Earth Planet. Sci. Lett. 224, 241 (2004); R. M. Wentzcovitch, T. Tsuchiya, and J. Tsuchiya, Proc. Nat. Acad. Sci. 103, 543 (2006).
- (7) R. M. Wentzcovitch, T. Tshuchiya, and J. Tsuchiya, Proc. Natl. Acad. Sci. USA, 103, 543 (2006).
- (8) E. Menéndez-Proupin and A. K. Singh, Phys. Rev. B 76, 054117 (2007).
- (9) P. Carrier, R. Wentzcovitch, and J. Tsuchiya, Phys. Rev. B 76, 064116 (2007); ibid., 76, 189901 (2007).
- (10) W. A. Bassett, J. Phys. - Condens. Matter 18, S921 (2006).
- (11) K. Hirose, J. Brodholt, T. Lay, and D. A. Yuen, Post-Perovskite: The Last Phase Change, Geophysical Monograph Series, 174, 350 pp. (2008).
- (12) G. Fiquet, D. Andrault, A. Dewaele, T. Charpin, M. Kunz, and D. Hausermann, Phys. Earth Planet. Inter. 105, 21 (1998); N. Funamori, T. Yagi, W. Utsumi, T. Kondo, T. Ushida, and M. Funamori, J. Geophys. Res. 101, 8257 (1996); N. Ross and R. Hazen, Phys. Chem. Miner. 16, 415 (1989); N. Ross and R. Hazen, Phys. Chem. Miner. 17, 228 (1990); W. Utsumi, N. Funamori, T. Yagi, E. Ito, T. Kikegawa, and O. Shimomura, Geophys. Res. Lett. 22, 1005 (1995); Y. Wang, D. Weidner, R. Liebermann, and Y. Zhao, Phys. Earth Planet. Inter. 83, 13 (1994).
- (13) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964); W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
- (14) R. M. Wentzcovitch, Phys. Rev. B 44, 2358 (1991).
- (15) S. Baroni, A. Dal Corso, S. de Gironcoli, P. Giannozzi, C. Cavazzoni, G. Ballabio, S. Scandolo, G. Chiarotti, P. Focher, A. Pasquarello, K. Laasonen, A. Trave, R. Car, N. Marzari, and A. Kokalj, http://www.pwscf.org/
- (16) S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
- (17) J. F. Nye, Physical Properties of Crystals, Their Representation by Tensors and Matrices (Clarendon Press, Oxford, 1985).
- (18) M. J. P. Musgrave, Crystal Acoustics, Introduction to the study of Elastic Waves and Vibrations in Crystals (Holden-Day, San Francisco, 1970).
- (19) W. B. Daniels and C. S. Smith, Phys. Rev. 111, 713 (1958).
- (20) A. K. Singh, H.-k. Mao, J. Shu, and R. J. Hemley, Phys. Rev. Lett. 80, 2157 (1998).
- (21) R. Hill, Proc. Phys. Soc. A 65, 349 (1952).

3000 K, 100 GPa | 3000 K, 120 GPa | |||

PV | PPV | PV | PPV | |

774.8 (0.0) | 933.4 (-3.3) | 844.4 (0.3) | 1069.5 (-2.8) | |

941.7 (0.0) | 756.0 (-2.1) | 1049.9 (0.5) | 846.8 (-1.9) | |

928.5 (0.0) | 949.0 (-2.9) | 1034.9 (0.5) | 1072.5 (-2.6) | |

287.2 (0.0) | 215.8 (-1.7) | 313.6 (0.1) | 286.6 (-1.4) | |

251.0 (0.0) | 164.4 (-1.1) | 265.6 (0.1) | 211.1 (-1.0) | |

248.4 (0.0) | 253.3 (-1.6) | 276.4 (0.1) | 314.9 (-1.2) | |

452.7 (0.0) | 376.7 (-1.3) | 520.4 (0.3) | 433.4 (-1.2) | |

373.9 (0.0) | 370.6 (-1.0) | 421.3 (0.2) | 413.3 (-0.9) | |

406.5 (0.0) | 434.3 (-1.1) | 455.7 (0.2) | 481.1 (-1.0) | |

567.9 (0.0) | 555.7 (-1.7) | 636.0 (0.3) | 627.2 (-1.5) | |

565.2 (0.0) | 553.3 (-1.7) | 632.2 (0.3) | 624.3 (-1.5) | |

562.4 (0.0) | 550.8 (-1.7) | 628.4 (0.3) | 621.5 (-1.5) | |

251.4 (0.0) | 223.8 (-1.2) | 273.3 (0.1) | 273.3 (-1.0) | |

249.2 (0.0) | 219.5 (-1.2) | 270.1 (0.1) | 268.6 (-1.0) | |

247.0 (0.0) | 215.2 (-1.2) | 267.0 (0.1) | 263.9 (-1.0) | |

5.04 (0.00) | 5.11 (-0.01) | 5.22 (0.00) | 5.29 (-0.01) | |

13.35 (0.00) | 12.87 (-0.01) | 13.79 (0.01) | 13.63 (-0.01) | |

7.03 (0.00) | 6.56 (-0.01) | 7.20 (0.00) | 7.13 (-0.01) | |

10.59 (0.00) | 10.41 (-0.01) | 11.01 (0.00) | 10.86 ( 0.00) | |

112.20 (0.02) | 108.33 (-0.10) | 121.21 (0.06) | 118.02 (-0.11) |