Two dimension spatial pattern formation in a coupled autocatalysis system

Accepted Manuscript Two dimension spatial pattern formation in a coupled autocatalysis system Zhang Li, Wei Sheng Chen PII: DOI: Reference:

S0307-904X(14)00191-7 http://dx.doi.org/10.1016/j.apm.2014.04.021 APM 9963

To appear in:

Appl. Math. Modelling

Received Date: Revised Date: Accepted Date:

9 March 2012 28 March 2014 7 April 2014

Please cite this article as: Z. Li, W. Sheng Chen, Two dimension spatial pattern formation in a coupled autocatalysis system, Appl. Math. Modelling (2014), doi: http://dx.doi.org/10.1016/j.apm.2014.04.021

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Two dimension spatial pattern formation in a coupled autocatalysis system 1 Zhang Li

2 3

Wei Sheng Chen

(School of Mathematics and Statitics Xidian University XI’AN ShaanXi PR. China 710071)

Abstract:

This paper addresses the cubic autocatalator kinetics modeling of

coupling via diffusion interchange of autocatalyst . By incorporating the effect of two identical cells, each governed by cubic autocatalator kinetics, considering the possibility of the spatiotemporal structures of two dimensional Turing patterns, a new model is proposed. Unlike previous models, the proposed model has two dimensional spatial variation. The equations and the local stability is obtained by linearizing about the spatially uniform solutions,firstly. It is shown that the necessary condition for the model undergoes bifurcation by using the singular purterbation theory. Next Landau constant and amplitude functions of two dimensional Turing patterns consisting of rhombic arrays of rectangles and hexagonal is obtain by singular perturbations theory. Finally, by the method of computer simulation of the model, we describe two different patterns. Keywords:

Turing pattern, nonlinear stability analysis, bifurcation,coupled sys-

tem,singular perturbation . MSC(2010): 35B20,35B35,35B36,35B40 1. Introduction Generally, most dynamics in chemical reaction are nonlinear. However, when systems approach thermodynamics steady state, their kinetics behavior are studied approximately by linear non-equilibrium thermodynamics. When systems are far away thermodynamics steady state, sometimes, nonlinear effects become principal factor of dynamics behavior, the nonlinear behavior couple by diffusive linear behavior. The coupling initially may give rise to the spontaneous appearance of the order and chaotic pattern.

Pattern formation in biological contexts has been studied by Nagorck et

al.[17] and Cruywagen and Murray [18]by coupling together a reaction diffusion sys-

Preprint submitted to Elsevier

April 23, 2014

tem governing the concentration of a pair of morphogens and a cell traction model of the elastic properties of the dermis. Each mechanism is capable of generating spatial patterns on its own, and the coupling leads to the formation of spatial patterns that are more complicated than those seen in either model separately and which are closer to those actually observed . Showalter and coworkers[19] have performed experiments and mathematically modeled the effect of coupling on traveling waves in chemical systems. In these experiments the reacting medium is the Belousov-Zhabotinskii (BZ) reaction in which an inhomogeneous medium is created by removing the ferroin (catalyst) from the bulk solution and immobilizing it in a Nafion ion-exchange membrane. This is then bathed in a solution of the remaining reactants and thus the reaction is restricted to the interfaces between the membrane and the bulk solution. Coupling occurs due to diffusive transport through the membrane of nonionic species such as HBrO2 and gives rise to the formation of new complex spatiotemporal structures. The feature common to nearly all isothermal oscillatory reactions is autocatalysis. Autocatalysis is shared not only by the Belousov-zhabotinskii (BZ)reaction, a wide range of halide-based oscillations, and the arsenite plus iodate reaction but by iodate-arsenous acid reaction. All mechanistic schemes include autocatalytic steps (or autocatalytic combinations of elementary steps). Forms of these schemes need to given a simpler prototypes in a logically ordered way. Singularly perturbed equations are often used as mathematical models describing processes in physics, chemical kinetics, and mathematical biology, and they often arise during investigation of applied problems of technology and engineering. As a class of elementary reaction that can replicate themselves and produce complex chemical morphogenesis and pattern,cubic autocatalysis has formed a great threat to human society and been shown to form the key steps in a number of more complex reaction schemes, see, for example, [1-3]. The model of cubic autocatalysis offer scientist in such far ranging fields as biology, physics, social science and mathematics, an important example of the phenomenon of critical behavior on the pattern and perturbed model. V.V. G and B.Y.D [4] investigate pattern formation in a fractional reaction -diffusion system. By the method of computer simulation of the model of excitable media with cubic nonlinearity they are able to show structure formation in the system with time and space fractional derivatives.In many epidemiological models, the corresponding incidence rate is nonlinear with respect to the numbers of suscepti-

2

ble and infective individuals. The mechanism of cubic autocatalysis is one of the main mechanism of forming nonlinear term of biological model.In [5],Yang et.al. include stochastic perturbations into SIR and SEIR epidemic models with saturated incidence and investigate their dynamics according to the basic reproduction number. In [6], certain two component reaction diffusion systems on a finite are known to possess mesa (box-like) steady state patterns in the singularly perturbed limit of small diffusivity for one of the two solution components. In [7], Stochastic partial differential equations are introduced for the continuum concentration fields of reaction diffusion systems.Spatially adaptive stochastic numerical methods are developed for approximation of the stochastic partial differential equations. Pattern formation in evolutionary dynamics has been studied as an example of phenomena inherent to many natural system. It is well known that reaction-diffusion models are used to study self-organization phenomena in physical, chemical and biological systems. In [8], two SIR epidemic models with different patterns of recruitment and difference in immunity are investigated. In [9], stable Turing patterns are presented in the simplest reaction diffusion system containing a single autocatalitic step in a continuously fed unstirred reactor. In [10], applying two types of Lyapunov functional techniques to an SIRS epidemic model with graded cure and incomplete recovery rates, Y.M. et.al. establish complete global dynamics of the model. In this paper, we consider the prototype chemical reaction scheme based on the cubic autocatalator which is proposed in [2] P →A A + 2B → 3B B→C

(rate k0 p),

(1)

(rate k1 ab2 ),

(2)

(rate k2 b),

(3)

where a and b are the concentrations of the species A and B respectively, and k1 , k2 are the rate constants(both reaction are assumed to be isothermal). Here we supposed that the reaction is taking place inside a closed system with the reactant A now being replenished by the slow decay of a precursor P , where p is the concentration of P and k0 is the rate constant. The reaction diffusion system has been discussed in some detail for the strip spatial variation patterns in one dimension system. The formation of spatial patterns is considered in [11], two types of one dimension patterns occur, standing wave patterns arising out of Hopf bifurcations, together with steady wave patterns arising out of pitchfork bifurcations.

In [12],we consider spatiotemporal structures and

bifurcation arising in two identical cells,which are governed by higher autocatalator 3

kinetics and coupled via diffusive interchange of autocatalyst. In the present paper, we consider a coupled system whereby two identical cells are coupled via the diffusion interchange of the autocatalyst B. In [13],we have considered one dimensional spatial diffusion model and given one strip pattern. Further, we consider the spatial variation to two dimensions and measure distance in these dimensions with the coordinate x and y. The vessel is taken to have length 1 and width 1 with impermeable boundaries at the ends x = 0, 1 and y = 0, 1.The equations are as in [11] but now the effect of spatial variation y included, namely  ∂a1 = D∇22 a1 + µ − a1 b21 ,  ∂t   ∂b 1 = Dδ∇22 b1 + a1 b21 − b1 + β(b2 − b1 ), ∂t (4) ∂a2 = D∇22 a2 + µ − a2 b22 ,   ∂t  ∂b 2 = Dδ∇22 b2 + a2 b22 − b2 + β(b1 − b2 ), ∂t 2

2

∂ ∂ Here, ∇22 = ∂x 2 + ∂y 2 , D and δ are constant diffusion coefficients and t is the time, and β is the dimensionless coupling parameter ,(a1 , b1 ) and (a2 , b2 ) are dimensionless concentration in cells 1 and 2.Consider the two dimensional domain spatial region defined by B = {(x, y)|0 < x < 1, 0 < y < 1}, whose rectangular boundary, we denote by ∂B. In addition to equations (4), we have the boundary conditions associated with the impermeable walls, that

∂bi ∂ai = 0(i = 1, 2), (x, y) ∈ ∂B (5) = ∂n ∂n and suitable initial conditions. It is assumed that the cells are sufficiently thin so that transverse diffusion across them can be considered to be instantaneous. Thus we have two identical regions, divided by a semipermeable membrane which allows the pass of autocatalyst B only, with the some reaction taking place in each region.We find that model which coupled by autocatalator have a unique spatial homogeneous steady state, a1 = a2 = µ? =

1 , b1 = b2 = µ, µ

µ > 0.

(6)

The steady state is stable provided that µ > 1, and unstable if µ < 1,for the unstirred system (D 6= 0, δ 6= 0). This paper is intended to develop an reaction diffusion model of two dimensional spatial diffusion that can effect on the stability of the steady state and pattern. By incorporating new singular perturbation method , we consider the stability of the steady state and give a necessary condition for pattern formation. The proposed model is proved to have the local behavior of two dimensional Turing patterns consisting of rhombic arrays of rectangles and hexagonal arrays of spots or nets,respectively. With the aid of new singular perturbations method, Landau constants and the stability of the equilibrium points of the amplitude equations, asymptotically solution is given. 2. The Linearized theory The equilibrium points (6) to our model system (4) represents a uniform steady state of spatially homogeneous exact solution . It was the stability of this solution to one dimensional 4

perturbations with which were studied by Merkin[13].In this section, we consider the linearized theory which was proposed in [16]. Here, we consider the temporal stability of the uniform state (6) for two dimensional spatial diffusion model. When subject to a disturbance of small amplitude with 0 < ε << 1, we suppose that the disturbance is imposed at t − 0 and is represented by the initial conditions 1 µ

+ ε¯ ai0 (x, y) µ + ε¯bi0 (x, y).

ai (x, y, 0) = bi (x, y, 0) =

(7)

wherei = 1, 2, a ¯i0 , ¯bi0 are bounded functions on B. We consider a solution of equations (4) in the form, ai =

1 + εai0 + O(ε2 ), µ

bi = µ + εbi0 + O(ε2 ), (i = 1, 2).

(8)

On substituting the above expressions into system (4), expanding and neglecting terms of O(ε2 ),we obtain the following linear systems governing the leading terms a10 , b10 , a20 , b20 ,at O(ε), 2 2 ∂a10 = D( ∂∂xa210 + ∂∂ya210 ) − µ2 a10 − 2b10 , ∂t ∂b10 ∂t ∂a20 ∂t ∂b20 ∂t

= = =

∂ 2 b10 2 ∂y 2 ) + µ a10 + b10 + 2 2 D( ∂∂xa220 + ∂∂ya220 ) − µ2 a20 − 2b20 , 2 ∂ 2 b20 2 Dδ( ∂∂xb20 2 + ∂y 2 ) + µ a20 + b20 + 2

Dδ( ∂∂xb10 2 +

β(b20 − b10 )

(9)

β(b10 − b20 ).

subject to zero-flux boundary conditions ∂bi0 ∂ai0 = 0, (i = 1, 2), = ∂n ∂n

(x, y) ∈ ∂B.

(10)

To solve this system of equations subject to the boundary conditions (10), we first define w(x, ~ y) to be the time independent solution of the spatial eigenvalue problem defined by 522 w ~ + k¯2 w ~ = 0,

(n · 5)w ~ = 0,

for (x, y) ∈ ∂B,

(11)

the eigenfunctions of which are ~ m,n cos(nπx) cos(mπy), k¯2 = π 2 (n2 + m2 ), w(x, ~ y) = C where n and m are integers, k¯ is called the wave-number and ω(ω =

2π ¯ ). k

1 ¯ k

is proportional to the wave length

From now on we shall refer to k in this paper as the wave-number, with finite domains

there is a discrete set of possible wave-numbers. Let w ~ k (x, y) be the eigenfunction corresponding to the wave-number k =

√

¯ Each eigenDk.

function w ~ k satisfies zero flux boundary conditions. Because the problem (9) is linear, we now look for solutions a1 (x, y, t), b1 (x, y, t) of (9) in the form  

 ai0 (x, y, t) bi0 (x, y, t)

=

∞ X k=0

 

 Ai0 k Bki0

5

 eσk t w ~ k (x, y), (i = 1, 2),

(12)

i0 where the constants Ai0 k , Bk are determined by the initial conditions and σk is the eigenvalue

which determines temporal growth. Substituting (12) into (8), we see that (12) have a non-trivial solution for Aik , Bik provided σ satisfies, det(M ) = 0. where     M =   



σk + kr2 + µ2

2

0

0

−µ2

σk + δkr2 − 1 + β

0

−β

0

0

σk + kr2 + µ2

0

2

−β

−µ

2 σk +

δkr2

   ,   

−1+β

√ with kr = rπ D(r = 0, 1, 2, · · · , ∞). On expanding this determinant, we find a dispersion relation which gives σk in terms of kr , µ, giving f (σk )g(σk ) = 0

(13)

with f (σk ) = σk2 + p1 (µ, kr )σk + q1 (µ, kr ), g(σk ) = σk2 + p2 (µ, kr )σk + q2 (µ, kr ), where p1 (µ, kr ) = (δ + 1)kr2 + µ2 − 1,

q1 (µ, kr ) = δkr4 + (δµ2 − 1)kr2 + µ2 ,

p2 (µ, kr ) = (δ + 1)kr2 + µ2 − 1 + 2β,

q2 (µ, kr ) = δkr4 + (δµ2 − 1 + 2β)kr2 + (1 + 2β)µ2 ,

We note that f (σk ) corresponds to the equation for the linear stability of the uncoupled system, g(σk ) is that of the coupled system and that the roots of these equations are given by q 1 ± (14) σik = (−pi ± p2i − 4qi ), (i = 1, 2) 2 The local stability of the steady state (6) is then determined by examining the behavior ± ± of σ1k and σ2k in the first quadrant of the (kr , µ) plane. So we discuss the neutral curves ± ± subjected to σki = 0 and Re(σki ) = 0 in the first quadrant of the (kr , µ) plane, which

are determined by f (σk ) = 0 and g(σk ) = 0. . The determination of the neutral curves is a necessary conditions to determine where the spatially uniform steady state (6) bifurcated from stable spatially nonuniform steady state. In [12], the equation of the neutral curves was examined in some detail and that discussion need not be repeated here. In [12], it was shown that the neutral curves,which we label L1 , corresponding to f (σk ) = 0 and L2 corresponding to g(σk ) = 0.That is L1 is given by   1 − (δ + 1)k 2 , 0 ≤ kr ≤ k ? , r r q µ2c1 (kr ) =  kr2 (1−δkr2 ) , ? kr ≤ kr ≤ 1δ . 1+δk2 r

6

(15)

± + ) = 0 and σk1 = 0 intersect where (kr? , µ? ) is the point on the (kr , µ) plane where the curves Re(σk1

and is given by

√

kr?2

=

1 δ2 [

µ?2

=

1 − (δ + 1)kr?2 ,

δ 2 + 1 − 1],

± The first of these branches 0 < kr ≤ kr? corresponds to Re(σk1 ) = 0, while the second corresponds q ± 1 to σk1 = 0. Note that L1 intersects the kr axis at kr = δ . ± ± The curves L2 on which σk2 = 0 and Re(σk2 ) = 0 are, respectively,   1 − 2β − (δ + 1)k 2 , 0 ≤ kr ≤ k¯r , r q µ2c2 (kr ) =  kr2 (1−2β−δkr2 ) , k¯ ≤ k ≤ 1−2β . r

1+δkr2 +2β

r

(16)

δ

± + where (k¯r , µ ¯) is the point where the curves Re(σk2 ) = 0 and σk2 = 0 intersect and is given by

p

2 k¯r

=

µ ¯2

2 = [1 − 2β − (δ + 1)k¯r ],

1 δ2 [

1 + δ 2 + 4βδ − (1 + 2βδ)],

We can see that L2 lies in the positive quadrant of the (kr , µ) plane only if 0 < β < 12 . A sketch of the neutral curves L1 and L2 is shown in Fig1. When the condition is satisfied, L2 lies strictly ± ± below L1 . Any point above the neutral curve L1 in Fig1 has Re(σki ) < 0 or σki < 0, (i − 1, 2),

whilst all points below it have two cases: ± ± 1.Any point above the neutral curve L2 has Re(σki ) < 0 but the point has at least Re(σk1 ) > 0. ± ± 2.The point below the neutral curve L2 has either Re(σk2 ) > 0, Re(σk1 )± > 0) or σk2 > 0 and ± σk1 > 0.

√ p=

2

µ

6

√ √ 1+2β( 2− 1+2β) δ

q= √ ( 2−1) √ δ

p1 =

√ √ ( 2− 1+2β) √ δ

1

1−2β δ √ 2−1 δ

q1 = 1 L1

M

1 δ

M = (kr? , µ?r ) N = (k¯r , µ ¯r )

À L2

1 − 2β N 0

p1 p

q1

q

kr

Fig.1.The neutral curves L1 and L2 with 0 < β < 12 ,

For a given value of µ, the occurrence of stable patterns requires that the spatially uniform steady state (6) should be stable in the absence of diffusion. This requires 7

± ± Re(σki (kr , µ)) < 0 or σki (kr , µ) < 0 for each r = 0, 1, 2 · · · ∞ and leads to the necessary

condition for pattern formation that is provided by   1, δ ≥ 1, 2 µ2 > µM c (δ, D) =  max{1, µ2 , µ2 }, 0 < δ < 1, 1 2

(17)

where µ21 µ22 N

= = =

2 2 kN (1−δkN ) , 2 1+δkN  2 2 k (1−δk  N +1 2 N +1 ) ; 1+δkN +1



0;

Integer part of{ π1

q

N≤

√ √1 π δD

− 1,

otherwise, 1 δλ }.

One branch of the neutral curves L1 and L2 will lead to Hopf bifurcation which can give rise to spatial and temporally varying periodic formation and others is Turing bifurcation which can lead to time independent spatially periodic formation. We mainly consider the spatially pattern formation. Therefore, for each µ > µM c (δ, D), the leading order terms (ui0 , vi0 ) in (8) decay exponentially as t → ∞, and any small disturbance imposed upon (6) is stable. However, for 0 < µ < µM c (δ, D), the terms in (8) are grow exponentially as t → ∞, and the system will diverge from the uniform state (6) when perturbed. Pattern formation may arise through a Turing instability for certain ranges of parameter values. A necessary condition for this to occur is that + ± the maximum on the σk1 = 0 curve should lie above the maximum on the Re(σk1 )=0 √ √ 2 ( 2− 1+2β) , 0 < β < 12 . The unstable curve, which gives us the condition that 0 < δ < 1−2β q 1 . When β > 12 , the curve L2 does not exist and modes are those which satisfy r < π1 Dδ

the only neutral modes which occur are identical to those in the uncouple system; this is called the strong coupling regime in [2]. When 0 < β < 12 , the system is called the weak coupling regime, the coupling gives rise to neutral modes which do not occur in the uncoupled system, corresponding to points on L2 . These all correspond to values of µ at which the uniform state is unstable, having lost stability at the larger value of µ corresponding to one of these modes on L1 . Any nonuniform steady states which bifurcate from these modes will be unstable (at least locally). Thus weak coupling gives rise to steady state solutions which do not occur in the uncoupled system. The most dangerous modes is kr2 =

1−δµ2 2δ .

8

3. Two dimensional pattern As we have already seen, pattern formation may arise through a Turing instability for certain ranges of parameter values. A necessary condition for this to occur is that the maximum on the σ+ = 0 curve should lie above the maximum on the Re(σ) = 0 curve, which give us the condition √ that 0 < δ < 3 − 2 2, 0 < β < 12 . Turing instability occurs at the largest value of µ corresponding to one of these neutral models. In order to investigate the possibility of occurrence for our model of those two dimensional patterns, we begin our discussion by putting µ = µc − µ1 ε2 ,

0 < ε << 1,

(18)

where µ1 is of O(1) , depending on whether the pattern form bifurcates initially into µ < µc (µ1 > 0) or µ > µc (µ1 < 0) and µc is given by µ2c =

kc2 (1 − δkc2 ) kc2 (1 − 2β − δkc2 ) 2 or µ = . c 1 + δkc2 1 + 2β + δkc2

(19)

where kc is the critical modes. We then proceed with the solution, as in [11], by looking for a solution by expanding about stationary states in powers of ε

ai ∼

1 µ

bi ∼

µ + εbi1 + ε2 bi2 + ε3 bi3 + · · · .

+ εai1 + ε2 ai2 + ε3 ai3 + · · · ,

(i = 1, 2)

(20)

On substitution of the (18-20) into (4) and equating like powers of ε, we obtain the recursive series of equations Lc [a11 , b11 , a21 , b21 ]> = 0,

(21)

Lc [a12 , b12 , a22 , b22 ]> = f2 (a11 , b11 , a21 , b21 ),

(22)

Lc [a13 , b13 , a23 , b23 ]> = f3 (a11 , b11 , a21 , b21 , a12 , b12 , a22 , b22 ),

(23)

with zero flux boundary conditions and where the operator Lc is given by  ∂ − D∆ + µ2c 2 0 0  ∂t  ∂  −µ2c 0 −β ∂t − Dδ∆ − 1 + β Lc =   ∂ 2  0 0 2 ∂t − D∆ + µc  0 where ∆ =

∂2 ∂x2

+

−µ2c

−β

∂ ∂t

    .   

− Dδ∆ − 1 + β

∂2 ∂y 2

3.1 A rhombic planform analysis In order to investigate the possibility of occurrence for our model of those rhombic type patterns, 9

we shall consider a solution of the leading term a1 , b1 of the form 

 a11

   b11    a21  b21



 A11

      B11 =     A21   B21



 C11

     σt  D  e cos(qc x) +  11     C21   D21

   σt  e cos(qc z),   

where z = cos(ψ)x + sin(ψ)y, ψ is characteristic angle of a rhombic pattern, qc =

(24)

kc √ D

is the critical

wave-number of linear stability theory. The local stability of the patterns solution may be determined by the introduction of the slow time scale τ = ε2 t with the amplitude Ai1 , Bi1 , Ci1 , Di1 (i = 1, 2). On substituting (24) into (21), the amplitude Ai1 , Bi1 and Ci1 , Di1 satisfy the same equation which is given by



 A11

   B11 Mc    A21  B21



 A11

      B  = −σ0  11     A21   B21

   ,   

(25)

where det(Mc − σ0 I) ≡ F (σ0 , kc , µc )G(σ0 , kc , µc ) = 0. Here

    Mc =    

(26) 

kc2 + µ2c

2

0

0

−µ2c

δkc2 − 1 + β

0

−β

0

0

kc2 + µ2c

2

0

−β

−µ2c

δkc2 − 1 + β

   .   

We may obtain σ0 = 0 or σ0 = −p < 0 with linearized theory, without general loss, the solution of (24) is given 

 a11

   b11    a21  b21



 d11

      e  = A(τ )  11     d21   e21



 d11

      e  cos(qc x) + B(τ )  11     d21   e21

    cos(qc z),   

σ0 = 0,

(27)

where d11 /e11 = d21 /e21 = −2/(kc2 + µ2c ), µc is given by (19). The amplitudes A and B are functions of the slow variable τ . To proceed, the equations at o(ε2 ) are most conveniently written

10

in vector form as



 a12

   b12 Lc    a22  b22



− µ1c b211 − 2µc a11 b11

    1 2   µc b11 + 2µc a11 b11 =     − µ1 b221 − 2µc a21 b21 c   1 2 µc b21 + 2µc a21 b21

    .   

(28)

We require only the particular integral part of the general solution of (27), since the homogeneous part can be absorbed into the leading order term by suitable modification of A and B. The appropriate solution of equations (28) can be written, after some algebra, as 

 a12 (x, y, t, τ )

   b12 (x, y, t, τ )    a22 (x, y, t, τ )  b22 (x, y, t, τ )     p22    q12  q22

 p11

      p21 =     q11   q21

 p12

where





 p13

      p  AB cos qc (x + z) +  23     q13   q23 

 p1i

   p2i    q1i  q2i

   2  [A cos(2qc x) + B 2 cos(2qc z)]+    

 p14

      p  AB cos qc (x − z) +  24     q14   q24



 R1i

      R 1  = [p  2i  Fi    R3i   R4i



   2  (A + B 2 )   

 R3i

      R  + q  4i     R1i   R2i

   ],   

with p=

1 2 1 2 e + 2µb d11 e11 , q = e + 2µb d21 e21 , µb 11 µb 21

Fi = 2(ξi2 − αi2 β 2 ),

ξi = θi αi + 2µ2b ,

R1i = αi β 2 + 2ξi − θi ξi ,

R2i = ξi (αi − µ2b ),

R3i = 2β(µ2b − αi ),

R4i = αi β(αi − µ2b ). (i = 1, 2, 3, 4)

and α(kb ) = kb2 + µ2b ,

θ(kb ) = δkb2 − 1 + β

α1 = α(2kb ), θ1 = θ(2kb ), √ √ α2 = α( 2 + 2 cos ψkb ), θ2 = θ( 2 + 2 cos ψkb ), √ √ α3 = α( 2 − 2 cos ψkb ), θ3 = θ( 2 − 2 cos ψkb ), α4 = α(0)

θ4 = θ(0).

, 11

(29)

At O(ε3 ) we again have a linear inhomogeneous problem with the same left hand sides as before, 

 a13

   b13 Lb    a23  b23





11 − ∂a ∂τ

    ∂b11   − ∂τ =   ∂a21   − ∂τ   21 − ∂b ∂τ



−( γ1 a11 + R1 I1 )

    1   γ a11 + R1 I1 +     −( γ1 a21 + R2 I2 )   1 γ a21 + R2 I2

       

(30)

where R1 I1 = 2µb (a11 b12 + a12 b11 ) + a11 b211 +

2 µb b11 b12

R2 I2 = 2µb (a21 b22 + a22 b21 ) + a21 b221 +

2 µb b21 b22

To obtain equations for the amplitudes A and B of the leading order solutions, we invoke the Fredham alternative for the solvability of equations(24). After some algebra, we find that this Fredham solvability condition requires that γA − A(D˜1 A2 + D˜2 B 2 ),

E dA dτ

=

E dB dτ

= γB − B(D˜2 A2 + D˜1 B 2 ),

(31)

where γ = −2µb µ1 , E > 0. Then evaluating the Landau constants D˜1 , D˜2 at σ0 = 0, it follows that , D˜1 =

d11 D1 + d21 D3 ,

D˜2 =

d11 D2 + d21 D4 .

(32)

where D1 =

d11 [ 34 e211 + (µb +

1 e11 µb d11 )(p21

+ 2p24 ) + µb de11 (p11 + 2p14 )] 11

D2 =

d11 [ 34 e211 + (µb +

1 e11 µb d11 )(p22

+ p23 + 2p24 ) + µb de11 (p12 + p13 + 2p14 )] 11

D3 =

d21 [ 34 e221 + (µb +

1 e21 µb d21 )(q21

+ 2q24 ) + µb de21 (q11 + 2q14 )] 21

D4 =

d21 [ 34 e221 + (µb +

1 e21 µb d21 )(q22

+ q23 + 2q24 ) + µb de21 (q12 + q13 + 2q14 )] 21

(33)

We now turn our attention to the rhombic planform amplitude equations of interest of (31) which possess the following equivalence classes of critical points (A0 , B0 ): (a) A0 = B0 = 0, (c) B02 =

γ , A0 = 0, D˜1

(b) A20 =

γ , B0 = 0, D˜1

(d) A20 = B02 =

γ . ˜ D1 + D˜2

Assuming that D˜1 , D˜1 + D˜2 > 0 and investigating the stability of the critical points we seek a solution of the amplitude equations of the form A ∼ A0 + εA1 eλτ + O(ε2 ), B ∼ B0 + εB1 eλτ + O(ε2 ).

(34)

Because the stability of (b) is the same as that of (c), we only investigate the stability of (a),(b),(d). 12

After substitution from (34) into (31), expanding and neglecting terms of O(ε2 ), we obtain the following the equation governing the eigenvalue of λ , (λ − Eζ)2 + Eb(λ − Eζ) + E 2 c = 0

(35)

where 2 b = (3D˜1 + D˜2 )(A20 + B02 ), c = (3D˜1 A20 + D˜2 B02 )(D˜2 A20 + 3D˜1 B02 ) − 4A20 B02 D˜2 .

We find that λ has the associated roots (a) λ1,2 = Eζ, (b) λ1 = −2Eζ, λ2 = Eζ(1 − (d) λ1 = −2Eζ, λ2 = 2Eζ(

D˜2 ), D˜1

D˜2 − D˜1 ). D˜1 + D˜2

They yield the stability criteria that (a) is stable for ζ < 0; (b), for ζ > 0, D˜1 < D˜2 ; and (d), for ζ > 0, D˜1 > D˜2 . The Landau constants D˜1 , D˜2 is determined by (33), D˜1 = f1 (θ, ψ), D˜2 = f2 (θ, ψ) where θ = Dδqc2 = δkc2 , kc is the critical wavenumber of linear stability theory, To simplify D˜1 , D˜2 further, we need the ratio

e11 d11

and

e21 d21 .

There are two case to consider,

depending on whether the bifurcation lies on L1 or L2 . When the bifurcation point lies on L1 , we find that e11 e21 −2 k 2 (1 − δkc2 ) = = 2 , µ2c = c 2 d11 d21 kc + µc 1 + δkc2 the same as for the uncoupled system, whereas when the bifurcation point lies on L2 , e11 e21 −2 kc2 (1 − 2β − δkc2 ) 2 = = 2 , µ = c d11 d21 kc + µ2c 1 + 2β + δkc2 The Landau constants is given as , D˜1 = D˜1 (x, θ)

=

3 2

−

x(2x−1) f1 (θ) 1−x g1 (θ)

D˜2 = D˜2 (x, θ)

=

3 2

−

x(2x−1) f3 (θ) 1−x [ g3 (θ)

− (2x − 1) fg22 (θ) (θ) − +

f4 (θ) g4 (θ) ]

13

2x−1 2β+3 1−x 1+2β

− (2x − 1)[ fg55 (θ) (θ) +

f6 (θ) g6 (θ) ]

−

2(2x−1) 2β+3 1−x 1+2β

where f1 (θ)

=

4(1 + x),

g1 (θ)

=

20θ − 3 + (12θ − 5)x,

f2 (θ)

=

(3 − 4θ)[(12θ − 5)x + 20θ − 3] + 2β[(−12θ + 3)x − 20θ + 15],

g2 (θ)

=

[20θ − 3 + (12θ − 5)x][20θ − 3 + 5 · 2β + (12θ − 5 + 3 · 2β)x],

f3 (θ)

= f3 (θ, s) = s(1 + x),

g3 (θ)

=

f4 (θ)

= f3 (θ, r),

g4 (θ)

=

f5 (θ)

= f5 (θ, s) = (3 − sθ)f3 (θ) + 2β(1 − sθ)[s + 1 + (s − 1)θ],

g5 (θ)

=

f6 (θ)

= f5 (θ, r),

g6 (θ)

=

g3 (θ, s) = (s + 1)sθ − s + 1 + [(s − 1)sθ − s − 1]x,

g3 (θ, r),

g5 (θ, s) = g3 (θ)[(s + 1)(sθ + 2β) − s + 1 + ((s − 1)(sθ + 2β) − s − 1)x],

g5 (θ, r),

withs = 2 + 2 cos ψ, r = 2 − 2 cos ψ. and x = θ (uncoupled), x = 2β + θ (coupled) Note that (a) and (b), as in the one dimensional analysis of Merkin, represent the homogeneous and striped states, respectively, while (d) can be identified with a rhombic pattern possessing characteristic angle ψ. Towards this end, we examine the signs of D˜1 and D˜2 , D˜1 + D˜2 , D˜2 − D˜1 for 0 < θ < 1 and 0 < ψ ≤ π. We first illustrate this procedure for a fixed value of ψ, namely ψ =

π 2,

by plotting

D˜1 + D˜2 and D˜2 − D˜1 , D˜1 , D˜2 versus θ in Fig. 2(a). Here the two θ − interval of stable square rhombic patterns are those which satisfied D˜1 + D˜2 > 0, D˜2 − D˜1 < 0. We observe that for the θ values between these intervals D˜1 > 0, D˜2 > D˜1 > 0 the system will exist stable stripes pattern. While for those to the left of the left-hand interval and to the right of the right-hand one no pattern is predicted by this analysis. Repeating this process for ψ = 2(b).

14

π 4

yields the companion plot of Fig.

30

25

20

15

10

5

0

−5

−10

−15

−20 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.7

0.8

0.9

1

0.7

0.8

0.9

1

Fig.2.(a) 140

120

100

80

60

40

20

0

−20

−40 0.2

0.3

0.4

0.5

0.6

Fig.2.(b) 140

120

100

80

60

40

20

0

−20

−40 0.2

0.3

0.4

0.5

0.6

Fig.2.(c)

15

30

25

20

15

10

5

0

−5

−10

−15

−20 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Fig.2.(d) Fig.2 For given ψ =

π 4

˜ D2, ˜ D2 ˜ − D1, ˜ D2 ˜ + D1 ˜ and β = 0 , the curves D1,

400

300

200

100

0

−100

−200

−300

−400 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Fig.3.(a) 600

400

200

0

−200

−400

−600

−800

−1000 0.2

0.3

0.4

0.5

0.6

Fig.3.(b)

16

0.7

0.8

0.9

1

400

300

200

100

0

−100

−200

−300

−400 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.6

0.7

0.8

0.9

Fig.3.(c) 400

300

200

100

0

−100

−200

−300

−400 0.1

0.2

0.3

0.4

0.5

Fig.3.(d) π 6

Fig.3 For given ψ =

˜ D2, ˜ D2 ˜ − D1, ˜ D2 ˜ + D1 ˜ and β = 0.6 , the curves D1,

Next we plot our Landau-curves for the fixed value of θ = 0.5 versus 0 ≤ ψ ≤ π in Fig. 4. 3000

2000

1000

0

−1000

−2000

−3000

−4000

0

0.5

1

1.5

2

2.5

3

3.5

˜ Fig.4 For given θ = 0.5,the curves D2 17

Restricting ourselves to the interval of interest 0 ≤ ψ ≤ π 3,

stable rhombic patterns flanking ψ =

π 2,

we see that there are two bands of

which have again been designated by shading, with no

pattern between these bands and stable stripes outside of them. This figure has been drawn for the extended interval

π 2

≤ ψ ≤ π in order to demonstrate graphically the symmetry about ψ =

π 2

characteristic of rhombic patterns. Repeating this procedure for other values of θ in the relevant interval, we observe from Figures 3 that there exist no stable rhombic patterns of characteristic angle ψ =

π 3

by virtue of the fact that lim D˜2 (ψ, θ) → −∞.

(36)

ψ→ π 3

This behavior occurs as a consequence of D˜2 containing terms proportional to the components of F3 , the denominator of which satisfies π F3 ( , kb , µb ) = 0, 3 where F3 (ψ, kb , µb ) = [(2δ(1 − cos ψ)kb2 − 1 + β)(2(1 − cos ψ)kb2 + µ2b ) + 2µ2b ]2 − [2(1 − cos ψ)kb2 + µ2b ]2 β 2 . These results of (36) are most easily deduced upon noting that π F ( , kb , µb ) = [(δkb2 − 1 + β)(kb2 + µ2b ) + 2µ2b ]2 − [kb2 + µ2b ]2 β 2 . 3 in conjunction with the equivalent relations µ2b =

k 2 (1 − 2β − δkb2 ) kb2 (1 − δkb2 ) , or µ2b = b 2 1 + 2β + δkb2 1 + δkb

respectively. We close this section with a morphological instability interpretation of the potentially stable critical points (b) and (d) of our amplitude equations when δ > 0 (µ < µc )relative to the Turing patterns under investigation. Then to lowest order the solution associated with these critical points is given by  a1 −

1 µ

   b1 − µ    a2 − µ1  b2 − µ where λc =







 d11

      e  ∼ A0  11     d21   e21

 d11

    ¶ µ  e11  2πx  cos + B0    λ c  d21    e21

  ¶ µ   cos 2πz ,  λc  

(37)

2π qc .

To follow the solutions away from their bifurcation points, we used the center finite difference method and alternating direction implicit method to give the contour plot 18

for function with A0 > 0 and B0 = 0 relevant to critical point (b)in the x − y plane of Fig. 5. Here the spatial variables are measured on units of λc . The initial conditions are randomly selected on the neutral curve. The colors scale is amplitude phase is measured on units of λc .Clear such alternating light and dark parallel bands produced by this critical point should be identified with a striped Turing pattern as anticipated above. In order to make an analogous interpretation of critical point (d), we consider the deviation function (37) with A0 = B0 > 0 and allow ψ to take on both the values employed in Fig. 2. We represent the contour plot for that function with ψ =

π 2

in Fig. 6. From the structure of the latter it is equally

clear that this critical point should be identified with a Turing pattern of square planform. Finally in a similar manner we generate the contour plot of Fig. 7. for ψ =

π 4

, which forms a family of

rectangles. −6

−6

x 10

x 10

1.5

1.5

1

1

0.5

0.5

0

0

−0.5

−0.5

−1

−1

−1.5

−1.5 −1

0

1

−1

0

1

−6

−6

x 10

x 10

−6

−6

x 10

x 10

1.5

1.5

1

1

0.5

0.5

0

0

−0.5

−0.5

−1

−1

−1.5

−1.5 −1

0

1

−1 −6

0

1 −6

x 10

x 10

Fig.5 The striped Turing pattern for critical point (b)

19

−9

−9

x 10

x 10

4

4

2

2

0

0

−2

−2

−4

−4 −2

0

2

−2

0

2

−9

−9

x 10

x 10

−9

−9

x 10

x 10

4

4

2

2

0

0

−2

−2

−4

−4 −2

0

2

−2

0

2

−9

−9

x 10

x 10

Fig.6 The Turing pattern of square planform for critical point (d) with parameter value ψ = −6

π 2

−6

x 10

x 10

1.5

1.5

1

1

0.5

0.5

0

0

−0.5

−0.5

−1

−1

−1.5

−1.5 −1

0

1

−1

0

1

−6

−6

x 10

x 10

−6

−6

x 10

x 10

1.5

1.5

1

1

0.5

0.5

0

0

−0.5

−0.5

−1

−1

−1.5

−1.5 −1

0

1

−1 −6

0

1 −6

x 10

x 10

Fig.7 The Turing pattern of rectangles planform for critical point (d) with parameter value ψ =

π 4

To demonstrate that the same thing happen with any ψ, we first put (37) for critical point (d)

20

in the form  a1 − µ1    b1 − µ    a2 − µ1 





 d11

      e  ∼ A0  11     d21   b2 − µ e21



 d11

  ·  µ ¶ µ ¶¸   e11 2πx 2πz  cos  + cos = 2A 0  λc λc   d21   e21

    cos ω1 cos ω2 ,   

(38)

where ω1 + ω2 =

2πx 2πz , ω1 − ω2 = , λc λc

(39)

or ω1 =

π π (x + z) = [(1 + cos ψ)x + sin(ψ)y], λc λc

(40)

ω2 =

π π (x − z) = [(1 − cos ψ)x − sin(ψ)y]. λc λc

(41)

and

From (39) to (41), we can then deduce that the intersecting level curves in the associated contour plot are two families of straight lines possessing slops of k1 =

−[1 + cos(ψ)] [1 − cos(ψ)] , k2 = sin(ψ) sin(ψ)

with k1 k2 = −1. Given that in general this is a rectangular planform we need to explain in what sense such a critical point can be identified with a rhombic pattern. To do so we form the quadrilateral depicted by dashed lines in Fig. 6. Its sides are each composed of two half diagonals collectively contained in the four light rectangles surrounding a dark one. Thus each side of that quadrilateral has the length of one of these diagonals λ = λc / sin(ψ). Each member of that family has aspect ratio W/L = tan(ψ/z), where W and L are its width and length, respectively, while ψ/z serves as its angle of inclination as well, an assertion most easily verified by the relation µ ¶ · ¸1 · ¸1 ψ 1 − cos(ψ) 2 (1 − cos(ψ))2 2 1 − cos(ψ) tan = = = = k2 . 2 2 1 + cos(ψ) 1 − cos (ψ) sin(ψ) Therefore we give such a Turing pattern as a rhombic array of rectangles. 3.2

A hexagonal planform analysis In order to investigate the possibility of occurence for our model of those hexagonal patterns,

21

we consider a solution of the leading term (a11 , a21 , b11 , b21 ) of the form 

 a11

   a21    b11  b21







A11 (t)

 A12 (t)

     √        A21 (t)   A22 (t)  =  [cos(qc x − φ1 (t)] +   cos[qc x − 3y − φ2 (t)]      2   B11 (t)   B12 (t)       B21 (t) B22 (t) 

 A13 (t)

  √    A23 (t)  x + 3y   + − φ3 (t)].  cos[qc 2  B13 (t)    B23 (t)

(42)

To determine the solvability conditions for these Landau constants most easily, we introduce the transformation  A11 (t)    A21 (t)    B11 (t) 





 A11

      A21 =     B11   B21 (t) B21



 A12 (t)

   σt e ,   



 A13 (t)



 A12

           A22 (t)   A23 (t)  1  A22  = =        B12 (t)   B13 (t)  2  B12      B22 (t) B23 (t) B22

   σt e ,   

(43)

φ1 (t) = φ2 (t) = φ3 (t) = 0. In this section, we shall consider parameter of bifurcation µ and the time derivatives which are accordingly replaced using

where qc =

kc √ δ

µ

=

µc − µ0 ε − µ1 ε2 ,

∂ ∂t

=

∂ ∂t

(44)

+ ε ∂τ∂ 1 + ε2 ∂τ∂ 2 ,

is the critical wavenumber of linear stability theory, µ0 is O(ε) and µ1 is O(ε2 ), µc

is given by (13), and the dependence of Ai , Bi (i = 1, 2) on the slow time scale τ1 ≡ εt, τ2 ≡ ε2 t appears explicitly. Substituting (42-44) into (21) and proceeding in exactly the same manner as in the last section, we take σ = σ0 , 

 a11

   b11    a21  b21 where

e11 d11

=



 e11

      d  = A(τ1 , τ2 )  11     e21   d21

−2 kc2 +µ2c ,



 e11

      d  cos(qc x) + B(τ1 , τ2 )  11     e21   d21

 √    cos qc x cos 3qc y ,  2 2  

A and B are function of τ1 , τ2 to be determined.

22

(45)

The equation at o(ε2 ) are most conveniently written in vector form 

 a12

   b12 Lc    a22  b22

 2µc µ0 a11 −

1 2 µc b11

 − 2µc a11 b11

      −2µc µ0 a11 + µ1c b211 + 2µc a11 b11 =     2µc µ0 a21 − µ1 b221 − 2µc a21 b21 c   −2µc µ0 a21 + µ1c b221 + 2µc a21 b21



      −      

∂a11 ∂τ1 ∂b11 ∂τ1 ∂a21 ∂τ1

    .   

(46)

∂b21 ∂τ1

This solution is only possible provided a solvability condition is satisfied, which determines the dependence of A and B on the slow time scale τ1 . After some algebra, we find that this solvability condition require that

  E dA dτ1  E dB dτ1

where γ0 =

1 4µc (1

=

2µc µ0 A − γ0 B 2 ,

=

2µc µ0 B − 4γ0 AB]

,

(47)

+ 2µ2c de11 ), E > 0. 11

The solution of equations (46) can be written 

 a12

   b12    a22  b22

    =   



 p11

   p21    q11  q12  p12    p22    q12 

  √  1 2  [ A cos(2qc x) + 1 B 2 cos(qc x) cos( 3qc y)] +  2 4     √  √   [AB cos( 3qc x ) cos( 3qc y ) + 1 B 2 cos( 3qc y)] +  2 2 4  

q22 

 p13

where

   p23    q13  q23

  1  (2A2 + B 2 ), 2  



 p1i

   p2i    q1i  q2i



 R1i

      R 1  = [p  2i  Fi    R3i   R4i



 R3i

      R  + q  4i     R1i   R2i

   ],   

with p=

1 2 1 2 e + 2µb d11 e11 , q = e + 2µb d21 e21 , µb 11 µb 21

23

Fi = 2(ξi2 − αi2 β 2 ),

ξi = θi αi + 2µ2b ,

R1i = αi β 2 + 2ξi − θi ξi ,

R2i = ξi (αi − µ2b ),

R3i = 2β(µ2b − αi ),

R4i = αi β(αi − µ2b ). (i = 1, 2, 3, 4)

and α(kb ) = kb2 + µ2b ,

θ(kb ) = δkb2 − 1 + β

α1 = α(2kb ), θ1 = θ(2kb ), √ √ α2 = α( 2 + 2 cos ψkb ), θ2 = θ( 2 + 2 cos ψkb ), √ √ α3 = α( 2 − 2 cos ψkb ), θ3 = θ( 2 − 2 cos ψkb ), α4 = α(0)

θ4 = θ(0).

, At O(ε3 ) we again have a linear inhomogeneous problem with the same left hand sides as before,         11 12 a13 − ∂a − ∂a R01 − R1 ∂τ2 ∂τ2              ∂b12    11  b13   − ∂b     − R + R 1  ∂τ2  ∂τ2  =   01  Lb  − + (48)      ∂a22    21   a23   − ∂a     − R − R 02 2 ∂τ2  ∂τ2        ∂b21 ∂b22 b23 − ∂τ2 − ∂τ2 R02 + R2 where R01

= (µ20 − 2µc µ1 )a11 − 2µb µ0 a12 − 2µ0 a11 b11

R02

= (µ20 − 2µc µ1 )a21 − 2µb µ0 a22 − 2µ0 a21 b21

R1

=

2 µb b11 b12

+ 2µb (a11 b12 + a12 b11 ) + a11 b211

R2

=

2 µb b21 b22

+ 2µb (a21 b22 + a22 b21 ) + a21 b221

We can determine the solvability conditions for these the amplitude phase function of A and B, dA E dτ 2

= (2µc µ1 − µ20 )A + 12 µ0 D0 B 2 − A(D˜1 A2 + D˜2 B 2 ),

dB E dτ 2

= (2µc µ1 − µ20 )B + 2µ0 D0 AB − B(2D2 A2 + 14 (D˜1 + 2D˜2 )B 2 ),

(49)

where E > 0 is constant, D˜1 and D˜2 are Landau constants in amplitude equations and are given respectively by, D˜1

=

2(p+q) e11 µc [( d11

+R41 + µ2c )( R214F + 1

R23 +R43 ) 2F3

+R31 + µ2c ( R114F + 1

R13 +R33 )] 2F3

+ 32 ,

D˜2

=

2(p+q) e11 µc [( d11

42 + µ2c )( R22F+R + 2

R23 +R43 ) F3

32 + µ2c ( R12F+R + 2

R13 +R33 )] F3

+ 2.

These amplitude phase equations to A and B are given by (47) and (49), dA dτ

=

E[σA − D˜0 B 2 − A(D˜1 A2 + D˜2 B 2 )],

dB dτ

=

E[σB − 4D˜0 AB − B(2D˜2 A2 + 14 (D˜1 + 2D˜2 )B 2 )],

(50)

where τ = τ1 +τ2 , σ = 2µc (µ1 +µ0 )−µ20 , D˜0 = γ −2µ0 D0 . Having determined the Landau constants of the amplitude equation (50), we now return to the six disturbance hexagonal planform amplitude 24

phase equations (47) and (49). In cataloguing the critical points of these amplitude phase equations and summarizing their orbital stability behavior, we shall employ the quantities, σ−1 =

2 2 2 −4D˜0 16D˜1 D˜0 32(D˜1 + D˜2 )D˜0 , σ1 = , σ2 = . D˜1 + 4D˜2 (2D˜2 − D˜1 )2 (2D˜2 − D˜1 )2

Then, these critical points can be obtained by a

:

A0 = B0 = 0;

b

:

A20 =

c±

:

2A0 = B0 = B0± =

d

:

A0 =

σ , B0 D˜1

= 0, (σ > 0);

−4D˜0 , B02 2D˜2 −D˜1

2[−2D˜0 ±

=

√

2

4D˜0 +(D˜1 +4D˜2 )σ] , (σ D˜1 +4D˜2

4(σ−σ1 ) , (σ D˜1 +2D˜2

> σ−1 );

> σ1 ),

where we are assuming D˜1 > 0, D˜1 + 4D˜2 > 0. In order to examine the stability of these critical points, we seek a solution of the amplitude phase equations (50) of the form A1 (t) =

A0 + εc1 ept + O(ε2 );

Aj (t)

=

B0 + εcj ept + O(ε2 )(j = 2, 3);

φi (t)

=

0 + εci+3 ept + O(ε2 )(i = 1, 2, 3),

with 0 < ε << 1. This sort of stability of pattern formation is meant in the sense of a family of solutions in the plane which may interchange with each other but do not grow or decay to a solution type from a different family. Such an interpretation depends upon the translation and rotational symmetries inherent to the amplitude phase equations, these invariance also limiting each equivalence class of critical points to a single member that must be considered explicitly. We obtained a is stable for σ < 0, while the stability behavior of b and c± which depends upon the signs of D˜0 and 2D˜2 − D˜1 . When b is equivalent to our original one dimensional periodic patterns exhibiting hexagonal symmetry. Finally, critical point d, which reduces to b for σ = σ1 and c± for σ = σ2 and hence is called a generalized cell by Kuske and Matkowsky[14], is not stable for any value of σ. For stable of these critical points, we obtained D˜0

2D˜2 − D˜1

stable structures

>0

≤0

c−1 for σ > σ−1

>0

<0

c−1 for σ−1 < σ < σ2 , b for σ > σ1

=0

<0

c± for σ > 0

=0

>0

b for σ > 0

<0

>0

c+ for σ−1 < σ < σ2 , b for σ > σ1

<0

≤0

c+ for σ > σ−1 25

Therefor, the stabilities of b and c± are determined by the sign of Landau constants D˜0 and 2D˜2 − D˜1 , we obtain in the case D˜0

=

γ(1−2x) 1+x

D˜1 = D˜1 (x, θ) =

3 2

D˜2 = D˜2 (x, θ) =

3 2

− 12 µ0 (γ > 0

is constant ),

−

x(2x−1) f1 (θ) 1−x g1 (θ)

f2 (θ) − (2x − 1) 2g − 2 (θ)

−

x(2x−1) f3 (θ) 1−x g3 (θ)

− (2x − 1) fg44 (θ) (θ) −

2x−1 2β+3 2−2x 1+2β

2x−1 2β+3 1−x 1+2β

, with θ = δkb2 . where f1 (θ)

=

4(1 + x),

g1 (θ)

=

20θ − 3 + (12θ − 5)x,

f2 (θ)

=

(3 − 4θ)g1 (θ) + 2β(1 − 4θ)(5 + 3x),

g2 (θ)

=

g1 (θ)[20θ − 3 + 5 · 2β + (12θ − 5 + 3 · 2β)x],

f3 (θ)

=

3(1 + x),

g3 (θ)

=

12θ − 2 + (6θ − 4)x,

f4 (θ)

=

(3 − 4θ)g3 (θ) + 2β(3 − 2θ)(2 + x),

g4 (θ)

=

g3 (θ)[6θ − 1 + 2 · 2β + (3θ − 2 + 2β)x].

withx = θ (uncoupled) and x = θ + 2β (coupled) The signs of D˜0 and 2D˜2 − D˜1 are shown in Fig.8(a)-(b). 8

6

4

2

0

−2

−4

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Fig.8(a) For different values of θ,the curves D˜0

26

3000

2500

2000

1500

1000

500

0

−500

−1000

−1500

−2000 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Fig.8(b) For different values of θ,the curves D˜2 − D˜1 Proceeding we note that to lowest order the iodide concentration associated with these critical points of b and c is given by  a1 −

1 µ

   b1 − µ    a2 − µ1  b2 − µ





 e11

      d  ∼ A0  11     e21   d21

    cos( 2πx ),  λc  

(51)

and  a1 −

1 µ

   b1 − µ    a2 − µ1  b2 − µ where λc =

2π qc , f (x, y)





 e11

      d  ∼ A±  11 0     e21   d21

πx = cos( 2πx λc ) + 2 cos( λc ) cos(

    f (x, y),   

(52)

√

3πy λc ).

We plot the function f (x, y) in Fig. 8,

where x, y are again being measured in units of λc . From these contour plots we see that (52) gives rise to hexagonal symmetry. Hence, recalling that the Turing patterns under consideration are classified by their light regions, we identify hexagonal arrays of nets or honeycombs with critical point c. Finally, we consider the type of spatial patterns we might expect from the solutions in (51) and (52). The solution of (51) give a striped Turing pattern the same as the description in the previous section and the solution of (52) gives the Turing patterns of hexagonal planform as in [15].

27

−6

−6

x 10

x 10

5

5

0

0

−5

−5 −4

−2

0

2

4

−4

−2

0

2

−6

4 −6

x 10

x 10

−6

−6

x 10

x 10

5

5

0

0

−5

−5 −4

−2

0

2

4

−4

−2

−6

0

2

4 −6

x 10

x 10

Fig.9 The Turing pattern of hexagonal symmetry for critical point (c) 4. Conclusions We have considered two dimensional spatiotemporal structures that can arise in two identical cells, each governed by cubic autocatalator kinetics, which are coupled via diffusive interchange of the autocatalysis B. We have shown that new stable spatially nonuniform structures, both steady (patterns) and oscillatory, arise purely through the coupling of the two cells and are not seen in the uncoupled system. However, these new structures are , in a sense, secondary as the primary bifurcations are the same as in the uncoupled system, with the coupling giving the possibility of further bifurcations. These initially lead to unstable solutions, which can then undergo further bifurcations to sable pattern forms. A major constraint in the present case is that both the coupled and uncoupled systems have the same single homogeneous stationary stage. Thus, bifurcations to spatially nonuniform solutions are from the same uniform state as in the uncoupled system. This is not case when the coupling between the cells is achieved via the reactant A. The case when the coupling is through autocatalyst B is significantly different from that when the coupling is through the reactant A. A necessary condition for the local temporal stability of homogeneous stationary stage is provided by   1, δ ≥ 1, 2 µ2 > µM c (δ, D) =  max{1, µ2 , µ2 }, 0 < δ < 1, 1

28

2

(17)

where µ21 µ22 N

= = =

2 2 kN (1−δkN ) , 2 1+δkN  2 2 k (1−δk  N +1 2 N +1 ) ; 1+δkN +1



0;

Integer part of{ π1

q

N≤

√ √1 π δD

− 1,

otherwise, 1 δλ }.

. On the other hand, We further deduced the possibility of occurrence of two dimensional Turing patterns consisting of rhombic arrays of rectangles and hexagonal by performing the appropriate weakly nonlinear stability analysis of the homogeneous solution and amplitude functions. This is not the case when the dimension of cells 1 and cells 2 is one dimension. REFERENCES [1] R. Hill, J. H.Merkin, Patten formation in a coupled cubic autocatalator system, IMA. J. Appl. Math.,1994 ,53, 265-322. [2] P.Gray and S.K.Scott,Autocatalytic reactions in the isothermal continuous stirred tank reactor:Oscillations and instabilities in the system A+2B → 3B, B → C, Chem.Eng.Sci., 1984,39,10871098. [3] A.Saul and K.Showalter, Oscillations and Travelling Waves in Chemical systems,ed. R.J. Field and M Burger, Wiley-Interscience, New York, 1985. [4] V.V. Gafiychuk and B.Yo. Datsko, Pattern formation in a fractional reaction diffusion system, Physica A 2006,365, 300-306 [5] Q.S. Yang, D.Q. Jiang, N.Z.Shi,C.Y. Ji, The ergodicity and extinction of stochastically perturbed SIR and SEIR epidemic models with saturated incidence, J. Math. Anal.Appl., 2012, 388, 248-271. [6] T.Kolokolnikov, M.J. Ward, J.Wei, Self replication of mesa patterns in a reaction diffusion systems, Physica D, 2007, 236, 104-122. [7] P.J.Atzberger, Spatially adaptive stochastic numerical methods for intrinstic fluctuations in reaction diffusion system, J. Comp. Physics, 2010, 229, 3474-3501 [8] Ji xuehui, Pei Yongzhen, Li Changguo, Two patterns of recruitment in an epidemic model with difference in immunity of individual, Nonlinear Anals. Real World Appli., 2010, 11, 2078-2090. [9] Dezso Hovvath, Agota Toth, Turing patterns in a single-step autocatalytic reaction, J.Chem.Soc.Faraday Trans., 1997, 93,4301-4303. [10] Y. Murroya, H.X.Li, T.Kuniya, Complete global analysis of an SIRS epidemic model with graded cure and incomplete recovery rates, J. Math. Anal.Appl., 2014, 410, 719-732. [11] D.J.Needham and J.H.Merkin, Pattern formation through reaction and diffusion in a simple pooledchemical system , Dynamics and Stability of Systems,1989,4,259-284. [12] Z. Li, L. SY, Bifurcation and Patterns Formation in a coupled higher autocatalator reaction diffusion system, Applied Mathematics and Mechanics, 2007, 28, 1235-1248. [13] L. Zhang, S.Y. Liu, stability and pattern formation in a coupled arbitrary order of autocatalysis system,Applied Mathematical Modeling, 2009, 33, 884-896. [14] Kuske and Matkowsky, On roll, square, and hexagonal cellular flames, European Appl.Math., 1994,5,65-93. [15] David.J.Wollkind and Laura E. Stephenson, Chemical Turing pattern formation analysis:comparison of theory with experiment, SIAM J.Appl.Math., 2000,61,387-431. [16] J.D.Murry,Mathematical Biology, Springer.2002 [17] B.N. Nagorcka, V.S. Manoranjian,J.D.Murry, Complex spatial patterns from tissue interactions:An illustrative model, J.Theoa.Biol., 1987, 128,359-374. [18] G.C.Cruyagen, J.D.Murry, On a tissue interaction model for skin pattern formation, J.Nonlinear Sci. 1992,2,217-240. [19] V.Gaspar,J.Maselko,K.Showalter, Transverse coupling of chemical waves, Chaos,1991,1,435-444.

29