J. Phys. Chem. 1985, 89, 118-128
118
obtains for the total hologram efficiency in the plane wave case qtotal = sin2 [ F G t / 2 ]
(A-12)
For Gaussian writing and reading beams This case yields qtotai
=
[
1-
7 1 sin (FGt)
(A-14)
For the limit of no depletion of the reacting material, it can be seen that the plane wave hologram achieves a maximum efficiency of 100% but the Gaussian beam hologram reaches a maximum of only 61%. This limit has been previously d i s c ~ s s e d . ' ~ The second case to consider is the general case where the photochemically active material may be used up. For the plane wave the integrations of eq A-8 and A-9 yield qtotal = sin2 [Fexp(-Gt)ll(Gt)]
(A-15)
which reaches a maximum in a time given by eq 14. The corresponding integral for the general Gaussian case must be evaluated numerically. To compare the values of t,, obtained for plane and Gaussian beams, it was assumed that the two beams had equal total power and that the Gaussian width is half of the plane wave width. This is equivalent to assuming that the peak Gaussian intensity was twice the average plane wave intensity. In this case the ratios of the time to reach maximum efficiency for the two beams was found to be
irrespective of the intensity of the writing beams (provided of course that the power densities were equal). This factor can be used to correct for the Gaussian nature of the interfering beams. For more complex beam profiles the correction factor will be different of course. The factor can, however, be obtained from the beam shape by numerical integrations of eq A-9. Registry No. Biacetyl, 43 1-03-8.
Statlonary Concentratlon Patterns in the Oregonator Model of the Belousov-Zhabotinskii Reaction Paul K. Beckert and Richard J. Field* Departments of Physics and Chemistry, University of Montana, Missoula, Montana 5981 2 (Received: July 12, 1984)
Stationary mosaic patterns in the concentrations of intermediateand catalyst species have been observed in various modifications of the oscillatory Belousov-Zhabotinskii (BZ) reaction when the chemical reagent is spread in a thin layer on a flat surface. They are somewhat similar in appearance to the familiar Zaikin-Zhabotinskii-Winfree traveling concentration waves, except that they are stationary. Concentration inhomogeneitiesin the observed stationary patterns repeat themselves over a length of 2 to 3 mm. The mechanism of this stationary pattern formation has not been clear. Interpretation as resulting only from the interaction of reaction and diffusion has been complicated by the obvious presence of convective and heterogeneous (e+, gas transfer across the surface of the reagent) effects. Furthermore, there has been only fragmentary theoretical evidence for reaction-diffusion supported stationary patterns in the Oregonator model of the BZ reaction. We now have found strong theoretical support for the formation of reaction-diffusion supported patterns in the Oregonator. Solutions have been obtained for the Oregonator partial differential equations that correspond to both small- and large-amplitude stationary concentration patterns. The small-amplitude patterns are calculated by a perturbation technique and have the profile of a cosine function. The large-amplitude patterns are calculated by numerical integration and have a spiked profile related to the bistability of the BZ reaction. Both types of patterns are boundary condition dependent. Their stability is investigated by numerical integration. The small-amplitudepatterns may be long-lived if the size of their container is very close to a multiple of a critical wavelength cosine function. However, they are ultimately unstable, decomposing usually to an oscillatory but spatially homogeneous state. In containers of size close to a critical length any small perturbation to an initially stationary, spatially homogeneous steady state tends to evolve toward a unique, quasi-stable, small-amplitude pattern. In one case, such a small-amplitude pattern evolves to a stable large-amplitudepattern. The largeamplitude patterns are shown to be stable by numerical integration. Multiple-peak patterns appear with concentration inhomogeneities large enough to be seen experimentally. A mathematical analysis is developed which allows the prediction of conditionsunder which large-amplitude patterns can appear. The calculated length scale of the large-amplitude patterns is very similar to that observed experimentally.
Some remarkable chemical reactions governed by nonlinear kinetic laws exhibit when well stirred bulk temporal oscillations1 in the concentrations of intermediate and catalyst species. The best known and understood of these reactions is the BelousovZhabotinskii (BZ) reaction, a metal-ion-catalyzed (Ce(IV)/Ce(111) or ferriin/ferroin are often used) oxidation and bromination of malonic acid in 1 M sulfuric a ~ i d . ~ Oscillations J are easily measured in [Br-] and the ratio [oxidized form]/[reduced form] of the metal-ion catalyst. The reactions leading to temporal oscillations in the BZ reaction may couple under unstirred conditions with diffusion of reactive
-
Department of Physics. Current Address: School of Oceanography, University of Washington, Seattle WA 98195. *Address correspondence to this author at the Department of Chemistry.
0022-3654/85/2089-0118$01.50/0
intermediates, mainly HBr02, Br-, and the oxidized form of the metal-ion catalyst, to yield traveling concentration inhomogeneitiesk6 with [Br-] much lower and [HBr02], as well as the ratio [oxidized form]/[reduced form] of the metal-ion catalyst, much (1) Field, R. J., Burger, M., Eds.; 'Oscillations and Traveling Waves in Chemical Systems"; Wiley-Interscience: New York, 1984. (2) Field, R. J.; Koros, E.; Noyes, R. M. J . Am. Chem. SOC.1972, 94, 8649. (3) Field, R. J. In 'Oscillations and Traveling Waves in Chemical Systems"; Field, R. J., Burger, M., Eds.;Wiley-Interscience: New York, 1984. (4) Field, R. J.; Noyes, R. M. J . Am. Chem. SOC.1974, 96, 2001. (5) Showalter, K.;Noyes, R. M.; Turner, H. J . Am. Chem. Soc. 1979,101, 7463. (6) Winfree, A. T. In 'Oscillations and Travrstng Waves in Chemical Systems":Field, R. J., Burger, M., Eds.;Wiley-InWence: New York, 1984.
0 1985 American Chemical Society
Concentration Patterns in the Oregonator Mode higher than in the bulk of the reagent. These traveling inhomogeneities can be qualitatively and to some extent quantitatively ~ n d e r s t o o d ~within - ’ ~ the framework of a simplified model”’ of the BZ reaction kinetics. It is found that the BZ traveling inhomogeneities are very similar in mathematical character to nerve impulses.12 Traveling BZ concentration inhomogeneities do not depend intimately on the size or shape of their container. They travel freely unless near a wall. However, there are several theoretical predictions that a chemical reagent with properties like the BZ reaction might also support stationary boundary-condition-dependent concentration patterns. Use of the word pattern here will imply a stationary spatial concentration inhomogeneity. Pioneering work on pattern formation in chemical systems was done by Turing13 working with a set of partial differential equations (PDEs) describing the reaction-diffusion properties of a hypothetical set of chemical reactions. He found solutions corresponding to periodic patterns of intermediate concentrations and suggested the phenomenon as a model of spontaneous pattern formation in biological systems. There is still considerable interest14 among biologists in this kind of model. Turing’s ideas were extended by Prigogine and c o - ~ o r k e r susing ’ ~ a derivative of Turing’s model called the Brusselator.16 Similar theoretical work using other models has been reported by Zhabotinskii and Zaikin” and by Kawcyznski and Zaikin.’* Experimental observation of pattern formation in the BZ reaction and related systems has been reported by Zhabotinskii and Zaikin,” S h ~ w a l t e r , and ’ ~ O r b h 2 0 Interpretation of these experiments is unfortunately clouded by the presence of convective and surface effects, and it is not clear that they are verification of patterns maintained only by the interaction of reaction and diffusion. The situation has been further complicated by the fact that there has been little theoretical basis in the Oregonator model of the BZ reaction to suggest stable reaction-diffusion patterns. We have now found such evidence in the form of solutions to the Oregonator PDEs showing both small- and large-amplitude patterns and whose stability can be tested by numerical integration. Small-amplitude patterns may be very long lived but are ultimately unstable. Large-amplitude patterns are stable. All of our analysis is done for one spatial dimension.
The Oregonator The Oregonator model and its relationship to the chemistry of the BZ reaction has been described in detail It reduces a very complex mechanism2g3of the BZ reaction to five steps involving only the three intermediates HBr02, Br-, and the oxidized form of the metal-ion catalyst. It has been successful in nearly quantitatively describing all of the experimentally observed behaviors of the BZ although some questions concerning its applicability to certain have been raised recently21-22
(7) Reusser, E.; Field, R. J. J. Am. Chem. SOC.1979, 101, 1063. (8) Field, R. J.; Troy, W. C. S I A M J . Appl. Math. 1979, 37, 561. (9) Field, R. J.; Noyes, R. M. J . Chem. Phys. 1974, 60, 1877. (10) Tyson, J. J. In “Oscillations and Traveling Waves in Chemical Systems”; Field, R. J., Burger, M., Eds.;Wiley-Interscience: New York, 1984. (11) Troy, W. C. In “Oscillations and Traveling Waves in Chemical Systems”; Field, R. J., Burger, M.,Eds.; Wiley-Interscience: New York, 1984. (12) Rinzel, J.; Miller, R. N. Marh. Biosci. 1980, 49, 27. (13) Turing, A. M. Phil. Trans. R. SOC.London, Ser. B 1952, 237, 37. (14) Winfree, A. T. “The Geometry of Biological Time”; Springer-Verlag: New York, 1980. (1 5 ) Nicolis, G.; Prigogine, I. ‘Self-organization in Nonequilibrium Systems”; Wiley-Interscience: New York, 1977. (16) Tyson, J. J. J . Chem. Phys. 1973, 58, 3919. (17) Zhabotinskii, A. M.; Zaikin, A. N. J . Theor. Biol. 1973, 40, 45. (18) Kawcyznski, A. L.; Zaikin, A. N. J. Non-Equilib. Thermodn. 1977, 2, 139. (19) Showalter, K. J. Chem. Phys. 1980, 73, 3735. (20) Orbln, M. J . A m . Chem. SOC.1980, 102, 4311. (21) Noszticzius, Z.; Farkas, H.; Schelly, Z. A. J . Chem. Phys. 1984.80, 6062. (22) Noyes, R. M. J . Chem. Phys. 1984,80, 6071.
The Journal of Physical Chemistry, Vol. 89, No. 1. 1985 119 modifications of the BZ reaction. The timespace behavior of the Oregonator is described in one spatial dimension by the scaled PDEs as follows: (ax/dt), = D’,.(d2~/dl2),
+ S(Y
- XY
+ x - qX2)
(az/at), = DL(a2z/d12), + W ( X - Z)
(1.1)
(1.3)
In eq 1.1-1.3 x = k l X / ( k 3 A ) ,y = k2Y/(k3A), z = k2k5Z/ ( k l k 3 A 2 )t, = ( k l k 3 ) 1 / 2 A qT ,= 2klk4/(k2k3),s = (k3/kl)1/2,w = ks/[(klk3)’/2A], Di = Di/[(klk3)1/2A],where 1 is the distance (mm), Tis the time (s), X , Y, Z , and A are, respectively, the molar concentrations of HBr02, Br-, oxidized form of the catalyst, and BrO), kl-ks are the mass-action rate parameters of the five steps in the model, f is a stoichiometric factor in step 5, and Di is the diffusion coefficient (mm2/s) of the ith species. All parameters are held constant here except D, and f which are varied in the ranges 0.5 -
a
s
Setting dx/dt = dy/dt = 0 in eq 2.1 and 2.2 yields the lines of zero derivative (zero rate of chemical change in time, or zero “chemical velocity”) given by
9
0
dx - 0 - y _ dt
=
x - 4x2 x-1
-
(3.1)
dv
co Q! -
N
a
0
-I
;
;
0 .o
4.0
8.0
DISTANCE IN MM Figure 1. Large-amplitude pattern f o r f = 2.5,Dz = 10 mmz/s, and L = 10 mm. This pattern grew under numerical integration from the initial waveform shown in Figure 4.
of existence of large-amplitude patterns. Our analysis is grounded in methods of phase-plane analysis that have proved useful25 in the study of the space-independent Oregonator ODEs (2.1-2.3) dx/dt = SO,- XY dy/dt = 1/s(-Y
+ x - qx2) - XY + fi)
dz/dt = W ( X
- Z)
(2.1) (2.2) (2.3)
These lines are called isoclines and divide the x-y plane into regions of positive and negative values of dx/dt and dyldt. Changes in z have been ignored even though z appears in eq 3.2 and dz/dt in eq 2.3. This is done because z changes slowly compared to x and y, and it is possible to consider z to be a parameter affecting only the relative location of the isoclines. Equations 3.1 and 3.2 are plotted in Figure 2 for z = 49.5.The circled arrows indicate the signs of dx/dt and dy/dt in various regions and predict the direction of time evolution there. There are three intersections of eq 3.1 and 3.2 at which dx/dt = dy/dt = 0. These are called nodes. The two outside nodes are locally stable. The inside node is infinitesimally unstable. The curved arrows in Figure 2 qualitatively indicate trajectories evolving from different initial values of x and y. All trajectories quickly approach one or the other of the stable nodes before z can change significantly. The trajectories cross the isoclines parallel to either the x or y axis, for either dx/dt or dy/dt will be zero at an isocline. The three crossings of the isoclines themselves enclose two connected lobes that will be most important in the following discussion. The limit cycle oscillations in the O r e g o n a t ~ r can ~ , ~ ~be 2. The system is essentially qualitatively t r a ~ k e d ~in ~ .Figure ~’ always on one or the other of the stable nodes. However, when it is on one stable node (say the lower right-hand one), z changes in such a way (increases) that this node moves toward the unstable node. (The actual movement of the nodes can be understood by moving the dy/dt = 0 isocline appropriately up or down relative to the fixed dx/dt = 0 isocline as z changes.) The first stable node eventually coalesces with the unstable node, and they both disappear. There is at this point only one crossing of the isoclines at the remaining stable node. The trajectory then jumps to that second stable node. Once there, z reverses direction (decreases)
It is helpful to first outline this phase plane view of the ODEs. ( 2 5 ) Troy, W. C., University of Pittsburgh, private communication.
( 2 6 ) Hasting, S.P.; Murray, J. D. SIAM J . Appl. Math. 1975,28,678. (27) Becker, P. K. M.S. Thesis, Department of Physics, University of Montana, 1983.
Concentration Patterns in the Oregonator Mode
The Journal of Physical Chemistry, Vol. 89, No. 1, 1985 121
1
Figure 3. A wind diagram in the x-y plane with two mappings. The isanemy lines are drawn for eq 3.1 and 3.2 with z = 49.5. The arrows now indicate the signs of d2x/dP and d2y/dP (the direction of the wind) rather than the direction of evolution of a mapping.
in such a way that the first stable node and the central stable node eventually reappear (the three intersections return), and, further, the second stable node starts to move toward the unstable node. Eventually the second stable node (in this case, the upper left one) coalesces with the unstable node, and the trajectory jumps back to the first node, and z again begins to increase. Oscillation occurs as this process repeats itself. This jumping back-and-forth between stable nodes as z changes illustrates the basic bistability of the BZ reaction and the Oregonator,28which will reappear in later analyses. If eq 2.1-2.3 had involved second rather than first derivatives with respect to time (e.g., d2x/dt2), the second derivatives could also be set equal to zero and the lines corresponding to eq 3.1 and 3.2 drawn in the x-y plane. However, they would now be lines of zero “acceleration” of “chemical velocity“ rather than lines of zero “chemical velocity”. The arrows in Figure 2 would not necessarily indicate the direction of motion but would indicate instead the direction of acceleration. Wind Diagrams. The basic ideas outlined above with time as the independent variable do not change much if space is the independent variable. Because we are searching for patterns that do not change with time, the time derivatives in eq 1.1-1.3 are set to zero, leading to a set of ordinary differential equations in space -D:(d2X/d12)
=
- XY + y - qx2)
-D’,,(dzy/d12) = 1/S(-Y -D:(d2Z/dP)
(4.1)
- XY + f i )
= W(X- Z)
(4.3)
The lines corresponding to eq 3.1 and 3.2 can still be drawn as shown in Figure 3, but the circled arrows are now reversed from their direction in Figure 2. It is necessary in order to make use of Figure 3 to map the x and y profiles of a spatial pattern (from I = 0 to I = L, where L is the length of the container) onto the x-y plane. For example, Figure 5a is a mapping of the concentration profile in Figure 1 created by plotting onto the x-y plane the values of x and y at incremental distances, Al, along the spatial axis. The size of AI does not matter so long as it is small enough to reproduce the (28) DeKepper, P.; Boissonade, J. In ‘Oscillations and Traveling Waves in Chemical Systems”;Field, R. J., Burger, M., Eds.; Wiley-Interscience: New York, 1984.
desired detail. The direction of a mapping is determined by the concentration gradients at a particular point in space. Any pattern will map onto a line of finite length as chemical concentrations and the length of the 1 axis must be finite. The use of mappings in the analysis of oscillatory dynamics has been extensively developed by Winfree.14 Mappings such as just described may be used in conjunction with diagrams like Figure 3 to characterize the possible form of patterns and what kinds of initial conditions may evolve to patterns under numerical integration. The following results are basically conjectures. They seem reasonable and work numerically but have not been proven formally. We start by recognizing that a first derivative over space (a gradient) is a kind of velocity stating how fast a concentration changes as space is traversed. A second derivative over space is a kind of acceleration. The arrows in Figure 3 do not indicate the direction of a mapping but instead indicate the direction of a “wind” that affects the direction a mapping will be accelerated in. We call diagrams like Figure 3 “wind diagrams”. The lines of zero second derivative will be called “isanemy lines”. It is important to remember that z does change slowly as x and y change across a pattern. This change in z causes the y-isanemy line eq 3.2 to move relative to the x-isanemy line eq 3.1. Equations 3.1 and 3.2 may have either one or three intersections, depending upon the value of z. Wind effects must be considered in conjunction with the boundary condition that all concentration gradients must be zero at both walls of the container. We observe that slightly away from a wall the second spatial derivative of a concentration (wind) must be opposed to the concentration gradient. Otherwise, the gradient could not reach zero at the wall. This leads to the result that all values of z in the profile of a pattern must be in a range where the isanemy lines have three intersections. Furthermore, the mapping of a stationary concentration pattern onto the appropriate wind diagram, as well as any initial waveform that can evolve into a stationary concentration pattern, must start in one lobe of the wind diagram and end in the other. To see the basis of the above statements we start at 1 = 0, where all concentration gradients are required to be zero, and follow a pattern across the 1 axis. It is assumed in this argument that the value of z at all points in space is such that the wind diagram has three intersections. It is apparent from the direction of the arrows in Figure 3 that all mappings starting with zero gradient and outside of both lobes will be blown away from the isanemy lines. The wind can never change sign without crossing an isanemy line, so such a mapping can never be decelerated to reach the boundary condition of zero gradients at the other side of the container. This is illustrated by line AB in Figure 3. All mappings of patterns must therefore begin inside one of the lobes of the wind diagram. Now consider a mapping starting with zero gradients inside of a lobe. The arrows show that it will be blown toward the other lobe (as in line SS’). Such a mapping will leave the first lobe, though necessarily heading generally toward the other lobe as it does so. The wind outside the lobes will act to push the mapping away from both lobes, and if the wind is strong enough, the mapping may continue on a route similar to line AB in Figure 3. If the mapping does enter the other lobe, it will find the wind blowing roughly in opposition to its direction of motion. It may then eventually reach zero gradients, the required boundary condition at the other wall. The mapping SS’ in Figure 3 could correspond to a stable pattern, as it starts in one lobe and ends in the other so that the wind near both walls blows in the opposite direction to the gradient. The easiest way to remain at zero velocity is to be in zero wind. Thus, we expect to find the termini of mappings of patterns close to the outside nodes. Indeed, because of the basic bistability of the reaction kinetics, we expect to find x and y usually close to one of the stable nodes for the value of z at that spatial point even in these diffusion-affected systems. The value of z changes as a mapping traverses a wind diagram, and the y-isanemy line moves relative to the x-isanemy line. It is possible to show with arguments similar to those above that, if the isanemy lines lose their three crossings as z changes, the
122 The Journal of Physical Chemistry, Vol. 89, No. 1, 1985
mapping will find itself in a wind such that it can never reverse its gradient and thus never can meet the second boundary condition. Thus, z must have a value at every point in space such that the isanemy lines have three intersections. For the Tysonlo rate parameters a n d f = 1.82 there are three intersections of the isanemy lines only for 2.33 < z < 125.4, a rather narrow range of the possible values of z. The spatial variation of z across a concentration pattern can be kept small by making D, large. This is the reason that the pattern in Figure 1 is stable. Large increases in D, and Dy do not induce patterns unless D, has a value such that there are always three intersections of the isanemy lines. A second requirement for stable patterns is that z must always be between the maximum and minimum values of x. To see this necessity rearrange eq 4.3 to d2z/d12 = - W ( X
Becker and Field
X
a
s
9 -
- z)/D’,
If z grows or starts (Le., is at a wall) above xmX,then dz/dl must be positive or zero and, because d2z/d12 is also positive (eq 5 ) , growth of z and dz/dl continues indefinitely. The boundary condition of zero gradient for z at the other wall never can be met. then d2z/dI2 Similarly, if z falls or starts (at a wall) below xmin, must be negative and dz/dl must be also negative or zero. The other boundary condition again never can be met. A final requirement for pattern formation is quite real but hard to express explicitly. This requirement is that neither the width of a pulse (as in Figure 1) nor the width of the container can be too narrow. There must be sufficient distance for gradients to fall to zero at the walls. A qualitative feeling for the distances required was obtained by numerical experience. The lengths required for a particular calculation were extrapolated from previous calculations. Application of the preceding ideas can be illustrated using D, = 10 mm2/s and f = 2.5. Initial conditions for a numerical integration are set up as follows. There are three intersections of the wind diagram only for 2.33 < z < 125.4. An intermediate value of z ((125.4 - 2.33)/2 = 61.5) is chosen and the concentration of z is set to that value at all spatial points. A pulse is created against the left wall by setting x and y equal to their value at the right node of the wind diagram for all points at the wall and an appropriate distance out into the container. The rest of space is set to the x and y values of the other node. The width of the pulse is adjusted by trial and error. These initial conditions are shown in Figure 4. The smoothed pattern in Figure 1 stabilized upon numerical integration. It is apparent looking at the z profile in Figure 1 that the width-of-container requirement results from the space needed for dz/dl to reach zero. Lengthening the container does not affect the shape of the pulse, and no other pulses form for these parameter values. The value of z changes very little across the container (Figure IC) for the very large and unrealistic value of D, used in this calculation. The resulting pattern is thus of little practical interest because it is patterns in z (usually a highly colored metal ion) that can be detected most easily experimentally. Figure 5b shows the three intersections in the wind diagram for the minimum value of z in Figure IC. The diagram based on the maximum z is essentially identical. The phase-plane mapping of Figure 1, a and b, is shown in Figure 5a. It stretches almost exactly between the outside nodes of Figure 5b as expected. The basic bistability of the Oregonator manifests itself again. A similar pattern is found withf= 2.5 for D, = 1 mm2/s, but none is found for D, = 0.1 mmz/s. However, forf= 0.6 patterns appear for values of D, as low as 10 X lW3 mmz/s. The one shown in Figure 6 stabilizes f o r f = 0.6 and D, = 50 X mm2/s, and initialization in the usual manner. Figure 7c shows the mapping of this pattern, which is more complex than the one in Figure 5 but equally as instructive. Figure 7a is the wind diagram for the maximum z and Figure 7b is the wind diagram for the minimum z (note shift in ordinate range in this plot). Point A in Figure 7c corresponds to the left node in the wind diagram for zmin. As z increases (in real space z is small on the right side of the container so motion to increase z is actually from right to left) the wind diagram changes from that in Figure 7b to that in Figure
9
>a 0
“ i l i a
0 -I
91 -
olo
I 4.0
8!0
L
DISTANCE IN MM Figure 4. Initial waveforms, (a) log x , (b) log y , (c) log z , that grew under numerical integration to the stable pattern in Figure 1. This initial condition was constructed by using the prescription developed in the text.
7a. The mapping follows the left node as it slowly rises (AB in Figure 7c). The values of x and y exchange magnitudes when the pulse is reached as the mapping quickly jumps from the left-hand to the right-hand node. The long part of the mapping (BC) actually corresponds to a very short distance in space. The jog in the right-hand side of the mapping corresponds to the movement of the right-hand node in Figure l a as dz/dl slowly approaches zero. Every point in the mapping is close to the right-hand or the left-hand node of a wind diagram corresponding to the value of z at that point in space. No numerical integration yielded a pattern that did not follow this rule. One apparent exception turned out to be a numerical artifact; unlike the other
The Journal of Physical Chemistry, Vol. 89, No. I , 1985 123
Concentration Patterns in the Oregonator Mode
9 N
X J
9
T I
0.0
I
I .o
I
LOG X 2'o
0
I
3.0
0
-
t. (3
s
0 0
9 I
'0.0
1
I .o
I
LOG
I
I
I
41 0
8!0
1 3.0
Y
X2'0
Figure 5. (a) Mapping into the x-y plane of the x and y waveforms in Figure 1. The numbers indicate the distance along the I axis a particular point on the mapping corresponds to. Note that the ends of the mapping match very closely the position of the outside nodes in the wind diagram. (b) Wind diagram in the x-y plane for eq 3.1 and 3.2 with log z = 1.86.
(v
N (3
0
-1
patterns described here, this one disappeared on a finer grid. The parameters 0,and f were manipulated to investigate what types of patterns could be found, always following the requirements derived in the previous sections. Most patterns looked much like those in Figures 1 and 6 . However, because all single-spike solutions meet the boundary condition of zero gradients, patterns can be pieced together in containers suitably long. For example, it is always possible to form a pattern with one spike of twice the normal width in the center of the container and another pattern with a spike against each wall (Figure 8). A most interesting situation occurs for the parameter values f = 0.6 and D, = 10 X mm2/s. This is the smallest value of D, for which a pattern was found. The important feature here is that an initial single pulse located against a wall repeats itself until the entire container is filled. In all other cases only a single peak is formed regardless of the size of the container. Figure 9 shows the evolution of a repetitive pattern from a single pulse. This type of pattern formation bears considerable resemblance to theories of s o m i t ~ g e n e s i in s ~embryo ~ development. Equally interesting is the connection between repetitive, large-amplitude solutions and the quasi-stationary solutions to be discussed next.
Quasi-Stationary, Small-Amplitude Solutions and the Transition to Large-Amplitude Solutions Perturbation Analysis. Following Turk@ and Prigogine**we now seek small-amplitude solutions to eq 4.1-4.3 expressible as a perturbation to the spatially homogeneous state. The infinitesimal stability properties of the spatially homogeneous steady state subject to diffusion have been investigated in two spatial dimensions by Greenbaum30 for parameter values for which the well-stirred steady state is stable. Real, positive eigenvalues were
(29) (a) Me, D. "An Introduction to Developmental Biology"; Wiley: New York, 1978. (b) Lipton, B. H.; Jacobsen, A. G. Den Biol. 1974,38, 73.
(30)Greenbaum, N. N. Proc. Model. Simul., Annu. Pittsburgh Conf., 14th 1983, 6, 989.
.. 0010
L
DISTANCE IN MM Figure 6. Stationary pattern obtained for D, = 50 X mm2/s,f = 0.6, and L = 10 mm. 1125 spatial grid points were used in the integration. (a) log x , (b) log y , (c) log z.
124
Becker and Field
The Journal of Physical Chemistry, Vol. 89, No. 1, 1985
>
9 -
c3
0
9
0
0
I
LOG X
A 9 0
S - I
' olo
1
0.4
I O!8
DISTANCE IN MM 1
I
-
00
\
3.d
I
1.0
LOG X
Figure 8. Large-amplitude pattern with pulses against both walls. The boundary condition of zero gradients is met at each wall and at the center of the container. f = 1, D, = Dy = 10 X mmz/s, D, = 1 X lo-) mm*/s, and L = 0.8 mm. Note that decreasing all of the diffusion coefficients simply shortens the length scale of the pattern.
vestigate in one spatial dimension the form and stability of such solutions. The perturbation analysis is standard. Details are given in the Appendix. The analysis leads to the conclusion that there are eigenvalue solutions to eq 4.1-4.3 corresponding to small-amplitude pattern formation. The form of these solutions for x is x = xo cx, cos (g,l) (6)
> c3
+
0 A
\
01
S I -
lO!O
0.0
1
I
1.0
210
3.0
LOG X Figure 7. Mapping and wind diagrams for the pattern shown in Figure 6. (a) Wind diagram for the maximum z in Figure 6 (log z = 2.613 at the left-hand side of the concentration profile). (b) Wind diagram for the minimum z in Figure 6 (log z = 1.153 at the right-hand side of the concentration profile). (c) Mapping of the x and y waveforms. The numbers refer to particular points on the 1 axis in Figure 6.
obtained, indicating the presence of solutions corresponding to pattern formation. Kopell and Howard3I worked under conditions where the well-stirred steady state is unstable and found evidence of time-dependent, spatially inhomogeneous solutions. We in~~
(31) Kopell, N.; Howard, L. N. Lect. Mark Life Sci. 1974, 7, 201.
with similar solutions for y and z . The quantity xo is the homogeneous steady-state value of x, C is an expendable parameter that determines the amplitude of the pattern, and g, is an eigenvalue with associated eigenvector element, x,. Such solutions exist only for containers of length L , = n r / g , with n = 1 , 2 , 3, etc. Small-amplitude solutions exist for the parameters used here (0, = D, = D, = 1 X lo-' mmz/s) only when 0.61 < f S 1.82. Requiring gm& = n r (n = 1) the critical & forf= 1.82 are (there are two eigenvalues) 0.41458025 and 0.58721092 mm. The solution for L, = 0.58721092 mm with C arbitrarily set to 4.8 X is shown in Figure 10. Stability of Small-Amplitude Solutions. Stability of smallamplitude patterns is checked by using them as the initial condition for a numerical integration of the original PDEs (eq 1 .l-1.3). All small-amplitude solutions are ultimately unstable to numerical integration, although some live for such a long time that we term them quasi-stable. This behavior stands in contrast to the apparent stability of similar solutions of the Brusselator modelI5 with the Dirichlet boundary condition. Time evolution of almost all small-amplitude solutions tested eventually leads to a state where large-amplitude, spatially ho-
Concentration Patterns in the Oregonator Mode
The Journal of Physical Chemistry, Vol. 89, No. 1 , I985
t
t
-
N
t Io
“0 X
X
(0
r!
N
I
I
1 :I
T=5000
I
1
1
r! N
t
r
1
11 0.0
8 .O
16.0
DISTANCE IN MM Figure 9. Growth of a single pulse to a multiple-peak pattern. The value of x is shown at t (scaled) = 0, 100, and 5000, and z is shown at t = 5000. Note that the waveform is asymmetric because the length of the container does not correspond to an integral number of peaks. f = 0.6, D, = 10 X 10” mmz/s, and L = 17 mm.
mogeneous oscillation in x , y , and z occurs. However, the time required for this homogeneity to occur depends upon how close the length of the container, L, is to L,, and upon the shape and
I
I
125
126 The Journal of Physical Chemistry, Vol. 89, No. 1, 1985
amplitude and shape close to but not identical with a small-amplitude solution with C = 4.8 X loy4. Lifetime calculations for both curves begin after their amplitudes exceed the 4.8 X level. The curves look similar to the eye at the starting point, yet there is a great difference in the lifetimes of the two waveforms. The one that grew to its lifetime starting point is quasi-stable, living about 4.5 X lo8 s, while the one started at its lifetime level lives only about 4000 s. Apparently a small difference in the shape of curves at a given amplitude has a large effect on stability. It is noteworthy that the addition of second-order terms to the perturbed solution does not significantly affect calculated lifetimes. Initial conditions with small C, even if very different from the cosine form (e.g., x = xo, y = yo C, z = zo C at all points in space), evolve quickly with L = L, to the quasi-stable cosine waveform corresponding to n = 1. The well-stirred steady state is unstable for the parameter values used above and temporal oscillation around that state are equivalent to the final state reached in the evolution of smallamplitude solutions. The question arises of whether a small-amplitude solution would be stable if the well-stirred steady state were stable. Such a state, where a small-amplitude solution also exists, can only be obtained by returning to the original FieldNoyeslo rate parameters. Even then, small-amplitude solutions turned out to be at most only quasi-stable, decaying to a spatially homogeneous steady state. Multiple-peak small-amplitude solutions are possible with D, = Dy = D, and containers of appropriate length, but they are shorter lived than the n = 1 solutions. Multiple-peak curves are unstable and tend to gain or lose a peak under numerical integration before flying off to a spatially homogeneous state. These multiple-peak small-amplitude solutions are important, however, because under some conditions they can grow into large-amplitude patterns. There is a small-amplitude solution when f = 0.6 and D, = 10 X mmz/s with L,values of 0.333 518 1 and 2.139000 mm. Any small perturbation of the homogeneous state in a container of L = L, evolves to a one-half cosine quasi-stationary waveform. Perturbations in a container of L = 2L, tend to evolve to a full cosine waveform. A container of L = 8L, produces a four-cycle waveform from a perturbation of appropriate cosine form and a three-, four-, or five-cycle waveform from a random perturbation. Figure 11 shows the rapid growth of a very small, four-cycle cosine perturbation to a four-cycle large-amplitude pattern that is similar to, but not identical with, the one in Figure 10 that grew from a single pulse. Ideal initial conditions were used to make this happen. However, because any small random perturbation grows to a small-amplitude waveform of roughly the correct number of peaks, it is likely that any small perturbation will grow to a large amplitude pattern. The special (still not well mm*/s understood) character of the f = 0.6 and D, = 10 X parameters, plus the similarity of L, for the small-amplitude solutions and the half-wavelengths of the large amplitude solution, creates a situation where large-amplitude patterns are almost certain to appear. The mechanism of their appearance would be random perturbations to the spatially homogeneous state.
+
+
Significance to Experiments Repetitive, stationary patterns similar in cross section to those in Figures 10 and 11 have been observed in the metal-ion-catalyzed BZ reaction by Zhabotinskii and Zaikin.17 They interpreted them as resulting from a diffusive instability caused by differences in diffusion coefficients of various intermediate species. They were able to model their results using a two-intermediate (variable) model with a cubic term in one differential equation. The diffusion coefficient of one intermediate was set equal to zero and that of the other intermediate (whose behavior was described by the cubic differential equation) was set to a normal value. Similar theoretical work was done by Kawcyznski and Zaikin'* using a three-variable model very similar in form to the Oregonator, but unrelated to any real chemistry. They allowed only one intermediate to diffuse and carried out a qualitative analysis similar to that presented here which indicated the existence of patterns supported only by a diffusive instability.
.
1
Becker and Field
T = 5000
T=5000
i, T=80
1 T=25
DISTANCE IN MM Figure 11. Growth of a small-amplitude, first-order perturbation solution to a multiple-peak, large-amplitude pattern. In this case the length of the container does correspond to an integral number of peaks, and the waveform is symmetrical. f = 0.6, D, = 10 X mm2/s, and L = 17.042. Scaled times are given.
Concentration Patterns in the Oregonator Mode
.
b c.
The Journal of Physical Chemistry, Vol. 89. No. I, 1985 121 be done by changing other parameters in the Oregonator in terms of eliminating the need for large differences in diffusion coefficients. The insights developed here will be most useful in searching for such sets of parameters. There is nothing implicit or explicit in our analysis which seems to require any difference in diffusion coefficient values for pattern formation. This could be of theoretical importance as most theory has been done under the assumption that there must be diffusion coefficient differences for pattern formationI5 to occur in reaction-diffusion systems. Acknowledgment. We thank Professor Richard J. Hayden for his interest in and contribution to this work and the University of Montana Computer Center for computing facilities. This work was partially supported by the National Science Foundation under Grant CHE-8023755.
Figure 12. Experimentallyobserved stationary concentration patterns in the system with initial conditions: [H2S0,J0= 0.189 M,[NaBrOJ, = 9.51 X IO-' M, [KBr], = 7.74 X M, [ferroinJo= 1.21 X lo-' M, and (C6H,(COOH),], = 8.27 X IO-*M. Temperature = 25 i 0.2 OC.
The photograph width is IO cm. Reproduced with permission from ref 19. Copyright The American Institute of Physics (1980). Later workers have observed patterns in bromate-containing ~ ~ a ferreagents related to the BZ reaction. S h ~ w a l t e rused roin-bromate system in which 4-cyclohexene-1,2dicarboxylicacid was used to absorb the bromine produced by the reaction. This system is not oscillatory. Orb5nZ0worked with an uncatalyzed bromate o s c i l l a t ~ r . ~ In . ~ ~both cases the observed phenomenon is remarkably similar mosaic patterns (Figure 12) that form spontaneously and remain stationary for several minutes. Presumably the limited lifetime of the patterns results from changes in reagent composition as chemical reaction occurs in a closed system. As in the case of the Zhabotinskii-Zaikin" patterns, the spatial scaling is of the order seen in our calculations, 2-3 mm. However, these authors probably correctly suggest that pattern formation does not result solely from a diffusive instability. They point to convective instability resulting from heat released by the chemical reactiod3 and/or gas transfer (bromine and/or oxygen) across the surface of the reagent. The potential importance of oxygen has been discussed by Winfree." The pattern found by Showalter did not depend upon contact with atmospheric oxygen while the one found by OrGn did. The kinetics of the BZ reaction do depend upon oxygen concentration.2 The effect is in the Brregeneration from brominated organic material, a step that is absent in the Showalter system. The effective value offin the BZ reaction is affected by oxygen.2 A periodic precipitation p h e n ~ m e n o n ~also ~ " ~has been mentioned.20 There seems to be a strong need to reinvestigate these systems now that a theoretical basis for reaction-diffusion based pattern formation in the BZ reaction exists. Some such experiments have recently been carried out by Wood and Ross.37 Experiments in gels where convection is eliminated and effective diffusion coefficients can be varied would seem to be a good place to start. While most of the large-amplitude patterns reported here formed with very large values of 0, compared to 0, and Dv, the most interesting (because of the formation of a repetitive pattern) appeared with the smallest D, used, 10 X mm2/s. While differences of diffusion coefficients of the order of a factor of ten are quite possible between small molecules and macromolecules, it is not likely that such a difference exists among species important in the BZ reaction. However, we have not yet explored what can
M.;Koros, E. J. Phys. Chem. 1978.82, 1672. (33) Gutkowicz-Kaufman, D.; Collins, M.A.; Ross, J. Phys. Ffuids 1979, 22, 1451. (34) Winfree, A. T. In 'Theoretical Chemistry: Advances and Perspectives"; Eyring, H., Henderson, D., Eds.; Academic Press: New Yotk, 1978; Vol. IV. (35) Flicker, M.;Ross,J. J . Chem. Phys. 1974,60, 3458. (36) Ortoleva, P.; Schmidt, S. In 'Oscillations and Traveling Waves in Chemical Systems"; Field, R. J., Burger, M.. Eds.; Wiley-Interscience: New York, 1984. (37) Wood, P. M.;Ross, J. J. Chem. Phys., in press. (32) Orban,
Appendix: Perturbation Analysis of Eq 4.1-4.3 We seek stationary, small-amplitude solutions to eq 4.1-4.3 expressible as a perturbation to the spatially homogeneous steady state. The analysis proceeds by substituting x = xo + x*, y = yo + y*, and z = zo + z* into eq 4.1-4.3. The quantity (xo,y@o) IS the spatially homogeneous solution given by
xo = Kl
-f- d + [(I -f - ~ 1 +) ~4 d f +
1)11'21/(29)
Yo = fxo/(l + xo) 20
(AI) (A2)
= xo
(A31
The quantities x*, y*, and z* are a first-order perturbation to the spatially homogeneous solution, and it is assumed that Ix*l
