A Metastable Fo-III Wedge in Cold Slabs Subducted to the Lower Part of the Mantle Transition Zone: A Hypothesis Based on First-Principles Simulations

The metastable olivine (Ol) wedge hypothesis assumes that Ol may exist as a metastable phase at the P conditions of the mantle transition zone (MTZ) and even deeper regions due to inhibition of the phase transitions from Ol to wadsleyite and ringwoodite caused by low T in the cold subducting slabs. It is commonly invoked to account for the stagnation of the descending slabs, deep focus earthquakes and other geophysical observations. In the last few years, several new structures with the forsterite (Fo) composition, namely Fo-II, Fo-III and Fo-IV, were either experimentally observed or theoretically predicted at very low T conditions. They may have important impacts on the metastable Ol wedge hypothesis. By performing first-principles calculations, we have systematically examined their crystallographic characteristics, elastic properties and dynamic stabilities from 0 to 100 GPa, and identified the Fo-III phase as the most likely metastable phase to occur in the cold slabs subducted to the depths equivalent to the lower part of the MTZ (below the ~600 km depth) and even the lower mantle. As disclosed by our theoretical simulations, the Fo-III phase is a post-spinel phase (space group Cmc21), has all cations in sixfold coordination at P < ~60 GPa, and shows dynamic stability for the entire P range from 0 to 100 GPa. Further, our static enthalpy calculations have suggested that the Fo-III phase may directly form from the Fo material at ~22 GPa (0 K), and our high-T phase relation calculations have located the Fo/Fo-III phase boundary at ~23.75 GPa (room T) with an averaged Clapeyron slope of ~−1.1 MPa/K for the T interval from 300 to 1800 K. All these calculated phase transition pressures are likely overestimated by ~3 GPa because of the GGA method used in this study. The discrepancy between our predicted phase transition P and the experimental observation (~58 GPa at 300 K) can be explained by slow reaction rate and short experimental durations. Taking into account the P-T conditions in the cold downgoing slabs, we therefore propose that the Fo-III phase, rather than the Ol, highly possibly occurs as the metastable phase in the cold slabs subducted to the P conditions of the lower part of the MTZ (below the ~600 km depth) and even the lower mantle. In addition, our calculation has showed that the Fo-III phase has higher bulk seismic velocity, and thus may make important contributions to the high seismic speeds observed in the cold slabs stagnated near the upper mantle-lower mantle boundary. Future seismic studies may discriminate the effects of the Fo-III phase and the low T. Surprisingly, the Fo-III phase will speed up, rather than slow down, the subducting process of the cold slabs, if it metastably forms from the Ol. In general, the Fo-III phase has a higher density than the warm MTZ, but has a lower density than the lower mantle, as suggested by our calculations.

However, some new high-P polymorphs, compositionally identical to but structurally different from Ol, have been discovered in the last decade or so, which thus introduces new variables to the metastable Ol wedge hypothesis. These new metastable high-P polymorphs were discovered in meteorites [35,36], in in-situ X-ray diffraction single-crystal compression experiments [37,38], in polycrystalline shock-wave compression experiments [39] and in first-principles simulation experiments [40]. Using micro-Raman spectroscopy and high-resolution transmission electron microscopy combined with molecular dynamics simulations, Van de Moortèle et al. [35] reported a new polymorph (ζ-(Mg, Fe) 2 SiO 4 ) in some shocked Ol grains of the Martian meteorites NWA 2737 and NWA 1950, but did not document the structural details. In their compression experiments with forsterite (Fo) single crystals at room T, Finkelstein et al. [37] discovered two new high-P polymorphs, forsterite-II at 50 GPa and forsterite-III at~58 GPa with the space group of P1 and Cmc2 1 , respectively (denoted as Fo-II and Fo-III hereafter in this study). Their meta-dynamics simulation showed a stepwise phase transition sequence from Fo to Fo-II and then to Fo-III via a shear mechanism. Using first-principles molecular dynamic calculations, Zhou et al. [40] found another polymorph with the P2 1 2 1 2 space group, which had lower enthalpy than Fo at P > 43. 46 GPa and 0 K (named as Fo-IV hereafter in this study). In their compression experiments with a San Carlos Ol single crystal, Santamaria-Perez et al. [38] observed some new "defect" peaks associating with the increase of the Si coordination and the amorphization process when P increased further, which was supported by their molecular dynamics simulations. While studying a heavily shocked meteorite, Tomioka and Okuchi [36] found a new high-P polymorph, the ε*-phase, occurring as nanoscale lamellas and having a topotaxial relationship with the host mineral Rwd. Consequently, they proposed that the ε*-phase could form from all other polymorphs via a shear transformation mechanism, and that such diffusion-free processes could complicate the high-P transformation of Ol and occur not only in shocked meteorites but also in lithospheres subducting into the deep Earth. In total, these studies reported five new Ol polymorphs, three of which had detailed structural information valuable for further first-principles calculations (Fo-II, Fo-III and Fo-IV). Because the phase transition of Ol to its high-P polymorphs may always involve several different processes such as nucleation and growth along the grain boundary where atomic diffusion plays a very important role, and lattice-coherent shear mechanism [1,3,6,22,26,[41][42][43], other kinetically favorable (lower activation/nucleation/growth barriers) phase transition paths where Ol transforms into some metastable polymorphs under high P and relatively low T, rather than directly transforms to the high-P stable phases Wds and Rwd (high activation barriers), are certainly possible. As such, these thermodynamically metastable phase transition paths in the inner cold part of a subducting slab could not complicate the phase transformation sequence of Ol only, but also shed new lights on the geodynamic process of the subducted materials in the deep Earth.
In this study, we have conducted extensive first-principles calculations to systematically examine the crystallographic characteristics, elastic properties and dynamic stabilities of these three newly discovered high-P polymorphs, namely Fo-II, Fo-III and Fo-IV, and their phases relations to other common mantle phase or phase assemblages (Fo, Wds, Rwd, akimotoite (Aki) + periclase (Pe), and bridgmanite (Brg) + Pe) at high P and high T. We did not consider the stishovite + Pe assemblage or the majorite phase, because the former is stable only at very low T (<125 K) for the composition Mg 2 SiO 4 and the latter is a low-P (<~22 GPa) and high-T (>1500 K) phase [44]. Our phonon calculation results show that the Fo-II and Fo-IV phases are dynamically unstable, but the Fo-III phase has a distinct dynamic stability for a vast P range. The calculated Fo/Fo-III phase boundary, with a transition pressure of 23.75 GPa at room T and a negative Clapeyron slope of~−1.1 MPa/K, is very similar to the Rwd/Brg + Pe phase boundary. This calculated phase transition pressure is quite different to the experimentally determined value (~58 GPa at room T; Finkelstein et al. [37]), which was presumably severely affected by a slow reaction rate of the Fo/Fo-III phase transition and short experimental durations. Further, combining the results of previous studies and our calculations, we infer that the Fo/Fo-III phase transition with a shear transformation mechanism is quite possible in the inner part of a slab during its fast subduction to the 660 km depth. Consequently, the Fo-III metastable phase may act as a 'bridge' phase linking the low-P Ol with the high-P stable phase or phase assemblages such as the Rwd, Brg + Pe and Aki + Pe, and play an important role in governing the dynamic interactions between the slabs and the ambient mantle.

Methods
All the computational work reported here was done on the Tianhe-2 Supercomputer System hosted by the National Supercomputer Center in Guangzhou, China. It was performed with the Vienna Ab-initio Simulation Package (VASP; Kresse and Furthmüller [45]) using density functional theory (DFT; Hohenberg and Kohn [46]; Kohn and Sham [47]). The projector augmented wave pseudopotentials [48] were used to describe the interactions between the ions and electrons. The electronic configurations used in our studied systems were Mg (3s 2 3p 0 ), Si (3s 2 3p 2 ), and O (2s 2 2p 4 ). The generalized gradient approximation (GGA; Perdew et al. [49]) with the PW91 functional [50], a technique well-known in reproducing phase relations (the primary interest of this study), was used to approximate the exchange-correlation interaction. In some particular cases, the local density approximation (LDA; Ceperley and Alder [51]; Perdew and Zunger [52]) method was also employed (see later discussion).
The static geometry optimization and ionic relaxation under different hydrostatic P were performed using the conjugate algorithm [53], and both the lattice parameters and the atomic positions were allowed to change during the optimization. All the initial structures used were adopted from previous studies (Fo: Hazen [54]; Wds: Horiuchi and Sawamoto [55]; Rwd: Hazen et al. [56]; Brg: Horiuchi et al. [57]; Aki: Horiuchi et al. [58]; Pe: Hazen [59]; Fo-II and Fo-III: Finkelstein et al. [37]; Fo-IV: Zhou et al. [40]). The planewave cutoff energy was set as 650 eV for all these structures for consistency, and the energy convergence criterion was set as 10 −8 eV. The k-point sampling grid [60] for each phase is listed in Table 1. The cutoff energy and k-point sampling grid were tested with the total energy converging to 1 meV/atom. This observation was in good agreement with Hernandez et al. [44], in which Fo, Wds, Rwd, Brg, Aki, Pe, stishovite and majorite were successfully simulated with a lower planewave cutoff energy of 500 eV. Likewise, it is believed that our higher planewave cutoff energy should be sufficient to converging the atomic forces, a requirement for the phonon calculations, which will be addressed below. Based on the optimized structures, the lattice vibrational properties were examined by calculating the dynamic matrix using the density functional perturbation theory (DFPT; Gonze and Lee [61]; Baroni et al. [62]) with a unit/primitive cell rather than a supercell, which could save a lot of computing time but lost little accuracy. Using the calculated lattice vibrational information at the Γ point in the first Brillouin zone, the phonon dispersion curve and corresponding density of state were calculated with the PHONOPY code [63], and were in turn used to calculate the thermodynamic quantities such as the equilibrium volume (V), isothermal bulk modulus (K T ) and its first pressure derivative (K T ), adiabatic bulk modulus (K S ), Gibbs free energy (G), thermal expansion coefficient (α) and thermal Grüneisen parameter (γ th ) at certain P and T with a quasi-harmonic approximation method (QHA; Born and Huang [64]; Carrier et al. [65]; Wentzcovitch et al. [66]). The q-point sampling grid for the QHA calculations for each phase was tested for the Gibbs free energy converging to 0.01 kJ/mol, and is listed in Table 1.

Crystallographic Characteristics of Fo, Fo-II, Fo-III and Fo-IV
The DFT optimized structures of the Fo, Fo-II, Fo-III and Fo-IV at 0 K and some selected P are plotted in Figure 1, using the VESTA program [67]. The lattice parameters of these phases at 0 K and 0 GPa are compared in Table 2. Their crystallographic details at some selected P are summarized in Supplementary Tables S1-S8. Essentially identical to the experimental results [54], all the O atoms in the Fo phase form an expanded and distorted hexagonally close-packed array stacked along the a direction whereas the Si atoms locate in the isolated tetrahedra which share corners with the [MgO 6 ] octahedra ( Figure 1a) results using the GGA technique [66,68], our optimized unit-cell parameters a, b, c, and V at 0 GPa and 0 K are slightly larger than those experimentally reported values (Hazen [54]; 1 atm and 296 K) bỹ 0.73%, 0.79%, 0.70% and 2.23%, respectively (Table 3).  [66,68], our optimized unit-cell parameters a, b, c, and V at 0 GPa and 0 K are slightly larger than those experimentally reported values (Hazen [54]; 1 atm and 296 K) by ~0.73%, 0.79%, 0.70% and 2.23%, respectively (Table 3). The Fo-II phase was experimentally determined as a triclinic phase (space group P1, Z = 4), with no crystallographic details refined from the single-crystal XRD data though (Finkelstein et al. [37]; ~52.4 GPa). Our optimized structure at 0 GPa shows that there are two types of Si atoms, Si(I)  The Fo-II phase was experimentally determined as a triclinic phase (space group P1, Z = 4), with no crystallographic details refined from the single-crystal XRD data though (Finkelstein et al. [37]; 52.4 GPa). Our optimized structure at 0 GPa shows that there are two types of Si atoms, Si(I) and Si(II), in equal proportions (see Supplementary Table S2). The Si(I) atoms locate in tetrahedral sites while the Si(II) atoms locate in octahedral sites (Figure 1b), with one [Si(I)O 4 ] tetrahedron pairing with one [Si(II)O 6 ] octahedron (via corner-sharing) to form one structural unit isolated to other Si-O polyhedra. As for the Mg atoms, they distribute on five crystallographically different sites (from Mg(I) to Mg(V) with the ratio of 2:2:1:2:1), and all form octahedra (see Supplementary Tables S2 and S3). P brings forth some changes in the crystal structure. Our optimized structure at 60 GPa (Figure 1c) shows that although the Si cations maintain as either fourfold-coordinated (50%) or sixfold-coordinated (50%), the Mg cations have divided into two groups, with one quarter of them in sevenfold coordination (forming distorted pentagonal bipyramid) and the rest in sixfold coordination, in good agreement with the simulated results of Finkelstein et al. [37]. The optimized lattice parameters at 52.4 GPa and 0 K are compared in Table 3 with the experimentally obtained ones at the same pressure but 298 K. The agreement between these two sets of lattice parameters is good, with our calculated values slightly larger (a by 0.58%, b by 0.56%, c by 0.31%, and V by 1.32%), a common phenomenon for the GGA method [66,68].
The Fo-III phase was experimentally determined as a non-centrosymmetric version of the Cmcm CaTi 2 O 4 -type structure [37], and its crystallographic details were well-refined by the single-crystal XRD method. The experimentally determined lattice parameters at 58.2 GPa (298 K) are compared in Table 3 with our calculated results at the same pressure (0 K). Clearly, our optimization overestimates the b parameter by 1.86%, but abnormally underestimates the a, c and V parameters by 1.02%, 2.15% and 1.28%, respectively. These unexpected minor discrepancies may be attributable to the experimentally observed high defect density in the Fo-III crystal under high P, which might have affected the single-crystal XRD result [37]. At 0 GPa, the optimized structure has only one type of Si atom and two types of Mg atom, Mg(I) and Mg(II). All the Si atoms locate in octahedral sites, forming slightly distorted [ It should be pointed out that all these chains run along the a direction, so that a distinct elastic anisotropy should be expected for the Fo-III phase. Compared to our optimized structure at 0 GPa, our optimized structure at 60 GPa shows that the coordination of the Mg(II) has become sevenfold, in good agreement with the results obtained at 58.2 GPa by Finkelstein et al. [37]. Higher P at least up to~100 GPa produces no further major structural change.
The Fo-IV phase was a theoretical discovery (space group P2 1 2 1 2; Zhou et al. [40]). Our optimized lattice parameters for 0 GPa are compared in Table 3 with the results of Zhou et al. [40]. Since we used the GGA method and Zhou et al. [40] used the LDA method, our calculated unit-cell parameters are all slightly larger than theirs (a by 1.28%, b by 1.30%, c by 1.51% and V by 4.15%), as expected. In our optimized structure at 0 GPa, there are one type of Si atom, and two types of Mg atom, Mg(I) and Mg(II), in equal proportions (see Supplementary Tables S6-S8 Figure 1f). Higher P leads to some structural variations: at 20 GPa, the coordination of the Mg(II) has changed from fourfold to fivefold; at 40 GPa, it has further changed to sixfold.
The P dependence of the crystallographic characteristics of the Fo, Fo-II, Fo-III and Fo-IV phases at 0 K has been characterized by two kinds of parameters: one is the average bond-length (d mean ) of the Mg-O and Si-O polyhedra, and the other is the effective coordination number (ECoN). Both parameters were calculated using the equations provided by Hoppe [69], as implemented in the VESTA code. The results are shown in . Figure 2 shows the results for the Fo phase. As P increases from 0 to 100  (dmean) of the Mg-O and Si-O polyhedra, and the other is the effective coordination number (ECoN). Both parameters were calculated using the equations provided by Hoppe [69], as implemented in the VESTA code. The results are shown in Figures 2-5. Figure 2 shows the results for the Fo phase. As P increases from 0 to 100 GPa, the dmean parameters for both the [SiO4] tetrahedron and the [MgO6] octahedron (solid curves) decrease smoothly, and the ECoN parameters remain nearly constant (dashed curves). Overall, the microstructural features of the Fo phase do not change suddenly.     Both parameters were calculated using the equations provided by Hoppe [69], as implemented in the VESTA code. The results are shown in Figures 2-5. Figure 2 shows the results for the Fo phase. As P increases from 0 to 100 GPa, the dmean parameters for both the [SiO4] tetrahedron and the [MgO6] octahedron (solid curves) decrease smoothly, and the ECoN parameters remain nearly constant (dashed curves). Overall, the microstructural features of the Fo phase do not change suddenly.

Compression Behavior of Fo, Fo-II, Fo-III and Fo-IV at 0 K
Our optimized volume data for the Fo, Fo-II, Fo-III and Fo-IV phases at different P and 0 K are shown in Figure 6. Using a least-squares method, the pressure-volume data have been fitted to the third-order Birch-Murnaghan equation of state (3rd BM-EoS; Birch [70]), where P is pressure (GPa), K T isothermal bulk modulus (GPa), K T the first pressure derivative of the K T , V 0 the volume at zero pressure, and V the volume at high pressure. The thus-obtained EoS parameters are listed in Table 4.
For the Fo phase, we have obtained V 0 = 295.96 Å 3 (fixed), K T = 127.3(5) GPa and K T = 3.96(3) (Figure 6a). The EoS of the Fo phase was determined by a range of techniques such as high-P single-crystal XRD [37,54,71,72], high-P powder XRD [73], Brillouin scattering spectroscopy [74], and theoretical simulation [75,76]; see Finkelstein et al. [37] for a summary of the literature. In general, the results reported by the earlier studies are in very good agreement with ours. In specific, the K T from Finkelstein et al. [37] is slightly larger than the present determination ( Figure 6a); their K T at 300 K (K T = 131.5(4) GPa, as obtained in this study by refitting their volume-pressure data with the 3rd BM-EoS) is~3.9% larger than ours at 0 K (K T = 126.6(1) GPa), on the basis of K T = 4. Two factors may have made their contributions to this phenomenon: (1) our theoretical calculation, done with the GGA method, is well known to lead to an underestimated K T [68,77]; (2) the experimental data, collected with the P medium of He which is unquestionably the best available P-transmitting medium but still develops some deviatoric stress at P >~40 GPa [78], are expected to generate slightly overestimated K T . If we consider the effect of T on the K T (0 K in the theoretical investigation versus 298 K in the compression experiments), the influences of these two factors become even more prominent.

Compression Behavior of Fo, Fo-II, Fo-III and Fo-IV at 0 K
Our optimized volume data for the Fo, Fo-II, Fo-III and Fo-IV phases at different P and 0 K are shown in Figure 6. Using a least-squares method, the pressure-volume data have been fitted to the third-order Birch-Murnaghan equation of state (3rd BM-EoS; Birch [70]), where P is pressure (GPa), KT isothermal bulk modulus (GPa), the first pressure derivative of the , V0 the volume at zero pressure, and V the volume at high pressure. The thus-obtained EoS parameters are listed in Table 4.   (3) 3.75(4) 4.11 (5) 3.94(4)    (2) 3.33 (2) 4.17 (6) 7.00 (11) c-axis c 0 (Å, fixed) 6.0188 5.7633 9.4266 5.0421 (9) 4.38 (2) 3.98 (2) For the Fo-II phase, we have obtained V 0 = 274.58 Å 3 (fixed), K T = 155.9(7) GPa and K T = 3.75(4) (Figure 6b), with the K T value potentially slightly underestimated because of the GGA method. Using the theoretically constrained 3rd BM-EoS for both phases, the volume of the Fo-II phase is thus~7.2% smaller than that of the Fo at 0 GPa, or~4.5% smaller at 52.4 GPa. The latter is in good agreement with the experimental observation of an~4.8% volume reduction associated with the Fo/Fo-II phase transition at~52.4 GPa and room T [37].
Our volume-pressure data of the Fo-III phase are compared in Figure 6c to the experimental data collected by Finkelstein et al. [37]. These two data sets are broadly compatible to each other. In detail, the experimental data, surprisingly, define a relatively large compressibility (or a relatively small K T ) for the Fo-III phase. Remarkably, Finkelstein et al. [37] reported a high density of defects in their Fo-III crystal, which might have resulted in the relatively large compressibility. Fitted with the 3rd BM-EoS, our volume-pressure data constrain the V 0 , K T and K T as 249.17 Å 3 (fixed), 184.9(8) GPa and 4.11(5), respectively ( Figure 6c); as usually, the K T value might have been slightly underestimated with the GGA method. Using the theoretically determined 3rd BM-EoS for both phases, the volume of the Fo-III phase is thus~15.8% smaller than that of the Fo at 0 GPa, or~10.4% smaller at 58.2 GPa. The latter is in general agreement with the~9% volume reduction extrapolated for the assumed Fo/Fo-III phase transition at~58.2 GPa and room T, taking into account the influence of the high defect density in the Fo-III crystal [37].
The Fo-III phase is a non-centrosymmetric version of the Cmcm CaTi 2 O 4 -type structure, or a post-spinel phase [37]. Compared on the basis of K T = 4, our K T value for the Fo-III phase (186.7 GPa) is only slightly larger than that of the Mg 2 SiO 4 -spinel, or the Rwd (184(1) GPa; see Liu et al. [79] for a summary of the literature), by~1.5%, which may have been severely affected by the GGA method used in our simulation. If we compare our GGA result for the Fo-III phase to the GGA result for the Rwd (169.3 GPa; Yu and Wentzcovitch [80]), a much larger difference in the K T values,~10.3%, is found. Similarly comparable K T values of the spinel and the CaTi 2 O 4 -type post-spinel phases with the Zn 2 TiO 4 composition, 162.9(5) versus 171(3) GPa (a 5.0% relative difference; K T = 4; the GGA method), were recently disclosed by Zhang et al. [81]. In comparison, Zhang et al. [82] found a much large difference of~31.3% in the K T values for the spinel (167.2(12) GPa with K T = 4; diamond-anvil cell compression experiments at room T) and the CaTi 2 O 4 -type post-spinel phases (219.6(26) GPa with K T = 4; the GGA method at zero T) with the Co 2 TiO 4 composition. There are however perplexing factors in the latter case, such as Co being a transition metal element and an intermediate phase For the Fo-IV phase, we have obtained V 0 = 260.00 Å 3 (fixed), K T = 163.4 GPa and K T = 3.94 ( Figure 6d). As expected, our V 0 value obtained with the GGA method is slightly overestimated whereas our K T value is slightly underestimated, compared to the values obtained with the LDA method by Zhou et al. [40]. Using our theoretically determined 3rd BM-EoS for the Fo-IV and Fo phases (Tables 3 and 4), the volume difference at 0 GPa is thus~12.2%, smaller than the volume difference between the Fo-III and Fo phases (~15.8%).
A linearized 3rd BM-EoS [83] was used to obtain the parameters of the equations of states for the crystallographic axes of the Fo, Fo-II, Fo-III and Fo-IV. The results are also listed in Table 4. Summarily, the most to least compressible axes, in order, are b > c > a for the Fo, a > c > b for the Fo-II, b > a > c for the Fo-III, and c > a > b for the Fo-IV. Our theoretically constrained anisotropic compressibility for the Fo is therefore well in line with the experimental results such as Downs et al. [71], Finkelstein et al. [37] and Zha et al. [74].

Phase Relations at High P and 0 K
With their nearly hydrostatic compression experiments and metadynamics simulations, Finkelstein et al. [37] claimed that the Fo-II and Fo-III phases are energetically metastable compared to their high-T counterparts (Fo, Wds and Rwd), and the Fo phase transforms to the Fo-II phase at 50 GPa and later to the Fo-III phase at~58 GPa. By performing ab initio molecular dynamics and DFT calculations (LDA), Zhou et al. [40] predicted that the Fo-IV phase is relatively metastable compared to the Fo phase at P <~22 GPa, compared to Wds at P <~35 GPa, and compared to Rwd at P <~43. 46 GPa, but has the lowest enthalpy among these phases at P >~43.46 GPa.
In this study, we have used the optimized structures and corresponding static energies (0 K and 0-100 GPa) for a large number of relevant individual phases, as listed in Table 1, to obtain the zero K enthalpies (H 0 ) of the Fo, Fo-II, Fo-III, Fo-IV, Wds, Rwd, Brg + Pe and Aki + Pe phase assemblages under the constraint of a bulk composition Mg 2 SiO 4 , and explored the phase relations among these phase assemblages. The results are summarized in Table 5 and shown in Figure 7. Table 5. P range of the thermodynamically stable phase assemblage for the Mg 2 SiO 4 composition and P for the polymorphic phase transition from Fo, Wds and Rwd to Fo-II, Fo-III and Fo-IV, based on our calculated enthalpies at 0 K (0-100 GPa; Figure 7). In terms of the enthalpies at 0 K, the stable phase assemblage for the bulk composition Mg 2 SiO 4 changes from the Fo to the Wds, Rwd, Aki + Pe, and eventually to the Brg + Pe, as P increases from 0 to 100 GPa (Table 5). This result implies that the Fo-II, Fo-III and Fo-IV phases are indeed thermodynamically metastable, in full agreement with Finkelstein et al. [37] and Zhou et al. [40]. Further, it is clear now that the Fo-IV phase does have a lower enthalpy than the Rwd at P > 55.5 GPa (the GGA result; Figure 7) or at P >~43.46 GPa (the LDA result; Zhou et al. [40]), but has higher enthalpy than the Fo-III, Aki + Pe and Brg + Pe phase assemblages (Figure 7). Zhou et al. [40] did not consider the phase assemblages Fo-III, Aki + Pe and Brg + Pe in their study, so that their simulation was incomplete.

Stable Phase Assemblage
As to the phase relations of the Fo and its polymorphs, Figure 7 shows that for the P interval of 22.1-24.9 GPa, the Fo-III phase has the minimum difference in enthalpies to the Fo phase, suggesting a possible formation of the Fo-III phase directly from the Fo phase. At P >~32.2 GPa, Figure 7 indicates that the Fo-II phase has the minimum difference in enthalpies, and might readily form from the Fo staring material due to the small energy barrier. At all P >~22.1 GPa, further, the Fo-III phase has the lowest enthalpy among the phases Fo, Fo-II, Fo-III and Fo-IV (Figure 7). Provided that the Fo-II phase does not break down to the phase assemblages such as the Brg + Pe and Aki + Pe, it eventually transforms to the Fo-III phase. These theoretically calculated phase transition sequences, Fo → Fo-III in the P interval of~22-24.9 GPa and Fo → Fo-II → Fo-III at P >~32.2 GPa, are broadly compatible with the experimental observations of Finkelstein et al. [37] and Newman et al. [39]. In the shockwave experiments of Newman et al. [39], the transition pressure for the Fo → Fo-III phase transition was constrained as~44 GPa, much higher than our predicted values (~22.1-24.9 GPa). This might largely reflect the fact of an extremely short pressurizing duration in the shockwave experiments. Similarly, the Fo → Fo-II → Fo-III phase transitions were observed at much higher P in the hydrostatic compression experiments with the diamond-anvil cell (at P >~50 GPa; Finkelstein et al. [37]). In this case, the low experimental T (300 K) and short experimental duration (mostly a few days) might have played important roles. Additionally, the Fo-IV phase has the minimum difference in enthalpies to the Fo phase in the P range from~24.9 to 32.2 GPa (Table 5), presumably suggesting that it might metastably form from the starting material Fo as well, but somehow either be kinetically hindered or missed in the compression experiments.
~50 GPa and later to the Fo-III phase at ~58 GPa. By performing ab initio molecular dynamics and DFT calculations (LDA), Zhou et al. [40] predicted that the Fo-IV phase is relatively metastable compared to the Fo phase at P < ~22 GPa, compared to Wds at P < ~35 GPa, and compared to Rwd at P < ~43. 46 GPa, but has the lowest enthalpy among these phases at P > ~43. 46 GPa.
In this study, we have used the optimized structures and corresponding static energies (0 K and 0-100 GPa) for a large number of relevant individual phases, as listed in Table 1, to obtain the zero K enthalpies (H0) of the Fo, Fo-II, Fo-III, Fo-IV, Wds, Rwd, Brg + Pe and Aki + Pe phase assemblages under the constraint of a bulk composition Mg2SiO4, and explored the phase relations among these phase assemblages. The results are summarized in Table 5 and shown in Figure 7. Table 5. P range of the thermodynamically stable phase assemblage for the Mg2SiO4 composition and P for the polymorphic phase transition from Fo, Wds and Rwd to Fo-II, Fo-III and Fo-IV, based on our calculated enthalpies at 0 K (0-100 GPa; Figure 7).  It should be interesting to examine what metastable phases might form from a different starting material of either Wds or Rwd. According to Figure 7 and Table 5, the Fo-III phase and Fo-IV phase might metastably form from the starting material Wds at 30.6 GPa and 41.1 GPa, respectively; the Fo-IV phase might further transform into the Fo-III phase with a lower enthalpy though. Kleppe et al. [84] hydrostatically compressed hydrous Wds up to~50 GPa, but did not observe any new framework Raman mode. In comparison, Figure 7 suggests that the Fo-III phase and Fo-IV phase might metastably form from the starting material Rwd at 36.2 GPa and 55.5 GPa, respectively; the Fo-IV might also further transform into the Fo-III phase. Interestingly, Kleppe et al. [85] hydrostatically compressed the Rwd phase to~56.5 GPa and observed a new framework Raman mode at~30 GPa (~802 cm −1 , as extrapolated to 1 bar). Presently it is unclear whether or not this new Raman mode is associated with any of the potential metastable polymorphs, Fo-III and Fo-IV.

Stable Phase Assemblage
Summarily, the Fo-III phase is by far the most stable metastable Mg 2 SiO 4 polymorph at P > 22.1 GPa on the ground of the zero K enthalpies of the Fo, Fo-II, Fo-III, Fo-IV, Wds, Rwd, Brg + Pe and Aki + Pe phase assemblages (Figure 7), and may play a pivotal role in the high-P phase relations of the composition Mg 2 SiO 4 at very low T.

Dynamic Stability of Fo, Fo-II, Fo-III and Fo-IV
For a dynamically (mechanically) stable crystal, all phonon frequencies (or eigenvalues) should be positive under the harmonic oscillator approximation [63,64]. In other words, any imaginary phonon frequency or negative eigenvalue means a dynamic instability. To examine the dynamic stabilities of the Fo, Fo-II, Fo-III and Fo-IV phases at different P within the harmonic oscillation approximation framework, we have calculated the phonon frequencies at the Γ point in the first Brillouin zone of these structures, and obtained the entire phonon dispersion curves using the PHONOPY code. In addition, the non-analytic term corrections, as outlined by Pick et al. [86], have been described and calculated using the method of Giannozzi et al. [87] and Gonze and Lee [61]. All the calculated dispersion curves are plotted and shown in Figures 8-11. Figure 8 shows that for the P range of 0-80 GPa the Fo structure is dynamically stable, in good agreement with previous theoretical result obtained at P up to 20 GPa [44,75]. At higher P such as 90 and 100 GPa (Figure 8j,k), the phonon dispersion curves display some imaginary phonon frequencies, suggesting a dynamic instability for the Fo structure at P >~85 GPa.

Dynamic Stability of Fo, Fo-II, Fo-III and Fo-IV
For a dynamically (mechanically) stable crystal, all phonon frequencies (or eigenvalues) should be positive under the harmonic oscillator approximation [63,64]. In other words, any imaginary phonon frequency or negative eigenvalue means a dynamic instability. To examine the dynamic stabilities of the Fo, Fo-II, Fo-III and Fo-IV phases at different P within the harmonic oscillation approximation framework, we have calculated the phonon frequencies at the Γ point in the first Brillouin zone of these structures, and obtained the entire phonon dispersion curves using the PHONOPY code. In addition, the non-analytic term corrections, as outlined by Pick et al. [86], have been described and calculated using the method of Giannozzi et al. [87] and Gonze and Lee [61]. All the calculated dispersion curves are plotted and shown in Figures 8-11. Figure 8 shows that for the P range of 0-80 GPa the Fo structure is dynamically stable, in good agreement with previous theoretical result obtained at P up to 20 GPa [44,75]. At higher P such as 90 and 100 GPa (Figure 8j,k), the phonon dispersion curves display some imaginary phonon frequencies, suggesting a dynamic instability for the Fo structure at P > ~85 GPa.  Figure 9 shows the phonon dispersion curves for the Fo-II phase from 0 to 100 GPa. At all these P, a large number of imaginary frequencies have been observed, suggesting a structural instability. This result is generally consistent with the experimental phenomenon that very preliminary XRD data could be obtained for the Fo-II phase at one single P only (~52.4 GPa; Finkelstein et al. [37]).  Figure 9 shows the phonon dispersion curves for the Fo-II phase from 0 to 100 GPa. At all these P, a large number of imaginary frequencies have been observed, suggesting a structural instability. This result is generally consistent with the experimental phenomenon that very preliminary XRD data could be obtained for the Fo-II phase at one single P only (~52.4 GPa; Finkelstein et al. [37]). As to the Fo-III phase, we have found no imaginary frequency for the investigated P range from 0 to 100 GPa (see Figure 10). In addition, we do not detect any significant vibrational mode As to the Fo-III phase, we have found no imaginary frequency for the investigated P range from 0 to 100 GPa (see Figure 10). In addition, we do not detect any significant vibrational mode softening at every high-symmetry point along the k-path, which means this phase will remain as dynamically stable under much higher P. softening at every high-symmetry point along the k-path, which means this phase will remain as dynamically stable under much higher P. Previous theoretical study showed that the Fo-IV phase was dynamically stable at 44 and 100 GPa (Figure 4 in Zhou et al. [40]). When our phonon dispersion curves ( Figure 11) are plotted using a more detailed and complete high-symmetry k-path in the Brillouin zone (Γ-X-S-Y-Γ-Z-U-R-T-Z|Y-T|U-X|S-R; Setyawan and Curtarolo [88]), however, a large number of imaginary frequencies emerge at all P. It follows that the Fo-IV phase is also dynamically instable in the studied P range. Previous theoretical study showed that the Fo-IV phase was dynamically stable at 44 and 100 GPa (Figure 4 in Zhou et al. [40]). When our phonon dispersion curves ( Figure 11) are plotted using a more detailed and complete high-symmetry k-path in the Brillouin zone (Γ-X-S-Y-Γ-Z-U-R-T-Z|Y-T|U-X|S-R; Setyawan and Curtarolo [88]), however, a large number of imaginary frequencies emerge at all P. It follows that the Fo-IV phase is also dynamically instable in the studied P range. Summarily, our extensive phonon dispersion data suggest that the Fo-III phase is the best candidate as a metastable phase forming from the Fo starting material, and the Fo-II phase is likely an intermediate transitional structure with no dynamic stability. Interestingly, recent shock-wave compression experiments with polycrystalline Fo staring materials performed at 44(3) and 73 (5) GPa have showed a direct transformation from the Fo phase to the Fo-III phase, without involving the Fo-II phase (Newman et al. [39]). In addition, the identity of the Fo-IV phase has not been established, either theoretically or experimentally. Nevertheless, the role of the Fo-IV phase is unlikely geologically important due to its much higher enthalpy than the Fo-III phase (Figure 7). Summarily, our extensive phonon dispersion data suggest that the Fo-III phase is the best candidate as a metastable phase forming from the Fo starting material, and the Fo-II phase is likely an intermediate transitional structure with no dynamic stability. Interestingly, recent shock-wave compression experiments with polycrystalline Fo staring materials performed at 44(3) and 73(5) GPa have showed a direct transformation from the Fo phase to the Fo-III phase, without involving the Fo-II phase (Newman et al. [39]). In addition, the identity of the Fo-IV phase has not been established, either theoretically or experimentally. Nevertheless, the role of the Fo-IV phase is unlikely geologically important due to its much higher enthalpy than the Fo-III phase (Figure 7).

Thermodynamic Properties of Fo-III
Using calculated phonon dispersion curves and corresponding vibrational densities of state, the thermodynamic properties of one phase can be calculated by the QHA method [44,65,66]. With the PHONOPY code, we have calculated some thermodynamic properties for the Fo, Wds, Rwd, Brg, Aki, Pe and Fo-III phases, such as the Gibbs free energies G, equilibrium volumes V, isothermal bulk moduli K T and their first pressure derivatives K T , adiabatic bulk moduli K S , bulk sound velocities v Φ , thermal expansion coefficients α, thermal Grüneisen parameters γ th , entropies S, isochoric heat capacities C V and isobaric heat capacities C P . Due to the dynamic instabilities of the Fo-II and Fo-IV phases, these properties cannot be derived.
Our calculated V, K T , K T , K S , v Φ , α, γ th , S, C V and C P at 1 atm and 300 K for the Fo, Wds, Rwd, Aki, Brg, and Pe phases are compared with the experimental results in Table 6. As expected, our results obtained by employing the GGA method show systematic deviations to the experimental results, slightly overestimating the V (maximally by 4.4%) and consequently underestimating the K T (maximally by 13.8%), K S (maximally by 14.5%) and v Φ (maximally by 5.5%). In addition, our calculated S, C V and C P closely match the experimentally determined values, with an utmost difference of~13.3%, 5.6% and 6.1%, respectively. Further, the largest deviation between our calculated results and the experimental determinations comes with the α (minimally by 9.4%), which is presumably largely due to the low accuracies in the experimental results. Overall, we believe that the thermodynamic properties of the Fo, Wds, Rwd, Aki, Brg, and Pe phases have been well simulated in this study. This claim is strongly supported by the excellent agreement between our values and those obtained by Hernandez et al. [44] using similar theoretic techniques (the GGA method) but with slightly different simulating parameters (e.g., the planewave cutoff energy of 500 eV): for the Fo, Wds, Rwd, Aki, Brg and Pe phases, the absolute relative difference is~1.9-19.9% for the α,~0.7-2.1% for the S,~0.1-2.2% for the C V , and~0.1-1.4% for the C P . The Fo-III phase was experimentally discovered lately [37,39]. More importantly, it formed from the Fo starting materials at very high P (P >~58 GPa in the hydrostatic compression experiments or P >~44(3) GPa in the dynamic shockwave experiments), and became amorphous as the experimental P decreased to ambient P. Consequently, nearly all of its thermodynamic properties are presently undetermined.
Our simulated thermodynamic properties of the Fo-III phase, as obtained by employing the GGA method, are listed in Table 6. To provide better constraints, these properties have been additionally calculated using the LDA method; this practice is deemed necessary because most thermodynamic properties of this phase cannot be experimentally measured, and because the LDA method has better performance in predicting mineral physical properties than the GGA method. The results obtained by these two methods for 1 atm and 300 K are compared in Table 6. The relative differences between these two sets of values are~6.1%, −18.2%, −17.8%, −5.2%, 17.4%, 4.2%, 11.5%, 4.2% and 4.6% for the V, K T , K S , v Φ , α, γ th , S, C V and C P , respectively. The true values of these thermodynamic properties are expected to be somewhere in between these two sets of values obtained by the GGA and LDA methods, presumably closer to the LDA results.
The P-V curves of the Fo-III phase calculated by both the GGA and the LDA methods at 0-100 GPa and 300-1500 K are plotted in Figure 12, along with the relevant parameters. For comparison, the experimental results from Finkelstein et al. [37] are also shown. Usually, the LDA results for the Fo-III phase are expected to be in better agreement with the experimentally observed volumes than the GGA results [44,66]. This is apparently the case for the Fo-III ( Figure 12). Comparing our LDA results with the experimental data, it is evident that the experimentally determined volumes have been surprisingly overestimated, by~1.5% on average, which might be attributable to the high defect density in the Fo-III crystals, as observed by Finkelstein et al. [37]. The Fo-III phase was experimentally discovered lately [37,39]. More importantly, it formed from the Fo starting materials at very high P (P > ~58 GPa in the hydrostatic compression experiments or P > ~44(3) GPa in the dynamic shockwave experiments), and became amorphous as the experimental P decreased to ambient P. Consequently, nearly all of its thermodynamic properties are presently undetermined.
Our simulated thermodynamic properties of the Fo-III phase, as obtained by employing the GGA method, are listed in Table 6. To provide better constraints, these properties have been additionally calculated using the LDA method; this practice is deemed necessary because most thermodynamic properties of this phase cannot be experimentally measured, and because the LDA method has better performance in predicting mineral physical properties than the GGA method. The results obtained by these two methods for 1 atm and 300 K are compared in Table 6. The relative differences between these two sets of values are ~6.1%, −18.2%, −17.8%, −5.2%, 17.4%, 4.2%, 11.5%, 4.2% and 4.6% for the V, KT, KS, vΦ, α, γth, S, CV and CP, respectively. The true values of these thermodynamic properties are expected to be somewhere in between these two sets of values obtained by the GGA and LDA methods, presumably closer to the LDA results.
The P-V curves of the Fo-III phase calculated by both the GGA and the LDA methods at 0-100 GPa and 300-1500 K are plotted in Figure 12, along with the relevant parameters. For comparison, the experimental results from Finkelstein et al. [37] are also shown. Usually, the LDA results for the Fo-III phase are expected to be in better agreement with the experimentally observed volumes than the GGA results [44,66]. This is apparently the case for the Fo-III ( Figure 12). Comparing our LDA results with the experimental data, it is evident that the experimentally determined volumes have been surprisingly overestimated, by ~1.5% on average, which might be attributable to the high defect density in the Fo-III crystals, as observed by Finkelstein et al. [37]. Our calculated α, γth, CV and CP for the Fo-III phase at 0-100 GPa and 300-1500 K, using both the GGA and the LDA methods, are fitted as polynomial functions of T at different P employing the equations used by Xiong et al. [91]. The only exception is for the γth calculated by the LDA method, in which case we have found that some small modification to the equation used is necessary in order to well reproduce our calculated values. These equations and the fitting constants are listed in Table 7 (the GGA results) and Table 8 (the LDA results). The quality of our data-fitting can be gauged by the average absolute deviation (AAD, %) of the theoretically calculated values and the fitted values. For the α, the AAD increases from ~0.10% at 0 GPa to 0.70% at 100 GPa (GGA) or from 0.13% at 0 GPa to 0.79% at 100 GPa (LDA). For the γth, the AAD increases from ~0.07% at 0 GPa to 0.22% at 100 GPa (GGA) or remains less than 0.03% for all P (LDA). For the CV, the AAD remains less than ~0.56% at all P (both GGA and LDA). For the CP, the AAD decreases from ~1.83% at 0 GPa to 0.18% at 100 GPa (GGA) or from 0.72% at 0 GPa to 0.21% at 100 GPa (LDA). The calculated α, γth, CV and CP for the Fo-III phase are shown in the Supplementary Figures S1-S4. Our calculated α, γ th , C V and C P for the Fo-III phase at 0-100 GPa and 300-1500 K, using both the GGA and the LDA methods, are fitted as polynomial functions of T at different P employing the equations used by Xiong et al. [91]. The only exception is for the γ th calculated by the LDA method, in which case we have found that some small modification to the equation used is necessary in order to well reproduce our calculated values. These equations and the fitting constants are listed in Table 7 (the GGA results) and Table 8 (the LDA results). The quality of our data-fitting can be gauged by the average absolute deviation (AAD, %) of the theoretically calculated values and the fitted values. For the α, the AAD increases from~0.10% at 0 GPa to 0.70% at 100 GPa (GGA) or from 0.13% at 0 GPa to 0.79% at 100 GPa (LDA). For the γ th , the AAD increases from~0.07% at 0 GPa to 0.22% at 100 GPa (GGA) or remains less than 0.03% for all P (LDA). For the C V , the AAD remains less than~0.56% at all P (both GGA and LDA). For the C P , the AAD decreases from~1.83% at 0 GPa to 0.18% at 100 GPa (GGA) or from 0.72% at 0 GPa to 0.21% at 100 GPa (LDA). The calculated α, γ th , C V and C P for the Fo-III phase are shown in the Supplementary Figures S1-S4. Table 7. Polynomial fitting results for α, γ th , C V and C P from 300 to 1500 K at different P (GPa), calculated using the GGA method.     Table 8. Polynomial fitting results for α, γ th , C V and C P from 300 to 1500 K at different P (GPa), calculated using the LDA method.

Phase Boundaries under High P and High T
The QHA method is a very effective approach to constrain the thermodynamic properties of solid phases at relatively low T [66]. When T is higher than some characteristic value (T C ), the phonon-phonon interactions become important, so that the QHA method cannot work properly. It is hence important to determine the T C of every phase in the study system if one aims at an accurate calculation of the phase relations. As suggested by Carrier et al. [65] and Wentzcovitch et al. [92], the T C of certain solid phase may be effectively evaluated by detecting the inflection point of the α-T profile. We have followed this suggestion in this study.
The T C for the Fo, Wds, Rwd, Aki, Brg, Pe and Fo-III phases has been evaluated from 0 to 30 GPa, with the result plotted in Figure 13. For every phase, the T C at 0 GPa is higher than 1000 K. The averaged slope dT C /dP is~26 K/GPa for Fo,~24 K/GPa for Wds,~23 K/GPa for Rwd,~23 K/GPa for Aki,~16 K/GPa for Brg,~18 K/GPa for Pe, and~24 K/GPa for Fo-III. Further, Figure 13 suggests that with the exceptions of the Fo-III phase and the Aki, all other phases in normal mantle may show various degrees of anharmonicity. As a result, the relevant phase relations calculated with the thermodynamic properties obtained by the QHA method may bear relatively large uncertainty.
In contrast, Figure 13 shows that similar phase relation calculation with identical techniques should be generally appropriate for the subducting slabs due to their much low T.

Phase Boundaries under High P and High T
The QHA method is a very effective approach to constrain the thermodynamic properties of solid phases at relatively low T [66]. When T is higher than some characteristic value (TC), the phonon-phonon interactions become important, so that the QHA method cannot work properly. It is hence important to determine the TC of every phase in the study system if one aims at an accurate calculation of the phase relations. As suggested by Carrier et al. [65] and Wentzcovitch et al. [92], the TC of certain solid phase may be effectively evaluated by detecting the inflection point of the α-T profile. We have followed this suggestion in this study.
The TC for the Fo, Wds, Rwd, Aki, Brg, Pe and Fo-III phases has been evaluated from 0 to 30 GPa, with the result plotted in Figure 13. For every phase, the TC at 0 GPa is higher than 1000 K. The averaged slope dTC/dP is ~26 K/GPa for Fo, ~24 K/GPa for Wds, ~23 K/GPa for Rwd, ~23 K/GPa for Aki, ~16 K/GPa for Brg, ~18 K/GPa for Pe, and ~24 K/GPa for Fo-III. Further, Figure 13 suggests that with the exceptions of the Fo-III phase and the Aki, all other phases in normal mantle may show various degrees of anharmonicity. As a result, the relevant phase relations calculated with the thermodynamic properties obtained by the QHA method may bear relatively large uncertainty. In contrast, Figure 13 shows that similar phase relation calculation with identical techniques should be generally appropriate for the subducting slabs due to their much low T. Figure 13. Calculated characteristic temperatures (TC) for Fo, Wds, Rwd, Aki, Brg, Pe and Fo-III from 0 to 30 GPa. The P-T profile of the cold slab or hot slab is sketched according to Thompson [93], while that of the normal mantle is from Anderson [94].
The phase relations of the bulk composition Mg2SiO4 at P conditions from the upper mantle to the uppermost part of the lower mantle have been calculated using the thermodynamic properties theoretically obtained for the Fo, Fo-III, Wds, Rwd, Aki, Brg, and Pe phase. The result is shown in Figure 14.
For the normal mantle temperatures [94], the Fo/Wds phase boundary has been located at ~15 GPa with an average Clapeyron slope (averaged value from 300 to 1800 K) of 2.2 MPa/K, in good agreement with previous experimental and theoretical results [66,[95][96][97]. On the contrary, the Wds/Rwd phase boundary locates at ~22.2 GPa and 1800 K with an average Clapeyron slope of 4.2 MPa/K, and the Rwd/Brg + Pe phase boundary locates at ~21.4 GPa and 1800 K with an average Clapeyron slope of −1.5 MPa/K, which obviously contradicts the well-accepted phase transition sequence, Wds → Rwd → Brg + Pe, of the normal mantle in the depth range of ~410-660 km. Indeed, the thus-obtained high-T phase relations possess some uncertainty.
At 300 K, the Fo phase will transform into the Fo-III phase at ~23.75 GPa if other stable phase transition paths are kinetically inhibited (i.e., Fo/Wds, Wds/Rwd and Rwd/Brg + Pe). Since the GGA method often overestimates a phase boundary by ~3 GPa [66], the actual phase transition P between the Fo and Fo-III at 300 K is probably ~21 GPa, which is significantly lower than the Rwd/Brg + Pe phase transition P of the ambient mantle (~24 GPa). Further, the averaged Clapeyron slope for the Fo/Fo-III phase transition has been estimated as −1.1 MPa/K for the T range of 300-1800 K, so that Figure 13. Calculated characteristic temperatures (T C ) for Fo, Wds, Rwd, Aki, Brg, Pe and Fo-III from 0 to 30 GPa. The P-T profile of the cold slab or hot slab is sketched according to Thompson [93], while that of the normal mantle is from Anderson [94].
The phase relations of the bulk composition Mg 2 SiO 4 at P conditions from the upper mantle to the uppermost part of the lower mantle have been calculated using the thermodynamic properties theoretically obtained for the Fo, Fo-III, Wds, Rwd, Aki, Brg, and Pe phase. The result is shown in Figure 14.
For the normal mantle temperatures [94], the Fo/Wds phase boundary has been located at 15 GPa with an average Clapeyron slope (averaged value from 300 to 1800 K) of 2.2 MPa/K, in good agreement with previous experimental and theoretical results [66,[95][96][97]. On the contrary, the Wds/Rwd phase boundary locates at~22.2 GPa and 1800 K with an average Clapeyron slope of 4.2 MPa/K, and the Rwd/Brg + Pe phase boundary locates at~21.4 GPa and 1800 K with an average Clapeyron slope of −1.5 MPa/K, which obviously contradicts the well-accepted phase transition sequence, Wds → Rwd → Brg + Pe, of the normal mantle in the depth range of~410-660 km. Indeed, the thus-obtained high-T phase relations possess some uncertainty.
At 300 K, the Fo phase will transform into the Fo-III phase at~23.75 GPa if other stable phase transition paths are kinetically inhibited (i.e., Fo/Wds, Wds/Rwd and Rwd/Brg + Pe). Since the GGA method often overestimates a phase boundary by~3 GPa [66], the actual phase transition P between the Fo and Fo-III at 300 K is probably~21 GPa, which is significantly lower than the Rwd/Brg + Pe phase transition P of the ambient mantle (~24 GPa). Further, the averaged Clapeyron slope for the Fo/Fo-III phase transition has been estimated as −1.1 MPa/K for the T range of 300-1800 K, so that the Fo/Fo-III phase transition P decreases as T increases: taking 900 K as an example, this phase transition P should be reduced by~0. 66 GPa. There are extra factors such as the iron content, stress state and water abundance, which may affect the phase transition P [38]. Additionally, it is worth to point out that this phase transition is endothermic (∆V = −4.86 cm 3 /mol and ∆H =~1.60 kJ/mol at 300 K), which may be beneficial for the slabs to maintain their low T. the Fo/Fo-III phase transition P decreases as T increases: taking 900 K as an example, this phase transition P should be reduced by ~0. 66 GPa. There are extra factors such as the iron content, stress state and water abundance, which may affect the phase transition P [38]. Additionally, it is worth to point out that this phase transition is endothermic (ΔV = −4.86 cm 3 /mol and ΔH = ~1.60 kJ/mol at 300 K), which may be beneficial for the slabs to maintain their low T. Figure 14. Calculated P-T boundaries for the Fo/Wds, Wds/Rwd, Rwd/Aki + Pe, Rwd/Brg + Pe, Aki/Brg and Fo/Fo-III reactions at 10-28 GPa and 300-1900 K. The P-T profile of the cold slab or hot slab is sketched according to Thompson [93], while that of the normal mantle is from Anderson [94]. All data used in the calculation are from our theoretical simulation with the GGA method. According to Wentzcovitch et al. [66], all these phase boundaries at T < ~1500 K are potentially overestimated by ~3 GPa because of the used GGA method, as indicated for the Fo/Fo-III phase boundary by the arrow.

Fo-III in the Subduction Zone
Our theoretical investigation shows that the Fo-III phase has a dynamically stable structure with very low enthalpy and Gibbs free energy, and therefore might form as a metastable phase from the starting materials of Fo, Wds and Rwd under some special circumstances (Figures 7 and  10). Indeed, it has been observed in two kinds of high-P experiments with the Fo staring material (Figure 15), long duration room-T hydrostatic diamond-anvil cell experiments at one end [37] and short duration high-T dynamic shockwave experiments at the other end [39]. Due to the long geological history and high T of the mantle material, the Fo-III phase should never exist in the normal mantle. In contrast, the Fo-III phase might metastably exist in the cold slabs subducted to the upper mantle-lower mantle boundary region and even beyond.
The existence of a metastable Fo-III wedge in the slabs is subject to several requirements. Firstly, the T of the slabs subducted beyond the 410 km depth should be low enough to inhibit the transformation from the Ol to the thermodynamically stable phase Wds and Rwd. According to some studies on the kinetics of the Fo/Wds and Fo/Rwd phase transitions [1,6,7,26], the highest T for negligible Fo/Wds and Fo/Rwd phase transitions lies in the range of ~800-1100 K, closely matching the T conditions in the cold slabs ( Figure 15). This was the exact basis on which the metastable Ol hypothesis was proposed in the past. Secondly, the Fo-III phase should be able to form directly from the Ol. Our theoretical results have suggested that this could take place at ~20 GPa for a slab T of ~900 K while the static compression experiments have demonstrated the phase transition sequence Fo → Fo-II → Fo-III in the P range of ~50-90 GPa. As discussed earlier, the Fo-II phase is not a dynamically stable phase, so that the phase transition sequence Fo → Fo-III seems more appropriate. As to the phase transition pressure, the static compression experiments at room T might have severely overestimated it due to possible slow reaction rate. Usually, a large P overstepping is required in order to activate a polymorphic phase transformation at a low experimental T like 300 K. Further, the short experimental duration in the diamond-anvil cell compression experiments (~10 4 to 10 5 s), much longer than in the dynamic compression process Figure 14. Calculated P-T boundaries for the Fo/Wds, Wds/Rwd, Rwd/Aki + Pe, Rwd/Brg + Pe, Aki/Brg and Fo/Fo-III reactions at 10-28 GPa and 300-1900 K. The P-T profile of the cold slab or hot slab is sketched according to Thompson [93], while that of the normal mantle is from Anderson [94]. All data used in the calculation are from our theoretical simulation with the GGA method. According to Wentzcovitch et al. [66], all these phase boundaries at T <~1500 K are potentially overestimated bỹ 3 GPa because of the used GGA method, as indicated for the Fo/Fo-III phase boundary by the arrow.

Fo-III in the Subduction Zone
Our theoretical investigation shows that the Fo-III phase has a dynamically stable structure with very low enthalpy and Gibbs free energy, and therefore might form as a metastable phase from the starting materials of Fo, Wds and Rwd under some special circumstances (Figures 7 and 10). Indeed, it has been observed in two kinds of high-P experiments with the Fo staring material (Figure 15), long duration room-T hydrostatic diamond-anvil cell experiments at one end [37] and short duration high-T dynamic shockwave experiments at the other end [39]. Due to the long geological history and high T of the mantle material, the Fo-III phase should never exist in the normal mantle. In contrast, the Fo-III phase might metastably exist in the cold slabs subducted to the upper mantle-lower mantle boundary region and even beyond.
The existence of a metastable Fo-III wedge in the slabs is subject to several requirements. Firstly, the T of the slabs subducted beyond the 410 km depth should be low enough to inhibit the transformation from the Ol to the thermodynamically stable phase Wds and Rwd. According to some studies on the kinetics of the Fo/Wds and Fo/Rwd phase transitions [1,6,7,26], the highest T for negligible Fo/Wds and Fo/Rwd phase transitions lies in the range of~800-1100 K, closely matching the T conditions in the cold slabs ( Figure 15). This was the exact basis on which the metastable Ol hypothesis was proposed in the past. Secondly, the Fo-III phase should be able to form directly from the Ol. Our theoretical results have suggested that this could take place at~20 GPa for a slab T of 900 K while the static compression experiments have demonstrated the phase transition sequence Fo → Fo-II → Fo-III in the P range of~50-90 GPa. As discussed earlier, the Fo-II phase is not a dynamically stable phase, so that the phase transition sequence Fo → Fo-III seems more appropriate. As to the phase transition pressure, the static compression experiments at room T might have severely overestimated it due to possible slow reaction rate. Usually, a large P overstepping is required in order to activate a polymorphic phase transformation at a low experimental T like 300 K. Further, the short experimental duration in the diamond-anvil cell compression experiments (~10 4 to 10 5 s), much longer than in the dynamic compression process (~10 −3 s) but extremely short compared to the slab-subduction process (generally > 10 12 s), should have led to severe overestimation of the phase transition pressure. It thus appears highly possible that in the cold slabs a direct metastable polymorphic phase transition from the Ol to the Fo-III phase happens at~20 GPa, or in the lower part of the MTZ. Thirdly, the Fo-III phase should not be readily transform to other polymorphic phases or break down to other phase assemblages such as the Brg + Pe. Our theoretical results have suggested that among the metastable Fo-II, Fo-III and Fo-IV phases, the Fo-III phase is the best polymorphic candidate for the pressures out of the P stability field of the Fo. The diamond-anvil cell compression experiments have demonstrated a vast stable/metastable P range for the Fo-III phase indeed. On the other hand, it is presently unclear how fast the Fo-III phase may break down to the phase assemblages like the Brg + Pe (Figure 7). Considering all these factors, we propose that the metastable Fo-III phase, rather than the metastable Fo phase, may exist in the cold slabs subducted to the lower part of the MTZ (below the~600 km depth). (~10 −3 s) but extremely short compared to the slab-subduction process (generally > 10 12 s), should have led to severe overestimation of the phase transition pressure. It thus appears highly possible that in the cold slabs a direct metastable polymorphic phase transition from the Ol to the Fo-III phase happens at ~20 GPa, or in the lower part of the MTZ. Thirdly, the Fo-III phase should not be readily transform to other polymorphic phases or break down to other phase assemblages such as the Brg + Pe. Our theoretical results have suggested that among the metastable Fo-II, Fo-III and Fo-IV phases, the Fo-III phase is the best polymorphic candidate for the pressures out of the P stability field of the Fo. The diamond-anvil cell compression experiments have demonstrated a vast stable/metastable P range for the Fo-III phase indeed. On the other hand, it is presently unclear how fast the Fo-III phase may break down to the phase assemblages like the Brg + Pe (Figure 7). Considering all these factors, we propose that the metastable Fo-III phase, rather than the metastable Fo phase, may exist in the cold slabs subducted to the lower part of the MTZ (below the ~600 km depth). Figure 15. Formation P-T conditions of metastable Fo-III. The P-T profile of the cold slab or hot slab is sketched according to Thompson [93], while that of the normal mantle is from Anderson [94]. The lowest temperatures for any possible phase transition from the Fo-rich Ol to its equilibrium phase or phase assemblage (from ~14 to ~24 GPa), as estimated in Sung and Burns [1], Wu et al. [26], Kirby et al. [6] and Devaux et al. [7], are shown as short lines 1, 2 (2a for hydrostatic condition and 2b for nonhydrostatic condition), 3 and 4, respectively. The potential P overestimate for the Fo/Fo-III phase boundary caused by the GGA method is indicated by the tiny red arrow.
Consequently, we have calculated the bulk sound velocity (vΦ) of the Fo and Fo-III phases at 900, 1200 and 1500 K, using the following equation, The results are compared in Figure 16 to the preliminary reference Earth model (PREM) of the mantle [98]. Clearly, the vΦ of the Fo phase at all these temperatures is much slower than the vΦ of the ambient mantle ( Figure 16a; the GGA result): taking the 900 K as an example, the difference between the bulk sound velocities vary from ~4.1% (450 km depth) to ~8.5% (720 km depth), which may be easily detected by seismic observations. Interestingly, the vΦ of the Fo-III phase at 900 K for the depth range of ~600-660 km is remarkably faster than the vΦ of the ambient mantle (Figure 16b): the difference between the bulk sound velocities is ~0.1% (GGA) or ~3.3% (LDA) at the 600 km depth, which should be similarly detected by the seismic methods. Indeed, fast seismic velocities Figure 15. Formation P-T conditions of metastable Fo-III. The P-T profile of the cold slab or hot slab is sketched according to Thompson [93], while that of the normal mantle is from Anderson [94]. The lowest temperatures for any possible phase transition from the Fo-rich Ol to its equilibrium phase or phase assemblage (from~14 to~24 GPa), as estimated in Sung and Burns [1], Wu et al. [26], Kirby et al. [6] and Devaux et al. [7], are shown as short lines 1, 2 (2a for hydrostatic condition and 2b for nonhydrostatic condition), 3 and 4, respectively. The potential P overestimate for the Fo/Fo-III phase boundary caused by the GGA method is indicated by the tiny red arrow.
Consequently, we have calculated the bulk sound velocity (v Φ ) of the Fo and Fo-III phases at 900, 1200 and 1500 K, using the following equation, The results are compared in Figure 16 to the preliminary reference Earth model (PREM) of the mantle [98]. Clearly, the v Φ of the Fo phase at all these temperatures is much slower than the v Φ of the ambient mantle ( Figure 16a; the GGA result): taking the 900 K as an example, the difference between the bulk sound velocities vary from~4.1% (450 km depth) to~8.5% (720 km depth), which may be easily detected by seismic observations. Interestingly, the v Φ of the Fo-III phase at 900 K for the depth range of~600-660 km is remarkably faster than the v Φ of the ambient mantle (Figure 16b): the difference between the bulk sound velocities is~0.1% (GGA) or~3.3% (LDA) at the 600 km depth, which should be similarly detected by the seismic methods. Indeed, fast seismic velocities are usually observed for the slabs subducted to the depth range of~600-660 km, but are conventionally explained by resorting to low T, the termination of the metastable Ol wedge (i.e., phase transition from the metastable Ol to the Wds or Rwd), or both. Our result has clearly pointed out a potentially important role of the Fo-III phase in producing these fast seismic velocity signals. At the depths below the 660 km seismic discontinuity, the v Φ of the Fo-III phase at 900 K seems slower than the v Φ of the ambient mantle. The P-wave and S-wave velocities of the Fo-III phase are out of the scope of this study, and is currently under investigation. are usually observed for the slabs subducted to the depth range of ~600-660 km, but are conventionally explained by resorting to low T, the termination of the metastable Ol wedge (i.e., phase transition from the metastable Ol to the Wds or Rwd), or both. Our result has clearly pointed out a potentially important role of the Fo-III phase in producing these fast seismic velocity signals. At the depths below the 660 km seismic discontinuity, the vΦ of the Fo-III phase at 900 K seems slower than the vΦ of the ambient mantle. The P-wave and S-wave velocities of the Fo-III phase are out of the scope of this study, and is currently under investigation. In practical subducted slabs, other factors might blur the picture. In the case of the metastable Ol phase below the 410 km depth, its vΦ-depth profiles in Figure 16a should be shifted upwards by ~3% on the ground that the LDA method better predicts the mineral physical properties than the GGA method (like in the case of the Fo-III phase shown in Figure 16b), but they would be moved downwards by ~1.6% due to the incorporation of ~10% fayalite in the Ol phase [99]. In addition, water in the Ol might elevate the vΦ-depth profiles by ~0.8% [100]. Moreover, the unknown T in the real slabs might bring forth the presumably largest uncertainty in the vΦ-depth profiles. Considering all these variables, the vΦ-depth profiles for the metastable Ol phase may lie below the PREM model, and could be detected by the seismic methods with confidence [16,17]. In the case of our proposed metastable Fo-III wedge below the ~600 km depth, these issues should be carefully evaluated in future.
The isothermal density-depth profiles of Fo and Fo-III at 900, 1200 and 1500 K are compared in Figure 17 to the PREM density model of the mantle [98]. If Fo metastably exists in the cold oceanic slab at the depths > ~410 km, it should provide strong positive buoyancy force at least partially suppressing further subduction of the slab: its density at 900 K (the GGA result) is lower than the ambient mantle by ~10.5-18.3% in the depth range of ~450-720 km (Figure 17a). Both a lower T in the slab (an unlikely and extremely low T of 600 K for example) and some substitution of Fe for the Mg in the Fo phase (the typical upper mantle Ol of (Mg0.89Fe0.11)2SiO4 for example) elevate the density-depth profile somewhat, but never alter the density contrast between the metastable Fo (or Ol) and the warm ambient mantle [101]. This observation is therefore in agreement with all existing studies such as Schmeling et al. [102], Tetzlaff and Schmeling [8], and Zhang et al. [101].
The case of the Fo-III phase is somehow more complicated (Figure 17b). If it metastably forms from a direct polymorphic phase transition from the Ol at ~20 GPa (or the ~600 km depth), its density will not be smaller than the ambient mantle, and it will not be able to provide any positive buoyancy force to stop the sinking oceanic slabs, taking into account the density-increasing effect of the Fe Figure 16. Calculated bulk sound velocity (v Φ )-depth profiles at 900, 1200 and 1500 K for Fo (a); and Fo-III (b). As a reference, the PREM bulk sound velocity model [98] is also sketched. The shown bulk sound velocity (v Φ ) at 450 km depth, 600 km depth or 720 km depth and 900 K is calculated The GGA results at different T are shown as red lines while the LDA results as blue lines. For the purpose of illustration, the bulk sound velocity differences at certain depth obtained by the GGA method and the LDA method are slightly horizontally shifted.
In practical subducted slabs, other factors might blur the picture. In the case of the metastable Ol phase below the 410 km depth, its v Φ -depth profiles in Figure 16a should be shifted upwards bỹ 3% on the ground that the LDA method better predicts the mineral physical properties than the GGA method (like in the case of the Fo-III phase shown in Figure 16b), but they would be moved downwards by~1.6% due to the incorporation of~10% fayalite in the Ol phase [99]. In addition, water in the Ol might elevate the v Φ -depth profiles by~0.8% [100]. Moreover, the unknown T in the real slabs might bring forth the presumably largest uncertainty in the v Φ -depth profiles. Considering all these variables, the v Φ -depth profiles for the metastable Ol phase may lie below the PREM model, and could be detected by the seismic methods with confidence [16,17]. In the case of our proposed metastable Fo-III wedge below the~600 km depth, these issues should be carefully evaluated in future.
The isothermal density-depth profiles of Fo and Fo-III at 900, 1200 and 1500 K are compared in Figure 17 to the PREM density model of the mantle [98]. If Fo metastably exists in the cold oceanic slab at the depths >~410 km, it should provide strong positive buoyancy force at least partially suppressing further subduction of the slab: its density at 900 K (the GGA result) is lower than the ambient mantle by~10.5-18.3% in the depth range of~450-720 km (Figure 17a). Both a lower T in the slab (an unlikely and extremely low T of 600 K for example) and some substitution of Fe for the Mg in the Fo phase (the typical upper mantle Ol of (Mg 0.89 Fe 0.11 ) 2 SiO 4 for example) elevate the density-depth profile somewhat, but never alter the density contrast between the metastable Fo (or Ol) and the warm ambient mantle [101]. This observation is therefore in agreement with all existing studies such as Schmeling et al. [102], Tetzlaff and Schmeling [8], and Zhang et al. [101].
The case of the Fo-III phase is somehow more complicated (Figure 17b). If it metastably forms from a direct polymorphic phase transition from the Ol at~20 GPa (or the~600 km depth), its density will not be smaller than the ambient mantle, and it will not be able to provide any positive buoyancy force to stop the sinking oceanic slabs, taking into account the density-increasing effect of the Fe substitution for the Mg in the practical mantle Ol (Figure 17b). But at P ≥~24 GPa and with the existence of Brg in the lower mantle, the role that the Fo-III plays is much similar to that of the metastable Fo, providing a positive buoyancy force to slow down the sinking process of the slabs. substitution for the Mg in the practical mantle Ol (Figure 17b). But at P ≥ ~24 GPa and with the existence of Brg in the lower mantle, the role that the Fo-III plays is much similar to that of the metastable Fo, providing a positive buoyancy force to slow down the sinking process of the slabs. Figure 17. Calculated isothermal density-depth profiles at 900, 1200 and 1500 K for Fo (a); and Fo-III (b). As a reference, the PREM density model [98] is also sketched. The shown density difference (Δρ) at 450 km depth, 600 km or 720 km depth and 900 K is calculated as Δρ = 100 × (ρ-Fo − ρ-PREM)/ρ-PREM, where ρ-Fo stands for ρ-Fo in (a), and ρ-Fo-III in (b). The GGA results at different T are shown as red lines while the LDA results as blue lines. For the purpose of illustration, the density differences at certain depth obtained by the GGA method and the LDA method are slightly horizontally shifted.
Moreover, it is rather intriguing to probe any potential relation between the Fo/Fo-III phase transition and the deep earthquakes in the subduction zones [6]. Our results have suggested that this phase transition might take place at likely depths (~600 km) with a shear mechanism, result in large volume reductions (~12.6% at 20 GPa and 900 K for example), and thus bear the ingredients necessary for the generation of the deep earthquakes.

Conclusions
With extensive first-principles calculations, we have systematically examined the crystallographic characteristics, elastic properties and dynamic stabilities of the Fo phase and its three newly discovered polymorphs, namely the Fo-II, Fo-III and Fo-IV phases, and probed the phase relations of these Mg2SiO4 polymorphs with other common mantle minerals/mineral assemblages such as the Wds, Rwd, Brg + Pe and Aki + Pe. In addition, we have theoretically calculated the thermodynamic properties of the Fo-III phase.
The major conclusions that we have reached are: (1) Among the metastable polymorphs of the Fo, Fo-II, Fo-III and Fo-IV phase at P > ~22 GPa, the Fo-III phase is the best candidate which can directly form from the Fo phase at zero K, in good agreement with the observation in the hydrostatic compression experiments of Finkelstein et al. [37]. Both the Fo-II phase and the Fo-IV phase are dynamically instable for the P range of 0-100 GPa, and thus have negligible influence on the subduction process of the slabs; (2) The Fo phase directly transforms to the Fo-III phase at pressures immediately beyond its P stability field, in good agreement with the shockwave experiments of Newman et al. [39]. With a small temperature dependence (averagely −1.1 MPa/K for the T range of 300-1500 K), the transition P for the Fo/Fo-III phase transition at 300 K is estimated as ~23.75 GPa using the GGA method, which is presumably underestimated by ~3 GPa [66]. This polymorphic phase transition occurs by intracrystalline shear deformation of the Fo structure at low T [37]; (3) The Fo-III phase likely forms in the cold slabs subducted to the lower part of the MTZ. As a result, the metastable Ol wedge hypothesis may be applicable to the upper part of the MTZ only while the metastable Fo-III phase should be important in the cold slabs in the lower part Figure 17. Calculated isothermal density-depth profiles at 900, 1200 and 1500 K for Fo (a); and Fo-III (b). As a reference, the PREM density model [98] is also sketched. The shown density difference (∆ρ) at 450 km depth, 600 km or 720 km depth and 900 K is calculated as ∆ρ = 100 × (ρ -Fo − ρ -PREM )/ρ -PREM , where ρ -Fo stands for ρ -Fo in (a), and ρ -Fo-III in (b). The GGA results at different T are shown as red lines while the LDA results as blue lines. For the purpose of illustration, the density differences at certain depth obtained by the GGA method and the LDA method are slightly horizontally shifted.
Moreover, it is rather intriguing to probe any potential relation between the Fo/Fo-III phase transition and the deep earthquakes in the subduction zones [6]. Our results have suggested that this phase transition might take place at likely depths (~600 km) with a shear mechanism, result in large volume reductions (~12.6% at 20 GPa and 900 K for example), and thus bear the ingredients necessary for the generation of the deep earthquakes.

Conclusions
With extensive first-principles calculations, we have systematically examined the crystallographic characteristics, elastic properties and dynamic stabilities of the Fo phase and its three newly discovered polymorphs, namely the Fo-II, Fo-III and Fo-IV phases, and probed the phase relations of these Mg 2 SiO 4 polymorphs with other common mantle minerals/mineral assemblages such as the Wds, Rwd, Brg + Pe and Aki + Pe. In addition, we have theoretically calculated the thermodynamic properties of the Fo-III phase.
The major conclusions that we have reached are: (1) Among the metastable polymorphs of the Fo, Fo-II, Fo-III and Fo-IV phase at P >~22 GPa, the Fo-III phase is the best candidate which can directly form from the Fo phase at zero K, in good agreement with the observation in the hydrostatic compression experiments of Finkelstein et al. [37]. Both the Fo-II phase and the Fo-IV phase are dynamically instable for the P range of 0-100 GPa, and thus have negligible influence on the subduction process of the slabs; (2) The Fo phase directly transforms to the Fo-III phase at pressures immediately beyond its P stability field, in good agreement with the shockwave experiments of Newman et al. [39]. With a small temperature dependence (averagely −1.1 MPa/K for the T range of 300-1500 K), the transition P for the Fo/Fo-III phase transition at 300 K is estimated as~23.75 GPa using the GGA method, which is presumably underestimated by~3 GPa [66]. This polymorphic phase transition occurs by intracrystalline shear deformation of the Fo structure at low T [37]; (3) The Fo-III phase likely forms in the cold slabs subducted to the lower part of the MTZ. As a result, the metastable Ol wedge hypothesis may be applicable to the upper part of the MTZ only while the metastable Fo-III phase should be important in the cold slabs in the lower part of the MTZ. Due to its high bulk seismic velocity, the Fo-III phase may make significant contribution to the fast seismic velocity signals usually detected for the cold slabs subducted into the MTZ, but conventionally assigned to low T, the termination of the metastable Ol wedge, or both; (4) In contrast to the role that the metastable Fo plays in governing the geodynamic process of the cold slabs subducted beyond the 410 km depth, the Fo-III phase is not able to provide a positive buoyance force to slow down the subduction process in the lower part of the MTZ, but is able to do so in the lower mantle.
Further, there are several limitations to our proposal of a significant role of the metastable Fo-III in the deep-subducted cold slabs, including the unknown effect of water and the absence of the kinetics data on the phase transition of the Fo-III phase to other phase assemblages such as the Brg + Pe. All these issues should be addressed in future.