Can photoemission tomography be useful for small, strongly-interacting adsorbate systems?

Molecular orbital tomography, also termed photoemission tomography, which considers the final state as a simple plane wave, has been very successful in describing the photoemisson distribution of large adsorbates on noble metal surfaces. Here, following a suggestion by Bradshaw and Woodruff (2015 New J. Phys. 17 013033), we consider a small and strongly-interacting system, benzene adsorbed on palladium (110), to consider the extent of the problems that can arise with the final state simplification. Our angle-resolved photoemission experiments, supported by density functional theory calculations, substantiate and refine the previously determined adsorption geometry and reveal an energetic splitting of the frontier π-orbital due to a symmetry breaking which has remained unnoticed before. We find that, despite the small size of benzene and the comparably strong interaction with palladium, the overall appearance of the photoemission angular distributions can basically be understood within a plane wave final state approximation and yields a deeper understanding of the electronic structure of the interface. There are, however, noticeable deviations between measured and simulated angular patterns which we ascribe to molecule-substrate interactions and effects beyond a plane-wave final state description.


Introduction
In recent years, photoemission tomography (PT) has evolved as a powerful technique to study the geometric and electronic structure of oriented layers of organic molecules [1,2]. As a combined experimental/theoretical approach, PT uses angle-resolved ultraviolet photoemission spectroscopy (ARUPS) over a wide angular range and seeks an interpretation of the photoelectron angular distribution, so-called momentum maps, in terms of the molecular orbital structure of the initial state. PT has found many interesting applications including the assignment of molecular orbital densities in momentum and real space [3][4][5][6][7][8], the deconvolution of spectra into individual orbital contributions beyond the limits of energy resolution [9][10][11][12], or the extraction of detailed geometric information [13][14][15][16][17][18]. These successes have been on relatively large molecules adsorbed on noble metal surfaces and the usefulness of PT to smaller molecules and particularly on more reactive surfaces has not been investigated.
Indeed, approximating the photoelectron's final state by a plane wave have been viewed very critically by many people in the community [19][20][21][22][23][24][25]. Recently, the plane-wave approximation lying at the very foundation of PT has been in particular questioned in an article in the New Journal of Physics by Bradshaw and Woodruff entitled 'Molecular orbital tomography for adsorbed molecules: is a correct description of the final state really unimportant?' [26]. In their paper, the authors highlight situations in which PT is likely to fail in an adequate Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. explanation of photoemission distribution, and also propose experiments to test the technique's limits of applicability. One such test involves studying small molecules such as benzene, which differs in two respects from the larger molecular systems so far investigated with PT. Firstly, larger molecules often have a quasiperiodicity in their orbital structure, which smaller molecules, such as benzene, do not have. Secondly, benzene's electron density tends to be commensurate with the substrate, unlike that of its larger oligomeric relatives. The first will cause the momentum maps, i.e. k-space distribution of the orbitals, to be less distinct, and the second could lead to scattering becoming dominant in the angular distribution.
Angle-resolved photoemission studies of benzene adsorbed on various metallic surfaces have already been performed by many groups in the past, thereby revealing the modification of benzene's electronic structure upon adsorption. Moreover, by applying symmetry arguments [27], the orientation of benzene on Pd(100) [28], Pd(110) [29], Ni(111) [30] and Ni(110) [31,32] has been inferred. More detailed information of the adsorption structure including the adsorption sites and heights has been obtained from photoelectron diffraction experiments for benzene on Ni(111) [33] and Ni(110) [34], and early scanning tunneling microscopy (STM) investigations have revealed adsorption sites of benzene on Ni(110) and Cu(110) [35]. With the advent of modern computer technology, the investigation of benzene on Al(111) [36] and on the Ni(111), Ni(100) and Ni(110) surfaces has also been supported by density functional theory (DFT) calculations [37]. More recently, the role of van der Waals interactions for the adsorption of benzene on metals came into the focus of theorists [38][39][40][41], and benzene has served as a computationally tractable model system for organic/metal interfaces allowing for the application of more accurate and highly sophisticated electronic structure methods [42][43][44][45].
In this work, we thus return to the classical benzene/Pd(110) system and present ARUPS experiments supported by density functional calculations. By making use of a toroidal electron spectrometer [46], we measure photoelectron angular distributions of molecular emissions over the entire polar and azimuthal angular ranges which allow us to go beyond simple symmetry arguments of the past and to test the theoretical predictions based on a plane wave final state. Our results are not only important for the further advancement of the photoemission tomography approach, but also serve as benchmarks for ab initio electronic structure methods in terms of adsorption site and electronic level alignment which we determine in our study.

Sample preparation and experimental setup
The experiment was conducted at room temperature under ultra-high vacuum conditions. The Pd(110) crystal was cleaned in a conventional way using several cycles of sputtering and annealing. It was then exposed to a background dose of 10 L of benzene at room temperature, which resulted in the c(4×2) low-energy electron diffraction structure of the benzene-saturated surface as reported earlier [29]. The photoemission experiments were conducted using the toroidal electron energy analyser [46] with an energy resolution of 200meV at the insertion device beamline of the Metrology Light Source (Physikalisch-Technische Bundesanstalt, Berlin, Germany). A p-polarized, 31 eV photon beam at 40°to the surface normal was used to excite photoelectrons, which were collected over a polar angular range of±85°in the plane of incidence. To collect the data in the full hemisphere of emission, the sample surface was rotated azimuthally over 180°in steps of 1°. The constant binding energy momentum maps were then obtained by slicing the three-dimensional datacubes of the photocurrent at the desired binding energies.

Computational details
All electronic structure calculations were performed using the Vienna ab initio simulation package [47][48][49] utilizing a repeated slab approach for simulating the surface. The metallic substrate has been modelled by seven palladium layers with a vacuum layer of about 2.7nm. To avoid spurious electric fields, a dipole layer [50] was inserted in the vacuum region. The interaction between electrons and nuclei is modelled using the projectoraugmented wave method [51] (PAW) with an energy cutoff of 400 eV. Partial occupancies close to the Fermi level are evaluated using the first-order Methfessel-Paxton [52] method with a smearing width of 0.2 eV. Geometry optimizations have been performed by using a damped second order equation of motion for updating the ionic degrees of freedom. The convergence criterion for the forces has been set to 0.1 eV nm −1 and the convergence criterion for the electronic energy has been set to 10 −8 eV. The exchange-correlation effects are described using the generalized gradient approximation [53] (GGA), and to account for van-der-Waals interactions, which are ill-described in standard GGA functionals, we employed the vdW-surf method [54], which includes screening of the substrate electrons in the determination of the vdW parameters, going effectively beyond the pairwise description, which is necessary for a reliable description of geometries and energies of inorganic-organic systems. For an efficient Brillouin-zone integration, the Monkhorst-Pack [55] scheme has been used for the k-point grid generation. For the benzene gas-phase, a 6×6×6 grid has been empoyed while for the adsorbed, freestanding and surface geometries a 12×12×1 grid has been used. After relaxing the structure, we have evaluated the density of states (DOS) by employing the Heyd-Scuseria-Ernzerhof (HSE) hybrid functional [56,57]. The momentum maps were simulated by calculating the ARUPS intensity I(k x , k y ) with the algorithm published in [3].

Adsorption geometry calculations
Based on the c(4×2) overlayer structure, we have determined the adsorption geometry of benzene on Pd(110). To this end we have considered in total 8 starting points for the geometry optimizations. Relating to the center of the benzene ring, we have considered on-top, short-bridge, long-bridge and hollow sites of the (110) surface. And for each adsorption site, we have tested two orientations of the benzene ring, namely one where the long diagonal of the hexagon connecting two opposite vertices (dashed line in figure 1) points along the close-packed Pd-rows, i.e. the [ ] 110 -direction, and one where this diagonal is directed across the Pd-rows. For the hollow site, these two orientations are depicted in panels (a) and (b) of figure 1, respectively.
In agreement with a recent DFT study of benzene adsorption on Pd(110) [39], for both orientations of the molecule, we have found the hollow site to be preferred followed by the long-bridge, short-bridge and on-top sites (see table 1 for a compilation of the adsorption energies). When comparing the two orientations of the molecule, we find the along configuration of figure 1(a) to be energetically favoured by about 0.1 eV based on the vdW-surf functional [54]. It is important to note that the low symmetry C s configuration shown in panel (a), in which the benzene ring is laterally shifted by about 50pm in [001] direction and tilted by about 8°, is energetically lower by 0.4 eV compared to the higher symmetry C 2v configuration with benzene centered above the hollow site and with the ring co-planar to the surface.
The benzene ring approaches quite close to the surface, indicating a quite strong interaction. We find carbon atoms at heights between 184and 196pm for the across orientation, and similar heights between 175and 204pm, for the along configuration with a somewhat larger spread owing to the tilt of the molecule. We have also analyzed how the internal geometry of the ring is modified upon adsorption. As indicated in figure 1, we observe an overall increase in the C-C bond length from the gas phase value of 140pm and distortions towards a more quinoidal structure as expected due to the reduction of symmetry from D 6h to the C 2v (across) or C s (along) point groups upon adsorption. Note that for the along orientation, the C-C bond length alternations (143-146 pm) are somewhat more pronounced than for the across configuration. As can be seen from the insets  Table 1. Calculated adsorption energies in eV of benzene/Pd(110) for four adsorption sites and two orientations of the benzene ring. Adsorption energies have been computed by employing the vdW-surf [54] scheme to account for van der Waals interactions. Numbers in parenthesis are values from the underlying GGA functional [53].

On top
Short bridge Long bridge Hollow in panels (a) and (b) of figure 1, we also observe an upward bending of the C-H bond as a result of chemical bond of C with Pd and the concomitant switching from sp 2 towards sp 3 hybridization. While for the across configuration, similar out-of-plane angles for all six C-H bonds are observed, the two C-H bonds pointing along the Pd-rows in the along orientation remain nearly parallel to the surface.

Density of states
For gas phase benzene, each of the six carbon atoms contributes one p z electron to the π electron molecular orbitals, which are delocalized over the aromatic ring. There are three occupied π orbitals which we have denoted as π 1 , π 2 and π 3 . According to the D 6h symmetry, π 1 and π 2 are degenerate and belong to the irreducible representation e 1g . They give rise to the highest occupied molecular orbital (HOMO) of benzene in the gas phase.
The HSE-value of the gas phase ionization potential is 6.9 eV. This can be seen from panel (a) of figure 2, where we have plotted the respective DOS. Note that the orbital structures of π 1 and π 2 are shown in figure 1. Clearly, the orientation of the nodal plane, but also the shape of the orbital density (triangular for π 1 ) distinguishes π 1 from π 2 . The energetically deepest π-orbital, π 3 with the a 2u symmetry label, is at about 10 eV below the vacuum level, while the topmost σ orbital, the doubly degenerate e 2g , is located in between π 1,2 and π 3 at about −9.3 eV. Next, we analyze how the geometry distortions and the intermolecular interactions in the( ) c 4 2 overlayer structure alter the electronic structure of the benzene ring upon adsorption. To this end, we take the relaxed adsorption geometry of figure 1(a) (hollow-site, along) and consider its freestanding layer with the adsorbed system's geometry frozen. Its DOS is depicted in figure 2(b). In addition to the total DOS (black line), we also evaluate the DOS projected onto the two topmost π-orbitals shown as red and turquoise lines, respectively. Clearly, we observe a splitting of roughly 0.1 eV which can be traced back to the c(4×2) structure and bond length and bond angle changes induced by the adsorption. Note that here π 1 ends up with a higher binding energy as compared to π 2 .
The impact of molecule/substrate interactions can be appraised from panel (c) of figure 2, where we plot the total DOS of the benzene/Pd(110) as well as its DOS projected onto π 1 and π 2 , respectively, which has been multiplied by a factor of 10 for better visibility. Overall the DOS is dominated by the Pd d-states in the energy range from −9.2 eV up to the Fermi level at −3.7 eV. The two π states of interest, π 1 and π 2 , have peaks located at −8.85 eV and −9.05 eV, respectively, which are broadened through strong molecule-substrate hybridizations, and which are significantly stabilized on adsorption with respect to the σ orbitals in agreement with experiment. Most interestingly, their energetic order has reversed compared to the free-standing layer. So the splitting has not only just doubled its value to 0.2 eV, but we observe a 0.3 eV stronger stabilization of π 2 compared to π 1 , leaving π 1 as the lowest binding energy molecular state on the surface. Let us now compare these theoretical findings to experimental ARUPS data as depicted in panel (d) of figure 2. In particular, we have measured the photoemission intensity for polar angles from θ=0°to θ=80°f or two emission planes along the principal azimuths of the substrate, thus along and perpendicular to the Pdrows, i.e. along the [ ] 110 and [001] directions of the substrate, respectively. Integrated intensities for these two azimuths are shown as red and turquoise lines in panel (d), respectively. Again, the spectra are dominated by emission from Pd d-states between the Fermi-level at −4.0 eV and about −7.5 eV. Thus compared to experiment, the theoretical d band width is significantly larger, indicating an error of the HSE functional presumably due to ill-described electron-correlations for the Pd d-states. Note, however, that the predicted decrease in work function to 3.7 eV is in reasonable agreement with the experimental value of 4.0 eV, and that, regarding changes in the work function upon adsorption, the theoretical value of ΔΦ calc =−1.3 eV agrees even better with the measured value of ΔΦ exp =−1.1 eV.
Comparison of the photoemission data with results from the clean Pd(110) surface clearly identifies the peaks at −8.3 and −8.6 eV as molecular emissions. We suggest that the peak at −8.3 eV originates from π 1 , while the lower lying emission at −8.6 eV arises from π 2 . Below we outline the arguments leading us to this assignment, which will be further confirmed by the photoemission momentums maps (MMs) discussed in the subsequent section.
But before we continue, it is important to note that a similar analysis of the DOS for the across configuration reveals practically no energetic splitting of the two HOMO π states, neither in the freestanding layer, nor for the benzene/Pd(110) interface, thereby contradicting the experimental findings. Thus, together with the slightly less favorable total energy for the across orientation, this provides stronger proof that benzene indeed prefers to adsorb along the row direction on Pd(110). Once the orientation of the molecule has been clarified, the assignment of the emissions at −8.3 eV and −8.6 eV to π 1 and π 2 , respectively, is straightforward. Inspection of figure 1(a) demonstrates that π 1 has a nodal plane across the Pd-rows, while π 2 has a nodal plane along the rows. Due to symmetry selection rules, no emissions are expected along the nodal planes, hence π 1 is expected to be observable along the [ ] 110 , while π 2 is visible in the perpendicular direction. We emphasize that not only the size of the splitting, 0.3 eV in experiment, agrees reasonably well with the theoretical value of 0.2 eV, but also the energetic order of π 1 and π 2 is reproduced correctly by the calculations.

Photoemission momentum maps
The analysis of the photoemission MMs over a large angular range further confirms the interpretation that the two molecular emissions at −8.3 and −8.6 eV arise from a breaking of the degeneracy of the topmost π state and that the assignment to orbitals π 1 and π 2 is indeed correct. Panel (a) of figure 3 presents simulated MMs of the gas phase π 1 and π 2 orbitals for benzene oriented along the [ ] 110 -direction as depicted in figure 1(a) by assuming a plane-wave final state [3]. Clearly the nodal structure of the real-space orbital distributions, but also the shape of the orbitals, for example the more triangular-like shape of the π 1 lobes, is transferred also to momentum space. Generally speaking, the features in k-space appear quite broad compared to previously reported MMs of longer molecules, such as oligophenyls [3,58] or oligoacenes [5,59], owing to the smaller size of one single benzene ring and its lack of quasi-periodicity in the orbital distribution. As a consequence, emissions from small molecules such as benzene are expected to be more difficult to pinpoint in ARUPS experiments due to their greater spread in k-space causing a weaker angular dependence.
Panel (b) of figure 3 demonstrates how the adsorption geometry affects the momentum maps of π 1 and π 2 . While the internal bond length changes have a negligible effect on the photoemission angular distribution, the 8°t ilt angle of the benzene ring leads to a shift of the lobes in the [001] direction. Due the presence of mirror domains on the sample surface, the momentum maps need to be symmetrized by adding up the simulated maps for +8°and −8°tilt angle. The resulting momentum maps shown in panel (c) appear again very similar to the gas phase geometry, however, with a slightly larger extension of the k-space features in the [001] direction. Let us now turn to the experimental momentum maps shown figure 3(d) corresponding to the emissions at −8.3 (top) and −8.6 eV (bottom), respectively. The MM at a lower binding energy clearly resembles π 1 with its main lobes along the [ ] 110 azimuth. The weaker features along the [001]-directions most likely originate from the energetic overlap with the higher binding energy feature ( figure 2(d)). Conversely, at the higher binding energy, the contribution of π 1 is still quite pronounced, while the fingerprint of a π 2 emission, intensity along [001], has clearly grown.

Discussion
Although the agreement between the simulated (figure 3(c)) and measured (figure 3(d)) momentum maps may be not as convincing as in previous photoemission tomography studies for larger molecules with more repeat units, see e.g. [3][4][5][6]8], the assignment of orbitals of a molecule as small as benzene is still possible. This allows a number of conclusions to be drawn.
First, we observe a clear energy splitting of the π 1 and π 2 states, which comprise the formerly degenerate HOMO of the gas phase benzene ring. This splitting has been overlooked in the early photoemission study of benzene/Pd(110) [29] and could be revealed in our work, despite worse energy resolution, because of (a) the availability of ARUPS data over a wide angular range and (b) the predictions of the angular emission patterns of π 1 and π 2 .
The energy splitting of π 1 and π 2 can also be used to determine the azimuthal orientation of benzene. The DFT total energy results alone do not give a compelling result for the orientation, because the total energy difference between the along and across alignment of only 0.1 eV is within the typical error bar expected from uncertainties in the exchange-correlation treatment particularly including van-der-Waals corrections. Only a combination of DFT and ARUPS leads to an unambiguous determination of the orientation. The most compelling evidence for the along orientation is thus the splitting of π 1 and π 2 which is only reproduced in DFT by this along orientation. Importantly, this splitting is rather insensitive to such intrinsic uncertainties of DFT and it can be understood by simple arguments with π 2 undergoing a strong bond-stabilization in this orientation. Inspecting the along adsorption structure of figure 1(a), we see that π 2 will interact with two adjacent Pd atoms of the close-packed Pd-rows in a bonding manner, while π 1 has a nodal plane between the Pd atoms. As a consequence, the interaction of π 2 with the substrate can be expected to be more pronounced.
The interaction with the substrate may also be responsible for the observed differences of the simulated and measured momentum maps of p 1 and π 2 . To start with, owing to the hybridization with the substrate, the molecular features are significantly broadened in energy, and as a consequence the signature of π 2 is still visible at the peak energy of π 1 and vice versa. However, there are also further deviations from the simulated momentum maps in terms of position, shape and relative intensities of the angular patterns. Overall, the intensity of π 2 is much smaller compared to π 1 . This may be related again to the fact that π 2 interacts more intense with the substrate, leading to a larger energy broadening and a concomitant reduction of spectral weight at a given constant binding energy. Also the differences in the k P position and the shape of the angular patterns may be-at least partly-ascribed to the hybridization with the substrate. The maximum of the MM of π 1 is located at nm −1 . Shifts of comparable magnitude have in fact been observed in our previous PT study of pentacene/Cu(110), which we have attributed to a spatial distortion of the orbital on account of the molecule-metal hybridization [60]. Moreover, the fact that the experimental MM of π 2 appears much sharper in the k [001] -direction than the simulated MM may aslo originate from molecule-substrate interaction. Not only does it lead to energetically broadened molecular resonances, but it also causes an enhanced intermolecular band dispersion which in turn results in substructures and sharper features in constant-binding energy momentum maps [60][61][62].
It would be advantageous to substantiate this interpretation by simulations of momentum maps for benzene interacting with Pd(110) rather than benzene in the gas phase, that is, by using the wave functions of the benzene/Pd(110) interface as intitial states in the simulation of the photoemission intensity. We have demonstrated that this may indeed improve the agreement with measured momentum maps [60][61][62]. For molecule/metal interfaces with considerable interactions, however, in the present case such simulations are hampered by the fact that the DFT electronic structure has a severe shortcoming. More precisely, when we take the Kohn-Sham eigenvalues as one-particle excitation energies, the π 1 and π 2 states are energetically overlapping with the Pd d-band, while experimentally they are clearly below the Pd d-band emissions, as can be seen from figure 2. As a consequence, the interacting strength of π 1 and π 2 with substrate states is significantly overestimated in the calculation, preventing a meaningful comparison of simulated momentum maps for the benzene/Pd(110) interface with experiment. Note that the DOS shown in figure 2 already results from a hybrid functional calculation, which in the case of benzene/Pd(110), provides no real improvement over GGA results. Here, the main problem seems to be related to ill-described correlations, resulting in a too large band width of the Pd d-band.
Let us finally return to the question of how important a correct description of the final state is for describing the photoelectron angular distribution [26]. Benzene on Pd(110) differs from the systems that have so far been investigated with the PT approach in at least two aspects. First, benzene (C 6 H 6 ) is considerably smaller than the molecules studied so far, such as pentacene (C 22 H 14 ), sexi-phenyl (C 36 H 26 ), perylene-tetracarboxylic dianhydride (PTCDA, C 24 H 8 O 6 ) or tetraphenyl-porphyrin (C 44 H 30 N 4 ), where the plane wave final state with the isolated molecules' initial states has yielded a very good description of the momentum maps. This may be due to an averaging effect in the final state substrate scattering as the intra-molecular orbital electron density in the larger molecules is not commensurate with the substrate. A second difference in the current system regards the notably more reactive Pd(110) surface compared to the coinage metal substrates used in previous PT work. Therefore, one might expect that the isolated molecule's orbitals may no longer be a good approximation for the initial state. Nevertheless, the plane wave final state description with an isolated molecule initial state still proves useful. In conjunction with the DFT results, it allows for an understanding of the global picture and even for some important conclusions to be drawn regarding the geometric (orientation of the molecule) and electronic structure (lifting of degeneracy and splitting π-states) of the adsorbate. This may be related to the experimental geometry in which we use p-polarized light and collect electrons in the half-plane where the angle between the polarization vector A and the emission direction k is small. Under these conditions, the plane wave final state approaches the independent atomic center method where spherical waves are used as final state [3,63]. Clearly for a better modeling of the PE distribution of strongly interacting adsorbates improvements in the treatment of both the initial and final state are required. We believe that for the present system, the major factor leading to the apparent breakdown of PT lies in the treatment of the initial state rather than in the final state. This is indicated by π 1 's distribution being reasonably described by PT with the isolated molecules initial state while the more strongly interacting π 2 is not.
To guide future theoretical developments in describing the initial state of strongly adsorbed systems and going beyond a plane wave final state approximation, more ARUPS data is needed for small-molecule/metal systems and other experimental geometries. In particular, as suggested by Bradshaw and Woodruff [26], experimental geometries where the polarization vector is held parallel to the emission direction (  A k) will be important to highlight errors arising from the plane wave final state approximation. Moreover, experiments for molecules of increasing length, such as linear polycyclic hydrocarbons, where benzene is the starting point, will be useful to shed light on the influence of final state scattering. Finally, experiments on substrates with varying interaction strength will serve as a benchmark for advanced electronic structure methods which are necessary for an accurate description of the initial state.
Concluding our paper with a philosophical comment, we stress that photoemission tomography in its present state using a plane wave final state is surely not 'correct', but it is, nevertheless, simple and predictive. We perhaps should not forget Popper's notion that 'Science may be described as the art of systematic over-simplification -the art of discerning what we may with advantage omit.' [64], or George Box's famous quote 'Essentially, all models are wrong, but some are useful.' [65].