Nucleation, domain growth, and fluctuations

in a bistable chemical system

Daniel Gruner and Raymond Kapral Chemical Physics Theory Group, Department of Chemistry, University of Toronto, Toronto, Ontario M5S lAI, Canada

Anna Lawniczak Department of Mathematics and Statistics, University of Guelph, Guelph, Ontario Nl G 2 WI, Canada

(Received 6 April 1993; accepted 17 May 1993) Phaseseparation and nucleation processesare investigated for a bistable chemical system. The study utilizes a reactive lattice-gas cellular automaton model to provide a mesoscopicdescription of the dynamics. Simulations of steady-statestructure, wave propagation, and critical nucleus size using this model are compared with results basedon the deterministic equationsof motion. The dynamic structure factor is computed for evolution from the unstable state and the effects of correlations are examined for early and late times. The study provides insight into these processesin a fluctuating, extendedmedium and also provides a test of the ability of the reactive lattice-gas method to describe the fluctuations in the system.

I. INTRODUCTION

Domain growth and nucleation in bistable systems with a nonconservedorder parameter have been studied often. In the early stage of the evolution, either from an unstable state or from a metastable state, molecular fluctuations play an important role and the standard description is in terms of a Langevin equation with form’

qbkt)

ar=-

63qqvrJ) 1 ) &$(r,t)

(1)

’ ‘(rJ)’

(2)

Here $(r,t) is the order parameter field, D is the diffusion coefficient, and &,t) is a Gaussian random force with white noise spectrum. In the late stage, once well-defined domains have developed, a macroscopic deterministic model is appropriate. The deterministic equation obtained by neglecting the noise term in Eq. ( 1)) is the timedependentGinzburg-Landau equation,

@ (r,t)

at

-b3(r,t)

+@(r,t) +v+ DV2qS(r,t)

where V( C#J) = -v~-&c$~+$$~ is the potential for this one-dimensional order-parameter system. If the stable states have the same stability ( Y=O) domain growth is controlled by the diffusive evolution of the curved interfaces that separatethe stable phasesand the characteristic domain size grows as a power law.lT2 Such phase separation dynamics can be observed in bistable chemically reacting systems. One of the moststudied systems of this type is the Schliigl model3 whose mechanism is 3938

J. Chem. Phys. 99 (5), 1 September

(4)

h 2X+ B+3X. k-2

The reaction-diffusion equation !@$Lk,a--k-,p(r,t)

+k2bp2(r,t)

---2p3(rJ)

+ DV2pW>,

with F the free energy functional

-=

4 As-X, k-1

1993

(5) (with proper resealing) has the same form as the timedependent Ginzburg-Landau equation. Here p(r,t> is the local concentration of speciesX. The mass-action rate law is the spatially-homogeneousversion of this equation obtained by neglecting the diffusion term. The Schlogl model has been studied often in various contexts. Its stochastic dynamics has been modeled using master, Fokker-Planck and Langevin equationsb7 Discrete models that mimic the reactive processesin the mechanism have also been used to investigate the dynamics.8’9 Numerous experimental studies of wave propagation in bistable chemical systems have been carried out; in the present context the work on the iodate-arsenousacid system is noteworthy since its overall kinetics is describedwell by cubic kinetics for a single chemical intermediate.” In this article we study phase separation dynamics, wave propagation, and nucleation in a lattice-gas cellular automaton model of the Schliigl reaction.’The automaton dynamics is a coarse-grained,mesoscopicmodel of the full reactive molecular dynamics and, as such, incorporates the effects of molecular fluctuations. The paper is organized as follows: In Sec. II we give a brief overview of the reactive lattice-gas cellular automaton method and present the details of the Schliigl automaton model. Wave propagation and phase separation dynamics starting from the unstable state are investigated in Sec. III. Here we examine the structure of the chemical wave front, its velocity, and the

0021-9606/93/99(5)/3938/8/$6.00

@ 1993 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

Gruner, Kapral, and Lawniczak:

Nucleation,

relative stability of the two phases in the spatiallydistributed medium. We quantitatively characterize the evolution from the unstable state by the nonequilibrium correlation function and dynamic structure factor. The critical nucleus size and spontaneous nucleation of the more-stable phase in a sea of the less-stable phase are also studied in this section. Section IV contains our conclusions. II. SCHLijGL

LA-l-l-ICE-GAS

= (PO R)“” C.

Assuming that the X particles are uniformly distributed on the lattice with average density p, the time rate of change of p is given in the mean field limit by

1) -p(t)

= 5 Wa)P(a$)

3939

coordination number of the lattice. Since this equation is a polynomial in the X species concentration, like the massaction rate law, P can be constructed so that mass-action rate law and automaton mean field rate law coincide. The elements of P are constructed to correspond to the steps in the Schliigl mechanism and satisfy the exclusion principle. The form of P we employ in this study is

AUTOMATON

A detailed description of the strategy for the construction of reactive lattice-gas automata was given earlier;’ here we outline the main features of the model and give details of the S&log1 model used in the present study. For extensions to multicomponent systems see Ref. 11. Molecules of the chemical species X are assumed to reside on the nodes of a square lattice L? labeled by the (now) discrete index r with discrete velocities i oriented along the lattice directions. The constrained chemical species, as well as the solvent in which the reactions take place, are not considered explicitly; instead, their effects are accounted for indirectly in the dynamics of the X species. We require that no two particles at a node have the same velocity; this exclusion principle restricts the maximum occupancy of a node to the coordination number of the lattice. Elastic collisions change the velocities of the particles and are described by a velocity randomization operator R that randomly rotates the velocity configuration independently at each node of the lattice. The rotation angles 0, n/2, n, 3~/2 have equal probabilities of l/4. Reactive collisions change the numbers of X particles as well as their velocities. The particle number changes are described by a chemical transformation operator C that prescribes the change due to reaction from a configuration with (r particles to a configuration with p particles. This change is specified by a reaction probability matrix P with elements P(ar 1p). Of course, the diagonal elements of this matrix give the probability that no reaction takes place. Following the change of speciesnumber, the velocity configuration is randomized to complete the reactive collision. Thus, the elastic and reactive collisions are described in the automaton by the product RO C. Particle motion is described by a propagation operator P which moves particles one lattice unit in the directions spectied by their velocities. The diffusion coefficient can be tuned by performing n propagation and randomization steps for each chemical transformation step, so the automaton rule takes the form @)

dp pp(t+

domain growth, and fluctuations

f

( l-;)“-a,

(6)

where the discrete-time difference is approximately equal to p for small reaction probabilities. Here m =4 is the

(7)

where

r;t-(a)=k1a,

rz+(a)=k2-

(mM-1)

a(a-l)b,

rl(a> =k-la,

(8) m2 (m-l)(m-2)

C(a)=k-2

da-l>(a-2>.

This reaction probability matrix differs from that in Ref. 9. Here we make sure that only those particle number changes that are specified by the Schliigl mechanism appear in P. This corresponds to a more microscopic view of the reacting system, and is more appropriate for the study of fluctuation effects since local particle number changes are governed by the steps in the reaction mechanism. Due to the presence of fluctuations the automaton dynamics is not mean field; however, provided the mean field description is not grossly violated it can be used as a guide to construct automaton models for particular parameter regimes. Furthermore, comparison of the full automaton dynamics with the mean field results allows one to determine the validity of the mean field rate law within the context of the automaton model. The calculations presented below were carried out on square NX N lattices with N= 5 12 and periodic boundary conditions. The length scale of the local inhomogeneities in the automaton simulations is determined by the relative magnitudes of diffusion and reaction. In this study we have in automaton units of lattice units taken n =6 (D=3/2 squared per time step). This choice eliminates most smallscale, short-time reactive recollision events which can lead to significant correlation corrections to the steady state average concentrations. In Ref. 9 (for a different reaction probability matrix) it was shown that the deviation between the steady state automaton and deterministic results decreasesas n increases. Ill. COMPETITION

BETWEEN STABLE STATES

The homogeneous steady states of Eq. (6) are given by the solutions of k,a-k-lp+k2bp2-k-2p3=0, and are sketched in Fig. 1 as a function of k-, for fixed values of the other parameters (cf. Fig. 1 caption). The two stable steady states will be called p1 and p2 while the unstable steady state will be called po. Since the one-dimensional order parameter dynamics derives from the potential

J. Chem. Phys., Vol. 99, No. 5, 1 September

1993

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

Gruner, Kapral, and Lawniczak:

3940

Nucleation,

-

0.012

0.010 0.0

0.014

0.016

0.018

0.020

~“-7

k,b

p3+q

k-2

p4,

the stability of the coexisting states can be deduced from the relative depths of the potential minima. In an infinite spatially-distributed system supercritical nuclei of the more-stable phase will always form and consume the lessstable phase;thus, the deterministic, less-stablestate is always metastable in the infinite system. Nevertheless, its lifetime can be extremely long, and in such a case fluctuation-induced transitions will occur rarely in simulations on finite systems. This is illustrated in Fig. 1, where the deterministic steady state solutions are compared with finite-duration automaton simulations. In these simulations the initial state was a random distribution of particles over the nodes of the lattice with mean concentrations correspondingto the deterministic fixed point values. The system was allowed to relax for 1000time steps and the mean concentration was determined by an averageover the lattice and over a time T=2000 steps, 4

T

P=T-’

x1 iv-2

c

rs3

(a)-

-i,

‘$T”

FIG. 2. Wave propagation in the S&log1 lattice-gas automaton, for systems at (a) k-,=0.0160; (b) k-,=0.0205. The lattice elements pictured here are the result of 2X2 coarse-granting on the original lattice.

FIG. 1. B&ability region for the Schlijgl model. The solid line represents the deterministic steady state concentrations, c= p/4 (expressed as particles per cell per velocity direction), while the points are obtained from the lattice-gas simulations. The error bars represent *one standard deviation. The other (fixed) system parameters are k,a=O.OOl, $c& =0.095, 16k-?=0.245.

k-I

.

0:022

k-1

V(p) = --k+zp+T

domain growth, and fluctuations

c Mr,tL

(10)

i=l

where N2 is the number of lattice sites, and qi(r,t) is a Boolean variable that is unity if the velocity direction i at node r at time t is occupied, and zero if it is unoccupied. On the long time scale of the simulation the averageconcentration p corresponds closely to the deterministic steady state concentrations. The small deviations in the interiors of the bistable domain and the larger deviations near its boundaries are due to local fluctuations in the X concentration. Fluctuation-induced transitions between steady state branches occur only at points very close to the edgesof the bistable domain on the time scale of our simulations.

A. Wave propagation

Wave propagation processes in a bistable spatiallydistributed system present a number of interesting features depending on the system parameter values. The phenomena are summarized by Harding and ROSS’~in a study of wave propagation in an optical bistable system; the effects of fluctuations are also examined using master and Fokker-Planck equations. To study simple wave propagation processeswe consider a class of initial conditions where each half of the lattice has averageconcentration equal to one of the deterministic stable states p(r)=

p1 for [r= (x,y):l
l
(11)

The wave front that is establishedbetween these two stable states will move in a direction determined by the relative stability of the two phases;the more-stable phase will consume the less-stablephase.The point of zero wave velocity correspondsto equistability of the two phases.This feature of wave front propagation has been taken as the definition of relative stability in spatially-distributed systems.3v’2 Examples of wave propagation for the parameter values indicated by arrows in Fig. 1 are shown in Fig. 2. The wave front moves with constant averagevelocity and the velocity as function of k-, is plotted in Fig. 3. The velocity is computed by averaging the concentration over y and determining the position of the wave front as a function of time [see Fig. 4(a)]. The point of equistability extracted from this graph is k-,=0.0195, which is close to the deterministic value of kwl =0.0197. The small difference presumably results from the same correlations that lead to small differences between the automaton simulations and the deterministic steady state values on the upper branch in the bistable domain. These results support the conclusion that the deterministic equistability condition determines the relative stability of the two phases in the inhomogeneous system.3”2 It is interesting to contrast the front propagation process for parameter values inside and outside the determin-

J. Chem. Phys., Vol. 99, No. 5, 1 September

1993

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

Gruner, Kapral, and Lawniczak: Nucleation, domain growth, and fluctuations

(a)“’ 1 0.5 , -‘-

I

I

100

200

0.4

n

0

300

400

500

300

400

500

X

0.016

0.018

0.020

0.022

FIG. 3. Wave velocity as a function of k-,. The points represent values measured from graphs such as Fig. 4(a). The solid line simply joins the points as a visual aid. c

istic bistable domain. Figure 4 shows the results of two automaton simulations starting from initial conditions ( 11) . Figure 4(a) is for parameters inside the bistable domain, and shows the expected advance of the wave front such that the more-stable phase consumes the less-stable phase; the front moves with constant velocity and shape (until collision betweenthe two counterpropagating fronts occurs). Figure 4(b) is for parameter values outside the deterministic bistable domain and shows the “crash” of the system directly to the stable state. Here, of course, the upper stable state does not exist and the mean concentration on that half of the lattice was taken to lie near the previously-stable steady state. Such behavior has been observed by Harding and Ross’* in their study of an optical bistable device. The structure of the interface can be determined from an examination of Fig. 2. Near equistability, evidence of long wavelength fluctuations transverseto the front can be seen [Fig. 2(b)]; such fluctuations are much less pronounced when the front propagates with higher speeds, away from the equistable point [Fig. 2(a)]. Roughness fluctuations, primarily normal to the interface, have been investigated for the Schliigl model.13The shape of the interfacial profile can be determined easily from Eq. (5) and for the casewhere the two phaseshave the same stability is given by13-” (12) wherepo=k2b/3k-,, and Q=[(k2b)2/3k-2-kk_l]1’2. Figure 5 comparesthis equation with the automaton results. It is clear from an examination of this figure that the automaton simulations for the averageinterfacial profile are in close accord with the predictions of the deterministic model. In contrast to results obtained for the Fisher equation,16 there is no discrepancy between the interfacial pro-

~; 100

0

200 2

FIG. 4. X density per velocity direction (c=p/4), plotted as a function of the column position in the lattice. c represents an average over the y direction and over two realizations. The different traces represent different times during propagation of the automaton, as indicated in the figures. See text for discussion. The initial conditions are given by Eq. ( 11) . (a) k-,=0.0160; (b) k...,=0.0220.

file obtained in the automaton simulation and the predictions of the one-dimensional deterministic model, which should be valid for a planar interface. B. Evolution

from the unstable

state

If the system is prepared in the unstable state, after a short time period, during which fluctuations play an important role, well-defined domains of the two stable phases form and grow. The long-time dynamics is driven by the curvature of the boundaries separatingthe stable phases.In an infinite system, although the average order parameter ( p - po) is zero, domains of arbitrarily large size exist in the system. In finite systems, different realizations of the evolution processlead to pure p1 or p2 phases(or m ixtures of these phases separatedby planar interfaces), but averaged over realizations the order parameter will again be zero. The domain evolution process is monitored by the nonequilibrium correlation function C(v)

=

(

N-2

r,sy

p(r+r’,t)pW,d

>

,

(13)

where the angle bracket denotes an averageover different realizations of the evolution process.Its Fourier transform is the dynamic structure factor S(k,t). Assuming that the domain size R(t) is the only characteristic length in the

J. Chem. Phys., Vol. 99, No. 5, 1 September

1993

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

Gruner, Kapral, and Lawniczak:

Nucleation,

domain growth, and fluctuations

0.4

0.3

0.2 C

C

0.2

0.1

0.1

0

0

100

300

200

400

500

0

I

I

I

400

800

1200

2

9 1600

I 2000

t

FIG. 5. Interphase profile at equistability. The ordinate represents column-averaged X density per cell per velocity direction, c=p/4. The abscissais the column position in the lattice. Initial conditions are given by Eq. (11). The thick line is given by Eq. (12), and the other superimposed traces represent the state of the automaton after 400,800, 1200, and 1600 iterations. Note the stability of the interphase profile.

system, and that its dynamics is governedby the interfacial curvature, R (t) - t “‘, the dynamic structure factor is predicted to have the scaling form’ S(k,t) -tF(kP),

in two dimensions. Here F is a universal scaling function. We have carried out automaton simulations of the phaseseparation dynamics following a “critical quench” to the unstable state. The system was started at the equistable point k-,=0.0195, and uniformly seededwith X particles so that the average concentration corresponded to the deterministic unstable steady state, pO. For this unstable steady state concentration we verified that a roughly equal number of realizations evolved to the two stable states (or mixtures of them). Figure 6 shows the evolution of the averageconcentration for 20 realizations and demonstrates this feature. The realizations for which the asymptotic concentration was not p1 or p2 were verified to correspond to situations where the two stable phases were separated by planar interfaces. An example of domain formation and growth is shown in Fig. 7 for one realization of the evolution process. At early times the boundaries are indistinct, while at later times sharp boundaries form and evolve slowly due to diffusive motion of the interfaces. In the realization shown in Fig. 7 the final configuration shows coexisting phases separatedby a nearly planar interface. The automaton results for S(k,t) are shown in Fig. 8, where log[S(k,t)/t] is plotted vs kt”2. From Eq. (14) this should be a universal function for large t and small k (but still large compared with the inverse system size) where the scaling assumption of a single length scale is expected to be valid. The plot shows results for several different times. Except for the shortest time (t= 600)) all results are in accord with the predictions of the scaling theory for the late-stage growth. These results may be compared with a calculation of the structure factor for deterministic dynam-

FIG. 6. Time evolution of the average lattice concentration, c, for 20 realizations of a critically-quenched system. Here k-1 =0.0195, corresponding to equistability, and the initial density is that of the deterministic unstable state [c(O)=O.lSS 9161.

its described by a coupled map lattice model for a system with cubic kinetics (cf. Fig. 11 of Ref. 15). Both the automaton and coupled map lattice calculations predict the same structure of the scaling relation. These calculations confhm that the correlations in the late-stage growth are properly described in the automaton. In addition, the results provide information on the early-stage growth, which is not given accurately by deterministic models, since they do not incorporate fluctuations which are crucial in this stage of the evolution. C. Nucleation

Nucleation processesprovide the mechanism for the evolution of the system from less-stable to more-stable states. While fluctuations are responsible for homogeneous nucleation processes,an estimate of the size of a critical nucleus is easily made on the basis of the deterministic equations. For a disk of the more-stable phase in a sea of the less-stablephase, the evolution of its radius is given by i (t) = v - D/R ( t) , reflecting the competition between the linear growth with velocity v and the tendency of the diffusive dynamics to reduce the curvature of the disk, which is given by the term D/R. Here we use the symbol R(t) again but now for the disk radius. The time evolution of the radius follows from the integration of this equation and is given by the parametric equation15 t=

R(t) --R(O) V

R,

(15)

where the critical radius is given by the condition that ii=0 and is R,= D/v. In order to investigate nucleation processesin the automaton we studied the evolution of an ensemble of disks with initial radius R(0) =R, of the more-stable phase in a sea of the less-stablephase, and monitored their eventual fates. For R,< R, the disk should shrink and disappear while if R0 > R, it will grow and consume the less-stable

J. Chem. Phys., Vol. 99, No. 5, 1 September

1993

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

Gruner, Kapral, and Lawniczak: Nucleation, domain growth, and fluctuations I

I

3943 I

I

I

I 600 1200”‘1800 2400 3000 -

0.011’. 0

0.2

0.4

0.6

0.8

1

1.2

.I .

1.4

&‘I2 FIG. 8. Dynamic structure factor at different times, as indicated in the figure. The traces represent the averageof 26 realizations of a critically quenchedsystem (for details see Fig. 6).

interesting to note the deviation of R, from that predicted by the deterministic theory. For k- 1=0.0160, the measured EC is somewhat larger than the deterministic value, while the converseis true for k-, =0.0205. These differencesreflect the effectsof fluctuations on theserather small critical nuclei. As a further test of the predictions of the deterministic theory, F ig. 9 presentsa comparison of the lattice gas results for R(t) at k-t =0.0160 with Eq. ( 15), for different values of R,. Note the close fit, especiallyfor large values of R, and for the measuredR,= 19. The observeddeviation from the deterministic prediction at small R (shorter times), is the result of the interplay of several effects; the effects of a transient tim e during which the disk initial

I

I

I

200 !r

FIG. 7. Phase separation dynamics of a critically-quenched system (k-i =0.0195). The times of the snapshotsare indicated in the figure. The lattice elementspictured here are the result of 2x2 coarse-grainingon the original lattice.

phase. O f course, in a fluctuating m e d ium these conclusions will be valid only in a statistical sense,since some realizations of the evolution process will lead to growth and others to shrinkage,evenfor R,-,> R, . W e have defined the automaton critical radius, R,, as that value of R. where half of the realizations lead to growth and half to shrinkage.W e have carried out such simulations for two valuesof k-, W e find 19<&<20 for k-,=0.0160, and 39
-I

R(t)

01 0

I

I

I

I

500

1000

1500

2000

t

FIG. 9. Radius as a function of time, for a system at k-,=0.0160. The system was prepared with a disk of the more-stable state, of radius R(0) =40, in a sea of the less-stablesteady state. The symbols represent the averageof two realizations of the lattice gas, while the curves are plots of Eq. (15) with different values of R,, as indicated in the figure and discussedin the text.

J. Chem. Phys., Vol. 99, No. 5, 1 September 1993

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

*

3944

Gruner, Kapral, and Lawniczak: Nucleation,

domain growth, and fluctuations , . . *., ! ..

‘, :,;

.,’

. : -.

.

.I

-;.

.

:

,,,.

,I

;.

.:..: ,.

,*

.

2

,,..

,’

,’

-11 .,.. ,,’ i ,,,A,.’ -

.

.

-,. . : ” .’ 1

0' 0

2000

I

I 4000

.--

6000'

“JsLv~: 8000

t

,.:I, ..q$:

FIG. 10. Radius squared as a function of time for a system at equistability (k-, =0.0195), prepared with a disk of radius R(0) = 150 of the highdensity steady state in a sea of the low-density steady state. The symbols in the figure represent the average of six realizations of the automaton, while the straight line represents a linear regression fit. See text for discussion.

condition evolves to form a smooth interface, the large curvature for small R, and the proportionately larger effects of fluctuations for small domains. In the case where the two steady states are equistable (i.e., v=O), the radius is predicted to evolve as k(t) =-D/R(t), which yields R2(t) =R2(0) -2Dt. In Fig. 10 we plot R2( t) vs t for the evolution of a disk of radius R (0) = 150 at equistability. A linear regressionfit of all the points produces a diffusion coefficient DZ 1.46, in very good agreementwith the actual diffusion coefficient used in the simulations ( D = 3/2). Since the probability of formation of a critical nucleus dependsexponentially on the size of the nucleus, except near the edges of the bistable domain where the critical nucleus size is small, a system prepared in the less-stable state will remain there for long periods of time. This is the reason for the observed steady state automaton structure displayed in Fig. 1. It is possible to observe nucleationinduced growth of the more-stable phase in the less-stable phase near the edges of the bistability domain. Figure 11 shows such a nucleation crash where supercritical nuclei form spontaneously and grow to consume the less-stable phase. This simulation was carried out for a parameter value of kml=0.0213, just inside the deterministic bistable domain, whose right edge is at kel=0.021 59. For km1 =0.0213 the critical radius is predicted to be R,z20 cells. This is sufficiently small so that fluctuations are important and the deterministic model breaks down in this region. IV. SUMMARY

The calculations presented in this paper demonstrate that reactive lattice-gas cellular automaton models are a useful tool for the investigation of chemical wave propagation processesin bistable chemical media. Such simulations allow one to study the structure of the propagating chem-

‘, : .

.

.,

,. . . .

‘.,. . .s,,, .’ . .. . : .,‘, ’ ,, I, . ., ., ‘1 ‘. %..,./,.; .z, *!,,, . . ..’ .* ‘:’ ., ,.‘.,..Jk$, ,:. i,r. .J’.,.. \ .. ,. t. -. 800

.:

-

.

,

,.’

:

I .

,,.I.

-

.. 4,. i

,P+q :;

y.;y,-:.: .‘, ,’ .: ,’

s-., ;

I . . -

~

- “,,‘,G . ‘I:., i,$$$ ;. .” ..:

7”.‘,, ‘. ,,‘,

:‘,,y

,““r:‘-,; .

. ’

., .

,’

:,’

,

FIG. 11. Nucleation crash of a system started in the less-stable steadystate, very close to the edge of the deterministic bistable domain. Here k- ,=0.0213, and the initial concentration was that of the deterministic upper stable state. The dark lattice sites correspond to the more-stable state, and the times are indicated in the figure. The lattice elements pictured here are the result of 2 x 2 coarse-graining on the original lattice.

ical wave front and its stability to internal fluctuations. We carried out direct simulations of the relative stability of stable phases in a spatially-distributed medium, and the observed phenomenology was shown to be in accord with the predictions of theory and experiments on an optical bistable device. Since fluctuations are incorporated in the model dynamics, nucleation-induced transitions between deterministic stable phasescould be investigated, and the

J. Chem. Phys., Vol. 99, No. 5, 1 September

1993

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

Gruner, Kapral, and Lawniczak:

Nucleation,

critical nucleus size could be determined and compared with the predictions of the deterministic theory. In purely deterministic models with only autocatalytic production steps, like that in the Schltigl mechanism, some anomalies are observedin the growth of nuclei in some dimensions.” The study of the effects of fluctuations on these irreversible models could prove interesting. Finally, domain formation and growth from the unstable state was examined. The nonequilibrium correlation function and corresponding dynamical structure factor conform to the predictions of scaling theory, confirming that the automaton correctly gives the diffusive dynamics of the curved interfaces separating the two phases,and that a single characteristic length enters in the description of the late-stage dynamics. Thus, the space-time fluctuations in the automaton are given correctly. We note that static correlations have also been investigated for another reactive lattice-gas automaton and yielded results in accord with the predictions of a fluctuating-hydrodynamics model. l8 All of these quantitative results support the use of the reactive lattice-gas cellular automaton method for the study of the reactive dynamics of spatially-distributed systems at the mesoscopiclevel. As such they provide a means to study the emergenceof self-organized structures in the presenceof fluctuations, a feature that is especially important near bifurcations, and when early-stage growth or noise-induced transition phenomena are under investigation. ACKNOWLEDGMENTS

This work was supported in part by a grant funded by the Network of Centers of Excellence program in association with the Natural Sciences and Engineering Research Council of Canada, and a grant from the Natural Sciences

domain growth, and fluctuations

3945

and Engineering Research Council of Canada. We would also like to thank Gary Doolen and the ACL at Los Alamos National Laboratory for making computer time available on the CM-200. ‘J. D. Gunton, M. San Miguel, and P. S. Sahni, in Phase Transitions and Critical Phenomena, edited by C. Domb and J. L. Lebowitz (Academic, New York, 1983), Vol. 8; J. D. Gunton and M. Droz, Introduction to the Theory of Metastable and Unstable States (Springer, New York,

1983). ‘S. M. Allen and J. W. Cahn, Acta Metall. 27, 1085 (1979). 3F. Schliigl, Z. Phys. 253, 147 (1972). 4A. Nitzan, P. Ortoleva, J. Deutch, and J. Ross, J. Chem. Phys. 61, 1056 (1974); A. Nitzan, P. Ortoleva, and J. Ross, Faraday Symp. Chem. Sot. 9, 241 (1974). ‘G. Dewel, D. Walgraef, and P. Borckmans, Z. Phys. B 28,235 (1977); Adv. Chem. Phys. 49, 311 (1982). 6M. E. Brachet and E. Tirapegui, Phys. Lett. A 81, 211 (1981). ‘G. Nicolis and M. Malek-Mansour,. J. Stat. Phys. 22, 495 ( 1980). *P. Grassberger, Z. Phys. B 47, 365 (1982). 9D. 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). loA. Saul and K. Showalter, in Oscillations and Traveling Waves in Chemical Sysfems, edited by R. J. Field and M. Burger (Wiley, New York, 1985), p. 419. “R. Kapral, A. Lawniczak, and P. Masiar, Phys. Rev. Lett. 66, 2539 (1991); J. Chem. Phys. 96,2762 (1992); D. Dab, J. P. Boon, and Y.-X. Li, Phys. Rev. Lett. 66, 2535 (1991). “R. H. Harding and J. Ross, J. Chem. Phys. 92, 1936 (1990); R. H. Harding, P. Paoli, and J. Ross, ibid. 92, 1947 ( 1990). 13F Schliigl and R. S. Berry, Phys. Rev. A 21, 2078 (1980); F. Schlogl, C: Escher, and R. S. Berry, ibid. 27, 2689 (1983). 14A. M. Albano, N. B. Abraham, D. E. Chyba, and M. Martelli, Am. J. Phys. 52, 161 (1983). I5 (a) G.-L. Oppo and R. Kapral, Phys. Rev. A 36, 5820 (1987); (b) R. Kapral and G.-L. Oppo, Physica D 23, 455 (1986). 16A. Lemarchand, A. Lesne, A. Perera, M. Moreau, and M. Mareschal (to be published). “D. J. Needham and J. H. Merkin, Nonlinearity 5, 413 (1992). “J. Weimar, D. Dab, J. P. Boon, and S. Succi, Europhys. Lett. 20, 627 (1992).

J. Chem. Phys., Vol. 99, No. 5, 1 September

1993

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