{'fraction_non_alphanumeric': 0.03306211832570727, 'fraction_numerical': 0.01643313689533, 'mean_word_length': 4.561188118811881, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quasiparticle spectra and excitons of organic molecules deposited on substrates: G0W0-BSE approach applied to benzene on graphene and metallic substrates

We present an alternative methodology for calculating the quasi-particle energy, energy loss, and optical spectra of a molecule deposited on graphene or a metallic substrate. To test the accuracy of the method it is first applied to the isolated benzene (C6H6) molecule. The quasiparticle energy levels and especially the energies of the benzene excitons (triplet, singlet, optically active and inactive) are in very good agreement with available experimental results. It is shown that the vicinity of the various substrates (pristine/doped graphene or (jellium) metal surface) reduces the quasiparticle HOMO-LUMO gap by an amount that slightly depends on the substrate type. This is consistent with the simple image theory predictions. It is even shown that the substrate does not change the energy of the excitons in the isolated molecule. We prove (in terms of simple image theory) that energies of the excitons are indeed influenced by two mechanisms which cancel each other. We demonstrate that the benzene singlet optically active (E1u) exciton couples to real electronic excitations in the substrate. This causes it substantial decay, such as {\Gamma} = 174 meV for pristine graphene and {\Gamma} = 362 meV for metal surfaces as the substrate. However, we find that doping graphene does not influence the E1u exciton decay rate.


I. INTRODUCTION
Nowadays, π-conjugated organic molecules and their derivative films 1 are increasingly used in many applications. They are often used in organic electronic devices such as field effect transistors 2 or organic transistors. 3 Also, their good charge mobility and small optical band gap makes these materials suitable for photovoltaic applications, such as in solar cells. 4 Moreover, the spatial localization of organic molecules allows the light absorbed by the molecule to be converted into substrate excitations such as surface plasmon or electron hole excitations. The latter could be of interest in biosensing applications. 5 This has spurred recent studies characterizing the formation via cyclization cascade reactions 6 and the charge transfer within combined donor-acceptor layers 7 of π-conjugated organic molecules on metal substrates. However, the basic building block of the most utilized organic molecules, such as aromatic hydrogen carbonates, is the benzene ring. This work is motivated by all these potential applications, and is focused on exploring the quasiparticle and optical properties of benzene deposited on semimetallic (pristine graphene) and various metallic (doped graphene and Ag(jellium) surface) substrates. Special attention is paid to examine the influence of the substrate on the molecular HOMO-LUMO gap, exciton plasmon interaction, efficiency of the molecule mediated light substrate plasmon conversion, and the decay of the molecular excitons to electron-hole excitations in the substrate.
In the formulation of the problem we shall use previous theories which are well established and have been tested by various spectroscopic experiments. [8][9][10] The quasiparticle properties of the deposited molecule will be investigated in the framework of Hedins's GW theory 8,11,12 while the optical properties will be investigated by solving the Bethe-Salpeter equation (BSE) whose practical application was first developed by Strinati, 13 and more recently by Louie et al. 8,[14][15][16] These methodologies have been successfully used to calculate the electronic HOMO-LUMO gaps 17 and optical gaps 18 of benzene on various substrates. However, the computational complexity of such calculations has limited their range of applicability, resulting in a need for simple models and benchmarking of the various levels of approximation. For example, recent studies have shown that quasiparticle corrections to energy levels are linearly correlated with the fractions of the levels' densities within the substrate, molecule, and vacuum, 19 with the quasiparticle gap renormalization proportional to the molecule's height above the substrate. 20 However, to address these issues, the theoretical method presented here has been modified. Specifically, the optical spectra is obtained directly from the imaginary part of the dynamical 4-point polarizability matrix L kl i j (ω), which is the solution of the matrix BSE. In the standard resonant approximation [14][15][16] the BSE reduces to an eigenvalue problem and the optical spectra is obtained in terms of BSE eigenvalues and eigenvectors. Moreover, we develop the methodology to examine dark excitons, i.e. excitons that cannot be excited by an external electromagnetic field. They can be seen as the energy loss of an external dipole driving the molecule, and its intensity also can be expressed in term of the imaginary part of the polarizability matrix L kl i j (ω). We also show that the equilibrium molecule/substrate separation is large enough so that their electronic densities do not overlap. As a consequence of this, the only modification which has to be done in the formulation, after the substrate is introduced, is to extend the bare Coulomb interaction propagator V by the substrate induced Coulomb interaction V + ∆W. Because of all these different approaches we have rewritten all previous expressions to arXiv:1309.1359v1 [cond-mat.mes-hall] 5 Sep 2013 make it clear where our approaches separate.
The theory is first developed generally and then applied to an isolated molecule and to a molecule deposited on various different substrates. To check the accuracy of our method, we compare the results obtained for isolated benzene with available experimental results. The calculated ionization energies and HOMO-LUMO gap for isolated benzene are in very good agreement with experimental results. [21][22][23] Also, the calculated energies of benzene excitons (dark, bright, triplet and singlet) agree remarkably well (within 100 meV) with experimental data. [24][25][26] We show that the introduction of the substrates reduces the HOMO-LUMO gap for about 2 eV, and that this reduction weakly depends on the type of the substrate. A somewhat surprising result is that energies of the benzene excitons are barely affected by the presence of the substrate. This is because the substrate reduces the HOMO-LUMO gap which reduces the exciton energy, but at the same time the substrate weakens the excited electron-hole interaction which increases the exciton energy. We find that these two effect almost exactly cancel and exciton energies remains practically unchanged. We also find that all these effects can be simply explained by applying image theory to the screening of the electron and the hole, as was also theoretically observed in Ref. 8.
The spectra of molecular excitations is obtained in such a way that an external probe (electromagnetic wave or dipole) can excite only excitations in the molecule but not in the substrate. This enables us to analyze the molecular spectra as spectra of driven/damped harmonic oscillator where the external probe is the driving force of frequency ω, exciton in the molecule is a harmonic oscillator of frequency ω 0 , and the substrate is the source of damping with damping constant Γ. We find that only the singlet E 1 1u exciton (bright exciton) decays when the molecule is in the vicinity of the substrate, while all other excitons fail to couple with the substrate and remain infinitely sharp. In the vicinity of graphene the E 1 1u exciton decays into a continuum of π − π * interband electron-hole excitations where Γ ≈ 176 meV. In the vicinity of the Ag(jellium) surface the E 1 1u exciton decays faster and Γ ≈ 362 meV. We also find that there are no extra peaks in any molecular spectra. This means that the excitons do not interact with the 2D plasmon in doped graphene, or with the surface plasmon on the Ag(jellium) surface.
In Sec. II we present the general methodology used to solve the BSE for the 4-point polarizability matrix L kl i j (ω) and explain how to use the imaginary part of L kl i j (ω) to obtain the optical absorption and energy loss spectra for an arbitrary system.
In Sec. III the developed formulation is applied to derive the quasiparticle spectra, optical absorption and energy loss spectra of an isolated benzene molecule. In order to calculate quasiparticle spectra, optical absorption and energy loss spectra of the molecule near a substrate, we demonstrate that we only need to replace the bare Coulomb interaction V with V + ∆W, where ∆W is the substrate induced Coulomb interaction.
In Sec. IV we present results where the developed theoretical formulation is first used to calculate the ionization energies and HOMO-LUMO gap in the isolated benzene molecule, and then the formulation is used to obtain the exciton energies and spectra of excitations in benzene deposited on various substrates. This is followed by concluding remarks in Sec. V.

II. FORMULATION OF THE PROBLEM
In this section we shall first present the general method we use to solve the Bethe-Salpeter (BS) equation for the 4-point polarizability L(r 1 , r 2 ; r 1 , r 2 , ω) which is also called the twoparticle correlation function. 13 We shall then present how to obtain the optical absorption and energy loss spectra for an arbitrary system from L.

A. General theory
In absorption experiments a photon creates an electron and a hole. In the lowest approximation we can consider them as two independent particles, which leads to infinitely long lived electron-hole pairs that can be described as a product of two one-particle Green's functions. However, in reality, the situation is much more complex. Because of the electronelectron interaction the excited electron and hole can interact with other excitations in the molecule or they can annihilate or interact mutually. These are all responsible for the creation of their bound states called excitons. Therefore, to give a realistic description of optical absorption phenomena we have to calculate the two-particle Green's function G 2 . However, since in the Dyson's expansion of G 2 there are two possible annihilations leading to independent electron-hole motion, one of them should be subtracted from G 2 . In this way, the quantity describing the propagation of the coupled electron-hole pair is defined as 13 where is the exact two particle Green's function and is the exact one particle Green's function. Each argument in (1) represents a four-vector, e.g. 1 ≡ (r 1 , t 1 ). The 4-point polarizability (1) satisfies a Dyson-like equation of the form [11][12][13]16 L(1, 2; 1 , 2 ) = L 0 (1, 2; 1 , 2 ) + larizability has the form and represents the Bethe-Salpeter kernel. Here represents the exact Hartree energy and Σ XC (3, 4) represents the exact exchange-correlation self-energy, where the fourvector 1 + ≡ (r 1 , t 1 + δ) as δ → 0 + . To calculate the exact kernel (6) we have to know the exact self-energy Σ XC . However, this is not possible, so we have to make some approximations. The most frequently used approximation (in order to determine the BSE kernel) is the static-GW approximation 16 where W represents the exact statically screened Coulomb interaction. If we assume that W weakly depends on G then The functional derivative in (6) can be performed analytically, and the BSE kernel becomes Here we note that this static approximation is only justified when the frequency of the characteristic collective modes in the system are high enough to instantly screen the charge density fluctuation caused by electron-electron or hole-hole scattering in the system. Even though this is not always justified, we will use the approximation (10) because it enables us to transform equation (4) into frequency (ω) space.
The Green's functions in (13) can be obtained by solving the Dyson equation G(r, r , ω) = G 0 (r, r , ω) where G 0 is the independent electron Green's function and the self-energy can be separated into the Hartree part plus the exchange-correlation part The Hartree term is as shown in Feynman diagrams in Fig. 3. The first quantity that needs to be obtained to solve this complicated set of equations is the propagator W of the screened Coulomb interaction, since it appears in the BSE kernel (12) and is essential for the GW approximation (17). It is the solution of the equation (18) where the response function χ can be obtained from the 4point polarizability by the coordinate annihilation χ(r 1 , r 2 , ω) = L(r 1 , r 2 ; r 1 , r 2 , ω).
The first step is to solve the Density Functional Theory Kohn-Sham (DFT-KS) equations for the system in order to obtain the Kohn-Sham (KS) orbitals ψ i (r) and energy levels ε i . Using these KS states we can construct the independent electron Green's function The second step is to determine L(r 1 , r 2 ; r 1 , r 2 , ω) within the random phase approximation (RPA). We begin by inserting (20) into (13) to obtain the noninteracting 4-point polarizability: where is the occupation factor and N is the number of the occupied states. We can assume a similar expansion for L L(r 1 , r 2 ; r 1 , r 2 , ω) = i jkl where ensures that contributions to L only come from transitions between empty and filled states. That is, the summation indices should satisfy the following conditions: After inserting (23) and (21) into (11) and using the fact that in RPA the second term in the BSE kernel (12) (containing W) should be ignored, we obtain a matrix equation for L L kl i j (ω) = L kl,0 i j (ω) + where The matrix of the BSE-Hartree kernel has the form where the bare Coulomb interaction matrix elements are shown in Fig. 5. Here we introduced the two-particle wave functions and the factor of two in (28) comes from spin. By solving equation (26) we obtain matrix elements L kl i j (ω), and after inserting them into (23) we obtain the 4-point polarizability L at the RPA level. Now the coordinate annihilation r 1 = r 1 and r 2 = r 2 gives us the response function (19) and by inserting it into (18) we obtain the propagator of the dynamically screened Coulomb interaction The third step the is G 0 W 0 approximation. After replacing the Green's function in (17) by G 0 given by (20), performing the ω integration in (17), and using the spectral representation of the time ordered W, the exchange-correlation self energy becomes: 12 where the spectral function is defined as Note that the Hartree term (16) is already included in the KS energies ε i . The first term in (32) is the screened exchange (SEX) term and the second term is the Coulomb hole correlation (COH) term, so this kind of GW approximation is usually called the COHSEX approximation. After inserting (32) into a Dyson equation (14) we obtain a matrix equation for the Green's function where the self-energy matrix elements are defined by In most cases the off diagonal matrix elements of Σ i j weakly influence the diagonal matrix elements of G ii (ω) 12 and therefore we may safely neglect them. In this case (34) becomes a simple scalar equation and the solution is the quasiparticle Green's function This represents the propagation of a quasiparticle (electron or hole) in state i. The exchange-correlation self-energy is likewise In (35) we had to subtract the KS exchange correlation contributions V XC i from the KS energy levels ε i because at this level of approximation these contributions are included in Σ XC i (ω). After inserting (32) into (36) and using expression (31), the exchange-correlation self energy becomes and The first term in (37) represents the bare and induced exchange (Fock) energy and the second term represents the Coulomb hole correlation energy. This corresponds to the polarization energy shift due to an extra electron or hole in the system. It should be noted here that even though for unoccupied states the KS-Hartree term is exact, for occupied states it contains an incorrect self-interaction term. However, the selfenergy Σ XC i (ω) also contains an incorrect self-interaction Fock term V ii ii , equal in amount and with opposite sign, so that these two terms exactly cancel. Therefore, the self-interaction Fock term V ii ii is useful and should not be extracted from Σ XC i (ω). The poles of equation (35) represent the new quasiparticle energies ε QP i . In the quasiparticle approach the solution of (40) is close to the real axis and we can expand equation (40) around ε i as and the solution is where we introduced the normalization factor In the quasiparticle approach the imaginary part of Σ XC i can be neglected and the new, renormalized quasiparticle Green's function becomes The fourth and final step is solving the full Bethe-Salpeter equation (11). First we insert the quasiparticle Green's functions (44) into (13) to obtain the noninteracting quasiparticle 4-point polarizability L 0 QP . This has the same form as (21) except that ε i is replaced by ε QP i . Then, by using L 0 QP and repeating the RPA scheme (23-31) we obtain a new 4-point polarizability L QP and dynamically screened Coulomb propagator W. Now we can solve the Bethe-Salpeter equation (11) with the full kernel (12). After inserting L 0 QP , W(ω = 0), and the expansion (23) into (11) and (12) the BS matrix equation becomes interacting quasiparticle 4-point polarizability has the form and the BSE kernel consists of two terms This approximation is usually called the Time Dependent Screened Hartree Fock Approximation (TDSHFA). The first term in (47) is a matrix of the BSE-Hartree kernel given by (28) and (29), and the second term is a matrix of the BSE-Fock kernel given by After using (31) the BSE-Fock kernel can be expressed in terms of the 4-point polarizability matrix as L QP the Fock term does not contain the spin factor 2 because it does not allow a spin flip. This will be discussed later in more detail. The total BSE kernel is shown in Feynman diagrams in Fig. 7. The first term is the BSE-Hartree or RPA term. This represents the Coulomb interaction between electron-hole creation and annihilation. The second term is the BSE-Fock term. This represents the interaction between electron and hole mediated by the bare Coulomb interaction. The third term is the induced BSE-Fock term. This represents the interaction between electron and hole mediated by the induced Coulomb interaction. The second and the third term together is also called the screened BSE-Fock term.
After solving equation (45) by using expansion (23) and coordinate annihilation r 1 = r 1 and r 2 = r 2 , we obtain the BS response function

B. Optical absorption spectra
Optical absorption spectra is usually calculated by solving the BSE in the resonant approximation. This means that the frequency dependent part of (23), i.e. L kl i j (ω), is represented as a sum over harmonic oscillators with frequencies ω S representing the frequencies of the excitations in system, and amplitudes A S representing their oscillator strengths. This reduces the BS equation (11) into a standard eigenvalue problem for ω S and A S . This can then be inserted into the absorption spectra formula. 15,16 In this subsection we shall describe an alternative approach in which the optical absorption spectra is obtained directly from L kl i j (ω). In an optical absorption experiment the incident electromagnetic wave couples to the electronic excitations in the system and is partially absorbed. In linear response theory the power at which the external electromagnetic energy is absorbed in the system can be obtained from the expression where Π is the current-current response function of the system, while E ext and A ext are external electric field and vector potential respectively. The symbol ⊗ denotes convolution in real space.
We shall assume that the incident electromagnetic field is a plane wave of unit amplitude where e is the polarization vector. If we also assume that the external scalar potential is Φ ext = 0, this implies that ∂t . If the wavelength λ is much larger than the dimension of the illuminated system or the crystal unit cell, the dipole approximation can be applied and the absorption power becomes In the Coulomb gauge (∇ · A = 0), there is an instantaneous interaction mediated by the Coulomb interaction V and a transversal interaction that is retarded and mediated by photons. In small systems such as a molecule, the interaction between charge/current fluctuations mediated by photons is negligible compared to the Coulomb interaction. This allows us to describe all interactions inside the molecule by the instantaneous Coulomb interaction V and the interaction of the molecule with the environment by both interactions. In this case this is only interactions with photons described by A ext .
As a result, the current-current response function can be expressed in terms of the response function (50), except that now charge vertices should be replaced by current vertices (shown as squares in Fig. 8 After inserting (54) into (53) the absorption power becomes: where the form factors are

C. Energy loss spectra
Since some of the excitation modes, so-called dark modes, cannot be excited with incident electromagnetic waves, and here we are interested in all types of excitations in the molecule, we have to design some alternative probe which is able to excite the dark modes as well. The simplest choice is an external time dependent charge distribution. The rate at which an external charge distribution is losing energy to excitations in the system is given by 27 where Φ ext (r, t) is the time dependent potential produced by the external charge, χ is the density-density response function of the system, and ⊗ denotes convolution in the real space. If we assume a simple oscillatory time dependence the power loss becomes: We note that equation (59) is the longitudinal equivalent of equation (53). After using the expression for the response function (50) and definition (30) the power loss can also be written in term of matrix elements L kl i j (ω) as where the form factors are

III. APPLICATION TO MOLECULAR SPECTROSCOPY
In this section we shall show how we may apply the general procedure described in Sec. II to calculate the optical absorption and energy loss spectra of benzene in gas phase, deposited on graphene, and adsorbed on a metallic substrate.
A. Spectroscopy of gaseous benzene

Numerical solution of the BSE
The first step is to determine the molecular ground state electronic structure. The benzene Kohn-Sham orbitals ψ i (r) and energy levels ε i are obtained by using the plane-wave self-consistent field density functional theory (DFT) code (PWscf), within the Quantum Espresso (QE) package, 28 using the Perdew-Wang GGA (PW91) exchange and correlation (xc)-functional. 29 We model the benzene molecule using a periodically repeated 22.845 a 0 × 22.845 a 0 × 22.845 a 0 unit cell. Since there is no intermolecular overlap, the ground state electronic density is calculated at the Γ point only. For carbon and hydrogen atoms we used GGA-based ultra soft pseudo potentials, 30 and found the energy spectrum to be convergent with a 30 Ry plane-wave cutoff.
The benzene molecule has 30 valence electrons, which corresponds to 15 doubly occupied valence orbitals. For the 4-point polarizability calculation we use 60 orbitals, i.e. 15 occupied and 45 unoccupied orbitals. In Sec. IV we will show that the excitation spectrum is mostly defined by transitions inside the π − π * complex, or between occupied states a 2u , e 1g , e 1g and unoccupied states e 2u , e 2u , b 2g . 10 The KS wave functions are periodic and can be Fourier expanded as where G are reciprocal vectors and Ω is the normalization volume. In this case the two particle wave functions defined by (30) are also periodic and can be expanded as where the Fourier coefficients C j i (G) = 1 Ω dre −iG · r φ j i (r), with use of expansion (62) and definition (30), can be expressed in terms of coefficients C i (G) as This transformation enables higher numerical efficiency in the calculation of the bare Coulomb interaction matrix elements V kl i j which are the most frequently used quantities throughout the calculation.
Since we study a single isolated benzene molecule, we have to exclude the effect on its polarizability due to the interaction with surrounding molecules in the lattice. This is accomplished using the truncated Coulomb interaction 31 where Θ is the Heaviside step function, and R C is the range of the Coulomb interactions, i.e. the radial cutoff. Since we choose the lattice constant L = 22.845 a 0 to be more then twice the range of the benzene molecule's density, choosing the radial cutoff to be R C = L/2 ensures that the charge fluctuations created within the molecule produce a field throughout the whole molecule but do not produce any field within the surrounding molecules. The definition (65) is very useful because the Coulomb interaction remains translationally invariant. This leads to a simple Fourier transform After using the definition (29), the Fourier transform of the Coulomb interaction (66) and the expansion (63), the bare Coulomb interaction matrix elements become: where Ω cell = L 3 is the unit cell volume.
After we obtain the KS spectra ε i and Coulomb matrix elements (67) we can perform the RPA scheme (26)(27)(28). The matrix of the noninteracting 4-point polarizability L kl,0 i j (27) can be obtained directly from the KS energies ε i , and the BSE-Hartree kernel (28) directly from the matrix elements (67). This gives us the 4-point polarizability matrix L kl i j (ω), which is needed to obtain W i j (ω) through (38) and finally the exchange-correlation self energy Σ XC i (ω) using (37). From Σ XC i (ω) we now obtain corrected quasiparticle energies ε QP i from (42) and the matrix of noninteracting quasiparticle 4point polarizability L kl,0 i j QP (ω) from (46). Using L kl,0 i j QP (ω) and repeating the RPA scheme (23-31) we obtain a new 4-point polarizability L kl i j QP (ω) and BSE-Fock kernel given by (49). Finally, using L kl,0 i j QP (ω), the BSE-Hartree kernel (28), and the BSE-Fock kernel (49) we can solve the Bethe-Salpeter matrix equation (45) for L kl i j (ω). The fact that we use 15 occupied and 45 unoccupied orbitals for the calculation means that the dimension of the Bethe-Salpeter kernel matrix is 1350 × 1350. However, this does not depend on the number of plane waves used in the expansion of the Coulomb interaction matrix elements (67), which is an important advantage of our method. This accelerates matrix calculations and at the same time allows us to perform very accurate calculation of the Coulomb matrix elements. For example, in expression (67) we use a 30 Ry energy cut off, which corresponds to an expansion over 35000 plane waves. The disadvantage of this method is that the dimension of the matrix V kl i j increases with the number of occupied states N. This means the method is not computationally efficient for very large molecules.

Determination of the energy loss and optical absorption spectra
To be able to detect all types of electronic modes in the molecule, we simulate two kind of experimental spectroscopic methods. First is an optical absorption experiment. This is simulated by the absorption of a plane wave of light illuminating the molecule, as shown schematically in Fig. 9a). Second is an energy loss experiment. This is simulated by the energy loss of an oscillating dipole placed close the molecule, as shown schematically in Fig. 9b).
The benzene absorption spectra is obtained using expression (55), where the form factors (56) after using the expansion (62) become Electronic modes are also excited by an external potential Φ ext (t) where the rate at which the probe is losing energy is given by expressions (60) and (61). After inserting the Fourier transform of the external potential in the definition of the form factors (61), using the Fourier expansion (63), and the fact that the Φ ext (r) is a real function, i.e. Φ ext (−q) = [Φ ext (q)] * , the form factors become: We model the external probe as a dipole with dipole moment p placed at position R from the center of the benzene molecule, as shown in Fig. 9b). In this case, the Fourier transform of the external potential has the explicit form

Singlet and triplet excitons
Besides the spatial symmetry of the molecular electronic excitations, which determines whether the excitation will be dark or bright, there are also two classes of solution of the BSE with respect to spin. If the spin-orbit interaction is negligible compared with the electron-hole interaction, as we assume here, then each quasiparticle state has an additional quantum number associated with spin, i.e. up ↑ or down ↓. This has a simple impact on BSE. If spins of an excited electron-hole pair are parallel (e.g. spins of states i and j in Fig. 7 are both ↑) then the final state Hartree interaction can either leave the spin configuration unchanged or flip both spins into the opposite direction (i.e. spins of states k and l are both ↓). However, they will remain parallel. On the other hand, the Fock interaction always leaves the spin configuration unchanged, and this is why there is no factor of two in the Fock kernel (49). The BSE kernel is then simply Ξ = Ξ H − Ξ F , and any excitons created in this way have a spin singlet configuration.
If the external perturbation instead creates an electronhole pair with anti-parallel spins (e.g. spins of states i and j are ↑ and ↓ respectively) then, because of the orthogonality, such a pair cannot be annihilated and the Hartree interaction is inactive. The Fock interaction, responsible for the mutual electron-electron and hole-hole scattering, survives and it does not change the initial spin configuration, i.e. the final spins are still anti-parallel. The BSE kernel then consists of the Fock term only, i.e. Ξ = −Ξ F and this type of exciton forms a spin triplet configuration. 15 Both spin classes of excitons will be investigated in Sec. III B.

B. Spectroscopy of the benzene deposited on graphene
We next investigate the energy levels, optical absorption and energy loss spectra of benzene deposited on a graphene substrate, as illustrated in Fig. 10. Graphene and benzene planes are chosen to lie parallel to the xy-plane, i.e. have a normal parallel to the z-axis. Ground state electronic and crystal structure is obtained by structural relaxation using the second version of the van der Waals density functional of Lee et al. 32 and the exchange functional (C09) developed by Cooper. 33 This combination of functionals gives good agreement with experimental data for simi- lar systems. 34 The initial, most favourable, geometry is taken from Ref. 35. For obtaining ground state electronic density we used a supercell with dimensions a × √ 3 2 a × a where a = 27.906 a 0 . We employed a plane wave basis set with ultrasoft pseudopotentials as implemented in QE. 28,30 The kinetic energy cutoff for the plane waves was 40 Ry, and 500 Ry for the density. We applied an 8 × 8 × 1 Monkhorst-Pack special k-point mesh to sample the Brillouin Zone. For the average equilibrium separation between benzene and graphene we obtain z 0 ≈ 6 a 0 . This is the same as that reported in Ref. 35. Because of the large separation, the electronic densities of these systems do not overlap, as can be clearly seen in Fig. 11. This shows the ground state electronic density in the xz-plane along the dashed line denoted in Fig. 10.
The fact that the electronic densities do not overlap simplifies the impact of the graphene to benzene energy spectra and response function significantly. More specifically, since there is no inter-system electron hopping, the only modification comes from the additional screening caused by polarization of the graphene. In other words, interactions between charge fluctuations in the benzene have to be additionally screened because of the polarization of the graphene. This simply means that the bare Coulomb interaction has to be modified as follows: and also shown in Feynman diagrams in Fig. 12. Here ∆W is the induced dynamically screened Coulomb interaction of the graphene. 36 Consequently, the matrix elements (29) have to be modified to where ∆W kl i j (ω) = Ω cell dr 1 dr 2 φ j i (r 1 )∆W(r 1 , r 2 , ω)φ k l (r 2 ). We note that the integration in (74) is performed over the benzene superlattice unit cell volume Ω cell = L 3 . By doing this we avoid the graphene mediated intermolecular interaction, i.e. the influence of the surrounding molecules through polarization of the graphene is completely excluded. However, this makes the numerical computation more demanding.
The modification (73) is exact at the RPA level. Namely, the benzene 4-point polarizability L obtained by solving equations (26)(27)(28), including the modification (73), represents an exact RPA 4-point polarizability screened by graphene. However, if we want to calculate the substrate renormalized molecular 4-point polarizability beyond RPA, the modification (72) is no longer sufficient.
For example, when graphene is absent, the electron-hole interaction can be mediated directly by V or indirectly via molecular polarization V LV, as shown by the BSE-Fock kernel in Fig. 7. Introducing graphene requires the replacing of V byW. This induces extra electron-hole interaction channels, such as interaction via graphene polarization ∆W and interaction via mixed molecular-graphene polarization V L∆W, ∆WLV and ∆WL∆W. However, the higher multiplicity processes, e.g. V L∆WLV or ∆WL∆WL∆W, etc. may be safely neglected. When the perturbed system is a small zero-dimensional object, such as a molecule, and the substrate is large higher dimensional object, such as a surface, the dominant substrate induced electron-hole interaction is via single substrate polarization ∆W, and all higher order processes are negligible. Since we use the modification (72), we actually do include higher contributions such as V L∆W, ∆WLV and ∆WL∆W. these multiple processes could influence the electron-hole interaction if the perturbed system is of the same dimensionality as the substrate. In that case, both systems can form coupled modes and the propagator of the screened interaction W can no longer be separated into propagators of individual screened interactions. 37,38 The impact of the modification (72) to the G 0 W 0 exchange correlation self-energy (37)(38)(39) is similar to its impact on the BSE-Fock term. After applying the modification (73) to (38), we neglect the above mentioned higher order processes. This does not influence the result significantly because we find that the dominant substrate induced modification of molecular self-energy comes from the single substrate polarization term ∆W, as shown in Ref. 8. After applying the modification (73) to the expressions (37)(38)(39), the bare exchange self en-ergy includes an additional exchange self energy term which contains single polarization ∆W It is very important to note that the modification (73) is unable to generate the substrate induced correlation self energy term which contains the single polarization ∆W. Such an induced correlation term may be defined as where may represent a significant correction to the self energy and therefore must be included by hand.
In conclusion, the introduction of the substrate requires modification of the bare Coulomb interaction defined as (72). However, at the same time, the self energy (37) should be corrected to include the induced correlation term of (76) and (77). Therefore, the only task is to calculate the matrix elements of the induced Coulomb interaction (74).
The first step is to perform a Fourier transform in the xyplane where ρ = (x, y), Q = (Q x , Q y ) is a two-dimensional wave vector and G are graphene reciprocal vectors in the xy-plane.
In Fig. 11 we see that the equilibrium benzene-graphene separation is z 0 ≈ 6 a 0 . In Ref. 36 it is shown that the centroid of the induced density (density induced by the external point charge) is at z im ≈ 2 a 0 from the graphene center, as shown in Fig. 10. This means that charge fluctuations in benzene feel like the "external" graphene field in the region z, z > z im . This is the region where the graphene induced density is zero. In this region the spatial part of the Fourier transform (78) has the simple form 27,36,39 From (79) we see that the exponential factor cuts the higher G components. Since the average benzene graphene separation is z 0 ≈ 6 a 0 , it is sufficient to keep only the G = 0 component. This has the consequence that ∆W becomes isotropic in ρ and Q space, and the Fourier transform of the graphene field in the benzene region can be written simply as The graphene electronic excitation propagator D(Q, ω) contains the intensities of all (collective and single particle) electronic excitations in graphene. The details of the calculation of the propagator D(Q, ω) can be found in Refs. 27, 36, and 39. Here we use the same parameters employed in the calculation of D(Q, ω), except that the response function χ 0 is calculated using a 201 × 201 × 1 Monkhorst-Pack special k-point mesh, in order to have a finer Q-point mesh.
After inserting (80) into the expansion (78) and then together with the expansion (63) in the definition of the matrix elements (74), we obtain where the form factors are defined as and The summation in (82) is over benzene superlattice recip- and n x , n y , n z ∈ Z. The two-dimensional Q integration in (81) is performed using a 61 × 61 rectangular mesh and the cutoff wave vector Q C = 0.3 a 0 .
In conclusion, the calculation of the quasiparticle spectra, optical absorption and energy loss spectra of the deposited benzene can be performed following the same recipe as for the isolated benzene (as described in Sec. III A) except that the Coulomb matrix elements have to be renormalized as (73) and calculated by using expressions (81-83). Also, the correlation self energy has to be corrected by the term (76).

IV. RESULTS AND DISCUSSION
In this section we use the formalism developed in Sec. III to calculate the quasiparticle properties, energies, and spectra of excitons in benzene deposited on various substrates. We also compare our results with available experimental data.

A. Quasiparticle properties of benzene on a substrate
The quasiparticle energies for gaseous benzene are calculated directly by using G 0 W 0 scheme (37-43), additionally corrected by (75) and (76) for the deposited benzene. As explained before, in order to get accurate energy shifts in gaseous benzene we use 45 unoccupied states, i.e. 60 benzene states in total. Table I shows a comparison between benzene ionization and affinity energies obtained experimentally and by using the G 0 W 0 scheme. The experimental ionization energies are taken from Refs. 21 and 22, while the electron affinity energy is taken from Ref. 23.  To compare our quasiparticle energies with experimental values we must first determine the exact vacuum level. However, at the DFT level, we were unable to obtain an accurate vacuum level. For this reason, we aligned the energy of the LUMO state 2e 2u with the experimental affinity energy of 1.12 eV, and shifted all other levels accordingly. As can be seen in Table I, the result of this procedure is that all quasiparticle energies, incorrect in the DFT calculations, are now in satisfactory agreement with the experimental ionization energies. We also see that the renormalization of the HOMO level is very weak while the LUMO state is shifted as much as 5 eV upward. In fact, its energy becomes positive, indicating that this state is unbound.
The benzene excitons are mostly composed of transitions between the occupied state 1a 2u , doubly degenerate occupied state 1e 1g , doubly degenerate unoccupied state 1e 2u and unoccupied state 1b 2g , which form the benzene π−π * complex. The molecular orbitals corresponding to these states are shown in Fig. 13. The energies of the doubly degenerate occupied states 1e 1g and doubly degenerate unoccupied states 2e 2u define the benzene HOMO-LUMO gap.
Table II contains our G 0 W 0 results for the HOMO-LUMO gap in benzene deposited on various substrates, as compared with available experimental and theoretical results. We can see that our approach for the HOMO-LUMO gap of benzene in gas phase is in excellent agreement with experimental results. In this way we are able to verify the accuracy of our approach.
Unfortunately, there are no available experimental results for benzene deposited on graphene or metal substrates. We find that all the substrates reduce the HOMO-LUMO gap approximately by the same amount, between 2 and 2.4 eV. This can be explained in terms of a simple image theory shift. This can be obtained from (75) and (76) by setting ω = ε i and the summation index j = i. Then the induced self-energy becomes: Since the static induced potential ∆W ii ii (ω = 0) is always negative, because positive/negative charges feel an attractive force from their negative/positive image charge. Thus, the degenerate HOMO states i = 14, 15 are pushed up while the degenerate LUMO states i = 16, 17 are pushed down, reducing the HOMO-LUMO gap. The gap obtained using this simple model deviates less than 10% from results obtained using the full expressions (75) and (76). This means that the dynamical effect only slightly corrects the simple image theory result (84). The same conclusion has been reached in Ref. 8 where the authors theoretically investigated the quasiparticle properties of benzene deposited on graphite. The benzene HOMO-LUMO gap is then reduced by 3.2 eV as shown in the last row of Table II. Considering the similar results obtained for pristine and doped graphene we expected the same result for graphite as well. Namely, the external charge is only able to induce charge in the surface region, i.e. charge in the first graphite monolayer, so for an external charge graphene should be the same as graphite. Surprisingly, it instead turns out that the gap is about 1 eV larger than it is for graphite. Perhaps the graphite effective image plane z 0 is shifted outward compared to graphene. This would strengthen the molecular orbital screening shift, and may explain the difference in the gap.

B. Excitons in gaseous benzene
To ensure that our methodology for determining the energy of molecular excitons is accurate, we first calculate the energy of excitons in gaseous benzene. This is easier to compare with the numerous experimental results which are available in the literature.
The energies of the excitons in benzene in gas phase are determined from the positions of the peaks in the optical absorption spectrum calculated from expressions (55) and (68). Here the 4-point polarizability L kl i j is obtained by solving the BSE (45-49). The incident electromagnetic wave is chosen to be x polarized, i.e. in (52) we set e =x. The energies of the excitons are also determined from the peaks in the energy loss spectrum calculated from expressions (60), (70) and (71). We chose an asymmetric external charge distribution so it can excite excitons of all symmetries. The dipole is placed in the molecular plane but shifted by 4 a 0 in the +x direction. The dipole is similarly polarized in x direction, i.e. in (71) we put p =x and R = 4.0x. Table III shows the energies of different excitons in gaseous benzene. To be consistent with available literature, we identified and labelled all excitons as shown in the first row of Table III. Empty space in the table means that the corresponding spectra does not contain a corresponding exciton peak. In other words, the exciton cannot be excited by a corresponding external driver.
From Table III  and singlet B 1 2u and E 1 2g excitons are dark excitons, which cannot be excited by an electromagnetic field. This division to bright and dark excitons is consistent with optical and energy loss measurements. [24][25][26] This means that our method simu- lates both classes of experiments well. In Table III we also see that the energies of all types of excitons, except for the dark E 1 2g exciton, which is overestimated by 0.75 eV, are in excellent agreement with the experimental data. Altogether, this suggests that the theoretical methodology we have developed works quite satisfactorily, and can be applied to molecules on substrates.

C. Excitons in benzene on a substrate
In this subsection we analyze the excitation spectra of benzene when the molecule is deposited on various substrates. In addition to the energy of excitons, special attention will be paid to the decay mechanism for excitons into real excitations within the substrate. In order to include these real excitations, we use full the dynamic BSE-Hartree kernel (28), (73) and (74), while the BSE-Fock kernel (49) is also renormalized according to (73) but remains static. Including the dynamical effects in the BSE-Hartree kernel causes the calculation of BSE to become very computationally demanding. However, we note that the transitions between three occupied (n = 11, 14, 15) and three unoccupied (n = 16, 17, 18) states forming the π − π * complex 24,26 participate most dominantly in forming all significant excitons in benzene. Therefore, we restrict our calculations to the transitions inside the π−π * complex. This reduces the dimension of the BSE kernel matrix to only 18 × 18. Table IV shows energies of the different excitons after depositing on various substrates. For the separation between the molecular plane and the graphene plane or Ag jellium edge, we take the equilibrium value z 0 ≈ 6 a 0 .
It is interesting to note that the exciton energies are very weakly affected by the presence of the substrates. Similarly, the optical gap of benzene was previously found to be only weakly dependent on the height above a metal substrate. 18 This seems to be a general property of the optical gap of weakly bound molecules on substrates. 20 There are three dominant factors that define the exciton en-ergy. First is the quasiparticle energy shift. This changes the HOMO-LUMO gap and therefore the zero order exciton energy. Second is the fluctuation-fluctuation interaction present in BSE-Hartree kernel. This increases the exciton energy. Third is the screened electron-hole interaction present in the BSE-Fock kernel. This reduces the exciton energy. In our case, the substrate reduces the quasiparticle HOMO-LUMO gap by more than 2 eV, as shown in Table II. This reduces the exciton energy, and it barely influences the fluctuation-fluctuation interaction in the molecule. On the other hand, the substrate significantly weakens the electronhole interaction, which increases the exciton energy. The latter can be explained by using simple image potential theory. When a substrate is present, the molecular electron interacts with the molecular hole but also with its negatively charged image. This reduces the attractive electron-hole interaction and increases the exciton energy. Therefore, there is a competition between the HOMO-LUMO gap reduction and weakening of the electron-hole interaction.
In our case, substrate induced electron(LUMO)hole(HOMO) interaction, which can be approximated by the matrix elements −∆W j j ii (ω = 0) > 0, almost exactly cancels the HOMO-LUMO gap reduction. Using (84), this can be approximated as where i = 14, 15 and j = 16, 17. Therefore, the exciton energies are indeed substantially affected by different mechanisms. However, these mechanisms cancel each other and the exciton energy remains almost unchanged. The theoretical model developed here allows us to analyze the molecule/substrate spectra in analogy to the spectra of a driven/damped harmonic oscillator. Namely, the calculation is performed in such a way that the external electromagnetic wave or dipole (driving force) can induce current or charge in the molecule (harmonic oscillator) but not in the substrate (damping source). However, the molecule interacts with the substrate and it can excite plasmons (leading to extra peaks in the molecular spectra) or electron-hole excitations (influencing final exciton width) in the substrate. The inverse exciton width represents the decay rate of the initially excited exciton.
Here we note that the molecule is a zero-dimensional object. This means there is no translational invariance within it, and Q is not a valid quantum number. This also means that the exciton at fixed frequency ω can decay into any of the electron-hole excitations with any momentum transfer Q. Since we use the static BSE-Fock kernel in our theoretical model, triplet exciton are always sharp peaks and cannot decay into substrate excitations. For this reason, we investigate the decay of singlet excitons.
To be able to distinguish a substrate induced exciton decay from intrinsic decay, we choose an intrinsic exciton decay which is very small, namely 1 meV. Since the interaction between the molecule and substrate is quite weak, the probability of transitions between the E 1u exciton of benzene and the substrate is well described by Fermi's Golden Rule 40 . As the FIG. 14: (Color online) Optical absorption and energy loss spectra of benzene in vacuum (black solid line), on pristine graphene (red dashed line), and in the vicinity of a Ag(jellium) surface (blue dashdotted line). For separation between molecular plane and graphene plane or jellium edge we take the equilibrium distance z 0 ≈ 6.0 a 0 . Black dotted lines are Lorentzian fits to the E 1u optical adsorption spectra.
substrates considered are in the wide-band limit, the final density of states is a Lorentzian distribution, so that the absorption spectra has the form Here H is the coupling matrix element between the molecule and substrate, ω 0 is the exciton energy, and τ ≈ /Γ is its lifetime. 41 Assuming |H| 2 is only weakly dependent on the energy, we may fit (86) to the calculated spectra to estimate the inverse lifetime of the exciton, i.e. the decay width Γ. Fig. 14 presents the optical and energy loss spectra of gaseous benzene (black solid line) and of benzene deposited on pristine graphene (red dashed line). Upper panels show the optical absorption spectra showing only those excitons that can be excited by the external electromagnetic wave. This is why we can notice the absence of triplet E 3 1u and E 3 2g excitons and singlet B 1 2u and E 1 2g excitons. Lower figures show the spectra of excitons excited by a dipole. These spectra show all types of dark and bright excitons. We see that the substrates do not generate any extra peak in the spectra. This means that the excitons do not interact with the π − π plasmons in pristine graphene. Also, as mentioned previously, we see that the exciton energies are almost unaffected by the substrate and that they are just slightly shifted toward higher energies. From energy loss spectra (lower panels of Fig. 14) we see that only the singlet E 1 1u excitons obtain a final width, while the dark B 1 2u and E 1 2g excitons remain sharp. By observing the low momentum/energy (Q ≤ 0.3 a 0 and ω < 10 eV) graphene spectra 39 we can see a wide interband π − π * electron-hole continuum and broad π − π * plasmon. So obviously the E 1 1u exciton decays into all types of π − π * excitations in graphene. This does not exclude the possibility that it can also interact with the weak graphene π−π * plasmon. By fitting E 1 1u optical spectra to a Lorentzian (86), we obtain Γ ≈ 174 meV as provided in Table IV. The spectra of benzene deposited on doped graphene (ε F = 1 eV) is almost identical to that for pristine graphene and is not shown here. This is expected because even though the doping induces extra intraband electron-hole continuum, it appears at energies (ω < 2 eV) 39 lower than the benzene exciton energies (ω 0 > 3 eV). There is also an absence of extra peaks in the spectra. This means that the excitons do not interact with the 2D plasmon of doped graphene. The decay width of the E 1 1u exciton is slightly reduced in the vicinity of doped graphene to Γ ≈ 162 meV. This could be because in the doped graphene the interband π − π * continuum is shifted to higher energies. This reduces their intensity at the exciton energy ω E 1 1u ≈ 7 eV. From this we can conclude that the excitons in larger molecules, as e.g. in terrylene C 30 H 16 where ω 0 ≤ 3 eV, would be more strongly influenced by graphene doping.
Fig. 14 also shows the optical and the energy loss spectra of benzene deposited on a Ag (jellium) surface (blue dasheddotted line). For the separation between molecular plane and jellium edge we chose z 0 ≈ 6 a 0 . We see that the singlet E 1 1u exciton significantly decays into excitations within the metal. By fitting its optical spectrum to a Lorentzian, we obtain for its width Γ ≈ 362 meV. This is because in the jellium metal there are many interband electron-hole channels into which it can decay. It is interesting to note that even in this case, when the phase space of the electron-hole excitation becomes very rich, the dark excitons (B 2u and E 2g in the lower panels of Fig. 14) remain sharp.
We conclude that only the bright exciton E 1 1u decays into real electron-hole excitations in the substrate. On the other hand, the dark excitons B 1 2u and E 1 2g do not interact with any type of real excitations in semi-metallic or metallic substrates and remains in well defined eigenmodes.

V. CONCLUSIONS
In this manuscript we presented how the molecular optical and energy loss spectra can be obtained directly from the dynamical 4-point polarizability matrix L kl i j (ω), which is the solution of the BS matrix equation. We demonstrated that the inclusion of substrates requires minimal intervention to the presented formulation. This implies that everywhere throughout the BSE-G 0 W 0 scheme the bare Coulomb interaction V can be replaced by the dynamically screened Coulomb interaction W(ω) = V + ∆W(ω), where ∆W(ω) is the substrate induced Coulomb interaction. This formulation has been successfully applied to the calculation of the quasiparticle energy levels and exciton energies in the isolated benzene molecule. The method has then been applied to calculate the electronic structure and excitations in benzene deposited on pristine and doped graphene and in benzene in the vicinity of a Ag(jellium) surface. It is shown that the substrates cause a reduction of the quasiparticle HOMO-LUMO gap (by about 2 eV), which weakly depends on the type of substrate. We have also shown that the energy of all excitons in the isolated molecule remain relatively unchanged when the molecule is deposited on a substrate.
By using an image theory based argument, we note that the exciton energies are under the influence of two mechanisms which tend to cancel each other out. The substrate reduces the quasiparticle HOMO-LUMO gap which reduce exciton energy. However, at the same time, the induced image-electron or image-hole weakens the electron-hole interaction which raises the exciton energy.
We pay special attention to the investigation of the interaction of different types of excitons with real electronic excitations in the substrate. It is noted that only the optically active E 1u exciton decays into the electron-hole excitations in the substrates. However, it does not couple to any plasmons in doped graphene or within the metallic surface.
Coupling to electronic excitation in the substrate causes a Lorentzian broadening of the E 1u exciton, whose width is Γ ≈ 174 meV for pristine graphene and Γ ≈ 362 meV for metal surfaces as substrate. We have also noticed that the exciton quenching could be tuned by graphene doping.
Although this effect is not observed in benzene, it should exist for larger π-conjugated complexes, such as terrylene C 30 H 16 . There the molecular excitons falls in the gap between upper and lower edges of the doped graphene intra and interband electron hole continuum respectively. Now that the developed formulation has been successfully tested, it has the potential to be applied to more computationally demanding and technologically interesting molecular systems.