Effect of cation configuration and solvation on the band positions of zinc ferrite (100)

Zinc ferrite ZnFe2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{2}$$\end{document}O4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{4}$$\end{document} belongs to the spinel-type ferrites that have been proposed as photocatalysts for water splitting. The electronic band gap and the band edge positions are of utmost importance for the efficiency of the photocatalytic processes. We, therefore, calculated the absolute band energies of the most stable surface of ZnFe2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{2}$$\end{document}O4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{4}$$\end{document}, the Zn-terminated (100) surface at self-consistent hybrid density functional theory level. The effect of Fe- and Zn-rich environments, cation exchange as antisite defects and implicit solvation on the band positions is investigated. Calculated flat band potentials of the pristine surface model ranges from -0.9\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-0.9$$\end{document} to -0.8\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-0.8$$\end{document} V against SHE in vacuum. For Zn-rich (Fe-rich) models this changes 0.3–0.9 (0.0–0.7) V against SHE. Fe-rich models are closest to the experimental range of reported flat band potentials. Solvent effects lower the calculated flat band potentials by up to 1.8 eV. The calculated band gaps range from 1.5 to 2.9 eV in agreement with previous theoretical work and experiment. Overall, our calculations confirm the experimentally observed low activity of ZnFe2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{2}$$\end{document}O4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{4}$$\end{document} and its dependence on preparation conditions.


Introduction
In the ongoing need and search for fuel sources that are neither fossil nor nuclear, hydrogen is considered as a promising energy carrier. One route to hydrogen from non-fossil sources is solar-powered water splitting. ZnFe 2 O 4 is a spineltype ferrite that has been proposed as photocatalyst material for water splitting [1,2]. Spinel ferrites MFe 2 O 4 crystallize in space group 227 ( Fd3m ) [1,3,4]. The distribution of Fe and M cations over octahedral and tetrahedral lattice sites is expressed by the formula [M 1−x Fe x ] tet [M x Fe 2−x ] oct O 4 with tet and oct denoting the 8a and 16d Wyckoff positions, 1 3 respectively. x represents the degree of inversion which can range between 0 and 1 corresponding to normal and fully inverse spinels. In general, ZnFe 2 O 4 is considered a normal spinel [5], which has been shown in a number of previous experimental [6][7][8][9][10][11][12] and theoretical studies [2,[13][14][15][16][17][18][19][20][21].
In silico methods allow the prediction of electronic band gaps and band edge positions, which are difficult to obtain experimentally. A suitable band gap ( ∼2 eV) is one requirement for water splitting. It allows to overcome the reaction enthalpy (1.23 eV) and the additional overpotentials of the associated redox reaction [22][23][24]. However, not only a suitable band gap, but also the positions of the highest occupied, HOCO, and lowest unoccupied band, LUCO -which approximate the valence band maximum, VBM, and conduction band minimum, CBM-relative to the redox potential of the respective donor/acceptor pair are essential for efficient photocatalysis [25]. For comparison of band positions, a common reference such as the vacuum or a standardized electrode potential has to be employed. Band positions of bulk structures have no physical meaning because of an undefined reference level of the one-particle levels in 3D periodic calculations [26]. This makes 2D-periodic surface calculations crucial for the investigation of suitable band positions.
Experimental information on the surface structure of spinels is scarce and few theoretical studies of the surfaces of ZnFe 2 O 4 exist [2,[27][28][29][30][31][32]. Low-index spinel surfaces are polar and belong to the unstable type 3 after Tasker's classification [33]. They can be stabilized by defect formation, by reconstruction, through doping, or by adsorption [34,35]. In previous works we investigated the stability of the lowindex surfaces of ZnFe 2 O 4 and electronic properties of the bulk-terminated pristine (100) surface [2,30]. Extending our previous work, we will focus on defect formation and adsorption on the Zn-terminated (100) surface [30]. ZnFe 2 O 4 is essentially a normal spinel with a degree of inversion x of 0.01 ≤ x ≤ 0.34 [36], even more than MgAl 2 O 4 with 0.11 ≤ x ≤ 0.44 [12], depending on the synthesis conditions. For the (100) surface of MgAl 2 O 4 , antisite defects have been observed experimentally [37]. In antisite defects M and Fe cations located on 8a and 16d sites, respectively, are interchanged. For ZnFe 2 O 4 surfaces, depth-depending Zn:Fe ratios were found experimentally [38], also indicating antisite formation. Cation exchange at the surface was also observed for ZnFe 2 O 4 thin films [39] and nanotubes [40].
Therefore we will consider local cation exchange and adsorption of Zn(OH) 2 and Fe(OH) 2 in our theoretical models to mimic these findings. Adsorbed OH species serve as representatives of water species in the inner Helmholtz layer to mimic the basic experimental conditions [38]. Aside from extending the models from our previous work [2] in this way, we will also consider the effect of implicit solvation via VASPSOL [41][42][43]. It is well known, that density functional theory (DFT) on the generalized gradient approximation (GGA) level of theory is insufficient to describe the electronic properties of open-shell transition metal compounds correctly. For the structure optimizations within VASP, we therefore utilized the fast, yet empirical GGA+U method [44], and electronic properties were calculated with a dielectric dependent self-consistent hybrid DFT method [2,20,45]. While experimental findings on ZnFe 2 O 4 surface structures have not been published to the best of the authors' knowledge, experimental flat band potentials of ZnFe 2 O 4 are available [38,[46][47][48][49][50][51] that can be compared with the calculated band positions.
First, we will describe the computational details. Then we will commence with discussing our model, the stability of different cation distributions, their effect on band gap and band positions as well as the effect of implicit solvation on the latter.

Computational details
We used a combined approach of CRYSTAL17 [52] (version 1.0.2) and VASP version 5.4 [53][54][55]. In VASP and CRYSTAL, cartesian coordinates are given in Å (10 −10 m). Therefore this non-SI unit will be employed throughout this work. As in our previous works on ZnFe 2 O 4 [2,20,30], we used the dielectric dependent self-consistent hybrid DFT method sc-PW1PW [45] based on the PW1PW global hybrid functional [56] to calculate the band energies. In CRYSTAL calculations, we used basis sets of triple-zeta quality with polarization quality optimized for solid-state calculations (pob-TZVP-rev2) [57]. This combination, denoted pob/ sc-PW1PW in the following, was employed to converge the model thickness in terms of the frontier band energies against the vacuum reference. For details on program settings the reader is referred to our last study [30]. Since bulk ZnFe 2 O 4 is reported to show antiferromagnetic behavior below 10 K [9], all surface models used antiparallel spin settings resulting in a net magnetic moment of zero. Up spins will be denoted as and down spins as . Having determined the necessary surface thickness, we used VASP to investigate the effects of implicit solvation via VASPSOL [41][42][43]. Surface models used 4 ×4× 1 k-points and the bulk was calculated with 4 ×4× 4 k-points. For the core electrons we employed the projector-augmented wave (PAW) [58,59] method using PBE.54 POTCAR files denoted Zn, Fe and O and a cutoff energy of 450 eV. Structure optimizations with hybrids in VASP are very costly. Therefore we utilized the more empirical but fast GGA+U method [44] combined with the PBE functional [60] for optimizations and then used sc-PW1PW for the electronic property calculations. We applied a uniform Hubbard-U of 4.5 eV for the Fe-d electrons following literature suggestions [13,15]. Since 1 3 VASP does not permit optimizing the bulk with fixed cell shape while relaxing the volume and atomic positions, we manually varied the lattice parameter a and optimized the atomic positions. To determine the optimized lattice parameter we used a second-order polynomial as a fit function. We used surface models optimized with PBE+U for sc-PW1PW single point calculations with and without implicit solvation, denoted as sc-PW1PW and sc-PW1PW/LSOL in the following. For the hybrid single points, PW91+U wavefunctions were used as a starting point [61]. In LSOL, water was chosen as solvent. The relative permittivity EB_K was set to 78.4, and the Debye length LAMBDA_D_K was set to the standard value 3.0 Å (0.3 nm). For surface models treated with VASP, the vacuum distance between two slabs has to be converged to prevent non-physical long-range interactions. 20 Å (2 nm) proved sufficient in the present case. Band positions from VASP calculations have to be corrected by the mean electrostatic potential in the vacuum center, which was obtained via the local potential without exchange contributions (LVHAR=T). Since we are using neutral models, the Fermi shift needed for one-electron levels obtained with LSOL is equal to this potential [41]. Furthermore, we set opposing spins for symmetry equivalent atoms, listed in Table S3 in the Supporting Information, to achieve perfectly antiferromagnetic behavior in the surface models.

Convergence tests
First, we checked the comparability between the electronic properties obtained with CRYSTAL and VASP. In Table S1 calculated band gaps and lattice parameters obtained with CRYSTAL-pob/sc-PW1PW and VASP-PBE+U as well as VASP-sc-PW1PW are given. Using the optimized CRYS-TAL-pob/sc-PW1PW bulk structure, a band gap of 3.04 eV is obtained with CRYSTAL-pob/sc-PW1PW, very close to the VASP-sc-PW1PW value of 3.10 eV. If the PBE+U optimized structure is used, VASP-sc-PW1PW gives a band gap of 3.13 eV, still in good agreement with the other results. This ties in well with our previous theoretical results [2,20] [67]. Using PBE+U(4.5 eV) the optimized lattice parameter, 8.519 Å (0.8519 nm) is larger than the dd-sc-PW1PW/pob value, but matches reported experimental values of 8.523 Å (0.8523 nm) [32].
Since CRYSTAL allows for cost effective hybrid optimizations, we used pob/sc-PW1PW to converge the band positions with the number of layers in the slab models. We use the same symmetric defect models and additional ghost atoms above the surface as in our previous study [30]. As shown in Figure S1 in the Supplementary Information, a surface model with 17 atomic layers or 4 ZnFe 2 O 4 layers is sufficient for converging the band positions within 0.1 eV. An error bar of ±0.1 eV will therefore be adapted for all band gaps and band positions of the surface models. To increase the number of bulk-like layers in the middle of the slab, we selected a larger model with 21 atomic layers despite its slightly larger deviations from the larger models. This model corresponds to 5 ZnFe 2 O 4 layers and 70 atoms. All surface calculations with VASP were performed with this model. Figure S2 shows that increasing the vacuum layer from 20 Å to 25 or 30 Å (2.0, 2.5, 3.0 nm) has no significant effect on the band positions. It also demonstrates a good match between CRYSTAL-pob/sc-PW1PW band positions and VASP-sc-PW1PW ones based on PBE+U structure optimizations. As expected, the band gaps calculated with PBE+U increase with increasing U, as demonstrated in Figure S3. However, it is not possible to match the PBE+U band positions and the center of the band gap with sc-PW1PW results. The PBE+U approach is therefore not further considered.

Surface models
As mentioned before, all low-index spinel surfaces, including the tetrahedrally terminated (100) plane, are polar and correspond to type 3 according to Tasker's classification [33]). We followed the approach of our previous works [2,30] for surface model generation, using defective but stoichiometric and symmetric surface models to eliminate artificial dipole moments. In depth discussion of model construction can be found in our previous work [30]. Aside from investigating a defective, stoichiometric and symmetric model, which represents an unrealistic structure in aqueous environment, we investigated the adsorption of Zn(OH) 2 at the Zn vacancy sites. Formally this corresponds to symmetric and defect-free surface models with adsorbed OH groups, a 1 3 situation which is more probable in basic solution. To investigate antisite-defects, we exchanged Zn and Fe cations close to the surface layers. To mimic the experimentally observed Fe-rich environments described by Warfsmann et al. [38], we also adsorbed Fe(OH) 2 . The model without adsorption is denoted as pristine (see Fig. 1). As an example for adsorption, the normal surface model with adsorbed Zn(OH) 2 is shown in Fig. 2. Not all cation-exchanged structures could be converged during the optimization, we therefore only present a small subset of all possible Fe/Zn configurations. The metal ion distribution and initial atomic spins of all models that could be optimized are listed in Table 1. Figures S4-S6 in the Supplementary Information depict the optimized surface models. All structures were visualized with VESTA [70].

Surface stability
In previous studies we calculated the energy required for the formation of an antisite defect in the bulk as ADFE=13 kJ/ mol [30]. Antisite defect formation in the bulk is equivalent to inversion. To compare inversion in the bulk and on the surface, we calculated the relative stability of the slab models with cation exchange and the respective normal model (see Table 2). Due to their differing stoichiometry it is not possible to compare the order of stability between Zn-rich (Zn(OH) 2 ) and Fe-rich (Fe(OH) 2 ) surface models. For the Zn(OH) 2 -anti2 model, we investigated the exchange of a Zn atom in the topmost layer with a fivefold coordinated Fe 3+ in the topmost Fe-containing layer. Surprisingly, this antisite defect is favored by 2 kJ/mol in comparison with the normal cation distribution. During the relaxation the tetrahedral Fe 3+ distorts from its initial position to increase the coordination number from 3 to 4. Fe 3+ also relaxes towards the surface, compare Fig. S4(b). When exchanging the same fivefold coordinated Fe 3+ with a lower-lying fourfold coordinated Zn, model Zn(OH) 2 -anti1, the exchange is unfavorable by 16 kJ/mol. This is similar to the calculated bulk value for antisite defect formation, 13 kJ/mol [30].
For the Fe-rich models with cation exchange that we could converge, antisites are less stable than the normal cation distribution. The ADFE for the exchange of a bulklike fourfold coordinated Zn 2+ with an adsorbed threefold coordinated Fe 2+ in Fe(OH) 2 -anti1 is 46 kJ/mol. Inversion of Fe and Zn in deeper layers in Fe(OH) 2 -anti2 results in further destabilization, ADFE = 66 kJ/mol. Exchange of the outermost Zn atom with a lower-lying Fe-atom (anti3 model) is also highly unfavorable with ADFE = 99 kJ/mol. Table 3 and Fig. 3 show the band gaps and band edge positions of our surface models with and without solvent interaction. As found in theoretical studies on other semiconductor surfaces [71], solvation leads to a positive shift of the band positions. The band gaps of the surface models are 0.5-0.2 eV smaller than the bulk band gap, which is also expected.

Electronic properties
On the surface, some atoms have a smaller coordination number in comparison with the bulk resulting in a reduced ligand field and a smaller HOCO-LUCO-splitting [71]. All calculated band gaps shown in Table 3 are therefore smaller than the bulk value (3.04 eV, obtained with CRYSTAL-sc-PW1PW, see Table S1). Generally, increasing the Zn content increases the band gap, while increasing the Fe content lowers the band gap with respect to a stoichiometric model. This trend has also been found in the experimental work by Warfsmann et al. [38]. While the anti3 model is a clear outlier with band gaps of 1.5±0.1 eV, the band gaps of the Fe-rich models are 1.5-2.0 eV without solvent correction and 1.9-2.2 eV with solvent correction. This lowering of the band gap can be explained by the projected densities of states (see Figs. S7-S10 in the Supplementary Information).  The adsorption of Fe(OH) 2 leads to mid-gap Fe-states that lower the band gap, independent of the cation distribution. In all surface models considered in this work, the LUCO/ CBM is dominated by Fe-states. For the pristine and Znrich models, the HOCO/VBM mainly consists of O-states. The investigated cation exchange of the Zn-rich models has a small effect on the band gaps and a moderate effect on the band positions. With 0.1 eV, the solvent effect on the band gaps of Zn-rich models is negligibly small. This ties in well with the findings of Granone et al. that the degree of inversion, and thus cation exchange, has no effect on the band gap. At variance, the band edge positions are strongly affected.
They increase by up to 1.8 eV for the Zn(OH) 2 -anti1 model, see Fig. 3. However, in solution the difference in band positions of the various Zn-rich antisite models is  In comparison with the pristine model, the effect of the solvent corrections on the band positions of the Zn-and Fe-rich models is much more pronounced. We attribute this to the terminating polar OHgroups which are directly exposed to the solvent.

Comparison with flat band potentials
The measured flat band potential E fb of a semiconductor is interpreted as CBM or VBM, depending whether it is of n-type or p-type, respectively. ZnFe 2 O 4 is an n-type semiconductor [48]. For n-type semiconductors, the flat band potential is often assumed to lie 0.1-0.2 eV below the conduction band [72]. This is a simplification, since the difference between the conduction band edge and the flat band potential depends on the charge carrier density and the effective mass of the electrons in the material [73]. Correspondingly, it also depends on the material's conductivity. For low-conductivity materials, the difference between the flat-band potential and the conductionband is significant [74], but it is negligible for doped or highly conductive materials. When comparing our calculated LUCO/CBM with experimental flat band potentials, listed in Table 4, we follow the approach of Arimi et al. [47] assuming E fb ≈ E CB . We only considered flat band p ri s ti n e + Z n (O H ) 2 − a n ti 1 − a n ti 2 + F e (O H ) 2 − a n ti 1 − a n ti 2 − a n ti 3   [46] 0.64 RHE (pH 13) −0.13 [38] Zn:Fe 0.52 0.75 RHE (pH 13.6) −0.05 [48] 0.84 RHE (pH 13.3) 0.06 [49] 0.82 RHE (pH 13.3) 0.04 [50] −0.72 NHE −0.72 [51] LB potentials from the work of Arimi et al. [47] that were taken from hematite-free samples identified via Raman measurements. Since our change in Fe:Zn ratio is overall small and constant, we solely used the flat band potential of the stoichiometric sample in the work of Warfsmann et al. [38]. All other experimental flat band potentials used for comparison are listed in Table 4. It also lists the ranges of the calculated LUCO positions, with individual LUCO positions listed in Table S2 in the SI. For the conversion of the calculated band energies vs. vacuum we used the standard hydrogen electrode, SHE, energy of 4.44 eV published by Trassati [75]. The SHE is also often referred to as normal hydrogen electrode, NHE [76,77]. Often potentials are given with respect to the reversible hydrogen electrode, RHE, which is dependent on the pH-value via E NHE =E RHE − 0.059 ⋅ pH [76,78]. However, they are not suitable for hydrogen production as they lie too high. This conclusion has also been drawn from experimental results [2]. Additionally, the photocatalytic performance of ZnFe 2 O 4 has been reported as below expectations [68,79]. When applying the solvent correction, the LUCO/CBM is lowered with respect to the SHE, however the corresponding HOCOs now lie too high for oxygen production. The LUCO positions of the Zn-rich models exceed the upper end of the experimental values without solvent correction at 0.3 to 0.9 V against SHE but are also too high for hydrogen production. The solvent correction evens out the effect of the cation exchange on the band positions and the HOCO/VBM positions could be promising for solar water splitting. Despite the deviation of the calculated LUCO values from the E fb range, they provide a possible explanation for the range of experimental flat band potentials, which may be due to different antisite defects, thus partial inversion. It has to be noted, that the conductivity of ZnFe 2 O 4 increases by two orders of magnitude when the degree of inversion exceeds 0.14 [69], also supporting an influence of cation exchange on the flat band potential. Regarding the deviation between the experimental and calculated values we identified two items for further investigation. First, we employed the dielectric constant of bulk water, whereas the dielectric constant of water at interfaces is significantly smaller [80]. Second, we explicitly adsorbed OH and no water molecules which could transfer electron density to the surface thus increasing the band energies [71]. This will be taken into account in future studies.

Summary and conclusion
We investigated the electronic properties of the (100) surface of normal ZnFe 2 O 4 , the effects of cation exchange via antisite formation, and of Zn-or Fe-rich conditions by Zn(OH) 2 and Fe(OH) 2 adsorption with self-consistent hybrid DFT and implicit solvent models. Both adsorption and solvent have an effect on the band gap and absolute band edge positions. Band edge positions of Zn-and Ferich models in vacuum reflect the experimentally reported flat band potentials, while the solvation correction deteriorates the agreement. The effect of the solvent correction on the band gap is not uniform and more pronounced for the more polar models containing OH-adsorption. Zn(OH) 2 adsorption increases the band gap in comparison to the pristine model and downshifts the band positions. For Zn-rich models, antisite formation in the topmost layer is slightly favored and in the bulk-like layers it is as endothermic as in the bulk. For Fe-rich models, antisite formation is thermodynamically unfavorable. Fe(OH) 2 adsorption mimicking Fe-rich conditions narrows the band gap due to the formation of iron mid gap states. Here the cation distribution has a smaller effect on the band position than for Zn-rich models. Overall, the cation distribution does influence the band positions and therefore, among other factors, also influences the flat band potential, which may explain the large scattering of experimental results.
Our modeling of explicit solvation through a small number of adsorbed OH-groups might not be a sufficient representation of the water adsorption of real ZnFe 2 O 4 surfaces. This will be addressed in future model studies where additional molecular water adsorption will be taken into account. Without solvent effects, neither the most stable Zn-rich model nor the most stable Fe-rich model have suitable band positions for solar hydrogen production. We therefore conclude that ZnFe 2 O 4 is not an active catalyst material in solar hydrogen production.

3
computer time with LUIS and the University of Bonn for providing computer time on the Bonna cluster.
Funding Open Access funding enabled and organized by Projekt DEAL.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.