Angle resolved photoemission from organic semiconductors: orbital imaging beyond the molecular orbital interpretation

Fascinating pictures that can be interpreted as showing molecular orbitals have been obtained with various imaging techniques. Among these, angle resolved photoemission spectroscopy (ARPES) has emerged as a particularly powerful method. Orbital images have been used to underline the physical credibility of the molecular orbital concept. However, from the theory of the photoemission process it is evident that imaging experiments do not show molecular orbitals, but Dyson orbitals. The latter are not eigenstates of a single-particle Hamiltonian and thus do not fit into the usual simple interpretation of electronic structure in terms of molecular orbitals. In a combined theoretical and experimental study we thus check whether a Dyson-orbital and a molecular-orbital based interpretation of ARPES lead to differences that are relevant on the experimentally observable scale. We discuss a scheme that allows for approximately calculating Dyson orbitals with moderate computational effort. Electronic relaxation is taken into account explicitly. The comparison reveals that while molecular orbitals are frequently good approximations to Dyson orbitals, a detailed understanding of photoemission intensities may require one to go beyond the molecular orbital picture. In particular we clearly observe signatures of the Dyson-orbital character for an adsorbed semiconductor molecule in ARPES spectra when these are recorded over a larger momentum range than in earlier experiments.

The imaging of molecular orbitals has attracted great attention in recent years in rather different areas of physics [1][2][3][4][5][6][7][8][9][10]. This interest is fueled by at least two motivations. First, images that in some sense reflect quantum mechanical probability densities are fascinating from a quantum mechanical, many-body point of view, because they have been interpreted to visualize 'a single electron' in a many-electron system, i.e. an object whose observability may be questioned on fundamental grounds. Second, the molecular orbital picture, although approximate, has proved invaluable for understanding many-particle systems and is frequently and successfully used in particular in the field of organic semiconductors. Having tools to visualize, e.g. the highest-occupied and lowest-unoccupied molecular orbitals (HOMO and LUMO, respectively) that are frequently employed for understanding charge transfer and interface effects helps our understanding of this important class of materials [3,7]. Angle resolved photoemission spectroscopy (ARPES) has emerged as one of the most powerful techniques in the field [3,[6][7][8][9][10][11][12] as it allows for studying large molecules in setups of technological relevance, e.g. adsorbed and interacting with a substrate. The question of what actually is visualized in molecular orbital (density) imaging experiments has been discussed from several perspectives [8,13,14] and the power of the molecular orbital interpretation of ARPES experiments has been impressively demonstrated. In this work, however, we take one step further and demonstrate in a combined theoretical and experimental study that in certain cases the interpretation of ARPES experiments for molecular semiconductors may require one to go beyond the single molecular orbital point of view. To demonstrate this, we review in section 2 the theory of the photoemission process which reveals the Dyson orbitals as the fundamental quantities for describing photoemission. In section 3 we present an approximate way of calculating Dyson orbitals that combines computational feasibility with predictive power for the qualitative features of ARPES intensity patterns. A decisive aspect of our approach is that it explicitly takes into account electronic relaxation effects, i.e. the changes in electronic structure that occur when one electron leaves the system. The latter effect is typically neglected in the molecular orbital interpretation of photoemission data. We discuss the approximations that are made in our approach and show how it can be implemented as a straightforward extension of existing electronic structure theory methods. Finally, we discuss how our calculated Dyson orbitals compare to orbital imaging data for the N 2 molecule and the molecular semiconductors pentacene, 3,4,9,10-perylene-tetracarboxylic-dianhydride (PTCDA), and 1,4,5,8-naphthalene-tetracarboxylic-dianhydride (NTCDA). We show that imaging experiments resolve the differences between molecular orbitals and Dyson orbitals. This demonstrates that even for organic semiconductor molecules for which the molecular orbital picture is frequently employed, a detailed understanding of ARPES experiments can require one to go beyond the molecular orbital point of view.

Photoemission and the Dyson orbitals
To set the ground for explaining our approach, we first review the basic equations from which the Dyson orbital emerges as the fundamental concept for interpreting vertical (in the Franck-Condon sense) photoemission intensities [15][16][17][18]. In doing so we explicitly point out the steps in which approximations are made.
According to Fermiʼs golden rule the probability for a transition from an initial N-particle state Ψ N in to a final state Ψ N f that is triggered by a photon of energy ω  of an electromagnetic field that is described by a vector potential A is given by Here, P is the momentum operator, and E in and E f are initial and final state energy, respectively. This matrix element describes photoemission, i.e. the photoemission intensity is proportional to → W in f , if the final state is one in which an electron has been detached. We presume that the Nelectron molecule is initially in its ground state Ψ N in,0 with energy E in,0 . In writing the final state one commonly makes a first approximation by assuming that the correlation between the emitted electron and the remaining ones can be neglected. Then, the final state can be approximated as an antisymmetrized product of the wavefunction γ x ( ) k of the liberated electron, with momentum  k and with σ = x r ( , ) denoting spatial and spin degree of freedom, and the wavefunction Ψ − − x x ( ,..., ) 1 , which is the Ith eigenfunction of the − N ( 1)-electron Hamiltonian of the ionized system. Note that for vertical photoemission, which we consider here, the (N − 1)-electron Hamiltonian contains the same nuclear coordinates as the N-electron Hamiltonian because nuclear relaxation is occuring on a significantly longer timescale. Thus, the final state wavefunction is written as  , and the second term is the kinetic energy of the emitted electron. Thus . Making a second approximation by assuming that A is constant in space, one can write the matrix element as In the second step we split the sum into the direct term i = j and the remainder ≠ i j. Photoemission results from the direct term because in this expression the perturbation only acts on the liberated electron. The term with ≠ i j is responsible for other effects, e.g. Auger transitions [16], and will not be considered in the following-which is the third approximation. The transition matrix element can then approximately be written as This naturally leads to the definition of the Ith Dyson orbital d I . It is the decisive concept for the present work. With the Dyson orbital, the many particle matrix element can be written in the form of a single-particle matrix element By defining = x x N and using the relation For further evaluation of the matrix element we need to specify the final state of the emitted electron, γ k . With the assumption that this is a plane-wave, i.e. assuming that the emitted electron is completely free and not affected by the remaining moleculeʼs potential (which is the fourth approximation), one obtains a transparent and charmingly simple result: inserting γ ∝ kr exp (i ) k into the matrix element and Fermiʼs golden rule, respectively, and letting the momentum operator act on the plane wave, the probability of finding an emitted electron with momentum  k becomes proportional to the Fourier transform of the Dyson orbital The plane waveʼs wave-vector, i.e. the emitted electronʼs momentum, follows from the energy conservation law Equation (10) shows that, apart from the light-polarization dependent scaling factor A k · , 2 the photoemission intensity is determined by the Fourier transform of the Dyson orbital. This opens a fascinating interpretation of angular resolved photoemission experiments: the ARPES intensity maps a (hemi)spherical cut [6] through the momentum representation of the Dyson orbital, or more precisely, through its absolute value, with the radius of the cut being determined by the kinetic energy of the emitted electron via equation (11). As the conditions for the validity of the four approximations that were mentioned in the course of the derivation are usually well fulfilled, ARPES experiments can be interpreted as visualizing Dyson orbitals. However, from the perspective of the above reviewed derivation the question 'When and why can molecular orbitals be observed?', which is underlying imaging experiments, can be phrased more precisely in the form 'It is plausible that Dyson orbitals can be observed, but when and why are Dyson orbitals close to molecular orbitals?' In the following we take a closer look at this question and give some answers.

Approximate Dyson orbitals from DFT
The most straightforward relation between Dyson orbitals and molecular orbitals is to identify one with the other. Within this approach, which we call the 'molecular orbital interpretation' or 'single orbital picture' of photoemission, the ith Dyson orbital (i.e. the Dyson orbital corresponding to the ith ionization potential) is identified with the molecular orbital φ x ( ) i corresponding to the ith eigenvalue i i In Hartree-Fock (HF) theory, Koopmans' theorem and its complete neglect of correlation and electronic relaxation build a formal bridge, but the accuracy of HF orbitals and eigenvalues in practice is very limited [19][20][21]. It has been argued that the identification of DFT orbitals with Dyson orbitals can be much more accurate because of the inclusion of correlation and the formal similarity between the Dyson equation and the Kohn-Sham (KS) or generalized KS equation [15,22,23]. Whether the KS or the generalized KS approach are to be preferred in the context of photoemission can be viewed from different perspectives. On the one hand, if one is interested in eigenvalues that well map ionization potentials and electron affinities, the generalized KS approach has benefits [21,23]. On the other hand it has been stressed that the KS orbitals are more physical than the generalized KS orbitals [20], and exact KS eigenvalues have been shown to be excellent approximations to the ionization potentials of upper valence electrons [19]. We will show below that in our present work it is not the difference between KS and generalized KS theory (which we both use here) that is decisive, but the quality of the approximation used for the exchange-correlation (xc) potential. Systems in which all orbitals are delocalized on similar length scales [24,25] can be well described by semi-local functionals [6,9]. For systems with a more complicated electronic structure, semi-local functionals and global hybrid functionals can lead to a distorted spectrum, whereas a self-interaction correction carried out in the KS approach [26] has been shown to yield KS orbitals that well reproduce ARPES experiments also in such situations [8]. However, despite these successes it must be kept in mind that even when the exact xc potential could be used, the identification between KS orbitals and Dyson orbitals remains approximate and cannot include exactly all the physical effects that play a role in ionization [23,27,28]. An inherent limitation of the single orbital picture of photoemission is the assumption that the only change in the electronic structure of a molecule upon photoionization is the depopulation of just one specific orbital. The state of the remaining, positively charged ion is assumed to be equal or at least very similar to the initial state except for the emptied orbital. This neglect of electronic relaxation, i.e. neglecting that the remaining electrons restructure once one electron has been removed, is to be questioned for typical molecular systems. In the following we explicitly go beyond the 'one-orbital emission and no relaxation' picture and study the first photoionization event, i.e. the transition from the ground state of the initial molecule to the ground state of the singly ionized system, from a many-particle perspective. We focus on the first transition because the first ionization potential is often the most relevant one in practical applications. In the single particle picture it corresponds to the view of 'an electron is emitted from the HOMO'. A formal reason for why it is natural to focus on the first transition is that the HOMO plays a special role in DFT, because its eigenvalue can rigorously be related to the ionization energy [29].
The starting point for our considerations is the above reviewed definition of the Dyson orbital as the overlap of the initial and final molecular ground-state wavefunctions, equations (6), (8). Using this definition for calculating the Dyson orbital of molecular semiconductors in practice requires one to face two challenges.
The first is that one needs to calculate the N-and ( − N 1)-electron wavefunctions, respectively. For molecules with many tens to hundreds of electrons, calculating an explicitly correlated wavefunction is not just a very demanding task-it is prohibitively expensive. Therefore, we here focus on the alternative approach of approximating the true wavefunction Ψ by a single Slater determinant Φ. This is possibly a very serious approximation. However, for the molecules that we are interested in it is typically justified because they are dominated not by static but by dynamical correlation and their wavefunction is predominantly of single-reference character. In other words, when written as a configuration interaction expansion, the fully correlated many-body wavefunction has one leading, dominating determinant with small corrections from many other determinants. The paradigm one-determinant approach of HF theory may still be inaccurate in such a situation. However, a KS Slater determinant, i.e. the Slater determinant built from the occupied orbitals that are obtained after converging a selfconsistent KS calculation, is often an acceptable approximation to the true wavefunction. While there is no guarantee, even if the exact xc functional were known, that the KS determinant approximates the true wavefunction well, experience shows that for systems that are not strongly correlated this is very often the case [30]. This finding can be rationalized by the observation that the KS determinant is density-optimal instead of energy optimal [31,32].
The second challenge is the calculation of the overlap integral between the N-and the − N ( 1)-electron wavefunction according to equation (6). We address it by expanding the N electron Slater determinant in minors [16] i.e. decompose it into a sum of Slater determinants build up by − N 1 orbitals. Each Slater determinant is multiplied by the single orbital φ i in, that depends on the remaining coordinate. Thereby, the Dyson orbital can be rewritten as This shows explicitly that the Dyson orbital does not correspond to a single-particle orbital, but is a coherent superposition of all N occupied orbitals that build the initial Slater determinant, weighted by w i . The weighting factors are given by the overlaps of ( − N 1)-particle determinants. The calculation of the w i is mostly a technical issue. We delineate the procedure in appendix A.
Central for our work is the choice of orbitals that we use to build up the initial and final state Slater determinants. The crudest approximation would be to use a simple Koopmans like picture where relaxation effects are completely neglected. The only change of the moleculeʼs electronic system upon ionization would then be that one initially occupied orbital is unoccupied after the ionization process. Hence, one would use the same set of orbitals to build the initial and final state Slater determinants. The expansion in minors of the initial state according to equation (13) would lead to one contribution where the single orbital is equal to the depopulated orbital and the − N 1 determinant from the corresponding expansion term would be equal to the unrelaxed ionized state. The overlap between the two latter determinants would therefore be one, whereas all other overlaps respectively weighting factors would be zero due to orthogonality. As a consequence the Dyson orbital would be equal to the emptied orbital and one would return to the single orbital picture.
Our aim is, however, to treat electronic relaxation explicitly in order to check how large its impact on vertical photoemission and orbital interpretation is, respectively. Therefore, we perform two separate self-consistent calculations, one for the ionized and one for the initial system, while keeping the moleculeʼs nuclear geometry fixed. In this way we find in addition to the initial orbitals a second set of orbitals-non-orthogonal to the first-representing the ground state electron density that minimizes the total energy of the positively charged molecule. Having two nonorthogonal orbital sets φ n in, and φ m f, in general leads to non-vanishing overlaps in equation (13), i.e. weighting factors w i that are not restricted to zero or one. Hence, relaxation effects can lead to more than one molecular orbital contributing to the Dyson orbital. How much the structure of the Dyson orbital differs from the one of a single molecular orbital therefore depends on how strongly electronic relaxation influences the ionization process, and on how large the differences between the different molecular orbitals are. A detailed description of the calculation of the overlaps with two different sets of orbitals is given in appendix A.

Photoemission from the N 2 molecule
As a first test for the relaxed Slater determinant approach (RSA) that we introduced in the previous section, we chose the nitrogen molecule because here the spatial structure of the Dyson orbital is known from higher harmonic generation (HHG) experiments [1]. As a first step we examine the influence that the choice of the DFT xc-functional has on the structure of the RSA Dyson orbital. Therefore, we have performed calculations with the widely used semilocal Perdew-Burke-Ernzerhof generalized-gradient-approximation (PBE-GGA) [33] and the Becke three-parameter Lee-Yang-Parr (B3LYP) [34][35][36] hybrid as a representative for a non-local generalized KS approach. In addition we used a self-interaction-correction (SIC) KS approach based on the generalized optimized effective potential (KS-SIC) [26,37]. In all cases the resulting RSA Dyson orbitals have an almost indistinguishable spatial structure, displayed in figure 1(a), and this is in very good agreement with the orbital that has been experimentally reconstructed by Itatani et al [1]. The Dyson orbital in each case is mainly dominated by the corresponding σ-HOMO, HOMO . All other orbital contributions are at least 30 times smaller.
It is quite instructive to investigate the situation with HF theory as the single-particle starting point. In contrast to the DFT results described above, HF theory leads to a degenerate highest eigenvalue which is associated with two π-shaped orbitals, see figure 1(c). This implies that the unrelaxed approximation, i.e. the molecular orbital interpretation in the HF approximation, leads to a Dyson orbital that is a superposition of these orbitals, in disagreement with the experiment. However, when one builds the final state from relaxed HF orbitals, i.e. uses a HF-based RSA, one finds a Dyson orbital that solely corresponds to the HF σ-HOMO-1, HOMO 1 (for a detailed explanation see appendix A). It is shown in figure 1(b) and, again, reflects all features of the orbital reconstructed by HHG. Therefore, a HF-based calculation leads to dramatically different results depending on whether electronic relaxation effects are taken into account or not. This failure of the unrelaxed HF approximation of course is not a major surprise since correlation effects are known to be important in multiple bonds [38], and HF theory for N 2 leads to the wrong energetic ordering of the HF HOMO and HOMO-1.
In summary, this first, small-molecule example demonstrates two things. First, treating relaxation effects can be decisive for a correct description of photoionization in terms of orbitals. Second, it points out a strength of the RSA. Equation (14) shows that the RSA Dyson orbital only depends on the overlaps w i . Therefore, as long as two different single particle calculations, e.g. a HF and a DFT one, lead to the same set of occupied orbitals, the corresponding DFT and HF RSA Dyson orbital will look the same, i.e. are independent of the internal energetic orbital ordering. The question of having the correct orbital ordering that can complicate ARPES interpretation based on molecular orbitals [6,8,9] thus becomes immaterial or at least mitigated.
Beginning with the simplest case, pentacene, we know that the molecular orbital ordering hardly depends on the used xc-functional [21], i.e. all functionals yield very similar HOMOs.
(The corresponding momentum maps are compared to each other in detail in [8].) It is a reassuring finding that it also does not matter which approach we use as a basis for the RSA calculation. For the above mentioned functionals and HF, the RSA Dyson orbital is equivalent to the HOMO

RSA HOMO
This simplicity in the analysis is due to pentaceneʼs simple electronic structure: with all frontier orbitals being π-orbitals delocalized on the same length scale, the ordering is not sensitive to the used xc-functional, and ionization hardly changes the spatial structure of the orbitals. The impressive resemblance of the pentacene HOMO momentum space map to the corresponding measured ARPES intensity has been discussed in detail in [6].
For the more complex molecule PTCDA the situation is quite similar. Again all approaches yield the same orbital, i.e. the obtained Dyson orbital is very close to the HOMO, and relation (15) holds in this case too. In appendix B, in particular figure A.1, we show in detail that the RSA Dyson orbital as well as the HOMO momentum space maps are in good agreement with experimental ARPES data.
However, the situation gets more involved for NTCDA, which is the main focus of the present investigation. In contrast to pentacene, the outer valence orbitals of NTCDA are alternately of πor σ-character. Their ordering is very much influenced by which xc-functional is used. An important aspect in explaining this finding is the self-interaction error whose magnitude strongly depends on the orbital character [24]. We see this effect in particular for the HOMO. Whereas B3LYP, HF and KS-SIC predict a π-HOMO, PBE shows a σ-HOMO. The two orbital types π and σ show distinct momentum space maps as displayed in figure 2. The π-orbital has its four main features at ≈ ± Å − k 1. . The σ-orbital on the other hand primarily has high intensities at higher k | |-values. This is a consequence of the nodal structure of the probability density in real space: the σ-orbitalʼs nodes are more closely spaced than those of the π-orbital.
Again we want to contrast the molecular orbital interpretation with the RSA, where PBE is particularly interesting, as it is the odd one out. With the HF N 2 case in mind our hope is that the RSA approach can compensate for shortcomings of the underlying single-particle theory. Indeed, the RSA based on PBE yields a Dyson orbital that differs from the PBE HOMO. It is a superposition of the PBE HOMO and the HOMO-1, in detail   photoemission process. However, for NTCDA the 'unification' is not as complete as for, e.g. the above discussed case of N 2 . This becomes visible when one builds the RSA Dyson orbital starting from a KS-SIC calculation. In this case the Dyson orbital is mostly determined by the

KS SIC HOMO
Also in the KS-SIC Dyson orbital there is more than one contribution, but the weighting factor for the HOMO is = w 0.98 HOMO , whereas the next contribution (HOMO-16) comes with a much smaller weight of = − w 0.06 HOMO 16 , and the intensity scales with the square of the weighting factors. Summarizing the results obtained for NTCDA, we can therefore say that for the momentum range ≈ Å , i.e. compared to previous experiments we extended the observation range to larger momenta. For PTCDA and NTCDA the influence of the Ag substrate on the orbital which is the HOMO of the isolated molecules is small [7,42]. The experimental pattern observed for the peak that traditionally has been associated with 'the HOMO' is depicted in figure 3(c). Comparing it to figure 2(a) shows without doubt that the measured intensity cannot be explained just by the PBE HOMO. The latter lacks the four peaks at ≈ ± Å − k 1. . From an experimental point of view there is no doubt that all of these features belong to the same ionization process, as their intensities occur uniformly in the same energy interval and well separated from other peaks.
Comparing to the RSA Dyson orbital in figure 3(a) shows that the PBE based RSA Dyson orbital clearly shows such features. In order to reveal them more clearly, we repeated the evaluation of the theoretical spectra and the experiment with horizontal light polarization (i.e. spolarization). In this way high k | |-values are more accentuated. In our setup the NTCDA monolayer is aligned such that the polarization is parallel to the long molecular axis (in our notation defined as the x-axis). On the theoretical side this brings the polarization factor of equation (10) into play. With the vector potential being parallel to the x-axis, the relation of the Dyson orbitals to the measured spectra is thus augmented by the smooth, here k x -dependent A k · x x 2 factor, and becomes As a consequence features at large k | |-values that are not well visible in figures 3(a) and (b) can now be compared much better to the experimental data.
In figure 4 we display the ARPES intensity as measured with horizontal light polarization. The experimental data has been overlaid with the contours of the PBE Dyson orbital on the left hand side and with the ones from the KS-SIC Dyson orbital on the right hand side. Both Dyson orbitals show features at large k | |. This in itself is a reassuring observation. There are, however, differences between the two calculations. First, in the PBE case the outer features emerge from the σ-orbital, whereas in KS-SIC they are purely attributed to the π-orbital. Second, there are differences in the relative intensities. For the PBE contour lines in figure 4(a) we see the trend that the intensity of the peaks increases with their position on the k x -axis so that the features at high k x are more apparent. This is in better agreement with the experiment than the relative intensities from the KS-SIC Dyson orbital. Third, however, there are also differences in the positions of the outer peaks, and in terms of these the RSA KS-SIC Dyson orbital is in better agreement with the experiment than the PBE one. For putting these findings into perspective one needs to be aware of additional issues. One of them is that higher k | |-values correspond to emission angles closer to the substrate. Thus, the interaction between the outgoing photo electron and the remaining system could be stronger and lead to larger inaccuracies within the plane wave approach. Another one is that for the high ∥ k -values the slight misalignment of the crystals surface normal with respect to the optical axis of the photoemission electron microscope leads to different k x -and k y -ranges. Additionally, this leads to distortions of the kimage that can lead to a mismatch of the k x -position of the emission peaks. We are therefore not able to make an unambiguous statement about which approach yields the better result. However, the overall picture that emerges from this study is consistent. It has been argued for a long time [43] that SIC eigenvalues can approximate relaxed ionization energies well. It thus appears natural that the SIC molecular orbitals can reflect effects which in another approach such as PBE must be included via explicit consideration of relaxation. The RSA Dyson orbital calculation for NTCDA amends the peaks at high k | |-values and as well captures the intensities at lower k | | that are completely missing in the molecular orbital interpretation of the PBE HOMO. Thus, the experimental data can be explained over the full range of measured electron momenta.

Conclusions
Motivated by ARPES experiments that can be interpreted as revealing orbital densities, we reviewed the theory of photoemission in which the Dyson orbitals, defined as the overlaps between N-and ( − N 1)-particle states, emerge as the relevant quantities. When Dyson orbitals are a priori approximated by single molecular orbitals the quality of the so predicted ARPES intensities can depend sensitively on which electronic structure method is used to calculate the molecular orbitals. We demonstrated that by approximating the initial and final molecular states by relaxed N-and ( − N 1)-electron Slater determinants, one achieves a description of the first photoemission process that is considerably less sensitive to the underlying single-particle picture and its orbital ordering. More importantly yet, this procedure allows for explicitly taking into account relaxation effects. ARPES experiments were done for the model organic semiconductors PTCDA and NTCDA covering a larger range of electron momenta than covered by previous experiments, and using variable light polarization. For PTCDA the RSA Dyson orbital approach leads to the same result as the molecular orbital approach and both are in agreement with the experimental data. For NTCDA the experiment reveals previously unobserved intensities at large electron momenta that cannot be explained based on molecular orbitals as obtained from, e.g. a usual GGA DFT calculation. The RSA Dyson orbitals show the additional intensities. Thus, ARPES experiments have reached an accuracy at which a simple molecular orbital interpretation may not be sufficient for explaining all of the observed features. The RSA Dyson orbitals discussed here provide a non-empirical and computationally relatively efficient approach for going beyond the molecular orbital interpretation of photoemission. Seen from a different angle, one can also argue that the present work gives further credibility to the concept of orbital imaging-as long as one keeps in mind that it is Dyson orbitals that are imaged.
The overlap of two Slater determinants with the same number of orbitals can be expressed as a determinant of overlap matrix elements between all single-particle orbitals This leads to the following expression for the weighting factor Looking again at the expansion for the initial ground state 3), respectively. If both orbital sets do not differ much, then the single particle overlaps will be close to either one or zero. Consequently, only the matrix element of the one reduced Slater determinant which contains the orbitals that are similar to the orbitals in the ionized determinant will be close to one. When there is only little relaxation the Dyson orbital is dominated by the depopulated orbital. Strong relaxation leads to two clearly different orbital sets and several non-vanishing weighting factors.
With this we can also elucidate why the HF N 2 Dyson orbital is determined by the σ-HOMO-1, in spite of the existence of the degenerate π-HOMOs: in the cation N 2 + the ordering of orbitals is inverted such that the σ-orbital lies about 14 eV above the two π-orbitals. According to equation (A.1) the weighting factors are determined by the overlap of the reduced Slater with the Slater determinant of the ionized system, each containing − N 1 orbitals. As the σ-orbital of N 2 + is not populated and also is spatially similar to the N 2 HOMO-1, the only combination of two determinants that can be close to one involves neither of the two σ-orbitals. This can happen only for the overlap corresponding to the weighting factor of the σ-HOMO-1 of N 2 . In a similar manner, the Dyson orbitals' structure of NTCDA and PTCDA can be explained. All approaches for PTCDA + yield a depopulated orbital that has a spatial structure that is close to the PTCDA HOMO. Therefore, the overlap argument strikes again and forces all weighting factors to be small, except for the one that corresponds to the HOMO. In the case of PBE NTCDA + the depopulated orbital is rather a superposition of the π-HOMO-1 and σ-HOMO of NTCDA, which is why the Dyson orbial also emerges as a superposition of these two orbitals.
Also note that the norm of the Dyson orbital is not one, but ranges between zero and one. It can be interpreted as a measure of the ionization probability.
Appendix B. Photoemission from the PTCDA/NTCDA molecule For the sake of completeness, we here discuss the PTCDA results in greater detail. As mentioned in section 5, all approaches lead to a similar intensity pattern for the first ionization of PTCDA, independent of the inclusion of relaxation and the used functional. Therefore, the Dyson orbital is in each case given by the corresponding π-HOMO. Contrary to the NTCDA case there are no additional orbital contributions. In . The horizontally polarized case, which should reveal any additional σ-orbital contributions if there are any, is well reproduced by the calculation as far as the peak positions are concerned. The trend is similar to the KS-SIC NTCDA result, i.e. the outer features are higher in intensity in the experiment than in the calculation. The magnitude of the differences, however, is smaller for PTCDA than for NTCDA. For PTCDA we thus come to the conclusion that the single molecular orbital approach is well justified, and relaxation effects are smaller than for NTCDA. Expressed in numbers, the difference between the relaxed and unrelaxed first ionization potential calulated with PBE for NTCDA is about 0.93 eV, whereas for PTCDA we find only a difference of 0.16 eV. This appears plausible, as a larger system with more electrons should be less affected by removing one charge.

Appendix C. Experimental and computational details
All experiments were performed under ultrahigh vacuum conditions at a base pressure of 2 × 10 −10 mbar. The NTCDA and PTCDA monolayer films were prepared in an attached preparation chamber by organic molecular beam deposition. Prior to deposition the Ag(110) substrate was cleaned by several sputter and annealing cycles. The cleanliness and surface quality was checked by low energy electron diffraction (LEED) and the absence of characteristic carbon 1 s signals in core level photoelectron spectra. The coverage and quality of the organic monolayers was monitored by their characteristic LEED patterns [39,[44][45][46] and the absence of any second layer signal in valence photoelectron spectra [40,47,48]. The ∥ k -dependent photoemission intensity distributions were recorded at the NanoESCA beamline at Elettra [49] using a FOCUS GmbH/Omicron NanoESCA photoelectron microscope [50,51]. With this setup the complete photoemission hemisphere can be measured at once, if the Fourier plane of the microscope is imaged with an activated transfer lens. A detailed description of the measurement and evaluation process is given in [40]. In contrast to previous investigations [6][7][8][9][10]40] we used different light polarizations (linear horizontal and vertical, and circular) to enhance and suppress photoemission intensities. To maximize the available ∥ k -range at manageable beam damage we used photon energies of around 50 eV and carefully checked degradation effects for each measurement.
The KS-SIC and PBE calculations for NTCDA, PTCDA, pentacene and N 2 were performed with the Bayreuth version [28] of the PARSEC real-space code, [52] where we employed a grid-spacing of 0.3 bohr (for N 2 0.2 bohr) and Troullier-Martins norm-conserving pseudopotentials [53].
The KS-SIC method has been explained in detail in [26,37], so we here just summarize very briefly its main features. Based on the self-interaction-correction functional of Perdew and Zunger [43] we construct-in contrast to traditional SIC approaches-a spatially local, multiplicative xc potential via the optimized-effective potential (OEP) equation [54]. This ensures that we stay within the realm of KS DFT with all its benefits. The direct construction of the OEP, which is possible but tedious for the SIC energy functional [55], was avoided by resorting to the generalized OEP equation and the corresponding generalized KLI approximation [26] with a complex-valued energy minimizing unitary orbital transformation [37].
For all HF and B3LYP calculations we used the QChem program package [56] and cc-PVTZ basis functions [57]. We explicitly checked and verified that the use of the pseudopotentials does not lead to any differences with respect to the all electron calculations in the PBE RSA Dyson orbitals.