Reactive dynamics in a multispecies Raymond

lattice-gas automaton

Kapral

Chemical Physics Theory Group, Department of Chemistry, University of Toronto, Toronto, Ontario MSS IA1 Canada

Anna Lawniczak Department of Mathematics and Statistics, Univeristy of Guelph, Guelph, Ontario NIG 2 WI Canada and CNLS, Los Alamos National Laboratory, Los Alamos, New Mexico 87545

Paul Masiar Chemical Physics Theory Group, Department of Chemistry, University of Toronto, Toronto Ontario MSS 1Al Canada and Ontario Centre for Large Scale Computation, University of Toronto, Toronto, Ontario MSS IA1 Canada

(Received 17 June 199 1; accepted 28 October 199 1) A multispecies reactive lattice-gas automaton model is constructed and used to study chemical oscillations and pattern formation processes in a spatially distributed two-dimensional medium. Both steady state and oscillatory dynamics are explored. Nonequilibrium spatial structures are also investigated. The automaton simulations show the formation of rings of chemical excitation, spiral waves, and Turing patterns. Since the automaton model treats the dynamics at a mesoscopic level, fluctuations are included and nonequilibrium spatial structures can be investigated at a deeper level than reaction-diffusion equation descriptions.

I. INTRODUCTION The full microscopic description of complex chemically reacting systems is a challenging problem that is difficult to solve by direct simulation of the equations of motion. This is especially true of multicomponent reacting systems constrained to lie far from equilibrium where many interesting phenomena occur on long distance and time scales. It is well known that in such circumstances interesting spatial and temporal structures can arise as a result of the interplay between the nonlinear reaction kinetics and the diffusion processesin the system.’ Examples include bistability, temporal oscillations, chemical waves, and Turing patterns. While these phenomena are usually described by reaction-diffusion equations, their existence is affected by molecular fluctuations in the system, especially near bifurcation points.’ If the goal is to attempt to obtain a molecular picture of these rather unusual spatiotemporal structures, it is desirable to seek a route other than direct molecular dynamics (MD) simulations of the equations of motion. There have been studies of simple models of far-from-equilibrium systems using full molecular dynamics,3 but such investigations have been necessarily limited in the number of particles and space and time scales which could be explored. The limitations of conventional MD techniques for the study of such systems has prompted our interest in alternative dynamical schemes. We present in this paper a lattice-gas cellular automaton (LGCA) description of multicomponent spatially extended reactive systems. While an outline of the method and illustrations of its utility were given in Ref. 4, we focus here on the full mathematical description and present a variety of simulation results to illustrate some directions for applications. Lattice-gas cellular automaton models’ have been applied with considerable success to the study of fluid flow problems in hydrodynamics. 6.7 The idea behind such mod2762

J. Chem. Phys. 96 (4), 15 February 1992

els is the construction of an artificial dynamics that captures many of the essential features of the real collision processes in reactive media, yet is simple enough to permit simulations on large systems composed of many particles. The reactive LGCA scheme described below satisfies these criteria and simulations of the dynamics of millions of particles for long times is a simple task. The automaton model occupies a position between full MD and macroscopic reactiondiffusion equation descriptions of the dynamics; it is a mesoscopic model of the dynamics. The LGCA is an alternative to other approximate MD methods, for instance, the Bird method’ which has been applied’ to the study of some of the phenomena described in this paper. It also permits a deeper level of description than other discrete methods such as “classical” cellular automata and coupled map lattices.” Finally, we note that a multispecies LGCA model has been constructed by Dab, Boon, and Li” and used to study spatiotemporal structures in a two-variable reaction model; their work complements this study. The outline of the paper is as follows: Section II begins with a brief qualitative description of the LGCA, followed by a mathematical formulation of the automaton dynamics. Sections III and IV present the derivation of discrete Boltzmann equations, continuity equations, and rate equations that follow from the LGCA equations of motion. In these sections the connection between given reaction-diffusion equations and specific LGCA models is established. This allows the construction of a LGCA dynamics that yields in the macroscopic limit a particular nonlinear chemical rate law. A two-variable reaction scheme, the Selkov model, is studied in Sec. V and an automaton corresponding to this reaction is constructed. The results of simulations of the Selkov LGCA are presented in Sec. VI. There we show that the automaton can produce chemical oscillations, chemical rings, and spiral waves characteristic of excitable media, and

0021-9606l92lO42762-15$06.00

0 1992 American Institute of Physics

Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

2763

Kapral, Lawniczak, and Masiar: Reactive dynamics in a lattice gas

Turing bifurcations. All of these phenomena are found only in multispecies reacting systems and have motivated the extension presented here. The conclusions are given in Sec. VII. II. LATTICE-GAS

CELLULAR AUTOMATON

MODEL

A. Overview The mathematical description of the multispecies reactive LGCA is rather technical in character and requires a somewhat cumbersome notation in order to insure generality. For this reason it is useful to give an overview of the principal features involved in its construction. We consider spatially extended systems composed of a number of species X,, 7 = l,...,n, which may undergo chemical reactions ofthe type a,X,

+ --- + a,X, -8,X,

+ -** +A&.

(2.1)

Basically, the problem is to construct an approximate dynamics for both the elastic and reactive collision processes that occur in the real system. To simplify the dynamical description space is made discrete by placing the system on some regular lattice with coordination number m; time is also taken to be a discrete variable. To complete the discrete phase space description, a discrete velocity which can point along any of the m lattice directions is associated with each particle. In order to have a computationally tractable scheme an exclusion principle is assumed where a given node of the lattice cannot be occupied by two particles of the same species with the same velocity.12 The collision processes are taken to be short ranged or local and occur at the nodes of the lattice. In the form of the model constructed here, only the dynamics of certain reactive species are followed and the solvent and possibly other species, which may be in excess or constrained, are taken into account through their effects on the dynamics of the species of interest. (This restriction is not an essential part of the model.) More specifically, elastic collisions of the X, species with the solvent (and each other) are modeled by random rotations of the velocities of the X,. Reactive collisions change the numbers and velocities of the particles of each species and in the automaton such collisions are effected by combining two operations: local reactions (2.1) at a node that occur with preassigned probabilites which are independent of the velocities, and random rotations. The net effect is a reactive event where both particle numbers and velocities are changed, just as in real reactive collisions. Finally, between such collisions the particles propagate freely and move from one lattice site to another. In more pictoral terms one way to visualize this multispecies LGCA model is as a “stack” of n “species lattices” with identical labeling of the nodes. Each species resides on its own lattice and the propagation and rotation steps are taken independently on each lattice but the dynamics on these lattices is coupled by the chemical reactions. While not a true molecular dynamics, the model captures the essential features of the full reactive system since the reactive and elastic collision processes mimic those in the real system. We now give a precise formulation of the rules that govern the dynamics of the system.

8. Mathematical

formulation

The goal of this subsection is the derivation of microscopic equations of motion for the system. These microscopic equations constitute the starting point for the derivation of chemical rate laws and reaction-diffusion equations. T’he links established between the microscopic automaton dynamics and the macroscopic equations provide a means for constructing automaton models that yield a given rate law or reaction-diffusion equation. The LGCA model for the chemically reacting system (2.1) can be described in formal terms: particles of species X, ,...,X,, move on a lattice 2 with coordination number m and periodic boundary conditions. Each node of the lattice is labeled by the discrete vector r and each particle has associated with it a discrete velocity ci, i = 1,...,m pointing in one of the m possible directions on the lattice. We assume that the ci are unit velocity vectors connecting the node to its nearest neighbors, where i increases counterclockwise and 1 corresponds to the positive direction of the x axis. The exclusion principle prevents two particles of the same species to be at the same node with the same velocity. Consequently, there may be 2” different configurations of particles at a node with the same velocity and 2”” different configurations of particles at a node. The dynamics of the reactive automaton can be neatly described by a set of microdynamical equations that are the analogs of Newton’s equations of motion for this system.5 The microdynamical equations are evolution equations for a set of Boolean random variables that describe the species type and discrete phase space coordinates of the particles. Let S, be a set of all m X n matrices s = [sir 1, (i = l,..., m, T = l,..., n) such that S=

s = [si7]:siT = 0 or 1, 1

Vi = l,..., m, O( i

s,(n

r=l

(2.2)

, I

andlet Si, l
(2.4)

The evolution of the system occurs at discrete time steps and the lattice gas automaton updating rule can be defined on the Boolean matrix field (2.3). At each node at each time step the updating rule consists of deterministic propagation followed by collisions, which may be either elastic or reactive. As in the one variable LGCA13 the evolution consists of a

J. Chem. Phys., Vol. 96, No. 4,15 February 1992 Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

2764

Kapral, Lawniczak, and Masiar: Reactive dynamics in a lattice gas

product ROCOP of three elementary operators: propagation P, chemical transformation C, and rotation R. In the propagation step each particle moves in the direction of its velocity from its node to a nearest-neighbor node. The velocity of each particle is conserved during this operation. Letting $denote the configuration in the ith direction after propagation we have P:$(k,r)

= qi(k - 1,r - q).

(2.5) The rotation operation can be described in the following way. At each node, independently of the others, the particle configuration on each species lattice is rotated through one of the possible angles 19~= 12rr/m, I = l,...,m with probabilityp,; the probabilities may be different for each species and they must be selected to preserve isotropy.i4 To explicitly construct the rotation operator we let {~lr(k,r):rd’, k = 1,2;. *}, be independent sequences of identically distributed, independent, Bernoulli-type random variables. Using these variables we may then construct new random variables (2.6a) and

7%

(2.6b)

with probabilities Pr(YE = 1) =Plr, (2.7) In Eqs. (2.6) and (2.7) we have dropped the dependence of 6 and y on k and r. We adopt this convention in the sequel when confusion is unlikely to occur. If 7: denotes a configuration of species X, after a rotation in the ith direction at a node then

R:rl: = $ tin, ,.r, where the addition in the index of 7 is modulo m. Since

+ 2 r"TC -vir I= 1

svir + AE(qT),

= I-

n

m (2.12) v 3ov “L,,)) I-d 0 r=l where 4, = fl, - ar. It is clear that if q, > 0 then particles of species X, are created while if qr < 0 then particles of species X, are destroyed. The chemical transformation step can be described formally in terms of the microdynamical variables vir. Let 55’ denote the set of all possible chemical transformations a -0, I.e., Ce = {(a,p):Vr

= l,..., n, O
O@,Gm}. (2.13)

For each input stoichiometric vector a let % (a) be the set of all possible output stoichiometric vectorsp different from a, I.e., (2.14)

5f (a) = @:(a,B)ECe,a#lS.

The reaction events in any realization of the evolution process are determined by combinations of random variables c gs such that

we may also write Eq. (2.8) in the form VP, =rlir

(2.11) 1 P(aP). B#U The off-diagonal elements of P with a7 >p, correspond to reactive transitions for which particles of species X, are destroyed while if ar
nap =

m-l ““-7=*-,g1

on the different species lattices and forms a central part of the reactive LGCA model. In this step the reactions (2.1) that the chemical species undergo are modeled by probabilistic rules for the increase or decrease of the number of particles. The numbers of the different species in the reactant and states can be specified by the vectors product a = (a, ,..., a,, ) and p = (fl, ,..., p, >, respectively, where O
@yp(k,r):i=

+vi+,,r)

(2.10)

which defines the rotation operator AiR,(q, ), that takes on the values { - l,O,l}. The propagation and rotation operations are intended to describe the free streaming and elastic collisions that occur in the real system; their net effect is to produce diffusive motions of the various species. Since the rotations are carried out independently on each species lattice and since the rotation probabilities need not be the same for all species, the different species may have different diffusion coefficients.14 If only propagation and rotation operations are considered the species lattices are uncoupled and the dynamics of each species is independent of that of other species. The chemical transformation C couples the dynamics

l,..., m, rdf,k

= 1,2;**}

(a$)& (2.15)

are independent sequences, independent of sequences of random variables defined in Eq. (2.6)) of identically distributed independent Bernoulli-type random variables. Denoting by ikc{l ,...,m} the velocity indices we define random variables PCi , ,..., iq ), for distinct sets of velocity indices {i, ,..., i,}, which govern the chemical transformations a -+fi r”p(i I,...,iq)

= fi j=l

g;@

J-J (1 ---g$ k # L...,q

x_Jr& j), (I - fcp”), (2-*6) where the prime on the product indicates that only distinct (distinct from each other and the i, ,...,i, > ik values are con-

J. Chem. Phys., Vol. 96, No. 4,15 February 1992

Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

Kapral, Lawniczak, and Masiar: Reactive dynamics in a lattice gas

sidered. The probability random variables are Pr[y@(i

distributions

associated with these

I,..., i,) = 1] =2,

velocity indices corresponding to velocity directions occupied by particles of X,

(2.17a)

Ua,d

,..., jq) =I)]

(2.18)

= Ci, ,..., ia,),

and by r(a) we denote one of the possible configurations of directions occupied by species X, for an initial stoichiometric vector a

and obviously = * -9.

(2.17b) “43 In terms of the random variables fP( i, ,...,i, ) we may write expressions for random variables r”,B (r,i) which describe all possibilities in the reaction a +p, where particles of species X, in velocity state i are produced ( + ) or destroyed ( - ), with particles of other species undergoing number changes specified by the stoichiometric coefficients of the reaction. The derivation of the explicit forms of these random variables is given in Appendix A. For a given input stoichiometric vector a and species X, we use r(a,r) to denote a subset of {l,...,m} of distinct Pr[ VT4

2765

r(a)

= {r(a,7):7=

l,..., n}.

(2.19) It is clear that the input vector a may be realized by many configurations of particles, i.e., by many r(a). We introduce the following products of Boolean fields labeled by an input stoichiometric vector a and a set I’(a) of occupied velocity directions: Q?q;r(a)]

= fi n VjT ,R,, (1 -vkr). (2.20) ?-= 1&I-(cz,r) If 7: denotes a configuration for a species r in the ith direction after a chemical transformation then

I

rl:= c mar>1

C Qah;rWl1(1-?iiT) 2 Y! (r,i)+ vi71- C yt!(74

r(a):l-(a.r)3i

t-c

@%

C

a:a,
Qahr(dl

[

- (a,r)

pE'I:

II

- (cr,T)

(2.21)

(1-“i),,~=,~~(~)i)+~i~[l-as~~~~~~(7,i)]], 1

I where

772 = vi7 -

% + (a,71 = CB:(a,PW,q, and

X

% - (a,71 = Cs:(a$k%,qr < 01, denote the sets of all possible output states p in the reaction a+fl in which particles of species X, are produced or destroyed, respectively. From Eq. (2.2 1) we can give a physical interpretation of the processes contributing to ~z. The two terms involving sums on a correspond to processes for which the number of particles of a species X, at a node decrease (expressed by the condition a, > 1) or increase (expressed by the condition a, cm). The first term in each curly bracket arises from reactive events while the second term corresponds to nonreactive events. The Q “[q;I?(a) 1 factor insures that the configuration at a node corresponds to an input stoichiometric vector a realized by particles located at r(a) velocity directions. For a such that ar ) 1 and IY(a,r) 3i Qa[q;r(a)]vi,

(2.22)

= Q”[q;r(a)l,

and for a such that ar .Cm and JY(a,r) 3i

Qah;rWl~iT

C

C

a:a,>i r(a):r(a.+)31

C y?! (r,i) px - ((I.T)

+ cz:a’ 7
3i



“[q’r(a)

xBvqaT) ff (rd.

I

(2.25)

The last two terms correspond to processes where the particle number of a species X, decreases or increases at a node, respectively. Formula (2.24) follows from Eq. (2.20) in Ref. 13 and the fact for any species u (2.26) where Iand Jare disjoint subsets of indices of { l,...,m) such that IU J = { l,...,m}. Finally we define a chemical transformation operator for a species X, , A$ ( q;afi) , by comparison of ?7C

=

Q%;r(a) 1

(2.23)

= 0.

Since Vi7

C

c

a:~~,>1t-(d:r(a,e3i

> 01,

QVWa)

formula (2.2 1) can be simplified to

I,

=

lir

+

C

A:

(q;ak%

(2.27)

(a,mw,a#B

(2.24) with Eq. (2.25). This operator has the property that for (aBW, a#B

J. Chem. Phys., Vol. 96,toNo. February 1992 Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject AIP4,15 license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

2766

Kapral, Lawniczak, and Masiar: Reactive dynamics in a lattice gas

Aic,(q;aB) =

0

if a, =fi,

-1 or 0 1 or 0

ifa,>p,.

the density of a species X, at time k (2.28)

if ar-cPT The cellular automaton rule consists of successive applications of the transformations described above, i.e., propagation followed by collision, RoCoP. The microdynamical equations for the automaton follow directly from the microdynamical equations for these transformations. From Eqs. (2.8), (2. lo), and (2.27) the result of the sequence of transformations RoC is

+T

I=

1 (n,Pk

&

.a#B

YX+ I.r(rl;aP)

=qir +AE(q) +AZ(q),

(2.29)

where we have defined

p,(k)

=N-’

(2.35)

C p,(k,rh r‘c?”

where Nis the number of nodes, and the local mean velocity of a species X,

u,(b)

=

[p,(k,r)] 0,

-‘X”=,ciNiT(k,r),

if p,(k,r)#O otherwise. (2.36)

Of course Ni,, pT, and u, depend on p, which we drop to simplify the notation. III. BOLTZMANN

EQUATION

By taking the expectation value of the LGCA microdynamical equation (2.3 1) , we can derive a discrete Boltzmann equation. Since the random variables lIr (k,r), 6 yD( k,r) are independent and also are independent of the evolution of the cellular automaton up to time k - 1, this implies that

Ah)=,z,caa,~=~s~A~+,,,(q;a~).

(2.30) 9 7 Using this result the microdynamical equations for the cellular automaton for a species X, take the form

(2.31) vi7 (kr) = 7; + A: (qp) + AZ ($‘I. This equation is a compact representation of the full microscopic dynamics of the automaton and can serve as a basis .for further theoretical analyses. The automaton rule can be cast in a somewhat more general form that is useful in simulations where the magnitudes of the diffusion coefficients need to be changed by large amounts. This is easily accomplished by performing different numbers of propagation and rotation steps on the various species lattices. The operator product PoRoCcan be generalized to fi 0*4 (2.32) U’,oR, )oC, r=l where the operator 0 with a + I superscript means I applications of the operator in its argument for every application of C, while a - 1 superscript means application of the operator in its argument every 1 th application of C. The evolution of the lattice gas can also be described in terms of a probability distribution P, [ k,s( * ) j at a time k with the initial distribution ,u. Since the random variables glr, 5 “governing the random rotations and chemical transformations, respectively, are independent, the entire Booleanfield{q(k); k = 0,1,2,...}definedinEq. (2.3) isaMarkov process. The Chapmann-Kolmogorov equation for Pp can be found in Ref. 13. The introduction of the probability distribution P,, and its corresponding expectation E,, enables one to define physically interesting average quantities; the mean population of a species X, in the ith direction N,,(b) =E, [vir(kr)], the local density of a species X,., p,(kr)

= 2 i=

1

Ni,(k,r),

(2.33)

NiT(k,r)

=Ng

+ Cc(N’)

+ 2 Pl7 I=

1

c

Ep A:+ 1.r(qP;crB).

In the above NE = Ni, (k - 1,r - ci) and for any matrix w=

[wir],

‘Ctw)

= 2 P[rwi+/,r - wirj I= 1

where addition is modulo m. If for all a, P(aa) = 1, i.e., no chemical transition takes place, each particle performs a random walk on its species lattice. These random walks are not independent since when two or more particles are simultaneously present at a node, they are rotated by the same angle in order to preserve the exclusion principle. Except for this weak dependence between trajectories due to the exclusion principle, all elastic collisions are independent and a kinetic Boltzmann equation for the average local number of particles can be derived. In the presence of reactive transformations collisions can no longer be assumed to be independent. For example, particles which have reacted to produce a new particle can collide reactively again with this new particle on nearby nodes. However, if the chemical transformation probabilities P( ap), a #/3 are sufficiently small, reactive transitions are rare events and we can again assume independence between collisions. In this approximation, we can factorize the expected value in Eq. (3.1), and the expected occupancy Ni, (k,r) satisfies the lattice Boltzmann equation Ni,(k,r)

=NR

+ CE(Np)

+ CiC,(N’),

(3.3)

where C,C,(N) = 5 pr+ C CF+,,,(N;a/?). I= 1 (w%4C,a#B

(2.34)

(3.1)

(a,8kv,o#8

(3.4)

Intheabove

J. Chem. Phys., Vol. 96, No. 4,15 February 1992 Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

2767

Kapral, Lawniczak, and Masiar: Reactive dynamics in a lattice gas

c C?+,.r(WB (Odvw~.a#B = ,,gmr(n):$&T)si Q”‘N’r(a)lbar (3.5) - ,:?>I N’r(a)ldarf r r(o):~,r)~iQa’

PT(O NiT(t,x) =-t

m equation (4.3 ) reduces to

the continuity dp, (t) -= dt

(4.4)

C (P,-a,VYaB (or.P)~%“,a#P

‘jj, ,zx(~IYLa~J

where b,, =

*.*,$o”m zP(afi), 1,” n 7

2 **. 2 &=++l B, = 0

(3.6)

X ( - l)‘-“‘m-[p,(t)‘.

and

Details concerning the derivation of Eqs. (3.6) and (3.7) are given in Appendix B. Since factorization in Eq. (3.3) does not hold exactly unless the system is in equilibrium and in a state strictly satisfying the Boltzmann hypothesis, Ni, as computed from these equations will not give the true expected occupancies of the process. Nevertheless, we suppose that under some conditions on the initial state analogous to Eq. (3.7) in Ref. 13, the automaton inherits the limiting behavior of Eq. (3.3) in some space-time scaling regimes. In the kinetic regime (3.3) converges to the continuous space-time Boltzmann equation JNi, (6x1 + Ci~Vlv,, (4x1 dt = CIR(N) + C~(N) = C,,(N).

(3.8)

Details concerning the derivation of this equation are omitted since they are similar to those leading to Eq. (3.8) in Ref. 13. IV. RATE LAW The results of Sec. III form the starting point for the derivation of the continuity equation and chemical rate law for the total densityp, of species X,. Summation of Eq. (3.8) over the velocity index i gives the continuity equation for local particle density. Since for each Z, I= l,...,m and (a$W, a#B 2 Cf+ ,,rW;aP) i= 1

= ,z, Cg(N;a;B),

The derivation of Eq. (4.5) is given in Appendix C. Note that the highest power of the density of a given species is limited by the coordination number of the lattice. If we identify coefficients of equal powers of the densities in the above equation with those in the phenomenological rate law we can establish relations between the chemical transformation probabilities and the rate coefficients in the chemical rate law. Using these considerations a set of LGCA rules consistent with a given macroscopic rate law can be constructed. V. SELKOV MODEL The Selkov model l5 is described by the following mechanism: 2 3Y, Y ;B. A 2x, X+2Y (5.1) k-3 k-2 k-1 We imagine that the concentrations of the A and B species are fixed by external flows of reagents and consider the dynamics of the intermediate species X and Y.16 The chemical rate law then takes the form of two coupled equations for the X and Y concentrations iIpx -=k,p, dt

(4.1)

4% -=a--xfx dt

Jp, (LX) + V-[u,(t,x)p,(t,x)] dt

= ig ,y&ZBc~(N;afi)*

d2iY _ dt-b-&,

(4.2)

then by Eqs. (2.36) and (3.4) the continuity equation for the local particle density p7 is

(4.3)

If the system is spatially homogeneous and in equilibrium velocity space, so that for each r and each i

in

--k-,p,

-kp,p:

i-k-&

dpy - k - 3Pb - k,p, -_ (5.2) -I- k,p,d k - 2~:. dt A phase diagram for this model was constructed earlier by Rehmus et al. ” and shows some of the interesting behavior that is found for the spatially homogeneous Selkov model. If we employ a scaling of variables similar to that of Ref. 17, the rate law takes the namely, k, t-+7 and pT = is,,/=, simpler form

and

$, Cz (WP) = 0,

(4.5)

- px& i- f@;, +pxj?;

-K/s;.

(5.3)

In our simulations we follow Ref. 17 and select x = 0.1 and = 1 and consider a and b as bifurcation parameters. With the scaling in Eq. (5.3) we can read off the bifurcation structure of the homogeneous system from Fig. 2 of Ref. 17. From this figure it is clear that even in the spatially homogeneous regime the Selkov model has a very rich bifurcation structure: There are regions in parameter space where the system possesses a single steady state (it may be either a node or a focus or exhibit excitability properties), a limit cycle, and bistability between two fixed points or between a limit cycle

K

J. Chem. Phys., Vol. 96, No. 4,15 February 1992 Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

Kapral, Lawniczak, and Masiar: Reactive dynamics in a lattice gas

2768

and a fixed point. We shall now apply the formalism of Sec. IV and construct a LGCA model for the Selkov reaction. In order to construct the Selkov LGCA we consider a square lattice (m = 4) I8 with periodic boundary conditions and write Eq. (4.5) for two species as (5.4a) where K& = ( - l)l+km-[-k xi

iP,(@(-1)-‘-j i=O j-0 (54b)

with P,(ij>

=

i i (& r’= 0 f = 0

(5.5)

-a,)P(ijli’j’).

The quantities K TXin Rq. (5.4) can be identified with the rate coefficients of any two-variable reaction model, and thus Rq. (5.4) provides the connection between the phenomenological rate coefficients and the LGCA reaction probabilities. From this comparison of coefficients we obtain expressions for the P, (ij) in terms of the rate coefficients in the Selkov model. We find P,(U) =$(j-

l)(j-

-#j-

2)k-, l)k,

P,(U) = - ycj-

l)(j-

-ik-,

+k,p,,

+-7-\(J-13

CEI+~-J-I

--+ e-s-q@‘+

ii +

t

t

:

~(1) = kg,,

~(12) = 16k-,,

~(2) = k-jpb,

~(13) = 8k,, p( 14) = (32/3)k,,

p(3) = k,, p(4) = 2k,, p(5) = %, p(6) = 4k, -k p(7)=k -17 ~(8)

2)k-,

+$(j-l)k, -jk, +k-,p,. (5.6) Equations (5.6) are a set of constraint conditions that the

Q@+-B

reaction probabilities P( ijlil”) must satisfy in order for the LGCA to yield the same macroscopic behavior as the Selkov model. Of course, there are an infinite number of LGCA models that are consistent with these constraints. This is to be expected since a set of macroscopic constraints cannot fully determine the microscopic dynamics. Rather one should regard the LGCA model as a means of specifying a microscopic dynamics that is consistent with a desired macroscopic rate law. The microscopic model should then be constructed to mimic as many features as possible of the true molecular dynamics of the system of interest. In this spirit, we restrict particle changes to increases or decreases by one particle in accord with the Selkov kinetics and reactive changes to those that occur in the Selkov mechanism. When these considerations are implemented the number of nonzero elements of the 25 X 25 reactive probability matrix is very small. The structure of the reaction probability matrix is given in Fig. 1. The reactive transition probabilities that appear in the figure are expressed in terms of the quantities p(k), (k = l,..., 22) which are defined below;

c10+0 3

@b, -

- 16k,,

p( 17) = 16k,, p( 18) = 64k-,

- 32k,,

- k,p,,

~(20) = 64k

- 48k,,

~(21) = 32k,f-

16km,,

~(22) = 64k, - 64ke,. ~(11) =jk,, Figure 1 and Eqs. (5.7) completely specify the reactive dynamics of this Selkov LGCA, but other equally valid probability matrices may be constructed. Finally we note that since the reactive probabilities must satisfy O(P( ijli’J’) < 1 there are additional restrictions on acceptable values of the rate coefficients that appear in the Selkov model (or any other model). VI. SELKOV LGCA SIMULATION

t

t

t

:

92-I + *7j12-l--* 2 tf 5 ,I2 &-I-*

I

/I

3 tf P /I2

+7$l3-~

t :

cey,2-1+

++I

/II

/I3

f

tr ? /I2 -P +e;213~+

t

t 4 + +a;;42 t tf

/I4

: ; t

0 /I2

+~;3~3-1+ to-..3 P

(5.7)

p( 19) = 24k,,

=2k-,,

p(9) =3k-,, p( 10 = 4k-,

~(15) = 64k-,, ~(16) = 64ke,

RESULTS

Below we present the results of simulations carried out on the Selkov LGCA described above in order to demonstrate that the reactive lattice-gas automaton can provide a new method for the study of spatiotemporal dynamics in chemically reacting media. We have selected parameter values in regions of the Selkov phase diagram where the system has very different kinds of behavior. A. Stable fixed point and limit cycle

FIG. 1. The nonzero reaction probabilities P(ijlff) are indicated by arrows from (0) + (I*J ) . The number on the arrow in the figure refers to the index k of the quantitiesp(k) defined in Eq. (5.7).

The simplest kind of attractor is a stable fixed point in a spatially homogeneous system. However, depending on the system’s parameters, the attracting fixed point may be a node or a focus. These two characteristic kinds of behavior for the Selkov LGCA are illustrated in Fig. 2 for a model where the diffusion coefficients of X and Y are equal. In Fig.

J. Chem. Phys., Vol. 96, No. 4‘15 February 1992 Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

2769

Kapral, Lawniczak, and Masiar: Reactive dynamics in a lattice gas

system state is spatially homogeneous on average in both cases. Spatially homogeneous stable limit cycle dynamics is observed in other regions of parameter space. Figure 3 shows the evolution to an attracting limit cycle starting from two different initial conditions, one inside the limit cycle and the other outside the cycle. This limit cycle has a strong relaxation character; the parameter values have been chosen to lie far from a Hopf bifurcation point. In all of the above results the effects of local fluctuations on the average species concentrations can be observed in the figures; e.g., the irregular trajectories in the montone decay to the fixed point noted above and the fact that the dynamics is perturbed about the limit cycle attractor. 6. Excitable ,132

(al

.I48

.I44

.I40

.I 36

cx

24,

.

,

.

I

I

I

I

a

I

I

.

,

*

I



I

3

I

8

,

Especially interesting dynamics occurs when the fixed point is excitable and perturbations beyond a threshold value cause the system to make a long excursion in phase space before return to the fixed point. During this long excursion the system is not susceptible to further perturbations; i.e., it is refractory. Small perturbations which do not exceed the threshold relax quickly back to the fixed point. If thresholdexceeding perturbations are introduced locally in an excitable medium, chemical waves of excitation can form and propagate through the reacting medium. The most frequently encountered chemical waves are produced by this mechanism; e.g., those in the Belousov-Zhabotinsky reaction are typically of this type. If a perturbation is introduced locally in an undisturbed excitable medium, a ring of excitation will

.l5L

0

2.0

4.0

61)

8.0

dynamics

8

II

1

x

I

z

1.1

1

.

I

-

I

8

I.

I



I



I

’ *

10.0

PIG. 2. (a) Phase plane plot of monotone decay to the fixed point from several initial conditions. System parameters: k,p. = 0.000 190 06, k-, =O.lk,, k, =0.005, k, =k-, =O.Ol, and k~,p,=0.00000559. System size: 256X 256. (b) Plot of oscillatory decay to the fixed point. The concentrations per velocity direction c, (r = xy), are plotted vs time. Sysk,p. =O.o001897, k-, =O.lk,, k, =O.oOl, parameters: tem k, =k_, =O.Ol, and k-,p,=O.OM)O41 11 and the system size is 512X 512. The time is in units of 104 time steps.

2 (a) we show in a phase plane plot that starting from several different initial conditions the decay is monotone close to the fixed point (apart from fluctuations) indicating the existence of a stable node. In contrast, in Fig. 2(b) we show the oscillatory decay to the fixed point characteristic of a stable focus. The system is close to the Hopf bifurcation point and shows the characteristic slow decay to the fixed point. The

Or,,,,,,,,,,,,,,,,.,.,.,,1 .I6

.20

.24

.28

.32

.36

CX

PIG. 3. Liiit cycle oscillation. The evolution from two initial conditions, one inside the limit cycle and the other outside indicated by heavy dots, is shown. System parameters: k,p, = O.COO189 7, k-, = O.lk,, k, = 0.001, k, = k-, = 0.01, and k_ )pb = O.M)o 025 30 and the system size is 256X 256. The time is in units of 104 time steps.

J. Chem. Phys., Vol. 96, No. 4,15 February 1992 Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

Kapral, Lawniczak, and Masiar: Reactive dynamics in a lattice gas

2770

propagate outward. The chemical wave has a refractory tail preventing back propagation of the wave. If two such chemical waves meet they annihilate. If the medium is disturbed or is inhomogeneous the rings of excitation may be broken and the free ends will curl to form spiral waves. These phenomena account for the structures of the most commonly observed chemical patterns.‘V19 It is possible to study these chemical waves using the Selkov LGCA by selecting reaction probabilities to correspond to homogeneous kinetics within the excitable region. Just as in real excitable media, we can initiate simple ring growth by introducing a local disk-shaped, threshold-exceeding perturbation in the LGCA. Figure 4 shows the evolution of a ring from this initial condition. Note that the system relaxes back to the fixed-point concentration after the chemical wave has passed. The wave is circular indicating that the particle fluctuations in conjunction with the



.

~ ,:

,,

‘..

;:

_

: _.:

.:

:.

.,_ ; .

:

:

.” ,

. ‘.

,,:. ) ” .

-;

I

“<.

, ,,‘.

.

,,::,

_; .‘. ; :.:, ,, -; :, ‘;,.

1 ‘. :,.;

.. . . .,

..,_

,

‘..,’ I:” ;, . .,‘. .: ‘,_ ‘,’ ;

i

:, :.



;. ..”

:,.‘. I, :..I

I

“_, ..’

.

-‘,

:,

..; ‘b.. ::. ,,::,.’ . ..-

.

.

,:y

:

propagation and rotation steps are strong enough to mask any dependence of the wave shape on the underlying square lattice. The chemical waves annihilate when they collide as in the real excitable system. The structure and velocities of the automaton chemical waves have been studied. Figure 5 shows the concentration profiles of the Y species as a function of the radial distance from the center of the ring for several different times. The propagating front and refractory tail are clearly evident in this figure. The chemical waves propagate with a constant velocity in this simulation so that curvature effects on the wave propagation speed are too small to be detected in the figure. Spiral wave formation in the LGCA can be induced by shearing the ring of excitation, simulating that process in the real excitable system. In Fig. 6 we observe that the free ends of the ring curl to form a pair of counter rotating spiral waves. Thus the automaton is able to generate the two major types of chemical waves that are typical of two-dimensional excitable media: rings and spirals. In the above simulations we introduced a disk-shaped perturbation. The radius of the perturbation must exceed a critical value (critical nucleus) for the perturbation to spread in the medium. One interesting problem to investigate is whether spontaneous fluctuations in the system can produce nuclei that exceed the critical radius and thus generate spatial structures. One way to increase the probability of forming such critical nuclei is to decrease the diffusion coefficients of the two species so that when the autocatalytic step in the reaction produces a local increase in the concentration of a species diffusion is unable to quickly smooth out this local concentration fluctuation. We have verified that this is

,.

-. .:

FIG. 4. Evolution of a ring of excitation in the excitable region. First panel shows the system state shortly after the introduction of the local diskshaped perturbation while the next panel shows the ring of excitation that has evolved from this initial condition. The concentration of Yis shown as gray shades in this and the following figures which depict spatial patterns. System parameters: k,p, = O.CQO083 85, k- , = 0.1/c,, k, = 0.0005, k2 = k- 2 = 0.01, and k- gb = O.ooO002 326 and the system size is 512x512.

I

0

80

I

I

120

I

I

I

I

160

I

I

I

200

FIG. 5. Radial concentration profiles for the ring growth in Fig. 4 for five equally spaced time intervals. The concentration of Yis plotted in arbitrary units vs the radial distance r. The plots are displaced along the concentration and time axes for clarity.

J. Chem. Phys., Vol. 96, No. 4,15 February 1992 Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

Kapral, Lawniczak, and Masiar: Reactive dynamics in a lattice gas

2771

one can observe the formation and growth of spontaneous nuclei giving rise to a complex dynamic spatial state that consists of annihilating wave fragments and newly formed nuclei and rings of excitation. I. .

.

;

C. Turing patterns ‘..

,.

,.‘. :

.’ .I

.

.: ..:. ;‘. I’ ‘.,

,. : .:,

;. . ..‘. ,.., . :

._.’ .. .

.‘.

FIG. 6. Spiral wave pair formation from “sheared” ring. The free ends curl to form a vortex-antivortex pair of spiral waves. System parameters are the same as Fig. 4 but the lattice size is 1024X 1024.

the case by carrying out simulations where the diffusion coefficients of both species were decreased by a factor of 3 by performing the propagation and rotation steps after every three chemical transformation steps. The results are shown in Fig. 7. In addition to the externally induced perturbation,

Turing bifurcations occur when a stable homogeneous steady state is destabilized by a difference in the diffusion coefficients of the two species giving rise to a steady inhomogeneous state. The technical conditions for the appearance of such bifurcations are well known*’ and pattern selection has been studied for such systems.2’ While these bifurcations were predicted by Turing** in 1952 and are speculated to be relevant for pattern formation processes in living systems, it is only recently that experimental evidence for Turing structures has been presented for chemical systems.23 Such experiments open up the possibility of detailed and controlled studies of these spatial structures24 and their study is a topic of considerable current interest.*’ The reactive LGCA method presented here provides a means for carrying out detailed theoretical and simulation investigations of such bifurcations. We have selected Selkov parameters to lie in the stable steady state region near a Hopf bifurcation point where the conditions for the onset of a Turing bifurcation are favorable.26 The analysis of the Selkov reaction-diffusion equa-

FIG. 7. Spontaneous nucleation of rings of excitation in the excitable region and spiral wave formation. In addition to ring growth from a disk-shaped perturbation, followed by shearing of the ring and spiral wave formation, waves of excitation spontaneously nucleate. In the last frame one sees that the spiral wavefronts have collided and “pinched off’ a wave fragment that can serve as the nucleus of another spiral wave pair. System parameters and lattice size same as Fig. 6 but the diffusion coefficients are one third of those in Figs. 4 and 6.

J. Chem. Phvs.. Vol. 96, No. 4, IF, Fnbruarv 1992 Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

2772

Kapral, Lawniczak, and Masiar: Reactive dynamics in a lattice gas

tion predicts that a Turing bifurcation for these parameters should occur when Ox/D,, > 24.*’ The evolution of the concentrations of the Xand Yspecies are shown in Fig. 2(b) for equal diffusion coefficients. One can observe the oscillatory decay to the fixed point characteristic of this region, with small tluctuations about the fixed point concentrations. The system is spatially homogeneous under these conditions as noted earlier. Figures 8 and 9 present the results of simulations for the same kinetic parameters but with unequal diffusion coefficients. We have selected D, = 250, so that in the activator-inhibitor language species X is the inhibitor and Y is the activator (see Murray, Ref. 1 and Ref. 25). The selected diffusion coefficient ratio is larger than that predicted from the deterministic reaction-diffusion equation for a Turing bifurcation so one expects pattern formation. One can seefrom Fig. 8 that the concentrations ofXand Y tend to steady state values, but the results in Fig. 9 show that the system is spatially inhomogeneous. The figure shows the inhomogeneous final state. The spatial structure persists over long times and is occasionally perturbed by concentration fluctuations in the system. In these simulations we used rule (2.32) with + i, = 5 and - Z,, = 5. If we let D be the diffusion coefficient for the automaton updating rule PRoC, then we have D, = 5D and D,, = D/5. Fluctuations are quite strong with this selection of diffusion coefficients and one can observe perturbations of the basic stationary inhomogeneous state, but its prominent features persist in spite of these rather strong fluctuations. As the magnitudes of the diffusion coefficients are increased, always keeping the ratio fixed, fluctuations play a smaller role. Of course, the characteristic wavelength of the pattern changes since it depends on

O.lD/i

o.24 k-----0.06--i 01 0

I

2.0

I

4.0

I

t

6.0

I a0

FIG. 9. Turing pattern. System parameters are those of Fig. 8.

D, and D,, individually. Simulations have been carried out that confirm these effects. It is also interesting to examine the bifurcations that arise from the spatially homogeneous oscillatory state as a result of diffusion coefficient differences. For parameter values slightly different from those used in the Turing simulation, the homogeneous oscillatory state is stable for equal diffusion coefficients. For these parameter values, if the diffusion coefficient ratio is changed to D, = 200, a steady inhomogeneous state results. Hence, both steady and oscillatory homogeneous states can be destabilized by diffusion coefficient differences to yield steady inhomogeneous spatial patterns.** Studies employing reactive lattice-gas automaton models will allow one to probe these pattern formation processes on a molecular level in order to obtain a deeper understanding of how molecular fluctuations influence their formation and growth. VII. CONCLUSION

I 10.0

FIG. 8. Unequal diffusion coefficients. The concentrations of c, and cv vs time are shown. System parameters are the same as Fig. 2(b) but D, = 5D and D,, = D/5, where D is the diffusion coefficient for the basic RoCoP rule. System size is 5 12 X 5 12. The time is in units of l@ time steps.

The multispecies reactive lattice-gas automaton employs a simple set of evolution equations that give rise to a dynamics like that of real reacting systems. By insuring that the essential characteristics of the microscopic reactive and elastic collision processes are correctly treated a near molecular-level description of complex chemically reacting fluids can be constructed. There are several advantages to this type of approach. Since the automaton describes particle motions and reactions, specific chemical reaction schemes can be modeled in contrast to “classical” cellular automaton models where the CA rule is a gross abstraction of the dynamics.” With the exception of the generation of random numbers for the rotation and chemical transformation probabilities all the dynamics is carried out at the bit level. Hence, problems associated with numerical errors and instabilities, as is the case for solutions of the reaction-diffusion equation,

J. Chem. Phys., Vol. 96, No. 4.15 February 1992

Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

Kapral, Lawniczak, and Masiar: Reactive dynamics in a lattice gas

do not arise. There is no increase in complexity if more complicated geometries or boundary conditions are included. This is perhaps one of the most important features of the reactive LGCA method since it permits the investigation of systems similar to those commonly studied in the laboratoryez9 For example, systems with flows imposed at the boundaries, like those in the recent investigations of Turing bifurcations,“*” spiral wave dynamics,” and turbulence in inhomogeneous systems composed of catalytic bead arrays3’ can be investigated. Apart from these features the LGCA model automatically incorporates fluctuations in the dynamics. Consequently, a variety of problems that cannot be described by the reaction-diffusion equation can be explored using this technique. These include the effects of fluctuations on bifurcations in far-from-equilibrium systems, the change in character offluctuations near to and far from equilibrium, fluctuation effects on chemical wave propagation processes, to name a few. In the brief sketch of possible applications in Sec. VI we noted that fluctuations can induce nucleation of rings of excitation in excitable media and can affect spatial patterns. Commonly, fluctuation effects in systems far from equilibrium are investigated using master equation methods.’ The LGCA method provides another way to study these phenomena. It is also complementary to interacting particleJ2 and random walk33 methods that have been used to study chemically reacting systems. The effects and nature of the reactive LGCA fluctuations have not been investigated fully in this paper since our main goal was the mathematical formulation of the automaton model and a demonstration of its potential utility. However, as stated above, the fact that fluctuations are present is perhaps one of the most important features of such reactive models and a few comments on how they arise in the dynamics are in order. If the LGCA fluctuations are to resemble those of the real system then every effort should be made to build automaton dynamics that is as close as possible to that of the real system. The elastic and reactive collisions that occur at a node should model the real collisions in the system. Of course, the collisions that occur at the nodes of the lattice should not be regarded as identical to those of the real system. Rather, they should be considered as effective collision processes, much like those that arise in approximate kinetic theory descriptions of fluids, such as the BGK kinetic equation,34 where the collision operator thermalizes particles at each “collision” with a predetermined frequency. Thus some coarse graining of the real dynamics is already implied in the automaton collision processes that occur at a node of the lattice. In our lattice-gas automaton the velocity is randomized in one time step; thus from a microscopic point of view it is reasonable to associate one time step with the velocity relaxation time t, in the system. Correspondingly, the distance a particle travels between such thermalizing collisions is the mean free path il, , and the lattice spacing can be associated with this length.35 Such physical considerations can aid in the actual construction of the reactive probability matrix since the constraint conditions corresponding to a given macroscopic rate

2773

law do not fully determine the elements of this matrix. The reactive events occurring at a node should mimic those of the actual reaction mechanism given the scaling described above. While rather bizarre choices of reactive probabilities may yield the same macroscopic rate law, there is no reason to expect such models to have fluctuations similar to those of the physical system if the microscopic processes are different. In summary, if an automaton dynamics can be constructed to correspond to this mesoscopic description of the real system one may hope that the fluctuations will be represented accurately. This is a subtle point and one that will be the subject of future study.36 In addition to these general considerations it worth remarking that the magnitude of the concentration fluctuations can be tuned by varying the number of propagation and rotation steps relative to the number of chemical transformation steps. Clearly, in the limit of an infinite number of propagation and rotation steps per reactive transformation the system is completely “stirred” when reaction occurs and no local concentration fluctuation effects arise from the depletion or increase of particles due to reaction. It is this competition that can lead to the breakdown of the macroscopic rate law for fast reactions. 35 These effects have not been fully explored but some information on their magnitude has been presented earlier for the Schiigl model.13 A complete investigation will require a careful consideration of the relative magnitudes of reactive versus diffusive time scales for systems of interest. The development of the automaton equations presented here form the starting point for the theoretical analysis of LGCA models. Because of the discrete nature of the velocity space, solutions of the discrete Boltzmann equation are possible and can yield insight into the nature of bifurcations in chemically reacting systems. Furthermore, the character of the correlations among the microscopic random fields and the breakdown of the factorization approximation used to obtain the discrete Boltzmann equation and chemical rate law can be investigated in the context of these LGCA models.

ACKNOWLEDGMENTS This work was supported in part by grants from the Natural Sciences and Engineering Research Council of Canada (R.K. and A.L.) and a grant funded by the Network of Centres of Excellence program in association with the Natural Sciences and Engineering Research Council of Canada (R.K.). We would like to thank the members of the CNLS, in particular S. Chen, G. D. Doolen, and Y. C. Lee, for useful discussions. A portion of the computations were carried on the Cray X-MP at the Ontario Centre for Large Scale Computation and we thank the members of this Centre for their assistance. We would like to thank J. Weimar for reading the original version of this manuscript and pointing out a number of errors to us. We are also grateful to B. Lawniczak for assistance in preparation of the manuscript.

Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to No. AIP4,license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp J. Chem. Phys., Vol. 96, 15 February 1992

Kapral, Lawniczak, and Masiar: Reactive dynamics in a lattice gas

2774

APPENDIX A In this Appendix we give the explicit expressions for the random variables l?E (7,~). We must introduce notation to denote the velocity state for the final stoichiometric configuration. First we consider processes where particles of species X, are produced in the reaction a -fi. The new X, particles can only occupy velocity directions that are not already occupied by existing X, particles. There are q, = fl, - aL, > 0 such directions, and we denote the set of velocity indices for this process by l? + (c&r)

= Cj, ,... J,,),

(Al)

which is a subset of the complement l?( a,r)’ of I(a,r) in Cl ,...,m}. In an exactly parallel fashion we consider processes where particles of X, are destroyed in the reaction a+fl. The number of X, particles that are destroyed is lql- 1 = a, - /I,, and the set of velocity indices where particles are destroyed is denoted by r - (c&r)

= Gil ,...jlq~l 1,

t.42)

which is a subset of r (a,r) . The other species at the node (apart from X, whose dynamics is of interest) may participate in the reaction a -fl and their velocity states must be specified. We denote by B(aP,r)

= bq,

> O,u#r),

(A3)

Yf (73) =

C c c q o=maS.r)r + (a/3,0)=D(cr8,d r - (c$?,K) *) X [ P(afl,r),P(ap,r),{i,jJ],

(A91

where the prime at the sum means that we take a sum over all distinct indicesj, such that &,l=Cj,

,... jq~~Cr(a,r)c-

{il.

(A101

The random variable vf (r,i) describes all the possibilities where the particles of species X, will be created, including always a particle in the ith velocity direction, and the particles of other species will be created or destroyed in a chemical transformation a --+fl with an initial configuration of particles l?(a) such that r(a,r) i9i. APPENDIX B The coefficients b,, and d,, are expected values of the following random variables: L

(Bl)

= Ep &%?a ~) Vf (rA

for a such that ar < m, and da, = E,

and D(aL?,r) = bq, cO,a#r}. (A4) all the species X, (a# r) for which particles are produced and destroyed, respectively. Then clearly P(afl,r) = Cr + (aB,a):aEB(aB,r)I, (A51 contains all the velocity indices where particles of species X,, U# r will be produced and rD(aB,r) = Cr - (&%a) :dCaP,rl, (A61 contains all the velocity indices where the particles of species X,, a# r will be destroyed. For each input stoichiometric vector a with an initial configuration of particles l?(a) such that r (a,r) 3 i and for each &Ce - (a,r) let us denote by l@ (r,i) the following sum of random variables: ff

each ,M5 + (a,r) we denote by y?! (r,i) the following sum of random variables:

(7;i) =

C f!(r,i), 032) BE%- (a.79 for a such that a, ) 1. From Eqs. (2.16) and (A9) it follows that

L =w& c+(n,r) dZ5.r) l-+:5n, (B3) where the prime at the sum means that we take a sum over all distinct indices j, such that CjJ=Cj,

,... &~ClW,rY--

Gil.

(B4)

Since for a such that a, < m, 55’+ (a,r) is a set of all output stoichiometric vectors 0 such that in a reaction a-/3, 8, > a,, then Eq. (B3) can be written as follows:

C c c c cr”” oEB(crS.r) ‘- + (cq3.0)=m-=@*er - (a5.x) %,I)

x [ rB(ap,r),rD(aa,r),{ij,,, I],

(A7)

where the prime on the sum means that we take a sum over all distinct indicesj, such that

ci,,,,}~{i,,...j,,,}Cr(a,r) - {il.

(A81

tB5) whereq, =fl,

-aa,,andbyEq.

b,, = 2 . . . 2 . . . i 5,=o &=a,+1

5,=0

The random variable f! (r,i) describes all the possibilities where the particles of species X, will be destroyed, including always a particle in the ith direction, and particles of the other species will be created or destroyed in the chemical transformation a-p with the initial configuration of particles r(a) such that l?(a,r) 3i. For each input stoichiometric vector a with an initial configuration of particles r(a) such that I (a,r) 3 i and for

Similarly we can computed,, (2.16) and (A7) da, =

2 &‘6 - (qr) xc

c

=maS.dr

(2.12)

8T - - ar P(Cra). m

(B6)

-a,

for a such that ar > 1. By Eqs.

2 + (crp.U)

c-pie, ndxa5?e r- (&,K) ci,,,,} n (aLo

J. Chem. Phys., Vol. 96, No. 4,15 February 1992 Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

(B7)

Kapral, Lawniczak, and Masiar: Reactive dynamics in a lattice gas

where the prime at the sum means that we take a sum over all . . . . distinct mdtces j,, such that

2775

velocity space, so that for each r and each i p,(t) =-=c m then Eq. (Cl ) becomes N,(t,x)

Cj,,,}~{j2,...j,,,}Cr(a,r)

-{il.

038)

Since for a such that a,) 1, %’- (a,r) is a set of all output stoichiometric vectors p such that in a reaction a-+fl, /3, < ar, then Eq. (B7) can be written as follows:

a,- I a7 - 1 ... c ... m CC a, -& - 1 P.=~ 51=o 5,=0

y

((3)

79

= i$l (a5j$ZpCWL3. t ,

(C3)

Also, Eq. (C2) implies that Eq. (2.20) becomes

da, = i

Q”[N,r(a)]

= fi cp(l -c,)“-“*

X=1

tB9) =

and by Eq. (2.12) da, =

J,

cq-- 1 2 . -. 2 -*- cm -a,-& 5, =o P,=O 5.=0 a,

40, (4x)

at

-

l)‘-“~~mPL;)m-‘p,o!

x

P(a;B).

(BIO)

APPENDIX C The continuity

(

,&

(C4) Since the number of different configurations I( a), such that I( a,r) 3 i, for a such that a, < m is given by the formula n m = I!Y5. t-1 , m +,Jd a, > and the number of different configurations I’( a), such that I (a,r) 3 i for a such that a, > 1 is given by the formula n,+ (qj)

equation for the local particle densityp,

is

+ V*[u,(t,x)p,(t,x)]

= J, (aa)z aZ5c: (N;afi)*

If the system is spatially homogeneous and in equilibrium

B w)~~.~+,YC~(N;~P)

= L,
m 1

n; (a,i> =%

fi m , K=l ( a, > then by Eqs. (3.5) and (C4) for each i

(Cl) in

(a,i)Q”[W(a)l

2.o.a,
- a,>W=,

- &:a,sldarnr ,” cf( I( c>

1 ,” cf( 1 -c,)“-=* ( I(> Using Eqs. (B6) and (BlO) for b,, and d,,, respectively, we obtain - La,>lL4’C=

z: (a~)E~r:.aZ5C~WiCZB) =L

m

~a:rr,4m~~=O~~~~~,=u,+1~~~~~~=O(~,

(C6)

(a,i)Q”[N;r(a)l

cK)m-=K

I

.

(C7)

-a,MaB

-2 ,:,,8,“1=o-*x;;~;.-.Z;n=o(P, -a,VY&)IK=’

Hi.=,

cP(l

cYp( 1 -c,)“-“”

-c,)“-=” ,

(C8)

which can be written compactly as (C9)

(&5)E~f;P#5 concentrationp,

of species XT, i.e., Eq. (4.5)

I ’See, for example, G. Nicolis and I. Prigogine, Self Orgunizafion in Nonlinear Systems (Wiley, New York, 1977); Oscillations and Traveling Waves in Chemical Systems, edited by R. J. Field and M. Burger (Wiley, New York, 1986); P. Berg&, Y. Pomeau, and C. Vidal, Order Within Chaos (Wiley, New York, 1984); J. D. Murray, MathematicalBiologV (Springer, New York, 1989).

‘The effects of fluctuations on such processes are often treated by a master equation approach. See, for instance, Nicolis and Prigogine, Ref. 1 and C. Vidal and H. Lemarchand, La R6action Crkatrice (Hermann, Paris, 1988). 3J. Boissonade Phys. Lett. A74,285 (1979); PhysicaA 113,607 (1982); P. Ortoleva and S. Yip, J. Chem. Phys. 65,2045 (1976).

J. Chem. Phys., Vol. 96, No. 4,15 February 1992 Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

2776

Kapral, Lawniczak, and Masiar: Reactive dynamics in a lattice gas

“R. Kapral, A. Lawniczak, and P. Masiar, Phys. Rev. Lett. 66, 2539 ( 1991); R. Kapral, A. Lawniczak, and P. Masiar, in Lectureson ThertnodynamicsandStatisticat Mechanics, edited by M. L. de Haro and C. Varea (World Scientific, NJ, 1991), p. 44. ‘U. Frisch, B. Hasslacher, and Y. Pomeau, Phys. Rev. Lett. 56, 1505 ( 1986); D. d’Hum%res, B. Hasslacher, P. Lallemand, Y. Pomeau, and J. P. Rivet, Complex Systems 1,648 ( 1987). ‘Lattice Gas Methodsfor Partial Dlrerential Equations, edited by G. Doolen (Addison-Wesley, New York, 1990); D. d’Humitres, Y. H. Quain, and P. Lallemand, in Discrete Kinetic Theory, Lattice Gas Dynamics and the Foundations of Hydrodynamics, edited by I. S. I. and R. Monaco (World Scientific, Singapore, 1989)) p. 102; see also articles in Lattice Gas Methods for PDE ‘s, edited by G. Doolen, Physica D 47, ( 199 1) . ‘There have also been studies of turbulent reactive flows: P. Calvin, P. Lallemand, Y. Pomeau, and G. Searby, J. Fluid Mech. 188,437 (1989); V. Zehnle and G. Searby, J. Phys. (Paris) 50, 1083, (1989); and to heat transfer and chemical reaction: 0. E. Sero-Guillaume and D. Bernardin, Eur. J. Mech. 9, 177 (1989). ‘G. A. Bird, Molecular Gas Dynamics (Clarendon, Oxford, 1976). 9 F. Baras, J. E. Pearson, and M. Malek Mansour, J. Chem. Phys. 93,5747 (1990). ‘OR. Kapral, J. Math. Chem. 6, 113 ( 1991). “D. Dab, J. P. Boon, and Y.-X. Li, Phys. Rev. Lett. 66,2535 (1991). 12The exclusion principle can be implemented in other ways, for example, by allowing only one particle, regardless of species type, to occupy a node with a given velocity. The manner in which the exclusion principle is applied can have consequences for the range of rate laws that can be investigated. For further discussion see Ref. 11. 13D. Dab, A. Lawniczak, J.-P. Boon, and R. Kapral, Phys. Rev. Lett. 64, 2462 ( 1990); A. Lawniczak, D. Dab, R. Kapral, and J. P. Boon, Physica D 47, 132 (1991). 14B. Chopard and M. Droz, J. Stat. Phys. 64, ( 1991). “E. E. Selkov, Eur. J. Biochem. 4, 79 ( 1968). ‘6This formulation, which does not explicitly consider the dynamics of the A and B species, suppresses fluctuations in these species. A more general treatment of this model where the dynamics of all four reacting species is followed will be given elsewhere. In this more general formulation the exclusion principle leads to the restriction of equal forward and reverse rate coefficients. “P. H. Richter, P. Rehmus, and J. Ross, Prog. Theor. Phys. 66, 385 (1981). ‘*With these rules on the square lattice there is a checkerboard invariant. See Ref. 13 for details.

“A. T. Winfree, SIAM Rev. 32, 1 ( 1990). r” See, for instance, Murray, Ref. 1. *’ D. Walgraef, G. Dewel, and P. Borckmans, Adv. Chem. Phys. 49, 3 11 (1982). “A. M. Turing, Philos. Trans. R. Sot. London, Ser. B 237,37 (1952). 23V Castets, E. Dulos, J. Boissonade, and P. DeKepper, Phys. Rev. L&t. 64,2953 (1990). “4. Ouyang and H. L. Swinney, Nature 352,610 (1991). ‘s I. Lengyel and I. Epstein, Science 251, 650 ( 199 1) . 26J Pearson and W. Horsthemke, J. Chem. Phys. 90, 1588 (1989); A. Arneodo, J. Elezgaray, J. Pearson, and T. Russo, Physica D 49,141 ( 1991). *‘The diffusion coefficient ratio at which a Turing bifurcation is observed can be reduced by considering parameters closer to the Hopf bifurcation point. We have not attempted such simulations where the relaxation times are longer since our intention was simply to demonstrate the feasibility of Turing bifurcation studies. *ax. Y. Li, Phys. Lett. A 147,204 (1990); A. Rovinsky, J. Phys. Chem. 94, 7261 (1990). 29Experimental situations like those in gel reactors and continuously fed unstirred reactors can be studied. See, W. Y. Tam, W. Horsthemke, Z. Noszticius, and H. Swinney, J. Chem. Phys. 88, 3395 (1988). “G. Skinner and H. L. Swinney, Physica D 48, 1 (1990). 3’J. Maselko, J. S. Reckley, and K. Showalter, J. Phys. Chem. 93, 2774 (1989). 32T. M. Liggett, Interacting ParticleSystems (Springer, New York, 1985); A. De Masi, P. A. Ferrari, and J. L. Lebowitz, J. Stat. Phys. 44, 589 (1986). 33See, for instance, D. Ben Avraham, Science 241, 1620 (1988); M. A. Burschka, C. A. Doering, and D. Ben Avraham, Phys. Rev. Lett. 63,700 (1989); B. J. West, R. Kopelman, and K. Lindenberg, J. Stat. Phys. 54, 1429 (1989). 34P. L. Bhatnagar, E. P. Gross, and M. Krook, Phys. Rev. 94,511 (1954). “In a dense fluid t, -ps and 1, - A. The chemical relaxation time rcchc,,, for our simulations is rchern- lO’t,, corresponding to a fast reaction. For such fast reactions our results show that macroscopic self-organization which is robust under internal fluctuations occurs on the scale of hundreds of angstroms. ‘60ne may also scale the observed LGCA properties, like the wave velocities and diffusion coefficients to correspond to a particular physical system and from this infer the length and time scales in the automaton discretization. The fluctuations will not necessarily be accurately given in this case without further coarse graining or other modifications in the dynamics.

J. Chem. Phys., Vol. 96, No. 4,15 February 1992

Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

Reactive dynamics in a multispecies lattice-gas ...

phenomena occur on long distance and time scales. It is well known that in such .... chemical rate laws and reaction-diffusion equations. T'he links established ...

2MB Sizes 1 Downloads 188 Views

Recommend Documents

selection and dispersal in a multispecies oak hybrid ...
Berkeley, California 94720 .... from Northern California are complex, suggesting introgres- .... veloped by the University of Montana, Numerical Terradyn-.

Reactive Distillation - Aussie Distiller
without attracting attention as a different class of operation. It was not until the. 1980s, thanks ...... of a formaldehyde water methanol mixture taken with an online technique with a 400 MHz NMR ...... An animation of the simulation can be viewed 

Reactive-short.pdf
logic. 1/2 of the bugs reported during a. product cycle exist in this code. --- quoted in Martin Odersky's. “Deprecating the Observer Pattern”. Page 3 of 206 ...

Reactive Distillation - Aussie Distiller
Schoenmakers and Bessling (Chapter 2) give an overview of the tools that are ...... cells have been coupled with spectroscopic analytics are described by ...

Reactive Distillation - Aussie Distiller
10.4.2 Continuation Analysis of Industrial Size Distillation Column 251. 10.5. MTBE and ...... SISER, a software product from the first EU-project. 2.3. Process ...... Hoffmann [11], good agreement was found between model prediction and experi-.

Multispecies grand-canonical models for networks with ...
Jan 13, 2006 - 3Dipartimento di Scienze Matematiche ed Informatiche, Università di Siena, Pian dei Mantellini 44, 53100 ... by the second-order correlations between vertex degrees, ... been recently shown 3 to display an additional type of.

Reactive Distillation - Aussie Distiller
tors and their colleagues with financial support from the Kompetenznetz Verfah- renstechnik Pro3 e.V., Germany, which is ... their support in collecting the manuscripts and in standardizing formats. Last but not least, we are very thankful to ......

Task Dynamics and Resource Dynamics in the ...
source-dynamic components that have not been considered traditionally as ... energy for the act (Aleshinsky, 1986; Bingham, 1988; Bobbert,. 1988; Van Ingen .... As an alternative explanation, Kugler and Turvey (1987) suggested that periods ...

Task Dynamics and Resource Dynamics in the ...
the Royal Society, London, Series B, 97, 155-167. Hinton, G. E. (1984). Parallel computations for controlling an arm. Journal of Motor Behavior, 16, 171-194.

Temporal dynamics of genetic variability in a mountain ...
*Département de biologie and Centre d'études nordiques, Université Laval, 1045 avenue de la Médecine, Québec, Canada G1V. 0A6, †Departamento de ... population monitoring and molecular genetic data from 123 offspring and their parents at. 28 mi

Firm Dynamics in Manufacturing and Services: a ...
Jan 15, 2007 - Phone: +39 06 4792 2824, Fax: +39 06 4792 3720. ... firms, even if not at the same pace, giving rise to a highly fragmented productive system.1 ..... in small scale business.10 The turbulence degree is also high, since ... distribution

Dynamics and Bifurcations in a Silicon Neuron
department of Electrical and Computer Engineering, Georgia Institute of. Technology ..... 5: Bifurcation in numerical integration: The bifurcation of the theoretical ...

Risk Premiums and Macroeconomic Dynamics in a ... - CiteSeerX
Jan 11, 2010 - problems generating suffi ciently large premiums and realistic real ... structure of the model based on real US data over the period 1947q1&2009q1. ... shows how this variation affects the predictive power of the price& dividend ratio

What Drives Housing Dynamics in China? A ... - Georgetown University
applied to study other shocks, such as fiscal, monetary, news or technology shocks. .... policies to either stimulate housing development or to prevent the housing ...

Spatially heterogeneous dynamics in a granular system ...
Dec 27, 2007 - system near jamming. A. R. Abate and D. J. Durian. University of Pennsylvania, Philadelphia, Pennsylvania 19104. (Received 23 August 2007; ...

Risk Premiums and Macroeconomic Dynamics in a ... - CiteSeerX
Jan 11, 2010 - T1+T3 !+D-+%-(" 1.72 3.46 0.95 1.33 0.97 0.98 1.00 0.86 &0.15. 3.95. &0.20 ...... volatility for stocks varies substantially over the business cycle.

Interest Rates and Housing Market Dynamics in a ...
Jul 21, 2017 - estimate the model using data on home listings from San Diego. .... willingness-to-pay for local amenities such as school quality, crime,.

Risk Premiums and Macroeconomic Dynamics in a ... - CiteSeerX
Jan 11, 2010 - Finally, Section 6 presents the results on the implied time variation ..... the index K. Dividends are defined as total income minus the wage bill (spot wage plus .... volatile compared to the data (0.98 versus 1.34), but the model ...