Charge-transfer states and optical transitions at the pentacene-TiO2 interface

Pentacene molecules have recently been observed to form a well-ordered monolayer on the (110) surface of rutile TiO2, with the molecules adsorbed lying flat, head to tail. With the geometry favorable for direct optical excitation and given its ordered character, this interface seems to provide an intriguing model to study charge-transfer excitations where the optically excited electrons and holes reside on different sides of the organic–inorganic interface. In this work, we theoretically investigate the structural and electronic properties of this system by means of ab initio calculations and compute its excitonic absorption spectrum. Molecular states appear in the band gap of the clean TiO2 surface, which enables charge-transfer excitations directly from the molecular HOMO to the TiO2 conduction band. The calculated optical spectrum shows a strong polarization dependence and displays excitonic resonances corresponding to the charge-transfer states, which could stimulate new experimental work on the optical response of this interface.


Introduction
Investigations on hybrid organic-inorganic systems [1][2][3][4][5][6][7][8] have opened up the possibility to combine the benefits of both constituents and devise new applications that outperform their either purely organic or inorganic counterparts. For example, the possibilities of tailoring the electronic band structure can be considerably extended compared to inorganic heterostructures, because the organic molecules can more flexibly adjust to the underlying lattice of the inorganic substrate [9,10]. The band structure of a hybrid system can be further tuned by modifying the properties of the organic material [1,10]. Such hybrid systems can also provide new types of semiconductor excitations, like the hybrid Frenkel-Wannier exciton [11] with a high oscillator strength combined with an enhancement of nonlinear optical response [2,12], or the hybrid charge-transfer exciton [5][6][7][8] where a bound electron-hole pair is formed with the electron and hole residing at different sides of the organic-inorganic interface of the system. These ideas have already been applied in photovoltaics [10], where dye-sensitized solar cells [3,[13][14][15][16][17] have been found to be promising in the search for efficient, low-cost, and environmentally friendly devices. The operational principles of solar cells are fundamentally dependent on charge-transfer excitations at the interfaces between the two constituents [5][6][7][8]. Even though charge-transfer excitations have been intensively studied during the last decades [5][6][7][8][18][19][20][21][22], the charge-transfer nature of the excited states [5,8,23] is not well understood. In order to enhance the performance of photovoltaic devices, it is especially important to characterize the mechanisms involved in the formation [5,23] and dissociation [5,6] of such excitons [8]. Therefore, it is desirable to study systems where the charge-transfer interface-excitation states are strongly identifiable, e.g., in the optical spectrum. Such model systems could reveal the fundamental character of the charge-transfer states through clearly visible features in the optical spectra.

Methods
To determine the geometry of the system, we optimize it using the SIESTA DFT code [37]. We use a slab geometry with a 1×6 supercell (with respect to the primitive surface unit cell) in the surface plane, and five layers of TiO 2 in the direction perpendicular to the surface. In total, we have 216 atoms in the unit cell and use 20Å of vacuum in order to minimize the interactions between periodically repeated slabs. To include van der Waals contributions to the energy and forces, we employ the optB88 functional [38]. Our double-ζ polarized (DZP) basis set of numerical atomic orbitals is generated using an energy-shift parameter of 100meV (see [37] for a detailed explanation of these parameters), and the core electrons are replaced by Troullier-Martins pseudopotentials [39]. Starting from the DFT results, we obtain the relevant band structure and wave functions to compute all the quantities that we need to set up and solve the BSE, as described below.

Excitonic resonances in optical spectra
An exciton consists of a Coulomb-bound electron-hole pair [29] that can be described, e.g., by the two-particle Green's function. Here, we analyze excitonic resonances in linear absorption by applying the standard BSE approach [27,40]. Technically, we solve an eigenvalue problem where nk  is the quasiparticle energy of an electron in band n with crystal momentum k  and f nk is the electronic occupation. A scissor operator can additionally be introduced for nk  without any loss of generality.
contains quasiparticle wave functions nk y as well as the bare (v) and screened (W) Coulomb interactions. The first term in(3) is the repulsive exchange contribution (with the factor of two appearing for the optically relevant singlet excitation) and the second is the attractive direct term; the one responsible for excitonic binding.
We construct the linear optical absorption from the imaginary part of the macroscopic dielectric tensor M ii e w ( ) in different Cartesian directions i, following the procedure in [27,28]. In practice, we make use of the representation Here V is the volume of the supercell (or quantization volume) related to the sampling of the Brillouin zone and S A A nm nm nm k k k , * º å l l l l ¢ ¢ defines the overlap between the (right) eigenvectors. In(4) and (5), are the dipole matrix elements with p î the momentum operator in Cartesian direction i and the 2 factor appears for a singlet exciton [40]. When non-local pseudopotentials are used to compute the wave functions a correction to the matrix elements of p î should also be included [41,42]. Note that the energies in the denominator in (6) should not be modified when a scissor operator is used.
In the Tamm-Dancoff approximation [27] (TDA) where particle-hole pairs are assumed to be uncoupled from the hole-particle pairs the solutions of(1) have the meaning of expansion coefficients of an exciton wave function where indices v and c denote occupied and unoccupied states, respectively.

Computational aspects of the spectrum calculations
We solve the equations above numerically using a set of newly developed codes that act as post-processing tools for the SIESTA code. In our numerical solutions we approximate the quasiparticle wave functions r nk y ( ) as the DFT Kohn-Sham solutions; an approximation commonly used in GW/BSE calculations. We take the band energies nk  to be the Kohn-Sham eigenvalues modified by a scissor operator to open up the band gap, as explained in detail later. In practice, the matrix elements of the DFT Hamiltonian are imported from SIESTA in a basis of local atomic orbitals and are rediagonalized to obtain wave functions and band energies for all relevant k-values. We modified the SIESTA code to export momentum matrix elements in the local basis, including the corrections due to the use of non-local pseudopotentials. These are subsequently imported and converted to the eigenstate representation to be used in (6).
The Coulomb integrals in(3) are computed by making use of a real space grid inside the unit cell and its corresponding plane wave expansion. For a given plane wave cutoff energy the number of grid points in real space as well as in reciprocal space is determined automatically and one can go between the two representations via Fourier transforms. First, the periodic part of the wave functions, which at this point are expanded in localized orbitals, are converted to the real space grid representation, then the electrostatic potential originating from each product state, e.g., v rv r rr r r d , is solved by going to reciprocal space where the Coulomb interaction is diagonal. After a back transform to real space, the potential is multiplied with the other product state on the grid and integrated, e.g., We use a 2Dtruncated Coulomb interaction [43,44], with the spatial cutoff set to half the cell in the normal direction. The Coulomb singularities are treated by replacing each diverging contribution with its numerical average over a small volume in k-space (defined by the original k-mesh) around the divergent point [44]. The averaging is performed using an auxiliary 10×10×10 Monkhorst-Pack k-mesh centered at, but excluding, the divergent point. A plane wave cutoff of 300 eV is used (differences to computed spectra with 600 eV are verified to be negligible).
The electronic part of the dielectric screening can, in principle, be approximately modeled ab initio using the random phase approximation. However, such a procedure would be numerically extremely costly because our unit cell is so large. Furthermore, it has been suggested [28,36] that contributions from phonon modes in bulk TiO 2 can be important for the effective dielectric function necessary to properly describe the lowest lying excitons. Rutile TiO 2 has a large difference between the purely electronic screening . A good agreement with TiO 2 experiments in bulk was found [36] by using ò=47.8 (now geometrically averaged over the directions) to compute excitonic absorption.
In our interface-dominated system the dielectric behavior will be strongly modified from that of the bulk. A simple estimate produces a dielectric constant (for q approaching zero) that is an average of bulk TiO 2 and vacuum [46], yielding when taking the best fit value above [36]. In [47] an image charge model was used to fit the dielectric constant for a molecule adsorbed above a TiO 2 (001) surface, with the best fit given by 2.76  = , which is lower than the estimate above. These limiting cases indicate the range of the expected ò as the screening changes locally across the TiO 2 -pentacene interface.
Since computing the screening microscopically is presently extremely challenging for such a large system, we model it by an effective dielectric constant over the reasonable range. Specifically, we explore three models where the dielectric constant is set to a low ( 1  ), middle ( 2  ) and high ( 3  ) value. For the low and high values, we choose the two averaged values of the TiO 2 bulk and vacuum discussed above, that is is chosen to give the average binding energy of 1  and 3  , assuming the binding energy to be 1  µ as is the case for a two-state system using(2).

Results and discussion
The geometry is optimized by keeping the middle TiO 2 trilayer frozen while letting all other atoms, including the molecule, relax. The lattice constant in the in-plane direction is fixed to the optimized bulk value obtained with the same functional and basis set (a=4.60 Å, b=2.98 Å), giving a surface 1×6 unit cell of a=6.51 Å and b=17.87 Å. The geometry is schematically shown in figure 1 where gray denotes Ti atoms, red O, blue C, and white H. The pentacene molecule is adsorbed over the surface, tilted inwards to the troughs formed by the Ti-O rows on the TiO 2 (110) surface with an angle of 24°, in agreement with the experimentally deduced 25° [26]. The closest C-Ti distance is 2.76Å and the closest C-O distance (coming from an oxygen row) is 2.86Å. Comparing the surface with and without the molecule shows that the surface bands are almost unaffected by the adsorption, except for the appearance of two flat bands in the band gap that are both situated below the Fermi level. Also the lowest conduction band comes down slightly, while still being almost degenerate with another band between the high-symmetry points X and M.
To assign the bands we take a look at the PDOS. The two gap states below the Fermi level mostly have contributions from the carbon atoms, that is, they belong to the pentacene molecule. The first unoccupied bands, however, belong mostly to the slab and have predominant contribution from Ti 3d orbitals. Indeed, this part of the band structure is very similar to the one coming from the bare surface, without the molecule, so the effect of the molecule on these bands seems to be small. The lowest unoccupied band is seen to be delocalized through the whole slab, so even though it has contributions from the surface atoms it is not a pure surface state. Higher up in energy, we see the contributions from the unoccupied pentacene states, centered at around 0.5 eV above the lowest surface conduction band, with the lowest weak feature at 0.25 eV above the same. The combination of well-separated molecular occupied states and and low-lying unoccupied states belonging to the substrate suggests the possibility of optical charge-transfer excitations from the molecule to the surface where the unoccupied molecular states are not expected to contribute much.

Band alignment
It is well-known that the band alignment can sometimes be a problem in DFT calculations for molecules on surfaces due to an unbalanced description of the different constituents. To investigate whether or not DFT qualitatively gives the correct band alignment for our system we look at a simple two-step procedure that has been successfully used to align molecular levels for surface adsorbates on metal and semiconductor surfaces [47][48][49]. In the first step individual shifts for the levels are computed in separate calculations for the bare substrate and the gas phase molecule. In the second step an image charge model is used to include the additional screening provided by the substrate that modifies the molecular electron removal and addition energies. In this way individual shifts for the substrate and molecular level are obtained, which can later be applied to the combined molecule-substrate DFT calculation to align the bands.
In order to obtain an unbiased description of the levels of the molecule and the surface we use (the negative of) the experimental ionization potentials and electron affinities as reference values. For the bare substrate, experimentally deduced values are in the range −8.5 to −8.7 eV for the highest valence band and −5.2 to −5.4 eV for the lowest conduction band [50].
)is the screened image charge, z is the height of the molecule above the surface and z 0 is the position of the image plane. To determine z 0 we fit E unocc D to match the exchange-correlation potential for the TiO 2 surface [47][48][49], giving us a plane 0.25 Å above the oxygen rows. We use the directionally averaged experimental  ¥ for bulk TiO 2 , that is 7.4  = , and consider the pentacene molecule to be located 2.8 Å above the surface, resulting in a shift of±1.08 eV. The image charge correction reduces the shift for the molecular levels to −0.92 eV and 0.87 eV for HOMO and LUMO, respectively.
Applying these shifts to the DFT band structure, as seen in figure 2, we see that the molecular HOMO is still in the band gap of the surface and the molecular LUMO is shifted up in energy in relation to the surface conduction band edge. In other words, our two-step procedure predicts that that the lowest excitations will be from the molecule to the surface, as the DFT band structure suggests, with little influence of the molecular LUMO that has a much higher energy. We also note that the band gap between the molecular HOMO and the conduction band edge has opened up from 0.55 to 2.02 eV(1.82 eV) where the value in parenthesis refers to the lower value used for the experimental reference. Furthermore, the corrected band gap shows a relative insensitivity to the dielectric constant that enters in the image charge formula; if the bulk 0  would instead be used (close to full screening) the value of the gap would be only 0.32 eV lower. In light of this, we feel justified in using the DFT band structure in the optical calculations, only opening up the band gap with a scissor operator.

Optical absorption and polarization dependence
All optical spectra are computed using one fully occupied band v 1 and five unoccupied bands c 1 -c 5 with a 55×19×1 Monkhorst-Pack k-mesh spanning the whole Brillouin zone (BZ). We use a scissor operator of 1.47 eV that gives a fundamental gap in our system (from the molecular state in the gap to the lowest conduction band) of 2.02 eV, in accordance with the higher value from our model for the band alignment. This estimate cannot be fully checked due to absence of experimental data. However, the exact value of the gap will not influence the relative positions of optical resonances within the TDA. Our results verify that this approximation gives visually identical spectra compared with the full solution of(1). For an easier interpretation, we discuss the TDA results in the remainder of this paper.
The computed linear absorption spectrum ( Im ii M e w µ [ ( )]) is shown in figure 3 for the three model dielectric functions, together with the non-interacting spectrum. The purple, blue, and orange areas correspond to response in x-, y-, and z-polarization direction, respectively. We observe that the Coulomb interaction changes all the spectra substantially, giving also a high sensitivity to the screening model used. In other words, the actual  should be directly deducible from absorption experiments.
While the quantitative spectra vary drastically as a function of  , some qualitative features remain  -independent. For example, the spectra exhibit a strong dependence on polarization direction, to a degree that different directions contain different resonances. This is seen most clearly for 1  and 2  , producing three separate excitonic resonances labeled A, B, and C. The lowest energy resonance A is clearly visible only for the z-polarized light while the B resonance is connected with the y polarization. The resonance C dominates the spectrum for the x and z polarizations, but practically vanishes for the y polarization. All resonances A, B, and C indicated in figure 3 originate from individual excitonic states. In the noninteracting system, by definition, we do not have excitonic resonances, but nevertheless even in this case the similar spectral regions corresponding to A, B, and C resonances are detectable, as shown in the lowest frame of figure 3.

Average hole and electron positions
For a given exciton, the averaged density of the electron and the hole follow from  almost identical distributions because they essentially originate from the same unoccupied band c 1 , although for different k-values and symmetries in the Brillouin zone, as will be discussed later in more detail. However, the C excitons mainly stem from the c 2 band, and although higher in energy, it has part of the density closer to the molecule and consequently has more overlap with the hole state. This partly explains the higher absorption peak in figure 3.

Band mixing in excitonic resonances
To characterize the spectral features seen in figure 3, we first separate the different contributions to the spectra by using only one conduction band at a time in the sum (4). Figure 5 shows this separation for the c 1 -c 5 band contributions in the polarization-averaged absorption for 1 3   -, as well as for the noninteracting case. In the noninteracting case, each eigenvector A l is connected to one conduction band only, but when the Coulomb interaction is included the transitions become mixed, i.e., A l will contain contributions from several c n . In this case, even when only one c n is considered in the sum of(4), the related density-change-matrix element v c i k n 1 r w ( ) will include contributions from other conduction bands than c n . Thus, the contributions of different bands to the spectra can be exactly disentangled by the used approach only in the noninteracting case, while the separation is only approximative when the Coulomb interaction is fully included. The Coulomb-induced mixing of transitions also results in some negative spectral features in the decomposition as seen in the interacting cases of figure 5. Nevertheless, the shown results give a good general picture of which c n bands the A, B and C resonances are most closely connected.
For the noninteracting case, we see that the A and B regions in the absorption originate from c 1 contributions while the C region comes predominantly from c 2 . The same holds for 3  which has a large enough screening for the bands not to mix in A l appreciably due to the Coulomb interaction. We observe that 1  and 2  yield similar trends, while making  smaller enhances mixing of transitions such that eventually all bands contribute to all excitonic resonances.

Excitonic states at the interface
To identify the k dependence of excitons one would optimally plot A nmk l and D nm i k as functions of k. However, when taken directly from our DFT-based results, both of these contributions may contain arbitrary phase jumps in k, due to the 'gauge-freedom' [52] problem, i.e., the freedom in choosing the phases of nk y . Additionally, as k is generally a state index rather than an actual wave-vector coordinate , there is no need of A nmk l and D nm i k being continuous functions of k; even without the phase jumps. Indeed, u r nk ( ) can become discontinuous in k at band intersections [53,54]. This complicates visualization, especially, when 'band-mixing' [55] is significant, i.e., when u u n m k k 2 á ñ ¢ | | | differs considerably from unity between k and k¢ in immediate vicinity but on opposite sides of the intersection, for all choices of n and m. Such band-mixings are expected to be important for large unit cells due to effects of band folding with respect to the bands for the bulk and/or primitive surface unit cells.
Since v 1 and c 1 are energetically isolated, they are not affected by the above band-mixing, yielding a continuous D c v i k 1 1 after matching the phases between k states. In fact, our diagonalization scheme produces ) is fulfilled. In this system, c 2 -c 5 bands exhibit so strong band mixing that following exciton through only one of these bands is not useful. Therefore, we investigate an effective exciton distribution f l for the C-resonance excitons of 1  and 3  , respectively. Each frame is centered at X k = ¢ that is denoted by a filled circle in the center of each frame. The A excitons are centered at X' with a predominantly even (real) s-like part in the Brillouin zone while having also an odd (imaginary) p-like part. The odd part contains from 40% ( 1  ) to 20% ( 3  ) of the weight in the total A vck l . The 1  B exciton (frames (g)-(i)) shows even symmetry (60% of the weight in A vck l ) and resembles a s-like state at M-point while having a sizable odd (40%) contribution located at X' (or at M). The B state changes drastically when ò is increased to 3  (figure 6(j)), because the even fraction is only 20% and A Re v c k ] has a d-like symmetry located at X'. At the same time, its p-like features centered at X' are strongly increased ( figure 6(k)). The C states experience strong band mixing, and we analyze only k f l in figures 6(m)-(n). These C excitons have even symmetry (more than 80% of the weight in A vck l ) and are predominantly s-like states centered at X'.
The major contribution to these states at different k comes from the c n -band u r  This detailed analysis of states yield the following conclusions: a resonances stem predominantly from c 1 -band s-like excitons located at the X' point and the C resonance from predominantly s-like c 2 -band states at X', independently of  . However, the B-resonance changes its nature from a predominantly c 1 -band s-like state at M-point to a c 1 p-like state at X' between 1  and 3  , respectively.

Combined symmetries of excitons and dipoles
Since we have a highly asymmetric unit cell and large portions of the BZ need to be considered, schemes based on symmetries of the carrier states to evaluate optical strengths of excitons are not convenient. Instead, we directly study the combined properties of A v c k [ ] in figure 6(f), we notice the similarity (especially for 1  ) between the symmetries and locations of the dipole and the exciton contributions. The found match makes the A resonances relatively well visible in the spectra for z polarization.
The same comparison as done above for the A resonance can be made for the B resonance in figures 6(g), (h), (j), and (k) with respect the dipoles of figure 7. This time, we find the match for y-polarization, explaining why the B resonance is optically active only for this particular polarization direction.  resonance are shifted in the BZ so that the position of X' and Γ points are exchanged. In addition, the x-and z-directional dipoles for C resonances are roughly a factor of seven and a factor of two stronger than the x-and z-directional dipoles related to c 1 band. These remarks explain why the C resonance is visible for x and z polarizations, but is optically inactive for the y polarization, and why the C resonances have a higher intensity compared to A resonances in figure 3.

Conclusions
We have modeled the pentacene overlayer of the rutile TiO 2 (110) surface by DFT and computed excitons and optical spectra by solving the Bethe-Salpeter equation (BSE) for a set of physically motivated system Hamiltonians. Flat bands coming from the pentacene molecule are found in the band gap of the pristine surface, and the BSE solutions reveal excitations with charge-transfer character from the molecule to the surface. The optical spectrum shows a large sensitivity to the polarization direction which we explain by investigating the exciton solutions and transition dipole matrix elements in the Brillouin Zone (BZ) assigning the three dominant transitions by their approximate symmetries and locations in the BZ. Many qualitative features of the excitons are weakly dependent on the dielectric screening chosen in the system Hamiltonian, although the energetic positions of the peaks are more sensitive to the screening. An experimental characterization of the optical spectrum of such a well-defined organic/TiO 2 system will certainly help to get a deeper understanding on the dynamics of formation and decay of charge-transfer excitations at this technologically relevant interface.
dielectric function M e that gives the optical spectrum. In the SBE, usually both the exchange and direct terms are screened, using the rationale of a background dielectric constant coming from inactive electrons that do not contribute directly in the optical process [29,36]. Since the screening in the SBE can be taken to be a parameter it can however also be chosen to coincide with the one of the BSE. In this case, the absorption coefficient [29,34] from SBE is proportional to the imaginary part of(4), i.e, to the results of figure 3. In summary, the BSE and the SBEs are two approaches to the same problem that give very similar (or identical, when using the same approximations) results in the linear regime, while the latter has applicability also for nonlinear phenomena, and the choice to use one or the other depends on the problem at hand.