Decay of metastable states: Mean relaxation time formulation

J. Chem. Phys. 110, 7133 (1999); https://doi.org/10.1063/1.478617 110, 7133 © 1999 American Institute of Physics. Decay of metastable states: Mean relaxation time formulation Cite as: J. Chem. Phys. 110, 7133 (1999); https://doi.org/10.1063/1.478617 Submitted: 08 December 1998 . Accepted: 15 January 1999 . Published Online: 07 April 1999 Alexander N. Drozdov, and J. Javier Brey ARTICLES YOU MAY BE INTERESTED IN First passage times for a tracer particle in single file diffusion and fractional Brownian motion The Journal of Chemical Physics 136, 175103 (2012); https://doi.org/10.1063/1.4707349 Generalized moment expansion for Brownian relaxation processes The Journal of Chemical Physics 82, 151 (1985); https://doi.org/10.1063/1.448788 The generalized Smoluchowski equation and non-Markovian dynamics The Journal of Chemical Physics 84, 5824 (1986); https://doi.org/10.1063/1.449892 Decay of metastable states: Mean relaxation time formulation Alexander N. Drozdova) Department of Applied Physics and Chemistry, University of Electro-Communications, 1-5-1 Chofugaoka, Chofu, Tokyo 182-8585, Japan J. Javier Brey Fı ´sica Teo ´rica, Universidad de Sevilla, Apartado de Correos 1065, Sevilla 41080, Spain ~Received 8 December 1998; accepted 15 January 1999! The mean relaxation time formalism introduced by Nadler and Schulten @J. Chem. Phys. 82, 151 ~1985!# in their generalized moment expansion method is extended to a general diffusion process in arbitrary dimensions. The utility of the approach is demonstrated by calculating analytically the rate of noise-induced transitions in a bistable system with an isolated transition point. The rate formula obtained summarizes in a uniform manner much of what had been done before in this field. Limitations of its validity are discussed and a perturbation procedure to systematically improve it is proposed. The validity of our theoretical predictions for the rate is confirmed by comparing with exact numerical results. © 1999 American Institute of Physics. @S0021-9606~99!50415-7# I. INTRODUCTION One of the longstanding problems in physical sciences has been the development of methods to calculate rate constants of equilibration process from a microscopic perspective. A small sample of this work can be found in Refs. 1–3. For an historical review of the field, see that of Landauer.4 Finally, general reviews of the present state of the art have recently been given by Mel’nikov5and Ha ¨nggi, Talkner, and Borkovec6~see also a collection of references in Refs. 7 and 8!. The phenomenon of escape from a locally stable state arises in a multitude of scientific contexts, but the main motivation to study it stems from its connection to chemical kinetics and the theory of diffusion in solids. Since the fundamental contribution of Kramers,9much work has been devoted to the analysis of stochastic models governed by the Fokker–Planck equation. It is a mesoscopic kinetic equation for the distribution function P(q,t) involving a deterministic drift vector Gand a diffusion tensor D. The former describes the deterministic path of the system, while the latter incorporates fluctuations away from this path. Here we deal with the Fokker–Planck equation of a generic type ~the summation rule over repeated indices is always implied, if not stated otherwise!, ] tP~q,t!5LP~q,t![ ] i@2Gi~q!1 ] jDij~q!#P~q,t!, ~1.1! supplemented by the natural boundary conditions nT•J~P!50 for qP ] R,n5normal to ] R,~1.2! where qT5(q1,...,qn)PR, ] Ris the boundary of R, and where we have introduced the Fokker–Planck operator L given by Eq. ~1.1!and the probability current reading Ji~P!5Gi~q!P~q,t!2 ] jDij~q!P~q,t!.~1.3! In what follows we assume that a positive stationary distribution P0(q)5P(q,t→`), the solution of LP050, ~1.4! exists. The quantity of interest is the escape rate of a system from one metastable state into another by crossing an intervening barrier. The flexibility of the above description makes Eq. ~1.1! very attractive for both theoreticians and experimentalists. Specific examples of this equation can be drawn from a vast amount of different fields ranging from nuclear physics to communication theory; accordingly, it has received a great deal of attention in recent years ~see, e.g., textbooks by Gardiner1and Risken2!. Since there is no unique way to determine the escape rate, various different methods have resulted. Here we mention specifically three rather general approaches to this problem. The calculation of the ratio of a stationary current at the top of the barrier to the population of the well, as originally proposed by Kramers,9is the method most frequently used in the past.1,2,9–16 An alternative derivation is based on the mean first passage time formalism.1,2,6,14,16–23 Within its scope, the escape rate is determined as the inverse of the mean time after which a stochastic trajectory starting within the well passes the stochastic separatrix for the first time. Finally, a more precise definition of the kinetic rate is adopted in the eigenmode expansion method.2,16,24–29 The separation of time scales which is typical in barrier crossing processes shows itself in a large gap in the spectrum of the Fokker–Planck operator separating one small nonzero eigenvalue from the rest of the spectrum. In such a case, one can easily recognize the smallest nonzero eigenvalue as the sum of forward and backward rates. Each of the three approaches has its own advantages and drawbacks. The Kramers method is more direct and simpler to use than the two other approaches, though the latter give some additional insight. Besides, the mean first passage time formalism reduces Eq. ~1.1!to the stationary backward Fokker– Planck equation supplemented by absorbing boundary con- a!Permanent address: Institute for High Temperatures, 13/19 Izhorskaya Street, 127412 Moscow, Russia. JOURNAL OF CHEMICAL PHYSICS VOLUME 110, NUMBER 15 15 APRIL 1999 71330021-9606/99/110(15)/7133/9/$15.00 © 1999 American Institute of Physics ditions. For a one-dimensional process, the respective Dirichlet problem is easily solved analytically in quadratures, giving an explicit result for the rate.1,2,6,19 Otherwise, its solution is far from straightforward because absorbing boundary conditions are known to present a special problem for systems with more than one degree of freedom. The eigenmode expansion method also reduces Eq. ~1.1!to a stationary problem but with the natural boundary conditions ~1.2!. Although the class of bistable processes for which the first nonzero eigenvalue may be calculated exactly is rather limited,2,30 the method is in a sense superior over the mean first passage time approach. From a mathematical point of view, diffusion problems of such a type are easier to solve since no external boundary conditions are required to determine the escape rate. In the present paper we propose an alternative method to obtain the rate which combines the principal advantages of the mentioned approaches. The basic idea is to define this quantity in terms of an integral time characteristic, a socalled mean relaxation time, introduced by Nadler and Schulter31 in their generalized moment expansion method ~for applications of the method, see also Ref. 32!. The remainder of the paper is organized as follows. In Sec. II the method of Nadler and Schulten is briefly reviewed and its connection to the method of eigenmode expansion is established. An approximate rate expression for a general twostate diffusion process is derived in Sec. III. The derivation is based on the assumption that the potential barrier is high enough compared to the noise strength so that any nonlinearity in the original Fokker–Planck operator can safely be neglected. In Sec. IV a perturbation approach to take into account the nonlinear contribution is outlined. In this way an improved rate expression is derived that involves the leading nonvanishing correction in terms of the inverse barrier height. In Sec. V, the accuracy of our theoretical predictions is tested by comparing with exact numerical rates in one- and two-dimensional potentials. The final section contains some concluding remarks. II. PRELIMINARIES A. Phenomenological description As a preliminary we briefly outline the problem of interest and a phenomenological approach to its solution. To model a two-state process, the system ~1.1!is assumed to have two domains of attraction RAand RBcorresponding to metastable states Aand B, respectively, and one transition attractor, referred to henceforth as a saddle S, which is unstable in the direction transverse to the stochastic separatrix ] V. Here we restrict our consideration to the case of a saddle point qSlocated at the origin, qS50. The generalization to an arbitrary unstable attractor is straightforward. It is clear that a detailed solution of the partial differential equation is a far from simple task not only in many dimensions but even in one dimension. To simplify the problem, one usually assumes that the dynamics of the two-state process obeys the phenomenological linear rate equations,33 A kB kAB~2.1! and ] tNA~t!52kANA~t!1kBNB~t!, ~2.2! ] tNB~t!5kANA~t!2kBNB~t!. Here NA(t) and NB(t) are the time-dependent nonequilibrium probabilities to find the system in, respectively, the A and the Bmetastable state, Ni~t!5 E Rdq s i~q!P~q,t!,~2.3! defined relative to a surface separating Aand B, s i~q!5 H 1, qPRi 0, otherwise, i5A,B,~2.4! while kAand kBare the rates of escaping from these states related to each other through the equilibrium fractions kANA e5kBNB e,Ni e5Ni~`!,i5A,B.~2.5! As the system of interest is assumed to be closed, NA1NB 51, the equations can be reduced to one equation for a function given by d N~t!5NA~t!2NA e.~2.6! This yields a single-exponential approximation of the dynamics, d N~t! d N~0!5exp~2kt!,~2.7! which depends upon the full equilibration rate, k5kA1kB. In this way the complicated nonstationary problem ~1.1!is reduced to a more simple problem of determining the equilibration rate. B. Generalized moment expansion method Next we briefly review the generalized moment expansion method of Nadler and Schulten.31 To this end, let us consider the equilibrium time correlation function of an observable M(t), C~t!5 ^ d M~0! d M~t! & 5 ^ d Mexp~tL1! d M & .~2.8! Without loss of generality we assume that C(t) is normalized, C(0)51. In the above L1is the backward Fokker– Planck operator L15Gi~q! ] i1Dij~q! ] ij 2,~2.9! d M5M2 ^ M & , while ^&denotes the average with respect to the stationary solution of the Fokker–Planck equation. The average has the properties of an inner product on the space of functions d M, i.e., ^ d M~0! d M~t! & 5 E RdqP0~q! d M~q!etL1 d M~q!. ~2.10! 7134 J. Chem. Phys., Vol. 110, No. 15, 15 April 1999 A. Drozdov and J. Brey The generalized moment expansion method starts with the Laplace-transformed correlation function C ˆ~ v !5 E 0 `dte2 v tC~t!5 ^ d M~ v 2L1!21 d M & ~2.11! and involves a two-point Pade ´-approximation around v 5`, C ˆ~ v !52 ( m50 ` t m~21/ v !m11,~2.12! t m5~21!m F dmC~t! dtm G t50 ,~2.13! and v 50, C ˆ~ v !5( m50 ` t 2m21~2 v !m,~2.14! t 2m2151 m! E 0 `dttmC~t!.~2.15! The expansion coefficients t mfor both m.0 and m,0 are evaluated from t m5~21!m ^ d M~L1!m d M & .~2.16! They contain information on the short- and long-time scales, respectively. Since it is impossible to determine and sum infinitely many terms in the series representations of Eqs. ~2.12!and ~2.14!, these must be truncated at some finite m 5iand m5j. The resulting approximation to C(t) is a superposition of n exponentials that exactly reproduce iderivatives of C(t)att50@Eq. ~2.13!# and jrelaxation moments @Eq. ~2.15!#, C~t!5( l51 n alexp~2 h lt!,~2.17! where i1j115 n and ( m 51 n a m ~ h m !m5 t m,2j<m<i.~2.18! The approximation includes in a balanced way both short- and long-time effects. However, the utility of the generalized moment expansion method appears to be restricted to onedimensional stochastic processes, because just in that case the moments of the backward Fokker–Planck operator L1 are available analytically.31,32 III. MEAN RELAXATION TIME FORMALISM To proceed further we note that by the Onsager regression hypothesis d N(t) decays to equilibrium in the same fashion as does the equilibrium time correlation function C~t!5 ^ d N~0! d N~t! & ^ d N~0! d N~0! & ,~3.1! i.e.,33 C~t!5 d N~t!/ d N~0!.~3.2! Consequently, by Eq. ~2.7!C(t);exp(2kt) at long times when the phenomenology is valid. Criteria for the validity of Eq. ~2.2!will be discussed in Sec. V of the present paper ~see also Ref. 34!. Here we only note that the dynamics of physical and chemical processes is often not known in detail. It is frequently the case that the largest relaxation time is the sole kinetic quantity obtainable from experiments. Therefore effective methods for its evaluation is one of the most fundamental problems in physics and chemistry. The method we outline below is an extension of the elegant formalism of Nadler and Schulten to a general twostate Fokker–Planck process. As we are interested in the long-time regime where the decay is unconditionally single exponential, we limit ourselves to n 51, i.e., C~t!5exp~2t/ t 21!.~3.3! In that event, the long-time rate is given by the inverse of the so-called mean relaxation time31 reading ~from here on, we shall drop the index 21 to keep the notation simple! t 5 E 0 `dtC~t!.~3.4! Interpretations of this quantity are the following. Should one think the decay be single exponential, as is the case for ~2.7!, then t 215k, and therefore kA5NB e t 21,kB5NA e t 21.~3.5! Furthermore, if one expands the Green function of the Fokker–Planck equation P~q,t u q0![etL d ~q2q0!5Pm 1~q0!Pm~q!e2lmt,m>0 ~3.6! in terms of the complete set of eigenfunctions ~below summation over the index iis not implied!, LPi52liPi,L1Pi 152liPi 1, E RdqPi 1~q!Pj~q!5 d ij,~3.7! one obtains C~t!5cmexp~2lmt!,m>1, ~3.8! and, accordingly, t 5cmlm 21~3.9! with ci5~NA eNB e!21 E RA dqPi~q! E RA dqPi 1~q!P0~q!.~3.10! Here the expansion coefficients satisfy 0<cm,1, Scm51. It is clear that the rate description is meaningful if 12c1!1, ~3.11! in which case l15 t 21.~3.12! The validity of the approximate relations ~3.3!and ~3.12! will be discussed in Sec. V of the present paper. 7135J. Chem. Phys., Vol. 110, No. 15, 15 April 1999 A. Drozdov and J. Brey On the other hand, by Eq. ~2.16!the formal expression for the mean relaxation time reads t 5~NA eNB e!21 E Rdq d N~q!P0~q!F~q!,~3.13! where d N(q)5 s A(q)2NA eand F(q) is determined from L1F~q!52 d N~q!,~3.14! on the space of functions orthogonal to P0, E RdqP0~q!F~q!50. ~3.15! In case the escape dynamics entails a one-dimensional Fokker–Planck process in the interval @a,b#, Eq. ~3.14!can be readily solved exactly to give31 t 5~NA eNB e!21 E a bdx@D~x!P0~x!#21 H E a xdy d N~y!P0~y! J 2. ~3.16! Yet another special case that allows the explicit expression for the mean relaxation time includes diffusive problems in a spherically symmetric potential, in a spherically symmetric domain. Other than that there is no closed form solution to Eq. ~3.14!. Here we deal with the problem in the limit of small diffusion coefficients ~or large barriers! Dij~q!5«D ¯ ij~q!,«!1, ~3.17! where the noise strength «measured in units of the barrier height is assumed to be small. In such a case, if the system starts within RAit will typically first approach the metastable state Aand stay within its neighborhood for a long-time T until an occasional fluctuation drives the system to the product state B. Accordingly, we may seek the function F(q)in the form F~q!5Tf~q!,~3.18! where the constant part Tis obtained by multiplying ~3.14! by P0and integrating over RA. This yields T52NA eNB e F E RA dqP0~q!L1f~q! G 21,~3.19! leading us finally to t 52 F E RA dqP0~q!L1f~q! G 21 E RA dqP0~q!f~q!. ~3.20! The above formula is formally exact in the sense that no approximation has been made to derive it. Unfortunately, it is still too complicated to use for quantitative calculations. Thus approximation schemes must be invoked. In order to construct an approximate solution we use the fact that Tis exponentially large, T;exp(«21), and hence the inhomogeneity 2 d N(q)/Tmay safely be neglected in Eq. ~3.14!, i.e., L1f50. ~3.21! For vanishing diffusion the resulting first-order partial differential equation has only solutions that are piecewise constant on the domains of attraction of the deterministic equations of motion. We choose these values as being NB eand 2NA e, such that the resulting function is orthogonal to the stationary solution. The presence of small diffusive terms in L1changes the behavior of fonly near the deterministic separatrix where the steplike behavior is smoothed out. Since the expression for the mean relaxation time, Eq. ~3.13!, contains the stationary distribution as a weight, only the region of the saddle point is of importance for sufficiently small «. Hence we may split the backward operator into a leading contribution L0 1 and a correction L1 1reading L15L0 11L1 1,~3.22! where L0 15Bi jqj ] i1D ˜ ij ] ij 2, Bi j5 ] jGi~q! u q5qS,D ˜ ij5Dij~qS!~3.23! describes the linear dynamics near the saddle and L1 15DGi~q! ] i1DDij~q! ] ij 2, DGi~q!5Gi~q!2Bi jqj,DDij~q!5Dij~q!2D ˜ ij ~3.24! the anharmonic correction. When writing Eq. ~3.24!we have taken into account the fact that for «→0 the saddle point coincides with a stationary point of the drift vector, i.e., G(qS)50. In passing we note that rescaling the coordinates qiby the square root of the noise strength renders the leading contribution L0 1independent of «and the correction L1 1pro- portional to the noise strength. Consequently, the nonlinearity L1 1may be considered as a small perturbation. A perturbation series for the function freads16,23,28 f5f01f11f21¯.~3.25! Neglecting L1 1completely we arrive at L0 1f050. ~3.26! Within this approximation the solution of Eq. ~3.21!reads f05N2 F N11 A 2/ p E 0 brdsexp~2s2/2! G ,~3.27! where N1and N2are constants of integration and b5 A m /Drr,r5w•q.~3.28! In the above, Drr denotes the r,r-component of the diffusion Drr5wTD ˜ w,~3.29! while m is the positive eigenvalue of the matrix Bwith the associated eigenvector w, i.e., Bi jwi5 m wj.~3.30! The chain of approximations made to construct ~3.27!constitutes the essence of the standard saddle point approximation method.11,15,20,27,28 One of the disadvantages of the method is that the resulting rate expression is asymptotic in 7136 J. Chem. Phys., Vol. 110, No. 15, 15 April 1999 A. Drozdov and J. Brey the barrier height. A way to improve it by taking into account finite-barrier corrections will be discussed in Sec. IV of the present paper. In the remainder of this section we restrict our discussion to the lowest order contribution f0. It is not hard to see that the constant N2involved in f0can be set to unity without loss of generality, N251. Next, the constant N1defined by Eq. ~3.15!is easily evaluated explicitly using the fact that P0has sharp maxima in the vicinity of the stable attractors A and B. There, the approximate solution f0is nearly constant. Consequently, the constant N1is given by N15~NB e2NA e! a 0.~3.31! Hereby the vector wwas assumed to be oriented so that ris positive at Aand negative at B. The numerator of Eq. ~3.20! can also be evaluated in the same manner to yield E RA dqP0~q!f~q!52 h 0NA eNB e.~3.32! Usually, one argues that the constants a 0and h 0appearing in the above equations are unity up to exponentially small corrections in the barrier height, and therefore can be neglected. We will see in the following that for low barrier heights these corrections also become important. Thus, using only the lowest order approximation f0, the mean relaxation time formalism leads to the following equilibration rate: k052~2 h 0NA eNB e!21 E RA dqP0~q!L1 1f0~q!,~3.33! where we have employed Eq. ~3.26!. The above formula reproduces various different rate expressions available in the literature for multidimensional Fokker–Planck processes.9–11,14,15,17,21 In order to illustrate this statement, let us explicitly evaluate the numerator of Eq. ~3.33!. After a partial integration the integral entering the numerator reads I5 E ] RA dSi@Ji~P0!f01P0Dij ] jf0#,~3.34! where dSdenotes the oriented surface element on ] RA.It should be noted here that the stationary probability current J(P0) does not in general vanish. As f0is practically constant away from the saddle, the surface integral over f0J(P0) vanishes. Furthermore, the main contribution to the surface integral over the second term comes from the neighborhood of the saddle where only the derivative with respect to r contributes while the other derivatives are negligible. Therefore, we may approximate the separatrix in the vicinity of that point by its tangential hyperplane r50, I' E r50dSrP0Drr ] rf052bDrr A 2/ p E r50dSrP0, ~3.35! where the minus sign arises from the opposite directions of dSiand ] r/ ] qi. Introducing a d -function into the integral on the right of ~3.35!we eventually arrive at k05bDrr h 0NA eNB e A 2 p E RdqP0~q! d ~w•q!.~3.36! Before proceeding further three remarks are in order. First, we note that Eqs. ~3.33!and ~3.36!do not necessarily coincide with each other. The difference may arise due to exponentially small corrections introduced by Eq. ~3.35!. Second, when applied to the Kramers problem and its generalization to time-dependent friction, Eq. ~3.33!agrees with rate formulas derived by Dekker25 and Talkner and Braun,21 respectively. Third, all the above-mentioned rate formulas presuppose the knowledge of the stationary solution. Since P0is positive, we may write P0~q!5 H E Rdqexp@2F~q!# J 21exp@2F~q!#,~3.37! where F~q!is a so-called generalized potential. The potential is easily obtained in closed form @by simple integration of Eq. ~1.4!# when the underlying dynamics obeys detailed balance.1,2 Otherwise, no general method exists for computing P0,35,36 and therefore the escape dynamics must be studied in each particular case separately.36,37 This makes establishing general properties of ka quite difficult task, which lies outside the scope of the present paper anyway. Here we assume that the generalized potential is known. In the weak noise limit, Eq. ~3.17!,F~q!is the Lyapunov function of the deterministic dynamics; consequently, it must be minimal on the stable attractors and maximal on the saddle of the deterministic dynamics. Assuming that all these attractors are point attractors and expanding F~q!in their neighborhood as F~q!5FM1 w ij Mqiqj1O~ u q u 3!,M5A,B,S,~3.38! allows one to perform all the integrals appearing in Eq. ~3.36!analytically. This leads in a straightforward way ~for more details, see the Appendix!to a rate formula of the standard form,14,15,20,27 kst5 m 2 p A u det w S u @ A det w Aexp~FA2FS! 1 A det w Bexp~FB2FS!#.~3.39! Some other approximate rate expressions available in the literature for particular Fokker–Planck processes9–11 follow from Eq. ~3.39!in a very natural fashion. Finally, to conclude this section we mention two drawbacks of the above steepest-descent approximation. First of all we note that using Eq. ~3.38!implies the differentiability of the generalized potential. The latter, however, is not a generic case. In many nontrivial examples the generalized potential shows singularities.35 Second, when deriving Eq. ~3.39!we have neglected in the series ~3.38!terms of order higher than u q u 2. This neglect is quite severe and introduces a large error that has to be compensated for by a large potential barrier. 7137J. Chem. Phys., Vol. 110, No. 15, 15 April 1999 A. Drozdov and J. Brey IV. FINITE-BARRIER CORRECTIONS As we already noted, the method used for approximately solving Eq. ~3.21!is asymptotic. It rests on the assumption that the barrier height measured in units of the noise intensity is sufficiently large not only in order that a rate description is meaningful but, additionally, that on the diffusional length scale at the saddle any nonlinearity can be neglected. This condition, though, is not always met in realistic barrier crossing processes. The problem of finite-barrier corrections has been studied actively in recent years and many different methods have been suggested for their evaluation.16,23,24,26–28,38,39 However, the utility of the results obtained so far is restricted to particular ~one- and twodimensional!Fokker–Planck equations. The aim of this section is to outline a general approach to the problem. An obvious way to improve the present solution is to take into account the nonlinearity L1 1@Eq. ~3.24!# fully neglected in the derivation of Eq. ~3.27!. Based on the series representation ~3.25!, a perturbation theory with respect to the nonlinearity L1 1can be performed, leading to a hierarchy of inhomogeneous equations of the form16,28 L0 1fm52L1 1fm21,m>1, ~4.1! with f0being determined by ~3.27!. Here we limit ourselves to first order in the perturbation L0 1f152L1 1f0.~4.2! Then, splitting off a Gaussian function from f1, f1~q!5b~2/ p !1/2 exp~2b2r2/2!h~q!,~4.3! and using the explicit expressions ~3.27!and ~3.24!for f0 and L1 1, respectively, we arrive at the following equation for the function h: ~L0 122rb2wiD ˜ ij ] j2 m !h5rb2wiDDijwj2DGiwi. ~4.4! The latter can be solved systematically by expanding h, DGi, and DDij, h5H1Hiqi1Hijqiqj1Hijkqiqjqk1O~ u q u 4!, DGi5Bi jkqjqk1Bi jkmqjqkqm1O~ u q u 4!,~4.5! DDij5Dij kqk1Dij kmqkqm1O~ u q u 3!, and equating like powers of qi. This results in a set of algebraic equations for the coefficients H,Hi,Hij, and Hijk reading ~2Bi j24b2wmD ˜ miwj2 m d ij!Hik5b2Dim kwiwmwj2Bi jkwi, H5~2/ m !D ˜ ijHij, ~4.6! ~3Bi j26b2wmD ˜ miwj2 m d ij!Hikp 5b2 2Dim kpwiwmwj21 6wiBi jkp , ~ m d ij12b2wmD ˜ miwj2Bi j!Hi56D ˜ miHmij . Once the first-order perturbation correction is known, the mean relaxation time formalism gives the following improved expression for the equilibration rate: k152~2 h 1NA eNB e!21 E RA dqP0L1 1f1,~4.7! where we have used L0 1f050 and the defining Eq. ~4.2!for f1. After a partial integration this rate can be rewritten in the form k15~ h 0/ h 1!k02~2 h 1NA eNB e!21 E ] RA dSiP0Dij ] jf1. ~4.8! Hereby we introduced the parameter h 1, h 15~2NA eNB e!21 E RA dqP0~f01f1! 5 h 01~2NA eNB e!21 E RA dqP0f1,~4.9! which is unity up to exponentially small corrections in the barrier height. As will be shown in the next section, these corrections are important for low barrier heights. V. APPLICATIONS A. A one-dimensional model Since exact results for the mean relaxation time are only available in one dimension, we first consider a typical model given by ] tP~x,t!5 ] x~x32x1D ] x!P~x,t!.~5.1! The dynamics is that of a Brownian particle moving in the symmetric bistable potential F~x!5~x221!2 4D,~5.2! in the large damping limit. The height of the potential barrier, DF5F(0)2F(61), is related to the noise strength D via DF5~4D!21.~5.3! Our aim is to test the validity of the rate description @Eqs. ~2.7!and ~3.12!#, as well as the utility of the various approximations discussed above for the escape rate. The mean relaxation time of the Fokker–Planck process ~5.1!was calculated in terms of Eq. ~3.16!by numerical quadratures. Numerically exact results for the first nonzero eigenvalue l1and the equilibrium time correlation function C(t) were obtained by a basis set method described earlier.29 In Fig. 1 these results are compared to the single exponential approximation C(t)5exp(2kt) taken with k5l1and k 51/ t . The sharp minima obtained with k51/ t mean that the corresponding approximation crosses the true correlation function at l1t51. The calculation clearly demonstrates that depending on the barrier height there are two different mechanisms determining the equilibration process. For low barrier heights (DF,5) the system belongs to the relaxational regime. In this regime the equilibration rate is not 7138 J. Chem. Phys., Vol. 110, No. 15, 15 April 1999 A. Drozdov and J. Brey very small, while the coefficient c1noticeably deviates from unity. Accordingly, the long-time behavior is determined by a set of low-lying eigenvalues. This indicates that the rate description as a whole loses its meaning for too low barriers. It is remarkable, however, that even for a low barrier (DF ;1) a large share of the system (c1;0.9) decays with one and the same rate l1and therefore the error made by using the single exponential approximation ~2.7!is relatively small ~just of order 10% as seen from Fig. 1!. On the contrary, at high barriers the system belongs to the rate regime. In that even the long-time behavior is characterized by a least nonvanishing eigenvalue that is separated from all larger ones by an exponentially large gap. Moreover, almost the whole system decays with the same rate, because c1'1 in this regime. The above observations are well illustrated by Fig. 2, which shows log10(12c1) as a function of the barrier height. The logarithmic plot clearly demonstrates that the difference 1 2c1decreases with DF exponentially so that already for DF*5 the dynamics is single exponential in the entire time domain. The same is found to be true for h 0, h 1, and l1 t . All these quantities also very rapidly approach unity with increasing barrier height ~the corresponding results for the product l1 t are not shown in Fig. 2, as they are indistinguishable from those for c1!. Finally, we compare in Fig. 3 the exact mean relaxation time @Eq. ~3.16!# to its zeroth- @Eq. ~3.33!# and first- @Eq. ~4.7!# order approximations. As we have seen from Fig. 2, for low barriers the numerator of the mean relaxation time is rather sensitive to the detailed shape of the function f(x); accordingly, the factors h 0and h 1in the denominator of Eqs. ~3.33!and ~4.7!must not be approximated by unity. Though the resulting corrections are ‘‘exponentially small’’ they do not much differ in magnitude from the ‘‘leading’’ algebraic corrections for low barrier heights. In particular, for DF51 the factors h 0and h 1are 0.80 and 0.71, respectively. Taking this fact into account reduces the error in the corresponding approximations by a factor of 3. With increasing barrier height the effect disappears. For DF*5, the inequalities 0.995< h 1, h 0,1 usually hold; consequently, the factors h 0and h 1can safely be neglected. One may thus conclude that in the rate regime, deviations between the numerically exact results and approximate rate expressions are essentially due to finite-barrier effects. As evidenced by Fig. 3, first-order corrections lead to a considerable improvement of the rate in this regime, reducing the error by factors, or even by orders of magnitude. The excellent agreement demonstrates the potential of the present approach. It is capable of describing rather accurately the largest relaxation time in a bistable potential over a broad range of barrier heights. Our comparison also includes results obtained with the standard rate expression ~3.39!. The latter is seen to be least accurate. It systematically overestimates the exact equilibration rate for all DF>2. FIG. 1. Logarithm of the relative error @~approximate!-~exact!#/~exact!in the equilibrium time correlation function @Eq. ~3.1!# for DF51 and 6 made by using the single exponential approximation, Eq. ~2.7!with k5l1~dashed lines!and k51/ t ~solid lines!. FIG. 2. Logarithm of the deviation of the different coefficients from unity. Solid line: y5c1with c1defined by Eq. ~3.10!; dashed line: y5 h 0, Eq. ~3.32!; dot-dashed line: y5 h 1, Eq. ~4.9!. FIG. 3. Relative errors in the inverse mean relaxation time made by using different approximate expressions. Dot-dashed line: standard rate formula kst , Eq. ~3.39!; dashed line: zeroth-order approximation k0, Eq. ~3.33!; open circles connected by dashed line: zeroth-order approximation evaluated with h 051, h 0k0; solid line: first-order approximation k1, Eq. ~4.7!; solid circles connected by solid line: first-order approximation evaluated with h 151, h 1k1. 7139J. Chem. Phys., Vol. 110, No. 15, 15 April 1999 A. Drozdov and J. Brey B. A double well coupled to a harmonic mode As a second example, we consider a model that has been introduced for studying charge-transfer reactions in polar solvents.40 It consists of a reactive coordinate x~describing the nonlinear system of interest!and a relaxational ~harmonic!mode ythat mimics a slowly relaxing polarization. In the limit of strong damping in both the x- and y-direction the modes undergo a diffusional motion with generally different diffusion constants. In dimensionless variables the dynamics is given by the Smoluchowski equation ] tP~x,y,t!5D~ ] xe2F ] xeF1 j] ye2F ] yeF!P~x,y,t!, ~5.4! with the potential of mean force F(x,y), reading F~x,y!51 D F U~x!1 g 2~y2x!2 G , U~x!51 4x421 2x2.~5.5! Here, g is the coupling constant, Dis a dimensionless temperature related to the barrier height by Eq. ~5.3!, and the anisotropy parameter j is the ratio of the damping constant in the direction of system and polarization coordinates xand y, respectively. Numerically exact results for the least nonvanishing eigenvalue of Eq. ~5.3!have been calculated in a previous paper.28 These are compared in Fig. 4 to the standard rate formula ~5.3!, the asymptotic mean relaxation time expression ~3.33!, and its value including the first-order corrections ~4.7!. As anticipated, the standard rate expression is the worst approximation to the least nonvanishing eigenvalue. It agrees relatively well with the exact results for large parameters j and DF. Large deviations, however, are found if both j and DF are small, i.e., in the case of a slow harmonic mode. The deviations are due to finite-barrier heights. It is seen that the zeroth-order mean relaxation time formula ~3.33!is in better agreement with numerical calculations than the standard rate expression for all values of j and DF. Finite-barrier corrections according to Eq. ~4.7!lead to a further considerable improvement of the rate. These corrections must include h 1for DF,5, in which case h 1cannot be approximated by unity. We have found that the factor h 1 strongly depends on the barrier height and is rather insensitive to j and g . Finally, we note that the errors of the various approximate rate expressions slowly decrease for a stronger coupling ~ g .1!and increase with decreasing g while the overall situation with varying g remains similar to that for g 51. VI. DISCUSSION In the present paper the mean relaxation time formalism of Nadler and Schulten31 that had been previously used for solving one-dimensional problems was extended to a general diffusion process in arbitrary dimensions. The utility of the approach has been demonstrated by calculating analytically the equilibration rate in a bistable system without manifest detailed balance. An asymptotic formula that recovers various different results already available in the literature for the equilibration rate is constructed by using the standard Kramers function. The derivation of this function is based on the assumption that the potential barrier is high enough compared to the noise strength so that any nonlinearity in the Fokker–Planck operator can be neglected. An improved function is constructed from a perturbation theory that gives the standard Kramers function as zeroth-order solution. In this way finite-barrier corrections to the rate are taken into account. The present formalism is as simple to implement as the Kramers method and still gives some additional insight. Its connection to the method of eigenmode expansion has been established. The validity of the theoretical predictions has been tested by comparing with exact numerical rates in one- and two-dimensional potentials. In all the cases considered we obtained excellent agreement between the theory and estimates of the rate from numerical calculations. This demonstrates the potential of the present approach. It is capable of describing rather accurately the largest relaxation time in a bistable system over a broad range of barrier heights. ACKNOWLEDGMENTS This work was supported in part by the Direccio ´n General de Ensen ˜ anza Superior e Investigacio ´n Cientı ´fica of FIG. 4. Same as in Fig. 3 but for the two-dimensional model, Eqs. ~5.4!and ~5.5!with g 51 and j 50.1, 1, and 10. Solid circles are for log10(12 h 1). 7140 J. Chem. Phys., Vol. 110, No. 15, 15 April 1999 A. Drozdov and J. Brey