Prediction of crystal structures and motifs in the Fe-Mg-O system at Earth's core pressures

Fe, Mg, and O are among the most abundant elements in terrestrial planets. While the behavior of the Fe-O, Mg-O, and Fe-Mg binary systems under pressure have been investigated, there are still very few studies of the Fe-Mg-O ternary system at relevant Earth's core and super-Earth's mantle pressures. Here, we use the adaptive genetic algorithm (AGA) to study ternary Fe$_x$Mg$_y$O$_z$ phases in a wide range of stoichiometries at 200 GPa and 350 GPa. We discovered three dynamically stable phases with stoichiometries FeMg$_2$O$_4$, Fe$_2$MgO$_4$, and FeMg$_3$O$_4$ with lower enthalpy than any known combination of Fe-Mg-O high-pressure compounds at 350 GPa. With the discovery of these phases, we construct the Fe-Mg-O ternary convex hull. We further clarify the composition- and pressure-dependence of structural motifs with the analysis of the AGA-found stable and metastable structures. Analysis of binary and ternary stable phases suggest that O, Mg, or both could stabilize a BCC iron alloy at inner core pressures.


I. INTRODUCTION
O, Fe, Si, Mg, Al, and Ca (CMAS+F) are the most abundant elements in terrestrial planets [1]. Among these planets, Earth provides essential general information, yet it is incompletely deciphered. All CMAS+F elements are lithophile (rock-loving) elements and are present in the Earth's rocky mantle and crust. Fe is the predominant element in the core and is a siderophile (metal-loving) element as well. Based on current knowledge, this classification is believed to be valid in the pressure and temperature (PT) range achieved in Earth's interior. Seismology and high-pressure data on iron shows that the Earth's core is ∼ 5-10 wt% [2] less dense than iron at expected conditions, i.e., ∼ 136-364 GPa and ∼ 4,000-6,500 K [3][4][5][6][7][8]. This indicates the presence of lighter elements in the core partitioned differently between its solid and liquid regions [9]. Extensive research has been carried out experimentally and computationally to shed light on the light elements' nature. Despite much progress, there is still a great deal of uncertainty [2,[10][11][12]. With every experimental or computational development, this question is revisited from a different angle. In particular, the development of materials discovery methods [13][14][15][16][17][18][19][20] has propelled the exploration of novel chemistries under pressure, which has fueled the debate on the possible nature of light elements in the core. Among the CMAS elements, the most likely light element candidates in the core are Si, and O. O is considered a required element today, a view that has evolved within the last decade [21][22][23][24]. The volatile elements S, C, and H are also regarded as likely candidates, but the abundances of C and H on Earth are still largely unconstrained. Mg, Al, and Ca have been branded lithophile elements up to core pressures. But the recent computational discovery of Mg-Fe compounds up to inner core pressures [25] suggests the possibility of Mg turning siderophile and its presence in the core. The formation of Fe-O [16] and Fe-Si [15] compounds with variable stoichiometry have been investigated using materials discovery methods. The theoretical prediction of pyrite-type FeO 2 [16] and its experimental confirmation [26,27] has been one of the greatest successes of this approach, which rarely explores the possibility of ternary compounds [17]. Given the present understanding that O is a required element in the outer core and should also exist in the inner core, we explore the possible formation of Fe-Mg-O compounds at typical core pressures of ∼ 350 GPa. The investigation of solid compounds provides the most critical information. Light elements must be present in both solid and liquid phases but more abundantly in the liquid phase. It is energetically more costly to accommodate these elements in the solid phase, a geochemical definition of incompatible elements. Therefore, the discovery of thermodynamically stable Fe-Mg-O solids is essential to investigate Mg's presence in the core.
Besides being essential for addressing Mg's presence in the Earth's core, the present study of Fe-Mg-O solids has significant ramifications for the mantle of terrestrial arXiv:2102.03402v1 [cond-mat.mtrl-sci] 5 Feb 2021 planets larger than Earth whose interiors can reach much higher pressures and temperatures. The B1-type isomorphous alloy (Mg 1−x Fe x )O, ferropericlase (x < 0.5) or magnesiumwüstite (x > 0.5) is the second most abundant phase of the Earth's lower mantle. Ferropericlase is the thermodynamically stable form of the Fe-Mg-O ternary compound up to ∼ 135 GPa and ∼ 4,000 K, i.e., CMB conditions in the Earth. The higher pressures and temperatures expected in Super-Earths, e.g., ∼ 4 TPa and ∼ 9,000 K at the CMB in a 20 M terrestrial planet [28], raises the possibility of other compounds and alloy structures with other compositions. The existence of other stable ternary phases under pressure may induce decomposition and recombination reactions between Fe-Mg-O, Fe-O, Mg-O, and Fe-Mg compounds under pressure, similar to what has been observed in the Si-Mg-O system [17,18]. Therefore, the current research can also provide the first glimpses on the essential (Mg 1−x Fe x )O alloy behavior in these planets. This paper is organized as follows. In the next section, we describe the computational method used in this work. Sec. III presents the crystal structure search results, the ternary convex hull, and the predominant structural motifs in low enthalpy structures. Section IV discusses some potential geophysical implications of these results. Conclusions are present in Section V.

II. METHODS
Crystal structures of Fe-Mg-O at high pressure were investigated using the adaptive genetic algorithm (AGA) [17,29]. This method integrates auxiliary interatomic potentials and ab initio calculations adaptively. The auxiliary interatomic potentials accelerate crystal structure searches in the genetic algorithm (GA) loop. At the same time, ab initio calculations are used to adapt the potentials after several GA generations to ensure accuracy. The structure searches were only constrained by the chemical composition, without any assumption on the Bravais lattice type, symmetry, atomic basis, or unit cell dimensions. In our AGA searches, the enthalpy was used as the selection criteria for optimizing the candidate pool. The candidate structure pool size in GA search is 128. At each GA generation, 32 new structures are generated from the parent structure pool via a mating procedure described in [30]. The structures in the pool were updated by keeping the lowest-energy 128 structures. The structure search with a given auxiliary interatomic potential sustained 600 consecutive GA generations. Then, 16 structures from the GA search were randomly selected for ab initio calculations to re-adjust the interatomic potential parameters for the next round of the GA search. This sequence of steps was repeated 40 times. In the AGA search in the Fe-Mg-O system, interatomic potentials based on the embedded-atom method (EAM) [31] were chosen as the auxiliary classical potential. In EAM, the total energy of an N-atom system is described by where ϕ (r ij ) is the pair term for atoms i and j at a distance r ij . F i (n i ) is the embedded term with electron density term n i = j =i ρ j (r ij ) at the site occupied by atom i. The fitting parameters in the EAM formula are chosen as follows: The parameters for Fe-Fe and Mg-Mg interactions were taken from the literature [32].
where D, α, r 0 are fitting parameters. The density function for O atoms are modeled by an exponentially decaying function, where α and β are fitting parameters. The form proposed by Benerjea and Smith [33] was used as the embedding function with fitting parameters F 0 , γ as For Fe and Mg, the parameters of the density function and embedding function were taken from ref. [32] as well. In the AGA scheme [17], the potential parameters were adjusted adaptively by fitting to the ab initio energies, forces, and stresses of selected structures. The fitting process was performed using the force-matching method with a stochastic simulated annealing algorithm implemented in the POTFIT code [34,35].
Ab initio calculations were carried out using the projector augmented wave (PAW) method [36] within density functional theory (DFT) as implemented in the VASP code [37,38]. The exchange and correlation energy are treated without the spin-polarized generalized gradient approximation (GGA) and parameterized by the Perdew-Burke-Ernzerhof formula (PBE) [39]. A plane-wave basis was used with a kinetic energy cutoff of 520 eV, and the convergence criterion for the total energy was set to 10 −5 eV. Monkhorst-Pack's sampling scheme [40] was adopted for Brillouin zone sampling with a k-point grid of 2π × 0.033Å −1 , and the unit cell lattice vectors (both the unit cell shape and size) are fully relaxed under fixed pressure (200 GPa and 350 GPa) together with the atomic coordinates until the force on each atom is less than 0.01 eV/Å. The phonon dispersions were computed with density functional perturbation theory (DFPT) implemented in the VASP code and the Phonopy software [41]. The formation enthalpy (H f ) of compound Fe x Mg y O z was calculated as   32 atoms per primitive cell). After the AGA search, we use the following method to determine the stability of compounds. For a ternary compound A x B y C z , we select three existing compounds A x1 B y1 C z1 , A x2 B y2 C z2 , and A x3 B y3 C z3 on the diagram; these can also be elementary or binary end-members. If A x B y C z can be written as If ∆H ≤ 0 for all scanned combinations of existing phases on the diagram, A x B y C z is determined as an energetic ground state. The energetic ground states form the convex hull as shown in Fig. 1(a). H d is introduced as the enthalpy above the convex hull to represent the relative stability on the phase diagram. By definition, all the ground-state phases have The AGA search found three new ternary ground state compounds at 350 GPa: FeMg 2 O 4 , Fe 2 MgO 4 , and FeMg 3 O 4 . They define the current Fe-Mg-O ternary phase diagram. The stability of these phases from 200 GPa and 350 GPa is shown in Fig. 1(b). This stability pressure range is computed by considering the relative stability of these phases against decomposition into all end-members (see Supplementary Fig. 1 for the stability range of all ground-state phases). We will discuss the construction of the phase diagram in the next section. We also identify low-enthalpy metastable structures such as FeMgO 3 with enthalpy very close to the convex hull (H d = 18 meV/atom). Here we first analyze these new ground states and low-enthalpy structures. Figure 2 shows the atomic structure, phonon dispersion, and electronic density of states for tetragonal FeMg 2 O 4 with space group I-42d. Fe and Mg atoms are coordinated with eight oxygen atoms to form a similar MO 8 (M for metal) polyhedra. Unlike a typical cubic polyhedron, this MO 8 consists only of triangular faces. These triangular faces form pentagonal caps and are similar to the Frank-Kasper polyhedra [42]. The Fe-and Mg-centered MO8 polyhedra pack in various edge-and face-sharing arrangements. This structure is the same as the I-42d -type Mg 2 SiO 4 found previously [17,43]. As shown in Fig. 1(b), this phase becomes the ground state at 349 GPa. Below this pressure, it decomposes into MgO and FeO 2 . As shown in Fig. 2(b), there are no imaginary phonon mode frequencies, which confirms this phase is dynamically stable. The electronic density of state in   Figure 3 shows the atomic structure, phonon dispersion, and electronic density of states for monoclinic Fe 2 MgO 4 with space group Cc. By inspecting the structure, we identify the same type of MgO 8 polyhedra found in I-42d FeMg 2 O 4 . However, the Fe-O polyhedra are more complicated. It contains two polyhedral types, one six-fold and one seven-fold coordinated as shown in Fig. 3(a). The FeO 6 is a highly distorted octahedron. The FeO 7 also shows the pentagonal cap similar to the MO 8 polyhedron found in I-42d FeMg 2 O 4 . As shown in Fig. 1(b), this phase's stability pressure range against the decomposition into MgO and Fe 2 O 3 starts at 325GPa. Phonon calculations confirm its dynamic stability. The electronic density of states also indicates this phase is metallic in Fig. 3(c).
The FeMg 3 O 4 Cmmm structure in Fig. 4 shows Fe-O and Mg-O octahedral building blocks. Visually it is similar to ferropericlase (Fe 1−x Mg x )O, which has a NaCltype (B1) structure. However, unlike the cubic structure and random Fe/Mg cation site occupancies of ferropericlase, FeMg 3 O 4 is an orthorhombic structure with ordered Fe/Mg site occupancies. The octahedra in the FeMg 3 O 4 structure are Jahn-Teller distorted because of the orthorhombic symmetry. This phase becomes stable against decomposition into FeO and MgO at ∼ 228 GPa. Phonon calculations in Fig. 4(b) also confirm its dynamical stability. Unlike the metallic Fe 2 MgO 4 and FeMg 2 O 4 phases, FeMg 3 O 4 Cmmm is a semiconductor.
Besides the ternary ground states, we also analyze a metastable FeMgO 3 Immm structure with H d = 18 meV/atom. This enthalpy difference is so small that the compound may become stable at high tempera-  tures. Structural analysis of FeMgO 3 Immm in Fig. 5 shows an interesting combination of various polyhedra. Iron shows three different oxygen coordination pulyhedra FeO 6 , FeO 7 and FeO 8 . The FeO 6 is an octahedron. The FeO 7 is a trigonal prism with an extra rectangular face capping neighbor. The FeO 8 is a cube. MgO shows one coordination polyhedron type, MgO 8 , not a cube but a triangular prism with two rectangular face capping oxygens. Such a combination of octahedra and prisms is similar to the building blocks in the complex   [44]. The appearance of cubic polyhedron is consistent with the observation of the CsCltype (B2) structure of FeO at Earth's core conditions [11]. This FeMgO 3 Immm structure shows an intermediate packing between NaCl-type FeO/MgO phase to the CsCl-type FeO/MgO phases. This phase is dynamically stable and metallic, as shown by the phonon dispersion and the electronic density of states in Fig. 5.

B. Construction of the Fe-Mg-O ternary convex hull
In this section, we discuss the construction of the ternary phase diagram and convex hull. In a binary sys-tem, the compositional space is one-dimensional so that the convex hull is a curve connecting the formation enthalpies of ground-state phases. In a ternary system, the compositional space is two-dimensional, and the convex hull consists of surface segments connecting the formation enthalpies of three stable phases, as shown in Fig. 1. In a discrete compositional space, these surface segments are triangles. Any new structure having formation enthalpy below this convex hull surface will be a new ground state. The convex hull surface needs to be reconstructed after the discovery of any new stable phase.
For binary references at 350 GPa, the ground-state phases of Fe-O, Fe-Mg, and Mg-O have been investigated in Refs. [16], [25] and [45], respectively (see Supplementary Fig. 2 for their crystal structures). Because these crystal structure searches have already covered the current study's pressure range, we do not perform new AGA searches for the binary phases. Still, we re-calculate the energetics of the previously found crystal structures. Ab initio calculations confirm these reported relative phase stabilities of Fe-O [16], Fe-Mg [25], and Mg-O [45]. By exploring an experimental database [46], we find two previously reported MgFe 2 O 4 stoichiometric compounds [47,48]. However, our calculations indicate these phases are metastable at Earth's core pressures. Based on these ground-state binary phases, we established the Fe-Mg-O ternary system's convex-hull shown in Fig. 6. At 200 GPa, all the AGA searched ternary compounds have relatively higher enthalpy than the elementary or binary ground-state references. Therefore no Fe-Mg-O stoichiometric phase can be a ternary ground state at 200 GPa. At 350 GPa three ternary phases become stable ground states. Detailed energetics and crystallographic information on these ground-state phases is given in Supplementary Tables S1 , Table S2 and Table S3.  Fig. 6 suggests a strong effect of pressure on these phases' stability, we now investigate how the structural motifs change under pressure. Because the current calculation does not include temperature effects on phase relations and at finite temperatures, the ground states may differ. We now focus on the ground-state phases and metastable phases with formation enthalpy within 0.8 eV/atom (∼9,000 K) above the convex hull. Phases in this energy range provide much better statistical information than just ground states. Here we employ a cluster alignment (CA) method [49][50][51] to analyze the motifs of Fe-centered and Mg-centered clusters in each structure. The clusters here are defined by a center atom and its first-shell neighbor atoms. The cluster alignment method identifies the similarity between an as-extracted cluster and the perfect template cluster. Here we use typ-ical motifs in the metallic alloys, including FCC-, BCCand HCP-type clusters as the templates. We also include common motifs in Fe-O binary compounds which are OCT (octahedron) and Prism (trigonal prism) [44]. The snapshots of these motifs are shown in Fig. 7. We first align the atomic cluster against these five motifs to check its similarity with the CA method. If the cluster does not match any of the above five polyhedra, it is marked as an "other" type. In Fig. 7 we show scatter-type plots of enthalpy differences between these structures and that of the convex-hull and their respective atomic volumes. We differentiate the cluster type and O or Fe concentrations with symbols and colors, respectively. Since we obtain similar results between the analysis of Fe-centered clusters and Mg-centered clusters, we only show the analysis of Fe-centered clusters in Fig. 7 and provide the analysis of Mg-centered clusters in Supplementary Fig. S3.
At 200 GPa, Figs. 7(a,c) indicate that crystal structures with octahedral type clusters and 50 mole% oxygen concentration generally have lower energy than the other structures. While there is no ternary ground state at 200 GPa, these structures are concentrated within an energy range very close to the convex hull. Therefore, it is likely that they become ground states at high temperatures. These structures are all similar to the NaCltype B1 structure and have variable Fe/Mg ratios and occupancies, which is essentially (Mg 1−x Fe x )O, i.e. ferropericlase (x F e < 0.5) or magnesiumwüstite (x F e > 0.5). This result is consistent with ferropericlase being a dom-inant phase not only in the Earth's lower mantle but at the higher pressures of Super-Earth's mantles. It's worth noting that the lowest enthalpy structure at 200 GPa is the same as the ground-state FeMg 3 O 4 Cmmm at 350 GPa shown in Fig. 4. Inspecting the higher energy range in Fig. 7(a) and (c), one finds that oxygen-rich structures generally have lower enthalpies than oxygen-poor structures. The oxygen-rich structures mainly contain octahedral clusters, while the oxygen-poor structures can have a greater variety of motifs, including FCC, BCC, HCPtype clusters. This is mainly because Fe and Mg start to alloy to form closely packed motifs under unsaturated oxygen conditions.
At 350 GPa, structures with 50% oxygen concentrations still have the lowest formation enthalpy. A few oxygen-rich phases become more stable and approach the convex hull energetically compared to 200 GPa. BCCtype clusters start to appear in these oxygen-rich structures at 350 GPa, indicating that the B2-type Fe-O clusters are favored at higher pressures over B1-type clusters at lower pressures. The situation with oxygen-poor structures at 350 GPa is similar to the one at 200 GPa. We note that at both 200GPa and 350GPa, several motifs ("other" type) that cannot be classified into the current simple cluster templates appear. Some of them are due to distortions, while some indeed form more complex clusters, e.g., the ones in Fig. 2 and Fig. 3.

IV. GEOPHYSICAL IMPLICATIONS
Our findings on the Fe-Mg-O system at core pressures appear to have some straightforward geophysical consequences. The Fe-rich side (right corner) of the ternary phase diagram in Fig. 1 Fig. 8. BCC-like Fe 2 Mg and hcp -Fe are likely to inter-alloy and form a eutectic system, with two coexisting solid phases for some composition-temperature ranges. Small Mg concentrations might produce hcp-like Fe 1−x Mg x alloys, but a BCC-like Fe 2+x Mg 1−x might precipitate and coexist beyond a certain concentration threshold. The situation is very similar for the Fe 2 O-Fe system. Therefore, the Fe-Mg-O system might contain Mg and O dissolved substitutionally in -Fe for small Mg and O concentrations, but beyond a certain concentration threshold BCC-like Fe 2+x+y (Mg 1/2−x O 1/2−y ) might precipitate. BCC-Fe can be stabilized at inner core pressures by alloying with S [52,53], and it has been argued, but not confirmed, that BCC iron could be stabilized at inner core conditions [54]. Therefore, the precipitation of BCC-like Fe 2+x+y (Mg 1/2−x O 1/2−y ) for non-negligible amounts of Mg, O, or both is not a surprising conclusion.
The ternary phases discovered in the O-rich side (left corner) of the phase diagram are relevant for the mantle of some Super-Earths. The absence of stable ternary phases at pressures lower than 228 GPa suggests that stable phases involving all three elements are solid-solutions of end-member phases with a small concentration of inter-alloying metals. For example, Fig. 7(a) shows that at 200 GPa, the low-energy structures are dominated by structures with octahedral coordination, with more Mg than Fe, and approximately 50% O, i.e., ferropericlase or B1-type (Mg 1−x Fe x )O. At 350 GPa, the oxygen-rich ternary phases FeMgO 3 , Fe 2 MgO 4 , FeMg 2 O 4 , and FeMg 3 O 4 emerge as ground states or low-enthalpy phases, besides the B1-type phase. One of them, I-42d FeMg 2 O 4 , has the same structure as I-42d Mg 2 SiO 4 , the stable silicate phase predicted to exist in the mantle of Super-Earths above 500 GPa [17,43]. Here emerges the possibility of an I-42dtype Mg 2 (Si 1−x Fe x )O 4 phase, with Fe substitutional in the Si site, or vice-versa, an unusual type of substitution in the Earth's mantle, unless as a coupled Mg-Si substitution. From the chemistry standpoint, the newly found phases at 350 GPa can all be viewed as combinations of binary end-members, e.g., FeMgO 3 as (MgO)(FeO 2 ), Fe 2 MgO 4 as (MgO)(Fe 2 O 3 ), FeMg 2 O 4 as (MgO) 2 (FeO 2 ), and FeMg 3 O 4 as (MgO) 3 (FeO). Such stable compositions suggest other stable stoichiometric phases might be found by exploring combinations of such end-member compounds, as seen in the Mg-Si-O system, i.e., (MgO) n (SiO 2 ) m phases [17,43]. Further AGA searches aiming at these complex compositions are needed to identify other possible ternary phases in the Mg-Fe-O system. Finally, O's greater intermixing with the metallic elements at 350 GPa suggests that Mg and O abundances might be non-negligible in the Earth's inner core. Also, core formation by Fe exsolution from the oxides might be a more complicated process during Super-Earths' core formation, or O and Mg might be more abundant light elements in Super-Earths' cores.

V. CONCLUSION
We use the AGA combined with ab initio calculations to identify high-pressure structures in the Fe-Mg-O system at 0 K across a wide range of stoichiometries. This procedure is a crucial preparatory stage for modeling the system at finite temperatures. At 350 GPa, we identify mechanically stable phases with FeMg 2 O 4 , Fe 2 MgO 4 and FeMg 3 O 4 compositions and one low-enthalpy phase with FeMgO 3 composition. These discoveries lead to the construction of the ternary phase diagram and convex hull at 350 GPa. While we have not found any groundstate stoichiometric ternary compound at 200 GPa, the metastable phases' analysis indicates that ferropericlaseor magnesiumwüstite-type phases with 50% oxygen are very close to the convex hull. Oxygen-rich phases are generally closer to the convex hull than the oxygen-poor phases at all pressures. Motif analyses show octahedral clusters are energetically favored at both pressures and BCC-type clusters start to appear in oxygen-rich phases at 350 GPa. In particular, the nature of iron-rich phases at 350 GPa indicates that Mg, O, or both simultaneously could stabilize a BCC-type iron alloy at inner-core pressures.