Fe-H binaries at ultrahigh pressures
We examined the lowest-enthalpy crystal structures of the Fe-H binaries at 300 GPa and then optimized them with GGA calculations. The results are consistent with previous studies19-21,27,28, where the phase of was found to be the most stable at core pressures. In this work, we have predicted extra Fe-H phases that would be stable when the H fugacity is high, such as P4/mmm-Fe3H5, C2/m-Fe3H11, and C2/c-FeH6. Our calculations also produced the stable pure phases of Fe(hcp-Fe) and H(Cmca-H), which agree with the available experimental29 and theoretical evidence30, respectively. The c/a ratio of the produced hcp-Fe phase is 1.58 at 0 K and 360 GPa, which is consistent with the previous research by Steinle-Neumann et al.31 and Vočadlo32. For each Fe-H compound, 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.19 and Sagatova et al.22 found metastable phases with compositions FexHy (x>y): Fe4H, Fe3H, and Fe2H. We confirmed that these phases were metastable at inner core pressures using the enthalpies of formation calculated here. As shown in Figure 1, for FexHy with x ≤ y, FeH, FeH3, Fe3H5, Fe3H11, and FeH6 were found to be stable compounds at core pressure while the phases of FeH7, FeH8, FeH9, and FeH10 are metastable as they fall outside of the convex hull21,22. C2/m-Fe3H11 is a new Fe-H phase found in our structure prediction. The crystal structures and cell parameters of these stable components are shown in Table S1 and Figure S1. The phonon dispersion curves of these stable phases are shown in Figure S2, and no imaginary phonon frequencies are found for any of these phases, which indicates their dynamical stability. In addition, these iron hydride compounds remain stable relative to pure Fe and pure H phases in the whole pressure range of the Earth's core (see Figure S3).
It is important to consider the underlying Fe sublattice 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 sublattice analysis, we found that the Fe in FeH and FeH3 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 Fe3H5 and Fe3H11, the Fe remains in a bcc lattice but with considerable distortion. The distortion of Fe3H5 sees one axis increasing by 20% relative to the other axes. In Fe3H11, 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 degrees. In FeH6, very major distortions are seen with α decreasing to 60 degrees and a structure that is somewhere between bcc and hcp ordering.
To discard 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 FexHy phases. As shown in Figure S5, H makes only small changes to the oxidation state of Fe. In pure Fe metal, this oxidation state is 0. For the hydrides, it ranges from 0.20-0.32 with no obvious pattern between H concentrations. 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 and phase distribution of Fe-H binaries
The previous section considered the stabilities of Fe-H binaries at static conditions, but the core is very hot (5000-6000 K). The thermal stabilities of the produced Fe-H binaries were determined at high pressures. The results are presented in Figure 2 which shows the Pm3m-FeH3 phase becomes less stable at high temperature, which indicates it should be metastable with the temperature increasing, while other Fe-H binary phases are still stable at high temperatures.
We also considered the anharmonic effect of H. The contribution of anharmonic terms to the free energy (Fanharm) 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 Å-3K-2 per atom33. We calculated the total Fanharm of hcp-Fe at 360 GPa and 5000 K is ~33 meV/Fe atom, which is increased by up to 16 meV/H atom due to the presence of H in the Fe-H binaries, which shows that the calculated Fanharm of FexHy is in the same order of magnitude as that of hcp-Fe. Therefore, the contribution of anharmonic terms is small enough that they do not change the relative thermal stabilities of Fe-H binaries.
The small energy differences of Fe-H binaries (similar to kBT) and the high mobility of H atoms (see below) means that no one Fe-H structure will be dominant at core temperatures and a mixture of different Fe-H structures will be present at core conditions. From the formation energy of each structure, we calculated their favorability using the Gibbs free energy distribution at a specific pressure and temperature:
where i are the different structures, and Ei is the formation energy of that structure per atom. The probability of any specific arrangement is then:
Calculated results are listed in Table 1 where it is shown that while remains the most stable Fe-H binary phase as previously predicted all of the predicted phases have a reasonable probability of occurring and that Fe-H in the core should be a mixture of all of these phases. Changing the temperature across a range of core temperatures (from 5000 to 6000 K) has little effect on the probability of the phases occurring and thus the distribution of H in the core should be insensitive to fluctuations of temperature in the core.
In conclusion, we find the is the most stable Fe-H binary phase and this does not vary with temperature but that the introduction of core temperatures allows a variety of Fe-H binary phases to form and that no individual binary Fe-H phase controls the properties of the core. Thus the properties of all stable Fe-H phases must be considered when determining the fate of H in the core.
Hydrogen concentration in the inner core
One or more light element-bearing components are needed to explain the density deficit of the Earth's core34. In addition to other light elements (Si, O, S, C, N, etc.), H has long been considered a potentially major light component in the core35. The H concentration of the core can be estimated precisely by comparing predicted seismic properties with real seismic observations. The value of H concentration in the Earth's inner core is important for its mineralogy but also due to the role of H during inner core crystallization. Furthermore, this value also can be used to constraint the partitioning of H between the inner and outer core.
Based on the calculated density of Fe-H binaries at relevant temperatures (5000-6000 K) (see Table S2), we estimated the density of H-bearing compositions in the core. To do this we determined the theoretical density differences between hcp-Fe and stable FexHy phases (FeH, FeH3, Fe3H5, Fe3H11, and FeH6) and then treated the compositional effect of H on the density as linear (Fe and FexHy are treated as an ideal solid solution). Table S3 shows the parameters of the Birch Murnaghan third-order equation of state for hcp-Fe and the stable FexHy phases. The density and Bulk modulus of hcp-Fe and Fm m-FeH phases are consistent with previous work19. The H concentration can then be estimated with:
where is inner core density from the PREM model, and are the computed T-Kelvin densities of Fe and Fe-H system, indicates the molar fraction of H in Fe-H system, and is the molar concentration of H.
Figure 3 shows the corresponding H concentration calculated by Equation 6. The matching concentrations vary between 0.11 and 0.34 mol (0.22-0.91 wt.%) if we use the individually different stable FexHy phases: when H is in the FeH phase, the concentration of H is about 0.34 mol (0.91 wt.%), the Fe3H5 phase is about 0.25 mol (0.59 wt.%), the FeH3 phase is about 0.19 mol (0.42 wt.%), the Fe3H11 phase is about 0.18 mol (0.39 wt.%), and the FeH6 phase is about 0.11 mol (0.22 wt.%). We conclude that temperatures have only a small effect on H concentration in the inner core (~ -0.3%, from 5000 K to 6000 K) while the effect of pressure is more visible (~ -2%, from 330 GPa to 360 GPa). Note that the concentration of H in the core decreases as the temperature and pressure increase except for FeH3, wherein the concentration of H increases slightly as the temperature increases (see Figure S6). This anomalous behavior is further evidence that perhaps FeH3 is metastable under inner-core conditions as also determined from the energy and shown in Figure 2.
The above results for H concentration were for single Fe-H phases but as shown above, H in the core should have a variety of structures. Using the phase distribution in Table 1, we can calculate the H concentration for a real phase distribution in the core. Using an ideal mixture of 30% of the FeH phase, 22% of the Fe3H5 phase, 11% of the FeH3 phase, 22% of the Fe3H11 phase, and 15% of the FeH6 phase, we calculated the total concentration of H that is possible in the inner core based solely on density is about 0.57 wt.%, which is equivalent to 3 times as much as ocean water. This result is higher than the recent estimate of 0.2-0.3 wt.% by Thompson et al.16, but is compatible with the 0.4-0.5 wt.% proposed by Bazhanova et al.19 and Tagawa et al.17, as well as in reasonable agreement with the experimental values of 0.6 wt.% (a ratio of H to Fe in the inner core, H/Fe = 0.34) observed by Okuchi11.
These calculations correspond to the concentration needed if H is the only light element in the core. Other light elements are likely to be present in the core and will also decrease the density. Therefore, we estimate the H concentration to be ~0.57 wt% as its upper limit. And how does so much H in the inner core affect its temperature? Hirose et al.36 have demonstrated that H limits carbon solution into liquid Fe and the melting temperature of Fe-H alloy was predicted to be lower than other binary Fe alloy systems at the CMB. Assume that the melting temperature of the Fe-H system has a linear relationship with the H content, based on the extrapolation of experimental data fitted by Simon-Glatzel equation37 and the H content (0.57 wt.%) in our simulation models, the melting temperatures of H-bearing components at the inner core pressure conditions are estimated, which is ~5300 K (see gray regions in Figure 4), indicating that H would induce a relatively low inner core temperature.
The seismic properties of solid iron hydrides
Light elements incorporated into Fe under extreme conditions can significantly alter its structure as shown by our structure prediction results where a conversion from hcp to fcc was observed and as suggested previously by other works38,39. A previous study by Caracas18 has investigated the effect of H on the seismic properties of solid Fe but they only considered the hcp structure of Fe-H phases and the consideration of temperature is not complete. Multiple Fe-H structures are important in the core and high-temperature effects on elasticity could be important. Therefore, the elastic properties of the different Fe-H alloy structures need to be investigated and should be distinguished from that of the pure hcp-Fe structure.
The static elastic constants of hcp-Fe, fcc FeH and FeH3 are shown in Table S4, and the derived seismic properties are listed in Table 2. Our hcp-Fe phase results are well consistent with previous work40,41. We find that H reduces the density of hcp-Fe significantly (each 1 wt.% H reduces the density by 7.6%), while H increase both the compressional and shear wave velocities (each 1 wt.% H increases the Vp and Vs by ~2.3% and ~1.5%, respectively). The value of Poisson's ratio (0.32-0.33) hardly changes with the different H concentrations or structures, and the Debye sound velocity is only a small increase (~2%) for each 1 wt.% of H.
Our results have a similar general trend to that observed by Caracas18 but they obtained a larger compressional wave velocities increase (~10% increase of Vp with 1 wt.% H). Alternatively, we predict that H reduces the elastic modulus of solid hcp-Fe, which is contrary to Caracas's conclusion but the wave velocities increase in both cases due to the reduction in density which is the largest control. The difference between our results and that of Caracas is likely due to the fact that different Fe-H structures exhibit different seismic properties. Therefore, determining the actual effect of H on velocities requires determining its distribution and structure. As explained above H can introduce a phase change to Fe but the predicted elastic effect of this is very small and the effect of H should be much larger than the effect of an hcp to fcc phase change in the Fe sublattice18,32.
These results were calculated statically, and the influence of temperature on the elastic and seismic properties also needs to be evaluated. We investigated the temperature effect on the most stable phase FeH by calculating elasticity with Molecular Dynamics (MD) simulations. The seismic properties of hcp-Fe and fcc-FeH at high-temperature (5000 K) are shown in Table 3. Compared with the static results, there is a reduction of Vp and Vs in FeH by ~8% and ~13%, respectively. However, the sound velocities of FeH at high temperatures still remain too high to account for the geophysical observations. Thus, Fe-H binaries alone are unlikely to provide a seismically suitable structure for the inner core.
Another interesting phenomenon is observed in our MD calculations. At the temperature of 360 GPa and 5000 K, we observe that while Fe-H is still solid (as observed in its RDF and RMSD plots, see Figure S7) the H atom shows a continual drift in time which is not seen at 3500 K. This is also observed in a manual inspection of the results as H atoms move between different Fe sites over time. This suggests that above ~3500 K the Fe remains in a solid sublattice but that the H sublattice converts to a reported superionic state42 where the H atoms can freely move. By determining the slope of the RMSD we determined the H diffusion coefficient to be 10-8-10-9 m2/s above the temperature of ~4000 K, which demonstrates that the H has significant mobility in Fe-H system. The superionic H may contribute to the conductivity of the core through the movement of charge carriers, but such an analysis is outside the scope of this work.
A recent study by Li et al.43 concluded that the Earth's core is a potentially large water reservoir by calculating the H or water partitioning between liquid Fe and silicate melts at core-mantle boundary (CMB) conditions. Moreover, previous studies15-17 have proposed that the water storage capacity of the outer core would be dozens of times that of ocean water. These studies all support that H is a vital core light element. However, with the incorporation of H and thus the wave velocities of Fe-H binaries remain too high to match the geophysical observations of the core even at core temperature conditions. A ternary combination Fe-X-H (X = O, S, Si, C, etc.) is needed to solve these dilemmas. More recently, Wang et al.44 reported that the superionic H in Fe-Si-H composition could simultaneously satisfy the density, sound velocities, and Poisson's ratio of the Earth's inner core, supporting that H is a significant component of the Earth's core. During the evolution of the early Earth, the first light element that Fe dissolved was H14; therefore, Fe-H-X ternaries were likely the first to form and the behavior of such ternaries, and how H affects the subsequent dissolution of light elements36, is critical in understanding the core’s composition and evolution.