Stability and anisotropy of (FexNi1−x)2O under high pressure and implications in Earth’s and super-Earths’ core

Oxygen is thought to be an important light element in Earth’s core but the amount of oxygen in Earth’s core remains elusive. In addition, iron-rich iron oxides are of great interest and significance in the field of geoscience and condensed matter physics. Here, static calculations based on density functional theory demonstrate that I4/mmm-Fe2O is dynamically and mechanically stable and becomes energetically favorable with respect to the assemblage of hcp-Fe and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R\bar{3}m$$\end{document}R3¯m-FeO above 270 GPa, which indicates that I4/mmm-Fe2O can be a strong candidate phase for stable iron-rich iron oxides at high pressure, perhaps even at high temperature. The elasticity and anisotropy of I4/mmm-(FexNi1−x)2O at high pressures are also determined. Based on these results, we have derived the upper limit of oxygen to be 4.3 wt% in Earth’s lower outer core. On the other hand, I4/mmm-(FexNi1−x)2O with high AV S is likely to exist in a super-Earth’s or an ocean planet’s solid core causing the locally seismic heterogeneity. Our results not only give some clues to explore and synthesize novel iron-rich iron oxides but also shed light on the fundamental information of oxygen in the planetary core.

In addition to oxygen-bearing iron alloys in Earth's core, iron oxides have a widespread occurrence in Earth's crust and mantle, which have a significant influence on the oxidation state and the phase balance of Earth's interior. Therefore, extensive experiments have been conducted at relevant P-T conditions of Earth's interior to investigate structural and physical properties of iron oxides 13,[16][17][18][19][20][21] . In particular, recent studies have successfully synthesized a series of iron oxides with new stoichiometries at high P-T conditions such as Fe 4 O 5 , Fe 5 O 6 , Fe 5 O 7 , Fe 7 O 9 and FeO 2 [21][22][23][24][25][26] . These findings have revealed a complex phase diagram of the Fe-O system at extreme conditions and indicated different scenarios in Earth's interior. For example, Fe 4 O 5 and Fe 5 O 6 as well as those well-known iron oxides are stable and some of them are likely to coexist from 10 GPa to 20 GPa resulting in a complicated oxygen buffer in Earth's transition zone 23,27,28 . However, it is obvious that these newly synthesized compounds are all oxygen-rich iron oxides with Fe/O < 1. Heretofore, investigations into iron-rich iron oxides with Fe/O > 1 are rare. Although some theoretical studies have been performed to explore iron-rich iron oxides under high pressure [29][30][31][32] , their structural and physical properties as well as thermodynamic stability with respect to other components in the Fe-O system are poorly understood.
Here we investigate the stability, elastic and seismic properties of (Fe x Ni 1−x ) 2 O under high pressure by first-principle calculations based on density functional theory (DFT). According to present results and previous data, we further discuss the existence of (Fe x Ni 1−x ) 2 O and its potential effect on geochemical and geophysical processes in Earth's and a super-Earth's interior.

Results
The static calculation results of Fe 2 O are displayed in Fig. 1(a) Supplementary  Fig. S2]. The phonon spectrum without any imaginary frequency implies that the I4/mmm structure is dynamically stable at 320 GPa.
In order to evaluate the relative stability of Fe 2 O versus FeO at pressures corresponding to Earth's core, the following chemical reaction is considered in the present study:  Fig. 1(b)]. These results are inconsistent with previously experimental observations, where B8 is the most stable structure above 100 GPa at moderate temperatures and will transform into B2 above 240 GPa at high temperatures [13][14][15]33,34 . However, our results are in good agreement with previous theoretical calculations 32 . The present calculations are performed at 0 K (static calculations) while previous experiments were conducted at high temperatures. In addition, the ideal chemical stoichiometric compound FeO is considered in our simulations but the non-stoichiometric Fe 1−x O was often used in experiments. For Fe the hcp structure is predicted to be stable at high pressures corresponding to Earth's core compared with the fcc structure [see Supplementary Table S1]. Based on these parameters, we have calculated the relative enthalpy of I4/mmm-Fe 2 O as a function of pressure compared with the assemblage of hcp-Fe and R m 3 -FeO [ Fig. 1(c)]. Figure 1(c) demonstrates that I4/mmm-Fe 2 O becomes energetically favorable above ~270 GPa and is likely to be stable at least up to 400 GPa. Additionally, the enthalpy difference between I4/mmm-Fe 2 O and the assemblage of hcp-Fe and R m 3 − FeO increases with increasing pressure indicating that I4/mmm-Fe 2 O becomes more stable upon compression.
The compression of volumes as a function of pressure displayed in Fig. S3 illustrates a marginal volume change across the chemical reaction (1). The volume reduction of reaction (1) is about 1.45% at 260 GPa and increases with increasing pressure. Previous experiments have detected B2-FeO under ultra-high pressure and temperature conditions 15 Table S1]. (It is worth mentioning that different arrangements of Ni in (Fe 0.5 Ni 0.5 ) 2 O were considered but they were unstable. When the atomic positions and unit-cell parameters were allowed to relax at each given volume, the I4/mmm structure would destroy.) The zero-pressure bulk modulus (K 0 ) of (Fe x Ni 1−x ) 2 O decreases with the Ni content [see Supplementary Fig. S4]. The elastic constants of (Fe x Ni 1−x ) 2 O increase monotonically with pressure [ Fig. 2 and see Supplementary Table S2]. These three components all present mechanically stable in the calculated pressure range supported by the Born-Huang criterion 35 . The substitution of Ni for Fe in (Fe x Ni 1−x ) 2 O tends to significantly decrease C 11 , C 33 , C 44 and C 66 and increase C 12 and C 13 . Figure 2 also reveals a moderate anisotropy in the axial incompressibility along a axis and c axis. The c axis is more compressible than a axis with C 11 > C 33 for Using the aforementioned C ij of (Fe x Ni 1−x ) 2 O, the adiabatic bulk and shear moduli (K S and G) at high pressures are calculated according to the Voigt-Reuss-Hill averages 36 [ Fig. 3(a) and see Supplementary Table S2]. The Ks and G of I4/mmm-(Fe x Ni 1−x ) 2 O increase monotonically with pressure and are all greater than PREM 37 . The substitution of Ni in (Fe x Ni 1−x ) 2 O slightly affects the K S but significantly reduces the G. Furthermore, the aggregate V P and V S of (Fe x Ni 1−x ) 2 O are obtained from their moduli and densities [ Fig. 3(b) and see Supplementary Table S2].  Table S2]. The AV P of Fe 2 O is 9.5% at 260 GPa and decreases to 8.2% at 400 GPa [ Fig. 3(c) and see Supplementary Table S3]. The AV P of (Fe 0.5 Ni 0.5 ) 2 O is similar to that of Fe 2 O while the AV P of Ni 2 O is about 1.7% larger than that of Fe 2 O at each pressure point. These results indicate that the small substitution of Ni for Fe will not change the AV P of (Fe x Ni 1−x ) 2 O significantly in the calculated pressure range. At 320 GPa, the fastest P wave (13.82 km/s) propagates along a axis direction and the slowest P wave (12.65 km/s) distributes in the (110) Fig. 3(d)]. The high AV S mainly results from the high anisotropy of V S2 , though AV S1 is, in fact not small. The layered I4/mmm structure consists of OX 8 blocks (X: Fe or Ni) along c axis contributing to a high AV S , which is similar to that of post-perovskite 40 .  32 . These structures are still unstable with respect to the dissociation, though they are very close to the convex hull. Furthermore, the P m 6 2 structure of Fe 3 O is found to be a mixture of phases consisting of Fe and Fe 2 O and the P4/nmm structure is found to consist of Fe, Fe 2 O and FeO. Therefore, combining previous data with our results, we propose that it is unlikely to form Fe 3 O or Fe 4 O by the chemical combination of Fe and Fe 2 O under Earth's core pressures. Our static calculations indicate that the I4/mmm-type Fe 2 O is a stable phase in iron-rich iron oxides at 0 K. The present static calculations have shown that the I4/mmm-type Fe 2 O is dynamically stable and has the smaller enthalpy and volume than those of the assemblage of Fe and FeO at high pressures, which all favor the formation of Fe 2 O by the combination of Fe and FeO. However, our high-pressure simulations are limited at 0 K in contrast to the high P-T conditions of Earth's core. As shown in Fig. 1(c), the enthalpy difference between the I4/mmm-type Fe 2 O and the Fe + FeO assemblage at high pressure of the ICB is of the order of 0.15 eV corresponding to ~1750 K. In terms of B2-FeO, which has been experimentally confirmed at high P-T conditions, the enthalpy difference is of the order of 0.35 eV (~4100 K). Thus it is likely that the I4/mmm-type Fe 2 O can be stable at high P-T conditions. On the other hand, the entropy is going to play an important role at high temperature but its effect on the relative stability of the target system is not considered in the present study. Previous ab initio simulations have found that when the entropic effect is included in the calculation, the Gibbs free energy of the reaction, 62Fe + FeO → Fe 63 O, can be substantially reduced by ~3 eV 31 . As for the present target system, the introduction of entropic effect may correspondingly favor the combination of Fe and FeO at high P-T conditions. But, on the contrary, it can have an opposite effect causing the reaction energetically unfavorable. That means it might change our conclusion significantly and invalidate our further interpretations for the deep Earth at ultra-high temperature. Further investigations on the present target system by either ab initio molecular dynamics methods or high P-T experiments are definitely required to verify the reaction, Fe + FeO → Fe 2 O. It is to be noted that Earth's core contains diverse light elements. The present of one light element can affect the relative stability and the amount of the other one in Earth's core 8,9 . Thus, the relative stability and physical properties of the iron-rich iron oxide or oxygen-bearing iron alloys should be investigated as a function of not only pressure and temperature but also compositions in order to comprehensively elucidate the nature of oxygen in Earth's core. But they are beyond the scope of the present study and we will work on them in the future.

Discussion and Implications
In terms of phase relations of the Fe-FeO system, Sherman 29 indicated there should be no solid solution between Fe and FeO under Earth's core conditions due to the aforementioned reasons. In contrast, later simulations implied that oxygen dissolved in iron might be stabilized at concentrations up to a few mol% under core conditions because of the significant entropic effect in the dilute solution 31 . Recent thermodynamic simulations suggested the feasibility of the ideal solution model to calculate the Fe-FeO liquid property under outer core conditions and yielded the eutectic compositions of Fe-7.2~9.1 wt% O 12 . Based on their thermodynamic simulations, they concluded that an overall oxygen-rich bulk outer core model should be excluded. In the present study, it is proposed that Fe 2 O can be a strong candidate phase for stable iron-rich iron oxides at high pressure, perhaps even at high temperature, which can potentially change phase relations of the Fe-FeO system at pressures corresponding to Earth's lower outer core.
The accurate measurement of the sound velocity of liquid iron alloys at high pressure is very challenging. Previous studies have discussed the sound velocity of the liquid phase can be addressed from the bulk sound velocity of the solid phase with identical compositions 7 . In addition, theoretical calculations have shown that the bulk sound velocity (V ϕ ) of the solid phase is comparable to V P of the liquid phase (the difference is about 2.5%) and both are marginally dependent of temperature at a given high pressure 41,42 . Indeed, recent experiments have presumed V ϕ of the solid FeS 2 is equivalent to V P of the liquid counterpart and thus estimated the sulfur content in Earth's outer core 43 . Furthermore, in the I4/mmm-type Fe 2 O, the coordination number of oxygen is eight and each iron is surrounded by four oxygen and nine iron, which is similar to the recent report of structural properties of the liquid oxygen-bearing iron alloys by first-principles molecular dynamics 44 . Thus, it is reasonable that the target system can be used to roughly derive the amount of oxygen in the lower outer core. We obtain V ϕ of (Fe x Ni 1−x ) 2 O as a function of pressure [ Fig. 4 and see Supplementary Table S2] and estimate the volume fraction of Fe 2 O by the following relation: where t is the volume fraction of Fe 2 O. We find a volume fraction, t ~ 35% [ Fig. 4]. Then, if oxygen is the only light element in Earth's lower outer core, the maximum possible oxygen content is ~4.3 wt%. However, thermoelastic parameters of (Fe x Ni 1−x ) 2 O at realistic outer core conditions are not well constrained, which certainly warrants further explorations. On the other hand, for a terrestrial super-Earth, unlike Earth, the metal core is mainly composed of solid iron alloys and only a small portion of the core is liquid because of its high mass. Therefore, there will be more oxygen in a super-Earth's solid core than that in Earth's inner core assuming that they have similar cosmochemical compositions. There also exist a type of so-called ocean planets, where a thick layer of ice covers the rocky materials and the center temperature can be much lower than that of super-Earths' 45,46 . Based on current results at 0 K, above 270 GPa, it is energetically favorable for the formation of the I4/mmm-type Fe 2 O. The enthalpy difference becomes larger upon further compression, which can stabilize the I4/mmm-type Fe 2 O even at high temperature. Thus, it is likely that I4/mmm-(Fe x Ni 1−x ) 2 O could exist in the solid core of a super-Earth or an ocean planet. In addition, Fig. 3(d) demonstrates that I4/mmm-(Fe x Ni 1−x ) 2 O has a high AV S in the calculated pressure range due to its layered structure. If (Fe x Ni 1−x ) 2 O can accumulate in the solid core locally, it will cause the locally seismic anisotropy in the solid core, which might be even observed by advanced instruments.

Conclusions
In conclusion, the high-pressure behavior of (Fe x Ni 1−x ) 2 O has been studied based on first-principle density functional calculations. The end member Fe 2 O, is predicted to undergo structural transitions from P6 3 /mmc to P m 3 1 at 212 GPa and further to I4/mmm at 266 GPa, consistent with previous simulations. The dynamically and mechanically stable I4/mmm-Fe 2 O becomes energetically favorable with respect to the assemblage of hcp-Fe and R m 3 -FeO above 270 GPa. These results indicate that I4/mmm-Fe 2 O can be a strong candidate phase for stable iron-rich iron oxides at high pressure, perhaps even at high temperature, which requires further investigations by ab initio molecular dynamics methods or high P-T experiments. The elastic and seismic properties of I4/mmm-(Fe x Ni 1−x ) 2 O at high pressures have been discussed in detail. Combining previous data with present results, we roughly estimate that if oxygen is the only light element in Earth's lower outer core, less than 4.3 wt% oxygen content is required to match the seismic observations. On the other hand, I4/mmm-(Fe x Ni 1−x ) 2 O may exist in a super-Earth's or an ocean planet's solid core. If (Fe x Ni 1−x ) 2 O can accumulate in its solid core locally, it is likely to cause the locally seismic heterogeneity because of its high AV S .

Methods
Three candidate structures (P6 3 /mmc, P m 3 1 and I4/mmm [see Supplementary Fig. S1(a)]) for Fe 2 O were considered in the present study and structural details could be found in Table S1. First-principle calculations were performed based on density functional theory with the projected augmented wave method (PAW) implemented in Vienna ab-initio simulation package (VASP) [47][48][49] . The Perdew-Burke-Ernzerhof (PBE) version of the generalized gradient approximations (GGA) was selected to treat the exchange correlation potential 50 . The kinetic energy cut-off was set to 1000 eV. The energy convergence criterion for the electronic self-consistent calculation was 10 −6 eV. The total energy difference was converged to 1×10 −5 eV/formula unit (f.u.) with respect to the energy cutoff and k-points, respectively. The force difference was converged to 1×10 −3 eV/Å (less than 0.1 GPa). The spin-polarization of iron of Fe 2 O with various structures was not included because the existence of magnetic moments under Earth's core conditions could be safely ruled out. For each crystalline phase, the atomic positions and unit-cell parameters were allowed to relax at each given volume to obtain the minimum total energy. Once the minimum total energies of each phase were obtained at different volumes, they were fitted to the third-order Birch-Murnaghan EoS 51,52 . Additionally, the enthalpy (H = E + PV) of each phase was compared with each other to identify the most stable structure at the given pressure.
The phonon dispersions were calculated using the Phonopy Code by the supercell method 53 and a 2×2×2 supercell was constructed for the I4/mmm-type Fe 2 O. Single crystal elastic constants of I4/mmm-(Fe x Ni 1−x ) 2 O were computed from the stress-strain relations (σ ij = c ijkl ε kl , where σ ij stands for the stresses, c ijkl for the elastic moduli and ε kl for the strains) 54 . For the I4/mmm symmetry, we derived six independent elastic constants (C 11 , C 12 , C 13 , C 33 , C 44 and C 66 ). We applied positive and negative strains of magnitude of 0.5% in order to accurately determine the stresses in the appropriate limit of zero strain.