Master equations and path-integral formulation of variational principles for reactions

GAVEAU, Bernard UFR de Mathématique 925 Université Pierre & Marie Curie 75252 Paris Cedex 05, FRANCE e-mail: gaveau@ccr.jussieu.fr

MOREAU, Michel Laboratoire des Physique Theoretique des Liquides Université Pierre & Marie Curie 75252 Paris Cedex 05, FRANCE e-mail: moreau@lptl.jussieu.fr

TÓTH, János Department of Analysis, Institute of Mathematics Budapest University of Technology and Economics Budapest, Egry J. u. 1. H-1521 HUNGARY e-mail: jtoth@math.bme.hu

April 2004

Abstract

The mesoscopic nonequlibrium thermodynamics of a reaction-diffusion system is described by the master equation. The information potential is defined as the logarithm of the stationary distribution. The Fokker–Planck approximation and the Wentzel–Kramers–Brillouin method give very different results. The information potential is shown to obey a Hamilton–Jacobi equation, and from this fact general properties of this potential are derived. The Hamilton–Jacobi equation is shown to have a unique regular solution. Using the path integral formulation of the Hamilton–Jacobi approximation of the master equation it is possible to calculate rate constants for the transition from one well to another one of the information potential and give estimates of mean exit times. In progress variables, the Hamilton– Jacobi equation has always a simple solution which is a state function if and only if there exists a thermodynamic equilibrium for the system. An inequality between energy and information dissipation is studied, and the notion of relative entropy is investigated. A specific two-variable system and systems with a single chemical species are investigated in detail, where all the defined relevant quantities can be calculated explicitly.

1 Reaction-diffusion systems

A uniform reactive system can be represented by a system of volume V at temperature T containing several reactive variables (or internal species) X_{m} (m = 1, 2, …, M) possibly in an inert solvent X_{0}. The species can react according to an arbitrary number of reversible reaction steps, labelled by r which can be represented in the form:

∑_{m=1}^{M} α_{r}^{m}X_{m} ⟺ ∑_{m=1}^{M} β_{r}^{m}X_{m} (r = 1, 2, …, R),  (1)

where the stoichiometric coefficients α_{r}^{m} and β_{r}^{m} can be nonnegative integers. This formulation can be extended to a nonhomogeneous system, if one assumes [5, page 169],[1], that the system is divided into cells small enough to be approximately homogeneous, exchanging molecules by linear rate. The molecules of a given species receive different labels in different cells, and diffusion is represented by the pseudo-chemical reaction X_{m} ⟺ X_{q}, X_{m} and X_{q} representing the same species in neighboring cells.

(1) however can only occur between molecules within the same cell.

From a stochastic point of view, the system can be described by the probability P(X_{1}, X_{2}, …, X_{M}) = P(X, t) of finding X_{m} particles of the mth species (m = 1, 2, …, M) at time t, where X_{m} denotes the number of molecules of species X_{m}.

We assume that the vector X = (X_{1}, X_{2}, …, X_{M}) changes according to a Markovian jump process.

2 Master equation and Fokker–Planck equation

2.1 Master equation

The system described above evolves by various subprocesses, and we call Wγ(X) the probability per unit time of the transition X ⟶ X + γ with given vectors of nonnegative integers γ. These vectors may be any of the elementary reaction vectors β_{r} − α_{r}, α_{r} − β_{r} (r = 1, 2, …, R). The state of the system at time t is described by the absolute probability distribution function P(X, t) whose time evolution is given by the master equation

\frac{∂P (X, t)}{∂t} = (ΛP)(X, t)  (2)

where Λ is the evolution operator [26, 27, 29, 39, 40, 41]

(ΛP)(X, t) := ∑ ; γ[Wγ(X − γ)P(X − γ, t) − Wγ(X)P(X, t)]  (3)

To calculate the large volume limit we introduce the concentration vector by x := X/V, and define the probability density function p through

P(B, t) = ∫_{B} p(x, t)dx  (4)

for all the measurable subsets B of the state space. Instead of the overall transition rates Wγ we can introduce the concentration dependent transition rates by wγ(x) := Wγ(X)/V. Using these notations the master equation can be rewritten for the function p as

\frac{∂p(x, t)}{∂t} = V ∑ ; γ[wγ(x − \frac{γ}{V} )p(x − \frac{γ}{V} , t) − wγ(x)p(x, t)].  (5)

This formalism can also describe reaction-diffusion processes approximately using the trick mentioned above.

2.2 Approximate Fokker–Planck equation

A standard approximation [41] of the master equation in the large volume limit is obtained by the Taylor expansion of the right hand side of Eq. (5) up to the terms of order O(1/V). This expansion gives the Fokker–Planck equation for p:

\frac{∂p(x, t)}{∂t} = − ∑_{m=1}^{M} \frac{∂}{∂x_{m}} (A_{m}p)(x, t) + \frac{1}{2V} ∑_{m,q=1}^{M} \frac{∂^{2}}{∂x_{m}∂x_{q}} (D_{mq}p)(x, t), (6)

where

A_{m}(x) := ∑ ; γ γ_{m}wγ(x), D_{mq}(x) := ∑ ; γ γ_{m}γ_{q}wγ(x) (m = 1, 2, …, M). (7)

The first quantity may obviously be called conditional expected velocity, and the second one conditional expected variance velocity [6]. It is well-known that keeping terms of higher order may lead to inconsistent results. Moreover, the present approximation does not respect the natural boundary conditions of the master equation. Both the master equation and the Fokker–Planck equation are approximately consistent in mean (and they are exactly consistent, if only processes with linear rates are present) with the usual deterministic equation

\frac{d\overline{x}}{dt} = A(\overline{x}(t))  (8)

for the average \overline{x} of x. The exact equation for the moments of the distribution obeying the master equation can also be derived from the equation for the generating function, see [5, p. 110] and [46, p. 129–133]. E.g. instead of (8) we have

\frac{d\overline{x}}{dt} = \overline{A(x(t))}.  (9)

2.3 Conservation laws and irreducibility

During a time interval [t, t + h] the change of the concentrations is

x(t + h) − x(t) = ∑ ; γ γηγ,  (10)

where the ηγ are (integer valued) random variables indexed corresponding to nonzero transition rates wγ. The stochastic process t ↦ x(t) will then satisfy some conservation laws. More specifically, for each point x_{0} of the concentration space ℝ^{M} let us define the stoichiometric compatibility class [5, 25]

S(x_{0}) := {x ∈ ℝ^{M}; x = x_{0} + ∑ ; γ γηγ; ηγ ∈ ℝ} ∩ (ℝ_{0}^{+})^{M},  (11)

where ξγ are integers corresponding to nonzero transition rates wγ. The whole space ℝ^{M} is a union of stoichiometric compatibility classes, and the process t ↦ x(t) stays in S(x_{0}) for all time (until the process is defined), if x(0) = x_{0}. We also suppose that each stoichiometric compatibility class contains exactly one stationary solution (which is not always so even in the deterministic case). Furthermore, we suppose that the stochastic model has a unique stationary solution for each x_{0}. This stationary distribution will be denoted by p(x|x_{0}), and it does not depend on the choice of the initial state so far as it remains in S(x_{0}). The master equation (5) induces another master equation in a smaller number of linearly independent variables on each S(x_{0}) (the other variables are linear functions of these variables), but the stationary distribution is supposed to be unique. The deterministic motion also evolves in S(x_{0}) what can easily seen by rewriting Eq. (8) into an integral equation form:

\overline{x}(t) = x_{0} + ∫_{0}^{t} ∑ ; γ γwγ(\overline{x}(s))ds = x_{0} + ∑ ; γ γ ∫_{0}^{t} wγ(\overline{x}(s))ds.  (12)

Similar remarks are also valid for the stochastic process associated to the Fokker– Planck equation. From now on it will be assumed that there are no linear conservation laws, a fact what can also be expressed that the master equation is irreducible. This assumption does not mean the restriction of generality because otherwise one can consider the master equation restricted to a stoichiometric compatibility class of the state space. We shall also assume that the zeros of the vector field (of the deterministic model) are isolated (although sometimes we shall consider the multiple zero case) on each set S(x_{0}).

3 Hamilton–Jacobi theories

Let us recall an approximation originally introduced by Kubo et al. ([28]; see also [34], and more recently [51, 52, 53]), giving physically more realistic results (see Section 4 and [11, 12]) than the usual Fokker–Planck approximation. Let us also mention that our basic reference is [36].

3.1 Hamilton–Jacobi theory for the master equation

The idea of the approximation is to write p(x, t) in the form of Wentzel–Kramers– Brillouin (WKB) expansion, valid for large V

p(x, t) = exp[−V Φ(x, t)][U_{0}(x, t) + \frac{1}{V} U_{1}(x, t) + …],  (13)

where Φ, U_{0}, U_{1} are unknown functions. Then, one uses the first term of this expansion in Eq. (5) and groups together terms containing decreasing powers of V. The highest order term in V is of the order O(V), and contains U_{0} as a factor. It is zero if and only if Φ satisfies the equation

\frac{∂Φ}{∂t} + ∑ ; γ wγ [exp (∇Φ · γ) − 1] = 0.  (14)

This is a Hamilton–Jacobi equation which can be integrated by the method of bicharacteristics [2, 32]. The next term of order O(1) in V in Eq. (5) when one uses the expansion (13) gives a first order linear equation for U_{0} which is a sort of transport equation:

\frac{∂U_{0}}{∂t} + ∑ ; γ [exp (∇Φ · γ) ∑_{m=1}^{M} \frac{∂}{∂x_{m}} (γ_{m}wγU_{0})]

+U_{0} ∑ ; γ exp (∇Φ · γ) (\frac{1}{2} ∑_{m,q=1}^{M} γ_{m}γ_{q}wγ \frac{∂^{2}Φ}{∂x_{m}∂x_{q}} ) = 0.  (15)

3.2 Hamilton–Jacobi theory for the Fokker–Planck equation

For large V, one can use the first term of the WKB-type expansion (here and above we leave including more terms for future investigations) for p in Eq. (6), i.e. to look for p in the form

p(x, t) = exp(−V Φ^{(FP)}(x, t))[U_{0}^{(FP)}(x, t) + \frac{1}{V} U_{1}^{(FP)}(x, t) + …].  (16)

Upon inserting the above expansion (16) into Eq. (6) and comparing the various powers of V one obtains equations for Φ^{(FP)} and U_{0}^{(FP)}:

\frac{∂Φ^{(FP)}}{∂t} + ∑_{m=1}^{M} A_{m} \frac{∂Φ^{(FP)}}{∂x_{m}} + \frac{1}{2} ∑_{m,q=1}^{M} D_{mq} \frac{∂Φ^{(FP)}}{∂x_{m}} \frac{∂Φ^{(FP)}}{∂x_{q}} = 0 (17)

\frac{∂U_{0}^{(FP)}}{∂t} + ∑_{m=1}^{M} \frac{∂}{∂x_{m}} (A_{m}U_{0}^{(FP)}) + ∑_{m,q=1}^{M} \frac{∂Φ^{(FP)}}{∂x_{m}} \frac{∂}{∂x_{q}} (D_{mq}U_{0}^{(FP)})

+\frac{1}{2} ( ∑_{m,q=1}^{M} D_{mq} \frac{∂^{2}Φ^{(FP)}}{∂x_{m}∂x_{q}} ) U_{0}^{(FP)} = 0. (18)

Eq. (17) is a standard Hamilton–Jacobi equation (with a standard Hamiltonian, quadratic in the momentum), whereas Eq. (18) is a kind of transport equation. Eqs.

(17) and (18) can be seen to be obtained from Eqs. (14) and (15) if one assumes that the derivatives ∂Φ/∂x_{m} are small, so that one is entitled to replace exp (∇Φ · γ)−1 by its Taylor expansion up to the second order. Then, Eq. (14) reduces immediately to (17) with A_{m} and D_{mq} as given in Eq. (7). This indicates that the Hamilton–Jacobi and transport equations (17)–(18) are less precise than the Hamilton–Jacobi and transport equations (14)–(15) directly deduced from the master equation. In fact, their range of validity is limited to neighborhoods of the stationary points (where ∇Φ is zero) of the action function Φ.

3.3 Stationary solutions

The stationary solution of the master equation or that of the Fokker–Planck equation will also be looked for in the form p(x) ∼ exp[−V Φ(x)][U_{0}(x) + \frac{1}{V} U_{1}(x) + …]. Then Φ satisfies the stationary form of the time dependent Hamilton–Jacobi equations (14) or (17) and U_{0} satisfies the stationary form of the transport equations (Eqs. (15) or (18)). In both cases, the stationary Hamilton–Jacobi equation can be written as

H(x, ∇Φ(x)) = 0,  (19)

where H is of the form

H^{(M)}(x, ξ) = ∑ ; γ wγ(x)[exp(ξ^{⊤}γ) − 1]  (20)

or

H^{(FP)}(x, ξ) = ξ^{⊤}A(x) + \frac{1}{2}ξ^{⊤}D(x)ξ,  (21)

for the master equation and for the Fokker–Planck equation, respectively, the variables ξ_{m} being the conjugate momenta of x_{m}. From Eq. (7) we deduce for small ξ

H^{(M)}(x, ξ) = H^{(FP)}(x, ξ) + O(|ξ|^{3}).  (22)

In the present work we shall only use the stationary equations, although there is some hope to extend these investigations to nonstationary cases. In some neighborhood of some stationary point of the vector field A it is always possible to study the stationary solution of the master equation by using the Hamiltonian originated in the Fokker– Planck approximation (Fokker–Planck Hamiltonian, see [14, pp. 7742–7743]).

3.4 The particular case of one chemical species

Let us consider the case of a single chemical species, i.e. when M = 1, and let us suppose that the system evolves only by transitions n ⟶ n ± 1. As before we go to the large volume limit introducing x := \frac{X}{V} , p(x) := V P(X), w_{±}(x) := \frac{1}{V} W_{±}(X), so that one obtains the exact master equation (5) and the Fokker–Planck equation (6) respectively, with A(x) = w_{+}(x) − w_{−}(x), D(x) = w_{+}(x) + w_{−}(x). The solution (13) of the master equation is ([11, 12, 34])

p(x) = U_{0}(x) exp[−V Φ(x|a)] = \frac{C}{\sqrt{w_{+}(x)w_{−}(x)}} exp (−V ∫_{a}^{x} log (\frac{w_{−}}{w_{+}} )) ,

where a is an arbitrary nonnegative value and C is a normalization constant. When w_{+} and w_{−} have no common zeros, p as given in (13) can be normalized with

C = [∫_{0}^{+∞} \frac{exp(−V ∫_{a}^{x} log (\frac{w_{−}}{w_{+}} ))}{\sqrt{w_{+}(x)w_{−}(x)}}dx]^{−1} .  (23)

The integral in (23) can be estimated by the saddle point method. The main contribution to the integral is obtained for a point x_{s} which is an absolute minimum of Φ(x|a). If exactly one such point exists, we have

p(x) ∼ \frac{C_{0}}{\sqrt{w_{+}(x)w_{−}(x)}} exp (−V ∫_{a}^{x} log (\frac{w_{−}}{w_{+}} )) ,  (24)

where C_{0} is now of the order O(1/\sqrt{V}), if the minimum of Φ is not degenerate. When there are several attracting points of A, say, x_{s}^{1}, x_{s}^{2}, …, x_{s}^{l}, …, one chooses x_{s}^{k} in such a way that for all l Φ(x_{s}^{(l)}|x_{s}^{(k)}) ≤ 0, and again one obtains (24).

Similarly, the stationary solution of the Fokker–Planck equation can be written as

p(x) ∼ U_{0}^{(FP)}(x) exp[−V Φ^{(FP)}(x|a)] = \frac{C}{D(x)} exp (2V ∫_{a}^{x} \frac{A}{D} ) , (25)

and is normalizable, provided D has no zero (or w_{+} and w_{−} have no common zero),

with C = [∫_{0}^{+∞} ^{\frac{exp(2V ∫_{a}^{x} \frac{A}{D})}{D(x)}}dx]^{−1} . Again, C can be evaluated by the saddle point

method, the main contribution coming form the absolute minimum x_{s} of Φ^{(FP)}(a), which is an attracting point of the vector field A: p(x) ∼ \frac{C_{0}}{\sqrt{V} D(x_{s})} exp (2V ∫_{x_{s}}^{x} \frac{A}{D} ) .

When w_{+} and w_{−} have a common zero, p given above is not normalizable, indicating that the expression in (13) for p is not valid. This is exactly the case of criticality [17].

3.5 Comparison of asymptotic results

The two approximations given by (24) and (25) are close to each other for small values of dΦ/dx, or when A(x) is small. In this case one has log(\frac{w_{−}}{w_{+}} ) ∼ −\frac{2A}{D}. The approximation (25) coming from the Fokker–Planck equation is exact near a zero of A. But when there are several zeros, the approximation fails, because the eigenvalues and the mean exit times calculated from the Fokker–Planck equation differ from the analogous quantities given by the master equation by an exponentially large factor [11, 12].

This indicates also that the limit theorems by Kurtz [31] are not applicable in the case when A has several zeros: they have been formulated for the case of detailed balanced systems which are known to have a single asymptotically stable deterministic stationary point. These theorems state that the stochastic process x associated with the master equation tend to the deterministic trajectory of A(x) (starting from the same initial point), and that the deviation is Gaussian, but these theorems are valid uniformly on finite time intervals. They cannot describe the situation for times like exp(kV). In particular, these theorems cannot describe chemically activated events like the passage over a potential barrier, and consider rate constants as given.

4 Construction, uniqueness and critical points of Φ

General basic properties of solution Φ to the Hamilton–Jacobi equation (19) given by (20) or (21) will be derived here and in the next section. The traditional method [2, 32] for constructing Φ does not work, but we shall still show a construction method and also prove the uniqueness of the smooth solution Φ. Finally, we also study the critical points of Φ. We shall treat H^{(M)} and H^{(FP)} together.

4.1 Lagrangians

With the Hamiltonian for the master equation and for the Fokker–Planck equation (20) and (21), respectively, the corresponding ẋ_{m} velocities are

ẋ_{m}(t) = \frac{∂H^{(M)}(x(t), ξ)}{∂ξ_{m}} = ∑ ; γ γ_{m}wγ(x(t)) exp(ξ^{⊤}γ),

ẋ_{m}(t) = \frac{∂H^{(FP)}(x(t), ξ)}{∂ξ_{m}} = A_{m}(x) + ∑_{q=1}^{M} D_{mq}(x(t))ξ_{q},

and the corresponding Lagrangians are

L^{(M)}(x, v) := ∑ ; γ wγ(x)[(ξ^{⊤}γ) exp(ξ^{⊤}γ) − exp(ξ^{⊤}γ) + 1],  (26)

L^{(FP)}(x, v) := \frac{1}{2}(v − A(x))^{⊤}D(x)^{−1}(v − A(x)).  (27)

In the Fokker–Planck case (27) ξ = D(x)^{−1}(v − A(x)) (with v = ẋ as usual), provided D(x) is nondegenerate, which is the case if we assume that there are no conservation laws. Under the hypothesis of the nonexistence of linear mass conservation laws L^{(M)} and L^{(FP)} are both nonnegative and both of them are zero if and only if ξ = 0. In the presence of conservation laws D(x) is necessarily degenerate.

Actually, (27) is the Gaussian form of the Onsager-Machlup principle in our case and a special (discrete) version of the Gyarmati principle for nonlinear conduction equations. See e.g. [4, 22, 23, 47].

4.2 Special paths

4.2.1 Deterministic paths

The deterministic paths \frac{d\overline{x}}{dt} = A(\overline{x}(t)) x(0) = 0 are obviously solutions of both Hamiltonian equations. Conversely, a path t ↦ (x(t), ξ(t)) which is a solution of the Hamiltonian equations, such that x(0) = 0, is the deterministic path, because of the uniqueness of paths under given initial conditions. Moreover, the Lagrangian is zero along a deterministic path, and conversely, and as a consequence, the variation of the action Φ along a deterministic path is zero.

4.2.2 Antideterministic paths

Let us assume now that Φ is a smooth solution of the Hamilton–Jacobi equation, and define ξ_{m}(x) := \frac{∂Φ}{∂x_{m}} . A solution s ↦ (x(s), ξ(s)) of the system

\frac{dx_{m}}{ds} = \frac{∂H(x, ξ)}{∂ξ_{m}} ξ_{m}(s) = ξ_{m}(x(s))  (28)

is a Hamiltonian path because

\frac{dξ_{m}(s)}{ds} = ∑ \frac{∂ξ_{m}}{∂x_{q}} \frac{dx_{q}}{ds} = ∑_{q} \frac{∂^{2}Φ}{∂x_{m}∂x_{q}} \frac{∂H(x, ∇Φ(x))}{∂ξ_{q}} = −\frac{∂H(x, ξ(x))}{∂x_{m}},

since H(x, ∇Φ(x)) = 0, therefore \frac{∂H}{∂x_{m}} = ∑_{q} \frac{∂H}{∂ξ_{q}} \frac{∂^{2}Φ}{∂x_{m}∂x_{q}} = 0. Moreover, Φ is increasing along such the paths because

\frac{dΦ}{ds} = ∑_{m} \frac{∂Φ}{∂x_{m}} \frac{dx_{m}}{ds} = ∑_{m} ξ_{m}(s)\frac{dx_{m}}{ds} = L ≥ 0.

The trajectories given by (28) will be called antideterministic paths (for reasons to follow).

4.3 Construction of Φ

The traditional method to construct a solution Φ of (19) (the action) is the following [2, 32]. One chooses a point x^{(0)} and considers a path (x(s), ξ(s)), a solution of the Hamiltonian system

ẋ_{m} = \frac{∂H}{∂ξ_{m}} ξ̇_{m} = −\frac{∂H}{∂x_{m}} ,  (29)

with the conditions

x(0) = x^{(0)}, ξ(0) = ξ^{(0)} x(t) = x H(x^{(0)}, ξ^{(0)}) = 0.  (30)

The unknowns are ξ^{(0)} and t, which will be implicitly fixed by condition (30). Then the function

Φ(x) := ∫_{0}^{t} ∑_{m=1}^{M} ξ_{m}dx_{m}  (31)

is the solution Φ(x|x^{(0)}) of (19) where in Eq. (31) the integral is taken along the path (x(s), ξ(s)), satisfying Eqs. (29)–(30). But the function Φ(x|x^{(0)}) is not differentiable at x^{(0)}, since ∇Φ(x|x^{(0)}) = ξ → ξ_{0} as x → x_{0}; but ξ_{0} depends on the path from x to x_{0}, so that ∇Φ is, in general, not defined at x_{0},. In our situation where p(x) ∼ U_{0}(x) exp[−V Φ(x)], we expect Φ(x) to be regular everywhere and to be peaked at a point x_{M} (at least), so that one cannot take for our Φ a function Φ(x|x_{0}) constructed by the traditional method as above.

We shall now describe the correct construction of Φ by a limiting process.

Let us consider a point x_{s} which is an attracting point of the vector field A. We take another point x^{(0)} and we construct the usual action Φ(x|x^{(0)}) using Hamiltonian paths starting from x^{(0)}, with energy 0. For x^{(0)} ≠ x_{s} this function is nontrivial and is not differentiable at x^{(0)}. The function Φ(x|x_{s}) is now defined as

Φ(x|x_{s}) := lim_{x^{(0)}→x_{s}} Φ(x|x^{(0)}).  (32)

The function Φ(x|x_{s}) is not simply the function Φ(x|x^{(0)}) taken at x^{(0)} = x_{s}. The reason is that if we choose x^{(0)} = x_{s}, ξ_{0} should vanish in order that H(x^{(0)}, ξ^{(0)}) = 0. This is obvious for H^{(FP)}, and can be seen to be also valid for H^{(M)} using the inequality e^{a} − 1 ≥ a. Therefore the trajectory never moves away from x^{(0)}, and the traditional action is 0. In general, Φ(x|x_{s}) is a nontrivial function, which is differentiable at x^{(0)} = x_{s}, has a strict minimum at x_{s}, which is a nondegenerate minimum, if x_{s} is a nondegenerate attracting point of the vector field A. The existence of the limit in (32) has been proven in [14, p. 7743–7745]. Moreover, we shall consider a given point x and a trajectory with energy 0 starting from x^{(0)} and arriving at x in a certain (unknown) time t. This trajectory has an initial momentum ξ which is a function ξ(x|x^{(0)}). In the Fokker–Planck case we have H^{(FP)}[x^{(0)}, ξ(x|x^{(0)})] = 0. If δ := x^{(0)} − x_{s}, then A(x^{(0)}) ∼ A(x_{s}) + A′(x_{s})δ ≃ A′(x_{s})δ, if δ → 0, then from the definition of H^{(FP)} we have

\frac{1}{2}ξ^{⊤}Dξ ≃ −A′(x_{s})δ · ξ.  (33)

Assuming that A′(x_{s}) ≠ 0 and D(x_{s}) is invertible we deduce from (33) that

|ξ(x|x^{(0)})| = O(δ),  (34)

showing that ξ(x|x^{(0)}) → 0 when x^{(0)} → x_{s}. From Eq. (34) and the definition of the velocity we get ẋ = A(x) + D(x)ξ. As we see that for the initial velocity |ẋ(0)| = O(δ) holds, the time needed to join x^{(0)} along the Hamiltonian trajectory of energy 0 will tend to infinity when x^{(0)} → x_{s}. It is precisely because the initial momentum is tending to zero that the limiting function Φ(x|x^{(0)}) will be differentiable at x_{s}, when x^{(0)} → x_{s}.

4.4 Uniqueness of Φ

If Φ is a smooth solution of the Hamilton–Jacobi equation (either H^{(M)} or H^{(FP)}), and x_{0} is a minimum of Φ, then the Taylor expansion of Φ at x_{0} is uniquely determined up to an additive constant [14, p. 7741]. In particular, if Φ is an analytic solution near x_{0}, it is unique. As a consequence, there exists at most one function Φ which is a global analytic solution of the Hamilton–Jacobi equation (up to an additive constant). Because of this fact we can determine the antideterministic path by using the analytic solution Φ of the Hamilton–Jacobi equation as in Paragraph 4.2.2, namely,

\frac{dx̃_{m}(s)}{ds} = \frac{∂H[x̃(s), ξ̃(s)]}{∂ξ_{m}}, ξ̃_{m}(s) = ξ_{m}[x̃(s)] ≡ \frac{∂Φ(x̃(s))}{∂x_{m}}.  (35)

4.5 Limiting behavior of trajectories

Consider the linear Fokker–Planck Hamiltonian with a constant positive matrix D and with a linear vector field A(x) := A^{0}x. Let (x(s), ξ(s)) be a trajectory with H = 0 such that x(0) = x^{(0)}, x(t) = x, where x^{(0)} and x are both nonzero. Then, if t → +∞, this trajectory has the following limit behavior:

1. for fixed s x(s) tends to the deterministic trajectory \overline{x}(s) starting from x^{(0)};

2. for fixed s x(t − s) tends to the antideterministic trajectory x̃(t − s) defined by Eq. (35) ending at x.

4.6 Critical points of Φ and zeros of the deterministic model (vector field)

The following facts are true [14, p. 7745]:

1. A nondegenerate critical point of Φ is a zero of the deterministic vector field A (both for H^{(M)} and H^{(FP)}). A nondegenerate minimum of Φ is a stable attracting point A.

2. Conversely, for H^{(FP)} the zeroes of A are critical points of Φ, and the stable attracting points of A are minima of Φ.

3. Conversely, for H^{(M)} the stable attracting points of A are minima of Φ.

5 Path integrals for the master equation and Fokker–Planck equation

5.1 The stochastic process of the Fokker–Planck equation

It is well known that there exists a stochastic process t ↦ x(t) associated to the Fokker–Planck equation (6). The weight of a trajectory is, up to a normalization factor,

exp(−V ∫_{0}^{t} L[x(s), ẋ(s)]ds),  (36)

where L is the associated Lagrangian defined in (27). In a small time interval [t, t + h] the state moves by Δ := x(t + h) − x(t) according to the Gaussian law (depending also on the actual state x):

\frac{(\frac{V}{2π} )d/2}{(det(D))^{1/2}} exp [−\frac{V}{2} (\frac{Δ}{h} − A(x))^{⊤}D(\frac{Δ}{h} − A(x))] .  (37)

Consider now a dual variable ξ of Δ and compute the Laplace transform L(x, ξ, h) of the Gaussian distribution (37) with respect to Δ. This is obviously a Gaussian integral in Δ, namely,

L(x, ξ, h) =
∫ \frac{(\frac{V}{2π} )d/2}{(det(D))^{1/2}} exp [−\frac{V}{2} (\frac{Δ}{h} − A(x))^{⊤}D(\frac{Δ}{h} − A(x))] exp(Δ.ξ)dΔ. (38)

The critical point which maximizes the exponent is \frac{1}{V} Dξ + A(x)h, and the critical value of the exponent becomes H^{(FP)}(x, ξ)h so that H^{(FP)}(x, ξ) = L(x, ξ, h)/h.

5.2 Path integral for the master equation

For the master equation Eq. (2) the stochastic process is a pure jump process (which sometimes specializes to a birth and death process), for which we cannot apply the method of Subsection 5.1. Nevertheless, we shall see that this method is valid in the large volume limit. First, we shall describe the jump process of the master equation. At time t, the stochastic process of the number of species is N(t). The time evolution is described as follows. During the time interval [t, t + h] the process stays in the state N(t), then a jump N(t) ⟶ N(t) + γ takes place with probability proportional to Wγ(N). The h length of the interval is exponentially distributed with parameter ∑γ Wγ(N). It is proved [15, p. 7753] that

L(ξ, N, h) = exp [V hH^{(M)}(\frac{N}{V} , ξ)] ,  (39)

where H^{(M)} is the Hamiltonian (20) associated to the master equation.

Inverting the Fourier transform of the probability q(x + Δ, h|x) of a transition x ⟶ x + Δ in time h, we have, using (39)

q(x + Δ, h|x) = ∫ exp(iξ.Δ) exp [V hH^{(M)}(x, −\frac{iξ}{V})] \frac{dξ}{(2π)^{s}}  (40)
= (\frac{V}{2π} )^{s} ∫ exp(V iξ.Δ) + hH^{(M)}(x, −iξ)dξ.

For large V, this is again estimated by the saddle point method, so that up to prefactors

q(x + Δ, h|x) ∼ exp [−V L^{(M)}(x, \frac{Δ}{h})h] ,  (41)

where L^{(M)} is the Lagrangian associated to H^{(M)} (and ẋ_{m} = ∂H^{(M)}/∂ξ_{m}.) This is similar to the Feynman path integral [7]. Equation (41) shows that the transition probability is in the form of a standard path integral. Namely, the weight of a trajectory is the exponential exp (−V ∫_{0}^{t} L(x(s), ẋ(s))ds) up to prefactors. We have proved [15, p. 7753] that Eq. (41) for the transition probability is consistent with the large volume asymptotics of the stationary state. In [14] we have seen that the stationary state is up to prefactors: p(x) ∼ exp(−V Φ(x)), where Φ is the solution of the Hamilton– Jacobi equation (19). Then, up to a prefactor one has ∫ q(x′, h|x)p(x) ∼ p(x′). Notice that Eq. (41) is also consistent, because we have proved in [14] that the Lagrangian L^{(M)} is positive.

5.3 Interpretation: Typical paths in the large volume limit

Equations (36) and (41) have the same structure. They assert that in the large volume limit the infinitesimal transition probability of an event x ⟶ x′ in a time interval of length h is exp (−V L[x, (x′ − x/h)h]) , where L is the Lagrangian of the corresponding theory. This means that a typical trajectory in the state space of concentrations will minimize the Lagrangian action and will be the first coordinate of a path [x(s), ξ(s)] satisfying the Hamiltonian equations of the corresponding Hamiltonian. This also implies that in the large volume limit the stochasticity is confined to the choice of initial momenta, ξ_{m}, which will be the cause of the evolution of the x_{m}.

6 A single chemical species:

exact and approximate solutions

6.1 The master equation and its adjoint for exit times

Let us again investigate in detail the case when we have a single chemical species. Let the number of particles be X, and V the volume. Suppose that system evolves only by transitions X ⟶ X ± 1 with probabilities W_{±}(X) per unit time. Let Ω := [aV, bV] be an interval of integers, and T_{Ω}(X) for X ∈ Ω the average of the first exit time of Ω, the stochastic process starting from X. It is well known [9] that T_{Ω}(X) satisfies the adjoint equation of the master equation with the right hand side −1, i.e.

W_{+}(X)[T_{Ω}(X + 1) − T_{Ω}(X)] + W_{−}(X)[T_{Ω}(X − 1) − T_{Ω}(X)] = −1
T_{Ω}(X) = 0 for X ∉ Ω.  (42)

Let us assume that the boundary point aV is absorbing and that bV is reflecting. One can find an exact formula for the solution of Eq. (42) [9]:

T_{Ω}(X) = ∑_{i=aV}^{X} φ(i) ∑_{j=i}^{bV} \frac{1}{W_{−}(j)φ(j)} with φ(i) := ∏_{l=i}^{bV} \frac{W_{+}(l)}{W_{−}(l)}.  (43)

Let us call x := \frac{X}{V} , so that for x ∈ [a, b] W_{±}(X) = V w_{±}(x), and then the continuum limit is:

T (x) = V ∫_{a}^{x} ∫_{y}^{b} \frac{dz}{w_{−}(z)} exp (V ∫_{y}^{z} log(\frac{w_{+}}{w_{−}})) dx′dy.  (44)

6.2 Approximate solutions (simple domains)

We assume that [a, b] contains a single attracting point x_{s} of the vector field A = w_{+} − w_{−}. In Eq. (44) the maximum of the argument in the exponential is obtained for y = a, z = x_{s}, and up to a prefactor

T_{Ω}(x) ∼ exp(V Φ(a|x_{s})),  (45)

where Φ(x|y) := ∫_{y}^{x} log(\frac{w_{−}}{w_{+}} ). We notice immediately that Φ(a|x_{s}) > 0, because w_{+} > w_{−} on [a, x_{s}] (and w_{+} < w_{−} on [x_{s}, b].) Then, T_{Ω}(x) does not depend on x in Eq. (45).

As a comparison, we could use the Fokker–Planck equation associated to the master equation. In that case we should solve

(A(x)\frac{∂}{∂x} + \frac{1}{2V} D(x)\frac{∂^{2}}{∂x^{2}} ) T_{Ω}^{(FP)}(x) = −1 T_{Ω}^{(FP)}(a) = 0 \frac{d}{dx}T_{Ω}^{(FP)}(b) = 0,

(46)

where A = w_{+} − w_{−}, D = w_{+} + w_{−}. As T_{Ω}^{(FP)}(x) = 2V ∫_{a}^{x} ∫_{y}^{b} ^{\frac{exp(2V ∫_{y}^{z} \frac{A}{D})}{D(z)}}dzdy, which can be approximated up to a prefactor by

T_{Ω}^{(FP)}(x) ∼ exp (2V ∫_{a}^{x_{s}} \frac{A}{D} ) .  (47)

If we compare the asymptotic results given by Eqs. (45) and (47), we see that \frac{T_{Ω}}{T_{Ω}^{(FP)}} ∼ exp \{V ∫_{a}^{x_{s}} [log (\frac{1+AD}{1−AD} ) − \frac{2A}{D}]\} , which is exponentially large with re-

spect to the volume.

We can also consider the eigenvalue problems

Λ^{∗}f_{1} = λ_{1}f_{1} R^{∗}f_{1}^{(FP)} = λ_{1}^{(FP)}f_{1}^{(FP)}

with the adjoint operator Λ^{∗} of the master equation, and the adjoint operator R^{∗} of the Fokker–Planck equation. Because these eigenvalues are λ_{1} ∼ \frac{1}{T_{Ω}} , we see that

\frac{λ_{1}}{λ_{1}^{(FP)}} ∼ exp(−V C),  (48)

where C is a certain positive constant.

This means that the Fokker–Planck approximation overestimates the velocity of the rate process, hence also the rate for crossing a barrier with respect to the more exact master equation.

6.3 Qualitative discussion of the typical trajectory

The typical trajectory of the stochastic process associated to the master equation or the Fokker–Planck equation can be described as follows. Essentially, the evolution t ↦ X(t) of the particle number is submitted to the ”force” exerted by the deterministic field A, which is attracting towards x_{s}, and by a diffusion force. The main contribution to the calculation of T_{Ω}(x) comes from trajectories joining point x ∈ Ω to point a in large times. By the results of Section 5 these trajectories are close to Hamiltonian trajectories (either for H^{(M)}, or for H^{(FP)}). For finite times, and for detailed balanced systems by Kurtz theorem [31] such a trajectory is, however, close to the deterministic trajectory, and so goes from x to a small neighborhood of x_{s}. At that point, Kurtz’s theorem is no longer applicable, and we must follow the reasoning outlined up to now. The trajectory finally returns to point a; following essentially the antideterministic trajectory (here this means the reverse of the deterministic trajectory from a to x_{s}). We shall see in Section 7 that this gives back Eqs. (45) or (48).

For example, if A(x) = −αx (α > 0), and D is a positive constant, the equation of motion for H^{(FP)}(x) = −αxp + \frac{1}{2}Dp^{2} are: ẋ = −αx + Dp ṗ = αp.

The trajectory which starts from x(0) > 0 at s = 0, and arrives at x at time t, is

x(s) = x(0)e^{−αs} + (x − x(0)e^{−αt})\frac{sh(αs)}{sh(αt)}.

When x(0)e^{−αt} < x, this trajectory starts at x(0), goes towards 0, reaches a minimum x_{min} ∼ 2\sqrt{xx(0)}e^{−αt/2} (so x_{m} tends to 0 when t → +∞) and goes towards x. In a sense, when t → +∞, the trajectory loses more and more time around 0 (see also [14, p. 7743] for a complete discussion). In the more complex case where Ω = [a, b] contains more than one zero of the vector field A (again with absorbing boundary condition at a and reflecting at b) equations (44) and (46) are still valid and exact. Depending on the position of x with respect to the extrema of Φ T_{Ω}(x) has different expressions which may or may not depend on x ([15, pp. 7755–7756]).

7 Exit times and rate constants

We shall now extend the previous results about the exit time to the multidimensional case. Let us consider the situation of Subsection 5.2 in the large volume limit, so that we define the state space as the space of concentrations x := X/V. We shall treat together the Fokker–Planck dynamics and the master dynamics. Let Ω be a certain domain in the space of concentrations and for a stochastic trajectory, starting from x ∈ Ω let t_{Ω} be the first exit time from Ω. The complementary of the probability distribution of t_{Ω} and the average exit time are

τ_{Ω}(x, t) = P(t_{Ω} > t|X(0) = x)  (49)
T_{Ω}(x, t) = ⟨t_{Ω}|X(0) = x⟩ = − ∫_{0}^{x} t\frac{∂τ_{Ω}(x, t)}{∂t}dt.  (50)

It is well known that [9]

\frac{∂τ_{Ω}(x, t)}{∂t} = Λ^{∗}τ_{Ω}(x, t) τ_{Ω}(x, 0) = 1 τ_{Ω}(x, t) = 0 x ∈ ∂Ω, (51)

where Λ^{∗} is the adjoint operator of the Fokker–Planck operator or of the master equation operator Λ.

Let us consider the eigenvalues λ_{n} of Λ in Ω with absorbing conditions on ∂Ω, ordered by decreasing order (they are negative), and let φ_{n} (resp. ϑ_{n}) be the corresponding eigenvectors of Λ (resp. Λ^{∗}): Λφ_{n} = λ_{n}φ_{n} Λ^{∗}ϑ_{n} = λ_{n}ϑ_{n}. The eigenvectors are normalized so that ∫_{Ω} ϑ_{n}(x)φ_{m}(x)dx = δ_{mn}. Then, δ(x − x_{0}) = ∑_{n} φ_{n}(x)ϑ_{n}(x_{0}).

Let p_{Ω}(t, x|x_{0}) be the probability density of the process x(t) without leaving Ω, i.e. ∫_{Ω_{1}} p_{Ω}(t, x|x_{0})dx is the probability of the event that the process is in the set Ω_{1} without having left the domain Ω in the interval [0, t], under the condition that it has started from x_{0} at time 0, i.e. ∫_{Ω_{1}} p_{Ω}(t, x|x_{0})dx = P(x(t) ∈ Ω_{1}, t ≤ t_{Ω}|x(0) = x_{0}), so that p_{Ω}(t, x|x_{0}) = ∑_{n≤1} e^{λ_{n}t}φ_{n}(x)ϑ_{n}(x_{0}), τ_{Ω}(x_{0}, t) = ∑_{n≤1} e^{λ_{n}t}ϑ_{n}(x_{0}) ∫_{Ω} φ_{n}, and finally,

T_{Ω}(x_{0}) = − ∑_{n} \frac{ϑ_{n}(x_{0})}{λ_{n}} ∫_{Ω} φ_{n}.  (52)

7.1 Estimation of relevant parameters for simple domains

Let us assume in the present subsection that Ω contains a single attracting point x_{s} of the deterministic vector field A and no other zero of A (except possibly on the boundary of Ω).

It can be proved [15, p. 7754–7755] that for large V

T_{Ω}(x) ∼ exp(V min{Φ(y|x_{s}); y ∈ ∂Ω}),

(up to a prefactor), and Φ(x|x_{s}) is the nonconstant regular solution of the Hamilton– Jacobi equation vanishing at x_{s}:

H^{(M)}(x, ∇Φ(x|x_{s})) = 0 Φ(x_{s}|x_{s}) = 0.  (53)

We have also shown above that the function Φ is unique. In particular, T_{Ω}(x) is independent from x (up to a prefactor). Clearly, the prefactor is such that T_{Ω}(x) should tend to zero for x approaching the boundary of Ω, but we are unable to give the form of this prefactor. Then, from Eq. (52) we see that ϑ_{1}(x) ∼ K, where K is an absolute constant (up to a prefactor). However, because of one of our normalization conditions ∫_{Ω} φ_{1}ϑ_{1} = 1, we see from Eq. (52) that

−λ_{1} ∼ \frac{1}{T_{Ω}(x)} ∼ exp(−V min\{Φ(y|x_{s}); y ∈ ∂Ω\})  (54)

and that τ_{Ω}(x, t) = exp(λ_{1}t).

The estimation of the first eigenvalue in Eq. (54) has been given by several authors for the Fokker–Planck dynamics [36, 50]. We remark that the proof in [36] is not logically consistent, although the result is correct (see [15, pp. 7754–7755]). If we allow the domain Ω to contain several critical points of Φ beyond a stable equilibrium point x_{s} then the mean exit time depends on the starting point. The full discussion of the mean exit times is given in [15, pp. 7756–7757].

8 Progress variables

8.1 Fundamental processes and progress variables

Now we suppose that we have two kinds of chemical species in the vessel of fixed volume V:

1. Internal species denoted by X_{m} (m = 1, 2, …, M) which are varying freely according to the chemical reactions in the vessel. The number of particles of species X_{m} is denoted by X_{m} and its concentration is x_{m} := X_{m}/V as above.

2. External species denoted by A_{l} (n = 1, 2, …, N), which are completely under control of external reservoirs.

The state of the system at time t is given by the vector x(t) which specifies the concentrations of the freely varying species. The concentrations a(t) may be functions of the state variables. The simplest situation is the case where the reservoirs maintain each a_{l} at a fixed concentration independent on t. We allow a general dependence in Sections 8–9, while in [18] the action of reservoirs has been discussed in more detail.

The species are reacting and diffusing in R subprocesses, of which we shall separately count the forward and backward steps from now on.

∑_{n=1}^{N} μ_{r}^{n}A_{n} + ∑_{m=1}^{M} α_{r}^{m}X_{m} ⟺ ∑_{n=1}^{N} ν_{r}^{n}A_{n} + ∑_{m=1}^{M} β_{r}^{m}X_{m},  (55)

and r labels the processes, as before. By convention, the forward process is from left to right, and the backward process is from right to left.

If from now on ξ_{r}(t) denotes the number of forward processes minus the number of backward processes of type r up to time t, then we have [31]

X(t) = X(0) + ∑_{α} γ_{r}^{m}ξ_{r}(t),  (56)

where γ_{r}^{m} is the difference between the stoichiometric numbers: γ_{r}^{m} := β_{r}^{m} − α_{r}^{m}. Let u_{r} be the number of processes per unit volume:

u_{r} := \frac{ξ_{r}}{V}  (57)

and call these numbers progress variables as usual. Finally, W_{r}^{±} denotes the probability per unit time that a forward or backward process of type α occurs in V. Because of Eq. (56) in rescaled variables we have

x(t) = x(0) + τ u(t),  (58)

and W_{α}^{±} can be considered as function of u, x(0), and a(t).

8.2 Dynamics of progress variables

The probability Q(γ, t) that at time t r_{α} processes of type α have occurred satisfies the following master equation

\frac{∂Q(γ, t)}{∂t} = ∑_{β} [W_{β}^{+}(γ − e_{β})Q(γ − e_{β}, t) + W_{β}^{−}(γ + e_{β})Q(γ + e_{β}, t)]
− (W_{β}^{+}(γ) + W_{β}^{−}(γ))Q(γ, t).  (59)

Let us now turn to the rescaled variables defined in (57) to define the rescaled probability densities q(u, t) and rates w_{α}^{±} by q(u, t) := V^{p}Q(γ, t) w_{α}^{±}(u) := \frac{W_{α}^{±}(γ)}{V}, and with these (59) may be rewritten as

\frac{1}{V} \frac{∂q(u, t)}{∂t} =

∑_{β} [w_{β}^{+}(u − e_{β})q(u − e_{β}, t) + w_{β}^{−}(u + e_{β})q(u + e_{β}, t)]
−(w_{β}^{+}(u) + w_{β}^{−}(u))q(u, t).  (60)

For large V we obtain an approximate Fokker–Planck equation from Eq. (60):

\frac{∂q(u,t)}{∂t} = − ∑_{β} \frac{∂(w_{β}^{+}(u)−w_{β}^{−}(u))q(u,t)}{∂u_{β}} + \frac{1}{2V} ∑_{β} \frac{∂^{2}(w_{β}^{+}(u)+w_{β}^{−}(u))q(u,t)}{∂u_{β}^{2}}. (61)

However, it has been shown that this equation does not give the correct results for stationary state and large time dynamics of the full master equation (see [12] and above, or [14]).

8.3 Approximate dynamics for large volumes

For large V again one may try a formal asymptotic expansion for q

q(u, t) = exp(−V ψ)(U_{0}(u, t) + \frac{1}{V} U_{1}(u, t) + …).  (62)

The variations of the prefactor U_{0} are usually negligible compared to the variations of the dominant exponential exp(−V ψ.) However, the prefactor becomes preponderant near criticality, when the argument of the exponential vanishes. On the other hand, the method of Kubo would make no sense if U_{0} was singular. For all these reasons, a further study of the prefactor is necessary for a sound mathematical foundation of the present formalism.

The expansion (62) is used either in the Fokker–Planck equation (61), or directly in the master equation (60). In both cases it is easily shown as above that the function ψ satisfies a Hamilton–Jacobi equation in the progress variables:

H′(u, ∇ψ(u)) = 0,  (63)

where H′ is either the master Hamiltonian H^{′(M)}, or the Fokker–Planck Hamiltonian H^{′(FP)},

H^{′(M)}(u, π) = ∑_{r}(w_{r}^{+} − w_{r}^{−})π_{r} + \frac{1}{2}(w_{r}^{+} + w_{r}^{−})π_{r}^{2},  (64)
H^{′(FP)}(u, π) = ∑_{r}[w_{r}^{+}(exp(π_{r}) − 1) + w_{r}^{−}(exp(−π_{r}) − 1)],  (65)

where π_{r} is the conjugate momentum of u_{r}. As above, we have H^{′(M)}(u, π) = H^{′(FP)}(u, π) + O(|π|^{3}).

We notice here that all the results proved previously for the Hamilton–Jacobi equation associated with the master equation are valid for the Hamilton–Jacobi equation in the progress variables.

8.4 Relation with the usual master equation

The usual master equation is an equation for functions of the state variable x. Although the formal derivation of the usual master equation and of its Fokker–Planck and Hamilton–Jacobi approximations are the same as in Subsection 8.3, their physical meaning is completely different, because u_{r} are not state variables. They are related to state variable by Eq. (58). Upon knowing u one can deduce the variation x(t) − x(0) of the state variable, but in general conversely not. Previously we have introduced the rate per unit volume wγ(x) of a transition x ⟶ x + γ/V. It is clear that

wγ(x) = ∑_{r:τ_{r}^{·}=γ} w_{r}^{+}(x) + ∑_{r:τ_{r}^{·}=γ} w_{r}^{−}(x).  (66)

The usual Kramers-Moyal expansion of the master expansion studied above and in textbooks [24, 26, 27, 39] and papers [12, 17] yields the usual Fokker–Planck equation in concentration variables,

\frac{∂p(x, t)}{∂t} = − ∑_{m=1}^{M} \frac{∂}{∂x_{m}} (A_{m}p)(x, t) + \frac{1}{2V} ∑_{m,q=1}^{M} \frac{∂^{2}}{∂x_{m}∂x_{q}} (D_{mq}p)(x, t),

where

A_{m}(x) := ∑ ; γ γ_{m}wγ(x) = ∑_{r} τ_{r}^{m}(w_{r}^{+}(x) − w_{r}^{−}(x)),  (67)
D_{mq}(x) := ∑ ; γ γ_{m}γ_{q}wγ(x) = ∑_{r} τ_{r}^{m}τ_{r}^{q}(w_{r}^{+}(x) − w_{r}^{−}(x)).  (68)

Finally, Eq. (58) shows that the evolution remains for all time in the subset

E(x_{0}) := {x ∈ ℝ^{M}; x = x_{0} + ∑_{r} τ_{r}^{m}u_{r}; u_{r} ∈ ℝ}.

It is clear that this set is the same as the set S(x_{0}) introduced in (11), because the set of 2R vectors γ_{r}^{·} (r = 1, 2, …, R) is exactly the same as the set of vectors γ(·, r) (r = 1, 2, …, R), and each subset E(x_{0}) carries a stationary probability distribution. We shall assume henceforth that we reduce the situation to a given set E(x_{0}), so that the dynamics is irreducible and has a unique stationary state in this set.

Notice that if d is the dimension of E(x_{0}) we can use d variables to parametrize E(x_{0}). The other variables x_{d+1}, x_{d+2}, …, x_{M} are still present (so that the chemical processes are the same) but they are linear functions of x_{1}, x_{2}, …, x_{d}.

If one introduces π := τ · ξ x := x(0) + τ · u, then (66) implies that the master Hamiltonian H^{′(M)} in the variables (u, π) reduces to the master Hamiltonian H^{(M)} in the variables (x, ξ) : H^{′(M)}(u, π) = H^{(M)}(x, ξ). In particular, if Φ is a solution of the Hamilton–Jacobi equation

H^{(M)}(x, ∇Φ(x)) = 0,  (69)

then the function ψ(u) := Φ(x(0) + τ · u) induces a solution of the Hamilton–Jacobi equation in the progress variable form:

H^{′(M)}(u, ∇ψ(u)) = 0,  (70)

because ∇ψ(u) = ∇Φ(x) · τ .

But conversely, a solution ψ of the Hamilton–Jacobi equation (70) does not necessarily produce a function in the state variables x and a fortiori does not define a solution of Eq. (69).

8.5 Free energy and rate constants in the unconstrained system

We consider now the vessel of volume V, in which the 2R processes take place, but now we switch off the exchange of external species A_{n} with the reservoirs (still maintaining the temperature T constant), so that the concentrations x_{m} and a_{n} vary freely according to the natural chemical processes (r = 1, 2, …, R) in the vessel. The state then will reach a thermal equilibrium. At thermal equilibrium, the probability distribution on the state space which consists now of the freely varying concentrations x_{m} and a_{n} is for large V : p_{eq}(x, a) ∼ U_{0} exp(−\frac{V F (x,a)}{k_{B}T}), where F is the free energy (of the state (x, a)) per unit volume, k_{B} is the Boltzmann constant, T is temperature, and U_{0} is a prefactor. At equilibrium, all processes r are supposed to satisfy detailed balance, as it is usual when dealing with realistic physical systems, written asymptotically for large V as w_{r}^{+}(x, a) exp(−\frac{V}{k_{B}T} F (x, a)) = w_{r}^{−}(x, a) exp(−\frac{V}{k_{B}T} F (x + \frac{τ _{r}^{·}}{V} , a + \frac{t_{r}^{·}}{V})), where t := ν − μ. Therefore,

k_{B}T log(\frac{w_{r}^{−}}{w_{r}^{+}} ) = \frac{∂F}{∂x} · τ + \frac{∂F}{∂a} · t.  (71)

For perfect gases or solutions one usually assumes that the rates are of the mass action form: w_{r}^{+}(x) = k_{r}^{+}x^{α{· \atop r}} w_{r}^{−}(x) = k_{r}^{−}x^{β{· \atop r}}, where k_{r}^{±} are temperature dependent rate coefficients. One can immediately check that the usual partial equilibrium form [13] of the free energy

F (x, a) = ∑_{m} F_{m}(x_{m}) + ∑_{n} F_{n}(a_{n})  (72)

(where F_{m} is the free energy of the ideal gas law at temperature T and concentrations x) does satisfy Eq. (71). In fact, the chemical potentials are

μ_{m} = \frac{∂F}{∂x_{m}} = \frac{dF_{m}}{dx_{m}} = k_{B}T log(x_{m}) + f_{m}(T),  (73)

and Eq. (71) reduces to

k_{B}T log(\frac{k_{r}^{+}}{k_{r}^{−}}) = ∑_{m} τ_{r}^{m}f_{m}(T) + ∑_{n} t_{r}^{n}f_{n}(T).  (74)

Here each f_{m} is calculated using the partition functions of the internal degrees of freedom of the species X_{m} and Eq. (74) is the usual expression for the equilibrium constant K_{r,eq} of the process r in term of the partition function of the internal degrees of freedom of the species appearing in the given process.

In many circumstances, like for imperfect gases or solutions, electrolytes, etc. one needs a more general formulation of the free energy F, not necessarily of the form (72). This is why we shall work with the most general form of the free energy.

From Eq. (71) it follows for all r

w_{r}^{+}(x, a) exp (\frac{1}{k_{B}T} \frac{∂F (x,a))}{∂u_{r}} − 1) + w_{r}^{−}(x, a) exp (−\frac{1}{k_{B}T} \frac{∂F (x,a))}{∂u_{r}} − 1) = 0, (75)

where \frac{∂F}{∂u_{r}} = ∑_{m} \frac{∂F}{∂x_{m}} τ_{r}^{m} + \frac{∂F}{∂u_{n}} t_{r}^{n}, so that \frac{1}{k_{B}T} F satisfies the Hamilton–Jacobi equation (63): H^{′(M)}(u, \frac{1}{k_{B}T} \frac{∂F}{∂u} ) = 0 with H^{′(M)} given by (65), where in w_{r}^{±} the variables are x = x(0) + τ · u a = a(0) + t · u.

9 Dissipation

9.1 Dissipation of information

From now on, we shall assume again that on a given set E(x_{0}) the state of the vessel reaches a stationary state p(x) ∼ U_{0} exp(−V Φ(x)) with the concentrations a being entirely controlled by the reservoirs and having fixed variations a(t).

The state x(t) macroscopically evolves according to the deterministic equations \frac{dx}{dt} = A(x) = τ · (w^{+}(x(t)) − w^{−}(x(t))). The value of the state function Φ evolves as

\frac{dΦ}{dt} = ∇Φ(x(t))ẋ(t) = ∑_{m,r}(w_{r}^{+}(x(t)) − w_{r}^{+}(x(t)))\frac{∂Φ}{∂u_{r}} .  (76)

We know from Subsection (8.5) that H^{′(M)}(u, ∇Φ(u)) = 0. But e^{x} − 1 ≥ x, so that 0 = ∑_{r} [w_{r}^{+}(exp(\frac{∂Φ}{∂u_{r}}) − 1) + w_{r}^{−}(exp(−\frac{∂Φ}{∂u_{r}}) − 1)] ≥ ∑_{r}(w_{r}^{+} − w_{r}^{+})\frac{∂Φ}{∂u_{r}} . As

a consequence, we obtain the inequality \frac{dΦ}{dt} ≤ 0, where equality is attained if and only if w_{r}^{+}∂Φ/∂u_{r} = w_{r}^{−}∂Φ/∂u_{r} for each r. In particular, if the deterministic state reaches a stationary point x^{(0)}, for which A(x^{(0)}) = 0, then dΦ/dt = 0, and either ∂Φ/∂u_{r} = 0, or w_{r}^{+} = w_{r}^{−} = 0.

The quantity dΦ/dt computed along a deterministic trajectory is always negative. It can be interpreted as a dissipation of information per unit time.

In fact, V Φ(x) can be considered as the average information which is obtained when the system is observed in the state x rather then being stochastically distributed with the stationary probability distribution p_{s}(x) ∼ exp(−V Φ(x)), which is the state of lowest information, when the system is coupled to the various reservoirs of heat and of chemical species A_{n}. Along a deterministic path Φ decreases with time while the state x tends to a deterministic stationary state which is a local minimum of Φ.

9.2 Dissipation of energy

The variation of free energy along the deterministic trajectory is

\frac{dF}{dt} = ∑_{i} \frac{∂F}{∂x_{i}} \frac{dx_{i}}{dt} + ∑_{i} \frac{∂F}{∂a_{l}} \frac{da_{l}}{dt} ,  (77)

or,

\frac{dF}{dt} = ∑_{i} \frac{∂F}{∂x_{i}} \frac{dx_{i}}{dt} + ∑_{i} \frac{∂F}{∂a_{l}} [\frac{da_{l}}{dt} ]_{c} − ∑_{i} \frac{∂F}{∂a_{l}} ([\frac{da_{l}}{dt} ]_{c} − \frac{da_{l}}{dt} ) ,  (78)

where [\frac{da_{l}}{dt} ]_{c} is the variation of a_{l} due to the chemical processes. It turns out [18][p. 684] that \frac{dF}{dt} −w ≤ 0, where w := − ∑_{l} m_{l} ([\frac{da_{l}}{dt} ]_{c} − \frac{da_{l}}{dt} ) , and the quantity \frac{dF}{dt} −w is the dissipation of energy in the system per unit time.

9.3 Inequality between the dissipation of information and that of energy

It is proved in Appendix A of [18] that the dissipation of information and the dissipation of energy satisfy the fundamental inequality \frac{1}{k_{B}T} (\frac{dF}{dt} − w) ≤ \frac{dΦ}{dt} ≤ 0, so that in absolute value the dissipation of information is always less than the dissipation of energy. Moreover, there is equality if and only if we have an equilibrium situation.

10 Acknowledgements

The present work has partly been supported by the National Science Foundation Hungary (Nos T037491 and T047132).

References

[1] Arnold, L.: On the consistency of the mathematical models of chemical reactions, In: Dynamics of synergetic systems, Ed. H. Haken, Springer, Berlin, 1980, pp. 107-118.

[2] Courant, R. and Hilbert, D.: Methods of mathematical physics, Interscience, New York, 1953.

[3] Cover, T. M. and Thomas, J. A.: Elements of information theory, Wiley, New York, 1991.

[4] Englman, R. and Yahalom, A.: Variational procedure for time dependent processes, Physical Review E, 69 (2004), 026120(10 pages).

[5] Érdi, P., Tóth, J.: Mathematical models of chemical reactions. Theory and applications of deterministic and stochastic models, Princeton University Press, Princeton, 1989.

[6] Érdi, P., Tóth, J.: Stochastic reaction kinetics—Nonequilibrium thermodynamics of the state space? Reaction Kinetics and Catalysis Letters 4 (1976), 81–85.

[7] Feynman, R. P. and Hibbs, A.: Quantum mechanics and path integrals, McGraw-Hill, New York, 1965.

[8] Forster, D.: Hydrodynamic fluctuations, broken symmetry, and correlation functions, Frontiers in physics, Vol. 47, Benjamin Cummings, Reading, MA, 1975.

[9] Gardiner, H.: Handbook of stochastic methods, Springer, Berlin, 1983.

[10] Gaveau, B., Martinás, K., Moreau, M., Tóth, J.: Entropy, extropy and information potential in stochastic systems far from equilibrium, Physica A 305 (2002), 445– 466.

[11] Gaveau, B.; Moreau, M. and Tóth, J.: Decay of the metastable state: different predictions between discrete and continuous models, Letters in Mathematical Physics 37 (1996), 285–292.

[12] Gaveau, B.; Moreau, M. and Tóth, J.: Master equation and Fokker–Planck equation: Comparison of entropy and rate constants, Letters in Mathematical Physics 40 (1997), 101–115.

[13] Gaveau, B.; Moreau, M. and Tóth, J.: Dissipation of energy and information in nonequilibrium reaction-diffusion systems, Physical Review E 58 (11) (1998), 5351–5354.

[14] Gaveau, B.; Moreau, M. and Tóth, J.: Variational nonequilibrium thermodynamics of reaction-diffusion systems. I. The information potential, Journal of Chemical Physics 111 (17) (1999), 7736–7747.

[15] Gaveau, B.; Moreau, M. and Tóth, J.: Variational nonequilibrium thermodynamics of reaction-diffusion systems. II. Path integrals, large fluctuations and rate constants, Journal of Chemical Physics 111 (17) (1999), 7748–7757.

[16] Gaveau, B.; Moreau, M. and Tóth, J.: Path-integrals from peV to TeV, (Florence, August 1998), Casalbroni R., ed., World Scientific, Singapore, 1999, pp. 52.

[17] Gaveau, B.; Moreau, M. and Tóth, J.: Information potential and transition to criticality for certain two-species chemical systems, Physica A 277 (2000), 455– 468.

[18] Gaveau, B.; Moreau, M. and Tóth, J.: Variational nonequilibrium thermodynamics of reaction-diffusion systems. III. Progress variables and dissipation of energy and information, Journal of Chemical Physics 115 (2) (2001), 680–690.

[19] Gaveau, B. and Schulman, L. S.: Master equation based formulation of nonequilibrium statistical mechanics, Journal of Mathematical Physics 37 (8) (1996), 3897–3932.

[20] Gaveau, B. and Schulman, L. S.: A general framework for non-equilibrium phenomena: The master equation and its formal consequences, Phys. Lett. A 229 (6) (1997), 347–353.

[21] Gaveau, B. and Schulman, L. S.: Theory of nonequilibrium first-order phase transitions for stochastic dynamics, Journal of Mathematical Physics 39 (3) (1998), 1517–1533.

[22] Gyarmati, I.: On the governing principle of dissipative processes and its extension to nonlinear problems, Annalen der Physik, 23 (1969), 353–378.

[23] Gyarmati, I.: Non-equilibrium thermodynamics. Field theory and variational principles, Springer Verlag, 1970.

[24] De Groot, S. R. and Mazur, P.: Nonequilibrium thermodynamics, Dover, New York, 1984.

[25] Horn, F. J. M. and Jackson, R.: General mass action kinetics, Archive for Rational Mechanics and Analysis 47 (1972), 81–116.

[26] Keizer, J.: Statistical thermodynamics of nonequilibrium processes, Springer, New York, 1987.

[27] Kubo, R.: Statistical mechanics, North Holland, Amsterdam, 1988.

[28] Kubo, R.; Matsuo, K. and Kitahara, K.: Fluctuation and relaxation of macrovariables, Journal of Statistical Physics, 9 (1973), 51–96.

[29] Kubo, R.; Toda, M. and Hashitsume, N. In: Statistical physics II, Springer, New York, 1986.

[30] Kullback, S.: Information theory and statistics, Wiley, New York, 1959.

[31] Kurtz, T. G.: The relationship between stochastic and deterministic models in chemical reactions, Journal of Chemical Physics 57 (1973), 2976–2978.

[32] Landau, L. D. and Lifsic, E. M.: Mechanics: Volume 1, Butterworth-Heinemann; 3rd edition, 1976.

[33] Landau, L. and Lifsic, E. M.: Statistical physics, Pergamon, Oxford, 1980.

[34] Lemarchand, H.: Asymptotic solution of the master equation near a nonequilibrium transition, Physica A 101 (1980), 518–534.

[35] Lemarchand, H. and Nicolis, G.: Stochastic analysis of symmetry-breaking bifurcations: master equation approach, Journal of Statistical Physics, 37 (1984), 609–629.

[36] Ludwig, D.: Persistence of dynamical systems under random perturbations, SIAM Rev. 17 (1975), 605–640.

[37] Moreau, M.: Letters in Mathematical Physics 1 (1975), 7.

[38] Moreau, M.: Note on the entropy production in a discrete Markov system, Journal of Mathematical Physics 19 (12) (1978), 2494–2498.

[39] Nicolis, G. and Prigogine, I.: Self-organization in nonequilibrium systems, Wiley, New York, 1977.

[40] Prigogine, I.: Nonequilibrium statistical mechanics, Wiley, New York, 1962.

[41] Prigogine, I. and Nicolis, G.: In: From theoretical physics to biology, Ed. M. Marois, Karger, Basel, 1973.

[42] Ruelle, D.: Statistical mechanics, Benjamin, New York, 1970.

[43] Schulman, L. S.: In: Finite size scaling and numerical simulation of statistical systems, Privman, V. ed., World Scientific, Singapore, 1990.

[44] Sieniutycz, S.: Some thermodynamic aspects of development and bistability in complex multistage systems, Open Sys. & Information Dyn. 7 (2000), 309–326.

[45] Sieniutycz, S.: Optimality of nonequilibrium systems and problems of statistical thermodynamics, International Journal of Heat and Mass Transfer, 45 (2002), 1545–1561.

[46] Tóth, J. and Érdi, P.: Indispensability of stochastic models, In: Nonlinear dynamics and exotic kinetic phenomena in chemical systems, (Gy. Bazsa ed.) Debrecen, Budapest, Gödöllő, 1992, pp. 117–143 (in Hungarian).

[47] Ván, P.: On the structure of the ’Governing Principle of Dissipative Processes’, Journal of Non-equilibrium Thermodynamics 21 (1) (1996) 17–29.

[48] Ván, P. and Muschik, W.: The structure of variational principles in nonequilibrium thermodynamics, Physical Review E 52 (4) (1995), 3584–3590.

[49] Ván, P. and Nyíri, B.: Hamilton formalism and variational principle construction, Annalen der Physik (Leipzig) 8 (1999), 331–354.

[50] Ventsel, A. and Freidlin, M.: On small random perturbation of dynamical systems, Russian Math. Surveys 25 (1970), 1–55.

[51] Vlad, M. O. and Ross, J.: Fluctuation-dissipation relations for chemical systems far from equilibrium, Journal of Chemical Physics 100 (10) (1994), 7268–7278.

[52] Vlad, M. O. and Ross, J.: Random paths and fluctuation-dissipation dynamics for one-variable chemical systems far from equilibrium, Journal of Chemical Physics 100 (10) (1994), 7279–7294.

[53] Vlad, M. O. and Ross, J.: Thermodynamic approach to nonequilibrium chemical fluctuations, Journal of Chemical Physics 100 (10) (1994), 7295–7309.

[54] Zheng, Q.; Ross, J.; Hunt, K. L. C. and Hunt, P. M.: Stationary solutions of the master equation for single and mulitintermediate autocatalytic chemical systems, Journal of Chemical Physics 96 (1) (1992), 630–641.