Energy Ordering of Molecular Orbitals.

Orbitals are invaluable in providing a model of bonding in molecules or between molecules and surfaces. Most present-day methods in computational chemistry begin by calculating the molecular orbitals of the system. To what extent have these mathematical objects analogues in the real world? To shed light on this intriguing question, we employ a photoemission tomography study on monolayers of 3,4,9,10-perylene-tetracarboxylic acid dianhydride (PTCDA) grown on three Ag surfaces. The characteristic photoelectron angular distribution enables us to assign individual molecular orbitals to the emission features. When comparing the resulting energy positions to density functional calculations, we observe deviations in the energy ordering. By performing complete active space calculations (CASSCF), we can explain the experimentally observed orbital ordering, suggesting the importance of static electron correlation beyond a (semi)local approximation. On the other hand, our results also show reality and robustness of the orbital concept, thereby making molecular orbitals accessible to experimental observations.

Introduction

Despite the fact that orbitals are not, strictly speaking, quantum mechanical observables,[1] a number of experimental techniques have evolved that enable orbital imaging.[2−7] While both the amplitude and the phase of the highest occupied molecular orbital (HOMO) of N2 could be recovered in three-dimensional space by high harmonics generated from intense femtosecond laser pulses, this tomographic reconstruction method seems to be applicable only for very simple molecules in the gas phase and only to the frontier occupied orbital.[3] On the other hand, scanning probe techniques applied to larger molecules adsorbed on surfaces have led to fascinating real space images of more complicated orbital structures with submolecular resolution.[4,8−1200] However, only frontier orbitals are accessible, and experiments require inserting an insulating decoupling layer between the metal substrate and the molecule and/or the use of functionalized tips.

It has been known for a long time that the angular dependence of the photoelectron emission from valence bands of molecular films contains rich information on the orbital structure of molecules which can be revealed by the angular resolved photoelectron spectroscopy (ARPES).[12,13] While standard angle-integrating photoemission experiments only reveal energy level positions of orbitals, the photoelectron angular distribution (PAD) is a fingerprint of the orbital structure, a momentum space image of the orbitals. More specifically, by approximating the final state of the emitted electron by a plane wave, the PAD can be interpreted as the Fourier transform of the initial state orbital,[5] or more precisely, its Dyson orbital.[14,15] Note that a Dyson orbital is the best approximation for the difference of the quantum mechanical states of a molecule and its cation, which is obtained by a simple direct one-electron ionization process.[1] This approach, termed photoemission tomography, has enabled the real-space reconstruction of molecular orbitals from ARPES data.[5−7,16] It also provides an orbital-by-orbital characterization of experimental spectra, not restricted to only the HOMO and LUMO.[14,17−19] Thereby, it yields detailed information on the energetic order and spatial structure of orbitals, which can be used as a most stringent test for ab initio electronic structure theory including density functional (DFT) calculations[20] as well as wave function-based approaches.

Results

Here, we focus on the PTCDA molecule which, as one of the most intensively studied organic semiconductors, serves as an experimentally and theoretically well-tested system.[21,22] Of particular interest for optoelectronic applications are its molecular π orbitals (MOs), which can be classified according to their symmetry and momentum distribution.[17] The frontier orbitals, i.e., the LUMO (A) and HOMO (B), as well the four next lower lying π-orbitals (C–F) are depicted in Figure , where the alphabetical labeling refers to the energetic order as obtained from DFT when employing a generalized gradient approximation (GGA)[23] for exchange-correlation effects.

Figure 1

Kohn–Sham orbitals at the GGA-level for PTCDA in the gas phase. The alphabetic labeling follows the orbital energy ordering of a GGA-DFT calculation, and the symmetry group of the orbitals are denoted.

Each MO shown in Figure exhibits characteristic PADs in accordance with its symmetry and nodal pattern.[6,14,16,17,24] In the so-called brick-wall (BW) phase of the PTCDA/Ag(110) monolayer, all molecules orient their long molecular axes along the substrate’s [001] direction, which makes the analysis of the data somewhat simpler. This has allowed us to unambiguously assign the observed emission at 1.9 eV below the Fermi level to the HOMO and a broader emission feature centered around 0.90 eV binding energy to the former LUMO occupied upon adsorption on silver due to charge donation from the metal.[17,24,25] Moreover, with our tomographic method that uses the energy and momentum dependence of ARPES data, we have been able to deconvolute a broad emission feature between 2.8 and 3.7 eV into individual contributions of the orbitals labeled C–F in Figure .[17] This data is reproduced for completeness in Figure a. It shows that, in contrast to the DFT prediction, a different ordering of the orbitals C–F must be deduced from experiment. For instance, orbital F turns out to have the lowest binding energy out of these four orbitals, while DFT predicts exactly the opposite.

Figure 2

Experimental projected density of states for orbitals C, D, E, and F resulting from an orbital tomography analysis of various PTCDA/Ag systems. Panel (a) corresponds to the brick-wall (BW) phase of PTCDA/Ag(110), (b) to the herringbone (HB) phase on Ag(110), (c) to the T-phase on Ag(100), and (d) to the herringbone structure of PTCDA/Ag(111). For each monolayer structure, the insets depict structural models as well as the theoretical maps of orbitals C–F used to deconvolute the spectra.

In our earlier work,[17] we have speculated that the reason for this discrepancy is rooted in molecule–substrate interactions. To check this assumption, we now scrutinize the orbital energies C–F on two other substrates with different molecule–substrate interaction strengths. Our choice fell on the (100) and (111) faces of Ag, because well-ordered monolayers are formed[26] and the adsorption heights of PTCDA have been well-characterized by X-ray standing waves.[27−29]

Judging from the measured adsorption heights of the carbon backbone ranging between 2.56 Å for Ag(110) over 2.81 Å for Ag(100) to 2.86 Å for Ag(111), different interaction strengths are indeed to be expected because also the local adsorption geometry varies significantly. Moreover, the alternative herringbone (HB) structure on Ag(110), panel b, which can be obtained by specific growth conditions, allows us to study the influence of intermolecular interactions.[30]

Our current ARPES experiments have been performed at the BESSY synchrotron using the toroidal electron energy analyzer[31] and focus on the orbitals C–F whose momentum space representations are depicted as the insets of Figure a. To evaluate our data for the other three monolayer structures b–d, the theoretical PAD maps that are used for deconvoluting the experimental data must take into account the multiple orientations of the molecules and possible mirror domains arising from the 4-fold and 3-fold symmetries of the Ag(100) and Ag(111) substrates, respectively. These maps are depicted in the insets of Figure b–d. The higher symmetry of the patterns for the orbitals C–F becomes evident. Nevertheless, they are still clearly distinguishable from each other, enabling a complete deconvolution into a projected density of states (PDOS).[17,18] The PDOS curves resulting from such a k-space deconvolution of experimental data are shown in Figure a–d for the PTCDA monolayers on Ag(110), Ag(100) and Ag(111), respectively. More data, including energy distribution cuts along several azimuths as well as the residues of the fits, are shown in Supporting Information Figure S1, which further emphasizes the reliability of our deconvolution procedure.

We have summarized the orbitals’ peak positions in the PDOS curves in Figure . Apart from an overall shift to higher binding energies when going from the most reactive Ag(110) to the least reactive Ag(111), which is explained by the different work functions of the substrates,[26] the most prominent finding is that the energetic order of the orbitals is always the same, namely F–D–C–E for ascending binding energy. This is surprising, because it is well-established that upon adsorption on metals large π-conjugated molecules usually undergo modifications of their electronic properties and/or geometry. Examples are geometric distortions, charge redistribution inside a molecule, and/or between the molecule and the metal, or lifting of degeneracies. The identical orbital ordering lets us deduce that it is neither molecule–substrate nor molecule–molecule interactions that are responsible for the deviation of the observed orbital ordering from the DFT prediction. For instance, the orbital ordering for the two phases on the Ag(110) substrate, i.e., the BW and HB phases, is identical, while the molecule–molecule interactions are quite different due to a distinct lateral arrangement of the molecules,[30] as can be seen from the insets in Figure . All these arguments lead us to the important conclusion that the experimentally observed energetic order reflects an intrinsic property of the PTCDA molecule.

Figure 3

Orbital energy positions of the PTCDA orbitals C, D, E, and F. Experimental values for PTCDA on the Ag(110), Ag(100), and Ag(111) surfaces are compared to theoretical predictions at various levels of sophistication (see text for details). As a reference level, the HOMO has been aligned to the experimental HOMO position of 1.9 eV for the BW phase on Ag(110).

Thus, we are left to explain the deviations between experiment and theory by questioning the predictions of DFT itself. First, we investigated whether the slight bend present in the adsorbed PTCDA[27−29] or the additional charge filling the LUMO upon adsorption may be responsible for the different orbital order (see block labeled DFT in Figure ). However, at the level of DFT, neither the geometric distortions (see Supporting Information for details) nor the additional electron results in a reversal of the computed orbital order. Calculating the energies of the system with half-occupied orbitals did not alter the order of the orbitals, either. Next, we have checked whether hybrid exchange-correlation functionals would lead to an improved orbital ordering. This has indeed been the case for the energetic positions of the anhydride σ-orbitals, which are more localized and therefore prone to self-interaction errors.[32] Here, we employ the range-separated hybrid functional HSE06[33] and an optimally tuned range-separated hybrid functional (OT-RSH).[22,34] The latter has been shown to be a method for obtaining the outer-valence quasiparticle excitation energies from a density functional- theory-based calculation, with an accuracy that is comparable to that of many-body perturbation theory within the GW approximation.[22,35−37] For both types of hybrid functionals, the energetic order of the orbitals remains essentially unaffected. Finally, we have also performed G0W0 calculations within many-body perturbation theory on top of our GGA result, which again led to the identical orbital order as the DFT starting point.

These results have motivated us to employ a computational method that is capable of incorporating static correlation effects, which are important when the ground state is well described only with more than one (nearly) degenerate Slater determinant. In such a situation, standard exchange-correlation functionals of DFT are known to underestimate the magnitude of the correlation energy.[38] To this end, we have performed wave function-based calculations using the complete active space self-consistent field (CASSCF) approach. In a CASSCF wave function, the occupied orbital space is divided into a set of inactive or closed-shell orbitals and a set of active orbitals. All inactive orbitals are doubly occupied in each Slater determinant. On the other hand, the active orbitals have varying occupations, and all possible Slater determinants are taken into account in the calculation. Due to the multiconfigurational nature of CASSCF wave functions, static correlations between orbitals are taken into account. For our CASSCF results, the orbitals are defined as natural orbitals that are computed from the first order density matrix γ(r,r′) of the many-electron wave function:The eigenvectors of γ(r,r′) are the natural orbitals ϕ(r), and the corresponding eigenvalues are the occupation numbers n. According to the extended Koopmans theorem,[39] orbital energies of natural orbitals are obtained from the eigenvalues of the Fock operator expressed in the natural orbital basis and can be rigorously interpreted as ionization potentials of the system. In our work, we approximate CASSCF orbital energies by the diagonal elements of the effective Fock matrix in the basis of the natural orbitals. This simplification is justified, because the computed orbital order turns out to be robust. Specifically, it is unaltered when ordering the orbitals (i) according to the occupation numbers n, or (ii) by using pseudocanonical orbitals, or (iii) when using state energy differences between the neutral molecule and the cation.

Results for PTCDA, both in the flat gas phase geometry as well as in the bent adsorption geometry (on the Ag(110) surface), for the neutral molecule as well as for the anion, and for different sizes of the active space are summarized in the bottom part of Figure . It is apparent that the order of orbitals C, D, E, and F is indeed altered compared to DFT. For instance, the CASSCF result for the neutral, flat PTCDA with 12 occupied and 12 unoccupied orbitals as active space (first row in the CASSCF block) moves orbital F to the low binding energy side and also pushes orbital C to higher binding energies. Compared to DFT, this leads to a much better agreement with the experimental observations. Enlarging the active space to 16 orbitals for both the occupied and unoccupied states shows that full convergence has not yet been achieved, since orbitals D and F as well as E and C change their places. It also suggests that static correlations among the orbitals C, D, E, and F may be significant.

We now analyze how the charge state of the molecule and its geometry affect the CASSCF results. Upon adsorption on the Ag surfaces, an overall charge transfer of 0.5–1.0 electron has been observed.[40] Our CASSCF calculations for the anion show a narrowing of the energy spread of the four orbitals, improving the agreement with experiment. The orbital order for the bent anion becomes close to the experimental one, importantly, orbital F is moved from the high to the low binding energy side of the spectrum. However, the comparison between the 13/12 (last but one row) and the 17/16 (last row) CAS also demonstrates that full convergence has not yet been reached, which, unfortunately, presents the limits of this computational approach. It should also be noted that for a one-to-one comparison between CASSCF results and experimental data, the influence of the substrate may have some effect. While the observed orbital ordering is independent of the crystal face and intermolecular arrangements, orbital energies may be affected by metallic screening, which is not captured in the CASSCF results.

An important question now is whether the strong static correlations evidenced from the sensitivity of our CASSCF results, in addition to influencing the orbital order, would alter the shape of the orbitals in a manner that would affect the experimental data deconvolution. This is, however, not the case, as illustrated in Figure . Here, we compare PADs of the orbitals C–F computed within GGA-DFT (top row) with the corresponding maps of the natural orbitals resulting from various CASSCF computations. Evidently, the momentum maps are quite robust, irrespective of whether they are derived from DFT or extracted as natural orbitals from the CASSCF calculations. This is because the overall symmetry and nodal patterns of orbitals are insensitive to the actual definition of the orbital (Kohn–Sham vs natural orbital), to the geometry of the molecule (flat vs bent), and the charging of the molecule (neutral vs anion). Clearly, subtle modifications in the PADs are discernible, however, these are small enough such that they do not affect the k-space deconvolution of the experimental data. We have cross-check our results for the BW phase on Ag(110) and found the same energy ordering when using CASSCF k-space maps for the photoemission tomography analysis as with GGA maps.

Figure 4

Theoretical photoelectron angular distributions (PADs) for the PTCDA orbitals C, D, E, and F ordered according to the binding energy from right to left. The top row depicts PADs of the Kohn–Sham orbitals using a GGA, while the remaining rows are based on natural orbitals of various CASSCF calculations.

Conclusions

In summary, we have demonstrated that photoemission tomography is able to provide an orbital-by-orbital characterization of organic molecules, thereby creating a compelling benchmark for ab initio electronic structure theory. For the case of the organic dye and semiconducting molecule PTCDA, the experimental data suggest an orbital ordering, which challenges density functional calculations, even if carried out with state-of-the-art exchange-correlation functionals. Despite inherent convergence problems, wave function-based complete active space calculations indicate the importance of static correlations in such molecules. Thus, we can understand why DFT, using a single-reference description, fails in predicting the correct orbital order in PTCDA, and presumably for similar types of molecules as well, no matter which state-of-the-art DFT method is used. This insight is based upon our CASSCF results demonstrating the importance of a multireference description. In more general terms, our findings substantiate the concept of orbitals as meaningful and experimentally observable objects and demonstrate the robustness of orbital patterns as observed in photoemission tomography. We find a good agreement between Dyson orbitals observed in photoemission tomography and natural orbitals computed from CASSCF, with respect to both energy and k-space distribution, while Kohn–Sham orbitals only agree in their k-space distribution but exhibit a different energy ordering. Our results have, however, also shown that the accurate treatment of correlations in such molecules is computationally challenging. Thus, for the future, it will be important to apply novel approaches such as density matrix renormalization group (DMRG) calculations or new DFT functionals developed for highly correlated systems to such systems.

Experimental and Computational Methods

PTCDA monolayer films were grown on Ag(100), Ag(110), and Ag(111) substrates as described in a previous publication.[26] Photoemission experiments were performed at BESSYII with 30 and 35 eV photons. An incidence angle of 40° with respect to the surface normal was used. The polarization direction was in the specular plane, which also contains the measured photoelectron trajectory. Photoelectrons were detected with a toroidal electron-energy analyzer,[31] which simultaneously records photoelectrons emitted with polar angles of −80° to +80° with respect to the surface normal. For our tomographic analysis of angle resolved photoemission spectroscopy, we collected ARPES data in the full half-space above the sample surface by stepwise rotating the sample azimuthally over 180°.

The DFT calculations were performed using the NWCHEM code,[41] utilizing a cc-pVTZ basis set. As an approximation for the exchange correlation potential, we employed a semi local functional,[23] the range-separated hybrid functional HSE06,[33] as well as an optimally tuned range-separated hybrid functional.[22,34] The CASSCF calculations were performed with the MOLPRO code.[42] Natural orbitals were obtained from the diagonalization of the first-order density matrix. The occupied orbitals and their energies were not altered by the change in basis set (6-31G*, cc-pVDZ, cc-pVTZ).