2 , 054604 ( 2018 ) Deep vs shallow nature of oxygen vacancies and consequent n-type carrier concentrations in transparent conducting oxides

J. Buckeridge,1,* C. R. A. Catlow,1 M. R. Farrow,1 A. J. Logsdail,2 D. O. Scanlon,1,3,4 T. W. Keal,5 P. Sherwood,5 S. M. Woodley,1 A. A. Sokol,1,† and A. Walsh6,‡ 1University College London, Kathleen Lonsdale Materials Chemistry, Department of Chemistry, 20 Gordon Street, London WC1H 0AJ, United Kingdom 2Cardiff Catalysis Institute, School of Chemistry, Cardiff University, Cardiff CF10 3AT, United Kingdom 3Thomas Young Centre, University College London, Gower Street, London WC1E 6BT, United Kingdom 4Diamond Light Source Ltd., Diamond House, Harwell Science and Innovation Campus, Didcot, Oxfordshire OX11 0DE, United Kingdom 5Scientific Computing Department, STFC, Daresbury Laboratory, Daresbury, Warrington, WA4 4AD, United Kingdom 6Department of Materials, Imperial College London, London SW7 2AZ, United Kingdom


I. INTRODUCTION
The development of materials that combine the properties of electronic conductivity and optical transparency has resulted in a broad range of technologies spanning capacitive sensing [1,2], energy applications [3], optoelectronics [4,5], photocatalysis and gas sensing [6,7], electrochromic devices [8][9][10], and beyond [11,12].Semiconducting oxides, known as transparent conducting oxides (TCOs), have proved to be among the most successful materials in this respect [13][14][15], particularly those based on In 2 O 3 , SnO 2 , and ZnO (and combinations of all three) [16][17][18].As a consequence of their electronic properties and ionic disorder, these wide-gap semiconductors tend to be either intrinsically n-type or easily n-dopable, while lacking any p-type activity, and oxygen deficient under reducing conditions (with ZnO remaining O-substoichiometric across the range of synthesis conditions) [19][20][21].Where the balance lies between the formation of stable free charge carriers and charged point defects, however, and the details of how these aspects combine to effect electronic conductivity and nonstoichiometry [22,23], is not well understood, despite many years of research into this topic.
Of the TCOs, those based on In 2 O 3 have been the most successful for applications, with Sn-doped In 2 O 3 , or ITO, forming the industry standard for touchscreen technology [17].Nominally undoped In 2 O 3 is an intrinsic n-type semiconductor [21].The link between this intrinsic conductivity and ionic disorder has been demonstrated through measurements performed at different oxygen partial pressures, P O 2 [30][31][32].Varying P O 2 will change the concentrations of oxygen vacancies, which may act as donors, and interstitials, which may compensate n-type carriers, in a manner that depends on the defect chemistry of the system.It has been observed that, in undoped In 2 O 3 samples, the conductivity increases significantly as P O 2 is reduced [62].Concurrently, the mobility decreases [21,63], which suggests that it is an increase in vacancy donor concentration, rather than a decrease of the concentration of compensating interstitials, that effects this increase in conductivity.Theoretical studies, however, vary widely in their predictions on the donor properties of vacancies, with some recent DFT studies predicting that vacancies are deep donors [39,42], which would seemingly contradict the experimental findings.One explanation that has been advanced is that, due to a reduction of formation energy, and consistent with observed surface electron densities, vacancies may occur in significant concentrations on the surface, where their ionization energies could be reduced [64][65][66][67][68][69].Such a reduction in ionization energy, however, has not been demonstrated to occur due to the difficulty in calculating the properties of charged defects on surfaces using standard supercell techniques.
SnO 2 is easily n-dopable [70,71], but undoped samples only gain significant carrier concentrations at elevated temperatures [33][34][35] (although Fonstad et al. [72] have shown concentrations of ∼10 17 cm −3 at room temperature).A clear increase in conductivity as P O 2 is decreased can be observed when temperature T > 1000 K, indicating the presence of a relatively deep donor [34].From these data, a donor activation energy of 0.15 eV in the dilute limit has been derived [72][73][74], or possibly at ∼0.3 eV [75], but the majority of computational studies indicates that the most likely donor, the oxygen vacancy, is deeper with the (2+/0) transition occurring at about 1 eV below the conduction band minimum (CBM) [20,41,43].
In this paper we determine the formation and ionization energies of oxygen vacancies in In 2 O 3 , SnO 2 , and ZnO using a hybrid quantum mechanical/molecular mechanical (QM/MM) embedded cluster approach, with a modified embedding po-tential formulated to calculate accurately the electronic properties.We describe the model in full and verify its accuracy through comparison with relevant calculations using more conventional, plane-wave based methods.In all cases we find that the defects are compact donors, but in agreement with Ágoston et al. [41], we find that vacancies in the bulk of In 2 O 3 are shallow donors, which can account for a significant proportion of the observed intrinsic n-type conductivity.For SnO 2 and ZnO we determine that the vacancies are deep, and, depending on the hybrid density functional used, determine transition levels either in excellent agreement with previous studies or slightly shallower (our results in these cases agree well with the experimentally determined ionization energies).We also calculate thermodynamic equilibrium concentrations based on our formation energies, resulting in n-type carrier concentrations in In 2 O 3 in good agreement with experiment and significant carrier concentrations, consistent with experiment, in SnO 2 and ZnO under certain conditions, despite the deep nature of vacancies in those cases.These conclusions can be drawn regardless of which hybrid functional we employ.Our results show that vacancies are crucial for intrinsic conductivity in In 2 O 3 and can contribute significantly in SnO 2 and ZnO, in contrast to previous studies.
The rest of the paper is structured as follows: In Sec.II we discuss our computational approach including a detailed description of the embedded cluster technique we employ, the defect reactions used to determine formation energies, the reference energy in our model and its application to determine charge carrier energies, and the method used to determine defect and carrier concentrations.In Sec.III we present our main findings, which is subdivided into: Sec.III A, where calculated ionization potentials are discussed; Sec.III B, where formation energies calculated using the different density functionals are presented and compared with other computational studies and experiment; Sec.III C, where we discuss the balance between compact and diffuse bound electron states and its significance to interpreting our results; and Sec.III D, where we present our calculated defect and carrier concentrations.Finally, in Sec.IV we summarize our results.

II. CALCULATIONS
We have employed the hybrid quantum mechanical/molecular mechanical (QM/MM) embedded cluster method [96] within the CHEMSHELL [97] package to calculate bulk and defect energies in the three TCOs, modeling each system in its ground state crystal structure (bixbyite for In 2 O 3 , rutile for SnO 2 , and wurzite for ZnO).In this approach, a defect (possibly charged) and its immediately surrounding region, typically of the order of 100 atoms (97 for In 2 O 3 , 89 for SnO 2 , and 86 for ZnO) [98], is treated using a QM level of theory, which is then embedded in a larger cluster (an outer shell of 9704 ions for In 2 O 3 , 10 357 for SnO 2 , and 10 460 for ZnO) treated at the MM level of theory.Embedding effective core potentials (ECPs) are placed at the interface between the QM and MM regions, to prevent spurious electron spillage, while the Madelung potential of the infinite solid is reproduced within the QM region and its immediate surroundings within a spherical active region (of radius 15 Å) by means of fitted point charges surrounding the combined QM/MM cluster.In this approach the defect is modeled at the true dilute limit.The advantages of this method over more commonly employed supercell approaches are: (i) a comprehensive account of polarization effects due to the presence of a charged defect, (ii) the lack of periodic image interactions that would have to be corrected for, and (iii) access to the vacuum level, which allows ionization energies to be determined with an absolute reference.A more detailed account of the approach can be found in Refs.[96,[99][100][101].We note that the method has previously been applied successfully to treat defects and band alignments in wide gap materials as well as systems of reduced dimensionality [23,27,[102][103][104].
For the QM part of the calculation we use density functional theory (DFT) with a triple zeta plus polarization Gaussian basis set for oxygens (Def2-TZVPP) [105] and a double zeta plus polarization set for In, Sn, and Zn cations (cc-pVDZ-PP) [106,107], including a 10 electron effective core potential (ECP) for Zn [107] and 28 electron ECPs for In and Sn [108], as implemented within the GAMESS-UK [109] software.To reduce the computational load, we have removed f functions from the oxygen basis set and some of the highly diffuse functions from the cation basis sets, which do not contribute to the bonding in these ionic solids.The resulting reduced basis set was carefully optimised to reproduce the total energy and eigenstates of the appropriate isolated ions.We have tested these reduced basis sets using a selection of calculated embedded cluster ionization energies and find that their inclusion changes the relative energies by less than 1%.For electron exchange and correlation, we have employed three hybrid functionals: PBE0 [110,111], which is frequently used in plane-wave basis calculations and can serve as a benchmark for our simulations (25% exact exchange); B97-2 (21% exact exchange) [112], which has been fitted to a broad range of thermochemical data, but is similar in form to PBE0 and HSE06 [113], another hybrid functional commonly in use in plane-wave basis calculations; and BB1k [114], which has been fitted to both thermochemical and kinetic data, resulting in a higher proportion of exact exchange (42%).
For the MM part of the simulation we have used pairwise interatomic potentials [115,116] that have been fitted to reproduce the structural, elastic, dynamical, and dielectric properties of the binary oxides.Details on the forcefields can be found in Ref. [25] for In 2 O 3 , Ref. [104] for SnO 2 , and Ref. [117] for ZnO.The GULP [118] software was employed for these calculations.
To treat the interface between the QM and MM regions, typically a large core ECP is placed on cation sites within a few Ångstrom of the edge of the QM region [96].In the present study, however, we found that such an approach resulted in erroneous ionization potentials for the bulk systems.Instead, a specially designed [102,104] local ECP was placed on cation sites within a range of 5 Å from the edge of the QM region.The ECP U p (r) has the form: where the parameters A i and Z i were fitted in order to minimize the gradients on the ions in the QM and interface region and the spread of deep core levels in the energy spectrum, according to the following procedure: An initial guess for the  [120].The final optimized parameters for each cation are given in Table I.The hybrid QM/MM method employed here includes a dielectric continuum approximation, meaning that, beyond a cutoff radius (in this study 15 Å), the MM ions are kept fixed.To account for the polarization of these fixed ions an a posteriori correction is applied using Jost's formula for the energy of polarization (E pol ) outside a sphere of radius R, in a dielectric medium of dielectric constant , due to the presence within the sphere of a charge Q [96]: When determining E pol for fully relaxed structures, we use the static dielectric constant ε in Eq. ( 2).Formation energies are determined (at the athermal limit) by calculating the enthalpy of the appropriate defect reaction: for oxygen rich conditions (high P O 2 , at or above normal atmospheric conditions), where 0 q 2 is the charge of the vacancy defect, q ∈ Z, and for oxygen poor conditions (extreme reducing conditions, with P O 2 practically equal to that in ultrahigh vacuum) [121], where M = In, Sn, or Zn and m = n = 1 for ZnO, m = 1, n = 2 for SnO 2 , and m = 2, n = 3 for In 2 O 3 .The energies of the reactants and products O 0 O and V q+ O are determined from the embedded cluster calculations, the O 2 energy from a molecular calculation using GAMESS-UK with the corresponding basis set and density functional, and the metal and binary oxide standard state energies are derived from experimental heats of formation/atomization (at the athermal limit) [122], as such energies are challenging to calculate in a consistent manner using the QM/MM approach.Using hybrid density functionals, however, one would expect the calculated heats of formation to lie within ∼10% of the experimental values [20,55], which implies that using experimental results will not alter the conclusions drawn from the computed formation energies significantly.
The electronic structure of a typical n-type TCO possesses a sharp contrast in curvature between the conduction and valence band edges; delocalized states constitute the disperse lower conduction band, while the upper valence bands are comprised of more localized states.A disadvantage of the hybrid QM/MM method is the difficulty in reproducing truly delocalized states, such as those necessary to determine electron affinities in n-type TCOs and hence the energy in the reactions above of an electron, e − (determining the energies of holes, however, is straightforward).As the defects studied here, however, all possess compact wave functions that are well confined within the QM region (see below), such electron energies can be approximated with a 'scissors'-like operation, by subtracting from the calculated ionization potential of the perfect system the experimental band gaps of 2.7 eV (In 2 O 3 ) [21], (3.6 eV SnO 2 ) [123], and 3.44 eV (ZnO) [124].One can then introduce the Fermi energy E F as a parameter, relative to the valence band maximum (VBM), of which each defect formation energy is a linear function with slope q, in order to elucidate likely defect properties at different doping conditions.The standard approach to presenting such results is as a graph of formation energy Thermal ionizations are calculated from the formation energies via the reactions: The value of E F which balances Eqs. ( 5), introduced via the e − terms, gives the transition levels from charge state to charge state for each defect.Equilibrium defect concentrations and the self-consistent Fermi energy have been calculated using the code SC-FERMI [125,126], which calculates the defect concentration [V q O ] according to: where N V O is the density of sites on which oxygen vacancies can form, g is the degeneracy of the defect state, k is Boltzmann's constant, and T is the temperature.Electron (n 0 ) and hole (p 0 ) concentrations are calculated as: where E g is the band gap, ρ(E) is the density of states, as the zero of the energy scale).As each concentration is a function of E F , either via the defect formation energy or the Fermi-Dirac function, the equilibrium values can be calculated self-consistently at T , given the constraint of overall charge neutrality [127,128].
The basis sets used in our QM/MM embedded cluster approach do not provide sufficient accuracy to reproduce the upper conduction bands.To calculate ρ(E) for each binary oxide, we have therefore used plane-wave DFT as implemented in the VASP [129][130][131][132] code, utilizing the PBE0 hybrid functional with the projector augmented wave method [133] to model the interaction between core and valence electrons (with twelve, thirteen, and fourteen valence electrons for Zn, Sn, and In atoms, respectively, and six for an oxygen atom).The unit cell structures were taken from the center of the appropriate relaxed QM cluster and a single point calculation was performed using an energy cutoff of 650 eV and 8 × 8 × 8, 12 × 12 × 16, and 16 × 16 × 12 Monkhorst-Pack [134] k-point meshes for bixbyite In 2 O 3 , rutile SnO 2 , and wurzite ZnO, respectively.These criteria allowed convergence in total energy up to 10 −4 eV per atom.The conduction bands were then shifted in energy using a 'scissors operator' in order to reproduce the experimental band gaps.We also performed the calculations using the GGA PBEsol functional, in order to test the effect of variation in ρ(E) with functional on the calculated carrier concentrations, finding that the results changed by less than 1%.As a further test, we included the valence band density of states from the QM/MM calculation, which also had an insignificant effect on the computed carrier concentrations.

III. RESULTS
In this section we first present the computed ionization potentials of the three oxides; next we present and discuss our calculated formation energies and thermal transition levels; then we turn to the subject of compact versus diffuse donor states, before finally giving our results on self-consistent Fermi energies and their consequences for the conductivity properties of the oxygen-deficient binary TCOs.

A. Ionization potentials
The calculated ionization potentials (IPs) are given in Table II, compared with results from previous experimental and computational studies.We note that, when determining the ionization potential from either experiment or from periodic slab-model calculations, the particular surface dipole present due to the crystal termination will affect the ease with which  electrons can be removed, hence the broad range of values observed [148].Our calculated values correspond to 'bulk ionization potentials,' as surfaces are not included in the embedded cluster model.Depending on the surface band bending, the bulk ionization potential can be above or below the surface IP, but would be expected to lie close (within about 1 eV).The agreement, therefore, between our calculations and others in the literature is good.Of the three functionals used, BB1k results in the highest IPs, then PBE0, then B97-2, a trend which reflects the different accounts of electron localization amongst the functionals.We find that the valence band of SnO 2 is deepest with respect to the vacuum level, as the IP is largest.Next is ZnO, then In 2 O 3 , but all three are relatively large, meaning that the VBM will be deep.Consequently, electron carriers will be stable, in contrast to hole carriers, which explains the observed difficulty in p-type doping these systems [149][150][151].

B. Formation energies
Our calculated formation energies and transition levels are shown in Fig. 1 for O-rich and O-poor conditions.As we demonstrate below, our results are in excellent agreement with previous calculations using a supercell approach, where applicable.We note here that the energies we calculate for O-poor conditions are slightly lower than those in many other studies; this difference originates in the values used for the heats of formation of the binary oxides that enter into Eq.( 4).As discussed in Sec.II, we have used experimental values for such heats of formation, while in supercell-based studies they are routinely calculated.In our discussion, the terminology 'shallow' with respect to electronic states is used for thermodynamic levels close to the relevant band edge or resonant within the bands; this term is employed to contrast those states with others deep in the gap.We now discuss each system in turn.

ZnO
For ZnO, we find that the V O is a deep donor of negative U type, having a thermal transition from the neutral to 2+ charge state [ (2 + /0)] deep in the band gap, in agreement with a variety of previous studies (see Fig. 1 and Table III).We calculate very similar values for the formation and transition energies using the PBE0 and B97-2 functionals, which agree well with supercell-based hybrid DFT studies and are slightly higher than those obtained using the generalized gradient approximation (GGA) [53,152].In particular, the results for E f (V 0 O ) in O-rich conditions and for (2 + /0) obtained using PBE0 are in very good agreement with those reported by Alkauskas and Pasquarello [53], who included 32% exact exchange in the hybrid functional rather than the standard 25%, while our calculated (2 + /0) is within ∼0.1 eV of that of Ágoston et al. [41], who also used PBE0.Oba et al. [51] and Frodason et al. [55] found that employing the HSE06 functional (with 37.5% exact exchange in order to reproduce the experimental band gap, which we denote modHSE06 in Table III) results in slightly higher formation energies and deeper transition levels, but overall our results are in reasonable agreement with theirs, which are also very similar to those obtained by Clark et al. [52], who replaced the local density approximation (LDA) exchange with a Thomas-Fermi screened exchange potential (sX functional).Clark et al. [52] also reported values calculated with the HSE06 functional, but these energies were for vacancies in the cubic system, hence we do not include them in Table III (they are nevertheless quite similar to the wurtzite results).The PBE0 and B97-2 O-rich formation energies also This study b Reference [53] c Reference [41] d Reference [51] e Reference [55] f Reference [52] g Reference [42] h References [48][49][50] agree very well with those obtained by Lany and Zunger [42], who used a corrected GGA (corrGGA) approach [153], but their study resulted in a transition level much deeper than our computed value, which also disagrees with the other studies shown in Table III [60,154].We note that, though we obtain a similar transition level, our formation energies are substantially lower than those of Janotti et al. [48][49][50], who used LDA with a Hubbard U parameter, corrected via an extrapolation scheme (corrLDA+U ), but their values also disagree with all other studies shown in Table III.
Although we find good agreement between our computed O-rich formation energies and the others shown in Table III (apart from Janotti et al., [48][49][50] see above), the energies we obtain for O-poor conditions are significantly lower.The difference between the O-rich and O-poor formation energies is the heat of formation of ZnO, H (ZnO), for which we take the experimental value at the athermal limit, H (ZnO) = −3.70eV [122].Calculating this property using the QM/MM approach is prohibitively challenging, but can routinely be determined using common plane-wave techniques (we have chosen to use experimental values as the BB1k and B97-2 functionals are not implemented in common plane-wave codes).The computed value can be less negative than the experimental value by as much as 1 eV, hence the discrepancy between our calculated formation energy in O-poor conditions and the others shown in Table III.There are experimental reports of formation energies in reducing conditions of 1 eV or above [36,88], but the O 2 partial pressure which these This study a Reference [152] b Reference [20] c Reference [41] d Reference [43] measurements were carried out under does not correspond to the extreme O-poor conditions our calculations represent.It is therefore difficult to ascertain the accuracy of our O-poor results in comparison with other DFT studies, but if the values we report are artificially low the conclusions we draw will not be significantly affected, a point to which we return below in Sec.III D. Using the BB1k functional, which includes 42% exact exchange, the computed E f (V 0 O ) is closer to that obtained from modHSE06 [51,55], as expected given that the modHSE06 results included 37.5% exact exchange.With BB1k, however, we find that (2 + /0), while still a deep level, is significantly shallower (at 0.25 eV below the CBM) than that obtained using either PBE0 or B97-2, or those obtained in previous studies, none of which employed a functional of this kind.Deep level transient spectroscopy (DLTS) has been used extensively to study defects in ZnO.There are several levels that are observed in many samples, particulary under suitable conditions for oxygen vacancy formation, within the ranges 0.27-0.31eV and 0.60-0.67eV [93][94][95][155][156][157][158][159][160][161][162].Our computed (2 + /0) using BB1k agrees well with the former DLTS signal and is a transition that we would expect to contribute significantly in the samples used experimentally.With BB1k, we calculate (+/0) = 0.64 eV; the metastable V + O state may form after ionization of the 2+ state, and could account for the DLTS signal within the 0.6-0.67 eV range.Furthermore, we also note that the 0.25 eV transition level is in good agreement with the activation energies (0.2-0.28 eV) derived from Hall measurements by Ziegler et al. [36].

Sn O 2
Similar to the case of ZnO, in SnO 2 we find that the V O is a negative U deep donor, in agreement with previous calculations (see Fig. 1 and Table IV).We find particularly good agreement when comparing our calculations using the PBE0 functional with the plane-wave supercell based study of Scanlon and Watson [20] (who used the same functional), as TABLE V. Calculated formation energies of neutral oxygen vacancies E f (V 0 O ) and the (2+/+), (+/0), and (2+/0) thermal transition levels [ (2 + /+), (+/0), and (2 + /0), referred to the conduction band minimum] in In 2 O 3 from this study using the PBE0, B97-2, and BB1k hybrid functionals, compared with previous calculations (Prev.Calc.) using a variety of functionals (see main text for details).Energies are given in eV.
This study b Reference [56] c Reference [42] d Reference [39] e References [57,58] f Reference [41] g Reference [60] h Reference [59] should be expected.There is some discrepancy between our calculated E f (V 0 O ) in O-poor conditions and theirs, which, as with the case of ZnO discussed above, is due to the difference in heat of formation used.We have used the experimental value of H (SnO 2 ) = −6.07eV [122], while Scanlon and Watson calculated it to be H (SnO 2 ) = −5.29 eV [20].Comparing with Ágoston et al. [41], who also used PBE0, we find our results computed using the BB1k functional to be closer to their values, but we note that their results also disagree somewhat with Scanlon and Watson [20], despite using very similar techniques.We note that, similarly to the case of ZnO, employing the GGA functional results in lower formation energies [152], but interestingly an earlier calculation using classical forcefields obtained E f (V 0 O ) = 3.65 eV in O-rich conditions [26], a result closer to those obtained using hybrid DFT.Our results disagree with Singh et al. [43], who used a corrected GGA+U (corrGGA+U ) approach and obtained higher formation energies and a deeper transition level, but their results also disagree with hybrid DFT calculations using supercell methods.
Using the BB1k functional results in a slightly higher formation energy and a shallower (but still deep) (2 + /0) level compared with those obtained using the other hybrid functionals.Mizokawa and Nakamura [75] measured conductivities and carrier concentrations as a function of temperature in SnO 2 samples under different external conditions, finding an activation energy of 0.15 eV that increased to 0.3 eV after O 2 adsorption, with an associated reduction in carrier density.One would expect a significant concentration of vacancies on the surface of SnO 2 , where the lower coordination would reduce the electron ionization energy.O 2 adsorption will passivate these vacancies, resulting in an increased contribution from vacancies in the bulk.The 0.15 eV level may then be associated with surface vacancies, and the 0.3 eV level is in good agreement with our result (2/ + 0) = 0.34 eV determined using BB1k.Indeed, in a more recent experimental study [163], where post-treatment resulted in an increase in activation energy and a reduction in intrinsic carrier concentration that was attributed to a reduction in hydrogen donor species, a deep level was determined at 0.28-0.35eV.The authors attributed this level to a deep acceptor, as it was assumed to be associated with N incorporation, but the dominant carrier type was not determined in their study.If, instead, this level can be attributed to V O , which becomes more dominant once the hydrogen donors are removed, then our calculated value is also in good agreement with this result.We note that a level at 0.28-0.35eV has been observed in many n-type samples [75,[164][165][166] (as has the ∼0.15 eV level) [73][74][75]167], and that, in a study on SnO 2 nanoparticles [168], a switching of activation energy from 0.11 to 0.35 eV was observed in the largest nanoparticle when changing from low to high temperature.

I n 2 O 3
In the cases of ZnO and SnO 2 , despite some discrepancies among different computational studies, there appears to be a consensus on the deep and negative U nature of the oxygen vacancy defect levels.For In 2 O 3 , the situation is different, as some studies claim the level is deep, while others claim that it is shallow.From our calculations (see Fig. 1), we find that the level is not of the negative U type, as we determine the + state to be stable for a small range of E F , and is shallow, as the (2+/+) and (+/0) transitions occur above the CBM.Comparing our results with previous calculations (see Table V), we find that the formation energies in the neutral state, similarly to the cases of ZnO and SnO 2 , are slightly higher than those obtained using GGA [56] and corrGGA [42] (apart from the study by Tanaka et al. [152] which resulted in formation energies much lower than any other calculations of which we are aware) and substantially lower than those determined with the corrLDA+U approach [57,58].We find, however, good agreement between our calculations and those obtained using GGA+U [39,56] and with those of Ágoston et al. [41], who used the HSE06 functional.
Here we can only compare the O-poor cases, but since the difference between the O-rich and O-poor formation energies is 1  3 H (In 2 O 3 ), any discrepancy between the computed formation energies at the O-poor limit arising from our use of the experimental heat of formation, H (In 2 O 3 ) = −9.60 eV [122], will be the difference in the computed heat of formation vs the experimental divided by a factor of three; we therefore expect better agreement between our O-poor results and other studies than we would for ZnO or SnO 2 .A positive U nature of the defect transition was also obtained by Ágoston et al. [56] and Liu et al. [56] using GGA+U approaches, in qualitative agreement with our work and in contrast to other studies.
Our finding that the vacancy introduces a resonance in the conduction band, and therefore acts as a shallow donor, is in agreement with Ágoston et al. [41], and at variance with Lany and Zunger and others [39,57,58,60], a point we return to in the next section.We note that our results agree well with those of Walsh and Scanlon [59], who studied the vacancy in the orthorhombic polymorph, o-In 2 O 3 , using hybrid DFT (HSE06) in a supercell approach, and found the level to be shallow.In o-In 2 O 3 , there are two inequivalent sites for oxygen vacancies, hence the two values of formation and ionization energy given in Table V.The shallow nature of the vacancy is supported by experimental studies on the variation in n-type conductivity and carrier concentrations with oxygen partial pressure [169], while the stability of V + O at accessible Fermi energies is consistent with the observed enhancement of ferromagnetism in Fe-and Ni-doped In 2 O 3 that occurs after vacuum annealing [170,171].Finally, our results demonstrate that the vacancy introduces a shallow level in bulk In 2 O 3 , which, if the defect contributes significantly to the observed intrinsic n-type conductivity, would explain that conductivity as a bulk rather than surface effect.

C. Compact versus diffuse defect states
It was pointed out by Lany and Zunger that the wave function associated with V 0 O in In 2 O 3 is well localized at the defect site [60].Indeed, we find well localized, 'compact' wave functions for V 0 O in all three systems, as demonstrated in Fig. 2, where we plot the molecular orbitals associated with the defect states introduced in the gap on formation of a neutral vacancy in each oxide system, as calculated using the BB1k functional (using the other two functionals gives nearly identical results).We disagree, however, with the conclusion that a calculated compact wave function cannot be associated with a shallow level.Compact donors have been described briefly by Stoneham [28], but their relative obscurity in comparison to the well-studied diffuse, effective-mass-like donor appears to be the source of some misinterpretations with regards to characterizing defects in TCOs (and other systems).
In the following, we use the term 'shallow' in the context of a compact donor to signify an occupied (donor) state close to but below the CBM, determined from its low vertical ionization energy [172].
One must consider the balance between different possible configurations of a donor defect in a given charge state.For example, in the case of V 0 O , the state can be compact, with both associated electrons localized within a ∼ 2 Å diameter spheroid, similar to those shown in Fig. 2, or can consist of an overall charge-neutral complex, combining a compact V + O with a weakly bound, diffuse, effective-mass-like electron, e − m * : , where both configurations represent distinct states of the defect [173].In Fig. 3 we depict the difference between V 0 O and V + O + e − m * but of course the argument could also be extended to include the complex V 2+ O + 2e − m * .Which configuration is most favorable depends on the relative energies.Here, as a first approximation, assuming isotropic and parabolic bands we can model e − m * using effective mass theory [175], i.e., the energy (in atomic units) of the electron is E CBM − m * /2ε 2 , where E CBM is the energy of the CBM, m * is the electron effective mass, and ε is the static dielectric constant.The balance to consider is then  within the embedded cluster approach, while the energies of the e − m * bound to make the defect overall charge neutral, can be added a posteriori.Thus, there is no contradiction when stating that, depending on the energetics, a compact donor may be shallow relative to the appropriate band edge [176].We note that this analysis can also be applied in supercell-based techniques.
In Table VI we give the calculated energies of V 0 O , V + O + e − m * , and V 2+ O + 2e − m * , determined at E F = E CBM for each oxide and density functional.We take computed values of ε from our interatomic forcefields, which are in good agreement with experiment, and experimental values of m * (in units of the electron mass m e ), averaging the diagonal tensor components in the noncubic cases of ZnO and SnO 2 : ε = 9.77 [117,177], m * = 0.22 m e [178,179] for ZnO; ε = 12.66 [104,180], m * = 0.26 m e [181] for SnO 2 ; ε = 9.05 [25,182], m * = 0.22 m e [68] for In 2 O 3 .The resulting binding energy for a single diffuse electron, taken relative to the CBM, is 0.032 eV for ZnO, 0.022 eV for SnO 2 , and 0.036 eV for In 2 O 3 , while the Bohr radii of the e − m * are 23 Å, 26 Å, and 22 Å, respectively.The results in Table VI differ from those given in Fig. 1 TABLE VI.Energy balance between the neutral vacancy modeled as a compact state with two localized, trapped electrons (V 0 O ), a compact state with one localized and one effective-mass-like electron (V + O + e − m * ), and a compact state with two effective-mass-like electrons (V 2+ O + 2e − m * ).The energies (in eV) are given for the Fermi energy at the CBM and are shown for each oxide and each density functional, in O-rich conditions.only in the added energies of the e − m * ; in the approach taken to determine the energies shown in Fig. 1 those electrons are assumed to be fully ionized.For both ZnO and SnO 2 , the neutral vacancy binding two localized electrons is most stable at the CBM, indicating that the defect is not a shallow donor.In contrast, we find for In 2 O 3 that the defects binding diffuse electrons are lower in energy than that binding localized electrons at the CBM.This conclusion matches that which we made when discussing the results presented in Fig. 1; here we are highlighting that the modeling of each charge state as compact is consistent with the possibility of the defect being a shallow donor.

D. Charge carrier and defect concentrations
We now turn to our calculated defect and excess carrier concentrations as a function of temperature, assuming thermodynamic equilibrium, derived from our computed formation energies and densities of electronic states.We consider the range of temperature 0 < T 1500 K, which encompasses all common synthesis temperatures of the three oxides and a majority of derived device operational temperatures.We continue to employ the three density functionals in this section in order to compare their effects and demonstrate that the main conclusions drawn from these results are largely independent of which hybrid functional we choose.The results are shown in Figs. 4, 5, and 6.We first discuss SnO 2 .
As shown above, we calculate (2 + /0) to lie substantially below the conduction band in SnO 2 .As a consequence we find that, in O-rich conditions where the formation energies are of the order of 4 eV, E F remains deep in the gap (between 1.8 and 1.4 eV below the CBM, see the insets on the right-hand side of Fig. 4).The carrier concentrations remain less than 10 15 cm −3 for T 1500 K, and with the PBE0 and B97-2 functionals we find that the system is semi-insulating, with the electron concentration n 0 equal to that of the hole concentration p 0 , and the oxygen vacancy concentration [V O ] an order of magnitude below both.With BB1k, where we calculated the transition level to be closer to the CBM, the vacancy formation energy at E F is lower than in the cases of the other functionals and [V O ] is correspondingly higher.Consequently, we find that there is an excess of electron carriers, at a concentration double that of [V O ] and over an order of magnitude greater than p 0 , but remaining below 10 15 cm −3 up to T = 1500 K.
In O-poor conditions, the formation energies are substantially lower, and we therefore find that E F is shifted closer to the CBM, varying from 1 to 0.5 eV below it as T is increased (see the insets on the left-hand side of Fig. 4), and that oxygen vacancies are found in significant concentrations, rising above 10 19 cm −3 for T > 1000 K (using the PBE0 and B97-2 functionals).With the PBE0 and B97-2 functionals, in this range of E F the dominant charge state for the V O is neutral, but at T > 1000 K we find that approximately 2% of the vacancies are thermally ionized, resulting in electron concentrations ranging up to 10 18 cm −3 at T = 1500 K.With the BB1k functional, as the transition level is shallower, a greater proportion of the vacancies are thermally ionized (at T = 1200 K, 27% are in the + and 65% are in the 2+ state), and we find that n 0 and [V O ] are close in value and increase to just above 10 18 cm −3 at T = 1500 K. From these results, we see that, despite such a deep transition level, in reducing conditions significant intrinsic charge carrier concentrations can arise at temperatures above 800 K, in good agreement with experiment [72][73][74].For T < 800 K, the concentrations we calculate are somewhat lower than experiment; the source of this discrepancy may be residual n-type activity from impurities or other defects which could dominate at reduced temperatures.
For ZnO, we calculated a similar (2 + /0) to that of SnO 2 , as can be seen in Fig. 1, although the results using PBE0 and B97-2 coincide more closely, and the formation energy of the neutral vacancy using BB1k is closer to those of the other functionals.Under O-rich conditions, the calculated concentrations shown in Fig. 5 are quite similar to those of SnO 2 .E F remains deep in the gap, more than 1 eV below the CBM in all cases, while, for the PBE0 and B97-2 functionals, the carrier concentrations remain below 10 15 cm −3 for T 1500 K and we have n 0 = p 0 with [V O ] < 10 13 cm −3 .Similarly to the case of SnO 2 , when using the BB1k functional, we calculate [V O ] to be higher than that computed using the other functionals as the transition level is shallower, meaning that the formation energy is lower at the self-consistent E F .This concentration of vacancies in the 2+ charge state results in an excess of electrons (which remains below 10 16 cm −3 for T 1500 K), with n 0 approximately double [V O ] and with a minority carrier concentration two orders of magnitude lower.
In O-poor conditions, however, the results differ from those of SnO 2 .As discussed in Sec.III B 1, due to our slightly lower formation energies and usage of the experimental value of H (ZnO), our calculated E f (V O ) in O-poor conditions are significantly lower than those of other DFT studies, at about 0.05 eV for PBE0 and B97-2 and 0.19 eV for BB1k (see Table III).Such a low formation energy would indicate high nonstoichiometry in extreme reducing conditions, for which there is some experimental evidence [183][184][185], but measurements to investigate the intrinsic electrical conductivity properties of ZnO tend to be performed under less extreme conditions [36,88,186].From such a low formation energy, we determine very high [V O ], of the order of 10 22 cm −3 for T > 400 K, using the PBE0 and B97-2 functionals.At such high concentrations, the assumption of noninteracting defects in the dilute limit breaks down, and we must treat our results with caution.For the BB1k functional, the formation energy is not quite so low, and we get vacancy concentrations an order of magnitude lower.We also find that n 0 rises above 10 19 cm −3 at T = 800 K, an intrinsic concentration that is about an order of magnitude greater than that determined by Halliburton et al. [78], who measured samples treated in Zn-rich conditions at elevated temperatures (with PBE0 and B97-2 n 0 = 10 17 cm −3 at T > 900 K and approaches 5 × 10 18 cm −3 at T = 1500 K, in good agreement with Ref. [78]).
If we take into account a rough estimate of the difference between the calculated enthalpy of formation and experiment (about 0.2 eV for the HSE06 functional) [55] and shift our formation energies correspondingly, the concentrations reduce by close to two orders of magnitude.Then, our results become qualitatively similar to SnO 2 , with intrinsic electron carrier concentrations of 10 16 -10 17 cm −3 for T > 1000 K. Regardless of which of the two values of oxygen chemical potential we take to represent O-poor conditions, we therefore conclude that ZnO will have high concentrations of vacancies in extreme reducing conditions and will have intrinsic carrier conditions of the order of 10 16 -10 19 cm −3 , despite having such a deep transition level, results that are consistent with experiment [36,78,183,184,186].
The transition level is significantly shallower in In 2 O 3 , and there is a small range of E F where V + O are stable, resulting in concentration values and trends with T that are somewhat different from those in the other oxides studied [187].In O-poor conditions, we find that E F varies between 0 and 0.5 eV above the CBM over the T range studied when using the PBE0 and B97-2 functionals (the BB1k functional results in the range 0.4 to 0.7 eV above the CBM).For PBE0 and B97-2, n 0 rises sharply at low T to above 10 18 cm −3 (at T = 100 K for B97-2 and T = 250 K for PBE0), then continues to rise but at a lower rate as ρ(E) above the CBM changes less rapidly as E F increases.At about T = 600 K, the charge state V + O becomes populated as well as V 2+ O , which leads to a greater increase in n 0 with T .That increase then slows down at higher T , where E F and T are high enough that a significant concentration of V 0 O is present.With BB1k, the dominant charge state is V 2+  O throughout the range of T studied, and we have n 0 of the order of 10 20 cm −3 , with [V 2+  O ] = n 0 /2.Such concentrations have been observed experimentally, typically when the samples are annealed in vacuum or N 2 or in low P O 2 environments [32,169,[188][189][190][191], but we note that in some cases lower concentrations, of the order of 10 18 -10 19 cm −3 , have been measured [32,62,63,191].Of course, we have not considered other intrinsic defects, such as oxygen interstitials or cation vacancies, which may be present and act as electron killers, resulting in a lower n 0 .
In O-rich conditions, the formation energies are over 2 eV higher than those in O-poor conditions, resulting in lower vacancy concentrations and an E F pushed into the band gap (somewhere between 0.5 and 1.2 eV below the CBM, depending on the functional used and T , see Fig. 6).We still find, however, that In 2 O 3 will be n-type, as [V 2+  O ] results in an n 0 that increases close to 10 17 cm −3 at higher T when using the PBE0 and B97-2 functionals (with the BB1k functional, as the formation energy is lower than the other functionals, n 0 is about an order of magnitude higher).When annealed in O 2 or in air, In 2 O 3 samples have displayed intrinsic carrier concentrations of this order, which persist even at lower T [29,32,62,63,190,191].Our results support the assignment of V O as the source of these carriers: Vacancies can be formed at elevated T , and a significant proportion can remain in the system at lower T due to kinetic barriers.We therefore conclude that, even in extreme O-rich conditions, due to the presence of V O In 2 O 3 will be intrinsically n-type, as is seen experimentally.
These results demonstrate the importance of taking into account the self-consistent Fermi energy when analyzing thermal transition levels and their effect on carrier concentrations.We have shown that, even in SnO 2 and ZnO, where the oxygen vacancy transition levels are deep in the gap, under certain conditions the materials can by intrinsically n-type, simply due to the presence of oxygen vacancies.

IV. CONCLUSIONS
We have found that oxygen vacancies form shallow, compact donors in bulk In 2 O 3 and are crucial for its intrinsic n-type conductivity.Oxygen vacancies in SnO 2 and ZnO form deep states, but nevertheless play an important role in the promotion of n-type carriers under certain conditions, a result that is evident from the analysis of the self-consistent Fermi energy and equilibrium defect and carrier concentrations.We have clarified that compact states can form shallow levels by studying the energy balance between different configurations of the vacancy donor in the three TCOs, finding that the vacancy binding effective-mass-like electrons, rather than localized electrons, is energetically favorable in In 2 O 3 , which accounts for the shallow nature of the defect state.We determine formation and ionization energies with the hybrid QM/MM embedded cluster approach that agree well with those obtained using a plane-wave supercell method and are in excellent agreement with such calculations when similar hybrid density functionals are employed to treat exchange and correlation.Our results agree well with experiment and demonstrate the key role played by oxygen vacancies in the conductivity properties of these TCOs.

FIG. 1 .
FIG. 1. Formation energies as a function of Fermi energy relative to the valence band maximum (VBM) of oxygen vacancies in (a) ZnO, (b) SnO 2 , and (c) In 2 O 3 , shown for O-rich and O-poor conditions.Results are given for three hybrid density functionals: BB1k (black), B97-2 (red), and PBE0 (blue).A vertical dashed line is placed at the position of the conduction band minimum for each system (derived from the experimental band gap; see the main text for details).The slopes indicate the charge states; a transition occurs where the slope of a line changes.For ZnO, the B97-2 result lies almost entirely below the PBE0 result.

FIG. 2 .
FIG. 2. The calculated molecular orbital associated with the defect level introduced by a neutral oxygen vacancy in (a) ZnO, (b) SnO 2 , and (c) In 2 O 3 , determined using the BB1k functional.Oxygen ions are represented by red, In by brown, Sn by darker gray, and Zn by lighter gray spheres.The molecular orbitals are indicated by yellow (positive) and purple (negative) isosurfaces of densities (in atomic units) 0.1, 0.08 and 0.06 (in the cases of In 2 O 3 and ZnO, the highest density isosurfaces are visible on the ions surrounding the vacancy).

FIG. 3 .
FIG. 3. Cartoon depicting a 'compact' neutral oxygen vacancy donor, with two electrons trapped in the vacancy, and a 'diffuse' donor, where one of the bound electrons is of the effective-mass type.r indicates the approximate extent of the wave function ψ, which is represented by the blue lines.The black line represents the shape of the electrostatic potential.

2 FIG. 4 .
FIG. 4. Calculated defect ([V O ], green line) and electron (n 0 , red line) and hole (p 0 , blue line) carrier concentrations as a function of temperature T for SnO 2 , determined using the PBE0, B97-2, and BB1k hybid density functionals, under O-poor and O-rich conditions.The insets show the self-consistent Fermi energy E F (black line) as a function of T , relative to the conduction band minimum.

FIG. 5 .
FIG. 5. Calculated defect ([V O ], green line) and electron (n 0 , red line) and hole (p 0 , blue line) carrier concentrations as a function of temperature T for ZnO, determined using the PBE0, B97-2, and BB1k hybid density functionals, under O-poor and O-rich conditions.The insets show the self-consistent Fermi energy E F (black line) as a function of T , relative to the conduction band minimum.

3 FIG. 6 .
FIG. 6. Calculated defect ([V O ], green line) and electron (n 0 , red line) and hole (p 0 , blue line) carrier concentrations as a function of temperature T for In 2 O 3 , determined using the PBE0, B97-2, and BB1k hybid density functionals, under O-poor and O-rich conditions.The insets show the self-consistent Fermi energy E F (black line) as a function of T , relative to the conduction band minimum.

TABLE I .
Embedding cation ECP parameters in atomic units [see Eq. (1) in the main text].

TABLE II .
Calculated ionization potentials (eV), shown in comparison with previous calculations and experiment.

TABLE III .
Calculated formation energies of neutral oxygen vacancies E f (V 0 O ) in O-rich and O-poor conditions and the (2+/0) thermal transition level [ (2 + /0), referred to the conduction band minimum] in ZnO from this study using the PBE0, B97-2, and BB1k hybrid functionals, compared with previous calculations (Prev.Calc.) using a variety of functionals (see main text for details).Energies are given in eV.

TABLE IV .
Calculated formation energies of neutral oxygen vacancies E f (V 0 O ) in O-rich and O-poor conditions and the (2+/0) thermal transition level [ (2 + /0), referred to the conduction band minimum] in SnO 2 from this study using the PBE0, B97-2, and BB1k hybrid functionals, compared with previous calculations (Prev.Calc.) using a variety of functionals (see text for details).Energies are given in eV.