Stimulated Raman signals at conical intersections: Ab initio surface hopping simulation protocol with direct propagation of the nuclear wave function

Femtosecond Stimulated Raman Spectroscopy (FSRS) signals that monitor the excited state conical intersections dynamics of acrolein are simulated. An e ff ective time dependent Hamiltonian for two C—H vibrational marker bands is constructed on the fly using a local mode expansion combined with a semi-classical surface hopping simulation protocol. The signals are obtained by a direct forward and backward propagation of the vibrational wave function on a numerical grid. Earlier work is extended to fully incorporate the anharmonicities and intermode couplings. C 2015 AIP Publishing LLC.


I. INTRODUCTION
The radiation-less decay of electronically excited states through Conical Intersections (CoIn) is an important deactivation pathway commonly found in molecular systems. In the vicinity of a CoIn, the Born-Oppenheimer approximation breaks down and the nuclear and electronic motions are strongly coupled. These effects can be probed, for example, by transient optical absorption, 1 photoelectron spectroscopy, 2 transient infra-red (IR), 3,4 and Raman spectra. 5 Vibrational marker bands can serve as a sensitive probe of the excited state photochemistry. A transient vibrational spectrum in an electronic valence state is obtained by first pumping the system with a actinic UV/VIS pulse, launching the nuclear dynamics. A subsequent IR or Raman pulse can then probe the vibrational modes after a delay, capturing the fingerprints of the photochemical processes. In a Femtosecond Stimulated Raman Spectroscopy (FSRS) 6-8 experiment, the (off-resonant) Raman probe consists of a picosecond (narrowband) pulse and a femtosecond (broadband) pulse. Spectroscopic signatures of CoIns are provided by means of vibrational spectroscopy of spectator modes such as CH vibrations. The mode frequencies and the corresponding Raman intensities are sensitive to the electronic state and can be utilized to probe the CoIns.
The simulation of transient vibrational spectra requires to take into account the dynamics of the vibrational modes and the respective transition matrix elements. However, since a full quantum simulation of the nuclear and electronic dynamics is usually not feasible, some approximations have been applied. Different levels of theory have been used to describe the time evolution of the molecular degrees of freedom. The coupling between the electronic states and the nuclear subspace is commonly accounted for by using a Surface Hopping (SH) protocol. [9][10][11] By treating the nuclear degrees of freedom classically, it can handle larger molecules. The vibrational modes of interest (marker bands) can then be described, e.g., by a normal mode approach 9,12 in a harmonic approximation or a) Electronic mail: mkowalew@uci.edu b) Electronic mail: smukamel@uci.edu in the basis of local modes 13 that take anharmonicities into account. The time evolution of the vibrational wave function can be described by expanding it in instantaneous eigenstates. This sum over states (SOSs) formalism ignores variations of the eigenvalues during each snapshot. This can be corrected by a direct numerical propagation of the nuclear wave function on a numerical grid. This level of theory, which is adopted here properly describes the full anharmonic evolution of the marker bands.
A semi-classical approach has recently been demonstrated, involving a Non-Adiabatic On-the-fly Molecular Dynamics (NA-O-MD) scheme combined with a normal mode tracking approach. 9 This approach looks at the change of the frequency and the polarizability transition matrix elements of the spectator modes during the dynamics through the conical intersections. The spectator modes were taken to be harmonic and the time integrated harmonic eigenvalues were used to generate the spectrum. Upon UV excitation, kinetic energy is converted into vibrational energy and the spectator modes could be excited into higher quantum states. FSRS spectra can thus be expected to show hot bands. C-H bonds have large anharmonicities making a description beyond the harmonic approximation neccesary. We present a higher level simulation protocol, which takes into account the full anharmonic potential of the spectator modes generated by NA-O-MD. The Raman signals are calculated by a direct propagation of the vibrational wave function on a numerical grid. 14 This description includes the excited vibrational states and their anharmonicities, as well as the inter-mode couplings between the spectator modes. The method is applied to the excited state dynamics of the acrolein molecule. The main reaction mechanism is a fast internal conversion taking the system back into electronic ground state. 15 The conjugated π-system formed by the C= =C and the C= =O group is involved in the non-adiabatic couplings between the electronic states.
The paper is organized as follows. Section II introduces the spectroscopic scheme and the signals. In Section III, we describe the propagation techniques followed by a description of the handling of the marker bands in Section IV. The results are presented in Secs. VII and VIII.

II. THE FSRS SIGNAL
In a FSRS experiment, an actinic pump pulse E 1 creates a vibrational wave packet in an excited electronic state ( Fig. 1) whose time evolution is then tracked by a Raman probe. The probe consists of a spectrally narrow picosecond pulse E 3 and a short femtosecond broadband pulse, and the a frequency dispersed probe transmission is recorded. For electronic offresonant transitions, the polarizability is used to describe the two-photon stimulated Raman interactions. The effective radiation-matter interaction Hamiltonian is thus given by where α n is the polarizability in the excited electronic state, and V e is the electronic transition dipole moment. FSRS is a six wave mixing process derived from the loop diagrams 14 shown in Fig. 1. The signal obtained by the frequency dispersed detection of the broadband Raman pulse (E 2 ) depends on the delay time T between the actinic excitation and the Raman probe, where ω r = ω − ω 3 is the Raman shift, and ℑ denotes the imaginary part. A Gaussian broad band pulse is assumed: where ω 2 is the center frequency of the broad band pulse and σ ω and its width. The monochromatic narrow band pulse is The ∆-dispersed signal entering Eq. (2) can be read from the diagrams in Fig. 1 (5)).S where ∆ enters through the broadband pulse as E 2 (ω + ∆) and spans the field modes that contribute to the signal. The retarded Green's function G(t 1 , represents the forward propagation of the field-free molecular Hamiltonian H, and θ(t) is the Heaviside step function. By assuming an impulsive excitation with respect to the actinic pulse E 1 , Eqs. (4) and (5) simplify tõ The multi-point correlation functions ⟨· · · ⟩ in Eqs. (6) and (7) may be calculated exactly by a direct numerical propagation of the joint electronic and nuclear wave function. Instead, we use a mean field approximation. Non-adiabatic surface hopping trajectories generate an effective time dependent vibrational Hamiltonian for the spectator modes. If the eigenvalues do not change during the pulse, the correlation functions in Eq. 6 and 7 can be expanded in eigenstates of the total system. 14 Here, we do not make this approximation and employ a full direct numerical propagation for the nuclear wave function corresponding to the time dependent effective Hamiltonian. In this level of theory, the electronic degrees of freedom are included on a mean field level through the time dependent trajectory but the nuclear evolution in the subspace of the marker bands is exact. At the heart of the calculation is the evaluation of the multi-point correlation functions from Eqs. (6) and (7). By assuming impulsive excitation (τ 1 = τ 5 = 0), the calculation reduces to a two-point correlation function that depends on t and τ 3 . The initial condition is the vibrational ground state in the electronic ground state at t = 0. For Eq. (6), the evaluation of the correlation function goes as follows: The dipole operator V † e projects the vibrational wave function into an excited electronic state. It is then propagated from time zero to time t. The polarizability operator α n acts on the time evolved nuclear wave function and represents the Raman interaction with the fields E 2 and E 3 . The wave function is then propagated backward from t to τ 3 and α n is applied the second time before the final backward propagation from τ 3 to time 0 takes place. Finally, the dipole operator V e projects the nuclear wave function back to the electronic ground state and the inner product with the initial wave function is formed. Equation (7) can be described similarly but it involves two forward and one backward propagation.

III. DIRECT PROPAGATION OF THE MARKER BANDS
The C-H stretch modes serve as spectators to probe the excited state dynamics. This follows the ideas formulated in Ref. 9. We use a mixed quantum-classical mean-field approach for the respective spectator modes. A NA-O-MD protocol is used to describe the time evolution of the nuclei and the nonadiabatic dynamics of the electronic states. The SH method 16 as implemented in Newton-X 17 is used to obtain the NA-O-MD, treating the nuclei classically while maintaining a wave function for the electronic states. This allows for the simulation of electronic state changes along the conical intersection. To reintroduce a quantum mechanical description of the selected vibrational modes, the mean-field model is applied on top of the trajectory. The potential V (q) of a set of selected modes {q} is calculated along the SH trajectory yielding a time dependent Hamiltonian This treats the active vibrational modes of the system quantum mechanically and all other modes as a classical bath. The exponential part of the propagator in G(t 1 ,t 2 ) is calculated by the Split Operator (SPO) scheme, 18 where is the nuclear kinetic energy operator with the reduced masses µ i . The initial wave function Ψ(q,t = 0) is the vibrational ground state of the electronic ground state Potential Energy Surface (PES) at time zero. To obtain the multi-point correlation functions (Eq. 6 and 7), Ψ(q,t = 0) is propagated forward and backward by G with the time dependent potential V (q,t) which corresponds to the electronic state at time t (according to the SH trajectory). The SPO method used to propagate the active modes using a grid in position space is practical for up to 3 modes. For larger systems, an approximate method like Multi Configuration Time Dependent Hartree (MCTDH) 19 or Hagedorn wave packets 20 must be used. The calculation of the potential V (q,t) by quantum chemistry ab initio methods becomes tedious in higher dimensions. A many body expansion of the potential, 21 keeping only the two-mode couplings, can be used to reduce the computational effort.
The simulation protocol may be summarized as follows: first, a trajectory is calculated based on the SH algorithm. Second, for every time step of the trajectory, the potential of the spectator modes is calculated by using the geometries x(t) of the trajectory and displacing them along the spectator modes q. The nuclear potential V (q,t) (Eq. (8)) is then calculated for the respective electronic state s(t) and is used to propagate the initial state according to Eq. (9). At this stage, the wave function is in the vibrational subspace. The electronic wave function is accounted for in the trajectory calculation in step one. The effect of the electronic properties on the vibrational wave function enters through the potential V (q,t) and the polarizability α of the excited electronic states. The procedure has to be repeated for a number of trajectories and ensemble averaging must be performed at the signal level.

IV. LOCAL MODE DESCRIPTION
The two C-H stretch modes shown in Fig. 2 are chosen as marker bands and the potential V (q,t) is calculated in the corresponding subspace. The normal mode tracking approach 9 cannot be used in a straightforward manner here, since not only the second derivatives but also the actual several displacements along a normal mode need to be calculated. However, the Hessian matrix only delivers reliable eigenvectors in the vicinity of a potential minimum. Explicitly calculating the potential in the subspace of the marker bands by means of displacement vectors delivers a more exact and reliable description of the modes and allows for larger displacements. It also directly includes anharmonicities and intermode couplings. A local mode approach is used instead of the normal mode tracking technique. 9,22 The construction of modes from bond lengths provides a unique measure for bonds which are not compromised during a course of a trajectory. This offers a more reliable reference than normal modes. Local modes provide a good approximation for the C-H stretch vibrations. Due to their small mass compared to the total mass of the molecule, their motions can be very well described as localized to a single bond, or a few distinct bonds (e.g., symmetric and asymmetric C-H stretch in CH 2 groups).
We employ an approximate local mode representation which avoids non-diagonal kinetic energy operators to simplify the calculations. Other possibilities to construct local modes are the G-Matrix 13,23,24 method. A convenient way to obtain a Cartesian type kinetic energy operator, as given in Eq. (10), is to create displacement vectors along the desired bonds and to orthogonalize them in internal, mass weighted coordinates. The local modes are constructed as follows: the Cartesian displacement vector along the bond is defined by where x i and x j are the Cartesian coordinates of the ith and jth atoms. Assuming that only local modes for terminal H-atoms should be constructed, the displaced geometry is given by a matrix B holding the displacements of the respective terminal atoms for the n selected modes in its columns. In the next step, the rotational and translational motion has to be projected out. This can either be done by applying the Eckart-Sayvetz conditions 24 directly to the displacements vectors or by utilizing the transformation matrices of a normal mode analysis. Here, we chose to re-utilize the normal mode analysis, where C is the transformation matrix from Cartesian displacements to normal modes 24 obtained from a normal mode analysis. Setting the low-frequency components of F (i.e., the first six rows) to zero corresponds to projecting out rotations and translation in the center of mass frame. The columns of the resulting matrix F ′ are then orthogonalized (i.e., by a Gram-Schmidt process) to create a set of displacement vectors which are free of kinetic couplings. Transforming F ′ back from normal mode representation to Cartesian displacements yields the approximate local modes, The columns b ′ i of B ′ hold the Cartesian displacements used to calculate the marker band potentials. The reduced mass µ i is readily obtained from the norm of b ′ i : µ i = ∥b ′ i ∥. The potentials V (q,t) entering Eq. (8) are then obtained by calculating the PES for different scalar displacements q 1 and q 2 of the vectors b 1 and b 2 .

V. EXCITED STATE POLARIZABILITIES
Calculating the correlation functions in Eqs. (6) and (7), requires the polarizabilities of the electronically excited states. The polarizability operator α n was obtained by using the coordinate dependent polarizability α(q) in the static zero frequency limit. This is a good approximation for off-resonant Raman processes and can be readily obtained from a quantum chemistry calculation. The tensor components of α r s are obtained by numerical derivatives of the dipole moments with respect to an electric field, r and s represent the Cartesian components x, y, z of the dipole moments µ of a particular electronic state and the electric fields. By directly using the polarizability α(q) along the mode coordinate q, no strict selection rules are enforced. This is in contrast to the harmonic oscillator model, where the slope of α(q) is used, only allowing ∆v = ±1 transitions. Instead, the matrix elements α i j = φ i (q) α(q) φ j (q) allow for Rayleigh scattering (i = j), and for all kinds of transitions between different eigenstates φ i (q). This gives rise to an intense elastic Rayleigh scattering component in the signal which can interfere with the vibrational bands in the numerical calculations.
To suppress this component at ω r = 0, the polarizabilities are shifted by their instantaneous expectation value to minimize the diagonal matrix elements α ii , For the simulations, α ′ was used in Eqs. (6) and (7).

VI. INSTANTANEOUS RAMAN SPECTRA
To interpret the features of the FSRS spectra, it is helpful to know the transition frequencies and transitions matrix elements of the marker bands. However, the eigenvalues of the C-H stretch modes are not static but do change over time. Moreover, the eigenvalues of the coupled and anharmonic two-mode system cannot be easily separated. The transition frequencies and transitions matrix elements can be analyzed instead in terms of an instantaneous spectrum. Here, we use two different types of spectra which allow to analyze the dynamical properties of the marker bands.
The instantaneous vibrational frequencies of the {q 1 , q 2 } two-mode system can be visualized by weighing the finite width peaks by the respective transition matrix elements, where φ n (q; t) are the instantaneous eigenfunctions of the Hamiltonian H(t) (Eq. (8)) for a certain time t and α k k (q,t) the polarizability tensor elements. The energy difference between two eigenstates is denoted by ∆ mn = E m − E n and the width of the Gaussian is σ s = 5 cm −1 , representing a homogeneous linewidth. Here, the 10 lowest eigenvalues E k and Alternatively we can examine the instantaneous Raman spectrum. Here, the initial wave function Ψ(q,t = 0) is propagated according to Eq. (9), and in addition to Eq. (16), the intensity is weighted by the population of the respective eigenstates, Rotational averaging is accounted for by using the isotropic part of the polarizability tensor.

VII. ACROLEIN NON-ADIABATIC TRAJECTORIES
Acrolein (or propenal, Fig. 2) is formed from glycerol when cooking oil is overheated. We choose acrolein to demonstrate the computational method on a molecular system of moderate size. The conjugated π-system composed from the C==C double bond and the carbonyl group (C==O) is involved in the lower lying electronic excitations 26 (see Fig. 3). Three excited states are accessible for a UV excitation, which are all connected with each other by CoIns and the S 3 state can decay

Transition
Calculated Experimental S 0 → S 1 (n → π * ) 3.6279 S 0 → S 2 (n → π * ) 6.8740 S 0 → S 3 (π → π * ) 7.5165 6.41(193 nm) radiation-less back to the electronic ground state (Fig. 5). The branching space of the CoIns between the states S 0 to S 3 is thus composed from the Molecular Orbitals (MOs) of the πsystem. Upon UV excitation into the S 3 state, the π-bonds are weakened by a π → π * excitation and a structural change of the carbon back bone is induced. The H-atoms serve as spectators and can thus be used as marker bands in the FSRS signal for the non-adiabatic excited state dynamics. The quantum chemistry calculations are carried out with MOLPRO 27 at the Complete Active Space Self Consistent Field (CASSCF) level of theory (sa4-CAS(6/5)/6-31+G * ) with 6 electrons in 5 active orbitals. All CoIns are optimized with the branching space algorithm implemented in MOLPRO. The occupied MOs in the electronic ground state are two π-orbitals for the CC and the CO bond as well as the oxygen lone pair. The Lowest Unoccupied MO (LUMO) corresponds to the first two π * -orbitals. It turns out that at the given level of theory, 4 states have to be taken in account. The S 2 and S 3 states are separated by ≈0.65 eV at the Franck-Condon point (Table I). Moreover, the S 2 and S 3 states are coupled by a conical intersection (see Fig. 4) which is in reach within the first 25 fs of the time evolution. The excitation energies are calculated at the CAS level of theory are given in Table I. A vertical, resonant UV excitation is most favorable for the S 0 → S 3 transition with a transition dipole moment of 0.18 a.u. as opposed to the S 0 → S 2 transition, which is much weaker (µ S2, S0 = 0.044 a.u.).
The nuclear dynamics triggered by the UV excitation is calculated by Tully's SH algorithm 16 (NEWTON-X 17 with a home made interface to MOLPRO). The dynamics of the nuclei is described by the classical equations of motion, while the electronic degrees of freedom are treated quantum mechanically (sa4-CAS(5/6)/6-31+G * ). The forces acting on the nuclei are derived from the gradient of a particular electronic state the trajectory is in at time T. The population P I of an electronic state S I is obtained by the average over all N trajectories, P is the electronic state of trajectory i. The initial state populated by the UV pulse is the electronic S 3 state. We assumed sudden excitation, thereby neglecting Franck-Condon factors. All four states are coupled to each other by CoIns and consecutive state jumps back to the ground state are observed in our trajectory simulations as well as under experimental conditions. 15 The time evolution for the first 500 fs of P I is shown in Fig. 5.
We focused on the S 2 /S 1 , and the S 1 /S 0 CoIn for the spectroscopic observations. The S 2 /S 1 CoIn strongly involves the C= =O bond (Fig. 6(a)) and thus the bonds C-H 3 and C-H 4 bonds are good candidates for observing spectroscopic signatures of this CoIn. The S 1 /S 0 CoIn affects the C= =C double bond (Fig. 6(b)) making the C-H 3 bond a good probe.

VIII. RESULTS
The high frequency normal modes, calculated for the S 0 state, are shown in Table II bond adjacent to the carbonyl bond (C-H 4 and C-H 3 ). We chose these as marker bands for the non-adiabatic dynamics. Both modes are highly localized to their respective C-H bonds in the normal mode picture and are thus suitable for the local mode description. The coordinates for the C-H 4 and C-H 3 will labeled as q 1 and q 2 , respectively. The selected modes are shown in Fig. 2. In the following, we show the direct propagation in the {q 1 , q 2 } mode system for a single trajectory as a basic proof of principle. The respective time evolution of the electronic energies and populations is shown in Fig. 7. At 147 fs, the molecule reaches the S 2 -S 1 CoIn and changes its state from S 2 to S 1 and back at 224 fs. Eventually, it reaches the CoIn at t = 355 fs and jumps back to the S 1 state.
The instantaneous spectra of the fundamental frequencies and hot bands corresponding to ∆v = 1 transitions resulting from Eq. tional excited states, can be observed around 3320 cm −1 and 3200 cm −1 . The anharmonicity of the C-H stretch vibrations is ≈100 cm −1 justifying a model beyond the harmonic approximation. Note that transitions between higher vibrational states have increasing matrix elements. The eigenvalues are significantly disturbed in the vicinity of the CoIn at 147 fs, 224 fs, and 355 fs. For comparison, the instantaneous Raman spectra are shown in Fig. 9. A comparison of the instantaneous spectrum in Fig. 9 with the instantaneous eigenvalues in Fig. 8 shows that the vibrational excited state populations are non-vanishing in the course of the trajectory.
The frequency dispersed FSRS spectrum (Eq. (2)) is shown in Fig. 10   A comparison of the FSRS and the spontaneous Raman spectrum (Fig. 11) demonstrates that a FSRS spectrum cannot simply be viewed as a simple snapshot Raman spectra with temporal broadening. 9,28 The numerical evaluation of Eqs. (6) and (7) slightly affects the line shape of the simulated spectra: the integration of the correlation function over finite time intervals behaves like a rectangular window function. Absorptive resonances thus appear as sinc-functions with minor side lobes, which can be seen, for example, in Fig. 10 at a delay of T = 50 fs. However, the basic features of the instantaneous spectrum qualitatively match the FSRS spectrum: the signal starts with a broad band centered at 3450 cm −1 , which represents the initial frequencies of the C-H stretch vibrations. Both modes, q 1 and q 2 are nearly degenerate after UV excitation in this setting and can thus appear as a single peak. At around 147 fs, there is a short lasting change in the frequency of the modes to lower frequencies (≈3100 cm −1 ) which is also present in the FSRS spectrum, however, with a more complicated shape. The passage through the CoIn back to the S 2 state (224 fs) is characterized by a drop in signal intensity. Between 300 and 400 fs, the frequencies of the modes split FIG. 12. FSRS spectrum of the coupled mode system {q 1 , q 2 }. ω 2 − ω 3 = −3500 cm −1 , σ ω = 3000 cm −1 . The time of passage through a CoIn is marked by the bold lines. and one transition is shifted towards 4000 cm −1 . This is less pronounced in the FSRS spectrum but visible. Together with the state change to the S 1 at 355 fs, the intensity of the band around 3500 cm −1 disappears.

IX. CONCLUSIONS
We have demonstrated the calculation of FSRS spectra by the combination of a semi-classical SH approach with the direct propagation of vibrational marker bands. The explicit treatment of the marker band nuclear potential takes into account anharmonicities, allowing larger displacement than in a harmonic description. This level of theory also includes the inter-mode coupling between the marker bands. Proper direct propagation of the nuclear wave function allows for a description of interference effects in situations where the transitions overlap in frequency. In principle, this technique would also allow for the inclusion of Fermi resonances. The combination of a SH hopping approach as a first step includes all nuclear degrees of freedom and the relevant electronic state reproducing the correct system dynamics. The potential obtained from the SH trajectories allows for a more explicit description of a few marker bands. The CoIn signatures appear in frequency shifts as well as intensity changes. Spectra for different pulse parameters demonstrate the lineshapes for different resolutions in the time domain.