Structural transition and re-emergence of iron’s total electron spin in (Mg,Fe)O at ultrahigh pressure

Within LDA+$U_{sc}$, both the IS and LS states of B2 (Mg_1-xFe_x)O can be obtained at ultrahigh pressure, while the HS state can only be obtained at $P\lesssim 0.29$ TPa. The Hubbard $U_{sc}$ of Fe in (Mg_0.875Fe_0.125)O at various volume/pressure are shown in Fig. 2. Note that B2 (Mg_0.875Fe_0.125)O is stabilized via rhombohedral distortion (as further discussed in Figs. 3 and 4), hence referred to as rB2 hereafter. At ultrahigh pressure ($0.5<P<1.8$ TPa), as shown in Fig. 2b, $U_{sc}$ is mainly affected by pressure (increasing with $P$ by $\sim$2.5 eV) and marginally affected by the Fe spin state and crystal structure (by $\sim$0.5 eV). In contrast, at $P<0.15$ TPa, Fe spin/valence state affects $U_{sc}$ by up to $\sim$2 eV, while pressure affects $U_{sc}$ by up to $\sim$0.5 eV [22, 49, 50, 52, 51, 53, 54].

Fe spin states in B1 and B2 (Mg_1-xFe_x)O exhibit distinct properties, due to the host minerals’ distinct crystal structures. In the B1 phase, Fe substitutes Mg in the 6-coordinate octahedral site (Fig. 1a, d), forming FeO₆ octahedra (Supplementary Fig. 1a and Note 1). In FeO₆ octahedra, the $t_{2g}$ orbitals have lower energy than the $e_{g}$ orbitals; the orbital configurations of HS and LS Fe²⁺ are $t_{2g}^{4}e_{g}^{2}$ and $t_{2g}^{6}e_{g}^{0}$, respectively (Supplementary Fig. 1b and Note 1). In the B2 phase, Fe substitutes Mg in the 8-coordinate site (Fig. 1b, e), forming FeO₈ polyhedra (Fig. 1c, f). For cubic B2 (Mg_1-xFe_x)O with $x=0.125$, after structural optimization, FeO₈ polyhedra remain cubic ($O_{h}$ symmetry), with all eight Fe-O bonds in the same length. Expectedly, the Fe-O bonds of the IS state are slightly longer than those of the LS state (Fig. 3a, f). Remarkably, cubic B2 (Mg_0.875Fe_0.125)O is dynamically unstable, regardless of the Fe spin state, as indicated by the soft phonon modes (negative phonon frequencies) shown in Fig. 3b, g. Nevertheless, the electronic structure of cubic B2 (Mg_0.875Fe_0.125)O still provides valuable insights. In FeO₈ cubes, the $t_{2g}$ orbitals have higher energy than the $e_{g}$ orbitals, and the orbital configurations of IS and LS Fe²⁺ are both $e_{g}^{4}t_{2g}^{2}$ (Fig. 3c, h). For the IS state, the $e_{g}$ orbitals are fully occupied, while the $t_{2g}$ orbitals are partially occupied by two spin-up electrons (Fig. 3c). A partially filled $t_{2g}$ band in the spin-up channel is thus formed, spanning across the Fermi energy (set as the reference, $0$ eV) from $-1.7$ to $1.3$ eV, as indicated by the density of states (DOS) shown in Fig. 3e. The $t_{2g}$ characteristic of the $t_{2g}$ band can be visualized via the integrated local density of states (ILDOS) over the energy intervals $-1.7<E<0$ and $0<E<1.3$ eV for the filled and empty $t_{2g}$ states, respectively (Fig. 3d). For the LS state, the $t_{2g}$ orbitals are partially occupied by one spin-up and one spin-down electron (Fig. 3h), forming a partially filled $t_{2g}$ band in the interval of $-1.2<E<1.6$ eV (Fig. 3j). The ILDOS of the filled and empty $t_{2g}$ states (Fig. 3i) resemble those of the IS state (Fig. 3d). For both spin states, the completely filled $e_{g}$ bands are embedded in the oxygen band spanning over $-20\lesssim E\lesssim-4$ eV (Fig. 3e, j). Evidently, the partially filled $t_{2g}$ band and thus the metallicity of cubic B2 (Mg_0.875Fe_0.125)O are the direct consequences of the cubic symmetry.

Depending on the Fe spin state, dynamically unstable cubic B2 (Mg_0.875Fe_0.125)O is stabilized via rhombohedral compression or elongation. As shown in Fig. 4a/f, the IS/LS state is rhombohedrally compressed/elongated, with shortened/stretched Fe-O bonds along the $[111]$ direction and rhombohedral angle $\alpha$ larger/smaller than $90^{\circ}$. The resultant rB2 structures for both spin states are dynamically stable with no soft phonon mode (Fig. 4b, g). With rhombohedral distortion, the FeO₈ cubes become FeO₈ dodecahedra ($D_{3d}$ symmetry), and the three $t_{2g}$ orbitals split into a singlet ($a_{1g}$), which is a $d_{z^{2}}$-like orbital along the [111] direction, and a doublet ($e_{g}^{\prime}$). For the IS state, with shortened Fe-O bonds along the $[111]$ direction, the $a_{1g}$ orbital has higher energy than the $e_{g}^{\prime}$ orbitals; the four spin-up electrons occupy the $e_{g}$ and $e_{g}^{\prime}$ orbitals, and the two spin-down electrons occupy the $e_{g}$ orbitals (Fig. 4c). With the splitting of $t_{2g}$ orbitals, an energy gap ($\sim$0.3 eV) is opened between the $e_{g}^{\prime}$ and $a_{1g}$ bands, as indicated by the DOS (Fig. 4e). The $e_{g}^{\prime}$ and $a_{1g}$ characteristics of these bands can be visualized via the ILDOS (Fig. 4d). For the LS state, with stretched Fe-O bonds along the $[111]$ direction, the $a_{1g}$ orbital has lower energy than the $e_{g}^{\prime}$ orbitals (Fig. 4h); the three spin-up and three spin-down electrons fully occupy the $e_{g}$ and $a_{1g}$ orbitals, with the $e_{g}^{\prime}$ orbitals left unoccupied, resulting in a gap ($\sim$2.0 eV) between the $a_{1g}$ and $e_{g}^{\prime}$ bands (Fig. 4j). The $a_{1g}$ and $e_{g}^{\prime}$ characteristics of these two bands can also be visualized via the ILDOS (Fig. 4i). For both spin states, the completely filled $e_{g}$ bands are embedded in the oxygen bands spanning over $-20\lesssim E\lesssim-4$ eV (Fig. 4e, j).

When the Fe concentration increases to $x=0.25$, a three-dimensional (3D) network of corner-sharing FeO₈ cubes is formed in B2 (Mg_0.75Fe_0.25)O (Fig. 1f), while for $x\leq 0.125$, the FeO₈ polyhedra are isolated/unconnected (Fig. 1c). For isolated FeO₈ polyhedra, rhombohedral distortion is allowed and favored, as observed in rB2 (Mg_0.875Fe_0.125)O (Fig. 4). In contrast, connectivity of the 3D FeO₈ network in B2 (Mg_0.75Fe_0.25)O suppresses the rhombohedral distortion and further stabilizes the cubic structure: Starting the structural optimization with rhombohedrally compressed/elongated rB2 IS/LS (Mg_0.75Fe_0.25)O, the crystal structure and FeO₈ polyhedra resume cubic symmetry within a few steps. Within LDA+$U_{sc}$, the LS state of cubic B2 (Mg_0.75Fe_0.25)O can be obtained throughout 0.2–1.8 TPa, while the IS state can only be obtained at $P<1.1$ TPa and is subject to magnetic collapse: The total magnetization ($M$) decreases from $2\mu_{B}$/Fe to $0$ in the region of $0.6<P<1.1$ TPa and vanishes at $P>1.1$ TPa (Fig. 5a). (Note: For $x=0.125$, the IS state retains $M=2\mu_{B}$/Fe up to $1.8$ TPa.) Regardless of the spin state and magnetization, cubic B2 (Mg_0.75Fe_0.25)O is dynamically stable with no soft phonon mode, even during the magnetic collapse (Fig. 5b–d). With FeO₈ cubes ($O_{h}$ symmetry), cubic B2 (Mg_0.75Fe_0.25)O and (Mg_0.875Fe_0.125)O have the same $3d$ orbital occupations and similar electronic structures. For cubic B2 IS (Mg_0.75Fe_0.25)O with $M=2\mu_{B}$/Fe (before the magnetic collapse), the spin-up $t_{2g}$ band spans across the Fermi energy ($0$ eV) while the spin-down $t_{2g}$ band lies above the Fermi energy (Fig. 5e), showing the same characteristic as cubic B2 IS (Mg_0.875Fe_0.125)O (Fig. 3e). Likewise, for cubic B2 LS (Mg_0.75Fe_0.25)O (Fig. 5g) and (Mg_0.875Fe_0.125)O (Fig. 3j), the $t_{2g}$ bands in both spin channels align, spanning across the Fermi energy. During the magnetic collapse ($0<M<$ 2$\mu_{B}$/Fe), the $t_{2g}$ bands in the spin-up and spin-down channels are shifted upward and downward, respectively (Fig. 5f).

To analyze the structural and spin transition of (Mg_1-xFe_x)O, we compute the equations of state (EoS) of the B1 LS, (r)B2 IS, and (r)B2 LS states (Supplementary Fig. 2 and Note 2); the relative enthalpies ($\Delta H$) of these states with respect to the (r)B2 IS state are plotted in Fig. 6. For $x=0.125$ (Fig. 6a), (Mg_0.875Fe_0.125)O transforms from the B1 LS state into the rB2 IS state at 0.642 TPa. Upon further compression, the rB2 IS state undergoes a spin transition to the rB2 LS state at 1.348 TPa. Throughout these transitions, (Mg_0.875Fe_0.125)O remains insulating (see Supplementary Fig. 1 and Note 1 for the insulating B1 LS state). Remarkably, in the simultaneous structural (B1–rB2) and spin (LS–IS) transition at 0.642 TPa, Fe’s total electron spin $S$ re-emerges from 0 to 1, opposite to the perception that $S$ decreases upon compression. Furthermore, the IS state is energetically favorable over a wide pressure range (0.612–1.348 TPa), despite that IS state has never been observed in the Earth’s lower-mantle minerals, including B1 (Mg_1-xFe_x)O, Fe-bearing MgSiO₃ bridgmanite and post-perovskite, ferromagnesite (Mg_1-xFe_x)CO₃, and the NAL phase [22, 49, 50, 52, 51, 53, 54]. For $x=0.25$ (Fig. 6b), a simultaneous structural, spin, and metal–insulator transition occurs at $0.539$ TPa, from the insulating B1 LS state to the metallic B2 IS state (notice that $S$ increases). Upon further compression, an IS–LS transition occurs in metallic B2 (Mg_0.75Fe_0.25)O at 0.855 TPa. From Fig. 6, effects of Fe concentration on the B1–(r)B2 transition pressure ($P^{B1/B2}_{T}$) can also be inferred. For Fe-free MgO ($x=0$), we find $P^{B1/B2}_{T}=0.535$ TPa (Supplementary Fig. 3 and Note 3), in agreement with other calculations [3, 4, 7, 8, 9, 10, 11, 12, 13, 14, 6, 5]. As $x$ increases, $P^{B1/B2}_{T}$ first increases to $0.642$ TPa at $x=0.125$ (Fig. 6a) and then decreases to $0.539$ TPa at $x=0.25$ (Fig. 6b), indicating a trend of decreasing $P^{B1/B2}_{T}$ in the region of $0.125\lesssim x\leq 1$, consistent with experiments: $P^{B1/B2}_{T}$ of FeO ($\sim$0.25 TPa) [47, 48] is much lower than that of MgO ($\sim$0.5 TPa) [1, 2, 3, 4]. [Note: (1) For $x=0.25$ (Fig. 6b), the EoS of the B2 IS state is fitted using the data points at $P<0.988$ TPa (where $M>0$). (2) To examine the robustness of the LDA+$U_{sc}$ results shown in Fig. 6, we perform extensive test calculations using various methods. Similar results are obtained, as shown in Supplementary Figs. 4, 5, and Note 4].

In the Earth’s mantle condition, Fe partitioning between B1 (Mg_1-xFe_x)O, Fe-bearing MgSiO₃ bridgmanite, and post-perovskite varies with pressure, temperature, and even the Fe valence/spin state [23, 21, 37, 38, 39]. Likewise, in exoplanet interiors, Fe concentration in B2 (Mg_1-xFe_x)O may vary with the depth or the interior region, due to the variation of Fe partitioning between B2 (Mg_1-xFe_x)O and other minerals phases, including post-perovskite and/or high-pressure silicates [15, 16]. Evident by comparing Fig. 6a and b, spin, structural, and metal–insulator transition of B2 (Mg_1-xFe_x)O can also be induced by the change of Fe concentration ($x$). In the depth/region with pressure of 0.642–0.855 TPa, if $x$ increases from 0.125 to 0.25, a simultaneous structural and metal–insulator transition occurs [insulating rB2 IS (Mg_0.875Fe_0.125)O $\rightarrow$ metallic B2 IS (Mg_0.75Fe_0.25)O]; in the depth/region with pressure of 0.855–1.348 TPa, if $x$ increases from 0.125 to 0.25, a simultaneous structural, spin, and metal–insulator transition occurs [insulating rB2 IS (Mg_0.875Fe_0.125)O $\rightarrow$ metallic B2 LS (Mg_0.75Fe_0.25)O]. Even in the depth/region of $P>1.348$ TPa, where only LS Fe²⁺ exists, if $x$ increases from $0.125$ to $0.25$, a simultaneous structural and metal–insulator transition occurs [insulating rB2 LS (Mg_0.875Fe_0.125)O $\rightarrow$ metallic B2 LS (Mg_0.75Fe_0.25)O]. On the other hand, if $x$ decreases from $0.25$ to $0.125$, the aforementioned transitions would be reversed. Based on the above analysis, metal–insulator transition is always included in the composition-induced transitions of the B2 phase, suggesting that variation of Fe partitioning can significantly change the electrical and thermal transport properties of exoplanet interiors.

At temperature $T\neq 0$, spin transition goes through a mixed-spin (MS) phase/state, in which different spin states coexist. For B2 (Mg_1-xFe_x)O, only the IS and LS states are relevant. Within the thermodynamic model detailed in Supplementary Note 5, the LS fraction ($n_{LS}$) in the MS phase can be written $n_{LS}=1/\left[1+3\exp(\Delta H/k_{B}Tx)\right]$, where $\Delta H\equiv H_{LS}-H_{IS}$, and the IS fraction $n_{IS}=1-n_{LS}$. Despite that lattice vibration is not considered, the results obtained from this approach have been shown in great agreement with room-temperature experiments [50, 22, 53, 54, 56]. In Fig. 7, the LS and IS fractions, compression curves, and bulk modulus of B2 (Mg_1-xFe_x)O at $T=300$ K are shown. For $x=0.125$, the IS–LS transition is smooth and spans over a wide pressure range (Fig. 7a), due to the small enthalpy difference ($\Delta H$) between the rB2 IS and LS states (Fig. 6a). The compression curves $V(P)$ of the MS, IS, and LS states are nearly the same (Fig. 7b); their difference is barely noticeable even by plotting the relative volume difference with respect to pure B2 MgO, namely, $(V-V_{\text{MgO}})/V_{\text{MgO}}$ (Fig. 7c). As a consequence, the bulk modulus $K\equiv-V\partial P/\partial V$ barely changes during the spin transition (Fig. 7d). For $x=0.25$, the spin transition is more abrupt (Fig. 7e), due to the larger $\Delta H$ between the B2 IS and LS states (Fig. 6b). With increased Fe concentration, the volume difference between the LS and IS states increases (Fig. 7f), resulting in prominent volume reduction ($\sim$0.5%) (Fig. 7g) and bulk modulus softening ($\sim$22%) (Fig. 7h). While the full elastic tensor ($C_{ij}$) and shear modulus ($G$) are not computed, the volume and elastic anomalies shown in Fig. 7g and h clearly indicate anomalous softening of bulk sound velocity $v_{\Phi}=\sqrt{K/\rho}$ ($\rho$: density) and compressional wave velocity $v_{P}=\sqrt{\left(K+\frac{4}{3}G\right)/\rho}$. Furthermore, within the phonon gas model, lattice thermal conductivity $\kappa\approx(1/3)C_{V}v_{P}^{2}\tau$ ($C_{V}$: heat capacity; $\tau$: average phonon scattering time) [57], suggesting that anomalous change of $v_{P}$ may play a role in the anomalous change of thermal conductivity. In the B1 phase, anomalous $v_{\Phi}$, $v_{P}$ [20, 28, 29, 30, 31, 32, 33] and $\kappa$ [35, 36] have all been observed. The volume/elastic anomalies accompanying the spin transition of the B2 phase may thus be a possible source of seismic and thermal anomalies in exoplanet interiors.

Experiments and computations have confirmed that HS–LS transitions in B1 (Mg_1-xFe_x)O and ferromagnesite (Mg_1-xFe_x)CO₃ are accompanied by anomalous changes of major thermal properties, including thermal expansivity, heat capacity, Grüneisen parameter, thermal conductivity, and thermoelasticity [20, 28, 29, 30, 31, 32, 33, 35, 36, 55]. As the temperature increases, the spin-transition pressure increases, and the transition is broadened [29, 55]. Remarkably, for (Mg_0.75Fe_0.25)CO₃, the anomalous increase of heat capacity retains its magnitude ($\sim$12%) without smearing out, even at high temperature [55]. Likewise, for the B2 phase, anomalous changes of thermal properties accompanying the spin transition can be expected. To investigate such anomalies at high $P$–$T$ conditions, vibrational free energy must be included. Thermal calculations are thus highly desirable and will be left for future studies.

In summary, we have investigated (Mg_1-xFe_x)O ($x\leq 0.25$) at ultrahigh pressure up to $1.8$ TPa via first-principles calculations. Our calculations indicate that Fe greatly affects the properties of (Mg_1-xFe_x)O. For $x=0.125$, insulating (Mg_0.875Fe_0.125)O undergoes a simultaneous structural and spin transition (B1 LS $\rightarrow$ rB2 IS) at $0.642$ TPa, followed by a spin transition (rB2 IS–LS) at $1.348$ TPa. For $x=0.25$, (Mg_0.75Fe_0.25)O undergoes a simultaneous structural, spin, and metal–insulator transition (insulating B1 LS $\rightarrow$ metallic B2 IS) at 0.539 TPa, followed by a spin transition (metallic B2 IS–LS) at 0.855 TPa. Remarkably, Fe’s total electron spin $S$ re-emerges from 0 to 1 in the B1–(r)B2 transition. In addition, structural, spin, and metal–insulator transitions of B2 (Mg_1-xFe_x)O can also be induced by the change of Fe concentration ($x$). These results suggest that Fe and spin transition may greatly affect planetary interiors over a vast pressure range, considering the anomalous changes of elastic, transport, and thermal properties accompanying the spin and/or metal–insulator transition.

In this work, all major calculations are performed using the Quantum ESPRESSO codes [58]. We use ultrasoft pseudopotentials (USPPs) generated with the Vanderbilt method [59]. The valence electron configurations for the generations are $2s^{2}2p^{6}3s^{2}3p^{0}3d^{0}$, $3s^{2}3p^{6}3d^{6.5}4s^{1}4p^{0}$, and $2s^{2}2p^{4}3d^{0}$ for Mg, Fe, and O, respectively; the cutoff radii are 1.4, 1,8, and 1.0 a.u. for Mg, Fe, and O, respectively. The aforementioned USPPs of Mg and O have been used in Refs. 15, 16, and the USPP of Fe has been used in Refs. 22, 49, 50, 52, 51, 53, 54. To properly treat the on-site Coulomb interaction of the Fe $3d$ electrons, we adopt the local density approximation + self-consistent Hubbard $U$ (LDA+$U_{sc}$) method, with the Hubbard $U$ parameters computed self-consistently [60, 61, 62, 63]. Briefly speaking, we start with an LDA+$U$ calculation with a trial $U$ (the “input $U_{in}$”) to obtain the desired spin state for (Mg_1-xFe_x)O. For this LDA+$U_{in}$ state, we compute the second derivative of the LDA energy with respect to the $3d$ electron occupation at the Fe site ($d^{2}E_{LDA}/dn^{2}$) via a density functional perturbation theory (DFPT) approach [63] implemented in Quantum ESPRESSO. This second derivative, $d^{2}E_{LDA}/dn^{2}$, is considered as the “output $U_{out}$” and will be used as $U_{in}$ in the next iteration. Such a procedure is repeated until self-consistency is achieved, namely, $U_{in}=U_{out}\equiv U_{sc}$. Phonon calculations are performed using the Phonopy code [64], which adopts the finite-displacement (frozen phonon) method. The third-order Birch–Murnaghan equation of state (3rd BM EoS) is used for the EoS fitting.

In this work, analysis for the IS–LS transition in (r)B2 (Mg_1-xFe_x)O at room temperature ($300$ K) is based on the thermodynamic model detailed in Supplementary Note 5. Plotted in Fig. 7, the IS/LS fractions and the EoS of the MS phase are given by Supplementary Eqs. 9 and 10, respectively.