Scholarly article on topic 'Nonequilibrium Fermi golden rule for electronic transitions through conical intersections'

Nonequilibrium Fermi golden rule for electronic transitions through conical intersections Academic research paper on "Physical sciences"

Share paper
Academic journal
The Journal of Chemical Physics
OECD Field of science

Academic research paper on topic "Nonequilibrium Fermi golden rule for electronic transitions through conical intersections"

A I ^^ The Journal of /

/All Chemical Physics .....Ij

Nonequilibrium Fermi golden rule for electronic transitions through conical intersections

Artur F. Izmaylov, David Mendive-Tapia, Michael J. Bearpark, Michael A. Robb, John C. Tully et al.

Citation: J. Chem. Phys. 135, 234106 (2011); doi: 10.1063/1.3667203 View online: View Table of Contents: Published by the American Institute of Physics.

Additional information on J. Chem. Phys.

Journal Homepage: Journal Information: Top downloads: Information for Authors:



Nonequilibrium Fermi golden rule for electronic transitions through conical intersections

Artur F. Izmaylov,1,a) David Mendive-Tapia,2 Michael J. Bearpark,2 Michael A. Robb,2 John C. Tully,1 and Michael J. Frisch3

1 Department of Chemistry, Yale University, New Haven, Connecticut 06520, USA 2Department of Chemistry, Imperial College London, London SW7 2AZ, United Kingdom 3Gaussian, Inc., Wallingford, Connecticut 06492, USA

(Received 20 September 2011; accepted 18 November 2011; published online 19 December 2011)

We consider photoinduced electronic transitions through conical intersections in large molecules. Starting from the linear vibronic model Hamiltonian and treating linear diabatic couplings within the second order cumulant expansion, we have developed a simple analytical expression for the time evolution of electronic populations at finite temperature. The derived expression can be seen as a nonequilibrium generalization of the Fermi golden rule due to a nonequilibrium character of the initial photoinduced nuclear distribution. All parameters in our model are obtained from electronic structure calculations followed by a diabatization procedure. The results of our model are found to agree well with those of quantum dynamics for a test set of systems: fulvene molecule, 2,6-bis(methylene) adamantyl cation, and its dimethyl derivative. © 2011 American Institute of Physics. [doi: 10.1063/1.3667203]


Improving performance of solar energy converting materials as well as obtaining insights in functional mechanisms of bio-molecules triggered by light absorbing chromophores would be impossible without considering radiationless transitions between different electronic states. One of the prevalent non-adiabatic features responsible for radiationless electronic transitions in large molecules of biological and technological significance is conical intersection (CI) topology.1'2 Conical intersections give rise to irreducibly quantum behavior which makes it impossible to use well-established molecular dynamics techniques based on the Born-Oppenheimer approximation. On the other hand' fully quantum treatment of electrons and nuclei is impractical for systems with a large number of degrees of freedom (DOF). Some intermediate solutions have been devised in the past decades, to name a few: mixed quantum-classical techniques3,4 (e.g., surface hopping), wave packet methods,5-7 semi-classical,8 and general path-integral based approaches.9,10 Although all these techniques alleviate burden of the full quantum consideration, they require numerical simulations and still can be quite computationally expensive owing to multiple potential energy surface (PES) calculations involved at every dynamical step. Also, the variety of time scales in large systems can easily make straightforward dynamical simulations incredibly long. In these circumstances, techniques that reduce the number of relevant DOF and parametrize PESs by simple model Hamiltonians amenable to analytical treatment are highly valuable not only because they make evaluation of non-adiabatic dynamics for large systems possible but also because these models pro-

vide qualitative insights to otherwise incredibly complex processes.

One of the simplest model Hamiltonian for non-adiabatic transitions through conical intersections is the linear vibronic coupling (LVC) Hamiltonian11

' A1 0 VW(p2+^2X2)/2 0 \

0 A2> 0 (p2+a$4)/2

ci xi xi di

a)Author to whom correspondence should be addressed. Electronic mail:

HLVC represents nuclear motion on two coupled diabatic PESs separated vertically in energy by A2 - A1 and parametrized by harmonic oscillators with frequencies coordinates and momenta p. Although both diabatic surfaces have the same set of frequencies, they differ by linear shifts df)xi and their adiabatic sets of frequencies can be different upon diag-onalization due to the off-diagonal couplings cx. In spite of its simplicity, HLVC parametrized with ab initio methods produced quite accurate results for vibronic spectra and dynamics in many molecules.1,12,13 However, the LVC Hamiltonian does not have an analytical solution and still would require quite involved quantum dynamical simulations to obtain population dynamics.14 In order to reduce the computational cost for large systems, one can consider only a restricted number of nuclear DOF which are the most relevant to the electronic transition. An ingenious and systematic approach based on a cumulant expansion of an autocorrelation function has been proposed recently for the LVC Hamiltonian.15 Of course, the reduction strategy is only possible, if relatively few effective coordinates affect non-adiabatic dynamics. Also, when some low frequency modes are involved in non-adiabatic transition, the temperature effects become important and the set of collective coordinates must be augmented with a bath which sets


135, 234106-1

© 2011 American Institute of Physics

FIG. 1. Photoinduced non-adiabatic transitions through (1) peaked and (2) sloped conical intersections. Red arrows are the ultra-fast laser excitations, purple arrows are the adiabatic oscillations on the donor surface and green arrows are the non-adiabatic transitions between the donor and acceptor surfaces.

the temperature and prevents unphysical recurrences related to a reduced number of collective coordinates. Such extension of the collective coordinates approach has been developed recently by Pereverzev and coworkers,16 and resulted in quite tractable reduced non-adiabatic dynamics within a space of electronic states and collective nuclear coordinates with a thermal bath representing the rest of the nuclear DOF.

In this work, we examine a possibility of an even simpler approach where we start with the LVC Hamiltonian and consider explicitly only dynamics of electronic populations, while all nuclear DOF are treated as a set of harmonic oscillators (bath) which stimulates electronic transitions and introduces the temperature. Within this framework, we study quite general photoinduced non-adiabatic transitions from an excited donor state through conical intersections to an acceptor state that can be either another excited state (Fig. 1, case 1) or an initial ground state (Fig. 1, case 2). Similar approaches are well known in the context of the electron and energy transfer processes, for example, the Marcus and Förster theories.17,18 In relatively rigid molecular environments where nuclear DOF are well described as a set of independent harmonic oscillators, the Marcus and Förster theories are classical limits of the Fermi golden rule (FGR) rate expressions that include thermal averages over states of harmonic oscillators. The FGR rate expression can be derived through a perturbative treatment of the off-diagonal couplings in the spin-boson Hamiltonian

Ai 0 Vv(P/2 0

0 A2) 0 P2 +«2x2)/2

'xidi}) Vc

Vc xidi

(Note that the only difference between HSB and HLVC is in the form of the off-diagonal couplings.) To be able to make first principles predictions with the FGR approach, one can obtain parameters of the spin-boson Hamiltonian from ab initio methods and use them in the FGR rate expression. This was successfully carried out in several studies for somewhat extended Hamiltonians including a Duschinsky rotation between sets of normal modes.19-21 In order to follow this strategy for photoinduced non-adiabatic transitions via conical intersections, one needs to account for the difference in non-adiabatic couplings in the spin-boson and LVC

Hamiltonians, and nonequilibrium character of the initial nuclear distribution involved in photoinduced processes. Both of these aspects were addressed in literature before but only separately. Back in the 1980s, Wagner considered linear electron-phonon couplings and derived the FGR rate expression for the equilibrium nuclear distributions.22 Recently, his results were rederived by Pereverzev and Bittner23 with further use in the context of the time-convolutionless master equation. Linear and exponential non-adiabatic couplings were also considered by Stuchebrukhov and coworkers in the context of generalizations of the Marcus theory for the electron transport in bio-molecules.24,25 As for nonequilibrium processes, theoretical models were mainly developed to treat photoinduced electron transport with the spin-boson Hamiltonian.26,27 Thus, in the current work we extend linear coupling non-adiabatic case to nonequilibrium distributions originating from a shifted ground state upon an ultra-fast photoexcitation (Fig. 1). Simplicity of the LVC model allows us to obtain electronic population of the donor state P(t) within a second order cumulant expansion as a two-dimensional integral of an analytical function F,

P (t ) = exp

f'dt' ft J0 J0

dt"F(tt"; Q,si, T)

where in addition to the LVC Hamiltonian parameters fl = Ak, a>i, df \ and ci we also introduce the shift coordinate si corresponding to the separation between minima of the ground and donor states, and the temperature T. An explicit mathematical form of F will be given by Eqs. (10)-(16), but its physical essence can be seen as an interplay of three rate-limiting factors: (1) Franck-Condon (FC) overlaps between vibrational states of the donor and acceptor diabats, (2) prefactors quadratic in couplings ci and in minima separation d(2) - d(1) originating from second order perturbational treatment of the coupling terms, and (3) phase oscillations due to a nonequilibrium character of the initial nuclear distribution. Finally, temperature can shift the balance between these three factors by changing populations of vibrational modes contributing to each of them. The regular time-independent FGR rate can be derived from Eq. (3) by starting with a stationary nuclear distribution and considering the infinite time limit. As will be shown in a couple of realistic examples, NFGR reproduces well the electron dynamics obtained in a more computationally demanding method involving propagation of frozen width wave packets. Hence, the formulated theory can be used for quick and computationally inexpensive estimates of the electronic population dynamics and its temperature dependence in large systems.

The rest of the paper is organized as follows. In Sec. II, we provide a short summary of the NFGR procedure, while more detailed derivation and discussions are given in the Appendices A and B. Section III reports details of numerical simulations and comparison of results of our approach with those of the frozen width wave packet propagation for three systems: fulvene molecule, 2,6-bis(methylene)-adamantyl, and 2-methylene-6-isopropylidene adamantyl radical cations. Section IV concludes with a discussion of

chemical insights obtained for test systems and proposes directions for future work.


Here, we present main steps of the NFGR method to obtain dynamics of electronic populations in photoinduced transitions similar to those depicted in Fig. 1.

• Adiabatic input: The minima of all involved states and the extremum of the CI seam are found using an appropriate electronic structure method. The Hessians of the donor and acceptor states are evaluated at the minimum of the donor state. The CI seam extremum is characterized by the gradients difference and non-adiabatic coupling vectors.

• Diabatization: Using the Hessians for the two states at the donor state minimum and characteristics of the CI seam minimum, the Koppel diabatization is performed (see an appendix of Ref. 28 for a detailed description). The diabatization routine provides a full second order Hamiltonian

! 2 , (1) Pi + Y,) xixj



2 , (2) Pi + Y) xixj

where xi and pi are coordinates of mass-weighted normal modes and their momenta. Hamiltonian (4) contains more parameters than is needed to construct the LVC Hamiltonian, and only those parameters which have counterparts in the LVC Hamiltonian are used. If curvatures yi(1) and y^ are different, a>i are chosen from the curvatures of the donor state (Ki(1)).29 Initial conditions: The shift parameters si are set as the difference between the donor state minimum and the ground state minimum. Rescaling: All parameters are rescaled according to

d(k) = - A-3/2

c , = ci/jlm,, i = SiM-3/2/V2,

Ak = Ak - V

^ 2ш}

This step is only a matter of notational convenience stemming from the use of the second quantization language in the derivation given in the Appendix B. 2D integration: The donor electronic population is calculated numerically as the two-dimensional integral

P (t ) = exp

fdt'f ' Jo Jo

dt"F(t', t")

F(t, t') = 2Re[/(t, t')ei(Ai-A2)(t-t')j, (10)

f (t,t') = {g(t-t') + [h(t-t')-u(t)][h(t -1') - u(t')]} x/pc(t - t')fs(t, -t'), (11)

h(t) = £ Cj ) + df)) + Cj ) - ))

x [(nj + 1)e- j - nje"a'],

g(t) = Y, Cj[(nj + 1)e-ita' + njém' ], (13)

i(t) = 2 CjSj cos (&>jt)

fs (t, t') = exJ 2i £ (df - dj2))sj [sin (j )

+ sin(«jt ')U,

/Pc(t) :

= A -I2j (d(1)-dj2))2(2nj + 1)+(dj1)-df))2[nje,fflj' +(n, +1)e-,fflj'^

and nj = (em'/T — 1) 1 are the Boltzmann populations of individual oscillators.

Three rate-limiting factors can be seen now in their explicit form: g and h2 functions have clear second order perturbational origins due to their proportionality to the coupling squares (cj). fFC introduces FC overlaps decaying exponentially with the square of the donor-acceptor distance (i?(1) — dj2))2. Finally, the terms u and fs include the shift vector components sj and represent effects of the nonequilibrium nuclear density on the perturbative prefactor [Eq. (14)] and the FC term [Eq. (15)].


All components of the adiabatic input to the NFGR procedure are obtained using electronic structure methods implemented in the Gaussian package.30 To keep computational expenses minimal but at the same time to obtain qualitatively correct topology of the electronic surfaces, all electronic structure calculations are done using the complete active space self consistent field (CASSCF) method with orbital state-averaging between intersecting surfaces and within the minimal STO-3G basis set. Quality of the basis set was not a matter of our particular concern because the main emphasis of this work is on testing the NFGR approach rather than providing results of chemical accuracy. The Koppel diabatization procedure was used as a part of the development

version of the mctdh package31 and the resulting parameters for the LVC Hamiltonian were employed in numerical integration of the F function in Eq. (9) using a simple two-dimensional trapezoidal scheme. To compare our results with more accurate methods, we also did quantum dynamical simulations employing the variational multi-configurational Gaussian (vMCG) wave packet method,7 implemented in a development version of the mctdh quantum dynamics program.31 The vMCG method can employ two types of nuclear wave-function expansions: single-set

f (k)(x,t) = £ Cf(i)xj(x, t),

and multi-set

f (k)(x,t) = £ Cjk)(t)xf(x,t),

where Xj (x, t) are the frozen width multidimensional gaussian wave packets (GWPs), and Cjk)(t) are the corresponding expansion coefficients. In the single-set expansion GWPs on two different electronic surface (k) are restricted to move identically, while in the multi-set expansion they are completely independent. Although the single-set expansion seams less flexible than its multi-set counterpart, depending on PESs of the system, in some cases, for a fixed number of GWPs, the single-set expansion can perform better than the multi-set expansion. Such situations arise when PESs of two states are quite different and initial electronic coupling is small. Then, GWPs on different surfaces can separate in space, and even when GWPs on one surface reach vicinity of CI, they will not be able to exchange population with GWPs on the other surface due to negligible spatial overlap. On the other hand, in the single-set basis this problem does not occur since for every GWP on an upper surface there is the same GWP right underneath, and both GWPs can always interact through interstate couplings. The vMCG code is also interfaced with the development version of the Gaussian package so that we can obtain quantum dynamics not only with the fixed LVC Hamiltonian but also with the fully quadratic Hamiltonian Eq. (4) evaluated on-the-fly at every propagation step. The latter Hamiltonian will be referred as the exact in the further discussions. The exact Hamiltonian within the vMCG method will allow us to assess quality of the LVC Hamiltonian separately from that of the subsequently introduced perturbation approximation in the NFGR approach. Although the developed NFGR formalism can be easily applied to any temperature, all comparisons with quantum dynamics results are done for the zero temperatures to which the vMCG method is currently restricted to.

To illustrate capabilities of our approach, we consider three systems: fulvene molecule, 2,6-bis(methylene) adamantyl (BMA), and 2-methylene-6-isopropylidene adamantyl (MIA) radical cations. This choice was motivated with the fact that the fulvene molecule and BMA cation were quite extensively studied in the past32-34 and their surface topologies are similar to those depicted in Fig. 1.

FIG. 2. Fulvene relaxation on the S1 excited state.

A. Fulvene

The ultraviolet spectrum of fulvene contains a relatively broad peak corresponding to the S0 ^ S1 transition.35 Moreover, deactivation of the S1 state proceeds without any fluorescence.35 These experimental observations were explained by finding areas of interstate crossings between S0 and S1 states.36 According to C2v symmetry, at the FC point, S0 and S1 states have A1 and B2 symmetries, respectively. Symmetry labels of adiabatic surfaces interchange after passing through the surface crossing, while diabats preserve their symmetry independent of geometry variations (at the FC point: A1 is the lower diabat and B2 is the higher diabat). Upon the S0 ^ S1 excitation, the main relaxation channel is elongation of the double CC bond connecting methylene group with the cyclopentadienyl ring (Fig. 2).32 This process initiates the radiationless transition S1 ^ S0 and engages several normal modes. All actively participating modes can be split into three categories: (1) coupling modes contributing to the non-adiabatic coupling vector (B2 symmetry), (2) tuning modes involved in the double bond elongation and contributing to the gradients difference vector (A1 symmetry), and (3) modes contributing to a torsional motion of the methylene group (A2 symmetry). Groups 1 and 2 are quite common in characterizing any system with CIs because these modes contribute to the branching space spanned by the non-adiabatic coupling and gradients difference vectors.37 Additionally, in fulvene, the torsional motion of the methylene group plays a decisive role in the CI character. Figure 3 schematically illustrates PESs involved in the radiationless transition. Here, we introduce explicitly only the torsional and the CC bond elongation coordinates, skipping the coupling coordinate but keeping in mind that every point of the intersecting seam is a CI point. As shown in Fig. 3, there are two types of relaxation paths: trajectories that do (dashed line) and do not (solid line) involve CH2 torsional displacements.32 The first type leads to sloped CIs and the second type gives rise to peaked CIs. Interestingly, the torsional coordinate has a negative curvature on the excited state as oppose to the positive curvature on the ground state. The negative curvature is quite small though, and full-dimensional simulations showed that the main relaxation path goes through the sloped intersection due to a quite steep slope along the double CC bond elongation coordinate at the FC point. The apparent shortcoming of the model LVC Hamiltonian is only one set of diabatic frequencies for both electronic states.

1. Nonequilibrium dynamics

To assess how different the LVC Hamiltonian38 dynamics is from that of the exact Hamiltonian, we simulated the vMCG wave packet dynamics with these two Hamiltoni-ans (Fig. 4). In both calculations, complete active space for

PIG. 3. Two types of conical intersections in fulvene.

CASSCF method was six electrons on six n orbitals. Previous investigations32 have found that the single-set expansion works better for this system and using 24 GWPs produce results which are very close to convergence with respect to population transfer. The two dynamics simulations are very similar for short times and agree only qualitatively for longer times. The main qualitative feature is a stepwise population decay due to a nonequilibrium character of the initial nuclear distribution on the upper B2 diabat and, as a consequence, multiple recrossings of the sloped CI. The main quantitative difference between dynamics on the exact and LVC surfaces is in the amount of population decay at every CI recrossing. Thus, we conclude that although the LVC Hamiltonian requires some extensions for quantitative agreement with the exact Hamiltonian dynamics, it reproduces very well qualitative features of non-adiabatic dynamics.

A comparison of the NFGR and vMCG dynamics within the LVC Hamiltonian (Fig. 5) reveals very good agreement

FIG. 4. vMCG dynamics with the exact (red) and LVC (blue) Hamiltonians. Both simulations were done with 24 single-set GWPs.

in the initial dynamics and in slopes of transitions as well as overall stepwise decay of the population. It is somewhat harder for NFGR to reproduce the amount of population dropped at each CI recrossing. We attribute this to considering only the forward population flow (donor to acceptor) and neglecting the backward population flow (acceptor to donor) in NFGR. Switching to the generalized master equation (GME) approach would be a possible strategy to introduce the backward population flow. However, it was found that a regular time-convolutionless version of the GME approach provides almost indistinguishable results from those of the NFGR method in fulvene. This means that the acceptor-donor population flow is negligible within the GME formalism. The reason for this underestimation is a GME assumption that nuclear density on the acceptor side can be taken in the same form as on the donor side. In fulvene, the energy difference between donor and acceptor minima is 2.7 eV, and therefore, upon the electronic transition the nuclear dis-


^ vMCG

FIG. 5. NFGR (black) and vMCG (blue, 24 single-set GWPs) population dynamics with the LVC Hamiltonian.

FIG. 6. Population dynamics of NFGR with the LVC Hamiltonian (black) and vMCG with the exact Hamiltonian (blue, 24 single-set GWPs). Red dashed lines illustrate slopes of the NFGR population stepwise decay. Each slope line is duplicated and shifted in parallel fashion to be compared with a corresponding vMCG decay slope.

tribution of the acceptor side is very different from that of the donor side. This difference is considered only perturba-tively in GME which is clearly insufficient for quantitative agreement in this case. One resolution of this issue would be to include the difference between nuclear distributions of the donor and the acceptor in some non-perturbative manner. The work on such a method is ongoing and results will be reported elsewhere. Nevertheless, it is quite amazing that such a simple method as NFGR can obtain the correct slopes of the population decay. Moreover, this is still true if we compare the NFGR results with those from the vMCG dynamics with the exact Hamiltonian (see Fig. 6).

2. Temperature dependence

Figure 7 shows temperature dependence of the nonequi-librium energy transfer process in fulvene. According to Eq. (11), there are three possible sources of the temperature dependence of the population dynamics: pre-exponential

t\ 0 K -

|r 500

I 1000 KM -


2000 KV "

°0 10 20 30 40 50

FIG. 7. Temperature dependence of nonequilibrium electronic population decay in NFGR for fulvene.

FIG. 8. The lowest energy coupling mode in fulvene.

Boltzmann populations of individual oscillators in the h and g functions [Eqs. (12)and(13)] as well as the Boltzmann populations in the exponent of the /FC function [Eq. (16)] that gives rise to the temperature dependence in the fully quantum analog of Marcus/Förster theory.17,39 Owing to symmetry restrictions in fulvene, all tuning modes have c ~ 0 and all coupling modes have df* ~ 0; hence, h function has a negligible contribution for any temperatures. By simply switching on and off the temperature dependence in the g function, we confirmed that the overall temperature dependence is dominated by the g prefactor in fulvene. Temperature increase creates higher Boltzmann populations of coupling modes and, thus, promotes the electronic transition. The lowest energy coupling mode in fulvene has a frequency of 902 cm-1 and corresponds to an in-plane ring deformation of B2 symmetry (see Fig. 8). Thus, the dynamics does not respond appreciably to the temperature increase until 1000 K when the lowest energy mode becomes activated (see Fig. 7).

B. 2,6-bis(methylene) adamantyl radical cation and its dimethyl derivative

The 2,6-bis(methylene) adamantyl and 2-methylene-6-isopropylidene adamantyl radical cations represent quite rigid molecular structures where intramolecular electron transport (IET) can occur between two unsaturated fragments separated by the adamantane cage (see Fig. 9).33,34 The energetic profile of this IET is depicted in Fig. 10. In contrast to fulvene, the adiabats have the same symmetry for all geometries 2B2 for ground and 2B1 for excited states, while the diabats change their symmetry after the crossing from 2B2 as a lower state to 2B1 as an excited state. In both systems, regions of conical intersections have been found between two C2v minima, Loc and Loc' (Fig. 9). In BMA, the minimum of the CI seam has D2d symmetry, and therefore, this system is an example of Jahn-Teller distortion where energy of the high symmetry

FIG. 9. Intramolecular electron transport in the 2,6-bis(methylene) [X = H] and 2-methylene-6-isopropylidene [X = CH3] adamantyl radical cations.

FIG. 10. Schematic adiabatic (3D) and diabatic (2D) energy profiles for BMA (both wells are equal) and MIA (the donor well is higher in energy).

FIG. 11. Equilibrium population dynamics in BMA: NFGR (solid black), vMCG 24 (solid red) and 47 (solid blue) multi-set GWPs with the LVC Hamiltonian, and vMCG 8 (dashed magenta) and 16 (dashed green) singleset GWPs with the exact Hamiltonian.

configuration (the minimum of the CI seam) is reduced by lowering the symmetry (the Loc and Loc' minima). The IET active modes in both systems are tuning modes of A1 symmetry and coupling modes of A2 symmetry. The tuning modes, as expected, mostly contribute to the antisymmetric stretching of the two n bonds, while the coupling modes involve an antisymmetric breathing of the rigid adamantane skeleton (see Fig. 10).33 In the adiabatic picture, the reaction coordinate is aligned mostly with the gradients difference vector until in the vicinity of the CI seam. Then, the system will try to overcome the minimum of the CI seam sideways and it is important for the non-adiabatic state coupling coordinate to be large enough because it is the one which creates a gap between electronic states. As it was found for both systems, the coupling vector norm is relatively small, and therefore, the resulting topology is similar to a simple curve crossing within (N — 1)-dimensional manifold, where N is the number of nuclear DOF. A small coupling vector results in relatively weak diabatic couplings between states and gives rise to the so-called "diabatic trapping" behavior when the system, after reaching the intersection of diabatic curves, keeps following the donor diabat rather than switching to the lower energy acceptor diabat. Small diabatic couplings in BMA and MIA make the NFGR perturbative treatment particularly attractive for these systems. Also, the quite rigid cage of adamantane makes the harmonic oscillator approximation very accurate for most of the modes. In addition, in BMA, the symmetry makes harmonic frequencies on different diabats exactly the same.

1. Equilibrium dynamics

Before considering more complex nonequilibrium dynamics in BMA and MIA systems, we would like to illustrate some differences in equilibrium (s = 0) IET processes for these systems. In both cases, to solve the CASSCF equations, we used the active space of three electrons within four n orbitals. Comparing the exact Hamiltonian vMCG dynamics

with that from the NFGR method reveals very good agreement but only for short times (see Figs. 11 and 12). Owing to presumably much better agreement between the LVC and exact Hamiltonians for these systems, we attribute discrepancy in the long time dynamics mainly to a finite number of GWPs used in the exact Hamiltonian vMCG dynamics. This assumption is numerically confirmed by the fact that the long time discrepancies can be postponed to longer times with increasing the nuclear basis set in both systems. However, systematic assessment of the LVC approximation for longer times would require further extensions of the nuclear basis sets which becomes prohibitively computationally expensive for the exact Hamiltonian vMCG calculations. Thus, we can only conclude that the LVC Hamiltonian introduces only very small deviations in dynamics for time scales where the employed nuclear basis sets are adequate. The NFGR rate expression uses analytic form for the trace element and, thus,

V 4 gwp

v * * 8 gwp

\ -""V ^^ s.*s

^^^ ------—

47 gwp\s. NFGR

29 gwp^N^

0 10 20 30 40

FIG. 12. Equilibrium population dynamics in MIA: NFGR (solid black), vMCG 29 (solid red) and 47 (solid blue) multi-set GWPs with the LVC Hamiltonian, and vMCG 4 (dashed green) and 8 (dashed dotted magenta) single-set GWPs with the exact Hamiltonian.

represents the complete basis set limit. In order to explore performance of the NFGR approach at longer times, we extended nuclear basis sets in our vMCG calculations to the multi-set expansion (which converges better for these systems) with increased number of GWPs. To keep the computational expenses feasible, we needed to switch to the LVC Hamiltonian38 and to reduce considered DOF to the IET active modes. This reduction procedure resulted in 23- and 28-dimensional LVC models for BMA and MIA, respectively. It also allowed us to carry out vMCG calculations with up to 47 GWPs on each diabatic surface and to illustrate very good agreement between NFGR and vMCG results even for longer times. In addition, Figs. 11 and 12 show that the electron transfer in MIA proceeds almost one order of magnitude faster than that in BMA. Further analysis of the rate components [Eqs. (10)-(16)] revealed the microscopic origins of this difference. First, terms containing products Cj (dj1-1 ± dj2)) have a negligible contribution due to a strict symmetry related separation of all modes to tuning (Cj ~ 0) and coupling

0) modes. Second, although the main contribution to

rates comes from j containing terms h(At)fFC(At), coupling strengths j Cj- in h(At)] for BMA and MIA have very similar values. Third, the FC prefactor fFC(At) is much larger for MIA than for BMA because of smaller reorganization of states £ j (dj:) - dj2))2 in MIA compared to BMA. Thus, the main effect of the methyl groups on IET kinetics is not the lowering of the Loc' energy minimum due to charge stabilization (Fig. 9) but rather reducing the separation between the Loc and Loc' minima in space. From the structural point of view, the stabilization of charge in the Loc' configuration of MIA can also be seen as making oscillators in two minima more similar and, thus, facilitating the non-adiabatic transition through increasing the FC overlaps between corresponding nuclear wave-functions.

2. Nonequilibrium dynamics

To investigate nonequilibrium dynamics, we consider photoexcitation from the minimum of the acceptor diabat (Loc' in Fig. 9) to the donor diabat (Fig. 10). Therefore, the shift coordinate corresponds to the gradients difference vector Sj = dj^ — dj2). For the vMCG method, the nonequilibrium case is even more basis set demanding than the equilibrium case: All our calculations with the exact Hamiltonian and small basis sets (up to 12 GWPs) provided quite inadequate results because of basis set limitations. Therefore, here we present only vMCG results with the LVC Hamiltonian for the reduced models where we can approach convergence with the basis set. Figures 13 and 14 illustrate the electron transfer dynamics in the NFGR and vMCG methods for BMA and MIA. NFGR performs very well in describing oscillations related to the nonequilibrium character of the initial distribution for quite extended time frame. Due to some differences in PESs, the MIA dynamics is even more demanding in terms of the nuclear basis than that of BMA. Figure 14 illustrates that basis deficiency can create dynamical artifacts in quite arbitrary places (cf. vMCG results with 29 and 40 GWPs), and it is highly valuable to have the complete basis set point

FIG. 13. Nonequilibrium population decay in BMA: NFGR (dashed black) and vMCG with the LVC Hamiltonian and 24 multi-set GWPs (solid red).

of view of the NFGR method, even though the latter has perturbative origins.

Interestingly, the stepwise decay that persists in fulvene and BMA almost completely disappears in MIA after 40 fs. This is a manifestation of the dephasing which takes place in MIA because of a larger number of DOF contributing to the shift coordinate. To see a connection between the population dynamics and the dephasing in a more transparent way, let us explore the evolution of the coherence along the shift coordinate

Trn [aspn(t)] =

J2- s2e

2 eitojt j "j

where as = Xj aj~sj/y Xj s2, and pn(t) represents dynamics of the nonequilibrium nuclear distribution (see the Appendix A). Since we consider the non-adiabatic transition as a perturbation, the pn(t) dynamics is unaffected by this transition and is governed solely by forces on the donor surface. Figure 15 shows dynamics of the real part of the coherence function for several choices of the shift coordinate.

\ \\ \\ \ N5. N X

\ \ v \ \ V>-- \ \ \ ^N \ 29 gWp^-

\ _ \40 gwp^-


20 40 60 80 100 120 140 t, fs

FIG. 14. Nonequilibrium population decay in MIA: NFGR (dashed black) and vMCG with the LVC Hamiltonian and 29 (solid red) and 40 (solid blue) multi-set GWPs.

FIG. 15. Coherence oscillations for several shift coordinates: (1) reaction j -

dom components (dashed black).

d. ) — iT. (solid blue), (2) single mode (solid red), and (3) ran-

FIG. 7. Lowest energy coupling modes in BMA and MIA.

modes involved in the shift coordinate, after 1 00 fs the coherence oscillations reappear in Fig. 1 5, reappearance of the related stepwise population decay can be also seen for longer times in Fig. 14.

Along with the nonequilibrium dynamics where Sj = j — j (reaction mode), we plot cases when only the largest single component of Sj is left, while all others were zeroed (single-mode shift), and when all Sj components were assigned random values in the same range as in real MIA. Corresponding population dynamics is given by Fig. 16. In general, the timescale of the large amplitude coherence oscillations in Fig. 15 matches the timescale of the stepwise population decay in Fig. 16. In the single mode case, the stepwise population decay persists without visible reduction in amplitude since there is no interference between different modes contributing to the shift coordinate. In the randomly assigned case, the population decay appears very much as in the equilibrium case because all nonequilibrium contributions cancel each other from the beginning. When we shift the nuclear distribution along the reaction mode, the situation becomes somewhat intermediate between the single mode and random cases. Initially, the decay of the coherence oscillations corresponds to the disappearance of the stepwise decay after 40 fs. Due to a still relatively small number of individual

3. Temperature dependence

Finally, we would like to discuss temperature dependence of non-adiabatic dynamics in BMA and MIA. Figures 18 and 19 show temperature dependence of the equilibrium IET process in the BMA and MIA cations. As in the fulvene case, among all temperature related contributions the g prefactor [Eq. (13)] dominates in both systems, and thus, temperature increase always promotes the transition by activating coupling modes. Temperature dependence in BMA mostly comes from 435 cm-1 coupling mode corresponding to a torsional motion of the methylene groups and involving deformations in the adamantane cage (see Fig. 17). In MIA, temperature effects appear at much lower temperatures than in BMA because the lowest coupling mode has frequency of only 32 cm-1, and has significant population at much lower temperatures. This mode corresponds to an antisymmetric rotational motion of the methyl groups (see Fig. 17) and, therefore, has hindered rotor character. The height of the rotational barrier is 723 cm-1 which should be enough to keep our harmonic oscillator consideration valid at least up to 500 K.

If we ignore the CI topology of the problem and assume some constant diabatic coupling instead of c\X\ terms,

FIG. 16. Nonequilibrium population decay in MIA for several shift coordinates: (1) reaction mode Sj = j1 — dj2 (solid blue), (2) single mode (dashed red), and (3) coinciding random components and equilibrium Sj = 0 results (dashed black).

FIG. 18. Temperature dependence of equilibrium population decay in NFGR for BMA.

FIG. 19. Temperature dependence of equilibrium population decay in NFGR for MIA.

then the temperature dependence of the rate expression originates solely from the fFC function. This assumption is equivalent to considering the fully quantum version of the Marcus/Förster theory.17'39 In BMA, raising the temperature increases the fFC defined rate, while in MIA the trend is the opposite because the minimum of the donor state is higher in energy than the minimum of the acceptor state (endother-mic reaction). Interestingly, adding the correct CI topology to the problem not only enhances the temperature dependence in BMA but it also reverses the rate response to the temperature in MIA.


We derived a nonequilibrium formalism for electronic transitions through conical intersections which is rooted in a very simple perturbative consideration. Starting with the LVC Hamiltonian [Eq. (1)], we treat all nuclear DOF as a harmonic bath that stimulates electronic transitions between the donor and acceptor states. Nonequilibrium initial conditions in nuclear DOF introduced in NFGR allows us to investigate non-adiabatic transitions induced by ultra-fast laser pulses. The main advantage of the NFGR approach is in ability to substitute the computationally expensive quantum dynamics step on the road from electronic structure calculations to characterization of the electronic transition process. Owing to computational simplicity of NFGR, we can shift the emphasis of the problem to the electronic structure part and to investigate large systems with more accurate ab initio methods including larger basis sets and dynamical correlation effects. Nuclear basis set completeness of the NFGR approach makes it also a good compliment to the frozen width wave packet techniques to asses their inevitable basis set incompleteness. Another advantage of the proposed approach is to quickly relate the structural features of the system to components of the rate expression. This capability can be used not only to obtain qualitative insights in electronic transitions in large nanoscale systems with the CI topology but also to perform the inverse design of such systems. As has been illustrated in BMA and IMA examples, temperature effects can be very important when low frequency modes have significant

contributions to the non-adiabatic coupling vector. Therefore, in order to obtain realistic estimates of the non-adiabatic rates, it is crucial to include temperature effects in simulations of large systems, and the NFGR approach does exactly that. Moreover, the temperature dependence through the simple Förster/Marcus exponential factor can be insufficient or even misleading in the presence of the conical intersection topology.

Using the precalculated LVC Hamiltonian and perturba-tional treatment of electron-nuclear dynamics makes computational cost of our approach negligible, but at the same time, these two features can be seen as weaknesses suggesting directions for further improvements of NFGR. As has been shown, the LVC Hamiltonian can provide quantitative agreement for short time dynamics and qualitative agreement for longer times. Interestingly, such a short time feature as the time scale of electronic transitions when nonequilibrium wave packet approaches states intercrossing is well reproduced even at longer times (Figs. 4 and 6). Without resorting to on-the-fly evaluation of the nuclear Hamiltonian, the LVC Hamiltonian can be somewhat extended by adding more flexibility through additional terms, for example the Duschin-sky rotation, to allow different surfaces to have different normal modes and frequencies. This extension will complicate the analytical function F in Eq. (3) but will preserve the integral form of the population expression. Also, from the computational point of view, an extended formalism will still require only two ab initio calculations: in the FC point and in the minimum of the CI seam. In a similar fashion we can include some anharmonic corrections, but one could envision that resorting to dynamically updated Hamiltonians will be necessary if anharmonic effects are strong. As for perturba-tional character of NFGR, it captures with quantitative agreement dynamics in systems with weak electronic couplings (BMA and IMA) and reproduces very well qualitative features of electron nonequilibrium dynamics in the strongly coupled case of fulvene. In the latter case, NFGR overestimates amount of population decay at each CI recrossing. Remaining in perturbative picture, this problem can be connected with the absence of the backward flow from the acceptor to donor. As we discussed, the main problem with conventional master equation approaches, which are usually employed to include the backward population flow, is their inability to treat highly nonequilibrium distributions appearing on the acceptor surface of fulvene. In such cases, one does need either to add the most nonequilibrium nuclear DOF to a set of DOF considered explicitly (e.g., approach taken in Ref. 16) or to modify conventional master equation approaches by including projectors on nonequilibrium states.


AFI appreciates stimulating discussions with Alexander Kandratsenka and Graham Worth. vMCG calculations were carried out at the Imperial College High Performance Computing Service and Yale University Faculty of Arts and Sciences High Performance Computing Center. AFI and DMT are grateful for support granted by Gaussian Inc. JCT is

grateful for support from the National Science Foundation, CHE-0909730.


The electronic population dynamics for the LVC Hamiltonian [Eq. (1)] can be written by projecting the full electron-nuclear density dynamics p(t) = e-lHLVC'p(0)elHLVC' on the electronic donor state ( 11 > ) and tracing out the nuclear DOF (Trn)

der population is

P(t) = Trn [(1|e-lHLVCtp(0)elHLVC111>]

We will use perturbation theory with respect to the electronic off-diagonal couplings од. However, first, in order to use the same harmonic oscillator states for both diabats, we remove on-diagonal linear terms df)xi in the LVC Hamiltonian with the polaron transformation (U)17 H = UJ,HLVC U, where U = , k enumerates electronic states (k = 1

donor and k = 2 acceptor), and zjk) = idf) /a?- are necessary geometric shifts. The polaron transformation results in the Hamiltonian where diagonal harmonic oscillators are not shifted anymore; instead, the off-diagonal terms are dressed by exponential operators

H = Ài|1)<1|-À2|2)<2| + (|1)<1| + |2)<2|)£ (p2+a\x})/2

|1)(2|e-£ izip'[J2cixi

+ |2>(1|e-(£<

V^ (!)

and the constant energy shifts undergo renormalization Дк = Ак j

"i i /(2 j. The dressing of the off-diagonal coupling corresponds to the introduction of the FC factors in the operator form. Following the conventional perturbation theory procedure, we define H = H0 + V, where H0 contains all operators with diagonal electronic projectors, and V incorporates both off-diagonal operators: V12 for the |1)(2| block and and V21 for the |2)(1| block. After the polaron transformation, the donor population dynamics [Eq. (A1)] becomes

P(t) = Trn [(1 \e—iHtUtp(0)UeiHt\1>] . (A3)

Expanding the evolution operator up to the second order and introducing the density factorization into the electronic (\ 1 >( 1 \) and nuclear (pn) parts give

P (t) = Trn [UlUfV )\1>p„(0)(1\U j2)t(t )\1>], (A4)

U ?\t ) = 1 -

fdt'f ' J0 Jo

dt"eiHat ' Ve


Note that the population in Eq. (A4) contains some parts of the overall fourth order corrections, while the pure second or-

P(2)(t) = 1 - 2Re f dt' f dt"(1|Trn[p„(t")

J0 J0 o-lH0(t '-t ")


V21] |1>.

Here, we used the trace invariance with respect to cyclic permutations. The last expression can also be connected to the time-independent FGR. To illustrate this, we rewrite P(2)(i) as

P(2)(t) = 1 - f dt'K(t'), h

K(t') = 2Re / dt"(1|Trn[p„(t")elH°(t'-t") x V12e-lHo(t'-t")V2l]|1>

is a time-dependent analog of a rate constant. Equation (A8) gives the NFGR time-dependent rate. The regular time-independent FGR rate can be derived from the NFGR rate under few additional assumptions: First, nuclear density pn(0) is taken as the Boltzmann density (pBZ), so that Pn(t") = pBz,40 and second, we consider times t which are longer than a characteristic decay time of the correlation function (1\Trn[pBZV12(i' — i")V21]\1> so that one can substitute t in the upper integral limit with infinity and obtain the time-independent rate. The latter condition requires the presence of many degrees of freedom, which is usually the case in condensed phase environments.

Population dynamics of Eq. (A7) becomes negative at longer times if the spectrum of nuclear DOF in the donor state overlaps with that of the acceptor state. A relatively easy way out of this problem is the so-called exponential resum-mation of the perturbation series, which is also known as the cumulant expansion.41 The main idea of this resummation is to do a perturbative expansion with preserving the exponential form of the evolution operator e—iHt. For a general exponential operator eA, the cumulant expansion of a thermal average is Tr[peA] = Tr[pBj], where first two terms in the expansion are B1 = Tr[pA] and B2 = Tr[pA2] — Tr[pA]2. Here, we substitute the thermal average of the exponential operator by the exponent of the thermal averages. Although, practically, we still need to truncate the cumulant expansion at some finite order; this expansion allows us to preserve the exponential structure in any finite order. In our case, Eq. (A3) can be rewritten as

P(t) = Trn [p„(0)(1|elHt11 >( 11e lHt|1>]

= eT,j Trn[pn(0)Bj]

and in the second order cumulant expansion, the donor population dynamics becomes

P (2)(t ) = e-/odi'K(i').

This expression is always positive and does not require any extra computational work compared to Eq. (A7). Equations (A7) and (A11) can also be seen as a consequence of the zeroth and first order kinetics for the donor population. To

make Eq. (A8) amenable to numerical calculations, next, we address the evaluation of the thermal averages involved in the time-dependent rate expression.


For better understanding of the physics involved in the time-dependent rate Eq. (A8), before considering a nonequilibrium situation with the linear couplings we will discuss two simpler cases: constant and linear couplings with the equilibrium nuclear density p„(t') = pBZ. In all cases, the time-dependent rates will be expressed as

K(t) = f dt'f(t,t')ei(Ä Jo

!-Ä2)(t-t ' )

where the electronic ei(Al A2)(t ) and nuclear f(t, t') contributions are separated.

a. Constant coupling with Boltzmann density In the constant diabatic coupling case, we have a fully quantum version of the Marcus and Förster theories with VC instead of XiCiX [see Eqs. (1 ) and (2)]. The nuclear contribution to the time-dependent rate in this case is

f (t, t') = V2cTrn[pBzeiHs(t-t)

x e^;^-^Fie-^B(t-t')e^jbf-zfp], (B2)

where HB = X j (f2 + x])/2 is the harmonic bath Hamiltonian. f(t, t) is the correlation function which can be seen as electronic transition promoting force originating from the harmonic bath fluctuations. In order to evaluate the trace in f(t, t), it is easier to work in the second-quantization notation

f (t, t') = VCTrn{pBze^i(t-t')-^(t-t'}, (B3)

ioj(t-t0), and

Oioj (t-t)

— a;e

'—j j j J

(k) = a(V„ -3/2д/2 are scaled displacements. Here, to

where Qk (t — t') = X j djC)(a)e d?) = —d(k)^ ; eliminate e±lHB(t-t ) factors of Eq. (B2), we used the follow' ing identity exajal(a\ a)e—xa]a = l(a^ex, ae—x) valid for any analytical function /.42 Further simplifications in Eq. (B3) are possible because all exponents are linear in creation and annihilation operators. For such exponential operators, we use the eAeB = eA + Be-[A,B]/2 identity43 followed by the Bloch theorem17 Tr [peA+B] = eTr[p(A+B)2]/2. With these two steps, evaluation of the final quadratic thermal average becomes relatively straightforward

f (t, t') = V2fFc(t - t'),

fpc(t -1')

= ei^j j

!))2[nje"J(t-,'4(nj +ï)e-i»J

where nj is the thermal Boltzmann population of the jth harmonic oscillator with the frequency The exponential decay in f(t, t') with the nuclear separation between minima of two electronic states (dj1) - dj2))2 can be related to the FC factors we introduced in operator form in Eq. (A2). Thermal averaging in Eq. (B3) makes various FC factors contribute

to f(t, t') according to their thermal weights given by the Boltzmann populations. In spite of two time arguments in f(t, t), it depends only on the time difference At = t - t, and thus, f(t, t') describes only one time-scale corresponding to the non-adiabatic transition. The time dependence of fFC(t - t) is relatively cumbersome in the fully quantum form. However, it is possible to simplify it, for example, in the zero temperature limit,17 and to show that the main contributions to the final rate expression [Eq. (B1)] come from resonance (iso-energetic) transitions between two states. In other words, vibrational levels of the donor which have iso-energetic vibrational levels of the acceptor surface contribute the most. Another limit where Eq. (B4) can be simplified is the high temperature classical limit; this is the limit where Eq. (B4) recovers the Marcus/Förster expression.17

b. Linear coupling with Boltzmann density The time-dependent rate for the linear coupling case is a correlation function with two additional linear terms [see Eq. (B3)]

f (t, t') = Trn{pBz^l(t-t,)L(t - t')e-Q2(t-'')ea2Le-^ },

where L(t — t') = Xj cj(ajelmj(t-t) + aje cj = cj/y!2mj. The analytical treatment of this correlation function uses thermal Wick's theorem44 that helps to reduce the linear operators L and L(t — t') to some time-dependent functions. We will skip the relatively lengthy and tedious application of Wick's theorem to Eq. (B6) and only state the final result here38

(B6) ioj(t-t 0) and

f (t,t') = [h(Ät )2 + g(Ät )]fpc(Ät ),

h(t) = Y,j J+J) + (j - j)

x [(n j + 1)e

g(t) = Ц j[(nj + 1)e-"aj + nje"Wj].

- nJeltwJ ],

Note that as opposed to the constant coupling case, the time dependence of the prefactor (h2 + g) modifies the iso-energetic (resonance) condition originating in /FC(At). This can be seen as inclusion of inelastic scattering of electrons: when an electronic state changes, there is an additional change of vibrational nuclear states encoded in Eq. (B7). We can further split the prefactor contribution into three parts: (1) the constant part of h, (2) the time-dependent part of h, and (3) the g function. The time-independent part of h ~ cj(j + dj2)) increases with the asymmetry of the electronic state minima with respect to the point of CI; here, the asymmetry is equivalent to the presence of the constant coupling VC = J2j Cj(dj1) + dj2)). The time-dependent part j(2)\

of h ~ cj (dj1) — dy2') naturally increases with the separation between the minima because of the coupling growth. The

dj' )2(2n j + 1)+(dj

contribution from the g function is a manifestation of a de-localized nature of quantum nuclei: even without any separation between the electronic state minima (¿?j1) = j = 0), the coupling is nonzero due to finite widths of nuclear wave-functions.

c. Linear coupling with shifted Boltzmann density In photoinduced non-adiabatic dynamics, we assume a vertical electronic transition (e.g., induced by femtosecond laser pulse) which puts the Boltzmann density distribution on the shifted donor state (see Fig. 1). This shifted density distribution can be written as ps = "(aj+sj)(ai +Si)/Trn[e-^j"(aj +sj)(ai +Si)], where sj are geometrical shifts between the ground and donor states. The correlation function we need to evaluate here is

f (t, t') = Trn{ps(t,)eQl(t-t')L(t - t')e-Q2(f-f)eQ2LeUl},

where ps (t) = "(aj e-"i' + sj )(a'e""'' + sj) /Trn[e-^" (aje ""i' +Sj)(aje""i'+ss)]. To proceed, we perform canonical transformations bj = aje-l"jt + Sj and bj = aje1"'1 + Sj that make the nonequilibrium density ps(t) algebraically equivalent to the previously considered Boltzmann density. This allows us to use the same operator contractions in thermal Wick's theorem as in the equilibrium case. Some differences in the rest of the correlation function expression only slightly increase the length of the derivation without requiring any additional algebraic tricks.38 Our key result is

f (t, t') = {g(At) + [h(At) - u(t)][h( At) - u(t')]}

x fs(t, -t/)fFc(t,t/), (B11)

u(t) = 2 CjSj cos ("jt), (B12)

fs (t,t') = exp

2l^ (dj1) - dj2)) sj [sin (j ) + sin (j ')]

Comparing with the equilibrium case, the inclusion of nonequilibrium effects manifests in two time modulated scalar products of the shift coordinate Sj with the coupling [cj, Eq. (B12)] and reaction [j — dj2), Eq. (B13)] coordinates. Scalar products account for the coupling growth when the shift is along the coupling modes and the coupling reduction (enhancement) and when the shift is in the opposite (same) direction with respect to the reaction coordinate (see Fig. 20). Time modulation of these scalar products comes from collective character of the shift coordinate where every contributing oscillator potentially has a different frequency. In cases when a large number of oscillators contribute to the shift coordinate, after a certain time, contributions from different oscillators will cancel each other and scalar product values will become negligible. This dephasing reduces Eq. (B11) to its equilibrium counterpart Eq. (B10). In cases when only a

FIG. 20. Schematic energy profiles illustrating the difference between excitations (1) further from and (2) closer to the intersection.

relatively small number of modes contribute to the shift vector, the scalar products can oscillate in time indefinitely producing oscillations in the rate. The time dependence of the scalar products also breaks the translational time invariance we had in all equilibrium cases. This is a manifestation of two time-scales present in the nonequilibrium case: first, the regular non-adiabatic transition time-scale, and second, the time-scale of the vibrational state phase oscillations related to the nonequilibrium character of the nuclear distribution.

1 Conical intersections: Electronic Structure, Dynamics and Spectroscopy, edited by W. Domcke, D. R. Yarkony, and H. Koppel (World Scientific, Singapore, 2004).

2B. G. Levine and T. J. Martinez, Annu. Rev. Phys. Chem. 58, 613 (2007). 3J. Tully, J. Chem. Phys. 93, 1061 (1990). 4R. Kapral and G. Ciccotti, J. Chem. Phys. 110, 8919 (1999). 5M. Ben-Nun, J. Quenneville, and T. Martinez, J. Phys. Chem. A 104, 5161 (2000).

6Multidimensional Quantum Dynamics, edited by H.-D. Meyer, F. Gatti, and G. A. Worth (Wiley-VCH, Weinheim, 2009).

7G. A. Worth, M. A. Robb, and I. Burghardt, Faraday Discuss. 127, 307 (2004).

8M. Thoss, W. H. Miller, and G. Stock, J. Chem. Phys. 112, 10282 (2000).

9A. Piryatinski, M. Stepanov, S. Tretiak, and V. Chernyak, Phys. Rev. Lett. 95, 223001 (2005).

10J. R. Schmidt and J. C. Tully, J. Chem. Phys. 127, 094103 (2007).

Atomic units are used throughout this work. 12H. Koppel, W. Domcke, and L. S. Cederbaum, Adv. Chem. Phys. 57, 59 (1984).

G. A. Worth and L. S. Cederbaum, Ann. Rev. Phys. Chem. 55, 127 (2004).

14H. Wang and M. Thoss, J. Chem. Phys. 119, 1289 (2003).

15L. S. Cederbaum, E. Gindensperger, and I. Burghardt, Phys. Rev. Lett. 94, 113003 (2005).

16A. Pereverzev, E. R. Bittner, and I. Burghardt, J. Chem. Phys. 131, 034104 (2009).

17A. Nitzan, Chemical Dynamics in Condensed Phases: Relaxation, Transfer, and Reactions in Condensed Molecular Systems (Oxford University Press, New York, 2006).

18R. Zwanzig, Nonequilibrium Statistical Mechanics (Oxford University Press, New York, 2001).

19Q. Peng, Y. Yi, Z. Shuai, and J. Shao, J. Chem. Phys. 126, 114302 (2007).

20R. Borrelli, M. D. Donato, and A. Peluso, J. Chem. Theory Comput. 3, 673 (2007).

R. Borrelli and A. Peluso, Phys. Chem. Chem. Phys. 13, 4420 (2011).

22M. Wagner, J. Phys. C 15, 5077 (1982). 23A. Pereverzev and E. R. Bittner, J. Chem. Phys. 125, 104906 (2006). 24I. Daizadeh, E. S. Medvedev, and A. A. Stuchebrukhov, Proc. Natl. Acad.

Sci. U.S.A. 94, 3703 (1997). 25E. S. Medvedev and A. A. Stuchebrukhov, J. Chem. Phys. 107, 3821 (1997).

R. D. Coalson, D. G. Evans, and A. Nitzan, J. Chem. Phys. 101, 436 (1994).

27M. Cho and R. J. Silbey, J. Chem. Phys. 103, 595 (1995). 28C. S. M. Allan, B. Lasorne, G. A. Worth, and M. A. Robb, J. Phys. Chem. A 114, 8713 (2010).

29This choice is one of the simplest among various judicious choices. As an ultimate solution for the parametrization problem, we consider building a rate expression for the full Hamiltonian Eq. (4) in future.

30M. J. Frisch, G. W. Trucks, H. B. Schlegel et al., gaussian development version, Revision H.01, Gaussian, Inc., Wallingford, CT, 2009.

31 G. A. Worth, M. H. Beck, A. Jackle, and H.-D. Meyer, mctdh development version 9.0, University of Heidelberg, Heidelberg, Germany, 2009.

32D. Mendive-Tapia, B. Lasorne, G. A. Worth, M. J. Bearpark, and M. A. Robb, Phys. Chem. Chem. Phys. 12, 15725 (2010).

33L. Blancafort, P. Hunt, and M. A. Robb, J. Am. Chem. Soc. 127, 3391 (2005).

34L. Blancafort, F. Jolibois, M. Olivucci, and M. A. Robb, J. Am. Chem. Soc. 123, 722 (2001).

35J. E. Kent, P. J. Harman, and M. F. O'Dwyer, J. Phys. Chem. 85, 2726 (1981).

36M. J. Bearpark, F. Bernardi, M. Olivucci, M. A. Robb, and B. R. Smith, J. Am. Chem. Soc. 118, 5254 (1996).

37D. R. Yarkony, Rev. Mod. Phys. 68, 985 (1996).

38See supplementary material at http://dx.doi.Org/10.1063/1.3667203 for the parameters of the LVC Hamiltonian and derivation details.

39P. F. Barbara, T. J. Meyer, and M. A. Ratner, J. Phys. Chem. 100, 13148 (1996).

40In the case when pn(0) is Boltzmann with respect to the donor surface, it is easy to see that pn(t) = Pbz: (1) pn(0) becomes pbz after the small polaron transformation, and (2) it stays Boltzmann because it commutes with H0 and, thus, does not acquire time dependence.

41S. Mukamel, Principles of Nonlinear Optical Spectroscopies (Oxford University Press, New York, 1995).

42W. H. Louisell, Quantum Statistical Properties of Radiation (Wiley, Canada, 1973).

43Note that [A, B] is a constant when A and B are bosonic operators linear in atand a.

44G. F. Giuliani and G. Vignale, Quantum Theory of the Electron Liquid (Cambridge University Press, Cambridge, 2005).