J . Phys. Chem. 1987, 91, 4606-4613
4606
cm3 molecule-’ s-I) is significantly larger than all but two of the previously reported values and is substantially below these two results which are 40% and 150% greater. N o evidence of reaction between SH and O2 was found. An upper bound for the rate constant was inferred which indicates that reaction with O2in the atmosphere is very slow. However, a smaller upper bound would be required to fully preclude a role for this atmospheric reaction. This upper limit determination reflects mainly the detection sensitivity for SH radicals of the
present technique and further investigation of this reaction is needed. Acknowledgment. We thank C. J. Howard for providing us with his results prior to publication. The research described was performed at the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. Registry No. SH, 13940-21-1; NO2, 10102-44-0; 02,7782-44-7.
luring Bifurcation and Stationary Patterns in the Ferroin-Catalyzed Belousov-Z habotinsky Reaction A. B. Rovinsky Institute of Biological Physics of the Academy of Sciences, Puschino, Moscow Region, 142292, USSR (Received: August 18, 1986; In Final Form: February 12, 1987)
A model of the bromate-ferroin-bromomalonicacid chemical system is studied. The model suggests that while the reaction running in a small well-mixed volume has a single, stable stationary state the same reagent spread in a thin layer may spontaneously form stationary periodic spatial patterns of either small or large amplitude, depending on conditions. The homogeneous state corresponding to the stable stationary state of an individual local element of the system may lose its stability because of diffusive coupling between spatially separated points. However, a situation is possible such that both the homogeneous and inhomogeneous(large amplitude) states are stable (bistability, despite the fact that the local system has the only stationary state). In this case, the periodic structure, being in contact with the homogeneous part of the reaction mixture, builds itself up to fill the whole vessel.
Many reaction-diffusion systems exhibit a rich variety of dynamic phenomena. The processes which are observed there can be divided into the two classes: traveling waves and stationary (Standing waves have also been f o ~ n d . ~ . ~ ) The traveling waves are often relaxational. That means that they are of large amplitude and characterized by several different time scales. Propagation of waves of this sort has the very remarkable feature of being in many respects almost independent of the actual mechanism of the underlying chemical or chemicallike process, Le., of the exact form of the kinetic nonlinearities and of the boundary conditions. Such relaxational systems can be described in terms of a rather simple model first suggested by Wiener and Ro~enblueth.~The principal feature of that model medium is that each of its local elements can be excited and immediately after the excitation it cannot be excited again during some refractory time, during which it relaxes to the original state. A most interesting class of standing patterns has been described by Turing.* H e discovered that a homogeneous medium of identical local elements, each being in the same stable steady state, can lose its stability and form a stationary spatialy periodic (1) Nicolis, G.; Prigogine, I. Self Organization in Non-Equilibrium Systems; Wiley-Interscience: New York, 1977. (2) Haken, H. Advanced Synergetics, Springer Series in Synergetics; Springer-Verlag: West Berlin, 1983. (3) Self-Organization. Auto-waues and Structures Far from Equilibrium, Krinsky, V. I., Ed.; Springer Series in Synerdetics; Springer-Verlag: West Berlin, 1984. (4) Non-equilibrium Dynamics in Chemical Systems, Pacault, A., Vidal, C., Eds.; Springer Series in Synergetics; Springer-Verlag: West Berlin, 1984. (5) Boiteux, A,; Ha,B. Ber. Bunsen-Ges.Phys. Chem. 1980,84,392-398. (6) Holmuhamedov, E. L. Eur. J . Biochem. 1986, 158, 543-546. (7) Wiener, N.; Rosenblueth, A. Arch. Inst. Cardiol. Mex. 1946, 16, 205-265. ( 8 ) Turing, A. Phil. Trans. R . SOC.London, Ser. B 1952, 2378, 31-12.
0022-3654/87/2091-4606$01.50/0
structure. The most striking point of the phenomenon is that it occurs despite the fact that the internal kinetic mechanism of the system opposes any deviation from the steady state, and diffusion opposes any deviation from homogeneity. The formation of Turing patterns is strongly dependent on the actual kinetic mechanism of the system and on the boundary conditions. Therefore studies of both waves and stationary patterns are equally important as they can provide a complementary information about the system. The best known example of an excitable reaction-diffusion system is the Belousov-Zhabotinsky reaction catalyzed by the ferroin/ferriin couple. A variety of wave patterns has been found and extensively investigated theree9 However, investigation of stationary structures is far less comprehensive and It is not yet clear whether hydrodynamic and surface effects always are involved in the formation of the patterns, although there is experimental evidence that in some cases they are (Showalter,” Agladze’*). Therefore, theoretical consideration of the phenomenon may provide a deeper insight into its nature and promote further experiments. The first step in this direction was made recently by Becker and Field.I3 They developed a “wind-diagram” method for a qualitative analysis of whether large-amplitude stationary patterns are possible. With this method they found sufficient conditions under which stationary patterns can be expected to develop in the Oregonator model and then demonstrated their existence computationally. Unfortunately, they found the stationary patterns ~
~~
~~~~
~
(9) Oscillations and Traveling Waves in Chemical Systems, Field, R. J., Burger, M., a s . ; Wiley-Interscience: New York, 1985. (10) Zhabotinsky, A. M.; Zaikin, A. N. J . Theor. Biol. 1973,40, 45-60. (1 1) Showalter, K. J. Chem. Phys. 1980, 73. 3735. (12) Agladze, K. I. An Investigation of Rotating Spiral Waues in a Chemical Actiue Medium; Akad. Sci. USSR: Puschino, 1983. (13) Becker, P. K.; Field, R. J. J . Phys. Chem. 1985, 89, 118-128
0 1987 American Chemical Societv
The Journal of Physical Chemistry, Vol. 91, No. 17, 1987 4607
Ferroin-Catalyzed Belousov-Zhabotinsky Reaction SCHEME I
H+
+ HBrO, + HBrO, + HBrOZ++ Br02' + HzO +
F e ( ~ h e n ) ~ ' + HBr02+ * Fe(phen),,+ 2HBr02 + HOBr
+ HBrOZ
+ HBr03
H+ + Br-
+ HBrOZ* 2HOBr H+ + Br- + HOBr Br2 + H 2 0 F=
H+ + Br-
+ HBr0,
H B r 0 2 + HOBr
(1)
(3) (4)
(5) (6)
(13)
(7)
F e ( ~ h e n ) , ~+ + CHBr(COOH)2 == Fe(phen):+ H+ + 'CBr(COOH)2 (8)
+
H20
+ 'CBr(COOH)z + H+ + Br- + 'COH(COOH)z
+ CHBr(COOH)2 CBr2(COOH)2+ HzO + CHBr(COOH)z * CBr2(COOH)z+ H+ + Br-
kIOBr Br,
CHBr(COOH)2
that set "Lo") and by Rovinsky and Z h a b o t i n ~ k y ' ~for, ~40 ~ OC. Field and FOrsterlingz4 elaborated recently a set of the rate constants close to those estimates. Their work also indicates that the protonated radicals HBrOZ+are hardly present (at least in a significant amount) in the system. This fact is very important for the cerium-catalyzed systems, but can be shown to have little importance for the ferroin-catalyzed system at least for sufficiently small catalyst concentrations (since oxybromine radicals are present in the ferroin-catalyzed system in very small concentrat i ~ n s ) . ' ~ Therefore, J~ Scheme I remains untouched here. Directly from Scheme I the following reduced model is derived:'9*20
(9)
Z -dz= x - a dr 1-z
where [ F e ( ~ h e n ) , ~ + ]Z = Cz
(10)
klA X = -x 2k4
[HBrO,] (1 1)
e=-
klA k4C
+ H 2 0+ CHOH(COOH)2 + H+ + Br(12)
only for the diffusion coefficient of the catalyst 10-10000 times larger than those of the other species. This value is rather unrealistic. It has been repeatedly pointed out, however, that the Oregonator, being an excellent model for qualitative studies of the BZ reaction, is not a perfect model for the quantitative description of the s y ~ t e m . ~ Therefore, . ' ~ ~ ~ Turing bifurcation and stationary pattern formation are studied here with a model more relevant for the quantitative description of the bromate-ferroin-bromomalonic acid s y ~ t e m . ' ~ J A ~ -more ~ ~ direct and informative approach is applied and more details are obtained. The domains of Turing instability are found in the space of experimental parameters. It is shown that stationary patterns can be formed through both a soft-mode instability of a uniform steady state and a hard-mode excitation. The latter means that the spatial system can exhibit bistability. In all cases the corresponding local system has a single, stable steady state, and the diffusion coefficients of all species are of the same order.
The Model It has been shown earlier'5*'8-20~22 that the ferroin-catalyzed Belousov-Zhabotinsky reaction can be quantitatively described by the modified FKN mechanism shown in Scheme I, with the set of rate constants first introduced by Rovinsky and Zhabot i n ~ k y ' and ~ , ~later ~ reestimated by TysonZ3for 25 O C (he called (14) Rovinsky, A. B.; Zhabotinsky, A. M. Teor: Eksp. Khim. 1978, 14, 182-1 9 1. (15) Rovinsky, A. B.; Zhabotinsky, A. M. Teor. Exsp. Khim. 1979, 15, 25-32. (16) Noszticzius, 2.;Farkas, H.;Shelly, Z. J . Chem. Phys. 1984, 80, 6062-6070. (17) Kuhnert, L.; Krug, H.-J.; Pohlmann, L. J . Phys. Chem. 1985, 89, 2022-2026. (18) Rovinsky, A. B. J. Phys. Chem. 1984, 88, 4-5. (19) Rovinsky, A. B.; Zhabotinsky, A. M. J . Phys. Chem. 1984, 88, 608 1-6084. (20) Rovinsky, A. B.; Zhabotinsky, A. B. In Fundamental Research in Chemical Kinetics, Shilov, A. E., Ed.; Gordon and Breach: London, 1986. (21) Rovinsky, A. B. J . Phys. Chem. 1986.90, 217-219. (22) Field, R. J.; K6ros, E.; Noyes, R. M. J . Am. Chem. SOC.1972, 94, 8649-8665. (23) Tyson, J. In Oscillations and Traueling Waves in Chemical Systems,
C = [ F e ( ~ h e n ) , ~ ++] [Fe(phen),,+] A = [HBrO,]
B = [CHBr(COOH),]
ho is the acidity function, q is the stoichiometric factor (q = Brproduction r a t e / F e ( ~ h e n ) ~reduction ~+ rate), and k*i are the rate constants of steps 1-12. The values of the rate constants are as folloWs:'9~~0 k l = 100 M-2 s-I k4 = 1.7 X lo4 M-' s-' k5 =
lo7 M-2 s-'
K* = 2 x 10-5 M s-1 kI2=
s-l
q = 0.5
Model 13 describes the local kinetics of the system. Consideration of diffusion yields
-dz=
dr
x- a -
Z
1 -z
+ 6A,z
Field, R. J., Burger, M., Eds.;Wiley-Interscience: New York, 1985; Chapter 3.
(24) Field, R. J.; F6rsterling, H.-D. J. Phys. Chem. 1986, 90, 5400-5407.
4608
Rovinsky
The Journal of Physical Chemistry, Vol. 91, No. 17, 1987
A
C
Figure 2. Null clines of system 13.
I Figure 1. Representation of the Turing bifurcation on the A-k plane.
where
--)
DZ
DX I DHBr02
DFe(phen)33+
6 = - DZ
p i = ri( kI2A2ho1
DX
k4C Dx
We assume hereafter Dx = below.
’I2
cm2/s. This model is studied
small-amplitude spatially periodic stationary solution, a largeamplitude periodic stationary solution, or a nonstationary solution with periodic or more complex temporal behavior is possible. Some examples of that sort will be given below. It should be stressed, however, that the loss of stability of a homogeneous steady state through the Turing bifurcation is not always a prerequisite for an onset of those patterns. A situation is possible under which both a homogeneous steady state and a periodic pattern of a large amplitude exist under the same conditions. Consider now more closely system 17 for two species and suppose
General Sketch of the Turing Bifurcation
If a local system li 24
=
(UI,
..., u,)
= f(u)
(15)
Let A = aij = implies
Stability of the corresponding local system
f = VI, ...,f,)
all
has a steady state
+ a22
0 u = uo
then the corresponding reaction-diffusion system u = f(u) + DAu (D is a diffusion matrix)
(16)
The eigenvalues of the L matrix are XI,2(k) =
has the same homogeneous steady state. The stability of system 15 is determined by the eigenvalues A,, ..., A, of the matrix A = df/du evaluated at the stationary p i n t . The state is stable if each of the eigenvalues has a negative real part. Stability of the uniform solution u ( x ) = uo of eq 17 is given by the eigenvalues of the matrix L = af/au - KD (18)
4A)’I2)/2
(19)
in a finite region. The solution u = uo is stable with respect to small perturbations if all eigenvalues of the L matrix have strictly negative real parts for all possible k . Note that the case k = 0 corresponds to the solution of local system 15. Let the function f be dependent on some parameter p which in our case is usually a concentration, diffusion coefficient, or reactor length. Figure 1 illustrates a possible situation when p changes. For p > pcrltthe solution u(x) = uo of eq 17 becomes unstable to the perturbation w oexp(ikx) with k close to ko. This kind of bifurcation was found by Turing.* H e showed that such an instability may result in a spatially periodic stationary pattern if Ab) is real or give rise to traveling waves (and, presumably, linear standing waves) if h(p) is complex. The latter is possible only in systems with no less than three essential components (variables). In this case Figure 1 should be considered as a plot of the real part of the eigenvalue. The linear theory of Turing bifurcation in multicomponent systems has been significantly advanced by Othmer and Scriven.2s The exact character of an appearing solution is determined by the nonlinear terms of the equation at the bifurcation point and by the boundary conditions. In a two-component system a
+ a22 - (Dl + D2)kZ A = A0 - (alID2 + a22D1)k2 + D,D2k4 u
= all
(22a) (22b)
is positive and the solution u = uo If A < 0, then one of the is unstable. That occurs when k lies in the interval
XI < k2 C X2
(23)
where x1.2 =
kiiD2
+ az2D1 f [ ( a 1 1 4+ a22D1)’ - 4 D i D ~ l ” ~ 2DlD2
Condition 23 can be obtained only if
a11D2 + a22D, > 0
(24)
Taking into account eq 20a, one concludes that a necessary condition for instability is that a l l and a22 must be of opposite sign. Therefore, if a , , > 0, then D2 > D , or if azz > 0 then Di > D2.26 Model Calculations
The matrix A for eq 14 is as follows: A = lajjl
~~
(25) Othmer, H.; Scriven, L. J . Theor. Biol. 1971, 32, 507-537.
(21)
where
where K = k2 is any positive number for an unlimited system or any eigenvalue of the problem
Au = K U
(-u f (a2 -
(17)
(26) Segel, L.; Jackson, J. L. J . Theor. Eiol. 1972, 37, 545-559
TI. re Journal of Physical Chemistry, Vol. 91, No. 17, 1987 4609
Ferroin-Catalyzed Belousov-Zhabotinsky Reaction
- I;
I
do2
do4
/
0.06
,
Figure 3. Regions of homogeneous (linearly stable) stationary states (shaded), bulk oscillations (bounded by for C = 3 X lo”, - - for C = and for C = 3 X and the Turing instability. The instability is possible on the part of the plane below the shaded region and beyond the bulk oscillations region. ho = 1 .
-
-.-.
Here zo and xo are the coordinates of the rest point of local system 13, which is the crossing point of the null clines x = 0 and i = 0 (Figure 2). Since u22is always negative, the Turing instability of the homogeneous steady state can occur only if a l l > 0, and consequently only if Dz > Dx It is easy to see that if a rest point of system 13 lies on the branch AB of the null cline x = 0 (then the point corresponds to a high ferroin concentration, “red state”) or branch.CD (a high ferriin concentration, “blue state”), then a l l < 0 and the homogeneous steady state of eq 14 is also always stable. If, however, the stationary point lies on the branch BC, then a,, > 0. There are two regions on this branch where all a22< 0 one is adjacent to point B, and the other is adjacent to point C. (However, a situation is possible when all a22< 0 on the whole branch BC.) If the stationary point is in one of these intervals the stationary state of system 13 is still stable, while that of system 14 can lose its stability at appropriate values of Dz/Dx Figure 3 shows the domains on the [bromate]-[bromomalonic acid] plane where the Turing bifurcation is potentially possible. These domains occupy the region of the plane below the line all(A,B) = 0 and beyond the oscillation region of system 13 (limited by the line a,,(A,B) a2,(A,B) = 0). The oscillatory region shrinks rapidly with decreasing total catalyst concentration, C, and the domain where the Turing instability can occur expands at its expense. Figure 4 shows the parts of the parameter plane where the homogeneous solution of eq 14 is unstable. U p to now the linear problem has been considered, and the results obtained are independent of the boundary conditions and even of the number of spatial dimensions. In order to get information about the state to which the system evolves after the loss of homogeneity, one should take into account these aspects as well as the nonlinearities. For this reason, further considerations will be confined to a one-dimensional vessel of finite size with impenetrable walls (Neumann no-flux boundary conditions). At this stage of our study, we can seek an answer to the question of whether a stable small-amplitude stationary (inhomogeneous) solution exists near the bifurcation point. This goal can be achieved by implementing a projection technique onto a so-called central manifold of eq 17.27*28This method makes use of the fact
+
+
+
(27) loss, G.; Joseph, D. Elementary Stability and Bifurcation Theory; Springer-Verlag: West Berlin, 1980.
I
I
1
a02
004
a06
I
ma
Figure 4. Regions where the homogeneous state is unstable. These regions lay below the shaded zone (absolutely stable steady states) and beyond the bulk oscillation domain (solid line) and are limited by the lines: for Dz/Dx = 1.5, - - - for Dz/Dx = 2, and for D,/Dx = 4.
-.-.
-.e---
that the behavior of a system near a bifurcation point is completely determined by evolution of some “normaln variables. These variables correspond to the eigenvalues of the L matrix crossing the imaginary axis. Therefore, the multicomponent partial differential equation (eq 17) is reduced to a single-component or two-component ordinary differential equation if, respectively, a real eigenvalue or complex conjugate pair crosses the imaginary axis. In our case, the answer is positive (that is, a small-amplitude periodic stationary solution is stable) if a coefficient of the cubic term of that equation (equivalent to Lyapunov’s G3 number) is negative. Appendix I describes how this number is calculated. Computations show that, in the region explored, G, is always negative and consequently the stable small-amplitude Turing structures exist. Further study involves direct numerical analysis of eq 14 with a finite-difference technique. The number of mesh points of the space interval was varied from 12 to 82. Most of the calculations were, however, checked with a more dense discretization, and that revealed that a further increase of the number of points does not change the results significantly. In the case considered, the stability of the computations was the most sensitive to the diffusive terms. Therefore the Crank-NicoLon implicit method29appeared to be the most efficient. Stationary solutions of eq 14 were also sought by using a Newtonian iteration technique. In that case, the space interval corresponding to a cycle of the periodic structure was divided into 12 subintervals. The problem thus was reduced to solving 24 algebraic nonlinear equations which give the stationary solution of the ODE system dr
t
(xj-1
dzi dr
- = xi - -a
Zi
1
- zj
+ xi+l - 2 x j ) / h 2 (14)
6 + -(zi-, + h2
Zj+l
- 2Zi)
where h is the space interval between two neighboring mesh points, (28) Hassard, B. D.; Kazarinoff, N. D.; Wan, Y.-H. Theory and Applications of Hopf Bifurcation; Cambridge University: London, 198 1. (29) McCracken, D. D.; Dorn, W. S. Numerical Methods and Fortran Programming, Wiley-Interscience: New York, 1965.
4610
The Journal of Physical Chemistry, Vol. 91, No. 17, 1987
Rovinsky
11 0
ob. O 0
L
.
I l3l
'
.ol.i
*
. .olg
. .0!6
m
Figure 5. A small-amplitudestationary pattern for A = 0.044, E = 0.57, c = 10-5, ho = 1, D,/D, = 1.5.
Figure 6. A large-amplitudestationary pattern for A = 0.02, E = 0.192, C= h, = 1, D,/D, = 25.
..., 12 and x1 = x2, z1 = z2, x l l = x12,z l l = z I 2due to
phenomena accompanying such a bifurcation will be described elsewhere.32
i = 1,
boundary conditions. Stability of the solutions was determined by calculating the eigenvalues of the Jacobian of the system using the QR method.30 Figures 5 and 6 show the stationary patterns evolving from the unstable homogeneous steady state. In these computations the homogeneous steady state was initially disturbed with a cosine function of amplitude about 5% of the steady-state values. In some cases, the amplitude of the resulting structure was unexpectedly large (Figure 6), despite the fact that near a bifurcation point a stable stationary solution of a small amplitude should exist. The situation becomes clearer looking at steady-state stability as a function of the ratio Dz/D, On these plots the maximum value of the x variable is taken as an ordinate (Figures 7 and 8). Figure 8 suggests a probable explanation. The small-amplitude solution is hardly detectable because it is either too close to the unstable branch or stable only over a very short interval of the parameter. Figure 8 demonstrates another interesting result. There is a parameter interval where both homogeneous and inhomogeneous stable stationary states exist. This conjecture has also been supported by direct computations using eq 14 (Figure 9). In this case, if a system of sufficient length rests in the homogeneous steady state, and a single pulse of a periodic pattern is placed at the wall of the container, the pulse reproduces itself until the whole vessel is filled with the periodic structure (Figure 10). A similar result was also reported by Becker and Field.I3 It is remarkable that exactly the same phenomenon was observed by Zhabotinsky and Zaikin in their experiments" (the corresponding illustrations are, however, presented in ref 31). Figure 8 also shows that a large-amplitude stationary structure may lose its stability through an oscillatory mode. Interesting (30) Wilkinson, J. H. The Algebraic Eigenualue Problem; Clarendon Press: Oxford, 1965. (31) Kawczinsky, A. L.; Zaikin, A. N. J . Non-Equilib. Thermodyn. 1977,
2. 139.
Discussion The principal results of the study can be summarized in the following statements. (a) A spatially homogeneous state of the Belousov-Zhabotinsky medium can become unstable at reasonably low ratios of diffusion coefficients of ferriin and bromous acid (slightly more than unity). However, the greater the ratio, the wider is the parameter region of Turing instability. This region also expands when the catalyst concentration falls. (b) A Turing bifurcation is more probable in the part of parameter space close to the "blue" boundary of the oscillatory domain of the local system (such a state, however, will be hard to maintain experimentally [reviewer's remark, Agladze (personal communication)]). (c) A large-amplitude stationary structure can arise as a result of the instability of the homogeneous state. The stability regions of the large-amplitude structure and of the homogeneous state often overlap (bistability). So, depending on initial conditions, the system can be found in either the homogeneous or inhomogeneous steady state for the same set of parameters. It can be transferred from one to the other by an appropriate perturbation. (d) Being in the bistability region, the system can "build up" the periodic elements at the expense of the homogeneous state. This phenomenon has been observed experimentally" and no explanation has been presented for it in terms of hydrodynamic flows or gas evolution. Therefore, this calculation seems to be strong evidence that reaction-diffusion stationary structures really exist in the Belousov-Zhabotinsky chemical medium. (e) Although diffusion-induced instability of the homogeneous state can occur only if the diffusion coefficients of the ferriin ion exceeds that of bromous acid, a large-amplitude stationary pattern ~~~
~
(32) Rovinsky, A B J Phys Chem , accepted for publication
The Journal of Physical Chemistry, Vol. 91, No. 1 7 , 1987 4611
Ferroin-Catalyzed Belousov-Zhabotinsky Reaction
'T s
I
I
I
1.4
is
4.45,
Figure 7. Steady-state bifurcation diagram for the parameter values corresponding to Figure 5 . Solid lines correspond to stable solutions; dashed lines, unstable solutions. The horizontal branch represents the homogeneous steady state.
L
I
I
I 0
-i I
tm
2 1
Figure 8. Steady-state bifurcation diagram for the parameters values corresponding to Figure 6. Solid lines represent stable solutions while the dashed lines represent unstable ones. The horizontal branch corresponds to the homogeneous stationary state. Inhomogeneous stationary solution (upper branch) loses its stability through Hopf bifurcation.
can exist even if DFe(phcn),3+ < DHBa2(which is more likely in the real system) (Figure 11). It was stated above that spontaneous symmetry breaking of a uniform steady state can only occur if diffusion coefficients of different species differ sufficiently. However, if the symmetry of the system is already broken with inhomogeneous initial conditions of appropriate value the resulting structure may be of any sort for any ratio of diffusion coefficients. In the case illustrated in Figure 11, this inhomogeneous structure C DHBa2. However, there is no reason is unstable for DFe(phcn)33+ to believe that a stable inhomogeneous solution is impossible if the ferriin diffusibility is less than that of bromous acid. More computations are needed to get a precise answer. As the Oregonator model used by Becker and Field13 is a rather general model it does not differentiate between the cerium-catalyzed system and ferroin-catalyzed system. The model used here is more specific and the results are applicable only to the systems with low-redox-potential catalyst such as F e ( ~ h e n or ) ~Fe(dipy),. Besides, the classical Oregonator does not involve the total catalyst concentration while the results presented here are strongly dependent on this parameter. These distinctions are important in view of possible experimental realization of the regimes. There is a problem concerning justification of low-order models like model 15. The-reduction of a high-order system of ordinary differential equations to a low-order set can be performed on the basis of singular-perturbation theory. Apparently, there is no such strong foundation for the reduction of high-order systems of PDE's. This fact does not constitute a difficulty for many questions concerning relaxation waves. It has already been pointed out that often the main properties of the relaxation autowaves are almost completely determined by the properties of the local elements of the medium. That means that the reduction can be made for the
Figure 9. Stationary pattern in the bistability region. The parameters are the same as for Figure 6 except for D z / D x = 3. The homogeneous steady state is also stable for this set of parameters.
local system (ODES). This is, however, not necessarily the case for diffusion-instability problems. Then the question arises of what are the conditions under which the same Turing bifurcation takes place in both the high-order system of P D E s and the reduced system. An attempt was made to consider this problem with respect to the third-order model preceding system 14. Appendix I1 shows that, for the homogeneous steady state of this three-component model lying on the branches AB or CD of Figure 2, a Turing instability is impossible just as for the second-order system 14. (The steady States of both systems are the same.) Computations (though not comprehensive) show that for the other homogeneous stationary states of this model, the Turing bifurcation can occur at approximately the same set of the parameters as in model 14. These results make us more-or-less sure that the consideration of the low-order is justified in this case. A firm mathematical foundation is still needed. however. Acknowledgment. The author is grateful to Prof. A. M. Zhabotinsky and Dr. Yu.A. Kuznetzov for useful discussions, and to Dr. I. A. Goryanin for permission to use and modify his Fortran program package for solving nonlinear algebraic equations and continuing the solutions.
Appendix I A restriction to the center manifold and the theory of stationary bifurcation is described in ref 28. Here this approach is loosely sketched and the algorithm for calculation of the G3number is given (the G3sign determines whether a small-amplitude stationary solution is stable near a Turing bifurcation point). System 14 is considered on a finite one-dimentional segment with Neumann no-flux boundary conditions. First consider a system of ODE'S x =f(x,p) (AI.1)
The Journal of Physical Chemistry, Vol. 91, No. 17, 1987
4612
Rovins kq
A
,
-
.
,
i
,
,
,
,
1
.
u
i
rm
d
tu-l
0
e
I
1
f
Figure 10. The reactor, which was initially in the homogeneous steady state, is filled with a periodic strucutre The time interval between part a-e is 20 s, f is the final pattern The parameters values are the same as for Figure 9 I
where p is a 1 vector directed along the eigenvector q corresponding to the zero eigenvalue A,@,) of the A matrix and w is an (n - 1) vector in the complementary subspace. Let A+ be the adjoint matrix to A (in the considered case A + is the transpose of A ) and q+ is the eigenvector of A+ corresponding to its zero eigenvalue (that is A+q+ = 0). Then the scalar product (q+,w) = 0. In the new variables, eq AI.2 takes the form
+ g(u,w) @ = B@)w + h(u,w) u = A,(p)u
1
I
'0
f 09 (DFeiphe"g/D"POg)
(AIS)
where the B matrix is real and has eigenvalues A j @ ) , j = 2, ..., n (A,@), A,@), ..., A,@) are the eigenvalues of A ) . Functions g and h are
Figure 11. Bifurcation diagram for the same set of the parameters as for Figure 6 with the exception that C = 1.8 X
where x is an n vector and p is a parameter. Let x0@) be a locus of stationary points of (AI.l) and p = p c be a critical value of ..., A,@,) the parameter at which the eigenvalues A,@,), A,@,), of the Jacobian matrix A = tlf/tlx(xo(p),p)obey: AI@,) = 0 and Re A,@,) < 0 for j = 2, ..., n. Thus, pc is a point of the bifurcation of the xo solution. The change of variables x(P) =
XOh)
+ Y@)
Expanding g and h in Taylor series 1 g(u,w) = -g,u2 2 1 h(u,w) = -hzu2 2
(AI.2)
1 + ugtw + -g3u3 + 0(luI4) 3!
(AI.7)
1 + uhlw + -ha& +O(~U~~,IW~~) 3!
+
yields j , = AY
+ FQ,p)
(AI.3)
where FQ,p) = f Q , p ) - Ay is at least a second-order function of y , and A is the Jacobian matrix. Make one more change of the variables y = u + w (AI.4)
Near the bifurcation point, p = pc p, A,@) = a p + O(p2),and p is a small parameter. Since all Re A, < 0 and IRe Ajl >> p (j = 2, ..., n), we can apply the methods of singular perturbation theory and exclude w from (AIS) with accuracy of the order of P: *
B' w = --h2ti2
2
(AI.8)
The Journal of Physical Chemistry, Vol. 91, No. 17, 1987 4613
Ferroin-Catalyzed Belousov-Zhabotinsky Reaction
predecessor of model 14. Put it in the unscaled form:
So (AI.5) becomes
X = klh&
Let a > 0 (Le. x ( p ) is stable to the left of pc and unstable to the right). Then the stationary solutions of (AI.1) not equal to x 0 b ) but close to it in the neighborhood of p c are determined by the stationary nonzero solutions of eq A1.9. In many cases (usually when there is some symmetry), g2 = 0. Then for p > 0 (p > p,) the small-amplitude stationary solution of (AI.9) exists until 1 G3 = g3 - j B ' g i h 2 < 0 (AI. 10) Formulas AI.2, AI.6, AI.7, and AI.10 provide the way for calculation of G3. If one considers a bifurcation of a homogeneous solution of the system. x = f ( x , p ) + DAx (AI.11)
'
- k4h#
KaBZ = 'h0(C - Z )
+ k12B - k 7 h d Y -
+ Aox z = g(x,z) + 6$z
:(
x(l
- x) - ( 2 q a1Z--z
g(x,z) = x
- a-
z=
2 kl ho2AXoC
K8B + 2klho2AXo
(AII.4)
The matrix of the linearized system AIL1 is shown in eq A I M .
(AII.5)
It is easy to see that for the steady state given by AII.2-AIL4 the only positive elements of the matrix are a23 and ~ 3 1 .System A.1 will undergo a Turing bifurcation if some real eigenvalue of the matrix
where f(x,z) =
(AII.3)
(15)
x =f(x,z)
+ DyAY
where X = [HBr02], Y = [Br], Z = [ F e ( ~ h e n ) ~ ~while + ] , the rest of the designations are the same as for system 14. At any stationary point X > k 7 A / k S . If the stationary point lies on branch AB (Figure 2), then approximately the stationary concentrations are
+
and q should be considered as an eigenfunction of the L operator. In this case the scalar product involves integration over the spatial region. If system AI.11 is considered on a one-dimensional segment with Neumann boundary conditions, then the eigenfunctions of the problem are cosine functions. It is not difficult to show that in this case both g2 and h2 are zeros. Then G3 = g3. Return now to system 15. Put it in the form
kshoXy
(AII.1)
( D is a diffusion matrix) in a bounded region, the same consideration is valid after the change A L = af/ax DA
-
- kSh&Y + k7hoAY + DXAX
+ 8)=) X+h
Z
1- z
Let
L=A-xD;
D=[:
Dx 0
(AII.6)
becomes zero at any positive x . In other words this happens when at a positive x P ( x ) = det(L) = 0 (AII.7)
In more detail P ( x ) = -60
- 6 , x - 62x2 - 6 3 x 3
where ko2 is such that det A = 0. Then up to a positive factor g3 is g3 = 428zzz923+ ql(fzzz923+ 3 f , , ~ 2 ~ 4+13f,,q12q2 + f x x ~ 1 3 ) where q=
[3
q+ =
[E:]
Appendix I1 Here we consider a model which is the three-component
We are considering the stationary state which is stable for the local system. So, according to the Hurwitz stability criterion a0 = det ( A ) < 0 It can be shown that until the stationary point remains on the branch AB a,, < 0 and a11a22 - a12a2,> 0 . Then the coefficients of the polynomial P ( x ) are all positive. Therefore it has no positive root. That means that the Turing bifurcation is impossible. The proof for the C D branch is similar. Registry No. BrOp-, 15541-45-4; bromomalonic acid, 600-31-7; ferroin, 14708-99-7.
