First-Principles Modeling of Polaron Formation in TiO 2 Polymorphs

: We present a computationally e ﬃ cient and predictive methodology for modeling the formation and properties of electron and hole polarons in solids. Through a nonempirical and self-consistent optimization of the fraction of Hartree − Fock exchange ( α ) in a hybrid functional, we ensure the generalized Koopmans ’ condition is satis ﬁ ed and self-interaction error is minimized. The approach is applied to model polaron formation in known stable and metastable phases of TiO 2 including anatase, rutile, brookite, TiO 2 (H), TiO 2 (R), and TiO 2 (B). Electron polarons are predicted to form in rutile, TiO 2 (H), and TiO 2 (R) (with trapping energies ranging from − 0.02 eV to − 0.35 eV). In rutile the electron localizes on a single Ti ion, whereas in TiO 2 (H) and TiO 2 (R) the electron is distributed across two neighboring Ti sites. Hole polarons are predicted to form in anatase, brookite, TiO 2 (H), TiO 2 (R), and TiO 2 (B) (with trapping energies ranging from − 0.16 eV to − 0.52 eV). In anatase, brookite, and TiO 2 (B) holes localize on a single O ion, whereas in TiO 2 (H) and TiO 2 (R) holes can also be distributed across two O sites. We ﬁ nd that the optimized α has a degree of transferability across the phases, with α = 0.115 describing all phases well. We also note the approach yields accurate band gaps, with anatase, rutile, and brookite within six percent of experimental values. We conclude our study with a comparison of the alignment of polaron charge transition levels


INTRODUCTION
Charge trapping in semiconductors and insulators has widespread interest in the fields of physics, chemistry, and materials science and controls performance for important applications, such as photocatalysis, superconductivity, tribocharging, magnetism, and microelectronics. 1−7 The trapping of electrons or holes may take place at pre-existing defects (such as vacancies or impurities) or even in the perfect lattice in some materials (socalled polaronic self-trapping). 8−13 Predictively modeling these effects is vital in order to provide a deeper theoretical understanding and to guide materials optimization for applications. Variants of Kohn−Sham density functional theory (DFT) 14,15 are often used to model charge trapping. 9,13 Although DFT is exact in principle, practical calculations require an approximation to the exchange-correlation (xc) potential. The self-interaction error (SIE) which plagues many xc-approximations, such as the local density approximation 16 or (semilocal) generalized gradient approximation, 17,18 results in the spurious delocalization of localized charges. 19−24 Other methods, such as DFT+U 25 and hybrid DFT, 26 aim to correct this behavior; however, they introduce additional parameters that are often fit to experimental data 27 limiting its predictive power. More sophisticated many-body methods that are in principle predictive, such as Møller−Plesset perturbation theory or the GW approximation, are extremely computationally demanding, especially for modeling polaronic charge trapping where full self-consistency for electrons and ions is essential. For these reasons, finding an inexpensive, predictive, and transferable method to model charge trapping in materials is urgently needed.
In this article, we use an inexpensive hybrid functional approach 28 utilizing the auxiliary density matrix method (ADMM) 29 to predictively model charge trapping in titanium dioxide (TiO 2 ). 30 The generalized Koopmans' condition−a known property of the exact functional−is used to parametrize the proportion of Hartree−Fock exchange in the hybrid functional (α), minimizing the effects of the SIE. 31−34 This leaves no parameters to be determined empirically allowing for a predictive approach. As a strongly correlated system capable of forming intrinsic polaronic defects, TiO 2 is an excellent material for testing this method. To illustrate the predictive capabilities of this approach, we investigate the possibility of polaron formation in a number of TiO 2 phases yet to be thoroughly studied in the literature. Alongside the well-known rutile and anatase phases, we consider the naturally occurring brookite and synthetically created TiO 2 (H), 35 TiO 2 (R), 36 and TiO 2 (B) 37 phases, all of which are phases that are stable or metastable in ambient conditions. Our findings suggest that SIE can be removed across all TiO 2 phases with a universal hybrid functional para-metrization. Furthermore, the use of Koopmans' condition yields predicted band gaps in excellent agreement with experimental results. We find that small (Holstein) polaron formation occurs in all phases of TiO 2 considered. The results for anatase and rutile are consistent with previous theoretical work and experimental evidence. The formation of multisite polarons is predicted in TiO 2 (H) and TiO 2 (R) due to the unique bonding present between anions and cations in the phases. These results suggest that hybrid functionals parametrized using Koopmans' condition enable predictions to be made with regards to intrinsic charge trapping. Furthermore, due to efficient evaluation of Hartree−Fock exchange terms permitted by the ADMM method the computational cost is only two to three times the cost of a standard DFT calculation. Therefore, the approach should find applications to model charge trapping at complex defects such as surfaces and interfaces relevant for technological applications but previously inaccessible to predictive modeling.
The structure of this article is as follows. We first briefly review previous experimental and theoretical work on polaronic charge trapping in TiO 2 and related materials in Section 2. Our methodology together with technical details of the hybrid functional implementation are described in Section 3. We then present our predictions for polaron formation in all TiO 2 phases in Section 4 before discussing the factors which may influence their accuracy and concluding.

BACKGROUND
It is known that the exact DFT functional should give a total energy that is piecewise linear with respect to the fractional number of electrons in a system 31 and has a slope corresponding to the frontier orbital eigenvalue ( Figure 1). 38,39 Approximations to xc such as the LDA and GGA display convex behavior due to residual SIE and give rise to spurious delocalization of charge. 22−24 The Hartree−Fock (HF) approximation on the other hand exhibits concave behavior and tends to overlocalize charge. One rationale behind hybrid functionals is that by mixing both of these opposing errors together, a near piecewise linear functional may be obtained. Similarly, the DFT+U approach introduces a correction that counters the convex error present with (semi)local xc-approximations. 40,41 Therefore, in principle, with the correct parametrization, both hybrid and DFT+U approaches can result in a functional which may be close to piecewise linear and therefore exhibit much reduced SIE.
The generalized Koopmans condition (GKC), 31−34 which must be satisfied for an exact Kohn−Sham or generalized Kohn−Sham functional, provides a convenient way to formulate the requirement of piecewise linearity where I(N) and E(N) are the total energy and ionization energy of an N electron system, and ε i (N) denotes the i th eigenvalue of an N electron system. We note that this condition only makes reference to the highest occupied molecular orbital (HOMO) and so only addresses piecewise linearity between N and N+1 electrons. Strictly speaking, this is a necessary but not sufficient condition; however, previous work has shown linearity is greatly improved when the condition is satisfied. 42−44 While application of the GKC to localized states (such as polarons) is straightforward using the above equation, its application to delocalized states has had less success. 45 A recent study has shown that Koopmans-compliant hybrid functionals can yield accurate electron densities in comparison to exact solutions for simple model systems. 46 An alternative approach to the GKC, which draws insights from the many-body perturbation theory, can also be used to parametrize hybrid functionals. This approach constrains hybrids to obey the condition α = ε ∞ −1 , where ε ∞ is the electronic contribution to the dielectric constant, which in turn ensures that hybrids have the correct form for the long-range tail of the Coulomb interaction. 47 If one is primarily concerned with removing localization errors associated with polarons, we believe it is more appropriate to ensure the GKC for the polaronic configuration is satisfied rather than any condition applied to properties of the bulk crystal.
DFT+U has been applied successfully to model polarons in a range of metal oxides by introducing a Hubbard-U correction on cations and, optionally, anions. 48, 49 Both the GKC and comparison to experimental data have been used to determine these Hubbard-U values. 48,49 One shortcoming of these approaches is, although a qualitatively correct description of the polarons is obtained, quantitative results are often incorrect. For example, the anatase band gap given by DFT+U (for reasonable parametrizations) is significantly underestimated from the experimental band gap of 3.2 eV. 50 A similar situation is found in other crystals, such as zinc oxide, ceria, and rutile-TiO 2 . 50 This underestimation naturally yields incorrect polaron charge transition levels (CTLs), hindering comparison to experimental data. The DFT+U approach also predicts a low (total) dielectric constant in titania, 51,52 meaning electron and ionic screening properties are underestimated. As finite size corrections are inversely proportional to the total dielectric constant ε, 53 an underestimated ε means that they become more of an issue.
Unlike DFT+U, hybrid functionals such as PBE0 or HSE06 generally yield well-behaved crystal properties. These functionals usually have at least one parameter, the fraction of exact exchange α, which determines the proportion of HF exchange used. Some hybrid functionals make a distinction between shorta n dl o n g -r a n g ee x c h a n g e( e . g . ,H S E 0 6 ) ,r e s u l t i n gi na n additional range-separation parameter ω. The parametrization of hybrid functionals usually comes in two varieties. A "one-sizefits-all" approach, where parameters are fixed by comparison to a large set of training data and then used for a variety of different systems. The performance of such hybrids usually varies

Journal of Chemical Theory and Computation
Article depending on the choice of system. An alternative approach is to fit hybrids for a specific system using, for example, the GKC or by fitting to experimental data. 42,44,54 Fitting parameters to experimental properties such as the band gap may give reasonable results; however, the approach is somewhat semiempirical. Furthermore, the predictive power of such an experimental fitting can also be questioned. In contrast, methods using the GKC to parametrize the hybrid functional have been shown to yield extremely accurate ionization and band gap energies. 42 However, the degree of transferability for these parameters remains somewhat unknown. A common issue regarding the use of hybrid functionals is the computational cost associated with them, usually around ten times the cost of a standard DFT calculation. To compensate for this, one is forced into using smaller supercells in periodic calculations. This can hinder the quality of a calculation as artificial interactions between periodically repeated defects tend to grow larger with smaller system sizes. 53 Therefore, while hybrid functionals can be superior to DFT+U, they come at a higher computational cost.
We now discuss some of the previous experimental and theoretical results for charge trapping in TiO 2 . The stability of electron polarons in rutile and anatase is a widely discussed question, with many differing reports from both experimental and theoretical work. While we only consider self-trapped electron and hole polarons in this work, numerous studies highlight the effects of reduced cation and anion species due to dopants, such as Nb, Li, and H, and vacancy defects. 55−61 Electron paramagnetic resonance (EPR), 62,63 infrared spectroscopy (IR), 64 and scanning tunneling microscopy (STM) 58 data suggest that small (Holstein) electron polarons can form in rutile localizing on titanium 3d orbitals. IR studies suggest that the polaron CTL is shallow, lying roughly 0.3 eV below the conduction band minimum (CBM). 64 Other data from deep level transient spectroscopy, STM, and optical absorption places the polaron CTL 0.5 eV, 65 0.7 eV, 58 and 0.8 eV 59 below the CBM, giving the CTL a range of 0.3−0.8 eV below the CBM. We note that oxygen vacancies give rise to defects known to lie across a similar, albeit slightly deeper, range of CTLs below the CBM. The experimental thermal activation energy associated with polarons in rutile is around 20−30 meV. 62 In rutile, the HSE06 hybrid functional has been shown to satisfy the GKC to within 0.1 eV, making it a suitable choice. 9,66,67 It predicts the electron polaron CTL is around 0.5 eV below the CBM, in good agreement with the 0.3−0.8 eV range, and a band gap of 3.4 eV, only a slight overestimate with respect to the experimental value of 3.0 eV. An alternative parametrization of the HSE functional (with a smaller range separation parameter, ω = 0.1 Å −1 ) predicts the polaron CTL is 0.8 eV below the CBM; 68 however, the linearity of this functional was not assessed. Another hybrid functional, B3LYP, significantly overestimates the band gap of rutile 69 resulting in much deeper defect levels than HSE06 and is known to not satisfy the GKC. 70 The polaron trapping energy E t , the difference between localized and delocalized electron solutions, provides a measure of polaron stability. This energy difference is predicted to be small with HSE06 and DFT+U (<0.1 eV), 58 and substantially larger with the RPA (0.6 eV), 71 although we note the RPA was done on a PBE0 geometric configuration. Polaron migration barriers calculated using HSE (with ω = 0.1 Å −1 ), the RPA, and DFT+U (with U = 10.0 eV) yield 0.03, 0.1, and 0.3 eV.
In defect-free anatase, electronic states are likely more delocalized than that of rutile as suggested by EPR data. 63,72 This information leads us to believe that, unlike rutile, excess electrons in bulk anatase do not form self-trapped (small) electron polarons but are much more diffuse. Further confirmation is obtained from STM experiments which find similar results in Nb doped samples. 58 This is also highlighted by conductivity measurements, which find Nb doped samples of anatase exhibit 3−4 orders of magnitude higher conductivity compared to rutile. 73,74 HSE06, 9 the random phase approximation (RPA), 71 and DFT+U (with U parametrized through the RPA) 58 all indicate that small electron polarons are unstable in anatase, further indicating why anatase has a higher electron mobility than rutile. On the other hand, B3LYP predicts that electron polarons in anatase are stable. 69,75−78 RPA results indicate that a large (Froḧlich) polaron forms when using large system sizes; 71 however, experiments indicate that this only occurs in the vicinity of a Nb defect. 58 The behavior of hole polarons in rutile and anatase is the opposite of the electron polaron case. Specifically, hole polarons can form in anatase but not rutile. HSE06 results predict that excess holes in anatase form self-trapped hole polarons (on oxygen 2p orbitals); however, hole polarons are not predicted to form in rutile (i.e., they are delocalized). HSE06 predictions suggest that hole polaron levels in anatase are a deep acceptor state, lying roughly 1.3 eV above the valence band minimum (VBM) with a trapping energy of 0.2 eV. EPR experiments suggest trapped holes can form in anatase. 79 The concentrations of trapped holes (O − species) were higher than those of trapped electrons (Ti 3+ species), primarily because most electrons were delocalized in the conduction band and hence EPR silent. In contrast, rutile samples showed larger EPR signals attributed to electron trapping than hole trapping, 80 in agreement with HSE06 results. Other EPR experiments somewhat confuse the picture as it is suggested that hole trapping in rutile is observed. 81 There have been very few studies of polaron formation in phases of TiO 2 other than rutile and anatase. This includes an experimental investigation into polaron self-trapping in brookite, 63 which finds that anatase and brookite have similar trapping properties. There has also been a theoretical study of polaron formation in using DFT+U which finds charge strongly localized on Ti sites in Li-doped TiO 2 (B). 82

METHODS
We use the hybrid DFT implementation within CP2K to perform all of our calculations. 28,29,83−85 Our Gaussian type orbitals are mapped onto CP2K's multigrid solver with a relative cutoff of 60 Ry. We use five multigrids, with the finest grid having a cutoff of 600 Ry. The primary oxygen and titanium basis sets we use are of double-and triple-ζ quality with valence and polarization exponents. 86 We use the Goedecker-Teter-Hutter (GTH) pseudopotentials available with CP2K. 87 As CP2K only samples the Γ-point in reciprocal space, we are required to converge crystal properties with respect to system size as opposed to the number of k-points. This is illustrated for the rutile phase in Supporting Information (SI) Figure S1 and Tables S1 and S2, where we find supercell dimensions of around 18 Å are well converged. We note that using larger system sizes results in band folding, which increases the number of bands sampled at each k-point in the Brillouin zone. This is equivalent to increasing the number of k-points sampled, meaning that the calculation remains accurate despite only sampling the Γ-point. Furthermore, using large system sizes instead of increasing the number of k-points reduces the size of the interaction between charged supercell images.

Article
The exchange part of the hybrid functional we employ (denoted PBE0-TC-LRC) has the following form  (2) where r 12 is the exchange integral distance, R c is the truncation (abbreviated to TC) parameter, E x HF,TC is truncated exact HF exchange, E x PBE,LRC is a long-range correction (abbreviated to LRC) for exchange, and E x PBE is PBE exchange. 28 In both the short-range (r 12 ≤ R c ) and long-range (r 12 > R c ) regions the full PBE correlation energy is used. We note that the value R c must be smaller than half the smallest supercell lattice vector to prevent exchange interactions between electrons in neighboring periodic images. This truncation is performed to reduce computational costs associated with HF exchange. However, as discussed in detail below we find predicted properties (e.g., lattice constants, band gaps, polaron trapping energies, and charge transition levels) converge rapidly with respect to increasing R c . For example, Table S3 (SI) shows that the lattice constants and band gap for rutile TiO 2 are well converged for R c = 6.0 Å.
In addition to HF exchange truncation in the functional, CP2K has other tools to alleviate computational cost. One such method is the auxiliary density matrix method (ADMM), 29 which approximates exchange integrals by mapping orbitals onto smaller, more localized, basis sets. These significantly improve the speed of HF exchange calculations, meaning the approach can be used to model larger systems with reasonable computational cost. In later sections we discuss variations in choices of ADMM basis sets to establish any error this may introduce in our calculations. CP2K also takes advantage of other exchange integral screening methods which reduce computational cost, all of which are shown in the example input file provided in the SI.
The procedure for determining the proportion of HF exchange for electron or hole polarons is as follows. We initially optimize cell vectors for a large defect-free supercell. As hybrid functional CP2K calculations use the Γ-point only, this must be done using the supercell rather than the unit cell. Following this a distortion is applied to the geometry, pushing ions within 2.2 Å of a chosen ion radially outward by 0.1 Å. This distortion creates a precursor potential well for localization of the polaron on a given ion. An additional electron (or hole) is then added to the supercell, and atoms are relaxed until the magnitude of forces on atoms is less than 0.01 eV/Å. The final calculation (vertically) removes the added electron (or hole) from the system, keeping ions in the polaronic configuration. The above calculations yield all the necessary energies to test the GKC (eq 1). We identify, with half-percent granularity, the range of α values which satisfy the condition to within 0.05 eV. We denote the degree to which the GKC is satisfied, also called the nonlinearity, as ξ: We set a limit of |ξ| ≤ 0.05 eV as physical properties of the crystal can change quickly with small changes in nonlinearity. We find that operating within this limit causes changes in the GKS gap of around 0.2 eV over the range −0.05 ≤ ξ ≤ 0.05 eV, as shown in Tables S6 and S7 (SI).

RESULTS
4.1. Anatase and Rutile. In this section we illustrate the application of the approach described above to model electron and hole polarons in anatase and rutile as test cases. The structures of the two crystals are shown in Figure 2. Initially we use the titanium FIT9 and oxygen pFIT3 (FIT9/pFIT3) ADMM basis sets (see the SI), although we test alternatives below to demonstrate they are of reasonable quality. We use 5 × 5 × 2 and 4 × 4 × 6 expansions of the anatase and rutile unit cells, yielding 600 and 576 atoms in the supercell. As mentioned previously, we demonstrate these supercell sizes are well converged in the SI.
In Figure 3 we tune α to satisfy the GKC for the electron and hole polarons in rutile and anatase. It can be seen that the anatase electron polaron case does not satisfy the GKC in the range 0.25 ≤ α ≤ 0.40. Attempting to reduce α to a lower value (α = 0.225) results in a delocalized conduction band state, suggesting that the electron polaron is not stable at the Koopmans-compliant value of α. Probing the point where the electron delocalizes closely, we find ξ = 0.29 eV at α = 0.230. A small move to α = 0.225 results in a completely delocalized electron, meaning that this method of constraining α predicts that even a metastable electron polaron configuration is not possible. We note that fully self-consistent optimization of the geometry for a given α is essential for accurate prediction of polaron properties. For example, if one fixes the polaron geometry to that obtained with α = 0.25, the electron polaron remains localized even when α is reduced to 0.10, contrary to the behavior seen when the geometry is self-consistently optimized.
Unlike the electron polaron case, the anatase hole polaron does satisfy the GKC at α = 0.105 with ξ ≈ 0 eV. The hole polaron spin density is illustrated in Figure 4, where it can be seen that the hole polaron is largely localized on a single oxygen ion. We define the trapping energy E t as the energy difference between localized and delocalized states, which for the hole polaron is −0.21 eV. While we would not advocate use of a non-

Journal of Chemical Theory and Computation
Article self-consistent geometry for prediction of polaron properties, we note that for deep traps, such as the hole polaron in anatase, there is little difference in polaron geometry between α = 0.25 and α = 0.105.
In rutile, we find electron polarons satisfy the GKC, while hole polarons do not, as shown in Figure 3. As mentioned in Section 2, this result is also in agreement with findings using HSE06 by Deaḱ et al. 9 However, unlike the anatase case, we find that a large number of α values satisfy the GKC. As the value of α decreases across the range where −0.05 ≤ ξ ≤ 0.05 eV, the electron polaron trapping energy increases from a value suggesting the polaron is mildly stable (E t = −0.07 eV at α = 0.14) to metastable (E t = 0.01 eV at α = 0.10). The localized solution cannot be found at α = 0.095. Ultimately, as the polaronic state is capable of satisfying the GKC we are inclined to believe that the polaron is stable. A low trapping energy suggests there is little difference energetically between localized and delocalized states, in agreement with HSE06 and DFT+U results. 9,58 The small trapping energy may explain why these polarons rapidly delocalize in EPR experiments when temperature is raised. 80 The spin density for the electron polaron is illustrated in Figure  4, where it is clear that, unlike the hole polaron in anatase, the electron polaron is much more diffuse.
In the region 0.100 < α < 0.140 the GKC is satisfied to less than 0.05 eV for the rutile electron polaron. Across this range, the polaron CTL moves from 0.22 to 0.65 eV below the CBM. We select the largest value of α satisfying |ξ| < 0.01 eV (α = 0.115) as the optimal value for evaluation of polaron properties. In this case the electron polaron CTL is predicted to lie 0.40 eV below the CBM. Experimental work puts the CTL at 0.3−0.8 eV below CBM as previously mentioned, which is in good agreement with our result. The prediction is also in good agreement with previous calculations performed using the HSE06 functional which places the polaron CTL at 0.5 eV below the CBM. 9 We draw attention to our previous discussion (in Section 2) that highlighted that TiO 2 band gaps from B3LYP and DFT+U are generally over-and underestimated, by around 1 eV. 50,69 In contrast, the HSE06 functional overestimates the gap only slightly by 0.2 eV. 9 Surprisingly, our procedure yields defect-free band gaps in both rutile and anatase that deviate only by around 6% from the experimental values of 3.0 and 3.2 eV (with α = 0.115 and α = 0.105) respectively. In both cases, the gap is underestimated by around 0.2 eV, as seen in Table 1. We stress that our results have not been optimized for the band gap. However, it seems that applying the GKC to these localized states results in a band gap very close to the experimental value. If we take the upper limit of our tolerance of |ξ|,wefind that the errors in the band gaps of rutile and anatase shrink further. We find that our lattice constants are also in good agreement with experiment, with errors of 0.3% (vector a) and 1.1% (vector c) in anatase and 0.7% (vector a) and 0.1% (vector c) in rutile. This is somewhat reassuring, as it suggests that the method is able to accurately describe both the polaronic and bulk properties of these crystals using only the GKC. This is also in line with results from finite systems, which report similar findings regarding gaps  The electron polaron, localized on a titanium 3d orbital, has a much more diffuse nature than the hole polaron, localized on an oxygen 2p orbital. Titanium and oxygen ions are represented by red and gray spheres respectively.

Journal of Chemical Theory and Computation
Article between the HOMO and the lowest unoccupied molecular orbital (LUMO). 42 The optimized values of α which satisfy the GKC for anatase and rutile are 0.105 and 0.115, which are very close together. In rutile, one could satisfy the GKC with |ξ| < 0.05 eV using 0.100 ≤ α ≤ 0.140, and for anatase one can use 0.095 ≤ α ≤ 0.115, meaning there is a large overlap of α values. This means that one could feasibly choose the same value of α to model both electron and hole polarons in these phases. This transferability in α values is somewhat unexpected; however, it allows for more than one phase to be described by the same functional. In later sections, we show that this transferability extends to other phases as well.
As the primary basis sets used are all of good quality, we expect that the ADMM basis will likely be a key source of inaccuracy in our calculations. To investigate the error introduced, we look at changes in bulk and polaron properties with differing titanium and oxygen basis sets. Table 1 illustrates that the choice of the ADMM basis set has little effect on both the trapping energy, local magnetic moment, and charge transition energy for both electron and hole polarons across both phases. We calculate magnetic moments via Bader analysis of the spin density. 88−90 The choice of basis set has no effect on the optimized value of α or the nonlinearity. Similarly, bulk quantities change very little with increasing the ADMM basis set quality. Using contracted basis sets, labeled with a lower case letter c, did not affect polaron or bulk properties. To this end, we find that using FIT9/pFIT3 provides the best balance between accuracy and computational efficiency.
In the preceding calculations we used an R c which was converged with respect to the predicted bulk lattice constant and band gap, resulting in the value R c = 6 Å. We now compare our parametrization with lower values of R c . By varying R c , the optimal value of α that satisfies the GKC (for the selected value of R c ) changes. As R c increases, the amount of exchange energy increases, resulting in a slightly more concave functional (see Figure 1). We have previously shown for rutile electron polarons that using R c = 6 Å with α = 0.115 satisfies the GKC with nonlinearity ξ < 0.05 eV. For the values R c = 4 Å and R c =2Åwe find that α = 0.115 and α = 0.2 satisfy the GKC to within 0.05 eV. The α parameter effectively remains the same when going from R c =6Åt oR c = 4 Å, suggesting that restricting exchange interactions over 4 Å yields similar results to 6 Å. However, for R c = 2 Å predicted polaron properties are very different. Taking R c = 6 Å yields a band gap of 2.8 eV, with the electron polaron CTL being 0.4 eV below the CBM. In contrast, R c = 2 Å yields a band gap of 3.1 eV and a polaron CTL that is roughly 0.9 eV below the CBM, much higher than that reported experimentally. This suggests that functionals with lower R c values have a tendency to predict more stable polaronic states. A similar conclusion can be drawn from the polaron trapping energy, where using R c =2Å and R c = 4 Å results in E t = −0.14 eV and E t = −0.02 eV respectively. Comparing the range-separation used in this work and the popular HSE functional, 27 HSE more smoothly weighs short-and long-range Coulomb interactions (using the complementary error and error functions). Short-range Coulomb interactions typically become negligible with a length of 2/ω. Values of ω are usually around 0.2 Å −1 , 94 meaning shortrange exchange is negligible beyond 10 Å.
The nonlinearity ξ for polaronic states also depends on the supercell size. We studied this for the following supercell sizes (number of atoms in brackets): 3 × 3 × 5 (270), 4 × 4 × 6 (576), and 5 × 5 × 8 (1200). Using the α optimized for the 576 atom system, nonlinearity was calculated for 270 and 1200 atom systems. We found ξ changed by 0.04 and 0.01 eV when moving from 270 to 576 and 576 to 1200 atom systems, respectively. The next reasonable (cubic) supercell size to test for nonlinearity is far too large (1944 atoms). However, as our chosen system size (576 atoms) has a nonlinearity only different by 0.01 eV to much larger systems (1200 atoms), we find that it provides a good balance between computational time and accuracy.
4.2. Brookite, TiO 2 (H), TiO 2 (R), and TiO 2 (B). Having benchmarked our approach for rutile and anatase, we now apply it to model polaron formation in four other phases of titanium . e As our approach predicts only anatase holes and rutile electrons are stable, only these two polarons are shown. The titanium basis FIT9 has 9 basis functions (FIT9), comprising of three s-, p-, and d-functions. This basis set has been optimized for TiO 2 and is shown in the SI. Titanium basis sets FIT10 and FIT11 have 10 and 11 basis functions, whilst the oxygen basis sets cFIT3, FIT3, cpFIT3, and pFIT3 have 4, 6, 5, and 7 basis functions. For each combination of the basis sets, we found the value of α which minimizes nonlinearity ξ. We found that the choice of basis set did not produce a difference in the optimized α. For anatase and rutile, we found the optimized α to be 0.105 and 0.115.

Journal of Chemical Theory and Computation
Article dioxide: brookite, TiO 2 (H), TiO 2 (R), and TiO 2 (B). Table S4 (SI) shows the bulk crystal properties of these phases, which each possess a number of lattice sites that are inequivalent by symmetry (see Figure 5). We denote lattice sites using shorthand notation: e.g. O a Br and O b Br label oxygen sites a and b in the brookite crystal. To avoid confusion, we abbreviate rutile and anatase to Ru and An, while TiO 2 (H), TiO 2 (R), and TiO 2 (B) are abbreviated as H, R, and B.
For all localized polarons, we determine the range of α values which satisfy |ξ| < 0.05 eV (shown in Figure 6). We find that across a large number of sites, there is an overlap in the range of acceptable α values. This means that it is feasible to model all phases and sites using the same α while still satisfying the GKC. We find α = 0.115 ensures |ξ| < 0.08 eV across all phases and sites. In the following subsections we describe the predicted properties of polarons in each phase in detail. We evaluate the polaron properties of each phase using an α that minimizes nonlinearity for that phase. This allows us to compare the CTLs of different polaron sites using the same functional (and hence the same lattice constant and band gap). These α values are shown in Table 2 alongside the predicted bulk properties. If instead, one uses α = 0.115 to model all phases, good quantitative agreement is seen. The relevant data comparing phase-specific values of α and α = 0.115 is shown in the SI.In Table 2 the stability of each of the phases relative to anatase is shown, where it can be seen that the energy differences between phases are small. For phases where data is available (rutile, anatase, and brookite), the energy differences are in line with previous (local) DFT calculations. 95,96 Diffusion quantum Monte Carlo results yield, ignoring the effects of lattice dynamics, a similar energetic ordering. 97,98

Polarons in Brookite.
We find that brookite is in many ways similar to anatase. Specifically, the GKC predicts the formation of hole polarons but not electron polarons. This is in good agreement with EPR data, which also suggests that the polaronic properties of the two crystals are quite similar. 63 The two inequivalent O sites in brookite, differing mainly by bond angles to neighboring Ti, are both capable of forming hole polarons. Both sites have (+/0) CTLs differing only by 0.02 eV at the optimal α for brookite (α = 0.105), as seen in Table 3. However, we find that their trapping energies differ by around 0.1 eV, making O a Br the most stable trapping site in brookite. The predicted band gap of brookite is 3.2 eV, consistent with the experimentally determined values of 3.1−3.4 eV. 99 4.2.2. Polarons in TiO 2 (H) and TiO 2 (R). We now discuss the phases TiO 2 (H) and TiO 2 (R) together, as our predictions suggest they exhibit similar charge trapping properties. TiO 2 (H) contains two unique oxygen sites and one unique titanium site. We find that the best α for all sites in this phase is α = 0.130, as seen in Table 2. A small hole polaron, similar to that formed in anatase, is predicted on one oxygen site (O b H ). On the second

Journal of Chemical Theory and Computation
Article oxygen site (O a H ), a multisite small polaron, or molecular polaron, is predicted to form previously unseen in titania. The multisite polaron, illustrated in Figure 7, is predicted to localize over two ions, with almost equal charge on both ions. Comparing the single-and multisite hole polarons, we find the (+/0) CTL from the multisite polaron lies deeper in the band gap (with respect to the VBM) by 0.1 eV. Furthermore, the multisite polaron is also more energetically stable, with a trapping energy 0.1 eV lower than the single-site polaron. On the only unique titanium site in TiO 2 (H), a multisite electron polaron can be found (illustrated in Figure 7). The (−/0) CTL of this polaron is twice as deep as the electron polaron found in rutile, lying 0.8 eV below the CBM.
The two oxygen (O a H ) ions on which the multisite hole polaron forms (see Figure 7) move closer together by 0.3 Å compared with their equilibrium position. In contrast, two nearby titanium ions (shared by the two oxygen ions on which the multisite polaron resides) move apart by 0.3 Å. The singlesite hole polaron (O b H ) results in an asymmetric distortion of nearby titanium ions, with bond lengths increasing by 0.1−0.3 Å from their defect-free positions. With geometries compared, it is clear that the nature of polarons of these two sites is different. Despite large changes in the geometry compared to the equilibrium positions, the multisite polaron remains more energetically stable than the single-site one. We find a similar behavior in the multisite electron polarons, where two titanium ions involved move together by 0.2 Å, while two surrounding oxygen ions are pushed apart by 0.2 Å.
Single-and multisite polarons are predicted to form in TiO 2 (R), on sites O a R and O b R respectively. The geometric distortions of both polarons are also quite similar to the TiO 2 (H) case. We note that in both phases multisite polarons are more stable than single-site polarons, as seen in Table 3.
One distinguishing factor between single-site and multisite polaron sites in the TiO 2 (H) and TiO 2 (R) phases is the bonding. The multisite electron polarons form over two cations ions that bond with two of the same anions shown in Figure 7. Anions which display the same double-sharing behavior also form multisite hole polarons. The multisite polarons reside on sites of the same symmetry. , is four-coordinated, and the GKC predicts a localized hole polaron is not possible on this site.

Polarons in
The electron polaron case in TiO 2 (B) is slightly more complicated than the hole polaron case. While the sites Ti a B and Ti b B are both capable of forming small polarons that satisfy the GKC to within 0.06 eV, the polarons begin to delocalize over the range 0 < ξ ≤ 0.05 eV. Small polarons on Ti a B and Ti b B delocalize at α = 0.145 and α = 0.140 respectively, with the more delocalized states having similar total energies to the small polarons. At α = 0.145, the small electron polaron on site Ti b B has a( −/0) CTL of 0.40 eV and a trapping energy −0.04 eV. As electron polarons are capable of satisfying the GKC over a small range, we find it plausible that they form in TiO 2 (B). However, the small trapping energies indicate that they may be less likely to form.
4.3. Comparison of Polarons Across Phases. The trapping energies, CTLs, and Hirshfeld spin density analysis of each polaron are shown in Table 3. These values are calculated using the ideal value of α for each phase (shown in

Journal of Chemical Theory and Computation
Article addition, multisite electron polarons (Ti H and Ti R ) are much more stable than the single-site electron polaron found in rutile. The deepest (−/0) and (+/0) CTLs correspond to multisite polarons in TiO 2 (H) and TiO 2 (R In contrast, CTLs from multisite polarons are distinctly split from those of small polarons, as seen in phases TiO 2 (H) and TiO 2 (R). The quantities given in Table 3 calculated from the transferable α (α = 0.115) are also shown in Table S5 (SI). As previously mentioned, the transferable and phase-specific α values are in good quantitative agreement. Further, for the range of α values determined for each polaron (given in Figure 6), we show the corresponding range of trapping energies, CTLs, and phase band gaps in Tables S6 and S7 (SI).
To compare polaron CTLs across different phases, the band alignment of the phases with respect to each other is required.
Here we align the VBM of all phases using the ionization energies obtained in a previous theoretical study, 100 which calculates these energies through a hybrid quantum mechanical/ molecular mechanical embedded cluster approach. This allows us to compare the calculated band gaps and polaron CTLs for all phases (Table 3) against a common vacuum level reference ( Figure 8). We find the deepest electron and hole states from the vacuum level are the Ti H and O b H polarons. In contrast, the shallowest electron and hole polaron sites are Ti R and O a B . It can be seen that, across phases, electrons and holes have an energy range which is separated by roughly 0.2 eV when aligned to a common vacuum level. Of the phases explored in this work rutile is the only phase of TiO 2 which cannot form hole polarons. A figure illustrating the alignment in Figure 8 for the transferable α is shown in Figure S2 (SI).

DISCUSSION
We now discuss some of the factors that may influence the accuracy of our predictions. As noted in Section 3 the use of the ADMM basis set allows for a vast speedup in calculation of exact exchange but also may introduce inaccuracies. However, we compared the bulk and polaronic properties of rutile and anatase with a variety of titanium and oxygen basis sets and found that the choice of the ADMM basis set has little influence. Our choice of pFIT3/FIT9 yields similar results to pFIT3/FIT11, suggesting that including polarization (f-)functions for titanium are largely unnecessary in these materials. Similarly, pFIT3/ FIT9 gives similar results to FIT3/FIT9, also suggesting that polarization (d-)functions for oxygen are not required. These tests suggest our choice of the pFIT3/FIT9 is sufficiently accurate for the prediction of polaron properties.
Another potential source of error is the truncation of exact exchange in the hybrid functional beyond a cutoff radius R c . Our approach was to converge predicted bulk properties with respect to R c to mimic the PBE0 functional but with reduced computational cost. This approach yielded R c = 6 Å. However, we also tested the effect of using smaller values of the cutoff radius. While results obtained using R c = 4 Å were not significantly different, reduction of the cutoff to 2 Å led to overstabilization polarons, absence of transferability between phases, and predicted polaron properties inconsistent with experiment. It is notable that empirical attempts to tune the amount of HF exchange in the PBE0 functional to give a reasonable description of the band structure of anatase TiO 2 arrived at α = 0.125. 101 This is very consistent with the value we have obtained from application of the GKC providing further support for the approach.
Our previous discussion (in Section 2) highlights that (small) electron polarons are stable in rutile but not in anatase. In addition, some evidence suggests hole polarons are more stable in anatase compared to rutile. This in agreement with both experiment 58,63,72 and other theory. 9,66,67,71 We have shown that the predicted band gaps of 2.82 and 2.94 eV for rutile and anatase are in very good agreement with the experimentally observed values of 3.0 and 3.2 eV. The rutile (−/0) electron polaron CTL is found to be 0.4 eV below the CBM in our work, which is in agreement with experimental data placing the level at 0.3, 0.5, and 0.7 eV. 58,64,65 We have shown that the rutile and anatase CTL data is also consistent with HSE06 findings. 9 In addition, the trapping energies of −0.02 eV and −0.20 eV for rutile and anatase are found to be in good agreement with HSE06. We found that anatase and brookite have very similar charge trapping properties, behavior which has been seen experimentally in ref 63. Specifically, the two phases have much more delocalized electron polarons compared to the rutile case. Experimentally, little information exists on polaronic trapping in the phases TiO 2 (H), TiO 2 (R), and TiO 2 (B). Therefore, our predictions are extremely valuable to guide experimental work.

CONCLUSION
We have demonstrated that hybrid functionals coupled with the generalized Koopmans' condition gives predictive capabilities with regards to localized charges in TiO 2 . We were able to assess the stability and properties of small electron and hole polarons in various TiO 2 phases. Furthermore, by fitting the fraction of exact exchange to reduce nonlinear behavior of energy with respect to the occupation number of these polaronic states, we were able to obtain a band gap differing from experimental results only around 6%. This suggests that our approach is not only capable of predicting polaron localization properties but also results in well-behaved bulk crystal properties. In addition to this, we showed that the fraction of exact exchange α is transferable between phases. Using α = 0.115 along with the PBE0-TC-LRC functional satisfies the GKC to within 0.08 eV for all phases. We have also explored parameters that influence the key results, such as the ADMM basis set, ensuring that they are of sufficient quality. Our work suggests that the approach we use is of reasonable quality, inexpensive, and transferable and should be  Table 2. For clarity, we show the subscript letter index only for each of the polaron sites. Electron and hole polaron energy ranges with respect to vacuum are labeled as e and h.

Journal of Chemical Theory and Computation
Article easily extendable to other metal-oxides capable of polaron formation.
Using this approach on the brookite, TiO 2 (H), TiO 2 (R), and TiO 2 (B) phases has shown that intrinsic polaron formation is not unique to rutile and anatase. Unlike the remaining phases, TiO 2 (H) and TiO 2 (R) display multisite electron and hole polarons, arising from the unique nature of bonds between anions and cations. TiO 2 (B) exhibits (small) hole polarons but has shallow electron polaron states. While theoretical exploration of brookite is scarce, EPR data suggests the phase capable of forming hole polarons but not electron polarons, in agreement with our own findings. Somewhat surprisingly, we find that rutile is the only phase of titania which is incapable of forming hole polarons. The techniques employed in this work have been shown to be consistent with experimental data where available and transferable across all of the TiO 2 phases considered. Furthermore, the computational efficiency of the approach means it can be applied to more complex systems, such as surfaces, interfaces, and doped systems, to provide greater insight into the nature of charge trapping in materials.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.8b00199.
Information on convergence, as well as predicted ranges of bulk and polaron properties; polaron data using the transferable hybrid parametrization for TiO 2 (PDF) Funding K.P.M. acknowledges support from EPSRC (EP/K003151/1, EP/P006051/1, and EP/P023843/1). This work made use of the facilities of Archer, the UK's national high-performance computing service, via our membership in the UK HPC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202/1). This work also made use of the facilities of N8 HPC Centre of Excellence, provided and funded by the N8 consortium and EPSRC (EP/K000225/1). The Centre is coordinated by the Universities of Leeds and Manchester.