Deterministic and Stochastic Asymptotic Behavior of Nonequilibrium

8756

J. Phys. Chem. B 1997, 101, 8756-8773

Deterministic and Stochastic Asymptotic Behavior of Nonequilibrium Chemical Systems with Time-Dependent Rate Coefficients† Marcel Ovidiu Vlad‡,§ and John Ross*,‡ Department of Chemistry, Stanford UniVersity, Stanford, California 94305-5080, and Center of Mathematical Statistics, Casa Academiei Romane, Calea 13 Septembrie 13, 76100 Bucharest, Romania ReceiVed: January 23, 1997; In Final Form: June 11, 1997X

We analyze the long time limit of the nonequilibrium solutions of a system of multivariable nonlinear kinetic equations with time-dependent rate coefficients, as achieved, for example, by temporal variation of the temperature. If the characteristic time scale attached to the change of rate coefficients is smaller than the relaxation time to equilibrium, then the system is constrained to evolve away, possibly far, from equilibrium. However, after a sufficiently large time the system forgets its past: in the long run all solutions of the kinetic equations tend toward a special (normal) solution which depends on the previous values of the rate coefficients, but it is independent of the initial state of the system. The normal solution may be very different from the equilibrium solution. The occurrence of this type of time-dependent normal regime for the deterministic kinetic equations is intimately connected to a similar behavior of the stochastic evolution equation of the system. In the long run all solutions of the stochastic equation for the state probability also tend toward a normal form which is independent of the initial preparation of the system. The logarithm of the normal form of the state probability is a Lyapunov function of the system of deterministic kinetic equations and plays the role of a generalized thermodynamic potential which may be used for developing a thermodynamic description of the chemical process. A Gibbsian ensemble description is introduced in terms of a multireplica stochastic master equation. The logarithm of the large time solution of the multireplica stochastic master equation is also a Lyapunov function, one for the stochastic evolution equations of the system. If the time dependence of the rate coefficients is generated by the variation of an external constraint, then the deterministic normal regime can be considered as representing the response of the system to an excitation. The coupling between the excitation and the system is nonlinear and multiplicative and generates interesting effects. A kinetic experiment is suggested for checking the existence of normal solutions for chemical systems with timedependent rate coefficients.

1. Introduction The classical treatments of chemical kinetics of ideal systems are based on the assumption that the rate coefficients are dependent on the state of the system defined by a suitable set of macroscopic variables (temperature, total pressure, external electric or magnetic fields, etc.), but are concentration- and timeindependent. There are however many processes for which the rate coefficients are time-dependent. One type of example arises when a chemical process is operated under externally controlled conditions, for instance under the influence of a periodic pressure field generated by a sound wave1,2 or under the influence of a given program of temperature variation.3,4 For example, in geochemistry, the quasiperiodical daily variations of temperature may have a strong effect on the time evolution of chemical reactions occurring in the atmosphere or in the ocean.5,6 A second type of situation arises when the time dependence of the rate coefficients is due to the reaction mechanism itself. For instance, for bimolecular diffusion-limited recombination reactions the rate coefficient has a time-dependent component which tends to zero for large times;7 in this case the time dependence of the rate coefficient leads only to small transient corrections for small times. An essentially different situation corresponds to a chemical process that takes place in a disordered system. In this case the energetic, structural, and †

This study is dedicated in honor of Daniel Kivelson, an outstanding scientist with the highest standards of scholarship. ‡ Stanford University. § Casa Academiei Romane. X Abstract published in AdVance ACS Abstracts, October 1, 1997.

S1089-5647(97)00306-4 CCC: $14.00

temporal disorder of the matrix in which the reaction process takes place may generate a dispersive kinetic regime for which the time dependence of the rate coefficients persists in the long run, and moreover is rate-determining for large times.8-11 Yet another situation corresponds to fast chemical processes for which there is a superposition between the time scale of the thermalization process describing the dynamics of nonreactive collisions and the time scale of the chemical processes.12-14 For rate processes with time-independent rate coefficients it has been shown that for large times a chemical process forgets its past; that is, the solutions of the kinetic equations tend for large time toward certain normal forms, which are independent of the initial preparation of the system.15,16 These normal forms may be one or many stable steady states, a stable limit cycle, or even a chaotic trajectory. The large time deterministic behavior of chemical systems is intimately related to the dynamics of concentration fluctuations.17-21 If the system is not subject to external constraints, the system eventually reaches the state of thermodynamic equilibrium for which the probability Peq of fluctuations is determined by Einstein’s fluctuation formula:

Peq ∝ exp(∆S/kB)

(1.1)

where ∆S is the variation of entropy due to fluctuations and kB is Boltzmann’s constant. Both the entropy of the system and the probability of fluctuations have extremum values for the state of thermodynamic equilibrium. Due to this extremum property, the entropy is, up to a constant additive factor, a © 1997 American Chemical Society

Nonequilibrium Chemical Systems

J. Phys. Chem. B, Vol. 101, No. 43, 1997 8757

Lyapunov function of the system of deterministic kinetic equations. A similar connection between the large time deterministic behavior of chemical processes and the fluctuation dynamics exists if, due to certain suitable constraints, the system cannot reach the state of thermodynamic equilibrium. By using different stochastic approaches, a number of researchers have shown that the large time behavior of chemical fluctuations is connected to the large time behavior of the solutions of the deterministic equations. A number of nonequilibrium analogs of the Einstein’s fluctuation formula (1.1) have been suggested, where the entropy S is replaced by various generalized nonequilibrium thermodynamic functions determined by the type of stochastic model considered. One theory has been suggested by Keizer18 (K) based on a set of phenomenological assumptions concerning the fluctuation dynamics and a set of fluctuationdissipation relations. A second approach has been suggested by Graham19 (G) based on the assumption that the fluctuations can be described by means of a Fokker-Planck equation. A third theory has been developed by Ross, Hunt, and Hunt20-24 (RHH) by assuming that the chemical fluctuations can be described by a chemical master equation. As far as we know, no similar thermodynamic and stochastic approaches have been used for the description of the large-time behavior of chemical systems with time-dependent rate coefficients. In this article we present a contribution to this subject, the investigation of the relationships between the large time deterministic and stochastic behavior of chemical systems with time-dependent rate coefficients. The comparative analysis of the stochastic and deterministic descriptions is used for developing a thermodynamic description of time-dependent chemical systems which is a direct generalization of the RHH approach for time-independent systems.20-24 The approach developed in this paper is the following. We start out from a homogeneous chemical master equation for a nonlinear chemical system and investigate the large time behavior of its solutions by generalizing a Lyapunov function suggested by Brey and Prados25 for the study of linear relaxation processes in glasses described by a master equation with discrete state variables and a finite number of states. We suggest a physical interpretation of this Lyapunov function by using an ensemble description of chemical fluctuations based on the Ramakrishnan’s approach to the theory of random point processes.26 The connection between the stochastic and deterministic descriptions is made by using an eikonal approximation27-30 for the solution of the master equation for large systems. Within the eikonal approximation the probability of fluctuations is determined by an equation of the type (1.1) where the entropy S is replaced by a chemical action which is the solution of an equation of the Hamilton-Jacobi type. This chemical action is a Lyapunov function of the deterministic kinetic equations. The two Lyapunov functions, for the master equation and for the kinetic equations, are used for deriving a thermodynamic description of time-dependent chemical systems. This thermodynamic description is possible due to the existence of the normal solutions of the deterministic and kinetic equations, that is, due to the fact that after enough large time the system forgets its past. 2. Stochastic Asymptotic Behavior We consider a closed chemical system described by a Markovian master equation

dPN(t) )

[WN′N(t) PN′(t) - WNN′(t) PN(t)] ∑ N′

(2.1)

where N is the composition vector of the system

N ) (N1, N2, ...) N1, N2, ... are the numbers of molecules of different chemicals present in the system, PN(t) is the state probability at time t and WNN′(t) is the transition rate from N to N′ at time t. We assume that there is no exchange of chemicals between the system and the environment; there is however an exchange of energy (for instance thermal or electromagnetic energy) which leads to an explicit time dependence of the transition rates. We assume that all states of the system are connected;31 that is, for any vectors N, N′, there is at least a path N f N1 f N2 f ... f N′ connecting the two states N and N′ for which the corresponding rates are different from zero. At this stochastic level of description it is not necessary to introduce detailed assumptions concerning the stoichiometry of the process. The stoichiometry, however, is very important for the description of the deterministic asymptotic behavior. The stoichiometric constraints necessary for the occurrence of the deterministic kinetic equations are discussed in section VIII. For chemical systems the number of vectors N, N′, ... is a possibly large but finite number. It follows that by attaching to each vector N a certain numerical label we can transcribe the master equation (2.1) in a matrix notation:

dP(t)/dt ) P(t) A(t)

(2.2)

P(t) ) (PN1, PN2, ...)

(2.3)

where

WNN′(t)|| ∑ N′

A(t) ) ||ANN′(t)|| ) ||(1 - δNN′)WNN′(t) - δNN′

(2.4)

are matrices with large but finite numbers of elements. According to these definitions, the probability vector is defined as a row vector rather than a column vector. Although rather unusual, this convention is mathematically consistent and leads to simpler equations. The formal solution of (2.2) may be represented in terms of a Green function,

∫t tA(t′) dt′)

ˆ exp( G(t|t0) ) T

(2.5)

0

which is the solution of the equation

dG(t|t0)/dt ) G(t|t0) A(t), G(t0|t0) ) I

(2.6)

P(t) ) P(t0) G(t|t0)

(2.7)

We have

Here T ˆ is the time ordering operator and I is the unity matrix. To investigate the asymptotic behavior of P(t), we consider pairs of different solutions Pu(t) and Pu′(t) of (2.2) where u, u′ ) 1, 2, .... These solutions correspond to different initial conditions Pu(t0) and Pu′(t0), respectively. We introduce the probability functional:

hu,u′ )

(u) (u′) P(u) ∑ N (t) ln[PN (t)/PN (t)] N

(2.8)

(u′) where P(u) N (t) and PN (t) are elements of Pu(t) and Pu′(t),

8758 J. Phys. Chem. B, Vol. 101, No. 43, 1997

Vlad and Ross

respectively. huu′ is an H function; that is, it fulfills the conditions

huu′ > 0, dhuu′/dt < 0, for Pu(t) * Pu′(t)

(2.9a,b)

huu′ ) 0, dhuu′/dt < 0, for Pu(t) ) Pu′(t)

(2.10a,b)

For the proof of (2.9)-(2.10) see Appendix A. For nonequilibrium steady states the stationary regime is usually investigated by considering the limit t f ∞. In our case, however, this limit cannot be used, because the normal solution which emerges in the limit t f ∞ is generally timedependent. An alternative approach is to maintain the current value of the time variable finite and to push the initial condition to minus infinity, t0 f -∞. The replacement of the limit t f ∞ with the limit t0 f -∞ should not be confused with the condition of time reversal, which is the cause of detailed balance in equilibrium statistical mechanics. The condition t0 f -∞ simply means that the process has been going on for a very long time, long enough to reach the normal solution. We emphasize that the master equation (2.1) does not obey the condition of time reversal. An important consequence of (2.9a,b) and (2.10a,b) is that for t0 f -∞ we have Pu(t) f Pu′(t). It follows that

Pu(t) f Pu′(t) as t0 f -∞

nan’s description of random point processes in terms of product densities.26 This procedure is general; it is not limited to the discrete systems with a finite number of states considered here. A more detailed account of this approach will be presented elsewhere. Rather than analyzing a single system we introduce an ensemble of a large number n of systems having different composition vectors N1, N2, .... Denoting by nN(t) the number of systems characterized by the composition vector N at time t, we have

n)

nN(t) ) constant ∑ ∀N

(3.1)

Although the total number of systems (replicas) in the ensemble is assumed constant, the numbers nN(t) are fluctuating quantities. For n f ∞ we expect that these fluctuations are negligible and that we can approximate the probability PN(t) by a relative frequency:

PN(t) = nN(t)/n as n f ∞

(3.2)

otherwise PN(t) should be interpreted as an average rather than an instantaneous relative frequency:

〈PN(t)〉 ) 〈nN(t)〉/n

(2.11)

(3.3)

As Pu(t) and Pu′(t) are arbitrary, it follows that for t0 f -∞ all solutions of the master equation (2.1) tend to the same vector:

Hence we need to analyze the problem of the fluctuations of the components of the vector

lim Pu(t) ) P*(t) as t0 f -∞

n ) (nN1, nN2, ...)

(2.12)

where the vector P*(t) is independent of the initial conditions Pu(t0), u ) 1, 2, ... corresponding to the different transient solutions Pu(t), u ) 1, 2, .... In the following we shall use the term “normal solution of the master equation” for P*(t). By examining the formal solution (2.7) of the master equation, we notice that (2.12) implies that

of the numbers of systems of all types. We introduce the probability P(n;t) that at time t the ensemble is formed of nN1 systems of type N1, nN2 systems of type N2, etc. The time evolution of this probability is given by a master equation similar to (2.1):

dP(n,t)

GNN′(t|t0)-∞) ) P* N(t) ) independent of N′ (2.13)

) dt

By making the choice (u′) P(u) N (t) ) PN(t); PN (t) ) P* N(t)

(2.14)

the H function (2.8) becomes

h(t) )

PN(t) ln[PN(t)/P*N(t)] ∑ N

(3.4)

{W(n′fn;t) P(n′,t) - W(nfn′;t) P(n,t)} ∑ ∀n′

(3.5)

where the rates W(nfn′;t) are equal to

W(nfn′;t) )

∑

∀N,N′*N

δn′(...,nN + 1,...,nN′ - 1,...)nNWNN′(t)

(3.6)

(2.15) and P(n,t) fulfills the initial condition

which leads to

P(n;t)t0) ) P0(n)

h > 0, dh/dt < 0 for PN(t) * P* N(t)

(2.16)

h ) 0, dh/dt ) 0 for PN(t) ) P* N(t)

(2.17)

Like (2.9) and (2.10), (2.16) and (2.17) also show that for t0 f -∞ any solution of the master equation (2.1) tends to P* N(t). The physical significance of this result is very simple: after a sufficiently large time, a Markovian system with connected states described by (2.1) forgets its past. 3. Physical Interpretation of the Stochastic H Function To clarify the physical significance of the stochastic H function, we introduce a statistical ensemble representation of the stochastic process described by the Markovian master equation (2.1). Such a representation stems from Ramakrish-

(3.7)

The general solution of (3.5) may be expressed in a closed form in terms of the Green function G(t|t0) of (2.1). We have (see Appendix B)

P(n;t) ) B

{ [ B

∑∑‚‚‚‚‚‚‚‚‚‚‚∑∏ ∏ β)1 u)1 n0

∀muβ g 0

nN0 β0!

]}

(GN0βNu(t|t0))muβ muβ!

P0(n0) (3.8)

where B is the number of possible composition vectors and the sums over muβ obey the constraints

nNu )

∑β muβ;

nN0 0β )

∑u muβ

(3.9)



For t0 f -∞ the Green function G(t|t0) tends toward a form independent of the initial composition vector N0, and the sums over n0 and muβ from (3.9) can be evaluated exactly. After some calculations we obtain

(Lyapunov) function of a stochastic equation, and (2) the asymptotic regime is generally not stationary but rather timedependent. 4. Variational Principle for the Most Probable Path

lim P(n,t) ) P*(n,t) ) independent of P0(n0)

t0f-∞

(3.10a)

where P*(n,t) is a multinomial distribution:

P*(n,t) )

n!

{P*N }n ∏ u

Nu

nN ! ∏ u

u

(3.10b)

u

Thus, the ensemble probability distribution P(n,t) also tends toward a “normal” form that is independent of the past of the system. For large values of n the factorials in (3.10b) can be approximated by means of Stirling’s formula, resulting in

gN ]-1/2 × ∏ u exp[-n∑gN ln(gN /P*N )] u

P*(n,t) = (2πn)-(B-1)/2[

u

u

u

u

as n f ∞ (3.11)

where

gN ) nN/n

gN ln[gN/P* ∑ N(t)] N

(3.13)

is the information gain (Kullback information32,33) obtained by observing the frequency distribution gN for a system whose theoretical probability distribution is given by the normal form P* N(t). In particular if gN is the same as a transient solution PN(t) of the master equation (2.1), that is, a solution that has not yet reached the normal form P*N(t), we have

h(t) ) T[gN)PN(t),P*N(t)]

∑i ν+ij Ai a ∑i ν-ij Ai,

(3.14)

where h(t) is the H function defined by (2.15). We arrive at the following physical interpretation of the H function (2.15): it is both a stochastic potential characterizing the asymptotic form P*(n,t) of the probability P(n,t) of the statistical ensemble fluctuations and the information gain obtained by observing a certain distribution gN ) PN(t) for a system whose assumed (a priori) probability distribution has the normal form P*N(t). Recalling an interpretation of the information gain as a distance between two distributions, it follows that h(t) is a measure of the distance of a solution PN(t) of the master equation (2.1) from its normal form P*N(t). This interpretation is similar to the one ascribed within the framework of the theory of Ross, Hunt, and Hunt20-24 for the exponent of the stationary probability distribution of a chemical reaction system. Similar interpretations have been suggested in nonequilibrium thermodynamics in connection with the solutions of Fokker-Planck equations.19 In comparison with these interpretations presented in the literature, our approach has two new features: (1) It provides an interpretation for the H-

j ) 1, 2, ...

(4.1)

where Ai, i ) 1, 2, ... are the notations for the different molecules making up the system and ν( ij are stoichiometric coefficients. We assume the existence of local equilibrium and the validity of the mass-action law. The forward and backward reaction rates r( j corresponding to the chemical reactions (4.1) are given by

(3.12)

is the relative frequency in the statistical ensemble of the systems with the composition vector N. From the definition (3.12) of gN we notice that the proportionality coefficient of the total number n of systems in the ensemble in the exponential term in (3.11)

T[gN,P* N(t)] )

To make a connection between the H function approach presented in the last two sections and the thermodynamic and stochastic theory of nonequilibrium processes of Ross, Hunt, and Hunt20-24 (RHH), we investigate the connections between the stochastic asymptotic regime described by the normal solution P* N(t) of the chemical master equation (2.1) and the large time behavior of the system described by the deterministic kinetic equations. The analogies with the RHH theory are further investigated in the following sections. We consider the thermodynamic limit for an ideal homogeneous chemical system in which the following reactions occur:

( r( j ) Vkj

∏i (Ni/V)ν ; ( ij

j ) 1, 2, ... as V f ∞ (4.2)

where Ni is the number of Ai molecules, V is the volume of the system, and k+ j , kj are forward and backward rate coefficients, respectively. By assuming that the stochastic as well as the deterministic behavior of the system is described by (4.2), it follows that the rates WNN′(t) from the stochastic master equation (2.1) are given by

WNN′(t) )

( ν δN′(N+f )r( ∑ j (N,t) ) ∑δN′(N+f )Vkj (t)∏(Ni/V) j;( j;( i

( ij

j

j

(4.3)

where

fj ) (fij)i)1,2,...

(4.4)

is the vector of the net stoichiometric coefficients + fij ) νij - νij , i ) 1, 2, ...

(4.5)

attached to the jth reaction. Usually in the RHH theory a reaction system is characterized by the concentrations of the components Ai:

ci ) Ni/V

(4.6)

However, the system considered here is closed so that the use of the reaction extents is more appropriate; the reaction extents allow us to get rid of the stoichiometric constraints characteristic of closed systems. The reaction extents depend on a reference value of the composition vector. As we are interested in the way the system forgets its past, we use as a reference composition vector the equilibrium composition vector Neq(t0) corresponding to a given time rather than the actual initial composition vector N(t0). By Neq(t0) we mean the composition vector that satisfies the conditions of chemical equilibrium (that is, of detailed balance) corresponding to the values of the rate eq coefficients k+ j , kj at time t. It follows that N (t0) is the


Vlad and Ross

solution of the equilibrium equations attached to the chemical reactions (4.1):

By combining (4.2) and (4.11)-(4.15), we obtain

ci(t) ) ai(t) +

eq eq Fj(t) ) r+ j (N (t),t) ) rj (N (t),t); j ) 1, 2, ... (4.7)

where Fj(t) is the common equilibrium value of the forward and the backward reaction rates corresponding to the jth reaction. By using (4.2) we get

∏i (ai)ν

Fj(t) ) k+ j (t)

+ ij

) kj (t)

Kj(t) )

∏i (ai)ν , ij

k+ j (t)/kj (t)

)

j ) 1, 2, ...

∏i (ai)

fij

(4.8) (4.9)

where

ai(t) ) Neq i (t)/V, i ) 1, 2, ...

(4.10)

are equilibrium concentrations and Kj(t) are equilibrium constants. Together with the conditions of mass conservation, (4.9) determines completely the functions ai(t). A theorem proved by Shear15 states that for ideal gases or solutions for a given ( set of rates k( 1 , k2 , ... there is only one equilibrium composition eq vector N (t0). The uniqueness of the equilibrium vector is necessary for the further development of the theory. As in our ( case the rates k( 1 , k2 , ... are time-dependent, at each moment we have a different set of equilibrium concentrations a1(t), a2(t), .... We emphasize that these equilibrium concentrations a1(t), a2(t), ... are only reference values which generally do not describe the actual state of the system. In terms of Neq(t0) we can introduce the extensive reaction extents Ξj(t), j ) 1, 2, ..., by means of the balance equations:

ci(t) )

Neq i (t0)

+

∑j fijΞj(t)

Γ( j )

∏i [ci(t)/ai(t)]ν

( ij

(4.13)

(4.14)

we come to

∑j fijRj(t),

i ) 1, 2, ...

( ij

(4.17a)

(4.17b)

We consider that the chemical system studied is very large and introduce a limit of the thermodynamic type, for which both the volume V of the system and the extensive reaction extents Ξj tend to infinity, but the intensive reaction extents ξj remain constant, and we express the state probability of the system in terms of the intensive reaction extents ξj:

PN(t) S P(ξ,t) dξ, as V f ∞ with ξ ) constant (4.18) By using a method suggested by Kubo, Matsuo, and Kitahara27,28 in this limit, we can write the master equation (2.1) in a form similar to a Schro¨dinger equation:

V-1

∂P(ξ,t) ) -H ˆ [t,ξ,V-1∇ξ]P(ξ,t) ∂t

{[

(4.11)

We assume that the chemical system is operated at constant volume so that the passage from the extensive extents Ξj(t), j ) 1, 2, ..., to the intensive extents ξj(t), j ) 1, 2, ..., and vice versa is straightforward. To express the reaction rates r( j in terms of the intensive reaction extents ξj(t), we introduce the equilibrium intensive extents Rj(t), j ) 1, 2, ..., which correspond to the equilibrium concentrations a1(t), a2(t), .... By using the convention

ai(t) ) ai(t0) +

∏i {1 + ∑j fij[ξj(t) - Rj(t)]}ν

(4.19a)

where the “Hamiltonian operator” H ˆ is given by

From (4.11)-(4.12) we have

Rj(t0) ) 0, j ) 1, 2, ...

)

a(t) ) [ai(t)]; R(t) ) [Rj(t)]; ξ(t) ) [ξ(t)]

(4.15)

∑j Fjj(t)

(

1 - exp -V-1

[ (

∂

∂

)]

∂ξj

)]

(Γ+ j (ξ,t)...) +

1 - exp +V-1 (Γj (ξ,t)...) ∂ξj

(4.12)

∑j fijξj

(4.16b)

where

We can also introduce the intensive extents ξj(t), j ) 1, 2, ...:

ci(t) ) ai(t0) +

(4.16a)

( r( j ) Fj(t)Γj [a(t),r(t),ξ(t)]

H ˆ [t,ξ,∇ξ]... )

ξj(t) ) Ξj(t)/V

∑j fij[ξj(t) - Rj(t)]

}

(4.19b)

By applying a variant of the Kubo-Matsuo-Kitahara approximation developed by Vlad and Ross,23 (4.19a) can be used for deriving a path integral representation for the state probability:

P(ξ,t) )

∫G(ξ,t|ξ0,t0) P(ξ0,t0) dξ0

(4.20)

where the Green function G(ξ,t|ξ0,t0) is equal to

G(ξ,t|ξ0,t0) )

∫∫D[ξ(t′)] exp{V∫t tL[ξ(t′),dξ(t′)/dt′,t′]dt′} 0

(4.21) Here D[ξ(t′)] is an integration measure in the space of functions ξ(t′); in the following we do not need its explicit expression. L[ξ(t′),dξ(t′)/dt′,t′] is a “chemical Lagrangian” given by

L[ξ(t′),dξ(t′)/dt′,t′] )

∑j Lj[ξ(t′),dξ(t′)/dt′,t′]

(4.22a)



where

Lj[ξ(t′),dξ(t′)/dt′,t′] ) dξj(t′) 1 ln dt′ 2F (t′) Γ-(ξ(t′),t′)

{

j

j

[[( ) ] ]} ] dξj(t′) 2 + dt′

1/2

4(Fj(t′))2Γ+ j (ξ(t′),t′) Γj (ξ(t′),t′)

[( )

The absolute maximum of the normal solution P*(ξ,t) of the master equation (2.1) corresponds to the deterministic normal path ξ*(t), which is the solution of the deterministic kinetic equations (4.25) and is independent of the initial state of the system. Since for t0 f -∞ the state probability P(ξ,t) tends toward the normal form P*(ξ,t), it follows that the transient optimal deterministic paths corresponding to the maximum of the state probability P(ξ,t), which also fulfill the deterministic kinetic equations (4.25), also tend toward the normal deterministic path ξ*(t).

-

dξj(t′) dt′

dξj(t′) 2 + 4(Fj(t′))2Γ+ j (ξ(t′),t′) Γj (ξ(t′),t′) dt′

+

1/2

5. Deterministic versus Average Paths

Fj(t′)[Γ+ j (ξ(t′),t′) + Γj (ξ(t′),t′)] (4.22b)

are individual chemical Lagrangians attached to the different reactions (4.1). The main steps of the derivation of (4.20)(4.22) are outlined in Appendix C. For a similar derivation for one-variable open systems with time-independent rate coefficients see Vlad and Ross.23 Equations 4.20-4.22 give a good approximation of the solutions of the master equations in the thermodynamic limit. The exact evaluation of the path integral (4.21) is impossible; however some conclusions can be drawn without calculating it. The chemical Lagrangian (4.22) has the following important property:

L[ξ(t′),dξ(t′)/dt′,t′] < 0 if for at least one j

The results derived in the preceding section lead to a variational principle for the most probable (deterministic) path:

-L[ξ(t′),dξ(t′)/dt′,t′] ) minimum

(5.1)

In particular, if t0 f -∞, (5.1) reduces to an equation for the most probable normal (deterministic) path:

-L[ξ*(t′),dξ*(t′)/dt′,t′] ) minimum

(5.2)

The aim of this section is to derive a similar variational principle for the average normal path corresponding to the normal stochastic solution P*(ξ,t) of the master equation (2.1). We note that the integrand in (4.27) may be interpreted as the probability density functional of a normal random path:

R*[ξ(t′)] D[ξ(t′)] ) D[ξ(t′)] ×

dξ(t′)/dt′ *

Fj(t′)[Γ+ j (ξ(t′),t′)

-

Γj (ξ(t′),t′)]

exp{V

(4.23)

and

∫-∞+∞L[ξ(t′),dξ(t′)/dt′,t′]dt′}

(5.3a)

with the normalization condition

L[ξ(t′),dξ(t′)/dt′,t′] ) 0

∫∫ R*[ξ(t′)] D[ξ(t′)] ) 1

if for all j dξ(t′)/dt′ ) Fj(t′)[Γ+ j (ξ(t′),t′) - Γj (ξ(t′),t′)]

(4.24)

(see Appendix D). It follows that the chemical Lagrangian L[ξ(t′),dξ(t′)/dt′,t′] has a maximum value for the deterministic path which obeys the kinetic equations dξ(t′)/dt′ ) Fj(t′)[Γ+ j (ξ(t′),t′) - Γj (ξ(t′),t′)], j ) 1, 2, ...

As we are interested in the evolution of a path from t0 f -∞ to t f +∞, we have replaced the superior integration limit in the exponent of (5.3a) by +∞. Similarly, in the normalization condition (5.3b) the path integral is taken with respect to all possible random trajectories ξ ) ξ(t′), -∞ < t′ < +∞, defined along the whole time axis. In terms of the probability density functional R*[ξ(t′)] D[ξ(t′)] we can define the cumulant generating functional:

(4.25) This maximum value of the chemical Lagrangian corresponds to an absolute maximum of the state probability P(ξ,t); thus the deterministic path is the most probable path. Another consequence of the path integral approach is that as t0 f -∞, the solutions of the kinetic equations (4.25) also tend toward a deterministic normal form:

ξ* ) ξ*(t)

(4.26)

where the vector ξ*(t) is independent of the initial state of the system. To prove that we express the normal solution of the master equation (2.1) as a path integral,

P*(ξ,t) )

∫∫ D[ξ(t′)] exp{V∫-∞t L[ξ(t′),dξ(t′)/dt′,t′]dt′} (4.27)

(5.3b)

∫-∞+∞[∑ηj(t′) ξj(t′) dt′]}〉} )

χ[η(t′)] ) ln{〈exp{

j

ln{

∫∫ R*[ξ(t′)] D[ξ(t′)] exp{∫-∞+∞[VL[ξ(t′),dξ(t′)/dt′,t′] + ∑j ηj(t′) ξj(t′)]dt′}}

(5.4)

where η(t) ) [ηj(t)], j ) 1, 2, ..., is a vectorial test function whose components ηj(t), j ) 1, 2, ..., have nonpositive real parts and the brackets 〈 〉 denote the averaging with respect to all possible paths ξ ) ξ(t′), -∞ < t′ < +∞. The different cumulants of the state vector ξ(t) may be computed by evaluating the functional derivatives of the cumulant generating functional (5.4) for η(t′) ) 0. We have

χ[η(t′))0] ) 0, 〈ξj(t1)〉 ) {δχ[η(t′))0]/δηj(t1)} 〈∆ξj1(t1) ∆ξj2(t2)〉 ) {δ2χ[η(t′))0]/δηj1(t1) δηj2(t2)}, etc.

we have

P(ξ,t) ) P*(ξ,t) as t0 f -∞

(4.28)

(5.5)


Vlad and Ross

In addition to the “theoretical” (a priori) probability density functional of a normal random path, R*[ξ(t′)] D[ξ(t′)], we introduce an “observed” (a posteriori) frequency density functional,

Π*[ξ(t′)] D[ξ(t′)] with

∫∫

{

}

Π*[ξ(t′)] D[ξ(t′)] (5.7) R*[ξ(t′)]

K[Π*[ξ(t′)]] is the information gain obtained as the result of the observation of the “a posteriori” frequency density functional Π*[ξ(t′)] D[ξ(t′)] for a system with a theoretical normal probability density functional R*[ξ(t′)] D[ξ(t′)]. To derive a variational principle for the average path, we are looking for the functional Π*[ξ(t′)] which extremizes the information gain K[Π*[ξ(t′)]] under the constraint that Π*[ξ(t′)] is normalized,

∫∫ Π*[ξ(t′)]D[ξ(t′)] ) 1

By combining (5.11)-(5.12), we notice that F[µ(t′)] is in fact a functional Legendre transform with changed sign of the cumulant generating functional χ[η(t′)]. It follows that

ηj(t) )

δ F[µ(t′)] , j ) 1, 2, ... δµj(t)

t ) arbitrary (5.9)

In particular, for η ) 0

δ F[µ(t′)] δµj(t)

|

) 0, j ) 1, 2, ...

From (5.15) and (5.16) we observe that the theoretical average path 〈ξ(t)〉 may be computed by extremizing the extremal information gain F[µ(t′)]:

F[µ(t′)] ) extremum for µ(t) ) 〈ξ(t)〉observed ) 〈ξ(t)〉 (5.18)

Λj1j2(t1,t2) )

|

δ2 F[µ(t′)] δµj1(t1) δµj2(t2)

∫-∞+∞∑ηj(t′) ξj(t′) dt′}

(5.10)

j

δµj2(t2)

〈∆ξj1(t1) ∆ξj2(t2)〉 )

|

δ2 χ[η(t′)] δηj1(t1) δηj2(t2)

where the Lagrange multipliers ηj(t′) are determined from the constraint (5.9). We obtain

δχ[η(t′)] , j ) 1, 2, ... δηj(t)

(5.19) µ(t′))〈η(t′)〉

(5.12)

δηj2(t2)

|

(5.20) η(t′))0

As the condition µ(t) ) 〈ξ(t)〉observed ) 〈ξ(t)〉 for F[µ(t′)] corresponds to the condition η(t′) ) 0 for χ[η(t′)], from (5.19)(5.20) we come to

∑j′ ∫Λj j′(t1,t′)〈∆ξj′(t′) ∆ξj (t2)〉dt′ ) δj j δ(t1 - t2) 1

where F[µ(t′)] is the extremal value of K[Π*[ξ(t′)]] consistent with the constraint that

) η(t′))0

δµj1[η(t′);t1]

ηj(t′) ξj(t′) dt′ (5.11)

j

µj(t) )

|

δηj1[µ(t′);t1]

F[µ(t′)] ) Kextremal[µ(t′)] )

∫-∞ ∑

) µ(t′))〈η(t′)〉

The cumulants of second order of the intensive reaction extents ξj(t), j ) 1, 2, ..., are equal to

Π*extremal[ξ(t′),µ(t′)] ) exp[-χ[η(t′)]] R*[ξ(t′)] ×

+∞

(5.17)

η(t′))0

By extremizing K[Π*[ξ(t′)]] with the constraints (5.8)-(5.9), we get the following expressions for the “a posteriori” frequency density functional Π*[ξ(t′)] D[ξ(t′)] and for the extremal value of the information gain K[Π*[ξ(t′)]]:

-χ[η(t′)] +

(5.16)

To see whether the extremum (5.18) is a minimum or a maximum, we evaluate the second functional derivative of F[µ(t′)] for µ(t) ) 〈ξ(t)〉observed ) 〈ξ(t)〉:

∫∫ ξ(t′) Π*[ξ(t′)] D[ξ(t′)];

exp{

q.e.d. (5.15)

(5.8)

and that the average path is given by a known function

µ(t) ) 〈ξ(t)〉observed )

∫-∞+∞∑ηj(t′) ξj(t′) dt′|η(t′))0 ) 0, j

Π*[ξ(t′)] D[ξ(t′)] ) 1

In terms of the “a posteriori” frequency density functional Π*[ξ(t′)] D[ξ(t′)] we define the information gain (Kullback information):

∫∫ Π*[ξ(t′)] ln

F[µ(t′))〈ξ(t′)〉] ) F[µ(t′)] ) -χ[η(t′)] +

(5.6)

K[Π*[ξ(t′)]] )

To prove this we combine (5.11)-(5.12) and (5.5)-(5.6), resulting in

2

12

(5.21)

It is the extremal value of the information gain obtained as a result of observing the average path µ(t) ) 〈ξ(t)〉observed. We expect that when the observed average path is the same as the theoretical one computed from the normal probability density functional R*[ξ(t′)] D[ξ(t′)], that is, when

Since by definition the cumulants of the second order 〈∆ξj1(t1) ∆ξj2(t2)〉 are non-negative, from (5.21) it follows that the same is true for the functional derivatives of second order of the functional F[µ(t′)]; that is, for the functions Λj1j2(t1,t2) defined by (5.19), and therefore the extremum of F[µ(t′)] is a minimum. From the considerations presented above it turns out that the average path corresponding to the normal stochastic solution P*(ξ,t) of the master equation (2.1) obeys a variational principle similar to (5.2) for the most probable normal path.

µ(t) ) 〈ξ(t)〉observed ) 〈ξ(t)〉

F[µ(t′)] ) minimum for µ(t′) ) 〈ξ(t)〉

µ(t) ) 〈ξ(t)〉observed

(5.13)

(5.14a)

6. Hamilton-Jacobi Approach

then this gain of information is equal to zero:

F[µ(t′) ) 〈ξ(t′)〉] ) 0

(5.22)

(5.14b)

In this section we develop a Hamilton-Jacobi approach for the study of the normal regime. This approach is used in the



next section for developing a thermodynamic descriprtion of the long time behavior of chemical systems with time-dependent rate coefficients. The method is based on the evaluation of the path integrals introduced in section 4 by means of an eikonal (Laplace) approximation. Both the transient and the normal solutions can be studied by using this method. We analyze here the normal solution only. We express the exponent in (4.27) in terms of a “chemical action”:

J*(ξ,t) )

∫-∞t L[ξ(t′),dξ(t′)/dt′,t′]dt′

(6.1)

and look for the paths ξj(t′) which maximize the chemical action J*(ξ,t); these paths obey the variational condition

∫-∞t L[ξ(t′),dξ(t′)/dt′,t′]dt′|ξ(t′))ξj(t′) ) 0

δJ*(ξ,t) ) δ

(6.2)

Equation 6.2 leads to the Euler-Lagrange equations

{[

]}

dξj(t′) ∂ ,t′ L ξj(t′), ∂ξjj dt′

)

d dt

{[

[[

]]

}

dξj(t′) ∂ ,t′ , L ξj(t′), dt′ dξj(t′) ∂ dt′ j ) 1, 2, ... (6.3)

]

where the overbar denotes the optimal paths, each one corresponding to a different initial condition. We approximate the chemical action J*(ξ,t) by the first three terms of its functional Taylor expansion around ξ(t′) ) ξj(t′):

j,t) + δJ*(ξj,t) + 1/2δ2J*(ξj,t) + ... = J*(ξ,t) ) J* extremal(ξ 1

2

Equations 4.27, 6.1, and 6.4 lead to the following expression for the normal probability density P*(ξ,t):

(6.5)

The pre-exponential proportionality factor in (6.4) is given, up to a proportionality constant, by the Gaussian path integral

∫∫ exp[V2 δ2J*(ξj,t)] D[ξ(t′)] which can be computed analytically. In the thermodynamic limit the main contribution to the normal state probability P*(ξ,t) comes from the exponential term (6.5), and the value of the path integral is practically constant. By inserting (6.5) into (4.20) and keeping the dominant terms in V, we get a Hamilton-Jacobi (HJ) equation for h,t): J* extremal(ξ

∂ J*(ξ,t) + H[t,ξ,∇ξJ*(ξ,t)] ) 0 ∂t

(6.6a)

∑j

(6.7)

dp(t)/dt ) -∇ξH[t,ξ(t),p(t)]

(6.8)

In particular, from (6.7) applied for

p ) ∇ξJ*(ξ,t) ) 0

(6.9)

we recover the deterministic kinetic equations (4.25). By using (6.1) we can express the chemical action J*(ξ,t) as a time integral of the chemical Lagrangian, L[ξ(t′),dξ(t′)/dt′,t′]: ξ(t) ) ξ ∫-∞t L[ξ(t′),dξ(t′)/dt′,t′]dt′|extremal

(6.10)

where the time integral is evaluated for the optimal paths obeying the Euler-Lagrange equations (6.23) or the Hamiltonian equations (6.7)-(6.8). Equation 6.10 is similar to the wellknown Onsager-Machlup relationship between the entropy and the thermodynamic Lagrangian for linear thermodynamic processes near equilibrium.35 The RHH theory20-24 suggests that stochastic potentials similar to J*(ξ,t) may be used for constructing H functions (i.e. Lyapunov functions) for the deterministic kinetic equations (4.25). By making an analogy with the case of open systems far from equilibrium with time-independent rate coefficients, we assume that such a Lyapunov function is given by the difference

∆J*(ξ,t) ) J*(ξ,t) - J*(ξ*,t)

(6.11)

where ξ* is the normal solution of the deterministic kinetic equations. In order that ∆J*(ξ,t) be a Lyapunov function, we should have

∆J*(ξ,t) < 0; d(∆J*(ξ,t))/dt > 0 for ξ * ξ* (6.12a,b) ∆J*(ξ,t) ) 0; d(∆J*(ξ,t))/dt ) 0 for ξ ) ξ* (6.13a,b)

where

H[t,ξ,p] )

dξ(t)/dt ) ∇pH[t,ξ(t),p(t)]

J*(ξ,t) )

J* j,t) + /2δ J*(ξj,t) (6.4) extremal(ξ

P*(ξ,t) ∝ exp{VJ* j,t)} extremal(ξ

intensive reaction extents. Since in all that follows all paths are extremal paths obeying the Euler-Lagrange equations (4.3), we now drop the overbars and the subscript “extremal”. Although less complicated than the master equations (2.1) and (4.27), the Hamilton-Jacobi equation (6.6a) is still very difficult to solve. The main advantage of the HJ equation (6.6a) is that it allows the development of an analogy between the chemical systems investigated in this article and the theory of Hamiltonian systems from analytical mechanics.34 The use of this analogy makes it possible to investigate some features of the stochastic and deterministic normal regimes without actually solving (6.6). From the theory of Hamilton-Jacobi equations it is known that solving (6.6.a) is equivalent to solving the Euler-Lagrange equations (6.3). An alternative approach is based on the integration of the characteristic system of (6.6a) (the Hamiltonian equations):

Fjj(t){[1 - exp(-p)]Γ+ j (ξ,t) + [1 - exp(+p)]Γj (ξ,t)} (6.6b)

By using the Hamilton-Jacobi equation (6.6), it is easy to prove that (6.12) and (6.13) are valid. Equations 6.12a and 6.13a are due to the fact that the normal deterministic solution ξ*(t′) is the most probable path corresponding to P*(ξ,t). Thus we have

P*(ξ,t) < P*(ξ*,t) for ξ * ξ*

(6.14)

is a Hamiltonian function attached to the Hamiltonian operator (4.19b) and

from which, by using (6.5), we come to (6.12a) and (6.13a).

p ) ∇ξJ*(ξ,t)

P*(ξ,t) ) P*(ξ*,t) for ξ ) ξ*

(6.6c)

is a generalized impulse vector conjugate to the vector ξ of the

(6.15)

The proof of (6.12b) and (6.13b) is a little bit more complicated;


Vlad and Ross

it proceeds in two steps. First we evaluate the derivatives of J*(ξ,t) with respect to ξj, j ) 1, 2, ..., and t for ξ ) ξ*. As P*(ξ,t) has a maximum for ξ ) ξ*, we have

[∇ξJ*(ξ,t)]* ) 0

(∂J*/∂t)* ) 0

(6.17)

Now we can evaluate the total time derivative of ∆J*(ξ,t) by assuming that the time evolution of the state vector ξ is given by the deterministic evolution equations (4.25):

) dt

∑j

( )( ) ∂J* dξj ∂ξj

)

+

FjΓ( ∑ j J;(

∂t

{ ( ) [ ( )]} 1(

∂J* ∂ξj

- exp -

∂J* ∂ξj

(6.18)

Here we have used the conditions (6.16) and (6.17) to evaluate the total time derivative of ∆J*(ξ,t), and (4.25) and (6.6) to evaluate the time derivative of the state vector ξ and the partial time derivative of ∆J*(ξ,t), respectively. Taking into account that for -∞ > a > +∞ we have

∆Ξj(t) ) Ξj(t) - Ξ*j

and taking into account that Fj and Γ( j are positive by definition from (6.18), we come to (6.12b) and (6.13b). The use of the Lyapunov function ∆J*(ξ,t) is not absolutely necessary. We have shown that the solutions of the kinetic equations (4.25) tend toward a normal form ξ*(t) by using the stochastic H theorem presented in section 2. However, the deterministic H theorem presented here is useful because it may serve as a basis for developing a generalized thermodynamic formalism for the time-dependent normal regime. This is the subject of the next section. 7. Thermodynamic Analysis In this section we try to develop a nonequilibrium thermodynamic formalism for chemical systems with time-dependent rate coefficients by starting from (6.5) for the probability of fluctuations P*(ξ,t) corresponding to the stochastic normal regime. The first step in the development of such a formalism is the evaluation of the moments of the fluctuating variables. For small fluctuations the cumulants of second order of the reaction extents ξ1, ξ2, ... are given by

)]

{ [( ) ]} { ∑( )

P*(ξ,t) = (2π)-R/2V-R/2 det

(7.1)

j1j2

(7.2)

is the displacement of the intensive extent ξj(t) of the jth reaction (4.1) from its deterministic normal value ξ* j(t). By expressing (7.1) in terms of the extensive reaction extents,

Ξj(t) ) Vξj(t); Ξ*j(t) ) Vξ* j(t)

-∂2J*

-1/2

∂ξu1 ∂ξu2

×

*

}

(7.6)

Equation 7.6 results by expanding the function J*(ξ,t) for t ) constant in a Taylor series around the normal solution ξ(t) and by keeping the first three terms:

J*(ξ,t) = J*(ξ*,t) +

1

∑

( ) ∂2J*

2 j1,j2 ∂ξj1 ∂ξj2

(ξj1 - ξ* j1)(ξj2 - ξ* j2)

*

(7.7a)

[( ) ] -∂2J* ∂ξj1 ∂ξj2

(7.7b)

*

should be positive definite. By examining (7.1), we notice that the positive-definiteness of the matrix M is consistent with the fact that the square form of the covariance matrix, M-1, should be also positive definite for nonvanishing and finite fluctuations. The thermodynamic conditions for the existence and stability of a certain state are usually expressed in terms of the total differentials of a certain extensive state function which is proportional to the logarithm of the probability of fluctuations. In our case such an extensive state function can be derived by taking the product of the chemical action J*(ξ,t) with the volume V of the system:

Ψ*(Ξ,V,t) ) VJ*(ξ)Ξ/V,V,t)

(7.8)

According to the definition (7.8) the function Ψ*(Ξ,V,t) is extensive; that is, it is a first-order homogeneous function of Ξj and V. To prove that, we note that in the Hamilton-Jacobi equation (6.6a) Ξj, j ) 1, 2, ..., and V enter in the form of the combinations ξj ) Ξj/V, j ) 1, 2, .... By taking this observation into account from (7.8), we get the following scaling equation of the Euler type:

-1

where R is the number of reactions (4.1) and

∆ξj(t) ) ξj(t) - ξ* j(t)

(7.5)

Equations 7.1-7.4 can be obtained by approximating (6.5) by a Gaussian:

M)

1 + a - exp(-a) ) 0 for a ) 0

[(

(7.4)

* j1j2

where

and

-∂2J* ∂ξu1 ∂ξu2

)]

-1

As the probability density P*(ξ,t) has a maximum for ξ(t) ) ξ*(t), it follows that for nonvanishing fluctuations the matrix

1 + a - exp(-a) > 0 for a * 0

〈∆ξj1(t) ∆ξj2(t)〉 ) V-R

[(

-∂2J* ∂ξu1 ∂ξu2

1 -∂2J* exp V (ξj - ξ* j1)(ξj2 - ξ* j2) 2 j1,j2 ∂ξj1 ∂ξj2 * 1

∂J*

dt

〈∆Ξj1(t) ∆Ξj2(t)〉 ) VR

(6.16)

Here and in the following the subscript * means that the derivatives are evaluated for ξ ) ξ*. From (6.6) and (6.16) we get

d∆J*

we have

(7.3)

λΨ*(Ξ,V,t) ) Ψ*(λΞ,λV,t), λ > 0

(7.9)

which describes the property of extensivity of the Ψ*(Ξ,V,t) function. Now a difficulty arises. Usually the thermodynamic formalism is used for the analysis of time-independent states. In our case, however, the normal solution ξ*(t) is time-dependent. From this point of view we distinguish two different cases: (a) The time scale for the regression of fluctuations, which is the same as the kinetic time scale, is smaller than the time scale for the variation of the normal solution. In this case the existence and stability of the normal solution ξ*(t) can be



discussed by evaluating the differentials of the potential Ψ*(Ξ,V,t) at a given time. (b) The time scale for the regression of fluctuations has the same order of magnitude as the time scale of variation of the normal solution ξ*(t). In this case, in order to discuss the existence and stability of the normal regime, we should differentiate the potential Ψ*(Ξ,V,t) with respect to the reaction extents and with respect to the time. (a) Fast Fluctuations. We have

dΨ*|t)constant )

∑j

pj dΞj

(7.10)

where

pj(ξ,t) )

∂Ψ* ∂J* ) , j ) 1, 2, ... ∂Ξj ∂ξj

(7.11)

( ) ( )

∂Ψ* ∂J* ) ) 0, j ) 1, 2, ... ∂Ξj * ∂ξj *

dξj

∑j pj dt

(7.13)

Equation 7.13 is similar to the expression for the entropy production in nonequilibrium thermodynamics. The condition for the existence of the deterministic normal regime is

dΨ*|t)constant ) 0 for ξ ) ξ*

(7.14a)

which corresponds to

Ψ ˙ *)0

(7.14b)

(7.15)

The stability of the normal regime is discussed in a similar way. The second differential of Ψ*(Ξ,V,t) at t ) constant is equal to

( )

∑ j ,j ∂ξ 1 2

j1

∂ξj2

Ξ

)V

(∂J*∂t )

(7.19)

ξ

is the rate of time variation of the stochastic potential Ψ*(Ξ,V,t) at constant Ξ. The rate Λ can be expressed in terms of the generalized reaction affinities pj. By using the Hamilton-Jacobi equation (6.6) and the expression of the Hamiltonian function H we come to

∑j Fj[Γ+j (1 - exp(-pj)) + Γ-j (1 - exp(+pj))]

Λ ) -V

For the normal regime the generalized reaction affinities pj are equal to zero and (7.20) leads to

Λ(ξ)ξ*) ) 0

(7.21)

(see also (4.17)). It follows that the condition for the existence of the normal regime is

dΨ* ) 0 for ξ ) ξ*

(7.22)

The rate of production of the stochastic potential Ψ*(Ξ,V,t) is equal to

Ψ ˙ *=

dΨ* dt

dξj

∑j pj dt + Λ ) V∑FjΓ( j [1 ( pj - exp(-pj)] j;(

)V

dξj1 dξj2

(7.16a)

(M)j1,j2 dξj ∑ j ,j

1

dξj2

(7.16b)

1 2

The stability condition is given by

d2 Ψ*(ξ*,t)|t)constant < 0

that is, it corresponds to a maximum of the potential Ψ*(Ξ,V,t); this is consistent with the fact that the matrix M is positive definite.

(7.24)

The second differential of the stochastic potential Ψ*(Ξ,V,t) is given by

d2 Ψ* )

∂2Ψ*

∑ j ,j ∂Ξ

j1

dΞj1 dΞj2 + 2

∂Ξj2

∂2Ψ*

∑j ∂Ξ ∂t

dΞj dt +

j

∂2Ψ* 2 dt ∂t2 (7.25)

By using the Hamilton-Jacobi equation (6.6) the time derivatives of Ψ*(Ξ,V,t) can be expressed in terms of the derivatives of Ψ*(Ξ,V,t) with respect to the extensive reaction extents ξj, j ) 1, 2, .... The calculations are very complicated in the general case. Fortunately for the analysis of stability it is sufficient to evaluate these derivatives for ξ(t) ) ξ*(t). After the repeated use of the Hamilton-Jacobi equation we obtain (see Appendix E)

( ) ( ) ∂2Ψ*

(7.17)

(7.23)

The last form of (7.23) has been derived by using the expression (7.20) for Λ and the deterministic kinetic equations (4.25). The expression for the dissipation function Ψ ˙ * has a structure similar to the relationship (6.18) for the time derivative of the Lyapunov function ∆J*. From (7.23) we have

1 2

In particular, in the vicinity of the normal solution ξ(t) ) ξ*(t) we have

d2 Ψ*(ξ*,t)|t)constant ) -V

(7.18)

Ψ ˙ *(ξ*ξ*) > 0; Ψ ˙ *(ξ)ξ*) ) 0

Ψ ˙ *>0

d2 Ψ*(ξ,t)|t)constant ) V

(∂Ψ* ∂t )

Λ)

Away from the normal regime we have

∂2J*

∑j pj dΞj + Λdt

where

(7.12)

(see also (6.16)). For pj * 0 the pj(ξ) functions are a measure of the departure of the system from the normal regime. The rate of production of Ψ*(Ξ,V,t) may be approximated by

dΨ* |t)constant ) V Ψ ˙ *= dt

dΨ* )

(7.20)

play the role of generalized reaction affinities. For the deterministic normal regime we have

pj(ξ)ξ*,t) )

(b) Slow Fluctuations. We have

∂Ξj ∂t

*

∂2Ψ* ∂t2

∑j′

)-

*

)V

( ) ( )

rj′

∂2Ψ*

∂Ξj ∂Ξj′ ∂2Ψ*

rjrj′ ∑ ∂Ξ ∂Ξ j,j′ j

(7.26)

*

j′ *

(7.27)


Vlad and Ross

where the net rates r1, r2, ... of (4.1) are given by + rj ) r+ j - rj ) Fj(Γj - Γj )

and of the probability diffusion coefficients

(7.28)

By inserting (7.26)-(7.27) into (7.25), we have

(M)j ,j (dξj ∑ j ,j

d2 Ψ*(ξ*,t) ) -V

1 2

1

1 + Dj ) (r+ + rj ) ) Fj(Γj + Γj ), j ) 1, 2, ... (7.36) 2 j Equations 4.21, 6.6, 7.19, and 7.35-7.36 lead to the following relationship among Λ, V, Dj, rj, and pj:

- rj1dt)(dξj2 - rj2dt)

1 2

(7.29)

We distinguish two different cases. (1) Stability. If at least one of the deterministic equations (4.7) is violated during a fluctuation, then there is at least one j for which dξj * rjdt. As M is positive definite from (7.29), we get the following condition of stability:

d2 Ψ*(ξ*,t) < 0

Λ)V

∑j

{

() } pj

2Dj tanh

r ) 2D tanh

d2 Ψ*(ξ*,t) ) 0 for dξj ) rjdt, j ) 1, 2, ...

(7.31)

By comparing the conditions for the existence and stability of the normal regime for fast and slow fluctuations we notice that they are practically equivalent to each other. The main difference is related to the existence of neutral stability for slow fluctuations. The extensivity property of the stochastic potential Ψ*(Ξ,V,t), expressed by (7.9), leads to a Gibbs-Duhem relationship for the generalized affinities pj, j ) 1, 2, .... By differentiating (7.9) with respect to λ and putting λ ) 1, we have

Ψ* )

∑j pjΞj + Vπ

(7.32)

(7.33)

is an intensive variable similar to the pressure from equilibrium thermodynamics. By differentiating (7.32) at V ) constant and comparing the result with (7.18), we come to

∑j Ξjdpj + Vdπ ) Λdt

(7.37)

Λ (2p) - V sinh(p)

(7.38)

Here the reaction label j ) 1 has been left out. As a final topic of this section we discuss the relationship between the thermodynamic description introduced above and the statistical ensemble description used in section 3. In the thermodynamic limit the normal solution (3.11) of the ensemble master equation (3.5) may be written as

P*(n,t) ∝ exp(Φ*) as n f ∞

(7.39)

where Φ*(ξ,t) is an “ensemble” stochastic potential given by

∫

Φ* ) n P(ξ,t) ln

[ ]

P*(ξ,t) dξ P(ξ,t)

(7.40)

The history-dependent probability density, P(ξ,t) dξ, may also be evaluated by means of the Hamilton-Jacobi approach used for the normal probability density P(ξ*,t) dξ. We have

P(ξ*,t) ∝ exp(VJ) as V f ∞

(7.41)

where the history-dependent chemical action J(ξ,t) obeys the same Hamilton-Jacobi equation as the normal chemical action J*(ξ,t),

where

π ) ∂Ψ*/∂V

- rj sinh(pj)

Equation 7.37 is similar to a fluctuation-dissipation relation for one-variable open chemical systems derived within the framework of RHH theory.23 In the particular case of a single reaction, (7.37) can be written in the classical form of a forceflux relationship:

(7.30)

if for at least one j, dξj * rjdt. Equation 7.30 is similar to (7.17) for fast fluctuations. (2) Neutral Stability. If the fluctuations from the normal regime obey the deterministic kinetic equations dξj ) rjdt for all j ) 1, 2, ..., then (7.29) leads to

2

(7.34)

∂ J(ξ,t) + H[t,ξ,∇ξJ(ξ,t)] ) 0 ∂t

(7.42)

with the difference that now the transient solution J(ξ,t) depends on the initial preparation of the system. After some algebraic transformations (7.40) for the grand canonical stochastic potential Φ*(ξ,t) can be rewritten as

(S*S )

Φ* ) n(〈Ψ*〉 - 〈Ψ〉) + n ln

(7.43)

Equation 7.34 is similar to the Gibbs-Duhem relation from equilibrium thermodynamics. The classical linear nonequilibrium thermodynamics is based on a set of “force-flux” relationships. In our case such relationships should establish a connection between the net reaction rates r1, r2, ... and the corresponding generalized reaction affinities pj, j ) 1, 2, .... Due to the explicit time dependence of the rate coefficients, the derivation of a complete set of force-flux relationships is not possible. An overall force-flux relationship can be however derived from the Hamilton-Jacobi equation (6.6). By expressing the forward and backward reaction rates r( j , j ) 1, 2, ... in terms of the net reaction rates rj, j ) 1, 2, ...,

and S and S* are two normalization factors corresponding to the probability densities P(ξ,t) dξ and P*(ξ,t) dξ, respectively:

+ rj ) r+ j - rj ) Fj(Γj - Γj ), j ) 1, 2, ...

S ) { exp[VJ(ξ,t)]dξ}-1

(7.35)

where the averages are computed by using the history-dependent probability density P(ξ,t) dξ:

〈...〉 )

∫...P(ξ,t) dξ

(7.44)

Ψ(ξ,t) is a stochastic potential corresponding to the transient chemical action J(ξ,t):

Ψ(ξ,t) ) VJ(ξ,t)

∫

(7.45)



and

where

∫

S* ) { exp[VJ*(ξ,t)]dξ}-1

(7.46)

∑i a (t )-∞) + f R(t)

λ(t) ) F(t)

Equation 7.43 shows that the relation between the grand canonical stochastic potential Φ*(ξ,t) and the canonical stochastic potential Ψ(ξ,t) is similar to the one between the Gibbs and Boltzmann’s entropies in equilibrium statistical mechanics. We conclude this section by examining the relationships between the thermodynamic formalism developed in this article and equilibrium and nonequilibrium thermodynamics. At first sight it seems that the conditions for the existence and stability of the normal solution derived in this article have nothing to do with classical thermodynamics because the corresponding equations depend explicitly on time. For example, in the generalized Gibbs-Duhem equation (7.34), the time t is taken as a variable in the space of macroscopic variables. In usual equilibrium and nonequilibrium thermodynamics, time is at most a nonthermodynamic parameter. Despite this difference, the analogies between our formalism and the usual thermodynamics are far from formal. By applying our approach in the particular care when the rate coefficients are time-independent but some constraints prevent the system from reaching equilibrium, our approach reduces to the thermodynamic and stochastic theory of Ross, Hunt, and Hunt. Moreover, if we remove the constraints, the system eventually reaches equilibrium and our description reduces to classical equilibrium thermodynamics. 8. Application to One-Reaction Systems In this section we apply the general theory developed in this article to the particular case of closed chemical systems with one reaction of the type

∑i ν+i Ai a ∑i ν-i Ai

i 0

has been evaluated by means of analytical or numerical integration and is thus known. If in addition we assume that the fluctuations of the reaction extent ξ(t) around the normal solution ξ*(t) of (8.2) are small, then we can also evaluate the normal chemical action J*(ξ,t) and the corresponding normal probability density P*(ξ,t). Under these circumstances we can linearize the stochastic evolution equations in terms of the displacement of the reaction extent ξ(t) around the normal value ξ*(t).

(8.3)

We assume however that the evolution of the deterministic component ξ*(t) of the reaction extent ξ(t) is generally nonlinear and is given by (8.2). After some calculations we come to a time-inhomogeneous Fokker-Planck equation for the normal probability density P*(ξ,t):

∂P*(ξ,t) ∂ ) {[λ(t)(ξ* - ξ) - r(ξ*,t)]P*(ξ,t)} + ∂t ∂ξ D(t) ∂2 P*(ξ,t) (8.4) V ∂ξ2

( )

i

+ 2 2 (νi ) - (νi )

∑i 2[a (t )-∞) + f R(t)]

D(t) ) F(t) + (ξ*(t) - R(t))

i 0

(8.6)

i

is a time-dependent probability diffusion coefficient. For this model with nonlinear deterministic kinetics and linearized fluctuations the Hamiltonian function H(ξ,p) is a square function of the generalized reaction affinity p:

H(ξ,p) ) [r(ξ*(t),t) + λ(t)(ξ - ξ*(t))]p - D(t)p2

(8.7)

and therefore the normal chemical action J*(ξ,t) is the solution of the Hamilton-Jacobi equation

∂ ∂ J*(ξ,t) + [r(ξ*(t),t) + λ(t)(ξ - ξ*(t))] J*(ξ,t) ∂t ∂ξ 2 ∂ D(t) J*(ξ,t) ) 0 (8.8) ∂ξ

(

(

)

)

Both the Fokker-Planck equation (8.4) and the HamiltonJacobi equation (8.8) can be solved analytically. In both cases the equations can be integrated along the characteristics by using a nonlinear transformation which allows us to get rid of the deterministic component of the motion; the detailed computations are not presented here. The results are

(8.1)

dξ*(t) ) r(ξ*(t),t) with dt r(ξ(t),t) ) F(t)[Γ+(ξ(t),t) - Γ-(ξ(t),t)] (8.2)

+ with fi ) νi - νi (8.5)

is a time-dependent regression rate of fluctuations and

P*(ξ,t) )

where the forward and backward rate coefficients k( are timedependent. We assume that the normal solution ξ*(t) of the deterministic kinetic equation,

∆ξ(t) ) ξ(t) - ξ*(t)

(fi)2

{

x

Vγ(t) Vγ(t) exp [ξ - ξ*(t)]2 4π 4

}

γ(t) [ξ - ξ*(t)]2 4

J*(ξ,t) ) -

(8.9) (8.10)

where

γ(t) ) {exp[2

∫-∞t λ(t′) dt′]}/{∫-∞t D(t′′) × t′′ exp[2∫-∞λ(t′) dt′] dt′′}

(8.11)

is a characteristic volume which depends on the time evolution of the rate coefficients. For illustrating these general results we consider the particular case where the chemical reaction (8.1) is a two-species bimolecular reversible reaction:

2A1 a 2A2

(8.12)

In this case the sum of the numbers N1,N2 of molecules of types A1 and A2 is constant:

N1(t) + N2(t) ) N1(t0) + N2(t0) ) NΣ ) constant

(8.13)

The theory developed in this paper makes the prediction that for systems characterized by the same total number N∑ of molecules all solutions of the kinetic equation (8.2) tend toward the same normal form ξ*(t) which is independent of the different possible initial conditions ξ(t0), ξ′(t0), .... We have checked this theoretical prediction of the theory through numerical integration of the deterministic kinetic equation (8.2) applied to the particular case of the bimolecular reaction (8.12) for


Vlad and Ross

Figure 1. Graphical representation of the relative values of the reaction extent for the bimolecular process (8.12). The curves shown correspond to a sinusoidal variation of the temperature and to different initial conditions for the reaction extent, ξ(0) ) 50, 110, 220, 260, and 284, respectively.

Figure 2. Graphical representation of the normal solution ξ*(t) of the deterministic kinetic equation (8.1) applied for the bimolecular reversible process (8.12). All transient solutions of (8.1), including the ones represented in Figure 1, tend toward this normal form.

constant N∑ and different initial values of the reaction extent ξ(t0), ξ′(t0), .... For circumventing the emergence of an equilibrium solution in the long run we have chosen a periodic variation in time of the equilibrium constant generated by the sinusoidal periodic variation of the temperature. For simplicity we have assumed that the reaction enthalpy change of (8.12) is temperature-independent:

∆Hr ) independent of T

(8.14)

and thus the time variation of the equilibrium constant of (8.12) is simply given by 0 K1(t) ) k+ 1 /k1 ) K1 exp[-∆Hr/[kBT(t)]]

(8.15)

where K01 is the value of the equilibrium constant in the limit of an infinite temperature and kB is Boltzmann’s constant. We also assume that the time scale for the heat exchange process between the system and its environment is at least an order of magnitude smaller than the characteristic chemical time scale, and thus the temperature of the system can be externally controlled. Figure 1 displays a number of kinetic curves corresponding to different initial values of the reaction extent. We notice that all transient solutions of the kinetic equation converge toward the same normal ξ*(t) solution, which depends only on the total number N∑ of molecules and on the time variation of the rate coefficients, but it is independent of the initial conditions ξ(t0), ξ′(t0), .... The time dependence of the normal solution ξ*(t) is shown in Figure 2. Figure 3 displays the dependence of the normal chemical action, J*(ξ,t), on the reaction extent as well as on time. We notice that J*(ξ,t) varies periodically in time. As expected, J*(ξ,t) is maximum for the case when the reaction extent is equal to the normal solution, ξ ) ξ*(t) displayed in Figure 2. It is interesting to compare the permanent normal regime generated by the periodic variation of temperature with the periodic regime resulting from the periodic variation of the influx concentrations in a continuous stirred tank reactor. We notice the existence of an important difference between these two systems. For the cyclical variation of concentration the driving force toward the normal regime is additive and the corresponding coupling is linear, whereas for the temperature variation the coupling is multiplicative and nonlinear. This nonlinear coupling determines the complicated shape of the normal solution displayed in Figure 2, which shows the existence

Figure 3. Graphical representation of the normal chemical action, J*(ξ,t), for the bimolecular reversible reaction (8.12) as a function of the reaction extent and of time. The corresponding surface displays periodicity with respect to the time variation, and it has the maximum value zero for the deterministic normal solution ξ ) ξ*(t), which varies in time, as shown in Figure 2.

of two local maxima per period, whereas the excitation function, which is assumed to be sinusoidal, has only one. A referee came up with an interesting example of a reversible chemical reaction that apparently contradicts the results concerning the eventual emergence of the normal solutions of deterministic kinetic equations for large times. The analysis of this second model is useful for understanding the mechanism of generating the normal solutions. Let us consider a first-order reversible reaction:

A1 a A2

(8.16)

and assume that the total concentration of A1 and A2 is equal to

Ctotal(t) ) C1(t) + C2(t) ) constant

(8.17)

The kinetic equation of the process is

dC1(t)/dt ) -C1(t)ktotal(t) + k1 (t)Ctotal

(8.18)

and has the following solution:

∫t tktotal(t′) dt′] + t t Ctotal∫t exp[-∫t′′ktotal(t′) dt′]k1 (t′′) dt′′

C1(t) ) C1(t0) exp[-

0

0

(8.19)



where ktotal(t) ) k+ 1 (t) + k1 (t)

(8.20)

is a total reaction rate. Equation 8.19 shows clearly how the normal solution may emerge. In a chemical system obeying the mass-action law the individual rate coefficients are positive or equal to zero, but never negative; only the difference between two rate coefficients may be negative. Moreover the approach developed in section 4 is based on the assumptions used in the derivation of Shear’s theorem; that is, in the absence of an external constraint the system evolves toward a state of chemical equilibrium for which all chemicals are present in finite concentrations different from zero. This is possible only if all rate coefficients are different from zero. In particular, for the chemical reaction (8.17) the individual rates and the total rate ktotal(t) ) k+ 1 (t) + k1 (t) are always positive, and as a result the factor in (8.19) depending on the initial condition tends toward zero for large times:

∫-∞t exp[-∫t′′t ktotal(t′) dt′]k-1 (t′′) dt′′

lim C1(t) ) Ctotal

t0f-∞

(8.21)

As pointed out by the above-mentioned referee, in the pathological situation when the individual rate constants have a sinusoidal shape, ( k( 1 (t) ) k1 * sin(ωt)

(8.22)

in (8.19) the term depending on the initial condition C1(t0) does not vanish and a normal solution does not exist. However, (8.22) is never fulfilled in a real system because an individual rate coefficient can never be negative. A physically realistic periodic dependence of time of the rate coefficients would be by a modified version of (22): ( ( ( ( k( 1 (t) ) k1 * sin(ωt) + κ1 *, with κ1 * g k1 * > 0

(8.23)

If (8.22) holds, the positivity of the constants κ( 1 * leads to the vanishing of the transient term in (8.19) in the limit t0 f -∞ and the system evolves toward a normal regime independent of the initial condition C1(t0). 9. Suggestion for an Experiment We have shown that under suitable circumstances the solutions of a system of kinetic equations with time-dependent rate coefficients may evolve toward a normal form which is independent of the initial state of the system. In this section we suggest an experiment for testing this theoretical prediction. In projecting such an experiment we must make sure that all constraints considered in the theory are fulfilled. A candidate for such an experiment is a closed but nonisolated batch reactor. The reactor is enclosed in a thermostat and its size must be much smaller than that of the thermostat. The reactions should not be too fast, so that the characteristic time scale of the heat transfer between the reactor and the thermostat is much smaller than the characteristic chemical time scale. Under these circumstances the time variation of the rate coefficient can be controlled by varying the temperature. The temperature variation is not the only way in which the time variation of the rate coefficients can be externally controlled. For instance for certain reactions the values of the rate coefficients depend on the variation of an external electromagnetic field. An experimental device based on such an effect may be easier to control than a chemical reactor enclosed in a thermostat.

In addition to the above-mentioned general constraints we must also take into account the stoichiometric constraints of the theory. Our theory is based on the description of the chemical process in terms of the reaction extents attached to the different couples of forward-backward reactions occurring in the system. We also asume that the initial state of the system is described in terms of the initial values of these reaction extents. The advantages of using the reaction extents as chemical variables is that they lead to a symmetric structure of the equations, which simplifies the theoretical description of the system. A disadvantage of using the reaction extents is that they fully determine the composition of the system only if the initial composition vector is known. Because of this, if the initial composition vector is unknown, then it may happen that two systems characterized by exactly the same reaction extents may have different compositions. It has been shown that two chemical systems characterized by the same reaction extents have the same composition provided that their reaction invariants are the same.36 By reaction invariants we mean the conservation relations generated by the stoichiometric constraints. For example in the particular case of the bimolecular reversible reaction (8.12) there is only one reaction invariant, given by (8.13), which expresses the conservation of the total number N∑ of molecules. In general, for many reactions, there are many reaction invariants:

∑u ujNj ) N∑(u) ) independent of time,

u ) 1, 2, ...

(9.1)

where uj are dimensionless numbers and N(u) ∑ are constant. These reaction invariants are first integrals of motion of the deterministic kinetic equations (4.25) in the same way as the energy or the total impulse are first integrals of the dynamic equations in classical mechanics. In the literature there are algebraic methods for determining the reaction invariants and the coefficients uj in terms of the stoichiometric coefficients of the elementary reactions (4.1).36 We emphasize that the specification of the chemical invariants is not equivalent to the specification of the initial conditions. Different initial conditions may correspond to the same reaction invariants. The systems characterized by the same initial condition tend toward the same normal reaction regime only if the values of the reaction invariants are identical for all systems. In the case of periodic variation either of temperature or of an external electromagnetic field the experiment should proceed in the following way. The same chemical process should be carried out many times with exactly the same external conditions. For different experiments the initial composition vectors should be different, but they all should correspond exactly to the same values of the reaction invariants (9.1). Under these circumstances all kinetic curves expressed in terms of the reaction extents should tend toward a normal form depending only on the reaction invariants N(u) ∑ and on the time variation of state parameters of the system, but independent of the initial reaction extents. For carrying out such an experiment it is not necessary to have full knowledge of the kinetics of the system. It is however necessary to have a detailed knowledge of the stoichiometry so that the experimenter can make sure that the reaction invariants have exactly the same values for the different tests. 10. Conclusions In this paper we have investigated conditions for which a chemical system with time-dependent rate coefficients may forget its past. We started our analysis by considering a

8770 J. Phys. Chem. B, Vol. 101, No. 43, 1997 stochastic Markovian master equation with time-dependent transition rates. On the basis of a statistical ensemble approach we have constructed a Lyapunov (H function) for this master equation and showed that its transient solutions tend toward a unique normal probability distribution which, even though generally time-dependent, is independent of the initial conditions. The physical meaning of the suggested H function has been clarified by using the theory of information gain. We emphasize that these results are not limited to chemical systems; they can be applied to any physical or chemical process described in terms of a Markovian master equation with timedependent rate coefficients. The general approach has been applied to a particular case of closed ideal chemical systems described by the stochastic version of the mass-action law. We have shown that for these chemical systems, in addition to the stochastic H theorem for the master equation, there is a second H theorem for the deterministic kinetic equations. These two H theorems have been used as a basis for developing a thermodynamic formalism for time-dependent processes. An experiment for checking the existence of the normal solutions of the deterministic kinetic equations has been suggested based on the variation of certain externally controlled parameters of the system such as the temperature or the intensities of an electromagnetic field. The structure of the theory developed in this article is hierarchical. The classical master equation provides a Lyapunov function for the deterministic kinetic equations. At a higher level the ensemble master equation describing the ensemble statistics provides a Lyapunov function for the classical master equation. This hierarchical structure of the theory leads to two types of thermodynamic descriptions, similar to the Boltzmann (system) and Gibbs (ensemble) descriptions in equilibrium statistical mechanics. Despite this formal analogy the detailed structure of these two levels of description is different. One particularity of the deterministic kinetic description is related to the stoichiometric constraints. At a deterministic level two chemical systems tend toward the same normal regime only if their reaction invariants are exactly the same. With respect to the reaction invariants a deterministic chemical system displays infinite memory in the sense that during a process the invariants remain the same. This type of memory effect has been recently investigated in a study concerning the frequency response of deterministic chemical systems.37 Some types of physical and chemical systems can evolve toward a time-dependent normal regime even though their rate parameters are time-independent. For example, most limit cycles investigated in chemical kinetics are generated by processes with time-independent rate coefficients. The results derived in this article are valid for any time-dependent normal regime, whether the corresponding rate coefficients are timedependent or not. Until now the research has been focused on the case where the time variation of the rate coefficients is generated by the external variation of certain constraints. For many chemical and physical processes, however, the time variation of the rate coefficients is generated by the underlying dynamics rather than by the controlled variation of a constraint. This is the case of chemical reactions and relaxation processes occurring in disordered systems such as the combustion of certain amorphous explosives, the relaxation in glasses, or certain recombination processes of active intermediates in radiochemistry. The analysis developed in this article can be extended to these processes, and work on this problem is in progress. An important generalization of the approach considered here corresponds to the case when the variation of the external

Vlad and Ross constraints is random. In this case, due to the multiplicative coupling with a random excitation, a pure deterministic description of the system is not possible anymore. The coupling between the internal and external fluctuations gives rise to stochastic intermittence which can be described in terms of a double stochastic fluctuation-dissipation relation. We have reported some preliminary results in the literature,38 and further work is planned. Acknowledgment. This research has been supported in part by the Air Force Office of Scientific Research. We thank the referees for their constructive criticism. Appendix A Equations 2.9 and 2.10 can be derived in two steps. The derivation of (2.9a) and (2.10a) is based on the well-known inequality

x ln(x) + 1 - x > 0 for x > 0, x * 1 x ln(x) + 1 - x ) 0 for x ) 1 (A.1) By summing in (2.1) over all possible composition vectors, it follows that, if for t ) t0 we have

PN(t)t0) ) 1 ∑ N

(A.2)

then

PN(t) ) ∑PN(t)t0) ) 1 ∑ N N

for any t g t0

(A.3)

That is, all solutions of the master equation (2.1) preserve the conservation of the normalization of the probability distribution. By using (A.3), we can rewrite (2.8) in the form

hu,u′ )

(u) (u′) (u) (u′) P(u′) ∑ N (t){[PN (t)/PN (t)]ln[PN (t)/PN (t)] + N (u′) 1 - P(u) N (t)/PN (t)} (A.4)

from which, by using (A.1) we get (2.9.a) and (2.10a). The derivation of (2.9b) and (2.10b) is a little bit more complicated. We evaluate the time derivative of the H function defined by (2.8), resulting in

d

( ) ∑ [ ( )

hu,u′ )

dt

d

(u) (u′) P(u) ∑ N (t) ln[PN (t)/PN (t)] + N dt

P(u) N (t)

N

1

d

P(u) N (t)

dt

P(u) N (t) -

1

( )]

P(u′) N (t)

d

P(u′) N (t)

dt

(A.5)

In (A.5) the time derivatives of the state probabilities can be computed from the master equation (2.1) applied for PN(t) ) (u′) P(u) N (t) and PN(t) ) PN (t), respectively. After some algebraic transformations involving the grouping of terms in two double sums followed by the interchange of the summation labels, (A.5) becomes

d

hu,u′ )

dt

∑ ∑ N N 2

1

WN2N1PN(u′)2

{ [ ] PN(u)1

PN(u)1 PN(u′)2

PN(u′)2

PN(u)1

}

ln + PN(u′)2 PN(u′)1 PN(u)2 PN(u)2 PN(u′)1

(A.6)



Now we use an algebraic inequality similar to (A.1):

system (B.5) in terms of G(t|t0). After lengthy calculations we obtain

a ln(b/a) + a - b < 0, for a, b > 0, a * b a ln(b/a) + a - b ) 0, for a, b > 0, a ) b (A.7)

L{Y;t} )

GN N (t|t0)YN }n ∑∏{∑ N

0 N0 ∀nN j

0 j i

For solving the ensemble master equation (3.5), we introduce the generating function:

L{Y;t} )

...{[YN ]n ∑ ∏ ∏ ∀n n n

N1

1

N1

where YN1, YN2, ... are complex variables attached to the different possible composition vectors N1, N2, .... By using a standard procedure in the theory of stochastic processes,26 we express the master equation (3.5) in terms of the characteristic functional L{Y;t}. We obtain a first-order partial differential equation in L{Y;t}:

∂

L{Y;t} )

∂t

∑ ∑ WN N (t)(YN N N *N 1 2

1

The state probability P(n;t) can be computed by evaluating the derivatives of the characteristic function L{Y;t}. From the definition (B.1) of the characteristic function L{Y;t} we notice that

2

1

∂ Y ) L{Y;t} (B.2) N 2 1 ∂YN1

L{Y;t)t0} ) L0{Y0)Y}

From (B.8) and (B.9) we obtain (3.8). Appendix C To give a path integral representation of the solutions of the chemical master equation (4.19), we introduce an infinitesimal propagator:

g(ξ,tfξ′,t+∆t) dξ′, with

(B.3) with the initial condition

...{[YN0 ]n ∑ ∏ ∏ ∀n n n 1

[YN0 2]nN2 ...}P0(n)

N1

(B.4)

N2

The characteristic system of the partial differential equation (B.2) is given by

YN1 )

dt

∆tf0 (C.1)

g(ξ,tfξ′,t) ) δ(ξ-ξ′)

N1

d

∫g(ξ,tfξ′,t+∆t) dξ′ ) 1,

This infinitesimal propagator is the solution of the equation

where

-

)

∂ V-1 g(ξ,tfξ′,t+∆t) ) -H ˆ [t,ξ,V-1∇ξ] g(ξ,tfξ′,t+∆t) ∂t (C.2)

with the initial condition

L0{Y0} )

(

1 dn L{Y;t}|∀YN )0 u nN1!nN2!... d(Y )nN1 d(Y )nN2... N1 N2 (B.9)

N2

∀|YN1|, |YN2|, ... e 1 (B.1)

P0(n0) )

L0{Y0)YG+(t|t0)} (B.8)

P(n;t) )

[YN2]nN2 ...}P(n,t),

i

i

By combining (A.6) and (A.7) we get (2.9b) and (2.10b).

Appendix B

0 N0 j

∑ [WN N (t)YN N *N 1 2

2

2

In the limit ∆t f 0 the solution of (C.1) and (C.2) can be represented as an inverse Fourier transform

∫-∞+∞...∫-∞+∞db × exp{-∆tV{H[ξ,ib,t] - ∑ibj(ξj - ξ′j)/∆t} + j

g(ξ,tfξ′,t+∆t) ) [V/2π]R

o(∆t)2}, ∆t f 0 (C.4)

- WN2N1YN1(t)] )

1

AN N (t)YN ∑ N 1 2

2

(B.5)

21

∫t tA+(t′) dt′)

(B.6)

0

By comparing (2.5) and (B.6), we notice that the Green function G(t|t0) can be expressed in terms of the Green function G(t|t0) of the master equation (2.1). We have

G(t|t0) ) {G+(t|t0)}-1

where

H[ξ,ib,t] )

where the matrix elements AN1N2(t) are given by (2.4). The Green function attached to (B.5) is

ˆ exp(G(t|t0) ) T

(C.3)

(B.7)

That is, G(t|t0) is the inverse of the transpose of G(t|t0), G+(t|t0). For evaluating the characteristic function L{Y;t} we integrate (B.2) along the characteristics. We consider the initial condition (B.3) and express the solution of the characteristic

∑j Fjj(t){[1 - exp(-ibj)]Γ+j (ξ,t) + [1 - exp(+ibj)]Γj (ξ,t)} (C.5)

is a Hamiltonian function corresponding to the operator (4.19b) and (b1, ..., bR) ) b are the Fourier variables conjugate to the reaction extents ξ1, ..., ξR. Equation C.4 is based on the assumption that for ∆t f 0 the reaction extents ξ1, ..., ξR and the differential operators ∂/∂ξ1, ..., ∂/∂ξR are approximately commutative. This assumption is commonly used in the literature in other physical and chemical contexts, and it is justified by the fact that in the Hamiltonian from (C.2) the derivatives ∂/∂ξ1, ..., ∂/∂ξR operate only to the left. By using this property, it is easy to check that the noncommutativity is important only for higher orders in ∆t. We take into account that our description of a chemical process is Markovian and consider a succession of a large number of very short time intervals. The propagator (Green function) G(ξ,t|ξ′,t′) of (4.19a)


Vlad and Ross

can be expressed as a product of infinitesimal propagators of the type (C.4).

∫ ∫ ∏{dξ(R) db(R)[V/2π]R} ×

G(ξ,t|ξ′,t′) ) lim ...

possible functions ξ(τ). Embedding the pre-exponential coefficient into the integration measure D[ξ(τ)], we obtain

G(ξ,t|ξ′,t′) )

nf∞

R

∫∫ D[ξ(τ)] exp{V∫t′tL[ξ(τ),dξ(τ)/dτ,τ]dτ}

n-1

(R) (R-1) )/∆t}} ) ∑ {H[ξ(R),ib(R),tR] - ∑j ib(R) j (ξj - ξj

exp{-∆tV

R)1

∫∫ Dt′t[ξ(τ),b(τ)] exp{-∆tV∫t′ {H[ξ(τ),ib(τ),τ] t

(C.12) where the chemical Lagrangian L[ξ(τ),dξ(τ)/dτ,τ] is given by a Legendre transformation with changed sign of the Hamiltonian H[ξ(τ),p(τ),τ]:

R

ibj(τ) dξj(τ)/dτ}dτ} ∑ j)1

(C.6)

L[ξ(τ),dξ(τ)/dτ,τ] ) -H[ξ(τ),p(τ),τ] + p(τ)‚dξ(τ)/dτ (C.13)

where Dt′t[ξ(τ),b(τ)] is a suitable integration measure over the space of functions ξ(τ), b(τ) and

tR ) t′ + R∆t; ∆t ) (t - t′)/n, ξ(n) ) ξ

R

∫t′t{H[ξ(τ),ib(τ),τ] - ∑ibj(τ) dξj(τ)/dτ}dτ} )

exp{-∆tV

j)1

∑R V∆τR{H[ξ(τR),ib(τR),τR] ibj(τR) dξj(τR)/dτR}} ∑ j)1

We notice that (C.12) is equivalent to the relationship (4.21) for the propagator G(ξ,t|ξ′,t′) used in section 4. The superposition relation (4.20), which expresses the general solution of the master equation (4.19a), is a consequence of the linearity of (4.19a) in the state probability. For evaluating the expression of the chemical Lagrangian L[ξ(τ),dξ(τ)/dτ,τ] we insert (C.5) for the Hamiltonian in (C.14) and compute the gradient with respect to p(τ). This operation leads to nonlinear R scalar eqations for the components of the vector p(τ). By solving these equations and inserting the resulting expressions into (C.13), we come to (4.22a) and (4.22b).

We introduce the dimensionless reaction rates

(C.8) zj ) rj/(2Dj)

and introduce the integration variables

(C.9)

(C.10)

j ) [(dξj/dt) - rj]/(2Dj) where rj and Dj are the net reaction rates and the probability diffusion coefficients defined by (7.35) and (7.36). Equation 4.22 for the chemical Lagrangian can be rewritten as

By applying the method of steepest descent, we come to

L[ξ(t′),dξ(t′)/dt′,t′] )

∫ ∏ R

{ db h (τR) exp{-V∆τR{H[ξ(τR),ib(τR)τR] -

∑j 2Djσ(zj,j)

(D.2)

where

R

∑ j)1

(D.1)

and the relative fluctuation rates

where pj(τR) are determined by the condition

∇p(τR)H[ξ(τR),p(τR),τR] ) dξ(τR)/dτR

(C.14)

Appendix D

R

ibhj(τR) ) ibj(τR) - pj(τR)

∇p(τ)H[ξ(τ),p(τ),τ] ) dξ(τ)/dτ

(C.7)

Next we introduce an assumption physically equivalent to (4.21). Following Kubo et al., in the thermodynamic limit we can evaluate the integrals with respect to b(R) in (C.6) by applying the method of steepest descent. We assume that the fluctuation paths in the (ξ,b) space are sufficiently smooth so that a good approximation in the limit n f ∞ in (C.6) corresponds to time intervals ∆τR which are small but large enough so that as V f ∞ the products V∆τR are very large. We write the exponential in (C.6) as

exp{-

with

ibj(τR) dξj(τR)/dτR}}} ≈

{exp{-V∆τR{H[ξ(τR),p(τR),τR] ∏ R R

∑ j)1

pj(τR) dξj(τR)/dτR}}} (C.11)

The pre-exponential proportionality coefficient in (C.11) can be expressed as the infinite-dimension limit of a Gaussian integral, which can be evaluated analytically. This proportionality coefficient can be included in the integration measure over the space of functions ξ(τ); its concrete form is of little importance, and to save space we do not give it here. By using (C.11) in the thermodynamic limit, we can express the propagator G(ξ,t|ξ′,t′) as a simplified path integral over all

σ(z,) ) -(z + ) ln

[x

]

1 + 2z + 2 + z + -1+ 1+z

x1 + 2z + 2

(D.3)

is a universal function which has the same structure for all reactions. In Appendix B of ref 23b it has been shown that for any real value of z the function σ(z,) has a unique maximum for ) 0 for which σ(z,) is equal to zero and that otherwise this function is always negative:

σ(z,)0) ) 0; σ(z,*0)

Deterministic and Stochastic Asymptotic Behavior of Nonequilibrium

Recommend Documents