Momentum space imaging of σ orbitals for chemical analysis

Tracing the modifications of molecules in surface chemical reactions benefits from the possibility to image their orbitals. While delocalized frontier orbitals with π character are imaged routinely with photoemission orbital tomography, they are not always sensitive to local chemical modifications, particularly the making and breaking of bonds at the molecular periphery. For such bonds, σ orbitals would be far more revealing. Here, we show that these orbitals can indeed be imaged in a remarkably broad energy range and that the plane wave approximation, an important ingredient of photoemission orbital tomography, is also well fulfilled for these orbitals. This makes photoemission orbital tomography a unique tool for the detailed analysis of surface chemical reactions. We demonstrate this by identifying the reaction product of a dehalogenation and cyclodehydrogenation reaction.


INTRODUCTION
In frontier orbital theory, orbital patterns are regarded as important for chemical empiricism-they determine molecular reactivities and reaction pathways (1,2). Orbitals are attributed as a reality of their own besides the total electron densities in which influential electronic structure theories are framed (3,4). This has spurred many successful efforts to image orbitals (5-7) despite subtleties of the single-electron orbital concept (8,9). Orbital imaging has also been applied to the investigation of chemical bonding (10)(11)(12). The imaging of orbital-derived local densities of states with lowtemperature scanning tunneling microscopy (STM) (6,10,11) is particularly attractive because it combines spectroscopic energy resolution with sub-angstrom spatial resolution, and the sensitivity to orbital patterns can be improved by functionalizing the STM tip (6,13). Furthermore, time resolution in STM-based orbital imaging has been achieved (14).
There is, however, a fundamental drawback of STM-based orbital imaging. The method is restricted to orbitals close to the chemical potential [Fermi energy (E F )], typically in an energy range E F ± 2 eV. In many cases, however, it would be advantageous to be able to identify, analyze, and image deeper-lying orbitals. A case in point are  orbitals. Because their electron density is concentrated along the atom-atom interconnections in molecules, they play an important role in the local bonding between neighboring atoms. The capability to measure  orbitals as fingerprints of local chemical structures is therefore critical.
Photoemission orbital tomography (POT) (7) is a recent technique for the orbital analysis and orbital imaging of molecules on surfaces. For example, valence spectra have been deconvolved modelfree into orbital projected densities of states (pDOS) (15,16), orbitals have been reconstructed from experimental data in two dimensions (2D) and 3D (17)(18)(19), and orbital patterns in momentum space have been recorded with femtosecond time resolution (20). While it has long been recognized that angle-resolved photoelectron spectroscopy (ARPES) can image orbitals in the molecular frame, the alignment of molecules in the gas phase is a challenge (21)(22)(23). In contrast, in POT, photoelectrons are collected in the half-space above the sample surface on which the adsorbed molecules are fixed in space (7,(17)(18)(19). The measured angular distribution of photoelectrons is governed by the spatial distribution of electrons in the initial state-the orbital. There is overwhelming evidence that, for  orbitals in conjugated molecules, this relation is a straightforward Fourier transform (24,25). This is valid if the final state of the photoelectron after photoemission can be represented by a plane wave (PW) [so-called plane wave approximation (PWA)].
Because it is based on ARPES, POT per se is not limited to a particular binding energy range if the photon energy is large enough. However, there is an important caveat. So far, the PWA has only been confirmed experimentally for  orbitals close to E F . For -conjugated organic molecules, these originate from linear combinations of only 2p z -orbitals on identical (all carbon) atoms. Under these circumstances, the PWA is equivalent to the independent atomic center (IAC) approximation (7,25) in which the photoelectron wave function is described in terms of unbound solutions of the atomic Schrödinger equation with energy E kin summed over all atoms of the molecule and taken in the asymptotic limit for infinite distances of the detector from the photoemitter (see Materials and Methods for details) (26). However, this equivalence relies on photoemission from only one type of atomic orbital. If, on the other hand, the photoemission occurs from several different atomic orbital types, for example, the carbon 2s, 2p x , and 2p y -orbitals of which  orbitals consist, then the PWA may become less accurate because the contributions of the different atomic orbitals do not any more assemble to the simple PW final state that essentially Fourier transforms the initial state orbital. Therefore, the applicability of the PWA to the final state of the photoemission needs to be questioned in this case.
This conceptual issue of POT for  orbitals is compounded by an experimental one. Because of the typically larger wave vectors of (the topmost)  orbitals, larger photon energies must be chosen, such that the characteristic emission patterns of these orbitals are still included in the photoelectron horizon. However, at the higher wave vectors and higher binding energies of the  orbitals, the intensity of substrate emissions is large, possibly causing a difficulty to distinguish orbital tomograms of the  orbitals from the background.
Here, we address by experiments the question of whether the angular distributions of  orbitals can be reliably measured and explained in terms of the PWA. If this was the case, then POT could straightforwardly be extended to the identification, analysis, and imaging of  orbitals and their chemical phenomenology. This would open a much more direct way to study chemical bonding than the commonly used x-ray photoelectron spectroscopy (XPS), which provides indirect information on local bonding by measuring the influence of these bonds on the binding energies and spectral shapes of atomic core levels. With -based POT, one would gain the possibility to access the bonding orbitals directly.
Our model is the on-surface dehalogenation and cyclodehydrogenation of 10,10′-dibromo-9,9′-bianthracene (DBBA). DBBA and similar molecules are often used to create carbon-based nanostructures by on-surface synthesis (27)(28)(29). The initially dehalogenated positions of the nonplanar precursor molecule DBBA not only can directly engage in an Ullmann-type C─C coupling (30,31) but also can become metalated or hydrogenated. Because the products of this reaction have an influence on the structural quality of the targeted nanostructures (28), the analytical discrimination between these competing reactions and products is important, and we used it as a test case for POT on  orbitals. The size of the molecular products suppresses photoelectron scattering by atoms of the metal substrate (32) because many atoms of the molecule are typically out of registry with the substrate, thus safeguarding the PWA as far as possible. Our model system is therefore a promising candidate to achieve our conceptual goal of proving momentum space imaging of  orbitals for chemical analysis and, at the same time, to address an actual scientific problem in the field of on-surface dehalogenation and cyclodehydrogenation reactions.

RESULTS
We deposited DBBA ( Fig. 1, 1) on Cu(110) (12,33). DBBA undergoes catalytic dehalogenation and cyclodehydrogenation reactions at 525 K (12,33), during which the adsorbate planarizes and the  conjugation expands to the entire fused carbon backbone, leading to the formation of a bisanthene-like species that is adsorbed flat on the surface-possible products in agreement with STM and XPS (12,33) are shown in Fig. 1. While a direct chemical bonding between the underlying substrate atoms and the carbon atoms in the 10,10′ positions ( Fig. 1, 2) (33), or even all carbon atoms along the zig-zag edge of bisanthene ( Fig. 1, 3), was excluded in a previous study based on the analysis of frontier  orbitals (12), the discrimination between bisanthene (Fig. 1, 4) and metalated bisanthene (Fig. 1, 5) has not yet been achieved. Several examples of metalmolecule complexes were reported in similar reactions (28,31,34), and hence, it is not clear how the valency of carbon atoms in question will be saturated. Both hydrogen and copper adatoms are abundantly available and mobile at elevated temperatures in the reaction environment on the surface.
In a first attempt to discern between products 4 and 5, we computed the adsorption structures and electronic properties of both species on Cu(110) by van der Waals-corrected density functional theory (DFT). We then calculated the pDOS by projecting the eigenstates of the DFT calculation of each molecule/substrate system onto the eigenstates (orbitals) of the corresponding isolated molecule (see Materials and Methods for details). We found that neither the pDOS in the energy window between the Fermi energy and the onset of the Cu d-band nor the theoretical k ∥ maps for the lowest unoccupied molecular orbital (LUMO) and highest occupied molecular orbital (HOMO), calculated by Fourier transforming the orbitals of the isolated molecules (see Materials and Methods for details), could discriminate between the two scenarios when compared to the experiment ( fig. S1). Note that upon adsorption on Cu(110), the LUMO of the bisanthene-like species becomes occupied by charge transfer from the metal (12), making it accessible to POT, as fig. S1 proves.
Crucially, Fig. 2A shows that one of the two uppermost  orbitals is strongly affected by the replacement of the hydrogen in the 10,10′ positions by Cu adatoms, leading to the expectation that products 4 and 5 could possibly be discriminated by measuring  orbitals.
To check whether  orbitals are at all detectable in POT, we recorded photoemission band maps I(E b , k), i.e., photoemission intensity I as a function of binding energy E b and wave vector k, in two azimuthal directions up to a binding energy of E b ≃ 13 eV (Fig. 2, B and C). We know from previous work (12,33) that the bisanthene-like product of the DBBA dehalogenation and cyclodehydrogenation is orientated with its zig-zag edges along the substrate [001] direction and with the armchair edges along [1 _ 1 0]. In addition to the frontier  orbitals, one can clearly identify in Fig. 2 (B and C) several molecular emissions below the d-band: a  band dispersing between E b ≃ 4 to 10 eV at |k ∥ | ≤ 1.5 Å −1 and a band between E b ≃ 5 to 13 eV at higher |k ∥ |. The large |k ∥ | of the latter suggest that it might be of  character.  1 0]). We therefore labeled the orbitals according to their type ( or ) and the number of nodal planes (0,1,2...) along k x and k y . For instance, the HOMO (2,3) exhibits two (three) nodal planes along k x (k y ), while the (7,3) orbital has seven and three nodes, respectively. (4,0) is the LUMO.
Next, we compared the theoretical k ∥ maps in Fig. 3 to experimental ones generated from the photoemission intensity data cube I(E b , k ∥ ) (movie S1) at fixed binding energies E b (see Materials and Methods). This comparison is well founded for  orbitals but still to be substantiated for  orbitals because of the uncertain applicability of the PWA for the latter. To search for individual orbitals in the experimental data cube, we used a momentum space deconvolution with theoretical k ∥ maps (figs. S3 to S7) (15). Details of the deconvolution are explained in Materials and Methods. Of the 15  orbitals and 27  orbitals in Fig. 3, we found 12  and 18  orbitals in the range from E b = 0 to 10 eV. Only nine  orbitals were not identified with certainty. The large number of identified  orbitals in excellent agreement with theory confirms that POT can be applied not only to p z -derived but also to s-, p x -, and p y -derived orbitals. Given the mathematical basis of the PWA (7, 26) as mentioned above, this is a remarkable finding for which no theoretical explanation has been given yet. Nevertheless, the experiment speaks for itself.

DISCUSSION
In the remainder of the paper, we focus on the two uppermost  orbitals just below the Cu d-band, using them to discriminate between products 4 and 5 of the surface reaction in Fig. 1. In Fig. 4A, the experimental k ∥ map at E b = 5.16 eV is displayed. Its pattern is complex and shows emissions from several states:  orbitals and  orbitals (circles), as well as high-E b tails of the Cu d-band (triangles). The strong emission lobes at |k ∥ | ≥ 2 Å −1 can be accounted for by the two uppermost  orbitals (7,3) (red) and (0,8) (blue), whose real space distributions are displayed in Fig. 2A. As can be discerned from their theoretical k ∥ maps (Fig. 4B, inset), the emissions at k ∥ = (0, ±2.89) Å −1 arise from (0, 8), while those at (±2.56, ±1.39) Å −1 can be from both (0,8) and (7,3). To quantitatively decompose the experimental k ∥ maps for varying E b , we again used the momentum space deconvolution with theoretical k ∥ maps (15) to extract the E b -dependent contributions a i of the component orbitals to the total POT intensity distribution I(E b , k ∥ ). Figure 4B shows the a i (E b ) for three orbitals: (0,8), (7,3), and (0,3). They have the meaning of an experimentally derived pDOS and can be compared to their theoretically calculated counterparts.
In Fig. 4 (C and E), we show simulated k ∥ maps for bisanthene (4) and metalated bisanthene (5) on Cu(110) based on the van der Waals-corrected DFT calculation of the complete molecule/substrate system. These maps have been obtained with a damped PW as the final state, which prevents the overrepresentation of bulk states from the substrate (35). For bisanthene, the six bright emission lobes at |k ∥ | ≥ 2 Å −1 (Fig. 4C) correctly reproduce the experimental pattern. In contrast, for metalated bisanthene, the four emission lobes of the (7,3) orbital dominate, while the emission from (0,8) is barely visible (Fig. 4E). This already leads us to conclude that bisanthene agrees better with the experimental data, and the metalation with Cu adatoms is unlikely.
This conclusion is further corroborated by comparing the calculated pDOS of (0,8) and (7,3) in Fig. 4 (D and F) to the corresponding experimental ones in Fig. 4B. For bisanthene, the calculated pDOS for both orbitals exhibit well-defined peaks with similar intensities and energies, therein resembling closely the experiment. In contrast, the calculated (0,8) pDOS for metalated bisanthene is spread out broadly in energy with its maxima well separated from that of  (7,3). This indicates that, in the van der Waals-corrected DFT calculation of the complete molecule/substrate system, the (0,8) orbital of metalated bisanthene hybridizes strongly with the substrate, presumably through the Cu adatoms, which are fully integrated into the lobe structure of (0,8) [but, notably, not of (7,3)] and which, upon adsorption, couple strongly to the metal surface. Because such a hybridization is not observed in the experimental pDOS of (0,8) (Fig. 4B), we can definitively rule out the metalated product 5. To cross-check, we repeated the deconvolution of the experimental data cube I(E b , k ∥ ) with theoretical k ∥ maps of free metalated bisanthene ( fig. S9). This delivers a very similar result to the experimental pDOS in Fig. 4B, which, however, is Binding energy (eV) inconsistent with the calculated pDOS for the metalated species in Fig. 4F. Thus, only bisanthene as a product of the surface reaction allows a consistent interpretation of our  orbital data. The analysis of the deep lying (0,3) orbital is in full accord with this result (see the Supplementary Materials).
In conclusion, we advanced photoemission orbital tomography, first by substantially extending the accessible binding energy range, particularly far below the d-band of the substrate, and then by successfully detecting  orbitals. In particular, we found that the PWA for the final state is also suited for  orbitals. For the product of the catalytic dehalogenation and cyclodehydrogenation of DBBA on Cu(110), we demonstrated that almost the complete orbital spectrum in a wide binding energy range is detectable. Because-depending on orbital shape, adsorption site, and registry with the sample-it is usually only a small number of orbitals that provide key information regarding chemical modifications, the ability to image the complete spectrum of valence orbitals with POT in an unrivaled binding energy range provides an invaluable advantage for the analysis of surface chemical reactions. Regarding DBBA/Cu(110), we found that the product of the catalytic dehalogenation and cyclodehydrogenation reactions is bisanthene. We anticipate that the absence of a metalation with surface adatoms at the dehalogenated positions of the precursor will aid the formation of high-quality nanostructures from this surface chemical reaction.

Sample preparation
Our experiments were performed in ultrahigh vacuum (~10 −10 mbar). The Cu(110) single crystal was cleaned by several cycles of sputtering by Ar + ions at 1 keV and subsequent annealing at 800 K. A film of the DBBA precursor (Sigma-Aldrich, CAS number 121848-75-7) was deposited by evaporation from a molecular evaporator (KENTAX GmbH) onto the crystal surface held at room temperature. Subsequently, the sample was annealed at 525 K to trigger the chemical reaction.

Photoemission experiments: Band maps and k ∥ maps
Photoemission experiments were conducted at the Metrology Light Source insertion device beamline of the Physikalisch-Technische Bundesanstalt (Berlin, Germany) (36). p-polarized ultraviolet light (photon energy = 57 eV) with an incidence angle of 40° to the surface normal was used. In this geometry, the A∥k condition (25), where A is the vector potential of the incident light and k is the wave vector of the photoelectrons, is approximately fulfilled for most molecular emissions in forward direction. The significance of this condition is described below. Two different types of photoemission experiments were conducted using the toroidal electron analyzer (37). To obtain the experimental band maps I(E b , k 1 _ 1 0 ) and I(E b , k 001 ), the photoemission intensity was recorded in the emission angle range from −85° to +85° along the [  function of kinetic energy E kin . The latter was converted to binding energy E b by fitting the experimental photoemission spectra around the chemical potential by the Fermi function to define the binding energy E b = 0 on the experimental E kin scale. Experimental k ∥ maps at a fixed binding energy E b were obtained from the 3D data cube I(E b k ∥ ), i.e., the intensity I of photoemission as a function of binding energy E b and parallel momentum k ∥ , recorded while rotating the sample around its normal in 1° steps. In this way, the full photoemission intensity distribution in the k ∥ plane perpendicular to the sample normal was recorded.

DFT calculations: Computational details
The geometry and electronic structure of the free and surfaceadsorbed molecules were calculated in the framework of DFT. We conducted density functional calculations for both gas phase molecules and for molecules adsorbed on the Cu(110) surface. For the former, we used the quantum chemistry package NWChem (38), while for the latter, we used a repeated slab approach and the Vienna Ab initio Simulation Package (VASP) (39)(40)(41) program.
The results of the gas phase calculations including the molecular orbitals and the corresponding theoretical k ∥ maps (see below) are available via a web-based database (42). In particular, the molecular orbitals of bisanthene (C 28 H 14 , 4) and metalated bisanthene (C 28 H 12 Cu 2 , 5) using a Perdew-Burke-Ernzerhof-generalized gradient approximation (PBE-GGA) exchange-correlation functional are accessible from this web interface using the database IDs 406 and 407, respectively. Real-space representations of bisanthene orbitals obtained with PBE-GGA with binding energies between 4 and 13 eV with respect to the vacuum level are displayed in fig. S2. In addition, we also computed the gas phase molecular orbitals of bisanthene using the range-separated hybrid functional HSE06 [ID = 480 in (42)] and the long-range corrected hybrid LRC-wPBEh [ID = 486 in (42)].
For the full molecule/Cu(110) systems, we applied the repeated slab approach. Here, the Cu(110) substrate was modeled with five atomic layers, a lattice parameter of a = 3.61 Å and a vacuum layer of at least 17 Å between the slabs to avoid spurious electric fields (43). The most favorable adsorption site for both molecules, bisanthene (C 28 H 14 , 4) and metalated bisanthene (C 28 H 12 Cu 2 , 5), was determined by testing several high-symmetry adsorption sites (hollow, top, short bridge, and long bridge) in a local geometry optimization approach, allowing all molecular degrees of freedom and the topmost two Cu layers to relax until forces were below 0.01 eV/Å. For these geometry optimizations, we used the GGA according to PBE (44) for the exchange-correlation potential with the D3 correction for van der Waals interactions (45). The projector augmented wave method (41,46) was used with a PW cutoff of 500 eV and a Monkhorst-Pack 3 × 3 × 1 grid of k-points with a first-order Methfessel-Paxton smearing of 0.2 eV. On the basis of the relaxed adsorption geometries, which were found to be the short-bridge position for bisanthene 4 and the hollow site for metalated bisanthene 5 (figs. S13 and S14), we analyzed the electronic structures in terms of the molecular orbital pDOS and simulated the photoelectron angular distribution (k ∥ maps) as described in detail below. Note that, for the pDOS, we also used the range-separated hybrid functional HSE06 (47,48) for comparison (see below).

DFT calculations: Molecular orbital pDOS
We analyzed the electronic structures of the molecule/metal interfaces by computing the molecular orbital pDOS (35). In addition to the full interface system with the Kohn-Sham states  nq , we performed a DFT calculation for a free-standing layer of molecules, cut out from the full interface, leading to the molecular orbital states ϕ iq . We then computed the molecular orbital pDOS ρ i by projecting the Kohn-Sham orbital of the full system  nq with Kohn-Sham energy ε nq onto a given molecular state ϕ iq and summing over the Brillouin zone and all states n of the full system In addition to the pDOS obtained from the PBE-GGA functional as shown in Fig. 4 (D and F), we also performed self-consistent DFT calculations using the range-separated hybrid functional HSE06 (47,48). The results are depicted in fig. S10 and fully support our conclusions drawn from the PBE-GGA results. Compared to the pDOS of the PBE-GGA calculation, the HSE06 functional leads to an overall shift to larger binding energies due to partial correction of self-interaction errors and thereby further improves the agreement with experiment, displayed in Fig. 4B. The pronounced difference between the two systems, C 28 H 14 /Cu(110) and C 28 H 12 Cu 2 /Cu(110), with regard to the two topmost  orbitals, (7,3) and (0,8), remains almost identical to the PBE-GGA calculation, thereby fully supporting our conclusions.

DFT calculations: Simulation of k ∥ maps
Our starting point for simulating the angle-resolved photoelectron intensity was the one-step model of photoemission in which the intensity I(E kin , k ∥ ) is given by Fermi's golden rule (49) Here, k ∥ = (k x , k y ) is the component of the photoemitted electron's wave vector parallel to the surface, which is related to the polar and azimuthal emission angles θ and ϕ by k x =|k| sin cos (3) k y =|k| sin sin (4) where k is the wave vector of the emitted electron, and E kin = ℏ 2 |k| 2 _ 2m , where ℏ is the reduced Planck constant and m is the electron mass. The photocurrent of Eq. 2 is given by a sum over all transitions from occupied initial states i described by wave functions  i to the final state  f characterized by the direction (θ and ϕ) and the kinetic energy of the photoemitted electron. The δ function ensures energy conservation, where  denotes the sample work function, E i denotes the binding energy of the initial state, and ℏω denotes the energy of the exciting photon. The transition matrix element is given in the dipole approximation, where p and A, respectively, denote the momentum operator of the electron and the vector potential of the exciting electromagnetic wave (50).
In the POT approach (25,50), the final state  f is approximated by a PW. Thereby, the photocurrent I i arising from one particular initial state i turns out to be proportional to the Fourier transform ˜  i (k) of the initial-state wave function multiplied by the polarization factor |A • k| 2 Regarding the applicability of the PW final-state approximation, it has been argued (7,51) that it can be expected to be valid if the following three conditions are fulfilled: (i)  orbital emissions from large planar molecules, (ii) an experimental geometry in which the angle between the polarization vector A and the direction of the emitted electron k is rather small (i.e., A∥k condition approximately fulfilled), and (iii) molecules consisting of many light atoms (H, C, N, and O). Experience on a variety of -conjugated molecules adsorbed on various substrates has shown (25,(52)(53)(54) that conditions (ii) and (iii) seem to be well satisfied within our experimental setup. Regarding the relaxation of condition (i) toward the inclusion of  orbitals, however, only little information has been gathered so far. Here, it should be noted that the arguments leading to (i) involve a comparison of the PW final-state approximation with the more sophisticated IAC approximation. In the IAC, the final state is treated as a coherent sum of atomic-like scattering states. However, when the initial state, written as a linear combination of atomic orbitals, consists of only a single atomic orbital type, as is the case for  orbitals, it can be shown that the IAC strictly reduces to the PW final-state description, because the atomic p z -orbital can be factored out of the respective summation [see the supplementary materials of (7)]. In contrast, for  orbitals, a mathematical equivalence of the IAC and PW final-state descriptions cannot be shown. This does not a priori preclude the applicability of the simple and intuitive PW approximation but asks for an experimental test.
Using the PW final-state approximation, we computed the expected (theoretical) photoemission k ∥ maps I i (k ∥ ) for an orbital i according to Eq. 5, both for gas phase molecules and for the full interacting systems of the molecules adsorbed on Cu(110). For gas phase orbitals, the construction of k ∥ maps can also be visualized in an intuitive manner ( fig. S11). Starting from the orbital in real space (x, y, z) ( fig. S11, A and D), first, the momentum space representation of the orbital ˜  (k x , k y , k z ) is computed by a 3D Fourier transform. Then, a hemispherical cut at |k|= √ _ 2m E kin / ℏ 2 ( fig. S11, B and E) and multiplication with the polarization factor according to Eq. 5 lead to the k ∥ map (fig. S11, C and F). As illustrated for the (0,3) (top row in fig. S11) and (7,3) (bottom row in fig. S11), the main features in  orbitals are typically located at somewhat larger k ∥ values than the ones of  orbitals, and, at the same time, the larger kinetic energies lead to a reduced photoemission intensity.
Photoemission k ∥ maps were also simulated for molecules adsorbed on Cu(110). Here, the starting point was Eq. 2, in which the initial states  i were taken as the Kohn-Sham orbitals of the interacting molecule/Cu(110) system, and the final state  f was again approximated by a PW that was, however, modified by an exponential damping factor with a damping constant of  = 0.5 Å −1 to mimic the finite mean free path of the photoemitted electrons. To achieve an acceptable momentum resolution in the resulting k ∥ maps, the sampling of the Brillouin zone was increased to a k-grid of 5 × 5 × 3. Further details about this approach have already been published elsewhere (35).

Deconvolution of the experimental data cube
In the orbital deconvolution procedure (15,50), we made use of the energy and momentum dependence of photoemission to deconvolve experimental data into individual orbital contributions. This provides an orbital-by-orbital characterization of the experimentally determined density states, i.e., an experimental pDOS that can be readily compared to a computed pDOS.
We obtained the experimental pDOS by a least squares minimization procedure. The experimental photoemission data can be viewed as a data cube I exp (E b , k ∥ ) = I exp (E b , k x , k y ): the photoemission intensity I exp measured as a function of the two momentum components parallel to the surface, k x and k y , and the binding energy E b (corresponding to the kinetic energy E kin as described above). The relation between k x and k y on one hand and emission angles (θ and ϕ) of the photoelectrons on the other follows from Eqs. 3 and 4. The deconvolution of the experimental data cube then consists in minimizing the squared differences between the experimental and theoretical k ∥ maps  2 ( a 1 , a 2 , … , a n ) = ∫ d k ∥ [ by adjusting the n weights a i (E b ) of all orbitals i with the theoretical k ∥ maps I i (k ∥ ) that are allowed to contribute to the measurement data. Because the minimization is performed for each binding energy E b separately, one thereby obtains an orbital pDOS given by the weight functions a i (E b .).
From Eq. 6, it is clear that the so-obtained experimental pDOS will also depend on the set of theoretical k ∥ maps I i (k ∥ ) that are used in the deconvolution procedure. It is therefore important to check how sensitive these theoretical k ∥ maps are with respect to the choice of the exchange-correlation functional. Fortunately, it turned out that these k ∥ maps are robust and remain almost unaffected by the particular choice for the exchange-correlation functional. This is illustrated in fig. S12 by computing the overlaps of bisanthene orbitals generated with the PBE-GGA functional [ID = 406 in (42)] with those obtained from the range-separated hybrid functional HSE06 [ID = 480 in (42)]. More precisely, fig. S12 shows the quantity in a gray-scale density map, where the darkest values indicate a maximum deviation of less than 9% from perfect overlap. It should be noted that, for most of the orbitals, the similarity of the PBE-GGA and HSE06 orbitals is even more pronounced and that the resulting theoretical k ∥ maps from these two functionals turn out to be virtually indistinguishable. We also note that a similar finding was made when comparing the PBE-GGA orbitals with those of the longrange corrected hybrid LRC-wPBEh [ID = 486 in (42)] or the global hybrid PBE0 [ID = 482 in (42)]. We are thus confident that the experimental orbital pDOS obtained from the photoemission data cube is indeed a meaningful and robust quantity.

SUPPLEMENTARY MATERIALS
Supplementary material for this article is available at https://science.org/doi/10.1126/ sciadv.abn0819 View/request a protocol for this paper from Bio-protocol.