Molybdenum Defect Complexes in Bismuth Vanadate
aa r X i v : . [ c ond - m a t . m t r l - s c i ] M a y Molybdenum Defect Complexes in Bismuth Vanadate
Enesio Marinho Jr ∗ and Cedric Rocha Le˜ao † Federal University of ABC (UFABC), 09210-580 Santo Andr´e, S˜ao Paulo, Brazil. (Dated: May 18, 2020)Monoclinic bismuth vanadate (BiVO ) is a promising n -type semiconductor for applications insunlight-driven water splitting. Several studies have shown that its photocatalytic activity is greatlyenhanced by high concentrations of Mo and W dopants. In the present work, we performed abinitio calculations to assess the most favorable relative position between Mo-related pairs in BiVO .Surprisingly, we verify that the lowest energy configuration for Mo V pairwise defects in BiVO occurson nearest-neighbor sites, despite the higher electrostatic repulsion and larger strain on the crystallattice. Similar results were observed for W V defect pairs in W-doped BiVO . We show that theorigin of this effect lies in a favorable hybridization between the atomic orbitals of the impurities thatis only verified when they are closest to each other, resulting in an enthalpy gain that overcomesthe repulsive components of the pair formation energy. As a consequence, Mo and/or W dopedBiVO are likely to present donor-donor defect complexes, which is an outcome that can be appliedin experimental approaches for improving the photocatalytic activity of these metal oxides. I. INTRODUCTION
Semiconductor metal oxides have been intensivelyinvestigated for photoanodes in photoelectrochemical(PEC) water splitting cells [1]. Promising materials forenergy conversion and storage through water splittingmust present chemical stability, relatively low cost, suit-able band edge positions, high optical absorption, longlived excitations and large carrier drift lengths [2, 3].Metal oxides are usually resilient to water corrosion andare low cost materials, but also present crucial intrinsiclimitations. For example, TiO with a wide band gaphas low efficiency in absorbing visible light, Fe O de-spite the moderate band-gap, has unfavorable band edgealignment relative to water’s electrolysis potentials [3, 4].Monoclinic scheelite-type bismuth vanadate (BiVO )has emerged as a promising complex metal oxide pho-toanode, since Kudo et al. [5] first reported its high visi-ble light photoactivity. BiVO has been estimated theo-retically to have a potential to harvest up to 11% of thesolar spectrum, delivering a photocurrent of 7.5 mA/cm ,with 9% solar-to-hydrogen conversion efficiency underAM 1.5 sunlight illumination [6]. This owes to BiVO moderate band gap, 2 . − . vs RHE[3, 9], requiring low applied external bias to drive PEC’shydrogen production.Despite these promising properties, the experimentalperformance of pristine BiVO as photoanodes is signifi-cantly limited by some key factors, such as short chargecarrier diffusion lengths and high electron-hole recombi-nation rate [3, 4, 10].Different experimental approaches have been adoptedto overcome these limiting factors of BiVO , including ∗ [email protected] † [email protected] crystal morphology control [11], heterojunctions [12, 13],tandem PEC devices [14], coupling with oxygen evolvingcatalysts [15, 16], and extrinsic doping [17–19].Doping is possibly the simplest attempt to improve theperformance of metal oxides in PECs by supplying ad-ditional free carriers. BiVO has been effectively n -typedoped with tungsten (W) and molybdenum (Mo) [20, 21].When Mo or W is embedded into V sites, even atmoderate concentrations, the photoelectrochemical per-formances of the doped BiVO electrodes has been no-ticeably enhanced [22]. Jeong et al. [23] have experimen-tally demonstrated that there exists an optimal dopingconcentration of Mo and W in BiVO that maximizesthe performance of the photoanode.The reported opti-mal concentrations are fairly high (8% and 10%, respec-tively), indicating that a comprehensive study of com-plexes involving these dopants in BiVO is important torationalize and potentially enhance their beneficial effectsfor PEC hydrogen production. Recently, Pakeetood etal. [24] have reported that Mo V defects do show a ten-dency to group with W V defects in co-doped material,giving rise to donor-donor complexes. The authors arguethat this is an unexpected trend owing to the repulsivecoulomb interaction between the donor defects.In the present work, we studied how Mo-defects in-teract with each other in a BiVO matrix. Using den-sity functional theory we investigated how the electronicproperties of BiVO with pairs of Mo substitutionaldopants depend on the relative position between theseimpurities. Applying the well established formalism offirst-principles calculations for point defects in solids,combined with analysis of the electronic properties, wealso observe that pairs of Mo-defects present lowest for-mation energy on nearest-neighbor sites. The forma-tion energy rises for intermediate distance between thedopants and then drops again as they get farther apart.We show that this surprising behavior can be rational-ized by competing effects between local lattice strain,electrostatic repulsion, and gain of enthalpic stabilitythrough hybridization of the electronic clouds of the de-fects achieved at short distances. II. COMPUTATIONAL DETAILS
Our ab initio calculations were based on the densityfunctional theory [25, 26], as implemented in Vienna abinitio simulation package ( vasp ) [27]. The projector aug-mented wave (PAW) method [28] were used to treat theelectron-ion interaction, and the exchange-correlation en-ergy was described by the generalized-gradient approxi-mation (GGA) as proposed by Perdew, Burke, and Ernz-erhof (PBE) [29].Structural optimizations were achieved using conju-gate gradient algorithm until the Hellmann-Feynmanforces on all atoms reach values lower than or equal to0 .
025 eV/˚A. Kohn-Sham orbitals were expanded into aplane-wave basis set with a cutoff energy of 500 eV. TheBi 5 d s p , V 4 s d , O 2 s p , and Mo 4 s p d s electronswere treated as valence electrons.We consider the base-centered monoclinic primitivecell containing 2 units of BiVO (12 atoms), with sym-metry described by the standard space group C2/ c . Theoptimized lattice parameters for the conventional unitcell were a = 7 .
325 ˚A, b = 11 .
765 ˚A, c = 5 .
179 ˚A, and β = 135 . ◦ , which are in good agreement with previoustheoretical and experimental reports [18, 30]. To studythe Mo-related point defects, we considered a 216-atomsupercell which was built up by a 3 × × primitive cell. The Brillouin zone was sam-pled using a Γ-centered 7 × × k − point grid for thestructural optimization of the BiVO primitive cell, fol-lowing the scheme proposed by Monkhorst and Pack [31],and a reduced Γ-centered 3 × × systems [20, 32]. For this reason, wemodeled the Mo-related point defects only as substitut-ing Mo in V sites (Mo V ).First, we studied a single Mo V defect in BiVO . Weanalyzed the electronic structure and also the equilibriumformation thermodynamics of this point defect, besidesthe study of the intrinsic point defects considering the Bi,V and O vacancies. Furthermore, we studied these sameproperties for pairs of Mo V defect into BiVO , examiningthe influence of the relative distances between substitu-tional V-sites in the electronic and structural properties,choosing substitutional sites 4 ˚A, 7 ˚A, and 10 ˚A far fromthe reference site (Fig. 1). III. RESULTS AND DISCUSSION
The GGA-PBE calculated band structure for theBiVO primitive cell is shown in Fig. 2a. The result Mo Bi V O d = 4 Å d = 7 Å d = 10 Å FIG. 1. Representation of the crystal structures of monoclinicBiVO doped with Mo. The green tetrahedrons highlight theV-sites selected for the substitution with Mo atoms. For pairof Mo V , we considered three relative distances of the substi-tutional V-sites: 4 ˚A, 7 ˚A, and 10 ˚A far apart from the firstV-site. yields a bandgap of 2.08 eV, which is 0.40 eV under-estimated relative to the experimental value [33]. Theband gap is indirect with CBM located at the R pointand the valence-band maximum (VBM) located in the G - L direction. The direct gap is about 0.15 eV largerthan the indirect gap, in excellent agreement with previ-ous experimental and theoretical reports [8, 9].Fig. 2b shows the density of states (DOS) of the BiVO primitive cell, projected onto the respective atomic or-bitals. Both Bi 6 s and O 2 p orbitals contribute to thehighest states of the valence band, whereas the V 3 d orbitals form the states of the bottom of the conduc-tion band. In Fig. 2c, the orbital-resolved DOS of the216-atom BiVO supercell containing a single Mo V isshown. The electronic contribution of orbital Mo 4 d oc-curs mainly in deep energy levels of the valence band, 7eV lower than the Fermi energy, and at the bottom of theconduction band. We do not notice significant changesin the shape of the valence and conduction band edges.Defect formation can be seen as a thermodynamic pro-cess in which atoms and electrons are exchanged betweenthe crystal containing point defects and the pristine onewhich acts as a particle reservoir. The analysis of pointdefect energies provides some fundamental informationabout the optoelectronic properties of semiconductors,such as if a certain defect tends to donate electrons to thesystem, accept electrons or remain neutral in the crystal.The occurrence of deep/shallow energy levels owing tothese point defects can also be assessed.The equilibrium defect concentration is calculated ac-cording to a Boltzmann distribution [34]: n eq N ≈ exp (cid:18) − ∆ H f k B T (cid:19) , (1) where n eq denotes the number of point defects at equi- D e n s i t y o f s t a t e s [ a r b . un i t s ] TDOS Bi 6 p Bi 6 s O 2 p V 3 d Mo 4 d c Energy [eV] (cid:1) (cid:0) (cid:2) (cid:3) (cid:4) (cid:5) (cid:6) (cid:7) a b D e n s i t y o f s t a t e s [ a r b . un i t s ] Energy [eV] (cid:8) (cid:9) (cid:10) (cid:11) (cid:12) (cid:13) TDOS Bi 6 p Bi 6 s O 2 p V 3 d G L R G X-6-5-4-3-2-10123456 E n e r g y [ e V ] k -point b b b G LRX
G: (0.0, 0.0, 0.0)L: (0.5, 0.5, 0.0)R: (0.0, 0.5, 0.0)X: (0.5, 0.0, 0.0)
FIG. 2. (a) Electronic band structure and (b) projected density of states of the BiVO primitive cell. Above the band structure,we depicted the selected high symmetry k -points in the Brillouin Zone. In (c), we presented the projected density of states ofthe 216-atom BiVO supercell with a single Mo V defect. The Fermi level was used as the zero-energy reference in each plot. librium, N is the total atomic sites involved in the defectformation, and ∆ H f is the enthalpy of formation.The formation enthalpy of a point defect in a chargestate q is given by∆ H f = ( E D − E p ) + X i n i µ i + q ( E v + µ e ) + E corr , (2)where E D and E p are the total energy of the supercellwith the defect and of the pristine supercell, respectively; µ i is the absolute value of the chemical potential of atom i ; n i is the number of such defect atoms added ( n <
0) orremoved ( n > q is the charge state of the defect; µ e is the chemical potential of the reservoir with which thesystem exchanges electrons or holes (Fermi level); and E v is the valence band maximum energy.Finally, the last term E corr corresponds to the finitesize correction, which should be included to remove spu-rious electrostatic interactions between the charged de-fect and its periodic images [35, 36] . Makov and Payne[37] (MP) described this correction energy, focusing oncubic cells, as follows: E MPcorr = q α M ǫL − πqQ ǫL , (3)where L = Ω − / is the linear supercell dimension (Ω isthe supercell volume), q is the defect charge state, ǫ isthe macroscopic dielectric constant of the medium, α M is the appropriate Madelung constant for the respectivesupercell geometry, and Q is the second radial momentof the localized charge distribution ρ c : Q = Z V sc r ρ c ( r ) d r . (4) We estimate the leading (first order) correction termin Eq. (3) for 216-atom monoclinic supercells by theEwald method, computing the Ewald energy ( E Ewald ) ofa point charge (H + ) placed into the supercell of interestand scaling the result by the calculated macroscopic di-electric constant. We obtained E Ewald = − .
454 eV, andthe calculated Madelung constant of BiVO supercell was α M = 2 .
8. Furthermore, while ǫ strictly is a tensor, inour calculations we employed the lowest ǫ diagonal ele-ment. For BiVO , our results of dielectric constants were ǫ xx = ǫ zz = 7 . ǫ and ǫ yy = 6 . ǫ , in good agree-ment with other theoretical works [18, 38]. Therefore,we adopted ǫ = 6 . ǫ in Eq. (3) aiming to apply theupper limit of this correction in the calculation of forma-tion energies.Lany and Zunger [39, 40] (LZ) have proposed to calcu-late the second radial moment in Eq. (4) considering thatcharge difference beyond the vicinity of the defect is pre-dominantly described by a delocalized screening chargeof density n s such that ρ c ≈ n s = q Ω (cid:18) − ǫ (cid:19) , (5)and therefore the second radial moment could be calcu-lated, substituting Eq. (5) into Eq. (4). For a generalgeometry, with a = b = c lattice parameters, we havethat Q = (1 / a + b + c ). Using this result in Eq.(4), the image charge correction yields: E LZcorr = (cid:20) − c sh (cid:18) − ǫ (cid:19)(cid:21) q α M ǫL ≡ (1 − f ) q α M ǫL , (6)in which the term c sh is the so-called shape factor. Forinstance, for a cubic cell c sh = π/ α ≈ .
369 [35, 40]. Forthe 216-atom monoclinic supercell of BiVO , we obtain c sh = 0 . , applying theLZ scheme: E LZcorr ≈ . q α M ǫL , (7)and therefore we verified that although the monoclinicsupercell is not approximately isotropic, the computedscaling factor (1 − f ) is in excellent agreement with(1 − f ) ≈ / Bi ),V (V V ) and O (V O ), and the extrinsic point defect con-sidered was the Mo substitutional on V site (Mo V ).Defect formation energies are conventionally definedwith respect to the chemical potential of the elemen-tal solid, as shown by Eq. (2). Numerical values of theatomic chemical potentials depend on the stoichiometricconditions under which the defects are created. In orderto assure the stable growth of the desired compound, wemust impose fundamental thermodynamic conditions toequilibrium chemical potentials, as detailed below [41]:( i .) To avoid precipitation of the atomic phases, thechemical potential of the atoms available to thecrystal growth (the so-called atomic chemical po-tential ) should be smaller than the chemical poten-tial of the respective elemental bulk or gas. Thatis: ∆ µ Bi,V,O ≡ µ Bi,V,O − µ bulk/gasBi,V,O ≤ . (8)( ii .) To maintain the thermodynamic stability of theBiVO crystal growth, the sum of the ∆ µ of thereacting elements must be equal to the heat of for-mation of BiVO :∆ H BiVO = ∆ µ Bi + ∆ µ V + 4∆ µ O , (9)where this heat of formation can be described asfollows:∆ H BiVO = µ bulkBiVO − (cid:2) µ bulkBi + µ bulkV + 4 µ gasO (cid:3) , (10)and each of these terms can be calculated by first-principles approach given the following definitionof chemical potential: µ bulk/gas = E total N formulas . (11)( iii .) The chemical potentials are further restricted byrequiring that other possible competing phases are not formed. In the present work, we consideredthe following competing phases: Bi O , V O andVO , as proposed in the Refs. [18] and [20]. Thus,we have that2∆ µ Bi + 3∆ µ O ≤ ∆ H Bi O , (12)∆ µ V + 2∆ µ O ≤ ∆ H VO , (13)2∆ µ V + 5∆ µ O ≤ ∆ H V O . (14)( iv .) Additional constraints must be posed to avoid thepossible formation of compounds containing Moand the forming elements of the host material. Wehave considered the following additional competingphases: MoO , and MoO :∆ µ Mo + 2∆ µ O ≤ ∆ H MoO , (15)∆ µ Mo + 3∆ µ O ≤ ∆ H MoO . (16)Applying the thermodynamic constraints described inEqs. (8)-(16), we achieved the accessible range for theatomic chemical potentials of Bi, V, and O, which is de-picted by the shaded area in the Fig. 3a. The formationenergies of the vacancies and Mo-related point defectsin BiVO were calculated using Eq. (2), considering theconditions “A” (O-poor condition) (Fig. 3b) and “B” (O-rich condition). Figs. 3b and 3c show that the energiesneeded to form V Bi and V V are higher than that to formV O , when E F is close to the top of the valence band(VBM). These differences tend to decrease considerablywhen E F is bordering the conduction band (CBM). Thisis in agreement with the fact that the highest valenceband states are constituted mostly from O 2 p orbitals,and the bottom of the conduction band is formed by V3 d states. In addition, V Bi and V V defects tend to benegatively charged throughout most of the allowed rangeFermi level can assume. V O defects tend to be positivelycharged. For this reason, V Bi and V V can be described ashole-producing acceptors, and V O are electron-producingdonors to the crystal.The position of the electronic state introduced by thedefects relative to the band edges of the host materialcan be estimated by the transition energy ǫ ( q/q ′ ), whichis defined as the Fermi energy at which the charge stateof a given defect spontaneously transforms from q to q ′ [41]. The transition energies of the analyzed point defectsin BiVO are shown in Table I. These results were ob-tained based on the formation energy curves representedin Figs. 3b and 3c. V Bi and V V defects are both deepacceptors, which means that their ionization energies aresignificantly above the thermal energy k B T . The pointdefects V O and Mo V are, in turn, shallow donors, whichmeans that they will be easily ionized and produce freeelectrons in the conduction band of the host material.These results are in agreement with other theoretical re-ports [18, 20] as well as experimental observations [23].As described by Eq. (1), the point defect concentrationvaries with the negative exponential of the formation en-ergy. That is, the lower is the formation energy of a given at distances of: d = 7 Å d = 10 Å d = 4 ÅBi-poor/ V-rich 0 1 2Fermi energy [eV]-4-202468101214 F o r m a t i o n e n e r g y [ e V ] VO V O B i V O (cid:14) μ O = 0 Bi O Condition "B" (cid:15) μ V [eV] (cid:16) (cid:17) (cid:18) (cid:19) (cid:20) (cid:21) a b Condition "A" [O-poor] (cid:22) μ Bi [eV] (cid:23)
12 0 (cid:24) (cid:25) (cid:26) (cid:27) (cid:28)
10 0 1 2Fermi energy [eV]-6-4-202 F o r m a t i o n e n e r g y [ e V ] -0.4 -0.6-0.200.2 F o r m a t i o n e n e r g y [ e V ] Bi-rich/ V-richBi-rich/ V-poor d Condition "A" [O-poor] e Condition "B" [O-rich] Eg PBE Eg PBE F o r m a t i o n e n e r g y [ e V ] c Condition "B" [O-rich] Eg PBE Eg PBE Mo V V O V Bi V V FIG. 3. Accessible range of chemical potentials that stabilize the formation of BiVO and the formation energy of point defectsunder two growth condition. (a) The stable chemical potential region of BiVO in the (∆ µ Bi , ∆ µ V ) plane with ∆ µ O = 0 eV(grey area). The conditions “A”(O-poor) and “B” (O-rich) were chosen to numerically represent chemical environments forcalculating the formation energy of the point defects. The curves in (b) and (c) show the formation energy of vacancies ofvanadium (V V ), bismuth (V Bi ), oxygen (V O ), and single substitutional Mo in V-site (Mo V ) under conditions “A” (O-poor)and “B” (O-rich), respectively. In (d) and (e), the plots show the formation energy of pairwise substitutional Mo V defectsconsidering three relative distances of the substitutional sites: 4 ˚A, 7 ˚A, and 10 ˚A, under conditions “A” (O-poor) and “B”(O-rich), respectively.TABLE I. Transition energies ǫ ( q/q ′ ) of Bi, V, and O vacancies, as well as single and pairwise Mo V substitutional defects inBiVO . The energy values were calculated relative to the CBM for the donors, and to the VBM for the acceptors.Acceptor point defect ǫ (0 / − ǫ ( − / − ǫ ( − / − ǫ ( − / − ǫ ( − / − Bi V ǫ (0 / + 1) ǫ (+1 / + 2) ǫ (+2 / + 3) ǫ (+3 / + 4) ǫ (+4 / + 5)V O V V − .
05 0.30 point defect, the higher will be the prevalence of this de-fect in the crystal. Thus, considering the expected PECproperties of the Mo-doped BiVO photoanodes, the op-timal thermodynamic condition in which we have highconcentration of shallow donor defects and low concentra-tion of deep acceptors is depicted by the condition “B”,with an O-poor crystal growing environment. The high-est concentration of donor defects in photoanodes intro-duces larger amounts of free electrons in the conduction band improving the photoelectrochemical performance.On the other hand, deep levels, known as trap statesare detrimental, since they capture photoexcited chargecarriers, facilitating electron–hole recombination throughShockley-Read-Hall (trap-assisted) recombination [42].We also calculated the formation energy of pairs of Mosubstitutional defects. Mo V -Mo V defect pairs are shallowdouble donors with low formation energies. To investi-gate the interaction between these defects and the result-ing variations in their electronic behavior we consideredthree relative distances ( d Mo-Mo ): 4 ˚A, 7 ˚A, and 10 ˚A,under the growth conditions “A” (O-poor) and “B” (O-rich) (Figs. 3d and 3e, respectively). Interestingly, themost stable pairwise defect configuration was achievedwith the Mo dopants being incorporated into the nearest-neighbor substitutional V-sites, with d Mo-Mo = 4 ˚A. Seo et al. [43] studied the interplay between the N sub-stitutional in O sites (N O ) and oxygen vacancies (V O )in BiVO , and their results have shown a energetic fa-vorable tendency to form N O -V O defect complexes inBiVO . In most of the band gap region, V O defect ex-hibits 2+ charge state, whereas N O defects present mostly1 − charge state. Hence, in this case the electrostatic at-traction of the charged point defects can favor the for-mation of N O -V O defect complex.For the Mo V -Mo V pair this result is surprising, sinceone would expect that the lattice strain induced by thedefects would be largest when they are on nearest neigh-bor sites, leading to a higher formation energy. Simi-larly, the electrostatic repulsion between the positivelycharged donors should drive them apart. We observe,however, that the formation energy for the defect com-plex is the lowest when they are on nearest-neighbor sites.For the intermediate distance configuration, it rises 58meV (7.5%) and then drops again 20 meV (-2.4%) forthe 10 ˚A separation. This hints to the occurrence ofsome energy reduction effect associated to hybridizationof the electronic clouds of the dopants in close proxim-ity, balancing out the larger repulsion and greater strainof this configuration. As a consequence, this implies atendency to form Mo V defect complexes.A similar behavior for Mo-W codeped BiVO has beenreported recently by Pakeetood et al. [24]. The authorsshow that Mo V -W V defect complex is more likely to formthan Mo V -V Bi or W V -V Bi , which are both donor-acceptordefects. To date, there is still a lack of physical under-standing about these outcomes [44].In order to understand this unexpected trend, we es-timate the local stress field through the difference of theaverage distance between atoms in the defective and pris-tine supercell. We used one of the Mo nucleus as the ref-erential, for the defective supercells, and the correspond-ing V-site as the referential for the pristine supercell.These distances were computed within a cutoff radius of10.10 ˚A, considering periodic boundary conditions. Thesum of the distances, for the defective and the pristinesupercell separately, was averaged by the total numberof atoms within the cutoff radius as follows: h d i = X r cutoff d atom /N atoms , (17)where d atom is the atomic spacing and N atoms is the totalnumber of atoms in each evaluated supercell. Therefore,the local strain was estimated by the difference between h d i of the defective ( h d i V ) and pristine ( h d i pristine ) su-percells: h ∆ d i = h d i V − h d i pristine . (18) Considering only the neutral charge Mo V defect pair,we obtained h ∆ d i equal to 0.072, 0.042, and 0.005˚A/atom for d Mo-Mo of 4 ˚A, 7 ˚A, and 10 ˚A, respectively.These results of h ∆ d i are exactly the same if we consid-ered the charged systems, with q = 2+. As expected,we notice that the increase in the separation betweenthe substitutional sites results in a decrease in the localstress. Therefore, lattice strain indeed favors Mo defectsfarther apart, independently from other effects. As dis-cussed above, we found that the lowest formation energy,however, happens when the two defects are incorporatedon neighboring V-sites.Our findings suggest that electronic effects are at playreducing the enthalpy of formation of the defect complexdespite the largest strain and electrostatic repulsion atshortest distance. To test this hypothesis, we first an-alyzed the charge density differences, which is shown inFig. 4 (upper panels). Since our purpose is to comparethe interaction between the Mo atoms in the three config-urations considered, we evaluated the charge distributionas follows:∆ ρ BiVO = ( ρ BiVO +2Mo V ) − ( ρ BiVO + ρ + ρ V ) , (19)in which ρ BiVO +2Mo V , ρ BiVO , ρ , and ρ V are thecharge densities of the BiVO with a pair of Mo V , thepristine BiVO , the two Mo atoms isolated in the simu-lation box, and of the BiVO supercell with the V vacan-cies (Mo incorporation sites), respectively. The ∆ ρ wascalculated for the BiVO supercell with a pair of Mo V defects at distances of 4 ˚A, 7 ˚A, and 10 ˚A from eachother. The results were plotted using vesta code [45].The excess negative charge due to the Mo V defects actsas a perturbation of the local charge density distribu-tions, being localized throughout the Mo-O tetrahedrons,with accumulation mainly around the Mo-O bond. Weobserve electron depletion close to the Mo nuclei. Fig. 4bsignalizes hybridization among the orbitals of the two Modefects, with no indication of nodal points. This is notobserved when the Mo defects are farther apart (Fig. 4cand Fig. 4d).In order to establish on firmer grounds whether elec-tronic hybridization is responsible for the observed ten-dency of Mo V to form stable pairs on neighboring sites,we performed a crystal orbital Hamiltonian population(COHP) analyses. The energy-resolved visualization ofchemical bonding in BiVO with a pair of Mo V was gen-erated using the crystal orbital Hamiltonian population(COHP), as implemented in the lobster code [46–49].This involves a transformation of the plane wave basisset used by vasp , to a localized basis set of Slater-typeorbitals (STO) [50]. The projected density of states isdefined as: PDOS i ( E ) = X n | c ni | δ ( E − E n ) , (20)where c ni are the coefficients of the molecular orbital ex-pansion regarding the linear combination of atomic or- -6-5-4-3-2-10123456 -3 -2 -1 0 1 2 3 E n e r g y [ e V ] − pCOHP [a. u.]Bi − VV − OBi − OMo − OMo − Bi a b c d -8-7-6-5-4-3-2-101234 -3 -2 -1 0 1 2 3 E n e r g y [ e V ] − pCOHP [a. u.] -8-7-6-5-4-3-2-101234 -3 -2 -1 0 1 2 3 E n e r g y [ e V ] − pCOHP [a. u.] -8-7-6-5-4-3-2-101234 -3 -2 -1 0 1 2 3 E n e r g y [ e V ] − pCOHP [a. u.]Mo − Mo FIG. 4. Crystal orbital hamiltonian populations of (a) BiVO pristine, or containing a double Mo V defect separated by (b) 4˚A, (c) 7 ˚A, and (d) 10 ˚A from each other. The charge density differences for each Mo-Mo relative distances are shown abovethe respective COHP plots. Blue and red volumes represent electron accumulation and depletion, respectively, with isosurfaceof 0 . e − ˚A − . The Fermi level was set to zero in the energy axes.TABLE II. Integrated COHP values (in eV) of the respective filled interactions up to the Fermi level for the Mo-relatedinteractions in BiVO with double Mo V defects in different substitutional-site distances d Mo − Mo .Interaction d Mo − Mo bitals ψ n = P i c ni φ i . Thus, the COHP is − COHP ij ( E ) = H ij X n c ni c ∗ nj δ ( E − E n ) , (21)where H ij is the matrix element h φ i | H | φ j i [50]. Theminus signal in the Eq. (21) is a mathematical artifactto describe bonding states as positive and anti-bondingstates as negative results of the COHP projection. Weused the following basis functions for the COHP calcu-lations: Bi 5 d s p , V 3 d s , O 2 s p , and Mo 4 s p d s ,and analyzed the interactions between Bi-V, Bi-O, V-O,Mo-Bi, Mo-O, and Mo-Mo.The COHP plot of pristine BiVO is shown in Fig. 4a.The curves reveal significant bonding interactions be-tween V and O for the filled electronic states, showing that the scheelite structure of the BiVO is stabilizedmostly by V-O bonds. On the other hand, we see bond-ing Bi-O interactions at energies around 4 eV bellow thevalence band maximum (VBM). We verify few antibond-ing states associated to Bi-O interactions near the VBM.Stoltzfus et al. [51] performed a study of the structureand bonding in different metal oxides, including BiVO ,focusing only on Bi O interactions. Our COHP resultsagree with their reports for this interaction. Our analysesalso show that the Bi-V interactions do not provide sig-nificant bonding/antibonding states in the valence band(Fig. 4a).In addition, we analyzed the integration of the Mo-related COHP up to the Fermi energy (ICOHP), whichcan provide useful quantification of Mo-related bondstrength [52, 53]. This furnishes a chemically intuitiveinterpretation of the electronic structure by comparingbond strength among the different elements in the struc-ture [54]. The results of ICOHP are shown in the Ta-ble II, in which the overall bonding and antibonding in-teractions up to Fermi level are represented by negativeand positive values, respectively. We observe a net bond-ing interaction of the Mo defects pairwise when they arein nearest-neighbor site. The bond strength of the Mo-Mo interaction vanishes when the defect pairs becomefarther apart. We also observe net anti-bonding statesbetween Mo and Bi atoms which remain on similar levelsin all cases. Therefore only in the nearest-neighbor con-figuration the net anti-bonding interactions are balancedout by Mo pairwise stabilizing interaction. IV. SUMMARY AND CONCLUSIONS
In summary, we investigated by first-principles calcu-lations the electronic nature and the interaction betweensubstitutional Mo defects in Mo-doped BiVO . Our anal-yses indicated that these shallow donors, that have beenshown to enhance the photoelectrochemical performanceof the material, display an unexpected inflection in their formation energy versus their relative distance curve.The most favorable configuration for Mo V defect pairs inBiVO is when they are on nearest-neighbor sites, despitethe electrostatic repulsion they exert on each other. Oursimulations also confirm that the lattice is more strainedwhen the defects are on nearest-neighbor sites. This isalso at odds with the observed lowest formation energyfor this configuration. The COHP analyses showed thatonly in the nearest-neighbor configuration bonding statesexist between the two Mo atoms, balancing out the elec-trostatic repulsion and lattice strain that would favorconfigurations with the defects farther apart. Our find-ings indicate that BiVO doped with Mo tends to formdefect clusters, pointing out to new research issues thatcan improve the photoactivity of these metal oxides. V. ACKNOWLEDGMENTS
This work was supported by the Brazilian FederalAgency for Support and Evaluation of Graduate Edu-cation (CAPES). We thank L´ıdia Carvalho Gomes andJuan Camilo Alvarez Quiceno for fruitful discussions.Computational resources were provided by the high per-formance computing center at UFABC. [1] A. Fujishima and K. Honda, Electrochemical pho-tolysis of water at a semiconductor electrode,Nature , 37 (1972).[2] J. K. Cooper, S. B. Scott, Y. Ling, J. Yang,S. Hao, Y. Li, F. M. Toma, M. Stutzmann, K. V.Lakshmi, and I. D. Sharp, Role of hydrogen indefining the n-type character of BiVO photoanodes,Chemistry of Materials , 5761 (2016).[3] Y. Yang, S. Niu, D. Han, T. Liu, G. Wang, andY. Li, Progress in developing metal oxide nano-materials for photoelectrochemical water splitting,Advanced Energy Materials , 1700555 (2017).[4] F. F. Abdi and S. P. Berglund, Recent develop-ments in complex metal oxide photoelectrodes,Journal of Physics D: Applied Physics , 193002 (2017).[5] A. Kudo, K. Ueda, H. Kato, and I. Mikami,Photocatalytic O evolution under visible light ir-radiation on BiVO in aqueous AgNO solution,Catalysis Letters , 229 (1998).[6] F. F. Abdi, N. Firet, and R. van de Krol, Efficient BiVO thin film photoanodes modified with cobalt phosphatecatalyst and W-doping, ChemCatChem , 490 (2013).[7] A. Kudo, K. Omori, and H. Kato, A novelaqueous process for preparation of crystal form-controlled and highly crystalline BiVO powderfrom layered vanadates at room temperature andits photocatalytic and photophysical properties,Journal of the American Chemical Society , 11459 (1999).[8] J. K. Cooper, S. Gul, F. M. Toma, L. Chen,Y. S. Liu, J. Guo, J. W. Ager, J. Yano,and I. D. Sharp, Indirect bandgap and opti-cal properties of monoclinic bismuth vanadate, The Journal of Physical Chemistry C , 2969 (2015).[9] A. Walsh, Y. Yan, M. N. Huda, M. M. Al-Jassim, andS. H. Wei, Band edge electronic structure of BiVO :elucidating the role of the Bi s and V d orbitals,Chemistry of Materials , 547 (2009).[10] K. Sivula and R. Van De Krol, Semiconductingmaterials for photoelectrochemical energy conversion,Nature Reviews Materials , 15010 (2016).[11] Y. Zhao, R. Li, L. Mu, and C. Li, Significance ofcrystal morphology controlling in semiconductor-basedphotocatalysis: a case study on BiVO photocatalyst,Crystal Growth & Design , 2923 (2017).[12] K. H. Ye, H. Li, D. Huang, S. Xiao, W. Qiu,M. Li, Y. Hu, W. Mai, H. Ji, and S. Yang, Enhanc-ing photoelectrochemical water splitting by combiningwork function tuning and heterojunction engineering,Nature Communications , 1 (2019).[13] Z. He, Y. Shi, C. Gao, L. Wen, J. Chen, and S. Song,BiOCl/BiVO p–n heterojunction with enhanced pho-tocatalytic activity under visible-light irradiation,The Journal of Physical Chemistry C , 389 (2014).[14] F. F. Abdi, L. Han, A. H. Smets, M. Ze-man, B. Dam, and R. Van De Krol, Efficient so-lar water splitting by enhanced charge separationin a bismuth vanadate-silicon tandem photoelectrode,Nature Communications , 2195 (2013).[15] T. W. Kim and K. S. Choi, Nanoporous BiVO photoan-odes with dual-layer oxygen evolution catalysts for solarwater splitting, Science , 990 (2014).[16] F. F. Abdi and R. van de Krol, Natureand light dependence of bulk recombina-tion in Co-Pi-catalyzed BiVO photoanodes, The Journal of Physical Chemistry C , 9398 (2012).[17] K. P. S. Parmar, H. J. Kang, A. Bist, P. Dua, J. S.Jang, and J. S. Lee, Photocatalytic and photoelectro-chemical water oxidation over metal-doped monoclinicBiVO photoanodes, ChemSusChem , 1926 (2012).[18] G. Wang, Y. Ling, X. Lu, F. Qian, Y. Tong,J. Z. Zhang, V. Lordi, C. Rocha Le˜ao, andY. Li, Computational and photoelectrochemi-cal study of hydrogenated bismuth vanadate,The Journal of Physical Chemistry C , 10957 (2013).[19] T. W. Kim, Y. Ping, G. A. Galli, and K. S. Choi, Simul-taneous enhancements in photon absorption and chargetransport of bismuth vanadate photoanodes for solar wa-ter splitting, Nature Communications , 8769 (2015).[20] W. J. Yin, S. H. Wei, M. M. Al-Jassim, J. Turner,and Y. Yan, Doping properties of monoclinic BiVO studied by first-principles density-functional theory,Physical Review B , 155102 (2011).[21] H. S. Park, K. E. Kweon, H. Ye, E. Paek, G. S.Hwang, and A. J. Bard, Factors in the metal dopingof BiVO4 for improved photoelectrocatalytic activityas studied by scanning electrochemical microscopyand first-principles density-functional calculation,The Journal of Physical Chemistry C , 17870 (2011).[22] W. Luo, J. Wang, X. Zhao, Z. Zhao, Z. Li,and Z. Zou, Formation energy and photoelectro-chemical properties of bivo 4 after doping at bi3+ or v 5+ sites with higher valence metal ions,Physical Chemistry Chemical Physics , 1006 (2013).[23] H. W. Jeong, T. H. Jeon, J. S. Jang, W. Choi, andH. Park, Strategic modification of BiVO for improv-ing photoelectrochemical water oxidation performance,The Journal of Physical Chemistry C , 9104 (2013).[24] P. Pakeetood, P. Reunchan, A. Boonchun, S. Limpi-jumnong, R. Munprom, R. Ahuja, and J. T-Thienprasert, Hybrid-functional study of native de-fects and W/Mo-doped in monoclinic-bismuth vanadate,The Journal of Physical Chemistry C , 14508 (2019).[25] P. Hohenberg and W. Kohn, Inhomogeneous electron gas,Physical Review , B864 (1964).[26] W. Kohn and L. J. Sham, Self-consistent equa-tions including exchange and correlation effects,Physical Review , A1133 (1965).[27] G. Kresse and J. Furthm¨uller, Efficient iterative schemesfor ab initio total-energy calculations using a plane-wavebasis set, Physical Review B , 11169 (1996).[28] P. E. Bl¨ochl, Projector augmented-wave method,Physical Review B , 17953 (1994).[29] J. P. Perdew, K. Burke, and M. Ernzerhof, Gen-eralized gradient approximation made simple,Physical Review Letters , 3865 (1996).[30] L. J. Cheng, C. J. Ping, and L. D. Yu, Crys-tal structure and optical observations of BiVO ,Acta Physica Sinica , 1053 (1983).[31] H. J. Monkhorst and J. D. Pack, Spe-cial points for brillouin-zone integrations,Physical Review B , 5188 (1976).[32] A. J. E. Rettie, H. C. Lee, L. G. Marshall, J. F.Lin, C. Capan, J. Lindemuth, J. S. McCloy,J. Zhou, A. J. Bard, and C. B. Mullins, Com-bined charge carrier transport and photoelectro-chemical characterization of BiVO single crys-tals: intrinsic behavior of a complex metal oxide,Journal of the American Chemical Society , 11389 (2013).[33] J. K. Cooper, S. Gul, F. M. Toma, L. Chen, P.-A. Glans,J. Guo, J. W. Ager, J. Yano, and I. D. Sharp, Electronicstructure of monoclinic bivo4, Chemistry of Materials ,5365 (2014).[34] P. A. Varotsos and K. D. Alexopoulos, Thermodynamicsof point defects and their relation with bulk properties (Elsevier Science, 2013).[35] H. P. Komsa, T. T. Rantala, and A. Pasquarello, Finite-size supercell correction schemes for charged defect cal-culations, Physical Review B , 045112 (2012).[36] C. Rocha Le˜ao and V. Lordi, Ab initio guided opti-mization of GaTe for radiation detection applications,Physical Review B , 165206 (2011).[37] G. Makov and M. C. Payne, Periodic bound-ary conditions in ab initio calculations,Physical Review B , 4014 (1995).[38] Z. Zhao, Z. Li, and Z. Zou, Electronic structure andoptical properties of monoclinic clinobisvanite BiVO ,Physical Chemistry Chemical Physics , 4746 (2011).[39] S. Lany and A. Zunger, Assessment of correction methodsfor the band-gap problem and for finite-size effects insupercell defect calculations: Case studies for zno andgaas, Physical Review B , 235104 (2008).[40] S. Lany and A. Zunger, Accurate prediction of defectproperties in density functional supercell calculations,Modelling and Simulation in Materials Science and Engineering , 084002 (2009).[41] C. Persson, Y. J. Zhao, S. Lany, and A. Zunger,n-Type doping of CuInSe and CuGaSe ,Physical Review B , 035211 (2005), 1504.06521.[42] J. S. Park, S. Kim, Z. Xie, and A. Walsh,Point defect engineering in thin-film solar cells,Nature Reviews Materials , 194 (2018).[43] H. Seo, Y. Ping, and G. Galli, Role of pointdefects in enhancing the conductivity of BiVO ,Chemistry of Materials , 7793 (2018).[44] For completeness, we also calculate the formation energyof W V defect pairs, considering the same relative dis-tances of the substitutional sites: 4 ˚A, 7 ˚A, and 10 ˚A.The W 5 p d s electrons were treated as valence elec-trons. We observe identical behavior between the defectpairs, with lowest formation energy for W V embedded innearest-neighbor V-sites. Regarding the neutral chargedBiVO with pairs of W V , for O-poor growth condition,we obtain ∆ H f = 0.68, 0.74, and 0.72 eV for relative dis-tances of 4 ˚A, 7 ˚A, and 10 ˚A, respectively. Similarly, forO-rich growth condition, we obtain ∆ H f = 1.39, 1.45,and 1.43 eV for 4 ˚A, 7 ˚A, and 10 ˚A, respectively.[45] K. Momma and F. Izumi, VESTA 3 for three-dimensionalvisualization of crystal, volumetric and morphology data,Journal of Applied Crystallography , 1272 (2011).[46] S. Maintz, V. L. Deringer, A. L. Tchougr´eeff,and R. Dronskowski, Lobster: A tool to ex-tract chemical bonding from plane-wave based dft,Journal of Computational Chemistry , 1030 (2016).[47] R. Dronskowski and P. E. Bloechl, Crystal or-bital Hamilton populations (COHP): energy-resolved visualization of chemical bonding insolids based on density-functional calculations,The Journal of Physical Chemistry , 8617 (1993).[48] V. L. Deringer, A. L. Tchougr´eeff, and R. Dron-skowski, Crystal orbital Hamilton population (COHP)analysis as projected from plane-wave basis sets,The Journal of Physical Chemistry A , 5461 (2011). [49] S. Maintz, V. L. Deringer, A. L. Tchougr´eeff,and R. Dronskowski, Analytic projection fromplane-wave and paw wavefunctions and appli-cation to chemical-bonding analysis in solids,Journal of Computational Chemistry , 2557 (2013).[50] S. Tao, I. Schmidt, G. Brocks, J. Jiang, I. Tranca,K. Meerholz, and S. Olthof, Absolute energy levelpositions in tin-and lead-based halide perovskites,Nature Communications , 2560 (2019).[51] M. W. Stoltzfus, P. M. Woodward, R. Seshadri, J. H.Klepeis, and B. Bursten, Structure and bonding inSnWO , PbWO , and BiVO : lone pairs vs inert pairs,Inorganic Chemistry , 3839 (2007).[52] D. K. Mann, J. Xu, N. E. Mordvinova, V. Yan-nello, Y. Ziouani, N. Gonz´alez-Ballesteros, J. P. Sousa, O. I. Lebedev, Y. V. Kolen’ko, and M. Sha-truk, Electrocatalytic water oxidation over AlFe B ,Chemical Science , 2796 (2019).[53] M. Khazaei, J. Wang, M. Estili, A. Ranjbar, S. Sue-hara, M. Arai, K. Esfarjani, and S. Yunoki, Novel MABphases and insights into their exfoliation into 2D MBenes,Nanoscale , 11305 (2019).[54] B. R. Ortiz, L. C. Gomes, J. R. Morey, M. Winiarski,M. Bordelon, J. S. Mangum, I. W. H. Oswald,J. A. Rodriguez-Rivera, J. R. Neilson, S. D.Wilson, E. Ertekin, T. M. McQueen, and E. S.Toberer, New kagome prototype materials: dis-covery of KV Sb , RbV Sb , and CsV Sb ,Physical Review Materials3