J . Phys. Chem. 1991, 95, 8706-8713
8706
bimolecular reactions with several possible product channels. A comparison with the rates of the isovalent BH species indicates that the AIH reactions are faster by a factor of 2-3. This may be due to the larger size of AIH and to its larger dipole moment. Since both O2and H 2 0 are present in the MOVPE environment, there may be interest in determining which reaction products may contribute to the contamination of AlGaAs films. The association reaction of AIH with D2 is very slow and occurs at the limits of our detection but offers a rough comparison, when paired with BH D2 results, to the group IVA isovalent reaction pair of C H D2 and SiH D2. In the case of C H and SiH, the reaction rate with D2 decreases by a factor of 5000. We observe the same trend in the group IIIA pair, BH and AIH. Begemann et al.36have SiH(v”=l) data that suggest that the difference in rates between C H and SiH is due to an increased barrier in the
+
+
+
association reaction. It is possible that the same effect occurs in the group IIIA pair. The reactions of AIH with saturated hydrocarbons are slower than our detection limit. In contrast, the reactions with unsaturated hydrocarbons are very rapid. In the case of C2H2,evidence from methyl-substituted analogues indicate that the AIH attacks the a bond. In the case of AIH + C2H4,the rate is slower with methyl substitution, which indicates a mechanism other than a addition may be present. Acknowledgment. We thank the Office of Naval Research for the funding of this research. Registry NO. AIH, 13967-22-1; H2, 1333-74-0 D2,7782-39-0; C3H8, 74-98-6; 02,7782-44-7; H20.7732-I 8-5; C2H2, 74-86-2; C H 3 m C H 3 , 503-17-3; CZH4, 74-85-1; (CH,),C=CHCH,, 51 3-35-9.
Canard Explosion and Excitation in a Model of the Belousov-Zhabotinsky Reaction Morten Brans* Mathematical Institute, The Technical University of Denmark, DK-2800 Lyngby, Denmark
and Kedma Bar-Eli Sackler Faculty of Exact Sciences, School of Chemistry, Tel-Aviv University, Ramat Aviv 69978, Israel (Received: February 5, 1991)
The experimental evidence for sudden formation of relaxation oscillations (hard transitions),or the sudden change of amplitude and period of oscillations,is reviewed. It is shown that, in addition to the well-known mechanisms for hard transitions, there is another way in which a very fast transition from small to large oscillations can occur. This mechanism termed canard explosion is analyzed in terms of a well-known chemical model, the two-uariable Oregonator. The theory of the canard explosion, which occurs close to a Hopf bifurcation, is analyzed and compared with computational results. The resulting difficulties of differentiating experimentally among the various transition mechanisms are discussed. The well-known phenomenon of excitation, i.e., a large deviation of a system (chemical, neural, or other) after perturbation from a stable steady state, is shown to be. closely related to the canard explosion. The common mathematical properties of the system, on which both phenomena depend, are discussed.
Introduction In many cases, chemical oscillations either start abruptly and go directly to large relaxation oscillations or change suddenly, within very small control parameter range, from small-amplitude oscillations to large ones. A complete citation of all the examples in the literature is impossible, and we shall give a few examples only. Thus in the BZ (Belousov-Zhabotinsky) reaction, using Mn2+ as catalyst and malic acid or a mixture of citric and malonic acids as organic substrates, Maselko’v2 has found, depending on experimental conditions, both small, increasing amplitude, oscillations and direct change of steady state to large relaxation oscillations (with and without hysteresis). Also in the BZ reaction conducted in a CSTR (continuous stirred tank reactor), Schmitz et al.’.‘ have found, as the flow rate changes, transitions from small sinusoidal oscillations to large relaxation ones. Similar results were obtained also by Hudson et al.5 In a system containing chlorite, bromate, and iodide ions in sulfuric acid, Alamgir and Epsteid have found a direct transition from steady state to large relaxation oscillations without hysteresis (1) Maselko, J. Chem. Phys. 1982, 67, 17.
(2) Maselko, J. Chem. Phys. 1983, 78,381.
(3) Schmitz, R.A.; Graziani, K. R.;Hudson, J. L. J . Chem. Phys. 1977, 67, 3040. (4) Graziani, K. R.;Hudson, J. L.;Schmitz, R. A. Chem. Eng. J . 1976, 12, 9. ( 5 ) Hudson, J. L.;Hart, M.; Marinko, D. J . Chem. Phys. 1979, 71, 1601. (6) Alamgir, M.: Epstein, 1. R. J. Am. Chem. SOC.1983, 105, 2500.
0022-3654/91/2095-8706$02.50/0
when reversing the system’s constraints. In a similar system, with bromide’ instead of iodide ion, the system remained excitable after the large oscillations disappeared as the flow rate was changed. Similar transitions occur in the permanganate hydrogen peroxide oscillator* and with the bromate sulfite ferrocyanide os~illator.~ In gas-phase reactions such as hydrogen or acetaldehyde oxidations,’*I2 temperature, pressure, and light emission pulses may appear directly and without hysteresis as a system’s constraint, such as ambient temperature, is changed. These phenomena occur also in closed systems in which the constraints are changing slowly in time. Thus Smoes” has found a sudden increase in amplitude and period in BZ reaction as sulfuric acid concentration is changed. Similarly, Burger and KorosI4 have found the BZ relaxation oscillations to appear as soon as bromomalonic acid reaches a critical concentration. (7) Orban, M.; Eptein, I. R. J . Phys. Chem. 1983, 87,3212. (8) Nagy, A.; Treindl, L. J . Phys. Chem. 1989, 93, 2807. (9) Edblom, E. C.; Yin Luo; Orban, M.; Kustin, K.; Eptein, 1. R.J . Phys. Chem. 1989, 93, 2722. (10) Gray, P.;Scott, S. K. In Oscillations and Traveling Waves in Chem-
ical Systems; Field, R.J., Burger, M., Eds.; John Wiley and Sons: New York, 1985; p 493. ( I 1) Gray, P.; Griffiths, J. F.; Hasko, S. M.; Lignola, P. G . Proc. R. Soc. London, Ser. A 1981, 374, 313. (12) Griffiths, J. F. In Oscillations and Traveling Waves in Chemical Sysfems; Field, R.J., Burger, M., Eds.; John Wiley and Sons: New York, 1985; p 529. ( I 3) Smoes, M. L. J . Phys. Chem. 1979, 71, 4669. (14) Burger, M.; Koros, E. J . Phys. Chem. 1980, 84, 496.
0 199 1 American Chemical Society
Canard Explosion and Excitation in the BZ Reaction
The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 8707
Air oxidation of benzaldehyde1s17 catalyzed by Co2+and Brshows a hard transition directly from steady state to relaxation oscillations. The possible transitions from stable steady state to oscillations are discussed by several One interpretation of direct transition from a steady state to large relaxation oscillations (hard transition) without hysteresis is that a limit cycle "collides" with a saddle node, forming a homoclinic orbit. Another interpretation may be that a subcritical Hopf bifurcation2' exists, but the hysteresis between the limit cycle and the steady state occurs within a very small region of constraints and thus cannot be observed experimentally. A sudden increase in amplitude and period can also occur, however, in the case of a supercritical Hopf bifurcation. In this case a small stable limit cycle emerges from the point of bifurcation, and thus one would usually expect to find slow changes in amplitude and period with no hysteresis between the steady state and oscillations. In some cases the period and amplitude can increase very quickly within a very small range of control parameter, and thus large oscillations can be suddenly created. This phenomenon is called a canard explosiofl and is well-known for the van der Pol equation. [The term 'canard" (French for 'duck") comes from the shape of the limit cycle at the transition in the case of van der Polzzcsystem.] It is thus possible for a small stable (unstable) limit cycle created at a super- (sub-) critical Hopf bifurcation to go, within a very small range of parameter, to large relaxation oscillations. A few reports appeared recently in the literature showing that canard explosion does indeed occur in systems of chemical interest. Thus, Gaspar and S h ~ w a l t e claim r ~ ~ a canard explosion in amplitude and period in a model for iodate oxidation of sulfite and ferrocyanide. Scott and T ~ m l i nas~well ~ as Gray and Thuraisingham" have reported the existence of canards in the Autocatalator*5 model. The Brusselator model has been discussed earlier from a different point of view by Baer and Erneux.22f Another related phenomenon is called e x c i t a t i ~ n . ~ *This ~~-~* is related to a system which is perturbed from its singly existing, stable steady state. While small perturbations will cause the system to return immediately to the steady state, slightly larger perturbations, beyond some well-defined threshold, will cause the system to revert to the steady state only after making a very large amplitude excursion, following quite closely the 'previously" existing limit cycle. Excitations have been treated theoretical-
-r-
/
d -
N
Figure 1. Schematic plot of the slow manifold S ynullin a w vs y coordinate system. Trajectories parallel to the y-axis when t = 0 are also shown. 1y224f-24*26.27*49 and e ~ p e r i m e n t a l l yin~ chemical ~~ systems and also
in neural systems by Rinzel et aLZ8 Also, FitzHughMhas treated this phenomenon in a model of nerve excitation (similar to the van der Pol oscillator). We shall refer to this work below. The excitability property of the medium, be it a chemical system, a neural one, or a cardiac muscle, is extremely important to the propagation of wavefronts in such media.37 This "excitation" should not be confused with the transition of a system from one stable steady state to another or from one stable steady state to a coexisting limit cycle (as occurs in a subcritical Hopf bifurcation) under the influence of sufficiently large perturbation^^^*^^ or by noise.47 The latter phenomenon requires the coexistence of two stable steady states or a steady state and a stable limit cycle, while we are concerned with systems that have only one stable steady state. The excitation occurs, therefore, before the Hopf bifurcation, where the steady state is globally stable. The purpose of this paper is thus 3-fold: (a) to give a simple mathematical model of the canard explosion; (b) to compute and show its behavior in the case of the two-dimensional Oregonator, while pointing out the experimental difficulties in observing this phenomenon; and (c) to show the very close relationship between the two phenomena canard explosion and excitation. Theory In this section, the theory will be sketched without going into the mathematical details; these will be given more fully in the Appendix. Further details can be found in the literature.22 We consider systems of the general form j = w - F(y,w)
(IS) Jensen, J. H. J . Am. Chem. Soc. 1983, 105, 2639. (16) Jensen, J. H.; Roelofs, M. G.; Wasserman, E. In Non Equilibrium Dynamics in Chemical Systems; Vidal, C., Pacault, A., Eds.; Springer Verlag: Berlin, 1984; p 35. (17) Roelofs, M. G.; Wasserman, E.; Jensen, J. H.; Nader, A. E. J. Am. Chem. Soc. 1983,105,6329. (18) Sheintuch, M.; Luss, D. In Non Equilibrium Dynamics in Chemical Systems; Vidal, C., Pacault, A., Eds.; Springer Verlag: Berlin, 1984; p 33. (19) Boissonade, J.; De Kepper, P. J . Chem. Phys. 1980,84, 501. (20) Bar-Eli, K.; Brans, M. J . Phys. Chem. 1990. 94, 7170. (21) (a) Hassard, B. D.; Kazarinoff. N. D.; Wan, V. H. Theory and Applications of Hopf Bifurcation; London Mathematical Society Lecture Notes Series No. 41; Cambridge University Press: Cambridge, 1981. (b) Marsden. J. E.; McCracken, M. The Hopf Bifurcation and its Applications; Springer Verlag: New York, 1976. (22) (a) Benoit. E.; Callot, J. L.; Diener, F.; Diener, M. Publication de PInsrirur de Recherche Marhematique Auancee; Strasbourg, 1980; Collect. Marh. 1981,32, 37. (b) Eckhaus, W. In Asymptotic Analysis II; Springer Lecture Notes in Mathematics No. 985; Verhulst, F., Ed.; Springer Verlag: New York, 1983. (c) Kaas Petemn, C.; Brans, M. In Proceedings of rhe Ilrh
IMACS World Congress on System Simularion and Scientific Computation; Wahlstrom, B., Ed.; North-Holland: Amsterdam, 1985; Vol. 4, p 285. (d) Brans, M. Math. Eng. Ind. 1988, 2, 51. (e) Brans, M. Mat-report 1989-22, Mathematical Institute, The Technical University of Denmark. (f) Baer, S. M.; Emeux. T. SIAM J . Appl. Math. 1986, 46, 721. (23) (a) Noszticzius, 2.;Wittmann, M.; Stirling, P. J . Chem. Phys. 1987, 86, 1922. (b) Finkcova, J.; Dolnik, M.; Hrudka, B.; Marek, M. J . Phys. Chem. 1990. 9 4 4 1 IO. (24) Ruoff, P.; Noyes, R. M. J . Chem. Phys. 1986, 84, 1413. (25) Ruoff, P. Phys. Lctr. 1982, 90, 76. (26) Field, R. J.; Noyes, R. M. Faraday Symp. Chem. SOC.1974, 9, 21. (27) Bar-Eli, K.; Noyes, R. M. J . Chem. Phys. 1987.86, 1927. (28) Rinzel, J.; Ermentrout, G. B. In Methods in Neural Modeling, Koch, C., Segev, I., Eds.; MIT Press: Cambridge, MA, 1989; p 136.
* = tGb,w,c,t)
(1)
(2)
where c and t are parameters. A main assumption is that the S = ynullisocline, defined by w - F(y,w) = 0, is a curve with a maximum M, located at (yo,wo),and minimum N, located at (yN, w N ) ,as shown in Figure 1. It is not very important that F is independent of the parameters, but it makes things a little simpler and suffices for our purpose. When t = 0, the ynullisocline consists entirely of steady states. The phase portrait is easily found and is depicted in Figure 1. One sees that the equilibria on the parts where the slope of S is positive (left and right branches of S) are stable, while they are unstable when the slope is negative (middle branch of S). From every initial point away from S, the trajectory will be a horizontal line going to (and away from) the portions of positive (negative) slopes of S. The ynullisocline is denoted the slow manifold S, also when t # 0. When t # 0, the steady state is located at the intersection between the isoclines S = ynulland wndl(found by solving G = 0). We assume that the steady state is at M when c = 0 and that there are no other steady states, i.e., no other intersections between yndl and wnUII.This assures that when t and c are small, there exists also exactly one steady state, located close to M . Some specific assumptions given in the Appendix2*' take care of this. In par(29) Ruoff, P.; Noyes, R. M. J . Phys. Chem. 1985,89, 1339. (30) Bar-Eli, K.; Geiseler, W. J . Phys. Chem. 1983, 87, 1352.
8708 The Journal of Physical Chemistry, Vol. 95, No. 22, 1991
N
I
N
Figure 2. Schematic curves in a w vs y coordinate system: S ynull, stable (Ms) and unstable (Mu) manifolds, flows, and typical trajectory (dashed line) for z # 0. (A) A < 0, Le., Msbelow Mu.(B) A > 0, Le., Ms above Mu.
ticular, we assume that the steady state moves from one side of M to the other as c passes through 0. Computing the Jacobian matrix and using a simple implicit function theorem argument, it can be shown22cthat a Hopf bifurcation occurs at a value of c close to 0. When e = 0, the slow manifold is inuariant. That is, for all initial conditions on S, the trajectory stays on S for all times. In fact, every point on S is a steady state when e = 0. In the general theory of differential equations, a number of results on the persistence of invariant manifolds are available: i.e., theorems on the existence of invariant manifolds near S even when e # 0. Thus, when e # 0, one can find on the interval (-m, yo] exactly one trajectory which stays close to S as t --. Furthermore, this trajectory is attracting. This trajectory is denoted Ms, the stable manifold, and its separation from S is of the order e. It is realized, of course, that, our system being of chemical nature, the (--, yo] interval means only all the positive values of y below yo, Le., 0 < y < yo. A completely analogous result holds on bN, m), Le, the right branch of S,but this is not very important for the present work. In a similar way one finds on the interval bo,Y N ] a trajectory, denoted Mu, with the following properties; it is the only trajectory that stays close (of the order e) to S as t increases, and it is repelling, just as the slow manifold is on this interval, when c = 0. This trajectory will be called the unstable manifold. Thus in Figure 2, two special trajectories are shown (in heavy lines): Ms, located near the left portion of S and toward which the trajectories are attracted, and Mu, passing through N, being near the middle part of S,and from which the trajectories are repelled. Both Ms and Mu continue beyond the maximum at M , but they do not have there their above-mentioned properties; namely, they are no more close to S and do not attract or repel, respectively. This is due to the folding of S at M. Phrased differently, it is not possible to find a global invariant manifold for the system. This is a persistence result on invariant manifolds: Note that both the existence and the stability properties that applied at e = 0 are preserved. However, the very special structure of the manifold is not preserved. At e = 0, it consists entirely of steady
-
Brans and Bar-Eli states-and thus of infinitely many trajectories-while it is simply one trajectory when e > 0. The term manifold may seem a little overwhelming for such a simple thing as a trajectory. The notation is inherited from higher dimensions, and thus we use this term here too. Now we have established the existence of two special trajectories, Msand Mu.Their separation a t point yo, defined as A = M s b 0 ) - MuCyo), has a significant influence on the dynamics. Let us look at a typical trajectory starting somewhere to the left of yo. It will be attracted to Msuntil it reaches y = yo; then it will come under the influence of Mu,from which it will be repelled. Now, if A < 0, Le., Mu is above Ms,the trajectory will be quickly pushed down and to the left and will come back to the vicinity of M sas the dashed line in Figure 2A. If, on the other hand, A > 0, Le., Mu is below Ms,the trajectory when reaching yo will go above Mu, as the dashed line in Figure 2B, go to the vicinity of the right branch of S,and move down near it until N, from which it will move quickly to the left branch of S, near Ms. In other words, A < 0 results in small trajectories, while A > 0 results in large trajectories. The border between the two types of trajectories will be when A = 0. At the point A = 0, the border between small and large responses, the two trajectories Msand Mu become one and the same, the canard manifold. Per definition, the value of c where this takes place is the canard point. Stated differently, the canard point is the value of c where it is possible to make a global reduction from two to one dimensions. There is a simple connection between the variation of A with c and the type of Hopf If the bifurcation is supercritical, i.e., the limit cycle is stable and coexists with an unstable steady state, A is negative at the bifurcation and changes sign to a positive value at a point where the steady state is unstable. If, on the other hand, the bifurcation is subcritical, A is positive at the bifurcation and changes sign to a negative value at a point where the steady state is stable. In the supercritical case, the attractor being a limit cycle, the border, as A passes through zero by changing the control parameter, will be between a small limit cycle and a large one. The small limit cycle will become suddenly a large one. This is termed canard explosion. In the chemical systems of interest, all the trajectories must be bounded. Thus, in the subcritical case, the small unstable limit cycle born at the Hopf bifurcation must be surrounded by a stable limit cycle. At the canard point, the small (unstable) limit cycle will experience an explosion and will merge with the large (stable) limit cycle in a looploop'8,20~22f bifurcation shortly thereafter. In this case the system will have a bistable region with both a stable steady state and a stable limit cycle. In a situation where the steady state is a global attractor, the positions of Msand Mu must be as depicted in Figure 2A. Then M u will serve as a separatrix between what should be regarded as small and large perturbations. In other words, a perturbation below M U will cause the system to revert quickly to the steady state (the dashed line in Figure 2A), while a perturbation above it will cause the system to revert to the steady state only after a large excursion near the right-hand side of S. M u thus serves as the excitation limit described above. Both the sudden increase in period and amplitude (canard explosion) and the excitation are related to one cause, namely, the behavior of the various manifolds in the system. Model and Computations The model is a t ~ o - v a r i a b l e ~ ~ *version ~ ' * ~ ' -of~ ~the Oregonator (31) Gaspar, V.; Showalter, K . J . Chem. Phys. 1988,88. 778. (32) Gaspar, V.; Peng, B.;Showalter, K. In Spatial Inhomogeneities and Transient Behaviour in Chemical Kinetics; Nicolis, G., Gray, P., Scott, S.K., Eds.; Manchester University Press: Manchester, U.K. 1990; pp 277-288. (33) Field, R. J.; Noyes, R. M. J . Chem. Phys. 1974, 60, 1877. (34) Tyson, J. J. J . Chem. Phys. 1977,66, 905. (35) (a) Gear, C. W. Numerical Initial Value Problems in Ordinary D~JerentialEquations; Prentice-Hall: Engelwood Cliffs, NJ, 1971; pp 209-229. (b) Hindmarsh, A. C. Gear: Ordinary Differential Equation Sysfem Soluer; UCID 30001 Rev. 3; Lawrence Livermore Laboratory, University of California: Livermore, CA, 1974.
The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 8709
Canard Explosion and Excitation in the BZ Reaction three-dimensional scheme of the BZ reaction, namely: k k l A Y - kzXY k3AX- 2k4XZ
+
TABLE I: Values of Rate Constants (3)
with dimensions
klA
1000 2 x 109 lo00 5 x 10'
k2
k3A
In the usual interpretation, A 3 [ B a r ] , X = [ H B a 2 ] , Y E[Br-1, Z = 2[Ce4+], and f is a stoichiometric factor relating the regeneration of Y to the consumption of 2. Making the Feudo-steady-state approximation for X, namely, assuming that X = 0, we obtain a two-dimensional system with Y and 2 as variables. We note that other possible reductions of the three-dimensional system to a two-dimensional one are given in the l i t e r a t ~ r e . ~ ' , ~ ~ The values of the various rate constants of the BZ reactions (for discussion of the various possible sets of rate constants see refs 38 and 48) are such that the obtained bifurcation will be subcritical, resulting in an exploding unstable limit cycle27~39 (compare Figure 3 of ref 27), with a small hysteresis region between stable steady state and oscillations. It seems very instructive, however, to study the case of an exploding stable limit cycle which occurs in the case of a supercritical bifurcation. This explosion can, a t least in principle, be observed experimentally. Since the parameter range in which the explosion occurs may be very small, as we shall see below, the experiments, as well as the computations, may turn out to be extremely difficult. Moreover, because of the smallness of the explosion region and its proximity to the bifurcation, it may be difficult to ascertain whether the bifurcation is super- or subcritical. The experimental as well as the computational results of the two cases may turn out to be very similar. The experimentalists as well as the computational chemists should be aware of these facts. In order to investigate the canard phenomenon with stable limit cycles, we have deliberately used a large value of klA so that the bifurcation will become supercritical, the limit cycle stable, and the explosion behavior manifested more clearly. The values of the various rate constants are given in Table I. The approximation made, namely, X = 0, is still justifiable as can be seen by comparing the full three-dimensional system with the approximated two-dimensional one, shown in Figure 3 (dashed line). These plots are very similar to the ones obtained with the small value of k l A used earlier,27J132indicating that the relaxation oscillations do not vary with large changes of this constant. According to the procedure outlined in the Appendix, the concentrations are transformed to nondimensional quantities, resulting in a set of two differential equations, namely:
9
dy/dT
- y - x@)y = c(x@)- Z ) fz
i = dz/dT
(6) (7)
(8)
- w)
(9)
CCfXb)
P
1000
4
2.5 X 0.001
e 'e
1
1
":::I 'F
8.4
z
7.8
L , , , , , , , ,
6.80
0.8
1.6
2.4
3.2
Y
Figure 3. Dashed line: Large limit cycle of the full three-dimensional Oregonator at f = 0.7147164889 (the value of the Hopf bifurcation). Full line: Large limit cycle of the two-dimensional system at f = 0.713654372198,Le., Af= 8.002X lod (just after the canard explosion).
1
8.6752
8'66881
8.6656
0.616
0.632
0.648
0.664
0.68
Y
j=w-y-x@)y kJ=
k4 k3
8.6624 0.6
with x@) given in eq 17 in the Appendix. Substituting w = fz, we obtain
dimensionless
which is a Litnard type system of the kind considered above. The two systems, namely, (6)-(7) and (8)-(9), are, of course, equivalent. In order to conform with previous work on the subject, (36) Guckenheimer, J.; Holmes, P. Nonlinear Oscillations, Dynamical Systems. and &furcations of Vector Fields: Springer Verlag: New York, 1983. (37) (a) Jahnke. W.; Skaggs. W. E.;Winfree, A. T. J. Phys. Chem. 1989, 93, 740. (b) Winfree, A. T. SIAM Rev. 1990, 32, 1. (c) Winfree, A. T.; Jahnke, W. J . Phys. Chem. 1989,93,2823. (d) Tyson, J. J. J . Chem. Phys. 1987, 84. 1359. (e) Keener, J. P.;Tyson, J. J. Physica 0 1986, 210, 307. (38) (a) Tyson. J. J. In Oscillarions and Traveling Waves in Chemical Systems; Field, R. J., Burger, M.,Eds.; John Wiley and Sons: New York, 1985; p 93. (b) Tyson, J. J. J . Phys. Chem. 1982,86,3006. (39) (a) Troy, W. C. In Oscillations and Traveling Waves in Chemical
Systems; Field, R. J., Burger, M..Eds.; John Wiley and Sons: New York, 1985; p 145. (b) Rinzel, J.; Troy, W. C . J. Chem. Phys. 1982,76,1775. (c) Janz, R. D.;Vanecek, D. J.; Field, R. J. J. Chem. Phys. 1980,73, 3132.
Figure 4. Small limit cycle of the two-dimensional system at f = 0.713654370698,i.e., Af = 8.0005 X lod Cjust before the canard ex-
plosion). ( 0 )Steady state.
we have used for the computational work, the form given in eqs 6 and 7, and thus, the results will be given in terms of y and z. Integration of the differential equations (6) and (7) has been performed by the Gear methodg5using error limits as small as 1X and time steps of a t most 0.001 of the period. Hassard's2' BIFORZ program has been used to calculate the system behavior near the bifurcation. Results By use of the data in Table I, various significant points of the system are calculated. These are summarized in Table 11, in dimensioned, nondimensioned, and differences from point M(yBwo) values. In particular, a supercritical Hopf bifurcation occurs at fH = 0.7 13646370198 (another, subcritical bifurcation occurs at f = 1.85285455836, but this will not concern us here) withy,, =
8710 The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 TABLE 11: Coordinates of "Important" Points point f
Y
Brans and Bar-Eli
W
2
X
Nondimensional Data M
N Hopf SSS'
canardb canardC
0.713587011970 0.71 358701 1970 0.713646370198 0.7 I2646370198 0.71 3654326247 0.71 3654370948
0.639851634248 1.35806262404 0.639899358548 0.639095176261 0.639905755149 0.639905791088
8.67767486947 7.12256806335 8.67695306202 8.68912035351 8.67685631898 8.67685577544
6.19227608096 5.08257206 188 6.19227607 6.19227008014 6.1 9227605027 6.19227605023
8.67767486947 2.7425 1670867 8.67695306202 8.689 12035351 8.67685631898 8.67685577544
Data with Point Mas the Origin M
N Hopf SSS"
canardb canard?
0 0 0.593582280004E-04 -0.940641 77 19978-03 0.673142770005E-04 0.673589780007E-04
0 0.718210989788 0.477243003871E-04 4.756457986789843 0.541209008880E-04 0.541568400578E-04
0 -1.555 10680613 -0.721807455250E-03 0.1 14454840314E-01 -0.818550495069E-03 -0.819094042470E-03
0
0
-1.l09704Ol908 -0.2386695996398-07 -0.600081591529E-05 -0.3069214926658-07 -0.307329628413E-07
-5.93515816081 -0.72 1807455363843 0.1 14454840314E-031 -0.81 8550495239E-03 -0.819094042413E-03
0.309613804048E-02 0.254128603094E-02 0.309613802855E-02 0.309613504007E-02 0.309613802513E-02 0.30961380251 1E 4 2
0.433883743474845 0.356128403l67E-02 0.433847653101 E-05 0.434456017676E-05 0.433842815949E-05 0.433842788772E-05
Dimensioned Data M N
Hopf SSS"
canardb canardC
0.713587011970 0.713587011970 0.713646370198 0.71 2646370198 0.71 3654326247 0.713654370948
0.3 19925817l24E-06 0.679031312018E-06 0.319949679274E-06 0.319547588131 E-06 0.319952877574E-06 0.3 19952895544E-06
0.433883743474E-02 0.356128403167E-02 0.433847653101E-02 0.434456017676E-02 0.433842815949E-02 0.433842788772E-02
"Stable steady state located at if=/= =fH = -0.001. bCanard point calculated by eqs IO and 48. CCanardpoint calculated from the explosion point in Table 111. TABLE III: Periods and Amplitudes of Limit Cycles with f > fH As Determid by Direct Integration f AfX IOd period" amplitude y 0.713646670198 0.3000 1.00923 0.4648800E-02 1.01441 0.713646870198 0.5000 0.571 36OOE-02 1.04577 0.103392OE-01 0.713647370198 1 .0000 1.06296 0.1 21 1360E-01 0.71 3648370198 2.0000 1.08545 0.1410720E-01 0.713649370198 3.0000 1.11570 0.1628400E-0I 0.713650370198 4.0000 0.1866900E-031 0.71 3651370198 5.0000 1. I 5720 0.2162780E-01 0.713652370198 6.0000 1.20150 1.28204 0.2556080E-01 0.713653370198 7.0000 0.2917320E-01 0.713653870198 7.5000 1.36190 0.3548460E-01 0.7 I3654270198 7.9000 1.5 1560 1.59520 0.3845980E-03 0.7 13654370198 8.0000 0.5078400E-01 0.71 3654370698 8.0005 1.93464 14.79870 0.713654371 198 8.0010 0.1950097E+01 0.7 I3654372198 8.0020 14.75820 0.1974892E+Ol 14.72350 0.713654375198 8.0050 O.l975266E+Ol 0.71 3654380198 8.0100 14.70310 O.I975538E+Ol 14.66970 0.713654400198 8.0300 O.I975625E+Ol 0.713654420198 8.0500 14.66130 0.1975916E+01 0.7 I3654470198 8.1000 14.65500 O.l975538E+Ol 0.7 I3654570198 8.2000 14.60940 0.1976119E+01 14.59220 0.1976455E+01 0.7 13654670198 8.3000 0.713654770198 8.4000 14.62440 O.l975656E+Ol 0.713654870198 8.5000 14.59280 0.1975943E+01 14.58400 0.1976453E+OI 0.7 13654970198 8.6000 0.7 13655070198 8.7000 14.57560 O.I976528E+Ol 14.58540 O.I976375E+Ol 0.713655370198 9.0000 0.713656370198 10.0000 14.52100 O.l976954E+Ol 14.51100 0.1977091E+01 0.7 I3657370198 11.oooo 0.713658370198 12.0000 14.49900 0.1976718E+01 14.48880 O.l976879E+OI 0.713659370198 13.0000 0.713660370198 14.0000 14.48500 O.I977059E+Ol 14.44 100 0.713746370198 100.0000 O.l958988E+Ol
amplitude 2 0.6820000E-03 0.840000E-03 0.1554000E-02 0.1838000E-02 0.217oooOE-02 0.2548000E-02 0.298oooOE-02 0.357oooOE-02 0.406oooOE-02 0.5278000E-02 0.6976000E-02 0.7876000E-02 0.1221 800E-01 0.1600104E+OI O.l600740E+Ol 0.1601462E+01 0.1601342E+O 1 0.1600450E+01 0.1601354E+01 0.160 1208E+OI 0. I601204E+OI 0.1600600E+01 0.1 601032E+01 0.1601502E+01 0.1601 388E+01 0.160 1348E+O 1 0.1600664E+01 0.1601 342E+OI 0. I601 370E+OI 0.1601 170E+01 0.1601354E+01 0.1601 352E+0 1 0.1600434E+01
#Period is given relative to its value at the bifurcation. 0.639899358548, z, = 8.67695306202, and w = 0.103860729, as shown in Table 11. As f increases above the given value of fH, a small limit cycle with period increasing from T = 2 s / w = 60.48 15 will be formed around the steady state. This small limit cycle is shown in Figure 4. The amplitude and period of this limit cycle increase in a&rd with the formula given by Hassard et aL2' As Af = f -fH increases further, the differences between the calculated period (7') and amplitude (by the formula) and those computed by direct integration increase, indicating increasing deviations from the linear formula2' as expected. Thus at Af =
0.3 X IO", TIT, = 1,009586 by the formula and 1.00923 by direct integration, a difference of 0.33%; at Af = 1 X IOd, however, this difference increases to 2.5% (Le., 10-fold). Further changes in the period T and amplitude are shown in Table 111. In particular, one observes the rather dramatic change occurring for Af in the range 8.0005 X I@ < Af < 8.001 X la". While before and after these limits the period and amplitude change rather moderately, a jump of 7.6-fold is observed in the period, and 38-fold and 131-fold jumps are observed in they and z amplitudes, respectively. A plot of any of these vs Af will simply appear as a step function. This change is in fact a continuous
The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 8711
Canard Explosion and Excitation in the BZ Reaction
0.001
TABLE I V Canard Point Calculated from hurtion 10 ~
order k 0 1 2
contribution fka 0.7135870119698 0.0673142773311 0.0713276575059
approximation to@ 0.7135870119698 0.7136543262471 0.7136543975748
relative err09 9.44E-05 6.26E-08 3.73E-08
afkandfc calculated by the procedure given in the Appendix, e q s 10 and 48. *Relative error is defined as (fc - f ~ ( ~ ~ m a r i ~ ) ) / f ~ ( ~where ~~~"~~, fc(numanc)= 0.713654370948 is determined from Table 111.
one, and only its size makes it look as if it is not. This "jump" is the canard explosion discussed above. It occurs at fc = fH + Af = 0.713654370948. This value should be compared to the one calculated from developed in the Appendix for the merge of the stable (Ms) and manifolds. The results, up to the second order, unstable (Mu) are shown in Table 1V. The relative error is seen to be -4 X and can be improved by adding more terms in eq 10. The large limit cycle is shown in Figure 3. The "hump" on the upper left goes above point M ,the nearly vertical portions are located very near the stable manifolds MS,and the nearly horizontal lines are fast transitions between them, exactly as described under Theory. Note also the scale difference between Figures 3 and 4 (SOX for the abscissa and 125X for the ordinate), to appreciate the magnitude of the "explosion". Note also the close proximity of the canard point to the biX lod), a difference that will not be measured furcation point (4 in most experiments. Thus the experiments may show, in this case, a sudden "creation" of large relaxation oscillation without hysteresis. The experiment will not be able to differentiate between this phenomenon and others described in the introduction and elsewhere.1J8-20 For values off < fH, the system is stable, and for all initial points, it will revert to the steady state. However, the system is excirable; namely, in a range of initial points it will revert to the steady state via a small inward spiral, while at other initial points, above a certain threshold, the returning to the steady state will nearly follow the previously stable limit cycle, Le., the trajectory will make a very large excursion before the final spiraling to the steady state, depending whether the initial point is below or above the unstable manifold Mu. Observations of thresholds of excitability date back to FitzHugh,4 who studied the Boenhofer-van der Pol oscillator. This is also a singular perturbation system, which can be analyzed with the methods proposed here. In particular, we see that the somewhat vague concept of "QTP separatrix", given in the above reference, is, in our terminology, the unstable manifold Mu. As can be seen in Figure 5, calculated for Af =f -fH = -0.001, i.e.,fF = 0.712646370198 (located to the left of M), the steady state IS stable. The figure shows plots of ynullwhich is identical with S,the slow manifold, described earlier, znull,and Mu,the unstable manifold. The origin in this figure is taken as the stable steady state, Le., the intersection of ynulland znUll.The character of the unstable manifold, namely, that it comes very near the unstable part of the slow manifold (the right-hand side of ynull), is clearly seen. Initial points below the line marked Mu will revert to the steady state in a small spiral, similar in size to the limit cycle in Figure 4;initial points above that line will revert to the steady state by first moving along a trajectory very similar to the limit cycle shown in Figure 3. The trajectory will first move toward the right-hand side of S,will follow it nearly to point N, and then will go over quickly to the vicinity of M sin order to follow it to the vicinity of the steady state to begin the inward spiraling. Plots of these two kinds of trajectories can be found in Figure 4 of ref 27. The line M u is the unstable manifold described above. It is calculated here in two ways: (a) by back-integration from point Nand (b) by finding the "borderline" between "small" and "large" trajectories of the system reverting to the stable steady state after a small perturbation. The two methods give the same results up
-0.001
4 -0.003
-0.005
-
-0:01
0.00
0.01
0.02
AY Figure 5. Plots of ynull,znull, and the unstable manifold Mu,for Af = f -fH = -0.001, where the steady state is stable and is taken as the origin of the figure. Points below Mu will revert to the steady state via "small" trajectories, while points above it will revert by "large" ones. Note the proximity of Mu to ynullon the right.
to at least 5 digits and thus are the same on the scale of Figure 5. The experiments of Noszticzius et al., Finkeova et al.,z3 and RuofP4 show this effect. A BZ system in a steady state is perturbed by adding Br- ions (or taking it out by adding Ag+ ions). After adding small amounts of Br- ions, the system returns immediately to the steady state, while after adding greater amounts of the same, beyond a certain, well-defined threshold, the system returns to the steady state after a huge change in the species concentrations. The border between the two types of behavior is very sharp, exactly as described here.
Discussion The paper deals with two seemingly unrelated phenomena namely, the canard explosion and excitation. While the latter is well documented in the literature, both in its own sake23-28and in relation to wave p r ~ p a g a t i o n ~in' , ~excitable ~ media, the former is nearly unknown in the chemical community. We show here that both phenomena have the same mathematical source, namely, the existence and mutual relationship between the stable and unstable manifolds. The canard explosion deals with the sudden increase of both period and amplitude of a limit cycle very near a Hopf bifurcation. Experimentally, it implies that one may encounter difficulties in determining whether a particular bifurcation is super- or subcritical. Both of them may result in seemingly direct transition to large relaxation oscillations (hard transitions), without hysteresis. We may note in passing that computational difficulties in determining the nature of the bifurcation (whether super- or subcritical) may also be traced to the Occurrence of a canard.20 The model discussed in this paper has some advantage over the Brussellator discussed by Baer and Erneux and mentioned in the introduction.**' The steady state near which the canard explodes in the latter model occurs at a point where the concentration of one species tends to zero while that of the other tends to infinity. In the present two-dimensional Oregonator, on the other hand, the explosion of the stable limit cycle occurs, as shown, at an acceptable concentration range. An example of the "explosion" occurring in the case of a subcritical Hopf bifurcation, where an unstable limit cycle explodes, can be. seen in Figure 3 of ref 27. The experimental and numerical difficulties in detecting it will be obviously multiplied. The experiments3" in which small, nearly sinusoidal, fast oscillations become, through a very tiny change in system's con-
8712 The Journal of Physical Chemistry, Vol. 95, No. 22, 1991
straints, slow large ones can also be explained by invoking a canard explosion in the system-the two types of oscillations existing on the two 'sides" of the 'explosion". The very fast increase in period during a canard may be mistaken for an approach to a homoclinic orbit, and one must be extremely careful in establishing these type of orbits.20 Since both the excitation and the explosion relate to the same mathematical properties of the system, one might expect to find one phenomenon when the other exists. Another phenomenon connected with the above is the bursting, in which large oscillations alternate either with small ones or with quiescent time. Various authors39 have obtained bursting solutions by letting the stoichiometric factorfdepend on one of the variables of the Oregonator (or an extended version of it). Whenfmoves slowly above and below a subcritical Hopf bifurcation, the solution will alternate between large oscillations and a quiescent time. We now observe that even if the bifurcation is a supercritical one, having a canard point nearby, the system can move between large and small oscillations similar to the ones obtained by Bar-Eli and Noyes," using the Field-KOrOs-Noye~'~ (FKN) mechanism. Finally, it is important to stress that nonbifurcation transitions such as the canard explosion are not limited to 2-D systems. First, in simulating the 3-Dmodel, exactly the same will occur, since there is a two-dimensional invariant surface, where the system has a canard explosion (as can be seen in Figure 3). Second, stable and unstable manifolds of higher dimension may exist in systems with more state variables, and they may also interchange position at some value. However, things will be more complicated, since, for instance, surfaces may intersect in many ways. To the authors' knowledge, no theory or even computational studies of such systems, from the canard point of view, exist.
Appendix In this Appendix we shall give some mathematical details to complete the discussion given in the text. This will be given while using the particular equations of the applied model. Transformation of Coordinates. The three original equations given in the text are X = k l A Y - k2XY
+ k3AX - 2kyuZ
L=fiJ-kIAY-
k&Y
Z = k,AX - ksZ
(11)
(12)
T
klAt
p =dy/dr z f i - y - x y i = dz/dT = (X - Z ) / P
y = Yk2/k3A z = Zk2kS/klAk3A d = klA/k3A p = klA/k5 4 klAkd/k2k3A
x = (2 x
d x / d r = (x
+ y - xy - 2qx2)/er
(14)
(40) (a) Bar-Eli, K.; Noyes, R. M. J. Chem. fhys. 1988,88, 3646. (b) Noyes, R. M.; Bar-Eli, K. In Spatial Inhomogeneities and Transient Behauiour in Chemical Kinetics; Nicolis, G.. Gray, P., Scott, S.K., Eds.; Manchester University Press: Manchester, U.K.. 1990, p 289. (41) Field, R. J.; KBrBs, E.; Noyes. R. M. 1.Am. Chem. Soc. 1972, 94, 8649. (42) Gaspar, V.;Showalter, K. J. Phys. Chem. 1990. 94, 4973. (43) Scott, S. K.; Tomlin, A. S.Philos. Tram. R. Soc. London 1990,332A,
106))~
+ y - xy - 2qx2 = 0, one gets
Assuming that x = 0, Le., x xcv) =
1 - Y + 1/(1 - Y ) 2 + 84Y (17)
44
to be inserted in eqs 15 and 16, which become now the relevant two-dimensional system. Defining, further, w = fz and e = l/p = we get p = w - y - x b ) y = w - Fcy) (18)
= eCfxb) - W ) = tGCy,w,c)
(19) which have the same form as the equations given in the text, with c = f - f o (wherefo is defined below). Note that FCy) is not a function of w and G(y,w,c) is not a function of E . This fact simplifies the deviations a little. The curve ynullwhich fulfils p = 0 is the slow manifold S = S(y,w) = w - y - x b ) y = 0 defined in the text and shown schematically in Figures 1 and 2. To the left of point M,and to the right of point N , aS/ay = -w'(y) C 0; thus it is an attracting invariant manifold. Between points M and N , asJay = -w%) > 0, making a repulsive invariant manifold there. When t = 0 the curve S consists entirely of steady states, and the trajectories are the horizontal lines shown in Figure 1. Tbe Steady State. The steady state is located at the intersection of ynulland wnull,Le., where the two equations W
ws
= Y s + x,cylvs
(20)
wss = f,x,cy) (21) are fulfilled. In other words, when ynull= wnull,the system is at steady state. Solving the two equations for ysl, one obtains
3f,
+ 24 + 1 f d(3f,+ 24 + 1)2- 8cf,
Y, =
+f,*) (22)
4
In our case only the minus sign point will be of interest (see below). Thus givenf, the values of y,, xsl, and of course w, are uniquely determined. If the steady state is to be located at point M,namely, at the maximum of ynull,the equation 1
+ XO + yo(dX/dY)o = 0
(23)
should be fulfilled too, where (yo,wo)is both steady-state and maximum M,calculated at the parameter value fo. Hopf Bifurcation. The Jacobian matrix resulting from eqs 18 and 19 has the form
where x' is given by
5.
(44) Gray, B. F.; Thuraisingham, R. A. J . Eng. Mat. 1989, 23, 283. (45) Gray, P.; Scott, S.K. Chem. Eng. Sci. 1983, 39. 29. (46) (a) FitzHugh, R. Biophys. J . 1961, 1, 445. (b) FitzHugh, R. Bull. Math. Biophys. 1955. 17, 251. (47) L'Heureux, 1.; Kapral, R.; Bar-Eli, K. J . Chem. fhys. 1989.91,4285. (48) (a) Field, R. J.; Fbsterling. H. D. J . Phys. Chem. 1986, 90,5400. (b) Bar-Eli, K.; Field, R. J. J . Phys. Chem. 1990, 94, 3660. (49) (a) Troy, W. C.; Field, R. J . S I A M J . Appl. Math. 1977, 32, 306. (b) Troy, W. C. J . Math. Anal. Appl. 1977, 58, 233.
(16)
y = (2 x 1 0 6 ) ~ = (2 x 1 0 3 ) ~ T = looof With the parameter values: e' = 1 p = 1000 q = 0.025
Using these, the following differential equations are obtained: 3
(15)
With the rate constants given in Table I, the relationship between the dimensionless coordinates and the real concentrations is
(13)
The following transformation of c o o r d i n a t e ~ ~ ~can f ' * be ~ ~done in order to bring the above equations into nondimensional form:
x = Xk2/klA
Brans and Bar-Eli
X'
= dx/dy = -1 /4q
+
y-l+4q
(24)
44d-
The condition for Hopf bifurcation to occur is that the trace of the Jacobian matrix will be zero, namely: 1 +x+yx'+c=O
or
(25)
The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 8713
Canard Explosion and Excitation in the BZ Reaction t=-1-X-yx’
(26)
From eq 22 one obtains two values of y,, resulting, of course, in two values for t for which the Hopf bifurcation occurs. When, as in our case, t = 0.001 is given, this equation provides a unique y, - yH, i.e., a point of Hopf bifurcation (using the minus sign of eq 22). Comparing eqs 23 and 26, it is seen that the Hopf bifurcation is located to the right of the maximum M, the deviation being of the order t. The bifurcation can be super- (sub) critical, depending whether the stability numbep a3is negative (positive). In the terminology of Hassard et al.2” this is equivalent to p2. For our system (eqs 18 and 19) this becomes fx’l - F” 160, = -F“’ + -F” f X ’I-
.
F’ -
= -3xt1- yx“‘+ (2x’+ yx’?
the common root can be canceled by division. This result, eq 36, is in agreement with the previous definition offo as the value of the parameter where point M is the steady state (see eq 21). Proceeding further to terms of order t2, we get W,WO’+ wIwI’= f l X ( Y ) - WI (37) which is solved by w2
a negative number, which makes the bifurcation a supercritical one, and the limit cycles emerging from it stable. The eigenvalues at the Hopf point will be those of the Jacobian matrix when its trace is zero, namely:
* d W
=
fixcy) - WI - WlWl’
(38)
WO’
(27)
=
WO(Y0)
fo=.cvo)
fx” - 2x’- yx“ f x ’ - 1 - x - yx’
= -10.9343331636
AH
This has a singularity at yo (defined as the abscissa of M), since dwo/dy = w,’ = 0 there. However, if the numerator is also 0, i.e.
(28)
a purely imaginary number. At the maximum point M , the steady state will be stable with
and again this is not defined unless (39) This procedure can be continued successively to determinefk and W&+l at Yo. Determine nowfi from the definition of eq 39. Since w,’(Yo) = 0, it follows from Taylor’s theorem that G(Y) exists such that wo’(Y) = NY)(Y - Yo) and analogously for the numerator of eq 35:
fay) - w o w L
(40)
=
*(Ym- Yo)
Thus
i.e., the real part of the eigenvalue is negative. The Canard Point. Consider the model equations, namely, eqs 18 and 19. We look for a solution that stays close to the slow manifold S defined by ynull= 0, i.e., eq 20. From eqs 18 and 19 we find that the trajectories fulfill dw (w - Y(1 + X(Y” = 4 w ) - w) (29)
Differentiating eqs 40 and 41, we get %(Yo) = WO”(Y0) and so from eq 42:
*(Yo) =fix’(Yo)
(43)
We seek a solution to eq 29 in the form w(Y) = w&)
(44)
+ tWl(Y) + €2W2(Y)+ ...
(30) We cannot expect to find such a solution for all values of the parameter f, and therefore we expand this also:
f = f o + tfi
(41)
(This result could have been obtained directly from eq 35 by using L‘Hbpital’s rule.) Differentiating eq 42, we get *’(y)G(Y)
+ t2f2 + ...
(31) Inserting eqs 30 and 3 1 in eq 29 and collecting terms of the same order in t, we get to order to:
Wl%) =
- *(Y)$’cy) (45)
@W2
and by differentiating eqs 40 and 41 twice, we find $’(Yo)
= f/ZWo”’(Yo)
*’(Yo) = !4(fo~”CVo)
- wo”(Ya)) (46)
and from eq 35 which is solved by
wo = Y ( 1 + x(Y))
(33) The derivatives of wo can be calculated from eq 33, and by inserting eqs 44 and 45 into eq 39, an expression forfi is finally obtained:
This is the slow manifold S as expected. Proceeding to order t l , we get (34) whose solution is (35)
