- Email: [email protected]
On the growth of nonequilibrium spatial correlations in a model reaction diffusion system
PHYSICA ELSEVIER
Physica A 231 (1996) 588-607
On the growth of nonequilibrium spatial correlations in a model reaction diffusion system Jerzy Gorecki ~,b,c, Kazuo Kitahara
c
Institute ~f Physical Chemistry, Polish Academy of Sciences and College ~f Science, Kasprzaka 44/52, PL-01-224 Warsaw, Poland b ICM, Banacha 2, PL-02-097 Warsaw, Poland c Department ¢~f'Applied Physics, Tokyo Institute ~['Technology, Meguro-ku, Tokyo 152, Japan
Received 13 December 1995
Abstract The mesoscopic description of a system with chemical reactions predicts that if the detailed balance condition is not satisfied then the nonequilibrium spatial correlations between concentrations of reactants may appear. The present work is concerned with the dynamics of their growth. We assume that the system is homogeneous at the beginning and study how fast the nonequilibrium correlations develop. The theory based on the master equation is compared with molecular dynamics simulations performed for a model system of "reacting" hard spheres. The qualitative characteristics of the growth of the spatial correlations, predicted by the theory based on a master equation, are confirmed in our simulations. However, this theory fails to give the correct quantitative description of this process.
1. Introduction It is known that nonequilibrium spatial correlations may appear in systems for which the detailed balance condition is not satisfied [1]. A few studies on this subject were concerned with the stationary correlations [2-5]. In this paper we focus our attention on the dynamics describing the growth of nonequilibrium correlations. We consider a system for which the average concentrations of reactants correspond to its stationary state. However, the reagents are randomly distributed in the system. In such a case the concentrations of reactants remain unchanged in time (they may fluctuate only), but the nonequilibrium correlations between reagents develop as reactions progress. A simple theory which describes the growth of fluctuations of correlations can be easily derived on the basis of a master equation. Its predictions are compared with molecular dynamics simulations, because this technique 0378-4371/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved PII S0378-4371 (96)00006-4
J. Gorecki, K. Kitahara/Physica A 231 (1996)588-607
589
allows us to study the correlations directly, by observation of positions of individual particles and measuring the distances between them. A closely related problem was concerned by Prigogine et al. [6] a few years ago, who observed how local correlations characteristic for a liquid develop in time if initially the particles form a triangular lattice. However, the theory describing correlations in a liquid is much more complex than this for correlations caused by the presence of chemical reactions and the authors of [6] did not attempt to give any quantitative description of the observed phenomenon, as we have done for the case of correlations originating from the presence of chemical reactions. 2. The model
The model which is considered below (compare also [7]) consists of the following reactions: kl
k2
X+X = X+Y, k-1
X + Y ~-- Y + Y . k-2
(1)
The same model was used in our recent paper [5] to study the stationary correlations of fluctuations in reagent's concentrations. From this work we know which is the form of the correlation function when time goes to infinity. Moreover, we have checked that the mesoscopic theory discussed below gives a good approximation for these correlations. The mesoscopic description of spatial correlations is based on the master equation for a spatially distributed system. Let us consider a system composed of cells characterized by an equal volume/2 and let Xi, Yi denote the number of molecules of X and Y in the ith cell. Information on the time evolution of the system can be extracted from the probability distribution P(..., Xi, Y/, ..., t) which describes the probability of finding Xi molecules of X and 1~ molecules of Y in the ith cell at the time t. The master equation for the system reads --P = dt
P +
P,
chem
(2)
diff
where
( )d d'7
P(...,Xi,Y}, .,t)= Z - \ - -{- 2~ k' 2l l x t i,x i chem
' '
1) + ~ - ~ X i Y i ~ X+i Y i
cells i
, y , i - - 1 ) ) P ( . .., y2k_2 - ~ 1 y i~
_/2kllx
Xi, Yi , ...,t)
+~--~-~( i+l)Xi+ ~.~( X i + I ) ( Y / - 1 ) +
(Xi-1)(Y/+I)+---~[,
)
i+1)I~
P(...,Xi+I,Y/-1
P(...,Xi-I,Yi+I
.... ,t) ...,t) (3)
590
J. Gorecki, K. Kitahara/PhysicaA 231 (1996) 5 8 8 ~ 0 7
and
d)
P(...,Xi, Yi,...,t)= ~-~,-(dfjXi +dfjYi)P(...,Xi, Yi,...,t) diff
cells i,j
+d~i(Xj + 1)P(..., Xj + 1, Yj,..., Xi - 1, Yi,..., t) +dfi(Yj + 1)P(...,Xj,Yj + 1, ...,Xi, Yi - 1, ...,t).
(4)
The factors 2kl, k_ l, ks and 2k_ 2 denote the probabilities that a corresponding reaction of (1) occurs within a time unit and d~"z3,dijYdescribe the probability of a diffusive jump of a particle of X or Y, respectively, between the ith and the jth cells. Calculating the average of Xi/f2 and Y//f2 one obtains the kinetic equations for concentrations of X and Y. In the following we shall denote these quantities by x ( t ) and y ( r ) respectively, where v instead of i denotes the position in our system. Considering the continuum limit, according to the standard approach [8], one may replace the diffusion term by the Laplacian of the appropriate concentration. Thus the averaging gives
dx -- klx2+k_lxy-k2xy+k_2y2+DV2x=f(x,y)+DV2x, dt dy d---t= klx~ - k-lxy + k2xy - k_2y 2 + DXY2y =g(x,y) + DU2y = -f(x,y) + DU2y,
(5)
where we assumed that the diffusion constants for X and Y are the same. It is clear that x + y = c(= const.) is a constraint of these kinetic equations. In the following we restrict our attention to a homogeneous, stationary state of the system (1). It is easy to prove that Eqs. (5) admit a single such state xs, y~; 0 _< x~, y~ < c, which is always stable. It may be shown [3] that if the detailed balance condition for reactions (1) is not satisfied in this state then the nonequilibrium spatial correlations between fluctuations of concentrations of reactants may appear. The stationary correlations have been studied in our recent paper [5]. Here we consider the question how do these correlations develop in time if at the beginning the reagents X and Y are perfectly mixed. The equations which describe the dynamics of the growth of spatial correlations of fluctuations in reactant's concentrations can be derived in a direct way from the master equation (F_x]. (2)). It is convenient to consider the correlation function of fluctuations of concentrations defined as
SNM(i,J,t):~2 ( (--N~ -- ns) (--~ --ms) >t
cells
(6) where M and N denote one of the reactants.
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
591
A tedious, but direct calculation leads to equations describing the dynamics of S y y (i, j, t), S r y ( i , j , t ) and S y y ( i , j , t ) . They are almost the same as those derived by Kitahara et al.[3], except of the fact that the time derivatives are taken into account,
-~sxxd
(r, r', t) = DV2,.Sxx(r,r',t) + D~72~,Sxx(r,r',t) - CxxS(v,r')
Of 2xx(r,r',t)+ +2 ~ ~.,y.
~
z.,y. (7a)
where C x x = klx, 2 + k2x,y8 + k_lx,y, + k_2ys 2,
.~Sxy(~.,r ' t) = D V 2 v S x y ( r , r ' , t ) + D V 2 , . , S x y ( r , r ' , t ) + C x x r ( r , r ') -
C~x(~,~',t)+
~
-
:vs,ys
xs,ys
xs,ys
2yy(7", r', t) ,
--~
(7b)
xs,ys
dsyr(r,
Of --2 ~
r ' , t ) = O V ~ , , C y v ( r , r', t) + DV2~,Syr(v,r',t) - Cxx6(r,v')
s r r (r, ~-',t) -
x,,y,
~
~..y~ (7c)
In the following we assume that the system is homogeneous and therefore the correlation functions should depend on the distance between r and v' only. Let us introduce new functions cr~, a~y and ~ruy defined as follows: ~(I
~ - ~' I,t) = ~ x x ( ~ , r ' , t )
-
xs~(r
-
~'),
~ ( I ~ - ~' I,t) ~xy(,,r',t), =
cruu(I r - r' I,t) = S r r ( r ,
r',t) - y,6(r
-
r').
(8)
In the functions O'nm defined above the fluctuations of concentration within a single cell are extracted. By calculating the Fourier transform of both sides of Eqs. (7) we obtain the following set of differential equations describing the dynamics of a,,,~ in the reciprocal space:
() Of
o'x~(q, t) + 2
(0,) ~y
O'~y(q,t)-C
-2D(27rq)2(rxx(q, t) ,
d
(9a)
( +(0~)
/0y,)~. ,,.)
~ryy(q,t)+C_2D(2rcq)2c~,y(q,t), :Vs~Ys
(9b)
J. Gorecki, K. Kitahara/Physica A 231 (1996)588-607
592
dayY(q't)=-2
-~x ~.,y. ~r~y(q, t) - 2
() NOf
.~,y. a y y ( q , t ) - C
-2D(27rq)2 ayy (q, t) ,
(9c)
where
Crnm(q, t) = f exp (--27rirq)crnm(r , t)dr J and
C = 2(klx~ - k - l x ~ y s ) .
(10)
If C equals to zero then the stationary (in a long time limit) functions a,,,~ converge to zero. It may be shown that this happens if the detailed balance condition is satisfied for reactions (1). This feature justifies calling the appearing correlations nonequilibrium. The stationary solution of Eqs. (9) is well known [1] and it reads
~(q, ~) = ~(q, ~) = -~(q,
~) =
C 1 2D n2 + (2~rq) 2 '
(11)
which in the space variables corresponds to C 1 a x x ( r , oc) = cryy(r, oo) = --axu (r, oo) -- - 87rD [ r I exp (-to I r ' I).
(12)
The constant n for reactions (1) is defined as follows:
= 1(2klx~
- ( k _ l - k : ) ( c - 2x~) + 2 k _ 2 ( c - x~))
(13)
and it is positive because the steady state of system (l) is a stable one. The relationship c%~([ r -
r ' I,t) + 2a~u( I r - r ' h t) + (ruu(I r -- r ' I,t) = 0
(14)
which comes out of Eqs. (9) is hold for any time t provided that it is satisfied by the initial condition. It may be written as
~xs(~,
~',t) - x~6(r - ~') + 2 Z x r ( , ~ , ,.', t) + ~ y y ( ~ , ~', t) - u86(~ -
~')
which expresses the fact that if the chemical identities of particles are neglected then no spatial correlations are present in the system considered. In molecular dynamics simulations performed at liquid densities, the equilibrium correlations caused by a hard sphere potential have to be extracted from the total correlations between reactants, in order to have relationship (15) satisfied.
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
593
For reactants which are perfectly mixed the initial condition reads
a z x ( q , t = 0 ) = a u y ( q , t = 0) = azy(q,t = O) = O.
(16)
All three independent solutions of the set of linear differential equations (9) can be easily found, however, because of the form of the source term and the initial condition, only the solution for which
crxz(q , t) = cr~y(q, t) = -~rzy(q, t) = -(r(q, t)
(17)
is satisfied develops in time and the amplitudes of the other solutions are equal to zero. Using condition (17) one may reduce each of Eqs. (9) to the form
C - 2D(27rq)2o'(q, t).
(18) The solution of Eq. (18) reads o'zx(q) = O'uu(q)
=
-O'zu(q)
C 1 C exp(-2D(x 2 + (27rq)2)t) + 2D ~2 + (2~q) ~ 2D ~2 + (2~q) ~
=
(19)
The inverse Fourier transform of Eq. (19) can be done analytically (compare [9]) leading to the expression C 1 = ~yy(r) = -O~v(r ) _ -87rD i r i exp ( - ~ 17" I)
~(r)
oo
1
t)fdq
C
q sin(27rrq) exp(-8DtTr2q 2) + (27rq) 2
2
o
C 1 C 1 -- - 8 7 r D l r ] e x p ( - ~ ]r ])+ 87r----~r----~(sinh i (~ Iv 1)+ ½exp(-n] r 1)
2 2v5-
/
2
(20)
where x
o/x,_0
Thus Eq. (20) gives us complete information on the time evolution of the correlations of fluctuations.
594
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
3. Molecular dynamics simulations of nonequilibrium spatial correlations in a system with chemical reactions The molecular dynamics technique applied in this paper is very similar to the one used in our previous paper [5] to simulate the stationary correlations of fluctuations. The technique is based on the concept of "reactive" hard spheres (see for example [10]). According to it all the reagents are represented by identical hard spheres (masses and diameters are the same) and their motion is ruled by a standard dynamics with all the collisions regarded as elastic. The additional parameter, which has no influence on the mechanical motion, labels the chemical properties of spheres. We consider both reactive and nonreactive collisions between spheres. A chemical reaction may occur if particles representing appropriate reagents collide. The probability of a reactive collision is usually lower than 1. In order to describe a stochastic character of chemical processes a random number is called for each collision which may lead to a reaction and its value is compared with the steric factor for the appropriate process. If the value is not larger than the steric factor the collision is regarded as a reactive one and the chemical identity parameters of the spheres involved are modified. In simulations we assume that the activation energies for all reactions (1) are equal to zero. The technique may be easily adopted for any activation energy of the reactions provided that the computation is long enough to give a large number of reactive collisions. However, it is well known [11,12] that every thermally activated chemical reaction creates a nonequilibrium state of the velocity distribution of the reagents involved and this effect is not taken into account by the theory based on the master equation presented in the previous chapter. The kinetic energy is the only form of energy in a system of hard spheres. As all reactive collisions are elastic from a mechanical point of view no kinetic energy is released or consumed when a reaction occurs and the system as a whole remains in the thermal equilibrium with respect to the translational motion. Its temperature is related to the average energy per sphere. It is very important from the computational point of view to maintain such equilibrium because it allows us to extend the size of simulations using a prerecorded equilibrium trajectory. A trajectory, which was calculated for a system characterized by the periodic boundary conditions, may be used to enlarge the size of simulations [5,13]. The original small system may be periodically expanded in any of the directions by any integer number of the box lengths. Of course, if a chemical identity of molecules is neglected then such expansion does not bring us any new information. However, for a multicomponent chemical system, in which the translational motion is not related to chemical identity, the situation is different. First, a different chemical composition may be initialized in various boxes by marking the equivalent (by periodicity) spheres in a different way. Secondly, a steric factor (if it is not equal to unity) differentiates the time evolution in various boxes, because a collision between the equivalent (by periodicity) objects may be reactive in one box and nonreactive in another. Moreover, the periodic expansion simplifies calculations of spatial correlations, because it is enough to calculate interesting intermolecular distances between spheres within a single box. Of course the correlations whose range extends over a single box may be affected by
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
595
the periodicity of our system. In the following we show the results which have been obtained by periodic extension of the following systems: (I) A low density system. The original simulations were done for N = 500 hard spheres placed in a cubic box with the side length d = 14.7cr, a being the diameter of a sphere. The packing fraction is r / = 0.0824. This system was expanded by 10 box lengths in each direction so the total number of molecules is 500000. (II) A medium density system (N = 1000 spheres, d = 14a, r / = 0.191). This system was expanded by 10 box lengths in each direction leading to 106 spheres. (III) A high density system (N = 1000 spheres, d = 12.Sir, ~ = 0.268). This system was also expanded by 10 box lengths in each direction and the total number of spheres is 106. The average time between successive collisions of a particle is equal to 1.29 ps, 0.405ps and 0.223 ps in systems I, II and III, respectively. In the following we give all values of time in these units. In all these systems reactants X and Y are represented by spheres with the mass (rn = 32 a.u.) and diameter (a = 5.~). The temperature equals 300 K. The basic set of steric factors (probabilities of a reactive collisions) for processes (1) is the same as used in [5], 81 = 0.8,
8_ 1 = 0.2,
S2 = 0.1,
S-2 = 0.8,
but we have also performed simulations for slower reactions with steric factors, s t = 0.4,
J---1 = 0.1,
I
s 2 = 0.05,
sl_2 = 0.4,
and s I/1 = 0.16,
.~,11 _1
= 0.04,
s It2 = 0.02,
811 -2
= 0.16.
For all these sets of steric factors the concentration of X at its stationary state equals to • ,
=
7)c.
At the beginning of simulations the chemical identities of spheres are assigned in a random way, such that the concentrations of X and Y correspond to their stationary values zs, ys. Because of the random distribution of reactants we may expect that the no spatial correlations of concentrations are present in the initial state. The simulations performed for the stationary correlations [4,5] indicate that they have a limited range, which is of the order of the mean free path. In order to save computer time and cover the most important part of correlation functions we have introduced a cutoff R = 1.8cr in the range of intermolecular distances. The interval of interesting distances [cr, R] is divided into a number of subintervals of the same length A (we have used A _-- 0.02~r or 0.05tr). Moreover, we have introduced interval of times [0, to], during which the system evolution is observed. This interval is divided into subintervals of the length At each. Simulations give us the average number of spheres representing the reactant N whose distance from a sphere representing reactant M lies within [~r + iA, tr + (i + 1)A] at the condition that the time variable belongs to the subinterval [jAt, (j + 1)At]. Let us denote this quantity a s GNM (i, .j). If in simulations the time variable exceeds to then the parameters describing the chemical identities of spheres
J. Gorecki,K. Kitahara/PhysicaA 231 (1996)588~507
596
are mixed in a random way in order to destroy existing correlations and a new reaction path starts. The results presented below represent the average over a few hundred reaction paths. The correlations between particles g(r, t) for the system as a whole (i.e. with neglected chemical identity of spheres) can be obtained by scaling G(i, j) by the value, calculated assuming a uniform distribution of spheres, g([cr + iA, ~r + (i + 1)A], [jAr, (j + 1)At])
G(i,j) = 47r((~r + (i + 1)A) 3 - (a + iA)Z)cAt '
(21)
where G(i, j) denote the average number of spheres within the space interval [~r + iA, cr + (i + 1)A] and the time interval [jAt, (j + 1)At], without consideration of their chemical identity. Of course g should be the same for all j, but because of limited statistical samples some fluctuations are possible. Similarly we can calculate the partial correlation functions g x x , gxY and g y y ,
g x x @ r + iA, o" + (i + 1)A], [jAr, (j + 1)At]) Gxx(i,j) = 47r((o" + (i + 1)A) 3 - (o" + iA)3)x,A, '
(22)
g x y ( [ a + iA, o" + (i + 1)A], [jAr, (j + 1)At]) G x y (i, j) = %r((a + (i + 1)A) 3 - (a + iA)3)x~At '
(23)
gyy([~r + iA, ~7+ (i + 1)A], [jAr, (j + 1)At]) ayy(i,j) = ~Tr((~r + (i + 1)A) a - (~r + iA)Z)y~A~
(24)
If the nonequilibrium spatial correlations of fluctuations of reagents concentration are present then the partial correlation functions are distinctly different from the total correlation function. However the functions obtained from Eqs. (22)-(24) contain information on both correlations introduced by the hard potential of spheres and on correlations related to the chemical processes. In order to compare the nonequilibrium spatial correlations observed in simulations with the theory (compare condition (15)), one needs to separate the equilibrium correlations introduced by the hard sphere potential from a partial correlation function. Having in mind that the accuracy of gNM/g is higher than gNM itself because fluctuations in samples of interatomic distances cancel out, it is convenient to calculate the chemical contribution to correlations of fluctuations in concentration in the following way: cr,,m([ r - r '
I , t ) = n , m , a0(I ,"
(gNM(I r -- r' I,[JAt, (j + 1)At]) _ 1) I) \ g--~r---r;]-,[jAt,(j÷ 1)At]) /
=
go([ r - r' l) \
- ~ ; - _ - ~ ; I-, [jAr, (j + 1)At])
'
where go (r) is the radial distribution function for the system of spheres at the given packing fraction (it may be calculated separately with high accuracy) and t denotes a point within
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
597
the selected interval of times, for example t = ½(2J + 1)At. Our previous study [5] has shown that the nonequilibrium correlations between concentrations of reagents calculated on the basis of Eq. (22) have the functional form given by Eq. (12) and the parameters x and C are given by Eq. (13) and Eq. (10), respectively. The rate constants may be calculated from steric factors and the collision frequency. D may be obtained from the expression valid for a hard sphere fluid characterized by the packing fraction r/, 1 D - 16
~r ~go~¢ +)Vmr/-~RT'
(25)
where the value of the radial distribution function at the point of contact of spheres go ((~+) can be calculated from the formula 1 + 4r/go(o"+) = (1 + 77+ 7/2 + r/a)
(27)
(1 - ~/)a
Alternatively one may get D directly from simulation data by measuring the average square of the particle shift as a function of time. However the decay constant n calculated using Eq. (13) is usually smaller (by a factor of 2) than the observed rate of decay of correlations at short distances (and this is the region we focus our attention on). On the other hand the "theoretical" value of a gives a good estimation of the decay of ¢rnm at large distances. Here, being interested in correlations appearing at short distances, we use the numerical fit for g rather than the result of Eq. (13). Studying the stationary nonequilibrium correlations in system (1) we have also found [5] that the value of the constant C predicted by the theory is usually much smaller than the real strength of fluctuations. We believe that it is associated with the excluded volume effect, which is not taken into account by the simple theory presented in the previous section. From the practical point of view it means that C should be rather fitted to the data than calculated from Eq. (10). A typical time evolution of the functions ¢,~,~ obtained from simulations is presented in Figs. la-c. These functions are plotted using dashed line for four selected moments of time (the dash length is an increasing function of time) and compared with the stationary correlation function o'nm (solid line). It can be noticed that short distance correlations grow first and the long distance correlations appear later. As the theory says, a~y has an opposite sign to the other functions ( ¢ ~ , (Txy). The absolute values of all functions ¢r~ (Fig. la), ¢~u (Fig. lb) and ~ruu (Fig. lc) are almost the same for all moments of time, as predicted by Eq. (20). This feature indicates the accuracy of simulations. An evolution of rr~u calculated from Eq. (20) for three first moments of time, for which the data in Figs. la-c were shown, is presented in Fig. ld. The result for the longest time is almost the same as for the second longest so we have not shown it. The parameters required by Eq. (20) were the same as used in [5] to fit the stationary correlation function to the form given by Eq. (13) (~ and C) or calculated from simulations (D). The basic feature of observed correlations dynamics is correctly reflected by the theory: both methods indicate that the short distance correlations are growing first (and this seems natural, as they are caused by chemical reactions, which occur when spheres are in touch)
J. Goreck~ K. Kitahara/ Physica A 231 (1996)588-607
598
(a) }
.
.
. . . .
X
•
. . -
.
.
.
,
.
-0.010
1
I:D C
:/
X cT -0,030 (b O u9 C__ O
~:
:
1.1
1.2
f
j
--
--
-0,050
© (D L L
0 (9
-0.070 1 0
1.5
1.4
1.5
1.6
distonce /sphere diometer
(b) >~ 0 . 0 6 0 ~D ID X
0.040 r,b (f] O -'~ 0 . 0 2 0 © L
0 (D
0.000
10
1.1
1.2
1.3
1.4
1.5
1.6
distonce /sphere diorneter Fig. 1. The comparison of time evolution of (a) axz/c 2, (b) ~x~/c 2 and (c) o'~y/c 2 obtained in molecular dynamics simulations with a/c 2 calculated using Eq. (20). The simulations were performed for the system I with the fastest reaction. The correlations are shown for three moments o f time: 0.116 (the dotted line), 0.271 (the short dashed line), 0.426 (the medium dashed line), 0.814 (the long dashed line). The solid line represents the stationary correlation function. (d) shows a/c 2 calculated from F_.q. (20) for first three moments of time. The following values were used to evaluate F-4. (17): D = 0.625cr2/ps, n = 2.4~r - 1 and C = 1 5 D .
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
599
(c) ....... ~............ . .................
i
>- - 0 . 0 1 0 %)
i ............
C ©
:
L
i
i.Z[/:.::.i-:.:-:.:-~:-:-:
~ [ . . . .
'
I-
r--
~
--
--
I
,---
)C
-0.030 © ffD oO C 0
--0.050
C~
o
©
-0.070
~llll~l~ll~lllll~rlllllllFii
1.0
1.1
1.5
1.2
distance
1.4
/sphere
1.5
1.6
diameter
(d) >- 0.060 K?
C 0
)< C 0.040
L) "" \t C 0 ~-~ 0 . 0 2 0
" '"" ".,
"'"%
x
@ o
0 0000 .0
1.
1.
1.4
1.5
distance /sphere diameter Fig. 1. Continued.
1.6
J. Gorecki,K. Kitahara/PhysicaA 231(1996)588-607
600
and then, as time progress, they spread out. However the comparison of Fig. 1b with Fig. 1d shows that the theory underestimates the short distance part of o, whereas it overestimates the growth of correlations at large distances. In order to introduce more detailed description of correlations dynamics let us introduce the scaled correlation function in the form
=
cr,~,~(r, t = oo) "
(28)
It follows from Eqs. (13) and (20) that u does not depend on reactants N and M chosen. As the time progresses u(r, t) changes from 0 to 1 for every r. Referring to Fig. 1 we expect that u(r, t) approaches unity for small r first and at large r later. The dynamics of u(r, t) is governed by the following equation:
1 0 u ( r , t ) _ 1 0 ( ~_~ ) (1) O Ot r 2 Or r2 u(r,t) + 2 x + +47rr exp(t~r)(1 -
u(r, t)),
~___~
,(r,t) (29)
which results directly from Eq. (7) and the form of the stationary correlations (Eq. (12)). As expected u(r, t) = 1 is the stationary state of Eq. (29). Using Eq. (29) seems to be more convenient than Eq. (20) because u(r, t) does not depend on the amplitude C as it scales out. u(r, t) as a function of distance for a few selected moments of time is presented in Figs. 26. Fig. 2a shows the growth of correlations for a system characterized by a low density (I) and fast reactions (si). Comparing the results of simulations (Fig. 2a) with theory (Fig. 2b) one may notice that for the first of moments of time considered (t = 0.116, the dotted line) our theory underestimates correlations at short distances, but overestimates them at large ones. The theory predicts a rapid growth of correlations in the interval t[0,116, 0.426], much faster than measured from simulations. According to the theory at t = 0.426 the nonequilibrium correlations of reactants are fully developed, whereas simulations show that for distances > 1.5~r they reached only half of their stationary values. Figs. 3 and 4 present results for systems at a medium density (II). If the reaction is fast (si) the differences between simulations (Fig. 3a) and Eq. (20) are enormous. For time t = 0.242 (solid line) the correlations are well developed according to simulations, but theory says that they are hardly visible. Their evolution observed in simulations remains an expanding front, whereas Eq. (20) predicts an uniform growth in the whole space. The results obtained for five times slower reactions (Figs. 4a and 4b; the set s~~) are in better agreement, but also in this case the rate of growth of correlations is underestimated by theory. The same kind of disagreement between simulations and theory may be observed in results for the most dense system (HI). The data for a fast reaction (si) are compared in Figs. 5a,b, for the slower one (s~) in Figs. 6a,b.
J. Gorecki, K. Kitahara/Physica A 231 (1996)588-607
601
(a) >~:~ 1.0 ©
L__
q--
q
x
'
- - -
L
_
_
@8
I
1 I
© q)
I --
x~06
I I
/
If)
t
_
0
O0.4 q) L
o 0.2 q)
o (p O@ 10
1.2
11
t.5
1.4
/.5
1.6
distonce /sphere di@meter
(b) 1.0
(f) C 0
0.8
© 0.6 L
© 0
g 0.4 o
",,
L)
.....
02
00
r
1.0
,
~
i
I
1.1
i
~
i
i
1.
2/
P
i
,
~
I
15
i
~
i~
I
1.4
i
f--i
i
15 q
1.
t
i
,
1 6
distonce /sphere diorneter Fig. 2. Comparison of (a) molecular dynamics simulations for the system 1 with rates si with (b) the theory. (a) The results of simulations are shown for four moments of time: 0.116 (the dotted line), 0.271 (the short dashed line), 0.426 (the long dashed line) and 0.814 (the solid line); (b) calculations are for O. 116 (the dotted line), 0.271 (the short dashed line) and 0.426 (the long dashed line). The values of parameters used are D = 0.625a2/ps and ~ = 2.40 "-1.
J. Gorecki, K. Kitahara/Physica A 231 (1996)588-607
602
(a) 1.0 C 0
08
b~ C 0
o0.4 -g k_
o 02 72)
©
@@
lo
1.1 1.2 15 distance /sphere
1.4 1.s diameter
(b} 1.0
tO
0.8
C 0 o
L_ L
0.6
o
"'.
(D
"'"
0.4 o
o
O3
0.2
" 4
0.0
l l ; J l l l l l l l
10
1.1
1.2
,1. 13
, ,,
1.
14,,,
~ I 1.5
,
, ,
, 1.6
distance /sphere diameter Fig. 3. The time evolution of the scaled correlation function v(r, t); comparison of (a) molecular dynamics simulations for the system II with rates s i with (b) the theory. (a) The results of simulations are shown for four moments of time: 0.054 (the dotted line), 0.104 (the dashed line) and 0.242 (the solid line); (b) calculations are for 0.242 (the solid line), 0.494 (the short dashed line) and 0.741 (the long dashed line). The values of parameters used are D = 0.229a2/ps and i¢ ---- 4 . 8 6 a -1 .
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
603
(a) x 1.0 C
x 08 C ©
© ~0.6 ...... : i ';
U3
i
©
-1 -i :
L i
4~
(30.4
~L~
LL
t _
L k~
i
© (.D
L_
I
'1 i_ I_
0.2
i -%
©
~ -!
-I
L ....
I_
. . . . . . . . . . .
0.0
10
. . . .
11~
. . . .
~
.... t
'
'
.=S ~ -2
-
_-
-
-
I .i2.... 1. 13 .... 1.4I ....
r q
1.
-
--
15 . . . . 1.6
distance /sphere diameter
(b) 1.0
0.8 (/3
o © ~0.6 L
o (.3 ~0.4 U
0.2
0.0
i-(-?-(?-£rr
.0
1.1
I
1.2
~FIilll~lllEillliS
1.3
1.4
1.5
1.6
distance /sphere diameter Fig. 4. The time evolution of the scaled correlation function u(r, t); comparison of (a) molecular dynamics simulations for the system 11 with rates s~' with (b) the theory. (a) The results of simulations are shown for four moments of time: 0.074 (the dotted line), 0.272 (the short dashed line), 0.519 (the long dashed line) and 1.21 (the solid line); (b) calculations are for 0.272 (the short dashed line), 0.519 (the long dashed line) and 1.21 (the solid line). The values of parameters used are D : 0.229a2/ps and n = 2.35(7 -1 .
604
7. Gorecki, K. Kitahara/Physica A 231 (1996)588-607
(a) )%) 1,0 C ©
X C © ©
0.8 t
i_~1 I
~06
'
(f)
1 I
--i i i
•
O
i
©0,4
L
L
(2) ~D
i
O2 U q) (D q)
,_
I
0@
1•0
11
12
distance
/sphere
114
I I5
diameter
(b) 1.0
(/) 0 8
-.--
0
'~"~)(~)0.6
"-""""
'"-.
~0.4
00
i
0
p ~ F
J
i
F
I
i
J
i
i
i
ITt
I
I ~
1.~ ~2 distance / s p h e r e
I
J
l
I
I
JTI
l
,
lq
1.s dTorneter
~
T T r
1.4
Fig. 5. The time evolution of the scaled correlation function t,(r, t); comparison of (a) molecular dynamics simulations for the system III with rates sl with (b) the theory. (a) The results of simulations are shown for four moments of time: 0.247 (the dotted line), 0.471 (the short dashed line), 0.919 (the long dashed line) and 1.10 (the solid line); (b) calculations are for 0.919 (the long dashed line), 1.10 (the solid line) and 1.35 (the short dashed line). The values of parameters used are D = 0.131~ 2 / p s and n = 7 . 9 a - 1 .
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
605
(a) x 1.0 C
X C © (L>
0.8
i-; I i
ii L], :
W)
i
,
t ~
C 0 i
©0.4
I- I
&) L_ m
o
i
,__
L
~
k_ 1
0.2 -o © © 0 (/3
"
, ~.. , 0.0
,
0
i
i
,
1 I1
i
]
I
7.].'~-, ............ q=..-~,-, T
i
I
12
~
J
i
i
L~ ] ~ ........~ - : ' =~i- -V..~- ~ ..,+- ]
I
1.,3
i
i
,
,
I
r - i-'-i 7]
- -
1.4
"'""ST"]"
~~
i
-6
distance /sphere diameter
(b) 1.0
u) C 0
0.8
© ~0.6 L
o t)
~0.4 I
©
k) u)
02
O0 10
1,1
1.2
1 5
1.4
1.5
1.6
distance /sphere diameter Fig. 6. The time evolution of the scaled correlation function v(r, t); comparison of (a) molecular dynamics simulations for the system III with rates s~ with (b) the theory. (a) The results of simulations are shown for four moments of time: 0.247 (the dotted line), 0.471 (the short dashed line), 0.919 (the long dashed line) and 1.10 (the solid line); (b) calculations are for 0.919 (the long dashed line), 1.10 (the solid line) and 1.35 (the dotted line). The values of parameters used are D = 0.131a2/ps and ~ = 5.6o "-1 .
606
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
4. Conclusions In this paper we have presented results of molecular dynamics simulations of the growth of nonequilibrium spatial correlations of reactants fluctuations in a system with chemical reactions. The considered systems exhibited such correlations because the detailed balance condition was not satisfied. Our simulations started from well stirred reactants, thus the correlations were restored during the progress of reactions. We observed that the correlations of fluctuations appear at short distances first and then spread out in space. We considered correlations at distances not exceeding a single molecular diameter, because in this range they are the most important. In this range of distances the growth of correlations is very fast and it is shorter than the average time for successive collisions of the same particle. These conclusions agree with observations of Prigogine et al. [6] who found that the radial distribution function (RDF) of a liquid is first reestablished at short distances. This process is also very fast and the first peak of RDF appears in less than a single collision time. We have also derived a simple theoretical description of the growth of nonequilibrium spatial correlations in a system with chemical reactions, based on the master equation for reaction-diffusion processes. It has already been shown that an analogous description of the stationary (independent on time) correlations is in a good agreement with molecular dynamics simulations [5] and with simulations based on the Bird technique [4]. This result is quite surprising because on one hand the theory is based on the large volume expansion, but on the other hand the range of correlations is of the order of the mean free path, which, in some of the systems studied [5], may be a fraction of the sphere diameter. Simulations show that this simple theory gives the qualitative description of correlations growth only. Both theory and simulations indicate that the correlation function approaches its stationary form at short distances first and then it spreads o u t - correlations at large distances approach their stationary values later. However, the quantitative description of the correlation dynamics is different. Molecular dynamics suggests front-like expansion of the stationary correlations, whereas the theory supports more uniform growth at all distances. We believe that the discrepancy comes mainly from the fact that the diffusive term in the master equation is replaced by the Laplace operator. It would be too naive to expect that a classical diffusive term may correctly describe processes which occur at times for which the motion of particles is dominated by a ballistic flow. In the case of chemical wave front propagation, generated by a fast reaction, it may be shown that the reaction-diffusion equation need to be associated with the equation describing the relaxation of a diffusive flow in order to get a correct wave speed [14]. The reactions considered in this paper were very fast and the diffusive flow relaxation seems also important.
Acknowledgements This work was supported by the grant KBN 2 P303 018 05 provided by the Polish State Committee for Scientific Research and by Grant-in-Aid in Priority Areas 07240106 from
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588--607
607
Ministry of Education, Culture and Science of Japan. One of the authors (J.G.) is grateful for JSPS support during his visit to Tokyo Institute of Technology. The authors would like to thank Professor U.M. Titulaer for helpful discussions. References [1] Y. Kuramoto, Prog. Theor. Phys. 52 (1974)711; G. Nicolis and M. Malek Mansour, Phys. Rev. A 29 (1984) 2845. [2] K. Kitahara, K. Seki and S. Suzuki, J. Phys. Soc. Japan 59 (1990) 2309. [3] K. Kitahara, K. Seki and S. Suzuki, in: Proceedings of Molecular Dynamical Processes, A. Blumen et al., eds. (World Scientific, Singapore, 1990)p. 403. [4] G. Nicolis, A. Amellal, G. Dupont and M. Mareschal, Journal of Molecular Liquids 41 (1989) 5. [5] J. Gorecki, K. Kitahara, K. Yoshikawa and I. Hanazaki, Physica A 211 (1994) 327. [6] I. Prigogine, E. Kestemont and M. Mareschal, in: Microscopic Simulations of Complex Flows, NATO ASI series, Vol. 236, M. Mareschal, ed. (Plenum Press, New York, 1990). [7] J.S. Turner, J. Phys. Chem. 81 (1977)2379. [8] C.W. Gardiner, Handbook of Stochastic Methods (Springer, Berlin, 1983). [9] I.S. Gradshteyn and I.M.Ryzhik, formula 3.954.1, Table of Integrals, Series and Products (Academic Press, New York, 1992). [10] J. Gorecki and J. Gryko, Computer Phys. Comm. 54 (1989) 245. [11] I. Prigogine and E. Xhrouet, Physica 15, (1949) 913; for the recent see J. Popielawski, in: The Dynamics of Systems with Chemical Reactions, J. Popielawski, ed. (World Scientific, Singapore, 1989), and references therein. [12] J. Gorecki and B.C. Eu, J. Chem. Phys. 97 (1992) 6695. [13] J. Gorecki, in: Far from Equilibrium Dynamics of Chemical Systems, J. Gorecki et al., eds. (World Scientific, Singapore, 1994); Molecular Physics Reports 10 (1995) 48. [14] J. Gorecki and J. N. Gorecka, Acta Physica Polonica B 26 (1995) 1177.
Physica A 231 (1996) 588-607
On the growth of nonequilibrium spatial correlations in a model reaction diffusion system Jerzy Gorecki ~,b,c, Kazuo Kitahara
c
Institute ~f Physical Chemistry, Polish Academy of Sciences and College ~f Science, Kasprzaka 44/52, PL-01-224 Warsaw, Poland b ICM, Banacha 2, PL-02-097 Warsaw, Poland c Department ¢~f'Applied Physics, Tokyo Institute ~['Technology, Meguro-ku, Tokyo 152, Japan
Received 13 December 1995
Abstract The mesoscopic description of a system with chemical reactions predicts that if the detailed balance condition is not satisfied then the nonequilibrium spatial correlations between concentrations of reactants may appear. The present work is concerned with the dynamics of their growth. We assume that the system is homogeneous at the beginning and study how fast the nonequilibrium correlations develop. The theory based on the master equation is compared with molecular dynamics simulations performed for a model system of "reacting" hard spheres. The qualitative characteristics of the growth of the spatial correlations, predicted by the theory based on a master equation, are confirmed in our simulations. However, this theory fails to give the correct quantitative description of this process.
1. Introduction It is known that nonequilibrium spatial correlations may appear in systems for which the detailed balance condition is not satisfied [1]. A few studies on this subject were concerned with the stationary correlations [2-5]. In this paper we focus our attention on the dynamics describing the growth of nonequilibrium correlations. We consider a system for which the average concentrations of reactants correspond to its stationary state. However, the reagents are randomly distributed in the system. In such a case the concentrations of reactants remain unchanged in time (they may fluctuate only), but the nonequilibrium correlations between reagents develop as reactions progress. A simple theory which describes the growth of fluctuations of correlations can be easily derived on the basis of a master equation. Its predictions are compared with molecular dynamics simulations, because this technique 0378-4371/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved PII S0378-4371 (96)00006-4
J. Gorecki, K. Kitahara/Physica A 231 (1996)588-607
589
allows us to study the correlations directly, by observation of positions of individual particles and measuring the distances between them. A closely related problem was concerned by Prigogine et al. [6] a few years ago, who observed how local correlations characteristic for a liquid develop in time if initially the particles form a triangular lattice. However, the theory describing correlations in a liquid is much more complex than this for correlations caused by the presence of chemical reactions and the authors of [6] did not attempt to give any quantitative description of the observed phenomenon, as we have done for the case of correlations originating from the presence of chemical reactions. 2. The model
The model which is considered below (compare also [7]) consists of the following reactions: kl
k2
X+X = X+Y, k-1
X + Y ~-- Y + Y . k-2
(1)
The same model was used in our recent paper [5] to study the stationary correlations of fluctuations in reagent's concentrations. From this work we know which is the form of the correlation function when time goes to infinity. Moreover, we have checked that the mesoscopic theory discussed below gives a good approximation for these correlations. The mesoscopic description of spatial correlations is based on the master equation for a spatially distributed system. Let us consider a system composed of cells characterized by an equal volume/2 and let Xi, Yi denote the number of molecules of X and Y in the ith cell. Information on the time evolution of the system can be extracted from the probability distribution P(..., Xi, Y/, ..., t) which describes the probability of finding Xi molecules of X and 1~ molecules of Y in the ith cell at the time t. The master equation for the system reads --P = dt
P +
P,
chem
(2)
diff
where
( )d d'7
P(...,Xi,Y}, .,t)= Z - \ - -{- 2~ k' 2l l x t i,x i chem
' '
1) + ~ - ~ X i Y i ~ X+i Y i
cells i
, y , i - - 1 ) ) P ( . .., y2k_2 - ~ 1 y i~
_/2kllx
Xi, Yi , ...,t)
+~--~-~( i+l)Xi+ ~.~( X i + I ) ( Y / - 1 ) +
(Xi-1)(Y/+I)+---~[,
)
i+1)I~
P(...,Xi+I,Y/-1
P(...,Xi-I,Yi+I
.... ,t) ...,t) (3)
590
J. Gorecki, K. Kitahara/PhysicaA 231 (1996) 5 8 8 ~ 0 7
and
d)
P(...,Xi, Yi,...,t)= ~-~,-(dfjXi +dfjYi)P(...,Xi, Yi,...,t) diff
cells i,j
+d~i(Xj + 1)P(..., Xj + 1, Yj,..., Xi - 1, Yi,..., t) +dfi(Yj + 1)P(...,Xj,Yj + 1, ...,Xi, Yi - 1, ...,t).
(4)
The factors 2kl, k_ l, ks and 2k_ 2 denote the probabilities that a corresponding reaction of (1) occurs within a time unit and d~"z3,dijYdescribe the probability of a diffusive jump of a particle of X or Y, respectively, between the ith and the jth cells. Calculating the average of Xi/f2 and Y//f2 one obtains the kinetic equations for concentrations of X and Y. In the following we shall denote these quantities by x ( t ) and y ( r ) respectively, where v instead of i denotes the position in our system. Considering the continuum limit, according to the standard approach [8], one may replace the diffusion term by the Laplacian of the appropriate concentration. Thus the averaging gives
dx -- klx2+k_lxy-k2xy+k_2y2+DV2x=f(x,y)+DV2x, dt dy d---t= klx~ - k-lxy + k2xy - k_2y 2 + DXY2y =g(x,y) + DU2y = -f(x,y) + DU2y,
(5)
where we assumed that the diffusion constants for X and Y are the same. It is clear that x + y = c(= const.) is a constraint of these kinetic equations. In the following we restrict our attention to a homogeneous, stationary state of the system (1). It is easy to prove that Eqs. (5) admit a single such state xs, y~; 0 _< x~, y~ < c, which is always stable. It may be shown [3] that if the detailed balance condition for reactions (1) is not satisfied in this state then the nonequilibrium spatial correlations between fluctuations of concentrations of reactants may appear. The stationary correlations have been studied in our recent paper [5]. Here we consider the question how do these correlations develop in time if at the beginning the reagents X and Y are perfectly mixed. The equations which describe the dynamics of the growth of spatial correlations of fluctuations in reactant's concentrations can be derived in a direct way from the master equation (F_x]. (2)). It is convenient to consider the correlation function of fluctuations of concentrations defined as
SNM(i,J,t):~2 ( (--N~ -- ns) (--~ --ms) >t
cells
(6) where M and N denote one of the reactants.
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
591
A tedious, but direct calculation leads to equations describing the dynamics of S y y (i, j, t), S r y ( i , j , t ) and S y y ( i , j , t ) . They are almost the same as those derived by Kitahara et al.[3], except of the fact that the time derivatives are taken into account,
-~sxxd
(r, r', t) = DV2,.Sxx(r,r',t) + D~72~,Sxx(r,r',t) - CxxS(v,r')
Of 2xx(r,r',t)+ +2 ~ ~.,y.
~
z.,y. (7a)
where C x x = klx, 2 + k2x,y8 + k_lx,y, + k_2ys 2,
.~Sxy(~.,r ' t) = D V 2 v S x y ( r , r ' , t ) + D V 2 , . , S x y ( r , r ' , t ) + C x x r ( r , r ') -
C~x(~,~',t)+
~
-
:vs,ys
xs,ys
xs,ys
2yy(7", r', t) ,
--~
(7b)
xs,ys
dsyr(r,
Of --2 ~
r ' , t ) = O V ~ , , C y v ( r , r', t) + DV2~,Syr(v,r',t) - Cxx6(r,v')
s r r (r, ~-',t) -
x,,y,
~
~..y~ (7c)
In the following we assume that the system is homogeneous and therefore the correlation functions should depend on the distance between r and v' only. Let us introduce new functions cr~, a~y and ~ruy defined as follows: ~(I
~ - ~' I,t) = ~ x x ( ~ , r ' , t )
-
xs~(r
-
~'),
~ ( I ~ - ~' I,t) ~xy(,,r',t), =
cruu(I r - r' I,t) = S r r ( r ,
r',t) - y,6(r
-
r').
(8)
In the functions O'nm defined above the fluctuations of concentration within a single cell are extracted. By calculating the Fourier transform of both sides of Eqs. (7) we obtain the following set of differential equations describing the dynamics of a,,,~ in the reciprocal space:
() Of
o'x~(q, t) + 2
(0,) ~y
O'~y(q,t)-C
-2D(27rq)2(rxx(q, t) ,
d
(9a)
( +(0~)
/0y,)~. ,,.)
~ryy(q,t)+C_2D(2rcq)2c~,y(q,t), :Vs~Ys
(9b)
J. Gorecki, K. Kitahara/Physica A 231 (1996)588-607
592
dayY(q't)=-2
-~x ~.,y. ~r~y(q, t) - 2
() NOf
.~,y. a y y ( q , t ) - C
-2D(27rq)2 ayy (q, t) ,
(9c)
where
Crnm(q, t) = f exp (--27rirq)crnm(r , t)dr J and
C = 2(klx~ - k - l x ~ y s ) .
(10)
If C equals to zero then the stationary (in a long time limit) functions a,,,~ converge to zero. It may be shown that this happens if the detailed balance condition is satisfied for reactions (1). This feature justifies calling the appearing correlations nonequilibrium. The stationary solution of Eqs. (9) is well known [1] and it reads
~(q, ~) = ~(q, ~) = -~(q,
~) =
C 1 2D n2 + (2~rq) 2 '
(11)
which in the space variables corresponds to C 1 a x x ( r , oc) = cryy(r, oo) = --axu (r, oo) -- - 87rD [ r I exp (-to I r ' I).
(12)
The constant n for reactions (1) is defined as follows:
= 1(2klx~
- ( k _ l - k : ) ( c - 2x~) + 2 k _ 2 ( c - x~))
(13)
and it is positive because the steady state of system (l) is a stable one. The relationship c%~([ r -
r ' I,t) + 2a~u( I r - r ' h t) + (ruu(I r -- r ' I,t) = 0
(14)
which comes out of Eqs. (9) is hold for any time t provided that it is satisfied by the initial condition. It may be written as
~xs(~,
~',t) - x~6(r - ~') + 2 Z x r ( , ~ , ,.', t) + ~ y y ( ~ , ~', t) - u86(~ -
~')
which expresses the fact that if the chemical identities of particles are neglected then no spatial correlations are present in the system considered. In molecular dynamics simulations performed at liquid densities, the equilibrium correlations caused by a hard sphere potential have to be extracted from the total correlations between reactants, in order to have relationship (15) satisfied.
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
593
For reactants which are perfectly mixed the initial condition reads
a z x ( q , t = 0 ) = a u y ( q , t = 0) = azy(q,t = O) = O.
(16)
All three independent solutions of the set of linear differential equations (9) can be easily found, however, because of the form of the source term and the initial condition, only the solution for which
crxz(q , t) = cr~y(q, t) = -~rzy(q, t) = -(r(q, t)
(17)
is satisfied develops in time and the amplitudes of the other solutions are equal to zero. Using condition (17) one may reduce each of Eqs. (9) to the form
C - 2D(27rq)2o'(q, t).
(18) The solution of Eq. (18) reads o'zx(q) = O'uu(q)
=
-O'zu(q)
C 1 C exp(-2D(x 2 + (27rq)2)t) + 2D ~2 + (2~q) ~ 2D ~2 + (2~q) ~
=
(19)
The inverse Fourier transform of Eq. (19) can be done analytically (compare [9]) leading to the expression C 1 = ~yy(r) = -O~v(r ) _ -87rD i r i exp ( - ~ 17" I)
~(r)
oo
1
t)fdq
C
q sin(27rrq) exp(-8DtTr2q 2) + (27rq) 2
2
o
C 1 C 1 -- - 8 7 r D l r ] e x p ( - ~ ]r ])+ 87r----~r----~(sinh i (~ Iv 1)+ ½exp(-n] r 1)
2 2v5-
/
2
(20)
where x
o/x,_0
Thus Eq. (20) gives us complete information on the time evolution of the correlations of fluctuations.
594
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
3. Molecular dynamics simulations of nonequilibrium spatial correlations in a system with chemical reactions The molecular dynamics technique applied in this paper is very similar to the one used in our previous paper [5] to simulate the stationary correlations of fluctuations. The technique is based on the concept of "reactive" hard spheres (see for example [10]). According to it all the reagents are represented by identical hard spheres (masses and diameters are the same) and their motion is ruled by a standard dynamics with all the collisions regarded as elastic. The additional parameter, which has no influence on the mechanical motion, labels the chemical properties of spheres. We consider both reactive and nonreactive collisions between spheres. A chemical reaction may occur if particles representing appropriate reagents collide. The probability of a reactive collision is usually lower than 1. In order to describe a stochastic character of chemical processes a random number is called for each collision which may lead to a reaction and its value is compared with the steric factor for the appropriate process. If the value is not larger than the steric factor the collision is regarded as a reactive one and the chemical identity parameters of the spheres involved are modified. In simulations we assume that the activation energies for all reactions (1) are equal to zero. The technique may be easily adopted for any activation energy of the reactions provided that the computation is long enough to give a large number of reactive collisions. However, it is well known [11,12] that every thermally activated chemical reaction creates a nonequilibrium state of the velocity distribution of the reagents involved and this effect is not taken into account by the theory based on the master equation presented in the previous chapter. The kinetic energy is the only form of energy in a system of hard spheres. As all reactive collisions are elastic from a mechanical point of view no kinetic energy is released or consumed when a reaction occurs and the system as a whole remains in the thermal equilibrium with respect to the translational motion. Its temperature is related to the average energy per sphere. It is very important from the computational point of view to maintain such equilibrium because it allows us to extend the size of simulations using a prerecorded equilibrium trajectory. A trajectory, which was calculated for a system characterized by the periodic boundary conditions, may be used to enlarge the size of simulations [5,13]. The original small system may be periodically expanded in any of the directions by any integer number of the box lengths. Of course, if a chemical identity of molecules is neglected then such expansion does not bring us any new information. However, for a multicomponent chemical system, in which the translational motion is not related to chemical identity, the situation is different. First, a different chemical composition may be initialized in various boxes by marking the equivalent (by periodicity) spheres in a different way. Secondly, a steric factor (if it is not equal to unity) differentiates the time evolution in various boxes, because a collision between the equivalent (by periodicity) objects may be reactive in one box and nonreactive in another. Moreover, the periodic expansion simplifies calculations of spatial correlations, because it is enough to calculate interesting intermolecular distances between spheres within a single box. Of course the correlations whose range extends over a single box may be affected by
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
595
the periodicity of our system. In the following we show the results which have been obtained by periodic extension of the following systems: (I) A low density system. The original simulations were done for N = 500 hard spheres placed in a cubic box with the side length d = 14.7cr, a being the diameter of a sphere. The packing fraction is r / = 0.0824. This system was expanded by 10 box lengths in each direction so the total number of molecules is 500000. (II) A medium density system (N = 1000 spheres, d = 14a, r / = 0.191). This system was expanded by 10 box lengths in each direction leading to 106 spheres. (III) A high density system (N = 1000 spheres, d = 12.Sir, ~ = 0.268). This system was also expanded by 10 box lengths in each direction and the total number of spheres is 106. The average time between successive collisions of a particle is equal to 1.29 ps, 0.405ps and 0.223 ps in systems I, II and III, respectively. In the following we give all values of time in these units. In all these systems reactants X and Y are represented by spheres with the mass (rn = 32 a.u.) and diameter (a = 5.~). The temperature equals 300 K. The basic set of steric factors (probabilities of a reactive collisions) for processes (1) is the same as used in [5], 81 = 0.8,
8_ 1 = 0.2,
S2 = 0.1,
S-2 = 0.8,
but we have also performed simulations for slower reactions with steric factors, s t = 0.4,
J---1 = 0.1,
I
s 2 = 0.05,
sl_2 = 0.4,
and s I/1 = 0.16,
.~,11 _1
= 0.04,
s It2 = 0.02,
811 -2
= 0.16.
For all these sets of steric factors the concentration of X at its stationary state equals to • ,
=
7)c.
At the beginning of simulations the chemical identities of spheres are assigned in a random way, such that the concentrations of X and Y correspond to their stationary values zs, ys. Because of the random distribution of reactants we may expect that the no spatial correlations of concentrations are present in the initial state. The simulations performed for the stationary correlations [4,5] indicate that they have a limited range, which is of the order of the mean free path. In order to save computer time and cover the most important part of correlation functions we have introduced a cutoff R = 1.8cr in the range of intermolecular distances. The interval of interesting distances [cr, R] is divided into a number of subintervals of the same length A (we have used A _-- 0.02~r or 0.05tr). Moreover, we have introduced interval of times [0, to], during which the system evolution is observed. This interval is divided into subintervals of the length At each. Simulations give us the average number of spheres representing the reactant N whose distance from a sphere representing reactant M lies within [~r + iA, tr + (i + 1)A] at the condition that the time variable belongs to the subinterval [jAt, (j + 1)At]. Let us denote this quantity a s GNM (i, .j). If in simulations the time variable exceeds to then the parameters describing the chemical identities of spheres
J. Gorecki,K. Kitahara/PhysicaA 231 (1996)588~507
596
are mixed in a random way in order to destroy existing correlations and a new reaction path starts. The results presented below represent the average over a few hundred reaction paths. The correlations between particles g(r, t) for the system as a whole (i.e. with neglected chemical identity of spheres) can be obtained by scaling G(i, j) by the value, calculated assuming a uniform distribution of spheres, g([cr + iA, ~r + (i + 1)A], [jAr, (j + 1)At])
G(i,j) = 47r((~r + (i + 1)A) 3 - (a + iA)Z)cAt '
(21)
where G(i, j) denote the average number of spheres within the space interval [~r + iA, cr + (i + 1)A] and the time interval [jAt, (j + 1)At], without consideration of their chemical identity. Of course g should be the same for all j, but because of limited statistical samples some fluctuations are possible. Similarly we can calculate the partial correlation functions g x x , gxY and g y y ,
g x x @ r + iA, o" + (i + 1)A], [jAr, (j + 1)At]) Gxx(i,j) = 47r((o" + (i + 1)A) 3 - (o" + iA)3)x,A, '
(22)
g x y ( [ a + iA, o" + (i + 1)A], [jAr, (j + 1)At]) G x y (i, j) = %r((a + (i + 1)A) 3 - (a + iA)3)x~At '
(23)
gyy([~r + iA, ~7+ (i + 1)A], [jAr, (j + 1)At]) ayy(i,j) = ~Tr((~r + (i + 1)A) a - (~r + iA)Z)y~A~
(24)
If the nonequilibrium spatial correlations of fluctuations of reagents concentration are present then the partial correlation functions are distinctly different from the total correlation function. However the functions obtained from Eqs. (22)-(24) contain information on both correlations introduced by the hard potential of spheres and on correlations related to the chemical processes. In order to compare the nonequilibrium spatial correlations observed in simulations with the theory (compare condition (15)), one needs to separate the equilibrium correlations introduced by the hard sphere potential from a partial correlation function. Having in mind that the accuracy of gNM/g is higher than gNM itself because fluctuations in samples of interatomic distances cancel out, it is convenient to calculate the chemical contribution to correlations of fluctuations in concentration in the following way: cr,,m([ r - r '
I , t ) = n , m , a0(I ,"
(gNM(I r -- r' I,[JAt, (j + 1)At]) _ 1) I) \ g--~r---r;]-,[jAt,(j÷ 1)At]) /
=
go([ r - r' l) \
- ~ ; - _ - ~ ; I-, [jAr, (j + 1)At])
'
where go (r) is the radial distribution function for the system of spheres at the given packing fraction (it may be calculated separately with high accuracy) and t denotes a point within
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
597
the selected interval of times, for example t = ½(2J + 1)At. Our previous study [5] has shown that the nonequilibrium correlations between concentrations of reagents calculated on the basis of Eq. (22) have the functional form given by Eq. (12) and the parameters x and C are given by Eq. (13) and Eq. (10), respectively. The rate constants may be calculated from steric factors and the collision frequency. D may be obtained from the expression valid for a hard sphere fluid characterized by the packing fraction r/, 1 D - 16
~r ~go~¢ +)Vmr/-~RT'
(25)
where the value of the radial distribution function at the point of contact of spheres go ((~+) can be calculated from the formula 1 + 4r/go(o"+) = (1 + 77+ 7/2 + r/a)
(27)
(1 - ~/)a
Alternatively one may get D directly from simulation data by measuring the average square of the particle shift as a function of time. However the decay constant n calculated using Eq. (13) is usually smaller (by a factor of 2) than the observed rate of decay of correlations at short distances (and this is the region we focus our attention on). On the other hand the "theoretical" value of a gives a good estimation of the decay of ¢rnm at large distances. Here, being interested in correlations appearing at short distances, we use the numerical fit for g rather than the result of Eq. (13). Studying the stationary nonequilibrium correlations in system (1) we have also found [5] that the value of the constant C predicted by the theory is usually much smaller than the real strength of fluctuations. We believe that it is associated with the excluded volume effect, which is not taken into account by the simple theory presented in the previous section. From the practical point of view it means that C should be rather fitted to the data than calculated from Eq. (10). A typical time evolution of the functions ¢,~,~ obtained from simulations is presented in Figs. la-c. These functions are plotted using dashed line for four selected moments of time (the dash length is an increasing function of time) and compared with the stationary correlation function o'nm (solid line). It can be noticed that short distance correlations grow first and the long distance correlations appear later. As the theory says, a~y has an opposite sign to the other functions ( ¢ ~ , (Txy). The absolute values of all functions ¢r~ (Fig. la), ¢~u (Fig. lb) and ~ruu (Fig. lc) are almost the same for all moments of time, as predicted by Eq. (20). This feature indicates the accuracy of simulations. An evolution of rr~u calculated from Eq. (20) for three first moments of time, for which the data in Figs. la-c were shown, is presented in Fig. ld. The result for the longest time is almost the same as for the second longest so we have not shown it. The parameters required by Eq. (20) were the same as used in [5] to fit the stationary correlation function to the form given by Eq. (13) (~ and C) or calculated from simulations (D). The basic feature of observed correlations dynamics is correctly reflected by the theory: both methods indicate that the short distance correlations are growing first (and this seems natural, as they are caused by chemical reactions, which occur when spheres are in touch)
J. Goreck~ K. Kitahara/ Physica A 231 (1996)588-607
598
(a) }
.
.
. . . .
X
•
. . -
.
.
.
,
.
-0.010
1
I:D C
:/
X cT -0,030 (b O u9 C__ O
~:
:
1.1
1.2
f
j
--
--
-0,050
© (D L L
0 (9
-0.070 1 0
1.5
1.4
1.5
1.6
distonce /sphere diometer
(b) >~ 0 . 0 6 0 ~D ID X
0.040 r,b (f] O -'~ 0 . 0 2 0 © L
0 (D
0.000
10
1.1
1.2
1.3
1.4
1.5
1.6
distonce /sphere diorneter Fig. 1. The comparison of time evolution of (a) axz/c 2, (b) ~x~/c 2 and (c) o'~y/c 2 obtained in molecular dynamics simulations with a/c 2 calculated using Eq. (20). The simulations were performed for the system I with the fastest reaction. The correlations are shown for three moments o f time: 0.116 (the dotted line), 0.271 (the short dashed line), 0.426 (the medium dashed line), 0.814 (the long dashed line). The solid line represents the stationary correlation function. (d) shows a/c 2 calculated from F_.q. (20) for first three moments of time. The following values were used to evaluate F-4. (17): D = 0.625cr2/ps, n = 2.4~r - 1 and C = 1 5 D .
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
599
(c) ....... ~............ . .................
i
>- - 0 . 0 1 0 %)
i ............
C ©
:
L
i
i.Z[/:.::.i-:.:-:.:-~:-:-:
~ [ . . . .
'
I-
r--
~
--
--
I
,---
)C
-0.030 © ffD oO C 0
--0.050
C~
o
©
-0.070
~llll~l~ll~lllll~rlllllllFii
1.0
1.1
1.5
1.2
distance
1.4
/sphere
1.5
1.6
diameter
(d) >- 0.060 K?
C 0
)< C 0.040
L) "" \t C 0 ~-~ 0 . 0 2 0
" '"" ".,
"'"%
x
@ o
0 0000 .0
1.
1.
1.4
1.5
distance /sphere diameter Fig. 1. Continued.
1.6
J. Gorecki,K. Kitahara/PhysicaA 231(1996)588-607
600
and then, as time progress, they spread out. However the comparison of Fig. 1b with Fig. 1d shows that the theory underestimates the short distance part of o, whereas it overestimates the growth of correlations at large distances. In order to introduce more detailed description of correlations dynamics let us introduce the scaled correlation function in the form
=
cr,~,~(r, t = oo) "
(28)
It follows from Eqs. (13) and (20) that u does not depend on reactants N and M chosen. As the time progresses u(r, t) changes from 0 to 1 for every r. Referring to Fig. 1 we expect that u(r, t) approaches unity for small r first and at large r later. The dynamics of u(r, t) is governed by the following equation:
1 0 u ( r , t ) _ 1 0 ( ~_~ ) (1) O Ot r 2 Or r2 u(r,t) + 2 x + +47rr exp(t~r)(1 -
u(r, t)),
~___~
,(r,t) (29)
which results directly from Eq. (7) and the form of the stationary correlations (Eq. (12)). As expected u(r, t) = 1 is the stationary state of Eq. (29). Using Eq. (29) seems to be more convenient than Eq. (20) because u(r, t) does not depend on the amplitude C as it scales out. u(r, t) as a function of distance for a few selected moments of time is presented in Figs. 26. Fig. 2a shows the growth of correlations for a system characterized by a low density (I) and fast reactions (si). Comparing the results of simulations (Fig. 2a) with theory (Fig. 2b) one may notice that for the first of moments of time considered (t = 0.116, the dotted line) our theory underestimates correlations at short distances, but overestimates them at large ones. The theory predicts a rapid growth of correlations in the interval t[0,116, 0.426], much faster than measured from simulations. According to the theory at t = 0.426 the nonequilibrium correlations of reactants are fully developed, whereas simulations show that for distances > 1.5~r they reached only half of their stationary values. Figs. 3 and 4 present results for systems at a medium density (II). If the reaction is fast (si) the differences between simulations (Fig. 3a) and Eq. (20) are enormous. For time t = 0.242 (solid line) the correlations are well developed according to simulations, but theory says that they are hardly visible. Their evolution observed in simulations remains an expanding front, whereas Eq. (20) predicts an uniform growth in the whole space. The results obtained for five times slower reactions (Figs. 4a and 4b; the set s~~) are in better agreement, but also in this case the rate of growth of correlations is underestimated by theory. The same kind of disagreement between simulations and theory may be observed in results for the most dense system (HI). The data for a fast reaction (si) are compared in Figs. 5a,b, for the slower one (s~) in Figs. 6a,b.
J. Gorecki, K. Kitahara/Physica A 231 (1996)588-607
601
(a) >~:~ 1.0 ©
L__
q--
q
x
'
- - -
L
_
_
@8
I
1 I
© q)
I --
x~06
I I
/
If)
t
_
0
O0.4 q) L
o 0.2 q)
o (p O@ 10
1.2
11
t.5
1.4
/.5
1.6
distonce /sphere di@meter
(b) 1.0
(f) C 0
0.8
© 0.6 L
© 0
g 0.4 o
",,
L)
.....
02
00
r
1.0
,
~
i
I
1.1
i
~
i
i
1.
2/
P
i
,
~
I
15
i
~
i~
I
1.4
i
f--i
i
15 q
1.
t
i
,
1 6
distonce /sphere diorneter Fig. 2. Comparison of (a) molecular dynamics simulations for the system 1 with rates si with (b) the theory. (a) The results of simulations are shown for four moments of time: 0.116 (the dotted line), 0.271 (the short dashed line), 0.426 (the long dashed line) and 0.814 (the solid line); (b) calculations are for O. 116 (the dotted line), 0.271 (the short dashed line) and 0.426 (the long dashed line). The values of parameters used are D = 0.625a2/ps and ~ = 2.40 "-1.
J. Gorecki, K. Kitahara/Physica A 231 (1996)588-607
602
(a) 1.0 C 0
08
b~ C 0
o0.4 -g k_
o 02 72)
©
@@
lo
1.1 1.2 15 distance /sphere
1.4 1.s diameter
(b} 1.0
tO
0.8
C 0 o
L_ L
0.6
o
"'.
(D
"'"
0.4 o
o
O3
0.2
" 4
0.0
l l ; J l l l l l l l
10
1.1
1.2
,1. 13
, ,,
1.
14,,,
~ I 1.5
,
, ,
, 1.6
distance /sphere diameter Fig. 3. The time evolution of the scaled correlation function v(r, t); comparison of (a) molecular dynamics simulations for the system II with rates s i with (b) the theory. (a) The results of simulations are shown for four moments of time: 0.054 (the dotted line), 0.104 (the dashed line) and 0.242 (the solid line); (b) calculations are for 0.242 (the solid line), 0.494 (the short dashed line) and 0.741 (the long dashed line). The values of parameters used are D = 0.229a2/ps and i¢ ---- 4 . 8 6 a -1 .
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
603
(a) x 1.0 C
x 08 C ©
© ~0.6 ...... : i ';
U3
i
©
-1 -i :
L i
4~
(30.4
~L~
LL
t _
L k~
i
© (.D
L_
I
'1 i_ I_
0.2
i -%
©
~ -!
-I
L ....
I_
. . . . . . . . . . .
0.0
10
. . . .
11~
. . . .
~
.... t
'
'
.=S ~ -2
-
_-
-
-
I .i2.... 1. 13 .... 1.4I ....
r q
1.
-
--
15 . . . . 1.6
distance /sphere diameter
(b) 1.0
0.8 (/3
o © ~0.6 L
o (.3 ~0.4 U
0.2
0.0
i-(-?-(?-£rr
.0
1.1
I
1.2
~FIilll~lllEillliS
1.3
1.4
1.5
1.6
distance /sphere diameter Fig. 4. The time evolution of the scaled correlation function u(r, t); comparison of (a) molecular dynamics simulations for the system 11 with rates s~' with (b) the theory. (a) The results of simulations are shown for four moments of time: 0.074 (the dotted line), 0.272 (the short dashed line), 0.519 (the long dashed line) and 1.21 (the solid line); (b) calculations are for 0.272 (the short dashed line), 0.519 (the long dashed line) and 1.21 (the solid line). The values of parameters used are D : 0.229a2/ps and n = 2.35(7 -1 .
604
7. Gorecki, K. Kitahara/Physica A 231 (1996)588-607
(a) )%) 1,0 C ©
X C © ©
0.8 t
i_~1 I
~06
'
(f)
1 I
--i i i
•
O
i
©0,4
L
L
(2) ~D
i
O2 U q) (D q)
,_
I
0@
1•0
11
12
distance
/sphere
114
I I5
diameter
(b) 1.0
(/) 0 8
-.--
0
'~"~)(~)0.6
"-""""
'"-.
~0.4
00
i
0
p ~ F
J
i
F
I
i
J
i
i
i
ITt
I
I ~
1.~ ~2 distance / s p h e r e
I
J
l
I
I
JTI
l
,
lq
1.s dTorneter
~
T T r
1.4
Fig. 5. The time evolution of the scaled correlation function t,(r, t); comparison of (a) molecular dynamics simulations for the system III with rates sl with (b) the theory. (a) The results of simulations are shown for four moments of time: 0.247 (the dotted line), 0.471 (the short dashed line), 0.919 (the long dashed line) and 1.10 (the solid line); (b) calculations are for 0.919 (the long dashed line), 1.10 (the solid line) and 1.35 (the short dashed line). The values of parameters used are D = 0.131~ 2 / p s and n = 7 . 9 a - 1 .
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
605
(a) x 1.0 C
X C © (L>
0.8
i-; I i
ii L], :
W)
i
,
t ~
C 0 i
©0.4
I- I
&) L_ m
o
i
,__
L
~
k_ 1
0.2 -o © © 0 (/3
"
, ~.. , 0.0
,
0
i
i
,
1 I1
i
]
I
7.].'~-, ............ q=..-~,-, T
i
I
12
~
J
i
i
L~ ] ~ ........~ - : ' =~i- -V..~- ~ ..,+- ]
I
1.,3
i
i
,
,
I
r - i-'-i 7]
- -
1.4
"'""ST"]"
~~
i
-6
distance /sphere diameter
(b) 1.0
u) C 0
0.8
© ~0.6 L
o t)
~0.4 I
©
k) u)
02
O0 10
1,1
1.2
1 5
1.4
1.5
1.6
distance /sphere diameter Fig. 6. The time evolution of the scaled correlation function v(r, t); comparison of (a) molecular dynamics simulations for the system III with rates s~ with (b) the theory. (a) The results of simulations are shown for four moments of time: 0.247 (the dotted line), 0.471 (the short dashed line), 0.919 (the long dashed line) and 1.10 (the solid line); (b) calculations are for 0.919 (the long dashed line), 1.10 (the solid line) and 1.35 (the dotted line). The values of parameters used are D = 0.131a2/ps and ~ = 5.6o "-1 .
606
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588-607
4. Conclusions In this paper we have presented results of molecular dynamics simulations of the growth of nonequilibrium spatial correlations of reactants fluctuations in a system with chemical reactions. The considered systems exhibited such correlations because the detailed balance condition was not satisfied. Our simulations started from well stirred reactants, thus the correlations were restored during the progress of reactions. We observed that the correlations of fluctuations appear at short distances first and then spread out in space. We considered correlations at distances not exceeding a single molecular diameter, because in this range they are the most important. In this range of distances the growth of correlations is very fast and it is shorter than the average time for successive collisions of the same particle. These conclusions agree with observations of Prigogine et al. [6] who found that the radial distribution function (RDF) of a liquid is first reestablished at short distances. This process is also very fast and the first peak of RDF appears in less than a single collision time. We have also derived a simple theoretical description of the growth of nonequilibrium spatial correlations in a system with chemical reactions, based on the master equation for reaction-diffusion processes. It has already been shown that an analogous description of the stationary (independent on time) correlations is in a good agreement with molecular dynamics simulations [5] and with simulations based on the Bird technique [4]. This result is quite surprising because on one hand the theory is based on the large volume expansion, but on the other hand the range of correlations is of the order of the mean free path, which, in some of the systems studied [5], may be a fraction of the sphere diameter. Simulations show that this simple theory gives the qualitative description of correlations growth only. Both theory and simulations indicate that the correlation function approaches its stationary form at short distances first and then it spreads o u t - correlations at large distances approach their stationary values later. However, the quantitative description of the correlation dynamics is different. Molecular dynamics suggests front-like expansion of the stationary correlations, whereas the theory supports more uniform growth at all distances. We believe that the discrepancy comes mainly from the fact that the diffusive term in the master equation is replaced by the Laplace operator. It would be too naive to expect that a classical diffusive term may correctly describe processes which occur at times for which the motion of particles is dominated by a ballistic flow. In the case of chemical wave front propagation, generated by a fast reaction, it may be shown that the reaction-diffusion equation need to be associated with the equation describing the relaxation of a diffusive flow in order to get a correct wave speed [14]. The reactions considered in this paper were very fast and the diffusive flow relaxation seems also important.
Acknowledgements This work was supported by the grant KBN 2 P303 018 05 provided by the Polish State Committee for Scientific Research and by Grant-in-Aid in Priority Areas 07240106 from
J. Gorecki, K. Kitahara/ Physica A 231 (1996)588--607
607
Ministry of Education, Culture and Science of Japan. One of the authors (J.G.) is grateful for JSPS support during his visit to Tokyo Institute of Technology. The authors would like to thank Professor U.M. Titulaer for helpful discussions. References [1] Y. Kuramoto, Prog. Theor. Phys. 52 (1974)711; G. Nicolis and M. Malek Mansour, Phys. Rev. A 29 (1984) 2845. [2] K. Kitahara, K. Seki and S. Suzuki, J. Phys. Soc. Japan 59 (1990) 2309. [3] K. Kitahara, K. Seki and S. Suzuki, in: Proceedings of Molecular Dynamical Processes, A. Blumen et al., eds. (World Scientific, Singapore, 1990)p. 403. [4] G. Nicolis, A. Amellal, G. Dupont and M. Mareschal, Journal of Molecular Liquids 41 (1989) 5. [5] J. Gorecki, K. Kitahara, K. Yoshikawa and I. Hanazaki, Physica A 211 (1994) 327. [6] I. Prigogine, E. Kestemont and M. Mareschal, in: Microscopic Simulations of Complex Flows, NATO ASI series, Vol. 236, M. Mareschal, ed. (Plenum Press, New York, 1990). [7] J.S. Turner, J. Phys. Chem. 81 (1977)2379. [8] C.W. Gardiner, Handbook of Stochastic Methods (Springer, Berlin, 1983). [9] I.S. Gradshteyn and I.M.Ryzhik, formula 3.954.1, Table of Integrals, Series and Products (Academic Press, New York, 1992). [10] J. Gorecki and J. Gryko, Computer Phys. Comm. 54 (1989) 245. [11] I. Prigogine and E. Xhrouet, Physica 15, (1949) 913; for the recent see J. Popielawski, in: The Dynamics of Systems with Chemical Reactions, J. Popielawski, ed. (World Scientific, Singapore, 1989), and references therein. [12] J. Gorecki and B.C. Eu, J. Chem. Phys. 97 (1992) 6695. [13] J. Gorecki, in: Far from Equilibrium Dynamics of Chemical Systems, J. Gorecki et al., eds. (World Scientific, Singapore, 1994); Molecular Physics Reports 10 (1995) 48. [14] J. Gorecki and J. N. Gorecka, Acta Physica Polonica B 26 (1995) 1177.