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

Abstract

The effect 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₂(ππ*)–S₁(nπ*) conical intersection in pyrazine (the system) which is bilinearly coupled to infinitely 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 Redfield theory [A. Kühl 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 effectively suppresses the irregular fluctuations of the electronic populations of the isolated two-mode conical intersection. While the population of the upper adiabatic electronic state (S₂) is very efficiently quenched by the system–bath coupling, the population of the diabatic ππ* 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 effective 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 first-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 field-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–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–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–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.¹⁶ In the case of a symmetry-allowed conical intersection, these vibrational modes can be classified as coupling and tuning modes.¹⁷ 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–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.³ Conical intersections, on the other hand, are multi-level systems with a complex internal dynamics. Even for minimal two-state 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,¹⁸ 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 Redfield 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–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³² 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–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.³⁸ An alternative formulation of the HEOM method is based on Fokker–Planck equations.³⁹ 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.³⁸ In practical applications of the HEOM method, it is convenient to model the bath and the system–bath coupling by the Drude spectral function⁷ or linear combination thereof.⁴⁰ The Drude bath is characterized by a system–bath coupling parameter λ and a memory relaxation time γ⁻¹. 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.⁴¹ 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.

II. Theoretical framework

A. Hamiltonian

We partition the total Hamiltonian into the system Hamiltonian, the bath Hamiltonian, and their coupling

H = H^(S) + H^(B) + H^(SB) (1)

We adopt a minimal model of a conical intersection which involves two excited electronic states |φ₁〉, |φ₂〉 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¹⁷

here, Q_j, P_j and Ω_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 = 1,2) are the vertical excitation energies and the intra-state electron-vibrational coupling constants, respectively, of the kth electronic state. λ_s denotes the inter-state 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_α,j, q_α,j and ω_α,j are the dimensionless momentum, coordinate, and frequency of the αth 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_α,j are the system–bath coupling constants. Alternative choices of system–bath coupling are discussed in ref. 20. The influence of the bath on the system dynamics is determined by the spectral density

The system–bath Hamiltonians specified as above can be used for deriving master equations in two different ways.^39–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–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〉, |n_t〉 for the coupling mode and the tuning mode, respectively, we can construct the matrix representation of the total Hamiltonian

where
with

In matrix form:

where

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 ρ(t), defined as the trace over the bath degrees of freedom of the full density matrix ρ_tot(t)

ρ(t) = tr_B{ρ_tot(t)} (12)

We assume instantaneous electronic excitation of the system from the electronic ground state |φ₀〉 to the diabatic state |φ₂〉 by an ideally short laser pulse. This corresponds to the preparation of the density matrix

ρ(0) = |φ₂〉|0〉〈0|〈φ₂| (13)

where |0〉 = |0_c〉|0_t〉 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,⁷

J_j(ω) = 2λ_jγ_jω/(ω² + γ_j²), j = c, t (14)

here, λ_j determines the coupling strength of the system modes (coupling mode and tuning mode) to the respective bath, and γ_j is the Drude decay constant, so that 1/γ_j is the bath correlation (memory) time. This choice of the spectral density facilitates the application of the HEOM method.³⁸

B. HEOM method

The HEOM method is a numerically exact method which allows the computation of the time evolution of the reduced density matrix ρ(t) beyond the Markovian approximation and for arbitrary system–bath coupling strength.^35–38

Using the matrix representation of the Hamiltonian (6) and adopting the high-temperature approximation (βħγ_j < 1, where β = (k_BT)⁻¹ is inverse temperature) which is fulfilled in our simulations, one obtains the following hierarchy of coupled equations of motion for the reduced density matrix

where l_c, l_t are non-negative integers. In eqn (15), the density operator with all indexes equal to zero, ρ₀₀(t), is the reduced density operator ρ(t). All other density operators ρ_lc,lt 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 [script L]^(s), and the relaxation operators Φ and Θ are defined as

Φ_j = iV^×_j (16)

where A°ρ = Aρ + ρA and A^×ρ = Aρ − ρA 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 Redfield theory does.²⁰

The system of eqn (15) represents an infinite number of coupled master equations. It can be shown that the condition for the termination of the hierarchy is N = l_c + l_t ≫ ω^(S)/min(γ₁, γ₂), where ω^(S) is a characteristic frequency of [script L]^(s).⁴⁷ In this case, the last term in the infinite hierarchy of eqn (15) can be replaced by

Low-temperature correction terms have to be included in eqn (15) when the condition βħγ_j < 1 is not satisfied (see ref. 37 and 38). The low-temperature correction results in the replacements

where ν₁ = 2π/βħ is the first bosonic Matsubara frequency.³⁶ Higher-order Matsubara terms can be included to further improve the accuracy of low-temperature correction.

C. Physical observables

Once the reduced density matrix ρ(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 defined as¹⁰

P^di_k(t) = tr_S{|φ_k〉〈φ_k|ρ(t)} (20)

The population of the adiabatic electronic states is given by

P^ad_k(t) = tr_S{|[small variant phi, Greek, tilde]_k〉〈[small variant phi, Greek, tilde]_k|ρ(t)} (21)

The adiabatic electronic states are related to the diabatic states via

where
and α(Q_c, Q_t) is the diabatic-to-adiabatic mixing angle. The adiabatic population P^ad₂ is the observable which directly reflects the radiationless electronic decay dynamics from the upper to the lower adiabatic energy surface. The diabatic electronic population P^di₂, on the other hand, is directly related to spectroscopic signals, e.g. to the total (frequency-integrated) spontaneous emission signal.⁴⁸ Note that the system modes Q_c and Q_t are traced out in the definition of P^ad₂ and P^di₂ 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 defined as¹⁰

〈Q_j〉_t = tr_S{Q_jρ(t)} (23)

〈P_j〉_t = tr_S{P_jρ(t)} (24)

For the nontotally symmetric coupling mode ν_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, defined as¹⁰

Under certain conditions, these probability densities are related to the time- and frequency-resolved fluorescence spectrum.⁴⁹

D. Parametrization of the Hamiltonian and computational details

The S₂[¹B_2u(ππ*)]–S₁[¹B_3u(nπ*)] conical intersection in pyrazine is a well-known example for ultrafast electronic relaxation via a conical intersection.^50–53 Various models of this conical intersection, comprising from 2 to 24 vibrational modes, have been constructed.^54–59 These models have served as testbeds for the benchmarking of novel methods for the theoretical description of electronically nonadiabatic dynamics.^60–72 The out-of-plane normal mode ν_10a of B_1g symmetry is the only normal mode which can couple the ¹B_2u and ¹B_3u electronic states in first order. The mode ν_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₂–S₁ conical intersection including the normal modes ν_10a and ν_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 ħω_c = 952 cm⁻¹ (2π/ω_c = 35 fs) and ħω_t = 597 cm⁻¹ (2π/ω_t = 56 fs). The intra-state electron-vibration coupling constants are κ₁ = −847 cm⁻¹ and κ₂ = 1202 cm⁻¹. The inter-state vibronic coupling constant is λ_s = 2110 cm⁻¹. The vertical excitation energies are E₁ = 31 800 cm⁻¹ and E₂ = 39 000 cm⁻¹. 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.

Fig. 1 Cuts through the adiabatic potential-energy surfaces of the electronic ground state (blue), S₁ state (green) and S₂ state (red) of the pyrazine model along the normal coordinates Q_t (a) and Q_c (b).

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 ρ(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.⁷³ 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 λ_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 ρ(t) up to 1000 fs on a single GPU requires from several hours to 10 hours, depending on the truncation number N. The matrix–matrix 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 harmonic-oscillator 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₂(t) of the diabatic ππ* state for different system–bath coupling strengths on a time scale of 1 ps. The time evolution of P^di₂(t) for the isolated system (λ_c = λ_t = 0, black line) exhibits an initial decay on a timescale of ≈30 fs, followed by a few quasi-periodic recurrences with the period of the tuning mode. Beyond about 400 fs, the ππ* population probability exhibits irregular fluctuations. The long-time average of P^di₂(t) is ≈0.4. For the system–bath coupling strengths considered here, the coupling to the bath does not affect the fast initial decay of P^di₂(t) which is driven by the strong nonadiabatic coupling at the conical intersection. On the other hand, the irregular fluctuations of P^di₂(t) for t ≥ 400 fs are wiped out by even very weak system–bath coupling (λ_c = λ_t = 2.1 cm⁻¹, magenta line). With increasing, but still weak, system–bath coupling (λ_c = λ_t = 10.6 cm⁻¹, red line), P^di₂(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₂(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₂(t) decreases with increasing system–bath coupling and is seen to converge to a small, but nonzero, value.

Fig. 2 Diabatic (a) and adiabatic (b) population probability of the upper electronic state of the pyrazine model for bath relaxation times γ_c⁻¹ = γ_t⁻¹ = 50 fs and different system–bath coupling strengths: λ_c = λ_t = 0 (black), 2.1 (magenta), 10.6 (red), 53.0 (blue), and 106.0 cm⁻¹ (green).

The time evolution of the population P^ad₂(t) of the adiabatic S₂ state is shown in Fig. 2(b). The initial decay of P^ad₂(t) of the isolated conical intersection (black line) occurs on a time scale of ≈20 fs and is somewhat faster than the initial decay of P^di₂(t). P^ad₂(t) approaches an asymptotic value of ≈0.2 at 1 ps. The system–bath coupling does not affect the short-time dynamics within the first ≈100 fs. As observed for the diabatic population probability, the irregular fluctuations of P^ad₂(t) are efficiently damped out by even very weak system–bath coupling (λ_c = λ_t = 2.1 cm⁻¹, magenta line). With increasing system–bath coupling, the population of the S₂ state quickly drops to zero, see Fig. 2(b). For the strongest system–bath coupling considered (λ_c = λ_t = 106.0 cm⁻¹), all recurrences of P^ad₂(t) are quenched (green line). It should be noted that the long-lasting coherent oscillations of P^di₂(t) observed for λ_c = λ_t ≥ 10.6 cm⁻¹ are absent in P^ad₂(t). Overall, P^ad₂(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₂ state to the S₁ state on sub-picosecond time scales for even weak system–bath coupling.

The expectation values of the position operator, 〈Q_t〉_t, and the momentum operator, 〈P_t〉_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), 〈Q_t〉_t and 〈P_t〉_t exhibit oscillations which are damped on a time scale of ≈300 fs. The mechanisms of this intrinsic vibrational dephasing of the conical intersection has been discussed in detail for the three-mode model of the S₂–S₁ conical intersection of pyrazine in ref. 55. A very weak system–bath coupling (λ_c = λ_t = 2.1 cm⁻¹ magenta line) effectively damps residual fluctuations of 〈Q_t〉_t and 〈P_t〉_t for t ≥ 500 fs. For stronger, but still weak, system–bath coupling (λ_c = λ_t = 10.6 cm⁻¹, red line), a significant 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₂(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 〈Q_t〉_t and 〈P_t〉_t are damped on time scales of a few hundred femtoseconds. Moreover, the long-time limit of 〈Q_t〉_t is seen to shift towards an asymptotic value, which corresponds to the minimum of the adiabatic S₁ PE surface, see Fig. 1. This reflects the flow of excess vibrational energy from the conical intersection into the environment.

Fig. 3 Expectation values of position (a) and momentum (b) operator of the tuning mode for bath relaxation times γ_c⁻¹ = γ_t⁻¹ = 50 fs and different system–bath coupling strengths: λ_c = λ_t = 0 (black), 2.1 (magenta), 10.6 (red), 53.0 (blue), and 106.0 cm⁻¹ (green).

Next we consider the effect of the bath relaxation time γ⁻¹ on the electronic populations and vibrational coherences. Fig. 4 shows the time-dependent population probabilities of the diabatic ππ* state (a) and the adiabatic S₂ state (b), calculated for a longer bath relaxation time (γ_c⁻¹ = γ_t⁻¹ = 166 fs) and three different system–bath coupling strengths. Compared to Fig. 2 (bath relaxation time γ_c⁻¹ = γ_t⁻¹ = 50 fs), the coherent beatings of P^di₂(t) last longer in Fig. 4(a). The adiabatic populations P^ad₂(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.

Fig. 4 Diabatic (a) and adiabatic (b) population probability of the upper electronic state of the pyrazine model for bath relaxation times γ_c⁻¹ = γ_t⁻¹ = 166 fs and different system–bath coupling strengths: λ_c = λ_t = 10.6 (red), 53.0 (blue), and 106.0 cm⁻¹ (green).

The corresponding expectation values of 〈Q_t〉_t and 〈P_t〉_t are shown in Fig. 5(a) and (b), respectively. The damping of the oscillations of 〈Q_t〉_t and 〈P_t〉_t due to system–bath coupling is less pronounced than in Fig. 3. The effect of dedamping of 〈Q_t〉_t and 〈P_t〉_t occurs for a higher value of the system–bath coupling strength, λ_c = λ_t = 53.0 cm⁻¹ (blue line in Fig. 5). For strong system–bath coupling (λ_c = λ_t = 106.0 cm⁻¹ (green line)), the oscillations of 〈Q_t〉_t and 〈P_t〉_t are completely quenched at ≈900 fs, while they were quenched already at ≈500 fs for the shorter bath relaxation time (Fig. 3).

Fig. 5 Expectation values of the position (a) and momentum (b) operator of the tuning mode for bath relaxation times γ_c⁻¹ = γ_t⁻¹ = 166 fs and different system–bath coupling strengths: λ_c = λ_t = 10.6 (red), 53.0 (blue), and 106.0 cm⁻¹ (green).

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) defined 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₂(Q_t,t) (upper panels) and P^ad₁(Q_t,t) (lower panels) of the tuning mode for zero system–bath coupling (left column) and system–bath coupling of intermediate strength (λ_c = λ_t = 53.0 cm⁻¹, right column). Let us consider the isolated conical intersection first (panels (a) and (c)). The localized wave packet prepared at t = 0 on the S₂ surface propagates to the right (towards negative Q_t) following the gradient of the S₂ surface (see Fig. 1(a)). Within half of a vibrational period, the wave packet reaches the conical intersection and partially transfers to the S₁ 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₂ → S₁ radiationless transition is incomplete. Beyond ≈200 fs, the structure of the probability density becomes more and more chaotic, especially on the S₁ 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₁ surface at t = 1 ps is peaked at Q_t ≈ 1.8 (panel (d)), which is not the case for the isolated conical intersection (panel (c)). This illustrates that the vibrational excess energy available after the S₂ → S₁ radiationless transition is effectively drained into the environment, allowing relaxation of the system to a stationary state.

Fig. 6 Time dependent probability density of the tuning mode in the S₂ (upper row) and S₁ (lower row) adiabatic electronic states for bath relaxation times γ_c⁻¹ = γ_t⁻¹ = 50 fs. Left column: λ_c = λ_t = 0; right column: λ_c = λ_t = 53.0 cm⁻¹.

The probability densities P^ad₂(Q_c,t) (upper panels) and P^ad₁(Q_c,t) (lower panels) of the coupling mode are shown in Fig. 7. As in Fig. 6, the left column shows the bath-free dynamics, while the right column displays the dynamics of the conical intersection coupled to a bath with λ_c = λ_t = 53.0 cm⁻¹. 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 = ±5 (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 ≈200 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₂(Q_c,t) and P^ad₁(Q_c,t). The peak of P^ad₂(Q_c,t) at Q_c = 0 rapidly disappears, while the peaks near the turning points decrease in height and move inwards, reflecting 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₁(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.

Fig. 7 Time dependent probability density of the coupling mode in the S₂ (upper row) and S₁ (lower row) adiabatic electronic states for bath relaxation times γ_c⁻¹ = γ_t⁻¹ = 50 fs. Left column: λ_c = λ_t = 0; right column: λ_c = λ_t = 53.0 cm⁻¹.

IV. Discussion

Considering the population dynamics of the diabatic ππ* 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 (λ_c = λ_t = 2.1 cm⁻¹), completely quenches the irregular fluctuations of P^di₂(t) beyond 400 fs (Fig. 2(a), magenta line). With somewhat larger, but still weak, system–bath coupling (λ_c = λ_t = 10.6 cm⁻¹), a qualitatively different behavior of P^di₂(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 (γ_c⁻¹ = γ_t⁻¹ = 166 fs), the coherent oscillations live longer.

It is interesting to compare these numerically exact results with earlier results obtained with multi-level Redfield theory for the same two-state two-mode conical intersection coupled in the same manner to an ohmic bath.¹⁹ In this earlier work, the secular approximation^4,24 to multi-level Redfield theory was employed. For a very weak system–bath coupling (η = 0.003 in ref. 19, λ_c = λ_t = 2.1 cm⁻¹ 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 Redfield theory gives correct results even with the secular approximation. For somewhat stronger system–bath coupling (η = 0.015 in ref. 19, λ_c = λ_t = 10.6 cm⁻¹ in Fig. 2(a)), on the other hand, secular Redfield theory fails dramatically. The pronounced coherent beatings of P^di₂(t) (red line in Fig. 2(a)) are completely absent in the simulations with secular Redfield theory.¹⁹ A later study performed with non-secular multi-level Redfield theory for a three-mode model of the S₂–S₁ conical intersection in pyrazine²⁰ reproduces the coherent oscillations of P^di₂(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.²⁶ for the two-mode model with a modified (Lindblad-type) Redfield theory. It can be concluded that non-secular multi-level Redfield 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₂(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₁ PE surface at Q_t ≈ 1.4 and thus away from the conical intersection at Q_t ≈ −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 (λ_c = λ_t = 2.1 cm⁻¹, magenta line) causes a damping of the oscillations of 〈Q_t〉_t and 〈P_t〉_t of the isolated conical intersection (black line), whereas a stronger system–bath coupling (λ_c = λ_t = 10.6 cm⁻¹) leads to a pronounced increase of the oscillations of 〈Q_t〉_t and 〈P_t〉_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 (λ_c = λ_t = 106.0 cm⁻¹, green line in Fig. 3), the wave packet becomes stationary at Q_t ≈ 2.0 at 500 fs. For a bath with longer memory time (γ_c⁻¹ = γ_t⁻¹ = 166 fs), the dedamping of 〈Q_t〉_t and 〈P_t〉_t is more pronounced at λ_c = λ_t = 53.0 cm⁻¹ 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₂(t) of the S₂ state. Again, a very weak system–bath coupling, λ_c = λ_t = 2.1 cm⁻¹, is remarkably efficient in suppressing the irregular fluctuations of P^ad₂(t) of the isolated conical intersection (Fig. 2(b), magenta line). With increasing system–bath coupling, the decay of P^ad₂(t) is significantly enhanced. For λ_c = λ_t ≥ 53.0 cm⁻¹, P^ad₂(t) decays within less than 200 fs. It is noteworthy that the long-lived coherent oscillations of P^di₂(t) are completely absent in P^ad₂(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, γ_c⁻¹ = γ_t⁻¹ = 166 fs, the quenching of P^ad₂(t) by the environment is less pronounced, see Fig. 2(b).

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 reflected in various level spacing distributions,¹⁷ 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 λ_c = λ_t = 53.0 cm⁻¹ 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₂–S₁ conical intersection in pyrazine,^20,55 presumably because the second tuning mode of the conical intersection (the ν₁ 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 ≈30 fs after excitation, is drained out of the system dynamics, resulting in a nearly stationary wave packet at ≈1 ps. Comparing panels (a) and (b) of Fig. 7, it can be seen that the ridge of the density surviving on the S₂ 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₂ surface are strongly suppressed and occur near the minimum geometry of the S₁ PE surface in the Q_c coordinate (see Fig. 1(b)), opposed to the isolated conical intersection, where the recurrences to the S₂ 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 first ≈30 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 specific 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 confirm 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 first 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 β-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 Redfield theory or similar models that are restricted to weak system–bath coupling and short bath relaxation time.^{19–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.⁴¹ In future extensions of the present work, the coupling of the material system with femtosecond pulsed laser fields will be included in a nonperturbative manner. While a nonperturbative treatment of the laser-matter interaction may not be a necessity for the laser fields 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 fields, which are required in the perturbative theory of nonlinear optics.⁷⁸ 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⁸² or femtosecond resonant Raman spectroscopy.⁸³ The combination of nonperturbative HEOM simulations of dissipative nonadiabatic dynamics with the calculation of spectroscopic signals directly from field-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.

Acknowledgements

Support from the Singapore National Research Foundation through the Competitive Research Programme (CRP) under Project No. NRF-CRP5-2009-04 is gratefully acknowledged. M.F.G. and W.D. acknowledge support from the Deutsche Forschungsgemeinschaft through a research grant and through the DFG-Cluster of Excellence “Munich-Centre for Advanced Photonics” (http://www.munich-photonics.de). M.F.G. thanks the members of the Zhao group for hospitality and numerous stimulating discussions. Fruitful discussions with Akihito Ishizaki are also thankfully acknowledged.

References

1. F. Bloch, Rev. Mod. Phys., 1957, 105, 1206.
2. A. G. Redfield, Adv. Magn. Opt. Reson., 1965, 1, 1.
3. G. W. Ford, M. Kac and P. Mazur, J. Math. Phys., 1965, 6, 504.
4. W. H. Louisell, Quantum Statistical Properties of Radiation, Wiley, New York, 1973.
5. A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fischer, A. Garg and W. Zwerger, Rev. Mod. Phys., 1987, 59, 1.
6. N. G. van Kampen, Stochastic Processes in Physics and Chemistry, North-Holland, Amsterdam, 1992.
7. U. Weiss, Quantum Dissipative Systems, World Scientific, Singapore, 1993.
8. W. Domcke, D. R. Yarkony and H. Köppel, Conical Intersections: Electronic Structure, Dynamics and Spectroscopy, World Scientific, Singapore, 2004.
9. W. Domcke, D. R. Yarkony and H. Köppel, Conical Intersections: Theory, Computation and Experiment, World Scientific, Singapore, 2011.
10. W. Domcke and G. Stock, Adv. Chem. Phys., 1997, 100, 1.
11. G. A. Worth and L. S. Cederbaum, Annu. Rev. Phys. Chem., 2004, 55, 127.
12. W. Domcke and D. R. Yarkony, Annu. Rev. Phys. Chem., 2012, 63, 325.
13. M. Klessinger and J. Michl, Excited States and Photochemistry of Organic Molecules, VCH, Weinheim, 1995.
14. M. A. Robb, F. Bernardi and M. Olivucci, Pure Appl. Chem., 1995, 67, 783.
15. D. R. Yarkony, Rev. Mod. Phys., 1996, 68, 985.
16. D. R. Yarkony, in Conical Intersections: Electronic Structure, Dynamics and Spectroscopy, World Scientific, Singapore, 2004, ch. 2, p. 41.
17. H. Köppel, W. Domcke and L. S. Cederbaum, Adv. Chem. Phys., 1984, 57, 59.
18. T. Gerdts and U. Manthe, Chem. Phys. Lett., 1998, 295, 167.
19. A. Kühl and W. Domcke, Chem. Phys., 2000, 259, 227.
20. A. Kühl and W. Domcke, J. Chem. Phys., 2002, 116, 263.
21. B. Balzer and G. Stock, J. Phys. Chem. A, 2004, 108, 6464.
22. B. Baltzer and G. Stock, Chem. Phys., 2005, 310, 33.
23. W. T. Pollard, A. K. Felts and R. A. Friesner, Adv. Chem. Phys., 1996, 93, 77.
24. K. Blum, Density Matrix Theory and Applications, Plenum Press, New York, 1981.
25. G. Katz, R. Kosloff and M. A. Ratner, Isr. J. Chem., 2004, 44, 53.
26. D. Gelman, G. Katz, R. Kosloff and M. Ratner, J. Chem. Phys., 2005, 123, 134112.
27. L. S. Cederbaum, E. Gindensperger and I. Burghardt, Phys. Rev. Lett., 2005, 94, 113003.
28. E. Gindensperger, I. Burghardt and L. S. Cederbaum, J. Chem. Phys., 2006, 124, 144103.
29. E. Gindensperger, I. Burghardt and L. S. Cederbaum, J. Chem. Phys., 2006, 124, 144104.
30. E. Gindensperger, H. Köppel and L. S. Cederbaum, J. Chem. Phys., 2007, 126, 034106.
31. E. Gindensperger and L. S. Cederbaum, J. Chem. Phys., 2007, 127, 124107.
32. I. Burghardt and J. T. Hynes, J. Phys. Chem. A, 2006, 110, 11411.
33. J. P. Malhado, R. Spezia and J. T. Hynes, J. Phys. Chem. A, 2011, 115, 3720.
34. J. P. Malhado and J. T. Hynes, J. Chem. Phys., 2012, 137, 22A543.
35. Y. Tanimura and R. Kubo, J. Phys. Soc. Jpn., 1989, 58, 101.
36. Y. Tanimura and P. G. Wolynes, Phys. Rev. A, 1991, 43, 4131.
37. A. Ishizaki and Y. Tanimura, J. Phys. Soc. Jpn., 2005, 74, 3131.
38. Y. Tanimura, J. Phys. Soc. Jpn., 2006, 75, 082001.
39. V. Chernyak and S. Mukamel, J. Chem. Phys., 1996, 105, 4565.
40. K.-B. Zhu, R.-X. Xu, H. Y. Zhang, J. Hu and Y. J. Yang, J. Phys. Chem. B, 2011, 115, 5678.
41. H.-G. Duan and M. Thorwart, J. Phys. Chem. Lett., 2016, 7, 382.
42. M. F. Gelin, D. Egorova and W. Domcke, J. Chem. Phys., 2012, 136, 034507.
43. A. W. Chin, J. Prior, R. Rosenbach, F. Caycedo-Soler, S. F. Huelga and M. B. Plenio, Nat. Phys., 2013, 9, 113.
44. A. Garg, J. N. Onuchic and V. Ambegaokar, J. Chem. Phys., 1985, 83, 4491.
45. Y. Tanimura, J. Chem. Phys., 2012, 137, 22A550.
46. M. F. Gelin, Y. Tanimura and W. Domcke, J. Chem. Phys., 2013, 139, 214302.
47. A. Ishizaki and G. R. Fleming, J. Chem. Phys., 2009, 130, 234111.
48. A. V. Pisliakov, M. F. Gelin and W. Domcke, J. Phys. Chem. A, 2003, 107, 2657.
49. D. Egorova, M. F. Gelin and W. Domcke, J. Chem. Phys., 2005, 122, 134504.
50. V. Stert, P. Farmanara and W. Radloff, J. Chem. Phys., 2000, 112, 4460.
51. T. Horio, T. Fuji, Y.-I. Suzuki and T. Suzuki, J. Am. Chem. Soc., 2009, 131, 10392.
52. Y.-I. Suzuki, T. Fuji, T. Horio and T. Suzuki, J. Chem. Phys., 2010, 132, 174302.
53. T. Suzuki, Molecules, 2014, 19, 2410.
54. R. Schneider and W. Domcke, Chem. Phys. Lett., 1988, 150, 235.
55. R. Schneider, W. Domcke and H. Köppel, J. Chem. Phys., 1990, 92, 1045.
56. C. Woywod, W. Domcke, A. L. Sobolewski and H. J. Werner, J. Chem. Phys., 1994, 100, 1400.
57. G. Stock, C. Woywod, W. Domcke, T. Swinney and B. S. Hudson, J. Chem. Phys., 1995, 103, 6851.
58. S. Krempl, H. Winterstetter and W. Domcke, J. Chem. Phys., 1995, 102, 6499.
59. A. Raab, G. A. Worth, H. D. Meyer and L. S. Cederbaum, J. Chem. Phys., 1999, 110, 936.
60. G. S. Latha and M. D. Prasad, J. Chem. Phys., 1996, 105, 2972.
61. U. Müller and G. Stock, J. Chem. Phys., 1997, 107, 6230.
62. G. A. Worth, H.-D. Meyer and L. S. Cederbaum, J. Chem. Phys., 1998, 109, 3518.
63. J. Y. Fang and S. Hammes-Schiffer, J. Phys. Chem. A, 1999, 103, 9399.
64. M. Thoss, W. H. Miller and G. Stock, J. Chem. Phys., 2000, 112, 10282.
65. M. Ben-Nun and T. J. Martínez, Adv. Chem. Phys., 2002, 121, 439.
66. D. V. Shalashilin and M. S. Child, J. Chem. Phys., 2004, 121, 3563.
67. P. Puzari, B. Sarkar and S. Adhikari, J. Chem. Phys., 2006, 125, 194316.
68. R. Martinazzo, K. H. Hughes, F. Martelli and I. Burghardt, Chem. Phys., 2010, 377, 21.
69. T. Grinev, M. Shapiro and P. Brumer, J. Chem. Phys., 2012, 137, 094302.
70. T. Zimmermann and J. Vanićek, J. Chem. Phys., 2014, 141, 134102.
71. M. Kanno, Y. Ito, N. Shimokura, S. Koseki, H. Kono and Y. Fujimura, Phys. Chem. Chem. Phys., 2015, 17, 2012.
72. U. Werner, R. Mitrić, T. Suzuki and V. Bonačić-Kouteckỳ, Chem. Phys., 2008, 349, 319.
73. http://arrayfire.com/ .
74. U. Manthe and H. Köppel, J. Chem. Phys., 1990, 93, 345.
75. W. Domcke, H. Köppel and L. S. Cederbaum, Mol. Phys., 1981, 43, 851.
76. M. Liebel, C. Schnedermann and P. Kukura, Phys. Rev. Lett., 2014, 112, 198302.
77. P. J. M. Johnson, A. Halpin, T. Morizumi, V. I. Prokhorenko, O. P. Ernst and R. J. D. Miller, Nat. Chem., 2015, 7, 980.
78. S. Mukamel, Principles of Nonlinear Optical Spectroscopy, Oxford University Press, New York, 1995.
79. L. Seidner, G. Stock and W. Domcke, J. Chem. Phys., 1995, 103, 3998.
80. M. F. Gelin, D. Egorova and W. Domcke, J. Phys. Chem. Lett., 2011, 2, 114.
81. M. F. Gelin, D. Egorova and W. Domcke, Acc. Chem. Res., 2009, 42, 1290.
82. J. Krčmář, M. F. Gelin and W. Domcke, J. Chem. Phys., 2015, 143, 074308.
83. B. J. Rao, M. F. Gelin and W. Domcke, J. Phys. Chem. A, 2016, 120, 3286.

---