A DFT study of the structures, stabilities and redox behaviour of the major surfaces of magnetite Fe3O4

The renewed interest in magnetite (Fe 3 O 4 ) as a major phase in diﬀerent types of catalysts has led us to study the oxidation–reduction behaviour of its most prominent surfaces. We have employed computer modelling techniques based on the density functional theory to calculate the geometries and surface free energies of a number of surfaces at different compositions, including the stoichiometric plane, and those with a deficiency or excess of oxygen atoms. The most stable surfaces are the (001) and (111), leading to a cubic Fe 3 O 4 crystal morphology with truncated corners under equilibrium conditions. The scanning tunnelling microscopy images of the different terminations of the (001) and (111) stoichiometric surfaces were calculated and compared with previous reports. Under reducing conditions, the creation of oxygen vacancies in the surface leads to the formation of reduced Fe species in the surface in the vicinity of the vacant oxygen. The (001) surface is slightly more prone to reduction than the (111), due to the higher stabilisation upon relaxation of the atoms around the oxygen vacancy, but molecular oxygen adsorbs preferentially at the (111) surface. In both oxidized surfaces, the oxygen atoms are located on bridge positions between two surface iron atoms, from which they attract electron density. The oxidised state is thermodynamically favourable with respect to the stoichiometric surfaces under ambient conditions, although not under the conditions when bulk Fe 3 O 4 is thermodynamically stable with respect to Fe 2 O 3 . This finding is important in the interpretation of the catalytic properties of Fe 3 O 4 due to the presence of oxidised species under experimental conditions.


Introduction
Magnetite, Fe 3 O 4 , is of significant importance as the main component of industrial catalysts in, for example, the dehydrogenation of ethylbenzene 1 which is used as the primary feedstock for the production of 85% of commercial styrene. 2,3 Fe 3 O 4 is also used as a catalyst for the water gas shift (WGS) reaction, where molecular hydrogen is formed from carbon monoxide and water, 4-6 the Fischer-Tropsch synthesis of hydrocarbons 7 and the Haber-Bosch process for the production of ammonia. [8][9][10][11] The high stability and catalytic activity as well as its low cost make Fe 3 O 4 the catalyst of choice for these heterogeneous processes. 12 Furthermore, Fe 3 O 4 is important in other applications, such as groundwater remediation, 13 and potentially in spintronic devices due to the conducting properties of only one channel of spins. 14,15 Above 120 K, Fe 3 O 4 crystallizes in the spinel structure 16 with space group Fd% 3m (cubic), 17 but when cooled below that temperature, it undergoes a phase transition known as the Verwey transition, where the space group changes to P2/m (monoclinic). 17 Thus, at room temperature, Fe 3 O 4 has the spinel face-centred cubic unit cell, on which we will focus in this paper. In this structure, the oxygen ions are regularly close packed along the [111] direction, separating layers of Fe ions, which appear in two different alternate arrangements. One is composed of Fe ions occupying two types of positions (octahedral (Fe B ) and tetrahedral (Fe A )) and the other one has only Fe B , shown in the scheme in Fig. 1 18,19 extra-terrestrial 20,21 and synthetic 22,23 Fe 3 O 4 crystals have been described by several authors. Among all the crystal habits in which this mineral has been found, the three most common are the octahedral morphology enclosed by (111) surfaces; a truncated octahedron by adding the (001) plane and as twinning on the (111) surface. 12,22 Zhao et al. synthesized Fe 3 O 4 under a systematic range of conditions using a polyol process, where the crystals obtained ranged from cubic, truncated octahedral to octahedral shapes, depending on pH. 22 The stacking sequence of the atomic planes perpendicular to the [001] direction can be represented as Fe A -(O-Fe B ) (atoms inside brackets are within the same layer), leading in principle, to two different bulk-like terminations for the (001) surface, which are all polar. There are also two possible non-dipolar reconstructions of this surface, i.e. when the slab is terminated by either 0.5 mono layers (ML) Fe A or 0.5 ML O-Fe B in both the top and bottom surface. Experimentally, this surface has been found to have a ffiffi ffi 2 p Â ffiffi ffi 2 p À Á R45 reconstruction for which different rationalizations have been given. Studies combining low-energy electron diffraction (LEED), X-ray photoelectron spectroscopy (XPS), X-ray photoelectron diffraction (XPD) and scanning tunnelling microscopy (STM), 24 as well as another work combining LEED with low-energy ion scattering (LEIS), 25 have suggested a (001) surface terminated by the reconstructed non-dipolar 0.5 ML of Fe A . On the other hand, a different study, combining STM, LEED, LEIS and XPS, has suggested a surface terminated by the reconstructed charge-compensated O-Fe B with one oxygen vacancy per unit cell. 26 Meanwhile, Voogt et al. 27 were unable to differentiate them based on reflection high-energy electron diffraction (RHEED) and LEED, suggesting as possible terminations: the reconstructed non-dipolar 0.5 ML Fe A layer or the reconstructed charge-compensated O-Fe B layer with oxygen vacancies or hydroxyl groups. 27 More recently, Parkinson et al. 28 have identified experimentally using STM and LEED images that the (001) surface terminations are temperature dependent. The 0.5 ML O-Fe B termination or one with wavelike structure and small defects, such as hydroxyl groups, is thermodynamically more stable at 923 K, while at lower temperatures (623 K) the surface terminated by 0.5 ML of Fe A is observed, although some point-defects may stabilise other terminations. 28 To date, most computational efforts have been concentrated on explaining the stability of the bulk-like dipolar O-Fe B termination, leaving largely ignored the reconstructed non-dipolar 0.5 ML Fe A termination. Pentcheva et al. 29 have studied the stability under varying redox conditions of one ideal and reconstructed stoichiometric (0.5 ML Fe A ) and several non-stoichiometric terminations using spin-polarised density functional theory (DFT) calculations. These authors found that the modified polar bulk-like O-Fe B termination was the most stable for the whole range of chemical potential, 29 which was validated by experimental X-rays diffraction (XRD) 29 and by the wavelike pattern along the [011] direction observed on experimental STM images. 30 Further studies by Parkinson et al. 31 of this surface termination using spinpolarized DFT + U calculations supported the Jahn-Teller distortion of this surface based on simulation of STM images of antiphase domain boundaries (APDBs). 31 In the [011] direction, Fe 3 O 4 is composed of alternating layers of (Fe A -Fe B -O) and (Fe B -O). After reconstruction, two non-dipolar terminations are possible: one terminated by the (Fe B -O) layer with 0.25 ML Fe A on top and another terminated by the (Fe A -Fe B -O) layer with 0.25 ML of Fe A vacancies. Single crystal studies carried out on this surface, involving the use of STM, LEED, scanning tunnelling spectroscopy (STS) and Auger electron spectroscopy (AES), have found a one-dimensional reconstruction along the [01% 1] direction, which was concluded not to have a simple bulk iron-oxide termination. 32,33 Subsequent studies supported by atomically resolved STM 34-36 suggested a model based on a surface terminated by a polar (Fe A -Fe B -O) bulk-like layer, in order to explain the atomic rows observed on the tops of ridges along the [01% 1] direction. However, the authors also left open the possibility of alternative models, including surface reconstruction, to interpret the STM images. 34 The (111) surface is the dominant cleavage plane of Fe 3 O 4 , and the stacking of the atomic layers perpendicular to this surface is Fe A1 -Fe B1 -Fe A2 -O 1 -Fe B2 -O 2 . All of the six possible different bulk-like surface terminations are dipolar. Only two reconstructions lead to non-dipolar terminations, i.e. 0.5 ML Fe B1 or 0.5 ML Fe B2 . Several possible terminations have been described from LEED and STM results: one dipolar plane showing close packed features (due to 0.75 ML of Fe B2 and 0.25 ML of O 2 over a close packed O 1 layer); 37 a reconstructed non-dipolar honeycomb plane (due to 0.5 ML of Fe B1 ), which was the most stable one; 37 a reconstructed dipolar 0.25 ML Fe A1 plane; 38 as well as a regular bulk-like dipolar Fe A1 termination and an intermediate case between the former two, which were obtained as a function of the annealing temperature. 39 Again, most of the computational studies have been directed towards the dipolar bulk-like terminations of the (111) surface. Martin et al. 40 used spin-polarized DFT calculations to study the dipolar non-stoichiometric bulk-like Fe B1 and Fe A1 terminations of the (111) surface and they found Fe B1 to be the more stable of the two, which they validated through comparison with experimental STM images. 40 Berdunov et al. 41 used DFT to study the dipolar non-stoichiometric bulk-like O 2 termination of the (111) surface which was also validated via comparison with experimental STM images. 41 Kiejna et al. 42 studied the non-stoichiometric bulk-like dipolar terminations of the (111) surface using DFT + U, and although they did not calculate the stoichiometric slab, they predicted the Fe A1 termination as the most stable for the whole range of chemical potential they considered. 42 Reduced surfaces of the Fe 3 O 4 (111) surface show superstructures with Fe 1Àx O(111) islands, 43 which makes the surface even richer in possible terminations.
Following the seminal work by Tasker 44 on the surface properties of ionic solids, in the present work we have used DFT + U to investigate the non-dipolar stoichiometric terminations of the low Miller index surfaces of Fe 3 O 4 , in order to complement previous experimental and computational studies. We report the equilibrium morphology of the crystals enclosed by stoichiometric non-dipolar surfaces and the factors that govern the redox properties of the most common surfaces, (001) and (111), which are also the most prominent surfaces of Fe 3 O 4 moieties. 45,46 We have also calculated the STM images of the different stoichiometric non-dipolar terminations of these surfaces to determine the most likely to appear in nanocrystals through comparisons with available experimental STM data. We have studied the redox processes by the systematic formation of single O vacancies and the adsorption of single O on the surface, as opposed to previous computational studies which have focused on bulk-like terminations and their modifications. This approach allows us to explore how these redox processes determine the surface properties by finely tuning the conditions of temperature and oxygen partial pressure on the stoichiometric non-dipolar surfaces.

Calculation details
We have used the Vienna Ab initio Simulation Package (VASP) [47][48][49][50] to carry out quantum mechanical calculations within the usual Kohn-Sham (KS) implementation of DFT. The Perdew-Burke-Ernzerhof (PBE) 51,52 version of the generalized gradient approximation (GGA) was employed as the exchange-correlation potential, together with the semiempirical method of Grimme 53 to model the long-range dispersion interactions. The core electrons, up to and including the 3p levels of Fe and the 1s of O, were frozen and their interaction with the valence electrons was described by the projector augmented wave (PAW) method. 54,55 KS valence states were expanded in a plane-wave basis set with a cutoff of 400 eV for the kinetic energy. An energy threshold-defining selfconsistency of the electron density was set to 10 À5 eV and the optimization of the structures was conducted via a conjugate gradients technique, which stops when the Hellmann-Feynman forces on all atoms are less than 0.01 eV Å À1 .
All calculations were spin-polarised, but spin-orbit coupling was not considered. Within the VASP code, it is possible to assign an initial spin population and orientation to each atom of the system, to converge to a particular spin configuration. Thus, the initial magnetic moments were set following a highspin ferrimagnetic structure, i.e. with opposite spins in the tetrahedral and octahedral sites, in agreement with experiment. 56,57 In order to describe the electronic and magnetic behaviour properly, an accurate treatment of the electron correlation in the localized d-Fe orbitals is crucial. Hence, we have used the Dudarev et al. 58 approach within the DFT + U 59 for improving the description of these localized states. This is a correction typically used where standard LDA and GGA functionals fail to describe the electronic structure properly. 60 The value for the on-site Coulomb interaction term in this study was U eff = 3.7 eV, which was obtained following the procedure described in ref. 15 but adjusted to a different DFT functional. The limitation of this approximation is the difficulty in choosing an adequate value for the U eff parameter, which is usually property dependent. [61][62][63] An alternative computational approach is the use of hybrid functionals, although in that case the calculated properties also depend on the fraction of the exact Hartree-Fock exchange that is added to the DFT exchange term, 60,64-67 and the calculations are too computationally expensive to be applied to the large surface models employed in this study.
Bulk calculations were performed using the rhombohedral primitive unit cell containing 14 atoms (Fe 6 O 8 ). Integrations in the reciprocal space were performed using a Monkhorst-Pack grid of 7 Â 7 Â 7 G-centred k-points, 68-70 which ensured electronic and ionic convergence. Test calculations with a higher number of k-points led to an energy difference smaller than 1 meV per cell. k-Point grids for the surface calculations were chosen in a such a way that a similar spacing of points in the reciprocal space was maintained.
Within this setup, we calculated a lattice constant for the bulk  15,57 following very closely the Néel model, 56 where the electronic configurations are e 2 " t 3 2" for Fe A and t 3 2g" t 1 2g# e 2 g" as well as t 3 2g" e 2 1.521 and À1.158 e À respectively. For further information, the density of states of the Fe 3 O 4 bulk phase are provided in Fig. S1 of the ESI. †

Surface models
In order to build slab models of the Fe 3 O 4 surfaces, two models are used in the literature to explain the reconstructions found in polar surface terminations: the electron-counting 72 and the dipole method. 44 Both models are based on the condition that the net surface charge or dipole perpendicular to the surface, respectively, must be zero. A surface structure satisfies the model of electron-counting (i.e., it is charge-or auto-compensated) if all the partially filled dangling bonds in the cations are empty and the partially filled dangling bonds in the anions are full. It assumes that the atomic orbitals are in the conduction or valence band respectively. To achieve this, the model postulates that a stable surface structure will be the one that (after reconstruction) is able to accommodate exactly all the electrons of the partially filled orbitals of the cations (in the conduction band) into the partially filled orbitals of the anions (valence band). However, the disadvantage of this approach is that this condition directly links to the conducting properties of the material under investigation. If the surface satisfies this model, it will be a semi-conductor; otherwise the remnant electrons will lead to a metallic surface.
The dipole method for reconstructing dipolar surfaces is a more robust option, at least with half-metallic materials like Fe 3 O 4 , 15 because it is not connected to the conducting properties of the structure. This method, pioneered by Tasker, 44 considers the crystal as a stack of planes, where three possibilities can arise. 44 In type 1, each plane has overall zero charge because it is composed of anions and cations in stoichiometric ratio which makes it non-dipolar, whereas in type 2 the stacking of three symmetrically charged planes cancels out the dipole moment perpendicular to the surface. In type 1 and 2, no reconstruction of the surface is needed because the repeat unit is non-dipolar perpendicular to the surface. However, in type 3 surfaces, alternating charged planes stack in a sequence that produces a dipole moment perpendicular to the surface, but the surface can be reconstructed through moving half of the ions with the same charge from the top to the bottom of the slab. This method also guarantees that the surface does not generate an electrical field within the crystal and therefore the potential felt at each ion site reaches the constant bulk value, a condition that is not necessarily satisfied by the electroncounting model.
All the surfaces in this study were created by cutting the geometry optimised bulk using the dipole method implemented in METADISE. 73 The resulting slabs are represented by keeping fixed the bottom atoms at their ab initio relaxed bulk positions to simulate the bulk phase of Fe 3 O 4 and by relaxing the rest of the slab during the optimization, giving a single relaxed surface. The slabs comprise 56 atoms (8 formula units), with 24 Fe and 32 O atoms. The Fe 3 O 4 (001), (011) and (111) surfaces were modelled with slabs having surface areas of 70.5, 99.7 and 61.1 Å 2 , respectively, and they were constructed of 9, 5 and 13 atomic layers, respectively. Fig. 2(c) provides a representation of their stacking sequence in each direction. For the (001) and (111) surfaces, the simulation slabs were symmetrical along the z axis. Top species in the (001) surface were (0.5 ML) Fe A atom and 2 Fe B and 4 O atoms (equivalent to 0.5 ML for each of the ions) for terminations A and B respectively, see Fig. 2. For the (111) surface, terminations A and B were terminated by half of the (Fe B ) 6 and (Fe B ) 2 bulk layers respectively. However, the simulation slabs for the (011) surface were asymmetrical along the z axis, with complementary top and bottom layers. Top layer of termination A was a single (0.25 ML) Fe A atom above a bulk-like O-Fe B layer, while its bottom layer was a bulk-like Fe A -Fe B -O layer with one (0.25 ML) Fe A vacancy. For termination B, top and bottom layers were the other way round.
In every simulation cell, a vacuum region of 12 Å was added perpendicular to the surface to avoid interaction between the periodic slabs. Different slab and vacuum thicknesses as well as numbers of relaxed layers were tested until convergence within 1 meV per cell was achieved. Since we are going to remove and add O atoms to the surfaces at one side of the slab only, we applied dipole corrections perpendicular to all surfaces in the calculations to enhance the electronic convergence. 74,75 We have used Bader analysis 76 in the implementation of Henkelman and co-workers [77][78][79] to analyse the charge transfer around the defects introduced in the stoichiometric surfaces. We have chosen this methodology for partitioning atomic charges, as it is based upon the charge density, which is insensitive to the metal oxidation state and the basis set used, unlike wavefunctionbased population schemes. [80][81][82] The DOS of the optimised surfaces is presented in Fig. S1 of the ESI. †

Calculation of surface energies
We have carried out energy minimisations of the (001), (011) and (111) slabs to obtain their surface energies. We derived the surface energy of the unrelaxed surface (g u ) from a single point calculation of the pristine symmetric stoichiometric slab before relaxation, via equation: where E slab,u is the total energy of the unrelaxed slab, E bulk is the energy of the bulk containing the same number of formula units as in the slab and A is the surface area of one side of the slab. During relaxation, the top surface was allowed to relax and the bottom one was kept fixed. As this slab model does not provide an isolated relaxed surface and both sides of the symmetric stoichiometric slabs are considered in the calculation of the energy, their surface energies (g r ) and (g u ), for the relaxed and unrelaxed sides respectively, are related by eqn (2), where E slab,r is the slab total energy after relaxation.
At this point it is also worth noting that eqn (1) and (2)   terminations A and B of (011), as the slabs are asymmetric and complementary. The cleavage energy (E clev = 2Ág r ) is thus related to the energy required to create both top and bottom surfaces of the slab. We have also calculated the degree of relaxation of each surface as a percentage (for (011) g r B E clev /2): The equilibrium morphology of a Fe 3 O 4 crystal is determined by the surface free energies and the related growth rates of the various surfaces, which provides a measure of the relative stabilities of the surfaces. The morphology is constructed according to Wulff's theorem, 83 where the distance from the centre of the particle to the surface is proportional to the surface energy. It is based on the Gibbs approach, 84 who proposed that under thermodynamic control the equilibrium form of a crystal should possess minimal total surface free energy for a given volume. Previous studies have shown 85,86 that using surface energies to calculate crystal morphologies provides good agreement with experiment as the difference in entropy between bulk and surface is small.

Redox processes of the (001) and (111) surfaces
We have also examined the redox properties of the most common Fe 3 O 4 surfaces, the (001) and (111), by removing or adding O atoms to form non-stoichiometric compositions of the top atomic layer. We have based the discussion of the stabilities of the nonstoichiometric terminations on the ab initio thermodynamics formalism 87 where the surface free energy (s) is calculated according to the equation: where Ds(T, p) is the difference between the surface energy of the stoichiometric surface and the surface free energy of the non-stoichiometric plane and G O is the excess of O ions at the top surface of the slab expressed in eqn (5) It is possible to express the chemical potential of molecular O 2 (m O ) in the gas phase as: Here the first term within the bracket is the DFT energy of the O 2 molecule. The second term is the difference in the Gibbs free energy per O 2 molecule at p 0 = 1 bar between 0 K and T, which in this study has been extracted from thermodynamic tables 88 to avoid its calculation in the gas phase. 89,90 The last term represents the change in free energy of the O 2 gas (assuming ideal gas behaviour) at constant temperature (T) when its partial pressure changes from p 0 to p.
We express m O with respect to half the energy of the O 2 molecule. The above convention makes m O a function of only experimental quantities. For consistency in the evaluation of the slab energies, we must subtract half of the energy of the O 2 molecule for each O atom in the slab. Expressing m O as described, it is possible to plot the surface free energies given by eqn (4) for different surface compositions as a function of m O , and discuss the redox behaviour of the surface. 89,90 Finally, for the calculation of the energy required to create an O atom vacancy or to add the atom on the surfaces, we need the energy of the O 2 molecule. However, it is known that GGA calculations fail in the description of the binding energy for this particular molecule, as is shown in the (over)binding of the O 2 molecule. 52 According to our calculations, the O 2 triplet ground state has an equilibrium bond length of 1.23 Å and a binding energy of À6.08 eV (with respect to triplet oxygen atoms), comparing well with previous computational studies. 62,91,92 However, this value lies 0.91 eV below the experimental binding energy (À5.17 eV). 93 Therefore, we have considered that half of the over-binding of the O 2 molecule, 0.46 eV, will be added to correct the adsorption or vacancy formation energies with respect to one O atom. The redox processes in the following sections are all reported with respect to the corrected value. A comparative analysis with uncorrected energy values is provided in the ESI. †

Calculation of scanning tunnelling microscopy (STM) images
The STM images were simulated according to the basic formulation of the Tersoff-Hamann approach 94 where the STM tip was approximated to an infinitely small point source. The tunnelling current between the surface and the tip in the STM experiments is proportional to the local density of the states (LDOS) integrated between the Fermi energy and the sample bias. We have used the program HIVE 95 for the production of our STM topographic images, where the DFT-based partial charge density was integrated from À2.5 eV to the Fermi energy. In the constant current mode, the tip of the STM is moved across the surface where its height varies to keep the charge density at a constant value, which is given by a constant LDOS. We map the simulated STM images by means of the heights as a function of the position of the tip over the surface. More details about the method can be found elsewhere. 96 3 Results and discussion

Stoichiometric surfaces
We have modelled different terminations of the three lowest Miller index surfaces of Fe 3 O 4 , shown in Fig. 2, whereas Table 1 summarises their surface energies before and after energy minimization. Before relaxation, the order of increasing surface energies, and therefore decreasing stability, is (001) o (111) o (011), which remains the same after relaxation. Note that this order was established by taking into account only the most stable termination (with lowest g) per surface, as these terminations would be the most likely to appear for each plane.  Table 2. The relaxation pattern of the top atomic layer of the surface slab agrees semiquantitatively with the B40% inward relaxation reported for the topmost layer of this termination based on LEIS analysis, 25 which is generally regarded to fit better than the more complex relaxation pattern reported before for this surface termination by Chambers et al. 24 Previous studies, purely theoretical 98 or combined with experiments, 99 98 or from a dipolar bulk-like Fe A terminated (001) surface. 99 However, we have not included dimers here as this lies outside the scope of the present study. The surface energy of termination B of the (001) surface is also reported in Table 1, but we do not consider this plane for further analysis because of its high surface energy, which makes it very unlikely to appear in the Fe 3 O 4 crystal morphology.
The stacking sequence of the Fe 3 O 4 (011) surface is shown in Fig. 2 and the vertical shifts of the ions towards the vacuum after energy minimization are listed in Table 2. One of the two lowest energy surface terminations, termination A, terminates with 0.25 ML of mono-coordinated Fe A at the surface followed by a bulk-like structure consisting of single rows of 4-coordinated Fe B ions shifted 25% in the [01% 1] direction and alternating with single rows of O ions with cubic packing. During energy minimisation, the protruding Fe A ions move 0.98 Å towards the bulk, thereby compressing the surface layer atoms underneath which move horizontally to accommodate this relaxation. Termination B has essentially the same relaxed surface energy as termination A but it differs in its structure. It is terminated with a bulk-like structure consisting of a single row of 4- 34 Finally, the bottom two panels of Fig. 2 represent the stacking sequence of the Fe 3 O 4 (111) surface terminations, while the vertical displacement of the ions in the surface regions during the optimisation are listed in Table 2. One of the two lowest energy terminations, termination A, contains 0.5 ML of 3-coordinated Fe B ions with c(2 Â 4) symmetry, occupying hexagonal close packed (hcp) hollow positions in the top layer. The next layer has a bulk-like structure consisting of rows of O ions along the [0% 11] direction with rhombohedral packing. The percentage relaxation experienced by this surface termination is the largest   22 We have derived a Wulff 83 crystal morphology of pristine Fe 3 O 4 using the lowest surface energies for each Miller index. Its calculated equilibrium morphology is then expressed as a cubic shape with truncated corners, Fig. 3(a). As expected, the (001) plane dominates the morphology, followed by the (111) surface truncating the corners of the cube. The (011) surface does not appear in the morphology of Fe 3 O 4 due to the mathematical relation between the energy of the surfaces and their position in the crystal. 92 Despite the (011) surface having a surface energy of the same order of magnitude as the others, it is not expressed in the Wulff construction due to competition with the (001) surface. The ratio between their surface energies: gð011Þ=gð001Þ ¼ 1: 43 4 ffiffi ffi 2 p , see Fig. 3(b), and, as shown in Fig. 3(d), the (011) surface would only become present in the crystal morphology if gð011Þ=gð001Þ o ffiffi ffi 2 p . There are many ways to modify the shape of nanoparticles, such as solvent, media, capping agents, temperature or viscosity, but the Wulff morphology shown in Fig. 3(a) expresses a particle produced in conditions of perfect thermodynamic equilibrium, vacuum and at 0 K. Nevertheless, our results compare well with the morphologies of crystals synthesised by Zhao et al., 22 who described the formation of Fe 3 O 4 under different pH conditions. They increased the OH À concentration and the resulting crystal shapes changed from cubic (or spherical -depending on other conditions) at low pH via truncated octahedral to octahedral at high pH values. All their crystals showed mainly the (001) and (111) surfaces but, in some cases, a little (011) surface was expressed due to certain conditions which may modify the surfaces' relative energies. The occasional appearance of the (011) surface is rationalised in terms of kinetically-controlled anisotropic growth of the Fe 3 O 4 nanoparticles. Zhao et al., 22 suggested that a high concentration of KOH in the solution can lead to selective adsorption of the hydroxyl anions to certain planes of the crystal, which slows down considerably their growth process. Therefore, the presence of these ions can affect the relative stabilities of the different crystal surfaces. The inversion of the nature of the inequality gð011Þ=gð001Þ 4 ffiffi ffi 2 p , which already lies close to ffiffi ffi 2 p , will cause the (011) plane to show up in the morphology.
3.1.2 Scanning tunnelling microscopy images simulation. From the optimised structures of the planes and terminations that are expressed in the morphology, i.e. termination A of (001) and terminations A and B of (111) surfaces, we have derived the topographical STM images. These images provide information about the spatial distribution of the valence band states in the vicinity of the Fermi energy (E F ). The above is particularly useful for systems where atoms (in our case O atoms) can be added or removed from many different positions on the surface. Modelled STM images may help to clarify experimental ones by direct comparison, for instance to identify between the two possible terminations of (111) surface, whose surface energies are very close, and also to validate the most stable termination of the (001) surface. The model also avoids any external perturbations, like the electric field of an STM experimental tip, which can influence the position of atomic species adsorbed on a surface. 96 The STM images in Fig. 4 are calculated on pristine Fe 3 O 4 (001) and (111) surfaces. Fig. 4(a) shows the STM image of the  The STM image of the Fe 3 O 4 (111) surface termination A is shown in Fig. 4(b) acquired at a density of 0.0055 e Å À3 and a distance of 1.50 Å to the tip. The image resolves the protruding Fe B as the brightest spots along rows in the [1% 10] direction and as dots in between these rows. The undulation of the rows is due to the 0.25 ML of O atoms in the 2nd layer that have moved towards the surface after minimization. From the modelled STM, we can even observe the rhombohedral packing of the sub-layer O ions.
The last STM image in Fig. 4(c) corresponds to termination B of the Fe 3 O 4 (111) surface obtained with a density of 0.0276 e Å À3 and the tip at 0.70 Å from the highest atom. The image acquired resolves the protruding Fe B as the brightest spots in the STM image with a p(1 Â 2) symmetry and the Fe A ions from the layer below which are always bonded to three O atoms immediately underneath. This atomic arrangement forms a pattern of incomplete hexagons (with Fe atom vacancies in one vertex of the imaginary hexagon) which can be seen as a thermally equilibrated structure with vacancies evenly distributed. Details of the layers further below are also visible in our STM image. Experimental studies of the Fe 3 O 4 (111) surface 37 have shown that among the two different terminations considered there, the one with 0.50 ML of Fe atoms is more stable than the one with 0.75 ML of Fe atoms and 0.25 ML of O atoms, agreeing well with our model of termination B of the (111) surface, whose simulated STM is shown in Fig. 4(c). The calculated vertical distance between the Fe A in the vertex of the hexagon and the O ion in its centre is 0.50 Å, which also agrees well with the value reported experimentally, 0.5 Å. 37 This experimental termination shows regions with full hexagons and others with incomplete hexagons (due to Fe vacancies). This atomic rearrangement may be a consequence of the high temperatures to which the surface was exposed.

Redox behaviour
We have studied the redox properties of the most stable terminations, A and B, of the Fe 3 O 4 (001) and (111) surfaces, respectively, by comparing the surface free energies corresponding to different O to Fe ratios at the surface. We maintained the number of Fe atoms in the slab as in the stoichiometric surface, but we modified the number of O atoms in the top layer by G (given by eqn (5)), as we were interested in studying the effect of different temperature and oxygen pressure on the stoichiometric non-dipolar surfaces. Because of the size of our supercells, and assuming that O atoms occupy bulk-like positions around the surface Fe atoms, 17 values of G are possible if we constrain the calculations to a maximum of one ML of adatoms or vacancies. However, due to the complexity of the (001) and (111) surfaces and in order to reduce the number of G to small values that reflect realistic m O , we have used only five values of G: G = 0 is the stoichiometric surface, G = +1, +2 are the partially oxidized surfaces, and G = À1, À2 are the partially reduced surfaces. We have represented all of them schematically in Fig. 5. 3.2.1 Reduction of the (001) surface. We discuss the first reduction process G = À1 by removing one O atom at the top surface of the slab, which leads to a 0.125 ML of O vacancies, with a vacancy formation energy (E vac ) calculated as where i takes values 0 and À1 in the first and second reduction respectively. At the surface there are three different types of O depending on the distance to the protruding 2-coordinated Fe A ion, see Fig. 5(a). Thus, the energy required to remove the first O from the surface is 2.60 eV for the atom furthest removed from this Fe A (see Fig. 5(a) for G = À1) and 3.28 eV for the one at intermediate distance. The vacancy created at the third type of O position has an even bigger energy and is therefore very unlikely. The comparison of these energies with the vacancy formation energy in the Fe 3 O 4 bulk (2.12 eV), suggests that under thermodynamic equilibrium any surface vacancies will migrate towards the bulk, a phenomenon that has also been observed to occur in another transition metal oxide VO 2 . 92 The tendency of the vacancy to migrate towards the bulk might seem contradictory with the fact that surface oxygen has a lower coordination number than bulk oxygen. However, this can be rationalized in terms of the oxygen vacancy-containing bulk material undergoing a different degree of relaxation than the oxygen vacancy in the surface slab, thereby driving the creation of the oxygen vacancy in the bulk. We proceed with the second reduction of Fe 3 O 4 (001) leading to G = À2. We removed an O located in the pristine row along the [110] direction, see Fig. 5(a) for G = À2, which is at intermediate distance to Fe A . This second vacancy is 3.23 eV less favourable than the previous state but it is just more likely than removing a more distant O ion from the row where the vacancy is now being created, 3.31 eV. This indicates that although the first vacancy is created preferentially in a position far away from the 2-coordinated Fe A , the second reduction might lead to a vacancy in the following O row at an intermediate distance from Fe A . As the energies to create the second vacancy in the two positions already described are within the DFT error, it might also be possible to find vacancies in the O positions further away from the row where the vacancy is now being created. We have characterised the G = À1 surface by means of a Bader analysis and compared the atomic charges with those on the pristine surface. The positive charge of the protruding Fe A ion was increased by a negligible amount (o0.05 e À ), where this small variation can be accounted for by the defect that was created at the farthest O location. The surface Fe B ions, however, are reduced, especially the ones closest to the vacancy with a variation in charge of 0.25-0.37 e À . This can be interpreted in terms of the number of O ions directly coordinated to the Fe B ions, see Fig. 5(a) for G = À1, where just over 80% of the electron density is transferred to the Fe B ions after removing the O atoms.
3.2.2 Reduction of the (111) surface. We have also explored different positions for the creation of the O (G = À1) vacancies in the Fe 3 O 4 (111) surface to find the lowest-energy configuration for this particular surface. We found that the process is thermodynamically even more unfavourable than on the (001) surface by 0.24 eV. The most likely vacancy is created in the centre of an incomplete Fe-hexagon, see Fig. 5 To remove an O atom coordinated to the 3-coordinated Fe A (opposite to the missing Fe in the incomplete hexagon) requires an energy of 3.56 eV. As in the (001) surface, any vacancy created in the (111) surface will be thermodynamically prone to migrate towards the bulk.
Creating a second vacancy among the atoms coordinated to the 3-coordinated Fe A , Fig. 5(b) for G = À2 costs 3.45 eV, which is less costly by 0.19 eV than removing the left O within the hexagon. These energies provide information about the consecutive reduction mechanism, where the first O vacancy is created in the centre of the incomplete Fe-hexagons and the next in one of the atom positions coordinated to the 3-coordinated Fe A .
The Bader analysis indicates that upon vacancy formation on the G = À1 surface, charge transfer on the Fe A is negligible and the protruding 3-coordinated Fe B is only slightly reduced. However, the Fe B ions whose charge is affected more are those in the 4th atomic layer (see Fig. 2) below the removed O atom to which they were previously directly coordinated. Altogether, the charge on those three Fe B is reduced by B0.89 e À , i.e. they have accepted 78.5% of the electron density previously held by the removed O.
3.2.3 Oxidation of the (001) surface. Another process we have studied is the surface oxidation by adsorption of one O atom leading to 0.125 ML of adatoms (G = +1). This process involves an adsorption energy per adatom derived from the equation, (where j takes values 0 and +1 in the first and second oxidation respectively). For j = 0, E ads is calculated at À1.87 eV. We considered that the O adatom is located similarly to the bulk structure, interacting simultaneously with the protruding 2-coordinated Fe A and one of the 5-coordinated surface Fe B , see Fig. 5 15 as is expected due to contraction of the top atomic surface layers after relaxation. The strongly exothermic adsorption suggests a favourable oxidation process, but as it is affected by m O , its evaluation requires a complete analysis of the gas partial pressure in equilibrium with the surface. Adding a second O atom (G = +2) is also an exothermic process, releasing 0.96 eV per adatom. The second O atom preferentially coordinates the protruding Fe A and a 5-coordinated Fe B , forming another O bridging structure collinear with the [% 110] direction, Fig. 5(a) for G = +2. As for G = +1, the top atomic contraction leads to short Fe-O distances, 1.85 Å. Another conformation for the second O adsorption is coordinating equivalent atoms but forming a V-shaped structure, leading to a weaker adsorption (E ads = À0.80 eV).
At this point, it is worth mentioning that although we started from the ideal terminations similar to the bulk when we added the first and second oxygen atom, this did not prevent them to relax to a different position. In fact, we can see in Fig. 5(a) for G = +1 (and +2), that after surface relaxation, the added oxygen has moved from its bulk site to another position, closer to the protruding Fe A . This finding agrees with the work of Reuter and Scheffler, 90 who found for RuO 2 (110) that terminations at positions different from the bulk can be important in non-stoichiometric compositions.
The Bader analysis on the density of the (G = +1) oxidised (001) surface shows the oxidation of the top layer Fe B by 0.60 e À while the protruding 2-coordinated Fe A ion only donates 0.04 e À to the newly added O atom. Hence, the O adatom gains 1.00 e À mainly from the surface metals whereas the charge of the surface anions (all of them more negative than the adatom) change by about 0.02 e À .
3.2.4 Oxidation of the (111) surface. On the (111) surface, the adsorption of one O atom (G = +1) led to the formation of a bridging structure between the protruding Fe B and one of the three closest 3-coordinated Fe A (see Fig. 5(b) for G = +1) releasing 3.00 eV. A less stable configuration is the one where the O-adatom is sitting atop the protruding Fe B providing an E ads of À2.04 eV.
The addition of a second O-adatom coordinating the protruding Fe B and one of the other two closest unoccupied 3-coordinated Fe A releases 2.30 eV; see the schematic representation in Fig. 5 Unlike the (001) surface, where the added oxygen atoms moved during the energy minimization, in the (111) we observed instead that the protruding Fe B ion moved from its bulk position after the addition of two defects (G = +2), see Fig. 5(b). This validates in our methodology the possibility of exploring nonbulk-like relaxed positions for any atom in the surface of our slab, as long as all non-equivalent bulk-like positions for the defects are carefully investigated.
The Bader analysis indicates that in the preferred structure for G = +1, the adatom gains 1.04 e À , where the charge of the other surface O atoms decreased as little as in the (001) surface. Amongst the two Fe ions coordinated to the added O atom, Fe B increases its charge by 0.08 e À , but Fe A by 0.39 e À . The charge on other surface Fe A and Fe B ions upon addition of the O atom changed by an average of 0.03 e À and À0.01 e À respectively.

Temperature and pressure effects
In this section, we discuss the thermodynamics of the redox processes at the (001) and (111) surfaces as a function of temperature and O 2 partial pressure in the gas phase. We express these macroscopic parameters by m O .
In Fig. 6(a), we have plotted m O in terms of temperature and the log p, along abscissas for easy comparison with the plots in Fig. 6(b) and (c). All the information used for the construction of Fig. 6(a) comes from experiment 88 and is independent from the calculations (see Computational methods section). Variations in T and p are necessary to modify the value of m O as required and to reflect different reducing or oxidising conditions. For example, the oxygen chemical potential is À0.3 eV (which is a typical oxidising value) at ambient conditions, i.e. at the intercept of T = 298.15 K and p = 0.21 bar, while more reducing conditions (lower values of m O ) can be achieved by increasing T while keeping the pressure constant (i.e. horizontal solid line in Fig. 6(a)).
The area between the two vertical dashed lines (m O from À3.13 to À2.44 eV) in Fig. 6 Fig. 6 as a function of G for inferring the intercepts of further oxidation/ reduction of these surfaces. We have limited this treatment to three defects, as a higher number can lead to the formation of At m O = À0.3 eV (ambient conditions), the (G = +3) oxidation of the Fe 3 O 4 (001) surface will take place, see Fig. 6(b). For the conditions where bulk Fe 3 O 4 is the most stable oxide, the (G = +3) oxidized (001) remains the most stable surface up to m O = À1.25 eV, from where the surface experiences a progressive reduction. In the early stages of this reduction, the unit cell loses two O atoms and remains so until m O = À1.85 eV. Beyond that chemical potential and until m O = À2.60 eV, which is just beyond the conditions in which the phase transition from Fe 2 O 3 to Fe 3 O 4 takes place, the most stable surface is the stoichiometric one. At lower values of m O = À3.00 eV (but still within the conditions in which Fe 3 O 4 is the most stable phase), the (G = À3) reduced surface is the favoured system, until reaching the conditions where the reduced bulk phase of FeO is the most stable oxide. A recent publication by Nie et al. 100 reports that the Fe 3 O 4 (001) surface is oxidised under exposure to 4 Â 10 À9 bar of oxygen at 923 K. They used low-energy electron microscopy (LEEM) and Raman spectroscopy to prove that Fe 3 O 4 grows at the expenses of Fe ions migrating from the bulk towards the surface. The Fe ion vacancies in the bulk, in turn, transform it into a-Fe 2 O 3 (hematite), which is the equilibrium iron oxide phase at the temperature and pressure of the experiment. Our results hence agree well with these experimental findings, although they are close to the limit in which the stoichiometric surface is the most stable one. The experimental conditions described above correspond to m O = À1.75 eV according to eqn (6), a value at which Fe 3 O 4 (001) surface is prone to suffer oxidation (see Fig. 6(b)), by adding one O atom per surface unit cell around the protruding Fe A . Our results, however, agree partially with those reported in a DFT study by Pentcheva et al. 29 as those authors found that the modified nonstoichiometric polar bulk-like Fe 3 O 4 (001) surface (Fe B -O layer) is the most stable under any chemical potential. However, the surface proposed by Pentcheva et al. is a generic oxidised (001) surface, created from a bulk-like termination, whereas our surface is gradually oxidised or reduced. However, regardless of terminations and reconstructions, we also predict our non-dipolar surface to be (G = +3 and +1) oxidized up to m O = À1.85 eV, but from this value of chemical potential onwards, our results predict a gradual reduction, which no longer agrees with the work by Pentcheva et al. as they predict the same oxidized surface for any value of m O .
The redox behaviour of the Fe 3 O 4 (111) surface is shown in Fig. 6(c). It indicates that the redox properties of the (111) surface are similar to the (001) surface, although the oxidized character extends to lower chemical potentials. The surface tends to be (G = +3) oxidised under the condition where Fe 2 O 3 is the most stable phase and up to m O = À2.45 eV, which is within the region where Fe 3 O 4 is the thermodynamically most stable iron oxide. From here, the surface loses two O atoms for a very short range of chemical potential, until m O = À2.95 eV, from where the surface loses a further two more O atoms, becoming now reduced up to m O = À3.25 eV. From this final value of the chemical potential, the surface becomes reduced (G = À3). In a previous DFT + U study, 42 Kiejna et al. studied the redox properties of the Fe 3 O 4 (111) surface. They only studied the non-stoichiometric dipolar bulk-like terminations and found that the Fe A1 surface, which corresponds with a generic oxidized one, is the most stable one up to m O = À2.6 eV, from which point their surfaces started to reduce gradually. Although we cannot make a direct comparison of our results due to the different terminations considered in both works, owing to our gradual redox processes, our results show the same trend.
Comparing both (G = +3) oxidised surfaces, the (111) is lower in energy than the (001) for the whole m O scale considered. Therefore the (111) surface will remain oxidised even at m O where Fe 3 O 4 (001) is not. We suggest then, from Fig. 6, that under the conditions in which the bulk Fe 3 O 4 is the most stable oxide, the phase transformation of the reduction of Fe 3 O 4 towards FeO will start from the (001) surface. On the other hand, Fe 3 O 4 oxidation to Fe 2 O 3 would take place initially on the (111).

Conclusions
In this paper, we have modelled three different surface orientations of Fe 3 O 4 crystals by using DFT methods within the GGA + U approximation. We have investigated the stabilities of their multiple (reconstructed) non-dipolar stoichiometric surface terminations and studied the redox properties of the most prominent surfaces. We have modified the redox conditions by creating O vacancies or adding O atoms to the most stable non-dipolar stoichiometric surface termination, under a wide range of chemical potentials, including ambient conditions and those conditions where bulk Fe 3 O 4 is the thermodynamically most stable oxide. In the initial stages of oxidation, the excess O atoms form bridging structures with the Fe ions at the surface, and in particular the Fe ions protruding from the surface. We found that some oxidised (non-stoichiometric) structures relaxed in such a way that it broke the bulk-like termination.
We conclude that the Fe-terminated (001) and (111) planes are the most stable Fe 3 O 4 surfaces, in agreement with previous experiments as shown by STM images. The equilibrium morphology of Fe 3 O 4 was found to be cubic with truncated corners, which means that (001) and (111) are the main surfaces exposed in the crystals. Although both (001) and (111) surfaces will be oxidized under ambient conditions, both surfaces suffer a gradual reduction, that starts at lower chemical potentials for the (001) surface including the stoichiometric plane.
The reduction of the (001) and (111) surfaces is thermodynamically favourable at the low end of the m O values in the region where Fe 3 O 4 is the most stable oxide. We found that, in both cases, the O vacancies are likely to migrate towards the bulk, thereby changing the phase structure.