Pyrene, a Test Case for Deep-Ultraviolet Molecular Photophysics

We determined the complete relaxation dynamics of pyrene in ethanol from the second bright state, employing experimental and theoretical broadband heterodyne detected transient grating and two-dimensional photon echo (2DPE) spectroscopy, using pulses with duration of 6 fs and covering a spectral range spanning from 250 to 300 nm. Multiple lifetimes are assigned to conical intersections through a cascade of electronic states, eventually leading to a rapid population of the lowest long-living excited state and subsequent slow vibrational cooling. The lineshapes in the 2DPE spectra indicate that the efficiency of the population transfer depends on the kinetic energy deposited into modes required to reach a sloped conical intersection, which mediates the decay to the lowest electronic state. The presented experimental–theoretical protocol paves the way for studies on deep-ultraviolet-absorbing biochromophores ubiquitous in genomic and proteic systems.

1 Experimental and theoretical methods

Experimental methods
This section describes the methods related to the experimental acquisition of the data presented in this paper.

Sample preparation
Pyrene was purchased from Radiant Dyes and used without further purification, with a purity of 95%. It was diluted in UV-quality ethanol (Carl Roth). The solutions were stirred with a glass-covered magnet for approximately 30 minutes at room temperature and then sonicated. Afterward the samples were filtered. Afterward the samples were filtered with 0.4 µm polypropylene syringe filters (Chromacol) and stored in a cool and dark place. The absorption spectra were measured using a UV2600 spectrophotometer (Shimadzu) in 1-mm quartz cells (Hellma), and additionally controlled in situ by measuring the transmission through a photodiode (UV-sensitive 2032 -Newport). Fluorescence spectra were measured with a Cary Eclipse fluorimeter, using a 1×1 cm fluorescence cell from Hellma. Excitation was 260 nm, excitation and emission slits were 20 and 1,5-2.0 mm respectively. Average time of acquisition was 2-5 seconds each wavelength.
The OD of the samples was adjusted to be 0.3-0.35 at maximum absorbance. In order to minimize light scattering and pulse broadening, the solutions were circulated through a home-built wire-guided jet [Pic15], with the thickness (path length) in the excited volume maintained at ≈ 250 µm. The excitation spot size was ≈ 50 µm and, with the jet pump speed used (3000-3500 rpm), a refreshment of excitation volume was guaranteed for the laser repetition rate of 1 kHz (see details in Ref. [Pic15]). The samples were excited with the UV-pulses having energies in range of 15 − 20nJ per beam.

Experimental setup
The experiments were conducted using an all-reflective 2DES-UV setup with passive phase stabilization based on a reflective diffractive optic beamsplitter (Holoeye) [Pro15], which allows us to perform both TA (transient absorption) and PE (photon echo) experiments with approximately 6-fs resolution. The collection of 2D spectra within a few hours is guaranteed by the long-term phase stability (≥ Λ/100) of the setup [Pro15]. A thin neutral density filter (OD = 2, BMW) in the local oscillator (LO) beam path created the probe pulse for the TA experiment and the heterodyning of the beam. The UV 1 pulses energies exciting the samples were 15 -20 nJ per beam. We also performed TA measurements without the LO-filter, in which case the energies of both pump and probe pulses were comparable. We did not detected a sufficient difference in the measured kinetics (both lifetimes and decay associated spectra -DAS) between the measurements with and without LO-filter. Because the LO-filter dispersion cause a chirp in the probe pulse, we phased the PE spectra with TA measurements done without the LO-filter. The broadband and coherent UV-pulses (6-7 fs FWHM) were produced using an achromatic frequency doubling of a VIS broadband spectrum generated from a home-built non collinear optical parametric generator pumped by a Coherent Elite USP laser system (Coherent) which delivers 800-nm laser pulses of 35 fs FWHM at the repetition rate of 1 kHz. The design of the achromatic doubling system, combined with a deformable mirror (Oko Technologies) for compensation of phase distortions in the UV-pulses, was similar to that described in [Bau04]. The temporal profiles of the UV-pulses were optimized using a genetic algorithm and measured by a frequency-resolved optical gating method (FROG) with a home-built setup [Pro15a], and directly in situ by recording the FROG-trace from a 150-mm fused silica cover slip. The temporal profiles of the pulses (after optimization) were characterized using a commercial software (Femtosecond Technologies).

Data collection
The frequency-resolved heterodyne photon echo signals were collected within a scan range of 128 fs with a step of 1 fs at fixed waiting times T; at each delay point 900-1800 single spectra were averaged (depending on the signal magnitude) to achieve high signal-to-noise ratio. Each photon echo scan took approximately 2-5 min (depending on averaging) per one delay point T. The pump-probe spectra were typically collected in a delay range of 5-10 ps with a step of 10-20 fs. At each delay point 3600 differential spectra (pump-on -pump-off) were averaged. To resolve the long-lived decay times, the delay scans were also performed in a 200 ps window with delay steps of 500 fs. In the present paper, we plot the differential spectra as measured, without normalization to the probe spectra. In order to achieve high spectral resolution necessary for resolving the oscillations, the delay step was set to 4 fs in heterodyne transient grating measurements and associated TAmeasurements.
The TA spectra were collected at the same time delays as the TG spectra and were individually used for phasing at each waiting time T. The contribution of the nonlinear solvent response which is more than 100 times higher than of the intrinsic signal and undesirable scattering effects from the sol-vent, prevented accurate measurements of intrinsic signal around zero delays which amplitude typically reaches ≈ 1.5mOD (corresponding to ∆T ≈ 40 counts) or, in absolute units ≈ 0.5%.

Data processing
The TA spectra were globally fitted to a multi-exponential decay model using a home-written software. Due to extremely strong contribution of the nonlinear solvent response and Rayleigh scattering we performed global analysis on the TA data starting 150 fs after excitation. Phasing of the 2D spectra is described in ref.
[Pro09] by taking into account the wavelength-dependent refraction index of the LO-filter material. The 2D spectra were then transformed to the time domain, filtered with a Gaussian-shaped kernel, and transformed back to the frequency domain. The aforementioned filtering lowers the resolution in 2D spectra to approximately 200 cm −1 .

Theoretical methods
This section describes the methods applied in the simulation of the TA and PE spectra.

Electronic structure computations
The characterization of the electronic structure of pyrene was performed using the complete active space self-consistent field (CASSCF) approach [Roo87], followed by second order perturbation theory (CASPT2) [And90] in its single-state (SS) flavor. The restricted active space (RAS) variation of the scheme (i.e. RASSCF/ RASPT2) [Mal90,Mal08] which reduces the number of simultaneous excitations was used thus allowing active spaces containing all valence π-orbitals. The electronic structure computations and geometry optimizations were performed within D 2h symmetry at the RASPT2/RASSCF(4, 8|0, 0|4, 8) level of theory thereby including all valence π-orbitals in the active space (see Figures S15-S18) and allowing for up to quadruple excitations. The computations were performed using the ab initio software package Molcas 8 [Aqu16] through its interface with COBRAMM [Wei18] allowing for parallelized numerical computation of the RASPT2 gradients. The stationary nature of the obtained critical points was verified by frequency computations (w/o symmetry). Due to the lack of analytical gradients for RASPT2/RASSCF in the quantum mechanical software of choice the frequency computations were performed at CASSCF(8,8) level, only after each geometry was first relaxed at this inferior level and it was verified that no significant changes occurred with respect to the RASPT2 optimization. The energies of the five lowest excited states of the system were recomputed at each critical point within their symmetry at SS-CASPT2/CASSCF(4, 8|0, 0|4, 16) level, thereby including further eight extravalence virtual orbitals for a better description of the dynamic correlation (see Figures S15-S18). The conical intersection optimizations were performed through a constrained energy minimization applying the projection approach by Bearpark at al. [Bea94], relying on numerical RASPT2 gradients and non-adiabatic couplings and implemented in COBRAMM[Wei18]. The CI optimizations were performed within C S symmetry at SS-CASPT2/SA-6-RASSCF(4, 8|0, 0|4, 7) level where the completely anti-bonding orbital was excluded from the active space to speed up computation. The ANO-L basis set contracted to C[3s,2p,1d]/H[2s1p] was used[Wid90] throughout and the Cholesky decomposition was adopted in order to speed up the evaluation of two-electron integrals [Aqu09].

Spectroscopy simulations
The third-order signal is simulated according to the experimental set-up in the rephasing (k I ) phase-matching direction in a box-CARS geometry with a fixed spatial orientation of the four pulses with wavevectors k 1 through k 4 , thereby scanning the time-interval t 1 between both pump pulses k 1 and k 2 from −t 1 to +t 1 . Thus, for positive delays t 1 (i.e. pump pulse k 1 arrives before pump pulse k 2 ) we compute the rephasing contribution to the nonlinear response, while for negative delays t 1 (i.e. pump pulse k 1 arrives after pump pulse k 2 ) we compute the non-rephasing contribution.
The nonlinear response can be expressed within the framework of cumulant expansion of Gaussian fluctuations (CGF) [Muk95, Muk04, Abr09]. It allows to calculate the shapes of electronic transition bands coupled to a bath for fluctuations with arbitrary time scales by using the formalism of line shape functions when the fluctuations follow Gaussian statistics. A detailed discussion of the working equations is provided in section 2. In summary, we used a bath with both slow and fast modes in the limit of decoupled populations and coherences (secular approximation) [Abr09, Abr07]. The slow bath modes are responsible for spectral diffusion during the three intervals t 1 , t 2 and t 3 . The Markovian approximation is applied to the fast modes. During t 1 and t 3 , these cause homogeneous line broadening, while during t 2 they induce population relaxation described by the Pauli master equation. Utilizing the outcome of the electronic structure calculations performed in this study, the acquired experimental data and the findings of previous TA experiments, we elaborated a sequential model that describes the population decay following excitation of the second bright ππ * state (see section 3.1). The Pauli master equation is solved for this model. We note that no explicit non-adiabatic dynamics was performed and the model is based on making reasonable hypothesis on the population state at different times, backed up by the static computational data.
The spectral lineshapes are simulated through composite line shape functions consisting of: a) a part responsible for the vibrational progression in the spectra described through the multidimensional uncoupled displaced harmonic oscillator (DHO) model [Muk95]); b) a part describing the homogeneous broadening due to coupling to a continuum of low-frequency modes represented by the line shape function of the semi-classical overdamped Brownian oscillator (OBO) [Muk95, Li94, But12]. In this work the OBO parameters were chosen so to reproduce the bandwidth in the linear absorption spectrum at room temperature (see section 3.3 and Figure 1 from main text). Within the DHO model the one-dimensional potential for the i − th state along each normal mode k is fully characterized by two parameters, frequency ω k , relative displacement d ik with respect to the GS equilibrium (see Figure S1). We applied three different protocols for obtaining the DHO parameters describing the various states involved in the relaxational dynamics. a) dynamics-based protocol: The parameters for describing the spectral dynamics of the GSB, SE and ESA associated with the S 4 state were obtained from an adiabatic SA-3-CASSCF(4,4)/6-31G* molecular dynamics simulation, following Newton's equations of motion for the nuclei in the second bright state, which is the third root with the utilized active space of frontier orbitals HOMO-1, HOMO, LUMO and LUMO+1. The nuclei were propagated for 250 fs with a time step of 0.5 fs and without initial kinetic energy (the so-called 0K trajectory) using the velocity Verlet algorithm implemented in our in-house molecular dynamics code COBRAMM interfaced with Molcas 8[Aqu16] for obtaining quantum-mechanical gradients[Alt07,Wei18]. Along the trajectory 100 roots were computed at a SS-CASPT2/SA-100-RASSCF(4, 8|0, 0|4, 8) level every 2 fs and the parameters ω k andd ik were obtained using a Fourier series to fit the temporal evolution of the electronic gaps E e (t)-E g (t) and E f (t)-E e (t). Further details are provided in sec. 3.2.1. b) normal mode analysis: The parameters for describing the dynamics of the dark 1B 3u state (not directly visible in the spectrum but affecting the spectral dynamics of the ESA), in which the population is eventually trapped, were obtained using a Cartesian coordinates projecting technique in which the geometrical deformations between CI(1B 1g /1B 2u /1B 3u ), the entry point on the 1B 3u PES, and the S 1 Figure S1: Displaced harmonic oscillator model for a three level system. minimum were projected on the normal modes of the system[Kur01]. Thereby we utilized ground state MP2 normal modes and frequencies.
For further details see section 3.2.2. c) gradient projection protocol: The parameters for describing the spectral dynamics of the ESA associated with the dark 1B 3u state, in which the population is eventually trapped, were obtained from electronic structure computations at the 1B 3u minimum. Energies, gradients and transition dipole moments were computed at CASPT2 level, ω k andd ik were obtained through projection techniques[Lee16, Fer12] (see section 3.2.2), utilizing ground state MP2 normal modes and frequencies. The protocol for obtaining the parameters is outlined in detail in section 3.
In particular, the gradient based approach used in earlier works reporting simulations of spectra probing in the visible [Nen18, Bor18] was abandoned in favour of a dynamics-based protocol for constructing the spectral densities of the deep-UV (i.e. up to 10 eV above the GS) excited states accessible (i.e. bright) from the initially populated second bright state of pyrene. This was necessary as due to the high state density in the deep-UV the calculation of numerical gradients, required by the projection protocol, becomes unreliable. The PP signal was computed as the marginal (with respect to the delay time t 1 ) of the two-dimensional Fourier transform (with respect to t 1 and t 3 ) of the nonlinear response of the system, which becomes equal to the third-order signal within the limit of temporally well separated ultrashort laser pulses (impulsive limit). Spectra simulations were performed with a development version of Spectron [Abr09]. The present model allows to incorporate (within certain limits) state-specific photo-induced absorptions and their spectral dynamics in photon echo and pump-probe spectroscopy simulations. We have recently applied this methodology to simulate with remarkable agreement to experiment time-resolved two-color UV-pump Vis-probe electronic spectrum of trans-azobenzene [Nen18].

Simulation protocol for nonlinear spectroscopy
The pump probe (PP) signal W (3) (t 2 , Ω 3 ) can be formally seen as the marginal of the Fourier transform of the third-order signal S (3) (t 1 , t 2 , t 3 ) [Muk95, Muk04, Abr09, Ham11, Abr07]. (1) In the following we will elaborate on the simulation protocol for S (3) (t 1 , t 2 , t 3 ), thus covering PP and two-dimensional (2D) nonlinear experiments. S (3) (t 1 , t 2 , t 3 ) is the convolution of the third order nonlinear response of the system R (3) (t 1 , t 2 , t 3 ) due to the interaction with three incident laser pulses with wave vectors k 1 , k 2 and k 3 with pulses E(r, t) = E(t)e ikr−iωt with central frequency ω and complex envelope E(t). Assuming temporally well separated ultrashort laser pulses (impulsive limit, i.e. the limit in which the pulse duration is less than a single vibrational oscillation period) the nonlinear reponse of the system R (3) (t 1 , t 2 , t 3 ) becomes equivalent to the third-order signal S (3) (t 1 , t 2 , t 3 ). In section 3.4 the influence of the pulse duration on the spectra is demonstrated.
can be written in terms of the perturbation of the system's density matrix by three external optical electric fields (the interaction with the optical filed being described in the impulsive limit through the coordinate dependent transition dipole moment operatorμ = ij µ ij i j ) and the field-free evolution of the density matrix during the intervals t 1 , t 2 and t 3 between the pulses: giving rise to eight independent contributions (Liouville pathways). ρ(t) is the time-dependent density matrix of the system, being in the ground state (GS) at equilibrium before the interaction with the first pulse (i.e. ρ(0) = g g ). For the electronic states i we consider the GS g, the manifold of singly excited states (ESs) e, e , accessible with the pump pulse (pair) and the manifold of higher lying states f , accessible through the probe pulse.
G(t)ρ(t) = e −iĤt ρ(0)e iĤt the retarded Green's function, describes the fieldfree system propagation governed by the vibronic HamiltonianĤ. The nonlinear response R (3) (t 1 , t 2 , t 3 ) depends on the spatial orientation of the wave vectors of the incident pulses, i.e. is emitted in a given phase-matched direction. By adjusting the experimental set-up (boxcar arrangement, collinear pump pulse arrangement) one can selectively detect sub-groups of Liouville pathways. For example, the rephasing nonlinear response R (3) is emitted in the rephasing phase-matching direction (k I = −k 1 + k 2 + k 3 ) and includes three possible interaction sequences of the incident electric fields with the system denoted as ground state bleach (GSB), excited state absorption (ESA) and stimulated emission (SE). Similarly, the nonrephasing nonlinear response can be obtained in the phase-matching direction k II = +k 1 − k 2 + k 3 . The definition of k I and k II phase-matching conditions in a boxcar arrangement is associated with a fixed order of the incident pulses (i.e. k 1 interacts first, k 2 second, k 3 third). By changing the temporal order of the first and second pulses one can record both k I and k II along the same phasematching direction. E.g. recording along the k I phase-matching direction by interchanging in time first and second pulses gives +k 2 − k 1 + k 3 , which is identical to a nonrephasing nonlinear response. In this particular experiment, the delay between both pump pulses is scanned from −τ to +τ along the phase-matching direction. Thus, for negative delay times τ the nonrephasing response is recorded, while for positive delay times τ the rephasing response is recorded. Finally, the 2D spectrum is obtained by Fourier transforming the entire signal. The simulations were performed so as to match the experimental conditions. Therefore, the resulting spectra are neither rephasing, nor nonrephasing, but rather the quasi-absorptive spectrum that could be also obtained by measuring rephasing and nonrephasing responses independently and plotting the real part of the sum of both representations. Such spectra are preferred for analysis as phase-twists are canceled to a large extent. Finally, the pump-probe response can be obtained either in the rephasing or in the nonrephasing directions, which become identical when t 1 = 0. Eq. 3 can be solved using second-order cumulant expansion within the framework of linear coupling of the system's energies to a Gaussian bath. This is known as the cumulant expansion of Gaussian fluctuations (CGF). It allows to calculate the shapes of electronic transition bands coupled to a bath for fluctuations with arbitrary time scales by using the formalism of line shape functions [Muk04, Abr07]. In the following we describe a bath with both fast and slow modes, i.e. modes that adjust themselves quickly and such that cause long-range memory effects. The Markovian approximation (memorylessness) is applied to the fast modes. During t 1 and t 3 , these cause homogeneous line broadening, while during t 2 they induce population relaxation.
In the limit of decoupled populations and coherences (secular approximation of the Green's function) the population relaxation is described by the Pauli master equation:ρ with K the rate matrix with elements K e e ,ee depicting the population transfer rate from state e into state e which is taken to be independent of the fluctuations of the slow modes. The solution of the differential equation is formally given by the population Green's function ρ e e (t) = e G e e ,ee (t)ρ ee (0) and the elements of the matrix G act as time-dependent weighting factors in the population pathways describing ESA and SE (see eq. 6). The slow bath modes are responsible for spectral diffusion during all three intervals t 1 , t 2 , and t 3 , causing correlations during the three intervals. The secular approximation allows to further partition the nonlinear response (eq. 4) into population (e = e during delay time t 2 ) and coherence contributions (e = e during delay time t 2 ). The population contributions then read: The coherence contributions read: ii e ge (0, t 1 +t 2 , t 1 +t 2 +t 3 , t 1 ) .
As there is only a single GS, no coherence contribution is observed for the GSB. In eqs. 6-7 ε i (i ∈ {g, e, f }) is the electronic contribution to the energy of the i-th ES (the energy of the GS ε g = 0 is used as a reference). ϕ(t 4 , t 3 , t 2 , t 1 ) are phase functions that describe the coupling to the bath and, hence, give rise to the time-dependent spectral line shapes [Abr07, Abr09].
where ξ e e = 1 − δ e e . The phase functions of the population contributions to the ESA and SE signals (ϕ ESA,i and ϕ SE,i ) have both a coherent (f cbe e ) part, the former comprising population-conserving Liouville pathways, the latter population-transfer pathways with the Markovian approximation for the fast modes. Both functions include correlations between t 1 , t 2 and t 3 . The coherent function f cba reads: cbe e reads: In the following we discard memory effects during the delay time t 1 in the incoherent part (all terms g be (t) and g ce (t) terms vanish). This implies that the dynamics initiated on the potential energy surface of state e (populated during the ES decay process) does not keep memory of the active modes in previously visited states. Hence, the dynamics during t 2 and t 3 are governed only by the Hamiltonian of the e state while the line shape during the delay time t 1 is dictated by the term g ee (t 1 ). We note that while this approximation affects the line shapes of the 2D spectra it does not affect the PP spectra where t 1 = 0, so that g be (t) and g ce (t) vanish by construction. g ij (t) in eqs. 9-10 is the line shape function, which is an integral transformation of the autocorrelation function of bath fluctuations C ij (t) [Muk95, Muk09]: In the case when the system is coupled to a set of high-frequency (intramolecular) modes and a continuum of low-frequency (intra-and intermolecular) modes the spectral density can be partitioned into: a) a weakly undamped or underdamped contribution due to the slowly decaying correlation function due to high-frequency modes, responsible for vibrational progressions in the spectra; b) an overdamped contribution due to the fast-decaying correlation function associated with the low-frequency modes responsible for the homogeneous broadening of the signals. The total line shape function then contains two parts as well [Val16]. The semi-classical overdamped regime is usually represented by the line shape function of the semi-classical Brownian oscillator (OBO) in the high-temperature limit [Muk95, Li94, But12]: In eq. 12 λ ij and Λ −1 are respectively the system-bath coupling strength and the fluctuation time scale. Undamped vibrations are described through the line shape function of the multidimensional uncoupled displaced harmonic oscillator (DHO)[Muk95]: with ω k the frequency andd ik the displacement of the i-th electronic potential along the k-th mass-weighted normal mode ( Figure S1). The state-specific displacementsd ik determine the magnitude of diagonal correlation functions C ii (t) and are related to spectroscopic parameters like the Huang-Rhys factors S ik or the reorganization energy λ ik through the following relations: 12 and can be used to compute the spectral densities either in the undamped where γ is a damping factor. Cross-correlation functions C ij (t) depend on the displacementsd ik andd jk associated with two transitions. Because negative displacements (with respect to the GS, Figure S1) are allowed, the corresponding spectral densities C ij (w) that can be formulated following eqs. 14 and 15 can have negative contributions[Nem10] (see next section).

Obtaining simulation parameters for Pyrene
The simulation of the pump-probe spectra for delay times t 2 up to 1000 fs within the framework described in the previous section requires the computation of a number of parameters to feed in the line shape functions of the overdamped Brownian and the undamped displaced harmonic oscillator models (eqs. 12-13). These include the strength λ ij and time scale Λ −1 of the system-bath coupling (12), as well as the normal modes, frequencies ω k and displacementsd ik (13). The computed line shape functions construct the phase functions (eq. 8). To obtain the nonlinear responses for the different Liouville pathways (eqs. 6-7) one has to further compute vertical transition energies ε i and dipole moments µ ij . Finally, the population dynamics must be incorporated by solving the differential equation 5. In the following we elaborate on how the individual parameters were obtained.

Solving the Pauli master equation: population dynamics
The system under study presents only one bright state under the envelope of the pump pulse (250 − 270 nm), S 4 . This is the second bright ππ * state of pyrene and the second state of B 3u symmetry in the FC region at the SS-CASPT2 / RASSCF(4, 8|0, 0|4, 16) / ANO-L[321,21] level within D 2h symmetry, hence labeled also 2B 3u . Due to the absence of further bright states under the pump pulse envelope the coherence contributions to the nonlinear response R k I ,SE,ii (eq. 7) in which the system is in an interstate coherence during the delay time t 2 vanish, leaving only the three population contributions (6). Utilizing the outcome of the analysis of the potential energy surfaces, molecular dynamics simulations, the current experimental data 13 2B 3u Figure S2: Population dynamics. and the findings of previous transient absorption experiments we elaborated a sequential model that describes the population decay following excitation of S 4 (2B 3u ) ( Figure S3). The computation at the FC (Frank-Condon) point shows the presence of another two electronic states in the vicinity of the bright 2B 3u state, denoted as 1B 1g (approx. 0.25 eV below 2B 3u ) and 2B 1g (approx. 0.40 eV above 2B 3u ), both spectroscopically dark and not directly accessible from the GS. Adiabatic molecular dynamics simulation performed in state 2B 3u demonstrate that the energy gaps are conserved throughout the simulation (see Figure S4). Due to the lack of SE signatures in the experimental PP spectrum and the proximity of the dark states (in particular 1B 1g ), we assume that the bright 2B 3u state is left before probing is initiated (i.e. earlier than 200 fs). Tentatively, we assign a lifetime of 100 fs to population transfer to the 1B 1g state. The 1B 1g state is characterized with a three-state conical intersection with the lowest bright S 2 state (labelled 1B 2u ) and the dark S 1 state (labelled 1B 3u ), ca. 0.1 eV above the local minimum, suggesting the potential splitting of the population in S 2 and S 1 . Previous experiments pumping in the near-UV and selectively exciting the S 2 state have revealed an efficient sub-100 fs relaxational dynamics to the S 1 state. Thus, we expect a rapid depopulation of the S 2 state. On the other hand, the dark S 1 state acts as a trap for the excited state population. Thus, we assign the 750 fs lifetime obtained from the global analysis to the effective population of S 1 . Within the time window of the simulation (1 ps) the population does not leave the S 1 state. Furthermore, excited state cooling is found to occur with 8 ps, thus the simulations reveal only "hot" population dynamics. Equation 16 shows the rate matrix corresponding to this model where the rates k 2B 3u = 1/100 fs −1 , k 1B 1g = 1/750 fs −1 were taken as the reciprocal lifetimes. The rate of the depopulation of the dark trapping 1B 3u state is set to k 1B 3u = 1/∞ fs −1 . Solving the Pauli master equation (eq. 5) 14 we obtain following expressions for the population dynamics: (17) Figure S3 shows the evolution of the population between 0 and 1 ps delay time starting in the 2B 3u state (ρ 2B 3u = 1). The population of the dark intermediate 1B 1g state ρ 1B 1g shows a peak around 230 fs. After 1 ps 70% of the population is the "hot" 1B 3u state. The density-matrix element ρ 2B 3u (t) k I ,SE,i (eq. 6) and causes the exponential decay of the signals associated with the 2B 3u state. Due to the vanishing transition dipole moment of the GS← 1B 1g transition and the lack of intense ESA signals in the probed spectral window the 1B 1g state can be seen as a truly "phantom" state which does not show in the spectra. Thus, population in the 1B 1g state does not give rise to any signal in the simulated spectra. The density-matix element ρ 1B 3u (t) multiply the incoherent part f (I) cbe e of the population contributions and lead to the build-up of the 1B 3u associated signals at times close to 1 ps. We note that due to the vanishing transition dipole moment of the GS← 1B 3u transition there is no incoherent contribution to the SE.

Electronic structure computations: constructing the displaced harmonic oscillator (DHO) model
The uncoupled DHO model ( Figure S1) neglects inharmonicities along the normal modes. Damping of the vibrations and internal vibrational redistribution are neglected as they are not expected to show in the simulated time window (0-1000 fs). Within this framework the one-dimensional potential for the i-th state along each normal mode is fully characterized by two parameters, frequency ω k and relative displacement with respect to the GS equilibriumd ik . Parameters are required for: a) the 2B 3u (S 4 at FC) and 1B 3u (S 1 at FC) states involved in the population dynamics (we remind that the intermediately populated 1B 1g state has no spectral signatures in the probed window); b) the higher lying ESs probed out of the S 4 state at early times (< 200 fs); c) the higher lying ESs probed out of the S 1 state at later times. In particular, the parameters for the S 4 state and of the ESs accessed during probing were obtained using a dynamics-based protocol (see sec.

Spectral densities from molecular dynamics
One way to obtain the frequencies ω k and displacements d ik required to feed the DHO model is by fitting the electronic gaps E e (t) − E g (t) (giving rise to GSB and SE signals) and E f (t) − E e (t) (giving rise to ESAs) computed along a single quantum-classical molecular dynamics trajectory without initial kinetic energy[Nen15], to the analytical expression for the classical timedependent fluctuation of the 1B 2u -S n energy gap (n belonging to a state from either the GS, g, or the f -manifold).
Eq. 18 allows extracting the mass-weighted displacement coefficients d ek and d f k , as well as the electronic transition energies ω eg and ω f g (the parameters d gk and ω gg describing the ground state are zero per definition). In the adopted DHO framework the spectral dynamics of the ESA is a function of the dynamics in the photoactive state, i.e. of the relative displacement of the higher lying excited states PES with respect to the PES of the photoactive state, and can induce positive or negative frequency correlations.
[Nem10] By definition the dynamics is restricted to the normal modes describing the molecular dynamics in the e-manifold (i.e. if d ek is zero, then so is d f k ). This is an implication of the missing state-specific modes. Furthermore, we emphasize that the extraction of the energy fluctuations of states from the f -manifold must be performed in a diabatic representation (i.e. by following the excited state wavefunction rather than the adiabatic root) in order to stay within the framework of the uncoupled harmonic oscillators as due to the high state density in the deep UV higher lying states intersect, thus rendering the resulting adiabatic potentials highly anharmonic (see Fig. S4). Although the derivation of the working equations is based on the Franck-Condon (FC) approximation, whereby transition dipole moments are independent of the nuclear coordinates, their fluctuations along the dynamics are accounted for in a semi-classical fashion by using the FC values for the first and second interaction with the pump pulse, while taking the values at time T for the interaction with the probe pulse and the local oscillator.

Spectral densities using gradient projection and normal mode
analysis.
An expression for the displacementsd ik can be obtained from the expression for the energy gap ∆E ij (q k ) between states i and j in the displaced Harmonic oscillator model.
expressed in mass-weighted normal mode coordinatesq. ε i − ε j is the adiabatic excitation energy between states i and j. When computed at the equilibrium of state j the gradient of the energy of state i equates the gradient of the energy gap where P is a matrix whose columns are the normal modes of the system ξ k expressed in normalized mass-weighted Cartesian coordinates Q, M is a diagonal matrix with the nuclear masses and f i is the Cartesian energy gradient. Taking the derivative with respect toq in eq. 19 and inserting the expression derived in eq. 20 gives [Lee16, Fer12] whered i andd j are k-dimensional arrays of the displacements along all modes on states i and j. An alternative to the gradient projection technique is the normal mode analysis [Kur01]. Starting with the definition of the normal mode matrix (eq. 20) one can express it in terms of finite differences as Thereby, ∆q becomes the array of displacementsd i along the normal modes of the system when the difference is taken with respect to the reference point on state i. Rearranging gives the final working equatioñ which allows to estimate how much every normal mode ξ k has to be displaced to connect two points in Cartesian coordinate space. In eq. 24 ∆Q is the difference in Cartesian coordinates between two geometries. With the normal mode analysis attention should be paid to global translational and rotational degrees of freedom which have to be removed prior to computing the difference in Cartesian coordinates. To this end we followed an iterative procedure relying on the vectors of inertia in order to minimize the distance in space between two geometries as outlined in ref.
[Kur01b]. Furthermore, it is paramount that both geometries used to compute ∆Q have been optimized at the same level of theory in order to avoid spurious contributions to the spectral densities

Electronic structure of the S 4 state and its ESA.
In order to simulate the third order nonlinear signal associated with the dynamics in the second bright state S 4 probed in the deep UV, knowledge of following quantities is required: i) the electronic structure of pyrene up to 80 kcm −1 (ca. 10 eV), thus covering not only the S 4 state (labeled e) but also the manifold of higher lying states (labeled f ) that fall within the envelope of the probe pulse (i.e. 35.5-40.0 kcm −1 above the S 4 state); ii) the Franck-Condon active modes of the S 4 state, source of the vibrational dynamics and reflected in a checkerboard pattern and coherent intensity beats in the 2D spectra; iii) the response of the states from the f -manifold to the vibrational dynamics in the S 4 state. ESA bands are an indirect probe of the coherent vibrational dynamics in the photoactive state manifesting in oscillations of signals intensities and/or spectral positions as a function of the delay time.
The frequency ω k and displacement d ik parameters used to construct the lineshape functions g ij (eq. 13) were extracted from a 250 fs adiabatic mixed quantum-classical molecular dynamics simulation treating electrons quantum-mechanically by solving the time-dependent Schroedinger equation, while applying classical Newtonian dynamics to the nuclei. The dynamics was initiated in the S 4 state at the Franck-Condon point without initial kinetic energy. The higher excited state manifold f was computed at the SS-CASPT2/SA-100-RASSCF(4, 8|0, 0|4, 8) level every 2 fs. At the Franck-Condon point we selected the states from the f manifold within the probed spectral window and with significant oscillator strength out of the S 4 state (i.e. the reference diabatic states, shown in green, blue, cyan, brown and magenta in Figure S4) and tracked the temporal evolution of the wavefunction associated with each reference state along the dynamics by searching for the adiabatic state with the greatest overlap with the reference. The energy profiles extracted for the five bright transitions are depicted in Figure S5 showing clear oscillatory dynamics properly captured by the fitting procedure (black dotted lines).
The spectral density of the S 4 state ( Figure S6,f) is dominated by the stretching of the central C-C bond (C 4 -C 10 ) with 1457 cm −1 frequency (see Figure S7). Several other C-C stretching and H-bending modes in the highfrequency region (> 1000 cm −1 ) are found to contribute as well, so as a 592 cm −1 breathing mode in the low-frequency region (see Figure S7). We note that the 1457 cm −1 mode is responsible for the checkerboard pattern in the 2D spectrum while the low-frequency modes give the fine-structure (visible at low temperatures or in crystals) and are one of the reasons for the broadening of the signals at room temperature. To assess the reliability of the dynamics simulation we computed the spectral  Table S5) shows a good agreement between all three methods, the difference between the dynamics-and projection-based protocols arising due to the limited resolution in frequency of the dynamics-based protocol (ca. 60 cm −1 ) due to the limited duration of the simulation. Figure S6,a-e show the spectral densities associated with the bright state from the f -manifold of pyrene under the probe envelope. Note, that asd j on the RHS of eq. 21 is different from zero (being the array of positive defined S 4 -specific displacementsd 2B 3u ) the displacements of the higher lying states d f can assume both positive and negative values (see also Figure S1). Specif- Figure S5: E f − E e (a-e) and E e − E g (f) energy gap dynamics and transition dipole moment fluctuations (gray) of the bright states D1-D5 from figure S4 along the 0K trajectory run in the S 4 state. Fits of the energy gaps following eq. 18 are shown in black. ically, in Figure S6 the spectral densities are constructed using displacements d f relative to the minimum S 4 min in order to highlight the differences between them. Table S4 instead lists the displacementsd f for all states with respect to the FC point as they enter in the line shape function eq. 9.
3.2.4 Electronic structure of the S 1 state and its ESA.
The FC point is the doorway for the excited wavepacket on the bright state S 4 . Its energy gradient or relative position in coordinate space with respect to the ES local minimum allows to extract the active modes governing the vibrational dynamics. Similarly, the three state CI(1B 1g /1B 2u /1B 3u ) is the doorway for the wavepacket on the 1B 3u state. Its analysis can help to extract the 1B 3u -specific active modes. The spectral density of the S 1 state ( Figure S8,a) is dominated by the C-C stretching mode with 1668 cm −1 frequency (see Figure S8c). Several H-bending modes in the high-frequency region (∼ 1000 cm −1 ) are found to contribute as well, together with two asymmetric backbone bending modes with frequencies 350 cm −1 and 539 cm −1 in the low-frequency region (see Figure S8). The spectral density itself was computed via normal mode analysis (eq. 24). We resort to normal mode analysis in this case as gradients generally show larger fluctuations in the crossing region due to strong wave function mixing. ∆Q was calculated as the difference between the geometries of the energy minimum on the S 1 state S 1 min and CI(1B 1g /1B 2u /1B 3u ) as The electronic spectrum at the S 1 minimum was calculated by averaging over 20 states in the A g and B 2g irreducible representations, respectively (i.e. SS-RASPT2/SA-20-RASSCF(4, 8|0, 0|4, 8)/ANO-L[321, 21]). Among all states only one state per symmetry (state 17 in symmetry A g and state 15 in symmetry B 2g ) fulfilled the requirements of having a large transition dipole moment from S 1 and absorbing in the 35.5-40.0 kcm −1 spectral window. Figures S8,b shows the spectral densities of these states obtained via the gradient projection method. The spectral densities are constructed using displacementsd f relative to the minimum S 1 min in order to appreciate the differences between them. Table S6 lists the displacementsd f for all states with respect to CI(1B 1g /1B 2u /1B 3u ).

Parameters for the overdamped Brownian oscillator (OBO) model
The parameters for the OBO model (eq. 12) were chosen so to reproduce the bandwidth of the linear absorption spectrum at room temperature (see Figure 1 in the main text) where inhomogeneous broadening contribution was neglected. Best fit was obtained at coupling strength λ = 300 cm −1 and time scale Λ = 1/50 fs −1 . For simplification, the same OBO parameters were used for all states.

Accounting for realistic pulses
Finite pulse duration is accounted for ad-hoc by convoluting the PP signal W (3) (t 2 , Ω 3 ) with a Gaussian function in the time domain A standard deviation σ = 3 fs was used (matching the experimental fullwidth-half-maximum of 6 fs). Figure S9,a-b show the PP signal prior and after the convolution. Despite the broadening of the bands, the high-frequency intensity beats (ca. 20 fs period) are still clearly visible. Furthermore, the signal was manipulated by applying a Gaussian mask in the frequency domain centered at 37.5 kcm −1 with a FWHM= 5 kcm −1 in order to account for the finite bandwidth of the probe pulse (see Figure 1 of the main paper). Figure S9,c shows the PP signal after applying the mask. Contributions below 250 nm (40 kcm −1 ) and above 280 nm (35 kcm −1 ) are cut off.

2D spectra
In the following we discuss in detail the simulated two-dimensional spectra. Note that each set of spectra in Figures S10-S13 has been normalized and absolute intensities should not be compared between figures.  . The off-diagonal contributions arise due to pump and probe pulse addressing a different vibrational level in the excited state. Weak SE contribution below 36.0 kcm −1 is also recognizable in the first 100 fs (labeled in red). It arises from de-exciting into a higher vibrational GS level (i.e. v0' → v1) with the probe pulse (see also Figure S11). At early times weak ESA contributions can be noticed above 38.0 kcm −1 (labeled in blue, see also Figure S11 and S12). After 200 fs the spectra begin to exhibit ESA features of the 1B 3u state populated with a time-constant of 750fs (see Figure S3). The "dark" 1B 3u state (no SE) is characterized by two intense ESA contributions in the probed window, a signal around 36.0-37.0 kcm −1 and a signal around 38.0-40.0 kcm −1 (see also Figure S13). These become clearly visible at later times (900fs and 1000fs). They overlap spectrally with the GSB signal and modulate its dynamics. Figure S11 shows the spectral dynamics during the very first 20fs after pumping. The majority of the population is in the second bright state S 4 and all signals arise from this state. The 20 fs temporal window was selected as the driving vibrational mode in the S 4 state is the 1456 cm −1 C-C stretching with a period of 23 fs. By inspecting the individual snapshots one can clearly see the vibrational dynamics of the GSB and SE contributions, with off-diagonal terms (e.g. 36.5/38.0 cm −1 due to pumping the v 0 → v 0 transition and probing the v 0 → v 1 transition), thus creating a checkerboard pattern. SE is clearly distinguished between 4 fs and 12 fs. The ESA contribution above 38 kcm −1 originates from several superimposed states exhibiting different oscillation patterns of energy and TDM (see Figure S5). Despite being weak the ESA modulates the dynamics of the GSB, in particular the overtones (probe at 38.5 kcm −1 and 40 kcm −1 ) which are partially or completely canceled. The ESA contributions to the spectrum during the first 20fs after pumping can be more clearly highlighted when working with cross-polarized pump and probe pulses, the reason being that all higher lying states contributing to the ESA have their transition dipole moments oriented orthogonally with respect to the GS → S 4 transition. As shown in Figure S12 with crosspolarized pulses the GSB/SE contribution is nearly completely covered by the ESA and only the most intense v0 → v0' transition remains visible.   Figure S13 shows the ESA contribution to the two-dimensional electronic spectrum close to 1 ps. By this time about 70% of the population has relaxed to the lowest dark 1B 3u state and there is virtually no population left in the initially pumped S 4 state. Thus, the observed signals arise exclusively out of the 1B 3u state. It is noticeable that the ESA spans over the entire spectral range and modulates the spectral dynamics of the GSB (see 500fs to 1000fs in Figure S10). Two intense signals form the ESA spectrum, one around 36.0-37.0 kcm −1 and another one between 38.0-40.0 kcm −1 . Both signals show a vibrational sub-structure giving rise to a time-dependent multi-peak pattern.  a root labeling follows the energetic order at the FC point in the A g irreducible representation in C S symmetry, accounting for both the electronic contribution f and the state-specific reorganization energy (eq. 14); b ε e/f is the electronic contribution to the energy of the i-th ES (GS energy ε g is set to zero); c spectrum origin (adiabatic transition energy) shifted to match the experimental value; d same displacement as for state S 4 due to instability of the gradient projection method. Table S2. Comparison of the parameters for simulating the signals associated with the S 4 state (GSB, SE) obtained through three different approaches: a) from single point computations at the SS-CASPT2/SA-100-RASSCF(4, 8|0, 0|4, 8) level along a 0K trajectory run at the SA-3-CASSCF(4,4) level; b) using gradient projection on the normal modes of the system obtained from a frequency computation at the equilibrium; c) using a projection of the difference between GS and S 4 equilibrium geometries in Cartesian coordinates on the normal modes obtained from a frequency computation at the equilibrium.

Spectral dynamics in the first 20 fs after pumping
Normal mode (MD) S 4 (MD) Normal mode (QM freq.) S 4 (gradient proj.) S 4 (Cart. coord. proj.) a ε e/f is the electronic contribution to the energy of the i-th ES (GS energy ε g is set to zero); b spectrum origin (adiabatic transition energy) shifted to match the experimental value. Table S3. Parameters for simulating the signals associated with the 1B 3u state (only ESA). Parameters for 1B 3u obtained through projection of the difference between the geometries of the 1B 3u minimum and the three-state CI(1B 1g /1B 2u ,1B 3u ) in Cartesian coordinates on the normal modes obtained from a frequency computation at the GS equilibrium. Parameters for the higher lying bright states obtained through projection of the excited state gradient, obtained at the SA-30-RASSCF(4, 8|0, 0|4, 8) level in D 2h symmetry, on the normal modes. 1668d 62 0.688 0.650 0.598 a root labeling follows the energetic order at the 1B 3u minimum in the A g irreducible representation in D 2h symmetry, accounting for both the electronic contribution f and the state-specific reorganization energy (eq. 14); b root labeling follows the energetic order at the 1B 3u minimum in the B 1g irreducible representation in D 2h symmetry, accounting for both the electronic contribution f and the state-specific reorganization energy (eq. 14); c ε e/f is the electronic contribution to the energy of the i-th ES (GS energy ε g is set to zero). Figure S14: Bond length changes during the optimization of conical intersections between various electronic states.

Rationalization of the structure of conical intersections
The following rationalization relies on state-specific orbital occupations and nodal structure analysis to predict geometrical (bond length) deformations that would lead to a conical intersection between electronic states. As seen in Table S5 at the equilibrium geometry of the second bright state S 4 (label 2B 3u ) it exhibits finite gaps with the lower states (0.25 eV to the dark 1B 1g state and 0.66 eV to the first bright S 2 state). We ask ourselves whether there exists an energetically low-lying crossing between both bright states that would avoid the intermediate population of the dark state. looking at the node pattern of the orbitals (Figures S15-S18) it becomes evident that this can be accomplished by elongating C 3 -C 4 , C 3 -C 7 , C 1 -C 2 and shortening C 4 -C 10 , C 7 -C 8 , C 2 -C 3 , as these deformations decrease the bonding and increase the anti-bonding interactions in the H-1 and L orbitals, thus effectively increasing their energies while simultaneously decreasing the antibonding and increasing the bonding interactions in the H and L+1 orbitals, thus effectively decreasing their energies. Performing a CI optimization (at dynamically correlated SS-CASPT2 level), we indeed observe the expected deformations ( Figure S14). However, due to the large energy gap (0.66 eV) it appears that very large deformations, associated with a insurmountable barrier, are required to reach a point of degeneracy. As seen in Table S6 at the equilibrium geometry of the dark state 1B 1g it exhibits finite gaps with the lower lying states (0.27 eV to the first bright state S 2 (1B 2u ) and 0.60 eV to the first dark state S 1 (1B 3u )). We ask ourselves what are geometrical deformations required to reach a CI with these states. The wavefunctions of the 1B 1g state are defined by the configuration state functions (CSFs) H→L+2 (0.73) and H-2→L (-0.30). The wavefunction of the first bright state S 2 (label 1B 2u ) is dominated by the H→L CSF (0.82), while that of the first dark state S 1 (label 1B 3u ) is a linear combination of the H→L+1 (0.59) and H-1→L (0.55). Comparing the electronic structure of the 1B 1g and S 2 state, it appears that in the 1B 2u state there is more electron density in the L orbital (occupation 0.91 vs. 0.32), while in the 1B 1g state there is more electron density in the L+2 orbital (occupation 0.80 vs. 0.09). Thus, any geometrical deformation that destabilizes L would destabilize 1B 2u more than 1B 1g . Conversely, any geometrical deformation that stabilizes L+2 would stabilize 1B 1g more than 1B 2u . By looking at the node pattern of the orbitals (Figures S15-S18) becomes evident that this can be accomplished by elongating C 4 -C 10 and C 2 -C 3 . Performing a CI optimization (at dynamically correlation SS-CASPT2 level), we indeed observe an increase of the C 4 -C 10 bond from 1.46Åto 1.51 Å ( Figure S14) which leads to the decrease of the energy gap. On the other hand, the electronic structure of the 1B 1g and S 1 (label 1B 3u ) offers more flexibility for geometrical deformations. Comparing the electronic structure of the 1B 1g and S 1 state, it appears that in the 1B 3u state there is more electron density in the H and L+1 orbitals (occupation 1.42/0.54 vs. 1.15/0.10), while in the 1B 1g state there is more electron density in the H-1 and L+2 orbital (occupation 1.89/0.80 vs. 1.51/0.08). Thus, any geometrical deformation that destabilizes H and L+1 would destabilize 1B 3u more than 1B 1g . Conversely, any geometrical deformation that stabilizes H-1 and L+2 would stabilize 1B 1g more than 1B 3u . By looking at the node pattern of the orbitals (Figures S15-S18) it becomes evident that this can be accomplished by elongating C 4 -C 10 , C 7 -C 8 and C 2 -C 3 , while shortening C 3 -C 4 , C 3 -C 7 , C 1 -C 2 . Thus, we have more geometrical parameters available to reach a CI between 1B 1g and S 1 . Indeed, a CI optimization circumvents the gap of 0.60 eV between both states, locating a CI just 0.1 eV above the 1B 1g minimum. As the elongation along C 4 -C 10 and C 2 -C 3 destabilize (albeit to a less extent) the S 2 , the found geometry turns out to be a three-state CI with 1B 1g , 1B 2u and 1B 3u lying within 0.07 eV.
7 Active space orbitals  9 Decay associated spectra from experimental data Global fitting of the experimental TA spectra show the presence of complex dynamics. Two different experiments were conducted on the sample sample system and with the same experimental conditions. A long scan study up to 200 ps with 500 fs steps, and a short scan study up to 10 ps with 20 fs steps. The global fitting was done in an iterative way by first letting the fit converge, then optimizing the results by fixing some of the life-times of one of the two sets of data, based on the results of the other. The short scan shows three time-scales, a 0.75 ps lifetime (blue line) is associated with the initial increase of the intensity of the ESA band (positive contributions around 255 nm) accompanied by the decrease of the bleach in the 272 nm region. It is followed by a vibrational cooling processes in the ES described with two time constants, 8 ps (green line), and 24 ps (red line) associated with a blue-shift of the ESA contribution and a decrease of the ESA intensity, respectively. The long scan study also showed three life-times, and the middle and longest life-times of the short scan study overlap well, within experimental errors with the two shortest life-times of the long scan study. The longest of the life-times was purposely fixed to infinity, as a residual signal is always present.     between 500 fs and 700 fs. Therefore the sensitivity of the frequencies is of approx. 167 cm −1 . The 333.6 cm −1 frequency was not considered as might be rising from some laboratory noise or other artifacts. The frequency of 667.1 cm −1 can be assigned to the 597 cm −1 detected in the FT of the residuals of the TA and HTG spectra, while the 1334 cm −1 and 1501 cm −1 , respectively to the 1238 cm −1 and 1413 cm −1 , although a mixing of the two is also possible, due to the relatively low resolution of the data presented here, which would explain the difference in the ratio of the intensities. The last frequency assigned is the 3002 cm −1 , which is clearly corresponding to the 2942 cm −1 detected in the FT of the residuals of the TA and HTG spectra. The frequency of 1001 cm −1 is present only in the peaks 3 and 4 and is currently under investigation.