Equilibrium Fermi's Golden Rule Charge Transfer Rate Constants in

Subscriber access provided by NEW YORK MED COLL

Article

Equilibrium Fermi’s Golden Rule Charge Transfer Rate Constants in the Condensed Phase: The Linearized Semiclassical Method vs. Classical Marcus Theory Xiang Sun, and Eitan Geva J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.5b08280 • Publication Date (Web): 09 Oct 2015 Downloaded from http://pubs.acs.org on October 13, 2015

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.

The Journal of Physical Chemistry A 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

The Journal of Physical Chemistry

Equilibrium Fermi’s Golden Rule Charge Transfer Rate Constants in the Condensed Phase: The Linearized Semiclassical Method vs. Classical Marcus Theory Xiang Sun and Eitan Geva∗ Department of Chemistry, University of Michigan, Ann Arbor, Michigan 48109-1055, USA E-mail: [email protected]

∗ To

whom correspondence should be addressed

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 In this paper, we present a comprehensive comparison between the linearized semiclassical expression for the equilibrium Fermi’s golden rule rate constant and the progression of more approximate expressions that lead to the classical Marcus expression. We do so within the context of the canonical Marcus model, where the donor and acceptor potential energy surface are parabolic and identical except for a shift in both the free energies and equilibrium geometries, and within the Condon region. The comparison is performed for two different spectral densities and over a wide range of frictions and temperatures, thereby providing a clear test for the validity, or lack thereof, of the more approximate expressions. We also comment on the computational cost and scaling associated with numerically calculating the linearized semiclassical expression for the rate constant and its dependence on the spectral density, temperature and friction.

I. Introduction Almost all chemical processes in molecular systems involve the redistribution of electronic charge. When the redistribution of electronic charge is significant, both in terms of the amount of charge transferred and the distance that is transferred across, one can identify different parts in the molecular system as donor and acceptor. Such charge transfer (CT) processes in donor-acceptor molecular systems are among the most studied chemical reactions. 1–23 They can occur spontaneously, as in redox reactions, 6,7,11,18,19 or can be induced by electron injection, as in electrochemical devices, or by photoexcitation, as in photochemical reactions and photovoltaic devices. 14,24–33 The rates of CT processes can vary on an extremely wide dynamical range and can be very sensitive to the underlying molecular structure, molecular environment, temperature, etc. Molecular level understanding of the mechanisms underlying CT processes and the factors that determine their rates is key for controlling and improving the efficiency of these reactions, and thereby advancing the technologies based on them. A fully quantum-mechanical expression for the CT rate constant can be obtained when the sys2 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


tem is assumed to start out with the nuclear degrees of freedom (DOF) at thermal equilibrium on the diabatic donor state potential energy surface (PES) and when the electronic coupling between the donor and acceptor electronic states can be assumed weak. The resulting equilibrium Fermi’s golden rule (FGR) expression for the CT rate constant puts it in terms of the Fourier transform of the fully quantum-mechanical real-time autocorrelation function of the electronic coupling at the frequency that corresponds to the free energy associated with the donor-to-acceptor CT reaction. 8,9,15,30,34–37 Importantly, the equilibrium FGR rate constant has the following appealing features. First, it is not restrictive with respect to the functional form of the diabatic donor and acceptor PESs. Second, it fully accounts for the quantum nature of the nuclear DOF. Third, it is explicitly dependent on the dynamics of the nuclear DOF via the electronic coupling real-time autocorrelation function. Fourth, it allows for non-Condon effects, which originate from possible explicit dependence of the electronic coupling on the nuclear DOF. The equilibrium FGR rate constant should be clearly distinguished from the commonly used classical Marcus expression for the CT rate constant. 2,38,39 Importantly, the equilibrium FGR rate constant is known to reduce to the classical Marcus rate constant only when the following conditions are met. First, the donor and acceptor diabatic PESs need to be described by identical parabolas, except for being shifted with respect to the corresponding equilibrium geometries and free energies. This condition implies that higher than second order terms in the cumulant expansion of the electronic coupling autocorrelation function can be assumed negligible. Second, the nuclear DOF need to be treated as classical. This condition requires that kB T ≫ h¯ ωmax , where

ωmax is the highest frequency of the Huang-Rhys-active modes. It should be noted that extensions that treat a few “inner-sphere” nuclear modes quantum-mechanically have also been proposed. 5 However, they require one to identify those modes ahead of time, and can only accommodate a few of them. Third, the electronic coupling autocorrelation function needs to decay on a time scale faster than that of the nuclear DOF. This transition-state-theory-like assumption implies that the CT rate constant only depends on the donor state equilibrium probability density of the nuclear DOF, and is therefore insensitive to nuclear dynamics. Fourth, the electronic coupling needs to be



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

assumed independent of the nuclear DOF (the Condon approximation). The popularity of the classical Marcus expression for the CT rate constant can be attributed to the fact that it is given by a simple analytical form in terms of three parameters—the electronic coupling coefficient, reorganization energy, and reaction free energy. As such, it offers a straightforward pathway for interpreting experimentally measured CT rate constants. However, the restrictive nature of the assumptions underlying the classical Marcus expression also makes it highly desirable to develop less restrictive methods for calculating the equilibrium FGR rate constant. 16 One example that demonstrates this need corresponds to photo-induced CT reactions, where the donor and acceptor states often correspond to a bright π − π ∗ state and a dark CT state, respectively, for which one expects significant deviations from the assumption of shifted but otherwise identical PESs. Another example corresponds to the inverted region, where the enhanced overlap between nuclear wavefunctions can make nuclear tunneling, as opposed to the classical activated process underlying classical Marcus theory, into the dominant CT pathway. 9,31–33,40 This scenario is particularly relevant to organic photovoltaic devices, where the rigid solid state environment and low dielectric constant of the organic host conspire to give rise to relatively small reorganization energies. 31–33 A third example corresponds to cases where the electronic coupling autocorrelation function decays on sufficiently slow time scales so as to make it sensitive to the underlying nuclear dynamics. 31 Lastly, the validity of the Condon approximation can become questionable, for example in cases where the transition between electronic states occurs through a conical intersection. 41,42 While the availability of fully quantum mechanical closed-form expressions for the equilibrium FGR CT rate constant when the diabatic PESs are harmonic and the electronic coupling terms are either constant, linear, or quadratic in the nuclear coordinates, makes these model systems into useful benchmarks, 34,41–43 it is highly desirable to develop strategies for calculating the equilibrium FGR CT rate constant under less restrictive conditions. For example, such a need can arise when dealing with CT that occurs in molecular liquids or solid state photovoltaic devices, which are most naturally described in terms of inherently anharmonic force fields. A common practice is to


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


map those systems onto one of the above mentioned models involving parabolic PESs. 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 that would allow for calculating FGR CT rate constants in systems described by anharmonic force fields in a direct manner and without first mapping them onto harmonic Hamiltonians. To this end, we propose to advance a general and direct strategy for calculating FGR CT rate constants for complex systems described by general anharmonic force fields, which would be based on the linearized semiclassical (LSC) method. Previous work provided numerous examples for the ability of the LSC method to account for a large number of quantum effects in complex condensedphase systems that, to the best of our knowledge, no other competing method can account for. 45–55 The LSC approximation for the FGR rate constant was originally derived in Ref. 48, where it was demonstrated on a few illustrative examples. Importantly, LSC is known to reproduce the exact quantum-mechanical equilibrium FGR CT rate constant for the canonical Marcus model, i.e. when the Condon approximation is valid and the donor and acceptor PESs are described by shifted parabolas. It should also be noted that we have previously observed that the LSC approximation remains accurate when the PESs are allowed to be different from each other. 48,56 In this paper, we present a comprehensive comparison between the LSC expression for the equilibrium FGR CT rate constant and a progression of more approximate expressions leading to the classical Marcus expression. The comparison is performed on a model system where the donor and acceptor PESs correspond to identical multi-dimensional parabolas that differ only with respect to their equilibrium geometries and energies. Results are presented for two different spectral densities and a wide range of temperatures and frictions. In addition, we also analyze the computational cost and scaling associated with numerically calculating the LSC approximation and its dependence on the spectral density, temperature, and friction. Importantly, the fact that the LSC approximation coincides with the fully quantum-mechanical result for this model system implies that it cannot be used to test the accuracy of the LSC approximation. However, rather than testing the accuracy of the LSC approximation, our goal 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 6 of 55

this paper is to provide a comprehensive comparison between the LSC expression and the series of approximations leading to the classical Marcus rate constant. The fact that for this model, the classical Marcus rate constant, Eq. (27), can be obtained under the aforementioned conditions makes this into a suitable model for this purpose. In addition, the fact that the exact fully quantummechanical equilibrium FGR rate constant is known in closed form also makes it a suitable model for examining the computational cost and scaling associated with numerically calculating the LSC approximation over a wide parameter range. The organization of this paper is as follows. The theory is outlined in section II. The model system is described in Sec. III. Results are presented and discussed in Sec. IV. Concluding remarks are given in Sec. V. Derivations of several key results are given in the appendices.

II. Theory A.

The equilibrium Fermi’s Golden Rule rate constant

We consider a two-state system with the following overall Hamiltonian:

Hˆ = Hˆ D |DihD| + Hˆ A|AihA| + Γˆ DA|DihA| + Γˆ AD|AihD| .

(1)

Here, D and A correspond to 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

(2)

ˆ = Rˆ 1 , ..., Rˆ N and Pˆ = Pˆ1 , ..., PˆN the mass-weighted nuclear coordinates and momenta. with R

Γˆ DA |DihA|+ Γˆ AD |AihD| is the electronic coupling, where Γˆ DA = Γˆ †AD are operators that only depend of the nuclear DOF and Γˆ DA → ΓDA = Const in the Condon approximation.


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


We assume that the system starts out in thermal equilibrium at the donor state:

ρˆ (0) = |Diρˆ Deq hD| ,

(3)

where (β = 1/kB T ),

ρˆ Deq =

exp(−β Hˆ D ) . Tr exp(−β Hˆ D )

(4)

It should be noted that photo-induced CT processes often start out with the system in a nonequilibrium state on the donor PES and that extensions of the theory to such cases are known. 41,57,58 However, here we choose to focus on situations where nuclear relaxation to equilibrium on the donor PES occurs on a time scale faster than that of CT, and for which the above assumption is justified. The main assumption underlying the equilibrium FGR rate constant for the transition from the donor state to the acceptor state, kD→A , is that the electronic coupling term, Hˆ I = Γˆ DA |DihA| + Γˆ AD |AihD|, constitutes a small perturbation compared to the uncoupled Hamiltonian, Hˆ 0 = Hˆ D |DihD|+ Hˆ A |AihA|. Thus, accounting for the effect of electronic coupling within the framework of secondorder perturbation theory combined with the complimentary assumption that the electronic cou−1 pling autocorrelation function, C(t) (see Eq. (6)) vanishes on a time scale much faster than kD→A ,

yields the following well known equilibrium FGR expression for kD→A :

kD→A

1 = 2 h¯

Z ∞

−∞

dtC(t) ,

(5)

where, i ˆ i eq HDt Γˆ DA exp − Hˆ At Γˆ AD C(t) = TrN ρˆ D exp h¯ h¯ i ˆ i ≡ exp HDt Γˆ DA exp − Hˆ At Γˆ AD . h¯ h¯ D

(6)

Here, TrN (· · · ) is the trace over the nuclear Hilbert space. Upon making the Condon approxima-



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

tion, Γˆ DA → ΓDA = Const, C(t) simplifies to give:

i ˆ i ˆ HDt exp − HAt C(t) → |ΓDA | TrN h¯ h¯ i ˆ i ≡ |ΓDA |2 exp HDt exp − Hˆ At . h¯ h¯ D 2

ρˆ Deq exp

(7)

B. The classical limit of the equilibrium Fermi’s Golden Rule rate constant— preliminary considerations The equilibrium FGR expression in Eq. (5) is given in terms of the fully quantum-mechanical real-time correlation function C(t), Eq. (6). The latter reflects the exact quantum dynamics of the nuclear DOF. However, from the point of view of computational feasibility and tractability, it would be desirable to treat the nuclear DOF as classical-like. Doing so in a self-consistent manner proves nontrivial, as demonstrated below. Naively, one may attempt to obtain the “classical limit” of the expression in Eq. (5) by replacing the quantum trace, TrN (· · ·), with a classical phase-space average,

R

dR0 dP0 (· · ·), the equi-

eq librium density operator, ρˆ D , with the corresponding equilibrium classical phase-space density,

ρDeq (R0 , P0 ), the quantum mechanical nuclear Hamiltonians, Hˆ D/A , by the corresponding classical Hamiltonians, HD/A (R, P), and the quantum-mechanical operators Γˆ AD/DA by the corresponding classical analogue, ΓAD/DA (R). Doing so would lead to the following approximation: 11

C(t) ≈

Z

eq dR0 dP0 ρD (R0 , P0 ) |ΓDA (R0 )|2 exp

i U (R0 )t h¯

,

(8)

where

U (R) = VD (R) −VA (R) .

(9)

is the vertical energy gap between donor and acceptor states at nuclear configuration R. However, the “classical limit” in Eq. (8) is not unique and thereby ill-defined. To demon8 ACS Paragon Plus Environment

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


strate this, we note that C(t) in Eq. (5) can be written in the following different looking and yet completely equivalent form: Zt i dτ Uˆ D (τ ) Γˆ DA (t) C(t) = TrN h¯ 0 Z i t ˆ exp+ − dτ UA (τ ) Γˆ AD . h¯ 0

ρˆ Deq exp−

(10)

Here, we have used the following identities: Zt i i ˆ i HDt = exp− Hˆ reft , dτ Uˆ D (τ ) exp exp h¯ h¯ 0 h¯ Z i ˆ i ˆ i t ˆ exp − HAt = exp − Hreft exp+ − dτ UA (τ ) , h¯ h¯ h¯ 0

(11)

where, exp± stands for positively-ordered (+) and negatively-ordered (−) exponentials, Uˆ D = ˆ = eiHˆ reft/¯h Ae ˆ −iHˆ reft/¯h , and Hˆ ref is an essentially arbitrary (except Hˆ D − Hˆ ref , Uˆ A = Hˆ A − Hˆ ref , A(t) for being time-independent and Hermitian) reference Hamiltonian. Taking the “classical limit” of the expression on the R.H.S. of Eq. (10), following the same procedure as that used to obtain Eq. (8), then leads to the following expression:

C(t) ≈

Z

eq dR0 dP0 ρD (R0 , P0 ) ΓAD (R0 ) ΓDA Rtref Zt i ref dτ U (Rτ ) . × exp h¯ 0

(12)

Here, Rtref is obtained by solving the classical equation of motion subject to the reference Hamiltonian, Hˆ ref , with the initial conditions (R0 , P0 ). Comparison of the two expressions for the “classical limit”, Eqs. (8) and (12), reveals that they are distinctly different, even though the corresponding fully quantum-mechanical expressions they originated from are guaranteed to coincide. More specifically, while Eq. (12) is affected by the dynamics of the nuclear DOF, Eq. (8) is not. Furthermore, the fact that the dynamics in Eq. (12) is dictated by an arbitrary reference Hamiltonian implies that it is not unique. For example,



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

if Href → HD , the dynamics captured by Eq. (12) would correspond to equilibrium fluctuations of eq

U (Rt ). This is because (R0 , P0 ) are sampled from ρD (R0 , P0 ), which is the classical equilibrium phase-space density that correspond to HD . However, if Href → HA , Eq. (12) would correspond to eq

eq

nonequilibrium dynamics associated with the relaxation process ρD (R0 , P0 ) → ρA (R∞ , P∞ ). The above mentioned difficulties associated with developing a classical limit of Eq. (5) can be traced back to the fact that the dynamics captured by Eq. (5) corresponds to the dephasing of electronic coherences, which in fact lacks a classical limit. To see this, it should be noted that since FGR is based on treating the electronic coupling perturbatively, the actual dynamics captured by Eq. (5) is dictated by the zeroth-order Hamiltonian, Hˆ 0 = Hˆ D |DihD| + Hˆ A|AihA|. Thus, given that the most general description of the initial state of the overall system is given by:

ρˆ (0) = |Diρˆ DD (0)hD| + |Aiρˆ AA(0)hA| + |Diρˆ DA(0)hA| + |Aiρˆ AD(0)hD| .

(13)

The density matrix of the overall system at time t would be given by: ˆ

ˆ

ρˆ (t) = e−iH0t/¯h ρˆ (0)eiH0t/¯h ˆ

ˆ

= e−iHDt/¯h |Diρˆ DD (0)hD|eiHDt/¯h ˆ

ˆ

+e−iHAt/¯h |Aiρˆ AA(0)hA|eiHAt/¯h ˆ

ˆ

ˆ

ˆ

+e−iHD t/¯h |Diρˆ DA (0)hA|eiHAt/¯h +e−iHAt/¯h |Aiρˆ AD (0)hD|eiHDt/¯h .

(14)

Hence, the electronic coherence at time t is given by: h i ˆ Dt/¯h −iHˆ A t/¯h i H e ρˆ AD (0) . hA|TrN [ρˆ (t)]|Di = TrN e

(15)

eq Thus, in the case where ρˆ AD (0) → ρˆ D , one can identify the dynamics captured by the integrand

of Eq. (6) with that of the dephasing of the electronic coherences. It should be noted that the substitution ρˆ AD (0) → ρˆ Deq should not be misinterpreted as implying that the system is in equilibrium. 10 ACS Paragon Plus Environment

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


eq What it implies instead is that the special case where ρˆ AD (0) → ρˆ D can be related to the coherence

dynamics captured by the FGR expression. The fact that the dynamics captured by C(t) in Eq. (6) corresponds to electronic dephasing is also the main reason for why predictions based on methods like fewest-switches surfacehopping, 20–22,59,60 Ehrenfest dynamics, 20 and ring-polymer molecular dynamics 19 can deviate significantly from the exact quantum-mechanical FGR rate constant (when it is known). More specifically, those discrepancies can be traced back to the inability of those methods to capture dynamical electronic coherences in regions of parameter space where C(t) decays on a time scale similar as that of the nuclear DOF, as captured by exp iHˆ Dt/¯h and exp −iHˆ At/¯h in Eq. (6).

In contrast, the classical Marcus expression for the rate constant relies on the assumption that C(t) decays on a time scale faster than that of the nuclear motions, and is therefore insensitive to how their dynamics is treated (see below). Thus, testing the ability of a method to capture the underlying nuclear dynamics cannot be based on using the classical Marcus expression as a benchmark and requires looking at regions of parameter space where the classical Marcus expression is not valid.

C.

The linearized semiclassical approximation

As opposed to the above mentioned ill-defined procedures for deriving a “classical limit” of what is essentially a non-classical quantity, the linearized semiclassical (LSC) approximation 48 provides a rigorous approach for obtaining a uniquely defined classical-like approximation for the fully quantum-mechanical C(t), Eq. (6). The deviation of the LSC approximation is detailed in Appendix A. It starts out by casting the fully quantum-mechanical expression for C(t), Eq. (6), in terms of a forward-backward real-time path integral. Here, the forward and backward paths originate from the forward and backward time evolution operators exp − h¯i Hˆ At and exp h¯i Hˆ Dt , respectively. The linearization approximation consists of truncating the power expansion of the forward-backward action in terms of the difference between the forward and backward paths after the first order. The linearization is justified by the expectation that rapidly oscillatory contributions 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

from forward-backward paths with large 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 resulting LSC approximation is given by (see Appendix A):

C

W-AV

Z 1 N eq dR0 dP0 Γˆ AD ρˆ D W (R0 , P0 ) (t) = 2π h¯ Zt i av av av × Γˆ DA W (Rt , Pt ) exp dτ U (Rτ ) , h¯ 0

(16)

where AW (R, P) =

Z

iZ·P/¯h

dZe

Z Z ˆ R − A R + 2 2

(17)

is the Wigner transform and (Rtav , Ptav ) are obtained via classical dynamics on the average PES, V (R) =

1 [VD (R) +VA (R)] , 2

(18)

starting with the initial conditions (R0 , P0 ). We label this approximation “W-AV” (Wigner (W) initial sampling with dynamics on average (AV) PES). Importantly, although Eq. (16) is similar to Eq. (12), there are two important differences between those two expressions: 1. In Eq. (16), the time evolution of (R, P) takes place on the rigorously derivable and uniquely defined average PES, V (R), as opposed to the arbitrary reference PES, Vref (R) in Eq. (12), which needs to be determined in an ad hoc manner. It should be noted that the dynamics does not take place on the donor or the acceptor PES. Also, since V (R) 6= VD (R), the dynamics captured by Eq. (16) corresponds to nonequilibrium relaxation rather than equilibrium fluctuations on the donor 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


eq 2. In Eq. (16), the Wigner transforms Γˆ AD ρˆ D W (R0 , P0 ) and Γˆ DA W (R0 , P0 ) replace the

corresponding classical expressions in Eq. (12). This implies non-classical initial sampling, thereby capturing the quantum behavior of the electronic coupling autocorrelation function on a short time scale. Importantly, one expects that capturing short time quantum effects in the time domain would allow for capturing high-frequency quantum effects in the frequency domain. In should be noted that high frequency in the context of CT reactions implies large reaction free energies, i.e. the inverted regime, where quantum effects are expected to be most pronounced (see below).

D. Intermediate approximations and the classical Marcus limit In this section, we consider the conditions under which the LSC expression for the equilibrium FGR rate constant reduces to the corresponding Marcus theory expression. To this end, we restrict ourselves to the case where the Condon approximation is valid, such that

C

W-AV

(t) =

|ΓDA |2 (2π h¯ )N

Z

eq dR0 dP0 ρD,W (R0 , P0 ) exp

Zt i av dτ U (Rτ ) . h¯ 0

(19)

Next, we define various intermediate approximate versions of the LSC expression whose accuracy we would later test on a benchmark model. The choice of approximations was motivated by the fact that using them would simplify and reduce the cost of the calculation. At the same time, one expects those approximations would also reduce the accuracy of the result. Our goal in this paper is to map the interplay between computational cost and accuracy over a wide range of parameter space. Initial sampling from the exact Wigner equilibrium phase-space density of a many-body anharmonic system is not possible in practice due to the sign problem associated with calculating the underlying multidimensional Fourier transform. Approximate methods for estimating the Wigner phase-space density in such systems have been proposed. 46,51,54,55 However, they typically require calculating the instantaneous normal modes 61 within an imaginary-time path-integral simulation



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

of the system. This requirement leads to a significant increase in computational cost. In comparison, initial sampling based on the corresponding classical equilibrium phase-space density is rather straightforward even for rather complex molecular systems. Thus, the first approximation eq

that we consider corresponds to replacing the Wigner phase-space density, ρD,W (R0 , P0 ), with the eq

corresponding classical phase-space density, ρD,Cl (R0 , P0 ). This yields:

C

C-AV

(t) =

|ΓDA |2 (2π h¯ )N

Z

eq dR0 dP0 ρD,Cl (R0 , P0 ) exp

Zt i av dτ U (Rτ ) . h¯ 0

(20)

We label this approximation “C-AV” (classical (C) sampling with dynamics on average (AV) PES). Even with classical initial sampling on the donor PES, having the dynamics take place on the average PES requires nonequilibrium MD simulations. Performing the dynamics on the donor PES would turn those into equilibrium MD simulations, thereby allowing one to take advantage of time reversal symmetry to enhance computational efficiency. Thus, the next approximation that we consider corresponds to performing both initial sampling and dynamics classically and on the donor PES: 59,60

C

C-D

(t) =

|ΓDA |2 (2π h¯ )N

Z

eq dR0 dP0 ρD,Cl (R0 , P0 ) exp

Zt i D U Rτ d τ . h¯ 0

(21)

We label the approximation “C-D” (classical (C) sampling with dynamics on donor (D) PES). It should be noted that one could also define a W-D approximation. However, we found that it did not provide new insights regarding the difference between Wigner and classical sampling, beyond those that could be obtained from comparing the W-AV and C-AV approximations. For this reason, the W-D approximation was not included in the comparison. The above mentioned approximations (W-AV, C-AV and C-D) call for simulating the dynamics of the nuclear DOF on the time scale defined by the lifetime of C(t). However, this requirement can be eliminated if this lifetime is shorter than the time scale of the nuclear DOF. Furthermore, in such a case the ambiguity regarding the PES on which the dynamics takes place is also eliminated,


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


since the FGR rate constant only depends on the donor equilibrium phase-space density:

C

W-0

(t) =

|ΓDA |2 (2π h¯ )N

Z

eq dR0 dP0 ρD,W (R0 , P0 ) exp

i U (R0 )t h¯

.

(22)

We label this approximation “W-0” (Wigner (W) sampling and independent (0) of choice of PES for dynamics). It should be noted that we choose to perform the initial sampling based on the Wigner phase space density in order to focus attention on the validity of assuming the abovementioned separation of time scales. Finally, the classical Marcus expression can be obtained from Eq. (22) by replacing the eq

Wigner phase-space density, ρD,W (R0 , P0 ), with the corresponding classical phase-space density, eq ρD,Cl (R0 , P0 ), and truncating the corresponding cumulant expansion after the second order. This

yields:

1 Cl 2 2 i − σD t C (t) = |ΓDA | exp hU iCl D h¯ 2¯h2 M

2

,

(23)

where

hU iCl D

=

1 (2π h¯ )

σDCl

N

Z

2

eq

dR0 dP0 ρD,Cl (R0 , P0 )U (R0 ) ,

2 Cl . = hU 2iCl − hU i D D

(24)

(25)

Substituting into Eq. (19) and performing the time integration in Eq. (5) explicitly, then yields:

M kD→A

|ΓDA |2 = h¯

s

2 # hU iCl D 2 exp − 2 . σDCl 2 σDCl "

2π

(26)

Further assuming that the donor and acceptor PESs are parabolic and identical except for a shift in equilibrium energy, ∆E (the donor-to-acceptor reaction free energy), and equilibrium geometry,



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

which gives rise to the reorganization energy Er , leads to the familiar classical Marcus expression: |ΓDA |2 M kD→A = h¯

r

# " π (∆E + Er )2 , exp − kB T Er 4kB T Er

(27)

where, 2 σDCl = −∆E − hU iCl Er = D . 2kB T

(28)

It should be noted that Eq. (28), and thereby Eq. (27), are only valid in the case where the donor and acceptor PESs can be given in terms of identical parabolas, whereas Eq. (26) does not rely explicitly on this assumption.

III. Model system The model system that we focus on in this paper consists of donor and acceptor PESs which are given in terms of identical multi-dimensional parabolas that differ only with respect to their equilibrium geometries and energies: N

Hˆ D = −∆E + ∑

j=1

N

Hˆ A =

∑

j=1

Pˆ j2 2

+

1 2

Pˆj2

1 + 2 2

N

∑ ω 2j Rˆ2j

,

j=1

N

eq 2 2 ˆ ω R − R ∑ j j j .

(29)

j=1

Here, {ω j | j = 1, ..., N} are the normal mode frequencies and ∆E is the donor-to-acceptor reaction free energy. In addition, we assume that Γˆ DA → ΓDA = Const. (the Condon approximation). It should be noted that for this model system the LSC equilibrium FGR rate constant coincides with the corresponding exact fully quantum-mechanical result, 48 which is known in closed form


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


(see Appendix C): ) N β ω h i ¯ j C(t) = CW-AV (t) = |ΓDA |2 exp − ∆Et − ∑ S j coth (1 − cos ω j t) + i sin ω j t h¯ 2 j=1 Z ∞ β h¯ ω 2 = |ΓDA | exp iωDA t − (1 − cos ω t) + i sin ω t . (30) d ω S(ω ) coth 2 0 (

Here, ωDA = −∆E/¯h, and S j =

eq 2 1 2¯h ω j (R j )

is the Huang-Rhys factor for the j-th mode, and S(ω )

is the Huang-Rhys density, which is given by:

S(ω ) ≡

1 2 ω j (Req j ) δ (ω − ω j ) . 2¯h ∑ j

(31)

It should be noted that the Huang-Rhys density is related to the spectral density, J(ω ), via the following identity:

S(ω ) =

4J(ω ) . πω 2

(32)

It should also be noted that C(t) = eiωDAt C1 (t), where C1 (t) is a smoothly decaying function (see Appendix C). Substituting back into Eq. (5), kD→A is then seen to be given by the Fourier transform, at frequency ωDA , of C1 (t). The intermediate approximations can also be calculated in closed form for this model and are given by:

C

C-AV

C

Z (t) = |ΓDA | exp iωDA t −

C-D

∞

2

2

(t) = |ΓDA | exp iωDAt −

0

2 (1 − cos ω t) + i sin ω t d ω S(ω ) β h¯ ω

Z ∞ 0

2 d ω S(ω ) (1 − cos ω t) + iω t β h¯ ω


,

,

(33)

(34)


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

W-0

2

(t) = |ΓDA | exp iωDAt −

Z ∞ 0

β h¯ ω d ω S(ω ) coth 2

Page 18 of 55

1 22 ω t + iω t 2

.

(35)

The fact that CW-0 (t) is a Gaussian function for this model, Eq. (35), makes it possible to perform the time integration analytically, so that the equilibrium FGR rate constant is given by the the Marcus-Levich formula: 16

W-0 kD→A

|ΓDA |2 = h¯

r

π (∆E + Er )2 exp − , a 4¯h2 a

(36)

where

a=

1R∞ 2 2 0 dω S(ω )ω coth

β h¯ ω 2

.

(37)

The classical Marcus rate constant can then be obtained from Eq. (36) by taking its high temperature limit, β h¯ ω ≪ 1, to obtain Eq. (27) with the reorganization energy given by:

Er = h¯

Z ∞

dω S(ω )ω .

(38)

0

IV. Results A.

Comparison of FGR rate constants

In this section, we compare the LSC rate constant to the aforementioned, more approximate, expressions that lead to the classical Marcus rate constant. We do so for two spectral densities. The first spectral density that we consider corresponds to the case where the nuclear modes form a narrow band: 43,48 (ω − ωop )2 1 π 3 exp − . J1 (ω ) = ηω 4 2σ 2 2¯h(2πσ 2 )1/2

(39)

Here, ωop is the central phonon frequency (such that 2π /ωop defines the typical nuclear time scale), σ = 0.1ωop is the phonon band width and η is the friction coefficient, which serves as a 18 ACS Paragon Plus Environment

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


measure of the strength of the coupling between the nuclear and electronic DOF. It should be noted that the LSC rate constant for this spectral density was reported in Ref. 48 for a specific choice of friction and temperature. The results reported below extend the analysis to a wide range of frictions and temperatures, and the comparison with those more approximate expressions for the FGR rate constant. The second spectral density that we consider corresponds to the case where the electronic DOF are coupled to a primary mode of frequency Ω, whose equilibrium geometry between the donor p and acceptor states is shifted by 2y0 = 2Er /Ω2, where Er is the reorganization energy. This primary mode is also bilinearly coupled to a bath of unshifted modes with the spectral density J0 (ω ) = ηω e−ω /ωc , where ωc ≫ Ω. It should be noted that 2π /Ω defines the typical nuclear time scale in this case. It should be noted that in this case, the friction coefficient, η , serves as a measure of the strength of the coupling between the primary mode and secondary modes. The fact that the primary and secondary modes are coupled implies that they are not the normal modes of this system. At the same time, the fact that they are bilinearly coupled implies that the Hamiltonian can be expressed in terms of normal modes. As is well known, the normal mode spectral density of this system can be obtained in closed form and is given by: 62

ηω 1 . J2 (ω ) = Ω2 Er 2 2 (Ω − ω 2 )2 + η 2 ω 2

(40)

It should also be noted that this choice of spectral density, with a relatively broad phonon band, is believed to be representative of CT reactions that take place in soft condensed phase environments, such as liquid solution and biopolymers. The FGR rate constants for J1 (ω ) are shown as a function of ωDA in Fig. 1, for different values −1 serve as energy and time units. Similar of temperature and friction. In this case, h¯ ωop and ωop

results for J2 (ω ) are shown in Figs. 2–6, for different values of Ω and Er . In this case, arbitrary energy/time units are used, since the results only depend on the relative values of Ω, Er , η and kB T . In all cases, we set ΓDA = 0.1. It should be noted that the condition kD→A τc ≪ 1, where



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

τc is the decay time of C(t), is satisfied for most the cases reported in Figs. 1–6. However, this condition, which is necessary in order to describe the dynamics by a rate constant, becomes harder to satisfy at the lowest temperatures, frictions and reaction free energies. In the very few instances where this is the case, lower values of ΓDA would be required in order for rate kinetics to be valid. However, since the dependence of the rate constant on ΓDA is through the multiplicative constant |ΓDA |2 , setting ΓDA to a lower values would not affect the general trends that we are concerned with here. Inspection of Figs. 1–6 reveals that generally speaking, the deviations between the exact/LSC expression and the various more approximate expressions decrease with increasing temperature and friction. This is consistent with the expectation that a classical description of the nuclear DOF becomes more valid and C(t) decays faster with increasing temperature and friction. The deviations between the exact/LSC expression and the approximate expressions are also seen to increase with increasing ωDA = −∆E/¯h. In this context, it is important to remember that the FGR rate constant is given in terms of the Fourier transform of C1 (t) at the frequency ωDA . Thus, one expects a classical treatment of the nuclear DOF to fail when the donor-to-acceptor reaction free energy becomes significantly larger than kB T (∼ 0.6 kcal/mol at room temperature). This is particularly relevant as one crosses over into the inverted region (|∆E| = h¯ ωDA > Er ), where the increased overlap between donor and acceptor vibrational wave functions is expected to enhance the efficiency of nuclear tunneling. In contrast, a purely classical treatment of the nuclear DOF treats the donor-to-acceptor reaction as an activated process. As a result, the LSC/exact results, which fully account for quantum nuclear effects such as tunneling, are seen to impose an upper bound on the rate constant. Accounting for quantum nuclear effects also leads to a much gentler decrease in the rate constant with increasing ωDA in comparison to the corresponding classical predictions. Furthermore, the temperature dependence of the exact/LSC rate constant is seen to be significantly weaker than that of the corresponding classical result, which corresponds to the Arrhenius dependence, typical of activated processes. Finally, the assumption that τc is much slower than the characteristic time


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


scale of nuclear dynamics is seen to break down rapidly with decreasing temperature, leading to predictions which are orders of magnitude smaller than the exact/LSC results. In comparison, the errors associated with treating the nuclear DOF classically grow more gradually with decreasing temperature. Thus, both W-0 and Marcus approximations, which bypass the need for accounting for the nuclear dynamics by assuming the above mentioned separation of time scales between τc and the time scale of nuclear dynamics, are in fact far more sensitive to decreasing temperatures and increasing frictions than approximations that account for it, even if they do so only approximately. The LSC approximation accounts for the quantum nature of the nuclear DOF in two ways: (1) Initial sampling of nuclear DOF based on the Wigner phase-space density, as opposed to the classical phase-space density; (2) Subsequent dynamics of nuclear DOF on the average PES, as opposed to the donor PES. While both seem to play a significant role at enhancing the rate constant, performing the dynamics on the average PES appears to be the more important of the two. In this context, it should be noted that while performing the dynamics on the average PES can be implemented in a relatively straightforward manner, initial sampling based on the Wigner phasespace density remains challenging in the case of a complex molecular system. 54,55 Next, we turn to the differences between the two spectral densities as reflected by the FGR rate constants that they give rise to. For J1 (ω ), the dependence of the exact/LSC rate constant on ωDA is oscillatory. The maxima, which correspond to multiples of ωop , represent 1-phonon, 2-phonon ... n-phonon resonances with the narrow phonon band. In contrast, J2 (ω ) is characterized by a broad phonon band, thereby giving rise to a smooth dependence of the rate constant on ωDA , which becomes monotonically decreasing in the inverted regime. The oscillatory character in the case of J1 (ω ) is seen to diminish with increasing temperature and friction. Interestingly, both the W-0 and Marcus approximations, which do not account for nuclear dynamics, are unable to capture the oscillatory behavior. In addition, the C-D approximation, which replaces the average PES with the donor PES is seen to shift the oscillatory pattern in a manner that depends on friction. The dependence of the results for J2 (ω ) on Ω and Er is also noteworthy. It should be noted



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

that the reorganization energy for J2 (ω ) is given by Er = 2Ω2 y20 , where 2y0 is the equilibrium geometry shift of the primary mode. Therefore, the reorganization energy being independent of friction implies that the classical Marcus rate constant is independent of η in this case. Thus, Figs. 2–4 show the effect of increasing Ω while fixing y0 fixed, while Figs. 2, 5, 6 show the effect of increasing y0 while keeping Ω fixed. Increasing the reorganization energy by either increasing Ω or y0 is seen to decrease the deviations between the exact/LSC and the other approximations. This can be attributed to the fact that quantum effects become prominent in the inverted region (¯hωDA > Er ). Thus, increasing Er extends the range of validity of approximations that rely on classical treatment of the nuclear DOF to higher ωDA . Finally, for the sake of comparison, we also show the results for J2 (ω ), for a choice of parameters used as a benchmark in several recent papers (see Fig. 7). 18,20,21,23 It should be noted that in Refs. 18,20,21,23 the parameters are given in terms of atomic units by Ω = 3.5 × 10−4 (77cm−1 ), Er = 2.39 × 10−2(5245cm−1 ), η = 1.2 × 10−3, kB T = 9.5 × 10−4(300K), and ΓDA = 5 × 10−5, which is equivalent to the choice Ω = 1, Er = 68.3, y0 = 5.84, η = 3.43, kB T = 2.71, ΓDA = 0.143 used here. As can be seen, the general trends observed in Figs. 2–6 are also observed for this choice of parameters.

B. Computational cost and scaling of LSC In this section, we consider the computational cost and scaling associated with calculating the FGR rate constant within LSC. To this end, we describe the nuclear DOF in terms of N normal modes whose frequencies are given by ω1 , . . . , ωN , such that ω j = j∆ω ( j = 1, 2, ...N). The equilibrium displacement of the j-th mode is then given by:

eq Rj

=

s

2¯hS(ω j )∆ω = ωj

s

8¯hJ(ω j )∆ω . πω 3j

(41)

For J1 (ω ), we used N = 250 and ∆ω = 0.002ωop . For J2 (ω ) we used N = 1000 and ∆ω = 0.015.


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


The numerical calculation of CW-AV (t) was carried out via a MC simulation,

C

W-AV

Zt 1 NMC sets i av U (Rτ ) dτ , (t) = exp h¯ 0 NMC (R ,P∑)∈ρ 0

0

(42)

W

where NMC is the number of trajectories averaged over. Initial sampling of (R0 , P0 ) was based on the phase space probability density N

ρW (R0 , P0 ) = ∏ ρ j (R j0 , Pj0) ,

(43)

j=1

where " β h¯ ω j β h¯ ω j 2 1 exp − tanh ρ j (R j0 , Pj0) = tanh π h¯ 2 2 h¯ ω j # 2 Pj0 ω 2j R2j0 × . + 2 2

(44)

Trajectories were propagated on the average PES, V (R), via the velocity-Verlet algorithm, with time step δ t = 0.002 for J1 (ω ) and δ t = 0.0005 for J2 (ω ). The numerically calculated LSC FGR rate constants are compared with the exact analytical result in Figs. 8 and 9 for J1 (ω ) and J2 (ω ), respectively. For J1 (ω ), the numerical results (Fig. 8) obtained with NMC = 5 × 106 are seen to be in excellent agreement with the exact analytical result. The agreement is still reasonable with NMC = 5 × 104 , but deteriorates somewhat in comparison to NMC = 5 × 106 , especially in cases where the FGR rate constant is small (corresponding to minima in the oscillatory resonance structure). It should be noted that the smaller rate constants will always be harder to converge since the same absolute error associated with numerically calculating the quantity of interest will give rise to a larger relative error when the value gets smaller. This trend can also be observed in the case of J2 (ω ) in Fig. 9, where the numerical result with NMC = 5 × 104 becomes unreliable due to noise at ωDA ∼ 10. Finally, we show the relative error averaged over the ωDA range considered here, (ωcut = 6ωop



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

for J1 (ω ) and ωcut = 10 for J2 (ω )),

δk =

1 #points

exact W-AV ωDA =ωcut k numerical (ωDA ) − kanalytical (ωDA )

∑

W-AV (ω ) kanalytical DA

ωDA =0

(45)

as function of NMC , for J1 (ω ), at different temperatures and frictions (see Fig. 10). The slower convergence with decreasing friction and temperature can be traced back to the, on average, smaller value of the corresponding rate constants. A similar trend is demonstrated in Fig. 11 for the case of J2 (ω ) with Ω = 1, Er = 2, kB T = 1 and η = 1.

V. Concluding Remarks Classical Marcus theory provides a very effective approach for fitting experiments and predicting general trends. However, obtaining a detailed molecular level understanding of the underlying dynamics requires one to move toward an approach that can start out with a molecular model, where the interactions would be described in terms of general anharmonic force fields. The approach should also account for non-Condon effects and the dynamics of the nuclear DOF, including quantum nuclear effects. The comprehensive comparison between the LSC method and approximations leading to the classical Marcus theory presented in this paper highlights the limitations of the two major assumptions underlying classical Marcus theory, namely treating the nuclear DOF as classical and assuming that their dynamics is slow enough so as not to matter. Both assumptions break down when |∆E| ≫ kB T . Furthermore, the predictions based on those assumptions were found to deteriorate rapidly once they cease to be valid, leading to rate constants that can be orders of magnitude smaller than the exact values. Overcoming those limitations call for a method that can account for the quantum nature of the nuclear DOF and provide a self-consistent description of their dynamics, at least on a short time scale. The LSC method does this, while still maintaining a semiclassical description of nuclear


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


DOF. Its ability to overcome those limitations is clearly demonstrated by the rather remarkable observation that its predictions for the canonical Marcus model (shifted but otherwise identical diabatic PESs within the Condon approximation) coincide with the exact fully quantum-mechanical result. However, the practical usefulness of the LSC method lies in the fact that it is rather straightforward to apply to complex molecular systems described by anharmonic force fields of one’s choice. It should be noted that the LSC method has been shown to be able to accurately capture a wide range of quantum effects in complex condensed-phase systems described by anharmonic force fields. 45–55 In the case of the equilibrium FGR expression for the CT rate constant considered here, this ability can be traced back to the use of Wignerized phase space densities and performing the dynamics on the average PES. These qualities turn the LSC method into an attractive platform for calculating CT dynamics in a manner that captures quantum effects, while not being restricted by assumptions like parabolic PESs, Condon approximation and neglecting dynamical nuclear effects. This is particularly important in the commonly encountered inverted regime, where the enhanced overlap between nuclear wave functions on the donor and acceptor PESs can turn nuclear tunneling, as opposed to the activated process underlying classical Marcus theory, into the dominant mechanism for CT. 9 Another advantage of the LSC method is that it can be extended to cases where the system starts out in a nonequilibrium state on the donor PES. 41,57,58 The work presented here can be extended in many directions, including testing the accuracy of LSC-based approaches in systems that involve a nonequilibrium initial state, non-Condon electronic coupling and anharmonic PESs, as well as going beyond the perturbative limit implied by FGR. Work on such extensions is underway and will be reported in future publications.

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 26 of 55

Appendix A Derivation of the linearized semiclassical approximation for the equilibrium Fermi’s Golden Rule rate constant In this appendix, we provide an outline of the derivation of the LSC approximation of C(t) (Eq. (6)). To this end, we use the position representation closure relations to put C(t) in the following form (one-dimensional notation is used for simplicity): i h eq iHˆ Dt/¯h ˆ −iHˆ A t/¯h ˆ C(t) =Trn e ΓDA e ΓAD ρˆ D =

Z

(A.1)

ˆ

− + − − iHDt/¯h − dx+ |xN i 0 dx0 dxN dxN hx0 |e ˆ

eq

+ + −iHAt/¯h + ˆ ˆ ˆ |x− i . hx− |x0 ihx+ N |ΓDA |xN ihxN |e D 0 0 |ΓAD ρ

(A.2)

Next, we express the forward and backward propagators in terms of discrete path integrals with N time slices separated by the time interval ε = t/N (the limit of N → ∞ will be applied at a later stage): −iHˆ At/¯h + hx+ |x0 i = N |e

Z

+ dx+ N−1 · · · dx1

N−1

ˆ

∏ hx+j+1|e−iHAε /¯h |x+j i j=0

m N/2 Z i + dx+ · · · dx+ e h¯ SA , = N−1 1 2π h¯ ε i

iHˆ Dt/¯h − hx− |xN i = 0 |e

=

Z

− dx− 1 · · · dxN−1

N−1

(A.3)

ˆ

∏ hx−j |eiHDε /¯h |x−j+1i j=0

mi 2π h¯ ε

N/2 Z

i −

− − h¯ SD . dx− 1 · · · dxN−1 e


(A.4)

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


− Here, the forward action, SA+, and backward action, SD , are defined as

SA+ =



N−1

m

∑ ε2

j=0



N−1

m − SD = ∑ ε 2 j=0

x+j+1 − x+j

ε x−j − x−j+1

ε

!2 !2



−VA (x+j ) ,

(A.5)



−VD (x−j ) .

(A.6)

The exact path-integral expression of C(t) is thus given by m N Z − − dx+ · · · dx+ C(t) = N dx0 · · · dxN 0 2π h¯ ε + − i eq hx+ |Γˆ ρˆ |x− ihx− |Γˆ |x+ ie h¯ (SA −SD ) . AD D

0

DA

N

0

N

(A.7)

We next change integration variables from {x+j , x−j } to {y j , z j }, which are given by

yj =

x+j + x−j 2

;

z j = x+j − x−j ,

( j = 0, . . ., N) .

(A.8)

It should be noted that the Jacobian of this transformation is unity: N−1 ∂ (x+ , x− ) j j

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

.

(A.9)

We then introduce the linearization approximation by expanding the difference between the forward and backward actions to first order in {z j }: − SA+ − SD ≈ε

N−1 h

∑

j=0

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

(A.10)

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

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

(A.11)

VD (y) +VA (y) . 2

(A.12)

V (y) =



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

Relabeling summation index, 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 =

j=1

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

,

(A.13)

j=1

we can then rewrite the forward-backward action: − SA+ − SD

N−1 h

i m ′ (y ) (2y − y − y ) −V j zj ∑ ε 2 j j+1 j−1 j=1 +U (y j ) + pN zN − p0 z0 .

=ε

(A.14)

To derive the last equality, we used the fact that in the limit of ε → 0 (N → ∞):

ε

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

(A.15) (A.16)

ε V (y0 )z0 −→ 0

(A.17)

ε U (y0) −→ 0

(A.18)

Substituting Eq. (A.14) back into Eq. (A.7), we obtain m N Z i C(t) ≈ dy0 dyN dz0 dzN exp (pN zN − p0 z0 ) h¯ 2π h¯ ε D E D zN E z0 zN z0 eq y0 + Γˆ AD ρˆ D y0 − yN − Γˆ DA yN + 2 2 " 2 #2Z Z N−1 i dy1 · · · dyN−1 exp ∑ εU (y j ) dz1 · · · dzN−1 h¯ j=1 ( ) i i N−1 h m ′ × exp ∑ ε ε 2 (2y j − y j+1 − y j−1) −V (y j ) z j . h¯ j=1


(A.19)

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


The integrations over z j , ( j = 1, . . ., N − 1) can be performed explicitly: Z ∞

i i hm ′ ε 2 (2y j − y j+1 − y j−1 ) −V (y j ) z j dz j exp h¯ ε −∞ h i 2π h¯ m ′ δ 2 (2y j − y j+1 − y j−1 ) −V (y j ) . = ε ε

(A.20)

The integration over y j , ( j = 1, . . . , N − 1) is facilitated by changing the integration variables from y1 , . . ., yN−1 to f1 , . . . , fN−1 , which are defined by

fj =

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

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

(A.21)

Taking advantage of the fact that the Jacobian of this transformation, |∂ y/∂ f | (the determinant of the (N − 1) × (N − 1) matrix whose (i, j)th element is ∂ y j /∂ f j ) satisfies the following identity (see Appendix B for proof): 1 m N−1 ∂ y 1 ∂ p0 lim ∂ f = m ∂ yt , N→∞ ε ε 2

(A.22)

the integral over y1 , . . . , yN−1 , z1 , . . ., zN−1 in Eq. (A.19) becomes

# i N−1 dy1 · · · dyN−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 " # N−1 N−1 i 2π h¯ ∂y exp = ∑ εU y j ( f j = 0)|V ε ∂ f { f j =0} h¯ j=1 " # N−1 N−1 ∂ p0 2π h¯ N−1 ε ε 2 i = . ∂ yN exp h¯ ∑ ε U y j |V ε m m j=1

2π h¯ ε

"

N−1 Z

(A.23)

It is important to point out that yt = yt (y0 , p0 ) follows a classical trajectory on the average PES, V , ′

since f j = 0 implies my¨ +V (y) = 0 (the classical equation of motion). Rt Now, in the limit of ε → 0 (N → ∞), ∑N−1 j=1 ε U y j |V → 0 dτ U (τ ). Recalling the definition of



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

ˆ Wigner transform of an operator A,

AW (y, p) =

Z

D z zE dzeizp/¯h y − Aˆ y + , 2 2

(A.24)

then lead to the LSC approximation for C(t):

C(t) ≈

Z 1 eq dy0 dp0 Γˆ AD ρˆ D W (y0 , p0 ) 2π h¯ Zt i ˆ dτ U (τ ) |V . × ΓDA W (yt , pt ) exp h¯ 0

(A.25)

Extending to the N-dimensional case, we obtain: Z 1 N dR0 dP0 Γˆ AD ρˆ Deq W (R0 , P0 ) C(t) ≈ 2π h¯ Zt i ˆ dτ U (τ ) |V . × ΓDA W (Rt , Pt ) exp h¯ 0

(A.26)

Appendix B Proof of identity Equation A.22 In this appendix, we provide the proof for the (purely classical) identity in Eq. (A.22). To this end, we start out by noting that the Jacobian of changing integration variables from y1 , . . . , yN−1 to f1 , . . ., fN−1 is given by: −1 ∂ f1 /∂ y1 · · · ∂ f1 /∂ fN−1 −1 ∂y ∂ f . . . = = .. .. .. ∂ f ∂y ∂ fN−1 /∂ y1 · · · ∂ fN−1 /∂ yN−1


(B.1)

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


The ( j, k)th element of the matrix

∂f ∂y

is given by

   − εm2 ,   ∂ fj  m ′′ = 2 ε 2 −V (y j ), ∂ yk      0,

Hence,

k = j±1

m N−1

ε2

det σN−1 ≡

m N−1

ε2

(B.2)

Other

  2 −1 0 · · · 0     −1 2 −1 · · · 0    N−1   ∂ f m .  0 −1 . . · · · 0  =   ∂y ε2   .. ..   .. . . 2 −1  .   0 0 0 −1 2  ′′ V (y1 )   ′′  V (y2 )  ε2  .. −  .  m  ′′ V (yN−2 )   ≡

.

k= j

gN−1 .

            ′′ V (yN−1 )

(B.3)

Now, consider the truncated j × j matrix σ j that consist of the first j rows and first j columns of σN−1 . It can be shown that this matrix satisfies the following recursion relation: 63 ε 2 ′′ g j+1 = 2 − V (y j+1 ) g j − g j−1 , m


(B.4)


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

since g j+1 =

0     . ..   σ j       0     −1 h i 2 ′′ 2 − εm V (y j+1 ) 0 · · · 0 −1 ε 2 ′′ = det σ j · 2 − V (y j+1 ) − det σ j−1 m 2 ε ′′ =g j · 2 − V (y j+1 ) − g j−1 . m 



(B.5)

(B.6) (B.7)

and

ε 2 ′′ g1 =2 − V (y1 ) m 2 ′′ ε 2 − V (y1 ) −1 m g2 = ε 2 ′′ −1 2 − m V (y2 ) ε 2 ′′ =g1 2 − V (y2 ) − 1, (g0 = 1) m

(B.8)

(B.9)

Rearranging the recursion relation, Eq. (B.4), we obtain: g j+1 + g j−1 − 2g j 1 ′′ = − V (y j+1 )g j . 2 ε m

(B.10)

In the limit ε → 0, Eq. (B.10) is equivalent to ϕ (t) = ε g j satisfying the following differential equation:

′′

d2 ϕ (t) V (t) =− ϕ (t) , 2 dt m


(B.11)

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


with initial conditions

ϕ (0) = ε g0 → 0 ε 2 ′′ dϕ (0) g1 − g0 =ε = 2 − V (y1 ) − 1 → 1 dt ε m

(B.12) (B.13)

In the next paragraph we show that m[∂ y(y0 , p0 ,t)/∂ p0 ] satisfies the same differential equation with the same initial conditions as ϕ (t), and that therefore m[∂ y(y0 , p0 ,t)/∂ p0] = ϕ (t). To this end, consider Lagrangian for particle moving on the average potential: 1 L = my˙2 −V (y) . 2

(B.14)

Next, define the Jacobian J(y0 , p0 ,t):

J(y0 , p0 ,t) ≡

∂ y(y0 , p0 ,t) . ∂ p0

(B.15)

Taking partial derivative with respect to p0 on both sides of Eular–Lagrange equation (see Ref. 63 Eqs. (12.33)–(12.36) on page 86), d dt

∂L ∂ y˙

−

∂L =0 , ∂y

(B.16)

then leads to the following equation: d dt

2 d ∂ L ∂ 2L ˙ ∂ 2L − J=0 . J + ∂ y˙2 dt ∂ y∂ y˙ ∂ y2

(B.17)

Substituting the Lagrangian from Eq. (B.14) into Eq. (B.17), we obtain an equation of motion for J: mJ¨+

∂ 2V (y) J=0 . ∂ y2


(B.18)


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

The initial conditions are

y(0) = a for all p ⇒ J(0) = 0 , p=

∂L ∂ y(0) ˙ 1 ˙ = 1 . = ⇒ J(0) = my(0) ˙ ⇒ ∂ y˙ ∂ p0 m m

(B.19) (B.20)

Comparing equation of motion and initial conditions of Jacobian J with those of ϕ (t), we observe that both obey the same differential equation of motion (Eq. (B.11) and Eq. (B.18)) and the same initial conditions except for a scale factor m for J. We can therefore conclude that

ϕ (t) = mJ(y0 , p0 ,t) .

(B.21)

1 ∂ p0 1 ∂ yt −1 1 1 = = = m ∂ yt m ∂ p0 mJ(y0 , p0 ,t) ϕ (t) 1 m N−1 ∂ y 1 1 = N−1 = = ∂ f . ε gN−1 ∂ f ε ε2 ε2 ε m ∂y

(B.22)

Finally,


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


Appendix C Proof that the LSC approximation for equilibrium FGR rate constant in the Condon case coincides with the exact quantum-mechanical result for a model consisting of shifted parabolic PESs Consider a model where the donor and acceptor Hamiltonians are given by (one-dimensional notation is used for simplicity): pˆ2 1 Hˆ D = + ω 2 xˆ2 − ∆E , 2 2 pˆ2 1 2 ˆ HA = + ω (xˆ − d)2 . 2 2

(C.1) (C.2)

Within the Condon approximation, the equilibrium FGR rate constant is given by

kD→A =

|ΓDA |2 h¯ 2

Z ∞

−∞

D i E i i ˆ ˆ dt e− h¯ ∆Et e h¯ HD′ t e− h¯ HAt

D

.

(C.3)

Here, pˆ2 1 2 2 + ω xˆ . Hˆ D′ =Hˆ D + ∆E = 2 2

(C.4)

E D i i ˆ ˆ The LSC expression for C1 (t) = e h¯ HD′ t e− h¯ HAt is given by (see Eq. (19)): D

C1 (t) =

Z

eq dx0 dp0 ρD,W (x0 , p0 ) exp

Zt i av δ U (xτ ) dτ , h¯ 0

(C.5)

where xtav is obtained by propagating from initial state (x0 , p0 ) for a time t on the average potential surface, 1 2 1 d 2 1 2 2 + ω d . V (x) = [VD′ (x) +VA (x)] = ω x − 2 2 2 8 35 ACS Paragon Plus Environment

(C.6)


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

Hence, xtav

p0 d d = + x0 − cos ω t + sin ω t . 2 2 ω

(C.7)

Let 1 δ U (x) =VD′ (x) −VA (x) = d ω 2 x − d 2 ω 2 . 2

(C.8)

Then, 1 δ U (xtav ) = d ω 2 xtav − d 2 ω 2 2 d p0 2 = dω cos ω t + sin ω t , x0 − 2 ω

(C.9)

and Z t 0

δ U (xτav ) dτ

d = d ω x0 − sin ω t + p0 d(1 − cos ω t) . 2

(C.10)

Finally, the initial Wigner phase-space probability density is given by: eq ρD,W (x0 , p0 ) =

β h¯ ω 1 tanh π h¯ 2 2 β h¯ ω tanh × exp − HD′ (x0 , p0 ) . h¯ ω 2

(C.11)

Substituting Eq. (C.10) and Eq. (C.11) into Eq. (C.5), we obtain ω d2 β h¯ ω coth (1 − cos ω t) + i sin ω t . C1 (t) = exp − 2¯h 2

(C.12)

which coincides with the exact quantum result. The extension for the N-dimensional case is trivial since the harmonic modes can always be assumed uncoupled in the normal mode representation.


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


References (1) Marcus, R. A.; Sutin, N. Electron Transfers in Chemistry and Biology. Biochim. Biophys. Acta 1985, 811, 265–322. (2) Marcus, R. A. Electron transfer reactions in chemistry. Theory and experiment. Rev. Mod. Phys. 1993, 65, 599–610. (3) Marcus, R. A. Electron transfer past and future. Adv. Chem. Phys. 1999, 106, 1–6. (4) Warshel, A.; Hwang, J. K. Simulation of the dynamics of electron transfer reactions in polar solvents: semiclassical trajectories and dispersed polaron approaches. J. Chem. Phys. 1986, 84, 4938–4957. (5) Jortner, J.; Bixon, M. Intramolecular vibrational excitations accompanying solvent-controlled electron transfer reactions. J. Chem. Phys. 1988, 88, 167–170. (6) Kuharski, R. A.; Bader, J. S.; Chandler, D.; Sprik, M.; Klein, M. L.; Impey, R. W. Molecular model for aqueous ferrous-ferric electron transfer. J. Chem. Phys. 1988, 89, 3248–3257. (7) Bader, J. S.; Kuharski, R.; Chandler, D. Role of nuclear tunneling in aqueous ferrous-ferric electron transfer. J. Chem. Phys. 1990, 93, 230–236. (8) Newton, M. D. Quantum chemical probes of electron-transfer kinetics: the nature of donoracceptor interactions. Chem. Rev. 1991, 91, 767–792. (9) Barbara, P. F.; Meyer, T. J.; Ratner, M. A. Contemporary Issues in Electron Transfer Research. J. Phys. Chem. 1996, 100, 13148–13168. (10) Topaler, M.; Makri, N. Path integral calculation of quantum nonadiabatic rates in model condensed phase reactions. J. Phys. Chem. 1996, 100, 4430–4436. (11) Chandler, D. In Classical and Quantum Dynamics in Condensed Phase Simulations;



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

Berne, B. J., Cicootti, G., Coker, D. F., Eds.; World Scientific: New Jersey, 1997; Chapter 2, pp 25–50. (12) Weiss, U. Quantum dissipative systems; World Scientific: London, 1993. (13) May, V.; Kühn, O. Charge and energy transfer dynamics in molecular systems; Willey-VCH: Berlin, 2000. (14) 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. (15) Nitzan, A. Chemical Dynamics in Condensed Phases; Oxford University Press: New York, 2006. (16) Jang, S.; Newton, M. D. Closed-Form Expressions of Quantum Electron Transfer Rate Based on the Stationary-Phase Approximation. J. Phys. Chem. B 2006, 110, 18996–19003. (17) Shi, Q.; Geva, E. A self-consistent treatment of electron transfer in the limit of strong friction via the mixed quantum-classical Liouville method. J. Chem. Phys. 2009, 131, 034511. (18) Hazra, A.; Soudackov, A. V.; Hammes-Schiffer, S. Role of Solvent Dynamics in Ultrafast Photoinduced Proton-Coupled Electron Transfer Reactions in Solution. J. Phys. Chem. B 2010, 114, 12319–12332. (19) Menzeleev, A. R.; Ananth, N.; Miller, T. F., III Direct simulation of electron transfer using ring polymer molecular dynamics: Comparison with semiclassical instanton theory and exact quantum methods. J. Chem. Phys. 2011, 135, 074106. (20) 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.


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


(21) Landry, B. R.; Subotnik, J. E. Communication: Standard surface hopping predicts incorrect scaling for Marcus’ golden-rule rate: The decoherence problem cannot be ignored. J. Chem. Phys. 2011, 135, 191101. (22) Landry, B. R.; Subotnik, J. E. How to recover Marcus theory with fewest switches surface hopping: Add just a touch of decoherence. J. Chem. Phys. 2012, 137, 22A513. (23) Huo, P.; Miller, T. F., III; Coker, D. F. Predictive partial linearized path integral simulation of condensed phase electron transfer dynamics. J. Chem. Phys. 2013, 139, 151103. (24) 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. (25) 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. (26) 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. (27) 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. (28) 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.



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

(29) 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. (30) Zhao, Y.; Liang, W. Charge transfer in organic molecules for solar cells: Theoretical perspective. Chem. Soc. Rev. 2012, 41, 1075–1087. (31) 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. (32) 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. (33) 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. (34) 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. (35) Jortner, J. Temperature dependent activation energy for electron transfer between biological molecules. J. Chem. Phys. 1976, 64, 4860–4867. (36) Mikkelsen, K. V.; Ratner, M. A. Electron tunneling in solid-state electron-transfer reactions. Chem. Rev. 1987, 87, 113–153. (37) Newton, M. D.; Sutin, N. Electron Transfer Reactions in Condensed Phases. Ann. Rev. Phys. Chem. 1984, 35, 437–480.


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


(38) Marcus, R. A. On the Theory of Oxidation-Reduction Reactions Involving Electron Transfer. I. J. Chem. Phys. 1956, 24, 966–978. (39) Marcus, R. A. Electrostatic Free Energy and Other Properties of States Having Nonequilibrium Polarization. I. J. Chem. Phys. 1956, 24, 979–989. (40) De Vault, D.; Chance, B. Studies of Photosynthesis Using a Pulsed Laser. Biophys. J. 1996, 6, 825–847. (41) Izmaylov, A. F.; Mendive-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. (42) 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. (43) Egorov, S. A.; Rabani, E.; Berne, B. J. Nonradiative relaxation processes in condensed phases: Quantum vs. classical baths. J. Chem. Phys. 1999, 110, 5238–5248. (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.



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

(48) 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. (49) 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. (50) Shi, Q.; Geva, E. A comparison between different semiclassical approximations for optical response functions in nonpolar liquid solutions. J. Chem. Phys. 2005, 122, 064506. (51) 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. (52) 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. (53) 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. (54) 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. (55) 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. (56) 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.


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


(57) 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. (58) Cho, M.; Silbey, R. J. Nonequilibrium photoinduced electron transfer. J. Chem. Phys. 1995, 103, 595–606. (59) Petit, A. S.; Subotnik, J. E. How to calculate linear absorption spectra with lifetime broadening using fewest switches surface hopping trajectories: A simple generalization of groundstate Kubo theory. J. Chem. Phys. 2014, 141, 014107. (60) Petit, A. S.; Subotnik, J. E. Appraisal of Surface Hopping as a Tool for Modeling Condensed Phase Linear Absorption Spectra. J. Chem. Theory Comput. 2015, 11, 4328–4341. (61) Stratt, R. M. The instantaneous normal modes of liquids. Acc. Chem. Res. 1995, 28, 201–207. (62) Garg, A.; Onuchic, J. N.; Ambegaokar, V. Effect of friction on electron transfer in biomolecules. J. Chem. Phys. 1985, 83, 4491–4503. (63) Schulman, L. S. Techniques and applications of path integration; Wiley: New York, 1981.



Figures and captions (Note: Figure 1–6 and Figure 8 should be two-column width format.)

η/ωop= 0.2

0

η/ωop= 1

η/ωop= 5 kBT ___ = 0.2

-5

log (k/ωop )

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

ℏωop

0

kBT ___ =1

-5

ℏωop

0

-10

kBT ___ =5

W-AV C-AV C-D W-0 Marcus

-5

0

5

0

5

ωDA / ωop

0

5

ℏωop

10

Figure 1. A comparison of FGR rate constants kD→A (on a semilog scale, with log base 10 hereinafter) as a function of ωDA = −∆E/¯h (with ∆E being the donor-to-acceptor reaction free energy), as obtained via different methods, at the indicated temperatures (T ) and frictions (η ), for the spectral density J1 (ω ), Eq. (39). It should be noted that the LSC result coincides with the fully quantum-mechanical result (both labeled W-AV here).


Page 45 of 55

Ω = 0.2

η = 0.5

0

η=1

Er = 0.08

η=5


-6

kBT = 0.2

0

log (k )

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


kBT = 1

-6 0

kBT = 5

-6 -12

0

5

10

0

5

10

0

5

10

15

ωDA Figure 2. A comparison of FGR rate constants kD→A (on a semilog scale) as a function of ωDA = −∆E/¯h (with ∆E being the donor-to-acceptor reaction free energy), as obtained via different methods, at the indicated temperatures (T ) and frictions (η ), for the spectral density J2 (ω ), Eq. (40) with Ω = 0.2 and Er = 0.08. It should be noted that the LSC result coincides with the fully quantum-mechanical result (both labeled W-AV here).



η = 0.5

0

η=1

Ω = 0.5 Er = 0.5

η=5


-5

kBT = 0.2

0

log (k )

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 46 of 55

kBT = 1

-5 0

kBT = 5

-5 -10

0

5

10

0

5

10

0

5

10

15



Page 47 of 55

η = 0.5

0

η=1

Ω=1 Er = 2

η=5

kBT = 0.2

-5 0

log (k )

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


kBT = 1

-5 0 W-AV C-AV C-D W-0 Marcus

-5 -10

0

5

10

0

5

10

0

5

10

kBT = 5 15

ωDA Figure 4. A comparison of FGR rate constants kD→A (on a semilog scale) as a function of ωDA = −∆E/¯h (with ∆E being the donor-to-acceptor reaction free energy), as obtained via different methods, at the indicated temperatures (T ) and frictions (η ), for the spectral density J2 (ω ), Eq. (40) with Ω = 1.0 and Er = 2. It should be noted that the LSC result coincides with the fully quantum-mechanical result (both labeled W-AV here).



η = 0.5

0

η=1

Ω = 0.2 Er = 0.8

η=5


-5

kBT = 0.2

0

log (k )

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

-5

kBT = 1

0 -5

kBT = 5

-10 0

5

10

0

5

10

0

5

10

15



Page 49 of 55

η = 0.5

0

η=1

Ω = 0.2 Er = 8

η=5

kBT = 0.2

-10 0

log (k )

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


kBT = 1

-10 0 W-AV C-AV C-D W-0 Marcus

-10 -20

0

20

0

20

0

20

kBT = 5 40

ωDA Figure 6. A comparison of FGR rate constants kD→A (on a semilog scale) as a function of ωDA = −∆E/¯h (with ∆E being the donor-to-acceptor reaction free energy), as obtained via different methods, at the indicated temperatures (T ) and frictions (η ), for the spectral density J2 (ω ), Eq. (40) with Ω = 0.2 and Er = 8. It should be noted that the LSC result coincides with the fully quantum-mechanical result (both labeled W-AV here).



0

-10

log (k )

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


-20

-30

0

100

ωDA

200

300

Figure 7. A comparison of FGR rate constants kD→A (on a semilog scale) as a function of ωDA = −∆E/¯h (with ∆E being the donor-to-acceptor reaction free energy), as obtained via different methods, for the spectral density J2 (ω ), Eq. (40), with Ω = ωc = 1, Er = 68.3, y0 = 5.84, η = 3.43, kB T = 2.71, ΓDA = 0.143. It should be noted that the LSC result coincides with the fully quantum-mechanical result (both labeled as W-AV here).


Page 51 of 55

η/ωop= 0.2

η/ωop= 1

η/ωop= 5

0 exact 5e4 5e6

-2

kBT ___ = 0.2

ℏωop

-4

log (k / ωop )

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 -2

kBT ___ =1

ℏωop

-4 0 -2

kBT ___ =5

ℏωop

-4 -6

0

2

4

0

2

4

0

2

4

6

ωDA / ωop Figure 8. A comparison of LSC FGR rate constants kD→A (on a semilog scale) as a function of ωDA = −∆E/¯h (with ∆E being the donor-to-acceptor reaction free energy), as obtained analytically (exact) and numerically with NMC = 5 × 104 and 5 × 106 , at the indicated temperatures (T ) and frictions (η ), for the spectral density J1 (ω ), Eq. (39).



exact 5e4

-2 log (k )

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

-4

-6

0

5 ωDA

10

Figure 9. A comparison of LSC FGR rate constants kD→A (on a semilog scale) as a function of ωDA = −∆E/¯h (with ∆E being the donor-to-acceptor reaction free energy), as obtained analytically (exact) and numerically with NMC = 5 × 104 , for the spectral density J2 (ω ), Eq. (40), with primary mode frequency Ω = 1 and Er = 2 at temperature kB T = 1 and friction η = 1.


Page 53 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


Figure 10. The relative error averaged over the ωDA range considered here as function of NMC , for J1 (ω ), at the indicated temperature (T ) and friction (η ).



0

log (error)

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

-1

-2

0.5

1 1.5 #Traj. x 106

2

2.5

Figure 11. The relative error averaged over the ωDA range considered here as function of NMC , for J2 (ω ) with primary mode frequency Ω = 1 and Er = 2 at temperature kB T = 1 and friction η = 1.


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 η = 0.5

0

η=1

η=5

LSC Marcus

-6

log (k )

Page 55 of 55

kBT = 0.2

0

kBT = 1

-6 0

kBT = 5

-6 -12

0

5

10

0

5

10

ωDA

0

5

10

15


Equilibrium Fermi's Golden Rule Charge Transfer Rate Constants in

Recommend Documents