Deterministic and stochastic theory of sustained oscillations in autocatalytic reaction systems

Physica 93A (1978) 114-137 (~) North-Holland Publishing Company

D E T E R M I N I S T I C AND S T O C H A S T I C THEORY OF SUSTAINED O S C I L L A T I O N S IN A U T O C A T A L Y T I C R E A C T I O N SYSTEMS

R. F E I S T E L and W. E B E L I N G

Sektion Physik, Wilhelm-Pieck-Universitiit Rostock, DDR

Received 18 August 1977 Revised 7 February 1978

A special type of reaction in a uniform s y s t e m with two reactants is considered which is a generalization of the first Lotka scheme. Using the P o i n c a r ~ - B e n d i x s o n t h e o r e m the analytical conditions for the existence of limit cycles are derived and some e x a m p l e s are treated numerically in the deterministic picture. The onset of limit cycles is considered as a second order phase transition. T h e master equations are formulated and a general analysis of limit cycle reactions in the stochastic picture is given. The fluctuations of phase and amplitude and the correlation f u n c t i o n s are discussed. Finally, Monte-Carlo solutions of the master equation are presented. The relations between the deterministic and the stochastic picture are discussed.

I. Introduction

The application of the qualitative methods of nonlinear dynamics ~'2) to nonlinear chemical reaction kinetics is one of the most interesting new fields in the theory of irreversible processes. First L o t k a and Volterra derived nonlinear oscillations for a special model of chemical and ecological systems. Since Prigogine introduced the concept of dissipative structures, nonlinear chemical oscillations have been at the center of i n t e r e s t S ) . Several theoretical models (Brusselator, Oregonator, etc.) were introduced for the qualitative interpretation of the spontaneous oscillations o b s e r v e d experimentally, e.g. for the f a m o u s Zhabotinsky-reaction as well as for biochemical oscillationsS'7). Since the real reactions which show oscillations are of e n o r m o u s complexity the theory is restricted up to now to simplifying models which represent only qualitatively the o b s e r v e d properties. H e r e a special generalization of the first L o t k a scheme is investigated which shows sustained oscillation for special values of the parameters. Three questions are at the center of our interest: (1) the kinetic transition from the stationary regime to the oscillatory regime; 114

SUSTAINED OSCILLATIONS IN REACTIVE SYSTEMS

115

(2) the analogy to phase transitions of second order; and (3) the stochastic description of this transition. In recent years many investigations have been devoted to the analogy of the kinetic transitions between different types of stationary states to thermodynamic phase transitionsS). Here the transition from stationary states to sustained oscillations is studied if the system parameters cross critical values. This investigation continues an earlier work which was devoted to the study of the generalized L o t k a - V o l t e r r a scheme in the deterministic 9) and the stochastic '°) picture.

2. Deterministic theory

2.1. Stability analysis We consider a chemical reactor with the constant input of raw material Y per unit time A0 which reacts autocatalytically to the intermediate X. The intermediate X reacts to a final product which is removed from the system

A0

Ai

,Y

<

'

X

~B(X) Here F ( X ) denotes the autocatalytic feedback. Assume that the equations of the chemical kinetics are given by dX dt = Y F ( X ) - B ( X ) - A,X, dY dt = A o - Y F ( X ) + B ( X ) ,

(1)

where Y F ( X ) describes the autocatalytic rate law of the formation reaction and B ( X ) the rate law of the decay reaction, X and Y being the concentrations of the intermediate and the raw material, respectively. The input rate A0 and the decay rate A, are assumed to be positive (A0 > 0, A, ~ 0). The rate functions F ( X ) and B ( X ) may be arbitrary continuous differentiable functions satisfying the conditions

F(X)>O

and

Ao+B(X)>~O,

F(0)/> 0

and

B(0) ~<0.

if X > 0 , (2)

116

R. FEISTEL AND W. EBELING

Introducing dimensionless time and concentrations by = A~t;

x = X A~ A---oo;

1 ~[

Ao\ f ( x ) = -~t t ~ x ~11,];

A~ Y = r --Aoo" b(x) = ~

1

B

(x Ao~,

(3)

A,]

w e g e t the f o l l o w i n g s y s t e m of t w o d i f f e r e n t i a l e q u a t i o n s : dx dr

--

--- ~t =

y/(x)

b(x) - x,

-

(4) dy= dr ~ = 1 - yf(x) + b(x). F o r f u r t h e r c a l c u l a t i o n s it is u s e f u l to define an a b b r e v i a t i o n g ( x ) =--1 + b ( x ) f(x)

if f ( x ) # O.

(5)

'

g ( x ) is t h e r a w m a t e r i a l b a l a n c e ; if y < g ( x ) , y is i n c r e a s i n g a n d v i c e v e r s a . So we get for f(x) # 0 another form of our differential equations: = [y - g ( x ) ] f ( x ) + 1 - x, = [g(x) - y]/(x).

(6)

It m a y b e o f i n t e r e s t to n o t e t h a t t h e e q u a t i o n s f o r t h e B r u s s e l a t o r a r e o b t a i n e d as a s p e c i a l c a s e o f eq. ( 4 ) i n t r o d u c i n g X = A x , Y = A y , f ( x ) = A 2 x 2, b ( x ) = B x - 1. T h e r e e x i s t s o n l y o n e e q u i l i b r i u m p o i n t o f the s y s t e m (4) or (6), r e s p e c t i v e l y , w h i c h is g i v e n b y x0 = 1;

y0 = g(1).

(7)

T h e e i g e n v a l u e s o f t h e J a c o b i a n o f (6) at this p o i n t are p = - ½(f + g ' f + 1) + X/~(f + g ' f + 1)2 _ f ,

(8)

w i t h f =- f ( x = 1) a n d g' =- d g ( x = 1)/dx. T h e s t e a d y s t a t e is s t a b l e if

- g ' < l + llf

(9)

a n d is o f o s c i l l a t o r y t y p e ( f o c u s ) if 1

2

holds. T h e t y p e o f e q u i l i b r i u m p o i n t (7) is r e p r e s e n t e d g r a p h i c a l l y in fig. 1. T h e s t a b i l i t y c o n d i t i o n m e a n s g e o m e t r i c a l l y t h a t t h e c u r v e ~ = 0 in the p h a s e p l a n e m u s t h a v e a s l o p e g r e a t e r t h a n - 1 at the p o i n t (7). In t h e limiting c a s e o f

g[

SUSTAINED OSCILLATIONS IN REACTIVE SYSTEMS

ll7

18

-g'= (1,1~')

16 i I l I

12.

,

uns~le node ' = 1 + 1If

10 ¸

/,/I I• //~///,,y///s/ // tocus /

6"

//~,,

/

/

stone focus

-g'=(1 -1~1 ~ //

/

/

oo'toNe

Fig. 1. Stability diagram of the equilibrium point (Xo, Yo). The values of g' and f at the equilibrium

point determine its stability properties in the presented way. For the models proposed by Lotka (f(x) = a3x; b(x)= 0) by Sel'kov (f(x)= a4x2; b ( x ) = 0) and for the model used in this paper (f(x) = a2 + 0.1.1"2; b(x) = 0) the set of points (g', f) for different values of the parameters a3, a4 or a2, respectively, form curves in the g'-f plane which are shown too. vanishing b a c k w a r d r e a c t i o n b ( 2 ) = 0 any a u t o c a t a l y t i c f e e d b a c k with the property f'>f(f+

l)

(ll)

m a k e s the equilibrium point unstable. As is easily seen, this condition involves that f ( x ) is at least piecewise c o n v e x in the interval [0, 1]. H e n c e , a strong positive (nonlinear) f e e d b a c k is a n e c e s s a r y condition f o r instability w h i c h is a prerequisite for sustained rotations a r o u n d this point. In special cases there is a s e c o n d attractor in the phase plane, and its stability analysis is of interest for our later considerations. If the condition f(0) = 0 is satisfied, we are able to give a t i m e - d e p e n d e n t solution of the differential e q u a t i o n s (4). F o r the initial conditions x ( 0 ) = 0, y(0) arbitrary, the f u n c t i o n s xO') = O;

y ( r ) = y(O) + "r

(12)

118

R. FEISTEL AND W. EBELING

are time-dependent solutions. These solutions c o r r e s p o n d to the situation where the concentration of the product is zero for all times and the concentration of the raw material increases unlimited with time. The stability analysis of the solution (12) yields that this trajectory is stable if f'(O) = 0 and, if f'(0) > 0, it b e c o m e s unstable for times greater than r0 =

1 + b'(0)

f'(0)

y(0).

(13)

An unlimited increase of the stock of raw material Y in the system seems to be physically unrealistic, therefore we shall restrict ourselves here to systems with f e e d b a c k functions fulfilling f(0) > 0 or b(0) < 0 or f'(0) > 0; b(0) = 0.

(14)

Then the trajectory (12), which leads to infinity, is not a solution or is unstable for large times. We mention that Sel'kov treated a reaction with f ( x ) = a4x 2, i.e. f(0) = 0 and f'(0) = 01~:2). 2.2. S u ~ c i e n t conditions f o r the existence o f a limit cycle The investigation of the stability map shown in fig. 1 d e m o n s t r a t e s clearly that the generalized L o t k a scheme studied here exhibits only transitions of stability-instability in the focus regime, i.e. f r o m a stable to an unstable focus. A quite similar picture is o b s e r v e d for the P r i g o g i n e - L e f e v e r - N i c o l i s reaction3). Mathematically this transition c o r r e s p o n d s to a so-called H o p f bifurcationZ), i.e. the two complex conjugate roots of the problem cross the imaginary axis for the bifurcation values of the p a r a m e t e r s given by - g ' = 1 + (I/f). This type of transition is called hard-mode instability as distinct from transitions in the node regime, i.e. from a stable to an unstable node, which are denoted as soft-mode instabilitiesS't3). In m a n y cases the H o p f bifurcation is connected with the a p p e a r a n c e of limit cycles. In order to investigate this point for the given reaction scheme we construct a c o m p a c t region in the phase space including the point (x0, Y0) given by eq. (7), which cannot be left by the representing point. This construction is different for the cases f(0) > 0 or b (0) < 0 and f(0) = 0, b(0) = 0. If f(0) > 0, we use the quadrangle P~P2P3P4 (see fig. 2a) with P1 = (xl, 0);

P2 = (Xl, Y2);

P3 = (x3, Y2);

P4 = (x4, 0),

(15)

where xl=O;

x3>l;

y2 >

sup {g(x)};

xE[xl.x3]

yz>l;

x4=x3+Y2 .

(16)

S U S T A I N E D O S C I L L A T I O N S IN R E A C T I V E S Y S T E M S

119

Y~ I

P,

i

_~

gO).

a)

Y

p~ b)

g(x~ In (X~/X,) i

i I

Y~

tin(~/x,) 1 "~(0)_ f'(o)

o. ~o.,).

\

J

',

'~a

Xa

Fig. 2. Construction of a c o m p a c t region used to apply the P o i n c a r 6 - B e n d i x s o n theorem: (a) the case f(0) > 0; (b) the case f(0) = 0, f'(0) > 0.

The lines P~P2 and P4Pt cannot be crossed by trajectories because the concentrations are non-negative numbers. Along P2P3 we find according to eqs. (6) # < 0, and along P3P4the relation ~ + # = 1 - x < 0 holds. This proof can be repeated in the same way in the case b ( 0 ) < 0 with xt < 1, 0 < xt < - b ( x t ) , because of :~ > 0 at Xl. If f(0) = 0 the situation is more complicated. We restrict ourselves to the case given by condition (14), f ' ( 0 ) > 0. The above picture remains unchanged except for the straight line along the y-axis, because x = 0 is now a trajectory going to infinity. P~ leaves the origin, away from both axes, PIP2 becomes a parallel x~ = const to the ordinate with sufficient small x~, and an appropriate

120

R. FEISTEL AND W. EBELING

curve defined by xg

y51(x) = f rl ++bb(r) (r)dr;

x E [ x t , x~]

(17)

x

connects P~ with a point P5 = (x~, 0) on the x-axis (see fig. 2b). To p r o v e the possibility of such a construction we formulate the conditions for P~, P2, P5 in a more precise way. Using the auxiliary quantity 3' inf {x l + b ( x ) ~ 0 < 3 " ~< x~0., x T b~-~J ~< 1,

(18)

we choose xt to be so small that the inequalities 0
xt+b(xO<3"[l+b(xt)];

1;

-lnxt>

xt + b(xO

3"f(xO

(19)

are fulfilled. Then it is possible to find an x~ satisfying

fx~ + b(x0] xl e x p { g ( x 0 } > x s > x t e x p ] yf(x,) {,"

0 < x5 < 1;

(20)

H e n c e for 0
1

(21)

and for x 5

f l + b(r) dr Yl = r + b(r)

(22)

Xl

the chain of inequalities xl + b(xO

f(xO

< y In xt x--2~< Yl <~ln XS <

(23)

is true. As fig. 2b demonstrates, the conditions (21) and (23) really guarantee the possibility of our construction for e v e r y b(x) and f(x). The definitions of P2, P3, P4 are given a b o v e [eqs. (15), (16)], but note that x~ has changed its value. N o w we have to p r o v e again that the velocity field at the boundary points to the inside. Along the lines P2P3 and P3P4 the situation remains unchanged. Along the parallel PiP2 to the ordinate, which is located completely " a b o v e " the curve ,~ = 0, one has .~ > 0, i.e. the flow is directed into the interior. For the curve PsP1 given by eq. (17) we define first the normal vector n pointing inside and the velocity vector v

+ b(x)'~, n = (~+b(x)/"

~ = (.~).

(24)

SUSTAINED OSCILLATIONS IN REACTIVE SYSTEMS

121

The scalar product following from eq. (4) is always positive n T. i)

:

y/'(x)(l - x) > 0,

(25)

since along the curve (except Ps) we have x < 1, y > 0. At P5 the velocity field is tangential to the curve ysl (17), but the stronger curvature leads the trajectory to the inside

,, = Y51(Xs)+

[1 + b(xs)l/(x~) [ X s + b(xs)] 3 (1 - x 0 >

Y~'dP~.

(26)

So we can conclude that the stated boundaries are either crossed by the trajectories f r o m outside to inside, or are not crossed at all. The internal equilibrium point (7) can be excluded by a small surrounding ellipse defined by the eigenvectors of the Jacobian of the system at this point. If the steady state is unstable and both eigenvalues have positive real parts, all trajectories leave the steady state and cross the ellipse. According to the P o i n c a r 6 - B e n d i x s o n theorem inside the region described a b o v e at least one stable limit cycle exists. We can summarize our results in the following theorem: if the f e e d b a c k function satisfies the conditions (1)

1//:(1) < - g ' ( 1 ) -

(2)

f(0)~0

1

(27)

and or

b(0)<0

or

f'(0)~0; b(x)>tO;b(O)=O,

(28)

then the reaction system exhibits at least one stable limit cycle. These requirements mean that the conditions for an unstable steady state (9) are satisfied and that there are no stable trajectories going to infinity due to the inequalities (14). 2.3. Special reaction models As f o r m e r investigations have shown, the backward reaction is not very important for the behavior of the system. Therefore, in what follows we shall restrict ourselves to the case b ( x ) =- 0 and discuss several special choices for the f e e d b a c k function f ( x ) . The simplest case discussed by L o t k a is f ( x ) = a3x. In this case a stable final state 5) always exists (see fig. 1). The case f ( x ) = a4x 2 is more complicated. Sel'kov could show that for certain values a 4 < 1 limit cycles exist11). But in this case / ' ( 0 ) - - - / ( 0 ) = 0 and the theorem given a b o v e cannot be applied. In recent papers S a m o y l e n k o and Sel'kov investigated models more general than (1) for e n z y m e reactions which exhibit

122

R. FEISTEL AND W. EBELING

also limit cycleslZ). Furthermore, we mention that a reaction model similar to ours was recently investigated by Ibanez et a1.14), which also shows a H o p f bifurcation connected with a transition to a limit cycle regime. In previous papers 9'~°'~5) the authors studied the case f ( x ) = a3x + a4x 2.

(29)

The application of the theorem derived above yields the existence of a limit cycle in the range a3>0;

(a3+a4)Z
(30)

Since / ' ( 0 ) = a 3 ~ - 0 , the solution (12) becomes unstable for large y, and consequently the corresponding time T = l/a3 may be considered according to (13) as an estimate for the magnitude of the time period of the limit cycle. Another example we want to discuss in more detail in this paper is given by the feedback function f ( x ) = a2 + a4x 2.

(31)

According to our theorem limit cycles exist in the region (a4 - - a2) > (a4 + az) 2.

(32)

Choosing a2 as the control parameter and fixing a4 = 0.1 we find a focus regime in the interval 0.0122< a 2 < 1.6822. Using the condition (32) we derive for the onset of limit cycles the critical value a~=0.07082 ( a z > a ~ means stability). In order to demonstrate the transition from a stable focus regime to a limit cycle regime, in figs. 3-5 three numerical examples are given. Near to the bifurcation point (az = 0.07, fig. 4) a limit cycle with a small amplitude is formed after a relatively long time. We note that in the critical point itself one observes a limit cycle with the amplitude zero. ~5) Far from the bifurcation point (a2 = 0.01, fig. 5) we find after a very short transition time a limit cycle with large amplitude with the typical form of a relaxation oscillation. Each period of these relaxation oscillations has three typical phases corresponding to the three sides of the triangle in the phase space. (1) Convergence to the trajectory (12) connected with a slowly increase of the raw material Y by import from the external reservoir. (2) After reaching the line p = 0 (y ~ g(0) = lla2) a fast autocatalytic reaction starts which consumes the raw material Y. The concentration of X increases quickly. (3) After exhausting the stock of raw material the autocatalytic reaction stops, the intermediate X decays rapidly, and the cycle starts again.

SUSTAINED OSCILLATIONS IN REACTIVE SYSTEMS

x ,AY

123

_N v

~,oflo,o

H

~1.

/I " \

L/

I ~1

/

I

\J

J/

I ~'x

0~

a,

2'o

i0

do

K0

Fig. 3. Numerical solution: stable focus near the threshold (connected points: a corresponding

stochastic realization). 2.4. The onset o f sustained oscillations as a second order phase transition Phase transitions of second order are characterized by two essential properties16): (1) the transition is connected with a change of s y m m e t r y of the system; and (2) in the transition point (;t-point) the two phases coincide. The concept of phase transitions has been applied to the kinetic transitions in nonlinear irreversible processes by m a n y authorsS'17-2°). The kinetic transitions which correspond to t h e r m o d y n a m i c phase transitions of second order are generally characterized by the feature that the system goes out of the way of external constraints by passing to a more effective physical mechanism. Lasers or catalytic reactions can be considered as examples of this kind of behavior. The analogy between the onset of sustained oscillations and second order phase transitions was first investigated by Drosdziok 2°) and Tomita et aj.13). It will be shown here that this analogy allows a proper understanding of the qualitative properties of our reaction model. We start with the following observations. (1) The onset of sustained oscillations is connected with a change of s y m m e t r y . In the stationary regime the system is invariant, for t ~oo, with

124

R. F E I S T E L A N D W. E B E L 1 N G

1,0

0~

L.u

~

v v

v

v

t

Fig. 4. Numerical solution: unstable focus near the threshold (connected points: a corresponding stochastic realization).

X,

20. i/~

iiI

IP t

/ 11

iI

l i It i

I

10

i

u

iI

iiI

iI I

lI

II i

/

i

ii

/I

i

ii¢ I

ii

/

II

lI iI

I

t

/

i i

i

ii

/

iI

i~(

iI

i

/

I

I

I I

•

1'0

20

.

30

Z,O

50

w

, .

60

.

"70

80

Fig. 5. Numerical solution: unstable node far from threshold.

SUSTAINED OSCILLATIONS IN REACTIVE SYSTEMS

125

respect to the transition t ' = t + C, with an arbitrary value of C; in the limit cycle regime the invariance is restricted to the transformation t ' = t + n T (n = integer, T = period). (2) In the transition point one observes a limit cycle with zero amplitudeS5). T h e r e f o r e both phases coincide in the A-point. The different mechanisms are given by the noncatalytic reaction [ ( x ) - a2 and the autocatalytic of second order f ( x ) ~ a4x 2. For large a2 the simple m e c h a n i s m dominates; the corresponding linear reaction equations have a stable equilibrium point. If a2 is sufficiently small the autocatalytic mechanism will govern the behavior of the system; and for small a4 limit cycles occur. As investigations of competitive reactions have showng), the appearing concentration of raw material reached is a measure of the effectiveness of the mechanism. We are able to define a very similar quantity with certain extreme properties in our model if the backward reaction b ( x ) is neglected. It is well known f r o m the stochastic description that the ensemble average of all quantities converges to a stable stationary value for t ~oo also if single systems carry out u n d a m p e d oscillations. Therefore, we average the kinetic equations (7) o v e r a suitable ensemble of initial conditions (or over a period in the final state, respectively) (Yc) = ( y f ( x ) ) - (x);

(9) = 1 - (y/(x)).

(33)

Stationarity leads to (x) = 1

(34)

independently of the occurrence of either stable states or limit cycles. Dividing the second equation by y before averaging we get d--~(ln y) = ( l / y ) - or(x)),

(35)

and in the final state ( l / y ) = 0r(x)).

(36)

For convex functions f ( x ) we have or(x)) >~f ( ( x ) )

(37)

and consequently ( l / y ) = or(x)) I> f ( ( x ) ) = f(Xo) = l/y0.

(38)

(The equality sign is valid iff the final state is a steady one.) The mean reciprocal value of the raw material concentration is in the limit cycle regime larger than the value at the (unstable) steady state (fig. 6). As a consequence of this statement the curve is continuous but not differentiable at the A-point and d e m o n s t r a t e s the onset of another reaction mechanism. A similar result

126

R. FEISTEL AND W. EBELING

0,3

~

~

timtt cycle

9,2

(V->

~tobl~

steQdy state 0,1

2.2

2/.

2'6 l ' ' threshold

218

3,0

-in A:.

Fig. 6. The mean reciprocal of the raw material concentration depending on the control parameter A: (dashed line above: slochastic result). w a s o b s e r v e d by M c N e i l et al. ~8) b y the a n a l y s i s of a o n e - d i m e n s i o n a l s y s t e m . T h e final s t a t e is a l w a y s c h a r a c t e r i z e d by a m a x i m u m of ( l / y } if f(x) is c o n v e x in the i n t e r v a l Xm~n
3. Stochastic theory

3.1. Master equation and phase fluctuations F o r s y s t e m s with large p a r t i c l e n u m b e r s the d e t e r m i n i s t i c p i c t u r e as u s e d in t h e f o r m e r s e c t i o n s is s u i t a b l e to d e s c r i b e t h e i r b e h a v i o r . This p i c t u r e c o r r e s p o n d s to the t h e r m o d y n a m i c limit. A d e s c r i p t i o n , w h i c h t a k e s into a c c o u n t f l u c t u a t i o n s a n d the d i s c r e t e c h a n g e of p a r t i c l e n u m b e r s , is given b y the m a s t e r e q u a t i o n t r e a t m e n t . D e n o t i n g the v e c t o r o f p a r t i c l e n u m b e r s of the d i f f e r e n t s p e c i e s b y r we w r i t e P(r; t) f o r the p r o b a b i l i t y of finding o n e s y s t e m o f an e n s e m b l e at t i m e t in the p o i n t r o f the p h a s e s p a c e . T h e p a r t i c l e n u m b e r v e c t o r r c h a n g e s w i t h

SUSTAINED OSCILLATIONS IN REACTIVE SYSTEMS

127

the probability W~(r)dt within dt by an elementary reaction in the channel i from r to r + Ari. The conservation of the number of systems is expressed by the continuity equation

00t P(r;

t) = ,~. {W~(r-

Ari)P(r- Ari; t ) - W~(r)P(r;

t)},

(39)

where the sum is to be extended over all elementary reactions. This is the so-called master equation or Pauli equation of the stochastic reaction system. As shown by several authors the solution of the master equation converges under rather general conditions asymptotically to a unique stationary distribution 21-24)

P(r; t)-->P°(r),

if t--->~.

(40)

For a limit cycle system p0 has the form of a probability crater. Owing to the H-theorem the mean particle numbers {r) converge to a time-independent value. This seems to be in contradiction to the deterministic picture. Similar problems have been investigated for the case of bistable systems25'26). The apparent different behavior of (r) in both pictures is due to the effect of the fluctuations. In the case of limit cycle systems the phase fluctuations are the most relevant fluctuations. Stratonovich zT) and Tomita et ai. j3) have treated the phase fluctuations of stochastic limit cycles with continuous variables. In our case of chemical limit cycles the discrete character of the steps may be essential especially for limit cycles with small amplitudes. For the investigation of this case we shall develop here a new method which is based on a waiting time analysis. First, we calculate the waiting time At(r) of the phase point at r. This is a first-passage-time problem and can be treated following Stratonovich27). We solve the master equation (39) under the initial and boundary conditions

P(r'; O) = 8r~'

P(r';

t) = 0,

if r' ¢ r.

(41)

Using (41) we find by a simple integration of eq. (39) unnormalized distribution

P(r';

t)=exp(-

~. W,.(r)t}Sr~,

(42)

and, furthermore, the waiting time distribution density 0 P(At)=o--~[1-P(r;At)]

= .~i W / ( r ) e x p { - A t .~i W/(r)}

(43)

128

R. F E I S T E L

AND

W. EBELING

with the mean and variance --

__

]1/2

I

Now let us consider a sequence r~, r, . . . . . rA of points forming a closed pathway enclosing the unstable stationary point which the system is able to pass. The mean rotation period is T --

rj

,

(45)

where the average has to be taken over all possible closed pathways around the stationary point. The variance of T is

The first term is due to the variance of the waiting times (44) and the second one describes the variance caused by the different pathways. After k rotations the variance becomes

[A T(k)] 2-- k A T : .

(47)

Now we define the number of rotations up to the phase balance A T ( n ) ~

T by

(48)

n =- ( T I A T ) 2.

If the amplitude of the limit cycle is small the Wi(r) in eq. (44) may be replaced by W~(r0). Then, according to eqs. (44) and (45), there follows for the mean rotation period T = (A)At(r0) = l

Wi(ro),

(49)

where I = (A) denotes the mean number of elementary reaction steps along the cycle. For the number of rotations up to the phase balance we find in the same approximation l2 n = l + ( A 2)-12"

(50)

Neglecting amplitude fluctuations we find n ~ I. We shall see later that this simple approximation overestimates the number of rotations up to the phase balance. Finally, we note that there exists another method for the stochastic in-

SUSTAINED OSCILLATIONS IN REACTIVE SYSTEMS

129

vestigation of limit cycles which starts from the correlation functions ~o(~') = {(r(7) - (r)) o (r(0) - (r))).

(51)

From the matrix of correlation functions the mean rotation frequency is obtained by ~3) 03= d-r-r~°xy

,=0"

For clockwise rotation 03 < 0 holds. Further information about oJ may be derived from the peaks of the spectral function z2) S(~o) = ~1. , f ~('r)

e i'~r

dr.

(53)

0

3.2.

Stochastic description o[ the reaction model

Let N be the particle number of the species X in a volume V and M the particle number of Y, respectively. The vector of particle numbers is (54) There are four reaction channels

AoV;

1) ~ Y

with

Arl=(+~)

and

Wl(r)=

2) Y ~ X

with

Arz= (+_I)

and

WE(r) = M E ( N ) ; (55)

3) X ~ Y

with

Ar3=(+ll)and

W3(r)=B(N)v;

4) X ~

with

At4=

W4(r) =

and

AIN.

Here F and B are the functions defined in eqs. (1) and (2). Some authors propose modifying such functions by replacing the powers X" by factorial products ~,~,u, , / . This makes allowance for the decrease of the number of free particles at the formation of n-particle complexes. We agree with this conception, but note that it becomes meaningless for large N. With (55) the master equation (39) reads

AoVP,N M-

130

R. FEISTEL AND W. EBELING +

+ AI(N + I ) P ( N [a0v+

+ 1, M; t)

MF(N)+ B(N) v+ A,N]P(N, M't).

(56)

The average without approximation yields

(57) dt

V

and for the steady state (N) V

A0 AI"

(58)

This is an exact result in agreement with eq. (34). In this model the mean value of the concentration of X at t ~oo is always equal to the steady state value X0 [eq. (6)] of the deterministic model. Neglecting higher order correlations, eq. (57) reduces to eq. (i). N o w we restrict ourselves to the special case B = 0 without backward reaction and multiply the master equation (56) by the harmonic sum

H(M)

= ~ 1

(59)

i=li"

Summation over all N and M yields

d_ dt (H(M))= AoV(M---~>-(F(N)>,

(60)

and at the steady state we get V

1

F

N

in analogy to eq. (38). In the region of stability far from the threshold the fluctuations are small and the mean value V I <~-T-~> ~ ~0

(62)

coincides with the reciprocal of the stationary raw material concentration. Near the A-point the fluctuations grow and so does (V[(M + I)). The phase

SUSTAINED OSCILLATIONS IN REACTIVE SYSTEMS

131

transition b e c o m e s " s o f t e n e d " in the stochastic picture (fig. 6), and a sharp A-point cannot be observed. In our example A0 = 1, A, = 1, V = 100 and F ( X ) = 0, 1X2+ A2 was chosen, therefore eq. (61) yields 100

( N ( N - 1))+ A2,

(63)

i.e. the value of (V/(M + 1)) can be calculated from the second m o m e n t of N. (This equation and (58) can be used as a simple test of whether the steady state is reached in numerical calculations.) For these p a r a m e t e r s we find a phase stability after n = 3 T V = 300T rotations by using the estimate n ~ I.

3.3. Monte-Carlo calculations f o r the reaction model and discussion Much information about the solution of the master equation (56) can be derived by the realization of the corresponding stochastic process on a computer. We approximate the probability P ( r ; oo) by the relative f r e q u e n c y of reaching the state r during a time T, where T is sufficiently large. For this purpose we take a simple time average until the corrections to the mean values b e c o m e small. The time required is much smaller than the time till the initial phase correlation is d a m p e d out. Of course the distribution created by this procedure is not yet smooth, but it reflects the essential properties of the real distribution. With V = 100 (one hundred particles correspond to the concentration 1) and A0 = 1, Al = 1, A4--- 0.1 the cases A2 = 1.0

(stable, far from the threshold),

A2 = 0.08

(stable, near the threshold),

A2 = 0.07

(unstable, near the threshold),

A2 = 0.01

(unstable, far from the threshold)

have been realized in comparison to the deterministic picture [figs. 7-10a (resp. b)]. The realizations were calculated up to about 200000 elementary reaction steps. The time behavior of two realizations is represented in figs. 3 and 4. First, one can see that the system-size expansion method of van K a m p e n 22) and Tomita 13) gives a qualitative correct result far from the threshold (figs. 7 and 10). The fluctuations to the secular motion are in general symmetrical. In the limit cycle we get a "probability crater" (limit gully), its height is

132

R. F E I S T E L A N D W. E B E L I N G

X

I o

>-

o6

~X

~D

d

SUSTAINED OSCILLATIONS IN REACTIVE SYSTEMS

133 ,,..,

I

LI~I I

,I ii' !11!111!

Z

kX

I I~ii:!

I !li II Ill~11 '

,§

2~

al

i

134

R. FEISTEL AND W. EBELING

determined by the local reaction velocity and the amplitude fluctuations. The way along the crest is practically the deterministic limit cycle (fig. 10). In the stable region a peak occurs, i.e. a 8-like distribution at the steady point (fig. 7). In the neighborhood of the threshold we see significantly increasing fluctuations (figs. 8 and 9), but the crossing of the critical values is not connected with a qualitative change. The single realizations exhibit an oscillating behavior on both sides of the threshold (figs. 3 and 4). The occurrence of limit cycles is reflected in the stochastic picture by the occurrence of m a c r o s c o p i c fluctuations around the (unstable) stationary point. Corresponding to the amplitude of the limit cycle the second m o m e n t s increase strongly over a few orders of magnitude behind the threshold (fig. 11), and so does ( V [ ( M + 1)) according to eq. (63). For ( - I n A2)> 3 the second m o m e n t of N is determined nearly exclusively by the size of the limit cycle (compare with the dotted drawn value of the deterministic calculation, which one gets by averaging over one period). In fig. I1 the measure of the phase diffusion n = ( T / A T ) 2 according to eqs. (45)-(48), the expression T E~ W~(ro), the mean number I of reaction steps per rotation, and the part of n caused by the waiting time variance [eqs. (46) and (44)] are presented. The curves of l and T E i W~(ro) are virtually identical.

10b /

(N~) -(N> 2

/

105 /

10~ ~ T ~ - I

I0 z

( Z~-tz )

n -(T)

134 E

10 °

i 4

0

5

6

7

-knA z

threshold

Fig. 11. Increasing fluctuations above the thresl~old and the phase stability measure n (see text).

SUSTAINED OSCILLATIONS IN REACTIVE SYSTEMS

135

which justifies the approximation (49) for this example completely [l was counted, T calculated from eq. (45)]. The agreement of I and T2/(E At:) is good with respect to the order nf magnitude, hence also the approximation

is confirmed. Furthermore, we see that the amplitude fluctuations are not meaningless for a system with a small particle number as in our case (which is also evident from figs. 9 and 10); on the other hand, they are responsible for the value of n [calculated from eqs. (44)-(48)] almost exclusively. This is demonstrated by the significant difference between n and T2[( 'y, At2), because the latter, according to eq. (46), is the value of n for omitted amplitude fluctuations. The " p e r i p h e r y " l of the limit cycle does not change essentially over the control parameter, in contrast to the "size" ((N 2 ) - (N)2) l/2 of the limit cycle. Obviously both quantities are not proportional. The explanation is that the fluctuations near the threshold often are perpendicular to the secular motion, whereas they are mostly tangential far from the threshold. Consequently, the limit cycle in the latter case is nearly triangular, one of the reaction channels (55) is significantly more probable than the others. Since the absolute value of n is about 102 in our examples, we have to wait for about 100 rotations (equivalent to nl--~ 106 reaction steps) for the phase balance. In the same manner the asymptotic part of the correlation function dies away. Fig. 12 shows the N - N autocorrelation function at A2 = 0.05, calculated with an ensemble of 400 systems. The first part is strongly d a m p e d corresponding to a fast "forgetting" of the " f o r m " of the initial distribution. Then only the phase correlation continues the oscillations. Therefore limit cycle systems can be considered as systems with a long-range order with respect to time. Finally, we briefly discuss the results for the different approximations for the rotation period. In the deterministic case we got a period of T -~ 16.7. The average T = (Ei At(rj-))~ 17.5 is in a good agreement with the peak of the spectral function (fig. 13) at to ~ 0.36 corresponding to T---17.5. The mean frequency 03 was calculated from eq. (52) to be 03 = - 0 . 4 , corresponding to T ~ 15.7 and a clockwise rotation sense. We conclude that the Monte-Carlo treatment of the master equation is consistent with the deterministic picture as well as with the estimates derived from the master equation, and yields also some new insights into the stochastic mechanism of limit cycle reactions.

136

R. FEISTEL

AND

W. EBELING

P i+% 6

L

2

3

-2

-L

5,o

1.0

Figs. 12 and (Monte-Carlo

13. The autocorrelation solution).

function

q NN and its cosine

transform

S,*$ for a limit cycle

References 1) A.A. Andronov, A.A. Witt and SE. Chaikin, Theorie der Schwingungen. Teil I und II (Berlin. 1965. 1969). 2) H. Margenau and H.M. Murphy (Eds.). Mathematik fiir Physik und Chemie. Bd. II (Leipzig, 1%6).

SUSTAINED OSCILLATIONS IN REACTIVE SYSTEMS

137

3) P. Glansdorff and I. Prigogine, Thermodynamic Theory of Structure, Stability and Fluctuations (London, 1971). 4) W. Ebeling, Strukturbildung bei irreversiblen Prozessen (Leipzig, 1976). 5) A.M. Zhabotinsky, Biofizika 9 (1964) 306. 6) H. Higgins, Ind. Eng. Chem. 59 (1967) 19. G. Nicolis and J. Portnow, Chem. Rev. 73 (1973) 365. U.F. Franck, Angew. Chem. Int. Ed. Engl. 17 (1978) 1. 7) I. Prigogine and G. Nicolis, J. Chem. Phys. 46 (1967) 3542. I. Prigogine and R. Lefever, J. Chem. Phys. 48 (1968) 1695. 8) F. SchliSgl, Z. Phys. 243 (1971) 303; 248 (1971) 466; 253 (1972) 147. A. Nitzan, P. Ortoleva, J. Deutsch and J. Ross, J. Chem. Phys. 61 (1974) 1056. H. Haken, Introduction to Synergetics (Springer, Berlin, 1977). 9) W. Ebeling and R. Feistel, Z. Phys. Chem. (Leipzig) 257 (1976) 705; Ann. Phys. (Leipzig) 34 (1977) 81. 10) W. Ebeling and R. Feistel, Proc. ICNO (Berlin, 1975). 11) E.E. Sel'kov, Europ. J. Biochem. 4 (1968) 79. 12) V.A. Samoylenko and E.E. Sel'kov, Biofizika (Moscow) 17 (1972) 862. 13) K. Tomita and H. Tomita, Progr. Theor. Phys. 51 (1974) 1731. K. Tomita, T. Ohta and H. Tomita, Progr. Theor. Phys. 52 (1974) 484. 14) J.L. Ibanez, V. Fairen and M.G. Velarde, Phys. Lett. ~laA (1976) 364. 15) W. Ebeling, R. Feistel and J. Schmelzer, Wiss. Z.W. Pieck Univ., Rostock, MNR, in press. J.L. Ibanez and M.G. Velarde, Studies Appl. Math., in press. 16) L.D. Landau and E.M. Lifschitz, Statistical Physics (Oxford, 1957). 17) G. Nicolis and I. Prigogine, Selforganization in Nonequilibrium Systems (New York, 1977). 18) K.J. McNeil and D.F. Walls, J. Stat. Phys. 10 (1974) 439. C.W. Gardiner et al.. J. Stat. Phys. 14 (1976) 307. 19) G. Czajkowski and W. Ebeling, J. Non-Equilib. Thermodyn. 2 (1977) I. 20) S. Drosdziok, Z. Phys. 261 (1973) 431. 21) F.R. Gantmacher, Matrizenrechnung, Bd.II. (Berlin, 1971). 22) H.G. van Kampen, Stochastic Processes in Physics, Lecture notes, Utrecht (1970). 23) A. Uhlmann, in: Stochastische Theorie Nichtlinearer Irreversibler Prozesse (Rostock, 1977). 24) F. SchliSgl, Z. Phys. B25 (1976) 411. 25) H.K. Janssen, Z. Phys. 270 (1974) 67. 26) G. Czajkowski, Acta Phys. Polonica, in press. 27) R.L. Stratonovich, Topics in the Theory of Random Noise, vols. I, II (New York, London, Paris, 1967).