Femtosecond dynamics of correlated many-body states in C₆₀ fullerenes

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–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–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 [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₆₀ 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₆₀ 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₆₀ 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–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₆₀ 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 with 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₆₀ fullerenes [20–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 density-functional 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₆₀ 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 highly-excited 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\times {10}^{10}$ W cm⁻². 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 figure 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 time-dependent electronic structure of intermediate states, which are 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 DOS of an anion deposited on a metal surface. Our experiments are compared to ab initio calculations.

Zoom In Zoom Out Reset image size

This paper is organized as follows. In section 2 many-body states below the C₆₀ ionization threshold are calculated as a guideline for the 2PPE experiments. Section 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 section 4.1 and discussed in terms of resonance-enhanced ionization and excited state lifetimes. Section 4.2 concentrates on the time-resolved experiments. A detailed theoretical analysis of the experimental data is given in section 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.

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 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] with 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 ${\phi }_{i}({\bf{r}})$ are confined (a sphere with a radius of 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 $| {{\rm{\Phi }}}_{\alpha }\rangle$ by

$$\begin{eqnarray}&&| {{\rm{\Phi }}}_{\alpha }\rangle =\displaystyle \sum _{i\in \mathrm{occ}}\displaystyle \sum _{j\in \mathrm{virt}}{A}_{{ij}}^{\alpha }{\hat{c}}_{j}^{\dagger }{\hat{c}}_{i}| {{\rm{\Phi }}}_{0}\rangle ,\end{eqnarray} \tag{ 1 }$$

where $| {{\rm{\Phi }}}_{0}\rangle$ denotes the determinant built by the ground-state KS orbitals and ${\hat{c}}_{i}$ (${\hat{c}}_{i}^{\dagger }$) is the annihilation (creation) operator with respect to the KS basis. The particle-hole amplitudes ${A}_{i,j}^{\alpha }$ are determined by Casida's equation based on linear response. The major approximation hereby is related to the exchange-correlation (xc) kernel ${f}_{\mathrm{xc}}({\bf{r}},{{\bf{r}}}^{\prime };\omega )$, defined as the functional derivative of the KS potential with respect to the density. We use the local-density approximation 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}}^{\alpha }$ with octopus, taking 75 occupied and 60 virtual orbitals into account, yielding well-converged results for excitation energies up to 10 eV. For testing purposes we also computed the binding energies of the virtual orbitals analogously to [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 figure 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 [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 [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 many-body state, which is identified as 6T_1u excitation [40]. The projection onto the KS basis is depicted in figure 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 figure 2. As clear from figures 1(b) and 2, the excitation G is predominantly composed of the transitions (i) from HOMO to virtual states above the first SAMOs with dominant angular momentum of L = 8 and L = 6, and (ii) from HOMO-1 to LUMO+2. The angular-momentum analysis of the individual orbitals (figure 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 HOMO. This is consistent with a full symmetry analysis of the initial orbitals and the ab initio calculations.

Zoom In Zoom Out Reset image size

In order to study electronic transitions into highly excited many-body states close to the ionization potential, fs pulses at 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 state-of-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 figure 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 $\approx 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 (BBO, 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 $\approx 2$ μJ.

Zoom In Zoom Out Reset image size

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 the 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₆₀ 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\times {10}^{10}$ W cm⁻², 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 C₆₀ powder (purity >99.5%) 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.

4.1. 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 $\approx 250$ laser shots for each excitation wavelength. The ${{\rm{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 figure 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 a resonance-enhanced two-photon photoionization at ${\lambda }_{\mathrm{exc}}=218\pm 0.5\,\mathrm{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}_{-3}^{+5}$ 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.

Zoom In Zoom Out Reset image size

4.2. Pump–probe delay dependence

Time-resolved mass spectrometry traces the excited state dynamics in neutral C₆₀ 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 figure 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 $\approx 3000$ laser shots per delay point. The time-dependent on-resonance ion signal shown in figure 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.

Zoom In Zoom Out Reset image size

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:

$$\begin{eqnarray}&&{S}_{\mathrm{ac}}(\tau )={\int }_{-\infty }^{+\infty }I(t-\tau )I(t){\rm{d}}{t},\end{eqnarray} \tag{ 2 }$$

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 ${\tau }_{\mathrm{vib}}$ the incoherent term is given by:

$$\begin{eqnarray}&&{S}_{\mathrm{inc}}(\tau )={\int }_{-\infty }^{+\infty }{S}_{\mathrm{ac}}(t-\tau )\exp (-| t| /{\tau }_{\mathrm{vib}}){\rm{d}}{t}.\end{eqnarray} \tag{ 3 }$$

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:

$$\begin{eqnarray}&&{S}_{\mathrm{tot}}(\tau )={a}_{\mathrm{ac}}{S}_{\mathrm{ac}}(\tau )+{a}_{\mathrm{inc}}{S}_{\mathrm{inc}}(\tau )+{a}_{\mathrm{bg}}\end{eqnarray} \tag{ 4 }$$

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 ${\tau }_{\mathrm{vib}}$ from the measurement, the experimental data is fit by the least squares method using expression (4). The background is set to ${a}_{\mathrm{bg}}=1$ and the laser pulses are assumed to be Gaussian with ${\tau }_{\mathrm{FWHM}}=100$ fs. Other variables, i.e., a_ac, a_inc, and ${\tau }_{\mathrm{vib}}$, are free fit parameters. The best fit curve (black line in the bottom graph of figure 5) yields a time constant ${\tau }_{\mathrm{vib}}=400\pm 100$ fs (95% confidence band) and an amplitude ratio ${a}_{\mathrm{ac}}\,:{a}_{\mathrm{inc}}\,:{a}_{\mathrm{bg}}=0.24\,:1.13\,:1$. The laser intensities in the interaction region in the present experiment are as low as $3.5\times {10}^{10}$ W cm⁻² which makes the contribution from the direct (nonlinear) two-photon process small. The ratio ${a}_{\mathrm{inc}}\,:{a}_{\mathrm{bg}}$ is close to 1 as expected for single photon ionization from an occupied transient state.

The observed exponential time constant ${\tau }_{\mathrm{vib}}=400\pm 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.

To 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 photovoltaics the Born–Oppenheimer (BO) approximation is an inherent limitation. 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–53] has been employed for the C₆₀ 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 section 2, i.e. electronic correlations can be taken into account only with great efforts.

5.1. 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₆₀ [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 figure 6. Our results compare very well with those tabulated in the literature, for instance table 6.2 in the book [48].

Zoom In Zoom Out Reset image size

As inferred from figure 6 the high-energy modes (which affect the electrons at 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 amplitudes ${\hat{Q}}_{\nu }$. On the KS level, the electron–vibron matrix elements are thus given by

$$\begin{eqnarray}&&{k}_{{ij}}^{\nu }=\langle {\phi }_{i}| {\left.\displaystyle \frac{\partial {v}^{\mathrm{KS}}}{\partial {Q}_{\nu }}| {\phi }_{j}\rangle \right|}_{{Q}_{\nu }=0}.\end{eqnarray} \tag{ 5 }$$

Details on the evaluation of equation (5) are provided in appendix A. Since we are opting for a model in the many-body basis of the excitations discussed in section 2, the matrix elements (5) are transformed according to equation (1) (see appendix A). We thus obtain the model Hamiltonian (atomic units are used throughout)

$$\begin{eqnarray}&&\hat{H}(t)={\hat{H}}_{\mathrm{el}}(t)+{\hat{H}}_{\mathrm{el}-\mathrm{vib}}+{\hat{H}}_{\mathrm{vib}},\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{el}}(t)=\displaystyle \sum _{\alpha }{E}_{\alpha }| {{\rm{\Phi }}}_{\alpha }\rangle \langle {{\rm{\Phi }}}_{\alpha }| +f(t)\displaystyle \sum _{\alpha \beta }{M}_{\alpha \beta }| {{\rm{\Phi }}}_{\alpha }\rangle \langle {{\rm{\Phi }}}_{\beta }| ,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{el}-\mathrm{vib}}=g\displaystyle \sum _{\alpha \beta }\displaystyle \sum _{\nu }{K}_{\alpha \beta }^{\nu }| {{\rm{\Phi }}}_{\alpha }\rangle \langle {{\rm{\Phi }}}_{\beta }| {\hat{Q}}_{\nu },\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{vib}}=\displaystyle \frac{1}{2}\displaystyle \sum _{\nu }\left(\displaystyle \frac{{\hat{P}}_{\nu }^{2}}{{M}_{\nu }}+{k}_{\nu }{\hat{Q}}_{\nu }^{2}\right).\end{eqnarray} \tag{ 9 }$$

Here, ${M}_{\alpha \beta }$ are the dipole transition matrix elements from our TDDFT calculations, while f(t) comprises the time-dependent fields. The prefactor g in equation (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 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 [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. The laser-induced dynamics is, besides the electron-vibron interaction, 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–61].

The time evolution of the occupation of the states depicted in the level scheme figure 1(a) is presented in figure 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 ${\rm{\Delta }}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\times {10}^{10}$ W cm⁻². The values for g and η are chosen to match the time-dependent pump-probe signal observed in the experiment (see section 5.3).

Zoom In Zoom Out Reset image size

As one can infer from figure 7, the population transfer between the ground state and the bright excited states is clearly not in the 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 $\approx 5.5\,\mathrm{eV}$ (see figure 7(b)). Closer inspection reveals that these DSs involve the excitation of the 3s and the 3p SAMOs; the respective weights are given in figure 7(a). This relaxation mechanism is the dominant consequence of the electron-vibron coupling. This behavior is expected, as dissipation pathways are generally preferred compared to bath-induced excitations. These thermalization processes are hence less pronounced and occur on a longer time scale.

5.2. Pump-probe dynamics

In order to compute the pump-probe signals from the dynamics of the density matrix ${\rho }_{\alpha \beta }(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

$$\begin{eqnarray}&&{N}_{{\bf{k}}}\propto \mathrm{Re}\displaystyle \sum _{\alpha \beta }{\int }_{-\infty }^{\infty }{\rm{d}}t{\int }_{-\infty }^{t}{\rm{d}}{t}^{\prime }{F}_{\mathrm{probe}}(t){F}_{\mathrm{probe}}({t}^{\prime }){{\rm{e}}}^{-{\rm{i}}({E}_{\alpha }+\omega -{E}_{\beta }^{+}-{\epsilon }_{{\bf{k}}})(t-{t}^{\prime })}| {M}_{\beta {\bf{k}},\alpha }{| }^{2}{\rho }_{\alpha \alpha }({t}^{\prime })\end{eqnarray} \tag{ 10 }$$

for the number ${N}_{{\bf{k}}}$ of released photoelectrons with momentum ${\bf{k}}$ and energy ${\epsilon }_{{\bf{k}}}$. The probe laser pulse is assumed as ${f}_{\mathrm{probe}}(t)={F}_{\mathrm{probe}}(t){{\rm{e}}}^{-{\rm{i}}\omega t}$ with the pulse envelop F_probe. The matrix element ${M}_{\beta {\bf{k}},\alpha }=\langle {{\rm{\Phi }}}_{\beta ,{\bf{k}}}| \hat{D}| {{\rm{\Phi }}}_{\alpha }\rangle$ ($\hat{D}$ denotes the dipole operator) describes the transition from the intermediate states $| {{\rm{\Phi }}}_{\alpha }\rangle$ to the final states $| {{\rm{\Phi }}}_{\beta ,{\bf{k}}}\rangle$ with one photoelectron and the ion state labeled by β (energy ${E}_{\beta }^{+}+{\epsilon }_{{\bf{k}}}$). Note that only the population of excited intermediate states and not the full density matrix enters equation (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 $| {\bf{k}}\rangle$, $| {{\rm{\Phi }}}_{\beta ,{\bf{k}}}\rangle$, we write the usual anti-symmetrized product ansatz

$$\begin{eqnarray}&&| {{\rm{\Phi }}}_{\beta ,{\bf{k}}}\rangle ={\hat{c}}_{{\bf{k}}}^{\dagger }| {{\rm{\Phi }}}_{\beta }^{+}\rangle ,\end{eqnarray} \tag{ 11 }$$

where $| {{\rm{\Phi }}}_{\beta }^{+}\rangle$ is an eigenstate of the ionized system with energy ${E}_{\beta }^{+}$. In this case the photoemission matrix element reduces to

$$\begin{eqnarray}&&\langle {{\rm{\Phi }}}_{\beta ,{\bf{k}}}| \hat{D}| {{\rm{\Phi }}}_{\alpha }\rangle =\langle {\bf{k}}| \hat{D}| {\phi }_{\alpha \beta }^{{\rm{D}}}\rangle .\end{eqnarray} \tag{ 12 }$$

Here, ${\phi }_{\alpha \beta }^{{\rm{D}}}({\bf{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

$$\begin{eqnarray}&&| {{\rm{\Phi }}}_{\beta =(m,\alpha )}^{+}\rangle \approx {\hat{c}}_{m}| {{\rm{\Phi }}}_{\alpha }\rangle ,\quad {E}_{\beta }^{+}\approx {E}_{\alpha }-{\epsilon }_{m}.\end{eqnarray} \tag{ 13 }$$

Here, ${\epsilon }_{m}$ stands for the KS eigenvalue of orbital ${\phi }_{m}({\bf{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 $| {\bf{k}}\rangle$ 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 equation (10) balances spectral resolution versus temporal resolution in terms of the convolution of the phase factor ${{\rm{e}}}^{-{\rm{i}}{\rm{\Delta }}{Et}}$, ${\rm{\Delta }}E={E}_{\alpha }+\omega -{E}_{\beta }^{+}-{\epsilon }_{{\bf{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 pump-probe signal $S(\tau )=\int {\rm{d}}{\bf{k}}\,{N}_{{\bf{k}}}$ simplifies to

$$\begin{eqnarray}&&S(\tau )\propto \displaystyle \sum _{\alpha \beta }{\int }_{-\infty }^{\infty }{\rm{d}}t{\int }_{-\infty }^{t}{\rm{d}}{t}^{\prime }\,{F}_{\mathrm{probe}}(t-\tau ){F}_{\mathrm{probe}}({t}^{\prime }-\tau ){P}_{\alpha \beta }({E}_{\alpha }+\omega ){\rho }_{\alpha \alpha }({t}^{\prime }),\end{eqnarray} \tag{ 14 }$$

where ${P}_{\alpha \beta }(\epsilon )=k\int {\rm{d}}{{\rm{\Omega }}}_{{\bf{k}}}| {M}_{\beta {\bf{k}},\alpha }{| }^{2}$ with $k=\sqrt{2\epsilon }$ is proportional to the energy-resolved ionization probability with respect to the initial state $| {{\rm{\Phi }}}_{\alpha }\rangle$ and final state $| {{\rm{\Phi }}}_{\beta }^{+}\rangle$. As in section 4.2, τ denotes the pump-probe delay.

5.3. Theory versus experiment

5.4. Steady-state versus time-dependent picture

Ultrashort pulses in the UV spectral range excite correlated many-body states of isolated C₆₀ molecules below the ionization continuum. The population is followed by subsequent UV pulses of the same wavelength that ionize the molecules. By recording the ${{\rm{C}}}_{60}^{+}$ ion yield as a function of time delay between the pump and probe pulses we observe an exponential decay with a time constant of ${\tau }_{\mathrm{vib}}=400\pm 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 ${\tau }_{\mathrm{el}}={10}_{-3}^{+5}$ 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.

This research was supported by the Deutsche Forschungsgemeinschaft through the excellence cluster 'The Hamburg Centre for Ultrafast Imaging (CUI)—Structure, Dynamics and Control of Matter at the Atomic Scale' (DFG-EXC1074), the collaborative research centres 'Light induced dynamics and control of correlated quantum systems' (SFB 925), 'Functionality of Oxide Interfaces' (SFB 762), the GrK 1355 and Grant Number PA 1698/1-1.

In this appendix we provide the details on how the first-order (Herzberg–Teller) electron-vibron coupling matrix elements, equation (5), were computed. After the calculation of the vibrational eigenmodes ν, one readily obtains the associated deformations of the fullerene cage. Denoting the collection of nuclear coordinates in equilibrium by $\{{{\bf{R}}}^{(0)}\}=({{\bf{R}}}_{1}^{(0)},\ldots ,{{\bf{R}}}_{60}^{(0)})$, the vibrational distortion is characterized by

$$\begin{eqnarray}&&\{{{\bf{R}}}^{(\nu )}(q)\}=({{\bf{R}}}_{1}^{(\nu )}(q),\ldots ,{{\bf{R}}}_{60}^{(\nu )}(q))\quad \mathrm{with}\quad {{\bf{R}}}_{m}^{(\nu )}(q)={{\bf{R}}}_{m}^{(0)}+q{{\bf{V}}}_{m}^{(\nu )},\end{eqnarray} \tag{ A.1 }$$

where ${{\bf{V}}}_{m}^{(\nu )}$ is the displacement eigenvector. Based on equation (A.1), we performed a DFT calculation for each vibrational mode, fixing the magnitude of the distortions at $\delta q=0.01$. Instead of taking derivatives of the KS potential, we employ the equivalent formulation

$$\begin{eqnarray}&&{k}_{{ij}}^{\nu }\simeq \displaystyle \frac{1}{\delta q}\langle {\phi }_{i}| {\hat{h}}_{\nu }-{\hat{h}}_{0}| {\phi }_{j}\rangle ,\quad {\hat{h}}_{\nu }=-\displaystyle \frac{1}{2}{{\rm{\nabla }}}^{2}+{v}^{\mathrm{KS}}({\bf{r}},\{{{\bf{R}}}^{(\nu )}(\delta q)\}),\end{eqnarray} \tag{ A.2 }$$

where we have approximated the derivative by finite differences. The KS Hamiltonian describing the molecule in equilibrium is denoted by ${\hat{h}}_{0}$, while ${\hat{h}}_{\nu }$ describes the distorted molecule. For the evaluation of equation (A.2) we insert a (approximate) completeness relation and obtain

$$\begin{eqnarray}{k}_{{ij}}^{\nu } & \simeq & \displaystyle \frac{1}{\delta q}[\langle {\phi }_{i}| {\tilde{\phi }}_{k}^{(\nu )}\rangle {\tilde{\epsilon }}_{k}^{(\nu )}\langle {\tilde{\phi }}_{k}^{(\nu )}| {\phi }_{j}\rangle -{\epsilon }_{i}{\delta }_{{ij}}],\\ {\hat{h}}_{\nu }| {\tilde{\phi }}_{k}^{(\nu )}\rangle & = & {\tilde{\epsilon }}_{k}^{(\nu )}| {\tilde{\phi }}_{k}^{(\nu )}\rangle ,{\hat{h}}_{0}| {\phi }_{i}\rangle ={\epsilon }_{i}| {\phi }_{i}\rangle .\end{eqnarray} \tag{ A.3 }$$

Note that the overlaps $\langle {\phi }_{i}| {\tilde{\phi }}_{k}^{(\nu )}\rangle$ incorporate the symmetry properties of the vibronically induced transitions.

The transformation into the many-body basis is accomplished by expressing the one-body coupling operator in second quantization, ${\hat{k}}^{\nu }={\sum }_{{ij}}{k}_{{ij}}^{\nu }{\hat{c}}_{i}^{\dagger }{\hat{c}}_{j}$, and evaluating $\langle {{\rm{\Phi }}}_{\alpha }| {\hat{k}}^{\nu }| {{\rm{\Phi }}}_{\beta }\rangle$ according to the algebra of the fermionic creation and annihilation operators.

For the derivation of the Lindblad master equation in the weak-coupling limit, we follow [58]. For the Hamiltonian (7)–(9) one obtains the following EOM for the electronic density matrix $\hat{\rho }(t)={\sum }_{\alpha \beta }{\rho }_{\alpha \beta }(t)| {{\rm{\Phi }}}_{\alpha }\rangle \langle {{\rm{\Phi }}}_{\beta }|$:

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\hat{\rho }(t) & = & -{\rm{i}}[{\hat{H}}_{\mathrm{el}}(t),\hat{\rho }(t)]\\ & & +\displaystyle \sum _{\alpha \beta {\alpha }^{\prime }{\beta }^{\prime }}{{\rm{\Gamma }}}_{\alpha \beta {\alpha }^{\prime }{\beta }^{\prime }}\left({\rho }_{\beta {\beta }^{\prime }}(t)| {{\rm{\Phi }}}_{\alpha }\rangle \langle {{\rm{\Phi }}}_{{\alpha }^{\prime }}| -\displaystyle \frac{1}{2}\{| {{\rm{\Phi }}}_{{\beta }^{\prime }}\rangle \langle {{\rm{\Phi }}}_{\beta }| \hat{\rho }(t)\}\right).\end{eqnarray} \tag{ B.1 }$$

The square (curly) brackets denote the commutator (anti-commutator). The vibronic bath enters into

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\alpha \beta {\alpha }^{\prime }{\beta }^{\prime }}=\displaystyle \sum _{\nu }{\gamma }_{\nu }({E}_{\beta }-{E}_{\alpha }){\delta }_{{E}_{\beta }-{E}_{\alpha },{E}_{{\beta }^{\prime }}-{E}_{{\alpha }^{\prime }}}{K}_{\alpha \beta }^{\nu }{K}_{{\alpha }^{\prime }{\beta }^{\prime }}^{\nu },\end{eqnarray} \tag{ B.2 }$$

where ${\gamma }_{\nu }(E)=({N}_{{\rm{B}}}(E)+1){A}_{\nu }(E)$. ${N}_{{\rm{B}}}(E)$ denotes the Bose distribution (displayed in figure 6) which accounts for the occupation of the vibronic modes for the given temperature. The vibrational frequencies determine the spectral function ${A}_{\nu }(E)=2\pi \delta (E-{{\rm{\Omega }}}_{\nu })$, which we replace by the smeared form

$$\begin{eqnarray}&&{A}_{\nu }(E)=\sqrt{2\pi /{\eta }^{2}}\exp [-{(E-{{\rm{\Omega }}}_{\nu })}^{2}/2{\eta }^{2}]\end{eqnarray} \tag{ B.3 }$$

to account for the finite lifetime of the vibrations.