To achieve high electrochemical or photoelectrochemical (PEC) efficiency is a fundamental objective for energy conversion and storage applications. [1][2][3] Transition metal oxides (TMOs) have been the most promising candidates for these applications, such as solar energy conversion and Li-ion batteries. 4,5 Chief among these TMOs is hematite (α-Fe 2 O 3 ), whose chemical abundance, remarkable stability, moderate bandgap, non-toxicity and potential for high solar-to-hydrogen efficiency make it one of the most promising materials to use in PEC devices, batteries, and supercapacitors. [6][7][8][9] However, the performance of Fe 2 O 3 based photoelectrodes is far below its theoretical limit, 6 which hinders their practical applications.
Hence, insights on how to improve properties important for PEC such as carrier conductance, light absorption, and surface reactivity are highly desired.
In particular, strategies to improve carrier conductivity of hematite are most critical for better PEC performance. This is because the carrier conduction of transition metal oxides such as hematite (α-Fe 2 O 3 ) is severely limited by the formation of localized carriers known as small polarons which yield low intrinsic mobility and concentration of carriers. [10][11][12] To overcome these limitations, several efforts have been made to dope hematite by tetravalent ions that yield improved photoelectrochemical performance of hematite photoelectrodes. [12][13][14][15][16][17][18][19][20][21] While the performance can be moderately improved via group IV and XIV dopants, the optimal doping concentration strongly varies with each individual dopant, [16][17][18][19][20][21] thereby requiring extensive experimental testing each time. 22 For example, researchers 17,18 have found that Ti-doped hematite photoanodes had the highest carrier density and photocurrent at a doping concentration of around 0.1%. Meanwhile, several works 19,20,23 have found that optimal PEC performance with Sn-doped hematite photoanodes was achieved at 3% Sn doping concentration. In all of these cases there is a direct correlation between optimizing carrier density and PEC performance; however, the mystery of extremely low optimal doping concentration for certain dopants remains elusive.
Two possible mechanisms could be responsible for the doping bottleneck: compensation by oppositely-charged defects or the clustering of dopants. However, the concentration of intrinsic p-type defects is expected to be negligible in Fe 2 O 3 , 24 which leaves a strong rationale for the clustering of dopants being the cause of low optimal doping concentration.
On the theoretical side, group IV and XIV dopants in Fe 2 O 3 have been previously investigated from first-principles, focusing on their electronic structure, formation energy, and polaron hopping barrier. [25][26][27][28] However, these studies cannot explain the low optimal doping concentration observed experimentally.
In this work, we will reveal the origin of the extremely low optimal doping concentration in Fe 2 O 3 through a joint theoretical and experimental study. We suggest a novel form of dopant clustering in polaronic oxides and conclude its critical role on determining carrier concentration. We begin by detailing our computational methodology, including our proposed model for disentangling the effects of dopant clustering. Next, the electronic structure of isolated and clustered Sn dopant formation is provided, which resembles an electric dipole and quadrupole, respectively. The binding energy of the clustered dopants as quadrupoles is computed to validate their thermodynamic stability. Then, the formation of the theoretically predicted Sn-Sn pairs are confirmed by experimental EXAFS, and their mechanistic origin is unraveled theoretically in terms of electrostatic, magnetostatic and strain effects.
Finally, carrier concentrations of Fe 2 O 3 with and without dopant clustering are computed to elucidate the underlying mechanism of the doping bottleneck. At the end, essential design principles are provided to yield higher conductivity in polaronic oxides for the advancement of energy conversion applications.
All Density Functional Theory (DFT) calculations were carried out using the open source plane-wave code QuantumESPRESSO 29 with ultrasoft pseudopotentials 30 and an effective Hubbard U 31 value of 4.3 eV for Fe 3d orbitals. 25,32 This U value is chosen for its ability to reproduce the bandgap of hematite (∼ 2.21 eV) but also has shown to capture physics of small polarons such as the polaron hopping barrier. 25,32 Plane-wave cutoff energies of 40 Ry and 240 Ry were used for wavefunctions and charge density, respectively. All calculations were performed with a 2 × 2 × 1 supercell (120 atoms) of the hexagonal unit cell with a 2 × 2 × 2 kpoint mesh for integration over the Brillouin zone. A 3 × 3 × 1 supercell was also tested to ensure convergence with supercell sizes (see supercell convergence in SI Table S1). The consistency between these supercell sizes also validates that the present calcula-tions are in the dilute limit and while the Sn at Fe concentrations of the aforementioned supercells are 2.08% and 0.93%, respectively, these concentrations do not yield interactions between dopant periodic images (even for systems with two dopants, see SI Figure S1). The actual concentrations of dopants are determined by evaluating charge neutrality directly from dopant formation energies at a synthesis condition as discussed later in this section. Finally, we note that we use the same U value for systems with dopants. While it is true that changing U will vary the bandgap (in this case the conduction band shifts due to U correction on Fe 3d), it has been shown that the formation energies computed with different U values were very similar and the ionization energies changed little when referenced to the free polaron level instead of the CBM for Fe 2 O 3 . 25
We will later demonstrate that single dopants resemble dipoles, while dopant-pairs resem-ble quadrupoles. Here, we demonstrate our analysis on the physical contributions to the quadrupole binding energies computed from first-principles as we will discuss later, by using an electrostatic, magnetostatic, and strain (EMS) model. In this model, the binding energy is obtained by separately computing electrostatic (∆ elec quad ), magnetostatic (∆ mag quad ), and strain (∆ strain quad ) contributions. For the electrostatic effect, we compute the electrostatic potential contribution to binding energy (∆ elec quad ) by taking the difference between the quadrupole (quad) and twice the dipole configuration (dipole):
Here k is the Coulomb constant, α is the Madelung constant, and r is the relative permeability (22.9 for hematite). 33 The summation goes over all polarons and dopants i and j, with relative charges q i and q j , and physical separation r ij .
Magnetic effects were computed using the Heisenberg Hamiltonian H spin = -1 2 ij J ij Ŝi • Ŝj , where J ij is the magnetic coupling between the spins of the i th and j th ion, and Ŝi is the spin of the i th ion. Here we use the magnetic exchange coupling constants computed from Ref. 34, which provided magnetic couplings for both the superexchange between two Fe(3+) or between two Fe(2+), as well as the double-exchange between Fe(3+) and Fe(2+). For high-spin Fe(3+) and Fe(2+), the value of Ŝi is 5/2 and 2, respectively, while the spin of the tetravalent dopants is zero (hence the magnetic interaction with these dopants is always zero). In this way, we can compute the magnetic contribution to the binding energy from the magnetic energy of the quadrupole system subtracted by two times of the dipole system:
In this work, we assume J ij between two Fe ions before and after doping are the same, as lattice distortions are generally small compared to Fe distances. With the above consideration and the non-magnetic nature of dopants, our computed magnetic energy is identical for all the dopants.
Finally, in order to compute the strain contribution to the binding energy, we evaluated the change in energy induced by lattice distortions before and after doping. As usual, taking this energy for the quadrupole system and subtracting twice the dipole system:
Here E strain X is the strain energy of system X, computed as a difference of total energy of the pristine system with its equilibrium geometry and with relaxed geometry from the doped system (first relax with dopants then substitute back Fe atoms to keep the same composition as pristine Fe 2 O 3 ). We note a similar approach was used in Ref. 35 to evaluate strain energies.
We computed the formation energy (E f q ) of each defect at a charge state q according to:
where E q (X) is the total energy of the defect system (X) with charge q, E prist is the total energy of the pristine system, µ i and ∆N i are the chemical potential and change in the number of atomic species i, and ε F is the electron chemical potential. A charged defect correction ∆ q was computed with techniques developed in Refs. 36,37 and implemented in the JDFTx code. 38 The chemical potentials were carefully evaluated against the stability of byproduct compounds as detailed in the elemental chemical potential section of the SI.
The corresponding charge transition levels of the defects were obtained from the value of ε F where the stable charge state transitions from q to q .
q|q =
The ionization energies are computed by referencing the CTLs to the free polaron state. [39][40][41][42] Namely in Fe 2 O 3 , it has been experimentally observed that photoexcited carriers relax on picosecond timescale to form small polarons, 11 which have been measured to form at energies ∼0.5 eV below the conduction band minimum (CBM). 10,12 Theoretically this free polaron level is computed as the charge transition level from q = 0 to q = -1 in the pristine system,
. By this method we obtain that the free polaron level is positioned at 0.497 eV below the CBM in excellent agreement with experimental observation. 10,12 From charge defect formation energies, charged defect concentration (c q ) can be computed as:
where g is the degeneracy factor accounting for the internal degrees of freedom of the point defect, k B is the Boltzmann factor, and T is temperature. Concentrations including intrinsic defects, extrinsic dopants and free electron polarons were computed by determine their charge neutrality condition 24,43,44 (additional details are provided in the defect concentration section of the SI). To best relate to experimental measurements of Fe 2 O 3 photoanodes, we computed the concentration first at a synthesis temperature of T S = 1073 K (800 • C is a common synthesis temperature 14,20 ), and then recomputed charge neutrality at normal operation (room) temperature T O = 300 K. The partial pressure of oxygen gas (p O 2 ) corresponding to atmospheric condition (1 atm) is used.
Substitutional doping by group IV and XIV elements was investigated theoretically by replacing a single Fe site by the dopant (X). Consistent with previous studies 25,27 and experimental observation, 14,15 we found this process yields the formation of small electron polarons corresponding to the identification of Fe(2+) after replacing Fe(3+) by the tetravalent dopant For the atomistic plots, gold=Fe, red=O, grey=Sn, and the yellow/blue (+/-) cloud is the isosurface of the polaron wavefunction (the isosurface level is 1% of the maximum). In the band structures, dark/light blue is spin up/down and ε F is the Fermi energy.
X(4+). (Note in this paper we use the notation X(i) to denote an ion X with valency i).
The electronic structure of the single Sn-doped system is shown in Figure 1b right panel, wherein the band structure exhibits a flat isolated occupied state in the gap corresponding to the small electron polaron (EP ) with tight spatial localization similar in size to the Fe-O bond lengths. Likewise, the projected density of states (PDOS) in Figure 1b shows a sharp isolated peak composed mostly by Fe 3d. The wavefunction of the EP is shown in Figure 1a with a clear d z 2 character. The EP forms at the Fe site nearest to the Sn dopant with a Sn-EP distance of 2.981 Å (d 0 in Figure 1a).
To investigate dopant-dopant interactions, we placed a second Sn dopant in the lattice.
All possible Sn-Sn pair configurations were tested, and the lowest energy configuration was clearly identified (see SI Figure S2, S3 and Table S2). The electronic structure of this corresponding configuration is shown in Figure 1c-d, which exhibit the formation of two EP states. We find the asymmetry of the local structure, which is a natural consequence of the corundum crystalline form, causes a noticeable energetic difference of 0.12 eV between EP 1 and EP 2 . Specifically, in Figure 1c, EP 1 has distances to the two adjacent Sn of The remaining group IV (Ti, Zr, Hf) and XIV (Si, Ge, Sn, Pb) dopants were also simulated in both single and pair dopant configurations with negligible differences in their electronic structure and polaron configurations from Sn (all electronic structures are presented in SI Figure S4-S9). Note for the present study, the configuration of two dopants is chosen to be the same for all dopants for the purpose of discussing chemical trends, as predicted by the case of Sn. It is possible that dopants may vary in their exact pair dopant configuration, for example see SI Table S3. This variation does not affect the main implications on car-rier concentration we conclude later. Most importantly, the stable configuration of a single tetravalent dopant (such as Sn) resembles an electric dipole where the Sn and EP represent positive and negative charge centers, respectively. In this way the system with two dopants resembles an electric quadrupole (two positive Sn centers and two negative EP centers).
Therefore, we will denote the single doped system as a first-order multipole (dipole) system and the pair doped systems as a second-order multipole (quadrupole) system. To examine the thermodynamic stability of dopant-polaron quadrupole, we studied their binding energy (∆ quad ) from two separate dipoles:
Here, E f (X) is the formation energy of the system with neutral dopants in a configuration X (e.g. quadrupole or dipole dopant system) following Eq. 4 at charge state q = 0. The quadrupole binding energy was evaluated for all group IV and XIV dopants considered in this study. We observed that the binding energy for all dopants is negative (∼ -0.1 to -0.2 eV), as shown in Table 1 (∆ DF T quad ), indicative of a strong tendency for dopants and polarons to aggregate. We note that we expect dopant clustering occurs during the cooling process from synthesis temperature (over one thousand K here) down to room temperature. At a synthesis condition, dopants will be all ionized and the binding of dopants into quadrupoles will not occur. Since all the binding energies are lower than kT at room temperature, it is expected that the quadrupoles are stable at room temperature. 2 for the 0.1% Sn and 1.0% Sn samples. For the 0.1% Sn sample (Figure 2a), the amplitudes of the further neighbor peaks are quite large and the data can be well fit (solid orange line) to the hematite structure, with a small expansion for the Sn-Fe pairs compared to hematite; roughly 0.1 for closer pairs but only 0.02 Å for Fe neighbors near 3.7 Å. This is the expected behavior around a substitutional dopant site when the dopant valence Sn(4+) is higher that the host valence Fe(3+), and this behavior has been observed in other similar situations. 45 The further neighbor Sn-O peaks are expected to contract very slightly, but because these small peaks overlap the larger Sn-Fe peaks, the pair-distances fluctuate too much. The first O shell, although split in hematite, collapses to a single peak with an average Sn-O distance of 2.05 Å, very close to the averaged first neighbor distance in hematite, 2.03 Å; this is a competition between a larger ionic radius for Sn(4+), and larger electrostatic force between
On the other hand, the EXAFS r-space plot for 1% Sn sample is quite different (Fig- ure 2b). The data up to 3 Å are very similar to that for 0.1% Sn -i.e. the phase of the real part of the Fourier transform, R(r), is the same. However in the range 3-3.8 Å, the phase changes dramatically and a dip develops in the amplitude near 3.3 Å which has the shape of an interference dip. It occurs close to the expected position for the Sn-Fe peak in an EXAFS plot (actual distance ∼3.7 Å: note that there is a calculable phase shift of peaks in r-space plots to lower r). These data can't be fit to a simple distorted hematite model and the shape of R(r) suggests that another peak is present (see SI Figure S11). Consequently, an additional peak, corresponding to one Sn-Sn pair with a distance close to 3.7 Å, was included in the fit.
The number of Sn-Fe pairs was correspondingly reduced from 6 to 5. This yielded the good fit (details in SI) shown in Figure 2b, and is a clear evidence that Sn-Sn pairs have formed.
Remarkably, this Sn-Sn pair distance (∼3.7 Å) matches the theoretically predicted distance of The plot for 0.1% Sn also shows a fit to the hematite structure; good agreement is obtained with a slight contraction of the Sn-O pairs and a slight expansion of the Sn-Fe pairs. At 1% Sn, the EXAFS changes significantly. Although the first two peaks are very similar, the region from 3-3.8 Å is quite different, particularly the shape of the phase (fast oscillating function), and a dip develops near 3.3 Å. These data cannot be fit to the hematite structure. The data suggest that there is another peak present; in the fit shown in part (b), one of the Fe neighbors at ∼3.7 Å is replaced with a Sn atom, forming a Sn-Sn pair. This leads to the excellent fit shown in (b). Fourier transform range, 3.5-13 Å-1 ; fit range in r-space, 1.1-4.2 for both plots. In both figures, the blue and gold bars at the bottom indicate the position of Sn-O and Sn-Fe peaks, respectively.
After confirming the existence of Sn-Sn clustering both theoretically and experimentally, we turn to investigate the mechanisms of their formation. Above the analogy was made between the single Sn-doped system and electric dipoles, so in order to probe this electrostatic interaction, we plotted the total energy of the single Sn-doped system as a function of Snpolaron distance, as shown in Figure 3a. The computed total energies were fitted to a Coulomb potential (-a/r + b) with the fitted values of a = 0.658 eV • AA and b = 0.230 eV, and a coefficient of determination (R 2 ) of 0.85. This validates a clear electrostatic attraction between the two bodies with opposite charges like a dipole. Furthermore, we find that b is close to the value of ionization energy of Sn (0.25 eV), as expected.
However, this Coulombic interaction competes with other factors. For example in Mo doped BiVO 4 , 46 strain causes short-range repulsion between the dopant and polaron, and dominates over the Coulombic attraction. 46,47 We find this was also reflected in the computed quadrupole binding energies (∆ quad ), as shown in Figure 3b, where we plotted them versus the ionic radius (R I ) of each dopant. Specifically, there is a roughly positive correlation between the ionic radius and the quadrupole binding energy. We attribute this to the compensatory size effects of the dopant and the polarons. Namely, the replacement of Fe(3+) with ionic radius of 64.5 pm by Fe(2+) with ionic radius of 78.0 pm, yields an expansion strain at the lattice site. This strain can be reduced by smaller radii dopants (e.g. Ti, Ge, Si) that will increase the magnitude of ∆ quad as shown in Figure 3b, or enhanced by larger radii dopants (e.g. Sn, Hf, Zr, Pb) that will decrease the magnitude of ∆ quad towards zero. We note that the case of Ti, an outlier in Figure 3b, possess stronger correlated 3d orbitals, which in turn exhibit stronger electron localization, may compensate local expansion from small polarons and lower its binding energy (similar to Si which intrinsically is smaller than the rest). Therefore, it is necessary to consider electrostatic and strain effects simultaneously, and also include magnetic effects (binding energies may also be modified by the antiferromagnetism of Fe 2 O 3 ). Hence, we propose a model of the quadrupole binding energy based on electrostatic, magnetostatic, and strain effects (abbreviated to EMS), in order to analyze the importance of each contribution:
Here ∆ elec quad , ∆ mag quad , ∆ strain quad , correspond to electrostatic, magnetostatic, and strain contributions to the quadrupole binding energy, respectively. The exact formulation for each component of the EMS model is detailed in the methods section (Eq. 1-3) and the results are summarized in Table 1. We evaluate how well this model reproduces DFT calculations by plotting them against each other in Figure 3c. The linear fitting shows adequate agreement between the simple EMS model and our exact DFT calculations (with a slope near unity and an R 2 value of 0.864), which justifies its use for interpreting the DFT binding energies.
As shown in Table 1, each of the three components contributes significantly to the overall quadrupole binding energy. The electrostatic interaction (∆ elec quad ) is typically the most dominant factor, and intuitively it is chiefly responsible for the attraction of dopants into the quadrupole configuration. Interestingly, we also found a non-negligible contribution from magnetostatic interactions (∆ mag quad , 42 meV) which further participates in the binding of quadrupoles. This effect is non-trivial but is an indirect consequence of placing non-magnetic dopants next to each other, which in-turn breaks fewer antiferromagnetic interactions and yields a lower energy configuration when dopant-pairs form. In contrast, the effect of strain (∆ strain quad ) typically mitigates the formation of dopant-pairs (increases the system's energy with clustering) due to the accumulation of lattice distortion. However, as aforementioned, this effect can be compensatory in the dopant cases with smaller ionic radii than Fe(3+), which can pack more efficiently next to the polarons with larger ionic radius as specified in Figure 3b.
Lastly, we discuss the effects of dopant clustering, which we will show to be responsible for the low optimal doping concentrations of Fe 2 O 3 observed experimentally. The computed thermodynamic charge transition levels (CTLs) are displayed in Figure 4a for both the dipole (single dopant) and quadrupole (pair dopants) systems. In quadrupole systems, there are two charge transition levels (gold and orange lines) which correspond to the ionization energies of two electron polarons (as shown in Figure 1c). Because quadrupoles yield consistently lower second CTLs, they possess very high second ionization energy ranging from 0.34 eV to 0.48 eV (the orange lines in Figure 4a, nearly doubled the ionization energies of corresponding single-doped systems denoted by the blue lines). Therefore, the ionization of both electrons after quadrupole formation is nearly impossible. On the other hand, the first ionization energies of quadrupoles compared to the ones of dipoles shift in a manner correlated with the ionic radius as shown in Figure 4b. In particular, dopants with a radius larger than the host site (Fe(3+)) have a lower first ionization after quadrupole formation (orange bars in Figure 4b) and vice versa.
To show how dopant clustering affects carrier concentration quantitatively, we compute the polaron concentrations with and without the effect of dopant clustering for three representative cases: Sn, Ge, and Ti in Figure 4c. Equilibrium polaron concentrations are computed following the defect formation energy and charge neutrality approach as detailed in the method section. 24,43,44 We then introduce clustering in the theoretical synthesis process, by considering the probability at which two dopants form into a quadrupole following a thermal Boltzmann distribution (details in the SI section on probability of quadrupole formation). As dopant-dopant pairs form during cooling process from synthesis temperature to room temperature, an intermediate temperature of T = 500 K is chosen as an approximate temperature where pairs start to form. We show the trends among different dopants remain the same for different temperatures (see SI Figure S13). In Figure 4c, we show that polaron concentration is reduced due to clustering (solid lines) compared to the case without clustering(dashed lines), and the magnitude of this reduction is closely related to their binding energy (∆ quad ) in Table 1. Specifically, dopants such as Ge and Ti with large binding energies show significant trapping of polarons due to clustering in the second and third panels of Figure 4c.
In order to avoid the adverse effects of clustering, Ti will be better suited to be doped at very low concentrations into hematite. This explains the experimentally observed very low optimal doping concentration of Ti in hematite (about 0.1%). 17,18 Furthermore, co-doping Ti with another dopant less prone to clustering may offer higher performance in hematite, which explains the success of recent co-doping strategies. [49][50][51] On the other hand, the small quadrupole binding energy of Sn makes it easier to dope Sn into hematite to higher concentrations without experiencing an immediate bottleneck. This explains why the optimal doping of Sn into hematite (3%) is an order of magnitude larger than Ti. 19,20,23 If neglecting the effect of clustering, Ti would be a better dopant than both Sn and Ge (see dashed lines in Figure 4c, which shows that Ti has predicted to contribute the highest carrier concentration without clustering). Therefore, to unlock the potential of Ti, it is necessary to mitigate its strong tendency of clustering. A strategy of co-doping Sn and Ti in hematite may offer the most effective strategy for maximizing performance of hematite photoanodes. Alternatively, co-doping of divalent dopants, such as Mg, has been shown to relieve lattice distortion and can also offer further improvement to PEC performance. 52 Ultimately, Ge is seen as the best dopant in group IV and XIV with the highest polaron concentration, which outperforms Sn and Ti in Figure 4c.
Finally, we remark that here we have focused on formations of dopant clustering at an early stage, which dominate at relatively low concentrations of doping in hematite. We also tested higher-order multipole clustering such as 'hexapole' formation in hematite (e.g. three Sn dopants with the three introduced EP , details in SI section on higher-order multipoles and SI Figure S14). We found that hexapoles also have a negative binding energy, and thus it is entirely possible that dopant clustering may grow even larger than the second-order multipoles we considered here. However, there will also be larger and more complicated strain effects and configurational entropy which can compensate binding energies of larger sized aggregations. Eventually, at even larger doping concentrations, the precipitation of different phases may occur, for example at 6% Sn-doping in hematite XRD shows SnO 2 secondary phases form. 20 Future theoretical work should provide further insights to PEC experiments by investigating the role of dopants and dopant clustering on polaron mobility via small polaron hopping 25,46,47,53 and optical absorption, 41,54,55 which together, along with carrier concentrations studied here, directly impact the photoconductivity of transition metal oxide based photoelectrodes.
This work discovers a new mechanism of dopant clustering via the aggregation of n-type dopants and electron polarons into dopant-pairs which resemble electric multipoles. These pairs are thermodynamically stable due to several contributions which we disentangle by a simple model involving three components: electronic, magnetic, and strain. Our model illuminates that binding occurs predominantly through electrostatic interactions but surprisingly is also mediated by magnetic interactions which together overcome strain to yield the consistently negative binding energies of tetravalent dopants in Fe 2 O 3 . EXAFS experiments confirm the existence of these Sn-Sn pairs which formed at 1% doping and have an identical interatomic distance compared to those predicted theoretically (∼ 3.7 Å).
The effect of doping with and without clustering on carrier (small electron polaron) concentration is carefully examined. We find doping in Fe 2 O 3 is limited by dopant clustering which traps electron polarons and severely lowers the carrier concentration with respect to doping concentration. This clustering is shown to be responsible for the doping bottleneck in hematite, where dopants such as Ti exhibit extremely low optimal doping concentration (i.e. 0.1%) for PEC application. Strategies to overcome this doping bottleneck are proposed;
specifically codoping with dopants that exhibit low binding energies for clustering (e.g. Sn-Ti codoping) is seen as an ultimate solution. Lastly, for single-type doping in group IV and XIV, we found Ge is the best dopant which can contribute the highest polaron concentrations even