Solid-liquid density and spin crossovers in (Mg, Fe)O system at deep mantle conditions

The low/ultralow-velocity zones in the Earth’s mantle can be explained by the presence of partial melting, critically depending on density contrast between the melt and surrounding solid mantle. Here, first-principles molecular dynamics simulations of (Mg, Fe) O ferropericlase in the solid and liquid states show that their densities increasingly approach each other as pressure increases. The isochemical density difference between them diminishes from 0.78 (±0.7) g/cm3 at zero pressure (3000 K) to 0.16 (±0.04) g/cm3 at 135 GPa (4000 K) for pure and alloyed compositions containing up to 25% iron. The simulations also predict a high-spin to low-spin transition of iron in the liquid ferropericlase gradually occurring over a pressure interval centered at 55 GPa (4000 K) accompanied by a density increase of 0.14 (±0.02) g/cm3. Temperature tends to widen the transition to higher pressure. The estimated iron partition coefficient between the solid and liquid ferropericlase varies from 0.3 to 0.6 over the pressure range of 23 to 135 GPa. Based on these results, an excess of as low as 5% iron dissolved in the liquid could cause the solid-liquid density crossover at conditions of the lowermost mantle.


Results
We consider two spin phases of (Mg, Fe)O liquid: The one with all Fe ions in low-spin (LS) is a non-magnetic phase and the other with all Fe ions in high-spin (HS) is a magnetic phase. The calculated pressure-volume curves of the LS ferropericlase liquid systematically lie below the corresponding HS curves (Fig. 1) so the low-spin phase shows higher density relative to the high-spin phase. For instance, the volume difference between two spin phases for 25% Fe concentration along the 3000 K isotherm decreases from about 0.6 cm 3 mol −1 at zero pressure to 0.3 cm 3 mol −1 at 100 GPa. They are somewhat smaller than those for solid ferropericlase (Fig. 1). In both cases, the spin associated volume differences are nearly insensitive to temperature and iron concentration.
Our FPMD simulations show that iron undergoes a pressure-induced high-spin to low-spin (electronic spin-pairing) transition in both solid and liquid ferropericlase (Fig. 2), as also found by many previous experiments 17,[27][28][29] and calculations 20,22,30,31 . The spin phase diagram can be calculated in terms of the fraction of low spin state (n) given 20 The predicted spin transition is sharp at 0 K occurring at 52 GPa so that all Fe ions are in the high-spin state (total spin S = 2) on low-pressure side and all in the low-spin state (S = 0) on high-pressure side. The transition pressure is bounded by previous calculations from the above 22 and below 20,31 . However, the transition at elevated temperatures becomes gradual exhibiting a pressure interval of coexisting LS and HS states with their respective proportions, n and (1 − n) determined by Equation 1. Temperature systematically widens the stability field of the mixed-spin phase  Our simulations also show that a gradual high-spin to low-spin transition occurs in liquid ferropericlase over a finite pressure interval about 55 GPa at 4000 K (Fig. 2, top). The spin crossover pressure is positively correlated with temperature: The transition pressure defined by n = 0.5 (that is, the middle point of the crossover) increases and the transition zone broadens with increasing temperature. Remarkably, the spin phase diagram of the liquid ferropericlase ( Fig. 2, top) almost matches with that of the solid ferropericlase (Fig. 2, bottom) along the 4000 K isotherm. The populations of LS states are comparable between the solid and liquid phases, being about 50% around 60 GPa. This correspondence actually exists over an extended pressure range, with low-spin populations being around 80 and 90% at 100 and 135 GPa, respectively. Our results thus suggest that liquid (Mg,Fe)O like its solid counterpart should be in a mixed-spin phase with iron ions perhaps mostly in the low-spin state at pressures corresponding to the deep parts of the lower mantle.
It is noted that the spin phase diagrams obtained here for both solid and liquid ferropericlase show notable differences with those from the recent FPMD calculations 22,30 . Interestingly, the discrepancies could be largely attributed to the use of the Hubbard U term ( Supplementary Fig. 2). Our GGA + U test simulations show that the inclusion of U shifts the HS-LS transition pressure in solid ferropericlase from 52 to 63 GPa at 0 K. The pressure shift becomes larger at higher temperature. The increased transition pressures corresponding to n = 0.5 at 0 and 4000 K are comparable with the previous GGA + U pressures 22 (Fig. 2, bottom). For ferropericlase liquid, the effects of U are so large that the pressure of co-existing HS and LS in equal proportions goes beyond 150 GPa at 6500 and 8000 K, again being consistent with the previous GGA + U prediction 30 . The difference in the Clapeyron slopes for the spin crossovers between the present and previous studies (Fig. 2) can be attributed to mainly the entropy factors ( Supplementary Fig. 3). Unlike the spin phase diagram, the U-term does not affect the density significantly as discussed later ( Supplementary Fig. 4).
The calculated pressure-volume-temperature (P − V − T) results ( Fig. 3) for all low spin iron-bearing liquids can be described with one single equation of state: The first term (Equation 2) represents the reference isotherm (T 0 = 3000 K) of pure MgO liquid taken to be the third order Birch-Murnaghan equation: V 0 (Å 3 ) = 30.0, K 0 (GPa) = 16.6, ′ = . K 6 14 0 . For all concentrations of low spin iron, the pressure-volume results almost match with those of the pure liquid. The calculated thermal contributions are the same within the computational uncertainty for pure and alloyed melts. They all increase considerably with compression thus requiring a strong volume-dependent coefficient B TH . This behavior is in sharp contrast to the case of the solid phase where the thermal pressure is weakly dependent on volume 13,14,32 . The thermal contributions in the liquid are systematically larger than those in the solid at most compressed volumes. . The value of n varies from 0 (i.e., 100% high-spin state) to 1 (i.e., 100% low-spin state). The black curves represent the n = 0.5 isolines, that is, the mid points of the spin crossovers in the P-T space. Also shown for solid ferropericlase are the experimental data (pluses 27 , asterisks 29 , horizontal white lines 28 ) and the previous GGA + U calculations (blue line) 30 along with the present GGA + U results (circles).
Scientific RepoRts | 6:37269 | DOI: 10.1038/srep37269 The effects of low-spin iron on pressure are small, compared to the high-spin iron effects ( Supplementary Fig. 5). The corresponding coefficient (B Fe ) in Equation 2 changes from slight negative to positive at high compression.
The calculated density-pressure isotherms for the pure and Fe-bearing liquids remain mostly parallel with each other (Fig. 4). So do the calculated solid density-pressure profiles. However, the isotherms are not parallel between the liquid and solid phases due to large difference in their compressibility. The liquid density increases more rapidly compared to the solid density with compression, and two densities approach each other as pressure increases. At the ambient pressure and 3000 K, the crystal-liquid density difference (∆ ρ ) is large (0.71 g/cm 3 for the pure phases). As pressure increases, the density gap decreases rapidly initially (0.25 g/cm 3 at 23 GPa and 3000 K) and then more gradually at high pressure (0.14 g/cm 3 at 135 GPa and 4000 K). Temperature widens the density gap somewhat because the liquid has larger thermal pressure (and hence larger thermal expansion) than the solid as our simulations have found. The gap widening at higher iron content has origin in the mass itself. For 25% Fe, the calculated values of ∆ ρ are 0.84, 0.30, and 0.18 g/cm 3 , respectively, at 0, 23, and 135 GPa. Our first-principles findings thus confirm that any iso-chemical solid-liquid density gap remains positive finite under conditions of the entire mantle.

Discussion
The density-pressure profiles of the liquid and solid phases for different (Mg, Fe)O compositions intersect with each other implying the possibility of solid-liquid density crossover (Fig. 4, Supplementary Fig. 6). Using these results, we can estimate the pressure at which the density of the liquid phase with a certain excess amount (∆ X) of iron dissolved in the liquid exceeds the solid density. This density inversion condition is defined by where X is the iron content of the solid (Mg 1−X Fe X )O phase. Corresponding to four iron concentrations that we simulated in this study, we can consider ∆ X values of 6.25, 12.5, 18.75 and 25%. Each value has more than one instance thereby making multiple pressure estimates possible, for instance, ∆ X of 6.25% has four comparable estimates. Thus derived ∆ X-pressure profiles along different isotherms (Inset of Fig. 4) show that the critical excess iron content (∆ X) of the liquid decreases rapidly initially with increasing pressure, dropping below 10% after 20 GPa and then below 5% after 100 GPa. A higher ∆ X value of 20% can cause the liquid-solid density crossover in ferropericlase at pressures as low as 10 GPa, as also found by the recent GGA + U calculations along the liquidus 30 . We now discuss the relevance of the estimated excess iron content (∆ X) for the density crossover to occur in ferropericlase by considering Fe partitioning between the solid and liquid states. The partition coefficient can be calculated from energetics of the following reaction: The equation 5 requires full knowledge of the Gibbs free energy, G = H − TS of all systems. The enthalpy H can be readily calculated from our simulations. The calculation of the entropy S is intensive, requiring the thermodynamics integration method 30,33 , which was not used in this study. However, it is likely that the vibrational and configuration contributions from different phases in the above reaction (Equation 4) will largely cancel out. Considering the same (low) spin state of all Fe atoms for the solid and liquid ferropericlase, we can also ignore the magnetic entropy. The net electronic contribution to the entropy calculated from the simulations turns out to be small and is included. Based on these arguments, it makes sense to use the enthalpy term only 34 , which was calculated as a function of pressure for different iron concentrations. The estimated value of ∆ G is 0.40 to 0.18 eV over the pressure range (23 to 135 GPa) of the lower mantle.
The following thermodynamic relation defines the iron partitioning between the solid and and liquid ferropericlase: Based on our calculated spin phase diagrams (Fig. 2), the solid and liquid phases likely contain low spin iron in similar amounts at any high pressure. Even if the spin transition zones may not match between the solid and liquid phases as found by the GGA + U calculations 30 , the transitions are very broad so each phase would still be in a mixed-spin state. This means that the spin transitions cannot impact the possible crystal-liquid density crossover significantly. By considering the pressure, temperature, and spin factors together for the (Mg, Fe)O system for several concentrations, our analysis shows that the excess iron content of ferropericlase melt for it to be denser than the solid phase can be as low as 5% at conditions of the lowermost mantle. Small amounts of excess iron (∆ X) mean that the dense melts may also see preferred partitioning of lighter incompatible elements like H, He, S, etc. The deepest part of the Earth's mantle perhaps existing as either a deep-seated liquid generated in contact with hot core or a residual dense magma ocean can be enriched in minor and trace elements as previously suggested 38,39 . Methods First principles molecular dynamics simulations of (Mg, Fe)O ferropericlase were performed using spin-polarized GGA 40 and projector augmented wave method 41 . Simulations based on the canonical (NVT) ensemble were performed over a wide volume (V) range resulting in pressure range of 0 to 140 GPa at temperature (T) from 3000 to 8000 K. The number (N) of atoms used were 64 for the small supercell with which most simulations were performed, compared to the large supercell (216 atoms) that was used to assess the finite system-size effects. Four iron concentrations including 6.25, 12.5, 18.75 and 25% were considered. For each composition, the initial atomic configuration was first melted at high temperature (6000 to 10000 K depending on the volume) and then quenched/cooled down to desired lower temperatures subsequently. The run durations varied from 10 to 100 ps (with a time step of 1 femtosecond used), depending on the volume and temperature condition of the simulation. Averages of the energy and pressure were computed by the blocking method 42 . The liquid state of the simulated system was assured by examining the mean square displacement and radial distribution function plots. The Pulay stress arising from the use of a finite cutoff of 400 eV at Γ point was added. Further details can be found elsewhere 5,13 .
For systems with correlated (d or f) electrons, a correction formally known as Hubbard U to the self-consistent solution is usually included in order to better describe the electronic structure features, e.g., opening up of band gap, better experimental agreement of the band width, particularly at the ambient conditions. The appropriateness of the U term at elevated conditions of pressure and temperature is not well founded, however. We have also performed a few GGA + U calculations for the comparison purpose. They showed that the effects of the U term 43 on the equation of state and density are small (Supplementary Fig. 1 and Fig. 4). The effects on the electronic density of states and the spin-crossovers are substantial ( Supplementary Figs 7-9), being mostly consistent with the recent GGA + U calculations 22,30 . As far as the density differences between the solid and liquid phases are concerned, they appear to largely be insensitive to the choice of GGA/GGA + U at all conditions ( Supplementary Fig. 4). So are the density differences between the high spin and low spin Fe-bearing ferropericlase.