J. Phys. Chem. B 2008, 112, 577-584
577
Theory of Non-Markovian Rate Processes† Ji-Hyun Kim and Sangyoub Lee* Department of Chemistry, Seoul National UniVersity, Seoul 151-747, South Korea ReceiVed: June 30, 2007; In Final Form: August 24, 2007
Starting from a multidimensional reaction kinetic equation with a general time evolution operator and a reaction sink function K, we derive formally exact expressions for the survival probability, reaction time distribution, and mean reaction time by using the projection operator technique. These rate expressions are given in the rational function form, with the irreducible memory function Ωirr as the key ingredient. This approach has an advantage over the direct perturbation approaches that use the reaction term as the small parameter, in that Ωirr has a structure that can be perturbatively treated with (K - 〈K〉) as the small parameter. The well-known Wilemski-Fixman-type rate expressions are reproduced as the zeroth-order approximation from the present theory. Practical methods for evaluating the formal rate expressions are presented, and the results calculated for a model of electron transfer in non-Debye solvents are compared with computer simulations. It is found that the present approach is very promising for the study of non-Markovian dispersive kinetics.
1. Introduction Many rate processes occurring in condensed-phase systems are severely retarded by the diffusive transport. The usual way for treating such rate processes is to use the general reaction kinetic equation given as1-3 ∂ ψ(Γ,t) ) [D(Γ) - K(A(Γ))]ψ(Γ,t) ∂t
(1.1)
Here, ψ(Γ,t) is the probability density function that the reaction system has not undergone reaction by time t with the phase space coordinates given by Γ. D denotes the multidimensional Fokker-Planck operator,4 and K is the reaction sink function. A(Γ) represents a set of variables used to define the reaction zone, a region in the phase space at which the reaction can occur. As shown by Weiss,5 using the reaction term as the small parameter, a systematic perturbation analysis can be carried out to obtain a formal expression for the survival probability. However, such a perturbative approach fails when the reaction strength becomes very large, unless the perturbation series is suitably re-summed. The well-known Wilemski-Fixman (WF) rate expression is one of such re-summed expression, although it was originally derived in a less systematic way.1 In this work, we first derive formally exact expressions for the survival probability, reaction time distribution, and mean reaction time by using the projection operator technique. The survival probability expression so obtained has a rational function form that does behave better in the case with large reaction strength than the conventional perturbative expression. In fact, it has a similar structure as that of the WF expression. The key ingredient of the rate expressions is the irreducible memory function.6-8 With its lowest-order approximation, the rate expressions of the present theory reduce to the WF-type expressions. The validity of such WF rate expressions requires that (i) the underlying stochastic process should be a stationary
Markovian process and (ii) the reaction system has to be effectively one-dimensional and the reaction sink function has the delta-function form. It should be noted that some of the rate expressions given in this work coincide with those obtained by Yang and Cao using a different procedure.9 We then present practical methods for evaluating the formal rate expressions. Since our main interest is in the influence of the diffusive transport on the overall reaction kinetics, rather than the details of reaction zone dynamics, we will restrict the calculation to the case in which the sink function has the deltafunction form so that the reaction zone is actually a surface in the Γ space. When the system maintains equilibrium distribution up to the reaction surface, the survival probability decays exponentially. However, as depletion develops near the reaction surface, the survival probability curve becomes dispersive and cannot be characterized by the mean reaction time alone. For dealing with this difficulty, we propose an approach to approximate the irreducible memory function in terms of a dynamical gauging function that measures the non-Markovianity. As a numerical example, we consider the escape from a twodimensional quadratic potential where the reaction coordinate is a non-Markovian Gaussian variable. This model problem has been used in describing the nonadiabatic outer-sphere electrontransfer reaction in non-Debye solvents, whose rate is controlled by the sluggish solvent motion.10 We test the accuracy of the theoretical results by comparing with the computer simulation results. We find that the present theory produces quite satisfactory results that go beyond the WF-type approximations. 2. Formal Rate Expressions We consider an irreversible reaction process that can be described by the reaction kinetic equation given in eq 1.1. For further mathematical development, it is convenient to introduce the adjoint equation to eq 1.1
†
Part of the “James T. (Casey) Hynes Festschrift”. * To whom correspondence should be addressed. E-mail: sangyoub@ snu.ac.kr.
∂ φ(Γ,t) ) -[L (Γ) + K(A(Γ))]φ(Γ,t) ∂t
10.1021/jp075099b CCC: $40.75 © 2008 American Chemical Society Published on Web 10/13/2007
(2.1)
578 J. Phys. Chem. B, Vol. 112, No. 2, 2008
Kim and Lee
where φ(Γ,t) ) ψ(Γ,t)/ψeq(Γ), with ψeq(Γ) denoting the equilibrium distribution, which satisfies Dψeq(Γ) ) 0, and L () -D †) is an operator defined by the relation4 Dψeqφ ) -ψeq L φ
(2.2)
L is self-adjoint with respect to the inner product defined by 〈A|B〉 ) ∫dΓψeq(Γ)A(Γ)B*(Γ), where the asterisk indicates a complex conjugate, such that 〈A|L B〉 ) 〈L A|B〉.11 In the present notation, the survival probability can be expressed as S(t) ) 〈1|φ〉. A formally exact expression for its Laplace transform Sˆ (s) )
∫
∞
0
P ) |1〉〈1|
∂ Q |φ〉 ) -Q (L + K)(P + Q )|φ〉 ∂t
(2.5)
where k ) 〈1|K|1〉. A formal solution to eq 2.5 is + K)
Q|φ0〉 t
Lirr ≡ Q[(L + K) - (L + K)|1〉〈1|(L + K)|1〉-1〈1|(L + K)] ) Q[(L + K) - K|1〉k-1〈1|K]
-t′Q( L + K)
Q(L + K)|1〉S(t - t′) (2.6)
where φ0 ≡ φ(Γ,0). In this section, we will assume that ψ(Γ,0) ) ψeq(Γ), so that φ0 is unity. In the next section, we will consider the more general case when the distribution along A(Γ) is arbitrary. Substituting eq 2.6 into eq 2.4 and multiplying by 〈1| to the resulting equation from the left, we can obtain the expression for Sˆ (s) as
M ˆ (s) ) Ω ˆ irr(s)[1 + k-1Ω ˆ irr(s)]-1
M ˆ (s) ) 〈1|K[s + Q(L + K)] QK|1〉
(2.11)
where the Laplace-transformed irreducible memory function Ω ˆ irr(s) is given by irr -1
) Q K|1〉
(2.12)
Substituting eq 2.11 into eq (2.7) gives
[
Sˆ (s) ) s +
k 1 + kΩ ˆ (s)
]
-1
(2.13)
with Ω ˆ (s) ) k-2Ω ˆ irr(s) denoting the normalized irreducible memory function. Weiss obtained a Markovian version of eq 2.13,5 and recently, Yang and Cao obtained eq 2.13 for non-Markovian end-to-end reaction dynamics of polymers.9 Taking k as the small parameter, they first obtained a perturbation series expression for the survival probability Sˆ (s) and then re-summed the series to a rational function as given in eq 2.13 since it is the form suggested by the attractive Wilemski-Fixman result. To relate Ω ˆ (s) with dynamic quantities in the absence of reaction, (s + L irr)-1 in eq 2.12 is expanded as (s + L
irr -1
)
)
(s + L )-1
∞
∑ {-Q [K - K|1〉k
-1
〈1|K](s + L )-1}n (2.14)
n)0
where we used the relation P L ) 0. The Ω ˆ (s) can be then expressed as the infinite series in k ∞
Sˆ (s) ) 1/[s + k - M ˆ (s)] -1
(2.10)
With eq 2.9, the memory function M ˆ (s) can be rewritten as
(2.3)
∂ P |φ〉 ) -P ( L + K)(P + Q )|φ〉 ) -|1〉kS(t) - P K Q |φ〉 ∂t (2.4)
0
is the irreducible operator of L + K defined by
Ω ˆ irr(s) ) 〈1|K(s + L
Applying P and Q () 1 - P ) to eq 2.1 from the left and noting that 〈1| L ) 〈0|, we obtain
∫ dt′e
irr
dte-stS(t)
can be obtained from eq 2.1 by using the standard projection operator technique, with the projector defined by
Q|φ〉 ) e-t Q ( L
where L
Ω ˆ (s) )
[s + Q(L + K)]-1 ) (s + QL)-1 [s + Q(L + K)]-1QK(s + QL)-1 (2.8) is not safe because it can result in negative S(t) for large values of k. As will be seen below, the irreducible memory function approach introduced by Cichocki and Hess in the theory of colloidal suspension6 provides a better approximation scheme for the present problem also. The operator identity in eq 2.8 is modified as [s + Q(L + K)]-1 ) (s + Lirr)-1 [s + Q(L + K)]-1[Q(L + K)|1〉k-1〈1|K](s + Lirr)-1 (2.9)
n n
n
(2.15)
n)0
(2.7)
The memory function M ˆ (s) may be expanded as an infiniteorder continued fraction, but a truncation keeping only a finite number of elements often leads to erroneous behavior of S(t). Another common approach using the iterative applications of the operator identity
∑ (-1) k Cˆ (s)
where C ˆ 0(s) ) χˆ 1(s) and n
C ˆ n(s) ) χˆ n+1(s) -
∑ χˆ (s)Cˆ j
n-j(s)
j)1
χˆ n(s) ) k-(n+1)〈1|K[(s + L )-1Q K]n|1〉
(2.16)
It can be noted that both χˆ n(s) and C ˆ n(s) are independent of the magnitude of k. We can relate χˆ n(s) to the sink-sink correlation functions as n-1
knχˆ n(s) ) N ˆ n+1(s) - s-1
∑k
χˆ j(s)N ˆ n-j(s)
j+1
(2.17)
j)0
with the (n+1)th-order sink correlation functions defined by N ˆ n+1(s) ) k-1〈1|K[(s + L )-1K]n|1〉.
(2.18)
Non-Markovian Rate Processes
J. Phys. Chem. B, Vol. 112, No. 2, 2008 579
The Laplace transform expression for the reaction time distribution function, F(t)() -dS(t)/dt), can be obtained from eq 2.13 as
{
k 1 Fˆ (s) ) s 1 + k[Ω ˆ (s) + s-1]
}
(2.19)
∞
∑ (-1) Yˆ (s) n
(2.20)
n
n)0
where n
Yˆ n(s) ) N ˆ n+2(s) -
∑ Nˆ
ˆ n-j(s). j+1(s)Y
j)1
The physical meaning of the sink-sink correlation function defined in eq 2.18 can be more easily understood in terms of its time domain expression Nn+1(t) ) k-1
∫ dτ ∫
τn-1
t
n-1 0
0
× 〈1|δ(A(0) - a0)
[
dτn-2 ‚‚‚
∫
τ2
0
n-1
]
n
dτ1
∫ ∏ da K(a ) m
m
m)0
∏ δ(A(τ ) - a ) δ(A(t) - a )|1〉 (2.21) j
j
j)1
n
In obtaining this expression, we have used the following relations A(t) ≡ A(e-tL Γ) ) e-tLA(Γ)
(2.22)
δ(A(t) - a) ) e-tLδ(A(Γ) - a)
(2.23)
e-tL
∏ B (Γ) ) ∏ e n
-tL
[∏
]
n
With eqs 2.15, 2.16, and 2.17, the quantity k[Ω ˆ (s) + s-1] can be rewritten as k[Ω ˆ (s) + s-1] )
has to be Markovian for the (n+1)-time joint PDF to be factorized as
Bn(Γ)
(2.24)
In eq 2.21, the integrand is just the (n+1)-time joint probability density function (PDF), P(an,t;‚‚‚;a1,τ1;a0,0). We can thus rewrite eq 2.21 as
P(an,t;‚‚‚;a1,τ1;a0,0) )
G(am,τm|am-1,τm-1) P(a0,0) (2.27)
m)1
where τ0 ) 0, τn ) t, and G(am,τm|am-1,τm-1) denote the conditional probability that A(τm) ) am, given that A(τm-1) ) am-1. For stationary processes, G(am,τm|am-1,τm-1) ) G(am,τm - τm-1|am-1,0), and the square-bracketed factor on the righthand side of eq 2.25 is factorized in the Laplace domain. Second, the reaction system has to be effectively one-dimensional, and the sink function K must have the form, κδ(A(Γ) - σ), where κ denotes the intrinsic sink strength on the reaction surface, A(Γ) ) σ. When the reaction zone has a finite width or when the reaction surface is anisotropic, this second requirement for the validity of the WF approximation cannot be satisfied, and much attention has been paid to such cases.1,12 In the present work, we address the failure of WF approximation due to nonMarkovianity and propose a practical approach to improve on the WF result. 3. Initial State Dependence of the Kinetics In this section, we restrict the theoretical development to model systems in which a global reaction coordinate A(Γ) can be defined and the sink function depends only on A(Γ). We then assume that the system is prepared at t ) 0 such that A(Γ) ) a0 but is equilibrated along the coordinates transverse to A(Γ). Then, in eq 2.1, φ(Γ,0) ) δ(A(Γ) - a0)/Peq(a0), with the reduced equilibrium distribution along A(Γ) given by Peq(a0) ) 〈1|δ(A(Γ) - a0)|1〉. It can thus be shown that the reaction time distribution is given by ∂ F(t|a0) ) - S(t|a0) ) ∂t
∫ daK(a)G (a,t|a )
where G (a,t|a0) denotes the conditional PDF that A(Γ(t)) ) a, given that A(Γ(0)) ) a0. Its Laplace transform expression is given by Gˆ (a,s|a0) ) Peq(a0)-1〈1|δ(A(Γ) - a0)(s + L + K)-1δ(A(Γ) - a)|1〉 (3.2) Using the operator identity
n
Nn+1(t) ) k-1
∫ ∏ da K(a ) m
[∫ dτ ∫ 0
m)0
τn-1
t
n-1
0
(s + L + K)-1 ) (s + L )-1
m
dτn-2 ‚‚‚
∫
τ2
0
]
dτ1P(an,t;‚‚‚;a1,τ1;a0,0) (2.25)
[
]
k 1 + kχˆ 1(s)
∞
∑ [-K(s + L )
-1 n
]
n)0
This shows clearly how the survival probability is related to the non-Markovian joint PDFs, which describe the probabilities of repeated visits to the reaction sink zone in the absence of reaction. The WF decoupling approximation is to assume that N ˆ n+1(s) )N ˆ 2(s)n or, in terms of χˆ n(s), χˆ n(s) ) χˆ 1(s)n. This approximation leads to Sˆ WF(s) ) s +
(3.1)
0
-1
(2.26)
We can now clarify the conditions for the validity of the WF approximation. First, the underlying stochastic process of A(t)
Fˆ (s|a0) can be expressed as Fˆ (s|a0) ) Peq(a0)-1
∞
∑ (-1) 〈1|δ(A(Γ) - a )[(s + L ) n
0
-1
K]n+1|1〉 (3.3)
n)0
Averaging eq 3.3 over Peq(a0) gives Fˆ (s) ) s-1k
∞
∑ (-1) Nˆ n
n+1(s)
(3.4)
n)0
This expression can also be obtained by expanding eq 2.19 with the aid of eq 2.20.
580 J. Phys. Chem. B, Vol. 112, No. 2, 2008
Kim and Lee
Multiplying and dividing the right-hand side of eq 3.3 by the same factor, 〈1|δ(A(Γ) - a0)(s + L )-1K|1〉, we obtain
which has no dependence on a0. Equation 3.13 leads to the WFtype expression for the reaction time distribution13
∞
Fˆ (s|a0) )
[∫ daK(a)Gˆ (a,s|a )] ∑ (-1) Nˆ n
0
Fˆ WF(s|a0) )
0 n+1(s|a0)
n)0
G ˆ (a,s|a0) ) Peq(a0)-1〈1|δ(A(Γ) - a0)(s + L )-1δ(A(Γ) - a)|1〉 (3.6)
Comparison of eq 3.11 with eq 3.14 shows that the memory function Ψ ˆ (s|a0) carries the initial state information due to nonMarkovianity. In the limit of κ f ∞, the reaction time distribution becomes the first passage time distribution (FPTD). From eq 3.11, we obtain
0 and N ˆ n+1 (s|a0) is defined as
〈1|δ(A(Γ) - a0)[(s + L )-1K]n+1|1〉 〈1|δ(A(Γ) - a0)(s + L )-1K|1〉
Fˆ FPTD(s|a0) )
∫ daK(a)Gˆ (a,s|a ) 0
1+Ψ ˆ (s|a0)
ˆ (s|a0) Ψ ˆ ∞(s|a0) ≡ lim κ-1Ψ κf∞
A practical method for evaluating Ψ ˆ ∞(s|a0) is to take it as a [M/M]-Pade´ approximant. For example, the [0/0]-Pade´ approximant of Ψ ˆ ∞(s|a0) yields
(3.8) Fˆ FPTD(s|a0) =
∞
∑ (-1) Yˆ (s|a ) n 0 n
κ-1N ˆ 02(s|a0)
) L
{∫
0 0 ˆ n-j (s|a0) j+1(s|a0)Y
(3.9)
K(A(Γ)) ) κδ(A(Γ) - σ)
κG ˆ (σ,s|a0) 1+Ψ ˆ (s|a0)
κnL
{∫ dτ ∫
τn
t
0
(3.11)
in eq 3.7 explicitly given by
0 N ˆ n+1 (s|a0) )
n
0
dτn-1 ‚‚‚
∫
τ2
0
}
dτ1G(σ,t;σ,τn;‚‚‚;σ,τ1|a0);s
G ˆ (σ,s|a0) (3.12) Here, L{‚‚‚;s} denotes the Laplace transformation; that is L{‚‚‚;s} ≡
∫
∞
0
dte-st(‚‚‚)
G(σ,t;σ,τn;‚‚‚;σ,τ1|a0) is the joint probability of observing the system at σ at times τ1, ‚‚‚, τn and t, provided that it started from a0 at time 0. For stationary Markovian processes, eq 3.12 reduces to 0 (s|a0) ) [κG ˆ (σ,s|σ)]n N ˆ n+1
Fˆ FPTD WF (s|a0) )
(3.10)
where κ denotes the intrinsic sink strength on the reaction surface, A(Γ) ) σ. Then, eq 3.8 reduces to
with
0
dτG(σ,t;σ,τ|a0);s
}
(3.16)
(3.13)
(3.17)
With this additional approximation, eq 3.16 yields the WF-type expression for the FPTD
From now on, we will assume the delta-function reaction sink
0 N ˆ n+1 (s|a0)
t
G(σ,t;σ,τ|a0) ) G(σ,t - τ|σ)G(σ,τ|a0)
j)1
Fˆ (s|a0) )
[G ˆ (σ,s|a0)]2
For stationary Markovian systems, we have
n
∑ Nˆ
G ˆ (σ,s|a0)
0
n)0
0 Yˆ 0n(s|a0) ) N ˆ n+2 (s|a0) -
(3.15)
Ψ ˆ ∞(s|a0)
where
where Ψ ˆ (s|a0) is given by Ψ ˆ (s|a0) )
G ˆ (σ,s|a0)
(3.7)
With the two equivalent expressions for Fˆ (s) in eqs 2.19 and 3.4, we can see that the expression for Fˆ (s|a0) in eq 3.5 can be rewritten as
Fˆ (s|a0) )
(3.14)
1 + κG ˆ (σ,s|σ)
(3.5)
Here, G ˆ (a,s|a0) is the reaction-free Green’s function
0 N ˆ n+1 (s|a0) )
κG ˆ (σ,s|a0)
G ˆ (σ,s|a0)
(3.18)
G ˆ (σ,s|σ)
Recently, Sokolov showed that the FPTD in non-Markovian systems can be calculated from the integral equation14 G(σ,t|a0) )
∫ dτF t
FPTD
0
(τ|a0)G(σ,t|σ,τ;a0) (a0 * σ)
(3.19)
where G(σ,t|σ,τ;a0) is the conditional probability of observing the system at σ at time t, provided that it visited σ earlier at time τ, starting from a0 at time 0; it is different from G(σ,t;σ,τ|a0). It is worth comparing his result with ours in eq 3.15. First, the time integral on the right-hand side of eq 3.19 cannot be deconvoluted analytically, in general. In contrast, eq 3.15 gives the FPTD explicitly but at the cost of having to evaluate the memory function Ψ ˆ ∞(s|a0) given in an infinite series form. Second, the WF-type approximation for the FPTD in eq 3.18 can also be obtained from eq 3.19 by assuming the stationary nature and Markovianity, that is G(σ,t|σ,τ;a0) = G(σ,t - τ|σ). The mean first passage time (MFPT), which is the first moment of the FPTD, can be obtained more easily from the survival probability expression rather than from eq 3.15. The initial-state-dependent survival probability Sˆ (s|a0) () s-1[1 Fˆ (s|a0)])is given from eq 3.11 as
{
Sˆ (s|a0) ) s +
κ[sG ˆ (σ,s|a0)]
}
1 + [Ψ ˆ (s|a0) - κG ˆ (σ,s|a0)]
-1
(3.20)
Non-Markovian Rate Processes
J. Phys. Chem. B, Vol. 112, No. 2, 2008 581
Taking the limits, κ f ∞ and s f 0 and noting that
that a higher-order Pade´ approximation should give a better result, but it requires rapidly increasing computing effort. The mean reaction time τrx is related to the survival probability as
ˆ (σ,s|a0) ) Peq(σ) lim sG sf0
we obtain the MFPT Ψ ˆ ∞(0|a0) -G ˆ (σ,0|a0)
τ(a0) ) lim Sˆ (0|a0) )
Peq(σ)
κf∞
ˆ (0) τrx ) Sˆ (0) ) k-1 + Ω (3.21) and its [1/1]-Pade´ approximant is
The leading divergent terms of both Ψ ˆ ∞(s|a0) and G ˆ (σ,s|a0) in the s f 0 limit are s-1Peq(σ) and thus cancelled out in the numerator on the right-hand side of eq 3.21. By averaging τ(a0) over Peq(a0), we obtain τ)
∫
∞
0
dt
[
Ψ∞(t)
Peq(σ)
-1
]
(3.22)
where Ψ∞(t) denotes Ψ∞(t|a0) averaged over Peq(a0). Comparison of eqs 3.15 and 3.18 indicates that the WF-type approximation for τ is given by τWF )
∫
∞
0
dt
[
G(σ,t|σ) -1 Peq(σ)
]
χˆ 1(0)(χˆ 2(0) - χˆ 1(0)2) + (χˆ 3(0)χˆ 1(0) - χˆ 2(0)2)k 1 (4.4) + k χˆ (0) - χˆ (0)2 + (χˆ (0) - 2χˆ (0)χˆ (0) + χˆ (0)3)k 2 1 3 1 2 1 This result is in contrast with the Markovian case. As shown clearly by Bicout and Szabo, the reaction and transport contributions to τrx are decoupled in the Markovian case.16 That is, τrx ) k-1 + χˆ 1(0), where χˆ 1(0) is explicitly given by eq 3.23 and is independent of k. However, in the non-Markovian case, the reaction and transport contributions are intermingled in Ω ˆ (0) due to memory effects. The MFPT τ ) lim τrx kf∞
4. Methods for Evaluating the MFPT and Survival Probability
Sˆ [1/1](s) ) χˆ 1(s)(χˆ 2(s) - χˆ 1(s)2) + (χˆ 3(s)χˆ 1(s) - χˆ 2(s)2)k 1 + k χˆ (s) - χˆ (s)2 + (χˆ (s) - 2χˆ (s)χˆ (s) + χˆ (s)3)k 2 1 3 1 2 1
is given by τ = τ[1/1] )
In this section, we consider practical methods for evaluating the formal rate expressions given in the previous sections that go beyond the WF-type approximations. Wilemski and Fixman proposed an iterative scheme to obtain higher-order expressions beyond the WF closure, but their iterative scheme does not satisfy the required condition, that is, the normalization of the reaction time distribution, Fˆ (0) ) 1.1 Weiss suggested a more systematic perturbative scheme, using k as the small parameter.5 However, such an approach does not work in the limit of k f ∞ and cannot be applied to the first passage problems. We will confine ourselves to the cases in which the system is initially prepared in equilibrium before the onset of reaction and the reaction sink function has the delta-function form in eq 3.10. As a first step beyond the WF approximation, we take the [1/1]-Pade´ approximation to Ω ˆ (s) in eq 2.15. Putting the resulting expression into eq 2.13, we obtain the Laplacetransformed survival probability as
{ [
τ[1/1] rx )
(3.23)
This expression for τWF was given by Pastor, Zwanzig, and Szabo.15
s+
(4.3)
]}
-1 -1
χˆ 3(0)χˆ 1(0) - χˆ 2(0)2 χˆ 3(0) - 2χˆ 1(0)χˆ 2(0) + χˆ 1(0)3
(4.5)
When the system maintains an equilibrium distribution up to the boundary of the reaction zone, the survival probability is exponential, and the MFPT is sufficient to characterize it. However, as depletion develops near the reaction zone, the survival probability displays a nonexponential behavior. The S(t) curve becomes more dispersive; that is, the nth-order moment of FPTD gets larger than the nth power of MFPT. Hence, we need the FPTD or survival probability itself. The survival probability expression in eq 4.1 is more difficult to evaluate than the MFPT. We therefore propose an efficient, though less rigorous, alternative approach to calculate the survival probability. The key dynamic quantity that determines the time dependence of the survival probability is the irreducible memory function Ω ˆ irr(s) defined in eq 2.12. Comparison of the exact Sˆ (s) in eq 2.13 and its WF approximation Sˆ WF(s) in eq 2.26 shows that Ω ˆ (s)[) k-2Ω ˆ irr(s)] is taken as χˆ 1(s) in the WF approximation. Indeed, the two functions have the common structure, which are reproduced below for easy reference Ω ˆ (s) ) k-2〈1|K(s + Lirr)-1QK|1〉
(4.6)
From eq 2.17, explicit expressions for χˆ n(s) are given in terms of the sink-sink correlation functions N ˆ n(s) in eq 2.25 as
χˆ 1(s) ) k-2〈1|K(s + L)-1QK|1〉
(4.7)
χˆ 1(s) ) k-1N ˆ 2(s) - s-1
The two functions behave identically at short times. In the limit s f ∞, we have
(4.1)
χˆ 2(s) ) k-2N ˆ 3(s) - 2s-1k-1N ˆ 2(s) + s-2 -3
-1 -2
-1 -2
χˆ 3(s) ) k N ˆ 4(s) - 2s k N ˆ 3(s) - s k N ˆ 2(s) + ˆ 2(s) - s-3 (4.2) 3s-2k-1N 2
Equation 4.1 shows that the WF approximation in eq 2.26 corresponds to using the [0/0]-Pade´ approximant. It is expected
Ω ˆ (s) = χˆ 1(s) ∼ s-1k-2〈1|KQK|1〉
(4.8)
This limiting behavior complies with the known behavior of the survival probability that Sˆ (s) ∼ [s + k]-1 for s f ∞ since Ω ˆ (s) ∝ s-1 , 1 in eq 2.13. It also tells that SWF(t) behaves correctly at short times, regardless of the magnitude of k.
582 J. Phys. Chem. B, Vol. 112, No. 2, 2008
Kim and Lee
Hence, it is worthy of examining the behavior of a dynamical gauging function defined by f(t) ≡ Ω(t)/χ1(t)
(4.9)
As just shown, f(t) becomes unity at short times. In the Markovian reaction system involving a highly localized sink, f(t) is unity since the WF approximation becomes exact in such cases. Hence, one can say that f(t) gauges the non-Markovian effects. From the variational theories of Doi17 and Portman and Wolynes,18,19 it is known that Sˆ (s) e Sˆ WF(s)
(4.10)
The mean force potential along the coordinate z can be determined from the relation Peq(z) ) 〈1|δ(x + y - z)|1〉 and has the quadratic form βV(z) ) z2/2〈z2〉, with 〈z2〉 ) 〈x2〉 + 〈y2〉. The reaction coordinate z is the solvent coordinate appearing in the usual electron-transfer theory. We thus have 〈z2〉 ) 2λkBT, with λ denoting the solvent reorganization energy and zc ) λ + ∆G, with ∆G denoting the total free energy change of the reaction. The inherent rate parameter κ is given by 2π|HRP|2/p, with HRP denoting the electronic coupling matrix element between reactant and product states. The time scale of the solvent fluctuations is characterized by the time correlation function φ(t)() 〈z|z(t)〉/〈z2〉). For non-Debye solvents, for example, alcohols and other hydrogen-bonding solvents, φ(t) is often assumed to be given by a biexponential form22
With eqs 2.13 and 2.26, this inequality implies that Ω(t) e χ1(t)
φ(t) ) cxe-t/τx + cye-t/τy (4.11)
Hence, for non-Markovian cases, f(t) decreases from unity as t gets larger. It seems natural to assume that f(t) approaches convexly to unity as t f 0, such that f(t) ∼ 1 - Rt˜β (R > 0, β > 1) for small dimensionless time ˜t. It is known that a reasonable approximation to the irreducible memory function often produces a quite good result.20 We will assume that the dynamical gauging function has the form f(t) )
1 (1 + ˜t β)R
(4.12)
Since f(t) behaves as ˜t -Rβ for large ˜t, there is a redundancy in determining the power-law exponent. Furthermore, we only have a single piece of information to determine the exponent. Therefore, we first set β to 2 and obtain the value of R from the equation τ(= τ[1/1]) ) lim Ω ˆ (0) ) kf∞
∫
∞
0
dtχ1(t)f(t)
(4.13)
χ1(t) can be obtained analytically for some tractable systems and may be evaluated very easily from computer simulations,21 as it is related to the probability of finding the system at the reaction zone at time t, provided that it was initially at the reaction zone. We will check the utility of eq 4.12 in the next section for a simple reaction model.
with cx + cy ) 1. Since 〈z2〉 ) 〈x2〉 + 〈y2〉 and 〈z|z(t)〉 ) 〈x|x(t)〉 + 〈y|y(t)〉, we have cx ) 〈x2〉/〈z2〉, cy ) 〈y2〉/〈z2〉, Dx ) 〈x2〉/τx, and Dy ) 〈y2〉/τy. The overall reaction time is given by eq 4.3 τrx ) k-1 + τR
(5.5)
Here, k ) 〈κδ(z - zc)〉 ) κPeq(zc) is the rate constant that would be observed if the solvent relaxed infinitely fast, and its expression was first obtained by Marcus;23 τR is the mean time required for reaching zc starting from equilibrium in the reactant potential well. In this work, we are interested in the case when the overall reaction rate is controlled by the activation rate in the reactant potential well. When k-1 , τR, τR can be equated to the MFPT τ, whose approximate expression is given by eq 4.5. Since the stochastic process under consideration is Gaussian, the sink-sink correlation functions in eq 2.25 can be expressed as k-nNn+1(t) )
∫ dτ ∫ ∫ dτ M
τn-1 n-1 0
t
0
τ2
0
1
dτn-2 ‚‚‚
n+1(t
- τn-1,τn-1 - τn-2, ...,τ2 - τ1,τ1)
Mn+1(t - τn-1,τn-1 - τn-2, ...,τ2 - τ1,τ1) ) P(zc,t;zc,τn-1;‚‚‚;zc,τ1;zc,0)
)
Peq(zc)n+1
5. Numerical Example To test the accuracy of the present theory, we consider a simple model problem, namely, the escape from a twodimensional quadratic potential. This model was used by Bicout and Szabo to describe the nonadiabatic electron-transfer dynamics occurring in a “biexponential” non-Debye solvent.10 The reaction kinetic equation we have to consider has the same form as eq 1.1, with Γ ) (x,y) and A(Γ) ) z(x,y) ) x + y. The evolution operator D(Γ) and the reaction sink function K(A(Γ)) are explicitly given by ∂ ∂ ∂ ∂ D(x,y) ) Dx e-βV eβV + Dy e-βV eβV ∂x ∂x ∂y ∂y
(5.1)
K(z) ) κδ(z - zc)
(5.2)
Here, β ) 1/kBT and βV(x,y) )
(5.4)
x2 y2 + 2 2〈x 〉 2〈y2〉
(5.3)
n+1
exp[(n + 1)γ] exp[-γ
∑ (A
-1 n+1
)ij]
i,j)1
(5.6)
xdet(An+1)
Here, γ ≡ V(zc)/kBT ) z2c /2〈z2〉 and the (n + 1) × (n + 1) normalized symmetric covariance matrix An+1 is given by
(
A n+1 ) φ(t - τn-3) 1 φ(t - τn-1) φ(t - τn-2) φ(τn-1 - τn-2) φ(τn-1 - τn-3) 1 φ(τn-2 - τn-3) 1 1
... ... ... ....
φ(t) φ(τn-1) φ(τn-2) φ(τn-3)
.
l 1
..
)
(5.7)
where only the upper triangular part is displayed for simplicity.
Non-Markovian Rate Processes
J. Phys. Chem. B, Vol. 112, No. 2, 2008 583
Figure 1. Comparison of the approximate theoretical MFPTs with the BD result. The open circles and the black solid circles represent the results of [0/0]- and [1/1]-Pade´ approximants to the exact MFPT, respectively; γ ) V(zc)/kBT, and cx is fixed to 0.5.
Then, {χˆ n(0);n ) 1,2,3} appearing in eq 4.5 can be calculated from eq 4.2 as
∫
χˆ 1(0) ) χˆ 2(0) ) χˆ 3(0) )
∫
∞
0
∫
∞
0
dt1
dt1
∫
∞
0
∫
∞
0
dt2
∞
0
dt1[M2(t1) - 1]
dt2[M3(t2,t1) - M2(t2) - M2(t1) + 1]
∫
∞
0
dt3[M4(t3,t2,t1) - M3(t3,t2) -
M3(t2,t1) - M2(t3)M2(t1) + M2(t3) + M2(t2) + M2(t1) - 1] (5.8)
with t1 ) τ1, t2 ) τ2 - τ1, and t3 ) t - τ2. To obtain the expression for χˆ 3(0), for example, we have changed the order of integrations as
∫ dt ∫ dτ ∫ ∞
0
t
0
2
τ2
0
dτ1 )
∫ dt ∫ dτ ∫ dτ ∞
∫
∞
∫
∞
∫ )∫
∞
)
0
)
0
)
0 ∞
0
t
t
0
1
0
dτ1
∫
∞
dτ1
∫
∞
∫ ∫
∞
dτ1 dt1
τ1
τ1
0 ∞
0
dt
2
τ1
∫
t
τ1
dτ2
∫
dτ2 ∞
τ2
dt
d(τ2 - τ1) dt2
∫
∞
0
∫
∞
0
d(t - τ2)
dt3.
We compare the approximate MFPT calculated from eqs 4.5, 5.6, 5.7, and 5.8 with the results obtained from Brownian dynamics (BD) simulations. We use the simulation algorithm that is detailed in the Appendix B of ref 10, and for each set of parameters, 100 000 BD trajectories are generated to calculate the statistical average. Figure 1 displays the ratios of the approximate theoretical MFPTs to the BD results. The variable for abscissa, η() τx/τy), is varied from 10 to 102, and cx is fixed to 0.5. The open circles represent the WF results, corresponding to [0/0]-Pade´ approximants to exact τ. The black solid circles represent the results of [1/1]-Pade´ approximants. One can see that τ[1/1] provides significant improvements on τ[0/0]() τWF) for all parameter values and that τ e τ[1/1] e τWF. With increasing non-Markovian effects (increasing η values for fixed cx), the deviation of τ[1/1]/ τ from unity gets larger. For larger γ, it is more difficult to
Figure 2. Comparison of the approximate theoretical MFPTs with the BD result. The open circles and the black solid circles represent the results of [0/0]- and [1/1]-Pade´ approximants to the exact MFPT, respectively; γ() V(zc)/kBT) is fixed to 4.
reach zc repeatedly. This means roughly that the contribution of higher-order terms in the series expression for τ becomes less important so that the low-order Pade´ approximation becomes accurate. However, the barrierless case with γ ) 0 is an exception. There may be some cancellation of errors due to symmetry, which renders τ[1/1] to be a good approximation. Figure 2 also displays the ratios of the approximate theoretical MFPTs to the BD results, but the variable for abscissa is now cx, which is another parameter controlling the non-Markovianity. The value of γ is fixed to 4. When cx gets closer to either zero or unity, φ(t) in eq 5.4 becomes almost a single exponential form, and thus, the system is in the near-Markovian limit. When cx has an intermediate value, the non-Markovian effects are pronounced. Hence, with fixed η, the τ[1/1]/τ versus cx curve has an approximately parabolic form. Let us then consider the survival probabilities when the reaction rate is controlled by the sluggish solvent dynamics (i.e., in the limit of k-1/τ f 0). Figures 3and 4 display the survival probabilities obtained from BD simulations (represented by open circles) and various level of theories. The cx is fixed to 0.5, and values of γ and η are varied, as shown in the legends. The relaxation time of the fast mode, τy, is used as the time unit. It is very striking that the method employing the approximate dynamical gauging function f(t) in eq 4.12 produces excellent results. When the value of R is determined from eq 4.13 using the exact τ calculated from the BD simulation, the method produces almost exact survival probability curves (represented by the solid curves), except when γ ) 0. This indicates that f(t) in the form of eq 4.12 is a good approximation to exact f(t). When the exact value of τ is not available, we have to use the approximate τ[1/1] to determine the parameter R. However, even with the approximate τ[1/1], the dynamical gauging function method produces the survival probability curves (represented by the dot-dashed curves) that are in far better agreement with the BD simulation curves than the WF approximation (represented by the long-dashed curves) or simple exponential curve, S(t) ) exp(-t/τ(BD)) (represented by the short-dashed curves), in all cases. When η ) 100, τ[1/1] is a poor approximation to the exact τ so that the S(t) curves generated by the dynamical gauging function method with τ[1/1] deviate considerably from the BD
584 J. Phys. Chem. B, Vol. 112, No. 2, 2008
Kim and Lee 6. Concluding remarks
Figure 3. Survival probabilities obtained from BD simulations and various levels of theories for (γ,η) values indicated with cx ) 0.5. The open circles represent the BD results. The solid and the dot-dashed curves represent the results obtained by employing the dynamical gauging function in eq 4.12, with R calculated using τ(BD) and τ[1/1], respectively; in (a) and (c), with η ) 10, τ(BD) = τ[1/1] so that the two curves are hardly distinguishable. The long-dashed curves are the WF results, and the short-dashed curves are the simple exponential curves, S(t) ) exp(-t/τ(BD)).
We have presented an irreducible memory function approach to non-Markovian rate processes. Exact, though formal, expressions have been given for the survival probability, reaction time distribution, and mean reaction time. For cases with an initial equilibrium distribution, they are eqs 2.13, 2.19, and 4.3, respectively, while for cases where the system is initially located at a particular location on the reaction coordinate and the sink function has the delta-function form, they are given by eqs 3.20, 3.11, and 3.21. The key ingredients in those expressions are the irreducible memory functions given by eqs 2.15 and 3.9 in the respective cases. With the lowest-order approximations for the irreducible memory functions, the WF-type rate expressions are reproduced. Practical methods for evaluating the mean reaction time and the survival probability have been given. The mean reaction time is expressed as the whole time integral of the irreducible memory function given as a series in the “equilibrium” rate constant k. It is found that the [1/1]-Pade´ approximant to the series gives quite satisfactory results for the model problem considered. When the survival probability is nonexponential due to non-Markovian dispersive kinetics, approximation of the irreducible memory function in terms of the dynamical gauging function is found to give impressively good results. Although we have only considered a model problem, in which the reaction coordinate is a Gaussian stochastic variable with a biexponential time correlation function, the present theory can be easily applied also when the time correlation of the stochastic variable is highly nonexponential (e.g., power-law decaying or stretched exponential) unless the MFPT diverges. Investigation of the general applicability of the present theory to various nonMarkovian systems is in progress. Acknowledgment. This work was supported by the Korea Research Foundation Grant funded by the Korean Government (MOEHRD) (KRF-2005-070-C00065). References and Notes
Figure 4. Survival probabilities obtained from BD simulations and various levels of theories for (γ,η) values indicated with cx ) 0.5. See the Figure 3 caption for detailed explanation.
simulation curves. However, when η ) 10, the S(t) curves generated by the dynamical gauging function method as well as the values of τ[1/1] are in very good agreement with the exact ones. Considering that the actual values of η are not so large for most non-Debye solvents (for example, the value of η for n-propyl alcohol at 20 °C is 7.75),22 the dynamical gauging function method with τ[1/1] is expected to be very useful. Moreover, if cx has a value larger or smaller than 0.5, the theoretical results should be in better agreement with the exact ones, as implied by Figure 2. In the barrierless case (γ ) 0), even with the exact τ (indeed, τ[1/1] itself does not deviate much from the exact τ when γ ) 0), the S(t) curves generated by the dynamical gauging function method deviate from the BD simulation results. It seems that a better approximation for the dynamic gauging function other than that given by eq 4.12 is needed.
(1) (a) Wilemski, G.; Fixman, M. J. Chem. Phys. 1974, 60, 866. (b) Wilemski, G.; Fixman, M. J. Chem. Phys. 1974, 60, 866. (2) Agmon, N.; Hopfield, J. J. J. Chem. Phys. 1983, 78, 6947. (3) Bicout, D. J.; Szabo, A. J. Chem. Phys. 1998, 108, 5491. (4) Zwanzig, R. J. Chem. Phys. 1974, 60, 2717. (5) Weiss, G. H. J. Chem. Phys. 1984, 80, 2880. (6) Cichocki, B.; Hess, W. Physica A 1987, 141, 475. (7) Kawasaki, K. J. Stat. Phys. 1997, 87, 981. (8) Wu, J.; Cao, J. J. Phys. Chem. B 2004, 108, 6796. (9) (a) Yang, S.; Cao, J. J. Chem. Phys. 2004, 121, 562. (b) Yang, S.; Cao, J. J. Chem. Phys. 2004, 121, 572. (10) Bicout, D. J.; Szabo, A. J. Chem. Phys. 1998, 109, 2325. (11) Shiwa, Y.; Kawasaki, K. J. Phys. C: Solid State Phys. 1982, 15, 5345. (12) Seki, K.; Barzykin, A. V.; Tachiya, M. J. Chem. Phys. 1999, 110, 7639. (13) Tang, J.; Marcus, R. A. Phys. ReV. Lett. 2005, 95, 107401. (14) Sokolov, I. M. Phys. ReV. Lett. 2003, 90, 080601. (15) Pastor, R. W.; Zwanzig, R.; Szabo, A. J. Chem. Phys. 1996, 105, 3878. (16) Bicout, D. J.; Szabo, A. J. Chem. Phys. 1997, 106, 10292. (17) Doi, M. Chem. Phys. 1975, 11, 107. (18) Portman, J. J.; Wolynes, P. G. J. Phys. Chem. A 1999, 103, 10602. (19) Portman, J. J. J. Chem. Phys. 2003, 118, 2381. (20) Na¨gele, G.; Baur, P. Physica A 1997, 245, 297. (21) Lee, J.; Yang, S.; Kim, J.; Lee, S. J. Chem. Phys. 2004, 120, 7564. (22) Hynes, J. T. J. Phys. Chem. 1986, 90, 3701. (23) (a) Marcus, R. A. Discuss. Faraday Soc. 1960, 29, 21. (b) Marcus, R. A. J. Chem. Phys. 1965, 43, 679. (c) Marcus, R. A.; Sutin, N. Biochem. Biophys. Acta 1985, 811, 265.
