18976
J. Phys. Chem. 1996, 100, 18976-18985
Stochastic Description of Temporal and Spatial Dynamics of Far-from-Equilibrium Reactions Raymond Kapral* and Xiao-Guang Wu Chemical Physics Theory Group, Department of Chemistry, UniVersity of Toronto, Toronto, Ontario, Canada M5S 1A1 ReceiVed: May 1, 1996; In Final Form: July 2, 1996X
The effects of internal molecular fluctuations on chemical oscillations and chaos are investigated using a lattice-gas model for the dynamics of spatially distributed reacting systems. For limit cycle oscillations, the manner in which temporal correlations develop long-range order as the system size increases is studied. For chaotic attractors the correlation functions indicate the existence of noisy periodic behavior for global concentration variables. These features are explored for the three-variable autocatalator model. The effects of fluctuations on transverse front instabilities in the cubic autocatalysis reaction are also studied. The front instability is investigated in the lattice-gas model, and fluctuations are shown to lead to diffusive roughening of the front below the instability point and, in addition, to modify the cellular structure above the instability point.
I. Introduction Fluctuations in chemical composition are an intrinsic feature of all reacting systems since collisions between molecular species may change their chemical identities. The study of such fluctuations near equilibrium in reactive systems has long played a central role in reaction rate theory. For example, the rate constant of a reaction is determined from the autocorrelation of the flux of the microscopic dynamical variable corresponding to the reaction coordinate.1 The nonequilibrium statistical mechanics of reactive systems close to equilibrium is well developed, and the molecular origin of macroscopic rate laws has been investigated. When a chemically reacting system is constrained to lie far from equilibrium complex attracting system states may arise; for example, limit cycle oscillations, chaotic dynamics, and a variety of spatially inhomogeneous states. The underlying probability distribution that characterizes these states is difficult to determine. As in systems near equilibrium, internal molecular fluctuations and their correlations form the basis of the dynamics that one observes on macroscopic scales. It is known that fluctuations can lead to unusual types of system response in this far-from-equilibrium regime.2 The principal tools that have been used to explore this regime are the master equation,3 in which the reaction is treated as birth and death process where chemical species are created and destroyed in local collision events and molecular motion is described by a random walk process, or other stochastic models.4-6 In this article we make use of a reactive lattice-gas model which implements a birth-death description of the reactive events and a random walk model of the diffusive motion, just as in the master equation approaches to this problem. This discrete space and time Markov chain model permits us to investigate complex system states of far-from-equilibrium reacting systems on large spatial domains. Consequently, we are able to investigate the spatiotemporal dynamics of systems whose mean field models exhibit oscillations or chaos or whose reaction-diffusion models may exhibit instabilities to inhomogeneous perturbations. Section II presents a sketch of the lattice-gas model that is utilized in this study. It is a model where no restriction is X
Abstract published in AdVance ACS Abstracts, November 1, 1996.
S0022-3654(96)01247-6 CCC: $12.00
imposed on the number of molecules of a species at a lattice node. The presentation is given in general terms so that the results may be directly applied to a variety of reactions, in particular the autocatalator and cubic autocatalysis reactions that are the focus of this study. The effects of fluctuations on limit cycle and chaotic oscillations in the autocatalator are the subject of section III. The chaotic attractor arises by a period-doubling cascade, and we examine the correlations of fluctuations in both the periodic and chaotic regimes as a function of system size. The propagation of chemical waves in these regimes is also studied. Transverse front instabilities in the cubic autocatalysis model are studied in section IV using the lattice-gas model. Here we investigate the nature of the propagating chemical front both below and above the instability point and discuss how the front dynamics can be described from a mesoscopic point of view. The conclusions are given in section V. II. Stochastic Lattice-Gas Dynamics A stochastic lattice-gas model will be used to investigate the role that internal molecular fluctuations play in determining the nature of the spatiotemporal dynamics of nontrivial attracting states of far-from-equilibrium reacting systems.7 As in all lattice-gas models of this type, the dynamics takes place on a lattice at discrete time steps. Reactions occur locally at the nodes of the lattice, and diffusion arises from the random walks that the molecules execute on the lattice.8 Since we shall study irreversible reaction schemes, we implement a lattice-gas model where there is no restriction on the number of molecules per node.8-11 We first formulate the model without exclusion in general terms that allows any reaction scheme to be studied; only in the applications will we restrict our considerations to irreversible kinetics. Consider a reaction scheme defined by the mechanism kj
νj‚X y\ z νjj‚X k -j
(1)
where X ) (X1,X2,...,Xn) is a vector of n chemical species, νj and vjj are vectors of stoichiometric coefficients for reactants and products, respectively, and kj and k-j are forward and reverse rate constants. The symbol j (j ) 0, 2, ..., r - 1) labels any of the r different reaction steps. © 1996 American Chemical Society
Temporal and Spatial Dynamics
J. Phys. Chem., Vol. 100, No. 49, 1996 18977
A microscopic reaction dynamics corresponding to this mechanism can be constructed in the following manner: We make space discrete and suppose that the dynamics takes place on some regular lattice.12 Suppose that at a given space point or node of the lattice there are r ) (R1, R2, ..., Rn) molecules of species τ ) 1, 2, ..., n. We place no restriction on the numbers of molecules that may reside at a node; only their average values are controlled by the dynamics.13 Given a configuration of molecules r at a node we wish to determine the probability that a transition to a configuration β ) (β1, β2, ..., βn) will take place as a result of local chemical reactions. We denote this probability by P(r,β). This transition probability may be described by a birth-death process3 as r
[
n
Rκ!
κ)1
(Rκ - νκj )!
P(r,β) ) h ∑ kj ∏ j)1
δβκ,Rκ+νjκj-νκj +
n
Rκ!
κ)1
(Rκ - νjκj)!
k-j ∏
]
δβκ,Rκ-νjκj + νκj (2)
This formula simply reflects the fact that reactions occur in accord with the reaction mechanism with probabilities determined by the rate constants and the numbers of molecules at a node. The delta function δβκ,Rκ+νjκj-νκj ensures that the postreaction configuration of molecules β is consistent with the reaction mechanism. In writing (2), we treat the forward and reverse reaction steps independently so that the prereaction configuration for a reverse reaction step characterized by k-j is r. The net rate of change of the number of molecules of species τ for an initial configuration r as a result of local reactions at a node is
n
r
n
dFτ(t)/dt ) ∑(νjjτ - νjτ) [kj ∏ Fκνκ - k-j ∏ Fκνjκ ] j
j)1
j
κ)1
where the density difference (Fτ(t + h) - Fτ(t))/h has been replaced by the time derivative with t ) hk. In addition to reaction, the molecules execute random walks on the lattice. More specifically, each molecule at a node is assigned a channel which is one of the links of the lattice connecting a node to its nearest neighbors. The particles need not only reside on the links of the lattice and channels, corresponding stop particles which do not propagate on a given time step may also be assigned.8 If the number of molecules at a node exceeds the number of channels, channels are multiply occupied. The channel configuration (collection of channel labels) is randomized, and particles on links are synchronously moved one lattice unit to neighboring nodes connected by the links. In order to be able to independently change the diffusion coefficients of the chemical species, each species may be assumed to reside on its own lattice, and the random walk process may be carried out mτ times for species τ for each reaction step. In the continuous time limit this random walk dynamics generates a diffusion process with diffusion equation
∂Fτ(r,t)/∂t ) Dτ∇2Fτ(r,t)
(9)
and diffusion coefficient Dτ ) mτ/(4h). Alternatively, one may carry out a diffusion step for each mτ reaction steps in which case Dτ ) 1/(mτ4h). Combining the reaction and diffusion steps the macroscopic dynamics, with neglect of correlations and fluctuations, is described by the reaction-diffusion equation r
n
n
∂Fτ(r,t)/∂t ) ∑(νjjτ - νjτ) [kj ∏ Fκνκ - k-j ∏ Fκνjκ ] + j
j
Pτ(r) ) ∑(βτ - Rτ)P(r,β)
(3)
β
(8)
κ)1
j)1
κ)1
κ)1
Dτ∇2Fτ(r,t) (10)
Using the expression (2) for P(r,β), we find In the following section we shall explore the microscopic dynamics generated by the stochastic lattice-gas model and determine the role that fluctuations play in modifying the dynamics predicted by either (8) or (10).
r
Pτ(r) ) h ∑ j)1
(νjjτ
-
[
n
Rκ!
κ-1
(Rκ - νκj !)
νjτ) kj ∏
]
n
Rκ!
κ)1
(Rκ - νjκj )!
- k-j ∏
(4)
Given the probability distribution of the prereaction configuration, p(r;G), which may depend on the mean density G ) (F1, F2, ..., Fn) of the chemical species, the rate of change of the mean density of species τ from discrete time step k to k + 1 is given by
Fτ(k+1) - Fτ(k) ) ∑ Pτ(r) p(r;G)
(5)
r
If diffusion is sufficiently fast compared to reaction, the distribution of molecules at a node will be Poisson, characterized by the local mean density
pτ(Rτ;Fτ) ) (FRτ τ/Rτ!) exp(-Fτ)
(6)
independently for each species, so that p(r;G) is the product distribution, n
p(r;G) ) ∏ pτ(Rτ;Fτ) τ)1
Substitution of (7) into (5) yields the mass action rate law
(7)
III. Oscillations and Chaos in the Autocatalator The study of the stochastic dynamics of systems whose mean field (mass action) rate laws show limit cycle or chaotic oscillations presents interesting questions that have begun to be explored over the past few years. For example, while a phase-space flow generated by the rate law in the limit cycle regime has a negative semidefinite Lyapunov spectrum, chaotic systems are characterized by the existence of at least one positive Lyapunov exponent. The exponential divergence of neighboring trajectories in phase space that follows from the existence of positive Lyapunov exponents has led to the suggestion that mean field models of reacting systems may break down in the chaotic regime.14 Also, it is shown that even in the limit cycle regime, within a period-doubling cascade to chaos, the response of the system to fluctuations can yield noisy attractors that resemble the chaotic attractor.2,15 The nature of global attractors in stochastic, spatially distributed systems with local chaotic dynamics has been discussed.16 It has been established that nontrivial global attractors can exist in spatially distributed systems with chaotic elements.17 Some aspects of these problems have been addressed recently. The validity of mean field rate laws for systems exhibiting deterministic chaos have been investigated using both reactive lattice-gas models18 and master equation approaches19 for the
18978 J. Phys. Chem., Vol. 100, No. 49, 1996
Kapral and Wu
Willamowski-Ro¨ssler model,20 a reversible, three-variable reaction model with quadratic kinetics that exhibits a period doubling cascade to a chaotic attractor. The autocatalor21 has also been used as a model to investigate some aspects of this problem.9 This model is irreversible and involves cubic reaction steps. The macroscopic rate law supports a globally attracting chaotic attractor that again arises by a period-doubling cascade for a suitable range of system parameters. Below we investigate the correlation of fluctuations of the species concentrations in both the limit cycle and chaotic regimes of this model. The autocatalator reaction mechanism is k0
k1
k2
A 98 X1, A + X3 98 X1 + X3, X1 + 2X2 98 3X2 (11) k3
k4
k5
X1 98 X2, X2 98 X3, X3 98 B which leads to the reaction-diffusion equation
∂F1/∂t ) κ0 + κ1F3 - κ2F1F22 - κ3F1 + D1∇2F1 ∂F2/∂t ) κ2F1F22 + κ3F1 - κ4F2 + D2∇2F2
(12)
∂F3/∂t ) κ4F2 - κ5F3 + D3∇2F3 Here the κj are proportional to the kj but incorporate any constant concentrations (A and B) that are fixed by the constraints. We shall consider only the case of equal diffusion coefficients, D1 ) D2 ) D3 ) D. Using (2), we may directly write an expression for the reaction probability matrix in the lattice-gas dynamics. The result is9
P(r,β) ) h{(κ0 + κ1R3)δβ1,R1+1δβ2,R2δβ3,R3 + [κ2R1R2(R2 - 1) + κ3R1]δβ1,R1-1δβ2,R2+1δβ3,R3 + κ4R2δβ1,R1δβ2,R2-1δβ3,R3+1 + κ5R3δβ1,R1δβ2,R2δβ3,R3-1} (13) The simulations were carried out on triangular lattices with N × N sites and periodic boundary conditions. The lattice size ranged from N ) 32 to N ) 512. The system parameters were fixed at κ0 ) 1.666 66, κ1 ) 0.153, κ3 ) 0.02, κ4 ) 4.0, and κ5 ) 1.0. The parameter κ2 was varied to change the state of the system. We consider initial conditions where the distribution of particles at a node is taken to Poisson with some specified average concentration. In the last subsection dealing with chemical waves we consider initial conditions which are inhomogeneous on long spatial scales. A. Dynamics of Global Concentrations. The global concentration is obtained by averaging the local particle numbers over all nodes in the lattice:
Fhτ(t) ) N-2
∑ Rτ(r,t)
(14)
r∈Lτ
where Rτ(r,t) is the number of molecules of species τ at node r at time t. We first consider small lattices where diffusion is effective in maintaining homogeneity over the entire lattice, and fluctuations play a significant role in the dynamics. The effects of fluctuations on the correlations of the global concentrations are determined by the structure of the temporal correlation function,
Cτ(t) ) 〈δFhτ(t) δFhτ(0)〉/〈(δFhτ(0))2〉
(15)
δFhτ(t) ) Fhτ(t) - 〈Fhτ(t)〉
(16)
where
and the angle brackets signify the average T
〈A〉 ) T-1 ∑ A(t′)
(17)
t′)1
We consider two system sizes, lattices with N ) 32 and N ) 64, and three parameter values, κ2 ) 12.5, 11.9, and 11.52, corresponding to deterministic period-1, period-2, and chaotic attractors, respectively. The calculations employ a time scale factor of h ) 0.001, and temporal averages were carried out over 2 × 106 time steps or T ) 2 × 103. The period of the period-1 limit cycle is tp1 ∼ 2.5, so the averages used in computing the correlation functions are over about 800 cycles of the oscillation. We have taken mτ ) 2, which gives a rough estimate of the diffusion length of lD ) (Dtc)1/2 ∼ 35 lattice sites. In fact, simulation shows that homogeneity is maintained over domains of more than twice this size. The correlation functions are shown in Figure 1. For each of the three attractors the correlation function for the deterministic dynamics is compared with that for the stochastic dynamics with N ) 32 and N ) 64. All of the correlation functions that we display are for the X1 species, C1(t), and we henceforth dispense with the subscript. The results in panel a pertain to the period-1 attractor. One observes that temporal coherence is almost completely destroyed in the N ) 32 system with significant structure persisting only to times t ∼ 6 where the oscillatory negative correlations disappear. In contrast, for the N ) 64 system there is a long-lived oscillatory structure which persists to about t ∼ 20. The results for period-2 display different behavior (panel b). The infinite-range, period-2 structure of the correlation function of the deterministic global concentration is substantially modified for both the N ) 32 and N ) 64 system sizes, which exhibit quite similar behavior. Correlations decay in an oscillatory fashion with oscillation period that is several times the cycle time on the attractor. The shorter time scale cycle time oscillations are superimposed on this longer time scale oscillatory decay. As expected, in the chaotic regime even the correlation functions for the deterministic dynamics exhibit decaying correlations with irregular structure (panel c). Once again the results for both the N ) 32 and N ) 64 system sizes are in close accord. Although the decay is somewhat more rapid than that in period-2 regime, the correlations for the stochastic dynamics present similar structure in both the period-2 and chaotic regimes. These results may be understood when the structures of the noisy and deterministic attractors are compared (cf. Figure 2). The period-1 regime was selected to lie far from the chaotic attractor while, of necessity, the period-2 attractor lies close to the chaotic regime since the period-doubling cascade comprises a small region of the κ2 parameter space (∼11.8 < κ2 < 12.0). As a consequence, internal noise arising from fluctuations of the chemical concentrations causes the system to explore the underlying manifold structure that gives rise to the chaotic attractor. The unstable fixed point of the deterministic system has stable and unstable manifolds. Fluctuations are very effective in causing the system to explore the unstable directions in phase space so that the noisy dynamics resembles that of the nearby deterministic chaotic attractor. Thus, noise causes the probability density to spread in phase space in a structured way so that the noisy period-2 attractor is geometrically similar to the chaotic attractor that appears at larger parameter values in the deterministic flow (cf. Figure 2e,f). However, the details of the probability density may differ in the stochastic and deterministic cases. For these small system sizes the fluctuation effects are quite striking. We have demonstrated for the wellstirred and spatially distributed Willamowski-Ro¨ssler model
Temporal and Spatial Dynamics
J. Phys. Chem., Vol. 100, No. 49, 1996 18979
Figure 2. Attractors in the period-1 ((a) and (b)), period-2 ((c) and (d)), and chaotic ((e) and (f)) regimes. In each set of panels the left attractor is for the deterministic system while the right attractor is for stochastic dynamics with N ) 64. All attractors are shown in projection on the (Fh1Fh2) plane, and both the ordinate and abscissa range from 0 to 6 for these variables. For clarity, only about 50 cycles of the trajectories are plotted.
Figure 1. Correlation functions for the period-1 (a), period-2 (b), and chaotic (c) regimes. In each panel the solid line corresponds to the results for the deterministic system while the heavy dashed line is for N ) 32 and the light dashed line is for N ) 64.
that as the system size increases the noisy attractor approaches that of the deterministic system. In this case the manner in which the system size must be increased in order to resolve additional periodic orbits in the period-doubling cascade has also been determined.18 Even for fairly small systems containing about ∼108 particles fluctuation effects on the global attractor in well-stirred systems are sufficiently small that external noise is likely responsible for the inability to resolve higher periodic orbits in experiments. Another measure of correlations in the system is provided by the variance of the concentration. Provided diffusion is sufficiently rapid the particle distribution at a node will be Poisson. This fact was used to determine the mean field rate law from the automaton dynamics. For a Poisson distribution the variance is equal to the mean. In Figure 3 the difference
Figure 3. Plot of V, the difference between the variance and the mean, as a function of time for the period-1 attractor at κ2 ) 12.5 on a 64 × 64 lattice.
between the variance and the mean for the global concentrations, V ) 〈(δFh1)2〉 - 〈Fh1〉, is plotted as a function of time for period-1 dynamics on a 64 × 64 lattice. The quantity V has small fluctuations around zero for the majority of the time but does show deviations with varying magnitude at periodically spaced
18980 J. Phys. Chem., Vol. 100, No. 49, 1996
Kapral and Wu
Figure 4. Correlation functions C(c) j) (solid line) and C(t) (dashed 1 (t;r line) in the (a) period-1 regime at κ2 ) 12.5 and (b) chaotic regime at κ2 ) 11.52. C(c)(t;rj) is shown for one value of jr lying at the center of a coarse-grained cell with L ) 32 on a 128 × 128 lattice.
time intervals. As is seen in panel a of Figure 2 the limit cycle oscillation is relaxational in character and during the fast portion of the cycle the the concentration changes rapidly. The correlations produced by the reaction are responsible for the deviations seen in the figure. B. Coarse-Grained Dynamics. We now consider the dynamics on larger lattices which we coarse grain into smaller cells with linear dimension L. The coarse-grained concentration field of species τ is defined as
Fhcτ(rj,t) ) L-2
∑
Rτ(r,t)
(18)
r∈Lτ(rj)
where Lτ(rj) is a domain of the lattice centered on jr containing L × L nodes. As in the preceding section, we may examine the decay of correlations of the coarse-grained density fields. We denote these correlation functions by a superscript “c” and argument jr to indicate that coarse-grained variables in cell jr are used in the definition (15): C(c)(t;rj). Figure 4 presents results for the correlation functions on lattices with N ) 128 which have been coarse grained into 16 cells with L ) 32. In each panel of this figure the correlation function for one value of jr is compared with that for the global concentration of X1 obtained from an average over the entire 128 × 128 lattice. The correlation functions for all values of jr show similar structure. Consider the results in panel a for κ2 ) 12.5 in the period-1 regime. Both C(c) and C exhibit much longer temporal coherence (in excess of 70 cycles of the oscillation) than that for simulations on the smaller 32 × 32 or 64 × 64 systems as can be seen by comparison with Figure 1. We also see that the correlations for the global concentration have smaller amplitude while displaying similar temporal coherence. Two effects come into play in these larger systems.
Roughly, the internal noise in the system scales as N-1 so the noisy attractor of the global concentration on the larger lattice should more closely approach that of the underlying period-1 oscillation. From a comparison of the results for N ) 32, 64, and 128 one sees how the oscillatory correlations grow with system size and approach the infinite-range correlations of the deterministic dynamics. The other effect that must be considered is the destruction of phase coherence on longer distance scales as a result of finite diffusion. This is signaled by the fact that the correlation function in cell jr has larger amplitude and consequently persists for longer times than that for the entire system. Although diffusion is able to maintain homogeneity over large spatial domains, concentration gradients do arise in the system, especially during the fast portions of the oscillation.9 As a result of partial phase averaging the global attractor occupies a smaller region of phase space, and consequently the correlation function of this global concentration has smaller oscillation amplitude. In the limit of a very large system the amplitude will shrink to zero. The fact that the temporal coherence in the 32 × 32 coarse-grained cells in the 128 × 128 system is substantially longer than that for an isolated 32 × 32 system implies that diffusion is able to maintain homogeneity over spatial domains larger than L ) 32. According to the scaling arguments for homogeneous systems where D is sufficiently large, the linear system size must be increased by approximately a factor of 7 in order to be able to resolve the structure at the next period doubling, period-2 in this case.18 The 128 × 128 lattice is too small to be able to resolve this structure, and as we have observed, the stochastic dynamics in the period-2 regime has many characteristics of the chaotic attractor.22 We next consider the correlations in the chaotic regime for κ2 ) 11.52 (cf. panel b of Figure 4). Several features of the correlations deserve comment. The correlations are much more highly structured than those for the smaller N ) 32 and N ) 64 systems: the incipient short-period oscillations that can be detected in the small-system correlation functions are very prominent in the N ) 128 lattice results. The appearance of nearly regular oscillations in the correlation function arises from the fact that in the chaotic regime noise has a large effect on the amplitudes of the chaotic osillations but is less effective in destroying the phase coherence that arises from cycling around the unstable fixed point that spawned the chaotic attractor. Thus, when the global or coarse-grained concentrations are considered, the amplitude fluctuations average to some nonzero value and the phase space trajectories have the appearance of a noisy limit cycle. For the small system sizes (N ) 32, 64) fluctuations are sufficiently strong to destroy much of the phase coherence, but already for N ) 128 the effects of fluctuations are more subtle and primarily affect the unstable directions of the flow which lead to the above-mentioned strong amplitude fluctuations with little loss of phase coherence. If the system size is increased, phase coherence will be lost as well, and the globally averaged attractor will appear as a noisy fixed point when plotted in phase space. Consequently, as the system size increases from small to large values, the global or coarse-grained attractors take on regular a limit-cycle-like structure which is eventually destroyed by loss of phase coherence in the system. We observe that the results for the global and coarse-grained correlations in Figure 4 are indistinguishable in contrast to the period-1 results for the same system size. The strong local instability due to the positive Lyapunov exponent leads to amplitude randomization on spatial scales smaller than L so that the coarse-grained cells are dynamically independent as far as this aspect of the dynamics is concerned, but phase coherence is largely maintained.9
Temporal and Spatial Dynamics
J. Phys. Chem., Vol. 100, No. 49, 1996 18981
Figure 5. Evolution of a circular wave in the period-2 regime at κ2 ) 11.9 on a 512 × 512 lattice. The concentration of X2 is shown as gray shades with high concentration black and low concentration white. The time scale parameter is h ) 0.0005, and moving from left to right and top to bottom, the frames are for k ) 1700, 1900, 2200, and 2400 time steps, respectively.
C. Chemical Waves. Thus far we have considered initial conditions which were spatially homogeneous apart from small uniformly distributed fluctuations. For the smallest system sizes studied above diffusion is able to maintain spatial homogeneity while for larger system sizes concentration gradients are generated by the dynamics leading to the destruction of phase coherence. Here we examine how chemical waves propagate in the fluctuating medium. A chemical wave may be initiated by introducing a diskshaped domain in the center of the lattice where the concentration field corresponds to a point on the limit cycle shifted in phase by approximately π from that of the uniform “sea” in which it resides. Depending on the character of the oscillation and its amplitude, the local phase difference will either be smoothed by diffusion or produce a propagating wave. An example of such wave propagation is shown in Figure 5. The system is in the period-2 regime and the lattice size 512 × 512. The evolution is shown for four different times covering a fraction of a global oscillation cycle. During this time period one observes propagation of the chemical wave as expected. However, already in the last two frames one sees spontaneous fluctuations which grow and give rise to local gradients in the system which interfere with the propagation of the circular wave. For later times the system is completely disorganized with the local concentrations varing erratically in space and time. This behavior is expected in view of the fact that period 2 lies near chaos and is influenced by fluctuations as described above. However, the system size is now large enough so that if the system were well stirred the period-2 orbit could be resolved. Thus, for periodic orbits near the chaotic regime local fluctua-
tions can have important effects on the spatiotemporal dynamics and wave propagation in small systems or systems with small diffusion coefficients. IV. Cubic Autocatalysis Fronts Propagating reaction-diffusion fronts can undergo transverse instabilities when the spatial dimension is greater than one. Such front instabilities are controlled by the reaction kinetics and the relative values of the diffusion coefficients of the reacting species. The cubic autocatalytic reaction
X1 + 2X2 f 3X2
(19)
has been studied extensively in this context.23-27 The corresponding reaction-diffusion equation is
∂F1/∂t ) -F1F22 + D1∇2F1 ∂F2/∂t ) F1F22 + D2∇2F2
(20)
On the basis of simple physical considerations, one expects that increasing the diffusion coefficient of the fuel X1 over that of the autocatalyst X2 will cause a planar front to become unstable provided the ratio δ ) D1/D2 is greter than a critical value δc. Consider a two-dimensional system infinitely extended along x with length Ny sites along y. Suppose species X1 initially fills the left half-plane and X2 the right half-plane. A chemical front will propagate from right to left as X2 consumes X1. Let φ(y,t) denote the profile of the front at time t. There have been a number of studies of the nature of the transverse front
18982 J. Phys. Chem., Vol. 100, No. 49, 1996
Kapral and Wu
Figure 6. Comparison of the computed front profile (diamond and long dashed line) for a system with Ny ) 64 with (23) (solid line) for δ ) 1. The computed results for a system with Ny ) 512 are shown as crosses and a short dashed line. System parameters are h ) 0.02 and D ) (0.08)-1.
analysis26 of (20). In addition, the nature of the front instability far beyond onset has been studied, and a type of biscale chaotic front dynamics has been described.27 In this section we consider a lattice-gas model of the front propagation in two spatial dimensions. The reaction probability matrix is trivial for this case and can be found from a specialization of (2) to the single irreversible reaction step (19). The result is
P(r,β) ) hR1R2(R2 - 1)δβ1,R1-1δβ2,R2+1
Figure 7. Front is shown for two times for a system with Ny ) 512. The concentration of X2 is shown as gray shades with high concentration black and low concentration white. The front propagates from right to left.
instability and the form of the profile, especially in two spatial dimensions. Close to the instability point the cubic autocatalysis reaction-diffusion equation may be reduced27 to the KuramotoSivashinsky equation28
∂φ(y,t) c ) ν∇2φ + |∇φ|2 + κ∇4φ ∂t 2
(21)
where c is the front speed and ν and κ are parameters that may be determined from (20). The front instability if signaled by a change in sign of ν < 0 giving rise to instability in the planar front. The critical value of δ has been determined to be δc ) 2.300 from a detailed analysis of the front dynamics in terms of the Kuramoto-Sivashinsky equation27 and using an eikonal
(22)
One may easily verify that this reaction probability matrix leads to the reaction-diffusion equation (20) using the method outlined in section 2. With this model we may examine the effects of local fluctuations arising from reaction and random particle motion on the front dynamics and bifurcations. The simulation results were obtained by considering the evolution on triangular lattices with Nx × Ny sites. Periodic boundary conditions were used along y, and dynamics was adjusted to simulate an infinite domain along x. The average front position was monitored and shifted backward by the number of units it had propagated along x. This was accomplished by deleting the appropriate number of rows from the end of the system behind the front and adding the same number to the beginning of the system ahead of the front. These newly added rows had X1 concentration corresponding to the initial state, F*1. In this way the propagation can be simulated for indefinite lengths of time. Diffusion was modeled as described in the previous section. For this model we have D ) 1/(4h), and D1 and D2 are multiples of this value. A. Front Dynamics below δc. Of course, in 1D no transverse instability is possible, but rather detailed studies of the dependence of the front velocity on the diffusion ratio have been carried out.23 In addition, the 1D dynamics forms the basis of analytical studies of the front in higher spatial dimensions. The reaction-diffusion equation (20) with equal diffusion coefficients (D1 ) D2 ) D) may be solved in 1D for the interfacial profile and front speed. The results are
F1(ξ) ) F* (1 + eVξ/D)-1
(23)
Temporal and Spatial Dynamics
J. Phys. Chem., Vol. 100, No. 49, 1996 18983
Figure 8. Front is shown for two times for a system with Ny ) 512 for two different system parameters both with δ ) 8: (a) h ) 0.04 and D1 ) 2D and D2 ) D/4 and (b) h ) 0.08 and D1 ) 4D and D2 ) D/2. Time increases from left to right. The coding is the same as in Figure 7. The double-headed arrow indicates the cellular length l*.
with V ) chF* (D)1/2 where limξf-∞ F1(ξ) ) F* and c ) 1/x2. The variable ξ describes displacements in the moving frame ξ ) x - Vk, where again k is the discrete time. From mass conservation it follows that F2(ξ) ) F* - F1(ξ). We have written these equations in units appropriate for comparison with the automation simulation. The relation between the automaton space and time scales and those in (20) are t ) k(h(F*)2) and r ) (D1/(F*)2)-1/2x where r and x are the continuous and discrete space variables, respectively. We have carried out 2D simulations of the cubic autocatalysis model for this equal diffusion case. Since no transverse instability occurs the propagating front will be planar if determined from the reaction-diffusion equation (20). In the lattice-gas model the front profile may be obtained from an average of the node particle numbers in the y direction. The
results of such a calculation for a system with small width (Ny ) 64) are shown in Figure 6 and compared with (23). The two results are in close accord as is the measured front speed (V ) 0.15 or c ∼ 1/x2. When Ny is larger fluctuations manifest themselves and give rise to transverse structure to the interface. For equal diffusion coefficients one has ν ) D and κ ) 0 for the KuramotoSivashinsky parameters. Internal fluctuations provide a source of noise, and we may account for this fact phenomenologically by appending a noise term η to the resulting equation:
∂φ c ) D∇2φ + |∇φ|2 + η ∂t 2
(24)
where η may be regarded as a Gaussian random variable with
18984 J. Phys. Chem., Vol. 100, No. 49, 1996
Kapral and Wu
Figure 9. Space-time plot of the front profile. Same parameters as Figure 8a.
Figure 10. Plot of the front velocity versus δ.
white noise spectrum. This is the Kardar-Parisi-Zhang (KPZ) equation.29,30 When there are no gradients, the front may be described by the Edwards-Wilkinson equation31 in which the nonlinear term is dropped. The noise induces roughness in the interfacial profile in the absence of instability that leads to an interfacial thickness, measured as deviations of the front from the mean interface position, that scales with Ny1/2. The interface for Ny ) 512 is shown in Figure 7 for two times and clearly displays the effects of diffusive roughening in this stable front regime. The profile obtained by averaging the concentration over y is shown in Figure 6. The thick interfacial profile for the large system length reflects the transverse structure of the front. B. Front Dynamics above δc. When δ exceeds δc, the planar front loses its stability, and a chaotic cellular front is
formed provided the system length Ny is large enough. As noted above, the instability is signaled by a change in sign of ν (ν < 0) so that long wavelength fluctuations of the front are unstable but short wavelength fluctuations are stabilized since κ > 0. The balance between these two effects leads to a characteristic length l* ) 2π/km, where km ) (-ν/2κ)1/2. Simulations of the front propagation have been performed for a variety of values of δ, and the front instability has been observed. By way of illustration the results for β ) 8 are shown in Figure 8a. This value of δ is far beyond the instability point where biscale chaos was observed in the reaction-diffusion simulations. However, the lattice-gas simulations were carried out for short times on relatively small systems (Ny ) 512) where the long-scale structure has not yet fully developed. The front instability is clearly evident in this figure (compare Figures 7 and 8a). The characteristic cellular length l* determined numerically from the reaction-diffusion equation simulation is indicated in this figure by a double-headed arrow. One can see the lattice-gas cellular structure on this scale. Figure 9 shows a space-time plot of the front profile extracted from the automaton simulations. The evolution of the cellular structure can be seen; for example, extrema in the profile coalesce and new extrema are born. Referring again to Figure 8a, one can see structure on scales small compared to l* due to fluctuations in the species concentrations. Fluctuations therefore modify the fine structure of the front and when these fluctuations are sufficiently strong they can lead to a destruction of the cellular pattern. An example of this is shown in Figure 8b where fluctuations have been enhanced by increasing the reaction rate. The cellular structure of the front is far less evident. The mean speed of the front was measured as a function of δ in the two-dimensional system where the front instability occurs. The results are presented in Figure 10. The front velocity shows the expected behavior, namely, that the velocity decreases as δ decreases. For a one-dimensional system Billingham and Needham23 have given analytical estimates of the dependence of the front velocity on D1 and D2. For δ . 1 and fixed D2 (this is how δ was varied in the simulations for δ
Temporal and Spatial Dynamics > 1) one has V ∼ δ-1/2. The velocity for the largest values of δ scales in this way, but the ratio of velocities approaches but does not equal the square root of the ratio of δ values. Of course, our calculations pertain to the 2D front with transverse structure and not to the 1D case for which the theoretical estimates were made. V. Conclusions The two case studies of how internal molecular fluctuations influence the dynamics of far-from-equilibrium reacting systems illustrated different features of the dynamics. In the autocatalator model the focus was on the oscillatory dynamics and the period doubling bifurcations leading to chaos. In the limit cycle regime, far removed from chaos, it was demonstrated how the correlations build long-range order as the system size increases leading, eventually, to temporal correlations on macroscopic scales. In the chaotic regime, or near to it, fluctuations have a more subtle effect. For example, in the period-2 regime the noisy attractor resembles the deterministic chaotic attractor for small system sizes. In the chaotic regime or near it, as the system size increases instability leads to random variations of the amplitude, but phase coherence is maintained. As a result the noisy chaotic attractor actually has the appearance of a noisy periodic orbit in phase space which is manifested in the oscillatory character of the correlation function. This effect of noise on the period doubling sequence was described earlier for stochastically forced one-dimensional maps16 and was observed in the autocatalator simulations.9 As the system size increases even further, phase coherence is destroyed in both the periodic and chaotic regimes. The second example considered the effects of fluctuations on a spatiotemporal front bifurcation in the cubic autocatalysis model. Here it was observed that fluctuations play a role, albeit a fairly obvious one, even when there is no transverse front instability: internal fluctuations provide a source of noise that can lead to a diffusively rough interface on long scales in the two-dimensional simulations. The front instability was observed when the diffusion coefficient of the fuel was sufficiently larger than that of the autocatalyst. Fluctuations were seen to modify the small scale structure of the chaotic cellular fronts and, for noise with large enough amplitude, to destroy or significantly modify the cellular structure. Mesoscopic studies of chemical instabilities and complex nonequilibrium dynamics provide a statistical mechanical basis that underlies commonly used reaction-diffusion descriptions. As such, they allow one to investigate the range of validity of macroscopic descriptions and study phenomena such as critical behavior near bifurcation points whose existence depends on fluctuations. Acknowledgment. This work was supported in part by a grant from the Natural Sciences and Engineering Research Council of Canada and by a Killam Research Fellowship (R.K.). References and Notes (1) Yamamoto, T. J. Chem. Phys. 1960, 33, 281. (2) Crutchfield, J. P.; Farmer, J. D.; Huberman, B. A. Phys. Rep. 1982, 92, 45. Kapral, R.; Schell, M.; Fraser, S. J. Phys. Chem. 1982, 86, 2205. (3) Gardiner, C. W. Handbook of Stochastic Processes; SpringerVerlag: New York, 1985. Nicolis, G.; Prigogine, I. Self-Organization in Non-Equilibrium Systems; Wiley: New York, 1977. van Kampen, N. G. Stochastic Processes in Physics and Chemistry; North-Holland: Amsterdam, 1981.
J. Phys. Chem., Vol. 100, No. 49, 1996 18985 (4) Examples of such studies using Fokker-Planck or Langevin models may be found in: Dewel, G.; Walgraef, D.; Borckmans, P. AdV. Chem. Phys. 1982, 49, 311. Nitzan, A.; Ortoleva, P.; Deutch, J.; Ross, J. J. Chem. Phys. 1974, 61, 1056. (5) Thermodynamic and stochastic formulations have been considered by: Peng, B.; Hunt, K. C.; Hunt, P. M.; Suarez, A.; Ross, J. J. Chem. Phys. 1995, 102, 4548. (6) It is also possible to investigate such questions using molecular dynamics. This necessitates working in the full position-velocity phase space of the system; the computations are computationally intensive, but progress in this direction is possible. See for instance: Ortoleva, P.; Yip, S. J. Chem. Phys. 1976, 65, 2045. Boissonade, J. Phys. Lett. A 1979, 74, 285. Gorecki, J. In Far-From-Equilibrium Dynamics of Chemical Systems; Gorecki, J., Cukrowski, A. S., Kawczynski, A. L., Nowakowski, B., Eds.; World Scientific: River Edge, NJ, 1994. (7) Such models may also be used to explore fluctuation effects near equilibrium where the detailed balance condition guarantees decay to equilibrium; however, our focus is on systems maintained far from equilibrium by flows of reagents so that detailed balance no longer applies, and complex attracting states may be obtained. (8) For a review of reactive lattice-gas models see: Boon, J.-P.; Dab, D.; Kapral, R.; Lawniczak, A. Phys. Rep. 1996, 55, 273. (9) Kapral, R.; Wu, X.-G. Physica D, in press. (10) Chopard, B.; Droz, M.; Cornell, S.; Frachebourg, L. In Cellular Automata for Astrophysics; Perdang, J., Lejeune, A., Eds.; World Scientific: Singapore, 1994. (11) Karapiperis, T.; Blankleider, B. Physica D 1994, 78, 30. (12) The restriction to regular lattices is not essential, and the dynamics on fractal lattices may be investigated. (13) Typically, in simulations the average number of molecules at a node is selected to be less than the coordination number of the lattice, or the number of channels, so that multiple occupancy of a channel rarely occurs. (14) Fox, R. F. Phys. ReV. A 1990, 41, 2960; 1990, 42, 1946. Fox, R. F.; Keizer, J. Phys. ReV. Lett. 1990, 64, 249; Phys. ReV. A 1991, 43, 1709. Keizer, J. E.; Fox, R. F. Phys. ReV. A 1992, 46, 3572. Keizer, J.; Tilden, J. J. Phys. Chem. 1989, 93, 2811. (15) Celarier, E.; Kapral, R. J. Chem. Phys. 1987, 86, 3366. (16) Bohr, T.; Grinstein, G.; He, Y.; Jayaprakash, C. Phys. ReV. Lett. 1987, 58, 2155. (17) Chate´, H.; Manneville, P. Europhys. Lett. 1992, 17, 291; Prog. Theor. Phys. 1992, 87, 1. (18) Wu, X.-G.; Kapral, R.; Wu, X.-G. Ann. N.Y. Acad. Sci. 1993, 706, 186; Phys. ReV. Lett. 1993, 70, 1940; J. Chem. Phys. 1994, 100, 5936; Phys. ReV. E 1994, 50, 3560. Kapral, R.; Wu, X.-G. In Chemical WaVes and Patterns; Kapral, R., Showalter, K., Eds.; Kluwer: Dordrecht, 1995; p 609. (19) Geysermans, P.; Nicolis, G. J. Chem. Phys. 1994, 99, 8964. Nicolis, G.; Baras, F.; Geysermans, P.; Peeters, P. In Chemical WaVes and Patterns; Kapral, R., Showalter, K., Eds.; Kluwer: Dordrecht, 1995; p 573. Baras, F.; Geysermans, P. NuoVo Cimento 1995, 17D, 709. (20) Willamowski, K.-D.; Ro¨ssler, O. E. Z. Naturforsch. A 1980, 35, 317. (21) Peng, B.; Scott, S. K.; Showalter, K. J. Phys. Chem. 1990, 94, 5243. (22) For the Willamowski-Ro¨ssler model the period-2 structure is resolved with N ) 300. (23) Billingham, J.; Needham, D. J. Philos. Trans. R. Soc. London, A 1991, 334, 1. (24) Scott, S. K.; Showalter, K. J. Phys. Chem. 1992, 96, 8702. Horvath, D.; Petrov, V.; Scott, S. K.; Showalter, K. J. Chem. Phys. 1993, 98, 6332. (25) Zhang, Z.; Falle, S. A. E. G. Proc. R. Soc. London, A 1994, 446, 1. (26) Milton, R. A.; Scott, S. K. J. Chem. Phys. 1995, 102, 5271. (27) Malevanets, A.; Kapral, R. Phys. ReV. E 1995, 52, 4724. (28) Kuramoto, Y.; Tsuzuki, T. Prog. Theor. Phys. 1976, 55, 356. Sivashinsky, G. I. Combust. Sci. Technol. 1977, 15, 137. (29) Kardar, M.; Parisi, G.; Zhang, Y.-C. Phys. ReV. Lett. 1986, 56, 889. (30) The deterministic Kuramoto-Sivashinsky equation may be reduced in 1 + 1 dimensions to the KPZ equation by projecting onto the long wavelength modes. This procedure renormalizes ν and gives rise to a noise term that accounts for the short wavelength structure [Yakhot, V. Phys. ReV. A 1981, 24, 642]. The source of noise in the present case arises from internal fluctuations in a region where there is no transverse instability. (31) Edwards, S. F.; Wilkinson, D. Proc. R. Soc. London, A 1982, 381, 17.
JP961247E
