Introduction

Molten iron oxide plays a crucial role in iron, steel, and copper production1,2 and has a major role in determining the behavior and properties of magmatic liquids3,4,5,6. Iron redox has previously been studied in multi-component silicate liquids and glasses3,7,8, where the oxidation states and coordination numbers have a direct effect on melt viscosity, density, phase stability, and heat capacity7. However, the assignment of the local Fe coordination environment and its relationship to the oxidation state is not straightforward. This is because iron redox depends on many factors, including oxygen partial pressure, temperature, and composition. Wilke7 has pointed out that the Fe coordination by oxygen in widely studied glass analogues is likely to differ from those found in high temperature melts. Consequently, the structural role of Fe is poorly understood due to the difficulty of studying liquids at high temperatures using direct methods9. Under oxidizing conditions, no suitable containers exist, and oxygen pressures of several atmospheres or supercooling are essential to stabilize molten Fe2O3. To our knowledge, the only structural measurements on a pure iron oxide melt are two 40-year-old X-ray diffraction studies by Waseda et al.10,11 at compositions close to FeO, and two X-ray absorption near edge studies12,13.

Structural interpretation of experimental data is complicated by the fact that both ferrous (Fe2+) and ferric (Fe3+) iron co-exist in melts, and both can be in a range of coordination environments. Fe2+ is the most abundant iron species in magmas and is usually 6-fold in crystalline minerals at ambient pressure. In the liquid state however, it has been debated if Fe2+ is primarily 5- or 6-fold coordinated by oxygen14,15 or 4-fold9,16. Historically Fe2+ was considered as being distributed between sites coordinated by either 4 or 6 oxygen atoms based on crystal structures16,17. However, it has also been suggested that Fe in the melt, similar to Al, can also be 5-fold coordinated by oxygen8,18. This has been supported by the fact that in some minerals Fe-O polyhedra exist with 5-short bonds5. A bimodal distribution between 4- and 6-fold Fe2+ in molten FeO would be in contrast to other mixed coordination oxide melts, such as Al2O3, TiO2 or MgO, in which liquids and glasses contain a single distribution of metal-oxygen bond distances, representing a mixture of 4, 5, and 6 coordinated units19,20. The coordination environment of Fe3+ is similarly unclear. In crystalline Fe3O4, ε-Fe2O3, and γ-Fe2O3, Fe atoms reside mainly in octahedral coordination with some tetrahedral sites. Contrastingly, in a soda-lime-silica glass Fe3+ has been considered to be tetrahedrally surrounded by oxygen atoms, along with minority higher coordination environments21, while measurements on Fe3+ in aluminosilicates indicate primarily 5- or 6-fold coordination in the liquid state15,16, implying considerable compositional dependence.

Here we have determined the structures of liquid FeOx (1.0 ≲ x ≲ 1.5) over a wide range of oxidation states and temperatures, using a combination of aerodynamic levitation with laser heating in various gas mixtures and high-energy X-ray diffraction. Levitation methods enable access to the oxidizing and metastable regions of the Fe-O phase diagram without interference from crucible materials, as well as the ability to supercool the liquid by the elimination of heterogeneous crystal nucleation caused by a container. Atomistic models have been obtained using Empirical Potential Structure Refinement (EPSR) modeling, see Table 1. The local Fe-O coordination number in FeOx melts is found to have majority populations of 4- and 5-fold polyhedra across a wide range of oxidation states, with structures that are very different from the crystalline phases. At the same time, the new data highlights gaps in our understanding of the structural role of iron, even in the pure FeOx system.

Table. 1 Parameters used for the Lennard-Jones reference potential in the EPSR simulations taken from the literature40 for Fe3+ and O2−, and modified for Fe2+.
Full size table

Results

The FeOx phase diagram

The locations of our liquid FeOx X-ray measurements within the FeO-Fe2O3 region of the phase diagram22 are shown in Fig. 1. At temperatures of ∼1673 K and below, the more oxidized liquids are deeply supercooled into the Wustite-Magnetite-Hematite phase fields. The x-ray data at ∼1973 K traverse the entire stable single-phase oxide liquid region. Our most reducing gas atmosphere (99CO:1CO2) is expected to bring the sample within the large miscibility gap between the liquid oxide and metal, however, buoyancy effects lead to measurement of the oxide melt portion only, see Supplementary Figs. S1 and S2 in the high energy X-ray experiments section for a detailed discussion.

Fig. 1: The FeO-Fe2O3 region of the equilibrium phase diagram22.
figure 1

The crystalline regions are in white, the single-phase liquid regions in grey and two liquid region in blue. Dashed lines are oxygen isobars labelled with log10 p(O2) (atm.). Circles and diamonds represent measurement conditions in this study and the star is the data point from Waseda et al.10. Our measurements in O2 gas are plotted to the right of the 100 = 1 atm. O2 line, because equation 12 of ref. 35 used to calculate Fe3+/ΣFe differs slightly from the most recent phase diagram22 (see uncertainty discussion and Supplementary Fig. S9). For measurements taken in the most reducing gas mixture (99CO:1CO2) the dome boundary between the oxide liquid and mixed metal + oxide region was used to determine the sample compositions (pink triangles).

Full size image

X-ray pair distribution functions

We have considered the X-ray pair distribution functions of FeOx melts, to be composed of six partial pair distribution functions: Fe2+-Fe2+, Fe2+-Fe3+, Fe3+-Fe3+, Fe2+-O, Fe3+-O, and O-O correlations. In the D(r) of Fig. 2, the intensity of the first peak representing the Fe-O bonds decreases from the most oxidizing to the most reducing atmosphere, and the peak position increases from 1.92(1) Å to 1.95(1) Å (Fig. 2). This indicates that on average the Fe2+-O bond lengths are longer (and weaker) than Fe3+-O. Concomitantly, the second peak in D(r) shifts from 3.34(3) Å to 3.15(2) Å with a decrease of intensity, which is attributed to changes in the relative amounts of Fe2+-Fe2+, Fe2+-Fe3+ and Fe3+-Fe3+ neighbors which each have different pair distances. In evaluating the local Fe-O coordination numbers, the method of choice of cut-off distance was typically found to shift the absolute values of nFeO by ±0.5, but did not affect the trend observed and plotted in Fig. 3. Here we define three criteria for the cut-off distance (i) the first minimum in gij(r) (ii) first minimum in r2gij(r) and (iii) a fixed value, where the Fe2+-O and Fe3+-O radial cut-off distances were chosen separately, see Supplementary Figs. S3–S5 in the extraction of coordination numbers section. The Fe-O coordination number was obtained using all three methods, and the average nFeO, was determined as a function of Fe3+/ΣFe, and is shown in Fig. 3a.

Fig. 2: High energy X-ray diffraction data on liquid FeOX.
figure 2

a Q-multiplied X-ray structure factors Q[S(Q)-1] and corresponding b differential distribution functions D(r) for liquid FeOx measured at ∼1673 K in different gases. Comparison is made to the digitized x-ray S(Q) from Waseda et al.10 and our Fourier transform of this data (dotted lines).

Full size image
Fig. 3: Coordination numbers from the empirical potential structure refinement fits.
figure 3

a Average Fe-O coordination number nFeO in molten FeOx at 1973 K, 1873 K, 1773 K, 1673 K, and 1573 K (colored symbols) from our EPSR models, as a function of Fe3+/ΣFe. nFeO was determined by the average of three cut-off criteria. The solid gray stars represent crystalline iron oxides. nFeO obtained by reanalyzing Waseda’s data10 is shown as the open orange star. nFeO in molten Fe2SiO4+x25 is indicated by the dashed line, as derived from peak fitting to X-ray diffraction data, which may lead to a smaller effective radial cut-off than used herein. b Average nFe2+O (open squares) and nFe3+O (solid circles) in molten FeOx at 1973 K (red) and 1673 K (blue). c Average O-Fe coordination number nOFe in molten FeOx at 1973 K (red circles) and 1673 K (black triangles) from our EPSR models. The estimated error bars are ±10 Fe3+/ΣFe (%) and ±0.25 on the nFeO or ±0.25(cFe/cO) for nOFe.

Full size image

EPSR modeling

In our EPSR models, which represent best fits of 3-dimensional atomic arrangements at the expected densities that reproduce our experimental data, nFeO is typically found to increase from ∼4.3(3) at 5.0% Fe3+ to 5.1(3) at 92% Fe3+. The values of nFeO as a function of Fe3+/ΣFe extracted from 1973 K to 1573 K are similar, indicating no temperature dependence within experimental error, except perhaps for the most reduced melts, where nFeO tends to increase slightly with temperature decreasing. The values of nFeO concur with previous findings on liquid iron oxide and silicate liquids and glasses obtained from XANES13,18, although linear extrapolations used for mixed valence FeOx to infer endmember nFeO appear to be inaccurate. Our results differ from early diffraction studies on compositions close to molten FeO, where Waseda et al.10 have reported an Fe-O peak position of 2.06 Å and nFeO = 5.5 at 1673 K using X-ray diffraction in a dry Argon atmosphere, and rFeO = 2.16 Å and nFeO = 5.7 at 1693 K using anomalous x-ray scattering11. A re-analysis of the Waseda et al.10 diffraction data through direct Fourier transformation of the digitized S(Q) and EPSR modeling reveals nFeO = 4.6(4), in line with our findings. In addition, we note that nFeO ≈ 5 for our most oxidized melt (92% Fe3+) approaches that found in the crystalline form of γ-Fe2O3, despite substantial differences in populations of 4-, 5- and 6-fold Fe-O polyhedra.

Discussion

Although we obtain Fe-O coordination numbers in reducing atmospheres that are substantially lower than the octahedral coordination found in Wustite, our measurements are in good agreement with first principles molecular dynamics (MD) calculations in related liquid systems (but none have yet been performed on FeOx). For example, it is found that nFe2+-O = 4.5 in FeMg3O4 (Holmström et al.20) and nFe2+-O = 4.9 in Fe2SiO4 (Ramo and Stixrude23). A similar average structure has been predicted for liquid MgO using ab initio MD19,24 (nMgO = 4.6) where Mg2+ has only a slightly smaller ionic radius than Fe2+. In the fayalite system, Alderman et al.25 found nFeO = 4.4 by x-ray diffraction, in reasonable agreement with Ramo and Stixrude23. Also values of 4.5 (Sun et al.26) and 4.7 (Guillot et al.27) have been predicted using classical MD simulations. Classical MD calculations of molten FeO by Seo et al.28 predict nFeO = 4.2, and are able to reproduce transport properties of the FeO-SiO2 system, although the model density is significantly lower than experiments would suggest29. To explore this further, we have carried out a series of classical MD simulations of molten FeO using seven published inter-atomic pair potential sets (see Table 2 and Supplementary Table S1 and Figs. S6–S8) and found the Guillot and Sator fayalite potential27 to be in the best (but not complete) agreement with the expected density and our experimental diffraction data, with nFe2+-O = 4.8, and four- and five-fold Fe-O polyhedral units dominating. However, no classical interatomic potential available for mixed valence iron oxides in the literature was found to reproduce the trend observed in Fig. 3.

Table. 2 A summary of Fe2+-O coordination numbers obtained from our MD simulations on FeO using various published pair potentials, compared to our EPSR model with the lowest value of Fe3+/ΣFe.
Full size table

A deeper understanding of the relative changes in nFeO can be obtained from an interrogation of the average coordination numbers of nFe2+O and nFe3+O. The coordination numbers of Fe2+-O and Fe3+-O in our EPSR models are similar to one another (and mostly invariant with temperature), similar to XANES findings in iron silicate glasses18, as both are found to increase with Fe3+/ΣFe, see Fig. 3b. Moreover, our EPSR models show that liquid FeOx primarily consists of 4-, 5- and 6-fold coordinated Fe2+ and Fe3+ atoms, and their relative populations vary with Fe3+/ΣFe (Fig. 4). No evidence was found for a bimodal FeO4 and FeO6 distribution in the liquid. Instead, a substantial amount of Fe3+ (38–49%) and Fe2+ (26–41%) are found in 5-fold coordination indicating that this is a dominant species in all FeOx melts. For ferrous iron, our EPSR models show Fe2+ predominantly occupies 4- and 5-coordinated sites, which is consistent with other divalent cations e.g. Ni2+, Fe2+, and Mg2+ in molten silicates30. For ferric iron, Fe3+, 4-fold and 5-fold polyhedra also dominate at low Fe3+/ΣFe, but for Fe3+/ΣFe >40% the number of 6-fold Fe3+ increases and the number of 4-fold Fe3+ decreases. This is opposite to the trend reported for binary silicate melts by Mysen15, which indicated nFe3+-O decreases from 6 to 4 with increasing Fe3+/ΣFe. In addition, it is a more complex scenario than that proposed by Hayashi et al.17 who only considered the existence of tetrahedral and octahedral Fe-O units in sodium-iron-silicate melts.

Fig. 4: The distribution of polyhedra in our EPSR models.
figure 4

a Fe2+Om b Fe3+On units (where m, n = 3, 4, 5, 6) and c OFep units (where p = 2, 3, 4, 5) from EPSR models at 1973 K (solid symbols) and 1673 K (open symbols), as a function of Fe3+/ΣFe (%) (Calculated using equation 12 of ref. 35). nFeO and nOFe were determined by the average of the three cut-off criteria. The clusters on the right correspond to representative OFe4, OFe3, and OFe5 units taken from the EPSR models.

Full size image

An alternative view of the change in liquid structure with redox can be obtained through the extraction of the oxygen-iron coordination number nOFe. In Fig. 3c, nOFe tends toward ~3.5 in pure O2 gas close to that found in γ-Fe2O3, but this increases only slightly to ~4.5 around FeO (compared to 6 in the crystalline phase). An analysis of the fraction of OFe coordination species with Fe3+/ΣFe in Fig. 4c in our EPSR models shows a dramatic effect. We find that the dominant OFe4 units account for approximately half the species in the melt and this remains relatively constant with Fe3+/ΣFe. For Fe3+ rich melts OFe3 units account for almost all of the remaining oxygen connectivity. As the availability of oxygen decreases and the percentage of Fe2+ increases, the OFe3 units are partially replaced by FeO5 units.

Finally, we consider the relationship between the presence of various local atomic structural units and the physical properties of liquid iron oxide slags. In binary ferrite melts it has been suggested that the atomic fraction of Fe3+O4 is the main factor in controlling the viscosity and the abundance of Fe2+ and Fe3+O6 units are main contributors to the electrical conductivity31. In addition, chemical reaction rates in the iron and steel making process are determined by concentrations and activities, where the activity is dependent upon bond strengths in different coordination environments17. However, a comparison of our diffraction data with a survey of classical MD simulations performed using various potentials in the literature has revealed the shortcomings of many interatomic potentials across the compositional range of FeOx melts, see Supplementary Fig. S7. Hence more sophisticated theoretical work and/or ab initio simulations are necessary to elucidate the precise ratios of 4-, 5- and 6-fold coordinated Fe2+-O and Fe3+-O in liquid FeOx, particularly in order to understand transport properties, and to address charge-transfer and small-polaron conductance phenomena. In fact, during the writing of this article we note that two ab initio molecular dynamics simulations on FeOx -bearing silicate melts have been published32,33. It is anticipated that our experimental pair distribution functions for molten iron oxides presented herein will provide ambient-pressure benchmarks to test models relevant to a wide range of industrial and natural processes.

Methods

Levitation and high-energy X-ray diffraction experiments

We formed powdered Fe2O3 (purity 99.995%) into roughly 3 mm diameter spheres using a laser hearth melter in the air. Weighed samples were floated in a 200–500 mL/min stream of gas, using a conical nozzle aerodynamic levitator in an enclosure25, and heated from above using a 400 W CO2 laser (Synrad Firestar i401). The droplet surface temperature was measured using 0.9 μm wavelength optical pyrometry (Chino IR-CAS) sighted onto the region of the sample where the X-ray measurements were made. A constant spectral emissivity of ελ = 0.92 was assumed for correcting the apparent temperatures and the accuracy of the pyrometer was ∼1% for temperatures in the range 1773–1973 K. The oxygen partial pressure, p(O2), experienced by the sample was controlled using different levitation gas mixtures, either 99CO:1CO2, 50CO:50CO2, CO2 (all CO:CO2 mixtures were 5% with a balance of 95%Ar), 99.9%Ar:0.1%O2 or pure O2, leading to a range of Fe3+/ΣFe at equilibrium. The p(O2) were calculated by standard methods34, based on the measured droplet surface temperatures and the Fe3+/ΣFe ratio was estimated using equation 12 from Jayasuriya et al.35, which has previously been found13 to be in accord with the molten iron oxide data of Darken and Gurry36. A detailed evaluation of the effect of errors in the calculated values of Fe3+/ΣFe is described in the uncertainty section and Supplementary Fig. S5.

High energy X-ray diffraction was performed at the Advanced Photon Source beamline 6-ID-D, (Argonne, IL, USA) using an incident X-ray beam of wavelength 0.123538 Å (100.36 keV). A Perkin Elmer XRD1621 area detector was placed ∼338 mm downstream of the sample and calibrated using a CeO2 powder standard (NIST) to provide a wide Q range of 0.5–25.0 Å−1. The x-ray beam was centered on the top 200 μm of the sample sphere to align with the laser beam and pyrometry measurement and avoid scattering from the levitator components. Before heating a sample, the chamber was purged with the selected process gas for 20 min. The levitated sample was then heated to a temperature just above its melting point, and held for 10–20 s to ensure complete melting. Diffraction patterns were collected for 40 s each in ∼100 K steps in a sequence of decreasing temperatures from 1973 K to at least 1673 K. Only scattered X-ray data from the top half of the Debye–Scherrer cone were used in the analysis to avoid significant attenuation corrections. The x-ray total structure factor S(Q) and the corresponding differential pair-distribution function D(r) = 4πρr[G(r) – 1] as a function of temperature and p(O2) were extracted. Details of the high-energy X-ray analysis used have previously been described by Shi et al.37 from samples prepared in a laser hearth melter38. Three-dimensional models describing the atomic structure were obtained using EPSR39. The Fe3+/ΣFe ratio was estimated using equation 12 from Jayasuriya et al.35, which was found to be in accord with the molten iron oxide data of Darken and Gurry36.

Atomistic modeling

Empirical potential structure refinement modeling

EPSR39 were used to obtain atomistic models of liquid FeOx based on the high-energy x-ray diffraction data. The EPSR simulations were performed on ∼6000 ions within a cubic box under periodic boundary conditions, and the temperature dependent densities fixed using the equations of Hara et al.29. A Lennard-Jones reference potential with Fe3+ and O2− parameters from Weigel et al.40 was used. For Fe2+, the reference potential parameters are given in Table 1, and were modified only slightly from those of Fe3+ by comparing initial Monte Carlo simulation predictions of the Fe-O bond length distributions to the diffraction data. The starting configuration was built according to the calculated Fe3+/ΣFe ratio35 and, following initial Monte Carlo equilibration, the empirical potential term was refined to improve agreement with scattering data, with the maximum allowed energy amplitude set to ∼30% of initial model energy. Once the goodness-of-fit parameter was minimized, structural data were collected over ensembles of at least 2000 configurations. The resulting EPSR fits were in excellent agreement with the high energy X-ray diffraction data shown in Fig. 2, see Fig. S1.

Classical molecular dynamics simulations

Classical molecular dynamics (MD) simulations with 0% Fe3+ were performed using the DL_POLY package41 on a system containing ∼6000 atoms. Several published potentials27,28,42,43,44,45 were tested. The results of Monte Carlo simulations under the EPSR reference potentials were used as starting configurations for the MD simulations, which were conducted either in the canonical (NVT) or the isothermal-isobaric (NPT) ensemble, using a Hoover-Nose thermostat or thermostat-barostat, at P = 1 atm and various temperatures. A time step of 1.0 fs was selected, with the first 25000 steps used for equilibration, and particle trajectories typically integrated over a further 50,000 steps (production runs) at each temperature point. The starting temperature was 2700 K and the systems were cooled stepwise in 100 K steps down to 300 K. In some cases (Fe2+-rich) crystallization was observed, and in all other cases glass formation occurred. The effect of Fe3+/ΣFe on the liquid structure was also investigated using the Teter and Guillot pair potentials27,42.