Nonequilibrium Fermi's Golden Rule Charge Transfer Rates via the

Subscriber access provided by SUNY DOWNSTATE

Article

Nonequilibrium Fermi's Golden Rule Charge Transfer Rates via the Linearized Semiclassical Method Xiang Sun, and Eitan Geva J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00236 • Publication Date (Web): 29 Apr 2016 Downloaded from http://pubs.acs.org on April 30, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Nonequilibrium Fermi’s Golden Rule Charge Transfer Rates via the Linearized Semiclassical Method Xiang Sun and Eitan Geva∗ Department of Chemistry, University of Michigan, Ann Arbor, Michigan 48109-1055, United States E-mail: [email protected] (Edited April 27, 2016)

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract The nonequilibrium Fermi’s golden rule describes the transition between a photoexcited bright donor electronic state and a dark acceptor electronic state, when the nuclear degrees of freedom start out in a nonequilibrium state. In this paper, we derive a new expression for the nonequilibrium Fermi’s golden rule within the framework of the linearized semiclassical approximation. The new expression opens the door to applications of the nonequilibrium Fermi’s golden rule to complex condensed-phase molecular systems described in terms of anharmonic force fields. We show that the linearized semiclassical expression for the nonequilibrium Fermi’s golden rule yields the exact fully quantum-mechanical result for the canonical Marcus model, where the coupling between donor and acceptor is assumed constant (the Condon approximation) and the donor and acceptor potential energy surfaces are parabolic and identical except for a shift in the equilibrium energy and geometry. For this model, we also present a comprehensive comparison between the linearized semiclassical expression and a hierarchy of more approximate expressions, in both normal and inverted regions, and over a wide range of initial nonequilibrium states, temperatures and frictions.

I. Introduction Charge transfer (CT) processes in donor-acceptor molecular systems are among the most commonly encountered nonradiative electronic relaxation processes. 1–5 They can occur spontaneously, as in redox reactions, or they can be induced by electron injection, as in electrochemical devices, or by photoexcitation, as in photochemical and photovoltaic devices. 6–24 In the case of photoexcitation from the ground electronic state to an excited bright donor state, the initial state of the nuclear degrees of freedom (DOF) is a nonequilibrium one. For example, in the case of impulsive excitation, the initial state of the nuclear DOF would correspond to equilibrium on the ground state potential energy surface (PES), which is typically significantly different from that corresponding to equilibrium on the excited donor PES. Furthermore, in the case of nonimpulsive photoexcitation, one expects the initial state of the nuclear DOF to be affected 2 ACS Paragon Plus Environment

Page 2 of 55

Page 3 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


by the laser pulse properties such as bandwidth and chirp. 25 A convenient way of describing the rate of photoinduced CT that starts out with a nonequilibrium initial state of the nuclear DOF, is via the nonequilibrium Fermi’s golden rule (NE-FGR). 26–28 Similarly to the equilibrium FGR (E-FGR), 1,5,5,13,29–35 the validity of the NE-FGR relies on the assumption of weak electronic coupling between donor and acceptor states, in the sense that the second-order time-dependent perturbation theory is valid. However, unlike E-FGR, NE-FGR does not describe the CT dynamics in terms of a rate constant, but allows for arbitrary nonequilibrium initial state of the nuclear DOF. It should also be noted that the NE-FGR is not limited to CT processes, and can also be applied to other types of electronic transitions, such as excitation energy transfer processes. 36 Importantly, the NE-FGR is based on assuming that the dynamics of the nuclear DOF is governed by quantum mechanics. However, a fully quantum-mechanical calculation of the NEFGR is only possible when the donor and acceptor PESs are parabolic. For those cases, the fully quantum-mechanical expression for the NE-FGR can be obtained in closed form and calculated in a feasible manner. 5,14–16,26,28,29,37–40 It should also be noted that the case of shifted, but otherwise identical, parabolic donor and acceptor PESs and constant donor-acceptor coupling (the Condon approximation) corresponds to the canonical Marcus model for CT. 41–43 While the availability of fully quantum-mechanical closed-form expressions for the NE-FGR in the case of parabolic donor and acceptor PESs make those model systems into useful benchmarks, it is highly desirable to develop strategies for calculating NE-FGR under less restrictive conditions. For example, such a need can arise when dealing with CT that occurs in liquid solution or solid state photovoltaic devices, which correspond to large, complex and disordered molecular systems, which are most naturally described in terms of inherently anharmonic force fields. A common practice is to assume that the PESs in those cases are effectively parabolic, so that the environment can be characterized by the corresponding spectral density. However, modeling inherently anharmonic systems as harmonic is nontrivial and of questionable general validity. 44 It would therefore be desirable to develop an alternative strategy for calculating NE-FGR CT rates



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

from molecular dynamics (MD) simulations based on force fields of one’s choice, so as to facilitate a direct relation between CT rates and the underlying molecular structure and dynamics. In this paper, we propose such a general and direct strategy for calculating NE-FGR CT rates for complex systems described by anharmonic force fields. This new strategy is based on the linearized semiclassical (LSC) approximation. 34,34,35,45–54 It starts out by casting the fully quantum-mechanical expression for the NE-FGR in terms of a forward-backward real-time path integral. The linearization approximation consists of truncation after the first order of the power expansion of the forward-backward action in terms of the difference between the forward and backward paths. The linearization is justified by the fact that rapidly oscillatory contributions from forward-backward paths with large forward-backward actions would be diminished by destructive interference. As a result, one expects the forward-backward path integral to be dominated by contributions from forward-backward paths with small forward-backward actions, which arise when the difference between the forward and backward paths is small, and for which the linearization approximation would be valid. The ability of LSC to account for quantum effects in complex condensed-phase systems has been established by numerous previous studies. 34,45–54 LSC is also known to reproduce the exact quantum-mechanical E-FGR CT rate constant when the coupling between donor and acceptor is assumed constant (the Condon approximation) and the donor and acceptor potential energy surface are parabolic and identical except for a shift in the equilibrium energy and geometry. 35 Furthermore, it was also previously observed that it remains accurate when the PESs are allowed to be different from each other. 34,55 Those observations suggest that the LSC method should also be able to capture at least some nuclear quantum effects within NE-FGR, which are known to be particularly important in the inverted regime. The remainder of this paper is organized as follows. The theory is outlined in Sec. II. The model system is described in Sec. III. Results are presented and discussed in Sec. IV. Finally, concluding remarks are provided in Sec. V. Derivations of several key results are detailed in Appendices A and B.


Page 4 of 55

Page 5 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


II. Theory A.

The Nonequilibrium Fermi’s Golden Rule

Consider a system with the overall Hamiltonian Hˆ = Hˆ 0 + Hˆ I , where the zeroth-order Hamiltonian, Hˆ 0 , and the perturbation, Hˆ I , are given by Hˆ 0 =Hˆ D |DihD| + Hˆ A|AihA|,

(1)

Hˆ I =Γˆ DA |DihA| + Γˆ AD |AihD|.

(2)

Here, D and A denote the diabatic donor and acceptor electronic states, respectively, and Hˆ D/A are the corresponding nuclear Hamiltonians Pˆ 2 ˆ , +VD/A R Hˆ D/A = 2

(3)

ˆ = Rˆ 1 , ..., Rˆ N and Pˆ = Pˆ1 , ..., PˆN are the mass-weighted nuclear coordinates and where R

momenta. The donor and acceptor states are coupled through the perturbation, Hˆ I , where we invoke the Condon approximation by assuming Γˆ DA → ΓDA = Γ∗AD = Const. The system is assumed to start out at the donor electronic state with the state of the nuclear DOF described by the density operator ρˆ D (0):

ρˆ (0) = |Diρˆ D (0)hD|.

(4)

The donor occupancy at time t > 0 is then given by h i ˆ ˆ PD (t) = Tr[ρˆ (t)|DihD|] = Tr e−iHt/¯h ρˆ (0)eiHt/¯h |DihD| ,

(5)

where Tr(· · · ) = TrN Tre (· · · ) denotes trace over both nuclear (N) and electronic (e) Hilbert spaces. NE-FGR is based on treating the electronic coupling term, Hˆ I , as a small perturbation, 5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 55

within the framework of second-order perturbation theory. It leads to the following approximate expression for PD (t), if back reaction is negligible: 26

PD (t) = 1 −

Z t 0

Zt ′ ′ k(t ) + · · · ≈ exp − dt k(t ) , ′

(6)

0

where, 2 Z t′ 1 k(t ) = dτ C(t ′, τ ) . 2Re h¯ 0 ′

(7)

Here, h i ˆ ˆ C(t ′, τ ) → |ΓDA |2 TrN ρˆ D (t ′ )e−iHA τ /¯h eiHD τ /¯h ,

(8)

where ˆ

′

ˆ

′

ρˆ D (t ′) = e−iHD t /¯h ρˆ D (0)eiHDt /¯h .

(9)

The E-FGR expression for PD (t), according to which

PD (t) = exp (−kD→At) , ˆ

eq

(10) ˆ

can be obtained from Eq. (7) when ρˆ D (0) → ρˆ D = e−β HD /TrN [e−β HD ] (β = 1/kB T ), and provided that t > τc , where τc is the lifetime of the following equilibrium correlation function, i h ˆ ˆ Ceq (τ ) = |ΓDA |2 TrN ρˆ Deq e−iHA τ /¯h eiHD τ /¯h .

(11)

Under those conditions, the E-FGR CT rate constant is give by: 35 2 Z ∞ 1 kD→A = 2Re dτ Ceq (τ ) . h¯ 0


(12)

Page 7 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


B. The Linearized Semiclassical Approximation for The Nonequilibrium Fermi’s Golden Rule The NE-FGR, Eq. (7), is based on treating the nuclear DOF fully quantum-mechanically. As such, its use is limited to the very few model systems where exact fully quantum-mechanical real time dynamics is accessible. Thus, from the point of view of computational feasibility and tractability, it is highly desirable to formulate an approximate NE-FGR which would be based on treating the nuclear DOF in a classical-like manner. However, doing so in a self-consistent manner is nontrivial since the τ -dependence of C(t ′, τ ), Eq. (8), lacks a well-defined classical limit. 35 The LSC approximation 34,35 provides a rigorous approach for obtaining a uniquely defined classicallike expression for the NE-FGR. The derivation of the LSC approximation for C(t ′, τ ) is detailed in Appendix A. In what follows, we restrict ourselves to outlining the main steps and underlying assumptions. We start out by writing C(t ′, τ ) in path-integral form (in what follows we adopt a onedimensional notation for simplicity): m N+n Z dx+ dx− dxt+′ dxt−′ dxt ′ −τ hx+ |ρˆ D (0)|x− C(t , τ ) = |ΓDA | 0 0 0 0i 2π h¯ ε Z i + + − − − dx+ 1 · · · dxN−1 dx1 · · · dxN−1 exp h SD − SD ¯ Z i ˜+ ˜− + + − − S − SD . dx˜1 · · · dxñ−1 dx˜1 · · · dxñ−1 exp h¯ A ′

2

(13)

Here, the forward and backward actions for the time periods 0 → t ′ , and t ′ → t ′ − τ (see Fig. 1) are



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 55

given by

+ SD =



N−1 j=0

− SD =

m

∑ ε2 

N−1

m

∑ ε2

j=0



n−1

m SÃ+ = ∑ ε  2 j=0  n−1 m − S˜D = ∑ ε 2 j=0

x+j+1 − x+j

!2

x−j+1 − x−j

!2

ε

ε

x˜+j+1 − x˜+j

!2

x˜−j+1 − x˜−j

!2

ε

ε



−VD (x+j ) ,

(14)



−VD (x−j ) ,

(15)



−VA (x˜+j ) ,

(16)



−VD (x˜−j ) ,

(17)

± ± ± with {0, ε , 2ε , . . ., N ε = t ′}, {0, ε , 2ε , . . ., nε = τ } (n ≤ N), x± 0 = x0 , xN = xt ′ , x˜0 = xt ′ −τ and ∓ x˜± n = xt ′ (the limit of ε → 0 will be applied at a later stage).

The linearization approximation relies on the assumption that the forward-backward path integral, Eq. (13), is dominated by forward-backward paths for which the forward and backward paths are close to each other. The linearization is implemented by truncating the power expansion of the forward-backward action in terms of the difference between the forward and backward paths at first order. To this end, we change integration variables in Eq. (13) from {x+j , x−j } to {y j , z j } and from {x˜+j , x˜−j } to {y˜ j , z˜ j }, such that yj = y˜ j =

x+j + x−j 2 x˜+j + x˜−j 2

;

z j = x+j − x−j ,

( j = 1, . . . , N − 1)

(18)

;

z˜ j = x˜+j − x˜−j ,

( j = 1, . . . , n − 1)

(19)

The linearized forward-backward path integrals are thus given by + − SD − SD ≈ε − SÃ+ − S˜D ≈ε

N−1 h

∑

j=1

i m ′ (2y − y − y ) −V (y ) j j+1 j−1 D j z j + pt ′ zt ′ − p0 z0 , ε2

n−1 nh

∑

j=1

i o m ′ (2 y ˜ − y ˜ − y ˜ ) −V ( y ˜ ) z ˜ +U ( y ˜ ) − p˜t ′ zt ′ . j j j j j+1 j−1 ε2 8 ACS Paragon Plus Environment

(20) (21)

Page 9 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Here, U (y) ˜ = VD (y) ˜ −VA (y) ˜ is the energy gap between the donor and acceptor PESs, and V (y) ˜ = [VD (y) ˜ +VA (y)]/2 ˜ is the average PES. Next, we change integration variables from {y j } to { f j } and from {y˜ j } to { f˜j }, such that m (2y j − y j+1 − y j−1 ) −VD′ (y j ), ( j = 1, . . ., N − 1), ε2 m ′ f˜j = 2 (2y˜ j − y˜ j+1 − y˜ j−1 ) −V (y˜ j ), ( j = 1, . . . , n − 1). ε

fj =

(22) (23)

The integrations over {z j }, ( j = 1, . . ., N − 1) and {˜z j }, ( j = 1, . . ., n − 1) can be performed explicitly:

i 2π h¯ dz j exp δ ( f j ), ε f jz j = h¯ ε Z 2π h¯ ˜ i ˜ ε f j z˜ j = d˜z j exp δ ( f j ). h¯ ε Z

(24) (25)

It should be noted that, Eqs. (24) and (25) imply that yt = yt (y0 , p0 ) in Eq. (13) follows a classicallike trajectory on the donor PES during 0 → t ′ , and the average PES during t ′ → t ′ − τ , since m (2y j − y j+1 − y j−1 ) −VD′ (y j ) = 0 2 ε m ′ f˜j = 2 (2y˜ j − y˜ j+1 − y˜ j−1 ) −V (y˜ j ) = 0 ε

fj =

d2 y(t) = −VD′ [y(t)], 2 dt d2 ε →0 ′ −−→ m 2 y(t) = −V [y(t)]. dt ε →0

−−→ m

(26) (27)

In addition, the following integration over zt ′ in Eq. (13) is given by Z

i dzt ′ exp (pt ′ − p˜t ′ )zt ′ h¯

= 2π h¯ δ (pt ′ − p˜t ′ ).

(28)

Substituting Eqs. (24), (25), (28), (20) and (21) into Eq. (13), we apply the following identities



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 55

to enable an initial value representation (see Appendix B in Ref. 35), 1 m N−1 ∂ y 1 ∂ p0 lim ∂ f =m ∂y ′ , ε →0 ε ε 2 t n−1 1 m ∂ y˜ = 1 ∂ p˜t ′ . lim ∂ f˜ m ∂ y˜ ′ ε →0 ε ε 2 t −τ

(29) (30)

It should be noted that the Jacobian |∂ y/∂ f | is the determinant of (N − 1) × (N − 1) matrix whose ˜ ∂ f˜| is the determinant of (n − 1) × (n − 1) matrix (i, j)th element is ∂ y j /∂ f j and the Jacobian |∂ y/ whose (i, j)th element is ∂ y˜ j /∂ f˜j . Finally, taking the limit of ε → 0, we arrive at the LSC expression for the electronic coupling correlation function, Z

CW-AV (t ′, τ ) ≈ |ΓDA |2 dy0 dp0 [ρˆ D (0)]W (y0 , p0 ) Z ′ i t −τ av × exp − dtU (yt ) . h¯ t ′

(31)

Here, [ρˆ D (0)]W (y0 , p0 ) is the Wigner transform of the operator ρˆ D (0), which is given by 1 [ρˆ D (0)]W (y, p) = 2π h¯

D i z zE dz exp pz y − Aˆ y + , h¯ 2 2

Z

(32)

and yt = yt (y0 , p0 ), pt = pt (y0 , p0 ) is obtained from forward in time classical dynamics on the donor PES during 0 → t ′ , followed by backward in time dynamics on the average PES during t ′ → t ′ − τ . We label this approximation “W-AV” (Wigner (W) initial sampling with dynamics on average (AV) PES). The LSC approximation in Eq. (31) can be straightforwardly generalized to the case of a multidimensional system (Nd denotes the dimensionality)

C

W-AV

′

2

(t , τ ) ≈ |ΓDA |

Z

dR0 dP0 [ρˆ D (0)]W (R0 , P0 ) Z ′ i t −τ av × exp − dtU (Rt ) . h¯ t ′


(33)

Page 11 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Here, [ρˆ D (0)]W (R, P) =

1 2π h¯

Nd Z

Z Z i dZ exp P·Z R − ρˆ D (0) R + h¯ 2 2

(34)

The new LSC-based expression for NE-FGR, Eq. (33), is the first main results of this paper. This equation gives rise to the following general prescription for calculating the NE-FGR CT rates: 1. Sample initial (R0 , P0 ) based on the Wigner phase-space density, [ρˆ D (0)]W (R0 , P0 ). 2. Propagate forward in time classically on the donor PES, VD (R), from time 0 to time t ′ , to obtain (RtD′ , PtD′ ). 3. Propagate backward in time classically on the average PES, V (R), from time t ′ to time t ′ − τ , av to obtain (Rav τ , Pτ ), and evaluate

R t ′ −τ t′

dtU (Rtav).

4. Repeat over a large number of trajectories until you obtain a converged average of i h R ′ exp − h¯i tt′ −τ dtU (Rtav) . It should be noted that despite being classical-like, LSC can capture some quantum effects by basing the initial sampling on the Wigner phase-space density, as opposed to the classical phasespace density, and by requiring that the dynamics during the time period t ′ → t ′ − τ is performed on the average PES. It should also be noted that the LSC prescription outlined above is based on a nonequilibrium MD simulation. One reason for this has to do with the fact that the initial phase-space density is different from the fully-classical equilibrium phase-space density. Another reason is that switching from the donor PES to the average PES at time t ′ , implies nonequilibrium dynamics during the time interval t ′ → t ′ − τ , even if the system reached equilibrium on the donor PES during the time interval 0 → t ′.

C.

Other Approximations

In this section, we consider a hierarchy of increasingly more approximate expressions for the NE-FGR rate constant. We focus on approximate expressions that would simplify the calculation and/or reduce the computational cost. 11 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 55

The first approximation to be considered is based on the fact that initial sampling based on the exact Wigner phase-space density is not possible for general anharmonic many-body systems. Approximate methods for estimating the Wigner phase-space density in such systems are available. 46,50,53,54 However, they too are computationally costly. Thus, the first approximation that we consider corresponds to replacing initial sampling based on the Wigner phase-space density, [ρˆ D (0)]W (R0 , P0 ) in Eq. (33), with the corresponding classical phase-space density, [ρD (0)]Cl (R0 , P0 ): CC-AV (t ′, τ ) =|ΓDA |2

Z

dR0 dP0 [ρD (0)]Cl (R0 , P0 ) Z ′ i t −τ av × exp − dtU (Rt ) . h¯ t ′

(35)

We label this approximation “C-AV” (classical (C) initial sampling with dynamics on the average (AV) PES). Calculating CC-AV (t ′ , τ ) still requires switching PESs from the donor PES during time interval 0 → t ′ to the average PES during time interval t ′ → t ′ − τ . Replacing the average PES during time interval t ′ → t ′ − τ with the donor PES would therefore simplify the calculation. Furthermore, comparison with the C-AV approximation would shed light on the importance of switching to the average PES during time interval t ′ → t ′ − τ . Thus, the second approximation that we consider involves dynamics on the donor PES during both time intervals:

C

C-D

′

(t , τ ) =|ΓDA |

2

Z

dR0 dP0 [ρD (0)]Cl (R0 , P0 ) Z ′ i t −τ D × exp − dtU (Rt ) . h¯ t ′

(36)

We label this approximation “C-D” (classical (C) sampling with dynamics on donor (D) PES). Another potentially useful approximation corresponds to the limit where C(t ′, τ ) vanishes with increasing τ on a time scale which is shorter than the time scale on which the nuclear DOF move. Importantly, being in this limit also implies that C(t ′, τ ) is independent of whether the average PES


Page 13 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


or donor PES are used during the time interval t ′ → t ′ − τ . The third and fourth approximations that we consider are based on this limit, and differ only with respect to whether the initial sampling is based on the Wigner or classical phase-space densities: CW-0 (t ′, τ ) =|ΓDA |2

Z

dR0 dP0 [ρˆ D (0)]W (R0 , P0 ) i × exp U (Rt ′ )τ . h¯

C

C-0

′

(t , τ ) =|ΓDA |

2

(37)

Z

dR0 dP0 [ρD (0)]Cl (R0 , P0 ) i × exp U (Rt ′ )τ . h¯

(38)

We label those approximations “W-0” and “C-0” (Wigner (W) or classical (C) initial sampling and independent (0) of choice of PES during t ′ → t ′ − τ .

III. Demonstrative Application To A Spin-Boson Model The spin-boson model is widely used for modeling CT. In accord with the physical picture underlying Marcus theory, it is commonly assumed that the donor and acceptor electronic states are coupled to a primary mode, which represents a collective solvation coordinate consisting of contributions from a large number of nuclear DOF. The PESs for the primary mode are assumed to be given by identical parabolas, except for a shift in equilibrium position and energy between donor and acceptor states. The primary mode is also assumed to be coupled to a thermal bath that consists of the remaining nuclear DOF, which are represented by independent secondary modes with parabolic PESs which are independent of the electronic sate. Importantly, the coupling between the primary and secondary modes is assumed bilinear. The second main result of this paper is that the LSC expression for NE-FGR coincides with the exact quantum-mechanical NE-FGR for the spin-boson model, in the case where the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 55

nonequilibrium initial state corresponds to a shifted donor equilibrium state and within the Condon approximation. The proof of this very appealing feature of LSC NE-FGR is detailed in Appendix B. The fact that LSC NE-FGR is exact for the above mentioned spin-boson model also implies that it cannot serve as a benchmark model for testing LSC. It can however serve as a benchmark for testing the accuracy of the hierarchy of approximations of the LSC approximation outlined in Sec. C, and thereby shed light on the role of different aspects of LSC. The specific model that we choose for this purpose assumes Ohmic spectral density with exponential cutoff for the bath of secondary modes. Such a model was originally proposed by Garg, Onuchic, and Ambegaokar (GOA), 56 and has been used extensively in recent years for benchmarking numerous approximate methods for describing quantum dynamics, including Fermi’s Golden Rule (FGR), 13,21,26,35,57 LSC, 35,58 mixed quantum-classical Liouville, 59–62 surface hopping, 63–65 path integral, 26,66,67 mean-field Ehrenfest, 64 quantum master equation, 67–73 noninteracting blip approximation (NIBA), 74–77 multiconfiguration time-dependent Hartree (MCTDH), 78 and hierarchical equation of motion (HEOM) 79–82 methods. The GOA model is based on the following Hamiltonian (operator designations with ˆ are dropped from this point on for the sake of simplicity): 56 Py2 1 h¯ ωDA σz + + Ω2 (y + y0 σz )2 HGOA =ΓDA σx + 2 2 " 2 # N−1 2 pα 1 2 cα y 2 +∑ + ωα xα + 2 . 2 2 ωα α =1

(39)

Here, σx = |DihA| + |AihD| and σz = |DihD| − |AihA|; ΓDA = ΓAD is the electronic coupling coefficient (constant within the Condon approximation); −¯hωDA is the donor-to-acceptor CT reaction free energy; y, Py and Ω are the primary mode mass-weighted coordinate, momentum and frequency, respectively; {xα , pα , ωα } are the mass-weighted coordinates, momenta and frequencies of the secondary bath modes, respectively; {cα } are the coupling coefficients between the secondary modes and the primary mode; 2y0 is the shift in equilibrium geometry of the primary


Page 15 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


mode between the donor and acceptor states, such that the reorganization energy is given by

Er = 2Ω2 y20 .

(40)

The GOA model also assumes that the secondary bath modes can be described by an Ohmic spectral density with exponential cutoff:

Jo (ω ) ≡

π N cα2 δ (ω − ωα ) = ηω e−ω /ωc . 2 α∑ ω =1 α

(41)

Here, η is the friction coefficient, which determines the coupling strength between the primary mode and secondary modes, and ωc is the cutoff frequency. The initial state is obtained by shifting the donor equilibrium state by s along the primary mode coordinate: " ( " 2 ##) N 2 Py2 1 2 p 1 y 1 c α α exp −β + Ω (y + y0 + s)2 + ∑ + ωα2 xα + 2 ρˆ D (0) = ZD 2 2 2 2 ωα α =1

(42)

where ZD is the donor partition function. The fact that the primary and secondary modes are harmonic and bilinearly coupled implies an equivalent form of the GOA Hamiltonian in terms of the corresponding normal modes. The latter are obtained by diagonalizing the Hessian matrix: 

      D=     

Ω2 +

c2 ∑α ωα2 α

c1

c1

c2

ω12 ω22

c2 .. . cN−1

···



cN−1       ,   ..  .   2 ωN−1

(43)

such that TT DT = diag(ω˜ 12 , . . ., ω˜ N2 ), with ω˜ j the normal mode frequencies. It should be noted that the normal modes are obtained via numerical diagonalization of the Hessian matrix rather than 15 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 55

by using the asymptotic effective spectral density, which is only known in closed form in the limit Ω ≪ ωc . 83 The donor and acceptor Hamiltonians in terms of normal modes are given by: N

HD = ∑

j=1 N

HA = ∑

j=1

Pj2 2 Pj2

+

1 2

1 + 2 2

N

∑ ω˜ 2j R2j + h¯ ωDA ,

j=1

(44)

N

∑ ω˜ 2j

j=1

eq

Rj −Rj

2

.

Here, {R j , Pj , ω˜ j } are the normal mode coordinates, momenta and frequencies, respectively, and eq

{R j } are the equilibrium donor-to-acceptor shift in equilibrium geometry along the normal mode coordinates:         

eq R1





1      − c1 eq  R2   ω12 T = 2y T  0  . ..  ..  .     c eq − ωN−1 RN 2

N−1

The initial state in terms of normal modes is given by: " ( N P2 1 1 j ρˆ D (0) = exp −β ∑ + ZD 2 j=1 2

         

(45)

N

∑ ω˜ 2j (R j + S j )2

j=1

,

#)

(46)

where, S j = TTj1 s .

(47)

The exact fully quantum-mechanical NE-FGR for the GOA model is known in closed form. As we show in Appendix B, the LSC NE-FGR coincides with the exact fully quantum-mechanical


Page 17 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


NE-FGR. Thus, CExact (t ′, τ ) = CW-AV (t ′ , τ ) = ( eq 2 " N ω ˜ (R j β h¯ ω˜ j j ) 2 coth (1 − cos(ω˜ j τ )) |ΓDA | exp iωDA τ − ∑ 2¯ h 2 j=1 ) # N iω ˜ j Req S j j sin(ω˜ j t ′ ) + sin(ω˜ j τ − ω˜ j t ′ ) + i sin(ω˜ j τ ) − ∑ h¯ j=1

(48)

The other approximations for this model can also be obtained in closed forms and are given by

N

2 ω˜ j (Req j )

2 C (t , τ ) = |ΓDA | exp iωDA τ − ∑ (1 − cos(ω˜ j τ )) + i sin(ω˜ j τ ) 2¯h β h¯ ω˜ j j=1 N iω ˜ j Req j Sj ′ ′ sin(ω˜ j t ) + sin(ω˜ j τ − ω˜ j t ) −∑ h¯ j=1 2 N ω ˜ j (Req 2 j ) 2 C-D ′ C (t , τ ) = |ΓDA | exp iωDA τ − ∑ (1 − cos(ω˜ j τ )) + iω˜ j τ 2¯h β h¯ ω˜ j j=1 N iω ˜ j Req j Sj ′ ′ sin(ω˜ j t ) + sin(ω˜ j τ − ω˜ j t ) −∑ h¯ j=1 2 N ω ˜ j (Req β h¯ ω˜ j ω˜ 2j τ 2 j ) 2 W-0 ′ coth + iω˜ j τ C (t , τ ) = |ΓDA | exp iωDA τ − ∑ 2¯h 2 2 j=1 N iω ˜ j Req j Sj ′ −∑ cos(ω˜ j t )ω˜ j τ h¯ j=1 2 N ω ˜ j (Req ω˜ j τ 2 j ) 2 M ′ C (t , τ ) = |ΓDA | exp iωDA τ − ∑ + i · ω˜ j τ 2¯h β h¯ j=1 N iω ˜ j Req j Sj ′ −∑ cos(ω˜ j t )ω˜ j τ h¯ j=1 C-AV

′

2

(49)

(50)

(51)

(52)

The results reported below were obtained using N = 200 normal modes with max frequency 15ωc . The primary mode parameters were set to Ω = 0.5ωc and y0 = 1.0 (Er = 0.5ωc ). The electronic coupling was set to ΓDA = 0.1.

The temperatures and frictions considered were

kB T /¯hωc = 0.2, 0.5, 1.0 and η = 0.5, 1.0, 5.0.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 55

IV. Results and Discussion Before presenting NE-FGR rates, it is useful to consider the E-FGR rate constants for the system under consideration. Fig. 2 shows the E-FGR rate constants as a function of ωDA , for η = 1.0 and kB T /¯hωc = 0.2, 0.5, 1.0. It should be noted that the E-FGR rate constant increases with increasing temperature at ωDA /ωc > 1. This trend can be traced back to the fact that CT is an activated process in this region. It should also be noted that the E-FGR rate constant decreases with increasing temperature at ωDA /ωc < 1. This trend can be traced back to the low barrier in this region, such that CT is dominated by the prefactor, which decreases with increasing temperature (for example, √ the prefactor in Marcus theory ∼ 1/ kB T ). Next, we consider the behavior of the NE-FGR explicitly time dependent rate coefficient, k(t ′) (see Eq. (7)). On the left panel of Fig. 3, we show the behavior of the exact k(t ′) for different values of ωDA , at kB T /¯hωc = 1.0 and η = 1.0, with the initial nonequilibrium state shifted by s = 1.0 (see Fig. 4). Following a transient nonequilibrium period, k(t ′) is seen to reach a plateau value which coincides with the corresponding E-FGR rate constant (see right panel in Fig. 3). The duration of the nonequilibrium transient time period is ∼ 10ωc−1 , and observed to be relatively insensitive to the value of ωDA . Further insight can be gained by focusing on the following two representative values of ωDA : (1) ωDA = 0, which corresponds to CT in the normal region, (2) ωDA = 2, which corresponds to CT in the inverted region. For each case we consider several different initial states corresponding to shifting the donor equilibrium state by different amounts, s, along the primary mode coordinate (see Fig. 4). Comparisons of the exact/LSC NE-FGR with the above mentioned hierarchy of approximations, in terms of both PD (t) and k(t ′ ), for different initial states, as well as different values of temperature and friction, are shown in Figs. 5–8 (normal region) and in Figs. 9–12 (inverted region). The following noteworthy observations emerge from a close inspection of Figs. 5–12: In all cases considered, the NE-FGR rate coefficients, k(t ′ ), reach a plateau following a


Page 19 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


nonequilibrium transient period, with the plateau value given by the corresponding E-FGR rate constant. This is the expected behavior in condensed phase systems when the nonequilibrium initial state equilibrates on the time scale of the population relaxation. This should be contrasted with molecules in the gas phase, where population relaxation dynamics can be dominated by nonequilibrium effects for its entire duration. 28 The deviation between E-FGR and NE-FGR population relaxation dynamics, PD (t), extend beyond the initial transient time period during which k(t ′ ) 6= kD→A and increases with decreasing temperature and friction, as well as with increasing the shift of the initial state from equilibrium. The deviations between the various approximations increase with decreasing temperature and friction. Among the approximations, W-0 and C-0 appear to deviate most significantly from the exact/LSC NE-FGR, which is indicative of an increased sensitivity to the τ -dependence of C(t ′, τ ). The C-D approximation is also observed to deviate significantly from the exact/LSC NE-FGR. This implies that when the NE-FGR is sensitive to the τ -dependence of C(t ′, τ ), it is important to account for the fact that the dynamics takes place on the average PES, as opposed to the donor PES. Finally, initial sampling based on the Wigner phase-space density as opposed to the classical one, was observed to be somewhat important, particularly for the lowest temperatures and frictions considered. During the initial transient time period, k(t ′) tends to be significantly larger than kD→A in the inverted region, whereas the opposite trend is observed in the normal region. As a result, the deviations in PD (t) between E-FGR and NE-FGR are significantly larger in the inverted region than in the normal region. This can be traced back to the fact that, in the inverted region, enhancing the nonequilibrium nature of the initial state by increasing |s|, also enhances its ability to penetrate regions of configuration space with significantly larger overlaps between donor and acceptor nuclear wave functions, during the time period 0 → t ′ . Enhancing the overlap between donor and acceptor nuclear wave functions increases the efficiency of nuclear tunneling, thereby giving rise to larger transient CT rates. It should be noted that this is the case even if the nonequilibrium state starts out centered at a region of configuration space that corresponds to lower overlap in



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 55

comparison to the equilibrium state (as in s = 1.0, 3.0, 5.0 in Fig. 4(b)). This is because the subsequent nonequilibrium oscillatory dynamics on the donor PES during the time period 0 → t ′ would take the system to regions of enhanced overlap after a relatively short time lag, as evidenced by the stepwise behavior of PD (t) within the NE-FGR. Furthermore, the larger |s|, the larger the overlaps accessible, and the larger the enhancement of k(t ′) over kD→A . In contrast, the main effect of the nonequilibrium oscillatory dynamics on the donor PES during the time period 0 → t ′ in the normal region, is to induce an enhancement of k(t ′) in the vicinity of the crossing point between the donor and acceptor PESs, and a reduction of k(t ′) away from it. While this too gives rise to a a stepwise behavior of PD (t), its physical origin is clearly rather different and classical in nature. Furthermore, because of the resonance nature of the enhancement at the crossing point, k(t ′) tends to be lower than kD→A except in the vicinity of the crossing point.

V. Concluding Remarks We proposed a new LSC-based strategy for calculating NE-FGR CT rates in complex condensed phase systems. The new strategy goes beyond our previously proposed LSC-based expression for the E-FGR rate constant, by being able to account for nonequilibrium initial states. We have also established that even though the new LSC NE-FGR strategy is based on classical trajectories, it is still able to reproduce the exact fully quantum-mechanical NE-FGR in the case where the coupling between donor and acceptor is assumed constant (the Condon approximation) and the donor and acceptor PESs are parabolic and identical, except for a shift in the equilibrium energy and geometry. The ability to reproduce the exact quantum-mechanical result for this important benchmark model, which is consistent with the physical picture underlying Marcus theory, suggests that the LSC NE-FGR should be able to account for quantum effects in complex condensed phase systems for which the fully quantum-mechanical result is not accessible. We emphasize that doing so does not require one to map the Hamiltonian of complex condensed phase systems onto spin-boson Hamiltonians. Rather, the new LSC NE-FGR strategy allows one to


Page 21 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


calculate CT rates directly from nonequilibrium classical MD trajectories as dictated by the force fields of one’s choice. We also presented a comprehensive comparison between LSC NE-FGR and a hierarchy of more approximate variants. The comparison was performed within the context of the GOA model, for which the LSC NE-FGR coincides with the exact fully quantum-mechanical NE-FGR. The comparison was also performed in both normal and inverted regions, and covered a wide range of nonequilibrium initial states, temperatures and frictions. Several interesting trends were observed regarding the robustness of various approximations as a function of friction and temperature in both the normal and inverted regions. The results demonstrate the importance of accounting for the dynamics of C(t ′, τ ) during the time interval t ′ → t ′ − τ , as well as the importance of doing so properly (for example, on the average PES as opposed to the donor PES). Extending the analysis to cases where the donor and acceptor PESs have different curvatures, 84 or are anharmonic, 18,85 would be highly desirable. Work on such extensions is underway in our group and will be reported separately in future publications. In conclusion, the new LSC NE-FGR strategy proposed herein opens the door to using the NE-FGR for calculating CT rates in complex condensed phase systems, such as photoinduced CT in molecular diads and triads in liquid solution 6,7,9,14,23 and CT in organic photovoltaic interfaces. 8,10,13,15,16 Work on such timely applications is underway and will be reported in future papers.

Appendix A Derivation of the linearized semiclassical approximation for the nonequilibrium Fermi’s golden rule charge transfer rate in the Condon case In this appendix, we detail the derivation of the LSC approximation of C(t ′, τ ) (Eq. 8). We start out by using the position representation closure relations to cast C(t ′, τ ) in the following form



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 55

(one-dimensional notation is used throughout for simplicity): i h i i ˆ ′ i ˆ ′ i ˆ ˆ C(t ′, τ ) =|ΓDA |2 TrN e h¯ HD τ e− h¯ HDt ρˆ D (0)e h¯ HD t e− h¯ HA τ =|ΓDA |

2

Z

ˆ

ˆ

′

− + − iHD τ /¯h + dx+ |xt ′ ihxt+′ |e−iHDt /¯h |x+ 0 dx0 dxt ′ dxt ′ dxt ′ −τ hxt ′ −τ |e 0i ˆ

′

ˆ

− iHDt /¯h − ˆ D (0)|x− |xt ′ ihxt−′ |e−iHA τ /¯h |xt ′ −τ i . hx+ 0 ihx0 |e 0 |ρ

(A.1)

Next, we express the forward and backward propagators in terms of discrete path integrals with N time slices for period t ′ and n time slices for period τ , such that t ′ = N ε , τ = nε (n ≤ N) (the limit ε → 0 will be applied at a later stage). + − − For the time period of 0 → t ′ , we set x+ N = xt ′ and xN = xt ′ and note that

m N Z i + 2 + + dx1 · · · dxN−1 exp S , h¯ D 2π h¯ ε i N Z 2 i − mi − − − iHˆ Dt ′ /¯h − dx1 · · · dxN−1 exp − SD . hx0 |e |xt ′ i = 2π h¯ ε h¯

ˆ ′ hxt+′ |e−iHDt /¯h |x+ 0i=

(A.2) (A.3)

− − + + − For the time period of t ′ − τ → t ′, we set x˜+ 0 = xt ′ −τ , xñ = xt ′ , x˜0 = xt ′ −τ , xñ = xt ′ and note

that m n Z i ˜+ 2 + + S , dx˜1 · · · dxñ−1 exp h¯ A 2π h¯ ε i n Z 2 mi i ˜− ˆ iHD τ /¯h + − − |xt ′ i = hxt ′ −τ |e dx˜1 · · · dxñ−1 exp − SD . 2π h¯ ε h¯

ˆ hxt−′ |e−iHA τ /¯h |xt ′ −τ i =


(A.4) (A.5)

Page 23 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


+ − − Here, the forward actions, SD and SÃ+ , and backward actions, SD and S˜D , are defined as

+ SD =

N−1



m

∑ ε2

j=0

x+j+1 − x+j

ε



−VD (x+j ) ,

 ! − − 2 x − x m j j+1 − SD = ∑ ε −VD (x−j ) , 2 ε j=0   ! + + 2 n−1 m x˜ j+1 − x˜ j SÃ+ = ∑ ε  −VA (x˜+j ) , 2 ε j=0   ! − − 2 n−1 x ˜ − x ˜ m j j+1 − S˜D = ∑ ε −VD (x˜−j ) . 2 ε j=0 N−1



!2

(A.6)

(A.7)

(A.8)

(A.9)

The exact path-integral expression for C(t ′, τ ) is therefore given by m N+n Z + ˆ − − + − dx+ C(t , τ ) =|ΓDA | D (0)|x0 i 0 dx0 dxt ′ dxt ′ dxt ′ −τ hx0 |ρ 2π h¯ ε Z i + − + + − − S − SD dx1 · · · dxN−1 dx1 · · · dxN−1 exp h¯ D Z i ˜+ ˜− + + − − S − SD . dx˜1 · · · dxñ−1 dx˜1 · · · dxñ−1 exp h¯ A ′

2

(A.10)

We next change integration variables from {x+j , x−j } to {y j , z j } for the time period of 0 → t ′ , and from {x˜+j , x˜−j } to {y˜ j , z˜ j } for the time period of t ′ − τ → t ′ , which are given by yj = y˜ j =

x+j + x−j 2 + x˜ j + x˜−j 2

;

z j = x+j − x−j ,

( j = 1, . . . , N − 1)

(A.11)

;

z˜ j = x˜+j − x˜−j ,

( j = 1, . . . , n − 1)

(A.12)

It should be noted that the Jacobians of these transformations are unity: N−1 ∂ (x+ , x− ) j j

n−1 ∂ (x˜+ , x˜− ) j j

∏ ∂ (y j , z j ) = ∏ ∂ (y˜ j , z˜ j ) = 1. j=1

j=1

(A.13)

The linearization approximation is introduced by expanding the difference between the forward



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 55

and backward actions to first order in {z j } and {˜z j }: + − SD − SD ≈ε − SÃ+ − S˜D ≈ε

N−1 h

∑

j=0

n−1 h

∑

j=0

i m ′ (y − y )(z − z ) −V (y )z j+1 j j+1 j D j j , ε2

(A.14)

i m ′ ( y ˜ − y ˜ )(˜ z − z ˜ ) +U ( y ˜ ) −V ( y ˜ )˜ z j j j j j , j+1 j+1 ε2

(A.15)

Here, U (y) and V (y) are the energy gap and average PES, respectively, which are defined by

U (y) = VD (y) −VA (y),

(A.16)

V (y) =

(A.17)

VD (y) +VA (y) . 2

Relabeling summation index for {y j , z j } (similarly for {y˜ j , z˜ j }) N−1

j+1→ j

∑ (y j+1 − y j )z j+1 −−−−→

j=0 N

N−1

∑ (y j − y j−1)z j = ∑ (y j − y j−1 )z j + (yN − yN−1 )zN ,

j=1

(A.18)

j=1

we can rewrite the forward-backward actions: N−1 h

i m ′ (2y − y − y ) −V (y ) ∑ ε 2 j j+1 j−1 D j z j + pt ′ zt ′ − p0 z0, j=1 i o n−1 nh m ′ − (2 y ˜ − y ˜ − y ˜ ) −V ( y ˜ ) z ˜ +U ( y ˜ ) − p˜t ′ zt ′ . SÃ+ − S˜D =ε ∑ j j+1 j−1 j j j ε2 j=1 + − SD − SD =ε

(A.19) (A.20)

In the above equalities, we used the fact that in the limit of ε → 0 (N → ∞, n → ∞), for {y j , z j },

ε

m (yN − yN−1 )zN −→ pN zN = pt ′ zt ′ ε2 m ε 2 (y1 − y0 )z0 −→ p0 z0 ε ′

ε V (y0 )z0 −→ 0


(A.21) (A.22) (A.23)

Page 25 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and for {y˜ j , z˜ j },

ε

m (yñ − yñ−1 )˜zn −→ pñ zñ = p˜t ′ z˜t ′ = − p˜t ′ zt ′ ε2 m ε 2 (y˜1 − y˜0 )˜z0 −→ p˜0 z˜0 = p˜t ′ −τ z˜t ′ −τ = 0 ε ′

(A.24) (A.25)

ε V (y˜0 )˜z0 −→ 0

(A.26)

ε U (y˜0 ) −→ 0

(A.27)

− Since the forward and backward paths coincide at xt ′ −τ = x˜+ 0 = x˜0 , we have y˜0 = xt ′ −τ , z˜t ′ −τ = − + − z˜0 = 0. In addition, notice that the fact that x˜+ n = xt ′ , and xñ = xt ′ implies z˜t ′ = −zt ′ .

Substituting Eq. (A.19) and Eq. (A.20) back into Eq. (A.10), we obtain m N+n C(t , τ ) ≈ |ΓDA | 2π h¯ ε D Z z0 E i z0 dy0 dz0 exp − p0 z0 y0 + ρˆ D (0) y0 − h¯ 2 2 Z i dy˜t ′ −τ dyt ′ dzt ′ exp (pt ′ − p˜t ′ )zt ′ h¯ #) ( " Z i N−1 h m i ε ∑ 2 (2y j − y j+1 − y j−1 ) −VD′ (y j ) z j dy1 · · · dyN−1 dz1 · · · dzN−1 exp h¯ j=1 ε " # Z Z i n−1 dy˜1 · · · dyñ−1 exp ε ∑ U (y˜ j ) d˜z1 · · · d˜zn−1 h¯ j=1 ( " #) i n−1 h i m ′ × exp ε ∑ 2 (2y˜ j − y˜ j+1 − y˜ j−1 ) −V (y˜ j ) z˜ j . (A.28) ε h¯ j=1 ′

2

The integrations over z j , ( j = 1, . . ., N −1) and z˜ j , ( j = 1, . . . n−1) can be performed explicitly: h h i i m i ′ dz j exp ε 2 (2y j − y j+1 − y j−1 ) −VD (y j ) z j ε h¯ h i m 2π h¯ ′ = δ 2 (2y j − y j+1 − y j−1 ) −VD (y j ) , ε ε h h Z i i i m ′ d˜z j exp ε 2 (2y˜ j − y˜ j+1 − y˜ j−1 ) −V (y˜ j ) z˜ j ε h¯ h i m 2π h¯ ′ δ 2 (2y˜ j − y˜ j+1 − y˜ j−1 ) −V (y˜ j ) . = ε ε Z


(A.29)

(A.30)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 55

The integrations over y j , ( j = 1, . . . , N −1) and y˜ j , ( j = 1, . . ., n −1) are facilitated by changing the integration variables from y1 , . . . , yN−1 to f1 , . . . , fN−1 and from y˜1 , . . . , yñ−1 to f˜1 , . . ., fñ−1 , which are defined as m (2y j − y j+1 − y j−1 ) −VD′ (y j ), ( j = 1, . . ., N − 1), ε2 m ′ f˜j = 2 (2y˜ j − y˜ j+1 − y˜ j−1 ) −V (y˜ j ), ( j = 1, . . . , n − 1). ε

fj =

(A.31) (A.32)

Taking advantage of the fact that the Jacobians of these transformation, |∂ y/∂ f | (the determinant of (N − 1) × (N − 1) matrix whose (i, j)th element is ∂ y j /∂ f j ) and |∂ y/ ˜ ∂ f˜| (the determinant of (n − 1) × (n − 1) matrix whose (i, j)th element is ∂ y˜ j /∂ f˜j ) satisfy the following identities: 35 1 m N−1 ∂ y 1 ∂ p0 lim ∂ f =m ∂y ′ , ε →0 ε ε 2 t n−1 1 m ∂ y˜ = 1 ∂ p˜t ′ , lim ∂ f˜ m ∂ y˜ ′ ε →0 ε ε 2 t −τ

(A.33) (A.34)

the integral over y1 , . . . , yN−1 , z1 , . . . , zN−1 and the integral over y˜1 , . . . , yñ−1 , z˜1 , . . ., zñ−1 in Eq. (A.28) become Z i N−1 h 2π h¯ N−1 m dy1 · · · dyN−1 ∏ δ 2 (2y j − y j+1 − y j−1 ) −VD′ (y j ) ε ε j=1 N−1 N−1 2 N−1 ∂ p0 ∂y π ε ε 2 2π h¯ h ¯ = = ∂y ′ ∂ f ε ε m m t { f j =0}


(A.35)

Page 27 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and

respectively.

"

# i n−1 dy˜1 · · · dyñ−1 exp ∑ εU (y˜ j ) h¯ j=1 i n−1 h m ′ × ∏ δ 2 (2y˜ j − y˜ j+1 − y˜ j−1 ) −V (y˜ j ) ε j=1 " # i n−1 2π h¯ n−1 ∂ y˜ exp = ∑ εU y˜ j ( f˜j = 0)|V ∂ f˜ ˜ h¯ j=1 ε { f j =0} # " n−1 n−1 ∂ p˜t ′ i 2π h¯ n−1 ε ε 2 ˜ = ∂ y˜ ′ exp h¯ ∑ ε U y˜ j ( f j = 0)|V , ε m m t −τ j=1

2π h¯ ε

n−1 Z

(A.36) (A.37)

It is important to point out that yt = yt (y0 , p0 ) follows a forward classical trajectory from time 0 to t ′ on the donor PES, VD , since f j = 0 implies my¨ +VD′ (y) = 0 (the classical equation of motion) ending up with phase-space point (yt ′ , pt ′ ). Furthermore, the integral over zt ′ in Eq. (A.28) gives rise to Z

i dzt ′ exp (pt ′ − p˜t ′ )zt ′ = 2π h¯ δ (pt ′ − p˜t ′ ) . h¯

(A.38)

This delta function δ (pt ′ − p˜t ′ ) implies that the dynamics backward in time over the time period t ′ → t ′ − τ on the average PES, V , should start out with the same momentum that the dynamics over the time period 0 → t ′ ended up with at time t ′. Alternatively, and equivalently, one can reverse the sign of momentum at time t ′ , pt ′ → −pt ′ , and propagate forward in time on the average PES from t ′ to time t ′ + τ . R t ′ −τ dτ U (yt |V ). Recalling the definition In the limit ε → 0 (N → ∞, n → ∞), ∑n−1 j=1 ε U y˜ j |V → t ′

of Wigner transform the density operator, 1 ρW (y, p) = 2π h¯

Z

D i z ˆ zE dz exp pz y − A y + , h¯ 2 2


(A.39)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 55

we arrive at the LSC approximation for C(t ′, τ ): C(t ′, τ ) ≈|ΓDA |2

Z

dy0 dp0 [ρˆ D (0)]W (y0 , p0 ) Z ′ i t −τ × exp − dtU (yt |V ) . h¯ t ′

(A.40)

Extending to the Nd -dimensional case, we obtain: ′

C(t , τ ) ≈|ΓDA |

2

Z

dR0 dP0 [ρˆ D (0)]W (R0 , P0 ) Z ′ i t −τ av × exp − dtU (Rt ) . h¯ t ′

(A.41)

Appendix B Proof that the LSC approximation for nonequilibrium FGR charge transfer rate in the Condon case coincides with the exact quantum-mechanical result for a model consisting of shifted parabolic PESs when the initial state is given by shifted donor equilibrium state. In this appendix, we provide the derivation of exact quantum-mechanical expression for C(t ′, τ ), and the LSC approximation for C(t ′, τ ), for the case where the donor and acceptor PES correspond to shifted parabolas and the coupling between them is within the Condon approximation. We also assume that the initial nonequilibrium state corresponds to a shifted donor equilibrium state. We show that for this model, the LSC C(t ′, τ ) coincides with the exact fully quantum-mechanical result. We start out by considering a model where the donor and acceptor Hamiltonians are given by


Page 29 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(one-dimensional notation is used throughout for simplicity): pˆ2 1 Hˆ D = + ω 2 xˆ2 + h¯ ωDA , 2 2 pˆ2 1 Hˆ A = + ω 2 (xˆ − d)2 . 2 2

(B.1) (B.2)

Within the Condon approximation, the NE-FGR rate is given by |ΓDA |2 k(t ) = 2 2Re h¯ ′

Z t′ 0

dτ eiωDA τ C1 (t ′ , τ ),

(B.3)

where h i i iH ˆ ′ t ′ − i Hˆ A τ Hˆ D′ τ − h¯i Hˆ D′ t ′ ˆ D h h h ¯ ¯ ¯ e ρD (0)e . e C1 (t , τ ) =TrN e ′

(B.4)

Here, pˆ2 1 2 2 Hˆ D′ = Hˆ D − h¯ ωDA = + ω xˆ . 2 2

(B.5)

Performing the trace within the position representation and inserting position representation closure relations, C1 (t ′, τ ), Eq. (B.4), can be put in the following form: ′

C1 (t , τ ) =

Z

i

ˆ

′

dx0 dx1 dx2 dx3 dx4 hx0 |ρˆ D (0)|x1ihx1 |e h¯ HD′ t |x2 i i

i

ˆ

i

ˆ

ˆ

′

hx2 |e− h¯ HA τ |x3 ihx3 |e h¯ HD′ τ |x4 ihx4 |e− h¯ HD′ t |x0 i.

(B.6)

The initial density operator, ρˆ D (0) is next assumed to be the equilibrium density operator of a system with the shifted donor Hamiltonian Hˆ g =

pˆ2 2

+ 21 ω 2 (xˆ + s)2 . The matrix element of the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 55

initial density matrix is given by 1/2 ω 2π h¯ sinh(β h¯ ω ) ω 2 2 ((x0 + s) + (x1 + s) ) cosh(β h¯ ω ) − 2(x0 + s)(x1 + s) , exp − 2¯h sinh(β h¯ ω )

1 hx0 |ρˆ D (0)|x1 i = ZD

(B.7)

where ZD is the partition function

ZD =

Z +∞ −∞

dxhx|ρˆ D (0)|xi =

1 2 tanh(β h¯ ω /2) sinh(β h¯ ω )

1/2

.

(B.8)

The four propagators are also known in closed form: 1/2 ω hx1 |e |x2 i = 2π i¯h sin(−ω t ′ ) 2 iω ′ 2 (B.9) exp (x + x2 ) cos(ω t ) − 2x1 x2 , 2¯h sin(−ω t ′ ) 1 1/2 ω −iHˆ A τ /¯h hx2 |e |x3 i = 2π i¯h sin(ωτ ) iω 2 2 ((x2 − d) + (x3 − d) ) cos(ωτ ) − 2(x2 − d)(x3 − d) , exp 2¯h sin(ωτ ) iHˆ D′ t ′ /¯h

(B.10)

1/2 ω hx3 |e |x4 i = 2π i¯h sin(−ωτ ) 2 iω 2 exp (x + x4 ) cos(ωτ ) − 2x3 x4 , 2¯h sin(−ωτ ) 3 1/2 ω −iHˆ D′ t ′ /¯h |x0 i = hx4 |e 2π i¯h sin(ω t ′ ) 2 iω ′ 2 exp (x + x4 ) cos(ω t ) − 2x0 x4 . 2¯h sin(ω t ′ ) 0 iHˆ D′ τ /¯h


(B.11)

(B.12)

Page 31 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Thus, C1 (t ′, τ ) can be written as a multidimensional Gaussian integral: ω 52 C1 (t , τ ) = 2π h¯ ω 52 = 2π h¯ ′

p

Z 2 tanh(β h¯ ω /2) 1 T T dx0 · · · dx4 exp − x Ax + b x + c sin(ω t ′ ) sin(ωτ ) 2 p 2 tanh(β h¯ ω /2) (2π )5/2 1 T −1 √ exp b A b+c , sin(ω t ′ ) sin(ωτ ) det A 2

(B.13) (B.14)

where the matrix A, vector b and constant c are defined as 

    ω A=  h¯     



coth(β hω )−i cot(ω t ′ )

−1 sinh(β h¯ ω )

0

−1 sinh(β h¯ ω )

coth(β h¯ ω )+i cot(ω t ′ )

−i sin(ω t ′ )

0

0

0

−i sin(ω t ′ )

i cot(ω t ′ )−i cot(ωτ )

i sin(ωτ )

0

0

0

i sin(ωτ )

0

−i sin(ωτ )

i sin(ω t ′ )

0

0

−i sin(ωτ )

i cot(ωτ )−i cot(ω t ′ )

¯

0

i sin(ω t ′ )

     ,     

(B.15) 

1−cosh(β hω ) s sinh(β h¯ ω¯ ) 1−cosh(β hω ) s sinh(β h¯ ω¯ )

    ω ωτ ) b=  id 1−cos( sin(ωτ ) h¯    1−cos(ωτ )  id sin(ωτ )  0



     ,     

ω 2 1 − cosh(β h¯ ω ) 2 1 − cos(ωτ ) . c= − id s sinh(β h¯ ω ) sin(ωτ ) h¯

(B.16)

(B.17)

Performing the Gaussian integral, we obtain the exact quantum-mechanical expression for



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 55

NE-FGR rate in the Condon case: (

ω d2 β h¯ ω C1 (t , τ ) = exp − coth (1 − cos(ωτ )) + i sin(ωτ ) 2¯h 2 ) iω ds − sin(ω t ′ ) + sin(ωτ − ω t ′ ) . h¯ ′

(B.18)

Next, we consider the LSC approximation for C1 (t ′, τ ), Eq. (A.40): ′

C1 (t , τ ) =

Z

Z ′ i t −τ dx0 dp0 [ρˆ D (0)]W (x0 , p0 ) exp − dtU (xt |V ) . h¯ t ′

(B.19)

Here, the Wigner phase-space density for initial sampling of the ground-state Hamiltonian Hˆ g is given by 2 p0 1 2 2 β h¯ ω β h¯ ω 1 2 exp − [ρˆ D (0)]W (x0 , p0 ) = tanh + ω (x0 + s) . (B.20) tanh π h¯ 2 2 2 2 h¯ ω With initial condition (x0 , p0 ), the system is first propagated on the donor PES, VD′ = 21 ω 2 x2 for time period 0 → t ′ . The position and momentum at time t ′ are given by: x(t ′ ) =x0 cos(ω t ′ ) +

p0 sin(ω t ′ ), ω

p(t ′ ) = − ω x0 sin(ω t ′) + p0 cos(ω t ′ ).

(B.21) (B.22)

Next, the system is propagated backward in time on the average PES, V from time t ′ to time t ′ − τ with the initial conditions given by (x(t ′ ), p(t ′)) as given by Eq. (B.22). This is equivalent to flipping the momentum ( p(t ˜ ′) = −p(t ′ ) and x(t ˜ ′) = x(t ′ )) and propagating forward in time from ˜ ′), p(t ˜ ′)). The final position is then given by time t ′ to time t ′ + τ , with initial condition (x(t d d p(t ˜ ′) ′ x(t ˜ + τ ) = + x(t cos(ωτ ) + ˜ )− sin(ωτ ). 2 2 ω ′


(B.23)

Page 33 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Recalling that the energy gap U (x) = VD′ −VA = d ω 2 x − 12 d 2 ω 2 , we obtain: d p(t ˜ ′) ′ x(t ˜ )− U (x(t ˜ + τ )) =d ω sin(ωτ ) . cos(ωτ ) + 2 ω ′

2

(B.24)

So the integral in Eq. (B.19) becomes Z t ′ −τ

ω d2 sin(ωτ ) + ω dx0 sin(ωτ − ω t ′ ) + sin(ω t ′ ) 2 t′ + d p0 − cos(ωτ − ω t ′ )) + cos(ω t ′ ) . −

dtU (xt |V ) = −

(B.25)

Substituting Eq. (B.25) and Eq. (B.20) into Eq. (B.19), we obtain the LSC approximation for C1 (t ′ , τ ): (

ω d2 β h¯ ω C1 (t , τ ) = exp − coth (1 − cos(ωτ )) + i sin(ωτ ) 2¯h 2 ) iω ds sin(ω t ′ ) + sin(ωτ − ω t ′ ) , − h¯ ′

(B.26)

which is identical to the exact quantum-mechanical result in Eq. (B.18). Extending to the Ndimensional case, we obtain: ′

C1 (t , τ ) = exp

(

N

−∑

2 ω˜ j (Req j )

j=1 N iω ˜ j Req j Sj

−∑

j=1

h¯

2¯h

"

β h¯ ω˜ j coth 2

(1 − cos(ω˜ j τ )) + i sin(ω˜ j τ )

sin(ω˜ j t ′ ) + sin(ω˜ j τ − ω˜ j t ′ )

)

.

# (B.27)

Acknowledgement This project was supported by the National Science Foundation through grant CHE-1464477.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 55

References (1) Barbara, P. F.; Meyer, T. J.; Ratner, M. A. Contemporary Issues in Electron Transfer Research. J. Phys. Chem. 1996, 100, 13148–13168. (2) Chandler, D. In Classical and Quantum Dynamics in Condensed Phase Simulations; Berne, B. J., Cicootti, G., Coker, D. F., Eds.; World Scientific: New Jersey, 1997; Chapter 2, pp 25–50. (3) Weiss, U. Quantum dissipative systems; World Scientific: London, 1993. (4) May, V.; Kühn, O. Charge and energy transfer dynamics in molecular systems; Willey-VCH: Berlin, 2000. (5) Nitzan, A. Chemical Dynamics in Condensed Phases; Oxford University Press: New York, 2006. (6) Liddell, P. A.; Kuciauskas, D.; Sumida, J. P.; Nash, B.; Nguyen, D.; Moore, A. L.; Moore, T. A.; Gust, D. Photoinduced charge separation and charge recombination to a triplet state in a carotene-porphyrin-fullerene triad. J. Am. Chem. Soc. 1997, 119, 1400–1405. (7) Liddell, P. A.; Kodis, G.; Moore, A. L.; Moore, T. A.; Gust, D. Photo switching of photoinduced electron transfer in a dithienylethene-porphyrin-fullerene triad molecule. J. Am. Chem. Soc. 2002, 124, 7668–7669. (8) Brédas, J.-L.; Beljonne, D.; Coropceanu, V.; Cornil, J. Charge-Transfer and Energy-Transfer Processes in π -Conjugated Oligomers and Polymers: A Molecular Picture. Chem. Rev. 2004, 104, 4971–5004. (9) Rizzi, A. C.; van Gastel, M.; Liddell, P. A.; Palacios, R. E.; Moore, G. F.; Kodis, G.; Moore, A. L.; Moore, T. A.; Gust, D.; Braslavsky, S. E. Entropic changes control the charge separation process in triads mimicking photosynthetic charge separation. J. Phys. Chem. A 2008, 112, 4215–4223. 34 ACS Paragon Plus Environment

Page 35 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(10) Tian, H.; Yu, Z.; Hagfeldt, A.; Kloo, L.; Sun, L. Organic Redox Couples and Organic Counter Electrode for Efficient Organic Dye-Sensitized Solar Cells. J. Am. Chem. Soc. 2011, 133, 9413–9422. (11) Mishra, A.; Fischer, M. K. R.; Bäuerle, P. Metal-Free Organic Dyes for Dye-Sensitized Solar Cells: From Structure: Property Relationships to Design Rules. Angew. Chem. Int. Ed. 2009, 48, 2474–2499. (12) Feldt, S. M.; Gibson, E. A.; Gabrielsson, E.; Sun, L.; Boschloo, G.; Hagfeldt, A. Design of Organic Dyes and Cobalt Polypyridine Redox Mediators for High-Efficiency Dye-Sensitized Solar Cells. J. Am. Chem. Soc. 2010, 132, 16714–16724. (13) Zhao, Y.; Liang, W. Charge transfer in organic molecules for solar cells: Theoretical perspective. Chem. Soc. Rev. 2012, 41, 1075–1087. (14) Lee, M. H.; Dunietz, B. D.; Geva, E. Calculation From First Principles of Intramolecular Golden-Rule Rate Constants for Photo-Induced Electron Transfer in Molecular DonorAcceptor Systems. J. Phys. Chem. C 2013, 117, 23391–23401. (15) Lee, M. H.; Dunietz, B. D.; Geva, E. Calculation from First Principles of Golden-Rule Rate Constants for Photo-Induced Subphthalocyanine/Fullerene Interfacial Charge Transfer and Recombination in Organic Photovoltaic Cells. J. Phys. Chem. C 2014, 118, 9780–9789. (16) Lee, M. H.; Dunietz, B. D.; Geva, E. Donor-to-Donor vs. Donor-to-Acceptor Interfacial Charge Transfer States in the Phthalocyanine-Fullerene Organic Photovoltaic System. J. Phys. Chem. Lett. 2014, 5, 3810–3816. (17) Baiz, C. R.; Kubarych, K. J. Ultrafast vibrational Stark-effect spectroscopy: exploring chargetransfer reactions by directly monitoring the solvation shell response. J. Am. Chem. Soc. 2010, 132, 12784–12785.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 55

(18) Strodel, B.; Stock, G. Quantum modeling of transient infrared spectra reflecting photoinduced electron- transfer dynamics. J. Chem. Phys. 2006, 124, 114105. (19) Bredenbeck, J.; Helbing, J.; Hamm, P. Labeling Vibrations by Light: Ultrafast Transient 2D-IR Spectroscopy Tracks Vibrational Modes during Photoinduced Charge Transfer. J. Am. Chem. Soc. 2004, 126, 990–1991. (20) Schmidtke, S. J.; Underwood, D. F.; Blank, D. A. Following the solvent directly during ultrafast excited state proton transfer. J. Am. Chem. Soc. 2004, 126, 8620–8621. (21) Xu, D.; Schulten, K. Coupling of protein motion to electron transfer in a photosynthetic reaction center: investigating the low temperature behavior in the framework of the spin– boson model. Chem. Phys. 1994, 182, 91–117. (22) Ishizaki, A.; Fleming, G. R. Quantum Coherence in Photosynthetic Light Harvesting. Annu. Rev. Condens. Matter Phys. 2012, 3, 333–361. (23) Manna, A. K.; Balamurugan, D.; Cheung, M. S.; Dunietz, B. D. Unraveling the Mechanism of Photoinduced Charge Transfer in Carotenoid–Porphyrin–C60 Molecular Triad. J. Phys. Chem. Lett. 2015, 6, 1231–1237. (24) Liang, K. K.; Lin, C.-K.; Chang, H.-C.; Villaeys, A. A.; Hayashi, M.; Lin, S. H. Calculation of the vibrationally non-relaxed photo-induced electron transfer rate constant in dye-sensitized solar cells. Phys. Chem. Chem. Phys. 2007, 9, 853–859. (25) McRobbie, P. L.; Geva, E. Coherent Control of Population Transfer via Linear Chirp in Liquid Solution: The Role of Motional Narrowing. J. Phys. Chem. A 2016, DOI: 10.1021/acs.jpca.5b09736. (26) Coalson, R. D.; Evans, D. G.; Nitzan, A. A nonequilibrium golden rule formula for electronic state populations in nonadiabatically coupled systems. J. Chem. Phys. 1994, 101, 436–448.


Page 37 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(27) Cho, M.; Silbey, R. J. Nonequilibrium photoinduced electron transfer. J. Chem. Phys. 1995, 103, 595–606. (28) Izmaylov, A. F.; Mndive-Tapia, D.; Bearpark, M. J.; Robb, M. A.; Tully, J. C.; Frisch, M. J. Nonequilibrium Fermi golden rule for electronic transitions through conical intersections. J. Chem. Phys. 2011, 135, 234106. (29) Kubo, R.; Toyozawa, Y. Application of the method of generating function to radiative and non-radiative transitions of a trapped electron in a crystal. Prog. Theoret. Phys. 1955, 13, 160–182. (30) Jortner, J. Temperature dependent activation energy for electron transfer between biological molecules. J. Chem. Phys. 1976, 64, 4860–4867. (31) Mikkelsen, K. V.; Ratner, M. A. Electron tunneling in solid-state electron-transfer reactions. Chem. Rev. 1987, 87, 113–153. (32) Newton, M. D.; Sutin, N. Electron Transfer Reactions in Condensed Phases. Ann. Rev. Phys. Chem. 1984, 35, 437–480. (33) Newton, M. D. Quantum chemical probes of electron-transfer kinetics: the nature of donoracceptor interactions. Chem. Rev. 1991, 91, 767–792. (34) Shi, Q.; Geva, E. Nonradiative electronic relaxation rate constants from approximations based on linearizing the path-integral forward-backward action. J. Phys. Chem. A 2004, 108, 6109– 6116. (35) Sun, X.; Geva, E. Equilibrium Fermi’s Golden Rule Charge Transfer Rate Constants in the Condensed Phase: The Linearized Semiclassical Method vs Classical Marcus Theory. J. Phys. Chem. A 2016, DOI: 10.1021/acs.jpca.5b08280. (36) Jang, S.; Jung, Y.; Silbey, R. J. Nonequilibrium generalization of Forster–Dexter theory for excitation energy transfer. Chem. Phys. 2002, 275, 319–332. 37 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 55

(37) Kestner, N. R.; Logan, J.; Jortner, J. Thermal electron transfer reactions in polar solvents. J. Phys. Chem. 1974, 78, 2148–2166. (38) Jortner, J.; Bixon, M. Intramolecular vibrational excitations accompanying solvent-controlled electron transfer reactions. J. Chem. Phys. 1988, 88, 167–170. (39) Leggett, A. J.; Chakravarty, S.; Dorsey, A. T.; Fisher, M. P. A.; Garg, A.; Zwerger, W. Dynamics of the dissipative two-state system. Rev. Mod. Phys. 1987, 59, 1–85. (40) Endicott, J. S.; Joubert-Doriol, L.; Izmaylov, A. F. A perturbative formalism for electronic transitions through conical intersections in a fully quadratic vibronic model. J. Chem. Phys. 2014, 141, 034104. (41) Marcus, R. A. On the Theory of Oxidation-Reduction Reactions Involving Electron Transfer. I. J. Chem. Phys. 1956, 24, 966–978. (42) Marcus, R. A. Electrostatic Free Energy and Other Properties of States Having Nonequilibrium Polarization. I. J. Chem. Phys. 1956, 24, 979–989. (43) Marcus, R. A. Electron transfer reactions in chemistry. Theory and experiment. Rev. Mod. Phys. 1993, 65, 599–610. (44) Gottwald, F.; Ivanov, S. D.; Kühn, O. Applicability of the Caldeira–Leggett Model to Vibrational Spectroscopy in Solution. J. Phys. Chem. Lett. 2015, 6, 2722–2727. (45) Shi, Q.; Geva, E. A semiclassical theory of vibrational energy relaxation in the condensed phase. J. Phys. Chem. A 2003, 107, 9059–9069. (46) Shi, Q.; Geva, E. Vibrational energy relaxation in liquid oxygen from a semiclassical molecular dynamics simulation. J. Phys. Chem. A 2003, 107, 9070–9078. (47) Shi, Q.; Geva, E. A semiclassical generalized quantum master equation for an arbitrary system-bath coupling. J. Chem. Phys. 2004, 120, 10647–10658. 38 ACS Paragon Plus Environment

Page 39 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(48) Shi, Q.; Geva, E. A derivation of the mixed quantum-classical Liouville equation from the influence functional formalism. J. Chem. Phys. 2004, 121, 3393–3404. (49) Shi, Q.; Geva, E. A comparison between different semiclassical approximations for optical response functions in nonpolar liquid solutions. J. Chem. Phys. 2005, 122, 064506. (50) Ka, B. J.; Shi, Q.; Geva, E. Vibrational energy relaxation rates via the linearized semiclassical approximations: Applications to neat diatomic liquids and atomic-diatomic liquid mixtures. J. Phys. Chem. A 2005, 109, 5527–5536. (51) Ka, B. J.; Geva, E. Vibrational energy relaxation of polyatomic molecules in liquid solution via the linearized semiclassical method. J. Phys. Chem. A 2006, 110, 9555–9567. (52) Shi, Q.; Geva, E. A comparison between different semiclassical approximations for optical response functions in nonpolar liquid solution II. The signature of excited state dynamics on two-dimensional spectra. J. Chem. Phys. 2008, 129, 124505. (53) Vazquez, F. X.; Navrotskaya, I.; Geva, E. Vibrational energy relaxation rate constants via the linearized semiclassical method without force derivatives. J. Phys. Chem. A 2010, 114, 5682–5688. (54) Vazquez, F. X.; Talapatra, S.; Geva, E. Vibrational Energy Relaxation in Liquid HCl and DCl via the Linearized Semiclassical Method: Electrostriction versus Quantum Delocalization. J. Phys. Chem. A 2011, 115, 9775–9781. (55) McRobbie, P. L.; Geva, E. A benchmark study of different methods for calculating one- and two-dimensional optical spectra. J. Phys. Chem. A 2009, 113, 10425–10434. (56) Garg, A.; Onuchic, J. N.; Ambegaokar, V. Effect of friction on electron transfer in biomolecules. J. Chem. Phys. 1985, 83, 4491–4503. (57) Ohta, Y.; Soudackov, A. V.; Hammes-Schiffer, S. Extended spin-boson model for nonadiabatic hydrogen tunneling in the condensed phase. J. Chem. Phys. 2006, 125, 144522. 39 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 55

(58) Huo, P.; Coker, D. F. Influence of environment induced correlated fluctuations in electronic coupling on coherent excitation energy transfer dynamics in model photosynthetic systems. J. Chem. Phys. 2012, 136, 115102. (59) Shiokawa, K.; Kapral, R. Emergence of quantum-classical dynamics in an open quantum environment. J. Chem. Phys. 2002, 117, 7852. (60) Rekik, N.; Hsieh, C.-Y.; Freedman, H.; Hanna, G. A mixed quantum-classical Liouville study of the population dynamics in a model photo-induced condensed phase electron transfer dynamics. J. Chem. Phys. 2013, 138, 144106. (61) Sergi, A.; MacKernan, D.; Ciccotti, G.; Kapral, R. Simulatimg quantum dynamics in classical environments. Theor. Chem. Acc. 2003, 110, 49–58. (62) Frantsuzov, P. A. Quantum consideration of electron transfer solvent control. J. Chem. Phys. 1999, 111, 2075–2085. (63) Jain, A.; Subotnik, J. E. Surface hopping, transition state theory, and decoherence. II. Thermal rate constants and detailed balance. J. Chem. Phys. 2015, 143, 134107. (64) Xie, W.; Bai, S.; Zhu, L.; Shi, Q. Calculation of Electron Transfer Rates Using Mixed Quantum Classical Approaches: Nonadiabatic Limit and Beyond. J. Phys. Chem. A 2013, 117, 6196–6204. (65) Landry, B. R.; Falk, M. J.; Subotnik, J. E. Communication: The correct interpretation of surface hopping trajectories: How to calculate electronic properties. J. Chem. Phys. 2013, 139, 211101. (66) Escher, J. C.; Ankerhold, J. Quantum dynamics of a two-level system in a structured environment: Numerical study beyond perturbation theory. Phys. Rev. A 2011, 83, 032122. (67) Evans, D. G.; Coalson, R. D. Incorporating backflow into a relaxation theory treatment of


Page 41 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the dynamics of nonequilibrium nonadiabatic transition processes. J. Chem. Phys. 1995, 102, 5658–5668. (68) Shi, Q.; Geva, E. A new approach to calculating the memory kernel of the generalized quantum master equation for an arbitrary system-bath coupling. J. Chem. Phys. 2003, 119, 12063–12076. (69) Zhang, M.-L.; Ka, B. J.; Geva, E. Nonequilibrium quantum dynamics in the condensed phase via the generalized quantum master equation. J. Chem. Phys. 2006, 125, 044106. (70) Gelin, M. F.; Egorova, D.; Domcke, W. Exact quantum master equation for a molecular aggregate coupled to a harmonic bath. Phys. Rev. E 2011, 84, 041139. (71) del Rey, M.; Chin, A. W.; Huelga, S. F.; Plenio, M. B. Exploiting Structured Environments for Efficient Energy Transfer: The Phonon Antenna Mechanism. J. Phys. Chem. Lett. 2013, 4, 903–907. (72) Kolli, A.; O’Reilly, E. J.; Scholes, G. D.; Olaya-Castro, A. The fundamental role of quantized vibrations in coherent light harvesting by cryptophyte algae. J. Chem. Phys. 2012, 137, 174109. (73) Kolli, A.; Nazir, A.; Olaya-Castro, A. Electronic excitation dynamics in multichromophoric systems described via a polaron-representation master equation. J. Chem. Phys. 2011, 135, 154112. (74) Vierheilig, C.; Bercioux, D.; Grifoni, M. Dynamics of a qubit coupled to a dissipative nonlinear quantum oscillator: An effective-bath approach. Phys. Rev. A 2011, 83, 012106. (75) Thoss, M.; Wang, H.; Miller, W. H. Generalized forward-backward initial value representation for the calculation of correlation functions in complex systems. J. Chem. Phys. 2001, 114, 9220–9235.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 42 of 55

(76) Goychuk, I. A.; Petrov, E. G.; May, V. Control of long-range electron transfer in dynamically disordered molecular systems by an external periodic field. J. Chem. Phys. 1997, 106, 4522– 4530. (77) Goychuk, I. A.; Petrov, E. G.; May, V. Bridge-assisted electron transfer driven by dichotomically fluctuating tunneling coupling. J. Chem. Phys. 1995, 103, 4937–4944. (78) Bonfanti, M.; Tantardini, G. F.; Hughes, K. H.; Martinazzo, R.; Burghardt, I. Compact MCTDH Wave Functions for High-Dimensional System-Bath Quantum Dynamics. J. Phys. Chem. A 2012, 116, 11406–11413. (79) Tanimura, Y. Reduced hierarchy equations of motion approach with Drude plus Brownian spectral distribution: Probing electron transfer processes by means of two-dimensional correlation spectroscopy. J. Chem. Phys. 2012, 137, 22A550. (80) Tanaka, M.; Tanimura, Y. Multistate electron transfer dynamics in the condensed phase: Exact calculations from the reduced hierarchy equations of motion approach. J. Chem. Phys. 2010, 132, 214502. (81) Tanaka, M.; Tanimura, Y. Quantum Dissipative Dynamics of Electron Transfer Reaction System: Nonperturbative Hierarchy Equations Approach. J. Phys. Soc. Jpn. 2009, 78, 073802. (82) Iles-Smith, J.; Lambert, N.; Nazir, A. Environmental dynamics, correlations, and the emergence of noncanonical equilibrium states in open quantum systems. Phys. Rev. A 2014, 90, 032114. (83) Sun, X.; Geva, E. Exact vs. Asymptotic Spectral Densities in the Garg-Onuchic-Ambegaokar Charge Transfer Model and its Effect on Fermi’s Golden Rule Rate Constants. J. Chem. Phys. 2016, 144, 044106.


Page 43 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(84) Freed, K. F. Influence of Frequency Shifts on Electron Transfer Processes. J. Phys. Chem. B 2003, 107, 10341–10343. (85) Yeganeh, S.; Ratner, M. A. Effects of anharmonicity on nonadiabatic electron transfer: A model. J. Chem. Phys. 2006, 124, 044108–044110.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 55

Figures and captions

x +0

x+t’ x t’-τ x-t’

x -0 0

t’-τ

t’ time

Figure 1. Illustration of the forward-backward path integrals underlying NE-FGR. The blue forward and backward segments correspond to the time intervals 0 → t ′ and t ′ → 0, respectively, and the red forward and backward segments correspond to the time interval t ′ − τ → t ′ and t ′ → t ′ − τ , respectively.


Page 45 of 55

-1

10

kBT=0.2 kBT=0.5

kD

A / ωc

-2

10

kBT=1

-3

10

-4

10

0

1

2

3

ωDA / ωc

4

Figure 2. Exact/LSC E-FGR rate constants, kD→A , as a functions of ωDA , at the indicated temperatures, as obtained for the GOA model (see text for model parameters).

0.03

ωDA=0.5

0.02

kD A / ω c

ωDA=0

k(tʹ) / ωc

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ωDA=1 ωDA=1.5

0.01 ωDA=2 ωDA=3

0

ωDA=4

0

5

10

15

20

0

1

ωc tʹ

2

ωDA / ωc

3

4

Figure 3. Exact/LSC NE-FGR rate coefficients, k(t ′), (left panel) and the corresponding E-FGR rate constants, kD→A , (right panel), at the indicated values of ωDA . The results were obtained for the GOA model, with s = 1.0, kB T = 1.0 and η = 1.0.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a)

Energy

5

D

Page 46 of 55

A

3 1

-1

Reaction Coordinate

5

D

(b)

Energy

3

-3

A

1

Reaction Coordinate Figure 4. The nonequilibrium initial states considered in this work in the normal (ωDA = 0.0, panel (a)) and inverted (ωDA = 2.0, panel (b)) regions. The donor and acceptor PESs are plotted along the GOA model’s primary mode coordinate with Ω = 0.5 and y0 = 1. The orange numbers (−3, −1, 1, 3, 5) correspond to the different values of s considered.


Page 47 of 55

ωDA = 0 1

η = 0.5

s=1

η=1

W-AV C-AV C-D W-0 C-0 eq

η=5

0.8

kBT ___ = 0.2

ℏωc

P(t) D

0.6 1

kBT ___ = 0.5

0.8

ℏωc

0.6 1

kBT ___ =1

0.8

ℏωc

0.6 0

10

0

10

ωc t

0

10

20

0.05

kBT ___ = 0.2

ℏωc

k(tʹ)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0 0.04

kBT ___ = 0.5 0.02

ℏωc

0.02

kBT ___ =1

0.01 0

ℏωc

10

0

10

0

10

20

ωc tʹ Figure 5. Comparison of NE-FGR donor populations, PD (t), (upper panel) and NE-FGR rate coefficients, k(t ′), (lower panel) as obtained via different methods and at different temperatures and frictions, for the GOA model with ωDA = 0 and s = 1.0. The dashed lines correspond to E-FGR.



ωDA = 0

η = 0.5

s = -1

η=1


η=5

1 0.8

kBT ___ = 0.2

ℏωc

P(t) D

0.6 1

kBT ___ = 0.5

0.8

ℏωc

0.6 1

kBT ___ =1

0.8

ℏωc

0.6 0

10

0

10

ωc t

0

10

20

0.05

kBT ___ = 0.2

ℏωc

0 0.04

k(tʹ)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 48 of 55

kBT ___ = 0.5

ℏωc

0.02 0.03

kBT ___ =1

0.02 0.01 0

ℏωc

10

0

10

0

10

20

ωc tʹ Figure 6. Comparison of NE-FGR donor populations, PD (t), (upper panel) and NE-FGR rate coefficients, k(t ′), (lower panel) as obtained via different methods and at different temperatures and frictions, for the GOA model with ωDA = 0 and s = −1.0. The dashed lines correspond to E-FGR.


Page 49 of 55

ωDA = 0

η = 0.5

s=3

η=1


η=5

1 0.8

kBT ___ = 0.2

ℏωc

P(t) D

0.6 1

kBT ___ = 0.5

0.8

ℏωc

0.6 1

kBT ___ =1

0.8

ℏωc

0.6 0

10

0

10

ωc t

0

10

20

0.05

kBT ___ = 0.2

ℏωc

k(tʹ)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0 0.04

kBT ___ = 0.5

ℏωc

0.02 0.03

kBT ___ =1

0.02 0.01 0

ℏωc

10

0

10

0

10

20

ωc tʹ Figure 7. Comparison of NE-FGR donor populations, PD (t), (upper panel) and NE-FGR rate coefficients, k(t ′), (lower panel) as obtained via different methods and at different temperatures and frictions, for the GOA model with ωDA = 0 and s = 3.0. The dashed lines correspond to E-FGR.



ωDA = 0 1

η = 0.5

s=5

η=1


η=5

kBT ___ = 0.2

0.8

ℏωc

P(t) D

0.6 1

kBT ___ = 0.5

0.8

ℏωc

0.6 1

kBT ___ =1

0.8

ℏωc

0.6 0

10

0

10

ωc t

0

10

20

0.04

kBT ___ = 0.2

ℏωc

0

k(tʹ)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 50 of 55

0.04

kBT ___ = 0.5

0.02

ℏωc

0

kBT ___ =1

0.02

ℏωc

0.01 0 0

10

0

10

0

10

20

ωc tʹ Figure 8. Comparison of NE-FGR donor populations, PD (t), (upper panel) and NE-FGR rate coefficients, k(t ′), (lower panel) as obtained via different methods and at different temperatures and frictions, for the GOA model with ωDA = 0 and s = 5.0. The dashed lines correspond to E-FGR. The dashed lines correspond to E-FGR.


Page 51 of 55

ωDA = 2 1

η = 0.5


s=1

η=1

η=5

kBT ___ = 0.2

ℏωc

P(t) D

0.95 1

kBT ___ = 0.5

0.95

ℏωc

1

kBT ___ =1

0.95

ℏωc

0.9 0

10

0

0.02

10

ωc t

0

10

20

0.01

kBT ___ = 0.2

ℏωc

0 0.02

k(tʹ)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.01

kBT ___ = 0.5

ℏωc

0

0.01

kBT ___ =1

ℏωc

0.005 0 0

10

0

10

0

10

20

ωc tʹ Figure 9. Comparison of NE-FGR donor populations, PD (t), (upper panel) and NE-FGR rate coefficients, k(t ′), (lower panel) as obtained via different methods and at different temperatures and frictions, for the GOA model with ωDA = 2.0 and s = 1.0. The dashed lines correspond to E-FGR.



ωDA = 2 1

η = 0.5

s=3

η=1


η=5

kBT ___ = 0.2

ℏωc

0.95

P(t) D

1

kBT ___ = 0.5

0.95

ℏωc

1

kBT ___ =1

ℏωc

0.9 0

10

0

10

ωc t

0

10

20

0.02

kBT ___ = 0.2

0

k(tʹ)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 52 of 55

ℏωc

0.02

kBT ___ = 0.5

0

ℏωc

0.02

kBT ___ =1

0.01 0 0

ℏωc

10

0

10

0

10

20

ωc tʹ Figure 10. Comparison of NE-FGR donor populations, PD (t), (upper panel) and NE-FGR rate coefficients, k(t ′), (lower panel) as obtained via different methods and at different temperatures and frictions, for the GOA model with ωDA = 2.0 and s = 3.0. The dashed lines correspond to E-FGR.


Page 53 of 55

ωDA = 2 1

η = 0.5


s = -3

η=1

η=5

kBT ___ = 0.2

0.95

ℏωc

P(t) D

1

kBT ___ = 0.5

0.95

ℏωc

0.9 1

kBT ___ =1

0.95

ℏωc

0.9 0.85 0

10

0

10

0

10

20

ωc t 0.02

kBT ___ = 0.2

0

k(tʹ)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ℏωc

0.02

kBT ___ = 0.5

ℏωc

0 0.02

kBT ___ =1

0.01 0 0

ℏωc

10

0

10

0

10

20

ωc tʹ Figure 11. Comparison of NE-FGR donor populations, PD (t), (upper panel) and NE-FGR rate coefficients, k(t ′), (lower panel) as obtained via different methods and at different temperatures and frictions, for the GOA model with ωDA = 2.0 and s = −3.0. The dashed lines correspond to E-FGR.



ωDA = 2 1

η = 0.5

s=5

η=1


η=5

kBT ___ = 0.2

0.95

ℏωc

P(t) D

1 0.95

kBT ___ = 0.5

ℏωc

0.9 1

0.95

kBT ___ =1

ℏωc

0.9 0.85 0

10

0

10

ωc t

0

10

20

0.04 0.02

kBT ___ = 0.2

0

k(tʹ)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 54 of 55

ℏωc

0.02

kBT ___ = 0.5

ℏωc

0 0.02

kBT ___ =1

ℏωc

0 0

10

0

10

ωc tʹ

0

10

20

Figure 12. Comparison of NE-FGR donor populations, PD (t), (upper panel) and NE-FGR rate coefficients, k(t ′), (lower panel) as obtained via different methods and at different temperatures and frictions, for the GOA model with ωDA = 2.0 and s = 5.0. The dashed lines correspond to E-FGR.


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical TOC Entry D Energy

Page 55 of 55

x +0

x+t’

A

x t’-τ

g

x-t’

x -0 Reaction Coordinate

0

t’-τ


t’ time

Nonequilibrium Fermi's Golden Rule Charge Transfer Rates via the

Recommend Documents