Femtosecond dynamics of correlated many-body states in C$_{60}$ fullerenes

Fullerene complexes may play a key role in the design of future molecular electronics and nanostructured devices with potential applications in light harvesting using organic solar cells. Charge and energy flow in these systems is mediated by many-body effects. We studied the structure and dynamics of laser-induced multi-electron excitations in isolated C$_{60}$ by two-photon photoionization as a function of excitation wavelength using a tunable fs UV laser and developed a corresponding theoretical framework on the basis of \emph{ab initio} calculations. The measured resonance line width gives direct information on the excited state lifetime. From the spectral deconvolution we derive a lower limit for purely electronic relaxation on the order of $\tau_\mathrm{el}=10^{+5}_{-3}$ fs. Energy dissipation towards nuclear degrees of freedom is studied in time-resolved techniques. The evaluation of the non-linear autocorrelation trace gives a characteristic time constant of $\tau_\mathrm{vib}=400\pm100$ fs for the exponential decay. In line with the experiment, the observed transient dynamics is explained theoretically by nonadiabatic (vibronic) couplings involving the correlated electronic, the nuclear degrees of freedom (accounting for the Herzberg-Teller coupling), and their interplay.


Introduction
Molecular junctions, molecular transistors and organic solar cells rely on charge transport channels with negligible energy dissipation during the carriers propagation time. In nanostructured materials and molecular complexes the characteristic timescale is determined by the long-range polarization interaction and by the formation and breaking of chemical bonds mediated by the electronic and nuclear motion. Transient structures and dynamics on the femto and sub-femtosecond timescale is the focus of ultrafast spectroscopy. Time-resolved experiments using femtosecond (fs) laser pulses unravel the dynamic response of promising materials that could serve for instance as molecular building blocks for organic photovoltaics. Polymer solar cells are commonly composed of a photoactive film of a conjugated polymer donor and a fullerene derivative acceptor [1,2,3], which makes use of the fullerenes' unique ability to form stable C − 60 anions. Electron correlation plays an important role in the formation of four bound states of the fullerene anion [4,5,6]. In fact, electronic correlations are responsible for the binding of the 2A g state, whereas the bindings of the states 2T 1u , 2T 2u and 2T 1g are less affected by electronic correlations (cf. Ref. [4] and further references therein).
With its special structure consisting of 174 nuclear degrees of freedom, 60 essentially equivalent delocalized π electrons, and 180 structure-defining localized σ electrons, neutral C 60 serves as a model for a large -but still finite -molecular system with many electronic and nuclear degrees of freedom. Because of the large charge conjugation, its finite "energy gap", and quantum confinement of electronic states, C 60 may be viewed as an interesting intermediate case between a molecule and a condensed matter system. In fact, applying solid-state concepts to the valence "Bloch electrons" on the C 60 sphere results in an "angular band structure" [7] from which other relevant quantities (such as plasma frequencies and group velocities) can be extracted. Photophysical studies of fullerenes using fs laser fields cover the whole range from atomic through molecular to solid state physics [8,9,10]. The molecular response is truly a multi-scale phenomenon. It ranges from attosecond dynamics in electronic excitation and ionization to statistical physics describing thermalization processes. So, light-induced processes in fullerenes cover more than 15 orders of magnitude in time [11].
Using low-temperature scanning tunneling microscopy (LT-STM) of C 60 molecules deposited on copper surfaces Feng et al. observed tunneling through electronic states that possess nearly atom-like character [12]. These "superatom" molecular orbitals (SAMOs), also discussed below, have a well-defined symmetry and can be characterized by the nodal structure (principle quantum numbers) n and angular momentum quantum numbers L [13]. In addition, these virtual states show a remarkable stability [14], i.e., their initial stage of decay proceeds substantially slower than other states (even LUMO or HOMO) which qualifies them as robust channel for hot electron transport. From a molecular physics point of view, SAMOs can be regarded as low-lying, mixed valence Rydberg states [15] that exhibit substantial electron density inside the hollow sphere. They form chemical bonds affected by hybridization when the system is excited optically or probed in STM in deposited nanostructures. The resulting free-electron bands in self-assembled one-dimensional wires or two-dimensional quantum wells are holding great promise for unique applications in molecular electronics [16], but also for new functionalities such as current carrying states and hence nanometer-sized magnetic field generators [17]. We note in passing that recently SAMO states have also been observed in planar, non-fullerene materials [18,19], as well as in isolated C 60 fullerenes [20,21,22], where they are energetically located below the known high-lying Rydberg states [23,24].
It is of fundamental importance for designing fullerene derivatives as building blocks for solid state chemistry to go beyond the characterization of static many-body electronic structure [25]. In particular for optimization and control of the charge flow and energy dissipation, rigorous dynamic studies on the fs time scale using ultrafast lasers are indispensable but, still in their infancy [26,27]. Additionally, time-dependent densityfunctional theory (TDDFT) calculations of the absorption and photoelectron spectra, accounting for full structural analysis, were performed [28,29] and constitute the basis for the further development presented below.
Accessing energy dissipation upon the excitation of correlated many-body states in gas phase C 60 became feasible recently which allows to connect the coherent quantum [28,30], classical and statistical mechanisms [31]. It is known from experiments with optical lasers that electron thermalization mediated by inelastic electron-electron collisions takes place on a time scale of 50−100 fs (see [11] and references therein). This is where the present experimental and theoretical work comes into play. Here, the objective is to reduce the complexity of laser-induced multi-photon processes by populating highlyexcited many-body states in a resonant one-photon transition in the ultraviolet (UV) spectral range at low laser peak intensity on the order of 3.5 × 10 10 W/cm 2 . This allows for a rather detailed probing of the correlated electron dynamics in highly excited states. The study is based on a resonance-enhanced multi-photon ionization (REMPI) scheme, i.e., two-photon photoemission (2PPE) spectroscopy as depicted in Fig. 1(a). The photoionization yield recorded for resonant excitation is enhanced as compared to an experiment performed in the off-resonance regime. Thereby, we trace the timedependent electronic structure of intermediate states free from any perturbation caused by metallic substrates affecting the energetics in LT-STM experiments. Furthermore, REMPI on gas phase fullerenes provides information on the neutral molecule whereas 2PPE and LT-STM essentially probe the binding energy and density of states of an anion deposited on a metal surface. Our experiments are compared to ab initio calculations. This paper is organized as follows. In Sec. 2 many-body states below the C 60 ionization threshold are calculated as guideline for the two-photon photoemission experiments. Sec. 3 describes some experimental details with a focus on the tunable fs laser system in the UV spectral range used for the time-resolved studies. Excitation energy dependent mass spectra are evaluated in Sec. 4.1 and discussed in terms of resonance-enhanced ionization and excited state lifetimes. Sec. 4.2 concentrates on the time-resolved experiments. A detailed theoretical analysis of the experimental data is given in Sec. 5 followed by a short summary and outlook. In Appendix A we provide details on how the electron-vibron coupling matrix elements were computed and the master equation in Lindblad form is derived in Appendix B.

Optical excitations
The first step of a REMPI or 2PPE experiment entails the calculation of excited states, which are typically more difficult to describe than ground-state properties. Often exact diagonalization (full configuration interaction) is not feasible for large systems in which case one may resort to only a few methods: TDDFT [32], equation-of-motion (EOM) quantum chemistry methods [33], and many-body perturbation theory (MBPT) based on a Green's function formulation [34]. For its reduced computational cost as compared to the other methods, we have employed the linear-response TDDFT approach (Casida's method) [35].
As a first step we calculated the Kohn-Sham (KS) orbitals using the octopus package [36]. A modified version of the asymptotically corrected functional by Leeuwen and Baerends [37], which was shown to considerably improve excited-state properties [38]. In order to account for a multitude of highly-excited states, we have chosen a relatively large box to which all KS states φ i (r) are confined (a sphere with radius 12Å with uniform grid spacing of 0.15Å). This ensures that higher virtual orbitals (including the SAMOs) are well represented. After converging the ground-state and computing a sufficient number of virtual orbitals, we computed the singly-excited (i.e., single particle-hole excitations) many-body excitations by Casida's method. Formally this procedure amounts to approximating the excited many-body states |Φ α by where |Φ 0 denotes the determinant built by the ground-state KS orbitals andĉ i (ĉ † i ) is the annihilation (creation) operator with respect to the KS basis. The particlehole amplitudes A α i,j are determined by Casida's equation based on linear response. The major approximation hereby is related to the exchange-correlation (xc) kernel f xc (r, r ; ω), defined as the functional derivative of the KS potential with respect to the density. We use the local-density approximation (LDA) for the xc kernel, as it is local and (within adiabatic TDDFT) frequency-independent. Casida's equation is thus transformed into an eigenvalue problem. We computed the Casida vectors A α ij with octopus, taking 75 occupied and 60 virtual orbitals into account, yielding wellconverged results for excitation energies up to 10 eV. For testing purposes we also computed the binding energies of the virtual orbitals analogously to Ref. [21]. We obtained very similar results for the low-lying states relevant for the present experiments.
The energies of the obtained many-body excitations are shown in Fig. 1(a), where we distinguish the states with vanishing dipole transition moment (these we refer to as dark states, DS) from the ground state GS, and optically accessible states (bright states, BS). The onset of the visible to UV (UV-vis) optical absorption is well documented (e.g., cf. Ref. [39] for a review and a comparison with previous experiments). Three distinct absorption peaks (labeled according to literature as C, E and G bands, respectively) are found in the UV-vis region. An overall agreement between our calculations of the optical absorption with the experimental results is found, though the excitation energies are slightly underestimated (which is typical for DFT calculations). Note, the specific experimental conditions may affect the spectral positions of the absorption peaks (as discussed in Ref. [39]). In particular, the present experiment probes the optical properties of isolated molecules by ultrafast pulses and thus potentially eliminates energy-loss channels (e.g. collisions and inter-or intra molecular decay) that might shift the absorption peaks to higher energy. For these reasons, we base the subsequent calculations on our DFT results without any adjustments. As detailed below, very good agreement to the present experiment is achieved, justifying this approach a posteriori. According to our calculations, excitation G represents one triply degenerate manybody state, which is identified as 6T 1u excitation [40]. The projection onto the KS basis is depicted in Fig. 1(b), where we illustrate the relative weight of the excitation from occupied (i) to virtual (j) orbitals by the thickness of the corresponding arrows. The relevant orbitals are also presented in Fig. 2. As becomes clear from Fig. 1(b) and Fig. 2, excitation G is predominantly composed of the transitions from (i) HOMO to virtual states above the first SAMOs with dominant angular momentum of L = 8 and L = 6, and (ii) HOMO-1 to LUMO+2. The angular-momentum analysis of the individual orbitals ( Fig. 2) clarifies why the transition to the many-body state associated to excitation G is optically allowed. Analogously one can conclude that neither the 3s- SAMO nor the 3p-SAMO can be populated in a direct single-photon dipole transition from the HOMO. This is consistent with a full symmetry analysis of the initial orbitals and the ab initio calculations.

Experimental setup
In order to study electronic transitions into highly excited many-body states close to the ionization potential, fs pulses at ultraviolet (UV) frequencies are required [41]. The present time-resolved mass-spectrometric study is based on second-order UV autocorrelation making use of a time-of-flight (TOF) mass spectrometer and stateof-the-art nonlinear optics. The laser setup for generating fs pulses in the range of 216-222 nm comprises a Ti:Sa laser system, an optical parametric amplifier (OPA) and several frequency mixing stages. The outline of the system is sketched in Fig. 3. The Ti:Sa laser (Amplitude Technologies) is the backbone of the overall generation scheme and provides 35 fs (FWHM) pulses with a pulse energy of 12 mJ behind the compressor at a repetition rate of 25 Hz and 800 nm central wavelength. This output is split into several arms. First, the beam is split in a 90:10 ratio. The more intense fraction is sent to a commercial OPA (Light Conversion, TOPAS-C + HE-TOPAS). The OPA is continuously tunable in the infrared spectral range (1140-3500 nm) and pumped by the Ti:Sa laser. For the subsequent frequency conversion in the UV the output of the OPA is tuned to 1200 nm. The 11 mJ of the 800 nm pump pulse are converted to ≈ 0.4 mJ at this wavelength. The low-intensity fraction of the Ti:Sa laser (10%) is equally split into two branches. One half is recombined with the 1200 nm output of the TOPAS in a β-barium borate (BBO, 0.2 mm thick) crystal to generate 480 nm light by sum frequency generation (SFG). The second half is frequency doubled in another BBO, and then together with the 480 nm beam directed to the third SFG stage (40 µm thick) to finally generate the 5.7 eV (218 nm) photons. The UV pulse energy was measured using a calibrated XUV photodiode and a pyro detector to be ≈ 2 µJ.
The output UV beam is directed into a vacuum chamber where it is split into two pulses by a reflective split-and-delay unit (SDU) in order to generate two synchronized pulse replicas. The complete nonlinear optical setup was simulated using the software package LAB2 [42] including dispersion induced by UV pulse propagation in air and through the 2 mm thick entrance window of the vacuum chamber. According to the calculation the 218 nm pulse duration in the interaction region is of the order of 100 fs FWHM with a spectral bandwidth of 2.8 nm. A coarse cross-correlation measurement performed between the 400 nm and 480 nm pulses of 150 fs FWHM supports the derived UV laser beam parameters.
The SDU consists of a Si split-mirror with one half mounted on a delay stage which can displace the mirror along its normal to set the time delay between the pump and probe pulses. The SDU is followed by a focussing mirror which spatially overlaps the two pulse replicas in the laser-sample interaction area. The laser beams are focused onto the C 60 molecular beam with a spherical mirror (f = 300 mm). Its reflectivity is above 80% in the 200-245 nm range. The beam waist in the interaction area is on the order of 150 µm. The maximum peak intensity reached in the experiments is approximately 3.5×10 10 W/cm 2 , which is derived from first principles based on Gaussian beam propagation, pulse energy, pulse duration and the far-field laser profile.
The molecular beam is produced by evaporation of gold grade C 60 powder in a resistively heated oven at 775 K. The UV laser beams are focused perpendicular to both, the effusive molecular beam and the TOF spectrometer axis. The ions created in the intersection volume are extracted by a static electric field (Wiley-McLaren configuration [43]), directed onto multichannel plates, and finally counted after amplification and discrimination by a digital oscilloscope. The mass resolution of the TOF spectrometer is 0.2% at M/q=720.

Excitation energy dependence
The photoionization signal recorded for resonant excitation is enhanced compared to an experiment performed off-resonance. The spectral width of the resonance yields information on the excited state lifetime. A UV wavelength scan was performed by tuning the IR wavelength of the OPA. The populated many-body state in the neutral molecule is subsequently ionized during the pulse duration of 100 fs. A full mass spectrum is accumulated over ≈ 250 laser shots for each excitation wavelength. The C + 60 yield was normalized to the relative pulse energy monitored by a photodiode. The C + 60 ion yield as a function of laser wavelength in the range of 216-222 nm is shown in Fig. 4. The cut-off at 216 nm corresponds to the OPA's lower wavelength limit of 1140 nm. Relatively large error bars at short wavelengths result from the corresponding low pulse energy and thus poor statistics. The C + 60 signal disappears for excitation wavelengths longer than 222 nm, thus representing the low-energy threshold of the resonance. The wavelength scan clearly indicates resonance-enhanced two-photon photoionization at λ exc = 218±0.5 nm. The observed width of the REMPI signal is of the order of 3.65 nm (94 meV) pointing towards ultrafast ionization within a lifetime that can be as short as 10 +5 −3 fs. The lifetime estimate is derived from the deconvolution of the observed resonance with a Gaussian laser pulse spectrum of 2.8 nm (FWHM) and a Lorentzian describing the homogeneous broadening.

Pump-probe delay dependence
Time-resolved mass spectrometry traces the excited state dynamics in neutral C 60 molecules directly in the time domain. The transient electronic structure is initiated by a UV 100 fs pump pulse and followed (probed) by a delayed pulse replica that photoionizes the molecule (see Fig. 1(a)). The single-shot mass spectra are taken at varying delay times between the UV pulses ranging from -60 fs to 900 fs with ≈ 3000 laser shots per delay point. The time-dependent on-resonance ion signal shown in Fig. 5 is derived by taking the average number of C + 60 counts for each delay and normalizing it to the relative pulse energy monitored as the UV stray light peak in the TOF spectrum. The pump-probe scan is repeated two times.
In the most general case the ion signal from a three-level system with a transient intermediate electronic state exposed to resonant excitation will result from two excitation pathways: direct two-photon photoionization from the ground state to the continuum and the ionization via the transient state (REMPI). Therefore, the total ion signal S tot can be described as a sum of three components [44]: the coherent term S ac (coherent artifact), the incoherent term S inc and a constant background a bg . The coherent artifact reflects the direct nonlinear ionization process and is proportional to the autocorrelation (AC) function of the laser pulse: where I(t) is the laser intensity and τ is the delay between the pump and probe pulses. The incoherent term S inc carries information about the population dynamics of the transient state. It is a convolution of the laser pulse AC with a symmetric decay function. In case of an exponential decay with a characteristic time constant τ vib the incoherent term is given by: The constant background a bg consists of contributions from each excitation pulse individually and is independent of the pump-probe delay. The total signal reads as: with a i being the relative amplitudes of the different components. In general the amplitudes have a ratio depending on the spatial overlap of the two pulses and the ionization pathway of the system. To extract τ vib from the measurement, the experimental data is fit by the least squares method using expression (4). The background is set to a bg = 1 and the laser pulses are assumed to be Gaussian with τ FWHM = 100 fs. Other variables, i.e., a ac , a inc and τ vib , are free fit parameters. The best fit curve (black line in the bottom graph of Fig. 5) yields a time constant τ vib = 400 ± 100 fs (95% confidence band) and an amplitude ratio a ac : a inc : a bg = 0.24 : 1.13 : 1. The laser intensities in the interaction region in the present experiment are as low as 3.5 × 10 10 W/cm 2 which makes the contribution from the direct (nonlinear) two-photon process small. The ratio a inc : a bg is close to 1 as expected for single photon ionization from an occupied transient state.
The observed exponential time constant τ vib = 400 ± 100 fs is significantly longer than the characteristic electron-electron interaction time derived from pump-probe spectroscopy [45] and pulse duration dependent studies [46,47] in the optical spectral range. It seems that electron thermalization mediated by inelastic electron-electron collisions on a time scale of 50−100 fs does not play a key role when high lying correlated many-body states are excited directly at rather low peak intensity. In the following a detailed theoretical description of the observed resonance and its time-dependent structure is discussed.

Simulations
In order to the describe the response of the molecule upon pulsed laser irradiation, the interplay between the electronic excitations and the vibrationally hot molecule has to be taken into account (for an overview on this topic we refer to the book [48] and further references therein). In full generality, an ab initio description for both the electronic and nuclear degrees of freedom and their coupling is not feasible currently. Hence, one has to rely on suitable approximations. In an Ehrenfest approach the electrons are described on the level of TDDFT and the nuclei are subject to classical equations of motion due to the forces exerted by the electron distribution. Despite the success of this molecular dynamics approach [49] for predicting vibrationally-assisted charge-transfer processes [50] in photovaltic molecules, the Born-Oppenheimer (BO) approximation is an inherent limitation hereby. For excited-state properties, where the nuclear and electronic dynamics are strongly mixed, BO-type molecular dynamics is not predictive. Molecular dynamics beyond the BO approximation [51,52,53] has been employed for the C 60 molecule [54]; merging such schemes with a treatment of the electronic excitations beyond the KS level is computationally too demanding for our system. Alternatively, one can treat the electrons in a single-particle atomic basis within a tight-binding model [55], removing the adiabaticity constraint with respect to the many-body states. Besides the inevitable empirical ingredient, such theory is also not directly compatible with the ab initio description of the many-body states in Sec. 2, i. e. electronic correlations can only be taken into account by great effort.

Initial laser-induced dynamics
In order to elucidate the laser-and vibration-induced dynamics we take a different angle. Since a considerable amount of energy is stored in thermally activated vibrations, which can only be transferred to the electronic subsystem in small portions, the vibrons can be treated as an effective heat bath for the electrons. A similar model has successfully been employed for incorporating the influence of the vibrations on charge-transfer processes in organic photovoltaic systems based on C 60 [56]. To construct an appropriate model for our case, several ingredients are required. For the vibrations we restrict ourselves to the harmonic approximation of the bottom of the BO surfaces. The vibronic eigenmodes along with their eigenfrequencies and reduced masses were computed using the Octopus code, as well. The resulting density of states (DOS) of the vibronic modes is shown in Fig. 6. Our results compare very well with those tabulated in the literature, for instance table 6.2 in the book [48].
As inferred from Fig. 6 the high-energy modes (which affect the electrons the most) are only weakly populated for T = 775 K. Therefore, the oscillations of the nuclei around their equilibrium positions can be considered small. Hence, the Herzberg-Teller (HT) expansion [57] of the full Hamiltonian, including electrons and nuclei, yields a reasonable description for both subsystems and their interaction. The first-order HT Hamiltonian amounts to approximating the electron-vibron coupling as linear in the mode amplitudeŝ Q ν . On the KS level, the electron-vibron matrix elements are thus given by Details on the evaluation of eq. (5) are provided in Appendix A. Since we are opting for a model in the many-body basis of the excitations discussed in Sec. 2, the matrix elements (5) are transformed according to eq. (1) (see Appendix A). We thus obtain the whereĤ Here, M αβ are the dipole transition matrix elements from our TDDFT calculations, while f (t) comprises the time-dependent fields. The prefactor g in eq. (8) is introduced as an overall scaling factor for the strength of the vibronic coupling. Ideally, g = 1 should be fixed; however, due to the perturbative description resulting from the HT expansion, the electron-vibron interaction might be underestimated. Hence, g is kept as a parameter. Instead of describing the dynamics of the full density matrix according to Hamiltonian (6) (which is a formidable task), we treat the vibrations, as explained above, as a heat bath. This allows to obtain a master equation for the density matrix in the electronic subspace only. Here we assume Markovian dynamics and thus employ the Lindblad master equation, following the standard derivation and formulation from Ref. [58]. More details are presented in Appendix B. This procedure requires an additional parameter: the vibronic broadening η, which corresponds to the lifetime of the vibrational modes. Note that the laser-induced dynamics is, besides the electronvibron interaction, basically treated on ab initio level endorsing the predictive power of the approach. Note that the electron-vibron coupling (8) includes, in principle, Jahn-Teller and Herzberg-Teller couplings, which where identified as the main mechanisms for vibrational coupling in fullerenes [59,60,61]. The time evolution of the occupation of the states depicted in the level scheme Fig. 1(a) is presented in Fig. 7(a). The driving pulse f (t) is chosen as a Gaussian pulse with a FWHM of 100 fs as in the experiment, while the central frequency is adjusted to the vertical excitation energy ∆E = 5.66 eV between ground state and the BSs corresponding to the G peak. The peak amplitude amounts to the intensity of 3.5 × 10 10 W/cm 2 . The values for g and η are chosen to match the time-dependent pump-probe signal observed in the experiment (see Subsec. 5.3).
As one can infer from Fig. 7, the population transfer between the ground state and the bright excited states is clearly not in a perturbative regime, as two Rabi cycles are apparent during the 100 fs UV pulse interaction. The relaxation dynamics, transferring part of the excitation to the dark states, takes place on two time scales: for short delays one can observe a rapid energy transfer, while for longer times the distribution thermalizes. The depletion dynamics of the laser-excited BS primarily takes place due to the coupling to two lower-energy states at ≈ 5.5 eV (see Fig. 7(b)). Closer inspection reveals that these DSs involve the excitation of the 3s and the 3p SAMOs; the respective weights are given in Fig. 7(a). This relaxation mechanism is the dominant consequence of the electron-vibron coupling. This behavior is expected, as dissipation pathways are generally preferred as compared to bath-induced excitations. These thermalization processes are hence less pronounced and occur on a longer time scale.

Pump-probe dynamics
In order to compute the pump-probe signals from the dynamics of the density matrix ρ αβ (t), as discussed above, an extension to the scattering states is required. However, a straightforward implementation of the Lindblad equation including both bound and unbound many-body states is not feasible. This is due to the large dimension of the Hamiltonian after the inevitable discretization of the continuum. Hence, we opt for a perturbation description which allows to compute the pump-probe dynamics from the time-dependent density matrix without incorporating the ionization dynamics explicitly. This is achieved by a method known from time-resolved photoelectron spectroscopy [62]. Adopting a straightforward derivation for the many-body case, one obtains for the number N k of released photoelectrons with momentum k and energy k . The probe laser pulse is assumed as f probe (t) = F probe (t)e −iωt with the pulse envelop F probe . The matrix element M βk,α = Φ β,k |D|Φ α (D denotes the dipole operator) describes the transition from the intermediate states |Φ α to the final states |Φ β,k with one photoelectron and the ion state labeled by β (energy E + β + k ). Note that only the population of excited intermediates states and not the full density matrix enters eq. (10). This is an approximation, which relies on the fact that pathway interferences play only a minor role for the considered two-step ionization process.
For the excited state with one electron in the continuum state |k , |Φ β,k , we write the usual anti-symmetrized product ansatz where |Φ + β is an eigenstate of the ionized system with energy E + β . In this case the photoemission matrix element reduces to Here, φ D αβ (r) stands for the corresponding Dyson orbital. As the sum over all excited states β is implied, a (computationally expensive) precise calculation of the Dyson orbitals can be omitted by approximating them by simple hole states, i.e., by assuming Here, m stands for the KS eigenvalue of orbital φ m (r). The matrix element (12) can thus be evaluated in terms of a superposition of matrix elements with respect to the KS orbitals. The scattering states |k were computed with respect to the spherically averaged KS potential [28], which is known to be an adequate approximation for angle-integrated quantities. Asymptotic corrections ensuring the r −1 behavior are incorporated by smoothly interpolating between the short-range and long-range regimes. Further orthogonalization with respect to the bound KS orbitals is performed.
To reflect the experimental situation, the integration over all photoelectron states has to be performed. Furthermore, we note that eq. (10) balances spectral resolution vs. temporal resolution in terms of the convolution of the phase factor e −i∆Et , ∆E = E α + ω − E + β − k , with the envelop function. Due to very long pulses as compared to one oscillation period, this convolution practically yields a Dirac δ-function with respect to the energy balance. Taking this into account, the total ionization pumpprobe signal S(τ ) = dk N k simplifies to where P αβ ( ) = k dΩ k |M βk,α | 2 with k = √ 2 is proportional to the energy-resolved ionization probability with respect to the initial state |Φ α and final state |Φ + β . As in Sec. 4.2, τ denotes the pump-probe delay.

Theory vs. experiment
The pump-probe signal based on the laser-driven and vibron-coupled dynamics of the density matrix can now be compared to the experiment. We remark that Eq. (14) assumes that the two pulses can be separated and does not account for the scenario of overlapping pulses. We thus limit the comparison between the theoretical calculations and the experiment, presented in Fig. 5, to the region of exponential decay for τ ≥ 300 fs. Furthermore, the dynamics presented in Fig. 7 indicates that the laser intensity exceeds the perturbative regime. However, even with stronger pulses, the bright states around 5.66 eV are the only accessible channels, as absorbing another photon leads to immediate ionization. This will, however, only affect the background signal. Therefore, we adjust the background S bg = S(τ → ∞) to the experiment. After solving the Lindblad master equation for various values of the parameters g and η, we found the best fit for g = 1.128 and η = 0.77 meV. The latter corresponds to a vibrational lifetime of ≈ 5.4 ps, which is in accordance with previous experiments [8]. The small deviation of the scaling factor g from unity underpins the predictive power of our treatment. Note that also g = 1 results in a decay dynamics which closely resembles, apart from the background, the experimental data, whereas varying g to smaller or larger values clearly deviates from the measurements.
We also calculated the ionization signal for varying photon energy ω and compared the resulting spectra in Fig. 4. To match the laser spectral bandwidth to the experiment, the obtained curves were convolved with a Gaussian having a FWHM of 2.8 nm. The theoretical spectra are centered around the vertical excitation energy, which is lower than what is observed in the experiment. This kind of underestimating of bandgaps and thus excitation energies is typical for most DFT calculations and for the C 60 molecular, in particular [40]. However, the theoretical and the experimental peak differ only by ≈ 10 meV. Generally, both the experimental as well as the theoretical spectra are considerably sharper as compared to optical absorption measurements. This is a major advantage of the current experimental setup: the dipole excitation dominates over possible loss channels, if the pulse strength is increased, which leads to a narrow resonance.

Steady-state vs. time-dependent picture
The interpretation of the present data on the relaxation dynamics requires the full picture including both, correlated electronic and nuclear degrees of freedom and their interactions [63,64]. It is known that highly excited electronic and vibrational C 60 states are strongly mixed [65,55]. In turn, relaxation channels open up that depopulate the electronically excited states by internal conversion close to conical intersections similar to electron-phonon coupling in solid state materials. We note that characteristic fragmentation patterns observed in TOF mass spectra using optical lasers as a function of pump-probe delay revealed (nonadiabatic) vibronic coupling on a time scale of τ vib = 200−300 fs [45], which is in good agreement with the present findings. Furthermore, nonadiabatic coupling of mixed valence and Rydberg states is ubiquitous in polyatomic molecules affecting potential energy surfaces, energy relaxation and dissociation dynamics [21]. Similar processes have been considered to evoke the loss of small neutral fragments from C 60 on a picosecond time scale [66].
The identification of vibronic coupling playing a key role in the electronic energy loss of correlated many-body states may open new vistas for optical control of chargetransport phenomena in smart materials containing these nanospheres. For instance, coherently induced radial symmetric "breathing" motion of the cage atoms strongly impacts the structure and dynamics of the molecule [67]. The carbon cage and the electron system start to exchange many eV in energy periodically on a sub-100 fs period timescale. The coherent oscillation prevails for several cycles [67], which might be interesting for novel ultrafast switching applications in molecular electronics.

Summary and Outlook
Ultrashort pulses in the ultraviolet spectral range excite correlated many-body states of isolated C 60 molecules below the ionization continuum. The population is followed by subsequent UV pulses of the same wavelength that ionize the molecules. By recording the C + 60 ion yield as a function of time delay between pump and probe pulses we observe an exponential decay with a time constant of τ vib = 400±100 fs which is explained by the coupling of electronic excitation to nuclear motion in the neutral molecule. The initial electronic relaxation can be as fast as τ el = 10 +5 −3 fs according to the evaluation of the resonance line width in single pulse experiments as a function of excitation wavelength. The experimental results are in good agreement with ab initio calculation of structure and dynamics including electronic correlation and vibronic coupling. However, the UV pulse duration of 100 fs did not allow to observe pure electron dynamics in time-resolved experiments. Future experimental work making use of shorter UV pulses shall reveal the predicted laser-driven Rabi oscillation and time-resolved transformation of the electronic orbitals, i.e., the coupling between different electronic states. density matrixρ(t) = αβ ρ αβ (t)|Φ α Φ β |: d dtρ (t) = −i Ĥ el (t),ρ(t) The square (curly) brackets denote the commutator (anti-commutator). The vibronic bath enters into where γ ν (E) = (N B (E) + 1)A ν (E). N B (E) denotes the Bose distribution (displayed in Fig. 6) which accounts for the occupation of the vibronic modes for the given temperature. The vibrational frequencies determine the spectral function A ν (E) = 2πδ(E − Ω ν ), which we replace by the smeared form to account for the finite lifetime of the vibrations.