CaSiO₃ perovskite at lower mantle pressures

1. Introduction

[2] The lower mantle is generally believed to consist mainly of (Mg,Fe)SiO₃ perovskite and (Mg,Fe)O magnesiowüstite, with CaSiO₃ perovskite, maybe, up to 7–8 volume % [Ringwood, 1975; O'Neill and Jeanloz, 1990; Ita and Stixrude, 1992; Hellfrich and Wood, 2001; Irifune, 1994] or 9 weight % [Kesson et al., 1998]. Despite its importance, there are many unanswered questions about the structure, the stability, the equation of state and the properties of CaSiO₃ under pressure and temperature, that complicates any attempt to model the lower mantle [Stacey and Isaak, 2001].

[3] At ambient conditions, CaSiO₃ forms the mineral wollastonite, which has a pyroxenoid structure. Under pressure, it transforms to walstromite, then to a mixture of larnite (Ca₂SiO₄) and titanite-structured CaSi₂O₅ [Gasparik et al., 1994; Swamy and Dubrovinsky, 1997; Shim et al., 2000a]. The perovskite structure appears at about 10 GPa and 1100 K, the phase boundary having a positive Clapeyron slope of about 2.7 MPa/K (as estimated from the phase diagrams published in the papers cited above). At lower mantle conditions CaSiO₃ has an ideal cubic perovskite structure, while at lower temperatures it is distorted. The small amplitude of the distortions is hardly observable by current high-temperature and high-pressure X-ray techniques and several orthorhombic and tetragonal structures have been proposed, based on experimental measurements [Shim et al., 2002; Kurashina et al., 2004, Ono et al., 2004] or theoretical calculations [Stixrude et al., 1996; Chizmeshya et al., 1996; Akber-Knutson et al., 2002; Magyari-Köpe et al., 2002]. CaSiO₃ perovskite is unquenchable at ambient conditions [e.g., Liu and Ringwood, 1975; Wang and Weidner, 1994].

[4] The equation of state (EOS) of the cubic perovskite phase of CaSiO₃ has been measured up to core-mantle boundary pressures by different groups [e.g., Mao et al., 1989; Tamai and Yagi, 1989; Wang and Weidner, 1994; Wang et al., 1996; Shim et al., 2000a, 2000b, 2002; Kurashina et al., 2004; Ono et al., 2004; Shieh et al., 2004]. The average of the published fits of the experimental results by third order Birch-Murnigham EOS, yields a unit cell volume, V₀ = 45.542 ų, bulk modulus, K₀, ranging from 232 to 288 GPa, and its pressure derivative, K′, within 3.9–4.5. In most of the experimental studies, K′ was fixed to 4.0. A variety of techniques were also employed in theoretical studies of CaSiO₃ [e.g., Wentzcovitch et al., 1995; Chizmeshya et al., 1996; Stixrude et al., 1996; Swamy and Dubrovinsky, 1997; Karki and Crain, 1998; Hama and Suito, 1998; Akber-Knutson et al., 2002], ranging from empirical interatomic potentials fitted on experimental thermodynamical data to pure first-principles calculations. The average of the theoretical results is similar to the one obtained for the experimental data: V₀ = 45.728 ų, K₀ = 246–305 GPa and K′ = 3.5–4.3. Many of these studies, both experimental and theoretical, focus on the behavior of the cubic modification of CaSiO₃ under pressure.

[5] In the present study we perform a detailed investigation of the major symmetry-allowed modifications of CaSiO₃, obtained as distortions from the parent cubic phase. We discuss the stability and the equation of state of these phases and demonstrate that the I4/mcm phase is the most likely stable static atomic configuration.

2. Computational Details

[6] We perform first-principles calculations based on density functional theory (DFT) with plane waves and pseudopotentials. We use Troullier-Martins pseudopotentials, previously tested for Ca [Karki and Wentzcovitch, 2003], and for Si and O [Wentzcovitch et al., 2004] with a 85 Rydberg (1 Rydberg = 13.605 eV) cutoff for the plane waves kinetic energy. We determine the crystal structure using damped variable-cell shape molecular dynamics [Wentzcovitch, 1991]. The cell parameters and atomic positions are fully optimized for initial configurations of various symmetries. The phonon band dispersions are calculated within perturbation DFT [Baroni et al., 2001].

[7] Special attention is paid to the use of equivalent sampling of the Brillouin zone for all the investigated structures. We base our calculation on a primitive 4x4x4 grid of special k points [Monkhorst and Pack, 1976] for the cubic structure. All the lower-symmetry structures have Brillouin zones larger than the cubic structure. The relation between these zones follows the group-subgroup relation that exists between the cubic and each of these structures. Consequently, we perform an accurate mapping of the k points from the cubic parent Brillouin zone to the lower-symmetry Brillouin zones. This guarantees the accuracy and the meaningfulness of the comparison between the different structures.

3. Various Structural Constructions

[8] The phonon dispersion relations of the Pmm cubic phase computed at several pressures show lattice instabilities in the R {1/2, 1/2, 1/2} and M {1/2, 1/2, 0} points of the reciprocal space (Figure 1). All the instabilities correspond to octahedral rotations. The instability in R increases with pressure, while the one in M slightly decreases with pressure. These results are consistent with a previous theoretical investigation [Stixrude et al., 1996]. We build the candidate low-temperature structures by superposing these unstable octahedral rotations into the Pmm cubic structure up to 160 GPa.

[9] The unstable phonon modes in R, three-fold degenerated, correspond to nearly rigid out-of-phase rotations of the SiO₆ octahedra around the cartesian axes. Following the notation proposed by Glazer [Glazer, 1972, 1975], these rotations correspond to a⁻, b⁻ or c⁻, if the rotation axis is parallel to [100], [010] or [001], respectively. The unstable phonon modes in M, one-fold degenerated, correspond to rigid in-phase rotations of the octahedra around the cartesian axis. Their notation is a⁺, b⁺ or c⁺ if the corresponding rotation axis is parallel to [100], [010] or [001], respectively. The absence of the octahedral rotation around a certain axis is denoted by a⁰, b⁰ or c⁰. Rotations along different cartesian axis, but with the same amplitude are denoted with the same letter, for example, a⁺a⁺c⁻ means in-phase rotations along the [100] and [010] cartesian axis, with the same amplitude and an out-of-phase rotation along the [001] cartesian axis with a different amplitude.

[10] Glazer [1972, 1975] obtained all the 26 possible lower-symmetry structures by adding up to three rotations, one around each cartesian axis. Later, Darlington [2002] reduced the number of structures (Table 1), by taking into account only those generated by rotations of same magnitude around equivalent axis. We consider in our analysis the structures generated by the latter approach.

Table 1. Derivation of the Low-Symmetry Structures and the Corresponding Glazer Notation

Space Group	Glazer Notation	Unstable Phonon Modes
Pmm	a0a0a0	-
I4/mcm	a0a0c−	R
Imma	a0b−b−	R,R
Rc	a−a−a−	R,R,R
P4/mbm	a0a0c+	M
I4/mmm	a0b+b+	M,M
Im	a+a+a+	M,M,M
P42/nmc	a+a+c−	M,M,R
Pnma	a−b+a−	R,R,M

4. Results and Discussion

[11] The properties of the lower-symmetry structures are very similar in terms of enthalpy, density and equation-of-state parameters. With the exception of the Rc structure (obtained from the addition of the three unstable R modes), all the other structures have lower energy than the cubic one.

[12] The relative energy differences between all phases are very small, as shown in Figure 2 and Table 2. They are computed based on fits with third-order polynomial expansions of the energy as a function of the compression f [Poirier, 2000]:

Table 2. Equations of State (EOS) and the Relative Energy Differences for All Considered Distortions of CaSiO₃ Perovskite^a

Structure	Pmm	I4/mcm	Imma	Rc	P4/mbm	I4/mmm	Im	P42/nmc	Pnma
BM3
V (Å3)	44.579	44.537	44.567	44.821	44.629	44.599	44.600	44.576	44.576
K (GPa)	250	249	249	247	247	250	250	248	249
K′	4.098	4.090	4.094	4.100	4.124	4.103	4.104	4.092	4.104
BM4
V (Å3)	44.588	44.547	44.576	44.832	44.641	44.609	44.610	44.566	44.588
K (GPa)	248	247	247	244	244	248	247	251	246
K′	4.206	4.213	4.218	4.236	4.261	4.219	4.229	3.977	4.248
K″	−0.002	−0.002	−0.002	−0.002	−0.002	−0.002	−0.002	−0.001	−0.002
Energy, meV/molecule
32.5 Å3	0	−102	−68	483	−10	−10	−5	−88	−68
43 Å3	0	−43	−32	357	10	−8	−8	−37	−32

• a The EOS parameters are derived from third- (BM3) and forth-order (BM4) Birch-Murnagham fits to the first-principles results. The uncertainty of the calculations is on the order of 20 meV/molecule.

[13] In general, the tetragonal structures are energetically preferred to the orthorhombic ones, the lowest in energy being the I4/mcm phase (obtained from one R mode). The energy difference between the cubic and the tetragonal I4/mcm structure is on the order of 43 meV/molecule at low pressure and 100 meV/molecule at high pressure. In order of increasing energy, after the I4/mcm phase we find the P4₂/nmc, Pnma and Imma (whose energy difference is on the order of the accuracy of the calculation), I4/mmm and Im (again with very small energy differences) and the cubic Pmm phase. Our results are different from those obtained considering only the octahedral tilt [Magyari-Köpe et al., 2002] due to the absence of further structural relaxation in the cited study. They are also different with respect to those obtained in a density-functional-based Variationally Induced Breathing model study [Akber-Knutson et al., 2002], who considered only Pm3m, I4/mcm, Pnma and P1 structures and where the orthorhombic one is energetically preferred at high pressure. The differences are perhaps not surprising given the subtlety of the problem and the more empirical approach used in the above study (where the electron density is obtained as a superposition of spherical ionic charge densities, thus not taking into account effects like polarization).

[14] In general our calculated c/a ratio is larger than 1.0, while the experimental measurements [Shim et al., 2002] showed a slightly smaller than 1.0 c/a ratio. With the exception of the P4/mbm structure, the c/a ratio exhibited by the tetragonal and the orthorhombic structures, is about 1.03, slightly increasing with pressure. The P4/mbm structure presents a c/a ratio on the order of 1.08 at 0 GPa, that decreases with pressure to about 1.07 at 160 GPa. The I4/mcm and the P4₂/nmc structures have c/a ratios about 1.02 at 0 GPa and about 1.03 at 160 GPa. The two orthorhombic structures, Imma and Pnma have the b/a ratio about 1.005 at low pressure and about 1.01 at high pressure, and the c/a ratio about 0.995 at all pressures. The I4/mmm structure, exhibits a c/a ratio of about 0.995 at all pressures. A possible alternative explanation to the c/a ratio lower than unity invokes distortions in octahedra instead of octahedral tilting [Shim et al., 2002], resulting in a tetragonal P4/mmm structure, but this has not been considered in our study.

[15] These differences in the c/a ratio between our theoretical results and the experimental data account also for differences in the diffraction peaks (when compared to Shim et al. [2002], Kurashina et al. [2004], and Ono et al. [2004]). Except for the P4/mbm theoretical structure, all the other have hardly observable splittings of the corresponding cubic (200) peak, and with very similar intensities.

[16] All the structures have very similar densities, ranging within about 0.02 g/cm³ at all pressures. The density of the cubic modification is 4.32 g/cm³ at 0GPa and 5.77 g/cm³ at 130 GPa. The Rc structure has the smallest density, while the I4/mcm structure has the largest density, 4.33 g/cm³ at 0GPa and 5.78 g/cm³ at 130 GPa. The density of CaSiO₃ is about 0.25 g/cm³ larger than PREM's [Dziewonski and Anderson, 1981], with a relatively constant difference at all pressures. Obviously, the temperature will decrease these differences.

[17] As expected from the variation of the energy and density with pressure, the pressure-volume relations are similar for all the studied phases. We fit third-order (BM3) and forth-order (BM4) Birch-Murnagham equations of state to first-principles results in all polymorphs. The results of the fits are summarized in Table 2. The values for the specific volume, V₀/Z, where Z is the number of molecules in the unit cell, range in a narrow interval, 44.5–46.0 ų. The bulk modulus within BM3 is about 250 GPa for all the structures, with the exception of Rc, and slightly smaller in BM4, where the spread of the value is somewhat larger: 246–251 GPa.

[18] This study has considered only 2 × 2 × 2 superstructures resulting from unstable phonons at M and R. They all have similar energies and it is anticipated that structures with larger primitive cells but similar octahedral rotation patterns will have similar energies as well. Since along any M-R line an entire optical branch is unstable, it is expected that different and larger superstructures result from other frozen commensurate and/or incommensurate phonons. The resulting structures will have tetragonal symmetry if the frozen phonons reside along the same MR line. Unstable phonons belonging to different MR lines will generate orthorhombic or tetragonal structures.

[19] This multitude of structural instabilities has two main consequences. From the structural point of view, the high temperature stabilization of the cubic structure may occur by two different processes: anharmonic fluctuations between these structures, or static disorder of octahedral rotations. At the moment, without knowledge of the energy barriers between these structures, this remains an open point. This question should be properly investigated by first principles molecular dynamics simulations.

[20] Second, since CaSiO₃ is permanently on the verge of structural transformations between these phases, it is likely to display giant responses to deviatoric or shear stresses. These stresses caused by seismic waves might induce phase transformations (actually symmetry transformations), regardless of temperature. Therefore, the seismic signature accompanying local increase in CaSiO₃ content should be a positive density anomaly and a simultaneous negative anomaly in V_S. The V_P anomaly is less certain due to the relatively large bulk modulus of CaSiO₃ at high pressure (125–145 GPa larger than PREM's at 135 GPa and T = 0 K). Temperature effects on the elasticity of CaSiO3 are expected to be quite dramatic due to these subtle phase transitions and molecular dynamics calculations are expected to provide the magnitude of all these effects. Although still controversial, anti-correlation between density and shear velocity appears to be detected by seismic tomography in places such as beneath the Central Pacific [Ishii and Tromp, 1999]. Excess iron [Humayun et al., 2004] is expected to be related to such anomaly, however, CaSiO3 might be associated with it as well.