Anomalous compressibility of ferropericlase throughout the iron spin crossover
See allHide authors and affiliations

Edited by Hokwang Mao, Carnegie Institution of Washington, Washington, DC, and approved March 17, 2009 (received for review November 29, 2008)
Abstract
The thermoelastic properties of ferropericlase Mg_{1−x}Fe_{x}O (x = 0.1875) throughout the iron hightolow spin crossover have been investigated by first principles at Earth's lower mantle conditions. This crossover has important consequences for elasticity such as an anomalous bulk modulus (K_{S}) reduction. At room temperature the anomaly is somewhat sharp in pressure but broadens with increasing temperature. Along a typical geotherm it occurs across most of the lower mantle with a more significant K_{S} reduction at ≈1,400–1,600 km depth. This anomaly might also cause a reduction in the effective activation energy for diffusion creep and lead to a viscosity minimum in the midlower mantle, in apparent agreement with results from inversion of data related with mantle convection and postglacial rebound.
Understanding of the Earth's lower mantle relies on indirect lines of evidence. Comparison of elastic properties extracted from seismic models with computed or measured elastic properties of candidate minerals at mantle conditions is a fruitful line of enquiry. For instance, it has shed light on the lower mantle composition (1–3) and on the nature of the D″ layer (4, 5). Such comparisons support the notion that the lower mantle consists primarily of ferrosilicate perovskite, Mg_{1−y}Fe_{y}SiO_{3}, and ferropericlase, Mg_{1−x}Fe_{x}O (hereafter, Pv and Fp, respectively). In contrast, evidence based on solar and chondritic abundances suggests a deep lower mantle chemical transition into a pure Pv composition at ≈1,000 km depth (6). A chemical transition with wide topography, gentle, and diffuse changes in elasticity and density is also supported by geodynamic modeling (7). The discovery of the spin crossover in Fp and Pv at lower mantle pressures (8, 9) introduces a new dramatic ingredient that demands a careful reexamination of these phases' elastic properties at appropriate conditions, the consequences for mantle elasticity, and reanalysis of lower mantle properties. This may, after all, support lower mantle models containing a chemical transition. Here, we show the effect of the spin crossover on the bulk modulus and bulk velocity of Fp at high temperatures. We also show the effect it should have on the bulk modulus of a homogeneous lower mantle with pyrolite composition and confirm and justify the origin of anomalies in the elasticity of Fp recently demonstrated at room temperature (10). We point out that such an elastic anomaly might alter the activation energy for diffusion creep (11, 12) in Fp, which might affect mantle viscosity.
Results and Discussions
The highspin (HS) to lowspin (LS) crossover (13) in ferrous iron in Fp has been detected by several techniques at room temperature (8, 10, 14–18) and recently up to 2,000 K (19). For typical mantle compositions the crossover may start as low as ≈35 GPa (18) and end as high as 75 GPa (8) at room temperature. The observed variations in the pressure range of the transition seem to be related to the variable degree of hydrostaticity in experiments. This pressure range broadens substantially with increasing temperature (19). This is actually a crossover that occurs continuously (20, 21) passing through a mixedspin (MS) state. Here, we extend the earlier thermodynamics formalism developed to investigate the spin crossover in Fp (21) by including the spinstatedependent vibrational properties. Equations of state do not have predictive quality unless they include vibrational effects. This is a particularly challenging task given the strongly correlated nature of this solid solution. We then obtain the high temperature compressibility and bulk velocity of Fp. We also address the potential effect the anomalous compressibility across the spin transition might have on the creep viscosity of Fp.
HighTemperature Properties.
Inclusion of the vibrational contribution to the free energy improves considerably agreement between the experimentally measured (18, 19) and our predicted pressure and temperaturedependent spin populations and compression curves at room temperature (see supporting information (SI) Figs. S1 and S2). A description of the calculation and a more detailed analysis of these results are given in SI Text (Figs. S1–S5 and Tables S1–S5). The quality of our predictions can also be tested by inspecting the thermal expansivity, α, shown in Fig. 1. At low or very high pressures α has normal behavior because the system remains, respectively, in pure HS or LS states. The normal thermal expansivity of the HS state is essentially the same as that of MgO (22, 23) as observed experimentally (25). Within the range of validity of the quasiharmonic approximation (QHA) (1, 24), the magnitude of α at 0 GPa agrees very well with measurements for other concentrations (22). Throughout the spin crossover α behaves anomalously. This type of anomaly has been measured in LaCoO_{3} perovskite, another system that undergoes 2 spin crossovers with increasing temperature at 0 GPa (26). The amplitude of the anomaly may be overestimated in our case because of the possibly narrower calculated crossover range. As discussed in the SI Text, the uniform distribution of iron in our calculations underestimates iron–iron interaction, decreases the pressure range of the crossover (18, 60), and increases the magnitude of the anomaly. Nevertheless, this effect is noticeable, could have significant consequences for the mantle geotherm, mantle dynamics, and temperatureinduced lateral heterogeneities. Ultimately a full geodynamic simulation with selfconsistent mineral physics parameters obtained onthefly will probably be necessary to answer these questions. Anomalies on several other thermodynamics quantities will be reported elsewhere (62).
The adiabatic bulk modulus, K_{S}, density, ρ, and bulk velocity, V_{Φ}, along several isotherms are shown in Fig. 2. Below 35 GPa, our calculated K_{S} and ρ are in excellent agreement with the adiabatic bulk modulus, K_{S}, and ρ measured at 300 K in the HS state for X_{Fe} = 0.17 (15). Inclusion of vibrational effects improves considerably the agreement with experiments. The remaining difference is consistent with the difference in iron concentration. There is a considerable reduction in K_{S} and V_{Φ} throughout the spin crossover that is consistent with the reduction in bulk modulus of Fp with X_{Fe} = 0.06 (10) shown in the same figure. The difference in the magnitude of the anomaly is also consistent with the difference in iron concentration, i.e., approximately a factor of 3. The magnitude of the anomaly is more noticeable at low temperatures: at 300 K the crossover pressure range is ≈36–48 GPa compared with the experimental one, ≈35–50 GPa (17), or ≈50–75 GPa (8, 10, 15–17, 19). Therefore, the difference between predictions and measurements are comparable to differences between experiments, but our results agree particularly well with the data of Fei et al. (18) (see Fig. S2).
Potential Effect of the Spin Crossover Transition in Fp on the Mantle Bulk Modulus.
The effect of the spin crossover in Fp along a typical geotherm (28) is shown in Fig. 3 A and B. The anomalies in K_{S} (25 ± 6%) and V_{φ} (15 ± 7%) predicted by a purely elastic model start at ≈40 GPa (≈1,000 km depth) and are most pronounced ≈70 ± 20 GPa (1,600 ± 400 km), but the crossover continues down to the coremantle boundary (CMB) pressure with a possible reentrance into the HS state because of the thermal boundary layer above the CMB (29). However, one should have in mind that anelastic effects might be enhanced in the region where the bulk modulus softens. In contrast, density increases smoothly throughout the entire pressure range of the lower mantle. The shaded areas correspond to possible values of these quantities due to uncertainties in the calculated static transition pressure and the narrower range of our transition pressure.
The net effect of the spin transition in Fp on the bulk modulus of a uniform aggregate with pyrolite composition (31) along a mantle geotherm (29) is shown of Fig. 3C. This comparison is made to elucidate and highlight an effect that may be quite subtle. We have adopted the first principles bulk modulus of Pv reported earlier by our group (1). The effect of iron on the bulk modulus of Pv without the effect of its own spin crossover was included as reported in ref. 32, K(x) = K_{0}(1 + bX_{Fe}), where b varies linearly between 0.079 and 0.044 from 0 GPa to 136 GPa, respectively. Experimentally, Pv's bulk modulus is not noticeably affected by the spin crossover (33). Theory predicts that LS ferrous iron will be displaced from the equilibrium HS site, and that the volume change will be quite insignificant throughout the crossover (34). At 0 GPa the aggregate consists of 80 wt % of Mg_{(1−x)}Fe_{x}SiO_{3}, with x = 0.12, and 20 wt % of Mg_{(1−y)}Fe_{y}O, with y = 0.1875. This translates into a monotonic increase in vol % of Pv in the lower mantle, from 79.6 vol % to 80.8 vol % from 23 GPa to 120 GPa. The bulk modulus of the aggregate was computed by using the Voigt–Reuss–Hill average (35). Compared with PREM's bulk modulus (K_{PREM}) (30), K_{pyr} shows a subtle undulation, i.e., a reduction of ≈4 ± 4%, which appears to be smoothed or cut through by PREM. The uncertainty in K_{pyr} is quite large and permits the signature of the spin crossover in Fp to fall within the uncertainty of global seismic constraints (27). The effect of the spin crossovers in Pv still needs to be better understood and more sensitive strategies need to be devised to identify the signature of this crossover in Fp, which is a subtle one at lower mantle conditions. Nevertheless, given that (i) our predictions should offer an upper bound value for these anomalies, (ii) that this is a purely elastic model (see next section), and (iii) the overall accuracy of the calculation, K_{PREM} does not appear to be inconsistent with K_{S} of a uniform pyrolite aggregate with a spin crossover in Fp along a typical adiabatic geotherm. The softening of the bulk modulus may be more noticeable in colder environments (slabs) and even less noticeable in hotter regions (superplumes).
Correlation Between Mantle Viscosity Structure and the Spin CrossOver in Fp.
The softening of K_{S} in Fp might also have an impact on mantle viscosity. Combination of a thermal convection model by using Newtonian viscous flow and seismic tomography data have implied the existence of a local minimum in mantle viscosity centered at ≈1,500 km (36, 37). We notice the proximity of the viscosity minimum and of the predicted anomaly in the bulk modulus of Fp in the mantle (Fig. 3A). As a relatively minor, weaker phase comprising ≤20 vol % of Earth's lower mantle, the influence of Fp on viscosity depends critically on its distribution. In a poorly mixed system, Fp grains will be isolated from one another by Pv grains, which have a viscosity ≈10^{3} times that of Fp far from the spin crossover (25). With Pv forming a loadbearing framework, the effect of Fp on viscosity will be modest. However, if phase separation occurs during largestrain deformation, Fp will markedly impact lower mantle viscosity. Recent shear deformation experiments on partially molten rocks, as well as on 2phase rocks in which the viscosities of the 2 phases are significantly different, demonstrate a profound segregation of the constituent phases (28, 38). Mineralogical segregation and compositional layering are also observed in highly strained, naturally deformed rocks (39, 61). Bands rich in Fp, separated by regions rich in Pv, are thus anticipated in a deforming lower mantle. Once phase separation occurs, strain localizes in the weak, Fprich layers causing a significant decrease in viscosity relative to the viscosity of a homogenously mixed, 2phase rock (41).
Here, we invoke an elastic strain energy model (ESEM) (11) for viscosity to estimate the potential impact of the bulk modulus softening on Fp's viscosity, η_{Fp}. A Newtonian subsolidus flow is assumed consistent with a diffusion creep deformation mode expected in the mantle and with the model used to infer lower mantle viscosity on the basis of convectionrelated and postglacial rebound data (36, 37). Fp's viscosity, η_{Fp}, is then: where G* is the extrinsic activation energy for the dominant deformation mechanism, i.e., ionic diffusion, and A is a constant. At a depth z, η_{Fp}(z), should be (12): where η_{Fp}(z_{0}) and G*(z_{0}) are Fp's viscosity and activation energy at a reference depth z_{0}, here assumed to be the top of the lower mantle.
The ESEM relates the activation energy for diffusion, G*(z_{0}), with the shear and bulk modulus of the system. The ionic diffusion process induces bond stretching and/or shearing depending on the diffusion path. As such, the diffusion barrier is related to different extents to shear and bulk modulus. This is usually described as a parameterized dependence on the pure shear and dilatational contributions, G*_{s}(z) and G*_{D}(z), to the activation energy, where δ is a free parameter. The other quantities are (12): with μ(z), K(z), and V(z) being shear and bulk moduli and volume, respectively. This model works well for metals, but the relationship between the diffusion barrier and the elastic moduli for ionic systems may not be this simple, even though there are indications that this model describes well the highpressure and hightemperature behavior of diffusion in MgO (42). Nevertheless, this model expresses a relationship that is very likely to exist in some similar form between viscosity and elastic moduli. Despite consistency between experimental data (42) and first principles results of migration barriers in MgO (43–45), similar investigations in Fp are still necessary to clarify this point. Much less is known about the shear modulus at this point. Room temperature measurements (10) have indicated that the shear modulus also softens throughout the spin crossover, but this has not been confirmed by theory or by more recent Brillouin scattering data (46). Therefore, the situation remains controversial and shear deformation may enhance or damp the bulk modulus related viscosity anomaly. Experimental data (42) and modeling (43–45) have suggested that G*(z_{0}) ≈ 300–330 kJ/mol at uppermost lower mantle conditions (z_{0} = 660 km, P = 23 GPa, T = 2,000 K). We then assume G*(z_{0}) = 315 kJ/mol (42). Even if the shear modulus were known with certainty at this point, one would still have to estimate δ to infer the impact of the shear modulus on Fp's viscosity.
The impact of the bulk modulus softening on Fp's viscosity predicted by a purely dilatational ESEM (δ = 0 in Eq. 11) is shown in Fig. 4 compared with the relative changes in lower mantle viscosity, η(z) (36, 37), with depth. All profiles have accentuated minima at ≈1,400–1,600 km. The decrease in Fp's viscosity near the CMB in our model is caused by the reentrance into the HS state owing to the thermal boundary layer (29) above the CMB, whereas the more drastic reduction in mantle viscosity beyond 2,000 km may be related with numerous additional factors (40), such as the approaching postperovskite transition or the temperature profile. It appears to depend also on the inversion model used to obtain the viscosity (36, 37).
The bulk modulus anomaly in Fp may not only affect the viscosity and dynamics of the mantle (49) but also its overall state (50) and properties. In general, it is anticipated that properties of Fp related with ionic diffusion, such as ionic conductivity, should improve in the MS state owing to its enhanced compressibility (anomalously “soft” bonds), even though ionic conductivity is not the prevailing electrical conduction mechanism at conditions explored so far (51–53). In contrast, such properties should deteriorate in the LS state compared with the HS state because of the reduction in lattice parameter. Heat (lattice) conductivity, instead, is expected to follow the opposite trend: it should be boosted in the LS state and damped in the MS state in comparison with the HS state. Seismic attenuation also results from an activated diffusion process. Experimentally, the attenuation Qfactor is related to viscosity as Q ∝ η^{β}, where 0.2 < β < 0.9 (54). Given this relation and the possibility of viscosity reduction by a factor of ≈10–50 in the middle of the crossover, one would normally expect Q^{−1} in Fp to increase by a factor of ≈10–30 in the same regime. However, it is unclear whether this relationship between Qfactor and viscosity holds for materials undergoing spin crossover, but in principle one should also expect enhanced attenuation throughout the crossover. Several aspects of the spin crossover in Fp still need to be investigated before its consequences are better understood.
Methods
Thermodynamics of the CrossOver Transition.
We treat Fp in the MS state as an ideal solid solution of HS and LS states. This approximation seems to be well justified by the concentrationindependent static spin transition pressure for concentrations up to x = 0.1875 (21). Therefore: where n = n(P,T) is the LS fraction, and V_{LS}, V_{HS}, K_{LS}, and K_{HS} are the equilibrium volume and isothermal compressibility of pure LS and HS states. Eq. 6 differs from the weighted average of the compressibilities by an additional term caused by the pressure dependence of n(P,T). This last term is ultimately responsible for the bulk modulus anomaly reported recently (10). According to Eqs. 5 and 6, the properties of Fp in the MS state may be determined from those of the LS and HS states plus the LS fraction, n(P,T), all of which must be computed by first principles.
In contrast to the previous thermodynamics treatment (21) we now include vibrational effects. It is impossible to address thermodynamics properties without them. The other approximations used in ref. 21 to compute n(P,T) are retained. They are:

The magnetic entropies are S_{HS}^{mag}(n) = k_{B}X_{Fe} ln[m(2S + 1)] and S_{LS}^{mag} = 0 for the HS and LS states, respectively. S and m are, respectively, the total spin quantum number and orbital degeneracies of the HS (S = 2 and m = 3) and LS (S = 0 and m = 1) states.

The HSLS configuration entropy is S_{conf} = −k_{B}X_{Fe} [n ln n + (1 − n) ln(1 − n)]. Fluctuations in n(P,T) are insignificant given the finite sample sizes. Because configurations are not expected to be static in this solid solution, this formula implicitly assumes the ergodic hypothesis, i.e., time and ensemble averages are equal.

3The Mg/Fe configuration entropy is insensitive to spin state.
n(P,T) is then obtained by minimizing the Gibbs free energy with respect to n. This leads to: where ΔG_{LS−HS}^{stat+vib}(P, T) is the difference between the static plus vibrational contributions to the free energy of the LS and HS states, X_{Fe} is the concentration of iron (here, 0.1875). Therefore, to obtain ΔG_{LS−HS}^{stat+vib}(P, T) and n(P,T) it is necessary first to obtain the vibrational spectrum and free energies of pure spin states within the QHA.
Vibrational Virtual Crystal Model (VVCM).
The thermal properties of Fp in pure spin states were computed by using the QHA (55) in which the Helmholtz free energy is given by: where U(V) is the volumedependent static total internal energy obtained by first principles and ω_{qj}(V) is the corresponding volumedependent phonon spectrum.
Current methodological limitations preclude a direct computation of the vibrational density of states (VDOS) of Fp within the first principles LDA+U approach (56). To circumvent this problem we developed a VVCM. The VC concept involves the replacement the atomic species forming the solid solution, in this case, magnesium and HS or LS irons, by an “average cation” that can reproduce the properties of the solid solution. Here, we develop a VC to compute only vibrational and thermodynamics properties. We are not aware of previous similar attempts in the literature. The development of successful VVCMs would be extremely useful to bypass the difficult problem of computing VDOS for numerous configurations involving hundreds of atoms representative of solid solutions, especially strongly correlated ones, so common in minerals.
The VVCMs corresponding to the pure HS and LS states consist of 2 atoms per cell in the rocksalt structure: oxygen and a virtual (cation) atomic species with a mass where M_{Mg} and M_{Fe} are, respectively, the atomic masses of magnesium and iron, with the latter being independent of the iron's spin state. The interactions of the VC cation in the solid are modeled for the purpose of computing vibrational and thermodynamics properties only.
The VVCMs are essentially periclase, MgO with modified interatomic force constants that reproduce the elastic constants of HS or LS Fp and cation masses as in Eq. 9. The force constants of periclase were previously computed by first principles and produce excellent phonon dispersions (57, 58). The force constants of the HS or LS VCs are obtained by matching the elastic constants extracted from the acoustic phonon dispersions (57, 58)^{†} close to k = 0 to the elastic constants of HS and LS calculated by first principles. There is a linear relationship between force constants D_{μν}(R^{ij}) and elastic constants C_{σταβ} (59): Here, Greek letters refer to Cartesian indices, C_{σταβ}(V) are the volumedependent elastic constants in cartesian notation, while D_{μν}^{ij} are the interatomic force constants between atoms i and j separated by R^{ij} when displaced in directions μ and ν, respectively. The sum in Eq. 10 is over all atomic pairs (i,j) and a_{μν,σταβ}^{ij}(V) are a set of volumedependent constants. Because of symmetry constrains, many of the a_{μν,σταβ}^{ij} constants vanish. Eq. 10 is a convergent summation since the force constants vanish rapidly with the interatomic distances. The convergence in Eq. 10 is guaranteed if the force constants vanish faster than 1/R^{5}, where R is the interatomic separation (59).
The force constants defined as are used to compute the phonon spectrum at each volume: and we need to obtain D_{μν}^{ij}(V) for HS and LS VVCMs. Fp and periclase have only 3 elastic constants, C_{11}, C_{12}, and C_{44} (Voigt notation) (35). We may modify 3 force constants of periclase independently to reproduce the static elastic constants of HS and LS Fp. We modified the 3 largest interatomic force constants of periclase, D_{xx}^{12} (Mg–O nearest neighbor longitudinal interaction), D_{xy}^{11} (Mg–Mg nearest magnesium interaction), and D_{xy}^{12} (the Mg–O nearest neighbor transverse interaction). All other force constants of MgO are at least 1–2 orders of magnitude smaller and remained unchanged. Changes in those force constants have only a minor effect on the elastic constants. More details of the VVCM developed here will be discussed elsewhere (62).
A comparison between the static bulk modulus obtained by fitting an equation of state to the energy versus volume relation in HS and LS Mg_{1−x}Fe_{x}O (x = 0.1875) and the bulk modulus obtained from the elastic constants of the respective VCs is shown in Fig. S3. The virtual crystals produce distinct vibrational density of states (VDOS) for periclase, HS, and LS Fp (see Fig. S4). The acoustic mode dispersions of the HS and LS VVCMs are precisely the same as those of HS and LS Fp. This ensures that thermodynamics calculations are carried out with the correct VDOS at low frequencies, which matter the most, and with a reasonably good weightaveraged VDOS at high frequencies as well. The VVCM should offer more accurate thermodynamics properties than a Debyelike model because of the more detailed structure of the VDOS. We considered carefully the pressure/temperature range of validity of the QHA. Full and dashed lines in all figures correspond to conditions within and outside its range of validity, respectively. The upper temperature limit of the QHA is adopted as the lowest temperature of the inflection points in the thermal expansivity of pure LS and HS Fp at every pressure (1, 24), i.e., ∂^{2}α/∂T^{2}_{P} = 0.
Acknowledgments
R.M.W. thanks the Japan Society for the Promotion of Science for support and the Department of Earth and Planetary Sciences of Tokyo Institute of Technology for hospitality during preparation of this manuscript. Calculations were performed with the Quantum ESPRESSO package at the Minnesota Supercomputing Institute and on the Big Red Cluster at Indiana University. This work was supported by National Science Foundation (NSF)/Division of Earth Sciences Grant 0635990, NSF/Division of Atmospheric Sciences Grant 0428774 (Virtual Laboratory for Earth and Planetary Materials), NSF/Division of Materials Research Grant 0325218 (Institute for the Theory of Advanced Materials in Information Technology), NSF/Division of Ocean Sciences, and University of MinnesotaMaterials Research Science and Engineering Center.
Footnotes
 ^{1}To whom correspondence should be addressed. Email: wentzcov{at}cems.umn.edu

Author contributions: R.M.W., J.F.J., Z.W., D.A.Y., and D.K. designed research; R.M.W., J.F.J., Z.W., and C.R.S.d.S. performed research; R.M.W., J.F.J., Z.W., and C.R.S.d.S. analyzed data; and R.M.W., J.F.J., Z.W., D.A.Y., and D.K. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/cgi/content/full/0812150106/DCSupplemental.

↵† The static elastic constants of Fp in pure spin states were computed by first principles by using the usual stress versus strain relations. The lower symmetry of the solid solution is restored by resymmetrization of the elastic constant tensor as described in ref. 32.

Freely available online through the PNAS open access option.
References
 ↵
 Wentzcovitch RM,
 et al.
 ↵
 ↵
 ↵
 Wentzcovitch RM,
 Tsuchiya T,
 Tshuchiya J
 ↵
 ↵
 van der Hilst RD,
 Bass JD,
 Matas J,
 Trampert J
 Williams Q,
 Knittle E
 ↵
 Kellogg LH,
 Hager BH,
 van der Hilst RD
 ↵
 Badro J,
 et al.
 ↵
 Badro J,
 et al.
 ↵
 Crowhurst JC,
 et al.
 ↵
 ↵
 Ellsworth K,
 Schubert G,
 Sammis CG
 ↵
 Gutlich P,
 Goodwin HA
 ↵
 ↵
 ↵
 Speziale S,
 et al.
 ↵
 Lin JF,
 et al.
 ↵
 Fei Y,
 et al.
 ↵
 Lin JF,
 et al.
 ↵
 Sturhahn W,
 Jackson JM,
 Lin JF
 ↵
 ↵
 Touloukian YS,
 et al.
 ↵
 van Westrenen W,
 et al.
 ↵
 Carrier P,
 Wentzcovitch RM,
 Tsuchiya J
 ↵
 Yamazaki D,
 Karato Si
 ↵
 Radaelli PG,
 Cheong SW
 ↵
 Masters G
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Holtzman BK,
 et al.
 ↵
 ↵
 ↵
 Holtzman BK,
 Kohlstedt DL,
 Morgan JP
 ↵
 ↵
 ↵
 Karki B,
 Khanduja G
 ↵
 ↵
 Wang J,
 SanchezValle C,
 Bass J
 ↵
 Grand SP,
 van der Hilst RD,
 Widyantoro S
 ↵
 ↵
 Boschi E,
 Ekstrom G,
 Morelli A
 Yuen DA,
 et al.
 ↵
 ↵
 ↵
 ↵
 ↵
 Karato SI
 ↵
 Wallace D
 ↵
 Cococcioni M,
 de Gironcoli S
 ↵
 ↵
 Karki BB,
 Wentzcovitch RM,
 de Gironcoli S,
 Baroni S
 ↵
 Ashcroft NW,
 Mermin ND
 ↵
 ↵
 Larkin LM,
 Zimmerman ME,
 Kohlstedt DL
 ↵
 Wu Z,
 et al.
Citation Manager Formats
Article Classifications
 Physical Sciences
 Geophysics