Spectral lineshapes in nonlinear electronic spectroscopy

We outline a computational approach for nonlinear electronic spectra, which accounts for the electronic energy fluctuations due to nuclear degrees of freedom and explicitly incorporates the fluctuations of higher excited states, induced by the dynamics in the photoactive state(s). This approach is based on mixed quantum-classical dynamics simulations. Tedious averaging over multiple trajectories is avoided by employing the linearly displaced Brownian harmonic oscillator to model the correlation functions. The present strategy couples accurate computations of the high-lying excited state manifold with dynamics simulations. The application is made to the two-dimensional electronic spectra of pyrene, a polycyclic aromatic hydrocarbon characterized by an ultrafast (few tens of femtoseconds) decay from the bright S 2 state to the dark S 1 state. The spectra for waiting times t 2 = 0 and t 2 = 1 ps demonstrate the ability of this approach to model electronic state fluctuations and realistic lineshapes. Comparison with experimental spectra [Krebs et al. , New Journal of Physics , 2013, 15 , 085016] shows excellent agreement and allows us to unambiguously assign the excited state absorption features.


Introduction
When studying ultrafast excited state (ES) dynamics by transient spectroscopy techniques one must consider dephasing processes which occur on a time scale comparable to the delay time between the pulses.The recorded spectra provide an averaged picture of the dynamics during the measurement (which can last several hundreds of femtoseconds).Nevertheless, the high temporal resolution (down to 5-10 femtoseconds) stems from the fact that different segments of the dynamics are averaged over depending on the delay times between the laser pulses which interact with the sample.[3][4] Another essential requirement for modeling nonlinear spectra of photoexcited molecules is the incorporation of higher lying ES.This is desirable for two reasons: first, absorptions to these states (excited state absorptions, ESA) dominate the spectra due to the high spectral density in the visible (Vis) and ultraviolet (UV) region; [5][6][7] second, ESA, together with the stimulated emission (SE) represent characteristic signatures of each photophysical or photochemical decay channel, and hence, their proper description will allow us to recognize electronic structural changes during the dynamics, as well as to disentangle dynamics occurring via competing channels. 8However, the higher ES also exhibits fluctuations coupled to the dynamics on the states on which the relevant photophysics and photochemistry occur (i.e. the photoactive states).An approach is desired which accounts for these fluctuations in the simulation of the nonlinear spectra.The accurate and complete computation of the higher ES manifold is challenging by itself and we demonstrated recently that a significant amount of preliminary work (related to calibration of the computational methods), as well as time and resource consuming computations is required to obtain useful results. 5,9,10The high cost of performing excited state computations drew our attention towards exploring approximate methods for describing photoinduced events.
In this report we present an efficient way of incorporating nuclear dynamics in nonlinear electronic spectroscopy within the framework of mixed quantum-classical dynamics.We use a mathematical model that relies on the cumulant expansion, describing the system's nonlinear response to the electric field, and utilizes a single trajectory as an input in combination with a Dipartimento di Chimica G. Ciamician, Universita `di Bologna, V. F. Selmi 2, 40126 Bologna, Italy.E-mail: artur.nenov@unibo.it,marco.garavelli@unibo.it;Tel: +39 0512099495 b Max Born Institute for Nonlinear Optics and Short Pulse Spectroscopy, This journal is © the Owner Societies 2015 state-of-the-art electronic structure computations to describe contributions due to coherent nuclear motion.We discuss different approximation levels with their limitations.We resort to pump-probe and two-dimensional (2D) electronic spectroscopy (third-order nonlinear techniques) with broadband probing in the Vis-to-UV region in order to assess the accuracy of the presented approach.Notably, 2D optical spectroscopy, which is now state of the art in the infrared regime, is becoming accessible in the Vis and UV region [11][12][13][14][15][16][17][18] and has recently been coupled to supercontinuum probing. 19The beauty of supercontinuum probing resides in its ability to resolve multiple ESA transitions and track their evolution over time.This is a crucial requirement for electronic spectroscopy, as we have shown in previous studies that the most characteristic ESA and SE signals can exhibit shifts of several thousands of cm À1 during the dynamics. 7,8ecently, Riedle and co-workers reported the 2D spectrum of pyrene, a polycyclic aromatic hydrocarbon (Fig. 1), obtained by pumping at the frequency of the bright electronic transition to the second excited state S 2 around 32 000 cm À1 and utilizing a supercontinuum probe in the region 16 000-38 000 cm À1 after a waiting time of 1 ps. 19The spectrum is characterized by a vibronic progression of ground state bleach (GSB) peaks and several well-resolved ESA bands.Using the introduced approach we simulate the 2D spectra of pyrene for waiting times t 2 = 0 and t 2 = 1 ps.
2 Single-trajectory based approach for 2D electronic spectroscopy Multiple field-matter interactions create higher order (nonlinear) responses.The experimental set-up for 2D electronic spectroscopy uses three pulses with wavevectors k 1 , k 2 and k 3 , which interact with the system.The third order nonlinear polarization emits a signal along a given phase-matched direction which is superimposed with a fourth pulse known as the local oscillator (LO).The third-order heterodyne signal 20,21 S ð3Þ t 1 ;t 2 ;t 3 ð Þ¼ involves three temporally separated interactions with the electric field E(r,t) = E(t)e ikrÀiwt with a central frequency w and a complex envelope E(t).S (3) (t 1 ,t 2 ,t 3 ) depends parametrically on the delay times t 1 , t 2 and t 3 between the field-matter interactions.The time-dependent signal is Fourier-transformed along t 1 and t 3 (i.e.delay times) to obtain the 2D spectrum for each waiting time t 2 .In the following we work in the semi-impulsive limit, in which the envelopes of the pulses are approximated by a d-function.In this limit the emitted signal becomes proportional to the nonlinear response 22 where r(0) is the equilibrium density matrix of the system on the basis of the nuclear amplitudes w i (t,Q), i.e. the expansion coefficients of the wavefunction C(t,Q,r) of the quantum-mechanical system formulated on the basis of the electronic wavefunctions F i (r;Q).The latter are parametrically dependent on the nuclear coordinates Q and can be formulated either on the adiabatic or on the diabatic basis.For the electronic states F i (r;Q) we consider the GS g, the manifold of singly ESs e, accessible through the pump pulse (pair), and the manifold of higher lying states f.The fieldfree propagation of the density matrix r(t) is governed by the retarded Green's function using the vibronic Hamiltonian T ˆnuc is the nuclear kinetic energy operator.V ij = hF i (r;Q)|H ˆel |F j (r;Q)i is an element of the electronic Hamiltonian matrix which becomes diagonal on the adiabatic electronic wavefunction basis, i.e.V ij = V i d ij , with the adiabatic electronic potential V i being the ith eigenvalue of the electronic Hamilton operator H ˆel at the nuclear postion Q. T ˆij = (2m) À1 d ij r Q is the non-adiabatic coupling term with elements d ij the nonadiabatic coupling vectors, 23 which are absent in the diabatic representation.Y(t) denotes the Heaviside step-function ensuring causality.
Before the interaction with the first pulse in eqn (2) r(t) is in the GS equilibrium (i.e.r(0) = |w g ihw g |).The interaction with the laser pulses enters eqn (2) through the coordinate dependent transition dipole moment m ¼ P ij w i j im ij ðtÞ w j which couples states i and j.The time argument t serves merely to indicate the time of interaction with the laser pulse.Following relations hold: t 1 = 0, t 2 = t 1 + t 1 , t 3 = t 2 + t 2 , and t = t 3 + t 3 .The expression in eqn (2) gives rise to eight different interaction sequences of the incident electric fields with the system, known as Feynman pathways.A suitably positioned detector can be used to select pathways subject to a phase-matching condition.Exemplarily, in the rephasing k I (k LO = Àk 1 + k 2 + k 3 ) phase-matching condition Before we proceed with re-writing eqn (2) within the framework of mixed quantum-classical dynamics, we briefly describe how nuclear and electronic dynamics are treated therein.The quantum system is expressed as a linear combination of adiabatic wavefunctions F i (r;Q) Note that the time-dependent expansion coefficients c i (t) in eqn (6), unlike the nuclear amplitudes w i (t,Q) (eqn (3)), have no explicit dependence on the nuclear coordinates Q.5][26] After inserting the wavefunction (eqn ( 6)) in the time-dependent Schro ¨dinger equation one obtains the equations of motion for the electronic degrees of freedom 27,28 i _ while the nuclear degrees of freedom are described using the classical Newtonian law of motion with F the force of the quantum-mechanical system on the classical bath.: Q in eqn (7) and Q ¨in eqn (8) denote the nuclear velocity and acceleration vectors.In a semi-classical, trajectory based Ansatz, the coupled eqn (7) and (8) are commonly solved applying the independent trajectory approximation, where each trajectory is driven by a single electronic potential (the force F in eqn ( 8) is then simply given by the gradient on a single electronic potential F i (r;Q)).The probability distribution is recovered in a stochastic manner by allowing the trajectories to change the surface on which they are evolving (i.e. a surface hopping event 29 ).Thereby, the decision of working either in the adiabatic or in the diabatic representation is taken based on the availability of non-adiabatic elements (necessary for adiabatic dynamics) or adiabatic-to-diabatic transformation matrices (necessary for diabatic dynamics).
The stochastic nature of surface hopping makes it meaningless within the single trajectory framework.Therefore, we are forced to restrict the dynamics on a single electronic state and to neglect inter-state couplings when relying on a single trajectory.This simplifies the electronic equations of motion (7)  to a set of independent equations Eqn ( 9) is written on a general basis, where E i can correspond either to an eigenvalue of the electronic Hamilton operator H ˆel (adiabatic representation) or to a diagonal element of the electronic Hamiltonian matrix V (diabatic representation).At that point we note that, while the outcome of the mixed quantum-classical dynamics simulations following eqn (7) does not depend on the choice of basis, this is no longer the case in the single trajectory framework utilizing eqn (9).Performing adiabatic dynamics, while neglecting the non-adiabatic coupling vectors d ij , allows for ultrafast changes of the ES wavefunction.On the other hand, diabatic dynamics in the absence of diabatic couplings V ij implies to follow the nature of the ES.In the following we will formulate the equations for the nonlinear response on the general basis.For pratical applications one has to estimate which representation poses a smaller approximation in the description of the ultrafast dynamics in the system of interest.This aspect is further discussed in the Computational methods, relating it to pyrene.Having worked out the basics for mixed-quantum classical dynamics in the single trajectory framework we now proceed to re-formulate eqn (2).We consider a system consisting of a single bright state e, which is vibronically coupled to one or more dark states e 0 , as is the situation in pyrene (see Results).The density matrix r(t) is propagated according to the coupled equations of motion for the electrons and nuclei.With eqn (9)  the propagation of the density matrix when the system is in a coherence state, i.e. r(t) = c i (t)c j * (t), simplifies to a function of the fluctuating electronic gap between the states i and j The choice of a reference Hamiltonian acting on the classical bath (eqn ( 8)) during a coherence propagation, as it is the case for the delay times t 1 and t 3 (Fig. 2), is not unique 2,20 as the bra and the ket side of the density matrix are subject to different electronic potentials.During coherence between the GS g and the singly excited manifold e, as well as between the singly ES manifold e and the higher ES manifold f we choose to propagate the nuclear degrees of freedom as if interacting with the singly excited electronic state e.Subotnik and co-workers proposed recently to use a swarm of trajectories propagated independently in the GS and in the ES, 3 thereby introducing two different reference Hamiltonians for either the bra or ket side and evaluating the coherence as the geometric average.The above approximation allows us to derive a simplified expression for each Feynman pathway (eqn ( 5)).Exemplarily, the response function (t 1 ,t 2 ,t 3 ) for the ESA under the rephasing phase-matching conditions reads R where the infinite ES lifetime (a consequence of eqn ( 9)) is corrected phenomenologically by introducing an effective ES lifetime T e .For the waiting time t 2 the system is in a population state r(t 2 ) = c e (t 2 )c e *(t 2 ), exhibiting coherent population dynamics modulated by T e .Eqn (11) (and its GSB and SE counterparts, This journal is © the Owner Societies 2015 which we do not explicitely formulate here) takes into account the electronic fluctuations due to nuclear motion in the ES (i.e.spectral diffusion) and can be directly employed in the simulation of 2D electronic signals by a 2D Fourier-transformation of the third order response function R (3) (t 1 ,t 2 ,t 3 ) subject to field envelopes with respect to time delays t 1 and t 3 (eqn (1)), followed by averaging over an ensemble of trajectories.On the single trajectory level the Fourier-transformation can introduce dispersive contributions in the spectra if the electronic gap fluctuations induce a shift of the mean transition frequency within the dephasing timescale. 30The dispersive features are (in part) averaged out in the ensemble averaged signal.This effect is particularly severe for large energetic fluctuations and ultrashort time-scale dynamics.
The cumulant expansion provides an exact solution to models with Gaussian statistics, i.e. the electronic degrees of freedom coupled to a harmonic bath.The Ansatz avoids tedious averaging over trajectories by introducing a functional form for the correlation function.By introducing the harmonic approximation, the linearly displaced Brownian harmonic oscillator model (BHO) thus represents an analytical solution of the coherence evolution for two (e.g.GS and ES) multidimensional harmonic oscillators. 3,20It provides an elegant alternative to extract relevant information of the spectral lineshape function g(t) from a single trajectory.The use of the BHO model is generally justified as the FC active modes that control the immediate relaxation after excitation are associated with bond relaxation deformations, whose GS and ES can be described to a good approximation through displaced harmonic potentials (short time approximation).ES specific structural deformations, which introduce anharmonicities in potentials are generally activated on a longer time scale.The coherence dynamics (10)) can be expressed within the BHO model as follows with o ij refers to the adiabatic transition between the GS and the ES (i.e.o eg ), while o k and D k describe the frequencies of the active vibrational modes and their amplitudes (i.e. the Franck-Condon coefficients).Upon integration eqn (13) yields coherent nuclear motion (underdamped limit) where we have omitted the temperature dependence in eqn ( 13 Using eqn (12) the expression for the response function The Condon approximation applies for m ge for the coherence time t 1 , while the coordinate dependence of m ef is considered explicitly through the dependence of the transition dipole moment on the dynamics for the delay time t 1 and the waiting time t 2 .Besides the coherent nuclear motion and the ES lifetime, the solvent induced fluctuations provide additional mechanisms of line broadening.While appearing naturally in the multi-trajectory description, the coherence decay due to The approximation of coherent superposition dynamics for the delay times t 1 and t 3 assumes that the wavepacket superpositions, created by the first pump pulse at time t 1 and the probe pulse at time t 3 , evolve coherently at two electronic potentials and nuclear decoherence effects are neglected.For each combination of delay time t 1 and waiting time t 2 the population (color coded and labeled), generated upon interaction with the second pump pulse at time t 2 , has propagated coherently on the electronic potentials, so that the third pulse probes a different segment of the higher ES manifold.The time-dependent fluctuation of the electronic gap is a function of (and hence contains information about) the active vibrational modes.In the semi-impulsive limit the signal intensity depends on the magnitude of the transition dipole moment at times t 1 , t 2 and t 3 .
This journal is © the Owner Societies 2015 Phys.Chem.Chem.Phys., 2015, 17, 30925--30936 | 30929 these fluctuations must be introduced phenomenologically in the single trajectory framework.When these fluctuations are small or very fast, the fast modulation limit applies, i.e. the lineshape function simplifies to a homogenous Lorentzian characterized by the pure dephasing rate g.The total homogenous dephasing contribution to the response function is thus given by G = g + 1/T e in eqn (15).Eqn (15) is applied in generating spectra based on a single trajectory simulation.It is particularly beneficial for simulating transient spectroscopy in the case of ultrafast dynamics for the delay times t 1 and t 2 , as it captures the spectral fluctuations and intensity modulation of the accessible (i.e.bright) states in the higher ES manifold f, which depend on the dynamics on the photoactive ES e.However, it must be noted that signals arising from other electronic states populated during the coherent dynamics are not incorporated in eqn (15), therefore the transient signal characterizes with an overall intensity decrease for a longer waiting time t 2 as a consequence of the finite ES lifetime.
The nonlinear response functions for GSB and SE can be formulated in a straight-forward way.The GSB reads We note that the present model does not include GSB recovery.The SE reads R In the Results section eqn ( 15)-( 17), together with their counterparts for the non-rephasing phase matching condition k II (k LO = + k 1 À k 2 + k 3 ), are used to simulate the absorptive twodimensional spectrum 20,21 of pyrene, as well as its marginal, the pump-probe spectrum, for waiting time t 2 = 0, thereby resolving the spectral features of the bright second ES S 2 .
Going beyond the restrictions posed by the introduced approach, namely, to provide the spectral lineshapes of ESA and SE associated only with the initially excited state, one faces the challenge of correlating non-equilibrium dynamics on different ESs subject to the inter-state couplings between them.In the special case when the fast initial dynamics in the bright state e due to the departure from the FC region is followed by slow dynamics associated with trapping of the ES population in dark states e 0 the non-linear response for waiting time t 2 longer than the ES lifetime T e of the bright state e will be dominated by the ESA features of the ES intermediates.This allows us to couple short-time ES dynamics, serving to describe the g-e coherence for t 1 , with a static description of the e 0 -f coherences for t 3 .In the static (snapshot) formulation of the nonlinear response it is assumed that the coherence dynamics is slower than the duration of the experiment and, therefore, spectral diffusion can be neglected.Therefore, one can use the energies obtained at the equilibrium geometry of the dark state to probe its ESA.The spectral lineshape at longer times is described solely through the homogenous dephasing rate.In principle population transfer among states could be incorporated either phenomenologically by a rate equation system 22 or by microscopic surface hopping dynamics. 2he coupling of fast and slow dynamics is particularly suitable for disentangling ES deactivation pathways through constructing the fingerprint spectra of the ES intermediates, which could be populated during the decay of the initially excited states.In the Results section we demonstrate how the idea can be applied in practice by coupling the short-time dynamics in the bright second ES of pyrene to the ESA spectrum of the local minimum on the dark first ES.The latter is known to be populated on a ultra-short time scale (T 2 E 80 fs).

Computational methods
The present study integrates different types of computations: ground and excited state optimizations, single point calculations, mixed quantum-classical dynamic simulations and the construction of linear and nonlinear electronic spectra.
The calculations aimed at the characterization of the electronic structure of pyrene have been performed in the framework of the complete active space self-consistent field (CASSCF) approach 31,32 augmented by second order perturbation theory (CASPT2). 33,34We also made use of the restricted active space (RAS) variation of the scheme (i.e.RASSCF/RASPT2), 35 which constructs the list of configuration state functions through a limited number of simultaneously excited electrons out of the occupied orbital set RAS1 and a limited number of simultaneous excitations into the virtual orbital set RAS3, while still allowing for all possible permutations of electrons among the orbitals in the RAS2 set.In the particular application the RAS2 set was left empty.We introduce the active space notation RAS(n,m|0,0|n,m), where n is the maximal number of simultaneously allowed holes(electrons) among the occupied(virtual) orbitals m in RAS1(RAS3).The optimization of the ground state was performed at the CAS(4,4) level of theory, including the highest occupied molecular orbital (HOMO, H), the lowest unoccupied molecular orbital (LUMO, L), H À 1 and L + 1 and the corresponding electrons.The two lowest excited states of the system (the so-called L a and L b states in Platt's notation 36 ) were optimized within D 2h symmetry, with the p-orbitals and electrons of pyrene (see ESI, † Fig. S1) distributed among the RAS1 and RAS3 subspaces and allowing up to quadruple excitations, i.e.RASSCF(4,8|0,0|4,8).The minimum on the dark state L b was subsequently characterized by means of single point calculations by averaging over 60 states (i.e.SA-60-RASSCF(4,8|0,0|4,8)/SS-RASPT2) to obtain the manifold of higher lying excited states accessible from the dark state and, thus, responsible for the visible ESA in the spectrum computed at longer waiting times.Within the RASPT2 calculations an imaginary level-shift correction of 0.2 au was used in order to avoid the presence of intruder states 37

and the core orbitals
This journal is © the Owner Societies 2015 were frozen.The CASPT2 standard zeroth-order Hamiltonian was used. 34For all computations the ANO-L basis set contracted to C[3s,2p,1d]/H[2s1p] was used. 38The Cholesky decomposition was adopted to speed up the evaluation of two-electron integrals. 39inimal active space SA-2-CASSCF(2,2)/6-31G* mixed quantumclassical dynamics simulations were conducted in the L a excited state for 250 fs with a time step of 0.5 fs, thereby employing zero initial starting velocity (i.e. a single 0K-trajectory).The small active space used in the dynamics simulations is necessary due to the inherent inability of CASSCF to describe the dark S 1 (L b ) and the bright S 2 (L a ) state at equal footing, which leads to the unphysical mixing of both states during the dynamics.Including solely the HOMO and the LUMO suppresses the description of the L b state which additionally requires orbitals HOMOÀ1 and LUMO+1, while considering the L a state which is characterized by a nearly pure H -L transition.Consequently, the minimal active space reveals solely the dynamics on the bright state and is blind for the state mixing in the L b /L a crossing region traversed during the simulation (see Results).Effectively, the single trajectory resembles more closely the diabatic picture, in which the nature of the ES is followed.We believe that for the particular study this quasidiabatic description is suitable for two reasons: (a) a crossing region between the bright and the dark states is encountered along the deformations dictated by the FC-active modes (see Results).As in the vicinity of a conical intersection the diabatic couplings are small, unlike the non-adiabatic couplings, the short-time dynamics are better described in the diabatic representation; (b) the effective ES lifetime T e , which we use, was obtained experimentally through fitting the decay of the transient ESA features at 550 nm, 40,41 which are associated with the electronic configuration of the bright (i.e. the diabatic) state.
The ability of the CASSCF(2,2) dynamics to reproduce the main dynamical features of the L a state of pyrene is discussed in detail in the Results section.Due to the overestimation of the energies at CASSCF(2,2) an ad hoc scaling of the electronic gap was performed.A linear constant factor of 0.65 was extracted out of averaging over calculations at 50 auxiliary points along the L a trajectory (i.e.structures at which the trend of the L a energy changes) computed at the SA-60-RAS(4,8|0,0|4,8)/SS-RASPT2/ANO-L(321,21) level.
The electronic structure computations and geometry optimizations were performed using the Molcas 7.9 software, 42 while the quantum-classical dynamics was operated through the velocity Verlet algorithm implemented in the Cobramm software 43 and interfaced using Molcas 7.9.
Employing the computed SS-RASPT2 energies and the SA-RASSCF transition dipole moments (TDMs), absorptive 2D electronic spectra were computed by implementing eqn ( 15)-( 17) within the sum-over-state (SOS) framework 44 in a development version of Spectron 2.7. 22In practice the g-e coherence is modelled using eqn (12), where o k and D k are extracted from the 0K-trajectory using a Fourier series to fit the temporal evolution of the electronic gap E e (t) À E g (0), while the adiabatic transition o eg is set to the electronic gap between the GS equilibrium g and the lowest energy point passed during the first oscillation of the wavepacket propagating in the ES e. Regarding the evaluation of the e-f coherence necessary to evaluate eqn (15) we assume that on an ultrashort time scale the wavepacket promoted with the probe pulse to a particular higher excited state will preserve its wavefunction character (diabatic representation).We selected at the FC point the states from the manifold of higher lying states f with significant oscillator strength out of the L a state within the probed spectral window (i.e. the reference diabatic states) and tracked the temporal evolution of the wavefunction associated with each reference state along the 50 auxiliary points by searching for the adiabatic state with the greatest overlap with the reference, thereby considering that only the CI-coefficients as the molecular orbitals change only slightly.The energy and TDM fluctuation profiles extracted in this way were subsequently fitted with Fourier series allowing us to compute the required quantities for any temporal resolution.We stress here that tracking diabatic states in that way is only an approximate and not always an unambiguous procedure.In the future we aim to interface a diabatization routine to the electronic structure computations, which will allow us to work directly in the diabatic representation, both for performing mixed quantum-classical dynamics simulations and for diabatic state tracking.
The effect of temporally finite pulses is not considered in this treatment, i.e. the spectral simulations are performed in the semi-impulsive limit.GSB and SE contributions appear as red peaks, while ESA appear as blue peaks in the 2D plots.A pure dephasing rate constant g corresponding to 200 cm À1 was used.

Results
As outlined in the Theoretical section, if the dephasing processes in the system occur on a time scale which is comparable to the duration of the experiment, the spectra reveal an averaged picture of the dynamics during the measurement.Therefore, the main photoinduced phenomena caused by the irradiation with the pump and probe laser pulses must be known in order to simulate linear and nonlinear spectroscopy.In the following, we characterize the photophysics and photochemistry of pyrene.
The GS minimum, hereafter denoted as 1 (gs) min , belongs to the D 2h point group.The FC point is characterized by two lowlying singlet excited states.Following Platt's nomenclature, 36 S 1 is the so-called L b state, whose wave function is described by the negative linear combination of the H -L + 1 and H À 1 -L configurations, while S 2 is referred to as the L a state, having a wave function dominated by the H -L configuration.The L b and L a vertical excitation energies at the FC point are equal to 3.42 and 3.77 eV, respectively (see ESI, † Table S1).The electronic transition from the ground state to the L b state is dark within D 2h symmetry, while the L a absorption is characterized with a high oscillator strength of 0.17.This implies that the interaction of the system with near-UV radiation causes the population of solely the L a state, from which consequently the main photoinduced events of the system take place.In order to study the L a photophysics, a SA-CASSCF(2,2)/ 6-31G* mixed quantum-classical dynamics simulation with zero starting velocity was conducted starting at the FC point (Fig. 3, green line).A legitimate question is whether such a limited computational level could describe the essential dynamics on the L a state with a sufficient accuracy.From experimental observations, it is known that after excitation into the L a state an ultrafast population transfer toward the L b state takes place at a time constant of 85 fs. 19As the time is too short to allow for vibrational energy redistribution, a conical intersection (CI) between the L a and the L b state must be encountered along the deformations dictated by the FC active modes.Therefore, the quality of the dynamics can be assessed by analyzing if pyrene explores regions of near degeneracy between the L a and the L b state.Full RASSCF(4,8|0,0|4,8)/RASPT2 computations at 50 auxiliary points along the trajectory, which provide the energies of both L b and L a demonstrate that both states swap order multiple times (marked in Fig. 3 with red vertical lines), with the state crossing being reached for the first time after less than B10 fs.In fact, the region of inverted state order is visited by the molecule ca.every 20 fs (8-10 fs, 26-28 fs, 44-46 fs, etc.).
The minimum on the L a surface in D 2h symmetry is found to be about 0.5 eV beneath the FC point and is characterized by the equilibration of the conjugated bonds in the periphery of pyrene (Fig. 4, 1 (L a ) min , see also ESI, † Table S1).The most pronounced bond variations are around 0.04-0.05Å.The analysis of several geometries in the region of inverted state order (see ESI, † Fig. S2 and S3) demonstrates that the crossings appear within 0.2-0.4eV above 1 (L a ) min , thus being below the FC limit of 3.77 eV.The crossings are characterized by a more pronounced displacement along the very same coordinates that connect the FC point to 1 (L a ) min , with bond variations up to 0.05-0.07Å. Accordingly, the CIs lie closer to the outer turning point along these coordinates.
Based on the above findings we postulate a simple onedimensional model that comprises the ultra fast photophysics of pyrene (Fig. 4).After irradiation the populated bright L a state undergoes a sub-10 fs bond length relaxation.The momentum accumulated in the FC active modes drives the system beyond the L a minimum to a turning point, which lies in the vicinity of the L a /L b crossing region.The wavepacket bifurcates and a part continues its dynamics on the L b surface, while the wavepacket that remains on the L a surface makes one turn in order to return to the crossing region ca.every 20 fs.The efficiency of the transfer could be roughly estimated to be 15% (assuming that within 85 fs the wavepacket reaches the CI region five times, a constant L a -L b population transfer rate and neglecting population transfer back to L a ).It is well known from earlier transient spectroscopy studies with picosecond resolution that the pyrene is trapped in the L b minimum (Fig. 4, 1 (L b ) min , see also ESI, † Table S1) before it eventually decays to the ground state through flourescence 40,45 (L b lifetime 339 ns 41 ).To summarize, the obtained results clearly demonstrate that the CASSCF(2,2) dynamics grasps the essential dynamics along the L a state responsible for reaching the CI region.
We now turn to the electronic spectroscopy of pyrene.First we present the linear absorption (LA) spectrum, which is obtained by a Fourier transformation of the first order response function given by R (1) The electronic energy gap fluctuations due to the lineshape function g(t) are obtained through fitting the dynamics with a Fourier series as outlined in the Computational methods.The homogeneous broadening G = g + 1/T e is affected by two  This journal is © the Owner Societies 2015 contributions, a lifetime broadening with a relaxation rate constant 1/T e = 1/85 fs (corresponds to a 125 cm À1 Lorentzian bandwidth), which accounts for the decay from the bright L a state to the dark L b state, 19,46 and a pure dephasing rate constant g which accounts for signal decay due to solute-solvent interactions.
The fit reveals that the initial relaxation dynamics of pyrene can be described through a few dominant vibrational modes, 400 cm À1 , 1200 cm À1 and 1600 cm À1 modes with Huang-Rhys factors 0.46, 0.44 and 0.41, superimposed with several less intense modes.These modes match well with the ones reported by Freidzon et al. to have the highest L a Huang-Rhys factors in the gas-phase (412 cm À1 (0.26), 1277 cm À1 (0.214), 1456 cm À1 (0.219), and 1700 cm À1 (0.279); Huang-Rhys factors in parentheses) 47 out of ab initio LA spectra generated within the Franck-Condon and Herzberg-Teller approximations.Freidzon et al. also report that similar modes are responsible for the vibronic profile of pyrene in water, based on computations of pyrene clustered by 44 waters, which represent its first solvation shell.
The vibrational modes with the highest Huang-Rhys factors contribute to the vibrational progression, which gives the spectral lineshape of the absorption spectrum.Fig. 5, left panel, compares the theoretical spectrum (red), based on the dynamics simulation in vacuo, and the experimental spectrum (black), obtained in aqueous solution. 48The rate constant g was set to 200 cm À1 to fit the bandwidth of the solvated spectrum.Both spectra show a vibrational progression with three absorption peaks with origin around 30 000 (29 000) cm À1 (i.e. the adiabatic transition energy o eg ).The vibrational bands are separated by B1400 cm À1 , with the 0-0 transition accumulating the highest oscillator strength.Fig. 5, right panel, compares the theoretical spectrum (red) and the experimental spectrum (black) recorded in a Ne-matrix at 6 K. 49 Thereby, the pure dephasing rate constant g was set to zero.Both spectra exhibit a finer vibronic structure, which reveals the progression due to the 400 cm À1 mode, not visible in the LA spectrum in solution.Notably, the theoretical spectrum characterizes with larger weights for the 0-2 overtones and a red-shift of the 0-0 transition of B2200 cm À1 , which is within the accuracy limits of the RASSCF/RASPT2 protocol and is probably due to the low level dynamics and the use of a uniform scaling factor.It should be noted as well that the pyrene absorption bands exhibit a red-shift of 1000 cm À1 in water. 48he LA spectrum of pyrene reveals the FC-active modes that govern the immediate excited state dynamics.However, changes in the electronic structure due to population transfer among excited states remain hidden.Ultrafast transient electronic spectroscopy provides the missing information.Each electronic state possesses characteristic bright transitions which give rise to ESA peaks in the transient spectra and can be used for state recognition.Pyrene is no exception to this rule, the bright L a and the dark L b state have essentially different electronic structures, therefore probing that the initial ultrafast dynamics on the L a state has different signatures than the subsequent dynamics in the L b state. 41,46To demonstrate this we generate the 2DUV spectra of pyrene and their marginals (i.e. the pump-probe spectra), for two different waiting time t 2 .For simulating the so called FC spectrum we assume that the pump (the second pump pulse in the case of 2D spectroscopy) and probe pulse arrive simultaneously at the sample, i.e. t 2 = 0.This spectrum should reveal the ESA transitions associated with the dynamics of the wavepacket on the L a surface.Another set of spectra is generated for a waiting time t 2 on the picosecond time scale and resolves the ESAs of the L b state.
The vertical SA-60-RASSCF(4,8|0,0|4,8)/SS-RASPT2 transitions at 50 auxiliary points along the scaled CAS(2,2) 0Ktrajectory, together with the associated TDMs, provide the required quantities for constructing the FC 2D spectrum of pyrene.We compare two different computational strategies.Assuming (unrealistically) that the time-scale for ES relaxation is much longer compared to the duration of the experiment we apply the snapshot approximation for simulating the coherence evolution for the delay times t 1 and t 3 and model the broadening of all optical transitions with a constant uniform Lorentzian linewidth (Fig. 6, left).On the other hand, the spectral lineshape in Fig. 6, right, incorporates the fluctuations of the electronic gap for the delay times t 1 and t 3 through the use of eqn (15).Prominent transitions are labeled (L a GSB/SE, peaks 1-4).Apart from the obvious inability of the static approach to reproduce the vibronic structure of the GSB, the comparison between both levels of theory clearly demonstrates that the static approximation reproduces the positions of the ESA only qualitatively, thereby neglecting the dynamically induced blue/red-shift and the peak broadening.We draw the reader's attention to two observations: first, while peaks 1, 3 and 4 blue-shift (by up to 4000 cm À1 ) when the fluctuation of the electronic gap is considered, peak 2 exhibits a pronounced red-shift of B2000 cm À1 ; second, peaks 1 and 4 whose wavefunctions are dominated by the negative (respectively positive) linear combination of the H À 4 -L + 1 and H À 1 -L + 4 configurations, both blue-shift, with the negative linear combination, which is dark in the FC region, gaining oscillator strength at the expense of the positive linear combination, which is covered by GSB/SE.The 2D spectrum in Fig. 6 recovers the vibrational structure of the GSB along o 1 and o 3 , with the third vibronic band (i.e. the 0-2 overtone) along o 3 being covered by a close lying intense ESA.Three characteristic ESAs appear in the spectrum at 18 500 cm À1 (peak 2), 23 500 cm À1 (peak 1) and 33 000 cm À1 (peak 4).ESA peak 2 is responsible for the intense absorption, observed in transient pump-probe spectra with high temporal resolution (30-60 fs) 19,46 around 17 000-18 000 cm À1 , and constitutes the most characteristic signature of the S 1 state in the Vis region.A weak absorption band at around 25 000 cm À1 , reported by Krebs et al., 41 matches the position of peak 1.To the best of our knowledge no transient data with sufficient temporal resolution that probes beyond 32 000 cm À1 has been reported in the literature, therefore it appears that intense absorption associated with peak 4 has not been detected yet.
The performed ab initio quantum-chemical calculations allow us to characterize the nature of the electronic excited state giving rise to the fingerprint ESA at 18 000 cm À1 (peak 2, see also Table 1).The dominant configurations can be described in the language of vibrational spectroscopy as the overtone of the GS -L a transition, i.e.H -L H -L + H -L. The doubly excited state appears red-shifted with respect to the sum of the single electron transitions out of the GS, a large anharmonicity of 1 eV is observed.Our computations reveal that the energy of the overtone transition decreases from the FC towards the 1 (L a ) min due to the larger stabilization of the doubly excited state in comparison to the L a state.This unique feature of overtones, documented by us for other systems as well, 8 causes peak 2 to red-shift by nearly 2000 cm À1 when dynamic fluctuations are considered in the simulation of the spectrum and explains the overestimation of the ESA transition in the static spectrum (21 000 cm À1 in Fig. 6, left) compared to the experimental findings (17 000 cm À1 41,46 ).
We now turn our attention to the spectrum at waiting time t 2 in the picosecond range.On the basis of the available literature on pyrene, it is evident that after 1 ps the excited state population has decayed to the dark L b state minimum 1 (L b ) min where it remains trapped on a time scale much longer than the duration of the experiment. 40,41,45Consequently, the appearing ESA signals are caused by electronic promotion out of the 1 (L b ) min .Employing SS-RASPT2 energies and the SA-RASSCF TDMs at the 1 (L b ) min we constructed the 2D spectrum for these conditions by coupling the short-time dynamics to the static description of the dark state intermediate (Fig. 7, left panel).We remind that this approach does not require information about the coherent dynamics for t 2 and should be applied only for generating the spectra of excited state intermediates, which are expected to trap the population for a sufficiently long time.
The simulated spectrum recovers the vibrational structure of the GSB, as well as two distinct ESA peaks at B21 000 cm À1 and B26 000 cm À1 (C, D in Fig. 7).Two less intensive peaks (A, B in Fig. 7) appear at around 18 000 cm À1 .Next, we elucidate the nature of excited states responsible for these ESA signals through analysis of the corresponding RASSCF wave functions.The latter are reported in Table 2. To underline that different states are revealed with respect to the FC spectrum (Fig. 6) we use letter-codes to label the ESA peaks in Fig. 7.The excited states associated with both most intense ESA peaks C and D can be regarded as combination bands of the lowest singly exited states of pyrene (L a , L b , B a , B b in Platt's notation 36 ).The excited state giving rise to peak C can be seen as L b + L a transition, while the state associated with peak D as L b + B a transition.As in the FC spectrum, the transitions exhibit pronounced anharmonicities and appear strongly red-shifted with respect to the Fig. 6 Theoretical quasi-absorptive 2D electronic spectrum of pyrene for a delay time t 2 = 0: (left) applying the snapshot approximation for the coherence evolution for t 1 and t 3 ; (right) applying eqn (15), thus explicitly including spectral diffusion for the coherence times t 1 and t 3 .The marginals of the 2D electronic spectra (i.e. the pump-probe spectra), obtained through integration over O 1 , are depicted on the side.Peak assignment is performed for the most characteristic transitions.The excited states associated with these transitions are given in Table 1.
Table 1 Vertical excitation energies DE (in eV and cm À1 ), transition dipole moments (in a.u.) out of the L a and wavefunctions for excited states responsible for the computed ESA peaks appearing in the spectra in Fig. 6 Label Riedle and co-workers reported recently an experimental 2D electronic spectrum of pyrene in methanol, recorded in a collinear pump-probe set-up utilizing a narrowband pumppulse pair, centered at around 32 000 cm À1 , i.e. at the second and third vibrational bands of the L a absorption and a supercontinuum probe pulse covering the Vis-to-UV spectral window (16 000-38 000 cm À119 ).A time delay t 2 equal to 1 ps was used.Such a result represents one of the few experimental examples employing Vis-UV broad pulses for 2D spectroscopy.The spectrum is shown in Fig. 7, right panel.Through the use of a narrowband pump''pulse pair the fundamental (0-0) was suppressed along o 1 , whereas it was recovered along o 3 thanks to the supercontinuum probing.The ESA peaks do not exhibit a vibrational structure along o 3 , but are noticeably elongated.The short lifetimes of the higher excited states or/and the large FC factors for the fundamental (0-0) transition could rationalize the broad and unstructured lineshapes.Knowledge of the energetic profiles of the higher excited states is necessary to shed more light on the matter.The comparison between theory and experiment shows an overall good agreement, which supports the validity of the applied approximations.On the basis of the ab initio quantum-chemical calculations it is now possible to unambiguously assign the bright ESA transitions C and D, while it can be spectulated that peaks A and B characterize the weak ESA at 20 000 cm À1 .

Conclusion
Bond length relaxations, the initial relaxation process in electronically excited molecular systems, are an immediate consequence of the altered electronic structure upon excitation and almost always occur on an ultrafast time scale.With the emerging laser technologies, approaching the sub-fs pulse duration, [50][51][52][53][54] the signatures of the relaxation dynamics, governed by FC-active modes, become visible.While these modes could be evaluated directly from ab initio computations using the time-(in)dependent Franck-Condon and Herzberg-Teller approximations or by fitting experimentally available high resolution LA spectra, the present approach, based on a mixed quantum-classical dynamics simulations, allows us to simultaneously record the fluctuations of the higher excited states as a function of the dynamics in the photoactive state(s).In that sense the present contribution is an improvement of our previously presented strategy 9 for incorporating higher excited states and is a further step towards the formulation of a computational protocol for the accurate computation of twodimensional electronic spectra.
The present example of the ultrafast decay of pyrene from the bright L a state to the dark L b state probed for waiting times t 2 = 0 and t 2 = 1 ps demonstrates the shortcomings of the static description and the ability of the proposed approach to model ESA fluctuations and realistic lineshapes.Without the explicit consideration of electronic fluctuations due to the nuclear degrees of freedom the vibronic structure and spectral shifts of the peaks are neglected.The presented strategy couples accurate Fig. 7 Theoretical (left) and experimental (right) quasi-absorptive 2D electronic spectra of pyrene for a delay time t 2 in the picosecond range (t 2 = 1 ps in the experiment) obtained through pumping at the frequency of the S 2 transition and supercontinuum probing in the Vis-to-UV region.In the simulation all vibrational bands of the GS -S 2 transitions are excited, while in the experiment the fundamental is suppressed.In the simulation it is assumed that probing occurs from the minimum on the dark S 1 1 (L b ) min , which is populated with a time constant of 85 fs and where the system is trapped on a time scale much longer than the experiment.The vertical transitions along o 3 are obtained from a single point SA-60-RASSCF(4,8|0,0|4,8)/SS-RASPT2 calculation at 1 (L b ) min .Peaks assignment is performed for the most characteristic transitions.The excited states associated with these transitions are given in Table 2.The experimental spectrum is taken from ref. 19.
) (high frequency limit: bh o k /2 4 1 with inverse temperature b = 1/k B T).The coefficients D k are related to the Huang-Rhys factors S k by the mean square displacement d of oscillator k (D k = 0.5d k 2 o k ; S k = 0.5d k 2 ).In practice, o k and D k are extracted from the ES trajectory (see Computational methods).

Fig. 2
Fig.2Schematic representation of the ES absorption signal in 2D electronic spectroscopy.The approximation of coherent superposition dynamics for the delay times t 1 and t 3 assumes that the wavepacket superpositions, created by the first pump pulse at time t 1 and the probe pulse at time t 3 , evolve coherently at two electronic potentials and nuclear decoherence effects are neglected.For each combination of delay time t 1 and waiting time t 2 the population (color coded and labeled), generated upon interaction with the second pump pulse at time t 2 , has propagated coherently on the electronic potentials, so that the third pulse probes a different segment of the higher ES manifold.The time-dependent fluctuation of the electronic gap is a function of (and hence contains information about) the active vibrational modes.In the semi-impulsive limit the signal intensity depends on the magnitude of the transition dipole moment at times t 1 , t 2 and t 3 .

Fig. 3
Fig.3Fluctuation of the L a -GS vertical gap along the scaled CAS(2,2)/ 6-31G* dynamics of pyrene on the bright L a state (green line).A constant scaling factor of 0.65 was derived by performing high level RAS(4,8|0,0|4,8)/ SS-CASPT2 computations at 50 geometries along the trajectory (8-10 fs, 26-28 fs, 44-46 fs, etc.).The energies of the L a and the L b states (obtained from the high level calculation) at these geometries are marked with red and blue dots, respectively.At several occasions along the trajectory (marked by red lines) a near-degeneracy and even inversion of the energies of the L a and L b states is observed indicating the vicinity of the crossing region.

Fig. 4
Fig. 4 Schematic representation of the photophysics of pyrene based on the 0K-trajectory dynamics augmented by high level electronic structure computations at relevant points.The dominant bond deformations are given, with the GS equilibrium geometry used as a reference.

Fig. 5
Fig. 5 Comparison of the experimental (black line) and theoretical (red line) linear absorption spectrum of pyrene recorded in aqueous solution at room temperature (left) and in a Ne-matrix at 6 K (right).The theoretical and experimental spectra are normalized with respect to the fundamental (0-0) transition.The spectral bandwidth was modeled through a lifetime broadening parameter of 125 cm À1 (corresponds to 85 fs lifetime of the L a state).Additionally, a dephasing rate constant of 200 cm À1 was used for the aqueous spectrum to model fluctuations due to solute-solvent interactions.The overtones of the three dominant modes which contribute to the vibrational progression are labeled.

Table 2
Vertical excitation energies DE (in eV and cm À1 ), transition dipole moments (in a.u.) out of the L b and wavefunctions for excited states responsible for the computed ESA peaks appearing in the spectrum in Fig.7