Research: Science and Education
W
The Illustration of Multistability
Guy Schmitz* Faculté des Sciences Appliquées, Université Libre de Bruxelles, C.P. 165, Av. F. Roosevelt 50, 1050 Bruxelles, Belgium;
[email protected] Ljiljana Kolar-Anic´ and Slobodan Ani c´ Faculty for Physical Chemistry, University of Belgrade, 11001 Belgrade, Yugoslavia ˇ ˇ Cupi c´ Zeljko Department of Catalysis and Chemical Engineering, IChTM, 11000 Belgrade, Yugoslavia
The phenomenon of multistability is usually described (1, 2) by discussing the steady-state solutions of some nonlinear process as given by an abstract mathematical model having a single variable x and evolving in accordance with the differential equation
dx = x 3 + µx + λ dt
(1)
where µ and λ are the parameters characterizing the dynamic state of the system. After the comment that the steady states are represented by the solutions of the cubic equation
x ss3 – µx ss – λ = 0
(2)
the reader can find the beautiful illustration of the possible solutions for xss (Fig. 1). The three solutions of eq 2 can be either all real or one real and two conjugate complex, depending on the parameters µ and λ. Equation 1 was proposed by Landau in a study of the stability of steady flow of fluids (1, p 13). We discuss hereafter an example of application in chemistry.
Along with this, one usually finds comments about the introduced steady states. They can be stable or unstable. A steady state is stable if a small perturbation of the system tends to decay. It is unstable if the perturbation tends to grow, displacing the system to another state (1–6 ). In the case of Figure 1b, when λ is close to zero, there are three steady states that are either stable or unstable. The states lying on the upper and the lower branches of the S-shaped curve are stable, whereas the ones lying on its intermediate part are unstable. If the initial steady state of the system belongs to the upper branch and λ increases, it will first follow this upper branch. At the point λ2 it will jump onto the lower branch. Then, if λ decreases, the system will stay on the lower branch until the point λ1 is reached, whereafter it will jump onto the upper branch. Hence, in the vicinity of λ = 0, the state of the system depends on its history and we have multistability. A mechanical illustration is given by the following sketch. The ball can come to rest either in position “a” or in position “b”.
b
a
xss
We have multistability when many stable steady states exist for a given set of the values of the parameters. These steady states depend on the values of the parameters. A bifurcation occurs when their number or stability changes as the value of a parameter changes. In Figure 1b we see two bifurcation points, λ1 and λ2. Figure 1c is a bifurcation diagram showing the region of the parameters where we have bistability. In chemistry, the illustration of multistability is usually performed on a short model, such as eqs M1–M2, with an autocatalytic step and only one intermediate species X.
(a)
0
λ
µ
A + 2X
(b) xss
3X
k2
(M1)
(c) λ
0
λ1 0
λ
λ2
X
k3
0
0
µ
Figure 1. Effect of parameters µ and λ on the steady-state solutions xss in the vicinity of the bifurcation point. (a) The folded surface F(µ,λ,xss) = 0. (b) The intersection with a plane µ = const > 0. (c) The region in the (µ,λ) plane of existence of three real solutions.
1502
k1
k4
B
(M2)
At equilibrium, the requirement of detailed balance of both reactions has to be satisfied. Then the equilibrium concentration xeq of intermediate species X is given by
x eq =
k 1 a eq k2
=
k 4 b eq k3
Journal of Chemical Education • Vol. 77 No. 11 November 2000 • JChemEd.chem.wisc.edu
(3)
Research: Science and Education
where aeq and beq denote the equilibrium concentrations of species A and B. In nonequilibrium stationary states, or steady states, the overall reaction is A → B or B → A. The species A and B are both reactants and products. They are called external species (7) and their concentrations a and b define the dynamic state of the system. On the other hand, X is an internal species not appearing in the overall stoichiometric equation. In a steady state, its concentration x is a function of the concentrations of the external species. In the case of the model M1–M2, the time evolution of x is given by the following equation.
dx dt
= k 1ax 2 – k 2 x 3 – k 3 x + k 4 b
(4)
Here, the steady state solutions can be discussed as those of eq 1, since the cubic equation k2xss3 – k1axss2 + k3xss – k4b = 0
(5)
can be transformed by using the substitutions
x ss = y ss +
µ=
k1a
k a –1 1 k2 3 k2
k3
k a λ= 2 1 27 k 2
3
+
(6)
3k 2 2
k1k3a 3k 22
(7)
–
k4b k2
into the following equation analogous to eq 2. y ss3 + µ y ss – λ = 0
(8)
(9)
Although the multistability phenomenon is well described by mathematical models similar to eq 1 or chemical models of the type M1–M2, teachers would feel better if they could illustrate multistability with a more realistic model of an actual reaction. To this end, we present here a model for the decomposition of hydrogen peroxide into water and oxygen in the presence of iodate and hydrogen ions, the Bray–Liebhafsky (BL) oscillatory reaction (8, 9): IO3, H+
2H2O2 → O2 + 2H2O
(D)
Survey of the Bray–Liebhafsky Oscillatory Reaction The simple global reaction (D) comprises a complex homogeneous oscillatory reaction (10–13 and references therein). Already in 1921, Bray proposed that the overall reaction proceeds through at least two complex reactions. One is the reduction of iodate by hydrogen peroxide (R) and the other is the oxidation of iodine by the same agent (O). 2IO3 + 2H+ + 5H2O2 → I2 + 5O2 + 6H2O
(R)
I2 + 5H2O2 → 2IO3 + 2H + 4H2O
(O)
+
The sum of reactions R and O gives the overall decomposition, D. Normally, their rates tend to become equal and we usually observe only the smooth decomposition. In a
narrow region of concentrations however, it is also possible that processes R and O are alternately dominant; the iodine concentration increases and decreases alternately and the reaction is periodic. Reactions R and O are themselves complex and involve numerous intermediates such as I, HIO, I2O, and HIO2. A model of these reactions ought to include several steps and at least be able to describe the oscillatory evolution of the overall process. At present, there is one such model with different variants (13–20). Although the version with eight steps (16 ) is the best one for describing the real process, the following simple model seems to be the most appropriate for pedagogical purposes. IO3 + I + 2H+ → ← HIO + HIO2
(R1,R1)
HIO2 + I + H → I2O + H2O
+
I2O + H2O → ← 2HIO +→ HIO + I + H ← I2 + H2O
(R2) (R3,R3) (R4,R4)
HIO + H2O2 → I + O2 + H+ + H2O
(R5)
I2O + H2O2 → HIO + HIO2
(R6)
The sum 2×(R1) + 2×(R2) + 2×(R3) + (R4) + 5×(R5) gives reaction R and the sum 2×(R1) + 3×(R2) + 2×(R3) + (R4) + 5×(R6) gives reaction O. Among the ten species, four are internal ones, namely I , HIO, I2O, and HIO2, and play the same role as X in the model M1–M2. During the oscillatory evolution, two species appearing in reactions R and O, namely I2 and H2O2, are external species. They play the same role as A and B in the model M1–M2. The concentrations of H+, IO3, and H2O can be taken to be constant (18). In this model, oxygen is a product without any influence on the time evolution of the reaction. The steady states of the system are given by the stationary solutions of four differential equations corresponding to the four internal species. Therefore, it seems that for an illustration of multistability in this system, we need four state variables and two parameters, which would involve a (4+2)-dimensional phase space. However, it suffices to consider a subspace consisting of the concentrations of three properly chosen species: the two external species defining the state of the system, H2O2 and I2, and one internal species, I. The four algebraic equations d[I]/dt = d[IOH]/dt = d[HIO2]/dt = d[I2O]/dt = 0 can be used to express the steady-state concentration of iodide [I]ss as a function of [H2O2] and [I2] for constant values of [IO3] and [H+]. The result is shown in Figure 2a. Multistability in the Model of the Bray–Liebhafsky Reaction After the above survey, it can be understood that the system is driven by process D, the sum of processes R and O. In an actual situation, these processes are always far from equilibrium, and their rates are complex functions of [H2O2] and [I2]. The hydrogen peroxide concentration is the classical control or bifurcation parameter, such as µ in eq 9. It is much larger than the concentration of iodine1 and is quasi constant in an open reactor. It decreases in a closed reactor but its variations during one oscillation are small enough that it can be regarded as constant. The iodine concentration is the second
JChemEd.chem.wisc.edu • Vol. 77 No. 11 November 2000 • Journal of Chemical Education
1503
Research: Science and Education
[I]ss
(a) F E G
[I2] F’
[H2O2]
E’
G’
(b)
(c)
[I]ss
[I 2]
E’
F’
G’
[I2]1 [I2]2
[I2]
[H2O2]
[I2] / (108 mol / dm3 )
Figure 2. Effect of H2O2 and I2 concentrations on the steady-state concentration of I in the vicinity of the bifurcation point F. The rate constants for the numerical calculations are taken from ref 16. (a) The folded surface F([H 2O2], [I 2], [I ] ss) = 0. (b) The intersection with a plane [H2O2] = const. (c) The region of existence of three real solutions.
3.0
2.5
2.0 0
30
60
Time / min [I] / (108 mol / dm3 )
control parameter, as λ in eq 9. It increases if process R is dominant and decreases if process O is dominant. In Figure 2a, one can see that at low H2O2 concentrations there is only one steady state. It is stable (15, 17 ). At higher concentrations, the surface is folded and three steady states exist. The two external ones are stable and the intermediate one is unstable; hence one encounters bistability. The point F denotes the transition from monostability to bistability. Figure 2c shows the corresponding region of existence of three real solutions. The intersection of the surface in Figure 2a with a plane corresponding to a given high concentration of H2O2 gives Figure 2b. The concentrations for which the rates of processes R and O are equal are given by the curve EFG in Figure 2a. Mathematically, it is deduced from the model R1–R6 by adding the condition d[I2]/dt = 0 to the other steady-state conditions. The part EF of this curve is on the stable part of the surface. For low hydrogen peroxide concentrations, the point representing the steady state of the system moves on the surface toward curve EF, on which one has the smooth catalysis. For high hydrogen peroxide concentrations, the curve FG is on the unstable part of the folded surface and cannot be reached. The system takes states either on the upper part or on the lower part. When it is in a steady state on the upper part, [I]ss is relatively large. Then reactions R1 and R4 take place from left to right. The dominant reaction is R and the iodine concentration increases. The point representing the steady state of the system moves toward the edge of the upper part of the folded surface where it falls on the lower part. There [I]ss is relatively low. Reactions R1 and R4 take place from right to left. The dominant reaction is O and the iodine concentration decreases. The point representing the steady state of the system moves in the opposite direction, reaches the edge of the lower part of the surface, and jumps onto the upper part. Its motion continues cycling around the unstable part of the surface. In an open reactor, the concentration of hydrogen peroxide is the same at each cycle. We get a limit cycle and sustained oscillations. In a closed reactor the concentration of hydrogen peroxide decreases slowly. The point representing the steady state of the system spirals around the unstable part of the folded surface. We get damped oscillations. When the system reaches the bifurcation point F the oscillations vanish. This simple explanation of the BL oscillations is illustrated in Figure 3. The time evolution of the concentrations was obtained by numerical integration of the differential equations associated with model R1–R6.
3
2
1 0
30
60
Time / min
Conclusion Nonlinear phenomena such as oscillating chemical reactions and chaos are illustrated by quite a number of papers in a very pedagogical way (21–24 ). However, we feel that the explanation of multistability, a basic problem for understanding the phenomena mentioned, is not yet offered in a form acceptable for students. We have illustrated it using a simplified model of a real chemical reaction. The possible steady states of the system are represented in a three-dimensional phase space, with the concentrations of three species that can be monitored experimentally plotted along the three coordinate axes.
1504
Figure 3. Numerical simulation of the BL reaction based on model R1–R6.
Acknowledgments ˇ Cupi´ ˇ c, and S. Ani´c gratefully acknowlLj. Kolar-Ani´c, Z. edge the partial support of the Fund for Science and Technology of Serbia and Federal Ministry for Development, Science and Environment of Yugoslavia. We also thank Prof. B. Mili´c for his comments.
Journal of Chemical Education • Vol. 77 No. 11 November 2000 • JChemEd.chem.wisc.edu
Research: Science and Education W
Supplemental Material
Supplemental material for this article is available in this issue of JCE Online. Note 1. In typical experiments, the concentration of hydrogen peroxide is about 0.1 mol/dm3 and the concentration of iodine is less than 10-3 mol/dm3.
Literature Cited 1. Drazin, P. G. Nonlinear Systems; Cambridge University Press: Cambridge, 1994. 2. Nicolis, G.; Prigogine, I. Exploring Complexity; Freeman: New York, 1989. 3. Gray, P.; Scott, S. K. Chemical Oscillations and Instabilities, Non-linear Chemical Kinetics; Oxford University Press: Oxford, 1990. 4. Gray, P.; Scott, S. K. J. Phys. Chem. 1985, 89, 22. 5. Gray, P. Proc. R. Soc. London, Ser. A 1988, 415, 1. 6. De Kepper, P.; Boissonade, J. In Oscillations and Traveling Waves in Chemical Systems; Field, R. J.; Burger, M. Eds.; Wiley Interscience: New York, 1985; p 223. 7. Clarke, B. L. In Advances in Chemical Physics, Vol. 43; Prigogine, I.; Rice, S. A., Eds.; Wiley: New York, 1980; p 1. 8. Bray, W. C. J. Am. Chem. Soc. 1921, 43, 1262.
9. Bray, W. C.; Liebhafsky, H. A. J. Am. Chem. Soc. 1931, 53, 38. 10. Noyes, R. M. J. Phys. Chem. 1990, 94, 4404. 11. Furrow, S. D. In Oscillations and Travelling Waves in Chemical Systems; Field, R. J.; Burger, M. Eds.; Wiley Interscience: New York, 1985; p 171. 12. Treindl, L.; Noyes, R. M. J. Phys. Chem. 1993, 97, 11354. 13. Schmitz, G. J. Chim. Phys. 1987, 84, 957. 14. Schmitz, G. J. Chim. Phys. 1991, 88, 15. 15. Kolar-Ani´c, Lj.; Schmitz, G. J. Chem. Soc., Faraday Trans. 1992, 88, 2343. 16. Kolar-Ani´c, Lj.; Miˇsljenovi´c, Dj.; Ani´c, S.; Nicolis, G. React. Kinet. Catal. Lett. 1995, 54, 35. 17. Kolar-Ani´c, Lj.; Vukeli´c, N.; Miˇsljenovi´c, Dj.; Ani´c, S. J. Serb. Chem. Soc. 1995, 60, 1005, 1187. ˇ c , Z.; ˇ Ani´c, S.; Schmitz, G. J. Chem. 18. Kolar-Ani´c , Lj.; Cupi´ Soc., Faraday Trans. 1997, 93, 2147. ´ S.; Kolar-Ani´c, Lj.; Körös, E. React. Kinet. Catal. Lett. 19. Ani c, 1997, 61, 111. ˇ Kolar-Ani´c , Lj. J. Chem. Phys. 1999, 110, 3951. ˇ c´ , Z.; 20. Cupi 21. Field, R. J.; Schneider, F. W. J. Chem. Educ. 1989, 66, 195. 22. Pojman, J. A.; Craven, R.; Leard, D. C. J. Chem. Educ. 1994, 71, 84. 23. Benini, O.; Cervellati, R.; Fetto, P. J. Chem. Educ. 1996, 73, 865. Strizhak, P.; Menzinger, M. J. Chem. Educ. 1996, 73, 868. 24. Copper, C. L.; Koubek, E. J. Chem. Educ. 1998, 75, 87.
JChemEd.chem.wisc.edu • Vol. 77 No. 11 November 2000 • Journal of Chemical Education
1505
