Polarization-induced renormalization of molecular levels at metallic and semiconducting surfaces
On the basis of first-principles GW calculations we study systematically how the electronic levels of a benzene molecule are renormalized by substrate polarization when physisorbed on different metallic and semiconducting surfaces. The polarization-induced reduction of the energy gap between occupied and unoccupied molecular levels is found to scale with the substrate density of states at the Fermi level (for metals) and substrate band gap (for semiconductors). These conclusions are further supported by GW calculations on simple lattice models. By expressing the electron self-energy in terms of the substrate’s joint density of states we relate the level shift to the surface electronic structure thus providing a microscopic explanation of the trends in the GW calculations. While image charge effects are not captured by semi-local and hybrid exchange-correlation functionals, we find that error cancellations lead to remarkably good agreement between the GW and Kohn-Sham energies for the occupied orbitals of the adsorbed molecule.
Solid-molecule interfaces are central to a number of important areas of physics and chemistry including heterogeneous catalysis, electrochemistry, molecular- and organic electronics, and scanning tunneling spectroscopy(1); (2); (3); (4). Most of our current understanding of level alignment at interfaces builds on effective single-particle descriptions such as the Kohn-Sham scheme of density functional theory (DFT)(5). Within such theories the energy levels of a molecule close to a surface are determined by hybridization, charge-transfer, and interface dipole fields – all properties of the static mean field potential defining the single-particle Hamiltonian. On the other hand, from photo-emission and electron transport measurements it is well known that the dynamic polarizability of the molecule’s local environment can have a large influence on the level positions(6); (7); (8); (9); (10). Such polarization effects, which are induced by changes in the charge state of the molecule, are not captured by available single-particle descriptions.
Many-body perturbation theory provides a systematic method to obtain the true single-particle excitations [sometimes referred to as addition/removal energies or quasiparticle (QP) energies] from the Green function of the system. In the GW approximation the electron self-energy is written as a product of the (non-interacting) Green function and a dynamically screened Coulomb interaction, (11); (12). It is instructive to compare this to the bare exchange self-energy given by , where is the unscreened Coulomb interaction. It is well known that the Hartree-Fock (HF) eigenvalues correspond to energy differences between the -particle groundstate and the unrelaxed -particle Slater determinants (Koopmans’ theorem). The effect of replacing with the screened and frequency dependent is two-fold: it introduces correlations into the many-body eigenstates, and it includes the response of the other electrons to the added electron/hole, i.e. relaxation effects. For a molecule at a surface, the latter effect is particularly important as it incorporates the attractive interaction between the added electron/hole and its induced image charge, into the QP spectrum.
Recent experiments on molecular charge transport have renewed the interest for theoretical modeling of polarization-induced level renormalization. First-principles GW calculations for a benzene molecule on graphite(13) as well as CO on NaCl/Ge(001)(14) have demonstrated significant reductions of the molecular energy gap due to image charge effects. Model GW calculations have been used to elucidate the qualitative features of the effect across different bonding regimes(15). Classical electrostatic models of various complexity have been developed to correct energy levels obtained from single-particle calculations.(16); (17); (18); (19).
In this work, we present a systematic study of image charge-induced renormalization at a range of different surfaces taking both a classical and quantum many-body viewpoint. We have performed DFT calculations with local density approximation (LDA) and hybrid (PBE0) exchange-correlation functionals as well as GW calculations for a benzene molecule weakly physisorbed on the metals Li, Al, Ti, Rh, Pt, and the semiconductors/insulators TiO, BaO, MgO, CaO, and NaCl. The results for the HOMO-LUMO gap of benzene are shown in Fig. 1. While LDA and PBE0 yields a substrate independent HOMO-LUMO gap, the GW gaps are reduced from the gas phase value by an amount which depends on the polarizability of the surface. For all systems, we find that the dependence of the QP gap on the distance to the surface can be described by a classical image charge model. However, the model parameters are sensitive to the microscopic details of the system and this limits the usefulness of the classical model in pratice. By evaluating the GW self-energy to second order we obtain a simple analytic expression which relates the level shift to the substrate’s joint density of states weighted by Coulomb interaction matrix elements. The model suggests that the HOMO-LUMO gap should scale with the substrate band gap (for semiconducting surfaces) and density of states at the Fermi level (for metallic surfaces). This trend is verified for the first-principles results and is further supported by GW calculations for simple lattice models. Finally, we analyze the deviation between the DFT and GW results in more detail. We find that the occupied Kohn-Sham levels obtained with LDA (PBE0) are in very good agreement with the GW results for benzene adsorbed on the metallic (semiconducting) surfaces, and we show that this is a result of significant error cancellation in the LDA/PBE0 approximations.
The paper is organized as follows. In Sec. II we outline the methodology used for the first-principles and model GW calculations. In Sec. III we investigate to what extent the first-principles GW results can be explained by a classical image charge model. In Sec. IV we derive a simple analytical expression for the polarization-induced level shift and show that it explains the main trends in both the first-principles as well as the model calculations. At the end of the section we analyze the description of occupied and unoccupied levels separately and discuss the effect of error cancellations in the DFT results. We conclude in Sec. V
ii.1 Ab-initio GW calculations
To model the solid-molecule interfaces we use a slab containing four atomic layers of the substrate in the experimentally most stable phase, and a benzene molecule lying flat above the surface followed by 12Å of vacuum. The benzene molecule is not relaxed on the surface but is fixed in its gas phase structure at a distance from the surface. An example of a supercell is shown in Fig. 1(a) for the case of benzene on NaCl(001). The number of atoms included in the supercell per atomic layer is 9 for Al, Rh, Pt, Ti; 12 for Li and TiO; and 16 for NaCl, MgO, CaO and BaO. This corresponds to distances between periodically repeated benzene molecules in the range 8.1 to 9.9 Å. All DFT calculations have been performed with the PWSCF code (27) which uses norm-conserving pseudopotentials (28). For exchange-correlation functionals we have used the local density approximation (29) as well as the PBE0 hybrid functional(30); (31). The Brillouin zone (BZ) was sampled on a 4x4x1 k-point mesh, and the wavefunctions were expanded with a cut-off energy of 40 Hartree.
In the GW(LDA) method one obtains the QP energies from the linearized QP equation
where and are LDA eigenstates and eigenvalues, and
The self-energy, , is evaluated non-selfconsistently from the single-particle Green function, i.e. , with . It is customary to use the random phase approximation for the screened interaction, i.e. with .
We have performed the GW calculations with the Yambo code (32) using the LDA wavefunctions and eigenvalues from the PWSCF calculations as input. The plasmon pole approximation has been applied with a frequency of 1 Hartree (the HOMO and LUMO energies of benzene change by less than 0.05 eV when the plasmon frequency is varied between 0.5 and 2.0 Hartrees). In the calculation of the self-energy we included a minimum of 200 empty states. We have checked that calculations are converged with respect to slab thickness, lateral supercell size, k-point mesh, all energy cut-offs, and that we reproduce the results previously reported in Ref. (13) for benzene on graphite at Å.
ii.2 Model GW calculations
In addition to the first-principles GW calculations, we have performed (self-consistent) GW calculations for two lattice models representing a metal-molecule and semiconductor-molecule interface, respectively. The model Hamiltonians contain three terms
describing the solid (metal or semiconductor), the molecule, and their mutual interaction, see Fig. 2. A metallic substrate is modeled by a semi-infinite tight-binding (TB) chain (we suppress the spin for notational simplicity),
A semiconducting substrate is modeled by
where refers to conduction and valence bands, respectively.
The molecule is represented by its HOMO and LUMO levels, i.e.
where e.g. , is the number operator of the HOMO level.
Finally, the interaction between the molecule and the terminal site(s) of the substrate TB chain(s) is described by
where is the number operator of the molecule. Note that since polarization of a semiconductor occurs via transitions between valence and conduction bands, only the interaction terms of the form given above contribute to the image charge effect (this will become clear in Sec. IV.1).
We set corresponding to a half filled band for the metal. We choose and so that the molecule contains exactly two electrons ( in the middle of the HOMO-LUMO gap). We consider the limit of zero hybridization between the solid and molecule so that interaction between the solid and molecule occurs only via the non-local . The model neglects interactions within the TB chain and between the molecule and interior TB sites (). These approximations are, however, not expected to influence the image charge physics described by the model in any qualitative way.
We obtain the Green function of the molecule from
where the Hartree potential due to has been absorbed in . The GW self-energy is calculated fully self-consistently using a recently developed GW scheme for quantum transport(20). The renormalized molecular QP levels are obtained as peaks in the spectral function .
Iii Classical Theory
In this section we investigate to what extent the GW results of Fig. 1 can be described by a classical image-charge model.
The electrostatic energy of a point charge, , located in vaccum at position above a polarizable medium filling the half-space , is given by (in atomic units)
The size of the image charge is , where is the relative dielectric constant of the medium(21). In 1973 Lang and Kohn showed that the energy of a classical point charge above a quantum jellium surface follows Eq. (8) with (corresponding to as expected for a perfect metal), with the image plane, , lying 0.5-0.9 Å outside the surface depending on the electron density(22). More recently, ab-initio GW calculations have found the same asymptotic form for the potential felt by an electron outside a metallic surface(21); (23); (24). From this it seems reasonable to conclude that the assymptotic position of the electronic levels of a molecule outside a surface would also follow the image potential of Eq. (8). This is, however, only true for the unoccupied levels whereas the occupied levels experience a shift in the opposite direction, i.e. the shift is upward in energy as the molecule approaches the surface. This is because the occupied levels represent the negative of the energy cost of removing an electron from the molecule. Similarly it has been found that the image potential leads to band gap narrowing at semiconductor-metal interfaces.(25); (26).
To test whether the gap reductions obtained in the GW calculations can be described by the classical image charge model we have fitted Eq. (8) to the calculated HOMO-LUMO gap for Å. In Fig. 2 we show the result of the fit for three systems (the fit is equally good for the other systems). The best-fit values for the effective image plane and the dielectric constant are given in table LABEL:table1.
As can be seen is generally smaller than the experimental optical dielectric constant of the bulk, . This is expected since the latter gives the long-range response of the bulk while probes the local response at the surface. Part of the discrepancy between and is clearly due to geometric effects. By taking the surface geometry into account, as done in Ref. (18), better estimates of can be produced from . On the other hand, electronic effects due to the local atomic structure of the surface cannot be captured by a classical model. For example, the BaO(111) surface is metallic due to surface states, and thus while . Similarly, impurities, defects, and surface roughness are expected to influence the local dielectric properties of the surface.
According to the classical image charge model all the molecular levels should experience the same shift (the sign of the shift being different for occupied and unoccupied levels). However, we have found that the best-fit values for and obtained by fitting the HOMO and LUMO levels speparately, are in general different – most notably for the metallic surfaces. This observation, which is discussed in more detail in Sec. IV.3, shows that the shape of molecular orbital also influences the size of the polarization-induced shift.
Iv Microscopic theory
In this section we first consider the GW self-energy for a molecule interacting with a surface to second-order in the electron-electron interaction. This leads to a simple microscopic model for image charge renormalization which relates the shift of molecular levels to the electronic structure of the surface, and explains general trends of the first-principles and model GW calculations. In the last section we consider the HOMO and LUMO levels separately and explain how error cancellations in semi-local exchange-correlation functionals can explain the surprisingly good agreement found between LDA eigenvalues and GW QP energies for the occupied levels of benzene on metallic surfaces.
iv.1 Second-order expansion
In quantum many-body theory, the effect of substrate polarization on the energy levels of a molecule enters the Green function via a self-energy operator. In general, the GW self-energy can be written symbolically as
where is the Green function of the non-interacting (Kohn-Sham) Hamiltonian, and is the polarization bubble. The first-order term, , is simply the static exchange potential while the remaining terms account for correlations and dynamic screening. In the following we consider the second-order term, explicitly. This corresponds to approximating the reponse of the substrate by its non-interacting response, .
For sufficiently large surface-molecule separations (Å) we can neglect hybridization effects, and the non-interacting eigenstates of the combined system can be taken as the eigenstates of the isolated molecule and surface. We denote these eigenstates by (’’ for adsorbate) and , respectively. To see how a given electronic level, , is renormalized by polarization processes in the substrate we consider the (time-ordered) matrix element , given by
The Feynman diagram corresponding to is shown in Fig. 4(a). The polarization and Coulomb matrices are given by
where we have defined the interaction strength,
Note that is simply the joint density of states (JDOS) of the substrate, shifted by , and weighted by the Coulomb matrix elements. The physically relevant retarded self-energy is readily obtained from Eq. (13)
where denotes the Cauchy principal value. Now, the renormalized QP energy can be obtained from the equation (neglecting off-diagonal terms)
A graphical solution to the QP equation is illustrated in Fig. 4(b,c) for the case of an occupied molecular level interacting with a metal or semiconductor surface, respectively.
From Eq. (14) it follows that the image charge effect does not broaden the molecular level because . We also note that the level shift is independent of the absolute value , and that the effect of changing the sign of is to change the sign of the level shift. These properties are all in line with the classical theory.
In the limit where varies little with and , is simply proportional to the shifted JDOS [the ”generic” cases illustrated in Figs. 4(b,c)]. In this case the level shift is simply determined by the form of the JDOS. For a metal, the JDOS raises linearly at with a slope given by the metal’s DOS at . This suggests that the level shift should increase with the substrate DOS at the Fermi level. For a semiconductor, the JDOS raises smoothly at , suggesting that the level shift should decrease with . In the following section we investigate these relations for the model and first-principles calculations. We mention that the second order approximation discussed above may not always provide a good description of the full GW self-energy. However, as we will show in the next section, it explains qualitatively the trends in GW calculations.
iv.2 Dependence of level shift on surface electronic structure
In Fig. 5 we show the HOMO and LUMO levels of the lattice models calculated with the HF and GW approximations. In all plots we vary one parameter of the model while keeping the remaining parameters fixed(34).
The upper panels refer to a metallic substrate and show the dependence on the levels on the interaction strength and the intra-chain hopping parameter . Note that the latter is inversely proportional to the projected density of states (DOS) of the terminal site evaluated at . The lower panels refer to a semiconducting substrate and show the dependence of the levels on and the substrate gap, . The HF eigenvalues are clearly independent of the non-local interaction between the molecule and substrate. This can be understood from Koopmans’ theorem which states that the HF eigenvalues do not include the electronic relaxations of the substrate induced by the extra electron/hole in the molecule. In contrast the GW levels vary in the way predicted by the simple model discussed in the previous section: The polarization-induced reduction of the HOMO-LUMO gap is stronger for larger as well as for larger substrate DOS at for the metals and smaller substrate band gap for the semiconductors. A more detailed discussion of level renormalization based on the lattice model for metallic substrates, including the case of strong metal-molecule hybridization, can be found in Ref. (15).
In Fig. 6(a,b) we plot the GW gaps from Fig. 1 versus the LDA band gap and DOS at for the semiconducting and metallic substrates, respectively. For the semiconductors the reduction of the HOMO-LUMO gap clearly correlates with . This indicates that the interaction strength, i.e. the matrix elements of Eq. (14), do not differ too much from one surface to another. For the metals, the HOMO-LUMO gap seems to scale with the metal’s DOS at . However, we note that Li(001) and BaO(111) deviate from the general trend followed by the other metals. This can be explained by the larger extend of the metallic wavefunctions of these systems into the vacuum region, which in turn leads to larger matrix elements. Indeed, Fig. 6(c) shows the average electron density evaluated in a plane lying Åabove the surface in the absence of the benzene molecule. The density outside the Li(001) and BaO(111) surfaces is significantly larger than for the other surfaces which on the other hand have quite similar densities.
iv.3 DFT eigenvalues and error cancellation
In Fig. 7 we plot the energies of the HOMO and LUMO levels of benzene at Å. For each surface, we have shifted the LDA, PBE0, and GW levels by the same amount so that the LDA HOMO is aligned with the HOMO in the gas phase. We note that the effect of substrate polarization is very similar for the GW HOMO and LUMO levels which are shifted up and down, respectively, by almost the same amount. This is indeed expected from the classical image charge model. Significant deviations from this trend are, however, seen for Li(001) and BaO(111). We ascribe this to the more extended nature of the metallic states on these surfaces which reduce the validity of the point charge approximation and can introduce differences between the and matrix elements.
Overall, the LDA and PBE0 eigenvalues for the HOMO are in better agreement with the GW QP energies than is the case for the LUMO. Moreover there is a general trend that the LDA eigenvalues come closer to the GW energies as we move from the insulating to the metallic surfaces. In fact, the LDA HOMO level is almost on top of the GW level on the metallic surfaces. This trend is clearly a result of sigificant error cancellation in the LDA. Indeed, it is well known that semilocal exchange-correlation functionals overestimate (underestimate) occupied (empty) molecular levels due to self-interaction effects. At the metallic surfaces this error is compensated by the missing image charge correction. PBE0 gives better estimates for the free molecule where it opens up the LDA HOMO-LUMO gap due to partial removal of self-interaction errors. In this case, the cancellation between the missing image charge effect and the remaining self-interaction error results in very good agreement between PBE0 and GW for the HOMO level on the semiconducting surfaces.
The cancellation between self-interaction errors and missing polarization effects will always be present in hybrid- and semilocal approximations. However, the relative size of the two contributions will in general depend on the shape of the molecule, its orientation with respect to the surface, the molecule-surface distance, and the type of substrate.
We have presented GW calculations for a benzene molecule physisorbed on different metallic and semiconducting surfaces. Upon physisorption the molecule’s HOMO-LUMO gap is reduced from its gas phase value due to dynamic polarization of the substrate. It was shown that a classical image charge model captures the qualitative features of the effect while the magnitude of the level shift is sensitive to the detailed atomic structure of the surface. In particular the presence of metallic mid-gap state at the surface of a semiconductor can have a large influence on the local response of the surface. Both local and hybrid exchange-correlation potentials fail to account for the polarization effects yielding Kohn-Sham eigenvalues of physisorbed benzene which are independent of the substrate. Nevertheless we found that a cancellation between self-interaction errors and missing image charge effects in the LDA leads to a very good agreement between LDA and GW energies for the occupied states of benzene on metallic surfaces. Similar conclusions were reached for the PBE0 energies on semiconducting substrates. Finally, we have derived a simple second-order approximation to the GW self-energy which expresses the polarization-induced shift of a molecular level in terms of the substrate’s joint density of states weighted by Coulomb interaction matrix elements. This model was used to explain general trends in the first-principles results, namely the scaling of the benzene’s HOMO-LUMO gap with the substrate density of states at (for metals) and the substrate band gap (for semiconductors).
Our results clearly demonstrates the importance of non-local correlations for the electronic levels at solid-molecule interfaces. We expect this to have important implications for the theoretical modelling of electron transport in organic and single-molecule devices.
KST and CR acknowledge support from the Danish Center for Scientific Computing. The Center for Atomic-scale Materials Design (CAMD) is sponsored by the Lundbeck Foundation. AR and JMGL acknowledge funding by the Spanish MEC (FIS2007-65702-C02-01), ”Grupos Consolidados UPV/EHU del Gobierno Vasco” (IT-319-07), e-I3 ETSF project (Contract Number 211956) and ”Red Española de Supercomputación”.
- J.K. Norskov, T. Bligaard, J. Rossmeisl, and C.H. Christensen, Nature Chemistry 1, 37 (2009)
- K. Moth-Poulsen and T. Bjornholm, Nature Nanotechnology 4, 551 (2009)
- A. Nitzan and M. A. Ratner, Science 300, 1384 (2003).
- W. A. Hofer, A. S. Foster, and A. L. Shluger, Rev. Mod. Phys. 75, 1287 (2003)
- F. Flores, J. Ortega and H. Vazquez Phys. Chem. Chem. Phys. 11, 8658 (2009)
- P. D. Johnson and S. L. Hulbert Phys. Rev. B 35, 9427 (1987)
- S. Kubatkin et al. Nature 425, 698 (2003)
- R. Hesper, L. H. Tjeng and G. A. Sawatzky, Euro. Phys. Lett. 40, 177 (1997)
- J. Repp et al. Phys. Rev. Lett. 94, 026803 (2005)
- A. Kahn, N. Koch, and W. Gao, J. Polym. Sci., B Polym. Phys. 41, 2529 (2003)
- F. Aryasetiawan and O. Gunnarsson, Rep. Prog. Phys. 61, 237 (1998)
- G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. 74, 601 (2002).
- J. B. Neaton, Mark S. Hybertsen, and Steven G. Louie Phys. Rev. Lett. 97, 216405 (2006)
- C. Freysoldt, P. Rinke, and M. Scheffler, Phys. Rev. Lett. 103, 056803 (2009)
- K. S. Thygesen and A. Rubio Phys. Rev. Lett. 102, 046802 (2009)
- S. Y. Quek et al., Nano Lett. 7, 3477 (2007)
- D. J. Mowbray, G. Jones, and K. S. Thygesen J. Chem. Phys. 128, 111103 (2008).
- K. Kaasbjerg and K. Flensberg, Nano Lett. 8, 3809 (2008)
- R. Stadler, V. Geskin, and J. Cornil Phys. Rev. B 79, 113408 (2009)
- K. S. Thygesen and A. Rubio, Phys. Rev. B 77, 115333 (2008)
- M. Rohlfing, N.-P. Wang, P. Krüger, and J. Pollmann, Phys. Rev. Lett. 91, 256802 (2003)
- N. Lang and W. Kohn, Phys. Rev. B 7, 3541 (1973)
- A. G. Eguiluz, M. Heinrichsmeier, A. Fleszar, and W. Hanke, Phys. Rev. Lett. 68, 1359 (1992)
- I. D. White, R. W. Godby, M. M. Rieger, and R. J. Needs, Phys. Rev. Lett. 80, 4265 (1998)
- J. C. Inkson, J. Phys. C, 6, 1350 (1973).
- J. P. A. Charlesworth, R. W. Godby, and R. J. Needs Phys. Rev. Lett. 70, 1685 (1993)
- S. Baroni et al., QUANTUM-ESPRESSO package, 2005 (http://www.quantum-espresso.org/).
- M. Fuchs and M. Scheffler, Comput. Phys. Commun. 119, 67 (1999)
- J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992)
- M. Ernzerhof and G. Scuseria, J. Chem. Phys. 110, 5029 (1999)
- C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999)
- A. Marini, C. Hogan, M. Grüning, D. Varsano Comp. Phys. Comm. In press (2009). http://www.yambo-code.org/
- For the metals: , , . For the semiconductor: , , , .