# Oxygen in the Earth’s core: a first principles study

## 1 Introduction

The Earth’s liquid outer core consists mainly of iron, but its density is about 10% too low to be pure iron (Birch, 1952), so that it must contain some light element. The nature of this element is still uncertain, and during the last 45 years the main candidates have been carbon (Birch, 1952; Clark, 1963; Urey, 1960; Wood, 1993), silicon (Birch, 1952; MacDonald and Knopoff, 1958; Ringwood, 1959, 1961, 1966), magnesium (Alder, 1966), sulphur (Clark, 1963; Urey, 1960; Birch, 1964; Mason, 1966; Murthy and Hall, 1970; Lewis, 1973), oxygen (Dubrovskiy and Pan’kov, 1972; Bullen, 1973; Ringwood, 1977), and hydrogen (Birch, 1952; Fukai and Akimoto, 1983; Suzuki et al., 1989). For a given light element, it is also uncertain what concentration is needed to explain the inferred density in the core. The arguments for and against each of the candidate light elements have been reviewed by Poirier (1994).

The aim of this paper is to use first-principles calculations to investigate the possibility that oxygen is the light element. First-principles calculations are well established as a reliable way of predicting the thermodynamic, structural and dynamical properties of solid and liquid materials, including liquid metals (Štich et al., 1989; Kresse and Furthmüller, 1993). We have recently reported calculations of this kind on pure liquid iron under core conditions (Vočadlo et al., 1997; de Wijs et al., 1998), which show that it is a simple close-packed liquid with a viscosity not much greater than that of many liquid metals at ambient pressure, contrary to some earlier suggestions (Secco, 1995). We have also used first-principles simulations to investigate a liquid iron-sulphur alloy under the same conditions (Alfè and Gillan, 1998a). We showed that the properties of the liquid are scarcely affected by the small sulphur concentration needed to explain the observed density.

The proposal that oxygen is the light element has a long and controversial history. Among the earliest proponents were Dubrovskiy and Pan’kov (1972) and Bullen (1973), the latter of whom suggested an outer core composition in the region of Fe<sub>2</sub>O (equivalent to 12.5 wt$`\%`$). Ringwood (1977) argued that oxygen should be seriously considered, and used seismic data to estimate the oxygen content as 28 mol percent (10 wt$`\%`$). However, it is not completely certain whether the Fe/O liquid is thermodynamically stable against phase separation under Earth’s core conditions at the composition that would be necessary to explain the density. It is known that the solubility of FeO in liquid Fe is very low ($`1`$ mol percent) near the melting temperature of pure iron (1811 K) at atmospheric pressure (Distin et al., 1971). However, the solubility increases rapidly with temperature, becoming $`6.5\%`$ at 2350 K (Fischer and Schumacher, 1978) and rising to the region of $`35`$ mol $`\%`$ at 2770 K (Ohtani and Ringwood, 1984). According to Ohtani and Ringwood (1984), an extrapolation of the available phase measurements would suggest that the region of immiscibility disappears entirely above $`3080`$ K at atmospheric pressure. It is also well established that the solubility of FeO in Fe increases with increasing pressure, and that the partial molar volume of FeO in liquid Fe is lower than that of pure liquid FeO itself. Ohtani et al. (1984) used their high pressure measurements on the solubility of FeO to suggest that the Fe/O system may show simple eutectic behaviour above a pressure of $`20`$ GPa, with no region of liquid immiscibility at any temperature. Subsequently, Ringwood and Hibberson (1990) showed by direct measurements that at 16 GPa addition of FeO to pure iron causes a depression of melting point, leading to a eutectic point at oxygen mole fraction of 28 $`\%`$ and a temperature of ca. 1940 K. Boehler’s (1992) measurements on the melting of Fe/O mixtures are consistent with these ideas. Experiments of Knittle and Jeanloz (1991) and Goarant et al. (1992) on the reaction between lower mantle material and molten iron at pressures above 70 GPa revealed that the liquid dissolves significantly amount of FeO. However, Sherman (1995) has recently used first-principles calculations on crystalline Fe/O phases to argue strongly against significant amounts of oxygen in the core. His calculations gave values for the enthalpy of formation of crystals of composition Fe<sub>3</sub>O and Fe<sub>4</sub>O starting from Fe in the hexagonal close packed structure and FeO in the NiAs structure. The Fe<sub>3</sub>O and Fe<sub>4</sub>O compositions were used to model substitutional and interstitial oxygen respectively. The enthalpies of formation were found to be so large that phase separation into FeO and Fe appears to be inevitable. However, it is not clear that Sherman’s results have any relevance to the outer core, since in the liquid phase oxygen does not have to be either substitutional or interstitial. Even for the solid phase it is not obvious that Sherman’s argument is robust, since Fe/O crystal structures other than those he studied might well give much lower enthalpies of formation.

We are mainly concerned in this paper with the Fe/O system in the liquid state. Our first-principles calculations, based on density functional theory and the pseudopotential method, will be used to address three questions: (a) Is the Fe/O liquid stable against phase separation under Earth’s core conditions? (b) If it is, what oxygen concentration is needed to reproduce the observed density? (c) At this concentration, do the structural and dynamical properties of the liquid differ appreciably from those of pure liquid iron at the same pressure and temperature? Our first-principles simulations of the liquid will provide strong evidence that it is thermodynamically stable, and that the observed density requires an oxygen concentration of $`2530`$ mol-percent. In studying the properties of the liquid, we shall be particularly concerned with the viscosity, since this is one of the most poorly determined properties of the outer core, with estimates from different experimental and theoretical methods spanning many orders of magnitude (Secco, 1995). We shall also investigate a number of other properties, including the nature of the short-range order, the atomic diffusion coefficients, and the electronic structure.

Although we are mainly interested in the liquid, we shall also present some results for the energetics of various crystalline forms of the Fe/O system. These crystal calculations serve two purposes: first, they demonstrate that our techniques are in complete agreement with those used by Sherman in predicting large enthalpies of formation for Fe<sub>3</sub>O and Fe<sub>4</sub>O in the structures he assumes; second, they demonstrate that there are crystal structures that give much lower formation enthalpies – an important fact in understanding how the Fe/O liquid can be stable against phase separation. It is not our intention to come to definite conclusions about the possible phase stability of the Fe/O solid solutions themselves, but our calculations suggest that their stability cannot be ruled out.

The paper is organised as follows. In section 2, we summarise the first-principles techniques on which the work is based. Section 3 presents our calculations on the energetics of Fe/O crystals. In section 4 we report our results on liquid Fe/O, including its structural properties, the evidence for its stability against phase separation, and its dynamical and electronic properties. The final sections present discussion and conclusions.

## 2 Methods

In first principles calculations, the solid or liquid is represented as a collection of ions and electrons, and for any given set of ionic positions the aim is to determine the total energy and the force on every ion by solving the Schrödinger equation. This is a formidable task if the number of atoms is large, but it was made feasible by the introduction of density functional theory (DFT) many years ago (Hohenberg and Kohn, 1964; Kohn and Sham, 1965; Jones and Gunnarsson, 1989; Parr and Yang, 1989). DFT treats electronic exchange and correlation in a way that allows the electrons to be described by single-particle wavefunctions, with the interaction between them accounted for by an effective potential. DFT can be applied in two ways: all-electron calculations, or pseudopotential calculations. The first approach includes such standard techniques as full-potential linearised augmented plane waves (FLAPW) and linearized muffin-tin orbitals (LMTO). In the pseudopotential approach, only valence electrons are explicitly treated, the effect of the core electrons being included by an effective interaction between the valence electrons and the cores. In both approaches, the accuracy with which the real material is described is governed by the approximation used for the electronic exchange-correlation energy. Until recently, the local density approximation (LDA) was the standard method. But in order to achieve the highest accuracy for transition metals it is essential to use an improved method known as the generalized gradient approximation (GGA) (Wang and Perdew, 1991). The present work is based on the pseudopotential approach and the GGA. A non-technical review of first-principles calculations based on the pseudopotential approach has been given recently by one of the authors (Gillan 1997).

There have already been extensive first-principles calculations on crystalline iron both at ambient pressure and at pressures going up to Earth’s core values (Stixrude et al., 1994; Söderlind et al., 1996). The calculations have been performed using different all-electron techniques and the pseudopotential technique, and a variety of properties have been studied, including the equilibrium volume, the elastic constants, the magnetic moment, the volume as a function of pressure and lattice vibration frequencies. The agreement between results obtained with different techniques is generally very close, and the agreement with experimental data is also good. Particularly relevant here is the recent comparison of the pseudopotential results for the pressure-dependent volume of hexagonal-close-packed iron up to core pressures with earlier all-electron results and with experimental measurements (Vočadlo et al., 1997).

Static first-principles calculations on crystals have been in routine use for many years. But to study liquids we need to do dynamical first-principles simulations, in which the calculated forces on the atoms are used to generate time evolution of the system, with every atom moving according to Newton’s equation of motion. This kind of first-principles molecular dynamics (FPMD) pioneered by Car and Parrinello (1985) has been extensively used to study liquid metals (Štich et al., 1989; Kresse and Furthmüller, 1993; Holender and Gillan, 1996; Kirchhoff et al. 1996a, 1996b), and it is known to give an accurate description of both structure and dynamics.

The present work was performed mainly with the VASP code (Vienna Ab initio Simulation Package) (Kresse and Furthmüller, 1996a, 1996b). As usual in pseudopotential work, the electron orbitals are represented using a plane-wave basis set, which includes all plane waves up to a specific energy cut-off. The electron-ion interaction is described by ultrasoft Vanderbilt pseudopotentials (Vanderbilt, 1990), which allow one to use a much smaller plane-wave cut-off while maintaining high accuracy. When we perform FPMD with VASP, the integration of the classical equation of motion is done using the Verlet algorithm (1967), and the ground-state search is performed at each time-step using an efficient iterative matrix diagonalisation scheme and a Pulay mixer (1980). This method differs from the original Car-Parrinello technique, which treated the electronic degrees of freedom as fictitious dynamical variables. In order to improve the efficiency of the dynamical simulation, the initial electronic charge density at each time step is extrapolated from the density at previous steps as described in our previous work (Alfè and Gillan, 1998a). In FPMD on metals, the discontinuity of occupation numbers at the Fermi level can cause technical difficulties. Since we are interested in high temperatures in our liquid simulations, the electronic levels are occupied according to Fermi statistics corresponding to the temperature of the simulation. This prescription also avoids problems with level crossing during the ground state search. The FPMD simulations are performed at constant temperature (rather than at constant energy), using the Nosé technique (1984).

The iron pseudopotential we use is the same as that used in our earlier work (Vočadlo et al., 1997; de Wijs et al., 1998; Alfè and Gillan, 1998a,b), and was constructed using an Ar core and a $`4s^13d^7`$ atomic reference configuration. The oxygen pseudopotential was constructed using a He core and the $`2s^22p^4`$ reference configuration. At Earth’s core pressures the distance between the atoms becomes so small that the Fe($`3p`$) orbitals respond significantly. The net effect is a small repulsion, which we determined from calculations on the h.c.p. crystal using a Ne core for the iron pseudopotential instead of Ar. We represent this repulsion by Fe-Fe and Fe-O pair potentials, as in our earlier work (Vočadlo et al., 1997; de Wijs et al., 1998). The resulting corrections to energy and forces are generally small. Non-linear core corrections (Louie et al., 1982) are included throughout the work.

For all the calculations to be reported, we use a plane-wave cut-off of 400 eV, which gives total energies converged within $`1020`$ meV/atom. In the calculations on crystals Brillouin-zone sampling is an important issue, and the sampling density we use is described in the next section. But for the liquid we use $`\mathrm{\Gamma }`$ point sampling, which experience suggests should be satisfactory. (We have done separate tests on pure liquid iron using 4 k-points, and we have found no detectable structural effects, while the average total energy difference with respect to the $`\mathrm{\Gamma }`$ point only calculations is of the order of 10 meV/atom, which is completely negligible for the purposes of the present work). The time step used in the dynamical simulations was 1 fs and we generally used a self-consistency threshold of $`1.5\times 10^7`$ eV/atom. With these prescriptions the drift of the Nosé constant of motion was less than $`60`$ K per ps.

## 3 Fe/O solid solutions

We have used the techniques described in the previous section to calculate the equilibrium properties and the enthalpy of formation $`\mathrm{\Delta }H`$ of various members of the Fe/O system, including those studied by Sherman (1995). Although Sherman’s calculations and ours are both based on DFT, the technical methods are completely different, since Sherman used the FLAPW method, whereas ours are based on the pseudopotential approach. One of our main aims is therefore to make detailed comparisons with his results in order to ensure that the two methods agree about the energetics of the systems. We shall also point out that there are Fe<sub>3</sub>O structures with enthalpies of formation much lower than those reported previously.

We begin by presenting our results for the equilibrium density $`\rho _0`$ (i.e. the density for which the pressure is zero), the bulk modulus $`K`$ at this density, and the pressure derivative $`K^{}dK/dP`$ for crystals of Fe, FeO, Fe<sub>3</sub>O and Fe<sub>4</sub>O. Pure Fe is in the $`ϵ`$ structure (hexagonal close packed); FeO is in the B8 structure (the NiAs structure); Fe<sub>3</sub>O is in the structure obtained from face-centred cubic Fe by replacing the atoms at the corners of the conventional cube by O atoms; and Fe<sub>4</sub>O is in the structure obtained from f.c.c. Fe by inserting an O atom at the centre of the conventional cube (see Fig. 1). We have done the calculations both spin restricted (the occupation numbers of every electronic orbital are equal for up and down spins, so that there are no magnetic moments) and spin unrestricted (the occupation numbers for up and down spins are allowed to vary independently). To sample the Brillouin zone (BZ) we have used Monkhorst and Pack grids (1976) using the sampling level that corresponds to 20 and 36 k-points in the irreducible wedge of the BZ respectively for the cubic (Fe<sub>3</sub>O and Fe<sub>4</sub>O) and the hexagonal (Fe($`ϵ`$) and FeO(B8)) crystals. Using these values, the total energies are converged within 5 meV/atom for the Fe<sub>3</sub>O, Fe<sub>4</sub>O and FeO(B8) structures, and within $`1015`$ meV/atom for the Fe($`ϵ`$) structure. We find that, except for Fe, all the structures are magnetic at low pressures, so that the system is stabilised if spin restriction is removed. This disagrees with the results of Sherman (1995), where a significant magnetic moment was found only for the FeO(B8) structure (a weak moment was found for Fe<sub>3</sub>O).

Since we regarded the disagreement concerning magnetic properties as disturbing, we repeated our calculations of the equilibrium magnetic moment of Fe<sub>3</sub>O and Fe<sub>4</sub>O using a completely independent electronic structure technique, namely the LMTO method (linearised muffin-tin orbitals), using the LMTO-46 code due to Krier et al. (1994). These calculations completely confirm the magnetic ordering in Fe<sub>3</sub>O and Fe<sub>4</sub>O and give numerical values for the magnetic moments that agree well with those given by our pseudopotential calculations. This suggests that the minimisation of the total energy with respect to magnetic moments may not have been fully under control in Sherman’s work.

For each system we have calculated the static internal energy $`E`$ for a series of volumes, and fitted the results to the Birch-Murnaghan equation of state:

$`E=E_0+{\displaystyle \frac{3}{2}}V_0K\left[{\displaystyle \frac{3}{4}}(1+2\xi )\left({\displaystyle \frac{V_0}{V}}\right)^{4/3}{\displaystyle \frac{\xi }{2}}\left({\displaystyle \frac{V_0}{V}}\right)^2{\displaystyle \frac{2}{3}}(1+\xi )\left({\displaystyle \frac{V_0}{V}}\right)^{2/3}+{\displaystyle \frac{1}{2}}\left(\xi +{\displaystyle \frac{3}{2}}\right)\right]`$ (1)   
$`\xi ={\displaystyle \frac{3}{4}}(4K^{}),`$   
where $`K`$ is the zero pressure bulk modulus, $`K^{}=(dK/dP)_{P=0}`$, $`E_0`$ is the equilibrium energy and $`V_0`$ the equilibrium volume. Our calculated values of $`\rho _0`$, $`K`$ and $`K^{}`$ are compared with Sherman’s results in Table 1, which also reports our calculated magnetic moments. The overall conclusion from the comparison is that the results agree well in the cases where magnetism is absent: pure Fe and spin-restricted Fe<sub>3</sub>O and Fe<sub>4</sub>O (Sherman’s calculations are effectively spin-restricted for Fe<sub>3</sub>O and Fe<sub>4</sub>O, since he found no moments). The equilibrium density agrees in those cases to better than $`2\%`$, and the values of $`K`$ and $`K^{}`$ to about $`10\%`$. However, our results show that magnetism has a strong effect on the equilibrium properties, so that it is important to treat it correctly. The case of FeO(B8) is problematic. Our spin-unrestricted calculations give a $`\rho _0`$ value in respectable agreement with that of Sherman, but the agreement is poor for $`K`$ and $`K^{}`$. It is clear that the accurate treatment of the volume dependence of magnetic moment is important in obtaining reliable values for these parameters, and our suspicion is that problems with moments may have affected Sherman’s values. However, we shall stress below that magnetic effects become unimportant at core pressures, so the most important feature of Table 1 is the good agreement between the two sets of calculations for the non-magnetic cases.

We turn now to the enthalpies of formation $`\mathrm{\Delta }H`$ of Fe<sub>3</sub>O and Fe<sub>4</sub>O, defined by:

$`\mathrm{\Delta }H(\mathrm{Fe}_3\mathrm{O})=H(\mathrm{Fe}_3\mathrm{O})H(\mathrm{FeO}(\mathrm{B8}))2H(\mathrm{Fe}(ϵ)),`$   
$`\mathrm{\Delta }H(\mathrm{Fe}_4\mathrm{O})=H(\mathrm{Fe}_4\mathrm{O})H(\mathrm{FeO}(\mathrm{B8}))3H(\mathrm{Fe}(ϵ)),`$ (2)   
where $`HE+PV`$ is the enthalpy per formula unit of each material. The total energy is taken directly from our Birch-Murnaghan fit, and the pressure $`P=dE/dV`$ is calculated from the derivative of the fitted form. We have done these calculations both spin restricted and spin-unrestricted, and our results are reported in the two panels of Fig. 2. In the spin-restricted panel we also show Sherman’s enthalpy results. The most important conclusion is that $`\mathrm{\Delta }H`$ for both Fe<sub>3</sub>O and Fe<sub>4</sub>O becomes very large at core pressures. In the range from 135 GPa (core-mantle boundary) to 330 GPa (inner-core boundary) $`\mathrm{\Delta }H`$ is at least 3 eV, which corresponds roughly to a temperature of $`3.5\times 10^4`$ K, so that it is exceedingly unlikely that Fe<sub>3</sub>O and Fe<sub>4</sub>O could be thermodynamically stable in the assumed structures.

Our spin unrestricted results indicate that the true values of $`\mathrm{\Delta }H`$ are considerably lower than Sherman’s results at low pressure, but at high pressure the differences between the magnetic and non-magnetic values of $`\mathrm{\Delta }H`$ become very small, so the conclusion is unaffected. The detailed agreement with Sherman’s values of $`\mathrm{\Delta }H`$ is only moderately good, but again this does not affect the conclusions about the very large size of $`\mathrm{\Delta }H`$.

We now want to ask whether the assumed crystal structures for Fe<sub>3</sub>O and Fe<sub>4</sub>O are actually the most stable. A glance at Wyckoff’s book Crystal Structures (1964) shows that compounds having the composition A<sub>3</sub>B crystallise in a bewildering variety of structures. We have picked some likely candidates and calculated their formation enthalpy. Most turn out to be unfavourable, with $`\mathrm{\Delta }H`$ values at least as great as those already reported in Fig. 2. However, we have discovered one that has a much lower value. This is the BiI<sub>3</sub> structure, which has a rhombohedral unit cell containing two formula units. Putting Fe<sub>3</sub>O into this structure and relaxing both the atomic positions and the shape of the cell, we end up with the triclinic structure shown in Fig. 1. To characterise the structure briefly, we note that at 300 GPa each oxygen is surrounded by 11 Fe neighbours at distances of between 1.76 and 2.57 Å, and two O atoms at distances of 2.02 and 2.34 Å (by contrast, in the cubic Fe<sub>3</sub>O structure at the same pressure, each oxygen has 12 Fe neighbours at 2.05 Å, and the nearest oxygen neighbours are at 2.9 Å). We have calculated the fully relaxed total energy of this structure at several volumes, and the resulting structural parameters are reported in Table 2. Spin unrestricted calculations show that the structure is weakly magnetic, but the moment and the energy stabilisation are so small that the effects can be ignored.

Fitting of the energies to the Birch-Murnaghan equation of state yields the enthalpy of formation shown in Fig. 2. Remarkably, $`\mathrm{\Delta }H`$ is very much lower than for the previous structures, and it decreases with increasing pressure. At the pressure of the inner core boundary it is only just over 1 eV. Since we arrived at this distorted BiI<sub>3</sub> structure in a rather haphazard way, it is quite likely that there are other Fe<sub>3</sub>O structures with even lower enthalpies. We cannot say at present whether there are Fe<sub>3</sub>O structures that are stable against phase separation under Earth’s core conditions, but it certainly does not look impossible. The low $`\mathrm{\Delta }H`$ for the BiI<sub>3</sub> structure will be highly relevant to our study of phase stability in the liquid Fe/O system.

## 4 The liquid

### 4.1 Thermodynamics

Our aim in choosing the thermodynamic parameters for our liquid simulations was to model liquid Fe/O near the thermodynamic state it would need to have to reproduce the known outer-core density at the inner core boundary (ICB). The temperature at this point is very uncertain, with estimates ranging from 4000 to 8000 K (Poirier, 1991). We took the value of 6000 K, which is intended to be a reasonable compromise. However, the density and the pressure are quite accurately known to be $`12000`$ kg m<sup>-3</sup> and $`330`$ GPa. This density is about $`10\%`$ lower than it would be if the core consisted of pure iron (Birch, 1952). The main problem in choosing thermodynamic parameters is that we do not know in advance the required oxygen concentration, so that a certain amount of trial and error is needed.

We started from our previous 64-atom simulation for pure liquid iron (Vočadlo et al., 1997), which had a mass density of $`13300`$ kg m<sup>-3</sup> and a calculated pressure of $`358\pm 6`$ GPa. Our first move was to hold the volume of the system fixed and to transmute the appropriate number of iron atoms into oxygen atoms to produce the density of $`12000`$ kg m<sup>-3</sup>. This resulted in a large reduction of the pressure, and we therefore reduced the cell volume to restore the original pressure. Naturally, this increased the density, and we therefore converted more iron atoms into oxygen to regain the density of $`12000`$ kg m<sup>-3</sup>. By repeating this cycle many times, one could in principle achieve the required density and pressure. But the calculations are very demanding, since at each state point one has to equilibrate the system and run it for long enough to obtain adequate statistics for the pressure, so that in practice a compromise between accuracy and computational effort is needed. After several iterations, we ended up with a simulation box containing 43 iron atoms and 21 oxygen atoms, i.e. mole fractions of $`x_{\mathrm{Fe}}0.67`$ and $`x_\mathrm{O}0.33`$. The resulting mass density of $`11600`$ kg m<sup>-3</sup> and pressure of $`342\pm 4`$ GPa are close to the known values at the ICB.

Since the mass density of $`11600`$ kg m<sup>-3</sup> is slightly below the known value at the ICB, it is likely that the concentration of $`x_\mathrm{O}0.33`$ is an overestimate. We have therefore taken a second thermodynamic state with a lower concentration. In order to facilitate comparisons with our calculations on crystalline Fe<sub>3</sub>O, we chose the value $`x_\mathrm{O}=0.25`$. This second simulation was performed on a system of 48 iron atoms and 16 oxygen atoms at the mass density of $`12200`$ kg m<sup>-3</sup>, and the resulting pressure was $`366\pm 8`$ GPa. We shall refer to the two simulations in the following as the ‘$`33\%`$ simulation’ and the ‘$`25\%`$ simulation’.

From the thermodynamic results just mentioned, we can estimate the oxygen concentration that would be needed to reproduce the known density and pressure at the ICB. Interpolating between the calculated density values and applying a small correction for the slightly different pressures in the two simulations, we estimate that the mole fraction $`x_\mathrm{O}=0.28`$ would reproduce the density $`12000`$ kg m<sup>-3</sup> at the ICB pressure.

In the next sections we describe the structural, dynamical and electronic-structure properties of the Fe/O liquid alloys.

### 4.2 Structure

We have simulated the $`33\%`$ system for 4.2 ps after 2 ps of equilibration. The structural properties of the system have been inspected by looking at the partial radial distribution functions (rdf), $`g_{\mathrm{FeFe}}(r)`$, $`g_{\mathrm{FeO}}(r)`$, and $`g_{\mathrm{OO}}(r)`$. The partial rdf’s are defined so that, sitting on an atom of the species $`\alpha `$, the probability of finding an atom of the species $`\beta `$ in the spherical shell $`(r,r+dr)`$ is $`4\pi r^2n_\beta g_{\alpha \beta }(r)dr`$, where $`n_\beta `$ is the number density of the species $`\beta `$ (the mole fraction of species $`\beta `$ times the total number of atoms per unit volume).

We have calculated averages of the rdf’s over different small time windows of the simulation and we find no meaningful differences between the windows. This confirms that the system is well equilibrated. In Fig. 3 we display the rdf’s calculated from the whole simulation. These show that the distance between neighbouring iron and oxygen atoms is significantly smaller than the iron-iron distance, the maximum of $`g_{\mathrm{FeO}}(r)`$ being at $`1.7`$ Å, while the maximum of $`g_{\mathrm{FeFe}}(r)`$ is at $`2.1`$ Å. It is interesting to notice that $`g_{\mathrm{OO}}(r)`$ has a first maximum at $`2.1`$ Å, which is much greater than the chemical bond length expected for O-O single or double bonds (1.47 Å and $`1.21`$ Å respectively). This is clear evidence that there is no covalent bonding between oxygen atoms. The presence of the O-O peak at $`2.1`$ Å indicates that oxygen atoms repel each other with an effective atomic diameter of $`2.1`$ Å. This fact shows that oxygen has two effective sizes in the liquid: a small one when it interacts with iron and a large one when it interacts with itself.

It is interesting to compare the structural properties of the alloy with those of pure liquid iron. In Fig. 4 we display the rdf calculated earlier for pure liquid iron at ICB conditions (Vočadlo et al., 1997) and the $`g_{\mathrm{FeFe}}`$ calculated here. The two are not very different, the only apparent effect being the broadening of the peak in the liquid alloy, which is probably due to the greater disorder in the alloy.

The integration of the first peak of the rdf’s provides a definition of the coordination number $`N_{\alpha \beta }^c`$ (the average number of neighbours of species $`\beta `$ surrounding an atom of species $`\alpha `$):

$$N_{\alpha \beta }^c=4\pi n_\beta _0^{r_{\alpha \beta }^c}r^2g_{\alpha \beta }(r)𝑑r,$$
(3)   
where $`r_{\alpha \beta }^c`$ is the position of the minimum after the first peak of $`g_{\alpha \beta }`$. We find the values $`N_{\mathrm{FeFe}}^c=11.0`$, $`N_{\mathrm{FeO}}^c=4.5`$, $`N_{\mathrm{OFe}}^c=9.2`$, and $`N_{\mathrm{OO}}^c=4.5`$. For comparison, the average coordination number found in our earlier simulation of pure liquid iron at ICB conditions was $`N_{\mathrm{FeFe}}^c=13.8`$ (Vočadlo et al., 1997). In interpreting these numbers, it is helpful to consider the coordination numbers that would be found if iron and oxygen atoms had exactly the same size and if atoms were packed in the same way as in pure liquid iron. In that case, the total number of neighbours of each iron atom, $`N_{\mathrm{FeFe}}^c+N_{\mathrm{FeO}}^c`$, would be the same as in pure iron, whereas in fact it is 15.5. This increase of coordination number is clearly due to the smaller size of oxygen, which allows more atoms to be fitted into the first shell of neighbours. On the other hand, the total number of neighbours of each oxygen atom, $`N_{\mathrm{OFe}}^c+N_{\mathrm{OO}}^c`$ is 13.7, which is almost the same as the coordination number in pure iron. We interpret this as the result of two competing effects. The smaller size of oxygen would lead to a smaller coordination number if all atoms in its shell of neighbours were iron. But since on average 4.5 of the neighbours are oxygen, which have a smaller size when interacting with iron atoms in the shell, the coordination number is increased again.

We note that the structure of the liquid is very different from that of the cubic Fe<sub>3</sub>O and Fe<sub>4</sub>O crystals discussed in section 3. Oxygen atoms in these crystals have respectively 12 and 6 iron neighbours. The coordination number of 9.2 in the liquid is roughly half way between the two. In the liquid, the radii from oxygen to iron and oxygen neighbours are equal to $`1.7`$ and $`2.1`$ Å respectively, whereas in the cubic Fe<sub>3</sub>O crystal at a similar pressure, the distances are 2.05 and 2.9 Å. On the other hand, in the BiI<sub>3</sub>-structure Fe<sub>3</sub>O, the O-Fe neighbour separation is spread over the range $`1.762.57`$ Å, and the O-O separation is in the range $`2.022.34`$ Å.

We now want to ask whether our simulated system is really in a single phase and whether we can detect any sign of phase separation. In studying this it is very helpful to calculate the static structure factors $`S_{\alpha \beta }(k)`$ defined by:

$$S_{\alpha \beta }(k)=\rho _\alpha ^{}(𝐤)\rho _\beta (𝐤),$$
(4)   
where $``$ denotes the thermal average (in practice evaluated as a time average). Here, $`\rho _\alpha (𝐤)`$ is the Fourier component of the number density of species $`\alpha `$ at wavevector $`𝐤`$, given by:

$$\rho _\alpha (𝐤)=N_\alpha ^{1/2}\underset{i=1}{\overset{N_\alpha }{}}\mathrm{exp}(i𝐤𝐫_{\alpha i}),$$
(5)   
where $`N_\alpha `$ is the number of atoms of species $`\alpha `$ and $`𝐫_{\alpha i}`$ is the position of the $`i`$th atom of this species. Phase separation is associated with fluctuations of the concentrations of the two species, and the structure factors give us quantitative information about the intensities of these fluctuations.

The connection between phase separation and structure factors can be made more precise. In the limit of zero wavevector, the structure factors of a liquid alloy can be rigorously expressed in terms of thermodynamic derivatives (Bhatia and Thornton, 1970):

$$\underset{k0}{lim}S_{\alpha \beta }(k)=\frac{k_BT}{(N_\alpha N_\beta )^{1/2}}\left(\frac{N_\alpha }{\mu _\beta }\right)_{V,T,\mu _\beta ^{}}.$$
(6)   
where $`\mu _\alpha `$ are the chemical potentials, and the notation indicates that the derivative is to be taken with the volume $`V`$, the temperature $`T`$ and all chemical potential except $`\mu _\beta `$ held fixed. But the condition for thermodynamic stability with respect to phase separation is

$$(\mu _\alpha /x_\beta )_{P,T}>0.$$
(7)   
At the consolute point (the point in the phase diagram at which phases start to separate) the derivatives $`(\mu _\alpha /x_\beta )_{P,T}`$ become zero. This implies that the matrix of derivatives $`(\mu _\alpha /N_\beta )_{V,T,N_\beta ^{}}`$ has vanishing eigenvalues, corresponding to variations of the numbers $`N_\alpha `$ that maintain the pressure constant at fixed volume. But the matrix $`(N_\alpha /\mu _\beta )_{V,T,\mu _\beta ^{}}`$ is the inverse of the matrix $`(\mu _\alpha /N_\beta )_{V,T,N_\beta ^{}}`$, so that when the latter becomes singular the former must acquire infinite eigenvalues. The consequence is that the values of the structure factors in the zero-wavevector limit must diverge if the system is unstable with respect to phase separation. This is also intuitively clear: as one passes from the miscible to the immiscible region, concentration fluctuations become ever larger, becoming of macroscopic size when the phases separate, and the increase in the fluctuations is reflected in the divergence of the quantities $`S_{\alpha \beta }(k0)`$.

Our calculated structure factors for the $`33\%`$ simulation are reported in Fig. 5. They have the form usually found in liquid alloys, with prominent peaks in $`S_{\mathrm{FeFe}}(k)`$ and $`S_{\mathrm{OO}}(k)`$ in the region $`k4`$ Å<sup>-1</sup> signalling the approximate spatial periodicity associated with the packing of the atoms. The significant feature for present purposes is the lack of any anomalous behaviour at small wavevectors. We recognize, of course, that because of the limited size of the repeating simulation cell, there is a lower limit to the wavevector that we can examine, which in the present case is 0.86 Å<sup>-1</sup>. But at least in the accessible region of wavevectors there is no indication of any tendency towards phase separation.

Before leaving the description of structure, we outline another method we have used to search for signs of phase separation. To explain this, let us imagine for a moment that the system had separated into phases of pure Fe and pure FeO. Then the Fe atoms in the Fe phase would have no oxygen neighbours, whereas the Fe atoms in the FeO phase might be expected to have 6 oxygen neighbours (we assume the FeO liquid to have a structure resembling that of crystalline FeO in the NiAs structure). On the other hand, in the unseparated Fe/O phase, the number of oxygen neighbours surrounding each Fe atom fluctuates around the value $`45`$ (see above). We can therefore distinguish between the two situations by studying the probability distribution for the number of oxygen neighbours surrounding Fe atoms. To do this, we use the cut-off distance $`r_{\alpha \beta }^c`$ defined above to decide when an atom of species $`\beta `$ counts as a neighbour of an atom of species $`\alpha `$, and we define the function $`P_{\alpha \beta }(n,r_{\alpha \beta }^c)`$ as the probability that an atom of species $`\alpha `$ has $`n`$ neighbours of species $`\beta `$. If there is a complete phase separation, we expect $`P_{\mathrm{FeO}}`$ to have peaks in the region of $`n=0`$ and $`n=8`$, but if there is no separation we expect a single peak in the region of $`n=45`$. Note that the rdfs contain less information than the $`P_{\alpha \beta }`$ functions, and cannot by themselves deliver the discrimination we need.

We present in Fig. 6 the function $`P_{\mathrm{FeO}}(n,r_{\mathrm{FeO}}^c)`$ calculated from our $`33\%`$ simulation with the cut-off $`r_{\mathrm{FeO}}^c=2.5`$ Å. This shows a single peak at $`n=4`$, and no sign of any structure at lower or higher values of $`n`$. This means that there is no indication whatever of any tendency towards phase separation.

### 4.3 Confirmation of phase stability

Our failure to find any evidence of phase separation strongly suggests that the Fe/O liquid is in fact stable. But it might be objected that a simulation lasting only $`46`$ ps does not give enough time for separation to occur, and that it would occur if the simulation were longer. To eliminate this possibility we have devised a method in which phase separation is artificially induced by an external force. We shall show that when the force is removed the phases spontaneously re-mix very rapidly. We performed this procedure on the $`25\%`$ simulation (it is not significant that the $`25\%`$ case was chosen for this, and we believe that the $`33\%`$ system would have behaved in the same way).

In our procedure, we notionally divide our cubic cell into two parts, consisting of the left region $`0<x<0.4`$ and the right region $`0.4<x<1.0`$ ($`x`$ is the coordinate along one of the edge directions in units of the cell length). The first step is to sweep all the oxygen atoms into the left region with an external force, so that this region contains something resembling FeO, while the right region contains pure Fe. To achieve this, we apply a constant force along the $`x`$-axis to all oxygen atoms lying in the right region. This force is in the positive $`x`$ direction for $`0.7<x<1.0`$ and in the negative direction for $`0.4<x<0.7`$. No force is applied to the iron atoms, and these are left free to redistribute themselves. Initially, the magnitude of the force was taken to be 1 eV/Å, but since this proved to be too weak it was increased to 3 eV/Å. After $`1`$ ps a complete phase separation was achieved, with all oxygen atoms in the left region, and we let the system equilibrate with the external force still present for a further 1 ps. We show in Fig. 7 a snapshot of one configuration taken from this period, which clearly shows the complete separation of phases. The external force was then switched off and the system was allowed to evolve for a further 2 ps. Remarkably, we found that after only 1 ps had elapsed the oxygen atoms became completely randomized throughout the cell, as can be seen in the snapshot shown in in the lower part of Fig. 7.

To characterise these events quantitatively, we use the probability function $`P_{\mathrm{FeO}}(n,r_{\mathrm{FeO}}^c)`$ described in section 4.2. Fig. 8 shows this function calculated by averaging over three short windows of 0.1 ps each, starting 0 ps , 0.5 ps and 1.0 ps after the external force was switched off. The first window shows a clear bimodal form, as would be expected for a two phase system. (We note that because our system is rather small, many atoms are near the boundary between the phases, so that the peaks in $`P_{\mathrm{FeO}}`$ are not as sharp as they would probably be in a larger system.) After 0.5 ps, the distribution has already become unimodal, and after 1 ps it is very similar to what we showed in Fig. 6 for the equilibrium $`33\%`$ simulation. The conclusion is clear: the system is not stable in the separated state and returns very quickly to the homogeneous state.

In section 3 we have used the calculations on Fe/O crystals to show that phase stability in high-temperature Fe/O systems might well be expected. In particular, we gave an example of an Fe<sub>3</sub>O solid structure whose enthalpy of formation is only $`1`$ eV. We have used our $`25\%`$ simulation to create another such structure. This was done simply by quenching the liquid at the rate of 3000 K ps<sup>-1</sup> until the atoms came to mechanical equilibrium in an amorphous structure (this is clearly a local minimum of the total energy function). We can regard this as a crystal with the Fe<sub>3</sub>O composition having an unusually large supercell. The calculated enthalpy of formation of this solid at the pressure of 290 GPa is reported in the left panel of Fig. 2, and we see that its stability is even greater that that of the BiI<sub>3</sub> structure proposed for Fe<sub>3</sub>O in section 3. This confirms the idea that even more stable structures may yet be found.

### 4.4 Dynamics

In studying the dynamical properties of the Fe/O liquid, our main concern is with the viscosity. However, we first give results for the atomic diffusion coefficients $`D_\alpha `$, which give a simple way of characterising the motion of the atoms. These are straightforwardly calculated from the mean square displacement of the atoms through the Einstein relation (Allen and Tildesley, 1987):

$$\frac{1}{N_\alpha }\underset{i=1}{\overset{N_\alpha }{}}|𝐫_{\alpha i}(t_0+t)𝐫_{\alpha i}(t_0)|^26D_\alpha t,\mathrm{as}t\mathrm{},$$
(8)   
where $`𝐫_{i\alpha }(t)`$ is the position of the $`i`$th atom of species $`\alpha `$ at time $`t`$, $`N_\alpha `$ has its usual meaning, and $``$ is the thermal average, in practice evaluated by averaging over time origins $`t_0`$. In studying the long time behaviour of the mean square displacement, it is convenient to define a time dependent diffusion coefficient $`D_\alpha (t)`$:

$$D_\alpha (t)=\frac{1}{6tN_\alpha }\underset{i=1}{\overset{N_\alpha }{}}|𝐫_{\alpha i}(t_0+t)𝐫_{\alpha i}(t_0)|^2,$$
(9)   
which has the property that

$$\underset{t\mathrm{}}{lim}D_\alpha (t)=D_\alpha .$$
(10)   
In Fig. 9 we display the iron and the oxygen diffusion coefficients calculated using Eq. (9). From this data we estimate $`D_{\mathrm{Fe}}0.8\times 10^8`$ m<sup>2</sup> s<sup>-1</sup> and $`D_\mathrm{O}1.0\times 10^8`$ m<sup>2</sup> s<sup>-1</sup>. These values should be compared with those obtained for pure liquid iron, $`D_{\mathrm{Fe}}0.40.5\times 10^8`$ m<sup>2</sup> s<sup>-1</sup> (Vočadlo et al., 1997; de Wijs, 1998), and those obtained for Fe and S in our Fe/S simulation $`D_{\mathrm{Fe}}0.40.6\times 10^8`$ m<sup>2</sup> s<sup>-1</sup>, and $`D_\mathrm{S}0.40.6\times 10^8`$ m<sup>2</sup> s<sup>-1</sup> (Alfè and Gillan, 1998a). This means that the two species of atoms in liquid Fe/O diffuse somewhat more rapidly than the atoms in both liquid iron and the liquid Fe/S alloy at the same pressure and temperature.

In the past, atomic diffusion coefficients have often been used to estimate the viscosity of liquids via the Stokes-Einstein relation, and this was the procedure used in our previous work on liquid Fe and Fe/S. Since the diffusion coefficients are larger in the present case, this procedure would lead us to expect a lower viscosity for liquid Fe/O. We have recently demonstrated (Alfè and Gillan, 1998b) that the viscosity can be more directly (and more rigorously) calculated in first-principles simulations using the Green-Kubo relations, i.e. the relations between transport coefficients and correlation functions involving fluxes of conserved quantities (Allen and Tildesley, 1987). The shear viscosity $`\eta `$ is given by:

$$\eta =\frac{V}{k_BT}_0^{\mathrm{}}𝑑tP_{xy}(t_0+t)P_{xy}(t_0),$$
(11)   
where $`V`$ is the volume of the system and $`P_{xy}`$ is the off-diagonal component of the stress tensor $`P_{\alpha \beta }`$ ($`\alpha `$ and $`\beta `$ are Cartesian components). The stress tensor is straightforwardly calculated, so that the stress autocorrelation function (SACF) $`P_{xy}(t_0+t)P_{xy}(t_0)`$ can also be obtained, but at first sight it might appear that very long simulations would be needed to gather adequate statistical sampling. However, we have recently shown that perfectly adequate $`\eta `$ values can be obtained with surprisingly short runs, and we have reported results for liquid aluminium and liquid Fe/S (Alfè and Gillan, 1998b).

In the left panel of Fig. 10 we display the the average of the five independent components of the traceless SACF divided by its value at $`t=0`$ which we denote by $`\varphi (t)`$. Since the traceless part of the stress tensor has zero average, $`\varphi (t)`$ goes to zero as $`t\mathrm{}`$. The statistical error on $`\varphi (t)`$ for all values of $`t`$ is $`5\%`$ of the value at $`t=0`$, and after $`0.20.3`$ ps the magnitude of $`\varphi (t)`$ falls below that error. In the right panel of Fig. 10 we display the integral $`_0^t𝑑t^{}\varphi (t^{})`$ of $`\varphi (t)`$ as a function of time. The limiting value of the integral for $`t\mathrm{}`$ is the shear viscosity. The error that one makes in evaluating that integral grows with time, since one integrates the noise together with $`\varphi (t)`$. We estimated the error in the integral as a function of time using the scatter of the SACF’s. Combining this estimate with an analytic expression for the error, we obtain the error estimate displayed in Fig. 10. From the point where $`\varphi (t)`$ falls below the noise one integrates only the latter, so one gains nothing by evaluating the integral beyond that point. If we assume that $`\varphi (t)`$ is zero above $`t0.3`$ ps, we obtain the value $`\eta =4.5\pm 1.0`$ mPa s. This value is roughly half the value reported earlier for liquid Fe (Vočadlo et al., 1997) and Fe/S (Alfè and Gillan, 1998a) at ICB conditions, and is not very much greater than the viscosity of typical liquid metals under ambient conditions; for example, the viscosity of liquid aluminium at atmospheric pressure 100 K above its melting point is 1.25 mPa s (Shimoji and Itami, 1986).

This result confirms our earlier conclusion (de Wijs et al., 1998) that the viscosity of the outer core is towards the lower end of the wide range of theoretical and experimental values reported in the literature.

### 4.5 Electronic structure

We have studied the electronic structure of our simulated Fe/O liquid in order to shed light on the nature of the bonding between the atoms and to help to interpret the structure discussed in section 4.2. The main tools used in this analysis are the electronic density of states (DOS) and the local density of states (LDOS). The DOS represents the number of electronic states per unit energy as a function of energy, while the LDOS is a projection of the DOS onto states of chosen angular momentum on atoms of chosen species. In performing this projection we took spherical regions of radius $`R`$ on the atoms, and in practice we chose $`R=0.6`$ Å for both iron and oxygen. This $`R`$ value is considerably smaller than half the interatomic distance in the liquid (see Fig. 3), so we expect to distinguish clearly between the electronic structures on different atoms. The results are not averaged over time but are calculated from the electronic energies and wavefunctions for selected time steps taken from the $`25\%`$ simulation.

The calculated DOS is shown in the upper panel of Fig. 11, and consists of four main features (energies are referred to the Fermi energy): two fairly narrow peaks at $`24`$ and $`11`$ eV; a large dominant peak spanning the range $`9`$ eV to $`3`$ eV; and a broad feature extending well above the Fermi energy. The LDOS shown in the lower panel of the Figure allows us to interpret these features. The peak at $`24`$ eV consists entirely of O$`(2s)`$ states. The peak at $`11`$ eV is mainly O$`(2p)`$, but with an appreciable contribution of Fe$`(3d)`$. The dominant peak from $`9`$ eV to $`3`$ eV consists mainly of Fe$`(3d)`$ states, but there is a significant peak associated with O$`(2p)`$ states just above the Fermi energy. The broad feature above the Fermi energy comes from both Fe$`(3d)`$ and O$`(2p)`$ states.

To see the implications of this structure for the interatomic bonding, we note that if the O$`(2p)`$ states were much lower in energy than the Fe$`(3d)`$ states, there would be very little hybridisation between the two kind of states, the O$`(2p)`$ levels would be completely filled, and the oxygen atoms would carry a net charge of $`2|e|`$. The Fe-O bonding would be ionic. On the other hand, if the O$`(2p)`$ states were at the same energy as Fe$`(3d)`$, we should expect strong hybridisation and little charge transfer, so that the bonding would be covalent. It is clear from Fig. 11 that the Fe-O bonding is intermediate between ionic and covalent. The O$`(2p)`$ states are below Fe$`(3d)`$, but not low enough to suppress hybridisation. There is a clear splitting of the O$`(2p)`$ levels into bonding and anti-bonding orbitals, though the bonding orbital clearly has much larger weight (in pure covalent bonding we should expect the weights to be equal). The peak in Fe$`(3d)`$ at $`11`$ eV is also clear evidence for hybridisation between O$`(2p)`$ and Fe$`(3d)`$. The implication is that there is a partial charge transfer from Fe to O, but not enough to give oxygen a charge of $`2|e|`$.

The bonding between Fe atoms is metallic. The partial filling of the Fe$`(3d)`$ levels gives the well known bonding mechanism emphasised by Friedel’s (1969) analysis of the cohesive and elastic properties of transition metal crystals.

Since the oxygen atoms carry partial charges and their $`2p`$ orbitals are almost full, no covalent bond is expected between them, and this is also clear from Fig. 3. To investigate the electronic state of oxygen in more detail, we have calculated the LDOS for oxygen atoms in different environments. We show in Fig. 12 the LDOS for two oxygen atoms denoted O<sub>a</sub> and O<sub>b</sub>, which have been chosen from the $`25\%`$ simulation so that O<sub>a</sub> has 10 Fe neighbours and one O neighbour while O<sub>b</sub> has seven Fe and four O neighbours. If there were any covalent bonding between O atoms, we should expect a larger bonding-antibonding splitting of the $`2p`$ states for the O<sub>b</sub> atoms than for O<sub>a</sub>, and we might expect a similar splitting (or at least a broadening) of the O$`(2s)`$ states for O<sub>b</sub>. The LDOS curves show neither of these effects. Instead, the main difference is the upward shift of the peaks for O<sub>b</sub> compared with O<sub>a</sub>. We believe this is direct evidence for the partial charge transfer: the valence electrons of O<sub>b</sub> feel a repulsive electrostatic potential due to the partial negative charges on the four oxygen neighbours, which raises their energy.

The main bonding mechanisms in the liquid are therefore the ionic-covalent Fe-O bond and the metallic Fe-Fe bond. Since the O$`(2p)`$ atomic orbitals are more compact than Fe$`(3d)`$ orbitals we expect the Fe-O distance to be shorter than the Fe-Fe distance, and this effect is clear from the rdf’s shown in Fig. 3. The partial charge transfer presumably also contributes to the shortening of this distance.

## 5 Discussion

Two of our main aims in this paper have been to probe the phase stability of liquid Fe/O under Earth’s core conditions, and to determine the oxygen concentration that would be needed to reproduce the known core density. In practice, these aims must be taken together: we want to know whether the alloy is thermodynamically stable at the appropriate concentration.

The results we have presented leave little doubt that if the liquid is stable then the mole fraction of oxygen must be in the region $`2530\%`$ (our best estimate is $`28\%`$), because anything much less than this would give a pressure that is too low at the known density. This is essentially the same as the value proposed many years ago by Ringwood (1977). In judging the robustness of this conclusion we recall some important facts: First, DFT electronic-structure methods of the type used here generally predict the density of materials at a given pressure to within a few percent. Particularly relevant here are recent DFT calculations (Stixrude et al., 1994; Söderlind et al., 1996; Vočadlo et al., 1997) on h.c.p. iron over the pressure range $`0350`$ GPa which are in excellent agreement with each other and with the experimental results. Similar comparisons for FeSi (Vočadlo et al., 1997) are also relevant. DFT calculations on oxides (including transition-metal oxides) generally predict the density with similar accuracy. A second important fact is that our earlier first principles simulations of pure liquid iron, based on exactly the same techniques, gave a prediction of the density at the pressure of the inner core boundary which is correct to $`2\%`$. (In fact the comparison was done the other way round: at the density of 13300 kg m<sup>-3</sup> and the temperature of 6000 K, our simulations gave a pressure of 358 GPa, compared with the value of 330 GPa estimated from experimental data.) Third, our estimate of the oxygen concentration is based on changes of pressure and density compared with our simulated pure iron. The calculations should be even more reliable for these differences then they are for the absolute values. We therefore believe that our value for the oxygen concentration required should be subject to an error of no more that $`5\%`$.

It is interesting to compare the liquid composition with that for the case of Fe/S. In our recent simulations (Alfè and Gillan, 1998a), we showed that liquid Fe/S is not far from reproducing the known pressure and density at the inner core boundary with a sulphur mole fraction of $`18\%`$. At this composition, the mass density of 12330 kg m<sup>-3</sup> gave a calculated pressure of $`349\pm 6`$ GPa. If we make a correction to bring the density to 12000 kg m<sup>-3</sup>, we find a sulphur mole fraction of $`23\%`$. This means that to achieve the required reduction in density we actually need a higher mole fraction of oxygen than of sulphur, even though the atomic mass of oxygen is only half that of sulphur. The reason is that the oxygen atom is smaller, a point to which we return below.

The question of phase stability is more complex. What seems certain is that earlier arguments against phase stability based on ab initio calculations of the energetics of Fe<sub>3</sub>O and Fe<sub>4</sub>O crystals (Sherman, 1995) do not really deliver the intended conclusion. This is not because the calculations were wrong. On the contrary, our calculations fully support their correctness. It is simply that the crystal structures assumed were not the most stable. We have presented an alternative structure for Fe<sub>3</sub>O which gives a formation enthalpy that is low enough ($`1`$ eV) to make phase stability under core conditions quite plausible.

We have used our simulations to probe the phase stability of the Fe/O liquid in the appropriate region of concentration, and all the indications are that it is stable. We therefore fully support the conclusion that has been drawn from high pressure experimental measurements (Ohtani et al., 1984; Ringwood and Hibberson, 1990) that liquid Fe and FeO are miscible under core conditions in the concentration region of interest. However, a word of caution is in order. The fact is that our simulated systems are rather small, and it is quite conceivable that a system that would be unstable in the thermodynamic limit could be stabilised by the artificial periodic boundary conditions in small simulation cells. Nevertheless, our calculations certainly provide strong support for the thermodynamic stability of liquid Fe/O under core conditions at the relevant concentration. In the end, we believe that the definitive theoretical approach to this question will be first principles free energy calculations on the appropriate solid and liquid phases. These would be computationally very demanding, but should certainly be feasible in the near future. (First principles free energy calculations have recently been used with success to calculate the melting properties of silicon and aluminium (Sugino and Car, 1995; de Wijs et al., 1998).)

The small size of oxygen is clear from our analysis of the liquid structure: the Fe-O nearest neighbour distance (the position of the first peak in the rdf) is only at $`1.7`$ Å, compared with $`2.1`$ Å for the Fe-Fe distance. We recall that in the Fe/S liquid, the Fe-S distance was $`1.95`$ Å (Alfè and Gillan, 1998a). An important feature of the liquid is that each oxygen has on average only 9 iron neighbours, whereas it has between 4 and 5 oxygen neighbours. This atomic environment of oxygen is very different from that produced by the cubic structure of Fe<sub>3</sub>O. By contrast in the BiI<sub>3</sub> structure of Fe<sub>3</sub>O oxygen has 11 iron neighbours and 2 oxygen neighbours. The dilemma in making crystal structures for Fe/O solid solutions at core pressures is that we want to achieve close packing because of the high pressure, but the atoms that are being packed have different sizes. It seems that the cubic structure of Fe<sub>3</sub>O is not a good solution. The BiI<sub>3</sub> structure is better, even though it means putting oxygen atoms next each other. There may be yet better ways.

Our analysis of the electronic structure of the liquid gives further insight. Here, the important feature is that Fe-O bonding is only partially ionic, with a substantial covalent contribution. The implication that there is partial electron transfer from Fe to O is relevant, because presumably if O carried a full ionic charge of $`2|e|`$ oxygen atoms would be more reluctant to become neighbours of each other. The electronic structure shows that there is no detectable covalent interaction between oxygens, so the fact that they become neighbours cannot be attributed to covalency.

Finally, we have studied the diffusion coefficient of iron and oxygen atoms and the viscosity of the liquid at $`33\%`$ composition. The finding is that both species diffuse more rapidly than in either pure liquid iron or the liquid Fe/S alloy: the diffusion coefficients have roughly twice the value that they have in pure Fe and Fe/S at the same pressure and temperature. We also find that the viscosity of the Fe/O liquid is about half what it is in those other systems. This means that if the major light element in the outer core was oxygen, our earlier conclusion (de Wijs et al., 1998) about the low viscosity of the outer core would be confirmed, and indeed strengthened.

## 6 Conclusions

We are led to the following conclusions: if the light impurity in the outer core is mainly oxygen, then its molar concentration would have to be $`28\%`$; in this region of concentration, we have strong evidence that the liquid is stable against phase separation; the proposed miscibility is not in conflict with the large formation enthalpies predicted by earlier ab initio calculations, because those calculations were based on assumed structures that are not the most stable; the proposed Fe/O liquid alloy has an even lower viscosity than that of pure Fe and the relevant Fe/S alloy under the same conditions.

## Acknowledgments

The work of DA is supported by NERC grant GST/O2/1454 to G. D. Price and M. J. Gillan. We thank the High Performance Computing Initiative for allocations of time on the Cray T3D and T3E at Edinburgh Parallel Computer Centre, these allocations being provided through the Minerals Physics Consortium (GST/02/1002) and the U.K. Car-Parrinello Consortium. We thank Dr. G. Kresse and Dr. G. de Wijs for valuable technical assistance, and Dr. L. Vočadlo and Dr. D. Sherman for useful discussions. We thank Dr. W. Temmerman for useful advice about the LMTO-46 code, and the Psi-k Network for making the code available.

## References

Alder, B. J., 1966. Is the mantle soluble in the core? J. Geophys. Res., 71: 4973-4979.

Alfè, D., and Gillan, M. J., 1998a. First principles simulations of liquid Fe-S under Earth’s core conditions. Phys. Rev. B, in press.

Alfè, D., and Gillan, M. J., 1998b. The first principles calculation of transport coefficient. Phys. Rev. Lett., submitted.

Allen, M. P., Tildesley, D. J., 1987. Computer Simulation of Liquids. Oxford University Press, Oxford.

Bhatia, A. B., Thornton, D. E., 1970. Structural aspects of the electrical resistivity of binary alloys. Phys. Rev. B, 2: 3004-3012.

Birch, F., 1952. Elasticity and composition of the Earth’s interior. J. Geophys. Res., 57: 227-286.

Birch, F., 1964. Density and composition of mantle and core. J. Geophys. Res., 69: 4377-4388.

Boehler, R., 1993. Temperature in the Earth’s core from melting-point measurements of iron at high static pressures. Nature. 363: 534-536.

Bullen, K. E., 1973. Cores of terrestrial planets. Nature, 243: 68-70.

Car, R., and Parrinello, M., 1985. Unified approach for molecular-dynamics and density-functional theory. Phys. Rev. Lett., 55: 2471-2474.

Clark, S. P., 1963. Variation of density in the Earth and the melting curve in the mantle. In: The Earth Sciences, edited by T. W. Donnelly, University of Chicago Press, Chicago.

Distin, P. A., Whiteway S., and Masson, C., 1971. Solubility of oxygen in liquid iron from 1785<sup>o</sup> to 1960<sup>o</sup>. A new technique for the study of slag-metal equilibria. Canad. Metallurg. Quart., 10: 13-18.

Dubrovskiy, V. A. and Pan’kov, V. L., 1972. On the composition of the Earth’s core. Acad. Sci. USSR Phys. Solid Earth, 7: 452-455.

Fischer, W. A., and Schumacher, J. F., 1978. Die Sättigungslöslichkeit von Reineisen an Saurstoff vom Schmelzpunkt bis 2046<sup>o</sup> C, ermittelt mit dem Schwebeschmelzverfahren. Arch. Eisenhüttenwes., 49: 431-435.

Friedel, J., 1969. Fundamental aspects in the bonding of transition metals. In: The Physics of Metals, ed. J. Ziman, Cambridge University Press, Cambridge, 494 pp.

Fukai, Y., and Akimoto, S., 1983. Hydrogen in the Earth’s core. Proc. Jpn. Acad., Ser. B, 59: 158-162.

Gillan, M. J., 1997. The virtual matter laboratory. Contemp. Phys., 38: 115-130.

Goarant, F., Guyot, F., Peyronneau, J., and Poirier, J.-P., 1992. High-pressure and high-temperature reactions between silicates and liquid iron alloys in the diamond anvil cell, studied by analytical electron microscopy. J. Geophys. Res., 97: 4477-4487.

Hohenberg, P., and Kohn, W., 1964. Inhomogeneous electron gas. Phys. Rev. 136: B864-B871.

Holender, J. M., and Gillan, M. J., 1996. Composition dependence of the structure and electronic properties of liquid Ga-Se alloys studied by ab initio molecular dynamics simulation. Phys. Rev. B, 53: 4399.

Jones, R. O., and Gunnarsson, O., 1989. The density functional formalism, its applications and prospects. Rev. Mod. Phys., 61: 689-746.

Kirchhoff, F., Holender, J. M., and Gillan, M. J., 1996a. The structure, dynamics and electronic structure of liquid Ag-Se alloys investigated by ab initio simulation. Phys. Rev. B, 54: 190-202.

Kirchhoff, F., Gillan, M. J., Holender, J. M., Kresse, G., and Hafner, J., 1996b. Structure and bonding in liquid Se. J. Phys. Condens. Matter, 8: 9353-9357.

Knittle, E., Jealnoz, R., 1991. Earth’s core-mantle boundary: results of experiments at high pressures and temperatures. Science, 251: 1438-1443.

Kohn, W., and Sham, L., 1965. Self-consistent equations including exchange and correlation effects. Phys. Rev. 140: A1133-A1138.

Kresse, G., and Hafner, J., 1993. Ab-initio molecular-dynamics for open-shell transition-metals. Phys. Rev. B, 48: 13115-13118.

Kresse, G., and Furthmüller, J., 1996a. Efficiency of ab-initio total-energy calculations for metals and semiconductors using a plane-wave basis-set. Comput. Mater. Sci., 6: 15-50.

Kresse, G., and Furthmüller, J., 1996b. Efficient iterative schemes for ab-initio total-energy calculations using a plane-wave basis-set. Phys. Rev. B, 54: 11169-11186.

Krier, G., Jepsen, O., Burkhardt, A., and Andersen, O. K., 1994. LMTO-46 code, Max-Planck-Institut für Festkörperforschung, Stuttgart.

Lewis, J. S., 1973. Chemistry of the planets. Annu. Rev. Phys. Chem., 24: 339-351.

Louie, S. G., Froyen, S., and Cohen, M. L., 1982. Non-linear ionic pseudopotentials in spin-density-functional calculations. Phys. Rev. B, 26: 1738-1742.

MacDonald, G. J. F., and Knopoff, L., 1958. The chemical composition of the outer core. J. Geophys., 1: 1751-1756.

Mao, H. K., Wu, Y., Chen, L. C., Shu, J. F., and Jephcoat, A. P., 1990. Static compression of iron to 300 GPa and Fe0.8Ni0.2 alloy to 260 GPa: implications for composition of the core. J. Geophys. Res., 95: 21737-21742.

Mason, B., 1966. Composition of the Earth. Nature, 211: 616-618.

Monkhorst, H. J., and Pack, J. D., 1976. Special points for Brillouin-zone integrations. Phys. Rev. B, 13: 5188-5192.

Murthy, V., and Hall, H. T., 1970. On the possible presence of sulfur in the Earth’s core. Phys. Earth Planet. Inter., 2: 276-282.

Nosé, S., 1984. A molecular dynamics method for simulations in the canonical ensemble. Molec. Phys., 52: 255-268.

Ohtani, E., and Ringwood, A. E., 1984 Composition of the core, I. Solubility of oxygen in molten iron at high temperatures. Earth Planet. Sci. Lett., 71: 85-93.

Ohtani, E., Ringwood, A. E. and Hibberson, W., 1984. Composition of the core, II. Effect of high pressure on solubility of FeO in molten iron. Earth Planet. Sci. Lett., 71: 94-103.

Parr, R. G., and Yang, W., 1989. Density-Functional Theory of Atoms and Molecules, Oxford University Press, Oxford.

Poirier, J.-P., 1991. Introduction to the Physics of the Earth’s Interior. Cambridge University Press, Cambridge.

Poirier, J.-P., 1994. Light elements in the Earth’s outer core: A critical review. Phys. Earth Planet. Interiors. 85: 319-337.

Pulay, P., 1980. Convergence acceleration of iterative sequences. The case of SCF iteration. Chem. Phys. Lett., 73: 393.

Ringwood, A. E., 1959. On the chemical evolution and densities of the planets. Geochem. J., 15: 257-283.

Ringwood, A. E., 1961. Silicon in the metal phase of enstatite chondrites and some geochemical implications. Geochim. Cosmochim. Acta, 25: 1-13.

Ringwood, A. E., 1966. Chemical evolution of the terrestrial planets. Geochim. Cosmochim. Acta, 30: 41-104.

Ringwood, A. E., 1977. On the composition of the core and implications for the origin of the Earth. Geochim. Cosmochim. Acta, 11: 111-135.

Ringwood, A. E., and Hibberson, W., 1990. The system Fe-FeO revisited. Phys. Chem. Minerals, 17: 313-319.

Secco, R. A. 1995. Viscosity of the Outer Core. In: Mineral Physics and Crystallography: A Handbook of Physical Constants, ed. T. J. Ahrens, Americal Geophysical Union, 218 pp.

Sherman, D. M., 1995. Stability of possible Fe-FeS and Fe-FeO alloy phases at high pressure and the composition of the Earth’s core. Earth Planet. Sci. Lett., 132: 87-98.

Shimoji, M., and Itami, T., 1986. Atomic Transport in Liquid Metals, Trans Tech Publications, Aedermannsdorf, p. 191.

Söderlind, P., Moriarty, J. A., and Wills, J. M., 1996. First principles theory of iron up to the Earth’s core pressures: Structural, vibrational, and elastic properties. Phys. Rev. B, 53: 14063-14072.

Štich, I., Car, R., and Parrinello, M., 1989. Bonding and disorder in liquid silicon. Phys. Rev. Lett., 63: 2240-2243.

Stixrude, L., Cohen, R. E., and Singh, D. J., 1994. Iron at high pressure: Linearized-augmented-plane-wave computation in the generalized-gradient approximation. Phys. Rev. B, 50: 6442-6445.

Sugino, O., and Car, R., 1995. Ab-initio molecular-dynamics study of first-order phase-transitions. Melting of silicon. Phys. Rev. Lett., 74: 1823-1826.

Suzuki, T., Akimoto S., and Yagi, T., 1989. Metal-silicate-water reaction under high pressure, I, Formation of metal hydride and implications for composition of the core and mantle. Phys. Earth Planet. Inter., 56: 377-388.

Urey, H. C., 1960. On the chemical evolution and densities of the planets. Geochim. Cosmochim. Acta, 18: 151-153.

Vanderbilt, D., 1990. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism. Phys. Rev. B, 41: 7892-7895.

Verlet, L., 1967. Computer ‘experiments’ on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. Phys. Rev., 159: 98-103.

Vočadlo, L., de Wijs, G. A., Kresse, G., Gillan, M. J., and Price, G. D., 1997. First principles calculations on crystalline and liquid iron at earth’s core conditions. Faraday Discuss., 106: 205-217.

Wang, Y., and Perdew, J., 1991. Correlation hole of the spin-polarized electron gas, with exact small-wave-vector and high-density scaling. Phys. Rev. B, 44: 13298-13307.

de Wijs, G. A., Kresse, G., and Gillan, M. J., 1998. First order phase transitions by first-principles free energy calculations: The melting of Al. Phys. Rev. B, 57: 8223-8234.

de Wijs, G. A., Kresse, G., Vočadlo, L., Dobson, D., Alfè, D., Gillan, M. J., Price, G. D., 1998. The viscosity of liquid iron at the physical conditions of the Earth’s core. Nature, 392: 805-807.

Wood, B. J. 1993. Carbon in the core. Earth Planet. Sci. Lett., 117: 593-607.

Wyckoff, R. W. G., 1964. Crystal Structures, 2nd edition, Vol. 2, Interscience, New York. 

|  |  | PP(spin-unrestricted) | PP(spin-restricted) | FLAPW |
| --- | --- | --- | --- | --- |
| Fe($`ϵ`$) | $`\rho _0`$ | 8910 | 8910 | 8780 |
|  | $`K`$ | 283 | 283 | 260 |
|  | $`K^{}`$ | 4.39 | 4.39 | 4.53 |
| FeO(B8) | $`\rho _0`$ | 5650 | 6980 | 5810 |
|  | $`K`$ | 92 | 258 | 173 |
|  | $`K^{}`$ | 4.96 | 4.4 | 2.93 |
|  | $`\mu `$ | 2.0 |  | not reported |
| Fe<sub>3</sub>O(cubic) | $`\rho _0`$ | 6500 | 7550 | 7420 |
|  | $`K`$ | 123 | 226 | 223 |
|  | $`K^{}`$ | 4.23 | 4.04 | 4.02 |
|  | $`\mu `$ | 2.34 |  | small |
| Fe<sub>4</sub>O | $`\rho _0`$ | 6880 | 7840 | 7830 |
|  | $`K`$ | 135 | 273 | 310 |
|  | $`K^{}`$ | 4.82 | 4.26 | 4.17 |
|  | $`\mu `$ | 2.0 |  |  |
| Fe<sub>3</sub>O(BiI<sub>3</sub>) | $`\rho _0`$ |  | 7930 |  |
|  | $`K`$ |  | 248 |  |
|  | $`K^{}`$ |  | 4.29 |  |

Table 1: Calculated equilibrium properties of crystals in the Fe/O system: equilibrium mass density $`\rho _0`$ (kg m<sup>-3</sup>), bulk modulus $`K`$ (GPa), the pressure derivative $`K^{}=dK/dP`$, and magnetic moment per atom $`\mu `$ (units of Bohr magneton). Results are given for the present spin unrestricted and restricted pseudopotential (PP) calculations and the FLAPW calculations of Sherman (1995). The structures of the cubic Fe<sub>3</sub>O, Fe<sub>4</sub>O and BiI<sub>3</sub>-structure Fe<sub>3</sub>O are shown in Fig. 1.

| Volume (Å<sup>3</sup>) | Pressure (GPa) | $`b/a`$ | $`c/a`$ | $`\alpha `$ | $`\beta `$ | $`\gamma `$ |
| --- | --- | --- | --- | --- | --- | --- |
| 22 | 467 | 1.01016 | 0.95443 | 123.76 | 90.70 | 79.04 |
| 24 | 325 | 1.00988 | 0.95322 | 123.82 | 90.88 | 78.79 |
| 26 | 226 | 1.01010 | 0.95164 | 123.85 | 91.03 | 78.55 |
| 28 | 156 | 1.00772 | 0.94984 | 123.89 | 91.26 | 78.28 |
| 30 | 105 | 1.00574 | 0.94874 | 123.85 | 91.39 | 78.07 |
| 32 | 67 | 0.99711 | 0.94465 | 123.82 | 91.79 | 77.69 |
| 34 | 40 | 0.99452 | 0.94303 | 123.76 | 91.97 | 77.37 |
| 36 | 19 | 0.99193 | 0.94258 | 123.57 | 91.98 | 76.95 |

Table 2: Cell parameters and pressure $`P`$ of Fe<sub>3</sub>O in the BiI<sub>3</sub> structure (see Fig. 1) calculated at a series of volumes (per Fe<sub>3</sub>O unit). The quantities $`a`$, $`b`$, $`c`$ are the magnitudes of the primitive translation vectors, and $`\alpha ,\beta ,\gamma `$ are the angles between the pairs $`(a,b),(a,c)`$ and $`(b,c)`$ respectively.

## Figure captions

Fig. 1: The crystal structures of cubic Fe<sub>3</sub>O (left), cubic Fe<sub>4</sub>O (centre) and the BiI<sub>3</sub> form of Fe<sub>3</sub>O (right) used to calculate the formation enthalpies of Fe/O solid solutions. Large and small spheres represent iron and oxygen respectively.

Fig. 2: Calculated enthalpies of formation (per formula unit) of solid solutions having compositions Fe<sub>3</sub>O and Fe<sub>4</sub>O. Left panel shows spin-restricted results from present work compared with FLAPW results of Sherman (1995); right panel shows present spin-unrestricted results. Key to style of curves: present cubic Fe<sub>3</sub>O  , cubic Fe<sub>3</sub>O of Sherman (1995)  , present cubic Fe<sub>4</sub>O      , cubic Fe<sub>4</sub>O of Sherman (1995)      , present Fe<sub>3</sub>O in the BiI<sub>3</sub> structure            . The isolated point shows formation enthalpy of the amorphous structure obtained by quenching the liquid (see Sec. 4).

Fig. 3: Radial distribution functions $`g_{\alpha \beta }(r)`$ obtained from simulation of liquid Fe/O at oxygen molar concentration of $`33\%`$.

Fig. 4: The iron-iron radial distribution function $`g_{\mathrm{FeFe}}(r)`$ from the present simulation of liquid Fe/O (oxygen molar concentration of $`33\%`$) compared with $`g_{\mathrm{FeFe}}(r)`$ from simulation of pure liquid iron at similar pressure and temperature (Vočadlo et al., 1997).

Fig. 5: Partial structure factors $`S_{\alpha \beta }(k)`$ calculated from simulation of liquid Fe/O at oxygen molar concentration of $`33\%`$.

Fig. 6: Probability distribution $`P_{\mathrm{FeO}}(n,r_{\mathrm{FeO}}^c)`$ for number $`n`$ of oxygen neighbours of an iron atom calculated from simulation of liquid Fe/O at oxygen molar concentration of $`33\%`$.

Fig. 7: Snapshots of the simulated liquid Fe/O system along three principal Cartesian axes. Top three panels show a configuration from the period when phase separation was artificially induced by application of an external force. Bottom three panels show a configuration from the later period after removal of the external force. Large and small spheres represent iron and oxygen respectively.

Fig. 8: Probability distribution $`P_{\mathrm{FeO}}(n,r_{\mathrm{FeO}}^c)`$ for number $`n`$ of oxygen neighbours of an iron atom calculated from simulation of liquid Fe/O at oxygen molar concentration of $`25\%`$. Results are average values for three windows of length 0.1 ps at times of 0 ps, 0.5 ps and 1 ps after removal of the external force used to induce phase separation.

Fig. 9: Time dependent diffusion coefficients $`D_\alpha (t)`$ for iron and oxygen calculated from the simulation of liquid Fe/O at oxygen molar concentration of $`33\%`$.

Fig. 10: Average over the five independent components of the autocorrelation function of the traceless stress tensor $`\varphi (t)`$ (left panel) and viscosity integral (solid curve) with its statistical error (dashed curve) (right panel).

Fig. 11: Electronic density of states (upper panel) and local densities of states (lower panel) calculated for liquid Fe/O at oxygen molar concentration of $`25\%`$. Energy is referred to the Fermi energy $`E_f`$.

Fig. 12: Local densities of states for two selected oxygen atoms taken from the simulation of liquid Fe/O at oxygen molar concentration of $`25\%`$. Atom O<sub>a</sub> has 1 oxygen neighbour and 10 iron neighbours; atom O<sub>b</sub> has 4 oxygen neighbours and 7 iron neighbours. Energy is referred to the Fermi energy $`E_f`$.

FIGURE 1

FIGURE 2

FIGURE 3

FIGURE 4

FIGURE 5

FIGURE 6

FIGURE 7

FIGURE 8

FIGURE 9

FIGURE 10

FIGURE 11

FIGURE 12