Iron Hydride in the Earth's Inner Core and Its Geophysical Implications

Hydrogen is potentially a key light element in the Earth's core. Determining the stability of iron hydride is essential for Earth's core mineralogy applications. We investigated the thermal stabilities of a range of Fe‐H binaries at core P‐T conditions. It is concluded that face‐centered cubic phase FeH is stable in the Earth's inner core. The high mobility of hydrogen in the Fe lattice suggests that hydrogen is transferred to a superionic state under the inner core conditions, where the superionic state transfer temperature of H in Fe fcc lattice is ∼500 K higher than that in hcp Fe system. The H concentration in the inner core is estimated to be ∼0.92 wt% to explain its density deficit, this value was further constrained to ∼0.21 wt% by matching the density jump at the inner‐core boundary. H alongside other light elements are required to account for the geophysical observations of the Earth's inner core.

However, there are still some controversies within the literature about the effect of H on the Earth's inner core. The compressional experiments on Fe-Si-H alloys conducted by Tagawa et al. (2016) suggested that 0.47 wt% H is required to match the density deficit of the inner core, which is higher than the recent estimate of the H content in the inner core (∼0.2-0.3 wt%) by Thompson et al. (2018). Bazhanova et al. (2012) predicted the structures of Fe-H binaries by variable-composition high pressure structure prediction methods and they speculated that the inner core would contain up to 0.4-0.5 wt% H. Alternatively, Caracas (2015) calculated the seismic properties of solid hcp Fe-H alloys, and claimed that H is not a vital core light element from a geophysical sense: solid Fe-H alloys could not reproduce the shear wave velocity of the inner core. More recent studies (e.g., He et al., 2022;Wang et al., 2021) reported that H in the hcp Fe lattice at inner core conditions is extremely diffusive, which is referred to as a "superionic state." To study the contents and influence of H on the properties of Fe in the Earth's core, the structure of iron hydride should be investigated first. The crystal structures and stabilities of Fm3m-FeH and Pm3m-FeH 3 (Bazhanova et al., 2012); the electronic properties, and superconductivity of I4/mmm-FeH 5 and C2/c-FeH 6 (Kvashnin et al., 2018); and the metallic nature of C2/c-FeH 6  have all been reported. A topological analysis and phase diagrams of Pm3m-FeH 3 , I4/mmm-FeH 5 , and Cmmm-FeH 6 at high P-T conditions have also been presented (Sagatova et al., 2020). However, these previous Fe-H binary structure prediction studies failed to fully consider the thermal stabilities of the structures at high pressures and temperatures. In particular, the potential existence of "superionic" H means that the extrapolation of low P-T data could lead to inaccurate stabilities and properties. The H budget of the inner core needs to be re-estimated at appropriate P-T conditions.
In this work, in order to further explore the possibility, structure, mobility, and concentration of H in the Earth's inner core, especially under high temperatures, we have employed evolutionary crystal structure prediction methods and density functional theory (DFT) (Kohn & Sham, 1965) calculations to examine the structural models of Fe-H binary at core pressure and temperature conditions. The influence of temperature on the stabilities of Fe-H binary has been simulated by using the PHONOPY package (Togo & Tanaka, 2015) within the quasi-harmonic approximation (QHA) framework. Molecular dynamics calculations are also performed to detect the state and mobility of H under core temperature conditions. The ionic conductivity of Fe-H alloy, as well as the H concentration in the Earth's inner core was determined, and its implications on the composition and evolution of the Earth's core are discussed.

Ab Initio Calculations
Fe-H binary structures were calculated by DFT (Kohn & Sham, 1965) within the generalized gradient approximation (GGA) (Perdew et al., 1996) and the projector-augmented wave (PAW) method (Blöchl, 1994). The calculations were performed in the Vienna Ab initio Simulation Package (VASP) code (Kresse & Joubert, 1999). We found an energy cutoff of 600 eV and a Monkhorst-Pack scheme (Monkhorst, 1976) with a k point grid of 2π × 0.05 Å −1 gave excellent stress tensors and structural energy convergence for the Fe-H systems. The variable-composition structure search has been conducted by the crystal structure prediction code USPEX (Lyakhov et al., 2013;Oganov et al., 2011), where the selected parameters have been demonstrated to have high reliability in many successful applications from single elemental solids to many-body compounds (Chen et al., 2015). The generated structures were then fully relaxed under the core pressure conditions.

QHA Calculations
Theoretical phonon dispersion and the Gibbs free energy of Fe-H systems were obtained by using the PHONOPY package (Togo & Tanaka, 2015) with quasi-harmonic approximation. The temperature and volume dependence of the Helmholtz free energy F(V, T) is given by: where U(V ) is internal energy, F el (V, T ) is electronic free energy, and the last term F vib (V, T ) is the vibrational free energy that can be calculated as: where N is the number of atoms, q is the wave vector, N q is the number of wave vectors, ω i (q) is the phonon frequency, k B is the Boltzmann constant, and ℏ is the reduced Plank constant. The contribution to the free energy from electronic excitations at different electronic temperatures F el (V, T) was included by using a Fermi smearing function when calculating the ground-state energy in the VASP code. The Gibbs free energy can be obtained from: Note that QHA neglects any anharmonic effects. However, a previous study has calculated the anharmonic contribution to free energy in an hcp-Fe system and found that it was less than 60 meV/atom at 6000 K (Alfè et al., 2001). In the Fe-H binary system, with the incorporation of H, the anharmonic contribution to the free energy may be higher, which will be estimated in the text.

MD Calculations
Molecular dynamics (MD) simulations were conducted to explore the state and mobility of hydrogen in the Fe lattice under core temperature conditions. MD calculations were run at the gamma point with a cutoff of 450 eV, the electron step energy converges to 10 −4 eV, and all runs were calculated with a cubic 2 × 2 × 2 unit cell of 64 atoms. By using the Nosé thermostat (Nosé, 1984), MD trajectories were implemented in the NVT ensemble for 5 picoseconds (ps), which was enough to reach reasonable convergence after 10 ps NPT simulations. We retrieved the mean-square displacements (MSDs) and Radial Distribution Function (RDF) from each MD simulation to ensure that the systems remained solid throughout.

Fe-H Binary at Ultrahigh Pressures
We examined the lowest-enthalpy crystal structures of the Fe-H binary at 300 GPa. The results are consistent with previous studies (Bazhanova et al., 2012;Kato et al., 2020;Kvashnin et al., 2018;Pépin et al., 2014;Zhang et al., 2018), where the phase of Fm3m-FeH was found to be the most stable at core pressures ( Figure 1). In this work, we have predicted extra Fe-H phases that would be stable when the H fugacity is high, such as Pm3m-FeH 3 , P 4 /mmm-Fe 3 H 5 , C2/m-Fe 3 H 11 , and C2/c-FeH 6 . Our static calculations also produced the stable pure phases of Fe (hcp-Fe) and H (Cmca-H) at high pressures, which agree with the available experimental (Tateno et al., 2010) and theoretical evidence (Pickard et al., 2012), respectively. The c/a ratio of the produced hcp-Fe phase is 1.58, which is consistent with the previous research by Steinle-Neumann et al. (2001) and Vočadlo (2007). For each Fe-H phase, we identified the most favorable crystal structure and calculated its formation enthalpy (∆H) from the elements at finite pressures. The ∆H of the most stable phases is defined as: These enthalpies, normalized per atom, are given in Figure 1. The stable phase is located on the convex hull with regard to a formation from basic elements. Previous research by Bazhanova et al. (2012) and Sagatova et al. (2020) found metastable structures with compositions Fe x H y (x > y): Fe 4 H, Fe 3 H, and Fe 2 H. We confirmed that these structures were metastable at inner core pressures using the enthalpies of formation calculated here. As shown in Figure 1 prediction. The crystal structures and cell parameters of these stable components are shown in Figure 2 and Table 1. No imaginary phonon frequencies were found in the phonon dispersion curves of these stable structures ( Figure  S1 in Supporting Information S1), indicating their dynamic stability. In addition, these iron hydrides remain stable relative to pure Fe and pure H phases in the whole pressure range of the Earth's core ( Figure S2 in Supporting Information S1).
It is important to consider the underlying Fe sub-lattice in these structures as the Fe lattice will be a major component in their elasticity. If light elements induced a considerable rearrangement of Fe, this should be seismically visible. By doing a sub-lattice analysis, we found that the Fe in FeH and FeH 3 remains in a perfect fcc arrangement, which indicates that the main structure of the H-bearing phase in the Earth's inner core should be fcc. In Fe 3 H 5 and Fe 3 H 11 , the Fe remains in a bcc lattice but with considerable distortion. The distortion of Fe 3 H 5 sees one axis increasing by 20% relative to the other axes. In Fe 3 H 11 , we see a slight tendency to form a monoclinic structure where the A-axis increases by 5% and the C axis increases by 6% relative to the B axis, and the angle α decreases to 80°. In FeH 6 , very major distortions are seen with α decreasing to 60° and a structure that is somewhere between bcc and hcp ordering. In the presence of large amounts of FeH 6 significant seismic distortions in iron should therefore be seen.
To examine the possibility that hydrogen makes very large changes to the electronic structure of the Fe, we calculated the Bader charge of Fe atoms in both hcp-Fe and Fe x H y phases. As shown in Figure 3, H makes only small changes to the oxidation state of Fe. In pure Fe metal, this oxidation state is 0 ( Figure S3 in Supporting Information S1). For the hydrides, it ranges from 0.20 to 0.32 with no obvious pattern between H concentrations, which is consistent with the previous work . These results indicate that the system is much closer electronically to behaving like Fe metal than an iron hydride (which would have formal charges of 1) even as the stoichiometry of the system approaches an iron hydride.

Thermal Stabilities of Fe-H Binary
The previous section considered the stabilities of Fe-H binaries at static conditions. Since the core is extremely hot, the thermal stabilities of the produced Fe-H binaries were determined at high temperatures. As shown in Figure 4, the Pm3m-FeH 3 structure exhibits less stability at high temperatures, which indicates that it becomes metastable as the temperature increases.   Table S1 in Supporting Information S1. We also considered the anharmonic effect of H. The contribution of anharmonic terms to the free energy (F anharm ) can be approximated by: where the volume dependence of α(V) is approximated by: with α 1 = 1.8 × 10 −9 eVK −2 and α 2 = 4.8 × 10 −10 eV Å −3 K −2 per atom (Alfè et al., 2001). We calculated the total F anharm of hcp-Fe at 360 GPa and 5000 K to be ∼33 meV/Fe atom, which is increased by up to 16 meV/H atom with the additional presence of H in an Fe-H binary. This shows that the calculated F anharm of Fe x H y is of the same order of magnitude as that of hcp-Fe. Therefore, the contribution of anharmonic terms to Fe-H energies is likely small enough that they do not change the relative thermal stabilities of Fe-H binaries. Considering that the H content of the Earth's core is likely to be small (less than 1 wt%), Fe-H binary phases with high H concentrations (e.g., Fe 3 H 5 , FeH 3 , Fe 3 H 11 , and FeH 6 ) are unlikely to form and will not be investigated further in this study. We will focus on the fcc-FeH phase and its properties and effects under the Earth's inner core conditions.

Superionic H in the Solid Inner Core
Recent studies have proposed that superionic H could soften the shear modulus of hcp Fe and this would allow hcp Fe with H to match geophysical observations of the Earth's core (He et al., 2022;Wang et al., 2021). To examine the possible presence and effect of superionic H in solid fcc Fe we have performed MD calculations on fcc phase FeH. The hydrogen atoms were observed to move throughout the supercell in these calculations. This was seen both through a constant increase in MSDs with simulation time ( Figure 5) and visual inspection of the MD runs where H was observed to move between different sites of the crystal over time.
The diffusion coefficients for H were calculated using the MSD of the ionic positions at the superionic state. The mean square displacement (MSD) is averaged over specific ions: where ⃖⃖ ⃗( ) is the displacement of the ith ion at time t, and N is the total number of specific ions in the system. The diffusion coefficient is defined as: where d is the dimension of the lattice on which ion hopping takes place. The value of D obtained at various temperatures was then fit with an Arrhenius equation: where A is a pre-exponential factor, ∆ is the activation enthalpy, k B is the Boltzmann constant, and T is the temperature. The ionic conductivity was calculated using the diffusion coefficients and the Nernst-Einstein equation: = fDc 2 (10)   Figure S4 in Supporting Information S1.
in which σ is the ionic conductivity, f is a numerical factor approximately equal to unity, D is the diffusion coefficient, c is the concentration of light elements defects, q is the electrical charge, k B is the Boltzmann constant, and T is the temperature.
We determined the H diffusion coefficient in fcc FeH to be ∼10 −8 -10 −9 m 2 /s (above a critical temperature of ∼3500 K) ( Table 2), demonstrating that H has significant mobility as fluids in the Earth's inner core, which is consistent with the previous study by Wang et al. (2021). The solid-superionic phase transition temperature for hydrogen in our fcc-FeH simulations is thus estimated to be ∼3500 K which is around ∼500 K higher than that of the previous study on hcp Fe-H system (He et al., 2022). This mobility or diffusion coefficient increases dramatically with temperature. Note that the diffusion coefficient may be subject to the size effect (Li et al., 2022), in which an analytical correction based on hydrodynamic theory (Yeh & Hummer, 2004) increases the diffusion coefficient of H by less than 10% at temperatures from 4500 to 7000 K. Therefore, the diffusion coefficient of H in this study might be slightly (<10%) underestimated.
The melting point of pure Fe at the inner-core boundary (ICB) pressure condition (∼330 GPa) has previously been determined to be ∼6200 K (Alfè et al., 2002;Sun et al., 2018) but the presence of hydrogen could reduce this melting point substantially (Hirose et al., 2019;Sakamaki et al., 2009). Our simulations of high temperatures up to 5000 K are likely beyond the melting points of FeH x (Kato et al., 2020) though this depends somewhat upon the H concentration. Our simulations are in a solid phase, however, as shown in Figure 5 where we observe that the MSDs of Fe are small and do not increase with simulation time even at 5000 K (and 360 GPa). This indicates that despite the liquid-like mobility nature of hydrogen the iron in fcc FeH remains as a solid phase. Due to the uncertainty of the FeH x melting points at the Earth's inner core, our simulations might be of a metastable solid state. Regardless of our main conclusions about the H diffusion coefficient remain since a similar H diffusion coefficient in liquid and solid hcp Fe phases has been proposed in previous studies (Umemoto & Hirose, 2015;Wang et al., 2021).
One obvious change that the presence of highly mobile H would induce in iron is an increase in conductivity. If this conductivity increase is sufficiently large, the concentration of hydrogen could be detected through it and accurately knowing the conductivity of the core is important for understanding the energy budget of the geodynamo and the heat flux of the Earth's core (De Koker et al., 2012;Xu et al., 2017). Previous calculations of this conductivity on hcp Fe-H alloys showed the contribution of superionic hydrogen to the overall conductivity (ionic conductivity) to be 2-3 orders of magnitude lower than the conductivity derived from the iron (electronic conductivity) and thus largely negligible (He et al., 2022). We determined an ionic conductivity in fcc FeH of the order of 10 2 -10 3 S m −1 from 3000 to 5000 K, which are of a similar magnitude to those determined on the hcp Fe-H system (He et al., 2022). Thus we conclude for both hcp and bcc hcp the presence of highly mobile hydrogen does not substantially alter the conductivity of the core. On the other hand, highly mobile hydrogen could affect the geodynamo in other ways by promoting convection and migration of core materials through the high mobility of light elements but this needs to be further investigated in the future.

Hydrogen Concentration in the Inner Core
One or more light element-bearing components are needed to explain the density deficit of the Earth's core (Anderson & Isaak, 2002). In addition to other light elements (Si, O, S, C, N, etc.), H has long been considered a potentially major light component in the core (Poirier, 1994). The melting point of Fe alloys should be a function of the concentration of light elements (He et al., 2022), and so the value of the H concentration in the Earth's inner core is not only important for understanding the core compositions but also the process of inner core crystallization. is also shown in Figure S5 in Supporting Information S1.

Table 2 The Calculated Hydrogen Diffusion Coefficients (D) in Fcc-FeH and
Hcp-FeH x at 360 GPa and Different Temperatures (3000,3500,4000,5000,and 6000 K) The H concentration needed to reproduce the density of the Earth's inner core was estimated by matching the observed density (PREM) (Dziewonski & Anderson, 1981) to the density of a mixture of hcp-Fe and the stable FeH iron hydride at the relevant pressures and temperatures. Figure 6 shows the compression curve of Fe and FeH phases as a function of pressure. Our calculated equation of state (EoS) for fcc-FeH is consistent with the experimental values given by Tagawa et al. (2022). The compositional effect of H on the density was treated as linear with Equation 11. Considering the Fe and Fe-H mixture, the molar concentration of H was determined with: where IC is inner core density from the PREM model, Fe and Fe−H are the computed densities of hcp-Fe and Fe-H iron hydride at relevant temperatures, indicates the molar fraction of H in Fe-H iron hydride, and is the molar concentration of H. Figure 7 shows the predicted H concentrations calculated by Equation 11 as the function of depth at 4000 and 5000 K. We found that temperature has only a small effect on H concentration in the inner core (∼1.6%, from 4000 to 5000 K) while the effect of pressure is more visible (∼3.3%, from 330 to 360 GPa). The concentration of H in the core decreases as the temperature and the pressure increase. The H concentration required to match the density deficit of the inner core is ∼0.34 mol (∼0.92 wt%). This result is very close to the recent estimate value by Tagawa et al. (2022) who proposed that the maximum H content of the inner core is about 0.8-0.9 wt% based on their laser-heated diamond-anvil cell (LHDAC) experiments.
The density jump at the ICB (Koper & Pyle, 2004;Tkalčić et al., 2009) can be used to further constrain the H budget in the inner core (Li et al., 2022). We used the effect of H on density in the solid determined from this study and the effect on the liquid from Li et al. (2022), and then directly compared the calculated density difference between solid and liquid states for pure Fe and Fe-H alloy under core conditions. Figure 8 shows the density jump between the solid and liquid as a function of X H /(X H + X Fe ). Using this method, the average H content in the inner core is further constrained to ∼0.21 wt%, which is equivalent to ∼1 ocean of water in the inner core. This result is consistent with many previous estimates (Caracas, 2015;Shibazaki et al., 2012;Thompson et al., 2018;Wang et al., 2021). Therefore, the upper limit of the H budget in the inner core can be significantly reduced from 0.92 wt% to 0.21 wt%.
It should be pointed out that the concentration of the light element in the Earth's core can be estimated more precisely by considering both the core's density and the elastic wave velocities, as well as the co-existence of other possible light elements. The calculations above are based on the assumption of H being the only light element in the core. Other light elements are most likely to be present in the core and will also reduce the Fe alloy's density. Hirose et al. (2019) demonstrated that H limits carbon solution into liquid Fe and that the melting temperature of Fe-H alloy was predicted to be lower than other binary Fe alloy systems at the core-mantle boundary (CMB). Thus the behavior of H in the core would affect other light elements partitioning into the Earth's interior during the later stages of its evolution. Kato et al. (2020) determined the crystal structure of FeH x (x ∼ 1) at high pressure and temperature in X-ray diffraction (XRD) measurements. They observed that FeH x forms fcc phase at the pressure of 43-137 GPa and temperature of ∼1000-2000 K. And the pressure would bend the dhcp-fcc  transition boundary back to low temperatures. In this work, we have predicted the stable Fe x H y phases and their thermal stabilities at the conditions of the Earth's core. Our results show that fcc phase FeH is stable in the inner core, which supports the result of Kato et al. (2020). The phase boundary of dhcp-fcc in the Fe-H system is still controversial, however. Isaev et al. (2007) computed the static free energy calculations and proposed that the dhcp-hcp transition occurs at 37 GPa, and then a second phase transition hcp-fcc happens at 83 GPa. Further compressional experiments suggested that the dhcp phase can be stable up to 136 GPa at ambient temperature (Pépin et al., 2014). Alternatively, a phase transition dhcp-fcc induced by temperature was observed up to 20 GPa in a multi-anvil apparatus (Ikuta et al., 2019;Sakamaki et al., 2009). The fcc phase appears to be increasingly stable over the dhcp phase with increasing pressure and temperature. Our calculations on enthalpy under static conditions show that the fcc phase becomes more stable at pressures above ∼42 GPa ( Figure S6 in Supporting Information S1). Thus it is unlikely that the dhcp phase is stable at inner core pressures and temperatures, both of which promote the stability of the fcc phase, and so it has not been explicitly examined in this study though we cannot rule out its existence at extreme pressures, temperatures and with the addition of hydrogen.

Discussion
A recent study by Li et al. (2020) concluded that the Earth's core is a potentially large water reservoir by calculating the H and water partitioning between liquid Fe and silicate melts at CMB conditions. Moreover, previous studies (Tagawa et al., 2016;Thompson et al., 2018;Umemoto & Hirose, 2015) have proposed that the water storage capacity of the outer core could be dozens of times that of ocean water. These studies all support the hypothesis that H is a vital core light element. However, with the incorporation of H, the seismic wave velocities of Fe-H binary remain too high to match the geophysical observations of the inner core even at core temperature conditions (He et al., 2022). A ternary Fe-X-H (X = C, O, Si, S, etc.) or a larger combination of elements is needed to solve these dilemmas. For instance, Wang et al. (2021) reported that superionic H in a Fe-Si-H composition could simultaneously satisfy the density, sound velocities, and Poisson's ratio of the Earth's inner core, supporting the conclusion that H can be a significant component of the Earth's core. Since Fe-Si and Fe-S have similar elastic properties, Fe-S-H and Fe-Si-S-H ideal mixing under inner-core conditions could also explain the seismic properties of the inner core Wang et al., 2021).
The question would then arise if the Earth's core is a potentially huge water reservoir (several times to dozens of times as much as oceans of water), when and where did the H or water come from, and what is the behavior or effect of H in the Earth's core? The early-Earth ingassing model suggests that nebular H was captured during Earth's accretion (Olson & Sharp, 2018) and a combination of ingassing of nebular H and chondritic water has been further proposed as an explanation of the low D/H (deuterium/protium) value of the Earth (Wu et al., 2018). Our estimated H concentration of the Earth's inner core also supports the early delivery of H or water from nebular materials and/or a chondritic source before or during core-mantle differentiation. A recent study reveals that H limits the solubility of C in liquid Fe, which suggests that H-rich/C-poor metals were likely added to the protocore in the late stages of core formation (Hirose et al., 2019). The melting temperature of Fe decreases remarkably with the incorporation of H, which would lead to a late crystallization of the inner core and indicates a young inner core age (<1 Gyr) (Gomi et al., 2018).
As for other light elements (e.g., C, O, Si, S) in the Earth's inner core, since most C and O preferentially partitions into the liquid outer core (Alfè et al., 2002;Li et al., 2019), the role of Si and S should be critical in determining the formation and properties of the inner core. The difference in silicon isotopic compositions between the bulk silicate Earth and chondrites could be explained by the partitioning of ∼5 wt% Si into the Earth's core (Georg et al., 2007). Based on the S partition coefficient (∼0.75) between the solid inner and liquid outer core (Alfè et al., 2002;Zhang et al., 2020), there could be ∼0.7-1.5 wt% S in the inner core according to the level of S abundance in the Earth's mantle and the S partitioning between metallic and silicate melts during core formation (Mcdonough & Sun, 1995;Rose-Weston et al., 2009). In addition, H is of critical importance in replicating the properties of Earth's inner core and may have been the first Fe-dissolved light element during the evolution of the early Earth (Iizuka-Oku et al., 2017). Fe-H-X ternary or H-bearing multicomponent systems were likely to form in the core in its subsequent evolutions. The behavior of such components, and how H affects the subsequent dissolution of other light elements, is thus critical for understanding the composition and evolution of the Earth's core.

Conclusions
Our study suggests that the Fe-H binary adopts numerous possible structures under core-like conditions, while the fcc structure is concluded to be a strong candidate for the H-bearing phase in the Earth's inner core. The high mobility of H in the solid Fe lattice at high temperatures indicates that H is transferred to a superionic state, where the H superionic state transfer temperature in Fe fcc lattice is ∼500 K higher than that in hcp Fe system. About 0.92 wt% of H is required to match the density deficit of the inner core, and this H concentration value is further constrained to ∼0.21 wt% by considering the density jump at ICB. H is a key light element for reducing the density and elastic modulus of Fe, but the wave velocities of the Fe-H binary still remain too high to account for the seismological observations of the inner core. Other light elements are, therefore, also required to match all the geophysical models.

Data Availability Statement
Supporting Information can be found in Supporting Information S1. Data sources for Figure 2 are given in Table 1 and Table S1 in Supporting Information S1. Other data presented in this work are available in the repository (https://doi.org/10.6084/m9.figshare.20765965.v1).