Controversies on the surface termination of α-Fe2O3 (0001) focus on its surface stoichiometry dependence on the oxygen chemical potential. Density functional theory (DFT) calculations applying the commonly accepted Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional to a strongly correlated system predict the best matching surface termination, but would produce a delocalization error, resulting in an inappropriate bandgap, and thus are not applicable for comprehensive hematite system studies. Besides, the widely applied PBE+U scheme cannot provide evidence for existence of some of the successfully synthesized stoichiometric α-Fe2O3 (0001) surfaces. Hence, a better scheme is needed for hematite DFT studies. This work investigates whether the strongly constrained and appropriately normed (SCAN) approximation reported by Perdew et al. could provide an improved result for the as-mentioned problem, and whether SCAN can be applied to hematite systems. By comparing the results calculated with the PBE, SCAN, PBE+U, and SCAN+U schemes, we find that SCAN and SCAN+U improves the description of the electronic structure of different stoichiometric α-Fe2O3 (0001) surfaces with respect to the PBE results, and that they give a consistent prediction of the surface terminations. Besides, the bulk lattice constants and the bulk density of states are also improved with the SCAN functional. This study provides a general characterization of the α-Fe2O3 (0001) surfaces and rationalizes how the SCAN approximation improves the results of hematite surface calculations.
Iron oxides are abundant on earth and hence deserve vigorous exploitation of their availability. Hematite (α-Fe2O3) 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,6 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 α-Fe2O3 (0001) surface is one of the most intractable problems remaining unsolved to date.15–17 As commonly accepted, the stoichiometry of α-Fe2O3 (0001) surfaces depends on the oxygen chemical potential μO.18 Only with an appropriate μO can the α-Fe2O3 (0001) surfaces be stable. The complicated relationship between the different stoichiometric α-Fe2O3 (0001) surfaces and the oxygen chemical potential has made the α-Fe2O3 (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 α-Fe2O3 (0001) surfaces has been confirmed, the monophase synthesis has not been well developed yet,15 so that the so-called biphase terminations are studied.19 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 α-Fe2O3 (0001) surface24–26 and the hematite/water interface properties,27 which influence photoelectrochemical water splitting.
The stoichiometry of α-Fe2O3 (0001) was once believed to include two stable structures.28 One is the Fe—O3—Fe—R termination, and the other is the O3—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 termination23 of the O=Fe—O3—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.29 Later, Lemire et al. presented evidence that such a double bond termination indeed exists for the α-Fe2O3 (0001) surface.30 By conducting IRAS, STM, and first-principle studies, Lemire et al. confirmed that the ferryl termination of α-Fe2O3 (0001) coexists with the Fe—O3—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.31 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 α-Fe2O3 (0001) properties. However, as is mentioned, the researchers studied the surface termination problem of the α-Fe2O3 (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 limits23 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 α-Fe2O3 (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.34 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.41 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),42 PBE+U, SCAN, and SCAN+U to the α-Fe2O3 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.
II. COMPUTATIONAL METHODOLOGY
Spin-polarized density functional theory (DFT) calculations are performed using the VASP 5.4.4 code43 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,44 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 3p6d74s1, where the 3p semicore electrons are also taken as valence electrons and 2s22p4 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 functional2 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−5 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 TN = 955 K, and slight canting happens above the Moring temperature TM = 260 K, which is however beyond the discussion of this work. On-site Coulomb corrections are applied with Ueff = 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 calculations21,45,46 with Ueff 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 Carter41 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 method47 and the Blöch corrected tetrahedron method48 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—O3—R, Fe—O3—Fe—R, O3—Fe—Fe—R, O2—Fe—Fe—R, O1—Fe—Fe—R, and O=Fe—O3—R (ferryl) terminations. The different stoichiometric α-Fe2O3 (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—O3—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,49
Eslab stands for the total energy of each slab, Ebulk for the total energy of the corresponding bulk structures, and A for the surface area of the slabs; NFe and NO stand for the number of Fe and O atoms in the slab; is the hematite chemical potential, i.e., the total energy of per Fe2O3 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 is the computed total energy of an oxygen molecule, and the second term is the binding energy correction, where Eb,exp and Eb,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 Fe2O3 would be reduced to Fe3O4,
refers to the total energy of per Fe3O4 formula unit.
The relative oxygen chemical potential μO,rel with respect to the oxygen molecule falls within the following interval:
III. RESULTS AND DISCUSSION
A. Bulk properties
By utilizing PBE, PBE+U, SCAN, SCAN+U, and HSE (12%), we computed the bulk lattice constants and the bandgaps for α-Fe2O3, 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%).2 SCAN and SCAN+U perform better on balancing the structural parameters and the bandgap than PBE and PBE+U.
|.||a = b (Å) .||c (Å) .||Ega (eV) .|
|PBE+U (U=4.0 eV)||5.088(+1.05%)||13.917(+1.21%)||2.16|
|SCAN+U (U=1.8 eV)||5.022(−0.26%)||13.740(−0.05%)||2.08|
|.||a = b (Å) .||c (Å) .||Ega (eV) .|
|PBE+U (U=4.0 eV)||5.088(+1.05%)||13.917(+1.21%)||2.16|
|SCAN+U (U=1.8 eV)||5.022(−0.26%)||13.740(−0.05%)||2.08|
The calculated bandgaps are derived from high symmetric k-mesh calculations.
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 α-Fe2O3 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” insulator53,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.2 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.
B. Stability of the Fe2O3 (0001) surface
Stability of surface terminations is one of the most fundamental problem for surface science, and the α-Fe2O3 (0001) surface is a typical example. To check the performance of SCAN and SCAN+U on determining the stability of various α-Fe2O3 (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—O3—Fe—R and O3—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—O3—R and Fe—O3—Fe—R results are possible terminations for the (0001) surface, and the inexistence of ferryl termination is inconsistent with experiments2 [see Fig. 5(b)]. SCAN and SCAN+U present similar judgment that the α-Fe2O3 (0001) surface terminates with either the Fe—O3—Fe—R or the O3—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 α-Fe2O3 (0001) surface,30 Zandi and Hamann demonstrated that the FeIV=O group is the key intermediate in the photoelectrochemical water oxidation reaction, implying that the intermediate is unstable.55 Note that ferryl and O3—Fe—Fe—R in the SCAN+U calculations have very close surface free energies at the oxygen chemical potential, where Fe—O3—Fe—R transforms to O3—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.,56 we calculated the domain for μOpoor to range from −1.45 eV to −1.19 eV according to
C. Electronic properties of the Fe2O3 (0001) surface
Three terminations, i.e., commonly accepted as stable Fe—O3—Fe—R, O3—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—O3—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 O3—Fe—Fe—R, an abnormal magnetic moment for the subsurface Fe (Fe2) is predicted, which may suggest a low-spin state of Fe2 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.58 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 Fe2O3 magnetization calculations, SCAN is better in computing the magnetic properties than PBE.
|.||Fe—O3—Fe—R .||O3—Fe—Fe—R .||O=Fe—O3—R .|
|.||Fe—O3—Fe—R .||O3—Fe—Fe—R .||O=Fe—O3—R .|
The experimental Bohr magnetic moment of Fe ions in hematite bulk is around 4.2 μB, as reported in Ref. 60.
Fe—O3—Fe—R is the most common termination of the α-Fe2O3 (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—O3—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 O3—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 Fe2. 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 α-Fe2O3 (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—O3—Fe—R and O3—Fe—Fe—R terminations dominate the (0001) surface. The surface free energies of the Fe—O3—Fe—R, O3—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 α-Fe2O3 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).