PBE-GGA predicts the B8↔B2 phase boundary of FeO at Earth’s core conditions

Significance FeO is a compound of great interest in materials physics and geophysics. It has a complex pressure–temperature (P,T) phase diagram, including subtle structural and magnetic transitions. Currently, it is challenging for prevailing theoretical approaches for strongly correlated materials, such as FeO at low P,Ts, to address anharmonic vibrations. We use density functional theory–based molecular dynamics and lattice dynamics calculations to demonstrate the dynamic and thermodynamic stability of an experimentally observed but theoretically unidentified B2-type phase of FeO at high P,Ts. We show that this phase is unstable at low Ts, but is stabilized at high P,Ts by anharmonic effects and thermal electronic excitations. This study establishes a theoretical framework for studying FeO and related Fe alloys at the Earth’s core conditions.

FeO is one of the major constituents of Earth and other terrestrial planets.It is not only the iron end-member of ferropericlase ((Mg 1−x Fe x )O), the second most abundant phase in the Earth's lower mantle, but also a significant alloying component in the Earth's core [1].As a classic "correlated" oxide, phase relations in FeO are also of great interest in condensed matter physics.Like ferropericlase, FeO undergoes a spin state change under pressure [2,3].FeO exhibits rich phenomenology at high pressures and temperatures (P, T ), e.g., polymorphic, magnetic, and insulating-metallic transitions.Such phase changes control fundamental material properties of the Earth, e.g., its thermal and electrical conductivity and magnetic susceptibility, to mention a few.
Here we perform ab initio calculations of the B8↔B2 boundary at ∼150-400 GPa, a relevant pressure range for the Earth's core, and T > 1000 K.Under such extreme conditions, anharmonicity is fundamental in determining dynamic and thermodynamic stabilities, especially for the B2 phase (e.g., the analogous bcc phase of elemental metals [25][26][27][28][29] and the B2 phase of binary compounds [30]).The high-spin (HS) to low-spin (LS) and insulator-to-metal transitions in FeO happen at ∼120 GPa [24,31], a typical mantle pressure.A recent density functional plus dynamical mean-field theory (DFT+DMFT) study [20] confirmed the LS and metallic (delocalized) electronic state in both B2 and B8 phases under core conditions.These observations suggest that DFT-based molecular dynamics (DFT-MD) could address the B8↔B2 phase competition.
We first address harmonic phonon dispersions at 0 K in both B2 and B8 phases using density-functional perturbation theory (DFPT) [32].DFT and DFPT calculations were performed using the PAW method [33] as implemented in Quantum ESPRESSO [34].We used the Perdew-Burke-Ernzerhof (PBE) [35] generalized gradient approximation (GGA) to compute the exchangecorrelation energy.Phonons were calculated on 4 × 4 × 4 and 4 × 4 × 2 q-meshes in the B2 and B8 structures, re-spectively.Figures.1(a) and 1(b) show harmonic phonon dispersions (dashed gray curves) for these phases.The B2 phase displays unstable modes around phonon wave vector q = X 0, 1  2 , 0 , while the B8 phase's phonons are all stable.The harmonic phonon instability in the B2 phase is not out of expectation.A similar soft mode behavior in the bcc phase of many elemental metals [25][26][27][28][29] drives the bcc to hcp phase transition at low temperatures.The B2 and B8 phases of FeO are analogous to the elemental bcc and hcp phases, respectively (see Fig. S1 of the Supplemental Materials (SM) [36]).Note that though B8 FeO shares the same P 6 3 /mmc space group with the elemental hcp phase, the atomic layers stacking in B8 ( differs from that in elemental hcp (AB AB).Imaginary frequencies at 0 K indicate that phonon-phonon interactions are critical in the dynamic stabilization of this structure, and entropic effects at high temperatures can stabilize it thermodynamically.Conventional harmonic phonons or quasiharmonic free energies cannot address the lattice dynamics and thermodynamic properties of this B2 phase.
We employ the phonon quasiparticle (PHQ) approach [37,38] to address the high-temperature anharmonicity of FeO under core P, T conditions.This approach assumes that a system with fully interacting phonons can be simplified as an effective system with noninteracting PHQs [39,40].Each PHQ can be described by two parameters, renormalized frequency, ω qs , and linewidth, Γ qs .A PHQ is numerically defined by mode-projected velocity autocorrelation function (VAF) [37,38], where • êi qs is the massweighted and qs-mode-projected velocity.M i , R i , and v i (i = 1, . . ., N ) are atomic mass, crystallographic atomic coordinate, and atomic velocity computed by DFT-MD of an N -atom supercell.êi qs is the harmonic phonon eigenvector obtained at an electronic temperature T el = T in the DFT-MD [29,39,40], where q should be commensurate with the supercell size and s labels the phonon branches at each q.For a well-defined PHQ, its VAF assumes an exponentially decaying cosine form A qs cos ( ω qs t) e −Γqst , where A qs is the initial oscillation amplitude.The well-defined VAF's power spectrum 2 assumes a Lorentzian line shape with a single peak at ω qs and a linewidth Γ qs .
To compute PHQs, ab initio molecular dynamics (AIMD) simulations were performed using the PAW PBE as implemented in VASP [41].The electronic temperature (T el ) was set the same as the ionic temperature (T ) using the Mermin functional [42,43].FeO was simulated with 128-atom supercells (4×4×4 for B2 and 4×4×2 for B8) with a Γ k-point sampling and a kinetic energy cutoff of 400 eV.The supercells are sufficiently large to converge the harmonic part of interatomic force constants.Thus they should also be sufficiently large to converge the anharmonic part since the anharmonic part of interatomic force constants has shorter ranges than the harmonic one [38,44].Simulations were conducted in the N V T ensemble on a series of volumes between 4.92 and 5.88 Å3 /atom and temperatures between 1000 and 4000 K controlled by Nosé thermostat [45,46].Each simulation ran for 50 ps, sufficiently long to converge PHQ parameters, with a time step of 1 fs.
Well-defined PHQs obtained on the AIMD-sampled q-mesh enable calculations of renormalized, i.e., anharmonic phonon dispersion using Fourier interpolation [37,38].Figures 1(a) and 1(b) show the anharmonic phonon dispersions obtained at T = 2000 K (solid orange curves) for the B2 and B8 phases, respectively.The unstable mode with imaginary harmonic frequency at q = X 0, 1  2 , 0 in the B2 phase stiffens drastically at high temperatures (see the gray square in Fig. 1(a) and Fig. S2 of the SM [36]).The anharmonic phonon dispersions are free of imaginary frequencies, indicating that B2 FeO is stabilized dynamically by phonon-phonon interactions at high temperatures [29,38,47].Frequency renormalization in the B8 phase is not as significant as in the B2, yet also not negligible.Thus, the following free energy calculations use temperature-dependent anharmonic phonon spectra (see Fig. S3 of the SM [36]) for both phases.The anharmonic phonon spectra were evaluated on much denser q-meshes (20 × 20 × 20 for B2 and 20 × 20 × 10 for B8) via Fourier interpolation [37,38] to approximate the thermodynamic limit.
We performed AIMD simulations for both B8 and B2 phases at the P, T conditions indicated in Fig. 2(a).The corresponding V, T conditions are shown in Fig. S4 of the SM [36].The B8 phase is dynamically stable at all studied P, T s, while the B2 phase shows structural and/or phononic instabilities at certain low P, T s.For instance, at T = 1000 K and P = ∼213 GPa (V = 5.88 Å3 /atom), the eight-fold coordinated Fe in the starting B2 structure transforms to a six-fold coordinated structure after thermal equilibration (see Fig. S5 of the SM [36]).By removing the B2 lattice constraints, a complete phase transition to B8 is realized.Compressed from ∼213 to ∼270 GPa (V = 5.52 Å3 /atom) at T = 1000 K, the structure no longer shows the B2→B8 transition.PHQs in the B2 phase are still not well-defined at this P, T .Fig. 2(b) shows the 1000 K VAF for the transverse acoustic mode at q = 1 4 , 1 4 , 1 4 (gray circle in Fig. 1(a)), which exhibits a pattern far distinct from an exponentially decaying cosine form.The corresponding power spectrum shows two peaks (Fig. 2(c)), indicating the breakdown of a well-defined B2 phase PHQ [29,38].This behavior signals the tendency of atoms to displace from the B2 equilibrium sites and the B2 structure to distort.In contrast, this mode is stable at 2000 K (∼275 GPa at the same V = 5.52 Å3 /atom), as indicated by the welldefined PHQ with an exponentially decaying VAF (Fig. 2(d)) and a single Lorentzian-shaped peak in the power spectrum (Fig. 2(e)).As seen in Fig. 2(a), higher P, T s systematically stabilize the B2 phase.
The extensive AIMD results indicated in Fig. 2(a) enable Gibbs free energy calculations on an equal footing for the B2 and B8 phases in a large P, T range.With Fourier interpolated anharmonic phonon spectra, the vibrational entropy can be obtained in the thermodynamic limit within the phonon gas model (PGM) [37,38,48], where . ω qs (T ) were obtained by fitting calculated ω qs 's at several temperatures and constant volume to a second-order polynomial in T [29,37,47].The Helmholtz free energy at constant volume can be obtained by integrating both the electronic and vibrational entropies [29,47], where T 0 is a reference temperature at which all PHQs are well-defined, E(T 0 ) and S el (T 0 ) are time-averaged internal (potential + kinetic) energy and electronic entropy obtained from AIMD at T 0 .The choice of T 0 does not change the resulting thermodynamics, so we set T 0 = 4000 K. S el (T ) at constant volume was computed as the ensemble average at temperatures shown in Figs.2(a) and S4 of the SM [36] and fit to a second-order polynomial in T [29].Contrary to S vib that suffers from finite-size effects and requires Fourier interpolation, S el and E are more insensitive to the simulation cell size [40], and AIMD ensemble averages converge well.
At the same V, T conditions, S vib is always larger for B2 than for B8 (Fig. 3(a)).The generally lower renormalized frequencies in B2 (see Fig. S3 of the SM [36]) contribute to its larger entropy and thermodynamic stability with respect to B8 at higher temperatures.S el varies nearly linearly with temperature (Fig. 3(b)), and a quadratic fitting for S el (T ) is sufficiently accurate [29].S el 's of both phases are similar but much smaller than S vib (see Figs. tions of state (EOS) were computed by fitting F (V ) to a third-order finite strain expansion.Pressures calculated as P = − ∂F ∂V T are shown in Fig. 3(d).At the same P, T conditions, the stable B2 phase volume is always smaller than the B8 volume, and the difference ∆(P V ) = (P V ) B8 − (P V ) B2 increases with pressure (see Fig. S6 of the SM [36]), which contribute to the enthalpic stabilization of B2 at high pressures.We thereby predict a negative Clapeyron slope ( dP dT = ∆S ∆V ) for the phase boundary.
The Gibbs free energy was calculated as G(P, T ) = F (V, T ) + P (V, T )V utilizing the fitted EOS.Comparing the Gibbs free energies of both phases (see Fig. S7 of the SM [36]), we obtain the B8↔B2 phase boundary shown in Fig. 4. The uncertainty in our free energy and phase boundary calculations are estimated as follows: Fig. S8 of the SM [36] shows the dependence of the Helmholtz free energy difference, ∆F , between the two phases on the k-mesh sampling.We consider the results on a much more computationally expensive 3×3×3 k-mesh in the 128-atom supercell to be fully converged since it differs by less than 1 meV/atom from the 2×2×2 k-mesh result.Hence, the difference between the Γ-point and 3 × 3 × 3 k-mesh results is ∼9 meV/atom.This is the adopted uncertainty in ∆F arising from k-mesh sampling.The uncertainty in ∆F arising from fluctuations in AIMD simulations is estimated by conducting five par-allel runs at a constant volume and 4000 K, which gives ∼0.1 meV/atom.Note that the dense q-mesh sampled in the entropy calculation (Eq.( 2)) mimics a sufficiently large supercell (16000 atoms) calculation in the thermodynamic limit.These combined effects give an uncertainty in ∆F of ∼9 meV/atom, and this value is passed to ∆G.This estimated uncertainty in the free energy difference compares well with the ∼10 meV/atom uncertainty reported in the calculation of the melting curve of iron [39] using the same PHQ approach.It is also similar to uncertainties in other free energy difference calculations using thermodynamic integration [49].The uncertainty in the P V term given by fitting the isothermal EOS at several volumes is a second-order effect, thus was disregarded here.The free energy uncertainty leads to an uncertainty of ±18 GPa in the transition pressure, shown as the shaded orange area in Fig. 4. The accuracy of our prediction, however, is better than the precision.The difference between the predicted and measured [15] phase boundaries is ∼5% of the pressure, which is a typical uncertainty by ab initio calculations (e.g., see [50]).The error bars in Fig. 4 label the reported experimental uncertainties in P, T [15].A recent estimation of the uncertainty in the experimental transition pressure [15] resulting from the choice of the Fe EOS [51] as a pressure scale is also ∼5% of the pressure [52] and is shown as the shaded gray area in Fig. 4. As such, the level of agreement between our predictions and measurements of this phase boundary is excellent.
In the temperature range shown, the B8↔B2 transition occurs at P > ∼240 GPa.The melting properties of FeO are beyond the scope of this study.At the innercore boundary (ICB) pressure, 329 GPa [15], the calculated transition temperature is 2490 K.Under such conditions, the B8→B2 transition is accompanied by a 1.5% density increase.Like the analogous hcp↔bcc transition in elemental metals [28,40,53], the Clapeyron slope of the B8↔B2 transition is also negative, −52 ± 5 MPa/K, which is in excellent agreement with that measured by experiments, −50 MPa/K [15].
In summary, we have investigated the B8↔B2 phase boundary of FeO at high P, T conditions of the Earth's core with ab initio calculations.We computed anharmonic free energies in the thermodynamic limit using phonon quasiparticle dispersions.The calculated phase boundary agrees with experimental observations [15] within uncertainties.The successful calculation of the B8↔B2 phase boundary demonstrates that the PBE-GGA functional describes well energy differences between the B8 and B2 phases of FeO at these P, T conditions.We might attribute this success to two factors: 1) FeO is metallic and non-magnetic at the relevant P, T s.A comparison between the electronic density of states calculated with PBE-GGA and DFT+DMFT [20]  P, T s, anharmonic effects on the B2 phase are dominant, and the vibrational entropy differences between the two phases are also well described by PBE-GGA.Therefore, the present results establish a theoretical framework for future predictive studies of the Fe-FeO system at these high P, T conditions, which is key to understanding the debated problem of oxygen partitioning between the liquid and the solid regions of the Earth's core and their density deficits.

FIG. 2 ., 1 4 , 1 4 (
FIG. 2. (a) P, T conditions covered by AIMD simulations.Hollow circles indicate conditions where the B2 structure is unstable and/or PHQs are not well-defined.Red circles indicate a dynamically stable B2 phase with well-defined PHQs.Blue crosses indicate a dynamically stable B8 phase.Pressures are ensemble averages from AIMD simulations.Modeprojected VAF of a TA mode at q = 1 4 , 1 4 , 1 4 (gray circle in Fig.1(a)) with a harmonic frequency of 304 cm −1 for B2 at V = 5.52 Å3 /atom and T = (b) 1000 and (d) 2000 K, respectively.(c) and (e) show the corresponding power spectra.

FIG. 3 .
FIG. 3. (a) Vibrational entropy S vib (V, T ) and (b) electronic entropy S el (V, T ) versus T at different V s.(c) Helmholtz free energy F (V, T ) and (d) pressure P (V, T ) versus V at different T s.B2 is shown by solid circles and solid curves, and B8 is shown by open circles and dashed curves.Circles indicate the V, T s sampled by the AIMD simulations.

FIG. 4 .
FIG.4.B8↔B2 phase boundary of FeO.The solid black curve is the calculated phase boundary, and the shaded orange area indicates the computational uncertainty.The x-ray diffraction measurements[15] for B2 (red circles with error bars) and B8 (blue squares with error bars) phases are shown for comparison.The dashed line is the experimental phase boundary[15], and the shaded gray area indicates a recent estimation of its uncertainty[52].The arrow indicates the inner core boundary (ICB) pressure.

Frequency (cm 1 ) 5 .
FIG. S2.Temperature-dependent renormalized frequency of the transverse phonon mode at q = X(0, 1 2 , 0) (gray square in Fig.1(a)) with imaginary harmonic frequency for B2 calculated at constant volume.All soft modes at different volumes acquire real renormalized frequencies and stiffen drastically with increasing temperature.
FIG. S9.Projected electronic density of states (DOS) for (a) B2 and (b) B8 FeO calculated by the PBE-GGA at T el = 1160 K and static P = 275 GPa.Contributions from Fe 3d states are shown by solid curves, and those from O 2p states are shown by shaded orange areas.Projected electronic DOS for (c) B2 and (d) B8 FeO calculated by the DFT+DMFT a at T el = 1160 K and relevant pressures are exhibited for comparison.a E. Greenberg, R. Nazarov, A. Landa, J. Ying, R. Q. Hood, B. Hen, R. Jeanloz, V. B. Prakapenka, V. V. Struzhkin, G. K. Rozenberg, and I. Leonov, Phase transitions and spin-state of iron in FeO at the conditions of earth's deep interior, arXiv:2004.00652.