New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Dynamical stability of FeH in the Earth's mantle and core regions

Edited by Hokwang Mao, Carnegie Institution of Washington, Washington, DC, and accepted March 14, 2007 (received for review November 1, 2006)
Abstract
The core extends from the depth of 2,900 km to the center of the Earth and is composed mainly of an ironrich alloy with nickel, with 10% of the mass comprised of lighter elements like hydrogen, but the exact composition is uncertain. We present a quantum mechanical firstprinciples study of the dynamical stability of FeH phases and their phonon densities of states at high pressure. Our freeenergy calculations reveal a phonondriven stabilization of dhcp FeH at low pressures, thus resolving the present contradiction between experimental observations and theoretical predictions. Calculations reveal a complex phase diagram for FeH under pressure with a dhcp → hcp → fcc sequence of structural transitions.
The solid inner core of our planet is believed to be composed of iron. The core density as determined from seismological data is, however, ≈2–5% lower than that of pure Fe at the corresponding conditions (≈300 GPa; ≈6,000 K). This finding gave rise to the suggestion that the inner core might also contain a certain amount of light elements, in particular, hydrogen. Xray diffraction experiments have shown that at room temperature and pressure >3.5 GPa Fe and H_{2} can react, forming nearly stoichiometric iron hydride (1–3). Moreover, FeH_{x} has been demonstrated to have an increased stability with respect to Fe, silicate, and water over increasing pressure–temperature conditions (4). Nuclear resonant experiments (5) have provided additional important information on the lowfrequency (LF) phonon dynamics of FeH_{x} up to ≈50 GPa. However, knowledge about the full phonon spectra of iron hydride at higher pressure is still lacking.
Because of a poor solubility of H in Fe at ambient conditions, the formation of iron hydride was considered to be unlikely until the first successful experimental synthesis of dhcp FeH was carried out (1, 6). The volume expansion was ≈17% compared with pure Fe. The dhcp structure of FeH was confirmed by neutron diffraction data (7) and also supported by Mössbauer experiments (8). In situ xray diffraction experiments showed that FeH remained in the dhcp structure up to at least 62 GPa (6), and none of the studies by Hirao et al. (3) and Saxena et al. (9) revealed any structural changes up to 75–80 GPa and 2,000 K.
The influence of hydrogen on the relative stability of iron phases and properties of iron hydride structures under high pressure have been studied theoretically (10, 11). In agreement with the results of neutron diffraction experiments (7) these studies show that hydrogen prefers to occupy octahedral interstitial positions in the closepacked Fe structures. For iron hydride, however, ab initio calculations (10, 11) predict the stability of the hcp structure rather than dhcp found in experiments (1, 3, 6–9). This discrepancy has to be resolved before one starts to seriously consider the probability of substantial hydrogen content within the core region of the Earth. To address this problem, we applied methods based on the first principles of quantum mechanics. Details about the performed calculations can be found in Methods.
Results and Discussion
The comparison of the total energies [E(V)] calculated for the dhcp, hcp, and fcc FeH structures (Fig. 1 a) show that hcp should be most stable up to ≈80 GPa (V = 10.22 Å ^{3} ) (throughout the article volumes and energies are given per formula unit FeH), where a hcp → fcc phase transition, promoted by a ferromagneticparamagnetic (FMPM) transition occurring in the fcc structure at ≈60 GPa, takes place. At ≈80 GPa a rather sharp FMPM transition also occurs in dhcp FeH. These results are in good agreement with previous ab initio calculations (10, 11).
Our freeenergy [F(V,T)] calculations at room temperature show that the mentioned discrepancy between the theoretically predicted (hcp) and experimentally observed (dhcp) FeH structures is resolved if the vibrational term [F _{vib}(V,T)] is accurately taken into account: F(V,T) = E(V) + F _{vib}(V,T). The inclusion of this term leads to a stabilization of the dchp phase with respect to hcp and fcc at low pressures (Fig. 1 b). In fact, this stabilization is already achieved if the zeropoint vibrations of hydrogen atoms are taken into account. According to our freeenergy calculations, there exist two phase transitions in the FeH system. At ≈37 GPa a dhcp → hcp phase transition takes place, and after the transition, the hcp phase remains stable up to ≈83 GPa, where a hcp → fcc transition occurs. Therefore, we find three pressure intervals where the dhcp (up to 37 GPa), hcp (37–83 GPa), and fcc (>83 GPa) FeH structures successively become stable (Fig. 2). Interestingly, based on experimental results, Hirao et al. (3) recently suggested that one can distinguish three pressure regions in the equations of state (EOS) of FeH: low (up to 30 GPa), intermediate (between 30 and 50 GPa), and high (>50 GPa). In Fig. 2 we compare the experimentally obtained EOS for iron hydride (2, 3) with those calculated for hcp, dhcp, and fcc FeH. We notice that the structures considered here are strictly stoichiometric and might be quite far from the composition of the phases obtained in real experiments. Deviations from stoichiometry can change the energy balance, shifting the phase transitions on the pressure and temperature scales.
Comparing the total and freeenergy results one might ask why in the low pressure limit the total energies and free energies indicate the existence of different phases, hcp and dhcp, respectively, but at high pressures the fcc phase is most stable according to both of them? To answer this question let us analyze the total energy and phonon contributions to the free energy for three volumes representing low (13.71 Å^{3}), intermediate (10.74 Å^{3}), and high (9.63 Å^{3}) pressure regimes (Fig. 3 Lower Inset). At low pressures (V = 13.71 Å^{3}) the phonon contributions to the free energy are found to be very similar for dhcp and fcc FeH [ΔF _{vib}(dhcp − fcc) ≈0.1 mRy at 300 K] but both are substantially lower than that for hcp FeH [ΔF _{vib}(hcp − dhcp) ≈2.2 mRy at 300 K]. Although the total energy difference between hcp and dhcp FeH speaks in favor of the former structure (Fig. 1 a), a large positive vibrational contribution makes the hcp phase less stable than dhcp in this pressure range. At intermediate pressures (37–83 GPa) the fcc phase has a lower vibrational energy [(ΔF _{vib}(dhcp − fcc) ≈1.22 mRy at 300 K] than hcp and dhcp (Fig. 3 Lower Inset). According to the total energies, however, the phase order in this pressure range should be E _{hcp} < E _{dhcp} < E _{fcc}. Because of an interplay between the electronic and phonon contributions to the free energy, hcp emerges as the most stable phase in the intermediate pressure range (37–83 GPa). In the high pressure limit (>83 GPa) the phonon energies for all three phases lie in a narrow interval (Fig. 3) that allows the electronic term to entirely determine the phase order and stability: E _{fcc} < E _{dhcp} < E _{hcp}. That is why we find fcc FeH to be most stable for these pressures also according to the free energies (Fig. 1 b). Thus, the stability of hcp (37–83 GPa) and fcc (>83 GPa) is determined mostly by the electronic term, whereas the stabilization of the dhcp phase at low pressures originates from a dominating phonon contribution to the free energy.
To further understand this behavior we have calculated an atomprojected partial phonon density of states, g _{i}(ω), and atomresolved phonon contributions to the free energy F _{vib} ^{i}(V,T) obtained from g _{i} (ω) where i is an atomic label and e _{nk} ^{i} is polarization vectors of atom i for a given mode and wave vector. In Fig. 3 we present the phonon density of states for the hydrogen and iron atoms for the three structures for V = 13.71 Å^{3}. The calculated spectra consist of two well separated (in frequency) parts: a LF part, which is mostly caused by the iron contribution with an admixture of hydrogen modes, and a highfrequency (HF) part, which is exclusively caused by the hydrogen optical modes. At low pressures the LF part is dominated by the iron contribution in all three structures. The application of pressure leads to a shift of the phonon spectra toward higher frequencies and an increase of hydrogen contribution to the LF modes.
Different from the fcc phonon spectrum, the HF parts of the spectra for hcp and dhcp have pronounced shoulders, extending to higher frequency in the case of the former phase (Fig. 3). The HF vibrational modes for hcp FeH bring a larger positive contribution to the free energy compared with the other structures that leads to its instability. Thus the decisive contribution for the phase stability comes from the optical H vibrations; therefore, the arrangement of the hydrogen sublattice is important.
The closepacked FeH structures have different atomic stacking sequences in the z direction (Fig. 2). In all of them a hydrogen atom has six hydrogen neighbors within the same atomic layer (n) perpendicular to the z axis. In addition, in fcc H also has six equally distant, but rather far situated, H neighbors in the (n − 1) and (n + 1) layers. In dhcp, H has three neighbors in the (n + 1) layer and one H atom in the (n − 1) layer at a short distance (≈c _{dhcp}/4) (Fig. 2). In hcp, H is squeezed between two hydrogens lying right above (n + 1 layer) and under (n − 1 layer) it at short distances (≈2c _{hcp}/4). Such a tight linear arrangement of hydrogen atoms in the hcp structure can lead to “unfavorable” HF optical vibrations. Indeed, detailed analysis of the calculated atomic vibrations shows that hydrogen vibrations along the z axis are responsible for the appearance of the large HF shoulder seen in the hcp spectrum that causes the instability of hcp FeH (Fig. 3).
To further compare our results with experimental data, we have calculated the pressure dependences of the Debye sound velocity (V _{D}) by integrating the calculated phonon spectra (12) (Fig. 4). The experimental Debye velocity have been extracted from the phonon densities of states obtained in nuclear resonant inelastic xray scattering experiments (5, 13). In the studied pressure interval the theoretically calculated velocities monotonically increase with pressure, except V _{D} for fcc FeH, which exhibits a sharp increase at ≈60 GPa where the ferromagneticparamagnetic transition takes place. The velocities for hcp FeH are higher than those for dhcp and fcc. The theoretical V _{D} values for dhcp FeH are higher than the experimental ones, which is not surprising considering the differences between the calculational setup and real experiments: the former deals with perfect monocrystalline FeH, whereas in the latter case we have polycrystals of unknown and most probably hydrogendeficient stoichiometry. Nevertheless, the calculated points follow the experimental ones rather well at least up to 20 GPa. At ≈37 GPa theory predicts a dhcp–hcp transition to occur; therefore, an abrupt increase of the Debye sound velocity in the vicinity of the transition pressure can be expected. Unfortunately, the available experimental points are rather few and scattered; therefore, more detailed experimental investigations are needed to verify our suggestions.
Our analysis of the temperature dependence of the FeH phase diagram showed that the pressure of the dhcp to hcp transition quite noticeably increased with temperature: 37 GPa (300 K), 45 GPa (700 K), and 50 GPa (1,000 K), whereas the pressure of the hcp to fcc transition changed insignificantly: 83 GPa (300 K) and 85 GPa (700 K). Therefore, the trend of increasing stability range for dhcp FeH is clearly indicated by our calculations.
Methods
The performed calculations are based on density functional theory (14) and density functional perturbation theory (15) as implemented in the Quantum ESPRESSO code (www.pwscf.org). The generalized gradient approximation in the PerdewBurkeErnzerhof parametrization (16) was chosen to treat the exchangecorrelation effects. The ultrasoft pseudopotentials in the form introduced by Rappe et al. (17) were used. The integration over the Brillouin zone was performed by using the special kpoints technique (18) with a broadening of 0.025 Ry according to the Marzari–Vanderbilt smearing scheme (19). We used 8 × 8 × 4, 16 × 16 × 16, and 12 × 12 × 8 kmeshes, for dhcp, fcc, and hcp FeH, respectively. Force constants were obtained by using eight dynamical matricies for all of the phases. The test calculations performed for fcc and hcp (dhcp) FeH using 16 and 21 dynamical matrices, respectively, have shown that such a dramatic increase in the number of qpoints led to changes in zero vibration energy <0.05 mRy. The phonon spectra with respect to the number of kpoints are within the accuracy of 0.01 THz. The electron wave functions were expanded in the plane wave (PW) basis set with a cutoff energy of 40 Ry, and PW with kinetic energy up to 600 Ry were used for describing the augmented charge. All of the calculations allowed for spin polarization. The ratio between the hexagonal aaxis and caxis (c/a) for the hcp and dhcp structures were optimized for each volume and found to be close to experiments (3) and only weakly dependent on pressure. For our phonon calculations we used (c/a) ^{hcp} = 1.63 and (c/a) ^{dhcp} = 3.26. We extended our calculations to nonzero temperatures by using the quasiharmonic approximation where the volume and temperature dependence of the free energy F(V,T) was calculated as a sum of the total energy, E(V), and phonon contribution, F _{vib} (V,T):F(V,T) = E(V) + F _{vib} (V,T). The last term is defined as: where k _{B} is the Boltzman constant, and ω _{nq} is the frequency of the phonon mode for wave vector q and volume V.
Acknowledgments
This work was supported by the Swedish Research Council (VR), FUTURA, European Science Foundation (EuroMinScI), the Swedish Foundation for Strategic Research (SSF), and the Swedish Energy Agency (STEM). The Swedish National Infrastructure for Computing provided computing resources.
Footnotes
 ^{¶}To whom correspondence should be addressed. Email: borje.johansson{at}fysik.uu.se

Author contributions: N.V.S., R.A., and B.J. designed research; E.I.I. performed research; E.I.I., N.V.S., R.A., Y.K.V., and B.J. analyzed data; and E.I.I., N.V.S., and R.A. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.
 Abbreviations:
 HF,
 high frequency;
 LF,
 low frequency.
 © 2007 by The National Academy of Sciences of the USA
References

↵
 Antonov VE ,
 Belash IT ,
 Degtyareva VF ,
 Ponyatovsky EG ,
 Shiryayev VI

↵
 Badding JV ,
 Mao HK ,
 Hemley RJ
 Syono Y ,
 Manghnani M
 ↵
 ↵
 ↵

↵
 Badding JV ,
 Hemley RJ ,
 Mao HK
 ↵

↵
 Schneider G ,
 Baier M ,
 Wordel R ,
 Wagner FE ,
 Antonov VE ,
 Ponyatovsky EG ,
 Kopilovskii Y ,
 Makarov E
 ↵
 ↵
 ↵

↵
 Ashcroft NW ,
 Mermin ND

↵
 Hu MY ,
 Sturhahn W ,
 Toekkner TS ,
 Mannheim PD ,
 Brown DE ,
 Zhao J ,
 Alp EE
 ↵
 ↵
 ↵
 ↵

↵
 Monkhorst HJ ,
 Pack JD
 ↵