Computational searches for iron oxides at high pressures

We have used density-functional-theory methods and the ab initio random structure searching (AIRSS) approach to predict stable structures and stoichiometries of mixtures of iron and oxygen at high pressures. Searching was performed for 12 different stoichiometries at pressures of 100, 350 and 500 GPa, which involved relaxing more than 32,000 structures. We find that Fe$_2$O$_3$ and FeO$_2$ are the only phases stable to decomposition at 100 GPa, while at 350 and 500 GPa several stoichiometries are found to be stable or very nearly stable. We report a new structure of Fe$_2$O$_3$ with $P2_12_12_1$ symmetry which is found to be more stable than the known Rh$_2$O$_3$(II) phase at pressures above $\sim$233 GPa. We also report two new structures of FeO, with $Pnma$ and $R\bar{3}m$ symmetries, which are found to be stable within the ranges 195-285 GPa and 285-500 GPa, respectively, and two new structures of Fe$_3$O$_4$ with $Pca2_1$ and $P2_1/c$ symmetries, which are found to be stable within the ranges 100-340 GPa and 340-500 GPa, respectively. Finally, we report two new structures of Fe$_4$O$_5$ with $P4_2/n$ and $P\bar{3}m1$ symmetries, which are found to be stable within the ranges 100-231 GPa and 231-500 GPa, respectively. Our new structures of Fe$_3$O$_4$ and Fe$_4$O$_5$ are found to have lower enthalpies than their known structures within their respective stable pressure ranges.


I. INTRODUCTION
Fe/O is an important planetary mineral and a significant component of the Earth's mantle, and it may also be important within the Earth's core.2][13][14] The identities and concentrations of light elements present in the core have yet to be established, although various arguments favour O, Si, C, S and H as the most likely candidates. 2 In an earlier paper we investigated the formation of iron carbides at high pressures, 3 and here we study iron oxides.O has been suggested as an important light element for explaining both the density jump at the Earth's inner core boundary and the observed seismic wave velocities. 4The possible presence of O within the core has stimulated a range of experimental [4][5][6][7] and theoretical 8,9 investigations of the Fe/O system, albeit only for a handful of stoichiometries.In this paper we present a broader investigation spanning 12 different stoichiometries at three different pressures.We search for stable structures at each stoichiometry and pressure, and construct convex hull diagrams to determine which phases are stable against decomposition.We also present results for transitions between solid phases up to pressures of 500 GPa.
Temperatures within the Earth extend to above 6000 K 10 .Atoms such as O could be present in the Earth's core as isolated impurities, as in a solid solution or, al-ternatively, Fe/O compounds could be formed.A solid solution has a much higher configurational entropy than an ordered Fe/O compound plus crystalline Fe, but the formation of crystalline Fe/O compounds could lead to more favorable bonding arrangements than are obtained in a solid solution.To help understand the conditions under which solid solutions and ordered compounds might be stable it is necessary to determine which compounds are energetically favorable, taking into account a wide range of stoichiometries.This is a formidable task in its own right, although it should be augmented by calculating finite temperature effects arising from the vibrational motion of the atoms, the configurational entropy of the solid solution, and the electronic entropy.We leave the estimation of finite temperature effects to subsequent work.
We have performed searches at 100, 350 and 500 GPa.The pressure at the bottom of the Earth's mantle is about 136 GPa, and our data at 100 GPa correspond to a pressure within the lower mantle.The pressures within the solid inner core are believed to be in the range 330-360 GPa, [11][12][13][14] and our results at 350 GPa are therefore representative of inner core pressures.The highest pressure at which we have performed searches of 500 GPa is above that found within the Earth, although such pressures occur in other planets, including exoplanets.

II. COMPUTATIONAL DETAILS
Determining the most stable structure of a material corresponds to finding the global minimum energy state arXiv:1508.05247v1[cond-mat.mtrl-sci]21 Aug 2015 within a high-dimensional space.A number of methods have been developed to tackle this problem, which have met with varying degrees of success.The method we use, which has been shown to be successful in determining structures that have subsequently been verified by experiments [15][16][17][18] , is called ab initio random structure searching (AIRSS). 19AIRSS has been applied to many systems at high pressures, including iron 20 and oxygen. 21ithin AIRSS the energy landscape is sampled by generating random structures.The structures are then relaxed to the local enthalpy minimum for some chosen pressure using the castep plane-wave density functional theory (DFT) code. 22This process is repeated until the lowest enthalpy structure is found several times or the available computational resources have been exhausted.
To improve the efficiency of the searches we can impose a number of biases.One bias we use is to reject initial structures in which the smallest inter-atomic separation is less than some chosen value.Structures which contain atoms that are very close together are extremely high in enthalpy and may even cause the geometry optimization procedure to fail.Since the number of local minima is known to increase exponentially with the system size 23 , an additional bias may be required for searches over larger unit cells.We use a bias that exploits the fact that low enthalpy structures tend to possess symmetry.We can create an initial structure with a particular symmetry by first selecting a "structural unit" containing A formula units (fu).This structural unit is chosen to be the lowest enthalpy structure found in an A fu search.A random space group containing B symmetry operations is then chosen and applied to the structural unit.This generates a larger initial structure containing A × B fu.In this investigation A and B range from 1-3 and 2-4, respectively.
Structures found to be metastable at one pressure may become stable at another.Expanding the enthalpy around the pressure p 0 at which a search is performed gives The second-order derivative in Eq. ( 1) is related to the bulk modulus and is computationally expensive to evaluate.A useful approximation can, however, often be obtained by neglecting the second-order term and using the fact that the first derivative of H with respect to pressure p at p 0 is equal to the volume v 0 , so that Since H(p 0 ) and v 0 are evaluated within the structure searching, Eq. ( 2) may be used to estimate the enthalpies of the phases at any pressure p.This allows us to determine approximate pressures at which phase transitions may occur, which we then refine by performing higher accuracy DFT calculations.We used ultrasoft pseudopotentials 24 to represent the cores of the Fe and O atoms.For the structure searching we used a fairly soft Fe pseudopotential in which only the 4s and 3d electrons were treated explicitly, while for O we treated the 2s and 2p electrons explicitly.For the higher quality calculations we used a harder Fe pseudopotential in which the 3s, 3p, 4s, and 3d electrons were treated explicitly.Both the harder Fe pseudopotential and the O pseudopotential have been successfully tested up to terapascal pressures. 20,21,25he final structures were obtained using a two-step procedure.We first used AIRSS in medium-quality calculations which were optimized for computational speed.The medium quality settings consisted of a k-point sampling grid of spacing of 2π×0.07Å−1 and a plane-wave cut-off energy of 490 eV.The lowest enthalpy structures for each stoichiometry were then further relaxed in a higher-quality calculation using a finer k-point sampling grid of spacing 2π×0.03Å−1 and a plane-wave cut-off energy of 1000 eV.In all of our searches we used the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) exchange correlation functional. 26he Fe/O system is known to be problematic for DFTbased methods because of magnetic and strong correlation effects.The presence of magnetism in FeO at high pressures has been documented, although there are some conflicting reports.Mössbauer spectroscopy has shown that a continuous transition from high-spin to low-spin states 27 occurs at 90-120 GPa, leading to an eventual collapse at 140 GPa and 300 K.This is in contrast to experiments performed using X-ray emission spectroscopy for FeO at 300 K, which have suggested that a magnetic state exists up to at least 143 GPa 28 , while a magnetic collapse in Fe 2 O 3 has been observed at 50 GPa. 29n addition, zero-temperature static lattice GGA DFT calculations have predicted a magnetic collapse of FeO at around 200 GPa. 30he existence of magnetic order under core conditions can be safely ruled out, as it would be destroyed by the high temperatures.We have not performed searches with spin-polarization because they are very costly, but to investigate the effects of magnetism on the enthalpies, we have re-relaxed our structures with high quality settings starting from a high spin state.If a relaxation did not achieve convergence we reduced the initial spin state and re-relaxed until convergence was achieved.We found that the small magnetic moments have a negligible effect on the enthalpies at the high pressures considered here (∼0.002 eV per atom at 100 GPa in the worst case).For these reasons we do not believe magnetism plays a significant role at the pressures we have investigated and hence we have not considered it further.
Similar arguments can be made about the effects of strong Coulomb correlations which arise in iron oxides from interactions involving the localized Fe d electrons.
The strength of such correlations is conventionally measured by the parameter U/W , where U is the Hubbard U energy and W is the d band width.Electrons tend to delocalize at high pressures, which has the effect of increasing W and reducing U . 30Under these conditions, where the correlations are less strong, GGA DFT is likely to be reasonably accurate. 31

III. CONVEX HULLS
We have plotted our energy data on convex hull diagrams to assess the stability of the phases with respect to decomposition.To construct the hulls we also need the enthalpies of the end members, pure Fe and pure O. Fe is predicted to adopt the hexagonal close packed (hcp)  structure over the pressure range of interest here, while for oxygen at 500 GPa we use the P 6 3 /mmc phase predicted by Sun et al. 21At pressures below 500 GPa we use the ζ phase. 32he enthalpy of formation per atom of a compound with respect to its elements is where N x and H x are the number of atoms and enthalpy per atom of element x, respectively, and H s is the enthalpy of the structure.A negative ∆H s indicates that the structure is stable with respect to decomposition into its elements, although it may be unstable with respect to decomposition into compounds of other stoichiometries.
Structures lying above the convex hull are unstable with respect to decomposition into other compounds and are therefore less likely to form.Structures lying on the convex hull are thermodynamically stable and are therefore more likely to form.The convex hull diagrams at 100, 350 and 500 GPa shown in Fig. 1 exhibit fairly smooth variations with composition.Figure 1(a) shows that the only stoichiometries that lie on the hull at 100 GPa are Fe 2 O 3 and FeO 2 .At 350 GPa, Fe 9 O, Fe 3 O, Fe 2 O, Fe 3 O 2 , FeO, Fe 2 O 3 and FeO 2 are found to lie on, or very close to the hull, and similar results are found at 500 GPa.Perhaps the most notable difference at 500 GPa is that FeO 4 is found to lie on the hull.FeO 2 and Fe 2 O 3 are predicted to be stable with respect to decomposition at all pressures investigated.Fe 3 O 4 , which is believed to be an important phase in the Earth's lower mantle 34 , is found to be marginally unstable with respect to decomposition at the pressures considered here.
Increasing pressure appears to stabilize Fe rich phases, an effect that is most noticeable between 100 and 350 GPa.There is also a stabilization of O rich phases with increasing pressure.For example, FeO 3 and FeO 4 can be seen to move considerably closer to the hull between 100 GPa and 500 GPa.The former also coincides with a transition from a low symmetry P 1 phase to a phase of Cmcm symmetry.
It is possible that structures that lie just above the convex hull correspond to metastable mixtures of phases rather than single phases.This is particularly likely for large unit cells.We have found several such cases in our results, and the fact that they lie just above the hull suggests that the interfacial energies between the different structures present are small.The P 6m2 structure of Fe 3 O is found to be a mixture of phases consisting of Fe and Fe 2 O while the P 4/nmm structure is found to consist of Fe, Fe 2 O and FeO phases.Finally, Fe 9 O (Amm2) is found to consist of Fe and Fe 2 O phases and the P m phase of Fe 3 O 2 (at 500 GPa) is found to consist of Fe 2 O and FeO phases.These mixtures of phases are indicated by square brackets in Table I.

IV. PHASE TRANSITIONS
We predict pressure-driven phase transitions for all of the stoichiometries considered in this study, with the exception of Fe 9 O.We have obtained accurate transition pressures by performing additional calculations at pressures around the values predicted using Eq.(1).Plots of relative the enthalpies as a function of pressure are shown in Fig. 4. Three of the stoichiometries considered, Fe 2 O 3 , FeO and Fe 3 O 4 , have already been investigated extensively and their properties, particularly at ambient pressure, are well known.However, most of the stoichiometries that we discuss here have not been studied before and hence little is known about their high pressure behaviour.Unless stated otherwise, all of the lowest enthalpy phases were found by searching and, to the best of our knowledge, have not been discovered before.
Fe 9 O.The most stable structure of Fe 9 O found has Amm2 symmetry, but it turns out to be a mixture of phases.
Fe 3 O.See Fig. 4(a).Between 100 and 350 GPa, the most stable structure of Fe 3 O of P 6m2 symmetry turns out to correspond to a mixture of phases rather than a single phase.At about 350 GPa we find a phase of P4/nmm symmetry to be the most stable, but it also turns out to be a mixture of phases.We also relaxed Fe 3 O between 100 and 500 GPa in the BiI 3 structure, which has previously been suggested as a low enthalpy phase. 8However, we found this structure to be between 0.3 eV and 2 eV per fu higher in enthalpy than our best structures found from searching.
Fe 2 O. See Fig. 4(b).We predict that Fe 2 O forms a structure of P 6 3 /mmc symmetry at 100 GPa.At about 180 GPa we find a transformation to a phase of P 3m1 symmetry, which is stable up to about 288 GPa, at which point we predict it to transform into a phase of I4/mmm symmetry.
Fe 3 O 2 .See Fig. 4(c).At 100 GPa we predict that Fe 3 O 2 will adopt a structure of P m symmetry.At about 127 GPa, we predict a phase transition to a structure of C2/m symmetry followed by a further transition to a phase of P m symmetry at about 486 GPa.
FeO. See Fig. 4(d).Previous DFT calculations for FeO have shown the stability of a phase of P 3 2 21 symmetry 35 between 65 GPa and 190 GPa.We also find this structure from searching and our results shown in Fig. 4(d) agree well with those previously reported.However, Oganov et al. 35 predict that this phase will transform to the distorted NaCl structure at 190 GPa.We instead find that the P 3 2 21 structure transforms to a new phase of P nma symmetry at about 195 GPa, which then undergoes a transformation to another new phase of R 3m symmetry at about 285 GPa.In addition to this, Fig. 4(d) shows another new phase of Cmcm symmetry, which we predict to be stable at pressures above 500 GPa.Our results suggest a somewhat different picture from that obtained from compression experiments up to 324 GPa 4 and at temperatures up to 4880 K.These differences may arise for a number of reasons.Firstly, experiments on FeO normally use a non-stoichiometric compound, Fe 1−x O, where x is typically about 0.05, whereas our calculations assume stoichiometric FeO.Secondly, we have performed static lattice calculations in contrast to real-world experiments which may yield phases stabilized by nuclear vibrational motion.Thirdly, kinetic barriers between structures could also influence the structures that are observed in experiments.Finally, the use of an approximate density functional will lead to further uncertainties in our results.It is worth remarking, however, that the PBE density functional, which obeys the uniform limit and gives a good account of the linear response of the electron gas to an external potential, 26 is likely to become more accurate as the pressure is increased.
A notable example of where our predictions for FeO differ from observations is at relatively low pressures where experiments (extrapolated to 0 K) have shown that FeO adopts the rhombohedrally distorted NaCl (rB1) structure, which is predicted to be stable up to about 100 GPa. 4 Above 100 GPa a phase transition to the NiAs (B8) structure is predicted.We find the inverse-B8 (iB8) structure to be more stable than B8 between about 100 GPa and 165 GPa.The stability of the iB8 structure has been predicted theoretically before 31,36 , although it has yet to be verified by experiment.Although our results correctly predict the B8 phase to be more stable than the rB1 phase at high pressures, our best structures obtained from searching are lower in enthalpy than the previously known phases across the entire pressure range considered.We are not aware of any experimental evi-    Relative Enthalpy per fu (eV) Relative Enthalpy per fu (eV) FIG. 4: (Color online).Plots of relative enthalpy against pressure for various stoichiometries.The reference structures appear as a horizontal line at 0 eV.Dashed lines represent structures known prior to this work.
dence suggesting that a P 3 2 21 phase of FeO is stable at 100 GPa, but we believe that this discrepancy may be explained by our arguments in the previous paragraph, and we note that earlier work suggested that the predicted stability of the P 3 2 21 phase may simply be an artefact of the PBE density functional. 35 4 O 5 .See Fig. 4(e).Compression experiments have identified a phase of Fe 4 O 5 with Cmcm symmetry to be the most stable between 5 and 30 GPa. 33 Our results show this phase to be unstable between 100 and GPa with respect to the new structures found in our searches.At 100 GPa we predict Fe 4 O 5 to adopt a structure with P 4 2 /n symmetry prior to a transition to a structure of P 3m1 symmetry at about 231 GPa, which we predict to be stable to at least 500 GPa.
Fe 3 O 4 .See Fig. 4(f).Under ambient conditions Fe 3 O 4 adopts an inverse spinel structure of F d 3m symmetry. 37etween ∼20 − 60 GPa, Fe 3 O 4 is observed to make a slow transformation to a new phase (sometimes denoted by h-Fe 3 O 4 ) with orthorhombic symmetry. 34,37The exact structure of h-Fe 3 O 4 has, however, been difficult to resolve. 34In light of this, a number of candidate structures have been proposed.Fei et al. proposed a structure of P bcm symmetry 38 which, although it agreed well with the available X-ray diffraction data, was later found to be inconsistent with Mössbauer spectroscopy data. 34It was later suggested that a structure of Cmcm symmetry (sometimes referred to using the non-standard Bbmm setting) is also consistent with the X-ray diffraction data and it currently seems to be favoured as the high pressure phase of Fe 3 O 4 . 39First principles calculations have shown the Cmcm phase to be more stable than the P bcm phase at pressures > 65 GPa. 37Our searches at 100 GPa found a new phase with P ca2 1 symmetry to be the lowest in enthalpy up to ∼ 340 GPa, above which a phase transition is found to occur to another new phase of P 2 1 /c symmetry.Figure 4(f) shows that the Cmcm phase is energetically unfavourable at the pressures considered.The P ca2 1 structure is found to be much more stable than the Cmcm phase, being ∼ 1.4 eV per fu lower in enthalpy at 100 GPa.The fact that the P ca2 1 phase is also of orthorhombic symmetry does not preclude it from being a candidate for h-Fe 3 O 4 .Ascertaining whether the P ca2 1 phase is indeed h-Fe 3 O 4 would require a more detailed investigation, which we leave to subsequent work.
Fe 2 O 3 .See Fig. 4(g).At ambient pressures the stable phase of Fe 2 O 3 is hematite, which is a corundum-type structure of R 3c symmetry. 40The high pressure phase of Fe 2 O 3 has been proposed to be either a Rh 2 O 3 (II)-type structure (of Pbcn symmetry) or a GaFeO 3 orthorhombic perovskite structure. 40Room temperature compression experiments have shown that hematite transforms to the Rh 2 O 3 (II)-type structure at about 50 GPa. 41However, compression experiments at higher temperatures (800-2500 K) predict a phase transition to occur between hematite and the high pressure phase at a much lower pressure of ∼ 26 GPa. 42Ono et al. have suggested that the observed transition at room temperature is between hematite and metastable Rh 2 O 3 (II). 40,42The perovskitetype phase is suggested to be kinetically inhibited and therefore only obtainable at high temperatures. 40,42This view is supported by hybrid-DFT calculations, which show that the Rh 2 O 3 (II) phase is metastable with respect to the perovskite phase at 50 GPa. 40At about 60 GPa and temperatures > 1200 K the high pressure phase is observed to transform into a CaIrO 3 -type (post-perovskite) structure. 40,43Hybrid-DFT calculations have predicted the anti-ferromagnetically ordered CaIrO 3 type structure to be the most stable between about 46 and 90 GPa with the low-spin Rh 2 O 3 (II) phase gradually becoming stable as the pressure approaches 90 GPa. 40It appears likely that this phase will become stable at around 100 GPa and hence become consistent with our results shown in Fig. 4(g).We found the (non-magnetic) Rh 2 O 3 (II) structure of space group Pbcn in our searches at 100 GPa.Our calculations find this structure to be more stable than both the CaIrO 3 and GaFeO 3 (non-magnetic) structures between 100 and about 233 GPa.We predict that the Rh 2 O 3 (II) phase transforms into a new phase of P 2 1 2 1 2 1 symmetry at 233 GPa.FeO 2 .See Fig. 4(h).At 100 GPa FeO 2 we predict a structure of P a 3 symmetry, see Fig. 2(d), which is stable across a wide pressure range of about 100-465 GPa.In this structure the Fe atoms form a face centred cubic (fcc) arrangement.At 465 GPa the P a 3 phase is predicted to transform into a phase of R 3m symmetry.
Fe 3 O 7 .See Fig. 4(i).Between about 100 and 160 GPa we predict a structure of P 2/m symmetry, and above 160 GPa we expect it to transform into a phase of I 43d symmetry.
FeO 3 .See Fig. 4(j).We find a structure of P 1 symmetry to be the most stable above 100 GPa, and we predict that it will transform into a different stucture of P 1 symmetry above 162 GPa.At about 395 GPa this structure is predicted to transform into a phase of Cmcm symmetry.
FeO 4 .See Fig. 4(k).We find a structure of F dd2 symmetry to be the most stable at 100 GPa, which is predicted to transform into a structure of P 2 1 /c symmetry at about 108 GPa.
A summary of the structures found is given in Table I, and the details of the structures are reported in the Supplemental Material. 44

V. DENSITY OF ELECTRONIC STATES
An analysis of the electronic densities of states of the phases was performed with the lindos code 45 which uses a linear extrapolation method. 46,47A Gaussian smearing of width 0.3 eV was applied to the density of states in all cases.Figure 3 shows the density of states for all of the structures found to lie on the convex hull.In all cases, with the exception of FeO 4 at 500 GPa, a significant density of states is observed at the Fermi energy, indicating that these structures are metallic.
The Fe 3s and Fe 3p states lie at, respectively, ∼87 and ∼54 eV below the Fermi energy in all Fe-O phases.We found the O 2s states to lie at ∼20 eV below the Fermi energy and almost all of the occupied electronic density of states around the Fermi energy arises from the Fe d bands, except in the O-rich structures.

VI. CONCLUSION
We have used DFT methods and the AIRSS technique to identify high pressure structures in the Fe/O system at zero temperature across a wide range of stoichiometries.This is a crucial preparatory stage for modelling the Fe/O system at finite temperatures.By constructing the convex hull at different pressures we have determined energetically favourable stoichiometries, some of which have not been considered before.Our results broaden the range of possible stable stoichiometries considerably from those studied previously.
We have found new structures of FeO and Fe 2 O 3 which are calculated to be stable against decomposition at 350 and 500 GPa.Structures with stoichiometries Fe 2 O 3 and FeO 2 are found to be stable against decomposition at each of the three pressures studied.Fe 9 O, Fe 3 O, Fe 2 O, Fe 3 O 2 , FeO, Fe 2 O 3 and FeO 2 are found to lie on, or very close to the hull at both 350 GPa and 500 GPa, although some of these structures are mixtures of phases, see Table I.At 500 GPa FeO 4 is also found to lie on the hull.
We have found structures of Fe 3 O 4 and Fe 4 O 5 at high pressures that are more stable than those predicted previously, although they are unstable to decomposition at the pressures studied here.
Increasing pressure tends to stabilize Fe/O compounds with respect to the separated elements over the range 100-500 GPa.Increasing pressure stabilizes Fe rich phases, an effect that is most noticeable between 100 and 350 GPa.There is also a stabilization of O rich phases with increasing pressure.
The energy reduction from compound formation in Fe/O at high pressures is much larger than in Fe/C.For example, at 350 GPa, the DFT results reported in Ref. 3 show that the minimum in the convex hull occurs at Fe 2 C with an enthalpy reduction, compared with separated Fe and C solids, of about 0.26 eV per atom, while in Fe/O we find an enthalpy reduction of about 1.76 eV per atom for FeO 2 .Our results show that Fe/O compounds are stable (or nearly stable) over a wider range of stoichiometries at 350 GPa than Fe/C compounds.
FIG. 2: (Color online).Structures found from searching that lie on the convex hull.The O atoms are shown in red and the Fe atoms are in blue.

FIG. 3 :
FIG.3:(Color online).Electronic densities of states for structures lying on the convex hull.Each plot shows the total density of states and its decomposition into contributions from Fe and O.The Fermi energies are at 0 eV.

TABLE I :
The lowest enthalpy structures found from searching at 100, 350 and 500 GPa.X is the fraction of oxygen atoms.Entries which show space groups within square brackets indicate mixtures of structures of phases with other stoichiometries, see text.N is the maximum number of formula units used in the searches.S P is the space group of the lowest enthalpy structure found at pressure P .In each case, the number of formula units in the primitive cell is given in round brackets.