Ab initio superionic-liquid phase diagram of Fe_1-xO_x under Earth’s inner core conditions

The key feature of the superionic Fe_1-xO_x alloy is the O diffusion in the crystalline lattice. Based on the ab initio MD simulation, we find O shows a similar mean squared displacement (MSD) in hcp, bcc, and liquid Fe phases at 323 GPa and 5500 K, a condition close to those at ICB (see Supplementary Information Fig. S3a). Thus, O is superionic in both hcp and bcc Fe under ICB conditions. Because ab initio MD simulation is highly limited by the time and length scales, we first use classical MD simulations to study the stability of superionic phases coexisting with liquid. Davies et al. [39] previously developed a Fe-O interatomic potential for IC conditions. While this potential can simulate the nucleation of Fe-O system [39], we find it failed to correctly describe the superionic state in Fe_1-xO_x. Here we developed a new Fe-O interatomic potential using the embedded-atom method [40](EAM) that accurately describes the superionic state in hcp and bcc lattices under ICB conditions (Supplementary Information Note 1). The MSD results from classical MD align qualitatively with ab initio simulations (Supplementary Information Fig. S3b).

To simulate the superionic hcp-liquid coexistence we constructed an hcp-liquid interface with 12,288 Fe atoms and randomly distributed O atoms for various O compositions. After a 2 ns MD simulation, O composition decreased in the hcp phase (Fig. 1a). The system reached equilibrium at around 100 ps, as indicated by the energy change in Fig. 1b. Extending the simulation to 2 ns showed that O atoms diffused throughout the simulation cell, providing sufficient data to compute their partitioning between hcp and liquid phases (Fig. 1c). The averaged O distribution suggests the O composition in the superionic phase ($x_{C}^{SI}$) is significantly lower than that in the liquid phase ($x_{C}^{L}$). Similar simulations across temperatures from 4600-6000 K at the same pressure yielded temperature-dependent $x_{C}^{SI}$ and $x_{C}^{L}$, representing the superionic hcp solidus and liquidus lines (Fig. 1d).

Based on the thermodynamic relations, when $Fe_{1-x}O_{x}$ superionic solution coexists in equilibrium with the liquid solution at the temperature of $T_{0}$, it satisfies the equilibrium condition,

where $G_{C}^{L}(x_{C_{0}}^{L},T_{0})$ and $G_{C}^{SI}(x_{C_{0}}^{SI},T_{0})$ are the absolute Gibbs free energies of the liquid and superionic phases with O compositions of $x_{C_{0}}^{L}$ and $x_{C_{0}}^{SI}$ at the temperature $T_{0}$, respectively (please refer to Supplementary Information Note 2 for detailed derivations). Since the MD simulation in Fig. 1 provides $x_{C_{0}}^{L}$ and $x_{C_{0}}^{SI}$ at various temperatures $T_{0}$, Eqn. (1) can be employed to compute the Gibbs free energy of the superionic phase $G_{C}^{SI}(x_{C_{0}}^{SI},T_{0})$, provided the liquid’s Gibbs free energy $G_{C}^{L}(x,T_{0})$ is known (see Supplementary Information Note 3).

A series of free energy calculations for a liquid solution is performed across the O composition range of 0-20 at.$\%$ with the compositions spaced equally, at various temperatures, as shown in Fig. 2a. We find these liquid’s free energy data can be well described by the regular solution model using the Redlich-Kister (RK) expansion [41],

where $G^{\text{Fe}}(T_{0})$ is the Gibbs free energy of pure Fe, $a$ and $L_{k}$ are fitting parameters. Only two RK terms ($k=0$ and $1$) are required to fit the liquid’s free energy data, achieving fitting errors of less than 0.1 meV/atom at all studied temperatures.

Based on Eqn. (1) and $G_{C}^{L}(x,T_{0})$, we can compute $G_{C}^{SI}(x_{C_{0}}^{SI},T_{0})$ for each $(x_{C_{0}}^{SI},x_{C_{0}}^{L},T_{0})$ combination obtained from MD simulations of superionic-liquid coexistence shown in Fig. 1. This results in sparse Gibbs free energy data for the superionic hcp state at a few temperatures, marked as solid circles in Fig. 2b. We then extend these data to a broader temperature range using the Gibbs-Helmholtz equation.

Here, $H(x_{0},T)$ is the temperature-dependent enthalpy with a specific O composition $x_{0}$ obtained from MD simulations. Using this approach, $G_{C}^{SI}(x_{C_{0}}^{SI},T)$ values are computed from $G_{C}^{SI}(x_{C_{0}}^{SI},T_{0})$ and $H(x_{C_{0}}^{SI},T)$ for different $x_{C_{0}}^{SI}$, and replotted as a function of compositions in Fig. 2b. We find that the free energy data of the superionic state can also be well-fitted by the RK expansion with only one RK term ($k=0$) to achieve a fitting error of less than 0.2 meV/atom.

With the absolute Gibbs free energy for both liquid and superionic solutions across various compositions and temperatures, the common tangent line approach can now provide the solidus and liquidus curves. We plot the relative Gibbs free energy using the 0$\%$ and 20$\%$ liquid free energy data as references for better visualization. Figure 2d shows the common tangent lines computed between $G_{C}^{\text{SI-hcp}}(x)$ and $G_{C}^{L}(x)$ curves at 5500 K, which resulted in intersections at $x_{C_{0}}^{\text{SI-hcp}}=0.88\%$ and $x_{C_{0}}^{L}=7.73\%$. These values are consistent with the equilibrium compositions of $x_{C_{0}}^{\text{SI-hcp}}=0.87\%$ and $x_{C_{0}}^{L}=7.70\%$ obtained from superionic hcp-liquid coexistence simulations under the same pressure and temperature conditions shown in Fig. 1c. More free energy data and their common tangent lines at other temperatures are shown in Supplementary Information Fig. S4. The superionic hcp solidus and liquidus curves computed from the free energy calculations are compared with those from MD simulations in Fig. 1d. Both methods result in a consistent superionic-liquid phase diagram, validating each other.

Because Eqns. (1-3) work for both superionic hcp and bcc phases, we repeated the calculations for the superionic bcc structure. The temperature-dependent free energies of superionic bcc are shown in Fig. 2c. The superionic bcc solidus and liquidus curves are shown in Supplementary Information Fig. S5. Figure 2d compares the relative free energy between superionic bcc and superionic hcp at 5500 K. More comparisons at different temperatures are shown in Fig. S4. These data suggest that the superionic bcc is metastable compared to the superionic hcp phase when the O composition is small. When the O composition is greater than $\sim$3 at.$\%$, the superionic bcc becomes more stable than superionic hcp in Fe_1-xO_x. However, the common tangent lines suggest when equilibrated with the liquid solution, the composition in the superionic phase is less than $\sim$1 at.$\%$. Thus, based on the classical potential, only superionic hcp can coexist with liquid at ICB.

We have obtained the absolute Gibbs free energy, $G_{C}$, of the liquid and superionic solutions for the classical system. Using the classical system as the reference state, we can perform TI to compute the Gibbs free energy of liquid and superionic phases for the ab initio system (noted as A system). The formulas for such TI calculation are derived in Supplementary Information Note 4. A series of TI simulations from the classical system to the ab initio system were performed for liquid, superionic hcp, and superionic bcc solutions at 5500 K, as shown in Fig. 3a. The TI path from C to A system is smooth and almost linear for all three phases, suggesting great similarity between the classical and ab initio systems [21]. As the ab initio MD simulations were performed with PAW8 potential ($3d$⁷4$s$¹ valence electrons), we further performed the FEP correction from PAW8 to PAW16 (3$s$²3$p$⁶3$d$⁷4$s$¹ valence electrons), as shown in Supplementary Information Fig. S8, which was demonstrated to be necessary for an accurate correction of the free energy data of Fe under Earth’s core conditions [21, 42]

Figure 3b shows the relative ab initio Gibbs free energy for liquid and superionic phases at 5500 K and 323 GPa. The regular solution model with the RK expansion of Eqn. (2) can also describe the compositional dependence of these free energies, providing fitting errors of less than 0.5 meV/atom. This error is larger than those in classical systems, mainly due to fewer data points and smaller length scales in ab initio calculations. Nevertheless, such a free energy error is sufficiently small for phase diagram calculations. In Fig. 3b, the ab initio free energy of the superionic bcc phase is higher than the superionic hcp phase when the O composition is small. This is consistent with the fact that pure Fe prefers the hcp phase under IC conditions [21, 22]. As the O composition increases, the free energy of the superionic hcp phase quickly increases. When the O composition is higher than $\sim$3 at.%, superionic bcc becomes more stable than superionic hcp. This trend is similar to the one found in the classical simulation. Thus, the O composition in Fe_1-xO_x changes the relative stability between the superionic bcc and hcp phases.

However, the common tangent line between liquid and superionic phases suggests that the superionic solution can only coexist with the liquid solution at small O compositions. Based on the $G_{A}^{\text{SI-hcp}}(x)$ and $G_{A}^{L}(x)$ curves, the common tangent line reveals ab initio solidus and liquidus points at 5500 K of $x_{A_{0}}^{\text{SI-hcp}}=0.47\%$ and $x_{A_{0}}^{L}=14.8\%$, respectively. Superionic hcp is more stable than superionic bcc at such a small O composition. We further employ the Gibbs-Helmholtz equation, i.e., Eqn. (3), to extend the free energy data to other temperatures for all three phases, as shown in Supplementary Information Fig. S9. The superionic hcp solidus and liquidus curves are computed using the common tangent approach. These data provide the ab initio phase diagram of superionic hcp and liquid Fe_1-xO_x at 323 GPa in Fig. 3c. The superionic solidus line shows small temperature dependence, while the liquidus line depends on temperature more strongly. This results in a partition coefficient $\left(D^{\text{SI/L}}=\frac{x_{A}^{\text{SI}}}{x_{A}^{L}}\right)$ that is strongly temperature-dependent.

The phase diagram in Fig. 3b indicates that the equilibrium O compositions in the solid IC and liquid outer core are highly correlated with the temperature at ICB. It provides a stronger constraint on the core’s composition and temperature than the one from the partition coefficient data alone. Based on the phase diagram, the densities of Fe_1-xO_x superionic and liquid solutions can be computed under the equilibrium conditions at different temperatures, as shown in Fig. 3d. In previous work [43], O was used to account for the large density difference between the solid and liquid core, i.e., the ”density jump”. By assuming O is the only light element and combining the partition coefficient with the density jump data, the O compositions in the solid core and liquid core were proposed to be $0.2\pm 0.1$ at.% and $8.0\pm 2.5$ at.% [43], respectively. We use the current density of superionic hcp and liquid Fe_1-xO_x solutions to match the density jump in the Preliminary Reference Earth Model [44]. It simultaneously constrains the O compositions as $0.35\pm 0.05$ at.% in the solid core and $9.7\pm 1.6$ at.% in the liquid core, and a temperature of $5790\pm 90$ K at ICB. Therefore, compared to the previous substitutional solid solution, the superionic phase nearly doubles O’s solubility in hcp Fe under Earth’s IC conditions. While the superionic bcc phase is metastable at this small O composition, it cannot be entirely excluded as we have seen that the bcc phase can be stabilized in Fe-Ni alloys [45]. If Ni is included, the Gibbs free energy of the superionic bcc phase is likely to be lowered in Fig. 3b. Thus, other elements’ effects must be included to fully resolve the core’s composition and structure.