Quasiparticle and Optical Properties of Rutile and Anatase TiO$_{2}$

Abstract

Quasiparticle excitation energies and optical properties of TiO2 in the rutile and anatase structures
are calculated using many-body perturbation theory methods.
Calculations are performed for a frozen crystal lattice; electron-phonon coupling is not explicitly considered.
In the GW method, several approximations are compared and it is found that inclusion of the full frequency
dependence as well as explicit treatment of the Ti semicore states are
essential for accurate calculation of the quasiparticle energy band gap.
The calculated quasiparticle energies are in good
agreement with available photoemission and inverse photoemission experiments.
The results of the GW calculations,
together with the calculated static screened Coulomb interaction,
are utilized in the Bethe-Salpeter equation to calculate
the dielectric function ϵ2(ω) for both the rutile and anatase structures.
The results are in good agreement with experimental observations,
particularly the onset of the main absorption features around 4 eV.
For comparison to low temperature optical absorption measurements that
resolve individual excitonic transitions in rutile, the low-lying discrete excitonic energy levels
are calculated with electronic screening only. The lowest energy exciton found in the energy
gap of rutile has a binding energy of 0.13 eV. In agreement with experiment,
it is not dipole allowed, but the calculated exciton energy
exceeds that measured in absorption experiments
by about 0.22 eV and the scale of the exciton binding energy is also too large.
The quasiparticle energy alignment of rutile is
calculated for non-polar (110) surfaces.
In the GW approximation, the valence band maximum is 7.8 eV
below the vacuum level, showing a small shift
from density functional theory results.

pacs:

Even after half a century of researchGrant (1959); Diebold (2003),
investigation of the fundamental properties of titanium dioxide (TiO2) crystal phases
remains important and fruitful, in part due to the
role they have in concepts to effectively utilize solar energy.
For example, TiO2 structures form the photoactive component in
heterogeneous photo-catalysts which, by absorbing energy
from the sunlight, degrade environmentally hazard materials Fox and Dulay (1993); Hoffmann et al. (1995)
and split water into H2 and O2Fujishima and Honda (1972).
Scintered anatase TiO2 nanoparticles provide the backbone for electron transport and
the substrate for organic chromophores in
the Grätzel photovoltaic solar cells Oregan and Gratzel (1991). In addition to that,
TiO2 has been widely used in various areas from optical
coatings to pigments Diebold (2003).
Fundamental to all of these applications are the relative alignments
of essential energy levels near the valence and conduction band edges
of TiO2 crystal phases and the corresponding optical transition energies.
If predictive computational methods are going to have impact
on the understanding and design of heterogeneous photocatalytic systems based on TiO2,
we must first establish that these methods can predict
the basic properties of the crystal phases, providing a coherent
framework for all the experimental facts.

Rutile and anatase are two common crystal structures
in which TiO2 is found.
In both phases, each Ti
atom in the crystal is surrounded by a slightly distorted octahedron
formed by six oxygen atoms. The distinct phases exhibit a different connection between
the distorted octahedra (TiO6). In the rutile phase
each octahedron shares two edges with its neighbors, while in the
anatase phase each octahedron shares four Pauling (1929).
In the rutile form, the crystal has a simple-tetragonal structure Abrahams and Bernstein (1971)
with a = b = 0.45936 nm and c = 0.29587 nm. The symmetry of the lattice
is described by the space group P42/mnm with the only internal
parameter u = 0.30479. In the anatase form the crystal structure
is body-centered tetragonal Horn et al. (1972); Fahmi et al. (1993) and
belongs to space group I41/amd. The three sides of the conventional
cell are a = b = 0.3784 nm and c = 0.9515 nm respectively. The internal
parameter u is 0.208 Horn et al. (1972); Fahmi et al. (1993).
The measurements quoted were done at room temperature;
the change in lattice parameters upon reducing the temperature to 15 K is
less than 0.001 nm Burdett et al. (1987).

Most of the early first-principles calculations of the properties of TiO2
were based on the local density approximation (LDA) in a Density Functional
Theory (DFT) based approach Hohenberg and Kohn (1964); Kohn and Sham (1965).
The crystal structures and ground-state
properties were accurately reproduced Glassford and
Chelikowsky (1992a); Mo and Ching (1995); Mikami et al. (2000); Asahi et al. (2000).
However, as has been more generally observed for semiconductors and insulators,
the energy gaps pertaining to optical properties were found to be too small.
The minimum energy gap in the LDA band structure underestimated
the band gap observed
in optical experiments Pascual et al. (1977); Tang et al. (1993)
by about 40% Mikami et al. (2000).
Calculations based on Hartree-Fock theory have been performed,
giving accurate structural properties for rutile and anatase, but with
a minimum energy gap that exceeded 10 eV Fahmi et al. (1993).
A hybrid approach, admixing a fraction of the bare exchange from Hartree-Fock,
also showed accurate structural properties for rutile Zhang et al. (2005),
with a band gap that is closer to experiment (3.4 eV) Muscat et al. (2001).
A more economical approach, approximately accounting for explicit
Coulomb interactions through a U parameter acting on the Ti 3d electrons
in a DFT+U approach,
overestimated the lattice parameters while still showing a band gap that
was smaller than experiment Nolan et al. (2008).

A direct approach to calculate
electronic excitation energies based on the Green’s function approach of many-body
perturbation theory (MBPT), specifically utilizing
the GW approximation for the electron self energy Hedin (1965),
has proven to be relatively accurate for a broad array of materials Hybertsen and Louie (1986); Aryasetiawan and Gunnarsson (1998); Aulbur et al. (1999); Onida et al. (2002).
Several calculations have been reported for TiO2
with different implementations of the GW method Oshikiri et al. (2003); van Schilfgaarde et al. (2006); Thulin and Guerra (2008); Mowbray et al. (2009).
However, the band gap was significantly overestimated in all
these reports.

The frequency dependent macroscopic dielectric function
probed in optical measurements has been extensively studied within the
framework of DFT for TiO2 crystal phases Glassford and
Chelikowsky (1992b, a); Mo and Ching (1995); Asahi et al. (2000).
Using the independent particle approach of
Ehrenreich and Cohen Ehrenreich and Cohen (1959),
the underestimate of the fundamental band gap immediately gives an error
in the optical threshold and
the overall shape of the dielectric function
calculated in this way was quite different from the experimental measurements
Cardona and Harbeke (1965); Hosaka et al. (1997).
In the framework of MBPT, the GW-based results for the electronic
excitation energies are input to a direct treatment of neutral excitations
through solution of the two-particle Bethe-Salpeter
equation (BSE), an approach that has provided
a satisfactory description of the optical properties of a number of systems Rohlfing and Louie (2000); Onida et al. (2002).
An application of the BSE approach for rutile and anatase TiO2
has only recently appeared in the literature Lawler et al. (2008).
The shape of the spectra are in much better agreement with experiment
as compared to the independent particle approach.

In this article, we critically assess the application of MBPT methods
to calculate the electronic and
optical excitations for TiO2 in the rutile and anatase crystal phases.
To treat the electronic excitation energies, we use the GW method
without self consistency or inclusion of vertex corrections.
Empirically, this approach is often relatively accurate,
although a full understanding of cancellations between self consistency
and vertex corrections remains
an open area of research Shirley (1996); Holm and von Barth (1998); Ku and Eguiluz (2002); Bruneval
et al. (2006a); Kotani et al. (2007); Shishkin et al. (2007); Shishkin and Kresse (2007).
Several different approximations to handle the frequency
dependence of the screened Coulomb interaction in the GW method
are compared. We find that use of plasmon pole models Hybertsen and Louie (1986); von der Linden and Horsch (1988)
results in significant overestimation of the band gap.
When the full frequency dependence of the screening is included,
together with explicit treatment of the Ti semicore states,
the calculated electronic excitation spectrum is found to
be in good agreement with photoemission and inverse photoemission spectra Tezuka et al. (1994); See and Bartynski (1994).
Interestingly, the calculated fundamental gap (3.34 eV and 3.56 eV for rutile and anatase respectively),
is still larger than the measured minimum gap from optical absorption
(3.03 eV for rutile Pascual et al. (1977) and estimated to be 3.3 eV for anatase Tang et al. (1993)).
As a first step towards application to heterogeneous photocatalytic systems,
we discuss the alignment of the valence and conduction band edges
at non-polar vacuum-solid (110) interface of TiO2.
We find that the GW method implies only very small corrections relative to
the LDA for the valence band
position, with the most of the band gap error going to shift the conduction band edge upwards.

To explore the role of electron-hole interactions and excitonic binding
energy, we have used the results from the GW based calculations
as input to the BSE approach.
Similar to the recent results of Lawler etal.Lawler et al. (2008),
the calculated frequency dependent dielectric function accurately
reproduces the main onset of absorption near 4 eV and gives
a good account of the frequency dependence for both rutile and anatase.
We also solve the BSE for the low-lying, bound exciton states for rutile.
The deepest exciton binding energy is calculated to be about 0.13 eV.
The dipole-forbidden character of the lowest exciton agrees
with low temperature measurements Pascual et al. (1977),
although the predicted exciton energy (3.25 eV) is still 0.22 eV
larger than experiments.
Also, the magnitude of the exciton binding energy is larger.
While the discrepancy for the exciton energy
could very reasonably be regarded as within the
expected errors of the MBPT methods used here, it may suggest
an important role for electron-phonon coupling
in screening and in further
renormalizing the energy gap in TiO2.

The rest of the article is organized as follows: In Sec. II, the
methodologies used in the DFT, electronic excitation and optical excitation calculations are
briefly summarized, the key approximations are discussed and
the numerical details are provided.
In Sec. III, we present the main results for the electronic
and optical excitations in rutile and anatase TiO2
and discuss them in comparison to available experiments.
Finally, we conclude the article
in Sec. IV with a short discussion, including the role of coupling to phonons.

ii.1 DFT Calculations

The LDA eigenvalues and eigenvectors of TiO2 are calculated
with a plane-wave basis set using norm-conserving pseudopotentials. Unless indicated
otherwise, the LDA calculations are carried out using the ABINIT packageGonze et al. (2002, 2005).
In TiO2 the Ti is nominally ionized to [Ti4+] and the
low lying conduction band states are of predominantly 3d character.
As we show below, artificially dividing the n=3 shell of Ti into frozen core (3s and 3p) and
valence (3d) contributions introduces a significant error to the energy band gap.
The pseudopotential of Ti which includes semicore electrons is generated
using the OPIUM package opi () in the Troullier-Martins scheme Troullier and Martins (1991)
with an initial configuration of (Ne)3s23p63d04s04p0.
The outermost five orbitals are included and the cutoff radii (in
Bohr) are 0.9, 0.9, 1.0, 0.9, and 0.9 respectively. Other
pseudopotentials are taken from the ABINIT pseudopotential database
generated using the FHI99PP package Fuchs and Scheffler (1999).

In all calculations,
the Perdew-Wang representation Perdew and Wang (1992) of Ceperly-Alder
exchange-correlation potential Ceperley and Alder (1980) is used.
When including the Ti semicore states,
a kinetic energy cutoff of 200 Ry is used to ensure the convergence
of the LDA results, as suggested by previous calculationsMikami et al. (2000); Vast et al. (2002).
To examine the accuracy of the pseudopotentials, we calculate the
optimized lattice constants for rutile,
finding a\leavevmode\nobreak= 4.5484 Å(4.5936 Å), c/a\leavevmode\nobreak= 0.6414 (0.64409)
and u\leavevmode\nobreak= 0.3040 (0.30479), agreeing with the experimental values
noted in parentheses Abrahams and Bernstein (1971) to the accuracy generally expected for LDA calculations.
We also compare our LDA calculations with results obtained using the VASP package Kresse and Furthmüller (1996, 1996)
with the recommended projector augmented-wave (PAW) pseudopotentials Blöchl (1994).
The difference between the two LDA calculations is within 0.5% for lattice parameters and less than
0.05 eV for bandgaps.
In the GW and BSE calculations described below,
the geometrical parameters of the unit cell for both rutile and anatase
phases are taken from
experimental measurementsAbrahams and Bernstein (1971); Cromer and Herrington (1955); Fahmi et al. (1993); Horn et al. (1972).

ii.2 GW Method

In MBPT, the evolution of the electrons in a material is described
by the one-particle Green’s function, with the effect of electron-electron interactions
represented by an electron self energy operator.
Well defined electronic excitations appear as peaks in the
corresponding spectral function.
Excitations with single particle character, namely quasiparticles, can be
obtained as solutions of a Schrodinger-like
equation Hedin and Lundqvist (1969)

(T+Vext+VH)ψn,k(r)+

∫dr′Σ(r,r′;En,k)ψn,k(r′)

=

En,kψn,k

(1)

where T is the kinetic energy, Vext is the external potential,
and VH is the average Hartree potential.
Σ is the self
energy of the electrons and the indices refer to Bloch states n, k. It includes all the exchange-correlation
effects contributed by surrounding electrons.
Since Σ is generally
non-Hermitian, En,k is complex with the real part giving the quasiparticle energy
and the imaginary part corresponding to the
width of the quasiparticle peak in the spectral function, i.e. the quasiparticle lifetime.

A practical approximation to calculate Σ
has proven to be the so-called GW approximation of Hedin Hedin (1965), in which the self energy Σ(r,r′;E)
is formally written as

Σ(r,r′;E)=

i2π∫dE′e−iδ+E′G(r,r′;E−E′)W(r,r′;E′),

(2)

Here G is the Green’s function of the electrons and W is the dynamically
screened Coulomb interaction determined by the inverse dielectric matrix ϵ−1(r,r′′;E), and δ+ is a positive infinitesimal time. The G and W in Eq. (2) refer
to the fully interacting Green’s function.
However, in practice,
using an initial LDA calculation to determine the screening through linear response calculations (not including the exchange-correlation
kernel) and to provide an initial, independent-particle Green’s function has often proven
to be sufficiently accurate.
There are specific examples where the LDA orbital character can be wrong, e.g. in some late 3d
transition metal compoundsAryasetiawan and Gunnarsson (1995); Faleev et al. (2004); Bruneval
et al. (2006b).
However, in TiO2, the Ti 3d is almost empty
and the valence band edge region is predominantly O 2p character
with minimal admixture of Ti 3d. The TiO2 case should be similar to the
vast majority of semiconductors and insulators in this regard. Also, the LDA wavefunctions are
sufficiently accurate
that a first order estimate of the self energy correction to the LDA eigenvalues is adequate.
The quasiparticle energy correction
ΔEn,k to a LDA orbital ϕn,k is obtained
through a reduced form of Eq. (1) as

ΔEn,k=Zn,k⟨ϕn,k|Σ(ELDAn,k)−VLDAxc|ϕn,k⟩,

(3)

where VLDAxc is the exchange-correlation potential and Zn,k
is the renormalization factor of the orbital defined as Zn,k=(1−∂Σ/∂E)−1|E=ELDAn,k.

The frequency dependence of the screened Coulomb interaction
(W) can often be addressed
using a generalized plasmon-pole (GPP) model Hybertsen and Louie (1986); von der Linden and Horsch (1988),
with substantial advantages in computational efficiency.
The GPP models have proven to be relatively accurate
for many semiconductors and insulators, including ionic crystals
such as LiCl Hybertsen and Louie (1985) and MgO Shirley (1998).
However, as discussed below, we find that use of the GPP
leads to a gap that is substantially too large for TiO2.
Several approaches to include the full frequency-dependent
dielectric matrix have been implemented and described
in the literature: (1) an analytical continuation method Rojas et al. (1995); Rieger et al. (1999),
(2) a direct method which carries out the integration in Eq. (2)
along the real axis Marini et al. (2002a, b); Shishkin and Kresse (2006), and
(3) a contour deformation (CD) method which deforms the integration
in Eq. (2) along the imaginary axis Lebègue et al. (2003).
We adopt the contour deformation method to carry out the calculations,
which is particularly efficient for evaluating self energy for
states near the gap region.

In the CD method, the correlation contribution Σc(r,r′;E)
of the self energy is written as the sum of two terms Lebègue et al. (2003); Bruneval (2005)

where Wp=W−Vcoul, Ef is the Fermi energy, η is
a small damping amplitude and Θ is a Heaviside function.
The first term in Eq. (II.2) comes from the integration along the
imaginary axis. As W is now smooth along the imaginary axis, a sparse
sampling of E is sufficient to converge the integration. The second
term is the residual contribution of poles near the real axis. It
is non-zero only while E>En,k>Ef or E<En,k<Ef.
For any E close to the Fermi surface, only Wp for |E−En,k|∼0
have non-vanished contributions to Σc in the second term. This makes the calculation
more computational efficient; W is a smooth function of
|E−En,k| around 0, due to the band gap,
and only relatively low frequencies need to be sampled.

We implement the contour deformation approach
based on a private branch of the YAMBO package Marini et al. (2009).
The integration along the imaginary axis in Eq. (II.2) is performed
with a non-uniform mesh of N points according to

E′′i=E0tan(i−12Nπ),i=1,2,...,N,

(5)

which maps the integration along the imaginary axis to an integration
on the [0, 1) interval. The energy E0 provides a scale for the overall
density
of the samples on the imaginary axis. Half of the mesh spans the energy scale
from zero up to E0
while the other half sample the higher energies. For TiO2, a mesh
of 50 points and an energy scale of 1 Ry were enough to keep the numerical
error of the integration within 1 meV. Wp on the real axis is
uniformly sampled with an energy increment of 0.1 eV and values between
mesh points are linearly interpolated.
The special case
in Eq. (II.2) for E→En,k must be handled with care.
A consistent treatment, that avoids the apparent divergence
and properly handles all the terms in Eq. (II.2),
is to add a small positive energy to E−En,k
(say δ=2.0×10−7 Ry) when necessary.
A very similar contour deformation approach has been implemented
in ABINIT Bruneval (2005), and we have carefully compared the
results for test cases. The ABINIT package has also been used for the
GPP model calculations.

For all the GW calculations, the energy cutoff
is 60 Ry for the evaluation of the bare Coulomb exchange contribution
Σx , and 20 Ry for the correlation part Σc.
A total of 160 bands are used for the calculation of both dielectric matrices
and self-energies. An unshifted 4×4×6 Monkhost-Pack (MP)
meshMonkhorst and Pack (1976) is used to sample the Brillouin zone (BZ)
of rutile, while for anatase an unshifted 4×4×4 MP mesh
is used.

A test of the convergence with respect to the number of bands included
is shown in Fig. 1 for the final full-frequency approach with
Ti semicore electrons treated explicitly as valence electrons in the pseudopotential.
In order to characterize the fully converged values, the data were
fit with two different empirical forms, E(N)=E0−b/N, and E(N)=E0−b⋅exp(−N/c).
We first check the validity of the fitting forms for the case of bulk silicon.
The exponential form closely represents the band edge and band gap
energies as a function of the number
of included conduction bands, yielding extrapolated results
in excellent agreement with those obtained via methods suggested by
Bruneval and GonzeBruneval and Gonze (2008). For rutile, the fitting
curves displayed in Fig. 1 are indistinguishable, but predict
slightly different N→∞ results for the absolute
shift of the valence band edge, as indicated by the horizontal dashed
lines in Fig. 1(b). In particular, the fit for the quasiparticle
energy gap indicates a converged quasiparticle energy of 3.37
eV for N=160. For the absolute shift of the valence band edge, the
convergence is somewhat slower, with extrapolated values of -0.12 eV and -0.31 eV. This suggests that
the valence band edge in the final results is 0.2 to 0.4 eV lower than the N=160 value.

Figure 1: Quasiparticle direct energy gap at Γ (a)
and energy of VBM (b) for rutile as a function
of the total number of bands kept in the full frequency dependent GW calculation.
Squares represent calculations with Ti semi-core electrons included explicitly
as valence electrons. The results are
fitted using two different functional forms described in the text
and displayed as solid lines that are indistinguishable in the figure.
However, in (b), they have different N→∞ limits,
displayed as dash-dot (1/N) and long dash (exp(−N/c)) lines in
the figure respectively.

ii.3 Bethe-Salpeter Equation and Optical Properties

A detailed description of the
BSE method is given in the literatureOnida et al. (2002); Rohlfing and Louie (2000).
We use the implementation in the public branch of the YAMBO package Marini et al. (2009).
With the Tamm-Dancoff approximation Fetter and Walecka (1971) and the use of a
static screened Coulomb interaction, the BSE is simplified to a generalized eigenvalue equation
Hvc,v′c′Av′c′s=EsAvcs.
The effective Hamiltonian
Hvc,v′c′ has an explicit form of

Hvc,v′c′=(Ec−Ev)δvv′δcc′+2¯Vvc,v′c′−¯Wvc,v′c′,

(6)

where the quasiparticle energies (taken from the GW calculations) enter on the diagonal,
¯Vvc,v′c′=⟨vc|VCoul|v′c′⟩,
¯Wvc,v′c′=⟨vv′|W(E=0)|cc′⟩
and the indices v,\leavevmode\nobreakc refer to the occupied valence and empty conduction
band states. For brevity, the explicit reference to Bloch
wavevector k for each state is suppressed.
The effective Hamiltonian
in Eq. (6) is explicitly written for the spin singlet excitations
which are directly probed in optical measurements.
For the spin triplet excitations, the effective Hamiltonian is modified
by dropping the so-called exchange term 2¯Vvc,v′c′.
In terms of Es and Avcs, the macroscopic dielectric
function ϵM(ω), including local field effects Rohlfing and Louie (2000); Onida et al. (2002), is expressed as

ϵM(ω)=1−limq→04πe2q2∑s∣∣∣∑vc⟨v∣∣e−iq⋅r∣∣c⟩Avcs∣∣∣2ω−Es+iη.

(7)

In YAMBO, we use the option to evaluate the response function recursively Benedict and Shirley (1999); Marini et al. (2009).
In order to study specific, low energy exciton states, we also
directly diagonalize the generalized eigenvalue equation.

Calculations of optical properties via BSE are more expensive computationally.
For both phases, the static dielectric matrices are calculated with
80 bands and a damping coefficient of 0.1 eV.
The dielectric function ϵ(ω) is calculated on an 8×8×12
MP mesh for rutile and on an 12×12×12 MP mesh for anatase.
For both cases, electron-hole
pairs within 15 eV are taken in to build up the BSE kernel, in which
the energy cutoff is 10 Ry for ¯V and 3.5 Ry for ¯W.
To calculate the excitonic
binding energy of rutile, we restrict the basis set
for the effective Hamiltonian to one conduction band and one valence band.
The exciton binding energy converges relatively slowly with BZ sampling.
The final results are reported based on a 12×12×18 MP mesh to
sample the BZ.
The energy cutoff is larger in this calculation for a more accurate representation of the BSE kernel, 14 Ry for ¯V and 6
Ry for ¯W.

Method

SC

VBM

CBM

EΓg

in PP

Vxc

Σx

Σc

Z

Vxc

Σx

Σc

Z

LDA

GW

GPP(HL)

n

-19.97

-24.59

3.94

0.82

-11.48

-3.57

-5.69

0.84

2.08

4.48

GPP(VDLH)

n

-19.97

-24.59

4.37

0.82

-11.48

-3.57

-6.35

0.85

2.08

3.60

FF(CD)

n

-19.97

-24.59

4.50

0.71

-11.47

-3.57

-5.78

0.73

2.08

3.73

GPP(HL)

y

-20.21

-24.67

4.09

0.83

-20.37

-12.13

-5.64

0.85

1.75

4.27

GPP(VDLH)

y

-20.21

-24.67

4.41

0.83

-20.37

-12.13

-5.99

0.85

1.75

3.70

FF(CD)

y

-20.21

-24.67

4.56

0.70

-20.37

-12.12

-5.88

0.72

1.75

3.38

Table 1: Analysis of the valence band maximum (VBM), the conduction band minimum (CBM)
and the direct energy gap EΓg
at the Γ point of the Brillouin zone calculated for rutile TiO2
using the GW method with several different approximations.
The methods refer to two different generalized plasmon-pole
(GPP) models and the full frequency-dependent (FF) approach described in the text.
The second column indicates whether Ti semicore states are explicitly included
as valence electrons.
For the VBM and CBM, the expectation value is shown for the
exchange-correlation potential Vxc in the LDA,
the bare exchange self energy (Σx),
the correlation part of the self energy (Σc)
and the renormalization factor (Z).
The band gap is shown in the LDA and for the GW method for each case.
Energies are given in eV and the renormalization factor Z is dimensionless.

iii.1 Electronic Excitation Energies in TiO2

The calculated electronic excitation energies in titanium oxides
are found to be sensitive to technical factors in the GW calculations.
We illustrate that here for the case of TiO2 in the rutile phase (Table 1).
First, explicit, self consistent treatment of the semicore electrons (3s and 3p) on the Ti
in the calculations for the solid affects the calculated energy gaps.
As discussed in the literature, although the 3s and 3p levels are well separated
from the 3d states energetically, there is significant spatial overlap Marini et al. (2002b).
The effect for TiO2
is already evident at the LDA level, where the gap is reduced by more than 0.3 eV
upon explicit inclusion of the semicore electrons
relative to freezing the semicore electrons in the pseudopotential.

Depending on the approximation used to treat the electron self energy
in the GW method, the net influence of the semicore electrons varies.
Focusing on the influence of the semicore electrons for the most accurate,
full-frequency method (FF), the change is greater than 0.3 eV.
A more detailed view of the contributions of
the LDA exchange-correlation potential Vxc, the exchange Σx
and the correlation Σc part of the self energy are also displayed
Table 1.
The valence band maximum (VBM) orbitals are largely composed of
oxygen p-states.
They are spatially separated from the Ti semicore electrons, so the changes are relatively small.
Explicit treatment of the semicore electrons as valence electrons only reduces
Vxc by 0.24 eV. The reduction in Vxc is partially compensated
by the decrease of Σx (0.08 eV) and Σc
(0.04 – 0.15 eV). Consequently the overall change of ΔEVBM
is less than 0.1 eV.
On the other hand, the conduction band minimum (CBM) orbital is
almost purely of Ti 3d character with substantial overlap to the semicore electrons.
The expectation value of (Vxc)
changes by about 9 eV. Of course, there is a corresponding, large change in the
pseudopotential between these two cases.
The bare exchange term in the electron self energy operator changes by a similar amount.
Their combined
contribution to ΔECBM is only about -0.34 eV.
The changes in the correlation part of the self energy Σc
of the CBM orbital is sensitive to the GW method. For the full-frequency
method, Σc is decreased by only 0.1 eV when semicore electrons
are included in the PP, so we find that most of the net effect on the band gap
comes from the difference between LDA and bare exchange interactions
with the semicore electrons.

In Table 1 the results of using different methods
to address the frequency dependence of W are shown.
These affect the correlation part of the electron self energy
Σc and the renormalization factor Z.
The results obtained with the full-frequency dependent
dielectric function, evaluated using the contour deformation method
are the reference results, designated FF(CD) in the table.
For comparison, results from two different generalized plasmon pole models are shown.
In the Hybertsen-Louie approachHybertsen and Louie (1986), designated GPP(HL), sum rules are applied
to each individual dielectric matrix element to develop a plasmon-pole model for its
frequency dependence.
As shown in the lower part of Table 1, when the semicore electrons of Ti were explicitly included in the pseudopotential, the calculated energy gap is
almost 1 eV too large.
In the approach of von der Linden and Horschvon der Linden and Horsch (1988), designated GPP(VDLH),
each dielectric matrix is first transformed to the basis of eigenpotentials
and then sum rules are applied to develop a plasmon pole model for each eigenpotential.
This model also overestimates the energy gap, but by a smaller amount.
The renormalization factor Z, in addition to entering the evaluation of the quasiparticle energies in Eq. (3),
indicates the degree of dynamical admixture of collective excitations into the quasiparticle.
Larger deviations from unity indicate more admixture.
Compared with the FF results, the value of Z indicates that the degree of admixture
is substantially underestimated by both GPP models.

We have also tested the GPP(HL) approximation for
anatase TiO2 as well as two other titanates,
SrTiO3 and BaTiO3.
In all three cases, the GPP(HL) approximation leads
to minimum band gaps that are too large by 0.7-0.8 eV.
A similar deviation for the renormalization constant, Z, is also observed.
These results suggest that the quantitative issues with the plasmon pole
model extend more generally to titanates.
In previous calculations of the loss function Vast et al. (2002)
and the finite wavevector dynamical scattering factor Gurtubay et al. (2004)
for rutile TiO2,
substantial structure is seen in the frequency dependence, well beyond
what could be easily accounted for by a single pole model.
These effects trace to an interplay between strong local field effects
and the Ti semicore p- to empty d-shell excitation.
However, further analysis of the frequency dependence of the
screening in Si and LiCl shows that deviations from a pole
model for a range of frequencies around the plasmon energy is not
sufficient to predict the performance of the GPP model as it is used
to evaluate the GW expression.
The dynamical screening at larger frequencies only enters in an integrated fashion,
resulting in substantial cancellations internally.
In the case of titanates, we find that the full frequency dependence is essential
for quantitative results.
Similar conclusions were drawn for the case of metallic Cu Marini et al. (2002b).

Rutile

Anatase

K-points

Γ

M

R

Γ

X

EValLDA

0.00

-1.06

-1.04

-0.48

-0.05

EValGW

0.00

-1.15

-1.12

-0.58

-0.06

ECondLDA

1.75

1.80

1.78

1.96

3.22

ECondGW

3.38

3.40

3.34

3.56

4.89

EDirectgap,GW

3.38

4.55

4.45

4.14

4.95

EIndirectgap,GW

3.34 (Γ→R)

3.56 (Δ→Γ)

Table 2: LDA and quasiparticle energy levels of rutile and anatase near the
Fermi surface at selected k-points of high symmetry. Corresponding
quasiparticle energy gaps are also displayed. The energy reference
is taken to be the valence band maximum for both LDA and GW results.
All energies are
in eV.

The quasiparticle energy levels of rutile and anatase from the highest valence band
and the lowest conduction band
at selected high symmetry k-points are displayed in Table 2.
While the energy gap of rutile is found to be a direct gap (at Γ)
in the LDA, our FF GW calculations indicate it as indirect (Γ→R).
However, the energy difference between the direct and indirect gap
is small. The energy of CBM at Γ is only 0.04 eV higher than
the energy at R.

The calculated quasiparticle energies can be directly compared
to spectroscopic measurements for electron removal or addition to the solid.
The calculated value of the quasiparticle energy gap, 3.34
eV, agrees well with electron spectroscopy measurements, photo-emission and
inverse photo-emission measurements Tezuka et al. (1994); See and Bartynski (1994).
In Fig. 2, the density
of states (DOS) of rutile derived from FF GW calculations is plotted
together with the experimental spectra measured
using x-ray photoemission and bremsstrahlung isochromat spectroscopyTezuka et al. (1994), which yielded a band gap of 3.3 ± 0.5 eV.
Overall, the shape of the
calculated DOS matches the shape of the experimental spectra, especially
around the band gap.
The experimental spectra measured using ultraviolet photoemission and inverse photoemission spectroscopy See and Bartynski (1994)
show very similar results with the inferred minimum energy gap about 0.2 eV smaller.

The LDA calculations show that the top of the valence band of anatase
lies in the Δ direction, somewhere about 0.88 times of the
distance from Γ to X.
The energy at this k-point is 0.05
eV higher than the energy of the VBM at X. Subtracting the energy difference
as a perturbation from the quasiparticle energy gap between X and
Γ, which is 3.62 eV from Table 2, the quasiparticle
energy gap of anatase is found to be 3.56 eV.
The photoemission data available for anatase show an overall occupied
band with the oxygen p-bands Sanjinés et al. (1994). That is similar to rutile.
To our knowledge, there is no inverse photoemission data available for anatase,
so the calculations can not be compared to a direct measurement of the quasiparticle energy gap.

iii.2 Optical Excitation Energies in TiO2

More precise measurements of the minimum energy gap
rely on optical absorption.
This introduces the extra complication that the
observed threshold for absorption will be altered
by interactions between the photoexcited electron and hole,
the formation of bound exciton states.
For rutile, the BSE calculation shows
a series of singlet bound excitonic states at Γ. The lowest
two show s-state symmetry in the electron-hole envelope. They have a binding energy
of 0.13 eV and 0.06 eV respectively. They are not dipole-allowed. The
third and forth are degenerate, with electron-hole envelope showing px,y symmetry,
and have a binding energy of 26 meV. They are weakly dipole allowed for
the electric vector perpendicular to the c-axis.
Together with the calculated direct quasiparticle gap energy from above,
we obtain the lowest energy singlet exciton at Γ with energy 3.25 eV
and the first dipole allowed singlet exciton with energy of 3.35 eV.
High resolution, low temperature optical absorption measurements
for rutile resolve several separate features Pascual et al. (1977).
A very weak, but sharp exciton feature at 3.031 eV is identified as
the 1s exciton which is electric quadrapole allowed. A stronger,
but still relatively weak, dipole allowed 2p exciton
feature starts at 3.034 eV. Finally, phonon-assisted features are also
identified that correspond to an indirect gap of 3.049 eV.

There are several important points of comparison.
First, the calculated lowest exciton energy at Γ
is about 0.22 eV higher than measured.
Broadly, this error is comparable to those encountered when
using the GW approximation for other semiconductors Aryasetiawan and Gunnarsson (1998); Aulbur et al. (1999); Onida et al. (2002).
However, it is important to be clear that the calculation
is performed for a frozen lattice with no account given
for electron-phonon interactions. In general, the electron-phonon
interaction will reduce the zero-temperature quasiparticle gap Cardona and Thewalt (2005).
Second, the present GW calculation gives the conduction band minimum
at the R point instead of being at Γ, as suggested by the optical measurements.
The energy differences are small; the calculated conduction band at R is 0.04 eV lower than at Γ.
In the analysis of the absorption data, the indirect gap is found to be 0.01 eV higher
than the direct gap Pascual et al. (1977), albeit including what ever role excitonic effects may have.
The difference between theory and experiment is too subtle to be resolved in these calculations, particularly
without the influence of electron-phonon interactions.
Third, the symmetry of the excitonic states from the calculation agrees
with the interpretation of the experiments.
However, the scale of the excitonic effects that we calculate using the static
dielectric matrix, and only including the electronic polarizability,
is substantially larger than suggested by the experiment.
Our calculated long-wavelength dielectric constant (∼ 8) is
slightly higher than the measured ϵ∞ (∼ 7) Cardona and Harbeke (1965), but similar to previous calculations Vast et al. (2002).
The lattice contribution is quite large, with ϵ0 = 111 Pascual et al. (1977).
This again points to the importance of the electron-phonon interaction.

The measured optical absorption in single crystal anatase at low temperature
does not resolve any significant structure Tang et al. (1995).
The energy dependence near the onset of absorption is consistent
with an Urbach tail. Analysis of the temperature dependence
leads to an estimate for the band gap for extended states of 3.42 eV Tang et al. (1995).
This exceeds the measured optical threshold in rutile by about 0.4 eV.
Since the measured absorption edge in anatase is featureless,
another way to characterize the absorption edge and
make comparison to rutile is to consider the energy
at which the low temperature absorption is 50 cm−1
for electric vector perpendicular to the c-axis.
In rutile, this occurs at 3.04 eV while it occurs at 3.30 eV in anatase Tang et al. (1995).
This suggests a more modest 0.26 eV difference between
the minimum energy gap of anatase and rutile.
Our calculated quasiparticle minimum energy gap in anatase is 3.56 eV, modestly
larger than the value deduced from the absorption measurements.
The calculated difference in gaps between anatase and rutile is 0.22 eV,
similar to the measured difference.

Figure 2: Density of states (DOS) of rutile derived from FF GW calculations
plotted with photoemission and inverse photoemission spectraTezuka et al. (1994).
The solid curve is the calculated DOS which is convoluted with a Gaussian
function of σ= 0.5 eV, while circles are photo- and inverse
photo-emission spectra.

So far, the calculated quasiparticle energies for TiO2 have been found
to agree well with electron spectroscopy,
but the minimum energy gaps, including excitonic effects, are larger compared to
the measured absorption threshold.
Furthermore, the strength of the calculated exciton binding energy
is larger than implied by the interpretation of the optical spectrum near threshold.
To get additional perspective, we calculate the
macroscopic dielectric function over a broad energy range,
including the correlations induced by electron-hole interactions.

In Fig. 3, we show the imaginary part of the dielectric
function of rutile for polarizations both perpendicular and parallel
to the tetragonal axis c. The solid curves are calculated from the BSE,
while the dashed curves are derived from optical reflectivity measurements
at room temperature Cardona and Harbeke (1965). For both polarizations, the
theoretical spectra are close to experiment up to about 6 eV.
In particular, the first strong peak at ∼ 4 eV for both
polarizations is reproduced very well by the BSE results.
Above 6 eV, the overall magnitude and prominence of the peaks
in theoretical spectra are distinct from experiment.
The ϵ2(ω) of anatase is displayed in Fig. 4,
where the experimental data Hosaka et al. (1997) were measured at 100
K. Similar to the rutile case, the theoretical calculations capture
the features around the onset of major absorption at 4 eV.
For perpendicular polarization, the calculated oscillator strength is systematically
too large starting at about 5 eV.
The calculated results are very similar
to the previous calculation for rutile TiO2Lawler et al. (2008).
For anatase, the calculated peak heights near 4 eV appear less intense in their spectra,
but this is largely due to their choice of a larger damping parameter, as is evident from
the broadening on the low energy side of their spectra.
In particular, we have analyzed the integrated oscillator strength (i.e. the contribution
to the f-sum rule) from the first peak in the anatase spectrum for parallel polarization.
We find that our oscillator
strength is essentially the same as theirs, but that both calculations show more
oscillator strength than is found in the experimental spectra by about 30%.

The systematic overestimation of oscillator strength at higher energy appears
to be a more general issue. For example, a recent BSE study for several alkaline earth metal monoxides shows some
similar excess oscillator strength at higher photon energies Schleife et al. (2009).
This may well trace to more fundamental assumptions in the methodology.
Two key issues are the use of the Tamm-Dancoff approximation, which has been
identified recently as the main source of errors in the calculations
of a confined system Grüning et al. (2009), and the assumption
of a statically screened Coulomb interaction Rohlfing and Louie (2000); Marini and Del Sole (2003).
Also, as noted by Lawler and coworkers, the f-sum rule for the oscillator strength
converges very slowly in the titanates and the experimental analysis that relies
on Kramers-Kronig transformations may have systematic errors as well Lawler et al. (2008).

Figure 3: ϵ2(ω) of rutile from 0 to 12 eV. Solid curves are
theoretical calculations with BSE, and dashed curves are experimental
results Cardona and Harbeke (1965) obtained at room temperature. In (a) the
direction of polarization is perpendicular to the tetragonal axis
c, and in (b) the direction of polarization is parallel to axis c.Figure 4: ϵ2(ω) of anatase from 0 to 12 eV. Solid curves are
theoretical calculations from BSE, and dashed curves are experimental
results obtained at 100K Hosaka et al. (1997). In (a) the direction
of polarization is perpendicular to the tetragonal axis c, and in
(b) the direction of polarization is parallel to axis c.

iii.3 External Energy Level Alignment: TiO2(110) Surface

The energetic position of the quasiparticle VBM with respect to the vacuum level
in the rutile phase for the (110) surface was calculated in several steps.
First at the LDA level, analogous to the determination of
work function for metals,
the electrostatic potential change between the bulk like region and the vacuum region
was determined Zhang et al. (1988); Baldereschi et al. (1988); Fall et al. (1999); Shaltaf et al. (2008).
Then in a second step, the bulk VBM position relative to the electrostatic potential
is determined.
Previous calculations of work functions based on metal slabs
based on the DFT Fermi energy
were in fairly good agreements with experimental measurements Fall et al. (1998, 2000).
In the third step, we will apply the perturbative correction to the LDA
exchange-correlation potential to determine the alignment
in the GW approximation Zhang et al. (1988); Shaltaf et al. (2008).

For semiconductors and insulators the Fermi level is separated from
the VBM, and the position of VBM varies as the surface structure changes.
To calculate the energetic alignment of rutile with respect to the
vacuum from the (110) surface, a stoichiometric slab of six layers
of Ti atoms is cut from the bulk structure with geometrical parameters
described in Sec. I. A vacuum space of the same thickness of the slab
is used to buffer the two surfaces of the slab. The surface is cleaved
as a (1×1), unreconstructed surface as observed in experiments Diebold (2003).
The k-point mesh is 2×4×1. After the slab is fully relaxed
with the atomic positions of the central two layers fixed, the electrostatic
potential inside the slab with respect to the vacuum is measured as
-12.65 eV.
The energy difference between the VBM (from the LDA calculation) and
the electrostatic potential inside the material is 5.05 eV, which
is determined in a separate bulk calculation to prevent the quantum
size effects Fall et al. (1998). Accordingly, the LDA value for the VBM
relative to vacuum is 7.60 eV. Previous LDA and gradient corrected
calculations for give 7.16 eV (LDA, three Ti layer slab) Vogtenhuber et al. (2002),
7.2 eV (PBE, 11 Ti layer slab) Kiejna et al. (2006) and 7.6 eV (PW91, 11 Ti layer slab) Kiejna et al. (2006).
In light of variations at the 0.2 eV level with number of layers in the slab model Kiejna et al. (2006),
the overall agreement is satisfactory.

The GW correction for the VBM calculated with FF model and 160 bands
is 0.07 eV. However, as noted in Sect. IIC, extrapolation to full
convergence with respect to the total number of bands will drive this
0.2 to 0.4 eV lower. We therefore
suggest a GW correction of -0.2 eV,
with about 0.1 eV uncertainty. The final prediction from GW for the
VBM alignment to vacuum at the clean rutile TiO2 (110) surface
is 7.8 eV.

In order to deduce the VBM alignment form experiment, two
results must be combined: (1) the work function which fixes the Fermi energy
relative to the vacuum; (2) the position of the VBM relative to the Fermi energy.
Experimental measurements of the work function of rutile from the
(110) surface Onda et al. (2004); Schierbaum et al. (1996) vary from 4.7 eV to
5.8 eV depending on the structure of the surface which is strongly influenced
by treatment (annealing, exposure to oxygen, etc.).
In addition, the position of the VBM relative to the Fermi energy also
depends on surface treatment Wendt et al. (2008).
Therefore, it is crucial to compare with data in which both values are
measured on the same sample.
To our knowledge, this is relatively rare.
Based on UPS measurements with (hν=21.2 eV), a work function of 5.2 eV
and a relative VBM position of 3 eV were measured for TiO2(110) Schierbaum et al. (1996),
which implies the VBM position relative to vacuum falls at 8.2 eV.
Similar measurements for TiO2(100) yield a work function of 4.9 eV and a relative VBM position of 3.1 eV respectively Xiong et al. (2007),
yielding the VBM position at 8.0 eV.
The difference in workfunction between these measurements is consistent with
a separate Auger Microprobe study of facet dependence Imanishi et al. (2007).

Based on this limited data set, the GW based prediction for the VBM
alignment is off by about 0.4 eV.
There are at least 0.1 eV uncertainties in both the theory and the
experiment. Since both the work function and position of Fermi level
are sensitive to the surface properties, the deviation between the
theoretical calculation and experimental measurements may well reflect
the complexity of the TiO2 surfaces. For example, recent studies
suggest that the commonly employed strategy of cleaning followed by
annealing in oxygen may not result in the ideal surface envisioned
(e.g. with no oxygen related defects) Wendt et al. (2008).
In particular, the physical origin of the widely observed defect states
around 1 eV below the CBM remains a point of vigorous discussionDi Valentin et al. (2006); Wendt et al. (2008).

We have presented a numerically well converged MBPT study of
the electronic and optical excitation energies in rutile and anatase
crystals of TiO2. The calculations are carried out in
the approximation of a frozen lattice, without consideration of electron-phonon coupling.
In most respects, the agreement with
experiment is well within the expected accuracy of this approach.
In particular, the calculated quasiparticle gap agrees with electron
spectroscopy measurements (photoemission and inverse photoemission),
the change in the minimum gap between rutile and anatase crystal
structures is reproduced,
and the main features of the optical spectrum agree with ellipsometry measurements.
The qualitative features of the zone center bound excitons calculated
for rutile agree well with low temperature absorption measurements.
However, the scale of the exciton binding energy is larger than that estimated from
experiment by about a factor of 10 and the calculated exciton energy is about
0.2 eV larger than measured.

The key theoretical assumptions in our application of MBPT include
evaluation of the electron self energy in the GW approximation without
iterating to self consistency or considering vertex corrections.
Self consistency would certainly increase the calculated energy gap
through the reduction in the screening
van Schilfgaarde et al. (2006); Shishkin and Kresse (2007).
Recent results for a set of other semiconductors and insulators shows that
approximate inclusion of vertex corrections in screening leads to
a partially compensating reduction of the calculated energy gapShishkin et al. (2007).
However, a fully
consistent approach
to include vertex corrections remains subject of current research in the field.
Because the Ti 3d electrons are almost completely ionized,
TiO2 should not be subject to the sorts of systematic errors
found in non-selfconsistent calculations for some late transition 3d metal
compoundsAryasetiawan and Gunnarsson (1995); Faleev et al. (2004); Bruneval
et al. (2006b)
The accurate results found for the key optical transition energies
contributing to the absorption (Figs. 3 and 4) support this,
and contrast to the case of Cu2O where non-selfconsistent
calculations showed substantial discrepenciesBruneval
et al. (2006b).
In the solution of the BSE, the calculated static (electronic) dielectric
matrix has been used and the equations were simplified through the Tamm-Dancoff approximation.
These approximations are part of the standard
MBPT treatment of optical spectra and
the low energy excitons are expected to be treated wellRohlfing and Louie (2000); Onida et al. (2002).
However, they may affect higher energy features in the spectraMarini and Del Sole (2003).

On physical grounds, we suggest that the most significant open issue
concerns the role of electron-phonon coupling.
In general, the electron-phonon self energy
will both lead to a smaller zero-temperature quasiparticle gap
and make a significant contribution to the temperature dependence of the energy gap Cardona and Thewalt (2005).
As noted in the text, there is a very large difference between
the electronic dielectric constant ϵ∞ and
the low frequency dielectric constant including lattice polarization ϵ0
for TiO2.
This suggests a relatively strong electron-phonon interaction
and there is a long standing debate over the polaronic character
of charge excitations in TiO2PNo ().

Based on the usual form of the Frohlich interaction, the dimensionless
coupling constant characterizing electronic coupling to the most important
polar optic mode for rutile (with mode energy about 0.1 eV)
is given by α=1.6√mb/me where mb is the
bare band mass and me is the free electron mass.
Using our DFT band dispersions to have estimates, the electron (hole)
band mass is about 0.6me (1.8me) along x or y and
1.6me (3.1me) along z. This suggests coupling constants
of α∼1−2 for electrons and α∼2−3 for holes which
would fall in the weak to intermediate coupling regime. Using the
usual weak coupling expression, the electron and hole renormalization
would be 0.1−0.2 eV and 0.2−0.3 eV respectively, both of which
act to reduce the quasiparticle energy gap. For anatase α∼1.6√mb/me,
essentially the same as rutile, based on the mode energy and dielectric constants
Gonzalez et al. (1997). The electron
(hole) band mass is about 0.4me (1.8me) along x or
y and 3.9me(1.0me) along z, suggesting slightly different coupling
constants of α∼1−3 for electrons and α∼2 for
holes with corresponding (weak coupling) electron and hole renormalization
of 0.1-0.3 eV and 0.2 eV respectively.
In weak coupling, the electron-phonon self energy is added to the results obtained
here based on the GW approximation.
For the analysis of the optical absorption edge, a more detailed calculation
is required because in the exciton-phonon coupling, the exciton is neutral and
the electron and hole distortions of the lattice
will partially cancel Schmitt-Rink et al. (1987); Rudin et al. (1990).
Taken together, if the large polaron regime is physically correct,
these rough estimates suggest that the effect of the electron-phonon coupling could
account for some of the differences between the present GW/BSE
results for the frozen lattice and experiment.
Firmer conclusions require a more extensive set of calculations,
beyond the scope of this article
Cardona and Thewalt (2005); Marini (2008); Bechstedt et al. (2005); Vidal et al. (2010).
In particular, it may be that a more complete consideration of self-consistency
and vertex corrections in the electron-electron contribution to the electron self energy will need
to be combined with an analysis of the electron-phonon contribution.
The two contributions should be treated
in a fully consistent theory.
In more empirical terms, an overstimate of the band gap based on
selfconsistent treatment of the electron-electron
interactions alone may be compensated by the electron-phonon contributions.

It may well be the case that the large polaron regime is not applicable for TiO2.
A recent THz spectroscopy study
of rutile gave a direct measurement of the electron scattering rate Hendry et al. (2004).
This data was analyzed with a Frohlich form for the electron-phonon interaction,
but regarding the coupling constant as a free parameter.
Using the Feynman approach Feynman (1955); Schultz (1959) to handle intermediate
to strong coupling, the analysis showed coupling constants for electrons
α∼4−6 depending on field orientation Hendry et al. (2004).
The inferred electron mobilities were consistent with earlier electron transport measurements Yagi et al. (1996).
These values suggest a substantially larger value for the
electron self energy of 0.4−0.7 eV Schultz (1959). An older estimate based on
small polaron theory also suggested 0.7 eV Eagles (1964). A recent DFT+U based study
suggested that an excess electron in rutile is indeed self trapped Deskins and Dupuis (2007). Although
the binding energy was not presented, the barrier for polaron hopping
was estimated to be 0.3 eV.
A similar study for an excess hole in rutile
suggested barriers of 0.5 - 0.6 eV Deskins and Dupuis (2009).
Taken together, if the small polaron regime is found to be physically relevant,
then the quasiparticle and excitonic energies will need to be fully reanalyzed.
For strong electron-phonon coupling a perturbative
approach to combine the electron-electron and electron-phonon self energies is no longer justified.
Furthermore, the electron-phonon coupling enters into the spectroscopic measurements
in distinct ways. The (inverse) photoemission and optical absorption would each need to be
properly analyzed.

Acknowledgements.

We thank Dr. A. Marini for access to the private
branch of the Yambo code.
W. K. thanks Dr. D. Prezzi for insightfull discussions on BSE.
Work performed under the auspices
of the U.S. Department of Energy under Contract No. DEAC02-
98CH1-886. This research utilized resources
at the New York Center for Computational Sciences
at Stony Brook University/Brookhaven National Laboratory
which is supported by the U.S. Department of Energy under
Contract No. DE-AC02-98CH10886 and by the State of New York.