Hyperbolic reaction-diffusion equations and ... - Science Direct

Bm ELSEVIER

Physica D 90 (1996) 119-153

Hyperbolic reaction-diffusion equations and irreversible thermodynamics: Cubic reversible reaction model Mazen AI-Ghoul, Byung Chan Eu 1 Department of Chemistry, McGill University, 801 Sherbrooke Street West, Montreal, Quebec H3A 2K6, Canada

Received 23 February 1995; revised 26 June 1995; accepted 7 July 1995 Communicated by Y. Kuramoto

Abstract In this paper we study chemical oscillations and wave formation in a cubic chemical reaction model (a modified Selkov model) by using hyperbolic partial differential equations, i.e., hyperbolic reaction-diffusion equations, instead of the usual parabolic reaction-diffusion equations. We have carried out a comparative investigation for the hyperbolic and parabolic reaction-diffusion equations in order to help determine which system of partial differential equations is more suitable for describing chemical oscillations and wave phenomena. It is shown that the hyperbolic system approaches the parabolic system as a parameter called the reaction-diffusion number increases to a large value. In contrast to the case of parabolic differential equations, the travelling wave speed has no singular behavior if the hyperbolic system is used. The hyperbolic system predicts that there is no necessity to take different diffusion coefficients to have the Turing instability unlike the parabolic system. Numerical solutions are presented for a one-dimensional case to show that the patterns change as the reaction-diffusion number is changed from the hyperbolic regime to the parabolic regime. The entropy productions are also calculated which show that, given the same boundary and initial conditions, different patterns have the same global entropy production, but the system gets organized to local structures of some particular frequencies and wave numbers of high entropy productions at the expense of energy and matter on the part of the rest of the frequency and wave number modes. Some modes appropriate to themselves most of energy and matter available to the system as a whole and get organized to local structures. It is also found that as the trajectories become chaotic the entropy production becomes relatively smaller than the ordered structures. The second law of thermodynamics, however, does not have direct control over such local organizations except for providing thermodynamically consistent evolution equations. In the case of a bistable system, the entropy production is found to make a rather abrupt change as the system makes transition from one stable state to another.

1. Introduction I n m a c r o s c o p i c p h e n o m e n a in b i o l o g i c a l , c h e m i c a l , a n d p h y s i c a l s y s t e m s , f o r m s a n d p a t t e r n s a r e u b i q u i t o u s l y c r e a t e d a n d m a i n t a i n e d for a s p a n o f t i m e a n d t h e n t r a n s f o r m e d to a n o t h e r o w i n g to s o m e l Also at the Physics Department, McGill University, Montreal, Quebec, Canada. 0167-2789/96/$15.00 (~) 1996 Elsevier Science B.V. All rights reserved SSDI 0167-2789(95)00231-6

120

M. AI-Ghoul, B.C. Eu

/ Physica D 90 (1996) 119-153

external or internal cause. As examples for such phenomena chemical oscillations and waves have been extensively studied experimentally and theoretically in recent years [1]. Since they are macroscopic p h e n o m e n a , continuum theory descriptions are believed to be appropriate and sought after for them. In the field of chemical oscillations and waves, parabolic reaction-diffusion equations have been used to study them theoretically. Since macroscopic processes and phenomena are believed to be subject to the thermodynamic laws, efforts have been made by Prigogine and his school [2,3] to frame a thermodynamic theory of such phenomena associated with dissipative structures. Their theory has been formulated under the local equilibrium hypothesis and linear thermodynamic force-flux relations for diffusion, heat conduction, and stress with the exception of chemical reactions. In the case of chemical reactions, the constitutive relations are generally nonlinear because of the mass action law. The linear thermodynamic force-flux relations necessarily give rise to a set of parabolic hydrodynamic equations, and the reaction-diffusion equations commonly used for studying chemical waves are typical examples. As a thermodynamic theory, the theory of dissipative structures under the local equilibrium hypothesis requires a generalization since there is no physically inevitable reason that the diffusion fluxes, heat fluxes, or stresses should be constant in time, especially, when the system is in the early stage of evolution, and should obey linear constitutive relations if the system is removed far from equilibrium or subjected to steep thermodynamic gradients. It is indeed possible to formulate a more general thermodynamic theory [4] by removing the local equilibrium hypothesis and generalize the linear thermodynamic force-flux relations to nonlinear relations. In such a generalized theory the evolution equations for macroscopic variables, such as the fluid velocity, concentrations, temperature, etc., can be hyperbolic but not parabolic. 2 In fact, the hyperbolicity of the evolution equations is a more desirable feature than the parabolicity since disturbances characterizing waves in macroscopic systems cannot propagate in infinite speed which parabolic partial differential equations for the evolution equations predict them to have [6]. This feature was already noticed by Nicolis and Prigogine in Ref. [2b], but no attempt has since been made, except in a couple of papers [7,8], to remedy the weakness of parabolic differential equations as mathematical representations for chemical waves in excitable media. It is true that as far as the mathematical structure of the differential equations in the long time is concerned, this aspect is not a great problem for many situations in practice, since if the steady state is stable the fluxes decay to their steady-state values in a sufficiently long time which is often the regime of time where experiments are done, and in that limit the aforementioned distinction between the types of evolution equations disappears because the time dependence is no longer a matter of interest. In fact, both parabolic and hyperbolic differential equations give rise to the same partial differential equations in space. In other words, if the p h e n o m e n o n of interest reduces to or requires a time-independent description, then the distinction of hyperbolic and parabolic systems disappears. However, if the system is nonlinear and evolves through transient states, despite the mathematical convergence of their structures in the longtime limit one can expect some subtle effects of the hyperbolicity that would not vanish and make the longtime behavior of the system different from the behavior of the system

2Second-order partial differential equations are classified into three different types according to the eigenvalues of the characteristic equation. Hyperbolic and parabolic types are two of them. To give simple examples, the well-known diffusion equation is a parabolic differential equation whereas the classical wave equation is a hyperbolic differential equation. The parabolic partial differential equations have only a first-order time derivative and thus the disturbance in the system instantaneously propagates to infinity (i.e., the speed of disturbance propagation is infinite), whereas the hyperbolic partial differential equations have a second-order time derivative and a finite speed of disturbance propagation. For discussions on the classification of second-order partial differential equations, see [5].

M. AI-Ghoul, B.C. Eu / Physica D 90 (1996) 119-153

121

described by a set of parabolic differential equations from the beginning. The reason is that the transient effects survive in time and affect the behavior of the solutions at later times because of the nonlinearity and the well-recognized sensitivity of nonlinear systems to their history. Motivated partly by these reasons and as an application of the generalized hydrodynamics, hyperbolic reaction-diffusion equations were studied in a previous paper by Cho et al. [8] on the irreversible Brusselator and their predictions were compared with those by the corresponding parabolic reactiondiffusion equations. In this paper, we continue comparative studies of hyperbolic and parabolic reaction-diffusion equations, this time, by using a model mimicking the nonlinearity of the modified form [9] of the Selkov model [10] for glycolysis and therefore reversible in chemical kinetics. We note that the Selkov model has been studied in the literature [11] by means of a lattice Boltzmann equation devised to yield parabolic reaction-diffusion equations. However, the Boltzmann equation generally gives rise to hyperbolic partial differential equations presented below, if the Maxwell-Grad moment method is used [4]. In addition to assessing the utility of hyperbolic reaction-diffusion equations, we examine the modes of energy-matter dissipation by the system and its possible connection with, or implications for, pattern formations when a particular pattern emerges from another as the system evolves under a given set of initial and boundary conditions. As will be shown later, a rather interesting picture emerges as to how the system utilizes energy or matter to organize itself into a pattern in space-time and maintain it. We hope to better understand eventually the role of thermodynamic principles, if they play a role at all, in pattern formation through this type of study. As in the previous paper [8] of this series, we assume that the system is free from convection, the temperature is uniform so that there is no heat flow, and there is no stress applied to the system. Furthermore, the fluid is incompressible. Especially, the first set of assumptions is liable to be broken in real systems, but they are made to reduce the number of evolution equations so that they become numerically more tractable. These assumptions are also implicit in almost all the works that deal with reaction-diffusion equations in the literature. They will have to be removed in more complete studies in the future. We consider the case of linear diffusion flux evolution equations since it is desired to make comparison with the results by the conventional parabolic reaction-diffusion equations. One must bear in mind that the diffusion flux evolution equations are not generally linear if there is a steep concentration gradient or diffusion is coupled with the stress in the system. Studies of these situations are reserved for future work. This paper is organized as follows. In Section 2 a set of evolution equations for concentrations, diffusion fluxes, and other macroscopic variables are presented in a general context in order to introduce the reader to the generalized thermodynamic theory [4] underlying the present study, The set is then reduced to the canonical forms applicable to the problem considered in this work. Presented together with the evolution equations is a rather brief listing of the basic equations in the attendant theory of irreversible thermodynamics which provides the thermodynamic basis to the model mimicking the modified Selkov model. This thermodynamic theory applies to other models of chemically reacting systems relevant to chemical oscillations and waves. This section includes the general formula for energy-matter dissipation. In Section 3, the evolution equations for the model assumed are presented and cast into dimensionless forms. These reduced evolution equations are examined of their mathematical features and their homogeneous and inhomogeneous steady states are studied together with the stability analyses of them in Section 4. In Section 5, numerical results are presented for the solutions of the evolution equations. The numerical results are organized so as to show the patterns emerging as the parameters are changed. We also present the results for energy-matter dissipation as various patterns evolve under given initial and boundary conditions. Section 6 is for discussion and concluding remarks.

M. Al-Ghoul, B.C. Eu / Physica D 90 (1996) 119-153

122

2. Hyperbolic reaction-diffusion equations and irreversible thermodynamics Chemical oscillations and waves are macroscopic phenomena and, as such, their description is made by means of the evolution equations of macroscopic variables which must be subjected to the thermodynamic laws. The generalized hydrodynamic theory [4] on which the present work is based is consistent with the thermodynamic laws. However, since the theory is not well known to the research workers in the field of chemical oscillations and waves, we would like to give, as a brief summary, a collection of evolution equations for macroscopic variables required to describe flow phenomena in a general context. The reader interested in the details of the fundamentals of this theory is referred to the literature on the subject [4]. Suppose the fluid consists of an r component reacting mixture. The generalized hydrodynamic equations for the system, consistent with the thermodynamic laws, are the partial differential equations governing the fluid density p, mass fractions ci of species i, velocity u, internal energy density ~, diffusion flux Ji, net heat flux Q'i = O . i - - fliJi given in terms of the heat flux Qi and the enthalpy/~i per mass of species i, shear stress//~ = P ~ - ½Tr P~, excess normal s t r e s s A i = Pi- Pi, etc., where Pi is the stress (pressure) tensor of species i and Pi is the hydrostatic pressure: Op Ot - - V . p u ,

(2.1)

pd,c i = -V'J i + AT(c),

(2.2)

pdtu = -V. P + pF ,

(2.3)

pd,~ = -V.Q

(2.4)

- P : Vu + ~ F~ .J~ , i=l

p d j i = - V . P i - pi(d,u - Fi) - J~ "Vu + A y ) ,

(2.5)

pdtO.; = - V . ql} Q) - (d,u - Fi)" (Pi - P i ~ ) - (Q'i - q~}Q)).~14 - J i d f l i - Pi "Vhi "~ AI Q) ,

(2.6)

pd,fI~ = - V . dt (s) - 2[(dtu - F i ) J i ] (2) - 2[(//~ + dis ) .Vu] (2) 7-i

(2.7)

pdtzi i

= - - V " Ib _ ~(B)

- ~2 (d,u - F i ) ' J ~

-

2pi[Vl/] (2) + AI s'

- ~2 (II i + A f t ) . Vu - p i d t ln(piv 5/3) - V ' ( J i P i / P i ) + A} B) .

(2.8)

In these equations, the various symbols are defined by Ji:Ji/p,

Ori=Q;/p,

l~li=l~i/p,

/~i~--Ai/~,

e:~ei, i=l

Q:~Qi, i=1

A i(Q) ' --i A (s) , and --i A -0.

q>l

k-I

(2.12a)

i-1 q>l

Here

Q~ : Q , - 12iJi + I~qiXqi.

(2.12b)

The --=c is called the calortropy production. This condition (2.12a) is a generalized form of the calortropy production obtained in Refs. [12,13]. The present generalized form includes a reactive contribution absent in the previous form. If we adopt the assumptions stated earlier, then there remain the evolution equations (2.2) and (2.5).


124

Furthermore, we linearize the kinematic terms Ys; and retain only the thermodynamic force term depending on the concentration gradient:

~Ji "~- -- (kB T / m i ) V P i

(2.13)

.

We also linearize the dissipation terms to the linear form: A~') = - ~

LJ,,

(2.14)

j=l

where L u are phenomenological coefficients obeying the Onsager reciprocal relations:

Lq = L j;.

(2.15)

The diffusion coefficients are defined in terms of L = (Lq):

Dq = ( k B T / m , ) ( L - ' ) q .

(2.16)

In addition to the aforementioned assumptions, the constitutive parameters Xj, are approximated by the linear form [4] (2.17)

Xji : -Ji/Pi : -Ji/ci • In these approximations, the calortropy production (2.12a) takes the form

I.~' ,.~¢ : - T

-1 ~ [Ji "V~lbi- I~Ji "~Xji~- 2 ~ilJikRk] ~- T-1 ~ Xji(~Ji ~- A I J ) ) i=1

k=l

>0-

-

(2.18)

i=1

The t e r m E i i~ji°VXji does not appear in the approximation neglecting higher-order fluxes other than Ji. Furthermore, if the solution is assumed to be ideal, then the chemical potentials are given by = tz°i(T, P) + k B T In p;. Consequently, by using this form of chemical potential /~; and approximations (2.13) and (2.17), we find

- J i "V~Jbi31-~ Xji~Ji = O . i=1

Moreover, the chemical reaction part of the calortropy production may be written in the form --chem i=I k=l

~iVikRk

= k B ~ (R~k*) - R~k-)) ln(R~k+) /R~-)) ,

(2.19)

k=l

where R~k÷) and R~k-) are the forward and reverse reaction rates of the kth reaction, respectively, for which we have used the mass action law. On use of (2.14) for AI s) in (2.18), we find the diffusion part of ~c to be in the form

i=1

XjiA}s): i=lj=l ~ ~ p:~LisJi'J s •

Therefore, the total calortropy production is given by the formula


=

-

ln(R~

/R k )+

k-=l

p,

LJ,.Jj>-O.

125

(2.20)

i~l j-I

In the approximations enumerated earlier, this form of calortropy production -=c completely coincides with the entropy production in the theory of linear irreversible processes [2a]. Therefore, we will henceforth use the more familiar but less general term, entropy production, for it in deference to the traditional terminology. (This does not mean that the calortropy production is the same as the entropy production in general.) This form of entropy production clearly consists of two parts: one is chemical and the other is diffusion. It will be examined numerically when various patterns emerge for the solution of the hyperbolic partial differential equation system under consideration. Under the assumptions enumerated earlier, the hyperbolic differential equations are ~p~=-~.J~+

( i = 1 . . . . . r),

~kRk

(2.21)

k-I

0

=

Z LiJj

j=l

(i = 1 , . . . ,

r).

(2.22)

This set of evolution equations will be used for a cubic reversible chemical reaction model.

3. Cubic reversible chemical reaction model

In the study of control and dissipation in chemical engines, Richter et al. [9b] proposed a modification of the Selkov model [10] which consists of three steps of chemical reactions k L

A.

"S, k

i k2

S + 2P.

"3P, k.. 2

k3

P.

"B k- 3

In these reactions A and B are kept at fixed concentrations and the intermediates S and P change in time. In this work, we take this chemical reaction system as a model and construct the evolution equations under the assumption that these species also depend on position. The reaction rates for the two intermediate species S and P are given by R s =ktp A - kips

R, =

2Ps0

-

- k2psp ~ + k_2p3p ,

-2P -

3Pp +

-3PB,

(3.1)

where k i and k_i represent the reaction rate constants for the forward and reverse reactions in the ith step, respectively. To make our equations simpler we assume that the phenomenological coefficients Lsp Lp, are equal to zero. We remark that this assumption is generally not valid for real systems, but it is taken here to make the equations simple to study. This assumption is also taken in virtually all works =

126

M. Al-Ghoul, B.C. Eu

/ Physica D 90 (1996) 119-153

on chemical waves in the literature. The evolution equations then are: Otp P = - ~ . J p

+ kzpsp 2 - k_zp 3 - k3p P + k_3p B ,

(3.2a)

Otp s = - V . j s

+ kip a - k_lp s - kzp 2 + k_zp 3 ,

(3.2b) (3.3a)

- LppJ e ,

OJ r = - ( k B T / m p ) V p v

(3.3b)

O,.Is = - (kB T~ m s)Vps - Ls~J s .

Here Ot = O/Ot. Since it is convenient to work with dimensionless equations, the following reduced variables are introduced: z=kat, X=

g =r/L

,

p p / ( k 3 / k 2 ) 1'2 ,

u = Je/Lk3(k3/k2)

Y=ps/(k3/k2)

O = Js/Lk3(k3/ke)

~/2 ,

A : (kl/k3)(k2/k3)l/2pA, K=k_2/k

2,

1'2 , 1/2 ,

B : (k_3/k3)(k2/k3)l/2pB,

R=k_l/k3.

With these reduced variables, the evolution equations take the forms 07X = -V~. u + B - X + X2y

- KX 3 ,

O,Y = -Ve. v + A - RY - X2y

O~u = - U r d f ( f f ) x V ~ g

+ KX 3 ,

(3.4b) (3.5a)

+ u),

O,V = --Nrd f - ' ( l ~ y v e Y

(3.4a)

(3.5b)

+ V) ,

where the reaction-diffusion number Nrd is defined by Nrd = [ k s T / ( m p m s ) l / 2 ] k 3 1 ( D s D p )

-1/2

(3.6)

and other dimensionless parameters are: f=

(3.7)

( m s D s / m p D p ) ~/2 ,

D x =-Dp = Op/(kaL2),

Dy =---Ds = Os/(k3LZ).

(3.8)

The diffusion coefficientsD e and D s are defined by the diagonal components of D# in (2.16). They are essentially the self-diffusion coefficients. The reaction-diffusion number gives a measure of relative time scales of the diffusion flux evolution to the chemical evolution of material species X and Y when the reduced time is reckoned in the relative scale of the rate constant k 3 and the mean diffusion constant to the sound speed; see (3.6). The characteristic determinant of this system of differential equations (3.4) and (3.5) is


-A 0 NrdDxf

0 -A 0

1 0 -A

0

0

Nrd b y / f

0

-h

127

1

0

(3.9)

=0

which gives the eigenvalues

/~1' )I'2' /~3' a4 =

+~(Nrdb.f), -+ V(NrdDy/f)

•

(3.9')

Since these eigenvalues are all positive and real, the system is seen to be hyperbolic. If the mean diffusion coefficient is of the order of 10 -9 m 2 S - I and the rate constant is of the order of l 0 1 2 S I, then the reaction-diffusion number Nrd is of the order of 10 2 at the room temperature. In this case, the transient time regime must be described by the hyperbolic wave equations, but in the iongtime regime the system behaves as if it is parabolic. It is helpful to put the equations in the form of wave equations to have better insight to them in relation to the classical theory of waves. To this end, if the first two equations (3.4a) and (3.4b) in the aforementioned set of evolution equations are differentiated with r and the flux evolution equations (3.5a) and (3.5b) are used, the following second-order partial differential equations in time and space are obtained:

aZ,z + N~dHa~Z- N~dDV~Z= NrdR(X,

r),

(3.10)

where Nrd = Nr~f, Z=

{q y

,

(3.11a)

0 ][

R = [(A B - X + X 2 y - K X 3 -RY-X2y+Kxg)f-2]=[O1

D=

[Do o

H =

LH,, ~

Oyf-2

=

[lo

f 2Jk 0

,-,xq H.,J'

f

ol -=k6,

/)y

B-X+X2y-KX3 1

2 LA - R Y - X 2 y

+

KX3J ==-f2R'"

(3.11b)

(3.11c) (3.11d)

with the definitions H,~ - 1 + (1 - 2 X Y + 3KX2)/]Vrd ,

H,v := --g2//~rd , Hv~ = ( 2 X Y - 3KX 2)/]Vrd , Hv v

= f - 2 + (R + X2)//~rd .

(3.12)

Eq. (3.10) is a set of coupled nonlinear telegraphist equations. This kind of wave equations also appears in the context of shear flows in a non-Newtonian fluid [4]. Since the eigenvalues in ( 3 . 9 ' ) - r o u g h l y

128


I Physica D 90 (1996) I19-153

speaking, of the inverse coefficient matrix to the spatial second derivative term in (3.10) - are the wave velocities (group velocity, more precisely) which are proportional to 3 ~ r d , we see that as ~ d or the wave speed increases to infinity, the coupled wave equations reduce to the parabolic reaction-diffusion equations

O,Z - DV2eZ = R'(X, Y),

(3.13)

for which use is made of H. = H(2qrd = ~) =f2 and

o f-2

f2 =

,

x

~=

o l~y

'

R'

=

LA - R Y - x 2 Y + K x

3

Eq. (3.13) clearly shows that the parabolic differential equation system is the infinite wave speed limit of the corresponding hyperbolic differential equation system. Since the disturbance cannot propagate at infinite speed, the hyperbolic system is physically more appropriate than the parabolic system, which can give rise to physically unrealistic situations. For example, if the medium is inflammable, our everyday experience indicates that the flame set off at a point (say, a boundary) cannot propagate instantly through the medium to infinity. But the parabolic reaction-diffusion equations would predict otherwise. In the case of one spatial dimension, the wave equations (3.10) can be cast into another form. Introduce the new variables a + and a - , a

+

=r+c~, m

where

a

(3.14)

=r-cl~,

~

c=[NrdW(OxDy)]

-I/2.

Then with the abbreviations al

and a 2 = 0 / & e - ,

the wave

equation takes the form

(3.15)

(I - f)(021 + O~)Z + 2(1 + f-')a,a2Z + Nr~H(a, + a2)Z = ~',dR(X, Y). T h e r e f o r e , in the case of f = 1, namely, m s D s = m p D p , the equation becomes

(3.16)

40102Z + 37rdX(a , + 02)Z = ~7,dR(X , Y ) .

This suggestive form of nonlinear wave equations will be studied in detail elsewhere. Before closing this section, we further consider the difference between the hyperbolic and parabolic systems with respect to the wave speed. For this discussion we take the A-to model [14,15] that has been used for calculating the wave speed for a parabolic system. In the ),-to model the reaction source term is given by

where A and w are real functions of r = (X 2 + y2)1/2. It must be noted that the reaction source term is still cubic in X and Y. The A(r) is assumed to have an isolated zero r0: A(r0) = 0 whereas to(ro) # 0. For 2 _~_

2

~

3The mean wave velocity c may be defined by 2/c 2 1/c x 1/cy. If Dx =

2

~)yf ,

-

~

1/2

then c = (N,~Dx) .


129

this model the source term FI in (3.10) is given by

R=

of2 , f-2

y =LR'.

(3.18)

T o m a k e the equations less c u m b e r s o m e we will assume f = 1 a n d / } x = b y = 1. In this case,

and

n =

- N,-d'

-- "~'-dl LXwx + Vk x Xwy + VkyJ '

(3.20)

where Ax = O A / a X , Wy = O w / O Y , etc. Since it is convenient to seek the polar representation of X and Y, we set

[x]

rrCOSOl =

Lsin 0 J "

T h e n the wave equation (3.10) with FI and H defined by (3.19) and (3.20), respectively, m a y be written in the following two equations: rrr

--

rO2r -- ]Wrd(r~ -- rO~)

+

(/Vrd -- A)r~ +

rO, w - rArr , = / ~ r d r A

,

rO~ + 2r~O~ - A?rd(2r~0¢ + rO~e ) -- r oo + r(/Vrd -- A)0~ - r~orr ~ = N r d r W ,

(3.21b)

where the subscripts ~- and ~ denote partial derivatives with respect to z and ~, respectively, and the subscript r denotes the derivative with respect to r. If there is a limit cycle, then r = a, a being a constant. F u r t h e r m o r e , since a travelling wave is looked for, we set 0 = o ' r - k~. Then, (3.21a) and (3.21b) reduce to the form k 2 - Or2//~r, 3 -- Ovt.O(O~)//~rd = A ( a ) ,

(3.22a)

[1 - A(a)/Nrd]o- = w ( a ) .

(3.22b)

Solving these equations for o- to calculate the dispersion relation, we obtain the travelling wave speed c f r o m the dispersion relation c=°/k=(lV~d/2)'/2

(

l_A(a)/k 2

~1,2

(3.23)

1 - A(ot)12NrJ

In the case of the parabolic system for the A-~o model, the travelling wave speed c is given by [15] c = ~r/k = w ( a ) / A ( a )

(3.24)

'/2 .

Thus, if the limit cycle approaches r 0, namely, if r = a ~ r 0, then A ( a ) ~ 0 and the wave speed diverges. In contrast to this, the wave speed of the hyperbolic system, (3.23), predicts that C'-*(/Vrd/2) 1/2

asr=a--->r

o .

(3.25)


130

/ Physica D 90 (1996) 119-153

If A = 7 - r 2, then r 0 = V~ and if there is a limit cycle of radius r 0 -- VY, then there is a travelling wave of critical wave number k c such that k c = O'(2//~rj)

1/2

(3.25')

which travels at the well-defined speed c given by (3.25). This result is not possible to obtain from the wave speed formula (3.24) predicted by the parabolic differential equations (3.13) in the A-w model since it diverges at the limit cycle. For this reason, if the parabolic equations are used for the h-o) model, travelling wave solutions can be examined only in the neighborhood of the limit cycle as was originally done by Koppel and Howard [15]. The hyperbolic system is free from such a difficulty. This is another distinguishing feature of hyperbolic and parabolic systems of differential equations.

.

Steady states

The steady states can be homogeneous or inhomogeneous. Since the properties of homogeneous and inhomogeneous steady states are significantly different, they will be studied separately. 4.1. H o m o g e n e o u s steady states

If the system is homogeneous in space, then V2~Zh = 0

(4.1)

and the wave equation becomes O~Zh + ~lrdHhO,Zh = l~lrdFl(Xh, Yh) ,

(4.2)

where Zh, Xh, and Yh denote the spatially homogeneous solutions and H h stands for H evaluated with such solutions. The homogeneous steady state solution is then defined by 0~Z0 = 0

(4.3)

or

R(X0, Y0) = 0.

(4.3')

Written out explicitly, this equation yields a pair of the algebraic equations B - X o + X2oYo - K X ~ = O, A - R Y o - XZYo + KX3o = 0.

(4.4)

This set may be rearranged to the equations Yo = ( a + B - X o ) / R ,

(4.5a)

K X 3 + a2XZo + a l X o + a 0 = 0 ,

(4.5b)

where a 2 = - ( a + B)/(1 - K R ) ,

(4.6a)


131

a 1 = R/(1 - K R ) ,

(4.6b)

a o = - R B / ( 1 - KR) = - B a 1 .

(4.6c)

T h e r e can be one or three real roots of the algebraic equation, depending on the values of the parameters. There are two parameters which are experimentally variable: A and B. Therefore, it is useful to investigate the domains of the variables A and B where the root structure changes. For this purpose it is convenient to use a new parameter set (B, C), where C = A + B, instead of the set (A, B). The discriminant of (4.5b) takes the following form:

f ( B , C) = 4a -2[B2 + 4AB/a

+ 4(A 2 - q3)/o~2] ,

(4.7)

where o~=al ,

A(C) = - a 2 C / 6 R + a3C3/27R 3 ,

q(C) = (a/3 )(aC2 /3R 2 - 1).

(4.8)

T h e r e f o r e , there is one real and two complex conjugate roots if f i B , C) > 0 which yields the conditions B + 2[,~(C) --- q3/2(C)l/a > 0

(4.9a)

B + 2[A(C) -+ q3/2(C)]/ot < 0 .

(4.9a')

or

T h e r e are three real roots if f i B , C) < 0 for which there holds the condition - 2 [ A ( C ) + q3/2(C)]/ot < B < -2[A(C) - q3/Z(C)]/ot.

(4.9b)

T h e r e are three real roots with two of them being equal if f i B , C) = 0 which is satisfied by B = -2[,X(C) - q3/2(C)]/a.

(4.9c)

In Fig. 1, the curve f(B, C ) = 0 is plotted in the (B, C) plane where the domains of positive and negative f ( B , C) are indicated: one real root for f ( B , C) > 0; three real roots, of which two are double roots, for f ( B , C ) = 0, that is, on the curve; three distinct real roots for f(B, C ) < O. T o analyze the stability of the homogeneous steady states in the time domain, we define the fluctuations of the solution from the steady state (X 0, Y0) and linearize the wave equations. Thus, with the definitions

x=X-Xo,

Y=Y-Yo,

we obtain

og~Zl +/~rdH00tZl =/VrdFloZl .

(4.10)

H e r e various symbols are defined as follows: Zt=[y ] ,

(4.11a)


132

1.30C

Physica D 90 (1996) 119-153

(a) Itablo

1.04C

Hopfbcanch iIm.a0t-'b~,lo

0,785

,,nst~lo

0 0.528 A I

H ~ 0270

mta~lo

0.013

i

i

(tO

120

i

i

180

270

k

Fig. 1. Stability phase diagram for h o m o g e n e o u s steady states. In the domain where f(B, C) > 0, there are one real root and two complex conjugate roots. In the domain where f(B, C) < 0, there are three distinct real roots. O n e root is unstable and the other two are stable; thus the system is bistable. O n the curve f(B, C) = 0, there are one single and a pair of double roots all of which are real.

[ - 1 + 2XoYo - 3KX2o Ro = L(-ZXoYo + 3 K X Z ) f -2

[-°x

X2o ] -(R - xZ)f-zJ '

1

H 0 = LHOrx ny°yJ '

(4.11b)

(4.11c)

with the definitions H x O x = l + ( 1 - 2XoY° + 3 K X o2) / N-r d , n x 0y = __X2/l~rd

z Hyox = (2XoYo _ 3KXo)/N~d o

Hyy = f

+

(R +

(4.11d)

Eq. (4.10) will be analyzed together with the inhomogeneous steady state.

4.2. Inhomogeneous steady states If the system is inhomogeneous in space, the steady state is defined by the condition 07Z t = 0

(4.12)

133


which yields the equation

-DV~Z,= R(At, Y,),

(4.13)

where Z,=

[Art] i..,

(4.14a)

,

R, = R ( X = A t , Y = Y , ) .

(4.14b)

At this steady state, the hyperbolic and parabolic systems of differential equations coincide. However, the temporal evolutions of the two systems (3.10) and (3.13) differ significantly since the manner in which the inhomogeneous steady state is asymptotically reached by the two systems can be qualitatively different. The evidence for this difference will be shown numerically in a later section. Eq. (4.13) describes steady-state Turing structures for the system of interest. Since it is useful to examine the linear stability of such structures, the evolution equation (4.13) is linearized with respect to the state of vanishing reaction rates, namely, (X0, Y0) defined earlier. Then, with the definition

x=At-X

o,

y=Yt-Yo,

(4.15)

Eq. (4.13) may be written in terms of new variables x and y as

-~7~Z, = MZ, + A(x, y) ,

(4.16)

where M = - D - 1 R0 ,

A(x, Y)

D

1[-

I

(4.17a)

Kx3 + (3KX0 - Y0)x2 - x2y - 2Xoxy ] _ [Kx3 + ( 3 K X ° _ Yo)X 2 + xZy + 2XoxY]f_2j,.

(4.17b)

Dropping the A term yields the linearized evolution equation V2~Zt= MZ t .

(4.18)

Combining this equation with the temporal equation (4.10) we obtain the linearized wave equation c3~Z! +/~rdHo0rZ/- NrdD~72~Z/-- -/VrdDMZt

(4.19)

which may be used to calculate the linearized waves and their dispersion relations. Solving this equation by the Fourier transform Z, = ~'~ ~ ~(k, ,0) exp[i(k • ~ - `or)] oJ k

(4.20)

we obtain the linear algebraic set ( - , 0 2 / - i`0NrdHo + k2NrdD + NrdDM)~(k, `0) = 0. The dispersion relation is given by the secular determinant

(4.21)

134

M. A l - G h o u l , B . C . Eu

/ Physica D 90 (1996) 1 1 9 - 1 5 3

d e t l ( - i w ) 2 1 - io~.AIraH0 + kZ/~rd D -I- NrdDM] = 0 .

(4.22)

This is a fourth o r d e r polynomial of z = - i w and k: P4(z,

k)

=

Z 4 + P z 3 + Q z 2 + Tz + S = 0

(4.22')

with the definitions Kxx = kZDxx + ( D n ) x x ,

(4.23)

gyy = k2Dyy + ( D M ) y y ,

p = N-r d ( n x 0x + n ~ y )

(4.24a)

O = l~rd(gxx + g y y ) + "lQ2 ¢Mo LIO o o " rd~.'" x x ' ' y y -- H x y H y x )

(4.24b)

,

T = Nr~[KxxHyy -2 o + KyyUxx o - H xoy ( D M ) y x _ H~x(DM)xy] ,

(4.24c)

S~

(4.24d)

--2 Nrd[KxxKyy - (OM)xy(DM)yx]

T o write these terms out in less formal forms, we define lx=l~ro,

(4.25a)

ly=fl~d/f z ,

e x = IVrdbx,

(4.26)

ey = N, d b r / f z ,

and with the abbreviations a I=I-2XOY

o + 3 K X 2,

a 2=-X

2,

b I=-3KX~+2XoY

o,

b 2 =

R + X2,

(4.27)

the following quantities: al = e x + % , ot 2 = ex(ly + b : ) + ey(l x + a , ) , Ol3 = b 2 E x l y -I- a l E y l x ,

/31 = (b 2 + a l ) ( l x + ly) + lxly + a l b 2 - a2bl , /3z = (alb2 - a2b,)(lx + ly) -~- lxly(b 2 ~- a , ) , /33 = (alb2 - a 2 b l ) l x l y , e = ex% .

(4.28)

T h e n , P, Q, etc. can be written as P=b2+al+lx+ly,

(4.24a')

Q = 41 k2 +/31 ,

(4.24b')


(4.24c')

T = ~2 k2 + f12, S = 8 k 4 + a3 k2

135

(4.24d')

+ f13"

For the polynomials to have all the roots with the negative imaginary part the following Hurwitz conditions [16] must be fulfilled: P > 0, PQPQT-

(4.29a) T>0,

(4.29b) (4.29c)

T 2 - P2S > 0 .

The first condition (4.29a) is independent of k and a function of B and C. The second condition T 2 - PZS are even functions of k. Since the polynomial P4(z, k) is real, it has two pairs of complex conjugate roots. The level curves of the Hurwitz conditions f2(k, C, B) = 0 and f3(k, C, B) = 0 for a given value of B are plotted in Fig. 2. Note that the first condition gives a point on the C axis. On the level curves the real parts of the roots vanish and the solution of the wave equation (4.19) becomes oscillatory in time. In the case of the parabolic system the linearized evolution equation takes the form

f2(k, C, B ) ~ P Q - T and the third condition ~ ( k , C, B ) =- P Q T -

Or Z, - DV ~ Z, = - [~MZ,

(4.30)

which gives rise to the secular equation z 2 + P ' z + Q' = 0,

(4.31)

where p, = a 1 + b 2 + (19 x + l~)y)k 2 = a 1 + b 2 + (lxly)-a(exlr + eylx)k 2 , Q,

= DxDy k 4 + (al/~y + ~ ~

b 2 D~ x)k

2+

alb2

_

bla2 = f 2 S

.

(4.32a) (4.32b)

The Hurwitz conditions for negative real roots are fi' ~ lxlyP' = (a~ + b2)lxly + (exly + sylx)k 2 > O, p,Q, ~p,sf2>O.

(4.33)

The level curves are plotted in Fig. 3 which are qualitatively different from the level curves for the hyperbolic system present in Fig. 2. There is only one level curve since the second condition Q ' ( k , C, B ) or S(k, B, C) is always positive for the parameter set chosen for C and B and hence does not appear in the figures. Comparison of the two figures indicates that the domain of instability is much larger for the hyperbolic system than the parabolic system which has a maximum k value beyond which instability does not occur whereas there is no such maximum value for the wave number in the case of the hyperbolic system. This difference in stability diagrams indicates that the dynamics of the two systems can be considerably different. Othmer [7] also observed that the stability diagram of a

136


/ Physica D 90 (1996) 119-153

0.06

~

0.05 0.04

m o.os 0.02

f(B,C)>O

.

0.01

.

f

f

(

B

,

C

)

=

O

f(B,C) 0. Therefore, if !~ = D I , then the stable steady states for which Re Ai < 0 cannot be destabilized by diffusion. This conclusion cannot be drawn if I~ ~ DI, namely, if /)x ~ / ) y as is well known [15,17,18]. In the case of the hyperbolic system under the assumption that /)x = by and Re Ai < 0, the secular determinant for (4.19) is det[-on

- Z V ~ r ~ H ° - ( O k z + z2)l[ = 0 ,

(4.36)

where z =--ito/VN~d. The presence of the term z'V/-~d H° prevents us from drawing a conclusion regarding the sign of Re z i similarly to the one drawn earlier for the parabolic system. To make this

138


point more definite, let us assume that V ~ d H ° = h I where h is positive or negative. Then, since Re Ai < 0 by assumption and k 2 -> 0, (4,36) implies that Re(z~ - zih ) < O.

(4.38)

Setting zi = ai + i/3~, we can express this condition in the form ½[h - (h e + 4/3~) 1/2] < ai ~. We have ascertained that the hyperbolic system becomes practically parabolic and gives the same results as the latter if Nra is increased to at least about 20 and beyond. (The maximum value studied was 50.) In the regime of large Nrd, although

!

a¸

C

b

J 10 17 26 33V~

\ 20

41

1'4 20

41~

Fig. 5. (a) T h e p o w e r s p e c t r u m o f t h e h y p e r b o l i c s y s t e m at N,d = 0.1. It i n d i c a t e s t h a t t h e r e a r e t h r e e f u n d a m e n t a l f r e q u e n c i e s in t h e c a s e o f b at ~ = 1, v2 = 17, a n d v3 = 43, a n d s u b h a r m o n i c s at v~ - 2 v 2 + v3 = 10, ~ - v2 = 26, 2 v 2 - v~ = 33, e t c . (b) T h e p o w e r s p e c t r u m f o r t h e c a s e o f t h e p a r a b o l i c s y s t e m ( t h e limit o f Nrd ~ ~ ) . It s h o w s t h a t t h e r e a r e t w o f u n d a m e n t a l f r e q u e n c i e s at v~ = 1 a n d v2 = 20. T h e o t h e r is a s u b h a r m o n i c : v1 + 2 v 2 = 41. (c) T h e p o w e r s p e c t r u m o f t h e h y p e r b o l i c s y s t e m at N,d = 0 . 0 1 . T h e f u n d a m e n t a l s a r e still d i s c e r n a b l e at Vl = 1, v2 = 20, a n d ~ = 41. H o w e v e r , t h e s p e c t r u m h a s b e c o m e b r o a d e r a n d d i f f u s e , indicating a chaotic motion.

140


/ Physica D 90 (1996) 119-153

2

1 >

0 c

-1 -2

-4

-2

0

2

4

10

0

20

30

40

CO

U 0.1 0.05

>

> 0.05

o

-0.05 0 0

0.05

0.1

-0.01-0.005

0

U

0.005 0.01

U

Fig. 6. (a) Trajectory of space-integrated (mean) fluxes u and v in the uv plane. This figure is in units of 10-19 for the abscissa and 10 -18 for the ordinate. Notice that their magnitudes are rather small. Nfd = 0.1. (b) The corresponding power spectrum of u indicating a chaotic motion. (c) Trajectory of space-integrated (mean) fluxes u and v in the uv plane in the case of N,o = 0.01. The fluxes do not oscillate and their motions are regular, whereas the concentrations X and Y fluctuate chaotically in this case.

2.2

1,8

1.4

1

stable state

~Ooo 0000

0

0

0

o O

0.6

O O O O

0.2 0,1

&

0,15

X Fig. 7. Trajectory of X and Y in the X Y plane in the case of the parameters in the stable steady state domain in Fig. 2a. The trajectory tends to the stable steady state. A = 0.2, B = 0.9, N,o = 0.1, /)x = 0.006, /)y = 0.0016, and f = 1.0.


141

0.2

0.18

x 0.16

o.~4[(A) o

0.5

a

0.5

b

(B)

0.5

0.4

0.4

0.3

0.3

x

x 0.2

0.2 !

0.1

0.1 ~

1

oi

0,5

0,5

C

d

0.8

0.5

0.7 0.4

0.6 0.5

0.3

X 0.4

X 0,2

0.3 0,2

0.1 0,1 0q

0

0.5

0 0

0.5

Fig. 8. (A) Comparison of travelling wave patterns predicted by the hyperbolic and parabolic systems. The solid curve is for the hyperbolic system whereas the curve of open circles is for the parabolic system. N~d = 0.1, Notice the sharp fronts in the case of the hyperbolic system. (B) Sequence of wave merging and splitting: a---~b---~c---~d. This process is repeated over a long time span.

142

9

10


/ Physica D 90 (1996) 119-153


143

the maximum amplitudes are a little reduced, the trajectories are similarly quasiperiodic in the X Y plane (Fig. 4c) as is the hyperbolic system, but the power spectrum in Fig. 5b shows that there are only two fundamental frequencies. Since this seems to suggest that the reaction-diffusion number plays a role in changing the number of fundamental frequencies, the value of Nrd is further decreased to Nrd = 0.01 SO that the hyperbolicity is amplified. It must be remarked here that for this case the initial conditions for fluxes were taken u(~, 0) = 0.1 and v(~, 0) = 0 to get physical solutions to the hyperbolic differential equations. Below this particular value of Nrd = 0.01 the present set of boundary and initial conditions failed to produce physical solutions. In any case, at Nrd = 0.01 the motion now begins to show a chaotic behavior and its power spectrum (Fig. 5c) shows a broad and diffuse structure. This trends probably will continue as the reaction-diffusion number gets smaller and the initial and boundary conditions are varied suitably. From this investigation the following picture seems to emerge for the particular case of parameters of the system: If there is no diffusion, there is a single frequency limit cycle. As diffusion is turned on at a very large value of Nra so that the disturbance propagates at infinite speed and thus the system is parabolic, there appears a two-fundamental frequency quasiperiodic m o t i o n - period doubling. As the reaction-diffusion number is reduced so that the disturbance propagates at a finite speed and the system thus becomes hyperbolic, the number of frequencies is increased to three or more and eventually to infinity, and the motion becomes chaotic. Unlike the parabolic system the hyperbolic system enables us to examine the evolution of diffusion fluxes in time. Intuitively, we know that the parabolic system of hydrodynamic equations emerges when the diffusion fluxes change rapidly compared with the concentrations. This picture is borne out when the reaction-diffusion number is sufficiently large. When the value of the reaction-diffusion number is in the range where the trajectories are not chaotic in the X Y plane (e.g., Nrd -->0.1), the diffusion fluxes fluctuate rapidly and their spatial mean values, namely, their spatial integrals, become very small in magnitude, indicating that they do not change over a significant span of time. This is shown in Fig. 6a and the power spectrum is shown in Fig. 6b. As evident from the figures in this regime of the reaction-diffusion number the diffusion fluxes may be regarded as stochastic variables. However, as the reaction-diffusion number is lowered well below Nrd = 0.1, say, to Nrd = 0.01 and the trajectories of X and Y become chaotic, the trajectories of u and v become regular as is shown in Fig. 6c which shows the trajectory of (u, v) integrated over space, namely, the spatial means of u and v. We thus see that the stochasticity is exchanged between the two sets of variables (X, Y) and (u, v), that is, the slow variables have become fast variables and the previously fast variables have become slow variables. This indicates that the hyperbolicity will become significant as the matter begins to diffuse on the time scale of change comparable to the concentration variations in the reacting species. This aspect of course is absent in the parabolic system since the fluxes are necessarily steady, namely, u = - D x V X and v = - / ) VY. Fig. 6d shows the trajectories in the (u, v) plane which are not integrated over space unlike those in Figs. 6a-c. If the parameters are chosen in the stable steady-state region of Fig. 2, then the trajectory in the X Y plane tends to a fixed point as shown in Fig. 7. There is a characteristic difference between the behaviors of the travelling waves predicted by the

Fig, 9. S p a c e - t i m e plot of concentration waves for X in the case of the paraboli c system for X = Y = 0.2 at the boundaries. (a) X, (b) Y, (c) u, and (d) v. T h e s e result coincide with those of the hyperbolic system at Nrd = 50. In the accompanying color code the intensity increases from left to right in the order of dark orange, yellow, green, light blue, dark blue, pink, and red. All the figures in color in this work conform to this color code. Fig, 10. S p a c e - t i m e plot of concentration waves for X in the case of the hyperbolic system for X = Y = 0.2 at the boundaries and NT,~ = 0 . 1 . (a) X, (b) Y, (c) u, and (d) v.


144

11

12

/ Physica D 90 (1996) 119-153


145

hyperbolic and parabolic systems. In the case of the former, the wave maintains a sharp front whereas the wave is rather diffuse and nonzero in amplitude throughout the space in the case of the latter, reflecting an instantaneous propagation of matter. This probably is the most noticeable from the experimental viewpoint when a travelling wave is involved. This characteristic difference is demonstrated in Fig. 8a where waves travelling from the boundaries to the center are shown. The solid curve denotes the prediction by the hyperbolic system at Nrd = 0.1 and the curve of open circles the prediction by the parabolic system. This kind of difference will be later exhibited in space-time for other values of the reaction diffusion number and for different boundary conditions in this work. The mid portions of these curves are not equal to zero because the concentrations are built up locally owing to the initial conditions chosen for the numerical study made here; see (5.1). It is also noticed that after the waves travel toward the middle and merge, they split and travel back toward the boundaries, and then get reflected and travel back toward the middle. This process is repeated. An example of this behavior is presented in Fig. 8b where a sequence of wave merging and splitting is shown. This behavior appears for all values of Nrd studied: from Nrd = 0.01 on to Nrd = 20 which is the maximum value of the reaction-diffusion number studied in this work under the initial and boundary conditions. It probably happens for higher values of the reaction-diffusion number since Petrov et al. [20] have observed the phenomenon with a system of one-dimensional parabolic reaction-diffusion equations for a cubic reaction model. The wave splitting behaviors are shown in space-time in Figs. 9a and b for the case of Nrd = 50. At this value of Nrd the hyperbolic system practically behaves as if it is parabolic and gives the results virtually coinciding with the latter in all qualitative aspects. In the figures in color, the amplitude of the wave increases in the following order of colors: dark orange, yellow, green, light blue, dark blue, pink, and red as shown in the color code strip in Fig. 9. In these figures the waves are plotted in space-time over 1000 units of 7 (reduced time); the abscissa is for time (7) whereas the ordinate is for space (~). Notice that there are more than one wave length discernable and the waves merge and split in the middle of ~ = [0, 1]. The hyperbolic system at Nrd = 0.1 predicts more complicated patterns of waves in space-time which are generally shorter in wavelength. These waves and the corresponding fluxes are shown in Figs. 10a-d. The wave merging and splitting phenomenon is still discernable at this value of Nrd. These patterns may be regarded as an example of one-dimensional patterns where the spots appear and disappear as time progresses. This behavior reminds us of the recent experimental results by Lee et al. [21] on birth and demise of spots in a reacting system. In Figs. l l a and l l b the waves are plotted in space-time in the case of Nrd = 0.01. The wave merging and splitting phenomenon is still noticeable in them. To study the effect of the boundary conditions, the hyperbolic system was solved for different boundary conditions X = 1.5 and Y = 2.5 which are quite different from those for the results in the previous figures. As shown in Figs. 12, the patterns have changed significantly. We observe that the nonstationary waves in the early stages (on the left edge in the figures) evolve to a steady peak in the center flanked by steady waves of long fixed wave lengths. These may be regarded as a pattern in space-time. They appear to satisfy the necessary requirements of Turing patterns that the pattern is stationary, symmetry is spontaneously broken, and the wavelength is intrinsic [22]. In particular, the central peak may be regarded as what is termed as the Turing hole in the literature [23,24]. For the

Fig. 11. S p a c e - t i m e plot of concentration waves for X in the case of the hyperbolic system for X = Y = 0.2 at the boundaries and Nr,L- 0 . 0 1 . (a) X, (b) Y, (c) u, and (d) v. Fig. 12. S p a c e - t i m e plot of the hyperbolic system at X = 1.5 and Y = 2.5 and N~a = 0.02. (a) X, (b) Y, (c) u, and (d) v.

146

13


/ Physica D 90 (1996) 119-153


147

1.57

1.56 T 1.55 M

1.54 1.53

1"520 0.5

J 1

L 1.5

2

2.5

3

315

4

4.5

N,d

Fig. 14. Space-integrated global entropy production vs. reaction-diffusion number. Since different patterns appear as N,, is increased, this figure indicates how the global entropy production changes with patterns.

same initial and boundary conditions chosen, the parabolic reaction-diffusion equations did not produce this kind of patterns where the Turing hole is flanked by waves; they produced only diffusive waves in space. In order to see the effect of the difference in the diffusion coefficients we have set them equal and solved the hyperbolic equations. The results are shown in Fig. 13A. As is shown, there still appear space-time patterns although they are quite different and less striking than those in Figs. 12a-d which are for different diffusion coefficients. The Turing theory [17] predicts that the diffusion coefficients must be different in order for patterns to appear, but the hyperbolic system indicates otherwise as we have discussed earlier in Section 4. Numerical evidence for this can be seen in Fig. 13A for which b , = b v = 0.006 and f = 0.1 were taken. We show how the patterns change as the parameter f is changed in Fig. 13B were b x = Dy = 0.006 and f = 1 are taken. In this case, a more structured pattern appears in the early time, but it is transformed to another stationary form of a different symmetry for the rest of time. These patterns are rather sensitive to the parameters taken.

5.2. Energy and matter dissipation in the case of a single steady state Macroscopic systems evolve to a dissipative structure at the expense of energy and matter on the part of the surroundings. The mode of energy and matter dissipation to create a structure therefore is of considerable interest. The theory of irreversible processes on which the present theory is based shows that the calortropy production must be positive semidefinite for the irreversible process involved. In the case of the present chemical system described by the hyperbolic reaction-diffusion equations under the assumptions stated earlier, this calortropy production coincides with the entropy production in the theory of linear irreversible processes as we have shown in Section 2. We now investigate how the

Fig. 13. (A) Space-time plot of concentration waves for X in the case of the hyperbolic system at X = 1.5 and Y = 2.5 and N , d = ( l . l , f = O . l , and b = D =0.006. (a) X, (b) Y, (c) u, and (d) v. (B) Same as in (A) except f o r f = l .


148

15

18

/ Physica D 90 (1996) 119-153

149


0.3 0 . 2 5 ~

T

1

13 0.2

0 . 1 5 ~ 0.10 200

400

600

800

1000

co Fig. 16. (a) Space-integrated entropy production corresponding to Fig. 12. (b) Its Fourier transform plotted in the (to, k) plane. The darkest shade corresponds to the lowest entropy production. e n t r o p y production (2.20) depends on the frequency and wave n u m b e r of a pattern evolved from a given set of initial and boundary conditions. Since we have noticed that the patterns change as the reaction-diffusion n u m b e r changes, we study the Nrd dependence of the entropy production. First, we have examined the global entropy production integrated over a period, namely, the integral of the e n t r o p y production in s p a c e - t i m e for a given set of initial and boundary conditions. These values turn out to be almost independent of Nrd over the entire range except for the region of small Nrd where it decreases significantly. This is the region where the motion gets chaotic and there is no recognizable pattern. It & interesting that a chaotic m o t i o n has a l o w e r entropy p r o d u c t i o n than a structured m o t i o n . This result is presented in Fig. 14. It suggests that the global e n e r g y - m a t t e r dissipation by the system is higher when there is an organized pattern than a chaotic state, but almost independent of the organized patterns formed in the system for a given initial and boundary conditions. That is, different patterns consume about the same amount of energy and matter globally. It suggests that the global e n e r g y m a t t e r dissipation by the system is independent of what is happening locally in the system for a given initial and boundary conditions. This is reasonable and what it should be in retrospect, since the system evolves to a local structure from a h o m o g e n e o u s state with the energy and matter provided by the surroundings, but the global consumption of energy and matter should be the same regardless of which structure is formed locally, if the boundary and initial conditions are the same. Next, we take the two-dimensional Fourier transform of the entropy production in s p a c e - t i m e and plot its logarithm in the ~o-k plane in the same color code as for the previous figures in color. The abscissa is for k and the ordinate is for o~. Figs. 1 5 a - d show how intense the entropy productions are for different modes (w, k); the low o~ modes have the highest intensities for all k and are thus favored by the system at the expense of energy and matter. Notice also that there are certain k values where patterns of high frequency ~o m o d e s are formed. Fig. 16 shows the space-integrated entropy production vs. ~- for the system and its

Fig. 15. Two-dimensional dissipation (entropy production) spectrum for the hyperbolic system at (a) Nr~ = 0.1, (b) Nr, = 2, (c) Nrd = 4, and (d) N,d = 50. Boundary conditions: X = 0.2 and Y = 0.2. /5x # / ) . Fig. 18. Space-time plot of the X concentration of the hyperbolic system in the bistable regime. Notice the transition region between the two stable states and the wavy behavior induction near the boundaries. A similar behavior was obtained in the case of parabolic system. Nrd =0.1 and the boundary conditions are: X=0.1, Y=0.1. (a) Y, (b) Y, (c) u, and (d) v.

150


/ Physica D 90 (1996) 119-153

Fourier transform for patterns in Fig. 12. The Fourier transform indicates there are many patches of high amplitude modes as in the case of Fig. 15a for a different set of boundary conditions. Since the high intensity of entropy production also means a high wave amplitude, these figures show that the system tends to create low wave-number or long wavelength waves (structures) at all frequencies at the expense of high wave-number waves (structures). There also appear islands of (to, k) modes that maintain structures appropriating energy and matter from the neighboring modes. The entropy production patterns for the hyperbolic and parabolic systems virtually coincide with each other as the reaction-diffusion number is increased, indicating that the hyperbolic system tends to the parabolic system as Nrd increases. Even though the entropy production computed and the reaction-diffusion equations strictly conform to the requirement of the second law, namely, is positive semidefinite, at least for this case studied, the second law of thermodynamics does not directly control creation and selection of a particular dissipative structure; it only controls the global aspect of the energy-matter consumption of the system. Pattern selection seems to be in the province of local dynamics dictated by the dynamical evolution equations- the reaction-diffusion equations in the present case.

(a)

stable 25

>-15

"'... stable

--> o~.:.".. t 2

i 1

X

(b)

>-

x

1

0 0

' 10

20

30

40

50

"E

60

70

80

90

100

Fig. ]7. (a) Phase space trajectory in the case of a bistable system. The trajectory moves from the lower corner to the upper left corner if the initial state was in the vicinity of the stable steady state at the lower corner. (b) Local behavior in the bistable region. The figure is the section at ~ = 0.0078.

151


a

lo

o

~

~

~

~

~ooo

2~

o3 Fig. 19. (a) Space-integrated entropy production in the case of the bistability. Notice the abrupt change between the two states. (b) Fourier transform of the entropy production plotted in the (w, k) plane. The darkest shade corresponds to the low intensity modes and thus low entropy production modes. 5.3. Bistable region We have also examined the region of bistability. The set of parameters is chosen to be C = 2.5, B = 0.005, and Nrd = 0.1 SO that the solution of the steady state equation gives rise to three steady states: X 0 = 2.2321, 0.00585, 0.03478. The first two are stable antt the third one is unstable. If the initial conditions are taken to be X(~, 0 ) = Y(~:, 0 ) = u(~, 0 ) = v(sc, 0 ) - - 0 . 1 and the boundary conditions are taken to be X ( 0 , t ) = X ( 1 , t ) = 0 . 1 and Y(0, t ) = Y(1, t ) = 0 . 1 at the boundaries, we notice that the system tends to the stable steady state at X 0 = 0.0585 to stay in its neighborhood for a short time and then moves to the vicinity of the unstable state at X 0 = 0.03478 which is very close to the stable state at X 0 = 0.00585, and then it gets kicked out to the neighborhood of the steady state at X0 = 2.2231. This is shown in panel (a) of Fig. 17 and the local behaviors of X and Y are presented in panel (b). Steady temporal oscillations take place near the boundaries as is shown in Figs. 18a-d. The integrated entropy production is computed for such process. It changes abruptly, as shown in Fig. 19, when the system makes transition from the unstable state to the stable state where oscillations take place near the boundaries. Such a Turing pattern is maintained in a higher dissipation state. Apparently, it costs a high uptake of energy and matter to maintain such a pattern. For the parabolic system no transition is

Fig. 20. Entropy production for a bistable case. Nrd = 0.1. This figure corresponds to Fig. 18.

152


/ Physica D 90 (1996) 119-153

induced and the system relaxes directly to the neighborhood of the stable steady state at X 0 = 0.0058547 instead of the other steady state X 0 = 2.2321. In Fig. 20 a three-dimensional plot of the entropy production is presented in space-time in the case of Nrd = 0.1. Near the boundaries where Figs. 18a-d indicate oscillatory structures the entropy production exhibits a corresponding oscillatory structure. The entropy production bursts to high peaks near the boundaries relatively early on in time and then settles down to a fairly regular structure as time progresses.

6. Discussion and concluding remarks In this paper we have presented some results of our ongoing study of chemical oscillations and waves by using hyperbolic partial differential equations- wave equations- instead of diffusion equations commonly used in the literature. The hyperbolic reaction-diffusion equations used follow directly from the irreversible thermodynamic theory of processes in systems removed far from equilibrium, and the present study is also an effort to apply the irreversible thermodynamic theory to chemical oscillations and wave phenomena. We have pointed out some formal differences between hyperbolic and parabolic systems of partial differential equations for systems of interest here and observed some numerical differences between them. Although hyperbolic differential equations are numerically more difficult to handle than the parabolic counterparts, there are conceptual and theoretical advantages to use hyperbolic differential equations as have been indicated in the text. In this work, we have also examined the entropy productions associated with wave phenomena in the hope of understanding the role of thermodynamic principles in pattern formations and chemical oscillations. We have observed some interesting features, which seem to have been obvious in retrospect. We do not have as yet a satisfactory understanding of thermodynamic reasons for pattern selections, if there is any ground for believing that there should be a thermodynamic reason at all. This line of question is completely open and should be a subject of more intense study in the future. The present thermodynamic theory is hoped to provide some theoretical framework for such study, and the present numerical solution approach gives reasons to believe that hyperbolic differential equations are reasonable evolution equations for chemical oscillations and waves. A study of a two-dimensional system is in progress at present. Many patterns including spirals are observed in this case. The results will be reported in the near future.

Acknowledgments This work has been supported in part by the Natural Sciences and Engineering Research Council of Canada and Fonds pour la Formation de Chercheurs et l'Aide h la Recherche of Quebec. We would like to express our thanks to Prof. David N. Harpp for his invaluable assistance in the production of the color figures.

References [1] Y. Kuramoto, ChemicalOscillations,Waves, and Turbulence (Springer, Berlin, 1984); F. Baras and D. Walgraef, eds., ChemicalDynamics:Experiments to MicroscopicSimulations, PhysicaA 188, Nos. 1-3 (1992).


153

[2] (a) I. Prigogine, Thermodynamics of Irreversible Processes (Interscience, New York, 1961); (b) G. Nicolis and I. Prigogine, Self-Organization in Nonequilibrium Systems (Wiley, New York, 1977). [3] P. Glansdorf and I. Prigogine, Thermodynamic Theory of Structure, Stability, and Fluctuations (Wiley, New York, 1971). [4] See, e.g., B.C. Eu, Kinetic Theory and Irreversible Thermodynamics (Wiley, New York, 1992). [5] S.G. Mikhlin, Mathematical Physics, An Advanced Course (North-Holland, Amsterdam, 1970). [6] I. Miiller and T. Ruggeri, Extended Thermodynamics (Springer, Berlin, 1993). [7] H.G. Othmer, J. Chem. Phys. 64 (1976) 460. [8] U.I. Cho and B.C. Eu, Physica D 68 (1993) 351. ]9] (a) Y. Termonia and J. Ross, Proc. Natl. Acad. Sci. USA 78 (1982) 2952, 3563. (b) P. Richter, P. Rehmus and J. Ross, Prog. Theor. Phys. 66 (1981) 385. [10] E.E. Selkov, Eur. J. Biochem. 4 (1968) 79. [11] S. Ponce Dawson, A. Lawniczak and R. Kapral, J. Chem. Phys. 100 (1994) 5211. [12] B.C. Eu, Phys. Rev. E 51 (1995) 768. [13] B.C. Eu, J. Chem. Phys. 102 (1995) 7169. [14] N. Koppel and L.N. Howard, Stud. Appl. Math. 42 (1973) 291; L.N. Howard and N. Koppel, Stud. Appl. Math. 56 (1977) 95. [15] J.D. Murray, Mathematical Biology (Springer, Berlin, 1989). [16] F.R. Gantmacher, Theorie des Matrices (Dunod, Paris, 1966) Vol. 2. [17] A. Turing, Phil. Trans. R. Soc. London Ser. B 327 (1952) 37. [18] J.E. Pearson and W. Horsthemke, J. Chem. Phys. 90 (1989) 1588. [19] C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zang, Spectral Methods in Fluid Dynamics (Springer, Berlin, 1988). [20] V. Petrov, S.K. Scott and K. Showalter, Phil. Trans. R. Soc. London A 347 (1994) 631. [21] K. Lee, W.D. McCormick, J.E. Pearson and H.L. Swinney, Nature 369 (1994) 215; Phys. Rev. E., to be published. [22] P. De Kepper, V. Castets, E. Dulos and J. Boissonade, Physica D 49 (1991) 161. [23] J.J. Perraud, A. De Wit, E. Dulos, P. De Kepper, G. Dewel and P. Borckmans, Phys. Rev. Lett. 71 (1993) 1272; V. Dufiet and J. Boissonade, J. Chem. Phys. 96 (1992) 664. [24] A. De Wit, G. Dewel and P. Borckmans, Phys. Rev. E 48 (1993) R4191.