Crystal structure, equation of state, and elasticity of phase H (MgSiO4H2) at Earth’s lower mantle pressures

Dense hydrous magnesium silicate (DHMS) phases play a crucial role in transporting water in to the Earth’s interior. A newly discovered DHMS, phase H (MgSiO4H2), is stable at Earth’s lower mantle, i.e., at pressures greater than 30 GPa. Here we report the crystal structure and elasticity of phase H and its evolution upon compression. Using first principles simulations, we have explored the relative energetics of the candidate crystal structures with ordered and disordered configurations of magnesium and silicon atoms in the octahedral sites. At conditions relevant to Earth’s lower mantle, it is likely that phase H is able to incorporate a significant amount of aluminum, which may enhance the thermodynamic stability of phase H. The sound wave velocities of phase H are ~2–4% smaller than those of isostructural δ-AlOOH. The shear wave impedance contrast due to the transformation of phase D to a mixture of phase H and stishovite at pressures relevant to the upper part of the lower mantle could partly explain the geophysical observations. The calculated elastic wave velocities and anisotropies indicate that phase H can be a source of significant seismic anisotropy in the lower mantle.

structure consisting of magnesium (Mg) and silicon (Si) atoms in an octahedral coordination and is similar to the crystal structure of δ -AlOOH 17,18 where aluminum (Al) atoms occur in octahedral sites. Owing to the possibility of the configurational disorder of the Mg and Si atoms in the octahedral sites, the space group symmetry of phase H was subsequently refined using the single-crystal X-ray diffraction method 19 .
Despite being an important candidate for transporting water into the deep Earth, the crystal structure and elasticity of phase H at high pressures are unknown. In this study, using first principles simulations, we report the crystal structure of phase H including the effect of ordering and disordering of Mg and Si atoms in the octahedral sites. We also determine the equation of state, full-elastic constant tensor, and the elastic anisotropy of phase H at pressures relevant to Earth's lower mantle pressure conditions.

Results
Energetics calculations, considering the effects of Mg and Si order-disorder reveal that the model-1 with P2/m space group symmetry has the lowest enthalpy across all pressures relevant for Earth's interior ( Fig. 1). Model-2 with space group symmetry of P2 1 2 1 2 show slightly greater enthalpy (~0.13 eV/Z, where Z = formula unit in the unit-cell) than model-1. The pressure-volume results for phase H are well described by the finite strain formulation of Birch-Murnaghan 20 (Fig. 2) with a zero pressure bulk modulus K O of 147.5 (± 6.8) GPa, the pressure derivative of the zero pressure bulk modulus ′ K O of 4.9 (± 0.2), and a zero pressure unit-cell volume V O of 58.9 (± 0.2) Å 3 . The predicted lattice parameters are in good agreement with the experimental results 18,19 (Fig. 2 O O at ~30 GPa. The hydrogen bond symmetrization has also been predicted for an isostructural dense hydrous phase such as δ -AlOOH phase using first principles simulations 21,22 and later confirmed by experiments 23,24 . The components of the full elastic constant tensor of phase H increase with pressure ( Table 2, Fig. 3). The principal elastic constants C 11 and C 22 , the off-diagonal elastic constants C 23 , C 16 , and C 26 , and the shear elastic constant C 66 show anomalous increase at ~30 GPa (Fig. 3). The pressure dependent anomalous behavior in the elastic constant is predicted for both the ordered (model-1) and disordered (model-2) phase H structures (Supplementary Table 2). The likely cause for such pressure dependent anomalous behavior in elasticity is the symmetrization of hydrogen bonds. Similar anomalous behavior in elasticity for the isostructural δ -AlOOH phase have also been explained by hydrogen bond symmetrization 21 . In the present study, although, the anomalous behavior occurs at ~30 GPa, it is possible that the quantum and thermal vibrational effects may affect the pressure where hydrogen bond symmetrizes and also the onset of anomalous behaviors 25 . At pressures beyond the hydrogen bond symmetrization, all the elastic constants stiffen steadily upon compression without further pressure dependent anomalous behavior. It is interesting to note that the predicted C 33 is smaller than C 11 Table 2). The model-2 type phase H is likely to be mechanically unstable at low pressure conditions since the shear elastic constant, C 55 of < 0 at ~0 GPa.  18 and (grey filled symbols) 19 . The bulk (K), shear (G) moduli, primary (V P ), and shear (V S ) velocities for phase H increase upon compression (Fig. 4). The bulk sound velocity (i.e., V Φ = (K/ρ ) 0.5 ) is in very good agreement with the recent study on shock wave experimental studies on DHMS phases 26 (Supplementary Fig. 1). The discontinuous behavior of K at pressures of ~20-30 GPa is related to the anomalous increase of the principal elastic constants, C 11 and C 22 , which is in turn related to the hydrogen bond symmetrization. The effect of Mg-Si order disorder do not show appreciable changes in the bulk and shear moduli throughout the range of pressures explored in this study (Fig. 4). Since phase H and δ -AlOOH are isostructural at conditions relevant to the lower mantle, these two phases could show appreciable solid solution through Mg + Si = 2Al substitutions. This has been observed in recent experimental studies 18,27 . We notice that in comparison to the aluminous end member δ -AlOOH phase, the primary (V P ), and shear (V S ) velocities are 3% and 4% slower for phase H (Fig. 4), respectively. However, the primary (V P ) sound velocity for the  volumetrically dominant lower mantle phase, bridgmanite, lies in between that of δ -AlOOH and phase H. The shear sound velocity (V S ) for bridgmanite lies in between that of δ -AlOOH and phase H for a pressure range of 30 to 60 GPa but is greater than both δ -AlOOH and phase H above 60 GPa (Fig. 4). Hence, lower velocities in the deep mantle could be due to a combination of iron enrichment, higher temperatures, and the presence of deeply subducted hydrous phases.

Discussions
The sound wave velocities vary as a function of the propagation direction. The anisotropy in the sound wave velocities could be calculated by solving the Christoffel's equations, det|c ijkl n j n l -ρ V 2 δ ik | = 0, where n, ρ , V, and δ ik are the propagation direction, density, wave velocity, and Kronecker delta, respectively 28 . The calculated shear wave (V S ) polarization anisotropy (AV S ) is defined as AV S = 100 × (V S1 -V S2 )/V S and the azimuthal anisotropy for primary wave (AV P ) is defined as A P = 100 × (V Pmax -V Pmin )/V P . The polarization anisotropy AV S increases with pressure, whereas the azimuthal anisotropy AV P remains unchanged beyond the hydrogen bond symmetrization pressure ~30 GPa. The azimuthal anisotropy AV P is reduced for the model-2 type disordered structured phase H (AV P ~ 18%) compared with that of model-1 (AV P ~ 32%) (Fig. 4, Supplementary Figs 2 and 3), it is still higher than that of bridgmanite (AV P ~ 12%) 29 , and post perovskite (AV P ~ 15%) 30 at lower mantle pressures of 100 GPa. There is no major change of the fast and slow directions of the velocities of phase H between 40 and 100 GPa. The core mantle boundary region is known to be seismically anisotropic and the horizontally polarized shear waves (V SH ) propagate faster than vertically polarized ones (V SV ) 31 . The calculated polarization anisotropy indicates that if the lattice preferred orientation of phase H is developed by aligning the c-axis vertically, the high polarization anisotropy with V SH > V SV could partly explain the observed seismic anisotropy at the bottom of lower mantle.
We also estimated the acoustic impedance contrast (Δ I/I) due to the decomposition of phase D to a mixture of phase H and stishovite at pressure conditions relevant for the upper part of the lower mantle (Fig. 5). Decomposition reaction involving the DHMS phases could partly explain observed seismic impedance contrast at a mid-mantle depths of ~1000 km 32 . If the impedance contrast is indeed partly related to the decomposition of phase D to phase H and stishovite, it would mean that a significant amount of water could be transported into the deep Earth through DHMS phases such as phase D and phase H down to lower mantle depths.

Methods
We used first principles calculations based on the density functional theory to predict the structure, equation of state, and elasticity of phase H. We used generalized gradient approximation (GGA) 33   the description of exchange-correlation functional. Norm-conserving pseudopotentials 34 have been employed to describe the ionic core potentials of silicon (Si), oxygen (O), and hydrogen (H), whereas the magnesium (Mg) pseudopotential is generated by the method of U. von Barth and R. Car 35 . The semi-core p-electrons are not included in the Mg psedopotential. These potentials were extensively tested in previous studies [35][36][37] . GGA has been successfully used in predicting high-pressure behavior of hydrous phases 38,39 and have been tested experimentally 40 . All structural parameters are fully relaxed at a static 0 K and 0-100 GPa by the damped variable cell shape molecular dynamics method implemented in the Quantum-Espresso codes 41 until residual forces become less than 1.0 × 10 −5 Ry/au. The electronic wave function is expanded in plane waves using a kinetic energy cutoff of 80 Ry. The irreducible Brillouin zone of the phase H structure is sampled on a 4 × 4 × 6 Monkhorst-Pack mesh 42 . In addition to the unit cell calculation of phase H, we also conducted the supercells in order to estimate the effect of disordering between Mg and Si onto the cell parameters. K-points in those supercells are sampled on a larger mesh in order to achieve the k-point sampling, which is equivalent to that for the unit cell in reciprocal space. The elastic constants are determined by using the stress-strain relations 43 . The magnitude of all applied strains was ± 1%. The linear relation was ensured for this strain range. Recent first principles simulations reported a fully optimized crystal structure of phase H 17 with a monoclinic symmetry (γ ~ 91° at 40 GPa) with space group P2/m that was slightly distorted from orthorhombic symmetry. Experimental studies reported that the structure of phase H is orthorhombic with space group P2 1 nm 18 and more recently based on single-crystal X-ray diffraction, a space group of Pnnm 19 (CaCl 2 type structure) has been proposed. Although, the polyhedral frameworks of crystal structures proposed by all these studies are similar, the Mg and Si were found to be disordered in the octahedral sites 19 . In order to mimic the effect of disorder of the Mg and Si in the octahedral sites, we used distinct structural models and evaluate their relative energetics (Fig. 1). Figure 1 shows the calculated model crystal structures. All these models are derivative of the previous first principles simulations with P2/m 17 space group symmetry. Model-1 corresponds to the ordered Mg and Si crystal structure (cell size a × b × c) 17 , Model-2 (a × b × 2c supercell) where Mg and Si atoms are alternatively arranged along the c-axis resulting in an orthorhombic space group symmetry P2 1 2 1 2, and Model-3 (2a × 2b × c supercell) where Mg and Si atoms are alternatively placed along the a-and b-axes with monoclinic space group symmetry P2/m (Fig. 1). We have fully optimized the cell parameters of these model structures at 0-100 GPa. We also calculate the full elastic constant tensor for model-1 and model-2 (Supplementary  Table 1). In order to compare the elasticity of phase H with that of the major lower mantle phase, we have also calculated the elastic constants of bridgmanite using the same pseudopotentials as those used for the calculation of phase H (Supplementary Fig. 1).