Electron small polarons and their transport in bismuth vanadate: a first principles study

Abstract

Relatively low electron mobility has been thought to be a key factor that limits the overall photocatalytic performance of BiVO₄, but the behavior of electrons has not been fully elucidated. We examine electron localization and transport in BiVO₄ using hybrid density functional theory calculations. An excess electron is found to remain largely localized on one V atom. The predicted hopping barrier for the small polaron is 0.35 eV (with inclusion of 15% Hartree–Fock exchange), and tends to increase almost linearly with lattice constant associated with pressure and/or temperature changes. We also examine the interaction between polarons, and discuss the possible concentration-dependence of electron mobility in BiVO₄.

---

1. Introduction

Bismuth vanadate (BiVO₄) is one of the most promising photoanode candidates for water splitting.^1–11 Its band gap E_g is sufficiently small (≈2.4 eV) to allow the absorption of visible light as compared to only ultraviolet light absorption by titanium oxide (TiO₂, E_g ≈ 3.0–3.2 eV), which is the most widely used. However, in practice, the solar-to-hydrogen (STH) conversion efficiency of BiVO₄ is far lower than the maximum theoretical value.¹² Previous experimental and theoretical studies have identified that the overall photoelectrochemical performance of BiVO₄ is largely limited by poor electron transport,^8,9 which facilitates the recombination of photogenerated electron–hole pairs. To overcome the drawback, many efforts have been made to find a robust and effective way to enhance the photoinduced charge separation, including doping^4,7–9 and heterostructure formation.^10,11 Although electron mobility appears to be a major bottleneck that we have to solve for realizing the desired efficiency of BiVO₄-based photoanodes, so far the nature of an excess electron and its transport properties in BiVO₄ have not been well studied.

In oxide semiconductors, excess charge carriers (electrons and holes) often remain localized.^13,14 The charge localization is commonly accompanied by lattice distortions, which results in a so-called polaron.^13–15 Density functional theory (DFT) is perhaps the most successful computational method for the study of semiconductor and oxide materials; however, in practice, the results rely on approximations for the exchange–correlation potential which accounts for the many-body electron–electron interaction. In particular, the microscopic description of localized charge states tends to depend strongly on the exchange–correlation (xc) energy functional employed. The (semi)local functionals of conventional DFT may fail to predict charge localization due to their inherent self-interaction error (which erroneously favors charge delocalization).^16,17 The self-interaction effects can be effectively canceled by the partial inclusion of non-local Hartree–Fock (HF) exchange;^16–18 the HF/DFT hybrid approach has been shown to be capable of predicting charge localization and transport in many metal oxide semiconductors.^19–21

In this work, we investigate the behavior of excess electrons in BiVO₄ using hybrid DFT calculations. We first identify the formation and structure of a small electron polaron, and evaluate its relative stability with respect to the completely delocalized (nonpolaronic) state by varying the fraction α of HF exchange. The α dependence of the polaron formation energy is also analyzed by decoupling it in terms of electronic energy gain and strain energy loss. Next, we show a variation in the predicted barrier for adiabatic electron hopping with α and discuss the underlying reasons. Finally, the possible effects of electron concentration and hydrostatic pressure on the mobility of small electron polarons are discussed.

2. Computational methods

Spin-polarized hybrid DFT calculations were performed using the Vienna Ab initio Simulation Package (VASP 5.2.2).²² The projector augmented wave (PAW) method with a planewave basis set (E_cut = 450 eV) was used to describe the interaction between the core and valence electrons.²³ We employed a PBE (Perdew–Burke–Ernzerhof)²⁴-based hybrid functional in which a fraction α of the PBE exchange is replaced with HF exchange. Here, the slowly decaying long-range part of the HF interaction was excluded following the Heyd–Scuseria–Ernzerhof (HSE) screened approach.²⁵ As such, the xc energy is given by

E^HSE_xc = αE^HF,SR_x(μ) + (1 − α)E^PBE,SR_x(μ) + E^PBE,LR_x(μ) + E^PBE_c,

where the range-separation parameter μ determines the partition of the short-range (SR) and long-range (LR) components. In our calculations, μ = 0.207 Å⁻¹ was used; this has already been demonstrated to be a reasonable compromise between accuracy and computational cost.²⁶

Bulk BiVO₄ was modeled using a supercell consisting of 36 BiVO₄ formula units with periodic boundary conditions in all three directions; for a charged system, a homogeneous countercharge was added to maintain overall charge neutrality. The atoms in the supercell were fully relaxed until the residual forces on all the constituent atoms become smaller than 0.02 eV Å⁻¹. The Brillouin zone was sampled using a single k-point (Γ). For electron-doped BiVO₄, we considered both localized (polaronic) and delocalized (nonpolaronic) electron states; to create a polaronic state, we initially applied a small perturbation around a selected VO₄ tetrahedron to break the lattice symmetry prior to structural relaxation, while a nonpolaronic state was obtained by relaxing the structure without such initial perturbation.

3. Results and discussion

Fig. 1(a) shows a schematic illustration of the tetragonal (scheelite) structure of BiVO₄ which consists of corner-sharing BiO₈ dodecahedra and VO₄ tetrahedra. As shown in the orbital-resolved density of states (DOS) in Fig. 1(a), the low-lying empty conduction bands are mainly composed of V 3d states with minor contributions from the Bi 6p and O 2p states. Note also the d orbital splitting into t₂ and e states by the tetrahedral crystal field; additionally, the compression in the z-direction causes a splitting in the e states^27,28 while leaving the d_z² levels at the bottom of the conduction band (CB).

Fig. 1 (left panels) Orbital projected density of states (DOS) for (a) pristine and (b) electron-doped tetragonal BiVO₄; the grey, green, blue, and red lines represent the Bi 6s, Bi 6p, V 3d, and O 2p states, respectively, and the energy zero is set at the Fermi level (E_F) which is indicated by the vertical dashed line. In addition, the V d_z² states (shaded in blue) are shown in the inset of (a), and the electron polaron state is indicated by an arrow in (b). (upper-right panel) Crystal structure representation of scheelite BiVO₄ with indication of BiO₈ dodecahedra (in purple) and VO₄ tetrahedra (in gray). (lower-right panel) Band-decomposed charge density of the localized electron state (with an isosurface value of 0.02 e Å⁻³); purple, silver, and red balls represent Bi, V, and O atoms, respectively, and the arrows indicate elongation of the V⁴⁺–O bonds.

As shown in Fig. 1(b), adding an electron to the BiVO₄ matrix creates a deep localized state within the band gap, which is evidently due to local trapping of the excess electron; most of the excess charge tends to reside on a V atom by preferentially filling up the empty d_z² states, and only a small fraction is on the surrounding four O atoms. The electron self-trapping is accompanied by significant local lattice distortions around the reduced V⁴⁺ site, forming a small polaron. It is worth noting that noticeable elongation of V⁴⁺–O bonds in the electron-doped case (with respect to V⁵⁺–O bonds in undoped BiVO₄) can be largely attributed to the filling of antibonding V 3d–O 2p states.

We also evaluated the relative stability of the polaronic state with respect to the completely delocalized (nonpolaronic) state (in which the excess electron spreads over entire V atoms and the Fermi level is shifted above the bottom of the CB, not shown here). The nonpolaronic state with no local lattice distortion has the nearly same configuration as the undoped structure, as expected.²⁹ As summarized in Fig. 2(a), the total energy difference between the polaronic and nonpolaronic states, referred to as the polaron formation energy (E_p),³⁰ varies with the fraction α of Hartree–Fock exchange added to the PBE-GGA functional; a positive E_p value indicates that the polaronic state is energetically more favorable than the non-polaronic state. The pure PBE functional (α = 0%) tends to erroneously favor electron delocalization, which is due to the incomplete cancellation of the Coulomb self-interaction.^16,17 The E_p is predicted to be 40 meV at α = 5%, and increases monotonically with increasing α; note that adding the HF exchange partially cancels the unphysical self-interaction in the Hartree approximation, which in turn makes the localized state energetically favorable.¹⁶

Fig. 2 Results from hybrid DFT as a function of the fraction α of Hartree–Fock (HF) exchange. (a) Polaron formation (E_p), electronic (E_e), and lattice strain (E_s) energies. The inset shows a graphical interpretation of E_p, E_e, and E_s, and D and L represent the delocalized and localized states, respectively. (b) Band gap (E_g) and localized polaron state location from the valence band (ΔE), as illustrated in the upper left corner. The relative polaron state location with respect to the band gap (ΔE/E_g) is also shown in the inset graph.

We would better understand the α dependence of E_p by decoupling it in terms of the electronic energy gain (E_e) due to lattice polarization and the strain energy loss (E_s) caused by local polaronic lattice distortions,^30,31 as illustrated in the inset of Fig. 2(a). For each case, E_s is obtained from the total energy difference between the undistorted lattice and the distorted lattice due to polaron formation but without the excess electron,³¹ and E_e is estimated by adding E_s and E_p. As shown in Fig. 2(a), while both E_e and E_s increase with α, the former increases more rapidly than the latter; apparently, the large increase in E_e is primarily attributed to the increased Coulomb interaction energy as a result of the increase of charge localization with inclusion of HF exchange.^16,17 A similar trend can be seen from DFT + U calculations;^30,32 a larger value of Hubbard U (analogous to α in hybrid DFT) leads to more localization of charge and in turn a greater E_p. Our calculations also show that the band gap (E_g) and the position of the localized polaron state (ΔE/E_g) are a strong function of α [Fig. 2(b)]. The E_g is predicted to be 2.5 eV at α = 5%, and almost linearly increases with increasing α, which is overestimated in comparison to the experimental value¹ of 2.34 eV while the pure PBE functional somewhat underestimates the band gap.³³ The relative position of the polaron state within the gap is found to be shifted closer to the valance band as α is increased, implying enhancement in the stability of the polaron which is also well demonstrated by the increase in E_p.

A small polaron may undergo thermally-activated hopping from site to site.¹³ The activation energy (E_a) for polaron hopping can be estimated by calculating the variation of potential energy as the polaronic lattice configuration is displaced from the initial to the final geometry; if the initial and final states are identical, the transition state may occur midway between the two local minima. In Fig. 3, we present the predicted E_a for electron hopping along the 〈101〉 (or 〈011〉) direction as a function of the HF exchange fraction α; another probable path in the 〈100〉 (or 〈010〉) direction was considered, but it turns out to be energetically less favorable (not reported here). E_a is found to monotonically increase from 0.10 eV to 0.46 eV as α increases from 5% to 25%. In addition, E_a decreases upon decreasing the supercell size, and the magnitude of the decrease becomes larger with decreasing α. These results can be explained based on the theory for small polaron transport, as described in the following.

Fig. 3 Predicted activation energies (E_a) for electron polaron hopping as a function of the fraction α of Hartree–Fock (HF) exchange for two different supercell sizes as depicted. The schematic diagram in the upper left corner illustrates the Marcus–Emin–Holstein–Austin–Mott theory for small polaron transfer between two adjacent sites; as detailed in the text, the E_a is determined by the reorganization energy (λ) which represents the energy needed to distort the initial-state structure to the final-state structure without electron transfer occurring and the electronic coupling between the initial and final states (H_AB). In this work, the initial (A) and final (B) states are V⁴⁺–V⁵⁺ and V⁵⁺–V⁴⁺, respectively. In the supercell diagrams, the band-decomposed charge densities of corresponding polaron states (with an isosurface value of 0.005 e Å⁻³) are also shown to demonstrate the supercell-size dependence of electron localization; purple, silver, and red balls represent Bi, V, and O atoms, respectively.

As illustrated in the inset of Fig. 3, for adiabatic electron transfer that follows the thermal hopping mechanism, the activation barrier can be estimated in terms of the reorganization energy (λ) and the electronic coupling between the initial and final states (H_AB),³⁴i.e., E_a = λ/4 − [|H_AB| + λ/2 − (λ²/4 + |H_AB|²)^1/2]. According to the Marcus–Emin–Holstein–Austin–Mott theories,^35,36 the initial and final states correspond to minima on their respective potential energy surfaces A and B which are assumed to be parabolic in shape. The steepness of the potential well is characterized by λ (which represents the energy needed to distort the initial-state structure to the final-state structure without electron transfer occurring). The height of the crossing point between the two potential energy surfaces is the diabatic activation energy, and the hopping barrier is reduced by the electronic interaction H_AB between the initial and final states while the localized electron undergoes migration in an adiabatic way.

The reorganization energy λ is directly related to the polaron formation energy E_p which tends to almost linearly increase with α.^34,36 Moreover, the more localized electron at a larger α may result in a smaller electronic coupling H_AB. Hence, it is expected that E_a will increase sublinearly as α is increased. In addition, a smaller supercell (which gives a shorter distance between the polaron and its periodic images) may result in an increase in H_AB, thereby leading to a further reduction in E_a. These hypotheses are well demonstrated and supported by our calculation results in Fig. 3.

In addition, for α < 20%, the substantial change in E_a with supercell size may suggest a possible dependence of electron mobility on its concentration in highly electron-doped BiVO₄, due to the overlap between the diffusing electron wavefunctions. If the polaron–polaron separation is sufficiently small, the polaron wavefunctions can overlap to become stabilized.³² As shown in the inset of Fig. 3, indeed, the excess electron appears to spread to a noticeable extent over neighboring V atoms particularly along the 〈101〉 and/or 〈100〉 directions, as the supercell size is decreased; for the small and large supercells, the distances between the polaron and its periodic images in the 〈100〉 direction are 10.3 Å and 15.4 Å, respectively. Clearly, the increased wavefunction overlap is largely responsible for the reduction of E_a from 0.35 eV (large supercell) to 0.29 eV (small supercell) at α = 15%; the concentration effect can be overestimated or underestimated, depending on the choice of α, such that its careful selection is necessary. A similar reduction in E_a due to the overlap of polaron orbitals has been also suggested for hole hopping in Li₂O₂ and Li₂CO₃.³² The orbital overlapping will diminish as the polaron concentration decreases, and appears to be unimportant when the polarons are separated by more than several atomic distances if they are strongly localized. Nonetheless, our theoretical study, warranting further experimental validations, highlights the possibility that the E_a of electron hopping could be concentration dependent in highly electron-doped BiVO₄.

Although an optimal α that gives an accurate description of electron localization and transport would hardly be determined, we could choose a reasonable value through comparison to experimental observations. In view of this, the good agreement between our predicted E_a (≈0.35 eV) at α = 15% for the low electron concentration case (large supercell) and the experimentally estimated value³⁷ of 0.29 eV in 0.3% W-doped BiVO₄ may suggest that α = 15% would offer a reasonable description of the behavior of small electron polarons in BiVO₄. Perhaps, we could also validate the selection of α via comparison to the spectroscopic characterization of electronic states. As shown in the inset of Fig. 2(a), the deep polaron level below the conduction band can be also detected using various optical measurements such as infrared and core-level X-ray spectroscopy and electron spin resonance spectroscopy;^31,38 therefore, we would like to call for experimentalists to investigate the formation and nature of localized polaronic states in electron-doped BiVO₄.

Finally, we examined the effect of hydrostatic pressure on the barrier E_a for electron hopping by varying the supercell volume (or lattice constant). Fig. 4 presents predicted E_a values at α = 15% as a function of the supercell volume normalized with respect to the equilibrium volume of V₀ = 2768 ų (=15.40 × 15.40 × 11.67 ų). E_a increases with volume (or lattice constant); that is, E_a increases (decreases) under tension (compression). The E_a variation of a few tens of meV is not trivial considering that electron mobility is exponentially dependent on E_a (μ = μ₀ exp(−E_a/k_BT), where μ₀ is the mobility prefactor¹³). The inset shows the calculated value of E_p, following the same trend as E_a; this indicates that polaron formation is facile under tensile conditions, which can be well understood given that the energy cost for the lattice distortions associated with a polaron is relatively less in comparison to the compressive case. The direct correlation between E_p and E_a may suggest that the increase of E_a with pressure is mainly attributed to an increase in λ, with a minor contribution from H_AB. We think that the extent of H_AB would be also affected to a certain degree by the changes in lattice constant, as the degree of electron localization is likely to vary with volume (see ρ_e variation in Fig. 4).

Fig. 4 Predicted activation energies (E_a) for electron polaron hopping as a function of normalized supercell with respect to its equilibrium value (V₀), i.e., ΔV = 100 × (V − V₀)/V₀; in the hybrid DFT calculations, 15% of HF exchange was included. The variations of the polaron formation energy (E_p) and the degree of electron localization around the VO₄ tetrahedron (ρ_e) are also shown in the inset.

The pressure dependence of E_a may raise a question regarding the temperature effect, as there is a direct relationship between temperature and volume, i.e., the volume (or lattice constant) increases with temperature. The measured volumetric thermal expansion coefficients for BiVO₄ are around 8 × 10⁻⁶ °C⁻¹ in the monoclinic phase and 50 × 10⁻⁶ °C⁻¹ in the tetragonal phase.³⁹ Taking those values, the volume expansion of BiVO₄ is estimated to be 0.2% (1.5%) as the temperature is increased by 100 °C in the monoclinic (tetragonal) phase, corresponding to an increase in E_a by about 2 (20) meV. This suggests that the temperature effect could be important if the sample temperature varies significantly.

4. Conclusions

We examined the behavior of excess electrons in BiVO₄ using hybrid functional calculations with the varying fractions α of HF exchange. Our calculation clearly shows that an excess electron remains largely localized on one V atom, accompanied by significant local lattice distortions around the reduced V⁴⁺ site. At α = 15%, the small polaron state is predicted to be about 0.5 eV more favorable than the completely delocalized state; the polaron formation energy E_p increases monotonically with increasing α, as a result of the larger increase in the electronic energy gain E_e relative to the strain energy loss E_s. The polaron hopping barrier E_a is found to increase monotonically from 0.10 eV to 0.46 eV as α increases from 5% to 25%; the predicted values are consistent with the recently reported experimental estimates (≈0.3 eV in 0.3% W-doped BiVO₄). Our study also demonstrates the possibility of overlap between the diffusing electron wavefunctions, leading to the concentration dependence of electron mobility in highly electron-doped BiVO₄ (where the polaron–polaron separation is sufficiently small). Moreover, the hopping barrier E_a tends to increase with lattice constant; however, under typical operating conditions for pressure and temperature, the variation of E_a appears to be small, but not negligible. The improved understanding could assist in identifying effective ways to enhance the electron mobility and thus photocatalytic performance of BiVO₄.

Acknowledgements

This work was partially supported by the Welch Foundation (F-1535), the National Science Foundation (DMR-1122603), and the Samsung Advanced Institute of Technology (SAIT). We also acknowledge the Texas Advanced Computing Center for the use of their computing resources.

References

1. S. Tokunaga, H. Kato and A. Kudo, Chem. Mater., 2001, 13, 4624.
2. J. Yu and A. Kudo, Adv. Funct. Mater., 2006, 16, 2163.
3. H. Luo, A. H. Mueller, T. M. McCleskey, A. K. Burrell, E. Bauer and Q. X. Jia, J. Phys. Chem. C, 2008, 112, 6099.
4. H. S. Park, K. E. Kweon, H. Ye, E. Paek, G. S. Hwang and A. J. Bard, J. Phys. Chem. C, 2011, 115, 17870.
5. Y. Park, K. J. McDonald and K. Choi, Chem. Soc. Rev., 2013, 42, 2321.
6. H. S. Park, K. C. Leonard and A. J. Bard, J. Phys. Chem. C, 2013, 117, 12093.
7. F. F. Abdi, L. Han, A. H. M. Smets, M. Zeman, B. Dam and R. van de Krol, Nat. Commun., 2013, 4, 1.
8. F. F. Abdi, N. Firet and R. van de Krol, ChemCatChem, 2013, 5, 490.
9. J. A. Seabold, K. Zhu and N. R. Neale, Phys. Chem. Chem. Phys., 2014, 16, 1121.
10. S. J. Hong, S. Lee, J. S. Jang and J. S. Lee, Energy Environ. Sci., 2011, 4, 1781.
11. J. Su, L. Guo, N. Bao and C. A. Grimes, Nano Lett., 2011, 11, 1928.
12. Z. Chen, T. F. Jaramillo, T. G. Deutsch, A. Kleiman-Shwarsctein, A. J. Forman, N. Gaillard, R. Garland, K. Takanabe, C. Heske, M. Sunkara, E. W. McFarland, K. Domen, E. L. Miller, J. A. Turner and H. N. Dinh, J. Mater. Res., 2010, 25, 3.
13. A. J. Bosman and H. J. van Daal, Adv. Phys., 1970, 19, 1.
14. A. L. Shluger and A. M. Stoneham, J. Phys.: Condens. Matter, 1993, 5, 3049.
15. D. Emin, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 48, 13691.
16. A. J. Cohen, P. Mori-Sánchez and W. Yang, Science, 2008, 321, 792.
17. P. Mori-Sánchez, A. J. Cohen and W. Yang, Phys. Rev. Lett., 2008, 100, 146401.
18. M. Ernzerhof and G. E. Scuseria, J. Chem. Phys., 1999, 110, 5029.
19. D. M. Ramo, A. L. Shluger, J. L. Gavartin and G. Bersuker, Phys. Rev. Lett., 2007, 99, 155504.
20. J. B. Varley, A. Janotti, C. Franchini and C. G. Van de Walle, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 081109.
21. T. Maxisch, F. Zhou and G. Ceder, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 104301.
22. G. Kresse and J. Furthmuller, VASP: The Guide, Vienna University of Technology, Vienna, Austria, 2001.
23. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
24. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
25. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207.
26. L. Schimka, J. Harl and G. Kresse, J. Chem. Phys., 2011, 134, 024116.
27. M. Itoh and T. Aoki, J. Phys.: Condens. Matter, 2010, 22, 045503.
28. V. V. Laguta, A. Vedda, D. Di Martino, M. Martini, M. Niki, E. Mihóková, J. Rosa and Y. Usuki, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 235108.
29. K. E. Kweon and G. S. Hwang, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 205202.
30. M. Setvin, C. Franchini, X. Hao, M. Schmid, A. Janotti, M. Kaltak, C. G. Van de Walle, G. Kresse, U. Diebold, 2014, arXiv:1401.7817 [cond-mat.str-el].
31. A. Janotti, C. Franchini, J. B. Varley, G. Kresse and C. G. Van de Walle, Phys. Status Solidi RRL, 2013, 7, 199.
32. J. M. Garcia-Lastra, J. S. G. Myrdal, R. Christensen, K. S. Thygesen and T. Vegge, J. Phys. Chem. C, 2013, 117, 5568.
33. K. E. Kweon and G. S. Hwang, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 165209.
34. J. Blumberger and K. P. Mckenna, Phys. Chem. Chem. Phys., 2013, 15, 2184.
35. D. Emin and T. Holstein, Ann. Phys., 1969, 53, 439; I. G. Austin and N. F. Mott, Adv. Phys., 2001, 50, 757.
36. N. A. Deskins and M. Dupuis, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 75, 195212.
37. A. J. E. Rettie, H. C. Lee, L. G. Marshall, J. Lin, C. Capan, J. Lindermuth, J. S. McCloy, J. Zhou, A. J. Bard and C. B. Mullins, J. Am. Chem. Soc., 2013, 135, 11389.
38. D. C. Cronemeyer, Phys. Rev., 1959, 113, 1222.
39. R. M. Hazen and J. W. E. Mariathasan, Science, 1982, 216, 991.

---