Benchmarking the Fundamental Electronic Properties of small TiO 2 Nanoclusters by GW and Coupled Cluster Theory Calculations

: We study the vertical and adiabatic ionization potentials and electron a ﬃ nities of bare and hydroxylated TiO 2 nanoclusters, as well as their fundamental gap and exciton binding energy values, to understand how the clusters ’ electronic properties change as a function of size and hydroxylation. In addition, we have employed a range of many-body methods; including G 0 W 0 , qs GW , EA/IP-EOM-CCSD, and DFT (B3LYP, PBE), to compare the performance and predictions of the di ﬀ erent classes of methods. We demonstrate that, for bare clusters, all many-body methods predict the same trend with cluster size. The highest occupied and lowest unoccupied DFT orbitals follow the same trends as the electron a ﬃ nity and ionization potentials predicted by the many-body methods, but are generally far too shallow and deep respectively in absolute terms. In contrast, the Δ DFT method is found to yield values in the correct energy window. However, its predictions depend upon the functional


■ INTRODUCTION
Titanium dioxide (TiO 2 ) nanostructures are key components of many 1−3 advanced technological products and processes, for example as photocatalysts for water and air purification 4 and photocatalytic water splitting, 5−11 as semiconductors in dyesensitized solar cells 12,13 and as photocatalytic and superhydrophilic coatings for self-cleaning windows. 14−16 All these applications depend on a combination of the wettability of TiO 2 by water or other polar solvents, its chemical and mechanical stability in the presence of water and illumination, and the specific optical and/or electronic properties of the TiO 2 nanostructures. TiO 2 is also the archetypal model system for oxide semiconductors where a large amount of fundamental experimental 17−19 and computational 20−59 studies on the optical, electronic, and photocatalytic properties of both bulk and nanostructured TiO 2 have been performed.
Key properties relevant for understanding the optical and electronic properties of TiO 2 nanostructures and materials in general 60 are (i) the fundamental gap (ΔE f , often referred to as the band gap for crystalline solids), (ii) the optical gap or absorption onset (ΔE 0 ), (iii) the exciton binding energy (EBE), and (iv) their ionization potential (IP) and electron affinity (EA); see Figures 1 and 2. The optical gap is the energy above which a system starts absorbing light and singlet excitons are formed, while the fundamental gap is the energy required to make pairs of free charge carriers. The exciton is a bound state of an excited electron in the conduction band and a hole in the valence band, in which the excited electron and hole are attracted to each other by their opposite charge. Triplet excitons can be formed from their singlet counterparts through intersystem crossing mediated by spin−orbit coupling but are not directly created upon light absorption due to the spin forbidden nature of the underlying transition. The exciton binding energy is the energy needed to separate the electron and hole (i.e., to convert them to free charge carriers). For systems that lack strict periodicity (i.e., clusters, polymers, and so on) the difference between the fundamental gap and optical gap is purely a measure of the extent to which singlet excitons are energetically stabilized with respect to free charge carriers due to the electrostatic interaction. Systems are inherently transparent to light below the optical gap and below the fundamental gap cannot carry a photocurrent. The energies of the valence and conduction bands relative to the vacuum level and the charge carriers in them correspond to the negative of the ionization potential (IP, the energetic cost of removing an electron and forming a cation) and electron affinity (EA, the energy released upon adding an electron and forming an anion). All these properties can in principle couple with phonons, i.e., nuclear degrees of freedom. For example, the experimental onset of light absorption will be slightly red-shifted with respect of the vertical optical gap due to contributions of (nonvertical) vibronic transitions. Similarly, while the lowest energy peak maxima in a nonvibrationally resolved photoelectron spectrum will approximately coincide with the vertical IP, the onset of electron emission is linked to the adiabatic IP, which includes the effect of nuclear relaxation after removal of an electron. Here, in analogy with the typical theoretical convention for solids, we will concentrate on vertical properties, but for selected clusters, we also explicitly explore the effect of nuclear relaxation.
The majority of computational studies on TiO 2 nanostructures are based on density functional theory (DFT) 61 and methods such as TD-DFT and GW, which are built on top of DFT. Specifically, the optical and/or fundamental gap within DFT is approximated by the difference between the highest occupied and lowest unoccupied (generalized) Kohn−Sham (KS) orbitals (Kohn−Sham (KS-) gap). 62 Other DFT studies 25,31,52,53,57 calculate the vertical and/or adiabatic electron affinity (EA) and ionization potential (IP) of nanostructures explicitly using a ΔSCF (ΔDFT) approach 63,64 and obtain the fundamental gap by taking the difference between the vertical EA and IP. Similar studies 26,31,44,45,51,56 use   . B3LYP/def2-TZVP optimized structures of the bare and hydroxylated TiO 2 nanoparticles studied. The (TiO 2 ) 2 trans, (TiO 2 ) 3 , (TiO 2 ) 4 , (TiO 2 ) 5 , (TiO 2 ) 6 , (TiO 2 ) 7 , and (TiO 2 ) 8 nanoclusters correspond to global minima candidate structures for that number of TiO 2 units, whereas (TiO 2 ) 2 cis, (TiO 2 ) 2 club, and (TiO 2 ) 3 alt represent higher lying minima. The hydroxylated (TiO 2 ) 2 (H 2 O) 2 cluster was manually constructed by adding -H and -OH groups to the bare (TiO 2 ) 2 global minimum candidate. Titanium, oxygen, and hydrogen atoms are shown as gray, red, and white spheres, respectively. time-dependent DFT 65 (TD-DFT) to calculate the optical gap as the lowest energy excited state of the nanostructure. Finally, there are studies 27,31,35,38,40,[44][45][46]48,55,58 that go beyond DFT either using many-body perturbation theory; Green's function based GW theory for the fundamental gap, 66−68 and GW in combination with solving the Bethe−Salpeter equation (BSE) for the optical gap, 69 as well as quantum chemistry methods. 70−73 The predominance of calculations based on (TD-)DFT is directly linked to the very favorable computational scaling with system size relative to the inherently more accurate alternatives, allowing for calculations on true nanosized systems. The drawback of (TD-)DFT, however, is the need in practice to use approximate exchange-correlation (XC) density functionals, as the functional form of the exact density functional is not known. The results obtained with (TD-)DFT are hence to a degree functionally dependent and ideally need to be validated by comparison with experimental or computational benchmark data.
Here we build further on our previous efforts 40, 44,45 to generate benchmark data for the optical gap and photoluminescence energy of small TiO 2 model nanoclusters using equation-of-motion coupled-cluster (EOM-CC) theory. We predict vertical and adiabatic −EA and −IP values and vertical fundamental gap values and exciton binding energies for the same model clusters using both Green's function based methods: G 0 W 0 , 66−68 quasi-particle self-consistent qsGW, 74−77 and coupled-cluster 70 (CC) theory based approaches: ΔCC and EA/IP-EOM-CC. 78 Such nanoclusters, see Figure 3, are arguably smaller than the nanostructures encountered experimentally in all but mass spectrometry type of experiments 17,18,79,80 but can be conveniently studied by both DFT and GW/CC methods. We also compare the results of different GW and CC methods, as well as their convergence with increasing basis-set size, and comment on the suitability of each to describe TiO 2 clusters. In the case of EA we also compare our G 0 W 0 and EA-EOM-CCSD predictions to experimental data. Finally, we consider the effect of hydroxylation on the fundamental gap and exciton binding energy of TiO 2 clusters.
GW and CC Methods. G 0 W 0 and qsGW. The central object in the GW formalism is the causal Green's function G. Formally it is obtained from the Dyson equation where H H is a noninteracting reference Hamiltonian, i.e., the Hartree Hamiltonian, and Σ is the so-called self-energy. In the GW approximation the self-energy is taken as the first-order expansion of the self-energy in terms of the screened Coulomb interaction W:Σ = GW. The zero-order expansion corresponds to Hartree−Fock (HF) theory. In general Σ is nonhermitian and energy-dependent. In the quasi-particle approximation the Green's function is expanded in terms of quasi-particles ψ QS with eigenvalues ε QS : The sum runs over all occupied and unoccupied single-particle states; μ is the Fermi energy. In this approach the Green's function hence keeps the analytical form of a noninteracting Green's function. In the G 0 W 0 approach, the ionization energies and electron affinities are obtained from the quasi-particle energies, ε QP , which are evaluated as a first-order perturbative correction to a set of noninteracting single-particle states. In practice these are usually the Kohn−Sham orbitals, ψ KS , with eigenvalues, ε KS : Although eq 3 is solved iteratively for ε QP , due to its perturbative nature the final ε QP values still depend on the starting point, that is, on the (arbitrary) set of eigenvalues ε KS and orbitals ψ KS . In this work we use the B3LYP XC functional 81,82 to generate the initial set of KS eigenvalues and orbitals for our G 0 W 0 calculations. 77 In qsGW self-consistency is reached at the level of the quasiparticle orbitals and energies. 74−76 In contrast to G 0 W 0 where only the diagonal matrix elements of the self-energy are needed, in qsGW also the off-diagonal matrix elements are required. To lift the nonhermiticity of Σ and obtain a single set of orthonormal quasi-particle states, the exact self-energy matrix elements are replaced by The solution of eq 1 is organized in an iterative scheme starting from the KS initialization. The QP orbitals of the (i + 1)th iteration ψ (i+1) (r) are expressed in the reference orbitals of the previous iteration, defining the coefficient matrix A: In this reference basis eq 1 becomes an eigenvalue problem: which yields updated ε n (i+1) and A n′n (i+1) , where the latter are used to construct new orbitals ψ n′ (i+1) (r). The final result of eq 6 has been shown to be independent of the starting point, 76 but both the stability of the iterative cycle and the rate of convergence can be greatly improved by using an optimal starting point. Hybrid functionals such as the B3LYP XC functional are known to give good starting points. Moreover, the actual implementation contains a linear mixing scheme in eq 6 to improve stability.
ΔCC and EA/IP-EOM-CC. In the EA/IP extension of the equation-of-motion CC approach 78 the kth wave functions of N ± 1 particle systems are represented in the form of the EOM-CC Ansatz where |Φ⟩ , T, and R(k) are the reference Slater determinant, typically from a previous Hartree−Fock calculation, the cluster operator for the N-electron system, and the excitation operator, respectively. In the rudimentary EA/IP-EOM-CCSD approximations the IP and EA excitation operators, R IP (k) and R EA (k) are represented as (9)

Journal of Chemical Theory and Computation
Article where a p and a + p are the annihilation and creation operators for the electron in pth spin−orbital and following the convention labeling occupied orbitals with i and j and unoccupied orbitals with a and b. The EOM-CCSD ionization potentials, electron affinities, and corresponding excitation operators are obtained by diagonalizing a similarity transformed Hamiltonian in normal-product form H̅ N (see ref 78 for details) in the N ± 1 particle subspaces spanned by single and double excitations; i.e., where the ω(k) values correspond either to ionization potentials or electron affinities. Contrasting Both Approaches. While both classes of methods have rather different starting points and are described in terms of different mathematical objects, they are also closely interrelated. For example, the fully interacting CC Green's function can in principle be obtained from the CC wave function and then be used to extract the CC self-energy Σ through the Dyson equations. Indeed, the EA/IP-EOM-CCSD methods used here are closely related to the Green's function− CCSD method 83−90 explicitly developed for this purpose, which utilizes the bivariational formulation of the CC formalism. The matrix elements of the CC Green function matrix, G pq , take the following form: where ⟨Φ|(1 + Λ) is the left eigenvector of the similarity transformed Hamiltonian H̅ = e −T He T and ̅ a q + and ̅ a p are the similarity transformed creation and annihilation operators a q + and a p ; i.e.,   Table 1). All G 0 W 0 calculations were started from a DFT calculation with the B3LYP XC functional. All values are in electronvolts.
The X p (ω) and Y q (ω) operators in eq 11 are the ionization potential/electron affinity type operators which assume singular values for ω values corresponding to the EA/IP-EOMCC eigenvalues ω k . In other words, all poles of the CCSD Green's function correspond to the eigenvalues of the IP/EA-EOM-CCSD equations.
Both classes of many-body methods can be systematically improved. By adding higher order excitations to the EA/IP-EOM-CC excitation operators, in analogy to the standard EOM-CC formulations, one can define a hierarchy of approximations converging to the full configuration interaction limit once all possible excitations are included. Similarly, GW can be improved upon by going beyond the zero-and firstorder terms of the self-energy expansion. In practice this knowledge is of limited use in predicting the relative accuracy of qsGW and EA/IP-EOM-CCSD for systems of interest. G 0 W 0 , being an approximation to qsGW, should be less accurate than qsGW, but the difference between the two methods can be minimized by choosing an optimal starting point. 76,91 Computational Details. The (vertical) ionization potential, electron affinity, fundamental gap, optical gap, and exciton binding energy were calculated for a series of (TiO 2 ) n clusters and two (TiO 2 ) n (H 2 O) m clusters (see Figure 3). The structures of the clusters discussed here were taken from our previous studies 40, 44,45 and were originally obtained through global optimization, 20,26,38,92 with the exception of the hydrated structures which were constructed manually by adding -H and -OH groups to the bare global minima candidate structures. 44 In our earlier works 40,44,45 the xyz coordinates of all the clusters were optimized using a combination of the B3LYP 81,82 hybrid XC functional and the def2-TZVP 93,94 basis set, using Turbomole 6.6. 95,96 As discussed in the Introduction, IP and EA of the clusters were calculated using a range of approaches: from taking the DFT KS orbital energies to using the ΔDFT, G 0 W 0 , qsGW, ΔCCSD, ΔCCSD(T), ΔCCSDT, and EA/IP-EOM-CCSD methods. The optical gap of selected clusters was approximated as the lowest vertical singlet excitation energy and calculated using EOM-CCSD. All the coupled cluster calculations, for reasons of computational tractability, employed the frozen core approximation where only the valence electrons are correlated. The (Δ)DFT and GW calculations were performed in Turbomole 6.6, 95,96 the qsGW calculations in a locally modified version of Turbomole 7.0, 76 while the coupled cluster calculations used NWChem 6.6. 97 All these calculations, finally, employed basis sets from the Ahlrichs 93,94,98 (def2-SVPD, def2-TZVPP, and def2-QZVPP) and Dunning 99,100 (aug-cc-pVTZ, aug-cc-pVQZ, and aug-cc-pV5Z) families of Gaussian basis sets.

■ RESULTS AND DISCUSSION
Effect of Basis-Set Size. We analyze the effect of basis-set size on the predictions of vertical −IP and −EA values by the ΔB3LYP, G 0 W 0 , and EA/IP-EOM-CCSD methods through a comparison of results for the TiO 2 molecule calculated with a range of basis sets from Ahlrichs and Dunning basis-sets families, where possible, extrapolating the results to the complete basis-set limit. Figure 4 shows the predicted vertical −IP and −EA values calculated with ΔB3LYP, G 0 W 0 and EA/ IP-EOM-CCSD and the Ahlrichs and Dunning basis-set families plotted against the number of basis functions in the basis set (Nbasis). Table S1 in the Supporting Information presents the same information in table form. Concentrating first on the results obtained with basis sets from the Ahlrichs family (filled symbols), the vertical −IP and −EA values are generally predicted to become more negative with increasing basis-set size, except for the case of ΔB3LYP, where −IP and −EA do not really change with basis-set size. Moreover, the change with basis set is generally larger for −IP than in −EA, at least in absolute terms. In addition, it appears that G 0 W 0 shows a much stronger dependence on the size of the basis set used than either EA/IP-EOM-CCSD or ΔB3LYP. Contrasting the results obtained with basis sets from the Ahlrichs and Dunning families (the open symbols), it seems that when comparing basis sets of the same cardinality (e.g., def2-TZVP and aug-cc-pVTZ), the Dunning basis sets always yield slightly more negative −IP and −EA values, in line with the fact that the latter always contain more primitive basis functions, as well as proper diffuse functions. Table 1 provides the extrapolations of the G 0 W 0 −IP and −EA values to the complete basis-set limit (CBS). Three different extrapolation methods were tried: (A) a linear extrapolation of −IP and −EA as a function of one over the number of basis functions (Nbasis) in each of the basis sets, i.e., a + b/Nbasis; (B) a linear extrapolation of −IP and −EA as a function of one over the cube of the cardinality (CN) of each of the basis sets, i.e., a + b/CN 3 ; and (C) an exponential extrapolation as a function of CN, i.e., a + b exp(c(CN)). All three extrapolation methods give similar results for −EA in the CBS limit, while for −IP in the CBS limit A yields a more negative value than either B or C. Regardless, the results of the extrapolation suggest that, as expected from the results in Figure 4, G 0 W 0 even for triple-ζ basis sets yields −EA values that deviate 0.5−0.6 eV and −IP values that differ by 0.2−0.3 eV from the CBS limit. Comparing these deviations to those found in a recent work studying a selection of 100 molecules with G 0 W 0 , 101 it appears that in the case of TiO 2 G 0 W 0 performs worse than average for −IP and better than average for −EA.
We cannot perform the same extrapolation for EA/IP-EOM-CCSD, due to linear dependency difficulties when solving the EA/IP-EOM-CCSD equations for the largest of the Dunning family basis sets. However, as discussed above, the data in Figure 4 and Table S1 show that the EA/IP-EOM-CCSD −IP and −EA values vary much less with increasing basis-set size than those obtained with G 0 W 0 , suggesting that the deviation between triple-ζ (and quadruple-ζ) basis set and CBS limit values for EA/IP-EOM-CCSD might be smaller than for G 0 W 0 . This weaker effect of basis-set size is in line with previous benchmark EA/IP-EOM-CCSD calculations for the C 2 and F 2 molecules. 78 For ΔB3LYP, not surprisingly, the results are not very sensitive to the basis-set size beyond triple-ζ basis sets.

Journal of Chemical Theory and Computation
Article For all methods, the effect of basis-set size on the predicted fundamental gap, discussed in more detail below, is much smaller than on the individual vertical −IP and −EA values. For example, for G 0 W 0 the fundamental gap changes only from 7.2 to 7.3 eV when going from the def2-SVPD to the aug-ccp-pV5Z basis set, where the latter, depending on the extrapolation method used, is approximately 0.1−0.2 eV from the extrapolated fundamental gap in the CBS limit of 7.4−7.5 eV. This much smaller change of the fundamental gap with basis-set size than in the constituting vertical −IP and −EA values is due to both displaying, as discussed above, large changes in the same negative direction.
For reasons of computational tractability, we limit ourselves in the remainder of this work to calculations using the def2-SVPD, def2-TZVPP and aug-cc-pVTZ basis sets. This might mean that our results for certain methods, e.g., G 0 W 0 , deviate in terms of absolute values from those at the basis-set limit, but we expect that trends will remain the same.
Comparing G 0 W 0 and qsGW. Table 2 compares the vertical −IP and −EA values of the TiO 2 molecule and (TiO 2 ) 2 and (TiO 2 ) 3 clusters obtained with DFT, G 0 W 0 and qsGW. Concentrating on the latter two methods, G 0 W 0 and qsGW appear to predict the same trends, which are discussed below in more detail, but qsGW predicts −IP values that are consistently deeper, by ∼1 eV, and −EA values that are consistently shallower, by ∼0.3 eV, than G 0 W 0 . This shift, at least for the −IP values, is similar in magnitude and direction to those previously observed when applying G 0 W 0 and qsGW to a range of simple organic and main-group molecules. 76 Use of aug-cc-pVTZ rather than def2-TZVPP in the G 0 W 0 calculations gives, in line with what was discussed above for the TiO 2 molecule, very similar results. The effect of adding diffuse functions hence appears small.
The −EA values of some of the small TiO 2 clusters were previously studied using G 0 W 0 by Marom and co-workers. 38 Our G 0 W 0 /B3LYP/def2-TZVPP and G 0 W 0 /B3LYP/aug-cc-pVTZ −EA values in Table 2 for the different isomers of (TiO 2 ) 2 and the (TiO 2 ) 3 global minimum candidate lie within 0.1 and 0.2 eV, respectively of Marom and co-workers' G 0 W 0 / PBEh values obtained using FHI-Aims and large (tier 4) numerical basis sets. The minor differences between our and their values are probably due to differences in the basis-set quality and the density functional, B3LYP vs PBEh, used to obtain the starting eigenvalues and eigenvectors.
As an aside, comparing the def2-TZVPP and def2-SVPD ΔB3LYP results in Tables 2 and S2, it is clear that while, as discussed above, the ΔB3LYP results in general have a weak dependency on the basis-set size, the predictions for the −IP values of the club isomer of (TiO 2 ) 2 and the (TiO 2 ) 3 global minimum candidate structure are significantly different. A comparison of the spin population for the cation structures in these cases suggests that this is the result of the SCF with the two basis sets converging to solutions with different distributions of the unpaired spin over the clusters.

Journal of Chemical Theory and Computation
Article (TiO 2 ) 2 dimer clusters, and the (TiO 2 ) 3 global minimum candidate cluster, as obtained with EA/IP-EOM-CCSD, ΔCCSD, ΔCCSD(T), and ΔCCSDT. The clusters seem to be dividable into two distinct sets: clusters for which EA/IP-EOM-CCSD and ΔCC methods predict similar values (i.e., the TiO 2 molecule and the (TiO 2 ) 2 trans cluster and (TiO 2 ) 3 global minimum candidate structure) and clusters for which EA/IP-EOM-CCSD and ΔCC methods give significantly different results for the IP (i.e., the (TiO 2 ) 2 cis and club clusters). The latter clusters are cases where the coupled-cluster calculation on the cation has large T1 amplitudes (e.g., 0.43 and 0.21 for the cis and club clusters, respectively, with CCSD/augcc-pVTZ). The large amplitudes suggest a significant multiconfigurational character of these cationic states and in turn that the N − 1 electron ROHF wave function might not be the best choice for a reference state. This in turn also might mean an unbalanced description of triples in CCSD(T) on this state. Since in EA/IP-EOM-CCSD the N electron RHF wave function of the neutral cluster is used instead as reference, such problems are circumvented and in such cases EA/IP-EOM-CCSD should yield more balanced results. The latter assessment is supported by the fact that the EA/IP-EOM-CCSD results show the same trends as G 0 W 0 and qsGW, while the ΔCCSD and ΔCCSD(T) results display a completely different trend. For the other clusters there are some minor differences, especially with respect to ΔCCSD(T), probably related to the absence of a description of triples in EA/IP-EOM-CCSD. The difference between def2-TZVPP and aug-cc-pVTZ results, and thus the effect of adding diffuse functions, finally appears to be again small.  38 ΔPBE results can be found in Figure S1 in the Supporting Information.  Tables 2, 3, and S2, G 0 W 0 , qsGW, and EA/IP-EOM-CCSD predict very similar trends for −IP and −EA with cluster size and between the different cluster isomers of (TiO 2 ) 2 , irrespective of the basis set used. The three many-body methods differ, however, in the predicted absolute values. qsGW consistently predicts the deepest −IP values and G 0 W 0 the most shallow. For −EA the differences are, just as in the case of the effect of basis-set size, smaller in absolute terms but larger in relative terms. ΔB3LYP predicts −IP and −EA values in roughly the same energy range as the three many-body methods, but displays a different ordering in the case of −IP, as well as suffers from the basis-set issues discussed above. KS-B3LYP predicts the same ordering as the many-body methods but yields much shallower −IP values and much deeper −EA values. While we only consider vertical potentials at this stage and hence cannot directly compare to experimental data, we expect similar observations to hold for adiabatic potentials, something that we consider in more detail below.
Any issue with describing −IP and −EA inherently carries over to the description of the fundamental gap by any given method combination. Not surprisingly thus all three manybody methods predict a similar nonmonotonic evolution of the fundamental gap with cluster size, where, in line with the error cancelation for the fundamental gap discussed above, the effect of the basis-set size on the calculated gap is very small for both G 0 W 0 and EA/IP-EOM-CCSD. qsGW consistently predicts the largest fundamental gap values while G 0 W 0 predicts the smallest with EA-IP-EOM-CCSD fundamental gap values lying between those of the two GW approaches. ΔB3LYP predicts fundamental gap values in the same energy range as the many-body methods but with a distinctly different pattern of maxima and minima in the fundamental gap for the smaller clusters, primarily due to its inability to reproduce the manybody method trend for −IP, discussed above. KS-B3LYP (see Figure S3 in the Supporting Information) in contrast predicts the same correct trend in the fundamental gap as the manybody methods but, not surprisingly, absolute gap values that are approximately twice as small as those found by G 0 W 0 and EA/ IP-EOM-CCSD.
Concentrating on the EA/IP-EOM-CCSD/def2-SVPD results, it appears that the clusters' IP values generally decrease with increasing cluster size for the global minima candidate structures. A similar trend can be seen for the −EA values, but here the undulations superimposed on this apparent trend are relatively larger than for the −IP values such that, over the size range considered, the (TiO 2 ) 3 minimum candidate structure has the deepest −EA value. The fundamental gap, finally, first appears to increase and then becomes more or less constant with increasing cluster size, but shows a dramatic dip for (TiO 2 ) 3 . This deep dip in the fundamental gap for (TiO 2 ) 3 appears to be the direct result of the cluster's relatively deep −EA value, which might be linked to the fact that the structure has an exposed three-coordinated titanium atom with a low onsite electrostatic potential; see below.
Considering the additional low-energy isomers studied for (TiO) 2 and (TiO) 3 , it is clear, see Table S2, that higher energy structures can have either more or less negative −IP values than the corresponding global minimum candidate structure, either more or less negative −EA values, and either larger or smaller fundamental gap values. It thus stands to reason that the trends in −EA, −IP, and the fundamental gap shown in Figures 5−7 are more likely the result of an underlying trend with cluster size in the types of structures that lie low in energy and the types of structural motives present in them than a direct effect of the cluster size per se.
Multiconfigurational Effects in EOM-CCSD. To illustrate the electronic structure of the N ± 1 electronic states corresponding to lowest IPs and EAs, in Table 4 we display  Figure S3 in the Supporting Information.

Journal of Chemical Theory and Computation
Article the leading excitations defining the EA/IP-EOM-CCSD IP and EA R(k) operators for (TiO 2 ) 1−5 (Table S3 in the Supporting Information shows the same information for larger clusters). For IPs the r H (= ̅ r H ), r H−1 (= − r H 1 ), and so on singly excited amplitudes correspond to the removal of an α (β) electron from HOMO and HOMO−1 orbitals, respectively. A similar convention has been assumed to characterize the EA states. It is interesting to notice that, for the TiO 2 molecule, where the IP (cationic) and EA (anionic) electronic states are dominated by determinants obtained from the N electron HF Slater determinant by removing/adding one electron from/to HOMO/LUMO orbitals, for larger (TiO 2 ) n clusters the structure of these states becomes increasingly more multiconfigurational. This departure from the single determinantal picture implied by Koopmans theorem is best illustrated by the example of the (TiO 2 ) 5 cluster where five orbital (or 10 spin− orbital) excitations contribute to the N − 1 electron coupledcluster wave function with weight factors greater than 0.1. Something similar was previously observed by us for the excited states of small TiO 2 clusters. 40 Hydroxylation. In the presence of water, the surface of TiO 2 particles will be (partly) hydroxylated. Water molecules will react with undercoordinated titanium and oxygen atoms on the surface to form surface hydroxyls. We studied the effect of hydroxylation by comparing results for the hydroxylated monomer (TiO 2 )(H 2 O) 2 , orthotitanic acid (Ti(OH) 4 ), with those for (TiO 2 ), as well as those of the dimer global minimum candidate structure (the trans dimer cluster) and its hydroxylated counterpart. As can be seen from a comparison of Tables 2, 3, S2, and 5, all methods predict that for both clusters' hydroxylation shifts the IP and −EA values downward and upward, respectively, and increases the size of the fundamental gap. However, the methods disagree on the exact nature and magnitude of the changes. Specifically, ΔB3LYP predicts only a small change in −IP and a large change in −EA for both clusters, while all other methods, including the use of the KS orbital energies, yield significant changes in both −IP and −EA for the monomer, and in the case of G 0 W 0 and EA-IP-EOM-CCSD significant changes in −IP and small changes in −EA for the dimer. Just as for the nonhydroxylated clusters the KS-B3LYP results again lie far away from the GW and CC results, while for ΔB3LYP the match with GW and CC results improves.
Exciton Binding Energy. When making the assumption that for all clusters in our study the lowest excited state is bright, i.e., not symmetry forbidden, we can then use (EA/IP)-EOM-CCSD to calculate their fundamental and optical gaps in a consistent approach and extract information about the exciton binding energy in these clusters. The result of such a comparison is shown in Table 6, which shows, as expected, that the fundamental gap is always predicted to be considerably larger than the optical gap. The exciton binding energy appears to decrease with increasing cluster size but not in a monotonous fashion.
The hydroxylated clusters, while having fundamental and optical gaps that are much larger than their nonhydroxylated

Journal of Chemical Theory and Computation
Article counterparts, have exciton binding energies that are comparable to the other clusters. While hydroxylation thus opens up the fundamental and optical gaps, it does not appear to significantly influence the screening of the excited electron−hole interaction and thus the exciton binding energy. The exciton binding energy values of both the bare and hydrated clusters, are approximately 3 orders of magnitude larger than that measured for the bulk. 102 Similar dramatic increases in the exciton binding energy, when going from the bulk to nanoscale, have been predicted 103−106 and experimentally observed 107 for other systems and are probably due to the effect of a combination of increased overlap between hole and electron and reduced dielectric screening of the interaction between them in such small systems.
Finally, Table S4 in the Supporting Information provides estimates from ΔCCSD of the excitation energy toward the lowest triplet exciton for selected clusters. As expected, and in line with Hund's rule, comparing these ΔCCSD estimates with the EOM-CCSD optical gap values in Table 6 shows that the lowest energy triplet exciton generally lies lower in energy than the lowest singlet exciton, on the order of hundreds of millielectronvolts and that hence the lowest energy triplet exciton typically is more strongly bound than its singlet counterpart.
Adiabatic Potentials. We approximate adiabatic potentials and vertical potentials for the relaxed anion and cation, respectively, by optimizing the anion and cation geometries, respectively, using a setup similar to that of the ground state geometries (i.e., B3LYP/def2-TZVP) and performing singlepoint ΔB3LYP/def2-TZVPP, G 0 W 0 /def2-TZVPP, and EA/IP-EOM-CCSD/def2-TZVPP calculations on the relaxed anion and cation geometries.Tables 7 and 8 below give vertical −EA values for the anions and vertical −IP values for the cations, as well as adiabatic −EA and −IP values for the TiO 2 molecule and the (TiO 2 ) 2 and (TiO 2 ) 3 clusters. For (TiO 2 ) 2 we consider the three different possible isomers, which is critical when comparing to experiment; see below. In all cases, we have employed no symmetry constraints and verified that the relaxed anion and cation geometries correspond to proper minima and where required have distorted the geometries along imaginary modes and reoptimized the geometries until minima were obtained. The cationic version of the club isomer of (TiO 2 ) 2 spontaneously interconverts into the trans isomer during optimization. All other clusters retain their approximate topology after optimization as anion or cation.
First, we compare our predicted −EA values to previous computational results from the literature. The G 0 W 0 vertical −EA values for the different (TiO 2 ) 2 isomers in Table 7 are similar to the corresponding G 0 W 0 values by Marom et al., 38 even if the match is worse than for the neutral cluster geometries. The EA/IP-EOM-CCSD adiabatic −EA values in Table 8 can be compared with the recent ΔCCSD(T) results of Dixon and co-workers. 58 Doing so we find an identical energy ordering of the predicted adiabatic −EA values for both approaches and a reasonable match between the explicit values (see the Supporting Information for a comparison of Dixon's and our ΔCCSD(T) results). In both cases some of the differences could be attributable to a different description of the anion potential energy surface and hence slightly different anion minimum energy geometries rather than the method used to calculate the electron affinity.   The G 0 W 0 and EA/IP-EOM-CCSD predictions of the vertical −EA and −IP for the optimized anionic and cationic clusters geometries in Table 7 display, as expected from the experience with both methods for neutral cluster geometries, similar trends, bar a minor disagreement about the relative −EA of the cis and club isomers. The observation that IP-EOM-CCSD consistently predicts deeper IP values than G 0 W 0 is also reproduced for the optimized cationic cluster geometries. More interestingly, ΔB3LYP appears to perform much better for the optimized anionic and cationic clusters geometries than for their neutral counterparts.
Comparison with vertical and adiabatic −EA values measured for size-selected anionic clusters by Wang's group 17,18 using a combination of mass spectrometry and photoelectron spectroscopy requires, in the case of the larger clusters, external input regarding which of the possible isomers is present. We follow here the assumption by Marom and co-workers 38 that the experimental cluster generation method employed preferentially generates the isomer with the most positive vertical EA (most negative vertical −EA) value for the corresponding neutral cluster. Based on both their and our results (see Tables  2 and 3) the relevant structure in the case of (TiO 2 ) 2 is then the club isomer, while for (TiO 2 ) 3 they predict that the neutral global minimum candidate structure has the most positive vertical EA. Focusing on these cluster structures that are most likely to accept an electron, we find a good fit between the experimentally measured vertical and adiabatic −EA values and the predictions by ΔB3LYP, G 0 W 0 , and EA/IP-EOM-CCSD.
Not surprisingly, relaxation results in deeper, i.e., more negative, vertical −EA values and more shallow, i.e., less negative, vertical −IP values. More interesting, G 0 W 0 and EA/ IP-EOM-CCSD predict very similar changes when comparing vertical potentials for the neutral and relaxed anion/cation structures. Also both for −EA and −IP the changes are roughly 2−8 times as large for the (TiO 2 ) 2 isomers than for the TiO 2 molecule, while the change in −IP values is 2−5 times as large as that in −EA. Holes thus appear to be considerably more strongly trapped than electrons, while it seems that larger clusters allow for more relaxation. The differences between the vertical −EA and −IP values of the neutral clusters and their adiabatic counterparts is considerably smaller and for selected cases even positive (−EA, TiO 2 molecule) or negative (−IP, cis isomer of (TiO 2 ) 2 ), respectively. It thus appears that any stabilization of the anionic/cationic state upon relaxation is largely compensated for by a destabilization of the neutral state. ΔB3LYP calculations yield a similar picture other than that holes are predicted to trap much stronger. As a result ΔB3LYP predicts, in contrast to ΔG 0 W 0 and IP-EOM-CCSD, adiabatic −IP values that are considerably shallower than their vertical counterparts for neutral structures.
Microscopic Picture. Considering the leading excitation contributions to the EA/IP-EOM-CCSD vertical −IP and −EA values in Table 4 for the different clusters and their HF orbitals (see Figures S4−S6 in the Supporting Information), it is clear that in all cases the anion involves the localization of an excess electron on the most exposed titanium atom with the smallest number of coordinating oxygen atoms. In some cases, multiple titanium atoms are involved, e.g., in the (TiO 2 ) 2 trans and (TiO 2 ) 4 global minimum candidate structures, but generally these atoms are symmetry equivalent. In all cases the relevant titanium atoms are those with the least negative onsite electrostatic potential (calculated using the formal ionic charges, equivalent of the Madelung potential in crystalline solids) suggesting a possible correlation between −VEA and cluster electrostatics. Such a localization makes sense, as simple electrostatics would suggest that the energetic penalty of localizing an excess electron on an atom should decrease when the onsite electostatic potential becomes less negative. Finally, as previously observed by Marom and co-workers, 38,55 the environment of these exposed titanium atoms, on which the excess electron localizes, can be reminiscent of that of Ti 3+ sites on TiO 2 surfaces. However, we should stress that all atoms in all clusters studied here have similar charges and there is no evidence for charge transfer. Indeed this is where semi-ionic materials such as TiO 2 108 differ from a semicovalent materials such as SiO 2 , where there can be charge transfer in the ground state of nanoclusters, e.g., so-called valence alternation pairs. 109 The case of the cations, and hence the −IP values of the clusters, is more complicated. First, the relevant orbitals (e.g., HOMO) are often very delocalized, involving multiple distinct types of oxygen atoms: for example, both the one-and twocoordinated oxygen atoms in the (TiO 2 ) 2 trans global minimum candidate structure; second, while a similar electrostatic argument can be made as for the anion that the hole should localize on the oxygen atom with the least positive onsite electrostatic potential, typically, a terminal 1-coordinated oxygen atom. There are also cases, e.g., the (TiO 2 ) 2 club structure, where the hole localizes on other oxygen atoms instead (the three two-coordinated oxygen atoms in the case of the club structure). The link between the −VIP value and the onsite electrostatic potential appears thus weaker than for −VEA. Visual inspection suggests that an additional contributing factor might be the number of atoms over which relevant orbitals are delocalized. The especially strong multiconfigurational character of the IP for the (TiO 2 ) 5 global minimum candidate structure, discussed above, appears linked to the existence of two near degenerate orbitals centered around the one-coordinated oxygen atoms on opposite sides of the cluster, which are not symmetry equivalent.
The effect of hydroxylation on −IP, −EA, and the fundamental gap can at least in part be explained by changes in the electrostatic environment of atoms. Hydroxylation increases the coordination number of undercoordinated atoms and increases the onsite electrostatic potential for such atoms and destabilizes the localization of electrons or holes on them. As previously observed in terms of the optical gap of TiO 2 nanoclusters, 44,45 hydroxylation averages out the electronic structure of clusters by "normalizing" the electrostatic environment of atoms.
It is tempting to blame the apparent failure of ΔB3LYP to predict exactly how −IP and −EA values vary between different bare TiO 2 clusters on the multiconfigurational nature of the cationic and anionic states of the clusters discussed above. The problem, however, with this explanation is that the use of the B3LYP KS orbital energies results in the successful recovery of the trends found for the many-body methods. Likely, the problem is more subtle and related, at least in part, to the localization of the hole and excess electron in the cation/anion and how this is described by DFT. This would explain why ΔB3LYP appears to perform better for relaxed anionic and cationic structures, in which by necessity the excess electron/ hole is unambiguously more localized. It also would be in line with the observation, discussed above, of the SCF for the cations of selected clusters converging to solutions with different spin distributions when changing the basis set and this resulting in substantially different predictions for −IP.
Finally, the notion that a good description of (de)localization is important is also supported by the fact that ΔDFT calculations using the PBE XC functional, which tends to give a more delocalized description, predict a different trend than that obtained with ΔB3LYP (see Figures S1−S3 in the Supporting Information).

■ CONCLUSIONS
We studied the ionization potentials, electron affinities, fundamental gaps, and exciton binding energies of small TiO 2 nanoclusters using a range of methods. We demonstrate that all many-body methods considered predict similar trends, but that G 0 W 0 appears rather sensitive to basis-set size. Moreover, we find that the difference between qsGW and G 0 W 0 predictions is generally larger than those between either of the GW approaches and EA/IP-EOM-CCSD. Specifically, in the case of −IP values, G 0 W 0 predicts the shallowest vertical potentials and qsGW the deepest, with the vertical −IP predicted by the EA/IP-EOM-CCSD method lying between the two but closer to the latter. Similarly, G 0 W 0 predicts the smallest fundamental gap values, qsGW the largest, and EA/IP-EOM-CCSD values roughly intermediate between both Green's function methods. For the vertical −EA values, the absolute difference in prediction between the three methods is smaller though not necessarily in relative terms. The B3LYP highest occupied and lowest unoccupied orbitals follow the same trends as those predicted by the many-body methods for −IP and −EA but are generally far too shallow and deep, respectively, in absolute terms. The total energy ΔB3LYP method, in contrast, yields values in the correct energy window but predicts different trends than the many-body methods in all electronic properties with cluster size. Assuming that the qsGW and IP/EA-EOM-CCSD methods are inherently the most accurate approaches among those studied, our results suggest that while ΔB3LYP allows one to predict the correct orders of magnitude of vertical −IP and −EA for bare TiO 2 clusters, the prediction of accurate trends in vertical −IP and −EA values appears to require a more advanced approach than ΔB3LYP and probably by extension ΔDFT in general. Conversely, G 0 W 0 , at least when using B3LYP orbitals as a starting point, allows one to accurately predict trends in −IP and −EA, but appears less reliable in terms of absolute values, especially in the case of −IP, although we cannot rule out this being in part due to different convergence rates with respect to basis-set size. Calculations of vertical potential values for the relaxed cation/ anion cluster geometries, as well as adiabatic potentials, generally display the same performance as the vertical potentials for the neutral cluster, other than an improved performance of ΔB3LYP, and in the case of −EA a good fit to experiment.
All electronic properties considered in this work are shown to vary nonmonotonically with cluster size while hydroxylation of undercoordinated atoms is demonstrated to be linked to an opening up of the fundamental gap and thus −IP moving deeper and −EA slightly shallower. We suggest that the change with cluster size, as well as upon hydroxylation, can be understood, at least in part, in terms of the on-site electrostatic potential and the number of atoms states are delocalized over. The exciton binding energy in these small clusters appears to decrease with increasing cluster size but again in a nonmonotonic fashion, while hydroxylation, in contrast, results in no real change in the exciton binding energy relative to the unhydroxylated bare clusters. Structural relaxation after removal/addition of an electron, i.e., trapping, is found to be consistently stronger for holes than electrons. The net effect on the predicted adiabatic potentials is, however, found to be relatively small because the stabilization of the cationic/anionic state is largely compensated for by a destabilization of the neutral state.