EQUILIBRIUM AND RECROSSINGS OF THE TRANSITION STATE: WHAT CAN BE LEARNED FROM DIFFUSION?

by

Z. Schuss^{1}

Abstract

.

We examine the assumptions and conclusion of (generalized) transition state theory (GTST) by considering the activation process in the diffusion (Langevin) limit. We find the asymptotic structure of the leading eigenfunctions and eigenvalues of the Fokker-Planck operator with a bistable potential, and hence the long time quasi equilibrium behavior of the phase space probability density function (pdf). Defining reactant and product as small neighborhoods Ω_{A} and Ω_{B} of the stable states A and B, respectively, we examine all possible recrossings of the transition state region (TSR) and find their contribution to the mean first passage time (MFPT) τ_{AB} from Ω_{A} to Ω_{B}. We show that the mean number of recrossings of the TSR is 1, hence τ_{AB} = 2τ_{AS}, where τ_{AS} is the MFPT from Ω_{A} to the stochastic separatrix S, which we use as a generalized transition state (GTS). The activation rate, that is, the rate at which trajectories arrive to Ω_{B} from Ω_{A}, is then shown to be given by κ = 1/2τ_{AS}, and in the limit of small noise is independent of the choice of Ω_{A} and Ω_{B}. We conclude that to obtain the correct rate in (G)TST (i) the quasi equilibrium density (qepdf) rather than the equilibrium density (epdf) has to be used, (ii) the qepdf contains a boundary layer near the stochastic separatrix, but otherwise the reactant qepdf≈epdf, and (iii) all recrossings of the (G)TS are accounted for if (G)TS=S, but not otherwise. We also consider the case of a single metastable state.

1Department of Mathematics, Tel-Aviv University, Ramat-Aviv, Tel-Aviv 69978, Israel

1. Introduction.

Transition state theory (TST) of activation and its generalizations (GTST) (see e.g., [1]) are based on the following assumptions (1) an equilibrium probability distribution of phase space trajectories, (2) that a certain configuration space surface, called the (generalized) transition state (G)TS, has the property that trajectories that cross the (G)TS from reactant to product never recross it. Because of (1) half the trajectories on the (G)TS cross in the product direction, and because of (2) they end up as product. It follows therefore that (3) the reaction rate κ is the outgoing half of the equilibrium phase space probability flux on the (G)TS, normalized by the reactant population. The underlying postulate which justifies the equilibrium assumption is that after a sufficiently long time the quasi equilibrium density of the reactant is sufficiently similar to the equilibrium density so that the outgoing half of the quasi equilibrium flux on the (G)TS, normalized by the quasi equilibrium reactant population, is well approximated by that of the equilibrium flux, normalized by the equilibrium reactant population.

Both theoretical [2] and experimental results [3] show that κ depends on the dissipation, contrary to the TST prediction. The failure of TST to show this dependence brought (1) and (2) into question. It has been widely recognized that both (1) and (2) are unrealistic, because activation is a non equilibrium process and because recrossings do occur. It is believed therefore that due to (2), (G)TST overestimates the rate [4]. There have been different attempts to modify the equilibrium density on the (G)TS to account for the non equilibrium nature of activation [5]-[7]. Another method for circumventing the equilibrium assumption is Kramers’ method of stationary flux, in which a source is placed at the bottom of the reactant well and an absorbing barrier is imposed sufficiently far outside it. Then the steady state normalized outgoing flux on the (G)TS is related to the rate. Attempts to account for recrossings of the (G)TS introduced a transmission factor k [8], and assumed a relation between k and the (mean) number of crossings ⟨n⟩ [3],[9],[10]. Also the choice of the (G)TS in (2) has been the subject of study and different authors proposed different surfaces, e.g., a surface through the saddle point which is perpendicular to the equipotential surfaces [4], the configuration of least probability [11], and a surface of minimal flux [12]-[14].

It seems that no satisfactory theoretical explanation of the failure of the equilibrium assumption has been offered, other than the incorrect predictions it implies. The methods used to modify the epdf in order to obtain the qepdf are mostly based on ad hoc assumptions about the behavior of the trajectories or of the pdf on the chosen (G)TS [5]-[7]. Also the choice of the (G)TS and the treatment of recrossings have not been adequately discussed in the literature.

In order to clarify some of these problems we consider them in the diffusion (Langevin, or Fokker-Planck) limit [2], [15]. The underlying assumption in this study is that the behavior of the phase space trajectories of the diffusion process approximates well that of the trajectories of a particle coupled to a bath of oscillators [16]-[18]. Such a limit corresponds to the Fokker-Planck approximation of the Liouville evolution equation for the probability density function (pdf) of the phase space trajectories [19]. The Fokker-Planck equation describes the evolution of the the pdf p(q, p, t) of the phase space trajectories of the particle, whereas the Liouville equation describes that of the joint pdf p(q, p, x, y, t) of the particle and the oscillators in multidimensional phase space. Thus,

p(q, p, t) = lim_{N→∞} ∫ ∫ p(q, p, x, y, t) dx dy,  (1.1)

where x = (x_{1}, …, x_{N}) and y = (y_{1}, …, y_{N}) are the phase space coordinates of the bath oscillators. In this approximation the behavior of the random trajectories governed by a (generalized) Langevin equation (GLE) [13], [20] approximates that of the ensemble of trajectories of he particle, if the initial states of the bath oscillators are chosen at random [17]. Therefore the study of the above mentioned problems in the diffusion limit may provide some insight into (G)TST and its variants.

We take advantage of the eigenfunction expansion of the nonequilibrium pdf p(q, p, t|q_{0}, p_{0})

p(q, p, t|q_{0}, p_{0}) = ∑_{n=0}^{∞} ϕ_{n}(q, p)ψ_{n}(q_{0}, p_{0})e^{−λ_{n}t},  (1.2)

where ϕ_{n} and ψ_{n} are the eigenfunctions of the Fokker-Planck operator (FPO) L and its adjoint L^{∗}, respectively, corresponding to the eigenvalue λ_{n}, and (q_{0}, p_{0}) is a point in the reactant part of phase space. We denote by D_{A} and ∂D_{A} the domain of attraction of the equilibrium state at the bottom of the reactant well and its boundary, respectively.

First we examine the equilibrium assumption (1) and its consequences. Since λ_{0} = 0 it follows from (1.2) that for times

\frac{1}{λ_{2}} ≪ t ≪ \frac{1}{λ_{1}}  (1.3)

the rate of change of the reactant population is λ_{1}, whereas for times

t ≫ \frac{1}{λ_{1}}  (1.4)

it vanishes. It is shown in Section 4 that

ϕ_{1}(q, p)ψ_{1}(q_{0}, p_{0}) ≫ ϕ_{0}(q, p)ψ_{0}(q_{0}, p_{0}) for (q, p) ∈ D_{A},  (1.5)

and ψ_{1}(q_{0}, p_{0}), ψ_{1}(q_{0}, p_{0}) are independent of (q_{0}, p_{0} in D_{A}. Thus for times (1.3) exp{−λ_{1}t} = O(1) so that the reactant pdf p(q, p, t|q_{0}, p_{0}) is given by

p(q, p, t|q_{0}, p_{0}) ∼ e^{−λ_{1}t}ϕ_{1}(q, p)ψ_{1}(q_{0}, p_{0}).  (1.6)

It follows that the qepdf (1.6) of the reactant is represented by ϕ_{1} and not by the equilibrium density ϕ_{0}. We conclude that the root cause of error in TST is the equilibrium assumption (1). There is however a certain degree of similarity between the equilibrium and the quasi equilibrium reactant densities. We show in Section 5 that

ϕ_{1}(q, p) ∼ ϕ_{0}(q, p)ψ_{1}(q, p),  (1.7)

in D_{A}, where ψ_{1} is a boundary layer function which connects smoothly to a large constant C_{A} in the interior of D_{A} and to \frac{1}{2}C_{A} on ∂D_{A}. This similarity explains the need for introducing corrections to the equilibrium density [5]-[7]. From the boundary layer structure (1.7) of ϕ_{1} we conclude that in order to obtain the qepdf (1.6) from the epdf a modification has to be introduced on ∂D_{A} and it must have the boundary layer form ψ_{1}.

Next we examine the problem of recrossings of the (G)TS. We show in Section 4 that in the limit of high barrier the stochastic separatrix S, the locus of points from which trajectories are equally likely to become reactant and product, is asymptotically ∂D_{A}. From the definition of S we find that average number of times ⟨n⟩, that a trajectory crosses a neighborhood of S before becoming a product is 2, so that k = 1. The calculation of the rate from the normalized flux on the (G)TS is based on the identity

λ_{1} = \frac{∫_{∂D} J(ϕ_{1}) · ν ds}{∫_{D} ϕ_{1} dq dp},  (1.8)

where J(ϕ_{1}) is the flux density corresponding to the quasi equilibrium density ϕ_{1}, D is the reactant region in phase space, and ν is the unit outer normal to ∂D. If ∂D = S is chosen, then it follows from (1.8) and (1.7) that the outgoing half of the quasi equilibrium flux on S, normalized by the reactant population, is the rate. We conclude that if the (G)TS is chosen to be the stochastic separatrix, then the relation (3) between the rate and the quasi equilibrium flux holds, but not otherwise.

Our conclusions about (G)TST can be summarized as follows. To obtain the correct rate in (G)TST (i) the qepdf the reactant is obtained from the epdf by a modification of the epdf, (ii) the modification has to have the form of a boundary layer at the stochastic separatrix S, and (iii) all recrossings are accounted for if (G)TS=S, but not otherwise. The structure of S in multidimensional phase space is needed only near the saddle point. In the high barrier limit it can be determined from the memory function in the GLE [21]. The determination of S in other cases has been discussed in [22], [23].

We find the quasi stationary behavior of the solution of the Fokker-Planck equation (FPE) by finding the asymptotic structure of the leading eigenvalues and eigenfunctions in the limit of high barrier. We consider both bistable and metastable dynamics, and discuss Kramers’ method of stationary flux. We find the average number of times a trajectory crosses a neighborhood of an equiprobable surface S_{p} before becoming a product (S_{1/2} ≡ S). Hence we find the relation between the flux on S and the rate. In Section 2 we formulate the problem and discuss the concept of escape. In Section 3 we discuss the relation between the MFPT, the absorption rate, the escape rate, and the eigenvalues of the FPO in bounded and unbounded domains. In Section 4 we find S and the asymptotic structure of the leading eigenfunctions. In Section 5 we consider the problem of recrossings and of their contribution to the escape rate. Finally, in section 6 we summarize the results and the conclusions of this paper.

2. Formulation.

We consider the system of n Itô stochastic differential equations

ẋ = b(x) + \sqrt{2ϵ} σ(x) ẇ,  (2.1)

where b(x) = (b^{1}(x), …, b^{n}(x)) is a smooth vector field in R^{n}, ẇ is a vector of k independent standard Gaussian white noises (k ≤ n), and σ(x) is an n × k noise matrix. The parameter ϵ is a measure of the noise intensity and is assumed to be small relative to other parameters of the problem such as the size of b(x) and σ(x), the sizes of the domains in which (2.1) is considered, and so on. The diffusion matrix a(y) ≡ {a^{ij}(y)} is given in terms of the noise matrix {σ^{ij}(y)} as

a(y) = σ(y)σ^{T}(y).  (2.2)

The noiseless dynamics

ẋ = b(x)  (2.3)

is assumed to have a finite attractor A with domain of attraction D_{A} whose boundary ∂D_{A} is a repeller (see fig.1). We assume that all domains mentioned below have boundaries with a piecewise continuous normal. When (2.1) is used to model an activated process by the Langevin dynamics of a particle in a potential well, activation occurs when a trajectory leaves a neighborhood Ω_{A} of A and does not return for a long period of time. Thus the activation rate κ is the escape rate from Ω_{A}. To be meaningful such a definition should be independent of the choice of Ω_{A}, of the location of ∂D, and of the noise strength ϵ, provided ϵ is small.

As noted in the Introduction, several definitions have been given for the GTS. Calculations of κ from the rate at which trajectories arrive at boundaries of various domains, have been used. One such calculation of κ is based on the rate at which trajectories are absorbed in the boundary ∂Ω_{A}, or equivalently, κ = 1/⟨τ(Ω_{A})⟩, where ⟨τ(Ω_{A})⟩ is the mean first passage time (MFPT) to ∂Ω_{A} [24], [25]. If Ω_{A} is replaced by the domain of attraction D_{A} of A, then for small ϵ, as noted in [26]-[28] trajectories arriving at ∂D_{A} are equally likely to return to Ω_{A} prior to absorption in the boundary ∂D of a larger domain D containing D_{A}, as they are to be absorbed before they return to Ω_{A}. Therefore the calculation of κ was modified in [22], [26], [29], to κ = 1/2⟨τ(D_{A})⟩. Another calculation identifies κ as the principal eigenvalue λ_{1} of the Fokker-Planck equation in the whole space [30] or in a specified domain D [31]. Not all such calculations are consistent with one another or with physically measured quantities. Thus for example, the normalized outgoing flux at the top of a potential barrier may not represent the rate, since trajectories arriving at the top of a potential barrier with outward pointing velocity may have a non-negligible probability of returning to the well.

The definition of an escape rate κ from Ω_{A} as the absorption rate κ_{abs}(D) in ∂D requires the calculation of the MFPT ⟨τ(D)⟩ from Ω_{A} to ∂D, or equivalently, of the principal eigenvalue λ_{1}(D) of the Fokker-Planck operator with absorbing boundary conditions on ∂D. The calculation of ⟨τ(D)⟩ and of λ_{1}(D) for multidimensional systems in domains D containing an attractor has a long history. Several types of domains have been considered, which we now classify. A domain D is classified as type I if D ⊂ D_{A} and the drift enters D. In this case ∂D is a noncharacteristic boundary. If D = D_{A}, we classify it as type II. Here the system does not drift into D_{A}, ∂D is a characteristic boundary, and is an invariant manifold for the drift. If D ⊃ D_{A} and the drift in D − D̄_{A} is attracted by ∂D, we classify D as type III (see the exact definition below).

The MFPT for domains of type I was calculated in [24], [25], [32] with increasingly more accurate results. Escape rates from domains of type II with drift derivable from a potential were calculated in [33]-[35] by a generalization of Kramers’ theory to higher dimensions. For general systems, not derivable from a potential, ⟨τ(D)⟩ was calculated in [36]-[39]. For domains of type III neither ⟨τ(D)⟩ nor λ_{1}(D) seem to have been calculated, other than for the one dimensional Kramers problem [2].

Definition. A domain D is of type III relative to (2.3) if

(i) D ⊃ D̄_{A}, (ii)

δ_{1} < dist(x, ∂D) < δ_{2} for all x ∈ ∂D_{A},  (2.4)

and δ_{1}, δ_{2} are positive constants, independent of ϵ, and (iii) trajectories of (2.3) which start in D − D̄_{A} reach ∂D in finite time.

An often encountered example of a domain of type III is the following. Let (2.3) be a bistable system with a finite attractor A and a second, more stable attractor B (possibly at infinity). If Ω_{B} is a neighborhood of B, inside the domain of attraction of B, then D ≡ R^{n} −Ω̄_{B} is a domain of type III relative to (2.3). This example corresponds, e.g., to a model of a chemical reaction. The noisy dynamics (2.1) models the motion of an atom bound by a stable chemical bond, corresponding to the attractor A. When the bond is broken due to the noise (e.g., to molecular collisions), a new more stable bond B is formed [2]. Another example of a domain of type III relative to (2.3) corresponds to the dynamics (2.1) with a single metastable state. In this case D contains D_{A} and trajectories which start in D −D̄_{A} are assumed to exit D in finite time and never return. This corresponds, e.g., to dissociation, where the metastable state A models the undissociated state of a molecule [40].

Because of the weak fluctuations, the trajectories of the system tend to stay in a small neighborhood Ω_{A} of A, with rare escapes to Ω_{B} or to an absorbing boundary ∂D, from which they do not return. In the case of bistability Ω_{B} usually corresponds to a more stable state of the system, so that the time spent in Ω_{B} is usually much longer than in Ω_{A}. Thus we may assume that trajectories arriving in Ω_{B} are absorbed there and never return to A. Therefore the bistable case can also be treated as metastable with absorption in ∂D.

The transition probability density function (tpdf) of trajectories of (2.1), p(y, t|x) ≡ P r{x(t) = y|x(0) = x}, satisfies the forward Kolmogorov (Fokker-Planck) equation [41]

\frac{∂p}{∂t} = L_{ϵ}p ≡ −∇·J,  (2.5)

where the probability current density J(x, y, t) is defined as

J^{i}(x, y, t) ≡ −ϵ ∑_{j=1}^{n} \frac{∂}{∂y^{j}} [a^{ij}(y)p(y, t|x)] + b^{i}(y)p(y, t|x).  (2.6)

The initial condition for (2.5) is given by

p(y, 0|x) = δ(x − y).  (2.7)

If trajectories of (2.1) are absorbed on the boundary ∂D of a given domain D, the tpdf satisfies the boundary condition

p(y, t|x) = 0 for x ∈ D, y ∈ ∂D.  (2.8)

It also satisfies the backward Kolmogorov equation

\frac{∂p}{∂t} = L_{ϵ}^{∗}p ≡ ϵ ∑_{i,j=1}^{n} a^{ij}(x)\frac{∂^{2}p}{∂x^{i}∂x^{j}} + ∑_{i=1}^{n} b^{i}(x)\frac{∂p}{∂x^{i}} for x ∈ D,  (2.9)

and

p = 0 for x ∈ ∂D, y ∈ D.  (2.10)

We assume that the tpdf p(y, t|x) can be represented by the eigenfunction expansion

p(y, t|x) = ∑_{n=1}^{∞} ψ_{n}(x)ϕ_{n}(y)e^{−λ_{n}(D)t},  (2.11)

where

L_{ϵ}ϕ_{n}(y) = −λ_{n}(D)ϕ_{n}(y) for y ∈ D,  (2.12)
ϕ_{n}(y) = 0 for y ∈ ∂D,  (2.13)

and

L_{ϵ}^{∗}ψ_{n}(x) = −λ_{n}(D)ψ_{n}(x) for x ∈ D,  (2.14)
ψ_{n}(x) = 0 for x ∈ ∂D.  (2.15)

This is the case for example, if D is a bounded domain, however we do not restrict our considerations to this case. The eigenvalues are ordered as 0 < λ_{1}(D) ≤ ℜe λ_{2}(D) ≤ ℜe λ_{3}(D) ≤ …, with λ_{1}(D) real [42]. We assume that the eigenfunctions are normalized so that

∫_{D} ϕ_{1}(y) dy = 1,  (2.16)

and are biorthogonal

∫_{D} ϕ_{m}(x)ψ_{n}(x) dx = δ_{mn}.  (2.17)

In Section 4 we calculate the probability P (x) of return from a point x to Ω_{A}, and in Section 6 we identify it asymptotically as the principal eigenfunction ψ_{1}(x). In Section 5 we calculate the principal eigenvalue λ_{1}(D) = κ = 1/⟨τ(D)⟩ for domains of type III, and prove

lim_{ϵ→0} \frac{λ_{1}(D)}{λ_{1}(D_{A})} = \frac{1}{2}.  (2.18)

3. The MFPT, the absorption rate, and the principal eigenvalue

Let D be a bounded domain which contains A. Then, under some mild positivity assumptions about a(x), the trajectories of the noisy system (2.1) which start in D, cross ∂D in finite (random) time τ with probability 1. Moreover, the MFPT

⟨τ(x)⟩ ≡ E(τ|x(0) = x)  (3.1)

from a point x in D to ∂D is finite and is independent of the behavior of the process at the boundary. These results may hold as well even if D is unbounded. Thus we may assume that ∂D is an absorbing boundary, so that any trajectory of (2.1) that reaches ∂D is instantaneously terminated. The MFPT is given by [32]

⟨τ(x)⟩ = ∫_{0}^{∞}∫_{D} p(y, t|x) dy dt,  (3.2)

or by

⟨τ(x)⟩ = ∫_{D} p(x, y) dy,  (3.3)

where p(x, y) is defined by

p(x, y) ≡ ∫_{0}^{∞} p(y, t|x) dt.  (3.4)

¿From the Fokker-Planck equation (2.5) and the initial condition (2.7) we find that

L_{ϵ}p(x, y) = −δ(x − y).  (3.5)

If δ(x − y) is replaced by f(y) in the initial condition (2.7), the pdf p(x, y) is independent of x, that is, p(x, y) = p(y), and satisfies the equation

L_{ϵ}p(y) = −f(y).  (3.6)

In particular, if we choose f(y) = ϕ_{1}(y), then p(y) = ϕ_{1}(y)/λ_{1}(D) and (3.3) implies that

⟨τ(D)⟩ ≡ ∫_{D} ϕ_{1}(x)⟨τ(x)⟩ dx = \frac{1}{λ_{1}(D)}.  (3.7)

The structure of the functions p(x, y) and p(y) has been discussed in [29]. Equation (3.5) with absorbing boundary conditions on ∂D represents Kramers’ method of stationary flux, in which a source is placed inside D and an absorbing barrier is imposed outside D_{A}. Equation (3.7) relates the solution of (3.5) with the escape rate. The relation of the out going flux on the (G)TS, as obtained from (3.5), and the rate is discussed in Sections 4 and 6.

Next we discuss the relationships between the absorption rate κ_{abs}(D) in ∂D, the MFPT ⟨τ(D)⟩, and the principal eigenvalue λ_{1}(D). The absorption rate constant κ_{abs}(D) represents the stationary relative rate of change of the population N(x, t) in D, given the initial condition (2.7), where

N(x, t) ≡ ∫_{D} p(y, t|x) dy.  (3.8)

The Fokker-Planck equation (2.5) implies that

\frac{∂N(x, t)}{∂t} = ∫_{D} \frac{∂p(y, t|x)}{∂t} dy = − ∫_{D} ∇_{y} · J(x, y, t) dy
= − ∫_{∂D} J(x, y, t) · ν(y) dS_{y} ≡ −F (x, t).  (3.9)

Recalling that J is the probability current density, we see that the rate of change of the population in D is the negative of the total probability flux on ∂D. It follows that the relative rate of change is given by

κ_{abs}(x, t) ≡ −\frac{∂N(x, t)/∂t}{N(x, t)} = \frac{F (x, t)}{N(x, t)}.  (3.10)

Using the eigenfunction expansion (2.11) we obtain

κ_{abs}(x, t) = \frac{∑_{n} λ_{n}(D)psi_{n}(x) ∫ ϕ_{n}(y) dy e^{−λ_{n}(D)t}}{∑_{n} ψ_{n}(x) ∫ ϕ_{n}(y) dy e^{−λ_{n}(D)t}} .  (3.11)

Obviously the relative rate of change (3.11) depends on both x and t, however its limit as t → ∞ is constant. Thus the concept of a rate constant for (2.1) is meaningful only in the limit t → ∞, since the relative rate of change becomes constant only after a sufficiently long time. Therefore the natural definition of the absorption rate constant κ_{abs}(D) is given as

κ_{abs}(D) ≡ lim_{t→∞} κ_{abs}(x, t),  (3.12)

which, in view of (3.11) is

κ_{abs}(D) = λ_{1}(D).  (3.13)

The corresponding notion of the MFPT must be defined in a manner consistent with the definition (3.13) of rate as a long time limit. After a long time t_{0} has elapsed, a large proportion of the trajectories have already been absorbed. Thus we have to consider only those trajectories that survived (were not absorbed) in ∂D at time t_{0}. The probability density of finding a trajectory at a particular point y, among the surviving trajectories in D at time t_{0}, is the conditional density

p_{C}(y, t_{0}|x) ≡ \frac{p(y, t_{0}|x)}{∫_{D} p(y, t_{0}|x) dy}.  (3.14)

Using the eigenfunction expansion (2.11) and the normalization (2.16), we obtain

p_{C}(y, t_{0}|x) = \frac{∑_{n} ψ_{n}(x)ϕ_{n}(y)e^{−λ_{n}(D)t_{0}}}{∑_{n} ψ_{n}(x) ∫_{D} ϕ_{n}(y) dye^{−λ_{n}(D)t_{0}}} → ϕ_{1}(y) as t_{0} → ∞.  (3.15)

Thus the conditional density p^{C}(y, t) of trajectories at time t+t_{0}, which were observed at time t_{0} ≫ 1 in D, is the solution of the Fokker-Planck equation (2.5) with the absorbing boundary condition (2.8), and the initial condition

p^{C}(y, 0) = ϕ_{1}(y).  (3.16)

Thus, for times 1/λ_{2} ≪ t ≪ 1/λ_{1} the qepdf is given by

p^{C}(y, t) = ϕ_{1}(y)e^{−λ_{1}(D)t}.  (3.17)

Therefore the MFPT for a trajectory, observed at a point y at a sufficiently large time t_{0}, must be calculated by (3.2), with the conditional density p^{C}(y, t) of surviving trajectories replacing the tpdf p(y, t|x). Now the concept of first passage time, which is consistent with the definition (3.12) of κ_{abs}(D), is the first passage time for the process (2.1) which starts with the initial density (3.16). The long time MFPT is therefore given by

⟨τ(D)⟩ = ∫_{0}^{∞}∫_{D} p^{C}(y, t) dy dt = \frac{1}{λ_{1}(D)}.  (3.18)

We see that, according to this definition, both ⟨τ(D)⟩ and κ_{abs}(D) are independent of initial conditions. Although p^{C}(y, 0) is not a density of the process x(t), defined by (2.1), it represents the shape of the long time tpdf with the decay e^{−κ_{abs}(D)t}.

Now we consider the FPE (2.5) in R^{n}. We assume that (2.3) is a bistable system with attractors A and B, whose domains of attraction, D_{A} and D_{B}, are separated by ∂D_{A}. We denote by ϕ_{i}, ψ_{i}, and λ_{i} the corresponding eigenfunctions and eigenvalues. We have λ_{0} = 0 and ϕ_{0} is the equilibrium pdf. We write (2.12) as

−∇ · J(ϕ_{n}) = −λ_{n}ϕ_{n},  (3.19)

where J(ϕ_{n}) is defined in (2.6) with p replaced by ϕ_{n}(y). Then for any domain Ω ⊂ R^{n}

λ_{n} = \frac{F_{n}}{N_{n}} ,  (3.20)

where

F_{n} ≡ ∫_{∂Ω} J(ϕ_{n}) · ν ds  (3.21)

and

N_{n} ≡ ∫_{Ω} ϕ_{n} dy.  (3.22)

In particular we choose Ω = D_{A}.

Next we describe the structure of the MFPT and its dependence on D. For domains of type I it was shown and II [25], [32], [43] that

⟨τ(D)⟩ = O(ϵ^{1/2}e^{Ψ̂(D)/ϵ}),  (3.23)

where Ψ̂(D) is a positive constant which depends on D and on the coefficients a^{ij}(x) and b^{i}(x) in D. Explicit expressions for Ψ̂(D) were given in [24], [25], [32] for various domains in terms of an eikonal function Ψ(x), which is a solution of the Hamilton-Jacobi type equation

∑_{i,j=1}^{n} a^{ij}(x)\frac{∂Ψ(x)}{∂x^{i}} \frac{∂Ψ(x)}{∂x^{j}} + ∑_{i=1}^{n} b^{i}(x)\frac{∂Ψ(x)}{∂x^{i}} = 0.  (3.24)

The constant Ψ̂(D) is given by

Ψ̂(D) = min_{x∈∂D} Ψ(x).  (3.25)

The function Ψ(x) is defined by the assumption that the principal eigenfunction ϕ_{0} has the WKB structure in D

ϕ_{0}(y) = K(y, ϵ)e^{−Ψ(y)/ϵ},  (3.26)

where K(y, ϵ) is a regular function of ϵ. The structure of Ψ(x) for bistable systems (2.3) is more complicated. If we assume that (2.3) has attractors A and B, with domains of attraction D_{A} and D_{B}, respectively, then a global smooth solution to (3.24) in R^{n} may not exist [44]. However locally, in D_{A} and D_{B}, (3.26) holds. In the case of detailed balance [41] (3.26) holds in R^{n} with K(y, ϵ) = 1 and Ψ(y) is the energy, so that ϕ_{0}(y) is the Boltzmann equilibrium density.

The function Ψ(x) decreases on the trajectories of (2.3) in D_{A}, so that Ψ̂(D_{1}) ≥ Ψ̂(D_{2}) if D_{1} and D_{2} are of type I and D_{1} ⊃ D_{2}, hence by (3.23)

\frac{⟨τ(D_{1})⟩}{⟨τ(D_{2})⟩} = O(e^{[Ψ̂(D_{1})−Ψ̂(D_{2})]/ϵ}).  (3.27)

If in addition δ_{1} ≤ dist(x, ∂D_{2}) ≤ δ_{2} for all x ∈ ∂D_{1}, where δ_{1} and δ_{2} are positive numbers independent of ϵ, then Ψ̂(D_{1}) > Ψ̂(D_{2}), so that

lim_{ϵ→0} \frac{⟨τ(D_{2})⟩}{⟨τ(D_{1})⟩} = 0.  (3.28)

This however, is not the case for domains of type III, since for such domains the limit (2.18) implies that ⟨τ(D)⟩ is to leading order independent of D.

We now discuss the relation between absorption and escape. We distinguish between two types of escape. Let Ω_{A} be a neighborhood of A in D_{A}. One type of escape is the event of reaching ∂Ω_{A} for the first time, regardless of the subsequent behavior of the trajectory. This type of escape describes several physical systems in which trajectories are terminated at ∂Ω_{A} [25], [45]. The rate of this type of escape is simply the absorption rate in ∂Ω_{A}. It involves neither multistability nor separation of time scales. A second type of escape involves multistability and separation of time scales. According to (3.23) trajectories spend time O(exp(Ψ̂(Ω_{A})/ϵ) in Ω_{A} before reaching ∂Ω_{A} (necessarily crossing it several times in rapid succession). If a trajectory reaches a distance O(1) beyond ∂Ω_{A} in R^{n} − Ω̄_{A}, it either returns to Ω_{A} or is absorbed in ∂D without returning to Ω_{A}. In the first instance, if the time of return is short relative to ⟨τ(Ω_{A})⟩, the trajectory is not considered to have escaped, whereas in the latter it is. This distinction means that relatively short departures do not represent transitions to a new physical state, e.g., going from reactant to product in a chemical reaction, whereas long departures do represent such transitions. The rate at which long departures occur can be defined in terms of an absorption rate in ∂D, where D is a domain of type III.

4. The probability of returns, ψ_{1}, ϕ_{1}, and λ_{1}

In this section we consider the probability P (x) of return from a point x to Ω_{A}, prior to absorption in ∂D, where D is a domain of type III, and its relation to the eigenfunctions ψ_{1} and ϕ_{1} in R^{n} and in domains of type II and type III. First we summarize the results of [27] for drifts which are continuous across ∂D_{A} (e.g., for smooth potential barriers), and then show that similar results hold for drifts whose normal component suffers a discontinuity across ∂D_{A} (e.g., for sharp potential barriers).

We denote by τ(∂Ω_{A}) (τ(∂D)) the first passage time to ∂Ω_{A} (∂D). The function

P (x) ≡ Pr{τ(∂Ω_{A}) < τ(∂D)|x_{ϵ}(0) = x}  (4.1)

is the probability that a trajectory starting at x will reach ∂Ω_{A} before ∂D. The stochastic separatrix S [22], [26]-[28], [46], [47] is defined as the locus of points x such that

P (x) = 1/2.  (4.2)

The probability P (x) is the solution of

L_{ϵ}^{∗}P (x) = 0 in D − Ω̄_{A}  (4.3)

with the boundary conditions

P (x) = 1 for x ∈ ∂Ω_{A}  (4.4)

and

P (x) = 0 for x ∈ ∂D.  (4.5)

(see e.g. [32]). In Section 6 we show that asymptotically (4.3)-(4.5) also define the principal eigenfunction ψ_{1}(x) of L_{ϵ}^{∗}.

The outer expansion of P (x) in D is found by assuming the regular expansion

P (x) = P^{0}(x) + o(1) as ϵ → 0.  (4.6)

In view of (3.23) the leading term P^{0}(x) must satisfy the reduced equation

∑_{i=1}^{n} b^{i}(x)\frac{∂P^{0}(x)}{∂x^{i}} = 0,  (4.7)

or equivalently,

\frac{dP^{0}(x(t))}{dt} = 0,  (4.8)

where x(t) is any trajectory of (2.3). It follows that P^{0}(x) = const. on every trajectory. By assumption, all trajectories outside D_{A} reach ∂D in finite time, so that the boundary condition (4.5) implies that

P^{0}(x) = 0 for x ∈ D − D̄_{A}.  (4.9)

Similarly, all trajectories in D_{A} reach ∂Ω_{A} in finite time, hence by (4.4)

P^{0}(x) = 1 for x ∈ D_{A}.  (4.10)

For continuous drifts the discontinuity in the outer solution is bridged by the uniform expansion [27], [28], given by

P (x) = \sqrt{\frac{1}{2π}} ∫_{−∞}^{χ(x)/\sqrt{ϵ}} e^{−z^{2}/2} dz.  (4.11)

The function χ(x) satisfies (see (4.3))

∑_{i=1}^{n} b^{i}(x)\frac{∂χ(x)}{∂x^{i}} = ∑_{i,j=1}^{n} a^{ij}(x) (χ(x)\frac{∂χ(x)}{∂x^{i}} \frac{∂χ(x)}{∂x^{j}} − ϵ\frac{∂^{2}χ(x)}{∂x^{i}∂x^{j}}) .  (4.12)

The boundary conditions (4.4) and (4.5) imply that

χ(x) > 0 for x ∈ D_{A} − Ω̄_{A}, outside a boundary layer  (4.13)

and

χ(x) < 0 for x ∈ D − D̄_{A}, outside a boundary layer.  (4.14)

The boundary conditions (4.4) and (4.5) are satisfied by (4.11) asymptotically. In view of (4.2) S is defined by χ(x) = 0. Now (4.13) and (4.14) imply that lim_{ϵ→0} S = ∂D_{A} (see the discussion of this argument in Section 6). Expanding χ(x) = χ^{0}(x) + o(1) as ϵ → 0 we obtain that on a trajectory x = x(t) of (2.3).

\frac{d}{dt}χ^{0}(x) = ∑_{i=1}^{n} b^{i}(x)\frac{∂χ^{0}(x)}{∂x^{i}} = ∑_{i,j=1}^{n} a^{ij}(x)χ^{0}(x)\frac{∂χ^{0}(x)}{∂x^{i}} \frac{∂χ^{0}(x)}{∂x^{j}}.  (4.15)

In view of (4.13), (4.14), and (4.15), χ^{0}(x) increases on the trajectories of (2.3) outside D_{A} and decreases inside D_{A}.

Next we consider drifts b(x) whose normal component b_{n}(x) suffers a discontinuity across ∂D_{A}. That is, we assume that ± lim_{ρ→0±} b_{n}(x) > 0. This is the case, e.g., if the potential barrier is sharp (see fig.1). The solution of (4.3) is still given by (4.11), however the local behavior of χ(x) near ∂D_{A} is not as given in [27]. The boundary layer to bridge the discontinuity of the outer solution across ∂D_{A} can be found in this case by introducing the local coordinates (ρ, s) near ∂D_{A}, where ρ(x) ≡ dist(x, ∂D_{A}) and s(x) = (s^{2}, …, s^{n}) are local coordinates in ∂D_{A}. Then we stretch ρ by setting

ζ ≡ \frac{ρ}{ϵ},  (4.16)

Q(ζ, s) ≡ P (x), and expanding

Q(ζ, s) ∼ Q^{0}(ζ, s) + ϵQ^{1}(ζ, s) + ….  (4.17)

Now (4.3) is to leading order

a(s)\frac{∂^{2}Q^{0}(ζ, s)}{∂ζ^{2}} + b_{n}^{±}(s)\frac{∂Q^{0}(ζ, s)}{∂ζ} = 0 for ±ζ > 0,  (4.18)

where

a(s) ≡ lim_{ρ→0} ∑_{i,j=1}^{n} a^{ij}(x)\frac{∂ρ(x)}{∂x^{i}} \frac{∂ρ(x)}{∂x^{j}}  (4.19)

and

b_{n}^{±}(s) ≡ lim_{ρ→0±} b(x) · ∇ρ(x).  (4.20)

Since ∂D_{A} is a repeller we have ±b_{n}^{±}(s) > 0. The matching conditions are

lim_{ρ→∞} Q^{0}(ρ, s) = 1, lim_{ρ→−∞} Q^{0}(ρ, s) = 0.  (4.21)

The solution is given by

Q^{0}(ζ, s) = 1 + \frac{b_{n}^{−}(s)}{b_{n}^{+}(s) − b_{n}^{−}(s)} exp\{−\frac{b_{n}^{+}(s)}{a(s)}ζ\} if ζ > 0  (4.22)

and

Q^{0}(ζ, s) = \frac{b_{n}^{+}(s)}{b_{n}^{+}(s) − b_{n}^{−}(s)} exp\{−\frac{b_{n}^{−}(s)}{a(s)}ζ\} if ζ < 0.  (4.23)

The stochastic separatrix S, given to leading order by Q^{0}(ζ, s) = 1/2, is found to be

ρ(s) = \frac{ϵa(s)}{b_{n}^{±}(s)} log \frac{±2b_{n}^{±}(s)}{b_{n}^{+}(s) − b_{n}^{−}(s)} if \frac{±2b_{n}^{±}(s)}{b_{n}^{+}(s) − b_{n}^{−}(s)} > 1.  (4.24)

Thus lim_{ϵ→0} S = {ρ(s) = 0} = ∂D_{A} in this case as well. Therefore (4.11) still holds, but with the local behavior (4.22) and (4.23) rather than that in the continuous case (see [27]).

Next we consider the eigenvalue problems (2.12) and (2.14) in R^{n} and in domains of type II and III. We we begin with a bistable system (2.3) in R^{n}, with attractors A and B. The domains of attraction D_{A} and D_{B} are separated by ∂D_{A}. We denote by Ψ_{A(B)}(y) the solution of (3.24) in D_{A(B)} and set ΔΨ_{A(B)} ≡ min_{∂D_{A}} Ψ_{A(B)}(y) − Ψ_{A(B)}(A(B)). In the case of bistable Langevin dynamics ΔΨ_{A(B)} is the barrier height of the A(B) well. If

ΔΨ_{A} < ΔΨ_{B},  (4.25)

then

λ_{1} = O(e^{−ΔΨ_{A}/ϵ}).  (4.26)

It follows that (2.14) with n = 1 is asymptotically the same as (4.3).

To find the expansion of ψ_{1} in R^{n} we begin with the outer expansion

ψ_{1} ∼ ψ_{1}^{0} + ϵψ_{1}^{1} + …  (4.27)

and find that

ψ_{1}^{0} ∼ C_{A(B)} in D_{A(B)},  (4.28)

where C_{A} and C_{B} are constants [32]. The boundary layer connecting the constants across ∂D_{A} is similar to P (x) in (4.11),

ψ_{1}(y) ∼ \frac{C_{A} − C_{B}}{\sqrt{2π}} ∫_{0}^{χ(y)/\sqrt{ϵ}} e^{−s^{2}/2} ds + \frac{C_{A} + C_{B}}{2} .  (4.29)

The orthogonality condition (2.17) and the WKB assumption (3.26) give

0 = ∫_{R^{n}} ϕ_{0}(x)ψ_{1}(x) dx = ∫_{R^{n}} ψ_{1}(x)K(x, ϵ)e^{−Ψ(x)/ϵ}
∼ (2πϵ)^{n/2} ∑_{i=A,B} C_{i}K(i, ϵ)H^{−1/2}(Ψ(i))e^{−Ψ(i)/ϵ} ≡ aC_{A} + bC_{B},  (4.30)

where H(Ψ(i)) is the Hessian of Ψ at i [32]. Assuming without loss of generality that min_{∂D_{A}} Ψ_{A}(x) = min_{∂D_{A}} Ψ_{B}(x) (this is the case, e.g., if detailed balance holds), (4.25) implies that

\frac{a}{b} = O(e^{−[Ψ(A)−Ψ(B)]/ϵ}) ≪ 1.  (4.31)

Hence C_{B} ≪ C_{A} so that

ψ_{1}(x) ∼ C_{A}P (x),  (4.32)

(see (4.11)). In particular we have ψ_{1} ∼ C_{A}/2 on ∂D_{A}. We normalize ψ_{1} below.

Next we determine ϕ_{1} in R^{n}. We consider two cases, K = 1 in (3.26) and K ≠ 1. For K ≠ 1 the structure of the principal eigenfunction ϕ_{1}(y) is a little more complicated than that for K = 1. The substitution

ϕ_{1}(y) = q(y)e^{−Ψ(y)/ϵ}  (4.33)

transforms (2.12) into

L̃_{ϵ}ϕ_{1}(y) = O(ϵ),  (4.34)

where L̃_{ϵ} is an operator of the form (2.9) with b(x) replaced by

b̃_{ϵ}^{i}(y) ≡ − (b^{i}(y) + ∑_{j=1}^{n} 2a^{ij}(y)\frac{∂Ψ(y)}{∂y^{j}} ) + O(ϵ).  (4.35)

It is easy to see from (3.24) that in appropriate local variables near the critical points of b(x), the dynamics (2.3) and

ẋ = b̃(x)  (4.36)

have the same equilibrium points with the same local stability properties. It follows that A is an attractor for both (2.3) and (4.36), and similarly saddle points are preserved. The domains of attraction D_{A} and D̃_{A} of (2.3) and (4.36), respectively, however may be different. In the case D_{A} = D̃_{A} the structure of q(y) is the same as that of ψ_{1}(y). This happens, for example, if (2.3) has no equilibrium points on ∂D_{A} (e.g., if D_{A} is a limit cycle in two-dimensions, see the analysis in [29, Section 4.2]), or if detailed balance holds. If D_{A} ≠ D̃_{A}, their intersection still contains a neighborhood Ω_{A} of A. Thus the outer expansion of q(y) is still C_{A} in Ω_{A}. From (4.33) it follows that ϕ_{1}(y) is sharply peaked at A, so that the averaging in (3.7) asymptotically gives

\frac{1}{λ_{1}(D)} = ⟨τ(D)⟩ ∼ ⟨τ(A)⟩.  (4.37)

Since ⟨τ(x)⟩ is asymptotically independent of x in D_{A} [32], the point A in (4.37) can be replaced by any point x in D_{A} outside a boundary layer near ∂D_{A}.

We assume K = 1, then [29]

ϕ_{1}(x) = ϕ_{0}(x)ψ_{1}(x).  (4.38)

Using the normalization condition (2.17)) with m = n = 1 we obtain from (4.30), (4.29), and (4.38)

aC_{A}^{2} + bC_{B}^{2} = 1,  (4.39)

and from (4.30) and (4.39)

C_{A} = \sqrt{\frac{b}{a}}, C_{B} = \sqrt{\frac{a}{b}}.  (4.40)

¿From the structure of the eigenfunctions we can determine the quasi stationary density in D_{A}. For times such that

\frac{1}{λ_{2}} ≪ t ≪ \frac{1}{λ_{1}}  (4.41)

the eigenfunction expansion (1.2) is given by

p(y, t|x) ∼ ϕ_{0}(y)ψ_{0}(x) + e^{−λ_{1}t}ϕ_{1}(y)ψ_{1}(x),  (4.42)

where exp{−λ_{1}t} = O(1). ¿From (4.40) it follows that ϕ_{0}(y) ≪ ϕ_{1}(y), hence in D_{A}

p(y, t|x) ∼ e^{−λ_{1}t}ϕ_{1}(y)ψ_{1}(x).  (4.43)

However for times t ≫ 1/λ_{1}

p(y, t|x) ∼ ϕ_{0}(y)ψ_{0}(x).  (4.44)

It follows that the rate of change of the population in D_{A} satisfies

−\frac{Ṅ}{N} ∼ λ_{1}  (4.45)

for times (4.41), however

−\frac{Ṅ}{N} ∼ 0  (4.46)

for times t ≫ 1/λ_{1}.

To calculate the flux F_{1} on ∂D_{A} we note that the drift b(x) is tangent to ∂D_{A}, so that the contribution of the second term on the right hand side of (2.6) to F_{1} vanishes. It follows that

∫_{∂D_{A}} J(ϕ_{1}) · ν ds = ∫_{∂D_{A}} ψ_{1}(x)J(ϕ_{0}) · ν ds + ∫_{∂D_{A}} ϕ_{0}(x)J(ψ_{1}) · ν ds
= \frac{C_{A}}{2} ∫_{∂D_{A}} J(ϕ_{0}) · ν ds + C_{A} ∫_{∂D_{A}} ϕ_{0}(x)ϵ ∑_{i,j} a^{ij}ψ_{1,j}(x)ν^{i} ds
∼ C_{A}\sqrt{\frac{ϵ}{2π}} ∫_{∂D_{A}} ϕ_{0} ∑_{i,j} a^{ij}χ_{j}(x)ν^{i} ds.  (4.47)

Hence, by (3.20),

λ_{1} = \frac{F_{1}}{N_{1}} ∼ \sqrt{\frac{ϵ}{2π}} ∫_{∂D_{A}} ϕ_{0} ∑_{i,j} a^{ij}χ_{j}ν^{i} ds/ ∫_{D_{A}} ϕ_{0} dx ≡ \frac{F̃_{1}}{N_{0}} .  (4.48)

Thus λ_{1} is the total flux F̃_{1} on ∂D_{A}, normalized by the equilibrium reactant population. Equations (4.48) and (4.46) indicate that the (G)TST assumption that the equilibrium flux is a good approximation to the quasi equilibrium flux is incorrect.

Next we consider the principal eigenvalue problem (2.12), (2.13) in domains of type II, that is, we choose D = D_{A}. If absorbing conditions are imposed on ∂D_{A}, then λ_{0} = 0 is no longer an eigenvalue, since no equilibrium density exists. The principal eigenfunction ψ_{1}^{D_{A}}, which satisfies the boundary condition ψ_{1}^{D_{A}} = 0 on ∂D_{A}, is given by

ψ_{1}^{D_{A}} = ψ_{1} − \frac{C_{A}}{2} ,  (4.49)

because for x ∈ ∂D_{A} (4.11) and (4.32) imply that ψ_{1} ∼ 1/2. Hence

ϕ_{1}^{D_{A}} ∼ ϕ_{0}(ψ_{1} − \frac{C_{A}}{2} ),  (4.50)

and consequently

λ_{1}(D_{A}) ∼ \frac{F_{1}(D_{A})}{N_{1}(D_{A})} = \frac{F̃_{1}}{\frac{1}{2}N_{0}} = 2λ_{1}.  (4.51)

It follows that the absorption rate in ∂D_{A} is twice the rate, that is, (2.18).

Finally, if D is a domain of type III, then asymptotically ψ_{1}^{D} ∼ ψ_{1}, and (4.51) holds with λ_{1} replaced by λ_{1}(D). The result (4.51) can be stated as

Theorem 1. If D is a domain of type III relative to (2.3), then

lim_{ϵ→0} \frac{λ_{1}(D)}{λ_{1}(D_{A})} = \frac{1}{2}.  (4.52)

Thus to leading order in ϵ, λ_{1}(D) is independent of D, as long as D is of type III. That is, changes in the shape of D do not affect λ_{1}(D) to leading order. In contrast, as D shrinks to D_{A}, λ_{1}(D) abruptly doubles as D nears D_{A}. As it shrinks yet further inside D_{A}, the quotient λ_{1}(D)/λ_{1}(D_{A}) increases exponentially fast in 1/ϵ.

A generalization of Theorem 1 to the eigenvalue problem in R^{n} can be stated as

Theorem 2. If (2.3) is a bistable system and D = R^{n}, then (4.52) holds.

Theorems 1 and 2 show that in the limit of high barrier the definitions of rate as an absorption rate of trajectories outside D_{A} or as the rate of change of the reactant population, or as twice the absorption rate on S, or as the stationary flux on an absorbing boundary ∂D with a source in D_{A}, all lead the same result.

We illustrate the difference between our non equilibrium considerations and the equilibrium considerations of (G)TST with an example. In the one dimensional Kramers problem [2] (2.1) is given in dimensionless phase space by

ẋ = y
ẏ = −γy − V′(x) + \sqrt{2γϵ} ẇ,  (4.53)

where V (x) is a bistable potential, γ is a dimensionless friction coefficient, and ϵ is dimensionless temperature (normalized by the barrier height). We denote the minima of V (x) by x_{A} and x_{B}, its local maximum by x_{C}, and set ω_{A(B)}^{2} ≡ V′′(x_{A(B)}), and ω_{C}^{2} ≡ −V′′(x_{C}). We assume V (x_{A}) < V (x_{B}) < V (x_{C}). The domain of attraction D_{A} of the attractor A ≡ (x_{A}, 0) in phase space is bounded by a separatrix ∂D_{A}, which passes through the saddle point C ≡ (x_{C}, 0). The direction of ∂D_{A} at C is that of the line y = −λ(x − x_{C}), where λ = [γ + \sqrt{γ^{2} + 4ω_{C}^{2}}]/2. Obviously, ϕ_{0} = exp(−E/ϵ), where E = \frac{1}{2}y^{2} + V (x). We have locally near C

χ(x, y) = \frac{ω_{C}^{2}}{\sqrt{γ(λ − γ)}}x + \sqrt{\frac{λ − γ}{γ}} y,  (4.54)

hence, according to (3.20),

λ_{1} ∼ \frac{ω_{A}}{4πω_{C}} (\sqrt{γ^{2} + 4ω_{C}^{2}} − γ) e^{−ΔV/ϵ},  (4.55)

where ΔV ≡ V (x_{C}) − V (x_{A}), which is Kramers’ result (see detailed calculations in [29], [38], and [48]).

In contrast, according to (G)TST, the rate is the normalized flux on the line x = x_{C}, y > 0,

κ = ∫_{0}^{∞} yϕ_{0} dy/ ∫_{−∞}^{x_{C}}∫_{−∞}^{∞} ϕ_{0} dx dy ∼ \frac{ω_{A}}{2π} e^{−ΔV/ϵ},  (4.56)

which is the TST result [49]. We note that if the outgoing flux on C is replaced by half the flux, calculated on any curve emanating from C to ∞, the result (4.56) is unchanged. Indeed, we have

J(ϕ_{0}) = (yϕ_{0}, −γϵϕ_{0,y} − [γy + V′(x)]ϕ_{0}) = (yϕ_{0}, −V′(x)ϕ_{0}),  (4.57)

so that

F_{0} = ∫_{C}^{∞} J(ϕ_{0}) · ν ds = ∫_{C}^{∞} ϕ_{0}(y dy + V′(x) ds) = ∫_{E_{C}}^{∞} e^{−E/ϵ} dE = ϵ e^{−ΔV/ϵ}, (4.58)

hence (4.56). We conclude that ϕ_{1} cannot be replaced by ϕ_{0} in the calculation of κ, regardless of the choice of the (G)TS.

5. Accounting for recrossings and the MFPT

A trajectory which starts in Ω_{A} and having crossed ∂D_{A} returns to Ω_{A} prior to reaching ∂D is said to recross ∂D_{A}. In the process of escaping from Ω_{A} to ∂D a trajectory may recross ∂D_{A} a random number of times. In this section we describe the escape process, calculate the mean number of times ⟨n⟩ that a trajectory crosses ∂D_{A} before it reaches ∂D, and show that for domains of type III the mean time ⟨τ(D)⟩ to reach ∂D (i.e., to escape) is asymptotically twice the MFPT ⟨τ(D_{A})⟩. We begin with a

Lemma. Let D be a domain of type III relative to (2.3). Then the MFPT from ∂D_{A} to ∂D, conditioned on reaching ∂D before ∂Ω_{A}, and the MFPT from ∂D_{A} to ∂Ω_{A}, conditioned on reaching ∂Ω_{A} before ∂D, are at most O(1/\sqrt{ϵ}) as ϵ → 0.

Proof. We prove only the second part of the lemma, since the first part follows from similar arguments. We denote by x^{∗}(t) the process x(t) in D−Ω̄_{A}, conditioned on the event {τ(∂Ω_{A}) < τ(D)}, that is, the trajectories of x^{∗}(t) consist of those trajectories of (2.1) which start in D − Ω̄_{A} and reach ∂Ω_{A} before ∂D_{A}.

The process x^{∗}(t) is a diffusion process with noise matrix σ(x) and drift

b^{∗}^{i}(x) ≡ b^{i}(x) + 2ϵ ∑_{j=1}^{n} a^{ij}(x)\frac{∂ log P (x)}{∂x^{j}}, (i = 1, …, n)  (5.1)

(see [50]), where P (x) is defined in (4.1) and is given by (4.11). We will show that the drift b^{∗}(x) is attracted to Ω_{A} everywhere in D − Ω̄_{A} and will estimate its component in the direction ”toward” Ω_{A}. To this end we change variables in the operator L_{ϵ}^{∗} in (2.9). We introduce the stretched variable

ξ ≡ \frac{χ(x)}{\sqrt{ϵ}}  (5.2)

to measure distance from ∂D_{A}. We assume that there exists a global set of variables s ≡ (s^{2}, …, s^{n}) on the hypersurfaces ξ = const. (see [39] for a more detailed discussion of this change of variables). In the variables (ξ, s) the operator L_{ϵ}^{∗} takes the form

L_{ϵ}^{∗}p = A(x)\frac{∂^{2}p}{∂ξ^{2}} + B(x)\frac{∂p}{∂ξ} + ∑_{k=2}^{n} C^{k}(x)\frac{∂p}{∂s^{k}} + O(\sqrt{ϵ}),  (5.3)

where

A(x) ≡ ∑_{i,j=1}^{n} a^{ij}(x)\frac{∂χ(x)}{∂x^{i}} \frac{∂χ(x)}{∂x^{j}} ≥ 0,  (5.4)
\sqrt{ϵ} B(x) ≡ ϵ ∑_{i,j=1}^{n} a^{ij}(x)\frac{∂^{2}χ(x)}{∂x^{i}∂x^{j}} + ∑_{i=1}^{n} b^{i}(x)\frac{∂χ(x)}{∂x^{i}},  (5.5)

and

C^{k}(x) ≡ ∑_{i=1}^{n} b^{i}(x)\frac{∂s^{k}}{∂x^{i}} .  (5.6)

The O(\sqrt{ϵ}) term in (5.3) contains mixed derivatives with respect to ζ and s^{k} and derivatives with respect to s^{k}. In view of (4.12) and (5.2),

B(x) = A(x)ξ.  (5.7)

The function P (x) in (5.1), given by (4.11), can be written in terms of the variables (ξ, s) as

P (ξ, s) = \sqrt{\frac{1}{2π}} ∫_{−∞}^{ξ} e^{−z^{2}/2} dz.  (5.8)

The components of b(x) in the ξ and s^{k} directions are denoted by (b_{ξ}, b_{s^{2}}, …, b_{s^{n}}), with

b_{ξ}(ξ, s) = A(x)ξ.  (5.9)

We define the domains D_{ξ} ≡ {ξ(x) > ξ} and note that if ξ_{1} < ξ_{2}, then D_{ξ_{1}} ⊃ D_{ξ_{2}}. In particular D_{0} = D_{A}. In view of (5.9), the drift b(x) points into D_{ξ} for ξ > 0, outside D_{ξ} for ξ < 0, and is tangent to ∂D_{0}. Equation (5.1) gives the ξ component of b^{∗}(x) as

b_{ξ}^{∗} = b_{ξ} + 2A(x)\frac{∂P (ξ, s)}{∂ξ} + 2 ∑_{k=2}^{n} A^{ξ,k}(x)\frac{∂P (ξ, s)}{∂s^{k}},  (5.10)

where A^{ξ,k}(x) is the coefficient of ∂^{2}P (ξ, s)/∂ξ∂s^{k} in (5.3). In view of (5.8) we have ∂P (ξ, s)/∂s^{k} = 0, so that

b_{ξ}^{∗}(ξ, s) = A(x) (ξ + 2\frac{∂ log P (ξ, s)}{∂ξ}) = A(x) (ξ + \frac{2e^{−ξ^{2}/2}}{∫_{−∞}^{ξ} e^{−z^{2}/2} dz} ) .  (5.11)

We denote by n(x) the inner unit normal to the surface ξ = const. Next we show that the drift b^{∗}(x) in D − Ω̄_{A} is attracted to ∂Ω_{A}. On ∂D_{A} we have ξ(x) = 0 and P (0, s) = 1/2. It follows from (5.11) that b_{ξ}^{∗}(0, s) = 2\sqrt{2/π} A(x), or in the original variables

b^{∗}(x) · n(x) = 2\sqrt{\frac{2ϵ}{π}}A(x) > 0.  (5.12)

We note that (5.12) holds in a boundary layer near ∂D_{A}. Thus the drift b^{∗}(x) on ∂D_{A} points into D_{A}. As ξ → ∞, that is, for x outside a small neighborhood of ∂D_{A} in D_{A}, the numerator in (5.11) decays exponentially, whereas the denominator converges to \sqrt{2π}. It follows that the second term in (5.1) decays exponentially so that the directions of b^{∗}(x) and b(x) coincide. Thus, by assumption, the drift b^{∗}(x) is attracted by ∂Ω_{A} in D_{A} − Ω̄_{A}. Next we consider b^{∗}(x) in D − D̄_{A}, outside a boundary layer, that is, we consider ξ → −∞. First we note that in this limit [51]

P (ξ, s) ∼ −\sqrt{\frac{1}{2π}} \frac{e^{−ξ^{2}/2}}{ξ} .  (5.13)

¿From (5.11) and (5.13) we obtain

b_{ξ}^{∗}(ξ, s) ∼ −A(x)ξ = −b_{ξ}(ξ, s),  (5.14)

where the equality is a restatement of (5.2). It follows that b^{∗}(x) points into D_{ξ} for all ξ < 0 and is therefore attracted to ∂D_{A} in D − D̄_{A}. Thus b^{∗}(x) is attracted to Ω_{A} everywhere in D − Ω̄_{A}. We see from (5.12) and (5.14) that the component b_{ξ}^{∗} is asymptotically identical to the component b_{ξ} inside D_{A} and to −b_{ξ} outside D_{A}, except for a boundary layer, where it is positive and O(\sqrt{ϵ}). That is, the conditional process x^{∗}(t) drifts toward Ω_{A} across the surfaces ∂D_{ξ} with normal speed which is either O(1) or O(\sqrt{ϵ}). Therefore the MFPT from any point in D − Ω̄_{A} to ∂Ω_{A}, conditioned on \{τ(∂Ω_{A}) < τ(D)\}, is at most O(1/\sqrt{ϵ}). This proves the lemma.

Theorem 3. If D is a domain of type III relative to (2.3), then

lim_{ϵ→0} \frac{⟨τ(D)⟩}{⟨τ(D_{A})⟩} = 2.  (5.15)

Proof. We partition the trajectories of (2.1), all of which eventually reach ∂D_{A}, into two sets,

S_{1,1} ≡ {trajectories which having reached ∂D_{A} for the first time,
reach ∂Ω_{A} before ∂D}  (5.16)

and

S_{1,2} ≡ {trajectories which having reached ∂D_{A} for the first time,
reach ∂D before ∂Ω_{A}}.  (5.17)

We denote by τ_{1} the first passage time to ∂D_{A}, that is ⟨τ_{1}⟩ = ⟨τ(D_{A})⟩. It follows from the lemma that

E[τ(D)|S_{1,2}] = ⟨τ(D_{A})⟩ + O(1/\sqrt{ϵ}) as ϵ → 0.  (5.18)

Also, as shown in Section 4,

P (S_{1,1}) ∼ P (S_{1,2}) ∼ \frac{1}{2}.  (5.19)

Next we partition S_{1,1} into two sets,

S_{2,1} ≡ {trajectories in S_{1,1} which having reached ∂Ω_{A} after the time τ_{1}, reach
∂D_{A} for the first time after τ_{1}, and then reach ∂Ω_{A} before ∂D} (5.20)

and

S_{2,2} ≡ {trajectories in S_{1,1} which having reached ∂Ω_{A} after the time τ_{1}, reach
∂D_{A} for the first time after τ_{1}, and then reach ∂D before ∂Ω_{A}}. (5.21)

Then, as above,

P (S_{2,1}) ∼ P (S_{2,2}) ∼ \frac{1}{2}P (S_{1,1}) ∼ \frac{1}{4}.  (5.22)

We denote by τ_{2} the first time after τ_{1}, that a trajectory in S_{1,1} returns to ∂D_{A} after reaching ∂Ω_{A}. By the lemma, the MFPT from ∂D_{A} to ∂Ω_{A} for trajectories in S_{2,1}, is at most O(1/\sqrt{ϵ}). The MFPT from ∂Ω_{A} to ∂D_{A} is ⟨τ(D_{A})⟩, since it is independent of the initial point of the trajectories in D_{A}, outside an ϵ-neighborhood of ∂D_{A} [32]. Thus by the lemma,

E[τ_{2} − τ_{1}|S_{2,1}] = ⟨τ(D_{A})⟩ + O(1/\sqrt{ϵ}).  (5.23)

Again, from the lemma

E[τ(D)|S_{2,2}] = 2[⟨τ(D_{A})⟩ + O(1/\sqrt{ϵ})] as ϵ → 0.  (5.24)

We proceed in a similar manner and obtain sets S_{n,1}, S_{n,2}, and times τ_{n}, such that

S_{n,1} ≡ {trajectories in S_{n−1,1} which having reached ∂Ω_{A} after the time τ_{n−1}, reach
∂D_{A} for the first time after τ_{n−1}, and then reach ∂Ω_{A} before ∂D} (5.25)

and

S_{n,2} ≡ {trajectories in S_{n−1,1} which having reached ∂Ω_{A} after the time τ_{n−1}, reach
∂D_{A} for the first time after τ_{n−1}, and then reach ∂D before ∂Ω_{A}}. (5.26)

Then, as above,

P (S_{n,1}) ∼ P (S_{n,2}) ∼ \frac{1}{2}P (S_{n−1,1}) ∼ \frac{1}{2^{n}} .  (5.27)

We denote by τ_{n} the first time after τ_{n−1} that a trajectory in S_{n−1,1} returns to ∂D_{A} after reaching ∂Ω_{A}. Again, by the lemma,

E[τ_{n} − τ_{n−1}|S_{n,1}] = ⟨τ(D_{A})⟩ + O(1/\sqrt{ϵ}),  (5.28)

and

E[τ(D)|S_{n,2}] = n[⟨τ(D_{A})⟩ + O(1/\sqrt{ϵ})], as ϵ → 0.  (5.29)

Since every trajectory reaches ∂D in finite time, having gone from ∂D_{A} to ∂Ω_{A} and back a finite number of times, the union of the sets S_{n,2}, (n = 1, 2, …) contains all trajectories. In addition, the sets S_{n,2} are mutually disjoint. It follows that

⟨τ(D)⟩ = ∑_{n=1}^{∞} E[τ_{n+1} − τ_{n}|S_{n,1}]P (S_{n,1}) = ∑_{n=1}^{∞} \frac{n}{2^{n}} [⟨τ(D_{A})⟩ + O(1/\sqrt{ϵ})]
= 2[⟨τ(D_{A})⟩ + O(1/\sqrt{ϵ})].  (5.30)

Since [32]

⟨τ(D_{A})⟩ = O(e^{Ψ̂(D_{A})/ϵ}),  (5.31)

the O(1/\sqrt{ϵ}) term in (5.30) is negligible relative to ⟨τ(D_{A})⟩, so that

⟨τ(D)⟩ ∼ 2⟨τ(D_{A})⟩,  (5.32)

hence (2.18).

Next we discuss the notion of recrossings and of the transmission coefficient k. If the transition state region (TSR) is chosen as a neighborhood of S, whose width is independent of ϵ, then, as is evident from the proof of Theorem 1, the probability that a trajectory leaves the TSR in the direction of A or B is independent of the choice of the TSR. We therefore choose the TSR to be the domain D − Ω_{A}.

Definition. A trajectory x(t) of (2.1) which starts in Ω_{A} is said to cross the TSR n times, if x(t) ∈ S_{n,2}. If a trajectory crosses the TSR n times, it is said to recross it n − 1 times.

Thus the number n of crossings of the TSR is a random variable. This definition is asymptotically independent of Ω_{A} and D (or Ω_{B}).

Corollary 1. If D is a domain of type III relative to (2.3), then the mean number of times a trajectory recrosses the TSR is 1.

Proof. From the definition and from the proof of Theorem 1 it is obvious that

⟨n⟩ = ∑_{n=1}^{∞} \frac{n}{2^{n}} = 2,  (5.33)

so that the mean number of recrossings is 1.

If the definition of the stochastic separatrix is modified to replace the probability 1/2 in (4.2) with a probability p, independent of ϵ, then according to (4.11), the resulting surface S_{p} will be in an ϵ-neighborhood of S. Then (2.18) is replaced by

τ_{AB} = \frac{1}{p}τ_{AS_{p}}.  (5.34)

It can be easily shown that τ_{AS_{p}}/p is asymptotically independent of p. Thus the convenient choice S = S_{1/2} does not limit the generality of the definition. As in the proof of Corollary 1, it can be shown that the mean number of crossings is 1/p.

Next we comment on the definition of the transmission coefficient k and its use. As mentioned in the Introduction, k was introduced to account for the recrossing trajectories which are neglected in (G)TST. If in the bottleneck case [52] the (G)TS≠ S, but passes through the saddle point, the probability of a return from the (G)TS strongly depends on the point where a trajectory hits the (G)TS. This is the source of the difficulty in calculating ⟨n⟩ in (G)TST, and of the overestimates of κ in the (G)TST results [3], [10]. If however the (G)TS is chosen as S, this difficulty is averted, as shown above.

6. Conclusions and discussion

Noise induced escapes from an attractor of a dynamical system account for a variety of physical phenomena, such as chemical reactions [2], [49], [53], nucleation [35], [54] the behavior of Josephson junction devices [55], [56], loss of lock in code tracking loops [36], [57], to name but a few. The noise in such systems describes internal or external fluctuations, brought about, e.g., by thermal vibrations, molecular collisions, noisy resistors, and so on. In many cases the strength of the fluctuations relative to the non fluctuating forces in such systems is small so that an escape is a rare event on the time scale of the dynamics. The escape rate determines physically measurable quantities, such as rates of molecular dissociation, isomerization, chemical reactions, the resistance of Josephson junction devices, the reliability of logic elements, the performance of tracking circuits, the stability of elastic structures, and so on.

In this paper we consider activation as a non equilibrium process. We define reactant and product as neighborhoods Ω_{A} and Ω_{B} of the stable states of the noiseless dynamics and show that in the limit of small noise this definition is independent of the choice of these neighborhoods. We show that the rate at which trajectories arrive from Ω_{A} to Ω_{B} for the first time is half the rate at which trajectories arrive from Ω_{A} to the stochastic separatrix S, and that this rate is also independent of these neighborhoods. We conclude that in order to obtain the correct rate in (G)TST (i) the qepdf has to be used rather than the epdf, (ii) the qepdf differs from the epdf in a boundary layer near S by a boundary layer factor, and (iii) all recrossings of the TSR are accounted for if (G)TS=S, but not otherwise. Thus no transmission coefficient is needed if (G)TS=S.

The need for (i)-(iii) in (G)TST is based on the assumption that the qepdf in multidimensional phase space is well approximated by the qepdf of the diffusion process in the sense of (1.1). We conclude that the modified (G)TST employed in [7] can be correct only if the chosen (G)TS is shown to be S, and if the modification of the epdf defined by the integral equation [7, (3.10] has the appropriate boundary layer structure.

The identification of S ∼ ∂D_{A} was based on the premise that ϵ was the only small parameter in the problem. This assumption is satisfied for example in Kramers’ model when the dissipation coefficient is not too small. In this case the above results and [38], [48] confirm Kramers’ result to be exact, not merely an upper bound. However in problems which involve more than one small parameter, then S ≠ ∂D_{A} in general. Indeed, S may differ considerably from ∂D_{A}. This is the case, e.g., in the extremely anisotropic overdamped motion of a two dimensional Brownian particle in a double well potential V (x, y). Here x and y are the coordinates along which the motion is fast and slow, respectively. If ∂^{2}V/∂x^{2} < 0 at the saddle point, then S may be completely unrelated to ∂D_{A}. This leads to many interesting phenomena such as the appearance of long non exponential transients, decay rates not determined by the height of the potential barrier, etc. [23], [58]-[60]. Also in the limit of small damping S ≠ ∂D_{A} [22]. In such cases the the factor 1/2 in the relation between the normalized flux on the GTS and κ should be used only if GTS=S.

In [47] an attempt was made to prove that lim_{ϵ→0} S = ∂D_{A}, for continuous drifts b(x), by employing the following argument. Using (4.11)-(4.15), in the limit ϵ → 0, χ^{0}(x) satisfies

∑_{i=1}^{n} b^{i}(x)\frac{∂χ^{0}(x)}{∂x^{i}} = 0 on S (on χ^{0}(x) = 0).  (6.1)

Hence lim_{ϵ→0} S is a characteristic surface, tangent to the vector field b(x). Since ∂D_{A} is a characteristic surface, it is concluded in [47] that lim_{ϵ→0} S = ∂D_{A}. However ∂D_{A} is not necessarily the only closed characteristic surface in D, even if there is only one attractor in D, and the above argument fails to choose S among the various closed characteristic surfaces. Therefore this argument is not complete. For example, the two-dimensional dynamics (in polar coordinates)

ṙ = −r(1 − r)^{2}(2 − r)
θ̇ = f(r, θ),  (6.2)

where f(r, θ) is a positive function, has a single attractor in the disk r < 2, centered at r = 0, and an unstable limit cycle at r = 1 (attracting from the outside and repelling from the inside). The second attractor is at infinity. Both circles r = 1 and r = 2 are closed characteristic surfaces, and the argument (6.1) does not identify which circle is S in the limit ϵ → 0. In fact, our approach shows that in the limit ϵ → 0, S is the circle r = 2, as follows from the argument of Section 4. Indeed, the outer solution to (4.3) is constant inside the disk r < 1. It is also constant in the annulus 1 < r < 2, since all trajectories in the annulus come arbitrarily close together near the attractor r = 1. It is easy to see that the solution of the boundary layer equation near r = 1 cannot match to two different constants on the two sides of r = 1. It follows that the outer solution is a single constant throughout the disk r < 2, so that lim_{ϵ→0} S = {r = 2}.

Acknowledgment. The authors thank E. Pollak for useful discussions of TST.

FIGURE CAPTION

Figure 1. The dynamics (2.3) is attracted to A and to ∂D and is repelled from ∂D_{A}. The domains Ω_{A}, D_{A}, and D are of types I, II, and III, respectively.

References

[1] P. Pechukas, in Dynamics of Molecular Collisions B, (W.H. Miller, Ed.), Plenum, N.Y. 1976.

[2] H.A. Kramers, Physica (Utrecht) 7, 284 (1940).

[3] G.R. Fleming and P.G. Wolynes, Phys. Today 43, 36 (1990).

[4] E. Wigner, Trans. Faraday Soc. 34, 29 (1938).

[5] M. Büttiker, Harris, and R. Landauer

[6] Melnikov and Meshkov

[7] E. Pollak, H. Grabert, and P. Hänngi, J. Chem. Phys. 91, 4073 (1989).

[8] H. Eyring, J. Chem. Phys. 3, 107, 1935.

[9] W.H. Miller, J. Chem. Phys. 61, 1283 (1974).

[10] H. Frauenfelder and P.G. Wolynes, Science 229, 337 (1985).

[11] M.G. Evans, Trans. Faraday Soc. 34, 49 (1938).

[12] J.C. Keck, Adv. Chem. Phys. 13, 85 (1967).

[13] E. Pollak, J. Chem. Phys. 85, 865 (1986).

[14] E. Pollak, J. Chem. Phys. 93, 1116 (1990).

[15] J.A. Christiansen, Z. Phys. Chem. B33, 145 (1936).

[16] R. Kubo, J. Phys. Soc. Japan 12, 570 (1957), Rep. Prog. Phys. 29, 255 (1966).

[17] G.W. Ford, M. Kac, and P. Mazur, J. Math. Phys. 6, 504 (1965).

[18] R. Zwanzig, J. Stat. Phys. 9, 215 (1973).

[19] S. Chandrasekhar, Rev. Mod. Phys. 15, 1 (1943).

[20] E. Pollak, S.C. Tucker, and B.J. Berne, Phys. Rev. Lett. 65, 1399 (1990).

[21] M.M. Dygas, B.J. Matkowsky, and Z. Schuss, SIAM J. Appl. Math. 46 (2), 265 (1986).

[22] M.M. Kłosek, B.J. Matkowsky and Z. Schuss, Berichte Der Bunsen-Gesellschaft Für Physikalische Chemie (1991).

[23] M.M. Kłosek, B.J. Matkowsky, A. Nitzan, M. Ratner, Z. Schuss, (in preparation).

[24] M.I. Freidlin and A.D. Wentzell, Random Perturbations of Dynamical Systems, Springer-Verlag, N.Y. 1984.

[25] D. Ludwig, SIAM Rev. 17 (4), 605, (1975).

[26] Z. Schuss and B.J. Matkowsky, SIAM J. Appl. Math. 36, 604 (1979).

[27] M. Mangel, D. Ludwig, SIAM J. Appl. Math. 33 (2), 256 (1977).

[28] M. Mangel, SIAM J. Appl. Math. 36 (3), 544 (1979).

[29] T. Naeh, M.M. Kłosek, B. J. Matkowsky, and Z. Schuss, SIAM J. Appl. Math. 50, 595 (1990).

[30] H. Risken, The Fokker-Planck Equation, Springer-Verlag, N.Y. 1984.

[31] B. Cartling, J. Chem. Phys. 87, 2638 (1987).

[32] Z. Schuss, Theory and Applications of Stochastic Differential Equations, J. Wiley & Sons, N.Y. 1980. (1975).

[33] H.C. Brinkman, Physica (Utrecht) 22, 149 (1956).

[34] R. Landauer and J.A. Swanson, Phys. Rev. 121, 1668 (1961).

[35] J.S. Langer, Ann. Phys. (N.Y.) 54, 258 (1969).

[36] B.Z. Bobrovsky and Z. Schuss, SIAM J. Appl. Math. 42, 174 (1982).

[37] B.J. Matkowsky and Z. Schuss, SIAM J. Appl. Math. 42, 822 (1982).

[38] B.J. Matkowsky, Z. Schuss, and C. Tier, J. Stat. Phys. 35 (3,4), 443 (1984).

[39] M.V. Day, SIAM J. Math. Anal. 20, 222 (1989), J. Math. Anal. Appl. 147, 134 (1990), preprint (1991).

[40] I. Oppenheim, K.E. Shuler, and G.H. Weiss, Stochastic Processes in Chemical Physics: The Master Equation, The MIT Press, Cambridge MA 1977.

[41] C.W. Gardiner, Handbook of Stochastic Methods, Springer Verlag, N.Y. 1982.

[42] M.H. Protter and H.F. Weinberger, Maximum Principles in Differential Equations, Prentice-Hall, Englewood Cliffs 1967.

[43] B.J. Matkowsky and Z. Schuss, SIAM J. Appl. Math. 33, 365 (1977).

[44] R. Graham and T. Tel, Phys. Rev. Lett. 52, 9 (1984) (see also Phys. Rev. A 33, 1322 (1986)).

[45] C.R. Doering, P.S. Hagan, C.D. Levermore, Phys. Rev. Lett. 59, 2129 (1987), A.W. Yu, G.P. Agrawal, R. Roy, J.Stat. Phys. 54, 1223 (1989).

[46] D. Ryter, J. Stat. Phys. 49, 751 (1987).

[47] D. Ryter, Physica 142A, 103 (1987).

[48] B.J. Matkowsky, Z. Schuss, and E. Ben-Jacob, SIAM J. Appl. Math. 42, 835 (1982).

[49] S. Glasstone, K.J. Laidler, H. Eyring, The Theory of Rate Processes, McGraw-Hill, N.Y. 1941.

[50] S. Karlin and H.M. Taylor, A Second Course in Stochastic Processes, Academic Press, N.Y. 1981.

[51] W. Feller, An Introduction to Probability Theory and its Applications I, Third Edition, J. Wiley & Sons, N.Y. 1968.

[52] M.M. Kreevoy and D.G. Truhlar, in Investigation of Rates and Mechanisms of Reactions, vol. 6, 4th ed. (C.F. Bernasconi, Ed.), J. Wiley, N.Y. 1986.

[53] S.H. Benson, The Foundations of Chemical Kinetics, McGraw-Hill, N.Y. 1960.

[54] H. Haken, Rev. Mod. Phys. 47, 67 (1975).

[55] E. Ben-Jacob, D.J. Bergman, B.J. Matkowsky, Z. Schuss, Phys. Rev. A, 26, 2805 (1982).

[56] E. Ben-Jacob, D.J. Bergman, Y. Imry, B.J. Matkowsky, Z. Schuss, J. Appl. Phys. 54, 6533 (1983).

[57] A.J. Viterbi, Principles of Coherent Communications, McGraw-Hill, N.Y. 1966.

[58] A.M. Berezhkovskii and V.Y. Zitserman, Physica 166A, 585 (1990).

[59] A.M. Berezhkovskii and V.Y. Zitserman, Chem. Phys. Lett. 172, 235 (1990).

[60] M.M. Kłosek, B.J. Matkowsky, A. Nitzan, M. Ratner, Z. Schuss, J. Chem. Phys. 95 (2) (1991).