2368
J . Phys. Chem. 1990, 94, 2368-2314
Coupling of Identical Chemical Oscillators K. Bar-Eli School of Chemistry, Raymond and Beverly Sackler Faculty of Exact Sciences, Tel- Aviv University, Tel- Aviv 69978, Israel (Received: April 14, 1989; I n Final Form: October 16, 1989)
Identical Belousov-Zhabotinskii oscillators, in a CSTR, modeled by the Field-Koros-Noyes (FKN) mechanism, are coupled in a diffusionlike manner. In addition to the obvious symmetric solutions, Le., solutions in which both CSTRs are oscillating in unison, or are in the same stable steady state, unsymmetric, broken-symmetry solutions may coexist under the same set of constraints. Thus, depending on constraints and initial conditions, the combined system can be in the following states: (a) stable symmetric steady state; (b) symmetric oscillations, when both cells oscillate in phase; (c) coexistence of symmetric and unsymmetric steady states; (d) coexistence of symmetric oscillations and unsymmetric stable steady state (broken symmetry); (e) coexistence of symmetric and unsymmetric oscillations. The latter differ from the former in phase, in amplitude, and in period. On the other hand, no unsymmetric oscillations were found to coexist with the symmetric steady state. All the initial conditions tried ended in either of the two, possible, stable states. The change of periods and amplitude of both types of oscillations are examined as a function of the system constraints, namely, concentrations and coupling rate.
Introduction The problem of coupling two chemical oscillators was first dealt with by Prigogine and Lefever,] who showed that two Brusselators when coupled in a diffusionlike manner can “break symmetry” and arrive at an unsymmetrical steady state. Bar-Eli2-5has shown that this is not unique to the Brusselator model, but that also other chemical oscillators can do so as well. The examined oscillators were (1) the simplest chemical oscillator6 resulting from the Noyes-Field-Thompson mechanism,’ (2) simple “Oregonator” model,s (3) the “Brusselator” model,] (4) an autocatalytic model by K ~ m a r(5) , ~ the Lotka-Volterra population model,I0 and ( 6 ) the Belousov-Zhabotinskii oscillator (BZ)” governed by the Field-Koros-Noyes mechanism (FKN).I2 All of them showed a coupling range in which the oscillations stopped, the symmetry was broken,’ and a stable steady state coexisted with the oscillations. Experimental verification of these results was given by Marek and ~ o - w o r k e r s ~and ~ - by ~ ~Bar-Eli and R e ~ v e n i .The ~ former authors used also coupled cells in a hexagonal configuration.14 Earlier, Tyson and Kauffman16 have used diffusionlike coupled Brusselators as a model for controlling the mitotic cycle in the acellular slime mold Physarum polycephalum under conditions where the two chemical species of the Brusselator diffuse at different rates. They obtained, in addition to the unsymmetric steady states, also, unsymmetric limit cycles, at the expense of the stability of the symmetric one. Comparatively little work has been done on these additional limit cycles. Field and Crowleyl’ have coupled two separated chemical (I) Lefevre, R.; Prigogine, I. J . Chem. Phys. 1968, 48, 1695. (2) Bar-Eli, K. J . Phys. Chem. 1984, 88, 3616. (3) Bar-Eli, K. J . Phys. Chem. 1984, 88. 6174. (4) Bar-Eli, K. Physica 1985, 1 4 0 , 242. (5) Bar-Eli, K.; Reuveni, S. J . Phys. Chem. 1985,89, 1329. (6) (a) Bar-Eli, K. In Nonlinear Phenomena in Chemical Dynamics; Vidal, C., Pacault, A.. Eds.; Springer Series in Synergetics 12; Springer-Verlag: Berlin, New York, 1981; pp 228-239. (b) Geiseler, W. J . Phys. Chem. 1982, 86, 4394. (c) Orban, M.; DeKepper, P.; Epstein, I. R. J. A m . Chem. Soc. 1982, 104, 2657. (d) Bar-Eli, K.; Geiseler, W. J . Phys. Chem. 1983, 87, 3769. (7) Noyes. R. M.: Field, R. J.: Thompson, R. C. J . Am. Chem. Soc. 1971, 93, 7315. (8) Field, R. J.; Noyes, R. M. J . Chem. Phys. 1974, 60, 187. (9) Kumar, V. R.; Jayaraman, V. K.; Kulkarni, B. D.; Doraiswamy, L. K. Chem. Eng. Sci. 1983, 38, 673. (IO) (a) Lotka, A. J . Phys. Chem. 1910, 14, 271. (b) Lotka, A. J . A m . Chem. SOC.1920, 42, 1595. (c) Volterra. V. Lecons sur la theorie Mathematique de la Lutre pour la Vie; Gauthier-Villars: Paris, 1931. ( 1 1 ) (a) Belousov, B. P. Sb. Ref. Radiat. Med. 1958, 145. (b) Zhabotinskii, A. M. Dokl. Akad. Nauk. S S S R 1969, 157, 392. (12) Field, R. J.; Koros, A.; Noyes, R. M. J . Am. Chem. Soc. 1972, 94, 8649. (13) Marek, M.; Stuchl, I . Biophys. Chem. 1975, 3241. (14) Stuch!, I . ; Marek, M. J . Chem. Phys. 1982, 77, 2956. (IS) Dolnik, M.; Marek, M. J . Phys. Chem. 1988, 92, 2452. (16) Tyson, J.: Kauffman, S. J . Math. B i d . 1975. I . 289
0022-3654/90/2094-2368$02.50/0
oscillators by coupling electrically the cerium ions via suitable voltage and current control. This coupling introduces logarithmic terms into the dynamical equations, unlike the linear coupling terms introduced in the previously described work. Recently, Hocker and Epstein’* coupled in a diffusionlike manner a two-variable model proposed by Boissonade and De Kepper.I9 The coupling rates were equal although the oscillators were not, Le., not with the same parameters. A variety of oscillations, including chaotic ones, and steady states were obtained. Chaotic oscillations were also obtained by coupling nonisothermal
oscillator^.^^ On the experimental side, Crowley and EpsteinZohave coupled two nonidentical BZ reactors and obtained experimental evidence relevant to the results obtained here. In this work, unlike the previous one, we are concerned with the coupling of two identical CSTRs (continuous flow stirred tank reactor), both containing BZ reactants and modeled by the FKN mechanism. In addition to the identity of the CSTRs, the coupling rates between them are also the same, thus introducing more symmetry into the equations. The FKN equations for the single CSTR were studied extensively, and in particular, De Kepper and Bar-EliZ1compared the experimental results to the theoretical predictions. Regions of monostability, bistability, and oscillations, simple and of the bursting type,22occurred in the same regions of the constraint space as predicted theoretically. The coupled CSTRs are examined (a) at constant equal concentrations as a function of the coupling rate, (b) at constant coupling rate when the concentration of malonic acid in the feed is changing, and (c) at constant coupling rate when the bromide ion concentration in the feed is changing. In most cases both coupled cells are in the oscillation region, but a few results are also given for other regions.
Model and Computations The model investigated is the Field-Koros-Noyes model: Br0,- + Br- + 2H’ s H B r 0 2 + HOBr k , = 2.1 M-3 s-’ k-, = 1 X lo4 M-l s-I
(FKN) (1)
(17) (a) Field, R. J.; Crowley, M. F. J . Phys. Chem. 1984, 88, 762. (b) Crowley, M. F.; Field, R. J. In Nonlinear Phenomena in Chemical Dynamics; Vidal, C., Pacault, A., Eds.; Springer Series in Synergetics 12; SpringerVerlag: Berlin, New York, 1981; p 147. (18) Hocker, C. G.; Epstein, I. R. J . Chem. Phys. 1989, 90, 3071. (19) Boissonade, J.; De Kepper, P. J. Phys. Chem. 1980, 84, 501. (20) Crowley, M. F.; Epstein, I. R. J . Phys. Chem. 1989, 93, 2496. (21) De Kepper, P.; Bar-Eli, K. J . Phys. Chem. 1983,87, 480. (22) (a) Bar-Eli, K.; Noyes, R. M. J . Chem. Phys. 1988, 88, 3646. (b) Roux, J . C.; Simoyi, R. H.; Swinney, H. L. Physica 1983, 8D, 257. (c) Maselko, J.; Swinney, H. L. J . Chem. Phys. 1986, 85, 6430. (d) Hudson, J . L.: Mankin. J. C. J . Chem. Phys. 1981, 7 4 6171.
63 1990 American Chemical Society
The Journal of Physical Chemistry, Vol. 94, No. 6, 1990 2369
Coupling of Identical Chemical Oscillators HBrO,
k2 = 2
+ Br- + H+ F! 2HOBr
Br03-
+ HBrO, + H +
k4 = 1 X IO4 M-, Ce3+ + BrO,'
k , = 6.5 Ce4+
k-, = 5
lo9 M-, s-I
X
X
lo-, M-I s-I
+ H20
2B1-0,'
k-, = 2
s-l
X
lo7 M-'
(3) s-l
+ H + e Ce4+ + HBrO, k-, = 2.4
lo5 M-2 s-l
X
IO7 M-l
(4) s-l
+ Br02' + H 2 0 e Ce3+ + Br0,- + 2H+
k6 = 9.6 M-' s-' X
lo7 M-'
Ce4+
+ MA
-
(5)
k-6 = 1.3 X IO4 M-3 s-I
2HBr0, e Br03-
k7 = 4
F!
X
(2)
+ HOBr + H+
k-7 = 2.1
s-l
gBr-
k6 = 0.53 M-' s-'
X
lo-', M-,
+ Ce3+ + products
(6) s-l 2211
(7)
Is,
"0
+
Results Single Cell. The results of computations and experiments of the FKN mechanism of a single CSTR are given in previous publications.21*22We shall review here the main results pertaining to the present work. A single CSTR has, depending on the constraints, regions of (1) single stable steady state with comparatively high concentration of bromide ions (with comparatively low concentration of ceric ions), SSII, (located at higher values of [Br-1, and not shown in Figure 1); (2) single stable steady state with low concentration of bromide ions (and high concentration of ceric ions), SSI (region D); (3) bistability region in which both stable steady states coexist (with hysteresis between them), region (23) Showalter, K.; Noyes, R. M.; Bar-Eli, K. J . Chem. Phys. 1978,69, 25 14. (24) (a) Tyson, J. J. In Oscillations and Traveling Waves in Chemical Systems; Field, R. J., Burger, M., Eds.;Wiley: New York, 1983. (b) Bar-Eli, K.: Ronkin, J. J . Phys. Chem. 1984.88, 2844. (c) Forsterling, H. D.: Lamberg, H.; Schreiber, H. Z . Naturforsch. 1980, 3 5 4 329. (25) (a) Gear, C. W. Numerical Initial Value Problems in Ordinary Differential Equations; Prentice-Hall: Englewood Cliffs, NJ, 1971; pp 209, 229. (b) Hindmarsh, A. C. Gear: Ordinary Differential Equation System Solver, UClD 30001 Rev. 3: Lawrence Livermore Laboratory, University of California: Livermore, CA, December 1974. (26) Hassard, B. D.; Kazarinoff, N. D.; Wan, Y. H. Theory and Applicarion of Hopf Bifurcation; Cambridge University Press: Cambridge, 198 1. (27) (a) Doedel, E. J. Cong. Num. 1981,30, 265 (Proc. 10th Manitoba
Conf. on Numer. Math. Computing). (b) Doedel, E. J. 'Auto: Software for Continuation and Bifurcation Problems in Ordinary Differential Equations". Applied Mathematics Department, California Institute of Technology, 1986.
-1 )
Lim I
40
g = 0.4615
A somewhat simplified version of it was used by Showalter et al.23 Some of the rate constants used were criticized lately,,, but in order to compare with previous data, we have kept the original rate constants as in ref 21. The dynamics of the system is determined by the set of equations C(i) = R(C(i) ko(Co(i)- C(i)) k,(C(j) - C(i)) (8) where i a n d j take the values 1 and 2. ko is the flow rate through the CSTR, k , is the coupling rate, and the function R of concentrations C is given by the chemical mass action rates as derived from the mechanism. The first two terms on the right-hand side specify the dynamics of a single CSTR, while the third one gives the flow of matter from one CSTR to the other. This flow is proportional to the difference in concentrations between the two cells and thus is of diffusionlike nature. The number of species in the above mechanism is nine, since water is taken as constant, out of which only five are independent. In the double system there are, therefore, 10 independent variables. The steady states, Le., C(i) = 0, were solved by Newton's method, their stability was determined by solving for the eigenvalues by the standard package (EISPACK), and the equations themselves were integrated by Gear's method.,, Program BIFORZ by Hassard et a1.26was used to test the behavior near the Hopf bifurcation points. Continuation package AUTO^^ by Doedel was used to continue the oscillations after the bifurcation.
+
I
120
80
160
[MA], /mM
Figure 1. Portion of [MA]o-[Br-]o subspace of constraints. Letters show various regions of steady states (marked by roman numerals). Roman numeral in parentheses: steady state is unstable. Regions E and F are oscillation regions. Solid line: hysteresis limit. Dashed line: Hopf bifurcation. Point S : standard conditions (see text). Arrows show directions of constraint change other than coupling rate. -4 4 r
-8.41'A' 0
'
100
'
'
200
'
I
300
'
'
400
'
I, 50
TIMEIS
Figure 2. Typical time plot of an oscillation of a single CSTR (or symmetric oscillation of coupled CSTRs). Shown is bromide ion concentration (log scale) vs time. Standard conditions.
A; (4) oscillation regions, in which there are either one or three unstable steady states, regions E and F; and (5) regions with three steady states, only one of which is stable, regions B and C. The various regions are shown in Figure 1. The oscillation region in the [Br-]o-[MA]o subspace of constraints is limited on the left and above, Le., at low values of [MA], and high values of [Br-I,, by Hopf bifurcations. The effect of these will be discussed below. In this work two identical CSTRs, both in the oscillation region, are coupled. We use as standard conditions, point S in Figure 1, the following constraints: [BrO3-Io= 30 mM, [Br-lo = 20 kM, [H+], = 1.5 M, [Ce3+Io= 3 mM, [MA], = 32.06 mM, and ko = 3.92 X s-l. This set of constraints position the system in the oscillation region E. At this point the system has one unstable steady state with two eigenvalues with positive real part. Relaxation oscillations occur with a period of 477 s and amplitude (ratio between maximum and minimum of bromide ion concentration) of 5250. Other constraints may result in more complex oscillations.22 The main features of these oscillations, shown in Figure 2, are slow increase of [Br-] from its minimum, at point A, t o point B, sharp increase to maximum (relative slope of 6) between points B and C, slow decrease to point D, and a very sharp fall t o minimum concentration (relative slope of 295) to complete the cycle. The slow increase and the slow decrease in concentration, take about 0.4 and 0.6 of the period, respectively.
2370 The Journal of Physical Chemistry, Vol. 94, No. 6,I990
Bar-Eli
TABLE I: Bifurcations at Standard Conditions (SeeText) symmetric solution
-
P
4 3“ 1.009028
k,ls-’
unsymmetric solution
-
3 2” 3.691 513
X
3 4 1 4.49037 X lo--’ 4-2 8.079681 X 2 - lb 5.314429
P
k,/s-‘
P k,/s-’
1-3 2.002702 X IO-) 2 oc 1.152978 X lo-) 1-3 3.695 705
-
-
3 4b 5.42435 X 0 2‘ 1.461 496
-
“Points where the unsymmetric solution bifurcates from the symmetric one. bTurningpoints. Between these values of coupling, the unsymmetric solution is stable (heavyarrows in Figure 3) -
,
10-7
i
I
1 UNSYM
-5.0-
I
-
-
N
-6.0-
-
.h
m
v
01 0
, 0 - 8 1 10-5
1
10-4
I
I
IO-*
-7.0-
IO-’
I00
IO’
k,/S-’
--
Figure 3. Steady-stateconcentrationsof bromide ion in the first CSTR (log scale) vs coupling rate (log scale). Arrows: points where p = 4 2. 2 or p = 3 1. Heavy arrows: Hopf bifurcations where p = 0 Dots: bifurcation point of unsymmetric solution from the symmetric one.
-
Standard conditions. The combined system is examined when both cells are at point S and the coupling rate is changed; then the coupling rate is fixed and point S is moved horizontally (changes in malonic acid concentration) and vertically (changes in bromide ion concentration). Changes in Coupling Rate: Steady States. The steady-state concentrations ([Br-1) of cell 1 are shown in Figure 3 vs the coupling rates together with the number p of the positive real part of the eigenvalues, which indicates the stability of the steady state (p = 0, stable: p > 0, unstable). Table I shows the change points i n p. The following features are observed: (a) The first is a symmetric solution in which the concentration of both cells are equal and are the same as that of the single cell. (This is the straight line with [Br-] = 9 X IO-* M for all k,). This solution is unstable for all values of coupling rates. The value of p for this line changes from 4, at low k , values, to 3, at medium ones, and finally to 2 at high coupling rates. This trend is easily explained as follows: at low coupling rates, the system is made essentially of two independent identical cells, each of which is unstable with two eigenvalues having positive real parts. The total number of eigenvalues with positive real parts is thereforep = 4. At very high coupling rates, the system behaves essentially as one “big” CSTR (identical with its two component parts) and thus has only two positive real eigenvalues p = 2. One should bear in mind that as the single cell has gone from region D to region E (Figure l), p has gone from 0 to 2 via a Hopf bifurcation, shown by the dashed line in Figure 1. At medium coupling rates, p goes from 4 to 3 and then from 3 to 2. As p changes by one unit, one eigenvalue becomes zero, and as is shown in ref 28, another solution, an unsymmetric one, must bifurcate from it. As we shall see later, the unsymmetric solution need not be connected to the symmetric solution and can form an kola separated from it.% (b) The second (28) loos, G.; Joseph, D. D. Elementary Stability and Bifurcation Theory; Springer: New York, 1981. (29) (a) Gray, P.; Scot, S.K. Chem. Eng. Sci. 1984,39, 1087. (b) Uppal, A . ; Ray. W. H.; Poore. A . B. Chem. Eng. Sei. 1976,31, 205.
Figure 4. Projection of phase-space trajectory on [Br-],-[Br-I2 subspace for unsymmetric oscillations at k, = 5 X s-I and standard conditions.
feature observed is two unsymmetric solutions with higher and lower bromide ion concentrations as shown in Figure 3. Because of the symmetry of the problem, if the first cell will be at the high concentration, the second one will be at the low concentration and vice versa. Thus the lower branch, with concentration very near the symmetric solution, depicts either the concentration of cell 2 or the concentration of cell 1 of the other solution. Following this unsymmetric branch, various changes of p are seen (Table s-I and k, = I). In particular, between k , = 11S 2 9 78 X 1.614 96 s-l, (at which coupling values Hopf bifurcations occur) p = 0; namely, all eigenvalues are negative and the unsymmetric solution becomes stable. Between these coupling values the system may exist in an unsymmetric stable steady state as was also shown experimentally in previous work.5~’3-i5*20 Efforts to determine whether the Hopf bifurcation is super or subcritical, by using the BIFORZ code,26failed. The error limits obtained were much greater than the values themselves. For this was not able to continue the oscillations beyond reason AUTO the bifurcation. Changes in Coupling Rate: Oscillations. At zero coupling rate, the two cells oscillate, as shown in Figure 2 , independently with arbitrary phase difference between them. The combined trajectory will be on a torus that is the direct product of the two independent limit cycles. As soon as the coupling is increased from zero, all the phase shifted limit cycles become unstable except two-a symmetric one with zero phase difference and an unsymmetric one with phase difference of 0.5. A typical phase trajectory of the unsymmetric oscillations is shown in Figure 4. A similar figure can be seen in the work of Mankin and Hudson.32 This is, of course a projection of the phase trajectory in the two indicated dimensions, namely, the Br- concentrations in the two cells. The symmetric phase point trajectory will be a straight diagonal (45’ slope) line in this subspace. The unsymmetric oscillations have a constant phase shift of 0.5 between the cells; Le., the maximum (minimum) of one cell appears half a cycle after the maximum (minimum) of the other. A trajectory starting with an arbitrary phase shift between the cells will develop into one of these limit cycles, namely, the symmetric or the unsymmetric one, depending on the initial phase, the phase difference, and the coupling rate. At very high coupling
The Journal of Physical Chemistry, Vol. 94, No. 6,1990 2371
Coupling of Identical Chemical Oscillators TABLE 11: Range of Phase Shifts Leading to Unsymmetric Oscillations
range of phase diff of 2nd cell" at IO4k,/s-' phase of 1st cell
0.2 0.073-0.935 0.231-0.922 0.098-0.940
0 min of [Br-Ib 0.43 max of [Br-IC 0.76 in midst of slow decr of [Br-Id
1 0.045-0.976 0.027-0.966 0.026-0.972
IO 0.1 12-0.987 0.027-0.953 0.016-0.979
5 0.080-0.987 0.016-0.978 0.016-0.979
" Initial point in this range will go over to unsymmetric oscillations, outside the range, to symmetric oscillations, at the indicated coupling rates. bZero phase is taken arbitrarily as the minimum of [Br-1, point A of Figure 2. CPointC in Figure 2. dBetween points C and D in Figure 2. TABLE 111: Unsymmetric Oscillations below Hopf at k, = 11.52978 X 104s-', Standard Conditions I 04k,/s-1 period/s amplitude" 104k,/s-' period/s amplitude" 11.5 11.4867 11.4634 11.4265 11.4 11.368 11.3 11.2 11.1 11.08 11.06 1 1.05 II 10.8 10.6 10.4 10.2
IO 9.8
7415 6245 5130 4285 3940 3628 3210 2860 2625 2600 2590 2575 2500 2250 2050 1900 1845 1825 1675
8910 8810 8910 8910 8810 8910 8810 8810 8910 8910 8910 8910 8910 8910 8910 9170 9170 9440 9720
9.5 9 8 7 6.31 5 4 3 2 1.5 1 0.9 0.8 0.7 0.6 0.5 0.4 0.2
1650 1550 1315 1209 1122 1030 934 843 754 683 606 585 570 555 540 528 516 494
10000 10000 10000 11220 11900 12000 11550 10900 9720 8660 7080 7080 6680 6310 6130 5960 5620 5700
"Amplitude of symmetric oscillations is 5250. rates, the system will, obviously, sustain only symmetric oscillations, regardless of the initial phase or the phase shift. Thus there is an intermediate region of phase shifts and coupling rates, which depends on the initial phase; starting inside this range, the system will go over to the unsymmetric oscillations. Starting outside this range, symmetric oscillations will ensue. The exact boundaries of this range are cumbersome to determine, and a few examples are given in Table 11. Starting with a small phase shift, the system will approach the symmetric oscillations, Le., with zero phase difference, at a rate that increases with the coupling (a rough "half-life" is estimated to be 8000/k, s). The period and amplitude of the symmetric oscillations remain unchanged at all coupling rates. On the other hand, the unsymmetric oscillations change drastically in period, amplitude and general shape as the coupling increases. Table I11 shows the change in amplitude and period as the unsymmetric oscillations approach the Hopf bifurcation at k, = 11.52978 X 10-4 s-'. Figure 5 shows a plot of log A vs log r, where A is the relative distance from the Hopf bifurcation, Le., l[kx- k,(Hopf)]l/k,(Hopf), and T is the oscillation period. Near the bifurcation, the plot of linear with a slope of 0.343, showing that the period increases sharply as the coupling approaches the bifurcation point, beyond which the system goes over to an unsymmetric stable steady state. This sharp increase in period and the difficulty of determining the type of the Hopf bifurcation may point to the existence of another bifurcation, possibly a global homoclinic bifurcation (saddle very near the Hopf, at least to within -0.1% of it. The amplitude (Table 111) increases from the same value as the symmetric oscillations, goes through a maximum at roughly k, = s-l, and then decreases slowly as shown in Figure 6, before the oscillations collapse to the steady state. The shape of the oscillations changes, too. At low coupling values, the oscillations resemble the symmetric or the single cell ones as shown in Figure 2, but from k , = 2 X IO4 s-l, the system (30) Guckenheimer, J.; Holmes, P. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields; Springer-Verlag: New York, 1983; p 319.
4.0
I
-0..
. \
*.'
!
-3.60
1
I
I
I
-2.78
-1.96
-1.14
-032
1"O
0.50
log A
Figure 5. Relative distance from the Hopf bifurcation point vs the period (log-log scale). A, unsymmetric oscillations,k, = 5 X IO4 s-I, change [Br-I,, till Hopf at 28.00047 ~LM; 0, symmetric oscillations, k, = 5 X lo4 s-I, change [Br-1, till Hopf at 28.00047 pM; a, unsymmetric oscillations, standard conditions, change coupling till Hopf at 11.529 78 X IO4 s-'; 0,unsymmetric oscillations, k, = 5 X IO4 s-', change [MA], till Hopf at 28.751 04 mM; 0 , symmetric oscillations, k , = 5 X IO4 s-', change [MA],, till Hopf at 25.758 12 mM. 1-4 . 2
u
2.2-
a
2 2.0-I
n
3 1.8W
1 1.6F a
$ 1.4-
I k,
x
104/S-'
Figure 6. Relative amplitude (to amplitude of symmetric oscillations) vs coupling rate. Standard conditions.
spends more and more time near points B and D (see Figure l), which are very near the respective steady states. Thus one cell spends a lot of time near one steady state, while the other cell spends this time near the other one, and then in a rather short time (-40 s), they exchange their concentrations as shown in Figure 7a. At close inspection, the cells make small but increasing amplitude oscillations in a period that corresponds to the imaginary part of the approaching Hopf bifurcation (see Figure 7a, insert). It should be noted that the unsymmetric oscillations exist even below the lower turning point of Figure 3 (Table I); thus their existence is independent of the steady states described above. s-I and k , = 1.46196 s-I, the Between k , = 11.529 7 8 X system can exist both in a stable unsymmetric steady state and in symmetric (single-cell-like) oscillations. As the coupling in-
Bar-Eli
2372 The Journal of Physical Chemistry, Vol. 94, No. 6, 1990 I /
”1
26 288 0830
0838 k
x
4
0846
lO4/8
P=O
I
I/
* 0,
0
b)
i
-48-
-\
-64r’\u1p\pI)
r
v
301
j
10
-
L-
1
10-2
1
10-1
100
-
k,/S-’ -8
II
0
- 9 6i
2000
’
4000‘$000
‘8ooo
TIME/S
Figure 7. (a) Unsymmetric oscillations at k , = 11.368 X IO4 sK1 (near the Hopf bifurcation at k , = 11.52978 X IO4 s-‘) and standard conditions Note the long times spent near the steady states, the comparatively short exchange times, and the small oscillations (with period, 12 s, appropriate to the imaginary part of the Hopf eigenvalue, 0.5268) just before the exchange. (Insert: magnified X 100.) (b) Unsymmetric oscillations at [Br-1, = 27.9 pM and k , = 5 X lo4 s-I (near the bifurcation at 28.00047 pM). Note that the system is mostly near the middle branch of the symmetric steady state. Short excursions go to the upper branch of the unsymmetric steady state (Figure 9).
creases and approaches the upper Hopf bifurcation, one needs smaller and smaller increments in order to remain on the unsymmetric stable steady state and to not go over to the symmetric oscillations, showing that the attraction basin of the stable branch decreases, while that of the symmetric oscillation increases. The approach to the steady state goes through small decreasing oscillations with frequency matching the imaginary part of the eigenvalue (Figure 7a). Beyond the upper Hopf, only the symmetric oscillations exist, and no other type of oscillations is observed. Changes with Malonic Acid Steady States. As the concentration of malonic acid decreases, the bifurcation points of the symmetric solution (Figure 3 and Table I) approach each other and coalesce at [MA], = 26.51250 mM and k, = 0.220 168 SKI. At this point p = 4 2 directly, without going through a p = 3 region. The unsymmetric steady state solution is “pinched” off from the symmetric one, and the “mushroomn of the unsymmetric solution seen in Figure 3 becomes an i ~ o l a , disconnected *~ from symmetric solution. At [MA], = 25.758 12 mM, the symmetric solution becomes stable with p = 0, at all coupling rates. At this point the single cell becomes stable, too.I6 The existence of the isola of the unsymmetric solution does not change its behavior. Thus the kola as well as the mushroom has a rather wide range of stability wherep = 0. This range decreases with increase of [MA], and finally disappears, as the upper and lower Hopf bifurcations of the unsymmetric solution (Figure 3) coalesce at [MA], = 67.0 mM, as shown in Figure 8. Inside the “plot” is the range of stability of the unsymmetric solution p = 0, outside p > 2, and it is unstable. The points of the creation of the isola (at [MA], = 26.51250 mM) and the “kink” in Figure 8 (at [MA], = 26.2929 mM) do not affect the system’s behavior. Symmetric stable steady states can therefore coexist with unsymmetric ones below [MA], = 25.758 12 mM in the range of coupling rates “inside” the plot of Figure 8. Thus at [MA], = 10 mM, starting with the unsymmetric unstable steady state near the Hopf (k, = 1.726 246 X IO4 s-]), and putting k, = 1.8 X IO4 s-I (JI = 0), the system ends in the unsymmetric stable steady state. Decreasing the k , to 1.72 X I 0-4 s-l 0, = 2), the system ends at
-
Figure 8. Hopf bifurcation of p = 0 2 of the unsymmetric solution as a function of malonic acid concentration and coupling rate. The system is stable, p = 0, “inside”and unstable “outside”. Insert: magnification of the “kink” on the left-hand side, at [MA] = 26.2929 mM and k , = 0.8404 SKIwhere p = 0 4 directly.
-
the symmetric steady state. The approach to the steady state, whether symmetric or not, is slower the nearer the system is to the Hopf bifurcation and goes through small oscillations of increasing amplitude with period appropriate to the imaginary part of the Hopf. In no case were unsymmetric oscillations obtained “outside” the plot of Figure 8 (p > 0 region below [MA], = 25.781 20 mM), and thus symmetric steady state and unsymmetric oscillations do not coexist. Changes with Malonic Acid: Oscillations. When the coupling is kept constant and the concentration of malonic acid is decreased, the unsymmetric oscillations will change their period and amplitude in a fashion similar to that described above. Thus the period will increase as shown in Figure 6, where it is plotted as a function of the relative distances from the Hopf bifurcation, Le., [MA], = 28.751 04 mM at k, = 5 X s-l. The final slope of the graph is seen to be the same as in the previous case, namely, 0.343. In both cases, the same line of Figure 8 is crossed; in the first case, it is crossed from left to right, keeping malonic acid concentration constant and changing the coupling. In the second case, it is crossed from top to bottom, keeping the coupling rate constant and decreasing malonic acid concentration. Near the bifurcation, the system spends more and more time near the approaching unsymmetric steady state, with a short excursion to change positions as described earlier (Figure 7a). Similar behavior is observed for the symmetric oscillationswhich exist at all coupling rates above the Hopf bifurcation at [MA], = 25.758 12 mM. Again log T vs log A final slope is 0.343. While the period of the symmetric oscillations increases with the decrease of MA, the amplitude decreases from 7940 at [MA], = 100 mM to 4570 at 25.78 mM. Similar results were obtained for a single
In both cases, before leaving the vicinity of the appropriate steady state, the system makes small, increasing amplitude oscillations with appropriate period to the imaginary part of the approaching Hopf. The symmetric oscillations will coexist with either unsymmetric oscillations or unsymmetric, broken-symmetry, stable steady state, depending whether the coupling is “outside” or “inside” the plot of Figure 8. Above [MA], = 67.0 mM, no stable state exists, but the existence of two kinds of oscillations continues. As the coupling increases, the symmetric solutions, namely, oscillations at high malonic acid concentrations and symmetric steady states at low ones, become more dominant and finally take over, as seen from Table IV: increasing the coupling from k , = 65 X IO4 to 70 X lo4 S-I takes the system from the unsymmetric solutions to the symmetric ones. Table IV shows also that, as the coupling increases, the unsymmetric oscillations split to large and small ones. The amplitude
The Journal of Physical Chemistry, Vol. 94, No. 6, 1990 2373
Coupling of Identical Chemical Oscillators TABLE IV: Oscillations at Various Coupling Rates and Otherwise Standard Conditions with [MAb = 0.1 M 10~k,/s-~ 1 5 10 15 20 30 40 50 60 65 70 100
amplitude large, small
period/s 156.2 188.3 212.5 228 238 253 263.8 27 1.5 278.8 282.5 148.2 147.9
8910
IO 000, 1.03 11 200, 1.73 11 200, 3.00
11 890, 4.47 11 890, 6.13 11 220, 12.59 10900, 17.78 10290, 25.11 10290, 30.73 7 940 7 940
TABLE V: Unsymmetric and Symmetric Oscillations with Various [Br-I, at k, = 5 X lo4 s-' 1Os[ Br-],/M
type of osc" unsym L unsym LIS1 unsym LIS1 unsym LIS1 unsym LIS1 unsym LIS1 unsym LIS' unsym LIS' unsym LIS' unsym LIS1 sym L sym L
"The symbol L"Smgives the number of large (L) and small (S) am-
1 2 2.2 2.3 2.5 2.7 2.15 2.76 2.77 2.78 2.785 2.79 2.792
plitude oscillations per cycle.
1 I 10-5 IO-^ [Br-]o/M
10-3
Figure 9. Steady-state concentrations of bromide ion vs bromide ion concentration in the flow, at constant coupling rate of k, = 5 X IO4 s-l. Note the S-shaped curve of the symmetric solution (solid line) and the unsymmetric one (dashed line) bifurcating from it at [Br-] = 25.375 96 pM. The Hopf bifurcation on the symmetric branch p = 0 2 is at
-
-
of the large ones is constant at 10000, while that of the small ones increases about 30-fold before the symmetric oscillations take over. Changes with Bromide Ion: Steady States. The main features of the steady states as a function of bromide ion concentration at constant coupling rate ( k , = 5 X lo-" s-l) are shown in Figure 9. The figure shows an S-shaped symmetric steady state with Hopf bifurcation p = 2 0 at [Bile = 28.00047 pM on the upper branch. All other steady states, both symmetric and unsymmetric (which branch from the upper symmetric branch at [Br-1, = 25.375 96 pM),have p > 0 and are unstable. The S-shaped curve of the symmetric solution is identical, of course, with that of the single cell, and the small region of the three unstable steady states (region F in Figure 1 ) can be easily seen.21v22 All solutions above the Hopf bifurcation, regardless whether the initial point is on a symmetric or unsymmetric trajectory, fall to the symmetric steady state. Again, the system cannot support coexistence of stable symmetric steady state and unsymmetric oscillations. Change with Bromide Ion: Oscillations. The behavior of the oscillations of both kinds near the Hopf follows the same pattern as discussed above, namely, a sharp increase of period when the Hopf bifurcation is approached. The results are summarized in Table V and Figure 5 . For both types of oscillations, the final slope of log T vs log A plots is 0.713, contrasted to 0.343 in the previous cases. It is interesting to note that the unsymmetric oscillations collapse to the symmetric steady state at the same point as do the symmetric ones. Thus the end of the oscillations is probably not related directly to the existence of the Hopf, but just to the existence of
-
1050 1060 1080 1138 1238 1693 1925 2525 3270 4030 5375 6510
12 230 12230 I 3 330 I 4 130 16790 26 600 31 620 47 300 66 800 69 180 69 180 75 860
Symmetric Oscillations 0.5 1 2 2.5 2.6 2.7 2.75 2.77 2.78 2.79 2.793 2.795
[Br-] = 28.00047 pM.
period/s amplitude Unsymmetric Oscillations 973 IO 000
424 439 478 53 1 560 845 1352 1930 2407 4330 5830 7820
4 470 4 570 5 500 6 680 7 500 15 800 47 300 70 790 75 860 75 860 1 5 860 75 860
an accessible stable state, and the existence of another bifurcation nearby is indicated. From Table V we see that the ratio of the amplitudes (unsymmetric to symmetric) decreases as the bifurcation is approached, although the shape of the oscillations is completely different as can be seen from Figure 7b, which shows the oscillations at k, = 5 X lo-" S-I and [Br-1, = 27.9 X lod M, Le., near the bifurcation. It is seen that the system spends most of its time near the upper (near the bifurcation) branch of the symmetric steady state. It then jumps for a short while to the upper branch of the unsymmetric state, returns to the symmetric state, then makes a large excursion before coming back to the symmetric state again. These oscillations are seen to be of a completely different nature than those shown earlier (Figure 7a), and thus it is no wonder that the approach to the steady states is different-Le., different slopes of Figure 5 . It should be noted that the symmetric oscillations for other values of malonic acid show exactly the same behavior as the single cell, Le., bursting oscillations of various types as described earlier.=
Discussion At first sight, the coupling between two identical CSTRs should not change anything in their behavior, and the combined system should behave exactly as the single cell does. These symmetric solutions do indeed exist and have the same stability. In other words, whenever the initial conditions are the same in both cells, the symmetric solution, whether steady state or oscillations, prevails. The many references cited in the Introduction corroborate, both experimentally and theoretically, that unequal concentrations, Le., "broken symmetry"' can be sustained in the two cells notwithstanding the equalizing effect of the diffusionlike coupling. The work of Bar-Eli and Geiseler3' showed that if the single CSTR system is bistable and can exist in either of two stable steady states (under the same set of constraints; region A of Figure l ) , then the coupled system can be either in states that are not too far from the original ones or, if the coupling is strong enough, in one state only. The first case is an unsymmetric, brokensymmetry solution, while the second one is a symmetric solution. These results were not surprising, since the original single CSTRs were not in the same steady states. (31) Bar-Eli, K.; Geiseler, W. J . Phys. Chem. 1981, 85, 3461. (32) Mankin, J. C.; Hudson, J. L. Chem. Eng. Sci. 1986, 41, 2651
J . Phys. Chem. 1990, 94, 2314-2311
2374
TABLE VI: Regions of Existence of Symmetric and Unsymmetric Oscillations and Stable Steady States“ coupling
low medium
high
low [Br-] high [Br-] low [MA] high [MA] low [MA] -high [MA] sym SS only sym osc + sym SS only unsym osc sym SS sym osc + sym SS only unsym SS unsym SS sym SS only sym osc only sym SS only
+
‘The exact limits of the various regions are given in the text. This table is only qualitative. Unlike the above case, in this report we deal with situations where the single CSTR has either one single steady state, stable (region D in Figure 1) or unstable (region E), or three unstable steady states (region F). The single CSTR will, therefore, oscillate or be in a single stable state. The report shows that, even in this situation, the obvious symmetric solution is not the only one, and an unexpected unsymmetric solution, oscillating as well as stable, can coexist with the symmetric ones. In Table VI the various regions of existence of the symmetric and unsymmetric solutions are qualitatively shown. It is seen that the system containing two identical CSTRs can, depending on constraints and initial conditions, be in the following states: (a) stable symmetric steady state, (b) symmetric oscillations, (c) coexistence of symmetric and unsymmetric steady states, (d) coexistence of symmetric oscillations and unsymmetric stable steady state (broken symmetry), and (e) coexistence of symmetric and unsymmetric oscillations. The latter differ from the former in phase, in amplitude, and in period. On the other hand, no unsymmetric oscillations were found to coexist with the symmetric steady state. All the initial conditions tried (even those with p = 2 ) ended in either of the two stable
states. While this paper was being written, the work of Crowley and Epstein” appeared. In that work two similar, but not identical, BZ oscillators were coupled and the results compared to calculations of two coupled Oregonators (without flow). Their experimental results corroborate the predictions given in this paper, in particular the coexistence of the symmetric and either unsymmetric oscillations or steady state and the lengthening of the period of the unsymmetric oscillations toward their “end” when they collapse to the unsymmetric steady state. Another feature of the oscillations, whether symmetric or unsymmetric, is that they “end” at a Hopf bifurcation, with a period increase obeying a power law with the distance to the Hopf. In all cases the final change of period T with the relative distance with different from the Hopf A obeys the relationship T = values for a. The same relationship between A and r holds for a single with similar a’s: a small a when [MA], is decreased and a large one when [Br-1, is increased. The particular exponent depends on the position of the bifurcation but not on the symmetry of the oscillation. The ‘‘Hopf‘ behavior, namely, small oscillations with period depending on the imaginary part of the vanishing real part eigenvalue, is observed only temporarily. This behavior is rather unexpected and indicates that something else may be hidden in the dynamics of the system. As has been mentioned earlier, we were unable to determine the nature of the Hopf bifurcations (super- or subcritical) nor were we able to continue beyond it by AUTO. In all cases the computation approached the Hopf to within -0.1%; it is therefore quite possible that the oscillations “end” not at the Hopf but at another such as saddle-loop, homoclinic, global b i f ~ r c a t i o n which , ~ ~ is located very near the Hopf. The possibility that the period does not go infinite, but just to a very high value near a supercritical Hopf bifurcation, cannot be excluded. It is quite possible that in the high dimensionality of the problem ( 2 X 5 = 10 dimensions) rather unexpected behaviors can be manifested.
Kinetics and Mechanisms of the C02 Laser Induced Decompositions of CFCI, and CF2CI2 R. N. Zitter,? D. F. Koster,**tT. K. Choudhury,f§ and A. Cantonitsl Department of Physics and Department of Chemistry and Biochemistry, Southern Illinois University, Carbondale, Illinois 62901 (Received: April 17, 1989; In Final Form: August 22, 1989)
The kinetics and reaction mechanisms of the decomposition of CFCI, and CF2CI2induced by a continuous C 0 2 laser are examined and compared with previous results for CF3CI. All three compounds clearly exhibit two temperature regimes, with similar reaction mechanisms in the low-temperature regimes and significantlydifferent mechanisms in the high-temperature regimes.
Introduction Fluorine-containing compounds have often been used for studies of CO, laser induced gas-phase chemical reactions because they typically have strongly absorbing infrared band frequencies matching those of CO, lasers. However, few of the reported investigations have been devoted to careful examination of decomposition kinetics and mechanisms, while standard pyrolytic techniques do not span the temperature range available with lasers. Furthermore, laser excitation can induce reactions far from gas cell walls to avoid the complications of wall reactions with fluorine and highly reactive radicals. ‘Department of Physics. *Department of Chemistry and Biochemistry. Present address: Department of Chemistry, The Catholic University of America, Washington, DC 20064. Present address: California Laboratories, Carlsbad, CA 92008.
We have reported’ a detailed study of the decomposition kinetics and mechanisms of CF3CI induced by a continuous (CW) C 0 2 laser. Similar studies of CFCI, and CF2Clzare presented here and form an interesting basis for comparison. For all three compounds, the data show two distinct temperature regimes, with reaction pathways that are dominant in one regime but not the other. The reaction mechanisms are similar in the “low”-temperature regimes but differ significantly at higher temperatures.
Experimental Details The techniques employed here are the same as those described elsewhere.’-3 Briefly, a laser beam of diameter 1.6 mm is absorbed ( I ) Zitter, R. N.; Koster, D. F.; Choudhury, T. K . J. Phys. Chem. 1989. 93, 3167. ( 2 ) Zitter, R. N.; Koster, D. F.; Cantoni, A,; Pleil, J. Chem. Phys. 1980, 46. 107.
OQ22-3654/9Q/2094-2314~02.50/0 0 1990 American Chemical Society
