# Multiorbital exciton formation in an organic semiconductor

Wiebke Bennecke  I. Physikalisches Institut, Georg-August-Universität
Göttingen, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany  Andreas
Windischbacher  Institute of Physics, University of Graz, NAWI Graz,
Universitätsplatz 5, 8010 Graz, Austria  David Schmitt  Jan Philipp Bange  I.
Physikalisches Institut, Georg-August-Universität Göttingen, Friedrich-Hund-
Platz 1, 37077 Göttingen, Germany  Ralf Hemm  Department of Physics and
Research Center OPTIMAS, University of Kaiserslautern, Erwin-Schrödinger-
Straße 46, 67663, Kaiserslautern, Germany  Christian S. Kern  Institute of
Physics, University of Graz, NAWI Graz, Universitätsplatz 5, 8010 Graz,
Austria  Gabriele D’Avino  Xavier Blase  Univ. Grenoble Alpes, CNRS, Inst
NEEL, F-38042 Grenoble, France  Daniel Steil  Sabine Steil  I. Physikalisches
Institut, Georg-August-Universität Göttingen, Friedrich-Hund-Platz 1, 37077
Göttingen, Germany  Martin Aeschlimann  Department of Physics and Research
Center OPTIMAS, University of Kaiserslautern, Erwin-Schrödinger-Straße 46,
67663, Kaiserslautern, Germany  Benjamin Stadtmüller  Department of Physics
and Research Center OPTIMAS, University of Kaiserslautern, Erwin-Schrödinger-
Straße 46, 67663, Kaiserslautern, Germany  Institute of Physics, Johannes
Gutenberg-University Mainz, 55128 Mainz, Germany  Marcel Reutzel  I.
Physikalisches Institut, Georg-August-Universität Göttingen, Friedrich-Hund-
Platz 1, 37077 Göttingen, Germany  Peter Puschnig  Institute of Physics,
University of Graz, NAWI Graz, Universitätsplatz 5, 8010 Graz, Austria  G. S.
Matthijs Jansen<EMAIL_ADDRESS>I. Physikalisches Institut,
Georg-August-Universität Göttingen, Friedrich-Hund-Platz 1, 37077 Göttingen,
Germany  Stefan Mathias<EMAIL_ADDRESS>I. Physikalisches
Institut, Georg-August-Universität Göttingen, Friedrich-Hund-Platz 1, 37077
Göttingen, Germany  International Center for Advanced Studies of Energy
Conversion (ICASEC), University of Göttingen, Göttingen, Germany

###### Abstract

Harnessing the optoelectronic response of organic semiconductors requires a
thorough understanding of the fundamental light-matter interaction that is
dominated by the excitation of correlated electron-hole pairs, i.e. excitons.
The nature of these excitons would be fully captured by knowing the quantum-
mechanical wavefunction, which, however, is difficult to access both
theoretically and experimentally. Here, we use femtosecond photoemission
orbital tomography in combination with many-body perturbation theory to gain
access to exciton wavefunctions in organic semiconductors. We find that the
coherent sum of multiple electron-hole pair contributions that typically make
up a single exciton can be experimentally evidenced by photoelectron
spectroscopy. For the prototypical organic semiconductor buckminsterfullerene
(C60), we show how to disentangle such multiorbital contributions and thereby
access key properties of the exciton wavefunctions including localization,
charge-transfer character, and ultrafast exciton formation and relaxation
dynamics.

##  I MAIN

Excitons, quasiparticles consisting of bound electron-hole pairs, are at the
heart of the optoelectronic response of all organic semiconductors, and
exciton formation and relaxation processes are largely responsible for energy
conversion and light harvesting applications in these materials. At the atomic
level, excitons are described by a two-particle correlated quantum-mechanical
wavefunction that includes both the excited electron and the remaining hole.
This wavefunction covers the complete shape of the exciton wave and thus
provides access to a number of critical exciton properties such as the orbital
character, the degree of (de)localization, the degree of charge separation,
and whether this involves charge transfer between molecules. Consequently, in
order to fully understand exciton dynamics and to exploit them in, e.g., an
organic solar cell, an accurate and complete measurement of the exciton
wavefunction would be ideal. Exemplary in this situation is the ongoing work
to understand the optoelectronic response of C60, a prototypical organic
semiconductor that is commonly used in organic solar cells ke_efficient_2015;
puente_santiago_tailoring_2020; yu_simplified_2020. Here, a topic of research
has been the optical absorption feature that occurs at 2.8 eV for multilayer
and other aggregated structures of C60 wang_aggregates_1993. Interestingly,
time- and angle-resolved photoelectron spectroscopy and optical absorption
spectroscopy studies have indirectly found that this optical transition
corresponds to the formation of charge-transfer excitons with significant
electron-hole separation stadtmuller_strong_2019; emmerich_ultrafast_2020;
hess_electroabsorption_1996; causa_femtosecond_2018; hahn_role_2016. Although
these hints are supported by time-dependent density functional theory
calculations that show the importance of delocalized excitations in C60
clusters mumthazmuhammed_impact_2022; habeebmokkath_delocalized_2021;
kobayashi_wannier-like_2020, quantitative measurements of the exciton
localization and charge separation have so far not been possible. Thus, the
C60 case highlights the need for a more direct experimental access to the
wavefunctions of the electron-hole pair excitations.

From an experimental point of view, our method of choice to access exciton
wavefunctions is time-resolved photoemission orbital tomography (tr-POT, see
Methods for the experimental realization used in this
work)puschnig_reconstruction_2009; jansen_efficient_2020;
wallauer_tracing_2021; neef_orbital-resolved_2022. In POT, the comparison with
density functional theory calculations (DFT) provides a direct connection
between photoemission data and the orbitals of the
electrons.puschnig_reconstruction_2009 The extension to the time-domain
promises valuable access also to the spatial information of excited electrons.
However, at least for organic semiconductors, it has not been explicitly
considered that photoemission of excitons requires the break-up of the two-
particle electron-hole pair and that only the photoemitted electron, but not
the hole, is directly detected. In fact, it is not clear to what extent
(tr-)POT can be reasonably used for the interpretation and analysis of such
strongly interacting correlated quasiparticles. Here we address this open
question and show how tr-POT can probe the exciton wavefunction in the example
system of a C60 multilayer.

Figure 1:  Schematic overview of time-resolved photoemission orbital
tomography of exciton states in C60. a, b, a femtosecond optical pulse (blue
pulse and blue arrow in (a) and (b), respectively) excites optically bright
excitons in a C60 film. The exciton electron-hole pairs are sketched in (a) as
correlated particles (shaded blue-red areas) with a blue sphere for the hole
and a red sphere for the electron. We probe the excitons in C60 with extreme
ultraviolet (EUV) pulses (purple pulse in (a), $h\nu=$ 26.5 eV) that break up
the electron-hole pairs and photoemit the corresponding electrons (red
spheres), of which we detect the kinetic energies and the momentum emission
patterns (yellow-green-colored disks in (a)). The optical excitation of
excitons in C60 is known to lead to the formation of a decay sequence of
singlet exciton states with varying charge-transfer
characterstadtmuller_strong_2019; emmerich_ultrafast_2020 (see (b), Si: ith
singlet excited state, S0: ground state). We are able to measure these exciton
dynamics and the corresponding orbital tomography momentum patterns by
adjusting the temporal delay between the optical excitation and the EUV probe
pulses.

We employ our recently developed setup for photoelectron momentum
microscopymedjanik_direct_2017; keunecke_time-resolved_2020; Keunecke20prb and
use ultrashort laser pulses to optically excite bright excitons in C60 thin
films that were deposited on Cu(111) (measurement temperature T $\approx$ 80
K; see Methods and Figure 1a). In the time-resolved photoemission experiment,
we detect the energy and momentum emission pattern of the photoemitted
electrons, which were initially part of the bound electron-hole pairs, i.e.
the excitons. Following the time-evolution of the photoelectron spectrum, we
can observe how the optically excited states relax to energetically lower-
lying dark exciton states with different localization and charge-transfer
characterstadtmuller_strong_2019; emmerich_ultrafast_2020 (Figure 1b and data
in Extended Fig. 6). In addition to the energy relaxation, we collect tr-POT
data, and investigate in how far these patterns can be used to access real-
space properties of the exciton wavefunctions. Specifically, we will address
two questions: which orbitals contribute to the formation of the excitons and
how this key information is imprinted in the energy- and momentum-resolved
photoemission spectra.

##  II Results & Discussion

###  II.1 The exciton spectrum of buckminsterfullerene C60

To lay the foundation for our study, we first discuss the theoretical
electronic properties of the C60 film. On top of a hybrid-functional DFT
ground state calculation, we obtain the exciton spectrum by employing the
many-body framework of $GW$ and Bethe-Salpeter-Equation ($GW$+BSE)
calculations (see Methods for full details). As shown in Fig. 2a, we model the
C60 low-temperature phase by the two symmetry-inequivalent C60 dimers 1-2 and
1-4, respectivelywang_orientational_2001, which are properly embedded to
account for polarization effects in the film. The calculated single-particle
energy levels are shown in Fig. 2b, where we group the electron removal and
electron addition energies into four bands, denoted as HOMO-1, HOMO, LUMO, and
LUMO+1 according to the orbitals of the parent orbitals of an gas-phase C60
molecule. These manifolds consist of 18, 10, 6, and 6 energy levels per dimer,
respectively, originating from the $g_{g}$+$h_{g}$, $h_{u}$, $t_{1u}$, and
$t_{1g}$ irreducible representations of the gas phase C60 orbitals
dresselhaus_c60_1996. We emphasize that the calculated $GW$ ionization levels
of HOMO and HOMO-1 of 6.7 eV and 8.1 eV are in excellent agreement with
experimental data for this C60 film (see SI and Ref. haag_signatures_2020).

Figure 2: Ab-initio calculation of the electronic structure and exciton
spectrum of C60 dimers in a crystalline multilayer sample. a, the unit cell
for a monolayer of C60, for which $GW$+BSE calculations for the dimers 1-2 and
1-4 were performed. b, electron addition/removal single-particle energies as
retrieved from the self-consistent $GW$ calculation. These energies directly
provide $\varepsilon_{v}$ in Eq. (2). c, results of the full $GW$+BSE
calculation, showing as a function of the exciton energy $\Omega$ from bottom
to top: the calculated optical absorption, the exciton band assignment S1 \-
S4, and the relative contributions to the exciton wavefunctions of different
electron-hole pair excitations $\phi_{v}(\bm{r}_{h})\chi_{c}(\bm{r}_{e})$.
Full details on the calculations are given in the Methods section. d, sketch
of the composition of the exciton wavefunction of the S1 \- S4 bands and their
expected photoemission signatures based on Eq. (2). In order to visualize the
contributing orbitals, blue holes and red electrons are assigned to the
single-particle states as shown in (b).

Building upon the $GW$ single-particle energies, we solve the Bethe-Salpeter
equation and compute the energies $\Omega_{m}$ of all correlated electron-hole
pairs (excitons). The resulting absorption spectrum (bottom panel of Fig. 2c)
agrees well with literaturewang_aggregates_1993. In addition, we obtain the
weights $X^{(m)}_{vc}$ on the specific electron-hole pairs that coherently
contribute to the $m$th exciton state, from which the exciton wavefunctions
are constructed in the Tamm-Dancoff approximation as follows:

$\psi_{m}(\bm{r}_{h},\bm{r}_{e})=\sum_{v,c}X^{(m)}_{vc}\phi_{v}^{*}(\bm{r}_{h})\chi_{c}(\bm{r}_{e}).$
(1)

This means that each exciton $\psi_{m}$ with energy $\Omega_{m}$ consists of a
weighted coherent sum of multiple electron-hole-transitions
$\phi_{v}(\bm{r}_{h})\chi_{c}(\bm{r}_{e})$ each containing one electron
orbital $\chi_{c}$ and one hole orbital $\phi_{v}$.

To gain more insight into the character of the excitons $\psi_{m}$, we
qualitatively classify them according to the most dominant orbital
contributions that are involved in the transitions. This is visualized in the
four sub-panels above the absorption spectrum in Fig. 2c. For a given exciton
energy $\Omega_{m}$, the black bars in each sub-panel show the partial
contribution $|X_{v,c}^{(m)}|^{2}$ of characteristic electron-hole transitions
$\phi_{v}(\bm{r}_{h})\chi_{c}(\bm{r}_{e})$ to a given exciton $\psi_{m}$.
Looking at individual sub-panels, we see first that characteristic electron-
hole transitions can belong to different excitons $\psi_{m}$ that have very
different exciton energies $\Omega_{m}$. For example, the blue panel in Fig.
2c shows the contributions of HOMO$\rightarrow$LUMO (abbreviated
H$\rightarrow$L) transitions as a function of exciton energy $\Omega_{m}$, and
we see that these transitions contribute to excitons that are spread in energy
over a scale of more than 1 eV (from $\Omega_{m}$ $\approx$ 1.7 eV to 3 eV).
This spread of H$\rightarrow$L contributions (and also H$-$n$\rightarrow$L+m
contributions) is caused by the fact that there are already many orbital
energies per dimer (cf. Fig. 2b) which combine to form excitons with different
degrees of localization and delocalization of the electrons and holes on one
or more molecules.

We now focus on four exciton bands of the C60 film, denoted as S1 \- S4, which
are centered around $\Omega_{S}\textsubscript{1}{}$,
$\Omega_{S}\textsubscript{2}{}$, $\Omega_{S}\textsubscript{3}{}$ and
$\Omega_{S}\textsubscript{4}{}$ at 1.9, 2.1, 2.8 and 3.6 eV, respectively (cf.
Ref. emmerich_ultrafast_2020). It is important to emphasize that each exciton
band S1 \- S4 arises from many individual excitons $\psi_{m}$ with similar
exciton energies $\Omega_{m}$ within the exciton band. Looking again at the
sub-panels, we see that the S1 and S2 exciton bands are made up of excitons
$\psi_{m}$ that are almost exclusively composed of transitions from
H$\rightarrow$L. On the other hand, the S3 shows in addition to
H$\rightarrow$L also significant contributions from H$\rightarrow$L+1
transitions (pink-dashed panel). The S4 exciton band can be characterized as
arising from H$\rightarrow$L+1 (pink-dashed panel) and H$-$1$\rightarrow$L
(orange-dash-dotted panel) as well as transitions from the HOMO to several
higher lying orbitals denoted as H$\rightarrow$L+n (yellow-dotted panel). We
emphasize that although orbitals from several different $GW$ energies
contribute, e.g., to an exciton in the S4 band, the exciton energy
$\Omega_{m}$ of each exciton $\psi_{m}$ has a single well-defined value.

###  II.2 Photoemission signatures of multiorbital contributions

In the following, we investigate whether these theoretically predicted
multiorbital characteristics of the excitons can also be probed
experimentally. As will be shown below, time-resolved photoemission
spectroscopy can indeed provide access not only to the dark exciton landscape
weinelt_dynamics_2004; dong_direct_2021; wallauer_momentum-resolved_2021;
Madeo20sci; schmitt_formation_2022, but also to the distinct orbital
contributions of exciton states. A key step in extracting this information
from the experimental data lies in a thorough comparison with simulations that
specifically consider the pump-probe photoemission process, a topic that has
recently attracted increased attention Popova2016; DeGiovannini2017;
hammon_pump-probe_2021; Reuner2023. Here, we rely on the formalism of Kern et
al.kern_exciton_tomog, which is based on a common Fermi’s golden rule approach
to photoemission Dauth2014. Assuming the exciton of Eq. (1) as the initial
state and applying the plane-wave final state approximation of POT, the
photoemission intensity of the exciton $\psi_{m}$ is formulated as

$I_{m}(E_{\mathrm{kin}},\bm{k})\propto\left|\bm{A}\bm{k}\right|^{2}\sum_{v}\left|\sum_{c}X^{(m)}_{vc}\mathcal{F}\left[\chi_{c}\right](\bm{k})\right|^{2}\times\delta\left(h\nu-
E_{\mathrm{kin}}-\varepsilon_{v}+\Omega_{m}\right).$ (2)

Here $\bm{A}$ is the vector potential of the incident light field,
$\mathcal{F}$ the Fourier transform, $\bm{k}$ the photoelectron momentum,
$h\nu$ the probe photon energy, $\varepsilon_{v}$ the $v$th ionization
potential, $\Omega_{m}$ the exciton energy, and $E_{kin}$ the energy of the
photoemitted electron. Note that $\varepsilon_{v}$ directly indicates the
final-state energy of the left-behind hole. In the context of our present
study, delving into Eq. (2) leads to two striking consequences, which we
discuss in the following.

First, we illustrate the consequences of the multiorbital character of the
exciton states on the photoelectron spectrum, and sketch in Fig. 2d the
typical single-particle energy level diagrams for the HOMO and LUMO states and
then indicate the contributing orbitals to the two-particle exciton state by
blue holes and red electrons in these states, respectively. For the S1 exciton
band (left panel), we already found that the main orbital contributions to the
band are of H$\rightarrow$L character (Fig. 2d, left, and cf. Fig. 2c, blue
panel). To determine the kinetic energy of the photoelectrons originating from
the exciton, we have to consider the correlated nature of the electron-hole
pair. The energy conservation expressed by the delta function in Eq. 2 (see
also Ref. weinelt_dynamics_2004; Zhu14jpcl) requires that the kinetic energy
of the photoelectron depends on the ionization energy of the involved HOMO
hole state $\varepsilon_{v}$ = $\varepsilon_{H}$ and the correlated electron-
hole pair energy $\Omega\approx\Omega_{S_{1}}$. Therefore, we expect to
measure a single photoelectron peak, as shown in the lower part of the left
panel of Fig. 2d. In the case of the S2 exciton the situation is similar,
since the main orbital contributions are also of H$\rightarrow$L character.
However, since the S2 exciton band has a different energy
$\Omega_{S\textsubscript{2}{}}$, the photoelectron peak is located at a
different kinetic energy with respect to the S1 peak.

In the case of the S3 exciton band, we find that in contrast to the S1 and S2
excitons not only H$\rightarrow$L, but also H$\rightarrow$L+1 transitions
contribute (Fig. 2d, middle panel, and cf. Fig 2c, blue and pink-dashed
panels, respectively). However, we still expect a single peak in the
photoemission, because the same hole states are involved for both transitions
(i.e., same $\varepsilon_{v}=\varepsilon_{H}$ in the sum in Eq. 2), and all
orbital contributions have the same exciton energy
$\approx\Omega_{S\textsubscript{3}{}}$, even though transitions with electrons
in energetically very different single-particle LUMO and LUMO+1 states
contribute. With other words, and somewhat counter-intuitively, the single-
particle energies of the electron orbitals (the LUMOs) contributing to the
exciton do not enter the energy conservation term in Eq. 2, and thus do not
affect the kinetic energy observed in the experiment.

Finally, for the S4 exciton band at $\Omega_{S}\textsubscript{4}{}=$ 3.6 eV,
we find three major contributions (Fig. 2d, right panel), where not only the
electrons but also the holes are distributed over two energetically different
levels, namely the HOMO (cf. pink-dashed and yellow-dotted panels in Fig.
2c,d) and the HOMO-1 (cf. orange-dash-dotted panels in Fig. 2c,d). Thus, there
are two different final states available for the hole, each with a different
binding energy. Consequently, the photoemission spectrum of S4 is expected to
exhibit a double-peak structure with intensity appearing $\approx$ 3.6 eV
above the HOMO kinetic energy $E_{H}$, and $\approx$ 3.6 eV above the HOMO-1
kinetic energy $E_{H-1}$, as illustrated in the right-most panel of Fig. 2d.
Relating this specifically to the single-particle picture of our $GW$
calculations, the two peaks are predicted to have a separation of
$\varepsilon_{{H-1}}$ – $\varepsilon_{{H}}$ = 8.1 – 6.7 = 1.4 eV.

In addition, Eq. (2) now also provides the theoretical framework for
interpreting momentum-resolved tr-POT data from excitons. Ground state POT can
be easily understood in terms of the Fourier transform $\mathcal{F}$ of
single-particle orbitals. A naive extension to excitons might imply an
incoherent, weighted sum of all LUMO orbitals $\chi_{c}$ contributing to the
exciton wavefunction. However, as Eq. (2) shows, such a simple picture proves
insufficient. Instead, the momentum pattern of the exciton wavefunction is
related to a coherent superposition of the electron orbitals $\chi_{c}$
weighted by the electron-hole coupling coefficients $X^{(m)}_{vc}$. The
implications of this finding are sketched in the $k_{x}$-$k_{y}$ plots in Fig.
2d and are most obvious for the S3 band. Here, the exciton is composed of
transitions with a common hole position, i.e., H$\rightarrow$L and
H$\rightarrow$L+1, leading to a coherent superposition of all 12 electron
orbitals from the LUMO and LUMO+1 in the momentum distribution. In summary,
multiple hole contributions can be identified in a multi-peak structure in the
photoemission spectrum, and multiple electron contributions will result in a
coherent sum of the electron orbitals that can be identified in the
corresponding energy-momentum patterns from tr-POT data.

###  II.3 Disentangling multiorbital contributions experimentally

These very strong predictions about multi-peaked photoemission spectra due to
multiorbital excitons can be directly verified in an experiment on C60 by
comparing spectra for resonant excitation of either the S3 or the S4 excitons
(cf. Fig. 2). The corresponding experimental data are shown in Fig. 3a and 3c,
respectively. Starting from the excitation of the S3 exciton band with $h\nu=$
2.9 eV photon energy (which is sufficiently resonant to excite the manifold of
exciton states that make up the S3 band around
$\Omega_{S\textsubscript{3}{}}=$ 2.8 eV), we can clearly identify the direct
excitation (at 0 fs delay) of the exciton S3 feature at an energy of E
$\approx$ 2.8 eV above the kinetic energy $E_{H}$ of the HOMO level. Shortly
after the excitation, additional photoemission intensity builds up at
$E-E_{H}\approx$ 2.0 eV and $\approx$ 1.7 eV, which is known to be caused by
relaxation to the S2 and S1 dark exciton statesemmerich_ultrafast_2020 and is
in good agreement with the theoretically predicted energies of $E-E_{H}$
$\approx$ 2.1 eV and $\approx$1.9 eV.

Figure 3: a-d, comparison of the time-resolved photoelectron spectra of
multilayer C60 for $h\nu=$ 2.9 eV excitation and $h\nu=$ 3.6 eV excitation
((a) and (c), respectively), both normalized and shifted in time to match the
intensity of the S3 signals (see Methods for full details on the data
analysis). As can be seen in the difference (b), for $h\nu=$ 3.6 eV pump we
observe an enhancement of the photoemission yield around $E-E_{H}\approx$ 3.6
eV as well as around $E-E_{H}\approx$ 2.2 eV. We attribute this signal to the
S4 exciton band, which has hole contributions stemming from both the HOMO and
the HOMO-1. To further quantify the signal of the S4 exciton, (d) shows energy
distribution curves for both measurements at early delays, showing the
enhancement in the $h\nu=$ 3.6 eV measurement.

Changing now the pump photon energy to $h\nu=$ 3.6 eV for direct excitation of
the S4 exciton band (Fig. 3c), two distinct peaks at $\approx$ 3.6 eV above
the HOMO and the HOMO-1 are expected from theory. While photoemission
intensity at $E-E_{H}\approx$ 3.6 eV above the HOMO level is readily visible
in Fig. 3c, the second feature at 3.6 eV above the HOMO-1 is expected at
$E-E_{H}\approx$ 2.2 eV above the HOMO level (corresponding to $E-E_{H-1}$
$\approx$ 3.6 eV) and thus almost degenerate with the aforementioned S2 dark
exciton band at about $E-E_{H}\approx$ 2.0 eV, which appears after the optical
excitation due to relaxation processes. Therefore, we need to pinpoint this
second H$-$1$\rightarrow$L contribution to the S4 exciton at the earliest time
of the excitation. Indeed, a closer look around 0 fs delay shows additional
photoemission intensity at about $E-E_{H}\approx$ 2.2 eV. Using difference
maps (Fig. 3b) and direct comparisons of energy-distribution-curves at
selected time-steps (Fig. 3d), we clearly find a double-peak structure
corresponding to the energy difference of $\approx$1.4 eV of the HOMO and
HOMO-1 levels. Thereby, we have shown that photoelectron spectroscopy, in
contrast to other techniques (e.g., absorption spectroscopy), is indeed able
to disentangle different orbital contributions of the excitons. In this way,
we have validated the theoretically predicted multi-peak structure of the
multiorbital exciton state that is implied by Eq. (2). We also see that the
photoelectron energies in the spectrum turn out to be sensitive probes of the
corresponding hole contributions of the correlated exciton states.

We note that the signature of the S3 excitons, even if not directly excited
with the light pulse in this measurement, is still visible and moreover with
significantly higher intensity than the multiorbital signals of the resonantly
excited S4 exciton band. This observation strongly suggests that there is a
very fast relaxation from the S4 exciton to the S3 exciton, with relaxation
times well below 50 fs (see Extended Fig. 6).

###  II.4 Time-resolved photoemission orbital tomography of exciton
wavefunctions

Based on the excellent agreement between the experiment and the $GW$+BSE
theoretical results, we are now ready to investigate to what extent the
momentum patterns from tr-POT data of excitons in organic semiconductors
contain information about the real-space spatial distribution of the exciton
wavefunction. In the experiment, we once again excite the S3 exciton band in
the C60 film with $h\nu=$ 2.9 eV pump energy, and we now use femtosecond tr-
POT to collect the momentum fingerprints of the directly excited S3 excitons
around 0 fs and the subsequently built-up dark S2 and S1 excitons that appear
in the exciton relaxation cascade in the C60 film (see Fig. 4a-c, where the
momentum map of the lowest energy S1 exciton band is plotted in (a), the S2 in
(b) and the highest energy S3 exciton band in (c); see Extended Fig. 6 for
time-resolved traces of the exciton formation and relaxation dynamics). We
note that the collection of these data required integration times of up to 70
hours, and that a measurement of the comparatively low-intensity S4 feature
when excited with $h\nu=$ 3.6 eV has not yet proved feasible. For the
interpretation of the collected POT momentum maps from the S1, S2, and S3
excitons, we also calculate the expected momentum fingerprints for the
wavefunctions obtained from the $GW$+BSE calculation for both dimers, each
rotated to all occurring orientations in the crystal. Finally, for the
theoretical momentum maps, we sum up the photoelectron intensities of each
electron-hole transition in an energy range of 200 meV centered on the exciton
band. The results are shown in Fig. 4d-f below the experimental data for
direct comparison.

Figure 4: a-f, comparison of the (a-c) experimental momentum maps acquired for
the three exciton bands observed in C60 with the (d-f) predicted momentum maps
retrieved from $GW$+BSE. Note that the center of the experimental maps could
not be analyzed due to a space-charge-induced background signal in this region
(gray area, see Methods). g-i, isosurfaces of the integrated electron
probability density (yellow) within the 1-2 dimer for fixed hole positions on
the bottom-left molecule (blue circle) of the dimer for the (g) S1, the (h)
S2, and the (i) S3 exciton bands.

First, we observe that the experimental momentum maps of the S1 and S2 states
are largely similar (Fig. 4a,b), showing six lobes centered at
$k_{\parallel}\approx$ 1.2 Å-1. These features, as well as the energy
splitting between S1 and S2 (cf. Fig. 2c), are accurately reproduced by the
$GW$+BSE prediction (Fig. 4d,e). Furthermore, also the $GW$+BSE calculation
shows very similar momentum maps for S1 and S2, suggesting a similar spatial
structure of the excitons. This is in contrast to a naive application of
static POT to the unoccupied orbitals of the DFT ground state of C60, which
does show a similar momentum map for the LUMO, but cannot explain a kinetic
energy difference in the photoemission signal, nor give any indication of
differences in the corresponding exciton wavefunctions. With this agreement
between experiment and theory, we now extract the spatial properties of the
$GW$+BSE exciton wavefunctions. To visualize the degree of charge-transfer of
these two-particle exciton wavefunctions $\psi_{m}(\bm{r_{h}},\bm{r_{e}})$, we
integrate the electron probability density over all possible hole positions
$\bm{r_{h}}$, considering only hole positions at one of the C60 molecules in
the dimer. This effectively fixes the hole contribution to a particular C60
molecule (blue circles in Fig. 4g-i indicate the boundary of considered hole
positions around one molecule, hole distribution not shown), and provides a
probability density for the electronic part of the exciton wavefunction in the
dimer, which we visualize by a yellow isosurface (see Fig. 4g-i). Obviously,
in the case of S1 and S2 (Fig. 4g,h), when the hole position is restricted to
one molecule of the dimer, the electronic part of the exciton wavefunction is
localized at the same molecule of the dimer. Our calculations thus suggest
that the S1 and S2 excitons are of Frenkel-like nature. Their energy
difference originates from different excitation symmetries possible for the
H$\rightarrow$L transition (namely $t_{1g}$, $t_{2g}$, and $g_{g}$ for the S1
and $h_{g}$ for the S2).kobayashi_wannier-like_2020

In contrast to the S1 and S2 excitons, the momentum map of the S3 band shows a
much more star-shaped POT fingerprint in both theory and experiment (Fig.
4c,f). This is to be expected, since the electronic part of the S3 excitons
contains not only contributions from the LUMO orbital, but also contributions
from the LUMO+1 orbital. Note, however, that we find the experimentally
observed star-shaped pattern to be only partially reproduced by the $GW$+BSE
calculation. An indication towards the cause of this discrepancy is found by
considering the electron-hole separation of the excitons making up the S3
band. Here, we find that the positions of the electron and the hole
contributions are strongly anticorrelated (Fig. 4i), with the electron
confined to the neighboring molecule of the dimer. In fact, the mean electron-
hole separation is as large as 7.6 Å, which is close to the core-to-core
distance of the C60 molecules. Although these theoretical results confirm the
previously-reported charge-transfer nature of the S3 excitons
stadtmuller_strong_2019; emmerich_ultrafast_2020, they also reflect the
limitations of the C60 dimer approach. Indeed, the dimer represents the
minimal model to account for an intermolecular exciton delocalization effect,
but it cannot fully account for dispersion effectshaag_signatures_2020 (cf.
Extended Fig. 5), which are required for a quantitative comparison with
experimental data. Besides the discrepancy in the S3 momentum map, this could
also be an explanation why the S2 in the present work is of Frenkel-like
nature, but could have charge-transfer character according to previous
studiesstadtmuller_strong_2019; emmerich_ultrafast_2020. However, future
developments will certainly allow scaling up of the cluster size in the
calculation, so that exciton wavefunctions with larger electron-hole
separation can be accurately described. Most importantly, we find that the
present dimer $GW$+BSE calculations are clearly suited to elucidate the
multiorbital character of the excitons, which is an indispensable prerequisite
for the correct interpretation of tr-POT data of excitons in organic
semiconductors.

##  III CONCLUSION

In conclusion, we have shown how the energy- and momentum-resolved
photoemission spectrum of excitons in an organic semiconductor depends on the
multiorbital nature of these excitons. By extending POT to fully-interacting
exciton states calculated in the framework of the Bethe-Salpeter equation, we
found that the energy of the photoemitted electron of the exciton
quasiparticle is determined by the position of the hole and the exciton energy
in combination with the probe photon energy. This leads to the prediction of
multiple peaks in the photoelectron spectrum, which we verify experimentally,
and allows disentangling the different orbital contributions, the wavefunction
localization, and the charge-transfer character. Similarly, the momentum
fingerprint provides access to the electron states that make up the exciton.
Most importantly, we introduce time-resolved photoemission orbital tomography
as a key technique for the study of exciton wavefunctions in organic
semiconductors.

##  IV Acknowledgements

This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) - 432680300/SFB 1456, project B01 and 217133147/SFB 1073,
projects B07 and B10. G.S.M.J. acknowledges financial support by the Alexander
von Humboldt Foundation. A.W., C.S.K., and P.P acknowledge support from the
Austrian Science Fund (FWF) project I 4145. The computational results
presented were achieved using the Vienna Scientific Cluster (VSC) and the
local high-performance resources of the University of Graz. R.H., M.A., and
B.S. acknowledge financial support by the DFG - 268565370/TRR 173, projects
B05 and A02. B.S. acknowledges further support by the Dynamics and Topology
Center funded by the State of Rhineland-Palatinate.

##  V Author Contributions

D.St., M.R., S.S., M.A., B.S., P.P., G.S.M.J. and S.M. conceived the research.
W.B., D.Sch. and J.P.B. carried out the time-resolved momentum microscopy
experiments. W.B. analyzed the data. W.B. and R.H. prepared the samples. A.W.,
C.S.K., G.D.A, X.B. and P.P. performed the calculations and analyzed the
theoretical results. All authors discussed the results. G.S.M.J. and S.M. were
responsible for the overall project direction and wrote the manuscript with
contributions from all co-authors.

## References

  * Ke  _et al._ (2015) Weijun Ke, Dewei Zhao, Corey R. Grice, Alexander J. Cimaroli, Jie Ge, Hong Tao, Hongwei Lei, Guojia Fang, and Yanfa Yan, “Efficient planar perovskite solar cells using room-temperature vacuum-processed C60 electron selective layers,” Journal of Materials Chemistry A 3, 17971–17976 (2015), publisher: The Royal Society of Chemistry.
  * Puente Santiago  _et al._ (2020) Alain R. Puente Santiago, Tianwei He, Oscar Eraso, Md Ariful Ahsan, Aruna N. Nair, Venkata S. N. Chava, Ting Zheng, Srikanth Pilla, Olivia Fernandez-Delgado, Aijun Du, Sreeprasad T. Sreenivasan, and Luis Echegoyen, “Tailoring the Interfacial Interactions of van der Waals 1T-MoS2/C60 Heterostructures for High-Performance Hydrogen Evolution Reaction Electrocatalysis,” Journal of the American Chemical Society 142, 17923–17927 (2020), publisher: American Chemical Society.
  * Yu  _et al._ (2020) Zhenhua Yu, Zhibin Yang, Zhenyi Ni, Yuchuan Shao, Bo Chen, Yuze Lin, Haotong Wei, Zhengshan J. Yu, Zachary Holman, and Jinsong Huang, “Simplified interconnection structure based on C60/SnO2-x for all-perovskite tandem solar cells,” Nature Energy 5, 657–665 (2020), number: 9 Publisher: Nature Publishing Group.
  * Wang  _et al._ (1993) Yi Ming Wang, Prashant V. Kamat, and L. K. Patterson, “Aggregates of fullerene C60 and C70 formed at the gas-water interface and in DMSO/water mixed solvents. A spectral study,” J. Phys. Chem. 97, 8793–8797 (1993).
  * Stadtmüller  _et al._ (2019) Benjamin Stadtmüller, Sebastian Emmerich, Dominik Jungkenn, Norman Haag, Markus Rollinger, Steffen Eich, Mahalingam Maniraj, Martin Aeschlimann, Mirko Cinchetti, and Stefan Mathias, “Strong modification of the transport level alignment in organic materials after optical excitation,” Nature Communications 10, 1470 (2019), number: 1 Publisher: Nature Publishing Group.
  * Emmerich  _et al._ (2020) Sebastian Emmerich, Sebastian Hedwig, Benito Arnoldi, Johannes Stöckl, Florian Haag, Ralf Hemm, Mirko Cinchetti, Stefan Mathias, Benjamin Stadtmüller, and Martin Aeschlimann, “Ultrafast Charge-Transfer Exciton Dynamics in C60 Thin Films,” The Journal of Physical Chemistry C 124, 23579–23587 (2020), publisher: American Chemical Society.
  * Hess  _et al._ (1996) Bret C. Hess, Don V. Bowersox, Shant H. Mardirosian, and Lawrence D. Unterberger, “Electroabsorption in C60 and C70. Third-order nonlinearity in molecular and solid states,” Chemical Physics Letters 248, 141–146 (1996).
  * Causa’  _et al._ (2018) Martina Causa’, Ivan Ramirez, Josue F. Martinez Hardigree, Moritz Riede, and Natalie Banerji, “Femtosecond Dynamics of Photoexcited C60 Films,” J. Phys. Chem. Lett. 9, 1885–1892 (2018), publisher: American Chemical Society.
  * Hahn  _et al._ (2016) Tobias Hahn, Steffen Tscheuschner, Christina Saller, Peter Strohriegl, Puttaraju Boregowda, Tushita Mukhopadhyay, Satish Patil, Dieter Neher, Heinz Bässler, and Anna Köhler, “Role of Intrinsic Photogeneration in Single Layer and Bilayer Solar Cells with C60 and PCBM,” J. Phys. Chem. C 120, 25083–25091 (2016), publisher: American Chemical Society.
  * Mumthaz Muhammed  _et al._ (2022) Mufasila Mumthaz Muhammed, Junais Habeeb Mokkath, and Ali J. Chamkha, “Impact of packing arrangement on the optical properties of C60 cluster aggregates,” Physical Chemistry Chemical Physics 24, 5946–5955 (2022), publisher: Royal Society of Chemistry.
  * Habeeb Mokkath (2021) Junais Habeeb Mokkath, “Delocalized exciton formation in C60 linear molecular aggregates,” Physical Chemistry Chemical Physics 23, 21901–21912 (2021), publisher: Royal Society of Chemistry.
  * Kobayashi  _et al._ (2020) Hajime Kobayashi, Shinnosuke Hattori, Raku Shirasawa, and Shigetaka Tomiya, “Wannier-Like Delocalized Exciton Generation in C60 Fullerene Clusters: A Density Functional Theory Study,” The Journal of Physical Chemistry C 124, 2379–2387 (2020), publisher: American Chemical Society.
  * Puschnig  _et al._ (2009) Peter Puschnig, Stephen Berkebile, Alexander J. Fleming, Georg Koller, Konstantin Emtsev, Thomas Seyller, John D. Riley, Claudia Ambrosch-Draxl, Falko P. Netzer, and Michael G. Ramsey, “Reconstruction of Molecular Orbital Densities from Photoemission Data,” Science 326, 702–706 (2009).
  * Jansen  _et al._ (2020) Matthijs Jansen, Marius Keunecke, Marten Düvel, Christina Möller, David Schmitt, Wiebke Bennecke, Jasmin Kappert, Daniel Steil, D. Russell Luke, Sabine Steil, and Stefan Mathias, “Efficient orbital imaging based on ultrafast momentum microscopy and sparsity-driven phase retrieval,” New Journal of Physics (2020), 10.1088/1367-2630/ab8aae.
  * Wallauer  _et al._ (2021a) R. Wallauer, M. Raths, K. Stallberg, L. Münster, D. Brandstetter, X. Yang, J. Güdde, P. Puschnig, S. Soubatch, C. Kumpf, F. C. Bocquet, F. S. Tautz, and U. Höfer, “Tracing orbital images on ultrafast time scales,” Science 371, 1056–1059 (2021a), publisher: American Association for the Advancement of Science.
  * Neef  _et al._ (2022) Alexander Neef, Samuel Beaulieu, Sebastian Hammer, Shuo Dong, Julian Maklar, Tommaso Pincelli, R. Patrick Xian, Martin Wolf, Laurenz Rettig, Jens Pflaum, and Ralph Ernstorfer, “Orbital-resolved Observation of Singlet Fission,” arXiv:2204.06824 [cond-mat] (2022), arXiv: 2204.06824.
  * Medjanik  _et al._ (2017) K. Medjanik, O. Fedchenko, S. Chernov, D. Kutnyakhov, M. Ellguth, A. Oelsner, B. Schönhense, T. R. F. Peixoto, P. Lutz, C.-H. Min, F. Reinert, S. Däster, Y. Acremann, J. Viefhaus, W. Wurth, H. J. Elmers, and G. Schönhense, “Direct 3D mapping of the Fermi surface and Fermi velocity,” Nature Materials 16, 615–621 (2017).
  * Keunecke  _et al._ (2020a) Marius Keunecke, Christina Möller, David Schmitt, Hendrik Nolte, G. S. Matthijs Jansen, Marcel Reutzel, Marie Gutberlet, Gyula Halasi, Daniel Steil, Sabine Steil, and Stefan Mathias, “Time-resolved momentum microscopy with a 1 MHz high-harmonic extreme ultraviolet beamline,” Review of Scientific Instruments 91, 063905 (2020a), publisher: American Institute of Physics.
  * Keunecke  _et al._ (2020b) Marius Keunecke, Marcel Reutzel, David Schmitt, Alexander Osterkorn, Tridev A. Mishra, Christina Möller, Wiebke Bennecke, G. S. Matthijs Jansen, Daniel Steil, Salvatore R. Manmana, Sabine Steil, Stefan Kehrein, and Stefan Mathias, “Electromagnetic dressing of the electron energy spectrum of Au(111) at high momenta,” Physical Review B 102, 161403 (2020b).
  * Wang  _et al._ (2001) Haiqian Wang, Changgan Zeng, Bing Wang, J. G. Hou, Qunxiang Li, and Jinlong Yang, “Orientational configurations of the C 60 molecules in the ( 2 × 2 ) superlattice on a solid C 60 (111) surface at low temperature,” Physical Review B 63, 085417 (2001).
  * Dresselhaus  _et al._ (1996) M.S. Dresselhaus, G. Dresselhaus, and P.C. Eklund, “Chapter 12 - electronic structure,” in  _Science of Fullerenes and Carbon Nanotubes_, edited by M.S. Dresselhaus, G. Dresselhaus, and P.C. Eklund (Academic Press, San Diego, 1996) pp. 413–463.
  * Haag  _et al._ (2020) Norman Haag, Daniel Lüftner, Florian Haag, Johannes Seidel, Leah L. Kelly, Giovanni Zamborlini, Matteo Jugovac, Vitaliy Feyer, Martin Aeschlimann, Peter Puschnig, Mirko Cinchetti, and Benjamin Stadtmüller, “Signatures of an atomic crystal in the band structure of a ${\mathrm{C}}_{60}$ thin film,” Phys. Rev. B 101, 165422 (2020), publisher: American Physical Society.
  * Weinelt  _et al._ (2004) Martin Weinelt, Michael Kutschera, Thomas Fauster, and Michael Rohlfing, “Dynamics of Exciton Formation at the Si(100) c(4 × 2) Surface,” Phys. Rev. Lett. 92, 126801 (2004), publisher: American Physical Society.
  * Dong  _et al._ (2021) Shuo Dong, Michele Puppin, Tommaso Pincelli, Samuel Beaulieu, Dominik Christiansen, Hannes Hübener, Christopher W. Nicholson, Rui Patrick Xian, Maciej Dendzik, Yunpei Deng, Yoav William Windsor, Malte Selig, Ermin Malic, Angel Rubio, Andreas Knorr, Martin Wolf, Laurenz Rettig, and Ralph Ernstorfer, “Direct measurement of key exciton properties: Energy, dynamics, and spatial distribution of the wave function,” Natural Sciences 1, e10010 (2021), type: Journal Article.
  * Wallauer  _et al._ (2021b) Robert Wallauer, Raul Perea-Causin, Lasse Münster, Sarah Zajusch, Samuel Brem, Jens Güdde, Katsumi Tanimura, Kai-Qiang Lin, Rupert Huber, Ermin Malic, and Ulrich Höfer, “Momentum-Resolved Observation of Exciton Formation Dynamics in Monolayer WS2,” Nano Lett. 21, 5867–5873 (2021b), publisher: American Chemical Society.
  * Madéo  _et al._ (2020) Julien Madéo, Michael K. L. Man, Chakradhar Sahoo, Marshall Campbell, Vivek Pareek, E. Laine Wong, Abdullah Al-Mahboob, Nicholas S. Chan, Arka Karmakar, Bala Murali Krishna Mariserla, Xiaoqin Li, Tony F. Heinz, Ting Cao, and Keshav M. Dani, “Directly visualizing the momentum-forbidden dark excitons and their dynamics in atomically thin semiconductors,” Science 370, 1199–1204 (2020).
  * Schmitt  _et al._ (2022) David Schmitt, Jan Philipp Bange, Wiebke Bennecke, AbdulAziz AlMutairi, Giuseppe Meneghini, Kenji Watanabe, Takashi Taniguchi, Daniel Steil, D. Russell Luke, R. Thomas Weitz, Sabine Steil, G. S. Matthijs Jansen, Samuel Brem, Ermin Malic, Stephan Hofmann, Marcel Reutzel, and Stefan Mathias, “Formation of moiré interlayer excitons in space and time,” Nature 608, 499–503 (2022).
  * Popova-Gorelova  _et al._ (2016) Daria Popova-Gorelova, Jochen Küpper, and Robin Santra, “Imaging electron dynamics with time- and angle-resolved photoelectron spectroscopy,” Phys. Rev. A 94, 013412 (2016).
  * De Giovannini  _et al._ (2017) Umberto De Giovannini, Hannes Hübener, and Angel Rubio, “A First-Principles Time-Dependent Density Functional Theory Framework for Spin and Time-Resolved Angular-Resolved Photoelectron Spectroscopy in Periodic Systems,” Journal of Chemical Theory and Computation 13, 265–273 (2017).
  * Hammon and Kümmel (2021) Sebastian Hammon and Stephan Kümmel, “Pump-probe photoemission simulated in real time: Revealing many-particle signatures,” Physical Review A 104, 012815 (2021).
  * Reuner and Popova-Gorelova (2023) Marvin Reuner and Daria Popova-Gorelova, “Attosecond imaging of photoinduced dynamics in molecules using time-resolved photoelectron momentum microscopy,” Phys. Rev. A 107, 023101 (2023).
  * (32) Christian Kern, Andreas Windischbacher, and Peter Puschnig, “Photoemission orbital tomography of excitons,” unpublished.
  * Dauth  _et al._ (2014) M. Dauth, M. Wiessner, V. Feyer, A. Schöll, P. Puschnig, F. Reinert, and S. Kümmel, “Angle resolved photoemission from organic semiconductors: orbital imaging beyond the molecular orbital interpretation,” New J. Phys. 16, 103005 (2014).
  * Zhu (2014) X. Y. Zhu, “How to draw energy level diagrams in excitonic solar cells,” The Journal of Physical Chemistry Letters , 2283–2288 (2014).
  * Schönhense  _et al._ (2018) B. Schönhense, K. Medjanik, O. Fedchenko, S. Chernov, M. Ellguth, D. Vasilyev, A. Oelsner, J. Viefhaus, D. Kutnyakhov, W. Wurth, H. J. Elmers, and G. Schönhense, “Multidimensional photoemission spectroscopy—the space-charge limit,” New Journal of Physics 20, 033004 (2018), publisher: IOP Publishing.
  * David  _et al._ (1991) Richard M. David, William I. F.and Ibberson, Judy C. Matthewman, Kosmas Prassides, T. John S. Dennis, Jonathan P. Hare, Harold W. Kroto, Roger Taylor, and David R. M. Walton, “Crystal structure and bonding of ordered c60,” Nature 353, 147–149 (1991).
  * David  _et al._ (1992) W. I. F. David, R. M. Ibberson, T. J. S. Dennis, J. P. Hare, and K. Prassides, “Structural phase transitions in the fullerene c60,” Europhysics Letters 18, 219 (1992).
  * Krishnan  _et al._ (1980) R. Krishnan, J. S. Binkley, R. Seeger, and J. A. Pople, “Self‐consistent molecular orbital methods. xx. a basis set for correlated wave functions,” The Journal of Chemical Physics 72, 650–654 (1980), https://doi.org/10.1063/1.438955 .
  * Frisch  _et al._ (1984) Michael J. Frisch, John A. Pople, and J. Stephen Binkley, “Self‐consistent molecular orbital methods 25. supplementary functions for gaussian basis sets,” The Journal of Chemical Physics 80, 3265–3269 (1984), https://doi.org/10.1063/1.447079 .
  * Adamo and Barone (1999) Carlo Adamo and Vincenzo Barone, “Toward reliable density functional methods without adjustable parameters: The pbe0 model,” The Journal of Chemical Physics 110, 6158–6170 (1999), https://doi.org/10.1063/1.478522 .
  * Grimme  _et al._ (2010) Stefan Grimme, Jens Antony, Stephan Ehrlich, and Helge Krieg, “A consistent and accurate ab initio parametrization of density functional dispersion correction (dft-d) for the 94 elements h-pu,” The Journal of Chemical Physics 132, 154104 (2010), https://doi.org/10.1063/1.3382344 .
  * Neese (2012) Frank Neese, “The ORCA program system,” Wiley Interdisciplinary Reviews: Computational Molecular Science 2, 73–78 (2012).
  * Neese (2022) Frank Neese, “Software update: The orca program system—version 5.0,” WIREs Computational Molecular Science 12, e1606 (2022), https://wires.onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1606 .
  * Jacquemin  _et al._ (2017) Denis Jacquemin, Ivan Duchemin, and Xavier Blase, “Is the Bethe–Salpeter Formalism Accurate for Excitation Energies? Comparisons with TD-DFT, CASPT2, and EOM-CCSD,” J. Phys. Chem. Lett. 8, 1524–1529 (2017).
  * D’Avino  _et al._ (2014) Gabriele D’Avino, Luca Muccioli, Claudio Zannoni, David Beljonne, and Zoltán G. Soos, “Electronic polarization in organic crystals: A comparative study of induced dipoles and intramolecular charge redistribution schemes,” Journal of Chemical Theory and Computation 10, 4959–4971 (2014), pMID: 26584380, https://doi.org/10.1021/ct500618w .
  * Li  _et al._ (2016) Jing Li, Gabriele D’Avino, Ivan Duchemin, David Beljonne, and Xavier Blase, “Combining the many-body gw formalism with classical polarizable models: Insights on the electronic structure of molecular solids,” The Journal of Physical Chemistry Letters 7, 2814–2820 (2016), pMID: 27388926, https://doi.org/10.1021/acs.jpclett.6b01302 .
  * Li  _et al._ (2017) Jing Li, Gabriele D’Avino, Anton Pershin, Denis Jacquemin, Ivan Duchemin, David Beljonne, and Xavier Blase, “Correlated electron-hole mechanism for molecular doping in organic semiconductors,” Phys. Rev. Materials 1, 025602 (2017).
  * Dai  _et al._ (1994) Sheng Dai, L. Mac Toth, G. D. Del Cul, and David H. Metcalf, “Ultraviolet‐visible absorption spectrum of C60 vapor and determination of the C60 vaporization enthalpy,” J. Chem. Phys. 101, 4470–4471 (1994), publisher: American Institute of Physics.
  * Hasegawa  _et al._ (1998) Shinji Hasegawa, Takayuki Miyamae, Kyuya Yakushi, Hiroo Inokuchi, Kazuhiko Seki, and Nobuo Ueno, “Origin of the photoemission intensity oscillation of ${\mathrm{c}}_{60}$,” Phys. Rev. B 58, 4927–4933 (1998).

##  VI Methods

###  VI.1 Femtosecond momentum microscopy of C60/Cu(111)

We apply full multidimensional time- and angle-resolved photoelectron
spectroscopy (tr-ARPES) to a multilayer C60 crystal evaporated onto Cu(111),
where the film thickness was such that no photoemission signature of the
underlying Cu(111) could be observed in our experiment. We verified the sample
quality by performing momentum microscopy of the occupied HOMO and HOMO-1
states simultaneously to the measurement of the excited states (see Extended
Fig. 5). Femtosecond exciton dynamics were induced using $\approx$100 fs,
$h\nu=$ 2.9 eV or $\approx$100 fs, $h\nu=$ 3.6 eV laser pulses derived from
the frequency-doubled output of a optical parametric amplifier. The exciton
dynamics were probed using our custom photoemission momentum microscope with a
500 kHz ultrafast 26.5 eV extreme ultraviolet (EUV) light source
keunecke_time-resolved_2020 that enables us to map the photoelectron momentum
distribution over the full photoemission horizon in a kinetic energy range
exceeding 6 eV and an overall time resolution of $\approx$100 fs (the EUV
pulse length is about 20 fs). The pump fluence was set to $90(10)$ $\mu$J/cm2
and $20(5)$ $\mu$J/cm2 for the $h\nu=$ 2.9 eV and the $h\nu=$ 3.6 eV
measurement, respectively. To prevent the free rotation of C60 molecules, we
cooled the sample down to $\approx$80 K wang_orientational_2001. In addition
to the resulting long-range periodic ordering of the C60 crystal, cooling was
also observed to prevent light-induced polymerization.

###  VI.2 Momentum microscopy data analysis

####  VI.2.1 Space-charge effects

In the momentum microscopy experiment, a balance has to be found between
sufficiently low pump and probe light intensities to avoid space-charge
effects, but also having sufficient intensity for the optical excitation
(pump) and reasonably short integration times (probe). Most of the time, for
our settings, small space-charge effects are present in the data, but do not
lead to strong distortions in the band-structure data and can be easily
corrected. Therefore, the first step in the data analysis was to subtract a
space-charge-induced delay- and momentum-dependent kinetic energy shift. For
this purpose, the central kinetic energy of the HOMO was determined for
normalization. To avoid the influence of the C60 crystal band structure
haag_signatures_2020, we fitted a two-dimensional Lorentzian
schonhense_multidimensional_2018 and shifted the kinetic energy distribution
accordingly, leading to the expected overall flat shape of the molecular
orbitals in the ARPES data.

####  VI.2.2 Replicas from the 13th harmonic

Although we use narrow-band multilayer mirrors to select the 11th harmonic at
$h\nu$ = 26.5 eV from our laser-based high-harmonic generation
spectrumkeunecke_time-resolved_2020, we observe subtle replicas of the HOMO,
HOMO-1, and HOMO-2 states in the unoccupied regime of the spectrum that are
caused by photoemission from the 13th harmonic at at $h\nu$ = 31.2 eV. To
quantify these replica signals, we fitted the static reference spectrum above
$E-E_{H}=1.2$ eV (i.e., in the unoccupied regime of the spectrum) with three
Gaussian-shaped peaks for the HOMO replicas and an exponential function to
account for residual photoemission intensity in the unoccupied regime that is
caused in this spectral region by the much stronger direct 11th harmonic one-
photon-photoemission from the HOMO state. After carrying out this fitting
routine, we are able to calculate clean 11th harmonic spectra (static and
time-resolved) via subtraction of the fitted 13th harmonic HOMO replicas. Note
that we only subtract the replica signals, but not the background signal that
is caused by one-photon-photoemission with the 11th harmonic from the HOMO
state, because this background is time-dependentstadtmuller_strong_2019, and
needs to be explicitly considered in the fitting procedure. The data shown in
Fig. 3 of the main text is processed in the way described above.

####  VI.2.3 Fitting procedure for the time-resolved data

From replica-free trPES data for $h\nu=$ 2.9 eV excitation, we determine the
amplitude $A_{i}$, kinetic energy $E_{i}$, and bandwidth $\Delta E_{i}$ for
the $i$th exciton signature using a global fitting approach. In particular, we
apply the model

$\begin{split}I(E,t)&=\frac{A_{S\textsubscript{3}{}}(t)}{\sqrt{2\pi}\Delta
E_{S\textsubscript{3}{}}}\mathrm{exp}\left[(E-E_{S\textsubscript{3}{}}(t))^{2}/\Delta
E_{S\textsubscript{3}{}}^{2}\right]\\\
&+\frac{A_{S\textsubscript{2}{}}(t)}{\sqrt{2\pi}\Delta
E_{S\textsubscript{2}{}}}\mathrm{exp}\left[(E-E_{S\textsubscript{2}{}})^{2}/\Delta
E_{S\textsubscript{2}{}}^{2}\right]\\\
&+\frac{A_{S\textsubscript{1}{}}(t)}{\sqrt{2\pi}\Delta
E_{S\textsubscript{1}{}}}\mathrm{exp}\left[(E-E_{S\textsubscript{1}{}})^{2}/\Delta
E_{S\textsubscript{1}{}}^{2}\right]\\\
&+A_{bg}(t)*\mathrm{exp}\left[-E/\tau\right].\end{split}$ (3)

Here, the last term is needed to account for the above-mentioned delay-
dependent photoemission intensity that is caused by a transient
renormalization of the HOMO state, as found in Ref. stadtmuller_strong_2019.

The fit results of this model applied to the $h\nu=$ 2.9 eV excitation and
momentum-integrated data are shown in Table 1, and Extended Fig. 6a for the
time-resolved exciton dynamics.

Exciton | Kin. Energy (eV) | Bandwidth (FWHM) (eV)
---|---|---
S3 (t=0 fs) | 2.768(2) | 0.606(4)
S2 | 1.978(2) | 0.406(3)
S1 | 1.667(1) | 0.362(2)
Table 1:  Peak parameters for $h\nu=$ 2.9 eV excitation, extracted using Eq.
3. Note that we give the full width at half maximum (FWHM) for the bandwidth.

For the measurement with $h\nu=$ 3.6 eV excitation, we account for the S4
exciton band by extending the model in Eq. 3 with a set of Gaussian peaks with
identical temporal evolution, given by

$\begin{split}I(E,t)=...+\frac{A_{S\textsubscript{4}{}}(t)}{2\sqrt{2\pi}\Delta
E_{S\textsubscript{4}{}}}(&\mathrm{exp}\left[(E-E_{S\textsubscript{4}{},upper})^{2}/\Delta
E_{S\textsubscript{4}{}}^{2}\right]\\\
+&\mathrm{exp}\left[(E-E_{S\textsubscript{4}{},lower})^{2}/\Delta
E_{S\textsubscript{4}{}}^{2}\right]),\end{split}$ (4)

which follows the same notation as Eq. 3. Here, we set
$E_{S\textsubscript{4}{},upper}$ to be close to 3.6 eV, and following the
$GW$+BSE calculation we set
$E_{S\textsubscript{4}{},upper}-E_{S\textsubscript{4}{},lower}=1.4$ eV.
Fitting this model to the momentum-integrated $h\nu=$ 3.6 eV excitation data,
we find $E_{S\textsubscript{4}{},upper}=3.59(1)$ eV, and for the FWHM of the
S4 we find $0.58(3)$ eV. The time-resolved amplitudes retrieved using this
model are shown in Extended Fig. 6b. Furthermore, this analysis was used in
Fig. 3 of the main text to subtract the exponential background
$A_{bg}(t)*\mathrm{exp}\left[-E/\tau\right]$ related to the transient
broadening of the HOMO state.

####  VI.2.4 Fitting procedure for the time- and momentum-resolved data

In order to analyze the time-resolved data also momentum-resolved and thereby
retrieve the momentum patterns that are shown in Fig. 1 and Fig. 4 in the main
text, we carry out the fitting routine separately for pixel-resolved energy-
distribution curves in the momentum distribution (1 pixel corresponds to
$\approx$ 0.02 Å-2). To optimize the signal-to-noise ratio, we apply a three-
fold rotational symmetrization and a mirror symmetrization that match with
crystal symmetry. Nevertheless, the signal-to-background ratio in a small
region around the center of the photoemission horizon remains below 1, and we
therefore exclude this region in our analysis (grey areas in Fig. 4 in the
main text). Also, for the momentum-resolved data, the replica HOMO background
signals due to the 13th harmonic amounts to 0, 1 or at most 2 counts in the
pixel-resolved (momentum-resolved) energy-distribution curves, and can
therefore not be fitted and subtracted accurately as described above for the
momentum-integrated data. As such, we need to ignore the HOMO replicas from
the 13th harmonic in the momentum-resolved analysis. We avoid overfitting of
the model in Eq. 3 by fixing the energy and bandwidth of the peaks in the
fitting routine to the parameters given in Table 1. Thus, the set of free
parameters in the momentum-resolved fitting procedure is limited to
$A_{S\textsubscript{3}{}}(k)$, $A_{S\textsubscript{2}{}}(k)$,
$A_{S\textsubscript{1}{}}(k)$ and $A_{bg}(k)$. This approach enables the
extraction of reliable momentum distributions also for the partially
overlapping energy distributions of the S1 and S2. The $1\sigma$ errors for
the full momentum maps are shown in Extended Fig. 7.

We note that the overall photoemission intensity of the S4 peak in the $h\nu=$
3.6 eV excitation data is comparably low due to the sub-50 fs decay to the
lower-energy S3 excitons (see Fig. 3 in main text, Extended Fig. 6b).
Furthermore, with $h\nu=$ 3.6 eV excitation, two-photon photoemission with 2
$*$ 3.6 eV = 7.2 eV is sufficient to overcome the work function, so that
space-charge effects could only be avoided by considerably reducing the
$h\nu=$ 3.6 eV pump intensity. Therefore, the signal-to-noise ratio in these
measurements was not sufficient for a momentum-resolved analysis of the S4
exciton data.

###  VI.3 Calculation of the C60 exciton spectrum

The _ab initio_ calculation of the exciton spectrum of the C60 film was
performed in two steps, using a $GW$+BSE approach. For the static electronic
structure, we perform calculations for two unique C60 dimers, which have been
extracted from the known structure of the molecular film david_structure_1991;
david_structure_1992 (see Fig. 2c, dimers 1-2 and 1-4 respectively). Starting
from Kohn-Sham orbitals and energies of a ground state DFT calculation
(6-311G*/PBE0+D3) krishnan_basisset_1980; frisch_basisset_1984;
adamo_pbe0_1999; grimme_d3_2010 using ORCA 5.0.1 neese_orca_2012;
neese_update5_2022, we employ the Fiesta code jacquemin_is_2017 to self-
consistently correct the molecular energy levels by quasi-particle self-energy
calculations with the $GW$ approximation. To account for polarization effects
beyond the molecular dimer, we embed the dimer cluster in a discrete
polarizable model using the MESCal program davino_mescal_2014; li_bse_2016;
li_correlated_2017. We found that mimicking 2 layers of the surrounding C60
film in such a way resulted in the convergence of the band gap within 0.1 eV
with a removal energy from the highest valence level of 6.65 eV. The close
agreement of this quasi-particle energy with the experimentally determined
work function of 6.5 eV gives us additional confidence in the choice of our
embedding environment. The calculated quasi-particle energy levels are shown
in Fig. 2a. Here the finite width of the black bars actually arises from
multiple energy levels forming bands on the energy axis. We characterize them
according to symmetry dresselhaus_c60_1996 as HOMO-1, HOMO, LUMO, and LUMO+1
bands, each consisting of 18, 10, 6 and 6 energy levels per dimer,
respectively. Note that we combine HOMO-1 is made up by states from two
different irreducible representations of the isolated gas phase C60 molecule
which are practically forming a single band.

Building upon the $GW$ energies, we compute neutral electron-hole excitations
by solving the Bethe-Salpeter equation beyond the Tamm-Dancoff approximation
(TDA). This yields the excitation energies $\Omega_{m}$ and the electron-hole
coupling coefficients $X^{(m)}_{vc}$ for a series of excitons labelled with
$m$. We first analyze the resulting optical absorption spectrum which is shown
in Fig. 2c as black solid line. It reveals a prominent absorption band around
$h\nu=$ 3.6 eV that is well-known from gas-phase
spectroscopydai_ultravioletvisible_1994. Secondly, the dimer calculation
reveals a strong optical absorption at $h\nu=$ 2.8 eV as well as a weakly
dipole-allowed transition at $h\nu=$ 2 eV. Both of these transitions are known
to only appear in aggregated phases of C60 wang_aggregates_1993, and cannot be
understood by considering only a single C60 molecule. We refer to the exciton
bands around $\Omega=$ 1.9, 2.1, 2.8 and 3.6 eV as S1, S2, S3 and S4,
respectively.

In line with our classification of the $GW$ energy levels, one can group the
composition of the excitons into four categories according to the contributing
quasi-particle energy levels. As visualized in Fig. 2c, this shows that S1 and
S2 are almost completely described by HOMO to LUMO transitions. On the other
hand, S3 is predicted to have a small contribution of HOMO-1 character,
however this contribution is too small to be reliably measured in our
experiment. Finally, for S4 we observe a clear and almost equal mixture of
HOMO-1 and HOMO contributions, which is also confirmed by our measurements.
Contributions from lower lying valence bands are negligibly small in the
studied energy window.

###  VI.4 Calculation of the exciton momentum maps

Based on our Kohn-Sham orbitals and BSE excitation coefficients, we calculate
theoretical momentum maps for each exciton according to Eq. 2 following the
derivation of Kern et al.kern_exciton_tomog. Note that for better readability,
Eqs. (1) and (2) are given within the TDA, however, can in general be extended
to include also de-excitation terms.kern_exciton_tomog In the present case, we
found the de-excitation contributions to be marginal (below 1%) without
affecting the appearance of the momentum maps and our interpretation. The
energy conservation term comprises the BSE excitation energies ($\Omega_{m}$),
the $GW$ quasi-particle energies for electron removal (i.e. ionization
potential $\varepsilon_{v}$), and the probe energy of $h\nu=$ 26.5 eV in
accordance with our experimental setup. Furthermore, we include an inner
potential to correct for the photoemission intensity variation of 3D molecules
along the moment vector component perpendicular to the surface. Here, we
choose a value of 12.5 eV, which has already been shown to match with
experimental C60 data hasegawa_potential_1998; haag_signatures_2020. (Note
that while considering the inner potential of the film is essential to
describe the ARPES fingerprint, a variation of the inner potential between 12
and 14 eV indicated no influence on the interpretation of our results.) As we
exploited the plane-wave approximation, the calculated photoemission intensity
is modulated by the momentum-dependent polarization factor
$\left|\bm{A}\cdot\bm{k}\right|$, which we modelled as p-polarized light
incoming with 68∘ to the surface normal according to experiments. To account
for the symmetry of the C60 film, the momentum maps were 3-fold rotated and
mirrored. Finally, application of Eq. 2 provides us with a 4D data set of
simulated photoemission intensity as a function of the excitation energy
($\Omega_{m}$), the kinetic energy ($E_{kin}$) and the momentum components
$k_{x}$ and $k_{y}$. Analogous to experiment, we referenced the kinetic energy
against the energy of the HOMO. The calculated ionization potential was
further used to set the photoemission horizon of the theoretical momentum
maps. Finally, to arrive at the theoretical momentum maps shown in Fig. 4, we
sum up the photoelectron intensities of each contributing electron-hole
transition in an kinetic energy range of 200 meV centered on the respective
exciton band.

##  VII Extended Figures

Figure 5: Static photoelectron momentum microscopy of the multilayer C60
sample, taken from the time-resolved data at -1000 fs. The high contrast and
line shape of the occupied molecular orbitals as shown in the energy-momentum
cut (a) and integrated energy distribution curve (b) confirm that an influence
of the Cu(111) substrate can be ignored. c, d: Furthermore, the clear
modulation of the HOMO (c) and HOMO-1 (d) momentum maps confirm the high
crystallinity of the multilayer at 80 Khaag_signatures_2020 Figure 6: a, b:
Exciton-resolved measurement of the femtosecond relaxation dynamics after
$h\nu=$2.9 and $h\nu=$3.6 eV excitation, respectively. The relative
intensities were acquired by fitting the model in Eq. 3 and 4 to the momentum-
integrated trPES data. Analysis of the time dependence indicates a sub 50 fs
lifetime of the S4 states. The comparably large S2 population in the $h\nu=$
3.6 eV measurement can potentially be explained by a relaxation of the S4 into
both the S3 and S2 states. Figure 7: Amplitude (a,b,c,d) and standard
deviation (e,f,g,h) of the extracted momentum maps of the S1, S2, S3 and the
exponential background ($A_{bg}(k)$), respectively. The data are scaled to the
mean value inside the circled area.