Dissipative dynamics at conical intersections: simulations with the hierarchy equations of motion method

The e ﬀ ect of a dissipative environment on the ultrafast nonadiabatic dynamics at conical intersections is analyzed for a two-state two-mode model chosen to represent the S 2 ( pp * ) – S 1 (n p * ) conical intersection in pyrazine (the system) which is bilinearly coupled to in ﬁ nitely many harmonic oscillators in thermal equilibrium (the bath). The system – bath coupling is modeled by the Drude spectral function. The equation of motion for the reduced density matrix of the system is solved numerically exactly with the hierarchy equation of motion method using graphics-processor-unit (GPU) technology. The simulations are valid for arbitrary strength of the system – bath coupling and arbitrary bath memory relaxation time. The present computational studies overcome the limitations of weak system – bath coupling and short memory relaxation time inherent in previous simulations based on multi-level Red ﬁ eld theory [A. K¨uhl and W. Domcke, J. Chem. Phys. 2002, 116 , 263]. Time evolutions of electronic state populations and time-dependent reduced probability densities of the coupling and tuning modes of the conical intersection have been obtained. It is found that even weak coupling to the bath e ﬀ ectively suppresses the irregular ﬂ uctuations of the electronic populations of the isolated two-mode conical intersection. While the population of the upper adiabatic electronic state (S 2 ) is very e ﬃ ciently quenched by the system – bath coupling, the population of the diabatic pp * electronic state exhibits long-lived oscillations driven by coherent motion of the tuning mode. Counterintuitively, the coupling to the bath can lead to an enhanced lifetime of the coherence of the tuning mode as a result of e ﬀ ective damping of the highly excited coupling mode, which reduces the strong mode – mode coupling inherent to the conical intersection. The present results extend previous studies of the dissipative dynamics at conical intersections to the nonperturbative regime of system – bath coupling. They pave the way for future ﬁ rst-principles simulations of femtosecond time-resolved four-wave-mixing spectra of chromophores in condensed phases which are nonperturbative in the system dynamics, the system – bath coupling as well as the ﬁ eld-matter coupling.


I. Introduction
The nonequilibrium dynamics of quantum systems which are interacting with a large (macroscopic) environment in thermal equilibrium is a problem with a long history in condensed-matter physics and chemistry, see, e.g., [1][2][3][4][5][6][7] and references therein. In the present work, we are addressing the question of how the complex nonadiabatic dynamics occurring at conical intersections of electronic potential-energy (PE) surfaces of polyatomic molecules [8][9][10][11] is affected by weak to moderate coupling to a condensed-phase environment. Conical intersections are hypersurfaces of exact degeneracy of adiabatic PE surfaces and are universally encountered in polyatomic molecules. [12][13][14][15] The simplest model of a conical intersection involves two electronic states and two vibrational modes which span the so-called branching space of the conical intersection. 16 In the case of a symmetry-allowed conical intersection, these vibrational modes can be classi-ed as coupling and tuning modes. 17 The singularity of the non-Born-Oppenheimer coupling elements at the seam of intersection and the pronounced local anharmonicity of the adiabatic PE surfaces cause exceptionally strong electronic inter-state couplings as well as vibrational mode-mode couplings. The interplay of these couplings gives rise to ultrafast radiationless relaxation processes following preparation of the upper electronic state by optical excitation. [8][9][10][11] In most applications explored in the physics literature, the quantum system of interest is simple and exhibits trivial dynamics when being isolated, such as a nuclear spin 1,2 or a harmonic oscillator. 3 Conical intersections, on the other hand, are multi-level systems with a complex internal dynamics. Even for minimal twostate two-mode models of conical intersections, the Hilbert space required for the description of the photoinduced dynamics can be quite large and the numerical diagonalization of the Hamiltonian can be challenging. Due to the intricate internal dynamics of conical intersections, the physics of energy and phase relaxation due to coupling to an environment can be much richer than the physics of a few-level system or a harmonic oscillator coupled to the same environment.
Early studies of the effects of an external bath on the photoinduced dynamics of two-mode and three-mode conical intersections were performed by Gerdts and Manthe, 18 Kühl and Domcke 19,20 and Baltzer and Stock, 21,22 assuming bilinear coupling of the active modes of the conical intersection to a harmonic bath with ohmic spectral density. The dissipative dynamics of these systems was numerically simulated with multi-level Redeld theory 23,24 which is valid within the limits of weak system-bath coupling and short memory relaxation time of the bath. Kosloff, Ratner and coworkers performed similar computational studies of conical intersections coupled to a bath, considering Lindblad-type relaxation operators or a bath of two-level systems. 25,26 Starting from multi-mode conical intersection Hamiltonians in the diabatic representation, Gindensperger, Burghardt and Cederbaum developed an effective-mode formalism which leads to a hierarchy of system-bath couplings. [27][28][29][30][31] In this approach, the short-time dynamics of the multi-mode conical intersection can be accurately reproduced with the inclusion of just a few effective bath modes. Burghardt and Hynes 32 and Malhado and Hynes 33,34 studied the effect of a dissipative solvent mode, described by a generalized Langevin equation, on the dynamics of a conical intersection related to retinal photoisomerization.
Herein, we employed the hierarchy equation of motion (HEOM) method developed by Tanimura and coworkers [35][36][37][38] to simulate the dynamics of a conical intersection in a dissipative environment. The HEOM method is a numerical technique for the solution of the Liouville equation of motion for the reduced density matrix of system-bath models. It is based on the construction and numerical solution of a hierarchy of master equations. 38 An alternative formulation of the HEOM method is based on Fokker-Planck equations. 39 The HEOM method is the most powerful method currently available for the numerically exact solution of master equations for a certain type of system-bath models. 38 In practical applications of the HEOM method, it is convenient to model the bath and the system-bath coupling by the Drude spectral function 7 or linear combination thereof. 40 The Drude bath is characterized by a system-bath coupling parameter l and a memory relaxation time g À1 . We will focus here on models with weak to moderate system-bath coupling, assuming that the strongly active modes of the conical intersection are included in the system Hamiltonian. Duan and Thorwart recently adopted an alternative strategy to simulate the dynamics of a two-state two-mode conical intersection coupled to a Drude bath. 41 By transforming the coupling and tuning modes of the conical intersection into the bath, they obtained an electronic two-level system which is coupled to a structured bath. The resulting master equation was solved with the HEOM method.

A. Hamiltonian
We partition the total Hamiltonian into the system Hamiltonian, the bath Hamiltonian, and their coupling We adopt a minimal model of a conical intersection which involves two excited electronic states |4 1 i, |4 2 i and two vibrational modes, the coupling mode and the tuning mode. Keeping only the leading (linear) terms in the electronic-vibrational couplings, the system Hamiltonian is written in the diabatic representation as 17 here, Q j , P j and U j are the dimensionless coordinate, momentum and frequency of the coupling mode (j ¼ c) and the tuning mode (j ¼ t), respectively. E k and k k (k ¼ 1,2) are the vertical excitation energies and the intra-state electron-vibrational coupling constants, respectively, of the kth electronic state. l s denotes the interstate coupling constant. We assume that the coupling mode Q c and the tuning mode Q t are bilinearly coupled to an external bath. Due to the different symmetry of the two system modes in the present model of a symmetry-allowed conical intersection, it makes sense to assume that each of the two modes is coupled to its own bath. The Hamiltonian of the environment is thus written as where p a,j , q a,j and u a,j are the dimensionless momentum, coordinate, and frequency of the ath oscillator of the bath of the coupling mode (j ¼ c) and the tuning mode (j ¼ t), respectively. We assume that the coupling and tuning vibrational modes are bilinearly coupled to the respective bath. The corresponding system-bath coupling Hamiltonian reads where c a,j are the system-bath coupling constants. Alternative choices of systembath coupling are discussed in ref. 20. The inuence of the bath on the system dynamics is determined by the spectral density The system-bath Hamiltonians specied as above can be used for deriving master equations in two different ways. [39][40][41][42][43] One can either retain the strongly coupled high-frequency vibrational modes in the system Hamiltonian or move the modes to the bath, resulting in a structured bath. The two descriptions are, in principle, equivalent and connected via a canonical transformation. [39][40][41][42][43][44] In the context of the HEOM method, the inclusion of all vibrational modes in the bath is the common choice. 40,41,45 In the present study we prefer to keep the strongly coupled high-frequency vibrational modes in the system Hamiltonian (cf. ref. 46). This way, we can treat fairly general system Hamiltonians (e.g., beyond the harmonic approximation for the system vibrational modes).
Introducing a basis set |n c i, |n t i for the coupling mode and the tuning mode, respectively, we can construct the matrix representation of the total Hamiltonian where In matrix form: The total (system + bath) dynamics is described by the Liouville-von Neumann equation The dynamics of the relevant system is described by the reduced density matrix r(t), dened as the trace over the bath degrees of freedom of the full density matrix r tot (t) We assume instantaneous electronic excitation of the system from the electronic ground state |4 0 i to the diabatic state |4 2 i by an ideally short laser pulse. This corresponds to the preparation of the density matrix where |0i ¼ |0 c i|0 t i denotes the vibrational ground state of the coupling and tuning modes in the electronic ground state. We model the bath by the Drude spectral density, 7 here, l j determines the coupling strength of the system modes (coupling mode and tuning mode) to the respective bath, and g j is the Drude decay constant, so that 1/g j is the bath correlation (memory) time. This choice of the spectral density facilitates the application of the HEOM method. 38

B. HEOM method
The HEOM method is a numerically exact method which allows the computation of the time evolution of the reduced density matrix r(t) beyond the Markovian approximation and for arbitrary system-bath coupling strength. [35][36][37][38]  Using the matrix representation of the Hamiltonian (6) and adopting the hightemperature approximation (bħg j < 1, where b ¼ (k B T) À1 is inverse temperature) which is fullled in our simulations, one obtains the following hierarchy of coupled equations of motion for the reduced density matrix v vt r lc;lt ðtÞ ¼ À À iL ðsÞ þ l c g c þ l t g t Á r lc;lt ðtÞ þ F c r lcþ1;lt ðtÞ þ F t r lc;ltþ1 ðtÞ þ l c g c Q c r lcÀ1;lt ðtÞ þ l t g t Q t r lc;ltÀ1 ðtÞ (15) where l c , l t are non-negative integers. In eqn (15), the density operator with all indexes equal to zero, r 00 (t), is the reduced density operator r(t). All other density operators r l c ,l t are auxiliary density operators which are introduced to describe nonperturbative and non-Markovian effects. The Liouvillian super operator corresponding to the system Hamiltonian H (S) is denoted by L (s) , and the relaxation operators F and Q are dened as where A r ¼ Ar + rA and A Â r ¼ Ar À rA for any operator A. Note that the HEOM method does not require the construction of the eigenstate representation for the system Hamiltonian, as multi-level Redeld theory does. 20 The system of eqn (15) represents an innite number of coupled master equations. It can be shown that the condition for the termination of the hierarchy 47 In this case, the last term in the innite hierarchy of eqn (15) can be replaced by v vt r lc;lt ðtÞ ¼ ÀiL ðsÞ r lc;lt ðtÞ Low-temperature correction terms have to be included in eqn (15) when the condition bħg j < 1 is not satised (see ref. 37 and 38). The low-temperature correction results in the replacements where n 1 ¼ 2p/bħ is the rst bosonic Matsubara frequency. 36 Higher-order Matsubara terms can be included to further improve the accuracy of low-temperature correction.

C. Physical observables
Once the reduced density matrix r(t) is given, the time evolution of system observables can straightforwardly be calculated. The most important observables for monitoring radiationless decay dynamics are the time-dependent populations of the electronic states. The population of the kth diabatic electronic state is dened as 10 The population of the adiabatic electronic states is given by The adiabatic electronic states are related to the diabatic states via where and a(Q c , Q t ) is the diabatic-to-adiabatic mixing angle. The adiabatic population P ad 2 is the observable which directly reects the radiationless electronic decay dynamics from the upper to the lower adiabatic energy surface. The diabatic electronic population P di 2 , on the other hand, is directly related to spectroscopic signals, e.g. to the total (frequency-integrated) spontaneous emission signal. 48 Note that the system modes Q c and Q t are traced out in the denition of P ad 2 and P di 2 and therefore can be considered as part of the bath. The coherence of vibrational motion of the system is monitored by the expectation values of the position and momentum operators of the system vibrational modes, which are dened as 10 For the nontotally symmetric coupling mode n c of the conical intersection, the expectation values of position and momentum operators vanish by symmetry. To visualize the coupled electronic-vibrational dynamics of the system, we consider the reduced probability densities of the coupling and tuning modes in the adiabatic electronic representation, dened as 10 Under certain conditions, these probability densities are related to the timeand frequency-resolved uorescence spectrum. 49

D. Parametrization of the Hamiltonian and computational details
The S 2 [ 1 B 2u (pp*)]-S 1 [ 1 B 3u (np*)] conical intersection in pyrazine is a well-known example for ultrafast electronic relaxation via a conical intersection. [50][51][52][53] Various models of this conical intersection, comprising from 2 to 24 vibrational modes, have been constructed. [54][55][56][57][58][59] These models have served as testbeds for the benchmarking of novel methods for the theoretical description of electronically nonadiabatic dynamics. [60][61][62][63][64][65][66][67][68][69][70][71][72] The out-of-plane normal mode n 10a of B 1g symmetry is the only normal mode which can couple the 1 B 2u and 1 B 3u electronic states in rst order. The mode n 6a of A 1g symmetry is the dominant tuning mode of this conical intersection. For the present purposes, a two-mode model of the S 2 -S 1 conical intersection including the normal modes n 10a and n 6a is suitable as a simple and generic model of a conical intersection.
The numerical values of the parameters of the model Hamiltonian H (S) of eqn (2) are taken from ref. 55. The frequencies of the coupling and tuning modes are ħu c ¼ 952 cm À1 (2p/u c ¼ 35 fs) and ħu t ¼ 597 cm À1 (2p/u t ¼ 56 fs). The intra-state electron-vibration coupling constants are k 1 ¼ À847 cm À1 and k 2 ¼ 1202 cm À1 . The inter-state vibronic coupling constant is l s ¼ 2110 cm À1 . The vertical excitation energies are E 1 ¼ 31 800 cm À1 and E 2 ¼ 39 000 cm À1 . Cuts through the potential-energy surfaces along the normal coordinates Q t and Q c are shown in Fig. 1. The bath temperature is set to 300 K.
The numerical calculation of the system dynamics is based on the representation of the Hamiltonian H (S) in a direct product basis of diabatic electronic states and harmonic-oscillator basis functions for the coupling and tuning modes. We have employed a basis set of 20 harmonic-oscillator functions for each mode, resulting in a dimension of 800 of the system Hamiltonian matrix H (S) . The numerical evaluation of r(t) via the HEOM method requires the propagation of a large number of coupled equations of motion for auxiliary matrices of large dimension, which is computationally demanding. We overcame the computational challenge by using a GPU (Graphic Processing Unit) implementation of the ArrayFire package which provides an easy-to-use API (application programming interface) and an array-based function set that facilitates GPU programming. 73 The numerical integration of the equations of motion for the auxiliary density matrices are performed on a NVIDIA Tesla M2050 GPU with 448 processors and 2 GB ECC-protected on-board memory. The HEOM eqn (15) has been solved via a fourth-order Runge-Kutta method with a time step of 0.2 fs. Depending on the value of system-bath coupling strength l j , j ¼ c, t, the truncation threshold N of the hierarchy was varied from 6 to 12, which was found to be sufficient for  obtaining converged results. With this efficient algorithm implemented on GPU hardware, the propagation of r(t) up to 1000 fs on a single GPU requires from several hours to 10 hours, depending on the truncation number N. The matrixmatrix multiplications in eqn (15) are carried out by the ArrayFire function matmul. The adiabatic populations P ad k are obtained by evaluating the expectation values of the projection operators on the lower and upper adiabatic electronic states in the diabatic representation, as described in ref. 74. The vibrational probability densities (25) are calculated by using a discrete variable representation (DVR) scheme. In this DVR scheme, a transformation between the harmonicoscillator basis and a grid in coordinate space is performed to represent operators in coordinate space (see ref. 75 for details).

III. Results
A. Electronic population probabilities and vibrational coherences Fig. 2(a) shows the time-dependent population P di 2 (t) of the diabatic pp* state for different system-bath coupling strengths on a time scale of 1 ps. The time evolution of P di 2 (t) for the isolated system (l c ¼ l t ¼ 0, black line) exhibits an initial decay on a timescale of z30 fs, followed by a few quasi-periodic recurrences with the period of the tuning mode. Beyond about 400 fs, the pp* population probability exhibits irregular uctuations. The long-time average of P di 2 (t) is z0.4. For the system-bath coupling strengths considered here, the coupling to the bath does not affect the fast initial decay of P di 2 (t) which is driven by the strong nonadiabatic coupling at the conical intersection. On the other hand, the irregular uctuations of P di 2 (t) for t $ 400 fs are wiped out by even very weak systembath coupling (l c ¼ l t ¼ 2.1 cm À1 , magenta line). With increasing, but still weak, system-bath coupling (l c ¼ l t ¼ 10.6 cm À1 , red line), P di 2 (t) exhibits, counterintuitively, a revival of the coherent oscillations driven by the tuning mode which extends now well beyond 1 ps. With further increasing system-bath coupling, the coherent oscillations of P di 2 (t) are damped out on a time scale of a few hundred femtoseconds, see Fig. 2(a), blue and green lines. The long-time limit of P di 2 (t) decreases with increasing system-bath coupling and is seen to converge to a small, but nonzero, value.
The time evolution of the population P ad 2 (t) of the adiabatic S 2 state is shown in Fig. 2(b). The initial decay of P ad 2 (t) of the isolated conical intersection (black line) occurs on a time scale of z20 fs and is somewhat faster than the initial decay of P di 2 (t). P ad 2 (t) approaches an asymptotic value of z0.2 at 1 ps. The system-bath coupling does not affect the short-time dynamics within the rst z100 fs. As observed for the diabatic population probability, the irregular uctuations of P ad 2 (t) are efficiently damped out by even very weak system-bath coupling (l c ¼ l t ¼ 2.1 cm À1 , magenta line). With increasing system-bath coupling, the population of the S 2 state quickly drops to zero, see Fig. 2(b). For the strongest system-bath coupling considered (l c ¼ l t ¼ 106.0 cm À1 ), all recurrences of P ad 2 (t) are quenched (green line). It should be noted that the long-lasting coherent oscillations of P di 2 (t) observed for l c ¼ l t $ 10.6 cm À1 are absent in P ad 2 (t). Overall, P ad 2 (t) reveals a bimodal relaxation dynamics. The ultrafast initial decay on a time scale of a few tens of femtoseconds originates from the nonadiabatic dynamics of the conical intersection itself. The coupling to the environment leads to a complete transfer of the population of the S 2 state to the S 1 state on sub-picosecond time scales for even weak system-bath coupling.
The expectation values of the position operator, hQ t i t , and the momentum operator, hP t i t , of the tuning mode are shown in Fig. 3(a) and (b), respectively. For the coupling mode, these expectation values vanish by symmetry. For the isolated conical intersection (black line), hQ t i t and hP t i t exhibit oscillations which are damped on a time scale of z300 fs. The mechanisms of this intrinsic vibrational dephasing of the conical intersection has been discussed in detail for the threemode model of the S 2 -S 1 conical intersection of pyrazine in ref. 55. A very weak system-bath coupling (l c ¼ l t ¼ 2.1 cm À1 magenta line) effectively damps residual uctuations of hQ t i t and hP t i t for t $ 500 fs. For stronger, but still weak, systembath coupling (l c ¼ l t ¼ 10.6 cm À1 , red line), a signicant enhancement of the amplitude of the tuning mode is observed. This "dedamping" of the tuning mode is the reason for the long-lived coherent oscillations of P di 2 (t) in Fig. 2(a). The mechanism of this effect will be discussed in Section IV. With further increasing system-bath coupling (blue and green lines in Fig. 3), the oscillations of hQ t i t and hP t i t are damped on time scales of a few hundred femtoseconds. Moreover, the long-time limit of hQ t i t is seen to shi towards an asymptotic value, which corresponds to the minimum of the adiabatic S 1 PE surface, see Fig. 1. This reects the ow of excess vibrational energy from the conical intersection into the environment.
Next we consider the effect of the bath relaxation time g À1 on the electronic populations and vibrational coherences. Fig. 4 shows the time-dependent population probabilities of the diabatic pp* state (a) and the adiabatic S 2 state (b), calculated for a longer bath relaxation time (g c À1 ¼ g t À1 ¼ 166 fs) and three different system-bath coupling strengths. Compared to Fig. 2 (bath relaxation time g c À1 ¼ g t À1 ¼ 50 fs), the coherent beatings of P di 2 (t) last longer in Fig. 4(a). The adiabatic populations P ad 2 (t) in Fig. 4(b) also decay more slowly than those in Fig. 2(b). The bath with longer relaxation time is less effective in draining coherence and energy from the dynamics of the system.
The corresponding expectation values of hQ t i t and hP t i t are shown in Fig. 5(a) and (b), respectively. The damping of the oscillations of hQ t i t and hP t i t due to system-bath coupling is less pronounced than in Fig. 3. The effect of dedamping of hQ t i t and hP t i t occurs for a higher value of the system-bath coupling strength, l c ¼ l t ¼ 53.0 cm À1 (blue line in Fig. 5). For strong system-bath coupling (l c ¼ l t ¼ 106.0 cm À1 (green line)), the oscillations of hQ t i t and hP t i t are completely quenched at z900 fs, while they were quenched already at z500 fs for the shorter bath relaxation time (Fig. 3).

B. Vibrational probability densities
The most detailed insight into the dynamics of the ultrafast internal-conversion process in the present model of a dissipative conical intersection is provided by the time-dependent probability densities P ad k (Q c ,t), P ad k (Q t ,t) dened in eqn (25). In what follows we call P ad k (Q j ,t), j ¼ c, t, probability densities of Q c and Q t , but it should be kept in mind that they are reduced probability densities, that is, the respective other vibrational mode is integrated out. Fig. 6 shows the probability densities P ad 2 (Q t ,t) (upper panels) and P ad 1 (Q t ,t) (lower panels) of the tuning mode for zero system-bath coupling (le column) and system-bath coupling of intermediate strength (l c ¼ l t ¼ 53.0 cm À1 , right column). Let us consider the isolated conical intersection rst (panels (a) and (c)). The localized wave packet prepared at t ¼ 0 on the S 2 surface propagates to the right (towards negative Q t ) following the gradient of the S 2 surface (see Fig. 1(a)). Within half of a vibrational period, the wave packet reaches the conical intersection and partially transfers to the S 1 surface, where it experiences a gradient in the opposite direction. This results in coherent wave-packet dynamics in the tuning mode at short times, see panels (a) and (c). For the isolated two-mode conical intersection, the S 2 / S 1 radiationless transition is incomplete. Beyond z200 fs, the structure of the probability density becomes more and more chaotic, especially on the S 1 surface. The coupling of the conical intersection to the environment results in a qualitatively different picture (panels (b) and (d)). The probability densities are smoothed and exhibit rather regular quasi-periodic behavior. The transfer of energy to the bath assists in the relaxation of the wave packet from the upper adiabatic surface to the lower adiabatic  surface. Recurrences of the probability density to the upper surface are strongly suppressed. For the dissipative conical intersection, the probability density of the tuning mode on the S 1 surface at t ¼ 1 ps is peaked at Q t z 1.8 (panel (d)), which is not the case for the isolated conical intersection (panel (c)). This illustrates that the vibrational excess energy available aer the S 2 / S 1 radiationless transition is effectively drained into the environment, allowing relaxation of the system to a stationary state. The probability densities P ad 2 (Q c ,t) (upper panels) and P ad 1 (Q c ,t) (lower panels) of the coupling mode are shown in Fig. 7. As in Fig. 6, the le column shows the bath-free dynamics, while the right column displays the dynamics of the conical intersection coupled to a bath with l c ¼ l t ¼ 53.0 cm À1 . For the isolated conical intersection (panels (a) and (c)) the probability density in Q c spreads very rapidly in a symmetric manner. The probability density on the upper surface exhibits three "ridges" at Q c ¼ 0 and Q c ¼ AE5 (the outer turning point in Q c ). The wave packet on the lower surface exhibits a symmetry-required node at Q c ¼ 0 and the width of the probability density shrinks with time. Beyond z200 fs, the probability densities on both surfaces are quite noisy. Panels (b) and (d) of Fig. 7 illustrate the effect of the system-bath coupling on the probability densities P ad 2 (Q c ,t) and P ad 1 (Q c ,t). The peak of P ad 2 (Q c ,t) at Q c ¼ 0 rapidly disappears, while the peaks near the turning points decrease in height and move inwards, reecting the loss of vibrational energy from the coupling mode. On the lower adiabatic surface, the width of the vibrational distribution narrows even more rapidly and the zero of P ad 1 (Q c ¼ 0, t ¼ 0) develops into a pronounced peak centered at Q c ¼ 0, see Fig. 7(d). As observed for the tuning mode, the probability densities of the coupling mode are much smoother than for the isolated two-mode conical intersection.

IV. Discussion
Considering the population dynamics of the diabatic pp* state, it is eye-catching how sensitive the electronic population dynamics is to the coupling of the system modes to an environment ( Fig. 2(a) and 4(a)). Even the weakest system-bath coupling considered here (l c ¼ l t ¼ 2.1 cm À1 ), completely quenches the irregular uctuations of P di 2 (t) beyond 400 fs ( Fig. 2(a), magenta line). With somewhat larger, but still weak, system-bath coupling (l c ¼ l t ¼ 10.6 cm À1 ), a qualitatively different behavior of P di 2 (t) is observed: the diabatic population probability exhibits completely regular coherent beatings with the frequency of the tuning mode ( Fig. 2(a), red line). With further increases in system-bath coupling, these coherent beatings are damped out on time scales of a few hundred femtoseconds. For a longer bath relaxation time (g c À1 ¼ g t À1 ¼ 166 fs), the coherent oscillations live longer. It is interesting to compare these numerically exact results with earlier results obtained with multi-level Redeld theory for the same two-state two-mode conical Fig. 6 Time dependent probability density of the tuning mode in the S 2 (upper row) and S 1 (lower row) adiabatic electronic states for bath relaxation times g c intersection coupled in the same manner to an ohmic bath. 19 In this earlier work, the secular approximation 4,24 to multi-level Redeld theory was employed. For a very weak system-bath coupling (h ¼ 0.003 in ref. 19, l c ¼ l t ¼ 2.1 cm À1 in Fig. 2(a)), the diabatic population probabilities are essentially identical. It can be concluded that for very weak system-bath coupling the spectral function of the bath does not matter and multi-level Redeld theory gives correct results even with the secular approximation. For somewhat stronger system-bath coupling (h ¼ 0.015 in ref. 19, l c ¼ l t ¼ 10.6 cm À1 in Fig. 2(a)), on the other hand, secular Redeld theory fails dramatically. The pronounced coherent beatings of P di 2 (t) (red line in Fig. 2(a)) are completely absent in the simulations with secular Redeld theory. 19 A later study performed with non-secular multi-level Redeld theory for a three-mode model of the S 2 -S 1 conical intersection in pyrazine 20 reproduces the coherent oscillations of P di 2 (t), albeit somewhat less pronounced and less regular than in Fig. 2(a), which may be a consequence of the additional tuning mode or the different spectral densities employed. A very similar result was obtained by Gelman et al. 26 for the two-mode model with a modied (Lindblad-type) Redeld theory. It can be concluded that non-secular multi-level Redeld theory gives a qualitatively correct description of the diabatic electronic population dynamics for weak coupling of a conical intersection to a bath (see also ref. 22). The mechanism of the dedamping of the tuning mode for intermediate system-bath coupling strength, which leads to the pronounced coherent beatings in P di 2 (t), was qualitatively explained in ref. 20. The coupling mode becomes particularly highly excited during the ultrafast internal conversion process at the conical intersection (see below) and is therefore highly susceptible to environmental damping. The efficient damping of the amplitude of the coupling mode reduces the intrinsic mode-mode coupling of the conical intersection, thus "liberating" the tuning mode to perform long-lived oscillations which in turn drive the diabatic electronic population dynamics. In addition, the center of the vibrational wave packet migrates towards the minimum of the S 1 PE surface at Q t z 1.4 and thus away from the conical intersection at Q t z À3.5 (see Fig. 2(a)) to a region where the PE surface is more harmonic than in the vicinity of the conical intersection. These phenomena are further illustrated by the dynamics of the expectation values of Q t and P t in Fig. 3 and 5. Fig. 3 shows that a very weak system-bath coupling (l c ¼ l t ¼ 2.1 cm À1 , magenta line) causes a damping of the oscillations of hQ t i t and hP t i t of the isolated conical intersection (black line), whereas a stronger system-bath coupling (l c ¼ l t ¼ 10.6 cm À1 ) leads to a pronounced increase of the oscillations of hQ t i t and hP t i t and a suppression of their damping (red line in Fig. 3), which is the dedamping effect discussed above. With further increase of the system-bath coupling, the coherent oscillation of the tuning mode is damped by the environment. For the strongest system-bath coupling considered here (l c ¼ l t ¼ 106.0 cm À1 , green line in Fig. 3), the wave packet becomes stationary at Q t z 2.0 at 500 fs. For a bath with longer memory time (g c À1 ¼ g t À1 ¼ 166 fs), the dedamping of hQ t i t and hP t i t is more pronounced at l c ¼ l t ¼ 53.0 cm À1 and the coherent oscillations of the tuning mode extend to longer times (blue line in Fig. 5).
Let us now consider the time evolution of the population probability P ad 2 (t) of the S 2 state. Again, a very weak system-bath coupling, l c ¼ l t ¼ 2.1 cm À1 , is remarkably efficient in suppressing the irregular uctuations of P ad 2 (t) of the isolated conical intersection (Fig. 2(b), magenta line). With increasing systembath coupling, the decay of P ad 2 (t) is signicantly enhanced. For l c ¼ l t $ 53.0 cm À1 , P ad 2 (t) decays within less than 200 fs. It is noteworthy that the long-lived coherent oscillations of P di 2 (t) are completely absent in P ad 2 (t). This observation underlines that the long-lived coherent wave-packet motion occurs exclusively on the lower adiabatic PE surface away from the conical intersection. With longer bath memory time, g c À1 ¼ g t À1 ¼ 166 fs, the quenching of P ad 2 (t) by the environment is less pronounced, see Fig. 2 The reduced probability densities of the coupling mode and the tuning mode in Fig. 6 and 7, respectively, provide additional insight. The probability densities in panels (a) and (c) clearly exhibit the irregular, chaotic character of the dynamics of the isolated two-state two-mode conical intersection. These irregularities are the time-domain signatures of the dense eigenvalue spectrum with incommensurable energy-level spacings which arise from the exceptionally strong nonadiabatic coupling at the conical intersection and the pronounced local anharmonicity of the adiabatic PE surfaces. These chaotic features of the dynamics at the conical intersection, which are reected in various level spacing distributions, 17 are highly sensitive to perturbations and are therefore eliminated by even very weak coupling to a dissipative environment. In contrast, the probability densities obtained for l c ¼ l t ¼ 53.0 cm À1 and shown in panels (b) and (d) of Fig. 6 and 7 are completely smooth and regular. The chaotic character of the dynamics was not observed in previous simulations for a three-mode model of the S 2 -S 1 conical intersection in pyrazine, 20,55 presumably because the second tuning mode of the conical intersection (the n 1 normal mode of pyrazine) already acts as a quencher of the chaotic dynamics.
The vibrational probability densities P ad k (Q c ,t) and P ad k (Q t ,t) of the dissipative conical intersection in Fig. 6(b) and (d) and 7(b) and (d) illustrate how the large excess energy in the tuning mode and the coupling mode, which is generated within z30 fs aer excitation, is drained out of the system dynamics, resulting in a nearly stationary wave packet at z1 ps. Comparing panels (a) and (b) of Fig. 7, it can be seen that the ridge of the density surviving on the S 2 surface at Q c ¼ 0 in panel (a) is turned into a valley by the coupling with the bath (panel (b)). The recurrences of the wave packet to the S 2 surface are strongly suppressed and occur near the minimum geometry of the S 1 PE surface in the Q c coordinate (see Fig. 1(b)), opposed to the isolated conical intersection, where the recurrences to the S 2 surface occur at Q c ¼ 0 and at the outer turning points of the motion in Q c , see Fig. 7(a).

V. Conclusions
We performed extensive simulations of the electronic and vibrational dynamics of a two-state two-mode conical intersection which is weakly to moderately coupled to a dissipative environment. The bath oscillators are bilinearly coupled to the coupling mode and the tuning mode of the conical intersection. The chaotic features of the intrinsic dynamics of the conical intersection, which are clearly visible in the electronic population dynamics and the reduced vibrational probability densities, are completely wiped out by even very weak system-bath coupling. While the ultrafast nonadiabatic dynamics within the rst z30 fs is unaffected by the coupling to the environment, the population probability of the upper adiabatic state is effectively quenched on sub-picosecond time scales by weak to moderate system-bath coupling. On the other hand, a counterintuitive dedamping of the dynamics of the tuning mode is observed for specic values of the system-bath coupling strength which depend on the bath memory time. This effect is a signature of the exceptionally pronounced nonseparability of the dynamics of the coupling and tuning modes of the conical intersection, to an extent that damping of the coupling mode liberates the dynamics of the tuning mode. Overall, the results conrm that the nonadiabatic wave-packet dynamics at a conical intersection is unusually sensitive to weak coupling to a thermal environment with the exception of the rst few tens of femtoseconds. Due to the complexity of the intrinsic dynamics of the conical intersection, the simulations exhibit phenomena which are quite different from the familiar scenarios of the optical Bloch equations or the Brownian oscillator. Our calculations show, in particular, that dissipative conical intersections not only can maintain vibrational wave-packets revealing the tuning modes, but also can create vibrational wavepackets featuring the coupling mode(s). Related effects have been observed in recent time and frequency resolved experiments on b-carotene and rhodopsin. 76,77 The present study extends earlier simulations of the dynamics of conical intersections coupled to an environment, which were based on multi-level Red-eld theory or similar models that are restricted to weak system-bath coupling and short bath relaxation time. [19][20][21][22]25,26 With the HEOM method, we were able to explore system-bath couplings in the nonperturbative regime. With this method, system-bath couplings of arbitrary strength can be investigated for complex systems such as conical intersections, as demonstrated recently also by Duan and Thorwart. 41 In future extensions of the present work, the coupling of the material system with femtosecond pulsed laser elds will be included in a nonperturbative manner. While a nonperturbative treatment of the laser-matter interaction may not be a necessity for the laser elds typically employed in spectroscopy, the inclusion of the laser-matter interaction in the equation of motion of the density matrix may be more straightforward and computationally efficient than the calculation of higher-order multi-time response functions and multiple nested integrals over the laser elds, which are required in the perturbative theory of nonlinear optics. 78 The calculation of femtosecond time-resolved four-wave-mixing signals from laser-driven equations of motion has been demonstrated, for example, for pump-probe and photon-echo spectroscopy, 46,79-81 two-dimensional electronic spectroscopy 82 or femtosecond resonant Raman spectroscopy. 83 The combination of nonperturbative HEOM simulations of dissipative nonadiabatic dynamics with the calculation of spectroscopic signals directly from eld-driven equations of motion will become a powerful tool for the analysis and/or prediction of state-of-the-art time and frequency resolved nonlinear spectra.