Interaction

of Turing and flow-induced

chemical instabilities

S. Ponce Dawson Theoretical Division and Center for Nonlinear Studies,Los Alamos National Laboratory, Los Alamos, New Mexico 87545

A. Lawniczak Department of Mathematics and Statistics, University of Guelph, Guelph, Ontario NIG 2WI, Canada and Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, New Mexico 87545

R. Kapral Chemical Physics Theory Group, Qepartment of Chemistq University of Toronto, Toronto, Ontario M5S IAI, Canada (Received 23 September 1993; accepted 14 December 1993) The interaction between the Turing instability and the instability induced by a differential flow is studied in the Selkov model. Both instabilities give rise to the formation of spatial patterns, and for a range of parameter values, these patterns can compete. The effect of anisotropic diffusion on the pattern formation process is investigated. Stripes with different orientations that travel with time and the suppression of patterns due to a competition of both instabilities are observed.

1. INTRODUCTION The interest in the study of pattern formation in reaction-diffusion systems or, more generally, in reactiontransport systems, has been increasing steadily since the seminal paper by Turing.’ In this work, Turing shows that a homogeneous steady state of a reaction-diffusion system may become unstable to inhomogeneous perturbations due to the interaction of the chemical kinetics and diffusion. As a result of this interaction, the system may evolve into an inhomogeneous steady state characterized by the presence of a spatial pattern. The simplest case in which this instability, called the Turing instability, takes place satisfies the following two conditions: (i) the kinetic system consists of two species-an activator (e.g. autocatalytic species) and an inhibitor-whose function is to prevent the autocatalytic growth; and (ii) the diffusion coefficient of the inhibitor is sufficiently greater than that of the activator. In this way, the activator can grow locally due to some random perturbation, but the lateral spread of its growth is prevented by the inhibitor. We call the first condition the kinetic condition and the second the transport condition. The role of the latter is to spatially “disengage” the reacting species. In the case of isotropic diffusion, the reaction-diffusion equations can be written generally as dP1 -D,V’p,=R,, at

IslsM,

(1)

where t is time and V2 is the Laplacian operator with respect to the spatial coordinate x. The number of species is M and pt(x,t) is the mass density of the species 1 at time t and position x. The diffusion coefficient is Dt , which we assume to be independent of any other variables of the problem. In the case of an inhibitor-activator system, M=2 and the following conditions are assumed to be satisfied: JRAIdpA>O; dR,ldpt
respectively. As noted before, the transport condition for the occurrence of the Turing instability is D,>D, . The formation of Turing patterns has been observed numerically and experimentally.’ It is also believed that the Turing instability might be one of the basic mechanisms of pattern formation in a variety of biological systems.3*4 In recent papers,5-7 it has been shown that the species can be spatially disengaged in another way, namely, by a differential bulk flow. The flow of the reactants at different flow rates, regardless of which is faster, can destabilize the homogeneous steady state of the system and lead to the formation of traveling-wave type patterns. This instability, called the differential-jlow-induced chemical instability (DIFICI), was studied analytically5*7 and numerically in one spatial dimension showing the formation of patterns.’ The patterns have also been observed in experiments.6 The analysis in Ref. 5 shows that in order for the DIFICI to occur in a two-species system, the kinetic condition is the same as for the Turing instability, namely, the system must be composed of an activator and an inhibitor. In systems involving more than two species, this condition is not necessary. The simplest system that can exhibit the differential-flow-induced chemical instability is described by the equations JPI -g= +A dtfv

R I? @A

- dX

-DAV2pA=RA,

where the inhibitor is completely immobilized and the activator flows with velocity u along the x direction and diffuses. This is the situation that was analyzed in Ref. 5. We study in this paper the interaction between the Turing instability and the differential-flow-induced chemical instability in two spatial dimensions for a system that flows through a porous medium and that can be formally described as a two-species system. We consider the general case in which diffusion is anisotropic. The introduction of this anisotropy allows, in certain cases, a natural separation of both

0 1994 American Institute of Physics 5211 00(7)/521 J. Chem. Phys. 100 April to 1994 Downloaded 21 (7), Apr12003 142.150.192.30.0021-9606/94/l Redistribution subjectl/8/$6.00 to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

Ponce Dawson, Lawniczak, and Kapral: Turing and flow-induced instabilities

5212

instabilities along two different spatial directions. Many real systems, especially biological ones, are anisotropic; thus, the study of the anisotropic case is more general and it can describe a wider variety of real situations. The organization of the paper is as follows: In Sec. II, we show how to obtain the basic set of equations starting from the standard hydrodynamic equations for a multicomponent reacting system. In Sec. III, we analyze the linear stability of the homogeneous solution in the presence of both unequal diffusion coefficients and differential flow. This analysis is relevant for an infinite system or for periodic boundary conditions. In Sec. IV, we present the results of a set of numerical simulations that illustrate the interaction of the Turing and the DIFICI instabilities in different situations. In Sec. V, we repeat the linear stability analysis for another choice of boundary conditions. Finally, the conclusions are summarized in Sec. VI.

II. REACTION-DIFFUSION-ADVECTION

EQUATIONS

We show in this section how to go from the hydrodynamic equations to a set of two coupled reaction-diffusion equations with flow. Rather than considering the most general situation, we focus on the case where a set of chemically reacting species, some of which undergo complex formation, diffuse and are advected by a flow field. Systems of this type are common in nature and the complex formation mechanism is believed to play an important role in Turing pattern formation in real systems.8 Consider a fluid composed of different reactants X, Y, S, and P, in a solvent, such that S reacts only with X to give the “inert” product P according to the law HX4P.

(4)

k-c

We assume that X, Y, and the solvent flow with the mean velocity u through a porous medium, and that the solvent is in excess, so that psol~>p1+p2, with pr and pz the densities of X and Y, respectively, and psol the density of the solvent. We also assume that the flow is incompressible, so that V*u=O and [d(p,,+pl+p2)]/(dt)~(dp,,)l(dt)=0, and that both S and P are immobile with respect to the porous medium. This last assumption together with Eq. (4) implies that the sum of the number densities of S and P, ns + n , is constant. Then, the evolution of the system can be descr&ed by equations for the densities p1 , p2, and pp and by the Navier-Stokes equations for the mean velocity u.~ Given that the moving species flow through a porous medium, we can replace the NavierStokes equation by Darcy’s law,” which states that the velocity u is proportional to the pressure gradient as long as this gradient is not too large. Darcy’s law can be written as

-vp=y u,

where i is the unit vector along the x direction and u is a constant, provided that p is practically uniform and constant. We concentrate now on a geometry with two spatial dimensions. Assuming that the medium in which the species diffuse is, in general, anisotropic, and that the anisotropy that each species feels is different, the basic equations can be written as JPl

dt

where p is the pressure, p=p1+p2+psol=psol, K is the hydrodynamic permeability, and Y is the kinematic viscosity. Thus, given a uniform gradient of pressure along some direction, which we can pick to be x, we can then assume that u=ui,

-Dlx

d2Pl -g-

=R,-k,p,p,+k-cp,,

(6)

2

%

+u 2

-Da

s

2

-Dzy s

=R2,

(7)

dp,--bm-k-cp,, at if we neglect cross-diffusion effects. In Eqs. (6)-(8), D,, and D,, are the diffusion coefficients of species X along the x and y directions, respectively, and D, and D,, are the corresponding diffusion coefficients for species Y. R, and R2 are the reaction terms associated with the chemical reactions other than Eq. (4), in which species X and Y are involved. We next assume that the time scale of the reaction (4) is short compared to the other time scales of the problem, so that pp immediately reaches the equilibrium value

kw, pp=

(9)

k-,

This is an adiabatic approximation in which pp is “slaved” to the other variables of the problem. Furthermore, we will as-. sume that species S is in excess, so that ps is uniform and almost constant throughout the evolution. Under these assumptions, if we add Eqs. (6) and (8) and replace pp using Eq. (9), we find a2Pl

+1 JPI ~dt+u~-D,x~-D

d2Pl =R,, ly ey’

(10)

where we have defined a={1 + [ (k,p,)lk-,I}. This reduction is the analog of that carried out by Lengyel and Epstein’ to rationalize the diffusion coefficient difference induced by complex formation of a chemical species. In our system, both the flow velocities and diffusion coefficients are renormalized. Then, the basic set of equations can be written as d2Pl lx -z

-DlY

d2Pl ayz =I?,,

2

9

(5)

d2Pl -g -D,y

dP1

+u x

+u z

-Dzx $

(11)

2

-D2,, s

=R2,

(12)

where rl=u/a, ~~=RII~, hl,=D~,la, adD,,=D,,la. SO we get a system in which species X and Y flow at different velocities along the x direction. We can always change to a frame of reference that travels with the mean velocity of X. Furthermore, we can also rescale the x and y coordinates in

J. Chem. Phys., Vol. 100, No. 7, 1 April 1994

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

Ponce Dawson, Lawniczak, and Kapral: Turing and flow-induced instabilities

a way such as to eliminate the anisotropy in the diffusion of one species, e.g., Y. In this way, we are left with the system of equations --dP1 at 3P2

dt

d2Pl g -D

DI,

dP2

+v x

d2Pl =R1, ly -p

-D2V2p2=R2,

(13)

(14)

where we have defined v =u - ri, the new y coordinate y = dm y, and the diffusion coefficient D2=Db and dropped the tildes to simplify the notation. This is the system that we shall study in this paper. As we have mentioned before, the anisotropy will allow us to identify, for certain choices of the parameters, excited wave vectors along the x direction driven by the differential-flow chemical instability, if v f0, and excited wave vectors along the y direction as driven by the regular Turing instability if D,,lD, is large enough. We shall concentrate on the particular set of reactions that constitute the Selkov model,” which is described by A;X, k-1

X+2Y;3Y, k-2

YSB. k-3

(15)

The Selkov model was originally derived from a study of oscillations in glycolysis. It is of biological interest and has recently been proposed as a simplified model for the study of symmetry breaking of chemical concentrations within the cell4 The Selkov model exhibits a very rich variety of behavior. In the spatially homogeneous case, its phase diagram has been studied in Ref. 12, showing the existence of regions with a unique steady state, but several distinguishable relaxation behaviors (exponential, oscillatory, and excitatory), limit cycles, bistability involving steady states or steady states and limit cycles, etc. In the inhomogeneous case, the Selkov model can undergo a Turing bifurcation and form spatial patterns. The Selkov model involves reactions among four species. Therefore, in principle, we should work with a set of four coupled reaction-diffusion equations. However, following Ref. 13, we shall suppose that the densities of A and B, pa and pb , are fixed by external flows and treat them simply as parameters. In this way, we are left with evolution equations for the densities p, of species X and pa of species Y of the form14 --dP1 at

J2P*

D,x z

5213

Ill. LINEAR STABILITY ANALYSIS Many different parameters appear in Eqs. (16) and (17) and each of them can be tuned to obtain different types of behavior. However, if we keep the values of the reaction rates fixed, then we can define three dimensionless parameters that characterize the evolution of the system. In fact, in all cases, we shall consider k 1pa = 2.656 673 X 10w3, k,=k-,=O.OlS 000, k,=6.650 OOOX1O-3, k-1 =6.650 OOOX10M4,and kv3pb=5.313 340X10w4. These are the actual dimensionless numbers that we use in our simulations, in which time is measured in time steps At, space is measured in units of the lattice spacing Ax, and mass is measured in units of the mass of species X (see the Appendix). For these parameter values, the system possesses a unique stationary and homogeneous solution, which in di1.331 141 and pz mensionless units is p1* = = 0.346 286. Moreover, it is possible to show that it is an activator-inhibitor system, with X the inhibitor and Y the activator. Then, the parameters that govern the evolution of and system are d,=D,,lDz, the d2=DlylD2, d3=v2/(D2k3). We can identify d, and d, as the parameters that control the Turing instability along the x and the y directions, respectively, and ds as the parameter that controls the differential-flow chemical instability. As mentioned earlier, the study of Ref. 5 was carried out for one spatial dimension, assuming that the inhibitor was immobilized. The equivalent of this situation would be, in our case, to set d, = d,= 0. We shall analyze how the asymptotic state changes for other choices of the parameters. Performing a linear stability analysis around the homogeneous stationary solution {p; ,pz}, we can find the real and the imaginary parts of the eigenvalues (+ and 7, respectively, as a function of the wave vector (k,,k,) of the perturbation. The real part (T provides the linear growth rate of the instability. Those modes with o-(k)>0 will grow in time. If all modes have negative a, then the homogeneous solution is stable. For each wave vector (k, ,k,,), there are actually two eigenvalues and only one of them can have a positive real part (labeled with the plus subscript). The real and the imaginary parts of these eigenvalues are given by c%=&

[v?! Tr(B)-+(Q-(vkx)2+{[Q-(vkx)2]2 +4(B,,

-B22)2(vk,)2}“2)“2],

(18)

and

d2Pl -D ly ay’

vk=&

[-~vkx+((vk,)2-Q+{[Q-(vkx)2]2

(16) +4(B11-B22)2(vkx)2}1’2)1’2],

ah ap2 dt+v -ax -D,V’p, =Rz=k-3pb-k3p2+k2p,p;-k-2p;.

where (17)

This is the set of equations which we study in this paper.

the

matrix

B

has

(19) elements

B 11

B12=(aRlM~P2)7 = [(aRM&dl -4,k,2-Dlyk;, B21=(aR2VG%d~and B22=[(~~2)l(~~2)1-D~(k~+k~),

Tr(B) is its trace, IB 1 is its determinant, and Q is defined by Q =(TrB)2-4jBj. With our choice of the reaction rates, it is

J. Chem. Phys., Vol.subject 100, No. 1 April 1994 Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution to 7, AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

Ponce Dawson, Lawniczak, and Kapral: Turing and flow-induced instabilities

5214

t

1.5.1o-4 u

u 0.0 0.05 05

2 0.0 0.05 4

05

FIG. 1. The linear growth rate of the instabilities (T as a function of the resealed wave vector Dzk for different choices of the parameters (a) d,=O, d*=O, and d,=2.50; (b) d,=O, d,=19, and d,=2.50; (c) d,=lO, d,=19, and d3=2.50; (d) d,=18, dz=19, and d,=2.50; (e) d,=l, d?=l, and d,=3.76; (f) d,=19, d,=19, and d,=O.

possible to show that the imaginary parts v cannot be zero if u is different from zero. For this reason, we obtain stationary patterns only if ds=O. We have plotted in Fig. 1 the real part of the eigenvalues a+ as a function of the “resealed” wave vectors k&. As we can see from Eq. (4), the growth rate cr is an even function of k, and k,. Therefore, the plot in the first quadrant (k,>O; k,>O) provides all the necessary information. Figure l(a) is the analog of the case discussed in Ref. 5, but in two spatial dimensions: d, = dz = 0 and d,=2.50. We can see that the wave vectors of the unstable modes are almost aligned along the direction of the flow (the x direction). Moreover, the most unstable mode lies along the x direction. We therefore expect a pattern composed of vertical stripes. Figure l(b) corresponds to d,=O, d,=19, and d3=2.50. This is an example in which we can identify the modes along the y direction as being driven by the regular Turing instability and the modes along the x direction as driven by the differentialflow chemical instability. We have chosen the parameters in a way such that maxk a(& = 0, k,,) = maxkx a(k;, f~,,= 0). From the linear stability analysis, we do not know a priori what kind of pattern will be stable. We shall see later that the numerical simulations eventually settle into a state characterized by the presence of vertical stripes. If we increase d2, while leaving d, constant, we favor the Turing instability along the y direction. Then, maxS o(/rx = 0, $,) will increase, while maxkx a(k,, ky = 0) will remain constant. Therefore, some type of “horizontal” pattern will be favored. If we increase d, instead, we favor the differential-flow-chemical

= 20,000

t = 40,000

t

= 80,000

FIG. 2. Snapshots of the density of species Y at times t=20 000, 40 000, and 80 000 for a case with d, =d,=O and d,=2.50. This is the extension to two spatial dimensions of the case discussed in Ref. 5, a case for which the homogeneous solution is unstable only due to the DIFICI. Since the most unstable modes have wave vectors aligned with the flow [see Fig. l(a)], the pattern is characterized by vertical stripes. In this and the following figures, the abscissa is the x direction (direction of the flow) and the ordinate is the y direction (normal to the flow).

instability. In this way, some type of “vertical” pattern will be more likely. If we increase d,, we favor the Turing instability along the x direction. However, the presence of the differential flow along this direction makes the prediction of the resulting pattern less obvious. In fact, the parameters of Fig. l(c) are d,=lO, d2=19, and d,=2.50. We see that there is no unstable mode with wave vector along the x direction. This situation can be understood intuitively in terms of the “disengagement” of the species. As we mentioned before, the basic effect of the differential flow is to disengage spatially the reacting species. The same effect is produced by a difference in the diffusion coefficients of the species. Thus, by increasing the diffusion coefficient of the X species in the x direction, we counteract the effect of the flow on species Y. The species X “catches up” with the species Y, and the growth rate of the modes along the x direction decreases. If we increase di even further, modes along this direction will eventually become unstable again. This is shown in Fig. l(d), which corresponds to d,=18, d,=19, and d3=2.50. Although it might appear that u depends only on Ikl in this case, the growth rate is actually not rotationally invariant. Therefore this case is not equivalent to the usual Turing instability without anisotropy or flow. We show the growth rate for the latter in Fig. l(f), which corresponds to dl=19, dz= 19, and d,=O. Finally, we point out that it is not necessary to have unequal diffusion coefficients in order to excite the differential-flow chemical instability, as one might anticipate from the fact that in Ref. 5 the case with an immobilized species was considered (d r = d, = 0). We plot in Fig. l(e) the growth rate for a case with d,=3.76 and d t =d2= 1, which means that both species diffuse at the same rate. We see that this figure is similar to Fig. l(a). IV. NUMERICAL

SIMULATIONS

In this section, we simulate numerically the set of Eqs. (16) and (17) in order to establish which patterns arise for the parameter values analyzed above. The numerical simulations were carried out using the lattice Boltzmann code of Ref. 15 with the modifications described in Appendix A. All the runs that we discuss in this section were made on a square spatial

J. Chem. Phys., Vol. 100, No. 7; 1 April 1994 Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

Ponce Dawson, Lawniczak, and Kapral: Turing and flow-induced instabilities

t = 15,000

t = 30,000

t

t = 60,000

= 10,000

5215

t = 60,000

t

= 70,000

FIG, 3. Snapshots of the density of species Y at times I= 15 000, 30 000, and 60 000 for a case with d,=O, d,= 19, and d,=2.50. In this case, there is a competition between the DIFICI and the Turing instability (times t = 15 000 and 30 OOO),and finally the system settles into a pattern characteristic of the DIFICI.

FIG. 5. Snapshots of the density of species Y at times t = 10 000, 60 000, and 70 000 for a case with d, = 10, d,=19, and d3=2.50. In this case, the enhancement of the diffusion coefficient of the X species suppresses the DIFICI giving rise to the formation of a horizontal pattern.

grid of 256 by 256 points with periodic boundary conditions and D,=O.O8 and we show in all figures snapshots of the density distribution of species Y at different times. By means of these simulations, we obtain results that go beyond the linear stability analysis of the previous section. We show in Fig. 2 the results at times t =20 000, 40 000, and 80000 for a case with d,=d,=O and d3=2.50 [the same as in Fig. l(a)]. This corresponds to the extension to two spatial dimensions of the case discussed in Ref. 5. Therefore, the Turing instability is absent and the homogeneous solution is unstable only due to the DIFICI. As we may see in Fig. l(a), the linear stability analysis predicts that the most unstable modes have wave vectors which are aligned with the flow. The final pattern obtained via the simulation is composed of vertical stripes that travel to the right, which is in agreement with the prediction of the linear stability analysis. Figure 3 corresponds to the parameters of Fig. l(b), d,=O, d2=19, and d3=2.50; and t=15000, 30000, and 60 000. This is a case for which we can associate wave vectors along the vertical axis with the Turing instability and wave vectors along the horizontal axis with the DIFICI. Furthermore, the parameters are such that we cannot decide from the linear stability analysis which pattern will be preferred. The simulation shows that there is a transient in which the pattern is composed of “spots” aligned along the vertical direction, indicating a competition between both instabilities. Finally, the system settles into a pattern formed of traveling vertical stripes which is characteristic of the differential-flow chemical instability. If we decrease the pa-

rameter d,, the situation is less favorable for the DIFICI, according to the linear stability analysis. We show in Fig. 4 the results of a simulation for times t = 10 000, 50 000, and 90 000, and parameters d,=O, d2=19, and d3=2.26, which only differ from those of Fig. 3 in the value of d3: d,(Fig. 4)
t

= 10,000

t

= 50,000

t

= 90,000 t

FIG. 4. Snapshots of the density of species Y at times t=lO 000, 50 000, and 90 000 for a case with d,=O, d,=19, and d3=2.26. In this case, after the competition between the DIFICI and the Turing instabilities (time t =50 000), the system settles into a pattern characteristic of the Turing instability.

= 10,000

t

= 40,000

t

= 90,000

FIG. 6. Snapshots of the density of species Y at times f = 10 000, 40 000, and 90 000 for a case with d,=O, d2= 19, and d,=O. This corresponds to a

casewhereonly the Turing instability(“along the y direction”) is present and the final pattern is formed of horizontal stripes.

J. Redistribution Chem. Phys., Vol. 100, No. 7, 1license April 1994 Downloaded 21 Apr 2003 to 142.150.192.30. subject to AIP or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

Ponce Dawson, Lawnkzak, and Kapral: Turing and flow-induced instabilities

5216 t

= 10,000

t

= 50,000

t

= 110,000

FIG. 7. Snapshots of the density of species Y at times r=lO 000, 50 000, and 110 000 for a case with d,=18, d,=19, and d,=2.50. In this case, the pattern shows the presence of spots, which are aligned along the diagonals. Though the linear growth rate is almost isotropic [see Fig. l(d)], the anisotropy is evident in the pattern.

clear that the asymptotic system state has been reached in this simulation. In this case, the pattern still shows the presence of spots, which are aligned along the diagonals and travel to the right. For other initial conditions, we found spots aligned along only one diagonal, while in the case that we show in Fig. 7, there are apparently different domains in which the spots are aligned along different directions.

to be convective for these boundary conditions, then our reaction-diffusion-flow system would be a candidate for the study of noise-sustained structures that have been considered in other contexts.‘6”7 We now determine whether the system of Eqs. (16) and (17) in one space dimension (x) with D,,=O is spatially unstable when we assume that it occupies the interval [O,a) and satisfies or(O) = pf and pa(O) = pz. As before, p: and p; are the values of the homogenous stationary solution. The system of linear equations that must be solved to determine how the perturbations fit = p1 - pr and j2 = p2 - pt evolve is similar to that discussed in Sec. III; however, the boundary conditions now require that the solution be of the form m

Pl’

a12

e~+w)pJXma+(k)

dk

u+-all

1

+eo~w)reLJx/2D2a~(k)

a12

c--all

P2=

V. THE EFFECT OF THE BOUNDARY CONDITIONS It is possible that the choice of boundary conditions can affect the differential-flow-induced chemical instability and destroy the patterns. The linear stability analysis described in Sec. III is appropriate for an infinite geometry if we let the wave vector take a continuum set of values, or to periodic boundary conditions if the wave vector takes a discrete set of values such that an integer number of wavelengths fits in the periodicity length. The simulations in Sec. IV were carried out for periodic boundary conditions and were consistent with the linear stability analysis of Sec. III. We can always imagine that periodic boundary conditions correspond to a small part of a large system where some type of moving pattern has been set up. Since the wavelength observed in the simulations is intrinsic and agrees with the linear stability analysis on the infinite line, this suggests that the boundary conditions do not impose the pattern. Furthermore, the patterns are qualitatively in accord with the experimental observations.6 These results suggest that although the patterns may be affected by the boundary conditions, they are not generated by them. The analysis and simulations of Sets. III and IV describe patterns that always travel due to the DIFICI. Therefore, we expect this instability to be “convective” or “spatial.“16 A spatial instability corresponds to a case in which any initial perturbation is damped as r+a at any fixed position in the original frame of reference, but there is some traveling frame of reference in which the perturbation grows. Thus, as the perturbation grows, and eventually spreads, both edges of it travel in the same direction, leaving behind the system in its unperturbed state. In this section, we repeat the linear stability analysis of Sec. III for another geometry and another choice of boundary conditions. For simplicity, we will only consider one spatial dimension, say x, and assume that Dlx=O, but the study can be extended to the more general case. If the DIFICI turns out

o I

sin(kx),

(20)

I

omdk[e u+(k)tevx/2D2a+(k) I

(21) where a+(k) and a-(k) are determined from the initial conditions, the matrix a has elements a 11= (~3?~)l(Jp~), and al2=(JRMdp2), a21 =(~R~Y(~PA a22 = ( 8R2)/( 8p2) calculated at the homogeneous solution {PT A$)

and 2

(22) In order to determine the long time behavior of this solution, we follow the method described in Ref. 16. Thus, we first calculate t_hevalue Re[a(k)] at k such that (dc)/(dk) =O. We find k=O and Re[a,(O)]O or Re(y-)>O. Furthermore, we see that Re(&)>O and Re(y,)>O for some real k’s. Therefore, the homogeneous solution is spatially unstable, also with this new set of boundary conditions.

J. Chem. Phys., Vol. 100, No. 7, 1 April 1994 Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

Ponce Dawson, Lawniczak, and Kapral: Turing and flow-induced instabilities

This means that the system, with these new boundary conditions, is unstable and that perturbations grow while traveling, leaving behind the system in the unperturbed state. Therefore, if we perform numerical simulations on a rather small domain with this new set of boundary conditions, we will see a traveling pattern if perturbations are applied continuously from the boundary at x=0. If not, the pattern will develop and then “escape” through the right boundary. If instead we use periodic boundary conditions, the pattern will be continuously reinjected back at x=0. The boundaries automatically perform this pumping. However, the basic features of the pattern formation processes are the same in both situations. There have been several studies of convectively unstable systems and the effect of noise at the boundaries in the development of structures (see, e.g., Ref. 17). Structures that arise in this way are called noise sustained and, from the above analysis, we see that DIFICI can give rise to noisesustained patterns. Therefore, it would be interesting to simulate a system with appropriate boundary conditions using the lattice gas technique which has intrinsic fluctuations. VI. CONCLUSIONS We studied the interaction of the Turing and the differential-flow chemical instabilities in a variety of cases through both linear stability analyses and numerical simulations with periodic boundary conditions. The effects that other boundary conditions can have on the DIFICI were also investigated. The basic features of the pattern formation processes can be understood on the basis of the analyses with periodic boundary conditions. Our simulations provide insight into the interaction between Turing and differentialflow-induced chemical instabilities. Our results should be relevant for physical systems where anisotropy and flow are involved. These features are commonly found in nature; for instance, many biological systems are highly anisotropic and are maintained out of equilibrium by flows of reagents. Though most biological processes require a more complicated description than the one provided by the equations we analyzed, the interaction of these two instabilities may play a role in biological pattern formation. We have seen that this interaction gives rise to different types of transient behaviors that reflect the competition between the instabilities. In almost all cases studied, the system eventually settled into a state characterized by a pattern composed of stripes which traveled in time in the direction of the flow. Thus far, we have investigated only a limited range of parameters, and other forms of spatiotemporal pattern formation processes are expected to occur as flow rates, diffusion, and kinetic coefficients are varied. Research on Turing and Turing-like instabilities and their interaction with other processes, such as Hopf bifurcaS tions, are currently being investigated.” All these studies show a very’ rich variety of behaviors, even in controlled laboratory experiments. The present study falls in this same general class and shows that Turing’s original scenario is not the only one that is likely to operate and that these more complex bifurcation scenarios might ultimately be more rel-

5217

evant than Turing’s original idea for a variety of biological and chemical pattern formation processes. ACKNOWLEDGMENTS This work was partly supported by the U.S. Department of Energy at Los Alamos National Laboratory and by grants from the Natural Sciences and Engineering Research Council of Canada (R.K. and A.L.). Part of this work was done during visits by A.L. at the Fields Institute for Research in Mathematical Sciences, Waterloo, Canada and the Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, NM. A.L. would like to thank both of these institutions for their hospitality and support. The simulations were done on computing resources located at the Advanced Computing Laboratory of Los Alamos National Laboratory. Use of these resources is greatly, acknowledged. We would also like to thank A. Rovinsky and M. Menzinger for discussions during the course of this work. APPENDIX We show in this appendix how to modify the latticeBoltzman method in Ref. 15 in order to simulate the set of Eqs. (16) and (17). In this case, the numerical scheme consists of solving a pair of dimensionless lattice Boltzmann equations of the form f,(x+ei

At+

l)-fs(x,4t)

= - f,(x,i,t)-f

tq(x,i,t) 7s

+k(x,i,t) b

’

s=l

and 2,

(Al) with f,(x,i,t) being the one-particle distribution function of species s at time t, position x, and velocity ei , labeled by the subscript i, and where t is measured in units of the time step At, space is in units of the grid spacing Ax, and R, is related to the right-hand side of Eqs. (16) and (17) by R,=R,/(m,At), with m, the mass of one particle of species s. In Eq. (Al), - rJ,’ [fJx,i,t) - f iq(x,i,t)] is the nonreactive collision term, for which we use the single relaxation time (rr) approximation to a local equilibrium distribution function f tq(x,i, t) that depends on time and space through the macroscopic quantities, in this case, the number density ~~,=(p,Ax~)/m,. The term [i,(x,i,t)]lb, with b equal to the total number of velocities at each grid point, takes into account the reactions. As usual, the space coordinate x and the velocity ei take on a discrete set of values, which are determined by the spatial lattice considered. Most of the numerical simulations described in this paper have been done on a square lattice, with four unit velocities, connecting first neighbors (16~4) and one zero velocity (i =O). We have also worked with a hexagonal lattice, six unit velocities (lsis6), connecting first neighbors, and one zero velocity (i=O). The scheme is basically the same as that described in Ref. 15, the only difference being the choice of the equilibrium distribution functions f eq(x,i, t) for each species. In order to simulate the anisotropic diffusion of species X, we choose

J. Chem. Phys., Vol. 100, No. 7, 1license April 1994 Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to AIP or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

Ponce Dawson, Lawniczak, and Kapral: Turing and flow-induced instabilities

5218

f qq(x,4~)=pi~,(xJ),

w

where pi’s are functions of the velocity direction i that satisfy ZFiipi= 1. In the case of the square lattice (b=S), we usep1=p3 andp2=p4, where e,=-e,=i and ez=-e4=j, and in the hexagonal case (b=7), we use p1 =p4 and el= -e4=i where and P2=P3=P5=P6r e,=-es=$i+fij) and e,=-e,=$(i--flj). In order to simulate the mean flow with constant and uniform velocity u of species Y, we choose f ;q(x,i,t) flz(x,t)

$+!jeixU+(f?ixd)2-+fi2,

= i g-i;‘,

for lSiG4,

for i=O,

$f

+u

2

-Lj2V2n2=R2,

where &=2p*(rr-l/2), Dz=5(75- l/2) in

&,=2p,(q-l/2), and square and case; &=3p2(r1-l/2), and D,,=(2Pl+P2)(~,-1/2), D2=$(r2--l/2) in the hexagonal case. Clearly, these are the dimensionless versions of Eqs. (16) and (17), from which we readily obtain the diffusion coefficients as a function of the parameters of the numerical code, namely, rS and pi. In the paper, we dropped the tildes for notational convenience, but all the results of the simulations, which are dimensionless, should be interpreted in the context of this appendix. the

W) in the case of a square lattice, and f ;q(%4t)

f+~eixv+5(eixii)2-~~2,

=

n2(Jbt)

I

3-i;“,

for l~i~6,

for i=O,

(A4) in the case of a hexagonal lattice, where V is the dimensionless velocity ii = u At/Ax. In all cases, the dimensionless satisfies density n, defined by n,(X,i,t)~Cif,(X,i,r) n,(xAt)

=c

i

f,(x,i,t)=

c

I

f zq(X,i,t)

WI

and is related to the densities pS of Eqs. (16) and (17) by ns=(psAx2)lm,.

In order to show how the solutions of Eqs. (Al) provide a solution of Eqs. (16) and (17), we perform a ChapmanEnskog expansion around the equilibrium distribution functions f, s f zq + f(,” assuming a certain scaling of the different terms. In this case: we need to assume the multiscaling a/(&)-c?, d/(&)--q R,--2, and U-E, where E is a small parameter. This multiscaling implies that the method is suitable to model solutions in the limit of long wavelength and low frequency. This also limits the values of rZ in Eq. (Al) and the size of the spatial domain that can be chosen. Once these expansions are introduced, summing the lattice Boltzmann equations over all velocities and keeping terms of up to order 2, we obtain the reaction-diffusion equations dnI -dt

_ d2nI _ Dlx --g -Do

d2nI ay” =li1

646)

‘A. M. Turing, Philos. Trans. R. Sot. B 237, 37 (1952). 2V. Castets, E. Dulos, J. Boissonade, and P. De Kepper, Phys. Rev Lett. 64, 2953 (1990); Q. Ouyang and H. Swinney, Chaos 1,411 (1991); F. Baras, J. E. Pearson, and M. Malek Mansour, J. Chem. Phys. 93,5747 (1990); A. Arneodo, J. Elzegaray, J. Pearson, and T. Russo, Phys. Status Solidi D 49, 141 (1991). 3G. Nicolis and I. Prigogine, Self Organization in Nonlinear Systems (Wiley, New York, 1977); J. D. Murray, Mathematical Biology (Springer, New York, 1989); H. Meinhardt, Models of Biological Pattern Formation (Academic, New York, 1982). 4B. Hasslacher, R. Kapral, and A. Lawniczak, Chaos 3, 7 (1993). 5A. Rovinsky and M. Menzinger, Phys. Rev. L&t. 69, 1193 (1992). 6A. B. Rovinsky and M. Menzinger, Phys. Rev. L&t. 70, 778 (1993). ‘A. B. Rovinsky and M. Menzinger (unpublished). ‘I. Lengyel and I. Epstein, Science 251, 650 (1991). ‘S. R. de Groot and P. Mazur, Non-Equilibrium Thermodynamics (Dover, New York, 1984). r”G . K . Batchelor An introduction to Fluid Dynamics (Cambridge University, Cambridge: 1967). r’E. E. Selkov, Eur. J. Biochem. 4, 79 (1968). r2P. H. Richter, P. Rehmus, and J. Ross, Prog. Theor. Phys. 66, 385 (1981). 13R. Kapral, A. Lawniczak, and P. Masiar, Phys. Rev. Lett. 66, 2539 (1991); J. Chem. Phys. 96, 2762 (1992). 141n order to be able to obtain Eqs. (16) and (17) in the four-variable case, following the arguments described in this section, we need to assume as a starting point that the fluid composed of species X, Y, A, and B flows with the divergence-free velocity II in the original frame of reference. “S. Ponce Dawson, S. Chen, and G. Doolen, J. Chem. Phys. 98, 1514 (1993). 16R. Deissler, J. Stat. Phys. 40, 371 (1985). “R. Deissler and K. Kaneko, Phys. Lett. 119, 397 (1987); R. Deissler and H. Brand, ibid. 130, 293 (1988); R. Deissler, J. Stat. Phys. 54, 1459 (1989). “J. I. Perraud, K. Agladze, E. Dulos, and P. De Kepper, Phys. Status Solidi A 188, 1 (1992).

J. Chem. Phys., Vol. 100, No. 7, 1 April 1994 Downloaded 21 Apr 2003 to 142.150.192.30. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp