J. Phys. Chem. B 2008, 112, 3735-3745
3735
Solvent Dynamics Effect in Condensed-Phase Electron-Transfer Reactions Jianjun Zhu* Department of Chemistry, Henna Normal UniVersity, Xinxian, People’s Republic of China, and Department of Chemistry, State UniVersity of New York at Stony Brook, Stony Brook, New York 11796
Yanbin Cheng Department of Chemistry, Henna Normal UniVersity, Xinxian, People’s Republic of China
Tong-Chun Bai Department of Chemistry, Suzhou UniVersity, Suzhou, People’s Republic of China
Yan Lu and Zhaorong Chang Department of Chemistry, Henna Normal UniVersity, Xinxian, People’s Republic of China
Dongqing Wei Department of Bioinformics, Shanghai Jiao Tong UniVersity, Shanghai, People’s Republic of China
George Stell Department of Chemistry, State UniVersity of New York at Stony Brook, Stony Brook, New York 11796 ReceiVed: September 22, 2007; In Final Form: January 6, 2008
Channel-based reaction-diffusion equations are solved analytically for two electron transfer (ET) models, where the fast inner-sphere coordinate leads to an ET reaction treated by Fermi’s golden rule, and the slow solvent coordinate moves via diffusion. The analytic solution has let us derive an ET rate constant that modifies the Marcus-Jortner formula by adding a constant R which we call a dynamic correction factor. The dynamic correction factor measures the effect of solvent friction. When the relaxation of solvent dynamics is fast, the dynamic correction can be neglected and the ET rate constant reduces to the traditional Marcus-Jortner formula. If the solvent dynamic relaxation is slow, the dynamic correction can be very large and the ET rate can be reduced by orders of magnitude. Using a generalized Zusman-Sumi-Marcus model as a starting point, we introduce two variants, GZSM-A and GZSM-B, where in model A, only one quantum mode is considered for inner-sphere motion and in model B, a classical mode for inner-sphere motion is added. By comparing the two models with experimental data, it is shown that model B is better than model A. For the solvents that have a relaxation time ranging between 0 and 5 ps, our result agrees fairly well with experimental data; for the solvents that have a relaxation time ranging between 5 and 40 ps, our result deviates from the experimental values. After introducing an adjustable scaling index in the effective time correlation function of the reaction coordinate, good agreement is achieved between the experiment and the theory for model B for all of the solvents studied in this paper.
I. Introduction Charge transfer, including electron transfer and proton transfer, in biological systems is a well-known topic. There are many experiments carried out in this field and there have been quite a few review articles on charge-transfer reactions in biomolecules.1-3 Our understanding of the electron transfer (ET) process is based on the Marcus-Levich-Jortner4-7 theory. The Marcus4-5 theory is classical, based on solvent reorganization in the presence of the external electrical field, and the LevichJortner6-7 theory is quantum mechanically based on Fermi’s golden rule. Marcus’ theory has led to an outer-sphere solventreorganization energy in the rate process, which collectively describes the contribution from the slow solvent modes. Fermi’s * Corresponding author.
golden-rule has led to a rate-process theory with contributions from all solvent modes. In the classical limit, Fermi’s goldenrule theory reduces to Marcus theory.3d Despite the complication of Fermi’s golden-rule theory, it has proven to be the most successful one in biological systems.3d,8 As the exact quantum calculation of the biological system is overwhelmingly complicated, the test of Fermi’s golden rule theory has been limited to fitting the experimental data with a few parameters.2,3,8 It is known that the classical Marcus theory deviates from experimental data seriously because of the neglect of highfrequency modes.2,3,8 The quantum theory catches the contribution from all modes but is too complicated to use. A hybrid ET picture to bridge these two limits has been proposed1,4b and employed by many people.7,9-18 The spirit of the hybrid ET theory is to conceptually divide the solvent into two groups,
10.1021/jp077637q CCC: $40.75 © 2008 American Chemical Society Published on Web 02/29/2008
3736 J. Phys. Chem. B, Vol. 112, No. 12, 2008
Zhu et al.
one consists of the inner-sphere solvents which will be called “ligands” in this paper and the other consists of the outer-sphere solvents. If the electron donor and accepter are in one molecule, all of the functional groups inside the molecule will be called “inner-sphere solvents”. The motion of the inner-sphere solvent could be very complicated. Harmonic-oscillator models for ligands have been usually used for theoretical studies, which has high-frequency modes and low-frequency modes for big molecules. The motion of the outer-sphere solvent is usually in slow modes. This hybrid ET picture has been incorporated into Fermi’s golden-rule theory by Jortner.7 The hybrid ET theory seems the most successful one so far in the application to the biological ET systems,2,3,8 where a high-frequency mode is usually used to characterize the inner-sphere ligand motion and a low-frequency mode is usually used to characterize the outersphere solvent motion. Another extension of the hybrid ET theory incorporates the solvent dynamic effect. Since the outer-sphere solvent is described collectively by one slow mode, its random nature has been neglected in the traditional Marcus theory. The reaction coordinate dynamics, or solvent dynamics, was first introduced by Zusman9 to ET theory. Since then, the solvent dynamic effect on electron transfer has received a lot of attention in the past three decades.10-18 Zusman’s theory is based on a diffusionreaction equation model, which has been a well-known model in the field of diffusion-limited reactions.19 In Zusman’s model, the reaction provides sinks to the diffusion of the solvent coordinate. Sumi-Marcus13 has enriched Zusman’s model by incorporating the hybrid ET concept, where the reaction part is induced by the fast ligand-vibration motion, and the diffusion part is from the slow solvent coordinate. In a series of papers, we have studied the solvent dynamic effect on ET reactions by solving the diffusion-reaction equation.14-15 The ZusmanSumi-Marcus model has been generalized to a variety of systems,14-18 and Fermis’ golden rule has been used to formulate the rate constant for the reaction part.14c This has led to the so-called generalized Zusman-Sumi-Marcus (GZSM) model.15 The solution of the GZSM model at high temperature has been given analytically, which leads correctly to both the SumiMarcus’ result and Zusman’s result in different limits. A dynamically corrected rate constant has been constructed from the solution of the GZSM model14-15
k)
ke 1+R
(1.1)
where we call R a “dynamic correction factor”. The R goes to 0 in the thermal-equilibrium limit and can be very large if the reaction coordinate is far away from an equilibrium distribution. ke is Marcus’s rate constant without considering the solvent dynamic effect. In contrast to the analytic approach, a numerical approach to solving Zusman-Sumi-Marcus model has also been employed.13,16-18 In this paper, we will continue our study of the solvent dynamic effect on electron-transfer reactions by focusing on biological systems. The reaction part will be fully treated quantum mechanically with the slow coordinate as a parameter. In the thermal-equilibrium limit, we will show that the thermalequilibrium rate constant is the same as that derived by Jorner et al.3a,7 We will split the GZSM model into two models called GZSM-A and GZSM-B, where in GZSM-A, only one quantum mode is considered for the inner-sphere motion and in GZSMB, one quantum mode and one classical mode are considered for the inner-sphere motion. We will also show that the
diffusion-reaction equation for GZSM-A can be solved exactly, and the diffusion-reaction equation for GZSM-B can be solved approximately. From the analytical solution of the diffusionreaction equation, we will show that eq 1.1 is still valid except that the dynamic correction factor R is replaced by its statistical average. Our goal is to provide a simple analytic formula that can be used to describe and explain the experimental data, where the influence of solvent dynamics on biological ET reactions can be known by inferring the value of R. There are many experiments16,20-21 and computations22,23 carried out for bio-molecular electron-transfer reactions in recent years. The match between theory and experiments is still a challenge because the theoretical models are either too simple or too computationally complicated. A better theory is needed to explain and correlate the experimental data. Mean first-passage time theory24 has been used to calculate the ET rate by a few groups13,16-18 where the diffusion-reaction equation has been solved numerically. Since mean first-passage time theory relies on numerical methods, it is not user-friendly, and it is not easy to factor out the solvent dynamic effect. Jortner and Bixon7c have developed an ET theory to consider the solvent dynamic effect in which the path-integral method and the Landau-Zener theory have been used. The Bixon-Jortner formulation has a clear analytic structure for the ET rate but seriously deviates from experiment.16 This is partly because of the complicated nature of the problem and partly because of the oversimplified dynamics for the reaction coordinate. The time correlation function of the reaction coordinate plays a crucial rule in the study of solvent dynamic effect.11-14 Hynes used a generalized Langevin equation to study the reactioncoordinate dynamics. He showed that the Langevin dynamics reduces to diffusion dynamics in the over-damped limit.11 The time correlation function of the reaction coordinate plays a crucial role in the analytic solution of the diffusion equations. Measuring the dynamics of the Born solvation energy has been a common practice24 in the experimental field, which has provided the needed information for the reaction-coordinate dynamics. It has been shown by us that the time correlation function of the reaction coordinate can be extracted from the dynamics of the Born salvation energy.14b We will incorporate the reaction coordinate dynamics in our formulation of ET theory. Multi-time-relaxation dynamics and multidimensional vibrational coordinates are introduced into ET theory by Zhu and Rasaiah14c and are deployed by computer simulation for Beaton30 by Rossky’s group.23 It is declared that intramolecular torsional dynamics also plays an important role in their calculation of the ET rate constant for Beaton 30. This paper is organized as follows: in section II, the freeenergy surfaces and the channel-based rate constants are discussed for our two specific models, model A and model B. In section III, the channel-based diffusion-reaction equation is solved by using a Green’s function method for both models, and the thermal equilibrium rate constant and the dynamically corrected rate constant are constructed. In section IV, numerical results are discussed, and a comparison between theory and experiment is given. II. Free-Energy Surfaces and Channel Rate Constant We will adopt the hybrid ET picture proposed by Marcus and the others,4,13-17 where x is used to represent the outersphere solvent coordinate, and q to represent the inner-sphere vibration coordinate. The explicit definition of the x coordinate
Solvent Dynamics Effect
J. Phys. Chem. B, Vol. 112, No. 12, 2008 3737
has been discussed by us in several places.14 Since the ligand vibration is faster than solvent reorientation, the fast coordinate q is separated from the slow coordinate x in the GZSM model. By using qi to represent the dimensionless coordinate for ith ligand and x to represent the polar solvent coordinate, the free energy surfaces can be written as14e,15 N
ai
(2.1a)
(qi - 1)2 + (x - x0)2/2 + ∆G0 ∑ i)1 2
(2.1b)
N
F2 )
bi
ai ) Mi(ωi ∆q0i )2 ) 2λi
(2.2a)
bi ) Mi(ω′i ∆q0i )2 ) 2λi′
(2.2b)
In eq 2.2, Mi is the reduced mass for ith ligand bond, ωi and ωi′ are the bond vibration frequencies corresponding to reactant and product state respectively, λiand λi′ are the bond reorganization energies for the forward and backward reactions, and ∆q0i , which has a dimension, is the shift of equilibrium vibration coordinate for bond i. In eq 2.1, x0 is the equilibrium value of x, which has been defined clearly by us previously for a dielectric continuum solvent.14 We have
4π c
∫
|P Bor(b,t) r -B Por,1eq(b)| r 2 db r
(2.3a)
In eq 2.3, B Por is the oriental polarization of the polar fluid,
x02 )
4π c
r -B Por,1eq(b)| r 2 db r ∫|PBor,2eq(b)
Fn |〈χm|χn〉|2 δ(Em(x) - En(x)) ∑ ∑ n)0 m)0
F1 )
Em(x) - En(x) ) mpωk′ - npωk - ∆G0(x)
N
F2 )
bi
(qi - 1) ∑ i)1 2
2
δ(x) )
+ ∆G (x)
(2.7)
1 2π
∫ dt exp(itx)
(2.8)
we have
|V12|2
R(x,k) )
2
p
(
∫ dt exp -
)
∆G0(x) it fk(t) p
(2.9)
where ∆G0(x) is the free energy gap defined in eq 3.4 and fk(t) is the Frank-Condon overlap integral for ligand k3c. So far, only one ligand bond k is considered in the derivation of eq 2.9. If all of the ligand bonds are considered, eq 2.9 is generalized to
R(x) )
|V12|2 2
p
(
∫ dt exp -
)
∆G0(x) it f(t) p
(2.10)
where f(t) is given by ref 3c with ωk ) ωk′ for the symmetric case N
f(t) )
fk(t) ) e-F ∑ Ak exp[∑ imkωkt] ∏ k k)1
(2.11)
{mK}
with
[ ]
(2.12)
Rk ) 2Sk[njk(njk + 1)]1/2
(2.13)
F ) Sk(2njk + 1)
(2.14)
Sk ) λk/pωk
(2.15)
Ak )
(2.4a)
0
(2.6)
and the delta function definition
ai
qi2 ∑ i)1 2
(2.5)
exp(-βEn) ) exp(-nβpωk)[1 - exp(-βpωk)] Q
eq
N
∞
In eq 2.6, Q is the partition function. Using the relation
(2.3b)
and B Por,i is the equilibrium value of B Por at state i. The outersphere solvent reorganization energy λo is given by λo ) x20/2. The justification of eq 2.3 for a dielectric continuum solvent has been given elsewhere.14a We will next discuss two hybrid models. a. Channel Rate Constants for Model A. We assume that all of the inner-sphere coordinates are fast compared with the outer-sphere coordinate x for model A. In the GZSM model, the reaction is caused by motion of the fast coordinates with the slow coordinate x seen as frozen. One can parametrically put the slow coordinate x into the free-energy gap. Subtracting x2/2, eq 2.1 can be rewritten as14c
∞
where V12 is the matrix coupling element, χn is the nuclei wave function for eigen state n in reactant well, χm is the nuclei wave function for eigen state m in the product well, and Fn is the population probability defined by
Fn )
where ∆G0 is the free energy gap for the reaction, and
x2(t) )
R(x,k) ) 2π|V12|2 p
qi2 + x2/2 ∑ i)1 2
F1 )
R corresponding to one ligand bond k with vibration frequencies ωk and ωk′ can be written as3,14c
(2.4b)
∏ k
njk + 1
Imk(Rk)
mk/2
njk
where Im(x) is the modified Bessel function.3c,7 Substituting eq 2.11 into eq 2.10, we have3c,7
where
∆G0(x) ) ∆G0 + λ0 - x x2λ0
(2.4c)
This definition of the solvent-coordinate-dependent free-energy gap distinguishes our work from other groups.16-19 For the free energy surfaces given in eq 2.4, Fermi’s golden rule rate constant
RS(x) ) |V12|2 e-F Ak 2 {mk} p
∑ ∫
dt exp
{∑ k
mkωkit +
∆G0(x)
}
it p
(2.16)
3738 J. Phys. Chem. B, Vol. 112, No. 12, 2008
Zhu et al.
Bm )
By using eq 2.8 for the delta function, we have
RS(x) )
2π|V12|2 px2λ0
e-F
(1 - ξ ) m![S(1 - ξ)] 2 1/2
∑ Akδ[x - xk}
(2.17)
{mk}
where the superscript S stands for symmetric and xk is defined by
∆G0 + λ0 + xk )
R (x) )
2π|V12|2 px2λ0
e-F
(2.18)
x2λ0
Amδ[x - xm} ) ∑ RSm(x) ∑ m m
(2.19)
where
xm )
∆G0 + λ0 + mpω
(2.20)
x2λ0
Am ) Im(R1)
[ ] nj1 + 1 nj1
m/2
(2.21)
and the channel rate constant RSm(x) is defined by
RSm(x)
)
2π|V12|2 px2λ0
e-FAmδ(x - xm)
(2.22a)
It is seen from eq 2.19 that the coordinate dependent rate constant derived from Fermis’ golden rule theory contains an infinite number of delta functions, which are usually called sinks in the diffusion-limited reactions. Now each sink provides a physical meaningful rate constant Rm(x) called the channel rate constant, which comes from the picture that a quantum transition from an initial state can go to any final state m. Each m state forms a particular channel for the quantum transition and couples to the slow solvent coordinate x. This channel-based picture has been first proposed by Jortner and Bixon.7c At high-temperature limit, eq 2.19 correctly reduce to the classical transition state result 14c
RS(x) )
2π|V12|2 px2πkBTλq
exp
{
}
[λq + ∆G0(x)]2 4kBTλq
(2.22b)
The Frank-Condon overlap integral in the asymmetric case with ωk! ) ωk′ has been worked out for one vibration mode previously,7b where the rate constant can be written as
RA(x) ) |V12|2 p2 with
e-S
Bm ∫ dt exp ∑ m
{
itmω′ +
|∑ m/2
[2S(1 - ξ)/ξ]-r
r)0
r!(m - 2r)!
|
2
(2.24a)
ξ ) (ω - ω′)/(ω + ω′)
(2.24b)
S ) Mωω′/p(ω + ω′)
(2.24c)
Using the definition of the delta function, we have
∑k mkpωk
Equation 2.17 and eq 2.18 give the Fermi’s golden-rule rate constant in the symmetric case with ωk ) ωk′ for the frozen coordinate x. If there is only one vibration mode, eq 2.17 reduces to S
m
∆G0(x)
}
it p
R (x) ) A
∑ m
RAm(x)
2π|V12|2 ) e-S p x2λ0
Bmδ[x - xm′} ∑ m
(2.25A)
where the channel rate constant RAm(x) is defined by
RAm(x) )
2π|V12|2 p x2λ0
e-SBmδ(x - x′m)
(2.25B)
and xm′ is the same as xm with ω replaced by ω′. It is easy to prove that Bm f Am when ωi ) ωi′. Although the classical transition state rate constant corresponding to eq 2.23 and eq 2.25 cannot be derived analytically, ref 15 has provided a numerical procedure to derive the classical transition state rate constant for the asymmetric case. Multisink-based R(x) has been used in reaction-diffusion models by other groups,19,21 but the sink constant and sink coordinate are derived here clearly without any ambiguity. Our quantum rate constant R(x) leads to the correct classical transition state limit which was not examined by other groups.19,21 b. Channel Rate Constants for Model B. For model B, we assume that the inner-sphere coordinates splits into a fast coordinate q and a slow coordinate y. This model was first introduced by Barbara’s group.16 The y coordinate is defined so that it is slow enough to be treated classically compared with the q coordinate, but fast enough to be treated as lying in thermal equilibrium distribution compared with the x coordinate. As in our previous treatment, we can parametrically put the slow coordinate y and x into the free energy gap,
a1 2 q 2
(2.26a)
a1 (q - 1)2 + ∆G0(y,x) 2
(2.26b)
F1 ) F2 )
where the new free energy gap is defined as
∆G0(y,x) ) ∆G0(x) + λy - y x2λy
(2.27)
and λy is the reorganization energy corresponding to the slow coordinate y. Applying Fermi’s golden rule to the q coordinate, we have
RS (x,y) )
2π|V12|2 p x2λy
e-F
Amδ[y - ym} ) ∑ RSm(y,x) ∑ m m (2.28)
where
(2.23) ym )
∆G0(x) + λy + mpω
x2λy
(2.29a)
Solvent Dynamics Effect
J. Phys. Chem. B, Vol. 112, No. 12, 2008 3739
As we have assumed that the probability distribution P(ym) for y coordinate always lies in thermal equilibrium, which is given by
P(ym) ) exp(-βym2/2)/x4kBT
(2.29b)
the coordinate dependent rate constant can be obtained by statistically averaging over the y coordinate
R (x) ) 〈R (x,y)〉 ) 2π|V12|2 λq (λq/pω)m S Rm(x) ) exp × pω m m! m px4πkBTλy exp[-βEm(x)] (2.30) S
The Smoluchowski operator has a time dependent diffusion coefficient D(t) which is related to the time correlation function of the reaction coordinate in the over damped limit by11,14b
D(t) ) -
1 d ln ∆(t) β dt
(3.4)
with
S
( )∑
∑
where the channel rate constant RSm(x) is defined by
RSm(x) ) 2π|V12|2 p x4πkBTλy
( )
exp -
λq (λq/pω)m exp[-βEm(x)] (2.30a) pω m!
and
Em )
(∆G0(x) + λy + mpω)2 4λy
∆(t) ) 〈x(t)x(0)〉/〈x2(0)〉
This relation between the time correlation function of the reaction coordinate and the diffusion constant has been discussed previously.11,14 It is seen that each quantum state m leads to a reaction channel for electron transfer. Equation 3.2 together with the channel-based rate constant RSm(x) constitutes our channelbased diffusion-reaction equation model. The advantage of this is that each channel-based diffusion-reaction equation can be solved exactly without introducing further complications.17 It should be pointed out that we have used the symmetric form RSm(x) for the reaction part in eq 3.2, but there is no reason that the asymmetric form RAm(x) cannot be used. The solution of the channel-based diffusion-reaction equation is given by a survival probability, which is defined by
(2.30b) Qm(t) )
By comparing eq 2.30a with eq 2.22a, it is seen that, in the expression of the channel rate constant, the delta function has been replaced by a Gaussian function in the channel rate constant. The advantage of model B over model A is that it adds on energy balance to reaction coordinate dynamics.
∫ dx Pm(x,t)
Q(t) )
V1(x) ) x2/2
(3.1a)
V2(x) ) (x - x0)2/2 + ∆G0 + mpω (m ) 0, 1, ...)
(3.1b)
where the x coordinate characterizes the nuclei motion for the classical system. This picture was first given by Jortner and Bixon7 and used by others.16,17 For each quantum state m, the channel-based reaction diffusion equation can be written as13,14
∂Pm(x,t) ) ∂t LPm (x,t) - RSm(x)Pm(x,t) (m ) 0, 1, 2, ..., N) (3.2) where we have neglected the backward reaction for simplicity, but one should keep in mind that it exists! In eq 3.2, Pm(x,t) is the probability to find the reaction system at coordinate x at time t for state m, and L is the generalized Smoluchowski operator corresponding to a harmonic potential well
L ) D(t)
[
d2 d dV +β dx dx d2x
]
(3.3)
(3.6)
The survival probability will be one if the reacting system is populated at reactant state and will be zero when time goes to infinity. The total survival probability can be written as
III. Solution of Channel-Based Reaction-Diffusion Equations a. The Channel-Based Reaction-Diffusion Equations. Since the reaction coordinate x characterizes the classical system and the coordinate q characterizes the quantum system, BornOpenheimer approximation tells us that the potential surfaces eq 2.1 are split into a series of parabolas when the quantum system lies in different eigen states where
(3.5)
Qm(t) ∏ m
(3.7)
Before we discuss a general solution, we will first discuss a static-state solution when the probability distribution is at thermal equilibrium. b. Thermal Equilibrium Rate Constant for Model A. The thermal equilibrium rate constant plays an important role in our theoretical formulation as it serves as a reference point when the effect of solvent dynamics can be neglected. We first derive the thermal equilibrium rate constant explicitly. If the relaxation of solvent polarization is very fast, the solvent can adjust itself instantaneously to the external field when electron-transfer takes place. In this limit, the probability distribution of the x coordinate does not depend on time, and thermal equilibrium distribution of x prevails. Solving eq 3.2, we have14a
Qm(t) ) exp(-km,eSt)
(3.8)
where km,eS is the thermal-equilibrium rate constant for channel m defined by
km,eS )
∫ dx RSm(x)Pmeq(x) )
1
x2πkBT
∫ dx RSm(x) e-βV(x)
(3.9)
where we have used a Boltzman distribution for the equilibrium probability distribution of x. Substituting eq 2.22a into eq 3.9 and carrying out the integral, we obtain the thermal equilibrium rate constant for channel m
3740 J. Phys. Chem. B, Vol. 112, No. 12, 2008
Zhu et al.
km,eS ) 2π|V12|2 p2 x4πkBTλ0
e-FAm exp[-βEm] (m ) 0, 1, 2, ...) (3.10)
ke )
km,e ∑ m
S
2π|V12|2
)
x4πkBTλ0
2
p
e-F
Am exp[-βEm] ∑ m
(3.11)
where
Em )
xm2 (∆G0 + λ0 + mpω1)2 ) 2 4λ0
(3.12)
If a further approximation is made that the vibration frequency is high, or the temperature is low, then βpω . 1, and we have the average quantum number nj1 ∼ 0, F ) S1, Am ) {S1m}/ {m!}, when eq 3.11 can be rewritten as3a,7
k eS )
km,eS ) ∑ m 2π|V12|2
( )∑
exp -
p x4πkBTλ0
∑ m
km,eA )
2π|V12|2
)
p x4πkBTλ
( )∑
exp -
λq
pω
m
(λq/pω)m m!
×
λq
pω
m
It is seen that the thermal equilibrium rate constants for model A and model B are very similar. Equation 3.18 reduces to eq 3.13 when the slow coordinates reduces to one. d. General Solution of Channel-Based Reaction-Diffusion Equation for Model A. The Green function method has been used extensively in the solution of diffusion equations.9,11b,14c We will use the Green-function method to find a general solution for the diffusion-reaction eq 3.2 for model A. In this section, we will limit our discussion to the symmetric case with one vibrational mode for a particular channel m, where RSm(x) is given by eq 2.22a. For channel m, the diffusion-reaction equation can be rewritten as
∂Pm(x,t) ) LPm (x,t) - RSm(x)Pm(x,t) ∂t
(3.19)
The Green function equation corresponding to the Smoluckowski operator L can be written as
(λq/pω)m exp[- βEm] (3.13) m!
dG(x,y,t) - LG(x,y,t) ) δ(x - y)δ(t) dt
(3.20)
The solution of the Green function is11,14
Equation 3.13 is the approximation that assumes one mode with a high and symmetric frequency, which has been used extensively in calculating the reaction rate in biological ET reactions.8 Inserting eq 2.25b into eq 3.9 and carrying out the integration, we have7a
k eA )
km,e ∑ m
S
exp[- βEm] (3.18)
which, by combining eq 3.7 and eq 3.8, leads to the well-known Fermi’s golden rule rate constant3a,7 S
ke ) S
2π|V12|
2
e-S
p x4πkBTλ0
Bm exp[-βEm′] ∑ m
(3.14)
where S is defined in eq 2.24c, Bm is defined in eq 2.24a, and Em′ is the same as Em except that ω is replaced by ω′. In the high vibration frequency approximation, βpω . 1, we have
[S(1 - ξ)]m Bm ) (1 - ξ2)1/2 m!
2π|V12| 2
p
(3.15)
2
x4πkBTλ
e-FAm exp[-βEm]
(3.16)
with
(∆G0 + λ + mpω)2 Em ) 4λ
[
x
β β (x - y∆(t)) exp 2 2 1 - ∆(t)2 2π(1 - ∆(t) )
2
]
(3.21)
where ∆ is the time correlation function of x coordinate. The Green function describes the propagation of the diffusion coordinate x which initially starts from a point source y. The general solution of eq 3.19 can be written as
Pm(x,t) ) Pmeq(x) -
∫0t dt′ ∫ dy RSm(y) Pm(y,t′) G(x|y,t - t′)
(3.22)
or in Laplace space
where ξ is defined by eq 2.24b. It is seen that eq 3.14 reduces to eq 3.11 when the vibration frequencies are the same (symmetric). The thermal equilibrium rate constants for model A are the same as those derived by Jortner and his co-workers.7 c. Thermal Equilibrium Rate Constant for Model B. To calculate this rate constant for model B, we substitute eq 2.30a into eq 3.9 and carry out the integral, when we have the thermal equilibrium rate constant for channel m
km,eS )
G(x,y,t) )
Pm(x,s) ) Pmeq(x)/s -
where the total reorganization energy of the classical system is defined by λ ) λo + λy. Summing over all channels, we obtain
(3.23)
where we have assumed an equilibrium distribution for the initial condition. Using eq 2.22a for RSm(y), we have
Pmeq(x) 2π|V12|2 -F e AmG(x|xm,s} Pm (xm,s) Pm(x,s) ) s p x2λ0 (3.24) where G(x|xm,s} is the Laplace transform of eq 3.21 with the point source replaced by the sink point xm. Integrating eq 3.24 with respect to the x coordinate, we obtain the survival probability
Qm(s) ) (3.17)
∫ dy RSm(y) Pm(y,s) G(x|y,s)
[
2π|V12|2 -F 1 1e AmPm (xm,s) s p 2λ
x
0
]
(3.25)
where F is given in eq 2.14, and λ0 is the outer-sphere solvent reorganization energy. Solving eq 3.24 and eq 3.25 together, we have
Solvent Dynamics Effect
J. Phys. Chem. B, Vol. 112, No. 12, 2008 3741
S km,e 1 Qm(s) ) - 2 s s [1 + a(s)]
(3.26)
where the kernel a(s) is defined by
a(s) )
2π|V12|2 p x2λ0
e-FAmG(xm|xm,s)
(3.27a)
Using the Green function defined in eq 3.21, we have the kernel a(s) in real space S a(t) ) km,e
[
x
2∆(t) 1 exp βEm 2 1 + ∆(t) 1 - ∆ (t)
]
∫
(3.28)
a(s) ) km,e /s + Rm
(3.29)
Inserting eq 3.29 into eq 3.26, we have
km,eS 1 Qm(s) ) - 2 s s [1 + R + kS /s] m m,e
(3.30)
Equation 3.30 leads to a single exponent decay
Qm(t) ) exp(- kmSt)
(3.31)
with a dynamically corrected channel rate constant defined by S km,e 1 + Rm
(3.32)
where we call Rm the “dynamic correction factor”, defined by
∫
2πV122 p x4πkBTλ0
{x
(
km ∑ m
S
)
S km,e
∑ m 1+R
(3.35) m
and
S
S dt[a(t) - km,e ])
(3.34)
where a total rate constant is defined as
k )
Take the s f 0 limit, we have
Rm )
Qm(t) ) exp(- ∑ kmSt) ) exp(- kSt) ∏ m m
(3.27b)
dt e-st[a(t) - km,eS] ) a(s) - km,eS/s
kmS )
Q(t) )
S
Equation 3.26 together with eq 3.27 is an exact solution. It is seen that the kernel a(t) contains two factors: one is the thermal equilibrium rate constant and the other is a function of ∆(t). When solvent dynamics is fast, ∆(t) f 0, the kernel reduces to the thermal equilibrium rate constant. To discuss the solvent dynamic effect on the ET rate constant, we need to quantitatively scale the solvent dynamic effect. Equation 3.27b suggests that this can be achieved by separating the rate constant from its equilibrium value. The requirement for the dynamic correction constant is that it will go to zero in the limit when the solvent dynamics is fast, and it will get larger and larger when solvent dynamics get slower and slower. In general, the dynamic correction constant is influenced by the time correlation function of the reaction coordinate. We will first discuss Zhu and Rasaiah’s14c,15 method to define a dynamic correction constant. Since the decay of Q(t) is not a single exponential, we will have to use approximations to extract a meaningful rate constant. Making use of the fact that a(t) f S km,e asymptotically when solvent dynamic relaxation is fast, the solvent dynamic correction can be extracted by introducing the following definition
Rm(s) )
It is seen that the “dynamic correction factor” is mainly influenced by the time correlation function of the reaction coordinate. If we consider all of the channels, the total survival probability can be written as
Am e-F e-βEm
∫
) }
2βEm∆ 1 exp -1 2 1+∆ 1-∆
∞
dt 0+ (3.33)
km,eS )
2π|V12|2 p x4πkBTλ0
( )
exp -
λq (λq/pω)m exp[-βEm] pω m! (3.36)
To simplify further, we might use a statistically averaged quantity R to replace each Rm by
R ) 〈Rm〉 )
1 e ∑ m
-βEm
e-βE Rm ∑ m m
(3.37)
Equation 3.35 can be rewritten as
kS )
keS 1+R
(3.38a)
where the thermal equilibrium rate constant keS is given by Jortner’s formula eq 3.13. If the ligand vibration is asymmetric, we have
keA k ) 1+R A
(3.38b)
where keA is given by eq 3.14. At this point, we have achieved our goal to characterize the solvent dynamic effect into a single constant R. When the solvent dynamics is very fast, R f 0, and when the solvent dynamics is very slow, R . 1. Equation 3.38 provides an easy and practical way to calculate the ET rate constant where the solvent dynamic effect is scaled out by the dynamic correction factor. e. General Solution of Channel-Base Reaction-Diffusion Equation for Model B. For model B, the channel rate constant km(x) is given by eq 2.30, which is a Gaussian. The change of km(x) from a delta function to a Gaussian function presents a challenge in solving eq 3.19. As a matter of fact, no analytical solution can be reached without approximations. The Green function method has been discussed in detail for the solution of diffusion equations with km(x) given by a Gaussian.9,11b,14c,23 By introducing a decoupling approximation13-15,27 and following exactly the same procedure as that of ref 14c, we have been able to derive an analytical solution for the channel survival probability in Laplace space
km,eS 1 Qm(s) ) - 2 s s [1 + a(s)]
(3.39)
S where km,e is given by eq 3.16, and a(s) is given is given in real space by
3742 J. Phys. Chem. B, Vol. 112, No. 12, 2008 S a(t) ) km,e
[
x
2η∆(t) 1 exp βEm 2 2 1 + η∆(t) 1 - η ∆ (t) Em )
Zhu et al.
]
(∆G0 + λ + mpω)2 4λ
(3.40a)
(3.40b)
where the reorganization energy fraction η is defined by η ) {λo}/{λ}. Comparing eq 3.40a with eq 3.27b, we may introduce an effective reaction coordinate time correlation function
∆A(t) ) η∆(t)
(3.40c)
This will allow us to rewrite the kernel a(s) in real space as S a(t) ) km,e
[
x
2∆A(t) 1 exp βEm 2 1 + ∆A(t) 1 - ∆A (t)
]
(3.40d)
It is clear at this point that the solution of model B is exactly the same as that given in model A except the time correlation function of the reaction coordinate is rescaled! It is seen that η plays a role of energy gating the dynamics. When λ f λ0, model B reduces to model A; if η f 0, the solvent-dynamics effect is turned off. Following the same procedure as given in the previous section and introducing the dynamic correction factor
Rm )
S ]) ∫ dt[a(t) - km,e S km,e
∫0+∞ dt
{x
) }
(
2βEm∆A(t) 1 exp -1 2 1 + ∆A(t) 1 - ∆A
(3.41)
we have
k ) S
km ∑ m
S
)
S km,e
∑ m 1+R
(3.42) m
If we replace Rm by its statistical average 〈R〉 as defined in eq 3.37, we have the dynamically corrected rate constant
k ) S
∑ m
kSm
)
keS 1+R
(3.43)
where kSe is given by eq 3.18. For the asymmetric case, we simply replace the superscript S by A. Equation 3.43 for model B looks the same as eq 3.38a for model A except the gradients are changed in model B where kSe is calculated from eq 3.18 and R is calculated from the rescaled time-correlation function. In summary, by solving the channel-based diffusion reaction equation with a full quantum mechanically based rate constant k(x), a dynamically corrected rate constant can be derived. The effect of solvent dynamics on bio-molecular ET reaction is measured by the dynamic correction factor. In this way, the effect of solvent dynamics can be scaled quantitatively by a single constant. Though our general solution has focused on the symmetric case, the results also apply to the asymmetric case. IV. Numerical Results and Discussion From our discussion in previous sections, it is seen that the rate processes in condensed phase ET reactions are very complicated. Its analysis spans quantum mechanics and nonequilibrium statistical mechanics. We have provided a concrete formula to calculate the ET rate constant in such complicated
Figure 1. Plots of dynamic correction factor against average relaxation time of solvents, where T ) 298 K, ∆G0 ) -1500 cm-1, hhω ) 1500 cm-1, V12 ) 400 cm-1, λq ) 8000 cm-1 λ ) 2400 cm-1, λy ) 1000 cm-1, λo ) 1400 cm-1. The average relaxation time 〈τ〉 for different solvent is given in Table 1. Line a: Jortner-Bixon’s formulation (see Apendix 2) marked by triangle. Line b,c: model A, where the full line is calculated from eq 3.35 and eq 3.38 marked by square, and the dashed line marked by triangle is the statistically averaged R calculated from eq 3.37. Line d,e: model B, where the full line is calculated from eq 3.43 marked with diamond, and the dashed line marked with triangle is the statistically averaged R calculated from eq 3.37. Line f: Zhu and Rasaiah’s formulation (see ref 14c) where the inner-sphere coordinates are treated classically.
systems. The overall ET rate is given in the form of eq 1.1 where the thermal-equilibrium rate constant is corrected by a dynamic correction factor. The thermal-equilibrium rate constant assumes that the Boltzman distribution prevails for the reaction coordinate in the reactant potential well, such that the equilibrium distribution is not altered when the reactant is depleted at the transition state. This is similar to the assumption of transition-state theory. When the equilibrium distribution for the reaction coordinate does not exist, the quantum ET rate process experiences a frictional effect that slows down the ET process.26 This slowdown mainly comes from the solvent dynamics. In our theory, the solvent dynamics is reflected in a constant Rm called the dynamic correction factor, for each channel. It is more convenient to scale the solvent dynamic effect by a single factor for the whole reaction. There are two ways to achieve this goal. One way is to calculate Rm for each channel and then calculate the overall ET rate constant by using eq 3.35. After kS or kA is derived, we use eq 3.38 to solve for R. Another way is to use a statistically averaged R for Rm, where R is calculated from eq 3.37. We show in Figure 1 that there is little distinction between these two methods. However, it is the statistically averaged R that has allowed us to write the overall ET rate constant in the form of eq 1.1! To illustrate how big the dynamic correction factor is, we have collected different solvents in which the relaxation dynamics of the Born solvation energy has been measured experimentally. The dynamic parameters from experiment are given in Table 1(see also Tables 2 and 3). As we have shown elsewhere,14b the reaction coordinate time correlation function for ET reactions is the same as that of the Born solvation energy. The dynamic parameters given in Table 1 will be taken as for ∆(t) which can be expressed explicitly as
∆(t) ) a1 exp(-t/τ1) + a2 exp(- t/τ2)
(4.1)
Solvent Dynamics Effect
J. Phys. Chem. B, Vol. 112, No. 12, 2008 3743
TABLE 1: Solvent Dynamic Relaxation Parameters for Various Solvent solvent
A1
τ1, ps
A2
τ2, ps
〈τ〉, ps
aceronitrilea acetonea dimethylsulfoxidea acetoneb acetoneb propylenecarbonatea ethylacetateb propylenecarbonateb benezonitrileb triacetinb methanola
0.73 0.47 0.57 0.65 0.81 0.46 0.83 0.69 0.55 0.92 0.40
0.27 0.31 0.33 0.7 0.9 0.43 1.5 1.1 3.6 3.4 1.16
0.27 0.53 0.43 0.35 0.19 0.54 0.17 0.31 0.45 0.08 0.60
1.05 0.99 2.3 2.3 7.0 4.1 7.0 5.8 4.7 >30 9.57
0.48 0.67 1.2 1.26 2.06 2.4 2.435 2.557 4.1 5.53 6.2
a Data from ref 14c. b Data from ref 16 b. Experimental data was fit into eq 4.1 for 〈τ〉.
TABLE 2: Parameters Used in Calculating ET Rate Constant for Betaine-30 in Different Solventsa solvent
T, K
∆G0
hhω
λq
λo
λy
acetonitrile acetone methylacetate ethylacetate benzonitrile GTA
298 298 298 298 298 313 308 303
11645 10954 10465 10363 11196 10557 10583 10609
1554 1763 1697 1600 1631 1800 1800 1800
1276 924 1016 952 1172 1143 1200 1258
3344 3004 2812 2828 3226 2962 2881 3069
1223 1223 1223 1223 1223 1223 1223 1223
a
〈τ〉, ps
λs
1921 0.5 1781 0.83 1589 1.66 1605 2.63 2003 4.76 1739 18 1658 28 1846 40
Data from ref 13b.
TABLE 3: ET Rate Constant for Betaine-30 Calculated by Using Different Modelsa solvent
T, K
‚tO Å
kobs
k(A)
k(B)
k(B)
k(Jortner)
acetonitrile acetone methylacetate ethylacetate benzonitrile GTA
298 298 298 298 298 313 308 303
0.5 0.83 1.66 2.63 4.76 18 28 40
2.0 1.4 0.77 0.67 0.47 0.40 0.36 0.34
2.82 1.38 0.98 0.82 0.42 0.11 0.06 0.05
2.36 (1) 0.87 (1) 0.52 (1) 0.31 (1) 0.24 (1) 0.01 (1) 0.003(1) 0.002(1)
2.36 (1) 1.40 (2) 0.85 (2) 0.64 (2) 0.52 (2) 0.38 (5) 0.35 (5) 0.33 (5)
0.45 0.20 0.06 0.03 0.01 0.0004 0.0002 0.00006
a Parameters used come from Table 2. Note all the calculations are done for a symmetric reaction where only a single frequency is used for the fast coordinate. The unit of 〈τ〉 is ps, and the unit for all the rate constant k is 1/ps. kobs: experimental value. k(A): result from model A. k(B): result from model B where the number of the scaling index is given in brackets. k(Jortner): result from Jortner-Bixon given in Appendix 1.
The average solvent relaxation time can be defined as
〈τ〉 )
∫0∞ ∆(t) dt ) a1τ1 + a2τ2
(4.2)
To compare our results with others, we have extracted the dynamic correction factor from the Jortner-Bixon’s formula and Zusman’s procedure in Appendix 1 and Appendix 2. We have found that the dynamic correction factor from Zusman’s method diverges. In Figure 1, R is plotted versus the average solvent relaxation time 〈τ〉 for different models. It is seen that the general pattern is that R increases with the increase of 〈τ〉. Jortner and Bixon’s formula is linearly proportional to the solvent relaxation time, and leads to the largest R. The result from model A is very close to the Jortner-Bixon result except slightly curved. The result from model B is much less than model A but greater than the Zhu-Rasaiah result,14c which is derived from the solution of the classical Sumi-Marcus model.13 Figure 1 also shows that the statistically averaged R is remarkably accurate in replacing the Rm in eq 3.37, eq 3.38,
Figure 2. Relationship between matrix coupling element and correct factor in methonal for different theoretical models with T ) 298 K, ∆G0 ) -1500 cm-1, hhω ) 1500 cm-1, λq ) 8000 cm-1 λ ) 2400 cm-1, λy ) 1000 cm-1, λo ) 1400 cm-1, A1 ) 0.4, A2 ) 0.6, τ1 ) 1.16, τ2 ) 9.57. a: Jortner-Bixon’s formulation (see Apendix 1). b,c: model A, where the full line is calculated from eq 3.35 and eq 3.38, and the dashed line is the statistically averaged R calculated from eq 3.37. d,e: model B, where the full line is calculated from eq 3.43, and the dashed line is the statistically averaged R calculated from eq 3.37. f: Zhu-Rasaiah’s formulation (see reference 14c) where the innersphere coordinates are treated classically.
and eq 3.43. In our calculation of R by using eq 3.27, we have found that typically there are three to six terms that give significant value in the statistical average. The general trend is that, the smaller the energy quantization, the less accurate is the agreement between eq 3.42 and eq 3.43. In another limit when the energy quantization of vibration turns to a continuum spectrum, the GZSM model collapses to the classical SumiMarcus model, and the replacement of eq 3.43 to eq 3.42 will be exact. The influence of the matrix coupling element on the dynamic correction factor is given in Figure 2. It is seen that R is essentially a parabola, which is expected. It is also seen that R from Bixon-Jortner formula lies on top, followed by our model A and model B, and R from Zhu-Rasaiah’s formula lies on the bottom. If we recall the fact that Zhu-Rasaiah’s previous work14c is based on the Sumi-Marcus model13 where the innersphere ET rate process is treated by classical transition-state theory, and our present work is based on Fermi’s golden-rule theory, where the inner-sphere ET rate process is treated quantum mechanically, then Figure 1 and Figure 2 tell us that the solvent friction effect is stronger in the quantum ET process than in the classical ET process. The same trend has also been seen for the adiabatic rate process.26 To compare our theory with experiment, we use the experimental data from Barbara’s group16 since all of the parameters we need are available. With the parameters given in Barbara’s paper, our calculations show that the thermal equilibrium rate constant, either kS in model A or kS in model B, deviates substantially from experimental data. The overall rate constants calculated from model A and from the Jortner-Bixon formula are also far away from the experimental value. The only model that survives is model B, where the theoretical rate constant agrees well with the experiment for solvents that have an average relaxation time less than 5 ps, and breaks down for a solvent like GTA. If we examine the plot of R in Figure 1, the dynamic correction factor becomes infinitely large as the solvent dynamics slows down. This obviously is not true for a GTA solvent.
3744 J. Phys. Chem. B, Vol. 112, No. 12, 2008
Zhu et al.
Figure 3. Plot of the inverse of the ET rate constant against the average solvent relaxation time for different models in different solvents. Black squares on line e are the experimental values. Data is from Table 3. a: Jortner-Bixon’s formulation (see Apendix 2). b: model A. c: model B where the scaling index n ) 1. d: Result from mean first passage time theory where data from ref 16b. e: model B with scaling index adjusted for different solvent.
To correct this, we need to modify our theoretical formulation in model B. The difference between our theoretical formulation for model B and model A lies in the fact that a slow reaction coordinate is introduced for the ligand motion in model B. This leads to a balance in the reorganization energies and to a modification to the dynamic correction factor. The modification to R is by eq 3.40c, which introduces an effective time correlation function of the reaction coordinate. The reorganization energy fraction η serves as the energy gating to the reaction coordinate dynamics. If η f 0, the solvent dynamic effect is also switched off. From this analysis, we might modify our definition of the effective time correlation function of the reaction coordinate by
∆A(t) ) ηn∆(t)
(4.3)
where n is an ad hoc scaling index that reduces the contribution of η in the effective time correlation function. This scaling index renormalizes the time correlation function and serves as a correction to the approximations we have introduced in solving model B. This scaling index allows us to bridge the gap between theory and experiments. As a matter of fact, we can reach good agreement with experimental data by using the scaling index n. It is seen from Figure 3 that our theoretical value of the rate constant matches quite well with the experimental value for model B even in the very slowly relaxed solvent GTA after introducing n. By using the experimental data from Barbara’s group, the conclusion reached by us is the same as that reached by them; that is, in our language, model B is better than model A. The essential physics behind this conclusion is that model B considered the reorganization motion of the in-sphere nuclei by providing a wider reaction window for quantum transition. The fast inner-sphere ligand reorganization competes with the slow outer-sphere solvent reorganization in ET reactions. For a protein molecule with many internal modes, the inner-sphere contribution to the total reorganization energy cannot be neglected. This essential physics is captured by model B and ignored by model A. Our analytic solution for model B shows that the inner-sphere reorganization provides an energy balance in solvent reorganization with a modification to the time
correlation function of the reaction coordinate. Though model B has assumed that the motion of the inner-sphere ligand is fast compared with the outer-sphere solvent, there is no reason to prevent them from switching their roles with the inner-sphere ligand motion slow and the outer-sphere solvent motion fast. This role switching remains to be observed in real systems. It is worthwhile comparing our method with that of the Barbara’s group. First, our solution is analytic, while their solution is numeric; second, we have built a rate constant analytically by compressing the dynamic correction into a single constant R, while Barbara’s group has used the mean firstpassage time theory to construct the rate constant. Last, we have solved the channel-based diffusion reaction equations. Instead, Barbara’s group has collapsed all of the channels into one equation. Since our formulation is analytical, it is user-friendly to the experimentalist if experimental data needs to be tested against a theory. ET in Beaton-30 has also been studies with computer simulation by Rossky’s group.26 They have reported that intramolecular torsional dynamics plays an important role in their calculation of the ET rate constant. Since the formalism used by us is different from theirs, there is no common ground to compare our result with their result. If the intramolecular torsional dynamics can be measured and separated from the outsphere solvent dynamics, two reaction coordinates will have to be used, and GZSM-B will have to be extended to have two diffusive coordinates that are in nonequilibrium distribution.14d We will leave that for future study. In conclusion, we have studied two GZSM models, model A and model B. By solving the channel-based reaction-diffusion equations, our analytic solution has allowed us to build a rate constant given in eq 1.1 where the solvent dynamic effect is compressed into a single factor R which we call a “dynamic correction factor”. The dynamic correction factor measures the effect of solvent friction on the biological ET rate process. If solvent dynamics is fast, the dynamic correction can be neglected, and the rate constant reduces to the thermal equilibrium rate constant, which is the traditional Marcus-Jortner formulation. By comparing the two models with a set of experimental data, we have shown that model B is generally better than model A, and good agreement between theory and experiment can be reached for model B after a scaling index is introduced in the effective time correlation function of the reaction coordinate. Acknowledgment. We would like to thank Professor James T. Hynes and Professor Jayendran C. Rasaiah for reading the manuscript. J. Zhu would like to thank Professor Zhao Xinshen sincerely for being his host at Beijing University. The work of J. Zhu was supported in part by National Science Foundation NIRT Award No. 29963. The work of G.S. was supported by the Division of Chemical Science, Office of Basic Energy Science, Office of Energy Research and U.S. Department of Energy. Appendix 1 In this appendix, we define an alternative dynamic correction factor to that found in Jortner and Bixon’s formulation.7c This will allow us to compare our theory with theirs. The JortnerBixon’s rate constant can be written as7c
k)
kSm
∑ m 1+R
(A1) m,JB
Solvent Dynamics Effect
J. Phys. Chem. B, Vol. 112, No. 12, 2008 3745
where kmS is defined by eq 3.36, and Rm,JB is defined by
Rm,JB )
4πV122τL (λq/pω)m exp(- λq/pω) pλo m!
not plotted Zusman’s dynamic correction factor in Figure 1 simply because it blows up. References and Notes
(A2)
where τL is longitudinal relaxation time for a Debye solvent. Though a detailed derivation is not given in Jortner-Bixon’s paper, its theoretical foundation can be traced to path-integral theory and Landau-Zener theory with an assumption that the solvent is a Debye solvent characterized by a single relaxation time. Similar definitions to eq A2 can also be found in other places.9,10,14 The key ingredient introduced in eq A2 is the τL that characterizes the solvent dynamics, although it is oversimplified. To make Jortner-Bixon’s formula useful in its application to real solvents, one has to replace τL by an average 〈τ〉 where once again an assumption is implied that eqs A1 and A2 can be applied to non-Debye solvents. In Figure 1, we have plotted Jortner-Bixon’s dynamic correction factor by using eqs A1 and eq A2 with τL replaced by 〈τ〉. Appendix 2 In this appendix, we will use Zusman’s procedure9 to derive a dynamic correction factor. By interpolating the long-time and short-time approximation, the Green function eq 3.21 can be approximated as
G(xm | xm,s) ) Pmeq(xm)/s + 〈τ〉/| xm |
(B1)
where 〈τ〉 ) -d∆(t)/dt is the mean relaxation time at t f 0. Inserting eq B1 into eq 3.27 and eq 3.26, we have
k)
kSm
∑ m 1+R
(B2) m,ZU
which leads to a dynamic correction constant for channel m
Rm,ZU )
2πV122 p x4πkBTλ0
Ame-F
〈τ〉 | xm |
(B3)
It is easy to see that the Zusman’s dynamic correction constant will diverge for a particular channel m when xm f 0. We have
(1) Marcus, R. A.; Sutin, N. Biochem. Biophys. Acta 1985, 811, 265. (2) Lewis, F. D.; Letsinger, R. L.; Wasielewski, M. R. Acc. Chem. Res. 2001, 34, 159. (3) (a) Bixon, M.; Jortner, J. AdV. Chem. Phys. 1999, 106, 35. (b) Anto, K.; Hynes, J. T. AdV. Chem. Phys. 2000, 106, 35. (c) Cukier, R. I.; Nocera, D. G. Annu. ReV. Phys. Chem. 1988, 49, 337. (d) Devault, D. Quantum Mechanical Tunneling in Biological Systems; Cambridge University Press: 1984. (4) (a) Marcus, R. A. J. Chem. Phys. 1956, 24, 966. (b) Marcus, R. A. Frady Discussion. Chem. Soc. 1960, 29, 21. (5) Hush, N. S. J. Chem. Phys. 1958, 28, 962. (6) Levich V. G. Physical Chemistry An AdVanced Treatise; Erying, H, Henderson, D. W., Jost, W., Eds.; Academic: New York, 1970. (7) (a) Freed, K. F.; Jortner, J. J. Chem. Phys. 1970, 52, 6272. (b) Ulstrup, J.; Jortner, J. J. Chem. Phys. 1975, 63, 4358. (c) Jortner, J.; Bixon, M. J. Chem. Phys. 1988, 88, 167. (8) Lewis, F. D.; Kalgutkar, R. S.; Wu, Y.; Liu, X.; Liu, J.; Hayes, R. T.; Miller, S. E.; Wasielewski, M. R. J. Am. Chem. Soc. 2000, 122, 12346. (9) Rusman, L. D. Chem. Phys. 1980, 49, 295. (10) Friedman, H. L.; Newton, M. D. Farady Discuss. Chem. Soc. 1982, 74, 73. (11) (a) Hynes, J. T. J. Phys. Chem. 1986, 90, 3701. (b) Fonseca, T. J. Chem. Phys. 1989, 91, 2869. (12) (a) Cukier, R. I. J. Chem. Phys. 1988, 88, 5594. (b) Zhu, J.; Cukier, R. I. J. Chem. Phys. 1993, 98, 5679. (c) Cukier, R. I.; Zhu, J. J. Phys. Chem B 1997, 101, 7180. (13) Sumi, H.; Marcus, R. A. J. Chem. Phys. 1986, 84, 4894. (14) (a) Zhu, J.; Rasaiah, J. C. J. Chem. Phys. 1991, 95, 3325. (b) Zhu, J.; Rasaiah, J. C. J. Chem. Phys. 1992, 96, 1435. (c) Zhu, J.; Rasaiah, J. C. J. Chem. Phys. 1994, 101, 9966. (d) Zhu, J.; Ma, R.; Lu, Y.; Stell, G. J. Chem. Phys. 2005, 123, 224505. (15) Zhu, J.; Wang, J.; Stell, G. J. Chem. Phys. 2006, 125, 164511. (16) (a) Barbara, P. F.; Jarzeba, W. AdV. Photochem. 1990, 15, 1. (b) Walker, G. C.; Akesson, E.; Johnson, A. E.; Levinger, N. E.; Barbara, P. F. J. Phys. Chem. 1992, 96, 3728. (c) Johnson, A. E.; Levinger, N. E.; Jarzeba, W.; Schleif, R. E.; Kliner, D. A. V.; Barbara, P. F. J. Chem. Phys. 1993, 176, 555. (17) Gayathri, N.; Bagchi, B. J. Phys. Chem. 1996, 100, 3056. (18) Zhao, Y.; Zhu, W. J. Chem. Phys. 2007, 126, 184105. (19) Szabo, A.; Lamm, G.; Weiss, G. H. J. Stat. Phys. 1984, 34, 225. (20) Hubig, S. M.; Bockman, T. M.; Kochi, J. K. J. Am. Chem. Soc. 1996, 118, 3842. (21) Matyushov, D. V. J. Phys. Chem. B 2006, 110, 10095. (22) Dun, X.; Li, X.; He, R.; Chen, X. J. Chem. Phys. 2005, 122, 084314. (23) Kim, H.; Hwang, H.; Rossky, P. J. J. Phys. Chem. A 2006, 110, 11223. (24) Szabo, A.; Schulten, K.; Schulten, Z. J. Chem. Phys. 1980, 72, 4350. (25) Barbara, P. F.; Jarzeba, W. Acc. Chem. Res. 1988, 21, 195. More references are cited in 14b. (26) Zhu, J.; Cukier, R. I. J. Chem. Phys. 1995, 102, 4123. (27) Distribution function of thermal equilibrium state is related to a wave function of ground state of harmonic oscillator. The average of two non-equilibrium states can not be carried out analytically without introducing an approximation. The decoupling approximation used here is a mathematical tool to simplify the integral that involves the average of states other than the ground state. Its use has been justified in certain limits.13,14 It does not have a physical meaning.
