Charge carrier transport dynamics in W/Mo-doped BiVO4: first principles-based mesoscale characterization

Viswanath Pasumarthi a, Taifeng Liu *b, Michel Dupuis *acd and Can Li de
aDepartment of Chemical and Biological Engineering, University at Buffalo, Buffalo NY, USA. E-mail:
bNational & Local Joint Engineering Research Center for Applied Technology of Hybrid Nanomaterials, Collaborative Innovation Center of Nano Functional Materials and Applications, Henan University, Kaifeng, 475004, China. E-mail:
cComputational and Data-Enabled Science and Engineering Program, University at Buffalo, Buffalo, NY, USA
dDalian National Laboratory for Clean Energy, Dalian Institute of Chemical Physics, Dalian, China
eState Key Laboratory of Catalysis, Dalian Institute of Chemical Physics, Dalian, China

Received 14th October 2018 , Accepted 21st November 2018

First published on 23rd November 2018

We present mesoscale characterization of carrier transport in W/Mo-doped ms-BiVO4 (BVO) to supplement our earlier study of electron and hole transport in bulk BVO. The mesoscale kinetic Monte Carlo (KMC) approach, supported by first principles-determined electron and hole hopping rates, captures the complex dynamics of electron carriers arising from light absorption and from metal doping. The computations and simulations support the observation that, for the doping level used in experiment, doping atoms do not affect significantly the transport dynamics of the charge carriers compared to stoichiometric BVO. W/Mo doping increases the electron carrier concentration and consequently the electrode conductivity. We used a density functional theory DFT + U method to characterize the electronic structure of doped BVO. Each W and Mo atom brings one more valence electron than a V atom. These excess electrons are mobile. We adopted the theories and methods of our earlier investigations on BVO. The DFT + U theory affords an accurate description of small electron polarons with electrons localized on W or Mo or V, as well as of the hopping barriers from W/Mo-to-V or V-to-V. Calculations on a 3 × 3 × 1 supercell of W/Mo-doped ms-BVO indicate that the excess electron from W or Mo does not reside on the doping atom. There is a shallow interaction region around the doping atom where the excess electron localized on a V atom is more stable than when localized on the doping atom and slightly more stable than when the excess electron is localized far away from the doping atom. This region is two nearest V atom-wide for W and three nearest atom-wide for Mo. The depths of these stabilization regions are very small, ∼−0.65 kBT (at 300 K) for W and ∼−0.73 kBT for Mo. Calculated V-to-V hopping barriers in the doped material are little affected within the region and not affected outside the region. Using our recently developed lattice-based kinetic Monte Carlo code adapted to account for doping atoms and their energy stabilization regions, we calculated the diffusivity of electrons for carrier density relevant to experiment as well as the conductivity. The stabilization regions have little effect on the diffusivity compared to the stoichiometric system because of the smallness of the stabilization. The diffusivity is found to decrease slightly with an increasing number of carriers, but the conductivity of a system with electron polarons arising from doping together with light absorption increases compared to that of the un-doped system. This work will set the foundation to study electron transport in gradient (W/Mo)-doped systems and other mixed phase systems.

I. Introduction

The holy-grail in efficient and cost-effective conversion of solar energy into electrical and chemical energy is solar energy-driven water splitting using semi-conductor-based photo-catalysts. Photoelectrochemical (PEC) cells are attractive due to their design whereby the direct contact between the semiconductor and the liquid electrolyte integrates efficiently the processes of light absorption and conversion of generated charge carriers. Such a design is much simpler and potentially more efficient than ‘integrated photovoltaic-electrolysis cells’.1,2 In any system intended for harvesting solar energy, including PEC cells, the semiconductor material plays the most crucial role. Semiconductor electrodes absorb sunlight and create excited charge carriers (electrons and holes). Upon separation, these charge carriers generate an electrochemical potential that, if suitable, can drive the reaction(s) at the surface of the electrode. Thus, the choice of the semiconductor material is critical in governing the performance of PEC cells. The characteristics of a viable photo-electrode are as follows: (i) it must generate a sufficient voltage to split water; (ii) it must have a suitable bandgap to harvest a significant portion of the solar spectrum while developing a sufficient photo-potential to drive the water splitting reactions; (iii) it must have well-positioned band edges to support the water splitting half-reactions; (iv) it must exhibit long-term stability under aqueous and oxidizing conditions; (v) it must require minimal kinetic over-potentials for charge transfer processes at the semiconductor–liquid junction (SCLJ); and lastly (vi) the semiconductor material must be made of abundant and cheap chemical elements.1–3

In the pursuit of ideal materials, semiconducting transition metal oxides are considered as potential candidates primarily owing to their chemical stability in aqueous environments. However, the poor semiconducting properties of binary oxides have resulted so far in small photocurrents compared to theoretical maxima.1 The limited efficiency in extracting the photo-generated charge carriers prompted a shift towards ternary and more complex multinary oxides. Among the large number of possible combinations for mixed metal oxides,4 bismuth vanadate BiVO4 (BVO) is one of the highest performing ternary metal oxide photo-anode materials reported to date.1 The most active polymorph of BVO has the monoclinic-scheelite crystal structure (ms-BVO) with impressive theoretical efficiency limits of photocurrent densities.5 The excellent activity of ms-BVO is a result of its well suited band structure for water splitting with an indirect band gap of ∼2.4–2.5 eV (ref. 6) and a favorable conduction band edge positioned close to the thermodynamic hydrogen evolution reaction (HER) potential.7 However, the major limitations that hinder this material from achieving theoretical efficiencies are rather poor electron–hole separation,7 poor electron transport,8 and poor water oxidation kinetics.9 In order to address these challenges, a wide variety of strategies have been implemented, which include morphology control, doping, and, finally, coupling with oxygen evolution catalysts.9

Of the above mentioned strategies, doping BVO with either Mo6+ (ref. 10 and 11) or W6+ (ref. 12 and 13) or simultaneous doping with Mo6+ and W6+ (ref. 8) at the V5+ sites14 was found to be effective in increasing the photo-electro-catalytic activity in terms of electronic conductivity. A general observation resulting from these doping studies is an increased carrier density as interpreted from Mott–Schottky plots.15 A consequence of the influx of excess mobile electrons from dopant sites is an increased carrier density n and thus an increase in electrical conductivity σ according to

σ = neμ(1)
where e is the elementary charge and μ is the carrier mobility.

Carrier transport in BVO and (W/Mo)-doped BVO has been the subject of many recent experimental and computational investigations, including from our group.16–23 We recently reported first principles-based mesoscale simulations of electron and hole transport in stoichiometric bulk BVO,16 whereby hopping activation energies were obtained from DFT + U calculations. This earlier work is illustrated in Fig. 1. We determined all possible hopping energies for electrons and for holes. We reported selected cluster calculations of the electronic coupling (VAB) underlying the small polaron theory that led us to identify the strong non-adiabatic character of the electron transport processes. We assigned values of ΔG* ∼ 0.37 eV for electron hops corresponding to a reorganization energy (in Marcus parlance) of λ ∼ 1.48 eV. Hole transport was shown to be bimodal, with fast but not transport-efficient hops among pairs of VO4 tetrahedra and slow but transport efficient hops across VO4 pairs. Nonetheless, the electron mobility is about one order of magnitude smaller than the hole mobility, so that the conductivity in BVO is gated by electron transport. We determined that electron hops have VAB electronic coupling (in Marcus/Holstein parlance) values of ∼600 cm−1, as confirmed by constrained DFT calculations,22 in between the regimes of diabatic and adiabatic hopping. Our value of ΔG* for electron polarons is in good accord with the experimental drift mobility value of ΔG* ∼ 0.300 eV.24 Lattice-based KMC showed that both electron and hole mobilities decrease slightly with increasing carrier density.

image file: c8ta09899a-f1.tif
Fig. 1 Crystalline structure of ms-BVO comprising 4 chemical structures per unit cell. MO4 (M = V/W/Mo) units are shown in polyhedra. The figure illustrates the migration of the excess electron from the dopant site (X = W/Mo) to a neighboring V site. A typical energy profile of electron hopping (through-space) between two V sites is illustrated in the inset. The intersection of potential energy curves was determined along the pathway interpolated between the initial and final participating polaron structures with their localized wavefunction states.

Abdi et al.20 and Han et al.21 showed that gradient doping of BVO with W atoms into n+–n homojunctions was a successful strategy to increase the carrier separation and BVO's overall conversion efficiency. Jovic et al.18,19 used soft X-ray absorption spectroscopy to establish the key characteristics of W/Mo-doping in ms-BVO: the excess electron from W/Mo does not reside on the doping atom, but rather as a small polaron localized on V sites. They also observed a small anisotropy in the electronic structure between the ab-plane and the c-axis. Our earlier study16 also yielded such a finding at the level of the dynamics. Overall, these authors concluded that “the improved conversion efficiency of W-doped BVO relative to the un-doped material is largely due to the increased carrier concentration in spite of increased carrier recombination rates”.

Park et al.8 applied the DFT + U approach to understand the effects of W or Mo doping on the electronic structure of BVO. They found that the excess electron from W or Mo localizes on V atoms and creates a localized state in the band gap. The small polaron formation involves a sizable lattice distortion around the reduced V4+ ion. They suggested that the excess electron residing at the V4+ site may undergo migration to an adjacent V5+ (or Mo6+ or W6+) site. Yin et al.25 studied the extrinsic doping properties of BVO and they found that substitutional Mo and W atoms on V sites are very shallow donors and have very low formation energies. Ding et al.26 explored the origin of the enhanced photocatalytic activity of Mo-doped BVO using DFT. They showed that the Mo atoms prefer to substitute the V atoms in the bulk phase, but for the (010) surface, Mo atoms prefer to substitute Bi atoms at the outermost layer. Overall it is accepted that doping with Mo or W in BVO brings excess electrons and these electrons form mobile small polarons.

Wu et al.22 and Zhang et al.23 provided further experimental evidence and computational understanding why W-doping and Mo-doping are the source of higher photocurrent. Several DFT-based levels of theory showed that the excess electron from W and Mo does not like the doping site but becomes a small polaron in the material. They calculated a hopping activation energy of ΔG* ∼ 0.25 eV and suggested that hopping barriers are reduced by the presence of doping atoms due to some strong electrostatic interactions. This value is somewhat smaller than our previously reported value of ∼0.37 eV (ref. 16) and also smaller than the experimental value of ∼0.30 eV.

In this paper, we report on the nature of the stability and mobility of electron polarons in ms-BVO when doped with W or Mo as substitutions at V sites. QM calculations confirm that the excess electron localization on W/Mo sites is energetically unfavorable. At the same time, localization of the excess electron becomes energetically favorable by a small amount on V sites extending up to ∼the 3rd nearest neighbor shell in the vicinity of the dopant site. It is widely accepted that slow carrier transport in complex metal oxides arises from the nature of carrier localization in the form of small polarons,3,27 so that the characterization of the electron transport in doped systems ought to account for the stabilization (trapping) near the dopant sites. Using data from first-principles calculations, we performed KMC simulations to characterize the nature of electron transport in the presence of extrinsic doping in ms-BVO. The paper is organized as follows: in Section II we describe the computational methods and details, first about the DFT calculations and then about the KMC implementation, in particular with regard to the treatment of substitutional defects. In Section III we present and discuss the results. We conclude in Section IV.

II. Computational methods

a. DFT calculations

We performed spin-polarized DFT calculations using the Vienna ab initio Simulation Package (VASP).28,29 The exchange correlation potential was described using the Perdew–Burke–Ernzerhof (PBE) functional within the generalized gradient approximation.30 The projector augmented wave method was applied to describe electron–ion interactions.31,32 We used the +U correction approach for the self-interaction correction,33 using the values of +Ueff = (UJ) = 5.0 eV (J = 1.0 eV) for V(3d), W(5d), and Mo(4d) states when performing charge localization calculations. The +Ueff values for the d orbitals of metal ions are taken from earlier work,16 affording localization of electron polarons on desired sites. The +U values for Mo(4d) and W(5d) were taken as the same as that for V(3d). We used a 3 × 3 × 1 supercell with 216 atoms. A plane-wave cutoff energy of 400 eV and a (1 × 1 × 2) k-point mesh were used for geometry optimization. The Γ-centered k-point mesh of the Brillouin zone sampling was set at (2 × 2 × 1) based on the Monkhorst–Pack scheme.34 The structures were relaxed until all forces on ions were less than 0.01 eV Å−1. The core electrons were represented by pseudo-potentials with 6, 11, 12, 12, and 5 valence electrons for plane-wave description of O, V, W, Mo and Bi, respectively. We used the conjugate gradient algorithm for ionic convergence with the energy difference cutoff set to 1.0 × 10−5 eV. For Mo or W single doping, one V atom was substituted by one Mo or W atom in the supercell, corresponding to a concentration of ∼2.8%. As for Mo and W co-doping, two V atoms far away from each other are substituted by Mo and W atoms in the supercell.

b. Polaron modeling

To describe polaron hopping, we made use of the Marcus/Holstein two-state model, of which several descriptions exist.35–40 The details of how to use DFT to obtain the important parameters of Marcus theory are published in several papers,36–38,41 including the definition of approximate reaction pathways. These papers also describe how to use the DFT data to obtain diffusion and mobility from Einstein's model of diffusion. In brief, the ‘initial’ and ‘final’ states in the two-state model of polaron hop correspond to localized spin densities on the ‘initial’ or ‘final’ sites, respectively, with the structures being optimized to give the lowest energy structure of the ‘localized’ polaron states. The spin density arises from one excess electron for electron hops and the removal of one electron from the supercell for hole hops. The approximate reaction pathway is generated by morphing linearly the ‘initial’ polaron stable structure and the ‘final’ polaron structure. DFT energies are calculated for several points along the pathway, possibly using the previous electron density as the initial guess of the DFT calculation to facilitate maintaining the localized spin character which we monitored. Such energy curves correspond to the ‘quasi-diabatic’ energy states when the localized spin character is maintained along the pathway and most importantly at the mid-point (crossing point). When the unpaired spin is shared by the ‘initial’ and ‘final’ sites at the mid-point, the curves correspond to the ‘adiabatic’ energy states. Calculated polaron curves for doped systems are shown in the ESI.

In our earlier work,16 we stated and showed that the electronic coupling elements VAB for all electron hops in bulk un-doped BVO are small, consistent with the non-adiabatic character of the Marcus curves that exhibit a cusp and a crossing region. We adopted the same approximation for electron hops in doped BVO. Dopant atoms were found to be essentially inaccessible as electron localization sites, with the free energy increases for electron localization on W/Mo sites being ∼26.12 kBT and ∼3.15 kBT, respectively, relative to their first shell V sites, so that hops to W/Mo sites could be taken out the list of KMC processes (for W) or treated with negligible VAB coupling (for Mo).

We note that, for organic semiconductors (like pentacene and other), dynamic disorders arising from thermal effects have been shown to lead to large fluctuations in structures and consequently in electronic coupling and that carrier transport can be significantly affected by such fluctuations in some cases.42 Unlike in these cases where the molecules in organic crystals are held together by weak interactions, bonding in inorganic crystals like BVO is much stronger and structural deformations due to temperature effects are expected to be much less significant than in organic crystals. It transpired also from our earlier work16 that the reorganization energy (in Marcus/Holstein parlance) is large (λ ∼ 1.48 eV based on ΔG* ∼ 0.37 eV). The magnitude of the electronic coupling VAB (∼600 cm−1 or less) is thus much smaller than the reorganization energy, and thermally-induced fluctuations of VAB are unlikely to have a significant qualitative effect on the predicted diffusivity and conductivity, in support of the applicability of the Marcus/Holstein small polaron hopping formalism.

c. KMC modeling

For an ensemble of polarons, each of which may undergo a number of hops, it is necessary to solve the coupled kinetics to describe the transport. To this end we implemented a lattice-based kinetic Monte Carlo simulation code.43–45 Lattice-based KMC is an approach that is particularly well suited for describing the dynamics of space-charge distribution in complex systems such as those considered here. In brief KMC involves creating a list of processes (eight possible nearest-neighbor hops for each electron polaron), each process having its own rate (eqn (2)). KMC uses a
image file: c8ta09899a-t1.tif(2)
stochastic algorithm to solve the coupled rate equations, leading to spatial and temporal characterization of the charge carrier distribution. Starting from a random distribution of a set of electron polarons across vanadium sites, the electrons ‘move’ stochastically in the presence of the other electrons during a simulation. The rate constants used at the start of the KMC simulations are taken from the DFT calculations described above. At each time step, the spatial distribution of the electrons affects differentially the vanadium site potentials, the polaron energies, and therefore the hopping rates. The site potentials are calculated through the Ewald summation approach because of the long-range nature of electrostatic interactions. The change in site relative energy for a polaron hopping from an ‘initial’ site to a ‘final’ site affects the activation energy for that hop, in accord with Marcus' expression in eqn (3) or the Bell–Evans–Polanyi principle,46,47 and, as a
image file: c8ta09899a-t2.tif(3)
consequence, its associated rate. As the space-charge distribution evolves, so do all the individual hopping rates. A detailed description of the program can be found in the Appendix. For a given concentration of charge carriers, the initial spatial distribution of the carriers and of the dopant sites is randomly generated. The diffusivity is calculated from the average of the mean square displacements (MSDs) of each electron (or hole), evaluated for 100 trajectories of 10 milliseconds each.

Special considerations are in order for doped systems. For un-doped ms-BVO, atoms are used as lattice sites carrying a charge equal to the full oxidation state of the constituent elements (Bi, V, and O have oxidation states +5, +3, and −2, respectively). When a dopant W/Mo is introduced as a substitution for a V atom, we learn from the QM calculations (see below) that the W/Mo substitution creates a region of small stabilization for electron localization right around the site of substitution, extending up to ∼3 nearest neighbor shells. In other words, the sites belonging to the 4th and higher nearest neighbor shells are treated as unaffected and belonging to the ‘bulk’ region, vis-à-vis the substitution site. At the same time, electron localization at the W/Mo site of substitution is made virtually impossible as observed in earlier studies.18,19,22 Accordingly the ‘coordination’ states of the W or Mo atoms are +5, with each dopant atom giving away an itinerant excess electron to the lattice. This observation made it possible to incorporate into the KMC model the effect of substitutions by W/Mo atoms by simply keeping the lattice site charges as +5 at the dopant sites, while adding one excess electron to the system. In this fashion the energy of any polaron distribution (calculated by the Ewald summation) would remain unchanged if it was not for electron polarons being localized on V sites inside the region of interaction of the dopant atom. The interaction region was then modeled by adding a (small) stabilization energy contribution to V sites inside the region. This is the quantity ΔG0 in the Marcus/Holstein expression of eqn (3). The modified activation energy ΔG* is given by eqn (3) and the altered rate is given by eqn (2). The magnitude of the stabilization energy was derived from DFT calculations for W/Mo doped supercells. For W dopants the W site was found to be very high in energy and the site was removed as a c site to which electron polarons might hop (in the list of possible KMC processes). For Mo dopants, the Mo site was given a relative energy of ∼+0.012 eV, giving it an ∼5% probability of holding an excess electron (see results below). All in all, doping was captured in the KMC model without modifying the electrostatics of the system lattice in terms of lattice site charges.

In this paper we studied the electron transport characteristics (mobility and conductivity) over a range of doping levels. Each data point in the figures represents the average quantity over 100 trajectories, with each trajectory initialized from a unique spatial distribution of dopant sites in the lattice. The standard error of mean is shown as the error bar in these figures. The lack of overlap between the error bars confirms the statistical significance of the trend. The observation of small magnitudes of uncertainties at each data point should not be misinterpreted as absolute independence of carrier transport characteristics over spatial distribution of dopant sites. The observation only confirms the independence of the characteristics from spatial distribution as long as the distribution offers uniformity. The finding that carrier transport can be controlled by engineering the spatial distribution of doping concentrations as in homo-junctions of gradient (W/Mo)-doped BVO will be the subject of a future publication.

III. Results and discussion

a. Electronic structure in doped BVO

We used DFT + U to calculate the relative energies of structures when the excess electron associated with the W/Mo atom was localized at the dopant site or at a nearest neighbor site or at a 2d nearest neighbor site or at a 3d nearest neighbor site, or further away. We used a 3 × 3 × 1 super cell, with 36 W (Mo) or V sites available for electron localization including the dopant site. The relative energies are displayed in Fig. 2. The calculations reveal a relative stability of charge localization in the nearest neighbor shells compared to the bulk of the lattice. The zero of energy corresponds to the excess electron localized beyond the 4th shell. When the electron is localized at the W site, its relative energy is ∼+0.660 eV (or ∼25.5 kBT with kBT = 0.0259 eV at 300 K), and in the case of Mo, the relative energy for the electron to be at the Mo site is ∼+0.079 eV (or ∼3.0 kBT). Since probabilities of energy states accessible to the electron localization are determined by Boltzmann distribution (eqn (4) and (5), we may conclude that the probability of accessing the energy state with the electron localized on the W site is essentially zero, while that with the electron localized on the Mo site is only ∼5%.
image file: c8ta09899a-t3.tif(4)
dEi = EiEbulk(5)

image file: c8ta09899a-f2.tif
Fig. 2 Relative DFT + U energies when the excess electron from W/Mo is localized within vanadium shells away from the dopant site.

For W, the electron localized on a nearest neighbor site (1st shell) has a relative energy of ∼−0.017 eV (or ∼−0.65 kBT); it is the most stable shell. For Mo, the most stable shell is the 2nd shell with a relative energy of ∼−0.015 eV (or ∼−0.60 kBT). The energy profile in Fig. 2 is semi-quantitatively similar to the one presented by Wu and Ping.22 The shape of the energy curve in Fig. 2 is suggestive of a region of stability, an interaction region with the stabilization energy tailing off to essentially zero beyond the 4th shell. The most stable shell for W-doped BVO is the first nearest neighbor shell while it is the second shell for Mo-doped BVO.

Note that the index in Fig. 2 is a shell index (0 is the dopant site, 1 is the nearest shell, 2 is the next nearest shell…). The distribution of stabilization energies in the doping region (on sites around the dopant site) cannot be represented in a similar fashion when using other parameters such as site distances, lattice directions, or connectivity between the dopant site and electron localization sites. Indeed there exist deviations in the distribution of stabilization energies for sites within shells. Fig. 2 arises from averaging the stabilization energies over sites within shells. For simplicity we defined a stabilization energy for each shell through three steps: (a) we calculated the Boltzmann population for each atomic site based on its DFT energy viaeqn (4); (b) we averaged the site populations over the sites within shells; (c) we defined, viaeqn (6), a shell energy (shown in Fig. 2) that yields the average population, giving rise to a shell probability viaeqn (7). Here N is the total number of metal sites in the model.

image file: c8ta09899a-t4.tif(6)
image file: c8ta09899a-t5.tif(7)
image file: c8ta09899a-t6.tif(8)

In this model, all sites within a shell are assumed to be equally accessible (they are equal energy states, i.e., degenerate states). Through fitting the data for a number of sites within shells (shell j degeneracy Ωj) up to j = 4, the shell j degeneracy Ωj in ms-BVO is reasonably well given by Ωj = 6.429j2 + 0.486j + 1.057. The lattice sites may be categorized into three regions based on the shell index j, such as the dopant site (j = 0), interaction region (0 < j < 4), and bulk region (j ≥ 4). The distribution of shell probabilities as derived from eqn (7) is shown in Fig. 3 for un-doped and doped systems within jmax = 7. While the degeneracy Ωj is responsible for the large shell probability in the bulk region, it can be seen in Fig. 3 that the shell probabilities corresponding to the interaction region are larger for a doped system than for an un-doped system. The larger probabilities indicate an increased accessibility of the electron localization states within the interaction region to the detriment of the dopant site and the ‘bulk’ sites. The nature of the trends shown in Fig. 3 remains unaffected with system size (increasing shell index maximum jmax); however, the difference in the bulk region probabilities between doped and un-doped systems becomes more insignificant, suggesting a shrinking influence of a single dopant with increasing size of the simulation supercell. At the same time, the kinetics in the bulk region are virtually unaffected since all the sites in the bulk are essentially degenerate.

image file: c8ta09899a-f3.tif
Fig. 3 Shell probability distribution around a dopant site in BVO. The increase in probability (highlighted with a yellow oval) for the shells around the dopant site (shell ‘0’) corresponds to the stabilization region.

Further discussion about shell probability distribution is included in the ESI. Excess electron localization sites upon Mo or W single doping and co-doping are illustrated in Fig. S1(a) and (b) for selected sites. In all cases, the excess electrons are well localized as confirmed by the spin density plots shown in Fig. S2. Upon W/Mo co-doping, there are two excess electrons. The DFT calculations shown in Table S1 indicate that the two electrons prefer to localize on two different V sites. A second configuration has one excess electron localized on a V site and the other electron on the Mo site with a probability of ∼3%. Interesting as well is the finding that the two electrons can also form a bi-polaron localized on one single V site, although with a low probability of ∼1%.

b. Hopping pathways in doped BVO

As indicated above, dopant atoms give rise to an interaction region around the doping site whereby the excess electron exhibits a small stabilization compared to the ‘bulk’ region. With Mo doping there is a small probability that the excess electron localizes at the Mo site. An excess electron never resides on a W site. In co-doping the stabilization energies are essentially additive.

We investigated the effect of the presence of doping atoms on activation barriers for hops in the neighborhood of dopant atoms. V-to-V hops further away from the dopant atoms than the 4th shell of V atoms do not feel the presence of the dopant atoms and the hopping barriers are unchanged. We selected four V atoms close to the dopant site and determined the Marcus/Holstein energy profiles for hops involving these atoms. The curves and the energies are presented in Fig. S3–S10 and Tables S3 to S10, with some of them involving Mo-to-V hops since Mo atoms have a finite probability to receive an excess electron. These curves complement the exhaustive set of curves reported in our earlier work.16 The overall outcome is that W and/or Mo doping does not affect the characteristics of hops outside the interaction region created by the dopant atom. Within the interaction region V sites have varying but always small stabilization energies that can very well be accounted for by the quantity ΔG0 in the Marcus/Holstein expression in eqn (3) or by the Bell/Evans/Polanyi rule46,47 on how they affect the activation barriers. In the same vein, changes in energy characteristics for hops involving a Mo atom can be accounted for in the same fashion. Considering this small to negligible probability associated with electron transfer processes involving dopant sites, we did not calculate the VAB for these processes. Besides, it is understood that changes in nuclear vibrational modes as a result of thermal fluctuations or crystallographic defects may modify the electron–phonon coupling.48 However, the nuclear vibrational modes in covalently-bound inorganic semiconductors are largely unaffected particularly at lower doping concentrations49 under purview in this study (doping concentrations up to 1%). Thus, we continued to use ΔG* ∼ 0.37 eV for electron hops to derive the reorganization energy (λ) from eqn (3) to be ∼1.48 eV.

c. Electron transport in doped BVO

We defined a model system for W- or Mo-doped BVO, based on the (9 × 9 × 4) supercell we used previously for undoped BVO.16 The number of dopant atoms was 0, 1, 2, 3, and 4 W/Mo atoms within the supercell, corresponding to dopant concentrations up to ∼0.309% (in the range of experimental levels). The number of excess electrons in addition to those introduced by the dopant atoms was kept constant at 4 electron polarons (arising from light absorption for example). The KMC data are shown in Fig. 4 for electron diffusivity for un-doped, singly-doped (with either W or Mo) and co-doped (doping ratio of W to Mo at 1[thin space (1/6-em)]:[thin space (1/6-em)]1) BVO. For the one electron polaron system, the diffusion coefficient (after correction for the difference in activation energy between computation and experiment) is only a factor of ∼5 smaller than the experimental value measured by drift mobility in the presence of an electric field. It can be seen that the average diffusion coefficient exhibits a decrease (cubic dependence) as a function of the electron polaron concentration. Indeed electrons in a multi-electron simulation are not completely free to drift. A move closer to other electrons is less favorable than a move away from other electrons.
image file: c8ta09899a-f4.tif
Fig. 4 Comparison of electron polaron mobility and conductivity for un-doped and doped BVO systems under the ‘3-shell’ implementation, referring to the number of V shells comprising the small stabilization region around each dopant atom.

When considering 3 nearest neighbor shells as the region of stabilization, every dopant site affects the stability of 97 sites including itself in the simulation supercell size of 9 × 9 × 4. Because of this we were able to explore doping levels only up to 0.386% (5 substitutions among 1296 V sites; here data up to only 4 substitutions are shown for easier comparison against co-doped simulations) if we assumed non-overlapping stabilization regions. In order to achieve larger doping levels, we repeated the simulations by considering 2 nearest neighbor shells as the region of stabilization and a total of 37 sites belonging to the region of influence (thus the ‘bulk’ region started from sites belonging to the 3rd nearest neighbor shell and beyond). This allowed us to explore doping levels up to 1% (13 substitutions among 1296 V sites). Mobilities and conductivities for the higher doping levels and still four additional excess electrons are shown in Fig. 5. Because the stabilization energies are small, we observed only a marginal difference in carrier mobilities between the 2-shell- and 3-shell-model implementations. The comparative data between the 2-shell and 3-shell models are presented in Fig. S11 in the ESI.

image file: c8ta09899a-f5.tif
Fig. 5 Comparison of electron polaron mobility and conductivity for un-doped and doped BVO systems under the ‘2-shell’ implementation, referring to the number of V shells comprising the small stabilization region around each dopant atom.

When comparing the electron mobilities between (singly- and co-doped) doped and un-doped BVO, we observe a marginal decrease (factor of ∼1.1) in the value of carrier mobility with an increase in the doping level. This observation is at odds with the interpretation of a factor of ∼20 decrease in carrier mobility upon introduction of 1% W doping.50 However our data support experiment in that the intrinsic mobility of electrons in doped BVO is smaller than that in un-doped BVO.

Furthermore, we performed co-doping simulations with a doping ratio of W to Mo of 1[thin space (1/6-em)]:[thin space (1/6-em)]1. For these simulations, we continued to ensure the non-interaction of any two doping regions to avoid complexities in the distribution of energies for electron localization when two or more doping regions overlap. When electron mobilities between singly-doped and co-doped BVO were compared, we noticed that the mean values of singly-doped mobilities (W and Mo) and co-doped mobilities are very comparable and within a factor of 0.99–1.01 from each other. However, owing to the shallow nature of the stabilization region, the electron mobilities of W/Mo singly-doped BVO are very close (differ by a maximum factor of only ∼1.006 for 12 substitutions). While simulations at higher doping concentrations are necessary for stronger evidence, current observations for co-doping may suggest a linear mixing of the individual contributions from singly-doped systems towards influencing electron transport.

IV. Conclusions

We presented mesoscale characterization of electron transport in W/Mo-doped BVO. The mesoscale KMC model, supported by first principles-determined electron hopping rates, captures the complex dynamics of electron carriers arising from light absorption and from metal doping. The computations and simulations support the observation that, for the doping level used in experiment, doping atoms do not affect significantly the transport dynamics of the charge carriers compared to stoichiometric BVO. W/Mo doping increases the electron carrier concentration and consequently the conductivity in the electrode. DFT calculations of W/Mo-doped BVO show a shallow interaction region around the doping atom where the excess electron localized on a V atom is much more stable than when localized on the doping atom, and slightly more stable than when the excess electron is localized far away from the doping atom. The calculated V-to-V hopping barriers are little affected within this interaction region and not affected outside this region. Mesoscale KMC simulations indicate that the stabilization regions have little effect on the diffusivity compared to the stoichiometric system because of the smallness of the stabilization. The diffusivity is found to decrease slightly with an increasing number of carriers, but the conductivity of a system with electron polarons arising from doping together with light absorption increases compared to that of the un-doped system. These findings are consistent with experimental observations. This work will set the foundation to study electron transport in gradient (W/Mo)-doped systems and other mixed phase systems.

V. Appendix

a. Program description

We developed a cross-platform object-oriented library of python modules and scripts called PyCT (Python-based Charge Transport) to perform lattice-based KMC simulations of charge carrier kinetics in crystalline semiconductor materials on mesoscopic length and time scales. The library is hosted as a private repository on GitHub. Interested users can contact one of us (VP) to get access to the repository. The workflow of the KMC code is depicted in Fig. 6. The code consists of two essential components, a ‘Simulation System’ and a ‘Simulation Engine’.
image file: c8ta09899a-f6.tif
Fig. 6 Workflow of the Python KMC code—PyCT.

The ‘Simulation System’ includes two object definitions: the first called ‘Material’ with all physical parameters related to the material of interest, including the system size, atomic coordinates, and charge hopping parameters (reorganization energy λ and electronic coupling VAB from Marcus/Holstein theory of the hopping model.51) The second object is ‘Neighbors’ that contains functions to generate the cell lists for hopping neighbors and atom pairs interacting within the short-range part of the Ewald summation.

The second essential component of the code, the ‘Simulation Engine’, includes a collection of objects related to the Variable Step Size Method (VSSM) implementation of KMC depicted in Fig. 7. The charge transfer reactions are modelled between neighboring sites via small polaron hopping steps.36 The charge transfer parameters for the elementary processes originate preferably from first-principles calculations or from experiment. The long-range coulombic interactions of polaron–polaron and polaron–lattice are computed by the Ewald summation technique. Since the KMC simulations are performed over a static lattice, we derived for computational efficiency an expression to compute directly the key quantity of interest, the change in free energy (ΔG0). This strategy reduces the computational scaling from O(N2) to O(N), where N is the total number of ions in the system. Instantaneous kinetic rates are computed using the activation energies as estimated from the change in (ΔG0) using the Brønsted–Evans–Polanyi relation52 or Marcus/Holstein expression. Since the spatial distribution of charges changes at each time step, the quantities (ΔG0) and the associated process rates k must be computed for all possible processes at each instant of time.

image file: c8ta09899a-f7.tif
Fig. 7 Schematics of the Variable Step Size Method (VSSM) algorithm for KMC simulations.

From eqn (3), the quantity ΔG0, the instantaneous energy ‘defect’ for a given hop, consists of two components: (1) one arising from the total electrostatic energy of the instantaneous charge distribution and calculated by the Ewald formalism and (2) the shift applied from the electron localization energetics due to doping atoms:

ΔG0 = ΔG0,Ewald + ΔG0,Shift(A1)

For a charge transferring from site k to site k′, the change in the system energy is given by

image file: c8ta09899a-t7.tif(A2)
image file: c8ta09899a-t8.tif(A3)

In doped-BVO systems, an interaction energy is given to electrons within the interaction region of a dopant atom. If the instantaneous energy stabilization contribution arising from the carriers within the interaction region of any dopant is denoted by δE(qi), upon a change of charge configuration associated with a KMC move (an electron may move into an interaction region or move out of an interaction region), the configuration stabilization energy ΔG0,Shift is given by

ΔG0,Shift = δE(qnewi) − δE(qoldi)(A4)

The quantity ΔG0 is computed within the function ‘get_process_rates’ that returns the kinetic rate of the process under consideration via the two equations

image file: c8ta09899a-t9.tif(A5)
image file: c8ta09899a-t10.tif(A6)

The choice of process to be enacted is made by generating a random number x1 within the function ‘do_kmc_steps’ as defined below:

image file: c8ta09899a-t11.tif(A7)
image file: c8ta09899a-t12.tif(A8)

The code is outlined in Fig. 8.

image file: c8ta09899a-f8.tif
Fig. 8 PyCT – KMC code outline; eqn (A1) to (A8) are shown next to the functions where they are implemented.

b. Further discussion on probability distribution

The nature of internal distribution of probabilities within the lattice may be quantified by segmenting between the pre-defined regions as presented in Table 1. The net increase in the probability for the doping region (dopant site + interaction region) upon doping is compensated for through a decrease in the probability for the bulk region.
image file: c8ta09899a-t13.tif(A9)
Table 1 Region-specific probabilities in un-doped and W/Mo-doped systems for jmax = 7
Dopant site Interaction region Bulk region
Undoped 0.001 0.104 0.895
W-Doped 0.000 (−0.001) 0.145 (+0.041) 0.855 (−0.040)
Mo-Doped 0.001 (−0.000) 0.148 (+0.044) 0.851 (−0.044)

The change in shell accessibility upon doping may also be noticed through the ratio of doped to un-doped probabilities (Pdoped/Pundoped) as defined in eqn (A9) with the quantity plotted against shell index j in Fig. 9. The zero or near-zero values of the ratio at the dopant site (shell 0) arise from the positive relative energy for electron localization at the dopant site, as discussed in the main text. In the interaction region the ratio is larger than 1.0, indicating increased stability for electron localization upon doping. The ratio peaks at j = 1 for W-doped and j = 2 for Mo-doped systems, consistent with the energetic distribution shown in the inset of Fig. 2. Owing to the zero stabilization energies of the ‘bulk’ region, the ratio remains constant at 0.828 for W-doped and 0.813 for Mo-doped systems. The decrease in the ratio for the bulk region from unity is a consequence of the probability distributions within the lattice as shown in Table 1. The constant value across the bulk shells suggests that the charge transport kinetics in this region are unaffected by the dopant site and continue to resemble that of an un-doped system.

image file: c8ta09899a-f9.tif
Fig. 9 (a) Ratio of doped to un-doped probabilities plotted against shell index j and (b) bulk value of probability ratio plotted against jmax.

The probability ratio for the bulk region, as shown in Fig. 9b shows an asymptotic behavior with increasing maximum shell index (jmax). Such quick convergence towards unity further reinforces the finiteness of the impact of the W/Mo-doping over the neighboring bulk in the ms-BVO lattice.

The glossary of the variables used in the above equations includes

qc: charge of the carrier (−1 for electron)

qoldi: charge configuration in the old state where the charge is located at site k.

qnewi: charge configuration in the old state where the charge is located at site k′.

α: Ewald summation parameter (α = 0.18 is used for computational efficiency).

image file: c8ta09899a-t14.tif: vector connecting ion i to ion j.

V: system volume.

[n with combining macron]: lattice vector in real space.

k: lattice vector in Fourier space.

δE(qi): energy shift to be applied corresponding to the charge configuration qi

Conflicts of interest

The authors declare no competing financial interest.


V. P. and T. L. contributed equally to this work. V. P., T. L., C. L., and M. D. designed the study and the calculations and contributed to the interpretation of the data. V. P., T. L. and, M. D. prepared the manuscript. We thank Dr Andrew Schultz for valuable discussions. The work was supported in part by the National Natural Science Foundation of China (grant # 21703054). M. D. also acknowledges funds from the University at Buffalo and support from the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Award Number(s) DE-SC0019086 (interpretation of the data and writing of the manuscript).


  1. K. Sivula and R. van de Krol, Semiconducting materials for photoelectrochemical energy conversion, Nat. Rev. Mater., 2016, 1, 15010 CrossRef CAS.
  2. M. G. Walter, et al., Solar Water Splitting Cells, Chem. Rev., 2010, 110(11), 6446–6473 CrossRef CAS PubMed.
  3. F. A. Fatwa and P. B. Sean, Recent developments in complex metal oxide photoelectrodes, J. Phys. D: Appl. Phys., 2017, 50(19), 193002 CrossRef.
  4. M. Woodhouse and B. A. Parkinson, Combinatorial approaches for the identification and optimization of oxide semiconductors for efficient solar photoelectrolysis, Chem. Soc. Rev., 2009, 38(1), 197–210 RSC.
  5. Y. Liang, et al., Highly Improved Quantum Efficiencies for Thin Film BiVO4 Photoanodes, J. Phys. Chem. C, 2011, 115(35), 17594–17598 CrossRef CAS.
  6. J. K. Cooper, et al., Indirect Bandgap and Optical Properties of Monoclinic Bismuth Vanadate, J. Phys. Chem. C, 2015, 119(6), 2969–2974 CrossRef CAS.
  7. T. W. Kim and K.-S. Choi, Nanoporous BiVO4 Photoanodes with Dual-Layer Oxygen Evolution Catalysts for Solar Water Splitting, Science, 2014, 343(6174), 990–994 CrossRef CAS PubMed.
  8. H. S. Park, et al., Factors in the Metal Doping of BiVO4 for Improved Photoelectrocatalytic Activity as Studied by Scanning Electrochemical Microscopy and First-Principles Density-Functional Calculation, J. Phys. Chem. C, 2011, 115(36), 17870–17879 CrossRef CAS.
  9. Y. Park, K. J. McDonald and K.-S. Choi, Progress in bismuth vanadate photoanodes for use in solar water oxidation, Chem. Soc. Rev., 2013, 42(6), 2321–2337 RSC.
  10. C. Le, et al., Mo-Doped BiVO4 Photoanodes Synthesized by Reactive Sputtering, ChemSusChem, 2015, 8(6), 1066–1071 CrossRef PubMed.
  11. W. Yao, H. Iwai and J. Ye, Effects of molybdenum substitution on the photocatalytic behavior of BiVO4, Dalton Trans., 2008,(11), 1426–1430 RSC.
  12. H. Ye, et al., Rapid Screening of BiVO4-Based Photocatalysts by Scanning Electrochemical Microscopy (SECM) and Studies of Their Photoelectrochemical Properties, J. Phys. Chem. C, 2010, 114(31), 13322–13328 CrossRef CAS.
  13. M. Li, L. Zhao and L. Guo, Preparation and photoelectrochemical study of BiVO4 thin films deposited by ultrasonic spray pyrolysis, Int. J. Hydrogen Energy, 2010, 35(13), 7127–7133 CrossRef CAS.
  14. W.-J. Yin, et al., Doping properties of monoclinic BiVO4 studied by first-principles density-functional theory, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83(15), 155102 CrossRef.
  15. D. K. Zhong, S. Choi and D. R. Gamelin, Near-Complete Suppression of Surface Recombination in Solar Photoelectrolysis by “Co-Pi” Catalyst-Modified W:BiVO4, J. Am. Chem. Soc., 2011, 133(45), 18370–18377 CrossRef CAS PubMed.
  16. T. Liu, et al., Bimodal hole transport in bulk BiVO4 from computation, J. Mater. Chem., 2018, 6, 3714–3723 RSC.
  17. T. Liu, et al., The nature of photogenerated charge separation among different crystal facets of BiVO4 studied by density functional theory, Phys. Chem. Chem. Phys., 2015, 17(36), 23503–23510 RSC.
  18. V. Jovic, et al., Soft X-ray spectroscopic studies of the electronic structure of M:BiVO4 (M = Mo, W) single crystals, J. Phys. Chem. A, 2015, 3(47), 23743–23753 CAS.
  19. V. Jovic, et al., A soft X-ray spectroscopic perspective of electron localization and transport in tungsten doped bismuth vanadate single crystals, Phys. Chem. Chem. Phys., 2016, 18(46), 31958–31965 RSC.
  20. F. F. Abdi, et al., Efficient solar water splitting by enhanced charge separation in a bismuth vanadate-silicon tandem photoelectrode, Nat. Commun., 2013, 4, 2195–2201 CrossRef PubMed.
  21. L. H. Han, et al., Efficient Water-Splitting Device Based on a Bismuth Vanadate Photoanode and Thin-Film Silicon Solar Cells, Chemsuschem, 2014, 7(10), 2832–2838 CrossRef CAS PubMed.
  22. F. Wu and Y. Ping, Combining Landau Zener theory and kinetic Monte Carlo sampling for small polaron mobility of doped BiVO4. Condensed-Mat.mtrl-sci, 2018, arXiv:1808.02507 Search PubMed.
  23. W. Zhang, et al., Unconventional relation between charge transport and photocurrent via boosting small polaron hopping for photoelectrochemical water splitting, ACS Energy Lett., 2018, 3, 2232–2239 CrossRef CAS.
  24. A. J. E. Rettie, et al., Anisotropic small-polaron hopping in W:BiVO4 single crystals, Appl. Phys. Lett., 2015, 106(2), 022106–022110 CrossRef.
  25. W.-J. Yin, et al., Doping properties of monoclinic BiVO4 studied by first-principles density-functional theory, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83(15), 155102 CrossRef.
  26. K. Ding, et al., Why the photocatalytic activity of Mo-doped BiVO 4 is enhanced: a comprehensive density functional study, Phys. Chem. Chem. Phys., 2014, 16(26), 13465–13476 RSC.
  27. A. J. E. Rettie, et al., Unravelling Small-Polaron Transport in Metal Oxide Photoelectrodes, J. Phys. Chem. Lett., 2016, 7(3), 471–479 CrossRef CAS PubMed.
  28. G. Kresse and J. Furthmuller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54(16), 11169–11186 CrossRef CAS.
  29. G. Kresse and J. Furthmuller, Efficiency of ab initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci., 1996, 6(1), 15–50 CrossRef CAS.
  30. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized gradient approximation made simple (vol 77, pg 3865, 1996), Phys. Rev. Lett., 1997, 78(7), 1396 CrossRef CAS.
  31. P. E. Blochl, Projector augmented-wave method, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50(24), 17953–17979 CrossRef.
  32. G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59(3), 1758–1775 CrossRef CAS.
  33. S. L. Dudarev, et al., Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57(3), 1505–1509 CrossRef CAS.
  34. H. J. Monkhorst and J. D. Pack, Special points for brillouin-zone integrations, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13(12), 5188–5192 CrossRef.
  35. K. M. Rosso and M. Dupuis, Electron transfer in environmental systems: a frontier for theoretical chemistry, Theor. Chem. Acc., 2006, 116(1–3), 124–136 Search PubMed.
  36. K. M. Rosso, D. M. A. Smith and M. Dupuis, An ab initio model of electron transport in hematite (α-Fe2O3) basal planes, J. Chem. Phys., 2003, 118(14), 6455–6466 CrossRef CAS.
  37. N. A. Deskins and M. Dupuis, Electron transport via polaron hopping in bulk TiO2: a density functional theory characterization, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 75, 195212–195221 CrossRef.
  38. N. A. Deskins and M. Dupuis, Intrinsic Hole Migration Rates in TiO2 from Density Functional Theory, J. Phys. Chem. C, 2009, 113, 346–358 CrossRef CAS.
  39. D. Emin and T. Holstein, Adiabatic Theory of Hall Mobility of Small Polaron in Hopping Regime, Bull. Am. Phys. Soc., 1968, 13(3), 464 Search PubMed.
  40. D. Emin and T. Holstein, Studies of Small Polaron Motion. 4. Adiabatic Theory of Hall Effects, Ann. Phys., 1969, 53(3), 439–520 Search PubMed.
  41. T. Liu, M. Dupuis and C. Li, Band Structure Engineering: Insights from Defects, Band Gap, and Electron Mobility, from Study of Magnesium Tantalate, J. Phys. Chem. C, 2016, 120(13), 6930–6937 CrossRef CAS.
  42. L. J. Wang, et al., Multiscale study of charge mobility of organic semiconductor with dynamic disorders, Phys. Chem. Chem. Phys., 2010, 12(13), 3309–3314 RSC.
  43. M. J. Hoffmann, S. Matera and K. Reuter, kmos: A lattice kinetic Monte Carlo framework, Comput. Phys. Commun., 2014, 185(7), 2138–2150 CrossRef CAS.
  44. M. Leetmaa and N. V. Skorodumova, KMCLib: a general framework for lattice kinetic Monte Carlo (KMC) simulations, Comput. Phys. Commun., 2014, 185(9), 2340–2349 CrossRef CAS.
  45. J. Yu, et al., Kinetic Monte Carlo Study of Ambipolar Lithium Ion and Electron-Polaron Diffusion into Nanostructured TiO2, J. Phys. Chem. Lett., 2012, 3(15), 2076–2081 CrossRef CAS.
  46. R. P. Bell, The theory of reactions involving proton transfers, Proc. R. Soc. London, Ser. A, 1936, 154, 414 CrossRef.
  47. M. G. Evans and M. Polanyi, J. Chem. Soc., Faraday Trans., 1936, 32, 1340 Search PubMed.
  48. L. Wang, et al., Multiscale study of charge mobility of organic semiconductor with dynamic disorders, Phys. Chem. Chem. Phys., 2010, 12(13), 3309–3314 RSC.
  49. V.-I. Merupo, et al., Structural and optical characterization of ball-milled copper-doped bismuth vanadium oxide (BiVO4), CrystEngComm, 2015, 17(17), 3366–3375 RSC.
  50. F. F. Abdi, et al., The Origin of Slow Carrier Transport in BiVO4 Thin Film Photoanodes: A Time-Resolved Microwave Conductivity Study, J. Phys. Chem. Lett., 2013, 4(16), 2752–2757 CrossRef CAS.
  51. R. A. Marcus and N. Sutin, Electron transfers in chemistry and biology, Biochim. Biophys. Acta, Rev. Bioenerg., 1985, 811(3), 265–322 CrossRef CAS.
  52. A. W. Hauser, et al., A systematic study on Pt based, subnanometer-sized alloy cluster catalysts for alkane dehydrogenation: effects of intermetallic interaction, Phys. Chem. Chem. Phys., 2016, 18(16), 10906–10917 RSC.


Electronic supplementary information (ESI) available: Free energy data of electron polarons localized on the sites involved in the doping region, corresponding spin density contours, interpolation energy curves for electron hops between sites involved in the doping region, comparison of electron mobility curves between 2-shell and 3-shell models, and the bar plot showing shell-wise residence frequency. See DOI: 10.1039/c8ta09899a
These authors contributed equally to this work.

This journal is © The Royal Society of Chemistry 2019