Improved description of hematite surfaces by the SCAN functional

Iron oxides are abundant on earth and hence deserve vigorous exploitation of their availability. Hematite (α-Fe₂O₃) has been applied in various domains on account of its outstanding features, including stability, nontoxicity, a narrow optical bandgap of around 2.0–2.2 eV,^1,2 and desirable cost efficiency. Diverse fields of research have been studying this material in the recent years. For industrial applications, hematite has been studied in domains including, but not limited to, ion-batteries,^3–5 biosensing,⁶ wastewater treatment,^7,8 and photoelectrochemistry,^9–14 attributing to its unique value to human society.

Despite the advantages of this material, there are still many problems that we have to solve to promote the understanding of hematite. The surface termination problem of the α-Fe₂O₃ (0001) surface is one of the most intractable problems remaining unsolved to date.^15–17 As commonly accepted, the stoichiometry of α-Fe₂O₃ (0001) surfaces depends on the oxygen chemical potential μ_O.¹⁸ Only with an appropriate μ_O can the α-Fe₂O₃ (0001) surfaces be stable. The complicated relationship between the different stoichiometric α-Fe₂O₃ (0001) surfaces and the oxygen chemical potential has made the α-Fe₂O₃ (0001) surfaces ambiguous. From the experimental point of view, the surface terminations of synthesized hematite are hard to determine due to various factors. For example, even though the existence of various stoichiometric α-Fe₂O₃ (0001) surfaces has been confirmed, the monophase synthesis has not been well developed yet,¹⁵ so that the so-called biphase terminations are studied.¹⁹ Theoretically, density functional theory (DFT) calculations with different schemes would produce divergent conclusions.^20–23 More efforts deserve to be paid to the surface termination problem, which is very fundamental for further investigation of hematite’s properties and applications. In particular, surface states strongly influence photoinduced reactions, such as photoelectrochemical water splitting. As mentioned in the present studies, the surface termination problem is very important for the molecule or cluster adsorption studies on the α-Fe₂O₃ (0001) surface^24–26 and the hematite/water interface properties,²⁷ which influence photoelectrochemical water splitting.

The stoichiometry of α-Fe₂O₃ (0001) was once believed to include two stable structures.²⁸ One is the Fe—O₃—Fe—R termination, and the other is the O₃—Fe—Fe—R termination (“R” stands for the repeating units accorded with the periodicity). Later in 2004, Bergermayer and Schweiger theoretically pointed out the existence of “ferrite-like” surface termination²³ of the O=Fe—O₃—R structure, and this is also called “ferryl termination.” Meanwhile, Rohrbash et al. predicted that the ferryl termination exists with Generalized Gradient Approximation (GGA), but not with GGA+U.²⁹ Later, Lemire et al. presented evidence that such a double bond termination indeed exists for the α-Fe₂O₃ (0001) surface.³⁰ By conducting IRAS, STM, and first-principle studies, Lemire et al. confirmed that the ferryl termination of α-Fe₂O₃ (0001) coexists with the Fe—O₃—Fe—R termination.

The inexistence of the ferryl termination in GGA+U implies deficiencies in the early DFT studies on hematite. Hematite is a typical transition metal system and is strongly correlated. Simple DFT schemes without corrections produce the delocalization error, which is famous for resulting in the underestimation of bandgaps.³¹ The most popular scheme to deal with this issue is to conduct an on-site Coulomb correction.^32,33 By applying a specified Hubbard U term on Fe atoms, the electronic structure of hematite bulk had been improved, and hence, it was believed that the on-site Coulomb correction would improve the prediction of the α-Fe₂O₃ (0001) properties. However, as is mentioned, the researchers studied the surface termination problem of the α-Fe₂O₃ (0001) with on-site Coulomb corrections, and acquired the result that the ferryl termination would not happen in the domain between the poor and rich limits²³ of the oxygen chemical potential. The inexistence of ferryl termination predicted with the on-site Coulomb correction is inconsistent with the experimental studies, implying that the on-site Coulomb correction may still be insufficient to address the surface termination problem. Hence to characterize theoretically the termination of the α-Fe₂O₃ (0001) surface, a more accurate description of the exchange-correlation functional is necessary.

In 2015, Sun and Perdew et al. reported that they had constructed a new meta-GGA that satisfies all the known possible exact constrains and appropriate norms, named as the Strongly Constrained and Appropriately Normed (SCAN) approximation.³⁴ Scientists have put significant efforts into testing of SCAN with various systems, confirming its good performance.^35–40 The descriptions of d-electrons with the SCAN approximation are still limited, because meta-GGAs are still avoiding full nonlocality. It is reported by Sai Gautam and Carter that SCAN underestimates the bandgaps of transition metal semiconductors and misestimates the ground state energies.⁴¹ Besides, they conducted SCAN+U studies, and found that SCAN+U is suitable for those systems.

In this work, we apply Perdew-Burke-Ernzerhof (PBE),⁴² PBE+U, SCAN, and SCAN+U to the α-Fe₂O₃ bulk and its different (0001) stoichiometric surfaces, and investigate the lattice constants, bulk density of states (DOS), surface termination, surface magnetic moments, and surface density of states. By comparing the results obtained with the different schemes, we find that SCAN and SCAN+U perform well in the hematite system calculations.

Spin-polarized density functional theory (DFT) calculations are performed using the VASP 5.4.4 code⁴³ with the projected augmented wave (PAW) method. The applied exchange-correlation functional is either the Generalized Gradient Approximation of the Perdew-Burke-Ernzerhof functional (GGA-PBE) or SCAN. The energy cutoff for the plane-wave basis is set to 520 eV for the full bulk relaxation,⁴⁴ and 400 eV for the other calculations with lattice constants fixed, including electronic structure calculations and slab relaxations. The pseudopotential of Fe atoms is constructed with the valence electron configuration of 3p⁶d⁷4s¹, where the 3p semicore electrons are also taken as valence electrons and 2s²2p⁴ for O atoms. For bulk calculations, the Monkhorst-Pack k-mesh is set to be 5 × 5 × 3, and for slab calculations, it is changed to 7 × 7 × 1. We also apply the Heyd-Scuseria-Ernzerhof (HSE) hybrid functional² with 12% Fock exchange to calculate [abbreviated as HSE (12%)] the bulk properties, where the k-mesh is set to be 4 × 4 × 2 and a corresponding coarse mesh of q-points with a reduce factor of 2. The energy convergence criterion of the electronic self-consistent field is 10⁻⁵ eV, and the structures are relaxed until the residual force on atoms is less than 0.02 eV/Å. It is common knowledge that the antiferromagnetic configuration is the most stable state for hematite below the Néel temperature T_N = 955 K, and slight canting happens above the Moring temperature T_M = 260 K, which is however beyond the discussion of this work. On-site Coulomb corrections are applied with U_eff = 4.0 eV (U = 5.0 eV and J = 1.0 eV) for Fe atoms in the GGA+U calculations and U = 1.8 eV for Fe atoms in the SCAN+U calculations. Both the previous PBE+U calculations^21,45,46 with U_eff of 4 eV and the SCAN+U calculations done here with U of 1.8 eV predicted bulk bandgaps, which agree well with the optical bandgap. This U term differs from the U of 2.7 eV used in the SCAN+U calculations by Sai Gautam and Carter⁴¹ chosen to match the fundamental bandgap. Gaussian smearing with a smearing width of 0.05 eV was applied for most structural relaxations. We decrease the smearing width to 0.02 eV for a few exceptions, in which the magnetic moments on some atoms are low. The Gaussian smearing method⁴⁷ and the Blöch corrected tetrahedron method⁴⁸ were applied for integrating the bulk and surface density of states, respectively. Applying the Blöch corrected tetrahedron method is necessary for the surface electronic structures, because the surface electronic structure is more complex than the bulk structure due to surface states.

The models are built as follows. The hexagonal bulk structure is shown in Fig. 1, where the lattice constant a is equal to b, and the angle included between the axes alongside is 120°. The spin state of this system is antiferromagnetic as is commonly accepted. Regarding the two Fe layers (named as an Fe buckle) separated by the O layers, the spin states alternate with the piling up of the atom layers along the c axis.

Various typical terminations are taken into consideration, listed as Fe—Fe—O₃—R, Fe—O₃—Fe—R, O₃—Fe—Fe—R, O₂—Fe—Fe—R, O₁—Fe—Fe—R, and O=Fe—O₃—R (ferryl) terminations. The different stoichiometric α-Fe₂O₃ (0001) surfaces are constructed by using symmetric (1 × 1) slabs with 16 Å vacuum space. Side views of these surface terminations are shown in Fig. 2. These surfaces are built with the lattice constants inherited from the corresponding bulk. All the slabs are built by cutting the bulk structure except ferryl, which is built by placing an O atom on the Fe—O₃—Fe—R model with the Fe—O bond vertical to the surface Fe atom. All the atoms in each termination are allowed to relax.

The surface free energy γ was calculated by following the framework proposed by Reuter and Scheffler,⁴⁹ 

E_slab stands for the total energy of each slab, E_bulk for the total energy of the corresponding bulk structures, and A for the surface area of the slabs; N_Fe and N_O stand for the number of Fe and O atoms in the slab; $μFe2O3$ is the hematite chemical potential, i.e., the total energy of per Fe₂O₃ formula unit, and μ_O is the oxygen chemical potential, which is variable in the range between two limits.

The rich limit μ_Orich exists, since the PBE functional overestimates the binding energy of the oxygen molecule. An energy correction is required to calculate the oxygen chemical potential. The rich limit of the oxygen chemical potential is determined as follows:

The first term $12EO2$ is the computed total energy of an oxygen molecule, and the second term is the binding energy correction, where E_b,exp and E_b,cal stand for the oxygen binding energy from experiments and calculations, respectively.

The poor limit of the oxygen chemical potential is defined as the chemical potential at which Fe₂O₃ would be reduced to Fe₃O₄,

$μFe3O4$ refers to the total energy of per Fe₃O₄ formula unit.

The relative oxygen chemical potential μ_O,rel with respect to the oxygen molecule falls within the following interval:

By utilizing PBE, PBE+U, SCAN, SCAN+U, and HSE (12%), we computed the bulk lattice constants and the bandgaps for α-Fe₂O₃, as shown in Table I. All the acquired lattice constants agree with the experimental data, with the errors less than 1.21%. Compared to PBE, PBE+U greatly improves the bandgap, but worsens the prediction of the lattice constants. SCAN gives overall better descriptions on both the lattice constants and the bandgap than PBE. SCAN+U further improves the two properties, which are very close to the experimental values, even better than HSE (12%).² SCAN and SCAN+U perform better on balancing the structural parameters and the bandgap than PBE and PBE+U.

The density of states (DOS) is computed with PBE, PBE+U, SCAN, SCAN+U, and HSE (12%), as shown in Figs. 3 and 4. PBE seriously underestimates the bandgap (0.58 eV) and predicts α-Fe₂O₃ as a “Mott-Hubbard” semiconductor, which is inconsistent with the “charge transfer” semiconductor suggested by the experiment.^51,52 SCAN increases the bandgap to 1.37 eV, which, however, is still smaller than the experimental optical bandgap, and it also reduces the hybridization between Fe 3d and O 2p orbitals. The bandgaps (2.16 eV and 2.08 eV) of PBE+U and SCAN+U schemes fall well within the experimental interval of 2.0–2.2 eV, and the “charge transfer” insulator^53,54 is correctly predicted. HSE (12%) gives a bandgap of 2.0 eV, which is close to the value of 1.95 eV reported by Pozun and Henkelman.² The profile of DOS computed with SCAN+U resembles that of HSE (12%), more than that of PBE+U. In a word, SCAN is better than PBE in describing both the structural and electronic properties, and SCAN+U gives better descriptions on the two properties than PBE+U.

Stability of surface terminations is one of the most fundamental problem for surface science, and the α-Fe₂O₃ (0001) surface is a typical example. To check the performance of SCAN and SCAN+U on determining the stability of various α-Fe₂O₃ (0001) terminations, we calculate the surface free energies for these terminations by using Eq. (1). Figure 5 shows the stoichiometry dependence of the surface free energies on the chemical potential μ_O. PBE predicts that Fe—O₃—Fe—R and O₃—Fe—Fe—R terminations dominate the (0001) surface, whereas the ferryl termination occurs in a very small window of the oxygen chemical, which is consistent with previous studies. PBE+U draws a different conclusion that Fe—Fe—O₃—R and Fe—O₃—Fe—R results are possible terminations for the (0001) surface, and the inexistence of ferryl termination is inconsistent with experiments² [see Fig. 5(b)]. SCAN and SCAN+U present similar judgment that the α-Fe₂O₃ (0001) surface terminates with either the Fe—O₃—Fe—R or the O₃—Fe—Fe—R termination, as shown in Figs. 5(c) and 5(d). Whether the ferryl termination exists is a subtle question. Although Lemire et al. reported that the ferryl termination formed on the α-Fe₂O₃ (0001) surface,³⁰ Zandi and Hamann demonstrated that the Fe^IV=O group is the key intermediate in the photoelectrochemical water oxidation reaction, implying that the intermediate is unstable.⁵⁵ Note that ferryl and O₃—Fe—Fe—R in the SCAN+U calculations have very close surface free energies at the oxygen chemical potential, where Fe—O₃—Fe—R transforms to O₃—Fe—Fe—R. These results suggest that the formation of the ferryl termination may strongly depend on the conditions in measurements or reactions.

Based on the experimental conditions reported by Ketteler et al.,⁵⁶ we calculated the domain for μ_Opoor to range from −1.45 eV to −1.19 eV according to

Three terminations, i.e., commonly accepted as stable Fe—O₃—Fe—R, O₃—Fe—Fe—R, and subtle ferryl, are selected to be further investigated for the magnetic moments and density of states. Local magnetic moments of the three terminations are shown in Table II. The local magnetic moments on Fe atoms of Fe—O₃—Fe—R termination are comparable to those in the bulk, indicating the high-spin state of the Fe atom. The magnetic moments converge quickly to the bulk magnetic moments with the Fe atom approaching the center of the slabs. PBE gives the smallest local magnetic moments, and SCAN significantly improves them. PBE+U and SCAN+U present the biggest magnetic moments. For O₃—Fe—Fe—R, an abnormal magnetic moment for the subsurface Fe (Fe₂) is predicted, which may suggest a low-spin state of Fe₂ and the disadvantage of PBE for surface calculations. Previous studies have shown that the application of the SCAN functional may produce erroneous results for some systems.^57,58 For example, PBE was reported to be more accurate than SCAN when conducting Fe magnetization calculations, but less accurate for Co, Ni systems.⁵⁸ In this work, PBE predicts the magnetization moment of around 3.6 µB on Fe ions, whereas the magnetic moments in the hematite bulk were reported to be around 4.0–4.6 µB,^21,59,60 which agree better with the SCAN results. That is, for Fe₂O₃ magnetization calculations, SCAN is better in computing the magnetic properties than PBE.

Fe—O₃—Fe—R is the most common termination of the α-Fe₂O₃ (0001) surface, whose stoichiometry is the same as that of the bulk hematite. The density of states of this termination is shown in Fig. 6. PBE predicts the Fe—O₃—Fe—R termination as a half-metal since the Fermi level crosses through the valence band (VB), whereas the other three functionals predict a semiconductor with surface energy states of Fe d characteristics existing in the bandgap. SCAN gives a qualitative consistent result with PBE+U and SCAN+U, except for the smaller energy band gaps

The density of states of O₃—Fe—Fe—R is shown in Fig. 7. All the four functionals predict the conductivity of this termination to be semiconductive. However, we have found that PBE gives an abnormal magnetic moment on Fe₂. Figure 7 is insufficient to uncover the origin of the phenomenon. To check the origin of the inconsistence between the PBE result and the other three functionals, we further investigate the state contribution of the surface atom layer, as shown in Fig. 8.

Some special peaks in the DOS are selected to check whether the inconsistency is caused by the narrow peak gaps, including two spin-up nearby peaks p1 and p2, a spin-up single peak p3, a spin-down single peak p4, and two spin-down nearby peaks p5 and p6. In the PBE result, p3 and p1-p2 peaks do not break asunder, and the spin-down p4 is located below the Fermi-level. The spin-up peak p3 and spin-down peak p4 below the Fermi-level imply the low-spin state of the surface Fe atoms in the PBE result, whereas the other results consistently predict the high-spin state. The SCAN, PBE+U, and SCAN+U results are more consistent when we investigate the position of the peaks. They all predict that the p1-p2 peaks are below the Fermi-level, and are split from the p3 peak. The p1-p2 peaks in the SCAN result are within the valence band (VB), whereas the peaks are apart from the VB for on-site Coulomb corrected results.

It is clear that the result given by PBE is totally different from those given by the other three functionals. The magnetization moments predicted by PBE is consistent with its DOS. Overall, it can be concluded that SCAN gives qualitatively consistent results with PBE+U and SCAN+U, whereas PBE does not.

The ferryl termination, as a special stoichiometric α-Fe₂O₃ (0001) surface, gives the density of states presented in Fig. 9. PBE predicts the termination to be conductive, with its Fermi-level crossing through the VB. SCAN predicts this termination to be semiconductive with a very small bandgap. The two functionals with the on-site Coulomb correction predict the termination to be semiconductive.

The DOS for the ferryl termination indicates that SCAN can give a qualitatively good description of the electronic structure of the (0001) surface, and that SCAN+U is better for such a strongly correlated system. The improvement by SCAN is moderate, and its insufficient widening of the bandgap can sometimes induce an equivocal consequence. A Coulomb U term applied to the SCAN calculations is hence still suggested for qualitative electronic structure analysis. The SCAN+U results are found to be consistent with the PBE+U results, but the value of U is smaller. Considering that SCAN performs better in the surface termination predictions and oxygen binding energy calculations (−1.34 eV for PBE and – 0.1 eV for SCAN), we suggest applying SCAN+U to strongly correlated systems for comprehensive studies.

In summary, we investigated hematite bulk properties, and surface stability and electronic properties for different stoichiometric (0001) surfaces with PBE, SCAN, PBE+U, and SCAN+U. The bulk bandgap and the local magnetic moment on Fe ions predicted by PBE are poor, compared to the experimental values. SCAN improves the description on the lattice constants, bulk density of states, magnetic moments, and surface density of states and always gives qualitatively consistent results with PBE+U and SCAN+U. By applying the on-site Coulomb corrections on Fe atoms in SCAN+U, all the above properties are improved further.

Although SCAN+U does not predict the existence of the ferryl termination, it gives a phase diagram, which is close to that of PBE. SCAN+U predicts that the Fe—O₃—Fe—R and O₃—Fe—Fe—R terminations dominate the (0001) surface. The surface free energies of the Fe—O₃—Fe—R, O₃—Fe—Fe—R, and ferryl terminations at the phase transition point are close to each other, suggesting the subtle role of ferryl in the surface stability. Overall, SCAN improves the description of the bulk and the (0001) surface properties of α-Fe₂O₃ over PBE, and SCAN+U is the best choice to determine these properties.

This work was supported by the National Natural Science Foundation of China (Grant No. 51876173), Natural Science Foundation of Jiangsu Province (Grant No. BK20190054), the Shaanxi Technical Innovation Guidance Project (Grant No. 2018HJCG-14), Natural Science Basic Research Plan in Shaanxi Province of China (Grant No. 2019JM-400), and China Fundamental Research Funds for the Central Universities. O.P. acknowledges funding of the U.S. Department of Energy (Grant No. DE-SC0014429).