Using particle swarm optimization with density functional theory, we identify the positions of hydrogen in a hypothetical Mg-end-member of phase egg (MgSiO4H2) and predict the most stable crystal structures with MgSiO4H2 stoichiometry at pressures between 0 and 300 GPa. The particle swarm optimization method consistently and systematically identifies phase H as the energetically most stable structure in the pressure range 10–300 GPa at 0 K. Phase Mg-egg has a slightly higher energy compared to phase H at all relevant pressures, such that the energy difference nearly plateaus at high pressures; however, the combined effects of temperature and chemical substitutions may decrease or even reverse the energy difference between the two structures. We find a new MgSiO4H2 phase with the P43212 space group that has topological similarities to phase Mg-egg and is energetically preferred to phase H at 0–10 GPa and 0 K. We compute the free energies for phase Mg-egg, phase P43212, and phase H at 0–30 GPa within the quasi-harmonic approximation and find that the effect of temperature is relatively small. At 1800 K, the stability field of phase P43212 relative to the other polymorphs increases to 0–14 GPa, while pure phase Mg-egg remains energetically unfavorable at all pressures. Simulated X-ray diffraction patterns and Raman spectra are provided for the three phases. Additionally, the crystallographic information for two metastable polymorphs with the P1 space group is provided. Our results have implications for the deep hydrogen cycle in that we identify two novel potential carrier phases for hydrogen in the mantles of terrestrial planets and assess their stability relative to phase H. We determine that further experimental and computational investigation of an extended compositional space remains necessary to establish the most stable dense hydrated silicate phases.

The water content of the Earth’s mantle is poorly constrained. In addition to its significance for planetary habitability when it resides in surface reservoirs, water that remains in the interior may lower the melting temperature of rocks (Gaetani and Grove 1998), enhance the rate of phase transitions (Kubo et al. 1998), change the position of phase boundaries (Wood 1995), lower the mantle’s viscosity (Mei and Kohlstedt 2000), affect mantle dynamics (Nakagawa et al. 2015), and facilitate element transport (Kogiso et al. 1997). Improving our knowledge of the budget and history of hydrogen depends on knowing its storage and transport mechanisms, as well as its original sources. Although most of the common hydrous phases observed at near-surface conditions decompose at mantle pressures and temperatures, there are several potential reservoir phases. In the upper mantle, “water” may be present at low concentrations in nominally anhydrous phases, such as olivine, garnet, and pyroxene and, in the transition zone, wadsleyite and ringwoodite (Gasparik 1993; Inoue et al. 1995; Smyth and Kawamoto 1997; Pearson et al. 2014). Mineralogists have also identified a range of so-called dense hydrous magnesium silicates, such as phase D [MgSi2O4(OH)2] and superhydrous B [Mg12Si4O19(OH)2], that may survive dehydration reactions in subduction zones and therefore serve as efficient transport vectors for water through the upper mantle and perhaps into the lower mantle. In addition to water recycled at subduction zones, if primordial water has persisted throughout Earth’s history, it must have also resided in some host phases. Numerical models have suggested that an early terrestrial global magma ocean may have retained substantial H2O in the silicate melt, given the solubility beneath a thick, CO2-dominated atmosphere (Gaillard and Scaillet 2014; Bower et al. 2019; Solomatova et al. 2021). Upon crystallization of the magma ocean, a portion of that dissolved water budget may have become incorporated into dense hydrous magnesium silicates. Thus, to improve our understanding of both the current state and evolution of Earth’s lower mantle, it is necessary to quantify its water budget and characterize the possible phases that may or may not exist at those pressures and temperatures. Note that one may equally well discuss the budget of hydrogen; given the ready availability of oxygen throughout the mantle and surface reservoirs of Earth, hydrogen and water are interchangeable when discussing inventories.

Several dense hydrous magnesium silicate phases have been suggested as possible carriers of hydrogen in the mantle. In subduction zones, serpentine (the product of lower-temperature, near-surface alteration of peridotite minerals) decomposes into a mixture of phase A, enstatite, and fluid. Subsequent reactions of hydrous magnesium silicate phases with increasing depth depend on the water content (Ohtani et al. 2000) and also, most likely, on temperature, chemistry, and oxidation state. In wet regions of the transition zone, the major hydrogen-bearing phases are likely hydrous wadsleyite, hydrous ringwoodite, and superhydrous phase B. At the top of the lower mantle, phase D becomes the dominant water carrier. It was previously assumed that all hydrous phases decompose beyond ~45 GPa (Shieh et al. 1998). However, it has recently been proposed that polymorphs of MgSiO4H2 may be important carriers of water in the lower mantle at pressures exceeding 45 GPa (Tsuchiya 2013; Nishi et al. 2014; Bindi et al. 2014; Panero and Caracas 2017, 2020). One such polymorph, phase H, topologically equivalent to δ-AlOOH, was experimentally found to exist at pressures of 35–60 GPa and temperatures below ~1500 K (Ohtani et al. 2014). In agreement with the experiments, a computational study predicts that phase H exists up to 60 GPa at 1000 K, decomposing into bridgmanite + ice VII at higher pressures (Tsuchiya and Umemoto 2019). Another new phase, Mg-bearing phase egg [note that phase egg sensu stricto is a polymorph of AlSiO3OH (eggleton et al. 1978)], was synthesized at 24 GPa and 1673 K, suggesting that it may also be a potential carrier of hydrogen at lower mantle conditions (Bindi et al. 2020). The Al-end-member of phase egg decomposes to δ-AlOOH + SiO2 stishovite below ~1500 K and to Al2SiO4(OH)2 phase D + Al2O3 corundum + SiO2 stishovite above ~1500 K at pressures of about 17–24 GPa (Fukuyama et al. 2017). However, the fate of Mg-bearing phase egg is yet to be determined. Here we explore the structure and relative stability of the Mg-end-member of phase egg in the context of various polymorphs of MgSiO4H2 using computational methods.

The Crystal Structure AnaLYsis by Particle Swarm Optimization (CALYPSO) package (Wang et al. 2012) seeks optimal structures by minimizing enthalpy during structural evolution by particle swarm optimization. Structural relaxation and enthalpy calculations are conducted with external optimization codes, which may be based either on density functional theory or interatomic pair potentials. In this study, the CALYPSO method was used for two purposes: (1) to determine the hydrogen positions in the hypothetical Mg-end-member of phase egg (Bindi et al. 2020) and (2) to search for any other energetically stable crystal structures with MgSiO4H2 stoichiometry at pressures between 0 and 300 GPa. Structures were optimized with density functional theory using the VASP package (Kresse and Furthmüller 1996). The first-generation structures were produced randomly; then, half of each subsequent generation was generated through particle swarm optimization and the other half was generated randomly. Each calculation consisted of about 40–60 generations with a population size of 30 structures per generation. To determine the hydrogen positions in the Mg-end-member of phase egg (hereafter, “phase Mg-egg”), the atomic positions of magnesium, silicon, and oxygen were fixed to the positions determined from single-crystal X-ray diffraction experiments (Bindi et al. 2020). The hydrogen atoms were then inserted into the structure either randomly or through particle swarm optimization, after which the structures were completely relaxed with VASP, allowing for all atoms to reach their equilibrium positions. The search was repeated eight times at pressures between 18 and 26 GPa to ensure reproducibility and self-consistency. To search for additional energetically stable structures with MgSiO4H2 stoichiometry, we conducted a full structure search with no pre-defined atomic positions at pressures of 0, 5, 10, 15, 20, 80, 135, 200, and 300 GPa, varying the number of formula units per unit cell from Z = 1 to Z = 4.

Ab initio calculations were performed using the projector-augmented wave (PAW) method (Blöchl 1994) implemented in VASP. The generalized gradient approximation (GGA) (Perdew et al. 1996) was used to approximate the exchange correlation terms. A plane-wave energy cut-off of 600 eV was used to ensure excellent convergence in volume and total energy (Online Materials1 Fig. OM1), and a k-point grid of <0.04 Å−1 was required to refine the transition pressures. The convergence criteria for electronic self-consistency and ionic relaxation loop were set to 10−5 and 10−4 eV, respectively. We ensured that forces acting on all relaxed atoms were less than 0.01 eV/Å. As discussed below, we found phase Mg-egg, phase H, a new low-energy structure with space group P43212, and two additional metastable but competitive candidate structures. These five structures were relaxed between 0 and 140 GPa to compare their relative stabilities. It is possible that, with certain cation substitutions, the competing candidate phases may be observed experimentally. Hence, all the studied structures, including those found to be metastable in the pure system, are provided in the Online Materials crystallographic information file (CIF1). X-ray diffraction patterns were simulated in the VESTA program (Momma and Izumi 2008).

The PAW method with the GGA approximation has successfully predicted the physical and elastic properties of a wide range of geological materials, including dense hydrous silicates (Li et al. 2006; Panero and Caracas 2017, 2020; Caracas and Panero 2017), silicate perovskites (Jung and Oganov 2005; Stixrude et al. 2007), pyroxenes (Yu et al. 2010), and carbonates (Oganov et al. 2008; Arapan and Ahuja 2010; Solomatova and Asimow 2017). The GGA approximation is known to slightly underestimate cohesive energies, resulting in the underestimation of bulk moduli and overestimation of volumes and phase transition pressures. While the experimentally synthesized phase H is characterized by disordered hydrogen atoms occupying half the 4g positions (Bindi et al. 2014), in our density functional theory calculations, the hydrogen positions are ordered, which reduces the symmetry from the experimentally determined orthorhombic Pnnm space group to the monoclinic Pm space group. In comparison, phase δ-AlOOH, the topologically equivalent structure to phase H, is characterized by ordered hydrogen at ambient pressure and undergoes an order-to-disorder transition at ~10 GPa, which results in an increase in symmetry from the P21nm to Pnnm space group and a change in compressional behavior (Sano-Furukawa et al. 2009, 2018; Kuribayashi et al. 2014; Ohira et al. 2019).

We compare the calculated lattice parameters, unit-cell volumes, and bond lengths of phase H to those experimentally measured by Bindi et al. (2014). Despite the difference in space groups between the computationally and experimentally studied phase H, we find that the GGA approximation overestimates lattice parameters by 0.2–1.4% and unit-cell volume by 1.8% (Online Materials1 Table OM1), while the O-H bond lengths are overestimated by about 3% (Online Materials1 Table OM2). We compare the equation of state parameters that result from fitting a Birch-Murnaghan equation of state to calculations covering the same pressure interval as in situ experimental volume data on phase H (Nishi et al. 2018), finding that the GGA method underestimates the bulk modulus by about 5–7% (Online Materials1 Table OM3). The higher degree of compressibility can be attributed to the lower symmetry of the structure and difference in hydrogen ordering. Although experiments have not yet constrained the hydrogen bond symmetrization pressure in phase H, previous simulations on AlOOH find that the PAW method with the GGA approximation overestimates the phase transition pressure of α-AlOOH to δ-AlOOH by <1 GPa (Li et al. 2006).

The method of lattice dynamics within the quasi-harmonic approximation was implemented using the finite-displacement method in the PHONOPY package (Togo and Tanaka 2015) to calculate the vibrational zero-point energies and thermal effects on the relative phase stabilities of the three most stable phases (i.e., phase Mg-egg, phase H, and the newly found P43212 phase) at 0–30 GPa. Supercells of 1×2×1 (64 atoms), 2×2×3 (96 atoms), and 2×2×1 (128 atoms) were created for phase Mg-egg, phase H and phase P43212, respectively, which proved to be sufficiently large for an energy convergence of about 10−4 to 10−5 eV/atom. Each supercell contains one atomic displacement with a magnitude of 0.01 Å and the number of supercells used for each phase depends on the space group (ranging between 16 and 48 calculations per phase at each pressure). Force constants were calculated for each displaced supercell using VASP. The convergence criterion for the ionic relaxation was decreased to 10−7 eV for the geometric relaxation of the undisplaced cells, and the convergence criterion for electronic self-consistency was decreased to 10−8 eV for all calculations related to the zero-point energy and thermal calculations. The resulting force constants were then used to calculate phonon-related properties (i.e., the zero-point energy and vibrational Helmholtz free energy at finite temperature), from which we derived the Gibbs free energies as a function of pressure and temperature.

The Raman spectra were computed as a function of pressure within the density functional perturbation theory (Baroni and Resta 1986; Baroni et al. 2001; Gonze et al. 2005), as implemented in ABINIT (Veithen et al. 2005; Gonze et al. 2002, 2009). We used the relaxed structures obtained with VASP, with norm-conserving pseudo-potentials and a k-point grid density of ~0.04 Å−1. The simulations use parameters similar to those employed to obtain the Raman spectra for minerals in the WURM database (Caracas and Bobocioiu 2011), which give reliable results for both anhydrous (McKeown et al. 2010) and hydrous minerals (Bobocioiu and Caracas 2014).

Structure searches

We identified the positions of the hydrogen atoms in the Mg-end-member of phase egg with the P21/n space group using constrained structure searches with the particle swarm optimization method (Fig. 1; Online Materials1 CIF). In phase Mg-egg (MgSiO4H2), half of the hydrogen positions are equivalent to those of Al-end-member phase egg (AlSiO4H) (eggleton et al. 1978; Schmidt et al. 1998; Schulze et al. 2018). We looked specifically at the question of whether the Mg-egg structure contains H2O groups; according to our results, it does not—every hydrogen atom is separately bonded with a single oxygen atom. With the full structure of phase Mg-egg identified, it is possible to compare the relative energy to its polymorph, phase H (Tsuchiya 2013; Nishi et al. 2014; Bindi et al. 2014). At 0 K, phase H is preferred to phase Mg-egg between 2.5 GPa and to at least 300 GPa, the maximum pressure explored with CALYPSO in this study.

Our structure search identified two energetically stable crystal structures with MgSiO4H2 stoichiometry in the pressure range from 0 to 300 GPa. At pressures between 0 and 10 GPa, the newly discovered structure with the P43212 space group is the energetically preferred phase (Fig. 2). The P43212 phase has topological similarities to P21/n Mg-egg. In Mg-egg, every two MgO6 polyhedra share edges, resulting in Mg2O10 dimers, which are all interconnected through their corners. Additionally, each Mg2O10 dimer shares four edges with two SiO6 octahedra and two corners with another two SiO6 octahedra. In the P43212 phase, the MgO6 octahedra share only corners with each other; every two MgO6 octahedra are at an acute angle of 56° (i.e., “dimer precursors”). Each MgO6 octahedron shares an edge with a SiO4 tetrahedron and three corners with three SiO4 tetrahedra. Evidently, the P43212 phase is a more open, less dense structure than Mg-egg. We provide simulated X-ray diffraction patterns of phase H, phase Mg-egg and phase P43212 at 0 GPa (Fig. 3).

Two additional triclinic phases with the P1 space group were discovered at low pressure. Although the energy differences are small, both triclinic structures have slightly higher energies than the P43212 phase at all pressures explored in this study. Due to the possibility that relative energies may change at high temperature and with chemical substitutions, we provide the structures of the two P1 phases, hereafter referred to as P1-a and P1-b (Fig. 1). The two P1 phases are characterized by corner-sharing SiO4 tetrahedra and MgO6 octahedra with H2O groups attached to the elongated axes (P1-a) or short axes (P1-b) of the MgO6 octahedra (Fig. 1).

At all pressures examined between 10 and 300 GPa, the particle swarm optimization method consistently and systematically identifies phase H as the energetically most stable structure. Phase Mg-egg has a slightly higher energy compared to phase H between 2.5 and 300 GPa; the energy difference between the two phases nearly plateaus at roughly 60 meV/atom at lower-mantle pressures (Fig. 2). This implies that the molar volumes of Mg-egg and phase H are similar and converge toward the same value, whereas P43212, P1-a, and P1-b are all much less dense. Hence, it is possible, given changes in relative energies due to temperature or chemistry, that phase Mg-egg could become energetically favorable relative to phase H at pressures of Earth’s lower mantle or the mantles of super-Earth exoplanets.

Standard density functional theory calculations do not include the kinetic energy corresponding to zero-point motion. The zero-point motion is larger for light elements such as hydrogen (Matsushita and Matsubara 1982; Natoli et al. 1993; Rivera et al. 2008), affecting the nature of hydrogen bonds and the structure of ice polymorphs (Benoit et al. 1998; Herrero and Ramírez 2011), such that zero-point energy may be significant for dense hydrous silicate phases (Tsuchiya et al. 2005). We investigated this effect for the MgSiO4H2 polymorphs, and we find that the zero-point energy corrections are small, e.g., ~1.5 meV/atom. The resulting Gibbs free energy differences are 0.5–14 meV/atom, which results in <1 GPa differences in phase transition pressures, and so do not affect the relative phase stabilities (Fig. 2 inset).

We applied the quasi-harmonic approximation to estimate the Gibbs free energies of phase H, phase Mg-egg, and phase P43212 along an isotherm at 1800 K, the approximate temperature of the deep upper mantle and transition zone. We find that the thermal effects are relatively small and do not change the sequence of stable MgSiO4H2 polymorphs with pressure. The stable energy crossover between phase P43212 and phase H shifts from 10 GPa (0 K) to 14 GPa (1800 K), while the metastable energy crossover between phase Mg-egg and phase H shifts from 2.5 GPa (0 K) to 9 GPa (1800 K). However, phase Mg-egg remains energetically disfavored relative to other polymorphs at all pressures and temperatures considered.

Equations of state

The pressure-volume results of our calculations for phase Mg-egg, phase P43212, and phase H were fitted with finite-strain equations of state (Fig. 4). The graph of reduced pressure (F) vs. Eulerian finite strain (f) reveals the level of Taylor expansion necessary to accurately fit equation-of-state data (Angel 2000). Both Mg-egg and P43212 phases show concave-down parabolic curves in F-f space that call for the use of the fourth-order Birch-Murnaghan equations of state (Fig. 4 inset). The situation for phase H is complex; the Pm structure of phase H undergoes a second-order phase transition to the P2/m structure at about 30 GPa due to H-O bond symmetrization (Tsuchiya and Mookherjee 2015; Lv et al. 2017). The phase transition is subtle in pressure-volume space but manifests as a clear break in slope at f = 0.05 in the F-f plot. The linearity of the low-pressure region warrants the use of a third-order Birch-Murnaghan equation of state between 0 and 25 GPa, while the concave-down behavior of the high-pressure region warrants the use of a fourth-order Birch-Murnaghan equation of state between 35 and 140 GPa. Phase H experimentally synthesized at 1273 K and studied at 300 K has Pnnm symmetry (Bindi et al. 2014), with disordered hydrogen occupying half the 4g positions and octahedral sites occupied equally by magnesium and silicon. On the other hand, in 0 K density functional theory studies, phase H becomes ordered, which reduces the symmetry to Pm. The lowest-energy ordering of the sites was chosen for the purpose of this study, such that the hydrogen positions correspond to the “hydrogen off-centered 1” (HOC1) arrangement in δ-AlOOH (Tsuchiya et al. 2002, 2008).

In the resulting monoclinic structure of phase H, the calculated results of density functional theory relaxation indicate a second-order transition at about 30 GPa from a Pm space group with asymmetric hydrogen positions to a P2/m structure with symmetrized hydrogen bonds, in agreement with previous density functional theory calculations (Tsuchiya 2013); the Pm and P2/m structures are otherwise topologically equivalent (Tsuchiya and Mookherjee 2015). Bond symmetrization has been experimentally observed in δ-AlOOH at ~18 GPa through neutron diffraction (Sano-Furukawa et al. 2018) and in phase D at ~40 GPa through anomalous axial and volumetric compression (Shinmei et al. 2008; Hushur et al. 2011). In phase D, bond symmetrization was experimentally observed at ~40 GPa through anomalous axial and volumetric compression (Shinmei et al. 2008; Hushur et al. 2011). Bond switching (O-H · · · O to O · · · H-O) was observed in the Al-egg phase (AlSiO3OH) at 14 GPa using a combination of IR spectroscopy and single-crystal X-ray diffraction (Liu et al. 2021). We do not observe bond switching in the Mg-egg phase; the additional hydrogen atoms prevent the occurrence of bond switching since the potential recipient oxygen atom is already bonded to a hydrogen atom. Meanwhile, the additional hydrogen atoms in the Mg-egg structure relative to the Al-egg structure do not experience bond switching due to the different local environment.

The symmetrization of the O-H bonds in phase H results in an increase in the zero-pressure bulk modulus K0T from 133(1) to 174(4) GPa and a decrease in the zero-pressure unit-cell volumes V0 from 7.41(1) to 7.25(2) Å3/atom (Table 1). The results of an X-ray diffraction study on phase H at 0 and 35–60 GPa are consistent with a change in compressibility at 30 GPa due to bond symmetrization (Nishi et al. 2018); however, more detailed experimental studies at pressures corresponding to the onset of bond symmetrization are needed. In comparison to phase H, we find that phase Mg-egg is more compressible with a bulk modulus of 113(1) GPa, consistent with convergence between the volumes of phases H and Mg-egg at high pressure and the leveling-off of the energy difference between them. Phase P43212 has the lowest bulk modulus, 76(1) GPa, and the largest initial volume, 9.08(1) Å3/atom, indicating that it is a low-pressure phase that is not expected to exist at lower-mantle pressures. Previous density functional theory calculations on phase H predict bulk moduli of 151.9 and 185.8 GPa and volumes of 7.34 and 7.20 Å3/atom for the asymmetric and symmetric hydrogen bond structures, respectively (Tsuchiya 2013). The differences in the fitted equations of state parameters for phase H between this study and Tsuchiya (2013) are consistent with expectations for the comparison of density functional theory calculations with core electron treated with norm-conserving pseudopotentials (Tsuchiya 2013) and with the PAW method (this study). We confirmed this by relaxing the reported unit cell of Tsuchiya (2013) with the PAW method, resulting in an increase of the unit-cell volume from 7.34 to 7.41 Å3/atom and an increase in the γ angle from 93.25° to 93.31°.

Raman vibrational modes

We provide the simulated Raman spectra of phase H, phase Mg-egg and phase P43212 at 0–40 GPa (Fig. 5). The spectra are computed for static structures, the main contribution of the anharmonicity being the broadening of the peaks. The vibrational frequencies of Raman bands are systematically related to features of the O-H bonds in the studied phases. The stretching frequencies of the O-H bonds depend on both O · · · O distance and H-O-O angle in the O-H · · · O bond environment (Hofmeister et al. 1999); an H-O-O angle of 0° corresponds to a linear O-H…O arrangement. Generally, larger O · · · O distances and H-O-O angles are expected to result in higher O-H stretching frequencies.

In the phase H structure, there are two types of O-H bonds, which are both oriented toward hydrogen-free oxygen atoms. The O-H · · · O bond environments are characterized by very small H-OO angles of 1.0° (H1-O3-O1) and 2.3° (H2-O2-O4) and moderate O · · · O distances of 2.58 Å (O3 · · · O1) and 2.51 Å (O2 · · · O4) at 0 GPa (Figs. 6a and 7). The restricted O-H environment in phase H results in low frequencies of the O-H stretching modes of about 2380 and 2425 cm−1 (Fig. 5a). By 20 GPa, the environment becomes even more restrictive with H-O-O angles of 0.2° (H1-O3-O1) and 0.7° (H2-O2-O4) and O · · · O distances of 2.41 Å (O3 · · · O1) and 2.39 Å (O2 · · · O4), decreasing the intensity of the O-H stretching modes and shifting them to 1570 and 1750 cm−1. By 30 GPa, the phase H structure achieves bond symmetrization and the OH-stretching modes shift to 1550 and 1710 cm−1 with intensities that are several orders of magnitude lower and so are not visible in Figure 5a. All calculated Raman vibrational modes and O-H · · · O bond environment information can be found in Online Materials1 datasets 1 and 2, respectively.

Phase Mg-egg is characterized by two types of O-H · · · O bond environments. The first type of O-H bonds are oriented toward a hydrogen-free oxygen atom with an H1-O1-O3 angle of 5.1° and an O1 · · · O3 distance of 2.45 Å, while the second type of O-H bonds are oriented toward a hydrogen-bonded oxygen atom with an H2-O4-O1 angle of 13.5° and an O1 · · · O4 distance of 2.89 Å at 0 GPa (Fig. 6b). The more restricted local environment of the first type of O-H bonds results in a lower frequency stretching vibration compared to the second type. The calculated Raman spectrum for phase Mg-egg at 0 GPa shows vibrational modes at ~2850 and ~3475 cm−1 (Fig. 5b), corresponding to the first and second type of O-H bonds, respectively. With increasing pressure, the O-H distances and H-O-O angles steadily decrease, resulting in a gradual shift of the O-H stretching modes to lower frequencies. Unlike phase H and phase Mg-egg, phase P43212 contains only one type of O-H bond. In phase P43212, the O-H bonds are all oriented toward hydrogen-free oxygen atoms at an H1-O2-O1 angle of 6.6° and an O1 · · · O2 distance of 2.75 Å at 0 GPa (Fig. 6c), resulting in intermediate vibrational frequencies at ~3000 cm−1 (Fig. 5c). With increasing pressure, the O1 · · · O2 distance and H1-O2-O1 angle decrease, shifting the O-H stretching modes to lower frequencies.

An experimental study on δ-AlOOH, the topologically equivalent structure to phase H, found that the O-H stretching modes exist between 2100 and 2700 cm−1 at 0 GPa (Ohtani et al. 2001). With ab initio calculations, Tsuchiya et al. (2008) found that the O-H stretching modes range from 2000 to 2700 cm−1 and that the number and frequency of the vibrational modes depend on the degree of hydrogen disorder and the O · · · O distances. It is also likely that the local environment of hydrogen would be affected by the chemistry of the mineral phase. Our predicted O-H vibrational frequencies of 2380–2425 cm−1 for Mg-end-member phase H exist within the range of frequencies observed in δ-AlOOH. Density functional perturbation theory calculations from Lv et al. (2017) on phase H using the CASTEP software package (Segall et al. 2002) produce very similar Raman spectra to those of this study (Online Materials1 Fig. OM2). The low O-H stretching vibrational frequencies observed in phase H and δ-AlOOH are characteristic of constrained O-H environments. Previous computational and experimental studies have reported similarly anomalous low O-H stretching frequencies in other hydrous phases: ~2000 cm−1 for ice X extrapolated to 25 GPa (Caracas 2008), ~2700 cm−1 for ice VII at 25 GPa (Hernandez and Caracas 2018), and 2500–3000 cm−1 for dimer acids at ambient pressure (Chen et al. 2017), where the O-H environments are highly constrained.

In this study, simulations based on the particle swarm optimization method, density functional theory, and the quasi-harmonic approximation were used to (1) identify the positions of hydrogen in phase Mg-egg, (2) identify additional candidate MgSiO4H2 structures, and (3) examine the effect of temperature on relative phase stabilities. Identifying the atomic coordinates of the hydrogen atoms in phase Mg-egg is fundamental to characterizing its structure and will facilitate future experimental and computational studies; despite the addition of H necessary to charge-balance the substitution of Mg for Al in phase egg, the structure does not develop H2O groups. We find that pure end-member Mg-egg does not have a stability field in the range of 0–300 GPa and 0–1800 K. With particle swarm optimization, we identified three other candidate MgSiO4H2 polymorphs: two structures with the P1 space group that is never energetically favored relative to other polymorphs and a structure with the P43212 space group that is preferred over phase H at 0–14 GPa at 1800 K. It is possible that phase P43212 may be stable at certain low-pressure conditions within the Earth’s mantle, particularly in subducting slabs just before reaching the transition zone. However, note that this study is restricted to the relative stability of phases with MgSiO4H2 stoichiometry only. In particular, we have not examined the effects of decomposition into less hydrous Mg-silicates and free H2O, due to the complex nature of the phase diagram of H2O and the sensitivity of computed physical and thermodynamic properties of H2O to various parameters of density functional theory models (e.g., representation of Van der Waals forces, flavor of pseudopotentials, the inclusion of zero-point motions, and inherent anharmonicity). A thorough density functional theory study of the decomposition of hydrous magnesium silicates into ice and anhydrous magnesium silicates has been performed (Tsuchiya 2013), predicting that phase H would be stable at 30–52 GPa at 0 K and 50–52 GPa at 1600 K (corresponding to a cold geotherm), decomposing above about 1900 K. Considering that phase Mg-egg and phase P43212 are both less stable than phase H above ~15 GPa, it is unlikely that the Mg-end-members of these phases would exist at any depth in the mantle. Furthermore, the relative stability of phase H may increase at high temperatures due to the configurational entropy from the disordering of Mg and Si on the octahedral sites.

Although pure phase Mg-egg is never favored relative to its polymorphs at the conditions explored in this study, chemical substitutions may decrease the energy difference between phase Mg-egg and phase H. Indeed, an (Mg,Al)-bearing phase egg solid solution was successfully synthesized and recovered from high-pressure experiments likely to have approached equilibrium. The observed synthetic composition was found to be Al0.65Mg0.35SiO4H1.35 (Bindi et al. 2020). The presence of aluminum in MgSiO4H2 may significantly alter the relative phase stabilities and their high-pressure behavior. It is thus possible that (Mg,Al)-bearing phase egg is thermodynamically favored over (Mg,Al)-bearing phase H or the assemblage of Mg-end-member phase H and Al-end-member phase egg under some natural conditions. Additional calculations are needed to examine the effect of aluminum substitution on the relative stability of the MgSiO4H2 polymorphs.

Dense hydrous magnesium silicate phases are important hydrogen carriers and/or reservoirs of hydrogen in the deep mantle, affecting the melting temperature of surrounding rocks, changing phase relations, and affecting mantle dynamics. The sequence of hydrous phases that exist in the mantle depends on the local geotherm, water content, and overall chemistry (Ohtani et al. 2000). Here we have examined the polymorphs of the end-member MgSiO4H2. We confirm that phase H is the preferred crystal structure with the MgSiO4H2 stoichiometry at pressures corresponding to the Earth’s transition zone and entire lower mantle. In fact, phase H continues to be the preferred polymorph of MgSiO4H2 to at least 300 GPa, indicating that it may be a candidate phase in water-bearing exoplanetary interiors if other hydrous magnesium silicate assemblages are not favored. Although we found that the Mg-end-member of phase egg is never energetically favorable relative to its polymorphs and phase P43212 is likely unstable relative to multiphase assemblages, it is possible that certain chemistries and water contents would favor them with respect to decomposition products and other dense hydrous magnesium silicate phases. Further investigations into the compositional space of dense hydrated silicates, e.g., the relative stabilities of solid solutions of phases egg, H, and P43212 along the MgSiO4H2-AlSiO3(OH) join, are necessary to resolve these questions.

This research was supported in part by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 681818–IMPACT to R.C.). We acknowledge access to the Irene supercomputer through the PRACE computing grant no. RA4947 and the PSMN center of ENS Lyon.

Deposit item AM-22-57937, Online Materials. Deposit items are free to all readers and found on the MSA website, via the specific issue’s Table of Contents (go to The CIF has been peer-reviewed by our Technical Editors.
Open access: Article available to all readers online.