feature article - WSU Math

---

J. Phys. Chem. 1985, 89, 22-32

22

FEATURE ARTICLE Sustained Oscillations and Other Exotic Patterns of Behavior in Isothermal Reactions P. Gray* and S . K. Scott School of Chemistry, University of Leeds, Leeds LS2 9JT, UK (Received: June 19, 1984)

- -

The aims of this Feature Article are to show how simple or complex isothermal autocatalytic reactions can give rise to oscillatory behavior by considering virtually the simplest possible example (A + 2B 3B; B C) under the simplest and most readily realizable of circumstances (cstr). Although this is only a two-dimensional scheme it vividly and clearly represents all the key aspects of many real systems, including the birth, growth, and extinction of stable oscillations. Moreover, it is a “robust” model and it does not collapse when reversibility and parallel alternative routes are added. In the recent past, mathematical schemes have usually been either too elaborate to be understood or have failed to satisfy such basic requirements as the principle of detailed balance. Yet other models have failed to conserve mass or have not represented stable oscillations. Such fatal flaws are absent not only from this prototype but also from its comparison (A + B 2B; B C saturating). The springs of chaotic behavior in deterministic chemical models are also touched on, but it should be noted that real cstr must fall short of perfect homogeneity and that there may still be a gap between expectation and experimental proof of chaos.

- -

I. Introduction Ever since the revival of interest in oscillatory behavior in isothermal systems, elaborate mechanisms and complex models1** have been the order of the day. Fundamental aspects have often been obscured. By contrast, the nonisothermal field has been well served ever since the mid-1950s by studies3s4 in depth of the simplest possible example of thermal feedback. These efforts have been rewarded: significant new propertiesS were still being found into the 1980s. The aim of this feature article is to perform the same service for isothermal systems with chemical feedback (isothermal autocatalyses) in the continuously flowing, well-stirred, “tank” reactor (cstr). There is no substitute for the hard and careful work on individual systems that is necessary to establish detailed mechanisms, but by showing the rich variety of behavior that basic models with very simple features can generate, we may curb excessive elaboration. I . Advantages of Open Systems. In closed systems there is no exchange of matter with the surroundings, and chemical reactions proceed to equilibrium. Closed systems are often very convenient experimentally. They offer economy of materials and ease of precise control. Because their final state is one of equilibrium, chemical thermodynamics is a powerful aid to interpretation. On the other hand, stationary states other than the final one are not in reach, multistability is impossible, and oscillatory behavior cannot be sustained. In open systems, where there is a continuous supply of fresh reactants and a matching removal of products, truly steady states are easily realized and oscillatory conditions can be maintained indefinitely. 2. Advantages of Good Stirring. Chemical reactions in open systems abound in nature (streams, fires, living organisms, single cells), and they are commonly associated with nonuniform concentration and temperature fields. Their mathematical treatment is correspondingly arduous, abounding in partial differential equations. The simplest of open systems to treat theoretically is one in which conditions are uniform, i.e. one with “perfect mixing”. It is also one of those that can be approximated experimentally, in (1) Gray, P. Ber. Bunsenges. Phys. Chem. 1980, 84, 309. (2) Gray, B. F. Spec. Per. Rep.: React. Kine?. 1975, 1 , 309. (3) van Heerden, C. Ind. Eng. Chem. 1953.45, 1242. (4) Uppal, A.: Ray, W. H.; Poore, A. B. Chem. Eng. Sci. 1974, 29, 967; 1976, 31, 205. (5) Balakotaiah, V.; Luss, D. Chem. Eng. Sci. 1983, 38, 1709; Chem. Eng. Commun. 1981, 13, 1 1 1 .

0022-3654/85/2089-0022$01.50/0

the continuously flowing, well-stirred, “tank” reactor (cstr): good stirring can come close to perfect mixing. The distribution of residence times is simple and is fully characterized by its mean value; partial differential equations are banished and replaced by ordinary ones. Finally, in stationary states, we normally have to deal only with algebraic relations. Since the pioneering work6 of Denbigh in the 1940s, reinforced in recent years by the Bordeaux school,’ solution chemists have turned increasingly to the cstr to study oscillatory reactions. Historically, these approaches led to the discovery of oscillatory behavior before that of multistability whereas the easiest route to the interpretation of these phenomena starts with the idea of multiple stationary states. Both phenomena, however, require some form of mechanistic feedback, and for both autocatalysis provides a simple chemical mechanism. 3. Autocatalysis as Isothermal Feedback. The feature common to nearly all isothermal oscillatory reactions is autocatalysis. (The obverse property of autoinhibition is considered briefly in section VII.) Autocatalysis is shared not only by the Belousov-Zhabotinskii reaction,”I0 a wide range of halide-based oscillators,” and the arsenite plus iodate reaction’* but by numerous enzyme systemsI3 and by chain-branching reactions in the gas phase.I4 All mechanistic schemes include autocatalytic steps (or autocatalytic combinations of elementary steps). However, even the simplified forms15J6of these schemes (Brusselator, Oregonator) are still complicated, and there is still a need to look at even simpler prototypes in a logically ordered way. Accordingly, we consider throughout steps with the overall stoichiometry A B. We may also take this as written to represent a reaction satisfying the rate law = ka. When the same change is catalyzed by a species Y

-

(6) Denbigh, K. G. Trans. Faraday SOC.1944, 40, 352; 1947, 43, 648; 1948, 44, 263, 479; J . Appl. Chem. 1951, 1, 227. (7) Vidal, C.; Pacault, A. “on-Linear Phenomena in Chemical Dynamics”; Springer-Verlag: New York, 1981; pp 10-14. (8) Belousov, B. P. Sb. Ref. Radiat. Med. 1959, 1 , 145. (9) Zhabotinskii, A. M. Dokl. Akad. Nauk. S S S R 1964,157, 392; Biof izika 1964, 9, 306. (10) Noyes, R. M.; Field, R. F.; Koros, E. J . Am. Chem. SOC.1972, 94, 1394. (11) Epstein, I. R.; Orban, M. ‘Oscillations and Travelling Waves in Chemical Systems”; Field, R. J., Burger, M., a s . ; Wiley: New York, 1984. (12) (a) Lintz, H.-G.; Weber, W . Chem. Eng. Sci. 1980, 35, 203. (b) Papsin, G . A,; Hanna, A.; Showalter, K. J. Phys. Chem. 1981, 85, 2575. (13) Aarons, L. J.; Gray, B. F. Chem SOC.Rev. 1976, 5, 359. (14) Semenov, N . N. “Chain Reaction”; Clarendon Press: Oxford, 1935. ( 1 5 ) Prigogine, I.; Lefever, R. J. Chem. Phys. 1968, 48, 1695. (16) Field, R. J.; Noyes, R. M . J . Chem. Phys. 1974, 60, 1877.

0 1985 American Chemical Society

Feature Article

The Journal of Physical Chemistry, Vol. 89, No. 1, 1985 23

-+

and has a rate satisfying the expression kay, we may represent this as A Y B Y. Autocatalysis corresponds to catalysis by the product B itself. Two exemplary span the whole range of behavior likely to be encountered:

+

rate = k,ab

quadratic autocatalysis cubic autocatalysis

rate = k,ab2

(la) (1b)

These may be represented by the “chemical” equations A+B-2B and A

+ 2B

-

(IC) 3B

(14

stract, but the algebra is not much more demanding than that involved in locating multiple stationary states. Sustained oscillations are possible in these two-dimensional models. Section V presents the qualitative, and some quantitative, aspects but does not dwell on the additional algebraic effort required. Chaos requires a third dimension. The circumstances in which there is a chaotic response to very small periodic perturbations (e.g. in pumping speeds or bath temperatures) are considered in section VI. Section VI1 shows how, when decay rates are not so simple as first order, even quadratic autocatalysis is sufficient to generate much of the above behavior. Finally, section VI11 considers some of the wider implications of these and related modern studies of similar systems.

Normally these will arise from a subscheme of elementary steps. For instance, the arsenite-iodate reaction produces iodide ions with an experimentally determined rate law of the formI2v2*

11. Cubic Autocatalysis in Open Systems: Kinetic Model and Mass Balance in a Well-Stirred Reactor

d[I-]/dt = ( k , + k2[I-]][I-][I03-][H+]2

Suppose there is a chemical reaction which can be represented by a first-order step

This does not, however, require us to believe that it requires a triple encounter as an elementary step. We return to questions of this kind in section VIII. It is important to note (i) that reversibility makes little difference to behavior until the equilibrium constant ([B]/[A]), becomes small and (ii) that less simple stoichiometry A + nB ( n + m)B can readily be incorporated if required.20*21 In a closed vessel, the rate of an autocatalytic reaction22passes through a maximum and product concentrations plotted against time lie on S-shaped curves. Such behavior is typical of a very wide range of systems indeed, ranging from many gas-solid reactions and solidsolid interconversions to the growth of biological populations constrained by a limited food supply and the spread of infectious diseases. Lotkaz3 provides a classical survey still attractive today. We shall concentrate mainly on the cubic rate law (lb), although in section VI1 we study simple or quadratic autocatalysis. Steps like (1 b) lie at the heart of many complex models. We shall see that it does not need a complicated pattern of auxiliary reactions to generate very exotic behavior indeed, despite much folklore to the contrary. 4. Arrangement of Text. The present article is concerned with studying the various patterns of behavior arising from the simplest model schemes based on reactions ICand Id. Where appropriate, we shall mention briefly experimental observations of particular features. The recent feature article by EpsteinZ4offers an admirable and up-to-date survey of the chemistry of some real systems: we set out a complementary account of the necessary analysis. Section I1 develops the concept of multistability and its relation to discontinuous jumps from one stationary state to another and the origin of hysteresis. These ideas are familiar to chemists in other contexts and can be vividly illustrated by means of a flow diagram. In section I1 we consider the simple situation where the catalyst B is indefinitely stable, undergoing no subsequent reaction. The composition of such systems requires only one variable for its complete description at all times. In section I11 we uncouple the concentrations of A and B by allowing the first-order decay B-C

-

to remove the catalyst. Two variables are now necessary to describe the evolution of the system. Section IV develops the techniques of local stability analysis. These ideas are more ab(17) w a y , P.; Scott, S . K. Chem. Eng. Sci. 1983,38,29. (18) Gray, P.; Scott, S. K. Chem. Eng. Sci. 1984,39, 1087. (19) Scott, S. K. Chem. Eng. Sci. 1983,38, 1701. (20) Lin, K. F. Chem. Eng. Sci. 1981,36,1447. (21) Scott, S . K. J . Chem. Soc., Faraday Trans. 2, in press. (22) Boudart, M. “Kinetics of Chemical Processes”; Prentice-Hall: Englewood Cliffs, NJ, 1968. (23) Loth, A. J. ‘Elements of Mathematical Biology”; Dover Publications: New York, 1956. (24) Epstein, I. R. J. Phys. Chem. 1984,88, 187.

A

-

rate = k3a

B

(2)

(we choose the subscript 3 for this step to remain consistent with work already p ~ b l i s h e d ’ ~ - ’ ~ Suppose ~ ~ ~ ) . also that this is accompanied by autocatalysis so that we can add either of the steps I C or Id. Then, if these are the only steps in the mechanism, the composition of the system at any time is fixed if we specify either the concentration of A or the concentration of B, because of the additional requirement [A10 = [AI + [BI

(3)

(assuming the system contained pure A initially), where [AIo is the concentration of A at t = 0. For this system, even if there was an initial nonzero concentration of B, the reactant and catalyst concentrations cannot vary independently. These are one-parameter systems. In many physical situations, however, the catalyst may not be indefinitely stable but will undergo further reaction, decomposition, or poisoning. The simplest decay of B is via a first-order process

B

+

C

rate = k2b

(4)

This adds another dimension to the system, as there are now three species concentrations and still only one relationship between them. Most importantly, the concentrations of A and B may now vary independently of each other-a requirement for oscillatory behavior. Almost all of the special features of these autocatalytic systems can be highlighted by studying the simplest case in the simplest circumstances. These correspond to the situation where the autocatalytic reaction swamps the uncatalyzed step. Accordingly, we shall study the system A

-

+ 2B B

3B

C

rate = klab2 rate = k2b

(5a) (5b)

Such schemes have zero rate when the concentration of B is zero (Le. at zero conversion). We may alter this by allowing a nonzero initial concentration of B, bo # 0. Before proceeding to the mass-conservation equations for this mechanism, we note that the product klao2has units of s-l. It is the inverse of a characteristic chemical time scale which plays an important role in determining the behavior observed under closed and open conditions. A second characteristic time, for the lifetime of the catalyst, is given by the inverse of the first-order rate constant k2 t2

= l/k2

(7)

(25) (a) Gray, P.; Scott, S . K. J. Phys. Chem. 1983,87,1835. (b) Gray, P.; Scott, S . K. Ber. Bunsenges. Phys. Chem. 1983,87,379. (c) Gray, B. F.; Gray, P.; Scott, S . K. J . Chem. Soc., Faraday Trans. I , in press.

Gray and Scott

The Journal of Physical Chemistry, Vol. 89, No. 1, 1985

24

Stable catalysts decay slowly and have long values of t2; unstable catalysts are characterized by short lifetimes. 1. Mass Balance Equations in Open Systems. For the kinetic mechanism ( 5 ) we may write down the two mass-balance equations for the concentrations of A and B:

where k f = l/tresis the inverse of the mean residence time. Note that we do not have to consider a third differential equation for the concentration of the product C as this is uniquely determined at all times by the concentrations of A and B in the reactor and the inlet conditions c = (a0 bo) - ( a b )

+

+

Equation 8 may be compactly recast in terms of a small number of dimensionless groups a = a/ao dimensionless reactant concentration @ = b/ao

dimensionless catalyst concentration

y. R is a cubic curve which passes through a minimum at the origin, corresponding to zero rate at zero conversion; it shows an inflection point at y = ]I3,reaches a maximum at y = 2 / 3 , and falls to zero again at y = 1, corresponding to complete conversion of A to B. The maximum value of R = 4/27 corresponds to a real maximum reaction rate

y = (ao - a ) / a o dimensionless extent of conversion

Po = bo/ao dimensionless inlet concentration of catalyst T,,,

= t,,,/t,h

= t2/tch

i2 T

Figure 1. Cubic autocatalysis, with infinitely stable catalyst B: (a) dependence of reaction curve R and flow line L on extent of conversion y = 1 - a with no inflow of B; (b) variation of stationary state with residence time corresponding to (a) showing extinction or washout; (c) dependence of R and L on y for nonzero inflow of B; (d) corresponding variation of ylr with ires showing ’ignition”, extinction, and hysteresis.

dimensionless residence time dimensionless catalyst lifetime dimensionless time

= l/tch

This yields da 1-a _ - -a@ + d7

Tres

The two dimensionless measures of the reactant concentration a and y are simply related by a=l-y

- -

In the limit of perfect catalyst stability ( k 2 0; t 2 system is reduced to one of a single variable with

p

= 1

+ Po - a = y + 8,

(10) m)

our (11)

For the general case of k2 > 0, y and @ are no longer equivalent measures (for Po = 0) nor even directly related. Nevertheless, an investigation of the extreme situation of an indefinitely stable catalyst, with simplification (1l), provides us with valuable insights which will be of use in understanding the more general case. 2. Stationary-State Extents of Conversion with Infinitely Srable Catalyst: Flow Diagrams. For systems in which the catalyst does not decay ( r 2 m) and for which there is no catalyst B in the inflow, we may use eq 11 to obtain a single mass-balance equation. In terms of the dimensionless extent of conversion y, this has the form

-

(note under these conditions d y l d r = -da/dr = d@/dr.) We denote the two terms on the right-hand side of (12) by R and L. These represent the rate of chemical conversion of A to B and the net rate of inflow of A, respectively. Figure l a shows how R and L vary with the extent of conversion

This finite bound on the rate distinguishes our model from some of the improperly formulated schemes26or “pool chemical” approaches in which the concentration of A is assumed constant and for which the reaction rate can become infinitely large. The flow line L is a straight line passing through the origin with a gradient equal to l/Tre,. For slow flow rates (long residence times), L is relatively flat; for fast flows, L is steep. Thus, the effect of varying the flow rate or residence time in our system is to vary the position of L relative to R (which is unaffected by changes in r,,). Stationary states of the reaction occur when the two rates R and L are in exact balance and d y l d r becomes zero. Such conditions are represented by the intersections of R and L on the flow diagram (Figure la). Any flow rate for which L intersects R three times corresponds to a system with three possible stationary states (multistability). We call these y ] , y2, and y3 in order of increasing extent of conversion. The lowest intersection is at the origin y 1 = 0. Here there is no reaction, and the outflow and inflow consist of pure A. The other two stationary states lie in the ranges 0

< y2 < y2 and ‘/z < y3 < 1

As the flow rate is increased, the flow line L gets steeper. This causes y2and -3 , to move closer together. For the special residence time T,,,

=4

1

i.e. -tre, = 1 / k l a o 2 4

R and L are tangential and y2 and y3 merge at the value ]I2. For any shorter value of r, yzand y3disappear and the system jumps to the only remaining intersection yIcorresponding to no reaction. This jump is known as extinction or washout. The dependence of the stationary-state extent of conversion on residence time swept out in this way is shown in Figure lb. ( 2 6 ) Karmann,

K.-P.: Hinze, J. J . Chem. Phys. 1980, 72, 5467.

The Journal of Physical Chemistry, Vol. 89, No. I, I985

Feature Article

3. Effect of Adding Catalyst to Inflow: Infinitely Stable Catalyst. If the concentration of B in the inflow is not zero, two new features arise. The mass-balance equation is now

The changes brought about by nonzero Bo occur only in the reaction rate term R . The rate of chemical conversion of A to B is now no longer zero for y = 0 and hence R does not pass through the origin but is "lifted up" to a finite rate (Figure IC). This also means that R and L, which are unaffected by the value of Po, no longer intersect at the origin. There may now therefore be two flow rates at which R and L become tangential, given by Between these values there are three intersections of the two curves and, hence, three stationary states. As the residence time is decreased toward the upper tangency, y2and y3 merge and extinction occurs, as before. At the lower tangency, corresponding to longer residence times, y1and y2 merge and the system jumps to y3. This is known as ignition. The variation of the residence time causes the stationary-state extent of conversion to sweep out the S-shaped curve shown in Figure Id. 4 . Stability of Stationary States: Relaxation Times and Slowing Down. The stability of a given stationary state may be tested by seeing if a small perturbation Ay = y - y,, grows (unstable stationary state) or decays to zero (stable stationary state). The time dependence of A y is given by

AT = A ~ o / ( l +

In general, the time taken for the perturbation to decrease from Ayo/ento Ayo/e"+' is given by 7*,+,

111. Systems with Catalyst Decay: Patterns of Stationary States (Isolas and Mushrooms)

Superficially, the change in the system brought about by allowing the catalyst to have a finite, rather than infinite, lifetime might be expected to be rather modest. In fact, profound alt e r a t i o n ~ are ' ~ made to the patterns of stationary states and their stability. We now have three chemical species in our scheme

-

+ 2B

0) as the time taken for Ay to decay to l / e of Ay,. This gives 7* = 2(e - l ) / A y ,

A where X and p reflect the dependence of dy1d.r on the extent of conversion

W 0 7 )

25

=

a+b+c sum of

(18)

concentrations = concentrations in inflow in reactor arising from the reaction stoichiometry. We thus have two independent concentrations in our system. For convenience we shall only consider situations in which there is no final stable product C added to the inflow, co = 0. I . Stationary-State Solutions and Flow Diagrams: Isolas for Systems with No Catalyst Inflow. A stationary state is achieved when the rates of change da/d7 and dP/d7 become simultaneously zero. We may add eq 9a and 9b to yield a relationship between CY,, and P,,, the dimensionless concentrations of A and B in the stationary state. For a system with no autocatalyst in the inflow, Po = 0, this yields

Comparing this expression with eq 1 1 , we see that P, is no longer equal to the dimensionless extent of conversion of the reactant but is reduced by the factor 1 7,,,/72. This difference reflects 7 0 = l / y 3 ( l - 2y3) for y3 the fact that some of the B formed from A has subsequently been converted to C. The longer the residence time, the greater the As we approach the extinction tangency, we have y2 y3 quotient 7,/72 and, hence, the higher the probability that a given 1/2. For these values, the relaxation times become infinitely large molecule will be converted to the final product C. as X tends to zero. In fact, infinite values cannot be a c h i e ~ e d , ~ ~ , ~ ' If we now substitute for @, from (19) into ( 9 ) ,we may obtain although there is a lengthening of 7. (Le. perturbations decay more for the stationary-state concentrations slowly). As X approaches zero, inequality 15 cannot be satisfied no matter how small Ayo is made. We must then consider both terms on the right-hand side of (14). For X = 0, and noting that when y2 = y3 = ll2,we replace (16) by p = 7r

= l / y 2 ( l - 2 y 2 ) for y2

+

--

(27) Heinrichs, M.; Schneider, F. W . J. Phys. Chem. 1981.85.21 12; Ber. Bumenges. Phys. Chem. 1980,84, 857. (28) Ganapathisubramanian, N.; Showalter, K.J. Phys. Chem. 1984.87, 1098; J . Chem. Phys. 1984,80, 4177.

R

=

L

26

Gray and Scott

The Journal of Physical Chemistry, Vol. 89, No. 1, I985

of the isola correspond to the two residence times for which L is tangential to R, i.e. for which

I

Figure 2. Cubic autocatalysis with catalyst of finite lifetime, no B in inflow: (a) dependence of R and L on stationary-state extent of conversion ysa;(b) variation of ysswith irefor unstable catalyst ( T < ~ 16), only the axis yy, = 0 represents intersections of R and L;(c) variation of ymwith rre for i2= 16 showing birth of an isola at yI = and 7,= 16; (d) variation of ysswith T~~ for stable catalyst (i2> 16) showing isola with two points of extinction.

We may represent the reaction rate and flow terms R and L as functions of ysson a flow diagram (Figure 2a). Again, R is a cubic curve passing through the origin yss = 0 and falling back to zero at the state of complete conversion yS = 1. The flow term L is still a linear function of ysspassing through the origin, but now the gradient is a more complicated function of the residence time: 2

gradient of L = dL/dy For the shortest residence times (highest flow rates) we may neglect the term rr,/rZ as small compared with unity. The gradient then is proportional to 1/rreS,and the flow line becomes less steep as the residence time decreases. At the longest residence times (lowest flow rates), 7,,/i2 is very much larger than unity, so the gradient becomes proportional to rIa. The flow line now becomes steeper as the residence time increases. This differs from the stable catalyst case studied in increases. section 11, where the gradient always decreases as ires In between these extremes of long and short residence times, the flow line achieves a minimum slope. This occurs when the mean residence time is equal to the characteristic lifetime (or half-life) of the catalyst, Le. when 7, = 72. The minimum gradient is given by It is the relative positions of this flow line with minimum slope and the line passing through the origin which is tangential to the reaction rate curve R that determine the pattern drawn out by the dependence of the stationary state as etc. on residence time. Systems with relatively unstable catalysts are characterized by low values of T ~ .For these the minimum gradient of L, 4/i2, is relatively high so the flow line is always steep. If L always lies above the tangent line (which has a gradient of I/J, R and L cannot have any intersections other than at the origin. The stationary state yss= 0 is unique for all residence times. This dependence of ysson T,,, is shown in Figure 2b, where there are no nonzero solutions. Under conditions where the catalyst is more stable, r2 becomes higher. The flow line of minimum gradient may now lie below the tangent line, and there will then be three stationary-state intersections over some range of residence time. The dependence of yss on ires for this situation is shown in Figure 2d. The two nonzero extents of conversion y2 and y3 lie on a closed curve or isola. There is both a lower bound and an upper bound on the residence times over which these solutions exist. The two ends

There are two such residence times, given by the roots of this quadratic, because L passes twice through tangency as T , , ~is increased: once as the gradient is decreasing before reaching its minimum and once as the gradient is increasing after its minimum. The borderline between parts b and d of Figure 2 w u r s when the loss line with minimum gradient exactly coincides with the tangent line. The flow line L may then just touch the reaction rate curve when rrep= i2.For this residence time, there is a second This represents the “birth” stationary-state solution, yss = of an isola as a point in the yss-~res plane and is shown in Figure 2c. The two ends of an isola are points of extinction of the nonzero extent of conversion. One of these accompanies a decrease in residence time (as found for the S-shaped curve in the previous section when the catalyst is infinitely stable), and the other accompanies an increase in T,,. 2. Stationary-State Patterns with Catalyst Znfow: Isolas and Mushrooms. If some catalyst B is included in the inflow, with a concentration Po = bo/ao,the relationship among CY,6, and y at the stationary state becomes

Substituting this result into 9, we obtain instead of eq 20 R

+

(1 - yss)(yss Po)’ =

I(1 + rres

):

2

yssp L (22)

We now no longer have yss = 0 as a solution. Again the effect of a nonzero concentration of autocatalyst in the inflow is to alter the reaction rate curve R, lifting it up in the neighborhood of the origin (Figure 3a). Thus, the lowest intersection corresponds to a nonzero extent of conversion. There are also now two lines passing through the origin which are tangential to some part of R: one corresponds to extinction and one to ignition. The flow line L still passes through a minimum gradient when rIes= i2.We must now consider the position of this line with respect to the two tangents (Figure 3b-f). If the catalyst is highly unstable, L is always steep and always lies above the upper extinction tangent. The only intersection is the one at low extents of conversion y,. This branch of solutions (Figure 3b) reaches a maximum when r,,, = i2. For the special case that the minimum slope of L corresponds to the upper tangent, an isola is born as a point (Figure 3c). With a slightly more stable autocatalyst, L may pass beneath the upper tangent but not reach the lower one. This leads to an isola (Figure 3d) which sits above the lowest (nonzero) branch of stationary states. The size of the kola increases as the catalyst lifetime becomes longer, as does the height of the maximum along the y, branch. A second special case arises when i2is such that L can just touch the lower ignition tangent. The bottom of the kola then just touches the lower branch of stationary states (Figure 3e). If the catalyst stability is further increased or if more catalyst is added to the inflow, the flow line may pass below the lower tangent. The lower two stationary states y1 and y2 then merge and disappear over some range of residence time. The resulting pattern (Figure 3f) is called a mushroom. There are now two extinction points and two points of ignition. Experimental observations of isola and mushroom dependences of the stationary-state extent of conversion on residence time have recently been made for the arsenite-iodate system by Ganapathisubramanian and Showalter.28 In order to provide a suitable removal step for the autocatalyst I-, these authors employed an additional flow of solvent through their reactor. (This is sufficient to enrich the pattern of stationary states but does not provide the

The Journal of Physical Chemistry, Vol. 89. No. 1 , 1985 27

Feature Article

equations (9) and two independent variables a and j3. This gives a total of four partial differential coefficients which determine the sign and character of the two exponents A+ and A- via the quadratic equation

" (3) aa d r

I i.e.

X2 - (Tr)X

+ det = 0

Depending on the signs and magnitudes of these partial derivatives, which are evaluated at the stationary state, A+ and Xmay be real or complex with positive or negative real parts. The time dependence of the perturbations is then described by

Figure 3. Cubic autocatalysis with catalyst of finite lifetime, nonzero inflow of B: (a) dependence of R and L on 7%;(b) unstable catalyst, no tangency of R and L; (c) birth of isola, minimum slope of flow line corresponds with upper (extinction) tangency; (d) isola, minimum slope of L lies between tangencies; (e) birth of mushroom, L touches lower tangent; (f) mushroom.

decoupling necessary to allow oscillations.) Good agreement between experiment and predictions based on a similar analysis to that described above was obtained by using known kinetic parameters. 3. Effects of Including the Uncatalyzed Step A B. Our simple model so far ignores the uncatalyzed conversion of A to B. We have in effect assumed that this always proceeds very much more slowly than the catalyzed step. The main effects of including the additional step A B rate = k,a

-

-

are qualitatively very similar to having a nonzero inlet concentration of the a u t o c a t a l y ~ t .Again ~ ~ ~ the lowest branch of stationary-state solutions is lifted up above the y = 0 axis. There is, however, one important change. It occurs at long residence times. Here dependences of yy,on 7, shown in Figures 2 and 3 all have a unique solution, and this solution tends to zero conversion as 71a becomes infinitely long. This suggests that, in the absence of an uncatalyzed component, there would be no net reaction even in a closed vessel. Hence, the system would not approach its state of chemical equilibrium. This remains the case in the cstr even if the concentration of autocatalyst in the inflow is nonzero. If, however, the uncatalyzed reaction is included, the lowest stationary state does not fall back to zero. Instead, this branch tends to complete conversion yss= 1 as 7,- tends to infinity, as required by equilibrium considerations. Nonetheless, the isola and mushroom patterns remain. They do not disappear, because the behavior at all but the longest residence times is only very slightly influenced by including the uncatalyzed step. IV. Character of Stationary States for Systems with Catalyst Decay: Bifurcation Diagram in Parameter Space and Birth of Sustained Oscillations In section 11.4 we analyzed the stability of the different stationary states to small perturbations. The determining quantity turned out to be the sign of a single coefficient X which was obtained by differentiation of the single governing equation with respect to the single variable y. For a system in which the catalyst decays we may examine the stabilities in an analogous way.29 We now have two governing

Aa = p 1 exp(X'7)

+ p 2 exp(X-r)

AS = q1 exp(X'7)

+ q2 exp(X-7)

(24)

where p l , p2, q l , and q2 are constants. 1. Systems with No Catalyst in Inflow, Po = 0. When the concentration of B in the inflow is zero, the lowest stationary state corresponds to no conversion (a1= 1, = 0, y1 = 0). Both the roots X+ and A- are real and negative. This means that the exponential terms in (24) decay monotonically to zero as 7 tends to infinity. A useful, alternative way of looking at this approach to the stationary state is to plot the variations of a and j3 against each other. In these a+ phase-plane diagrams the stationary states are represented by points, known as singular points or singularities, with coordinates a,,and p,,. There may be one or three of these depending on the number of solutions of eq 20. The variations of a and p in time now correspond to motions or trajectories across the plane. The monotonic decay in time to aI and PI, governed by real, negative values of A+ and X-, appears as a monotonic movement across the a-p phase plane to the point of a = y1 = 1, j3 = PI = 0. This monotonic approach characterizes the singularity as a stable node (Figure 4a). The middle solution (a2,p2, or y2) is always such that Xf and A- are both real and always have opposite signs. Although one of the experimental terms decays as T increases, the other, with a positive value for A, increases with time. The stationary state is unstable, and the system cannot reside there. The motion across the a-j3 phase plane corresponding to these values of X is shown in Figure 4b. All but two trajectories diverge from the singularity which is known as a saddle point. (The trajectories have the same configuration as the contours of the saddle displayed by the transition state on a potential energy surface in the theory of chemical reactions.) The sign and character of the roots A+ and A- for the highest stationary-state extent of conversion (y3,a3,f13) vary considerably. At the shortest residence times at which the isola exists, XI and X2 are real and both negative, so we have a stable node on the fast-flow side. As T, is increased, the roots may become complex conjugates with negative real parts. When this occurs, Aa and A0 decay via a series of damped oscillatory overshoots and undershoots. This sort of response characterizes a stable focus (sf), as the variation of j3 with respect to a appears as a spiral (Figure 4c). For longer residence times, the roots are still complex, but the real parts may change from negative to positive. The perturbations now grow in time, with a spiraling away from the stationary state in the a-/3 phase plane (Figure 4d). This is known as an unstable focus (uf). The roots finally become real but remain positive before the right-hand extinction of the isola. The perturbations now grow monotonically (exponentially) in time. The stationary state is an unstable node (un; Figure 4e). (29) Andronov, A. A,; Vitt, A. A,; Khakin, S. E. "Theory of Oscillators"; Pergamon: Oxford, 1966.

28

Gray and Scott

The Journal of Physical Chemistry, Vol. 89, No. I , 1985

I I ’

,

I

(r)

Figure 4. Different characters of singularities (stationary states) in a-,3 phase plane: (a) stable node, sn; (b) saddle point, sp; (c) stable focus, st; (d) unstable focus, uf; (e) unstable node, un.

2. Systems with Catalyst in the Inflow, Po > 0. Including some catalyst in the inflow not only leads to a wider range of possible patterns for the dependence of the stationary states on the residence time but to a wider range of characters displayed by the lowest of the three solutions. For small nonzero inlet concentrations of @, a1may change from a stable node to a stable focus for some range of 7ra. As Po gets larger, the lowest solution may also become unstable (uf or un). This instability is only found for mushroom patterns and on the branch of y1at long residence times. The y1 branch at short residence times (to the left of the mushroom in Figure 30, the so-called flow branch, is always stable. The middle stationary state yz is still always a saddle point. The uppermost branch shows the same variety of character and stability displayed when there is no catalyst in the inflow. 3. Dependence of Stationary-State Patterns and Character on the Parameters r2 and Po. What pattern is displayed by the dependence of the stationary-state concentrations on residence times (isola, mushroom, etc.) depends upon the dimensionless values of the catalyst lifetime 7 2 and the inlet concentration Po. We may subdivide each of these patterns according to whether their difference branches show changes from stable to unstable solutions as rre8is varied. Thus, we find nine different possible case^.^*^^^ These are presented in Figure Sa. Each of these corresponds to different ranges of values of r2and P,,. These ranges can be portrayed as areas on a 72+O plane as shown in Figure 5b. Thus, if Bo = 0.1 and 72 = 40, the system will display pattern viii-a mushroom which has unstable portions along both the uppermost branch y3and the lowest branch yl as well as the saddle point y2.

0.1

b

0

(30) Escher, C. E.; Ross, J., private communication. (31) Boissonade, J.; De Kepper, P. J . Phys. Chem. 1980, 84, 501.

1,15

0.1

Figure 5. (a) Nine possible patterns for the dependence of the stationary-state extent of conversion on residence time: solid lines represent stable solutions (sn or sf), dashed lines give unstable states (un,uf, or sp), and the symbol represents point of Hopf bifurcation; the middle branch is always a saddle point. (b) Values of catalyst inflow concentration Po and lifetime T~ at which the above patterns are found.

that developed for nonisothermal reactions in the chemical engineering literature. 1 . Birth of Limit Cycles. If the concentrations a and P of the reactant A and autocatalyst B do oscillate regularly, then the variation of (Y with respect to P in the phase plane appears as a closed loop or limit cycle. (This is most easily seen by drawing out the changes in (Y and /3 with time as the system moves round the limit cycle in Figure 6a in the direction of the arrows. The resulting concentration histories are as shown in Figure 6b.) Such limit cycles originate at points where a stable focus changes to an unstable focus. This is known as a point of Hopf bifurcation.3293 The mathematical condition for this is that the real parts of the roots X+ and X- (see section IV) pass through zero. This provides an analytical criterion which can readily be applied. When there is no autocatalyst in the inflow (Po = 0), the residence time at which Hopf bifurcation occurs is given by

V. Sustained Oscillations in Autocatalytic Reactions: Limit Cycles Systems that have a stable focal character show a damped oscillatory approach (or “ringing”) to their final stationary state. There are a number of chemical reactions, however, which display long-lived oscillatory sequences that become apparently indefinitely sustained when carried out in open systems. Some previous discussions of such oscillatory behavior have been based on the ”cross-shaped” diagrams of De K e ~ p e r ,and ~ ] some are presented in terms of the switching between two different branches of the S-shaped curves. The present model allows a clear and rigorous analysis to be employed without either oversimiplification or loss of physical significance. This approach mirrors

0.05

7,,

=

T*,,~

N

=

j/4723/2(1

+ (1 - 472-1/2)1/2)2

723/2(1- 2 7 ~ - ’ -/ ~272-l - ...)

There is no such change possible for systems in which the half-life of the catalyst is too short, such that 72 < 16 (Le. 7 2 ‘ 1 2 < 4). The initial growth in size of the limit cycle is proportional to the square root of the difference between T~~ and 7*,,. The complete analytical expression can also be derived,I8at the expense of some algebraic effort, for our simple model. The “amplitude” of the corresponding oscillations also grows in this way and can (32) Hopf, E. Ber. Verh. Saechs. Akad. Wiss. Leipzig,Math.-Phys. KI. 1942, 94, 3.

(33) Hassard, B. D.; Kazarinoff, N. D.; Wan, Y.-H. “Theory and Applications of Hopf Bifurcation”; Cambridge University Press: Cambridge, England, 198 1.

Feature Article

The Journal of Physical Chemistry, Vol. 89, No. 1 , 1985 29

Figure 6. (a) Limit cycle in a+ phase plane. (b) Sustained oscillations corresponding to motion around the limit cycle.

be represented as a projection onto the response (cu,)-constraint ( T , ~ diagram ) (Figure 7) as an envelope originating at T*,. This representation is extremely compact and conveys a livelier feeling for the growth of the oscillations than that given by Uppal et al.;4 they display only the upper limit of the extent of conversion or the switching between branches of stationary states, which would require a different diagram for each residence time. 2. Stability of Limit Cycles. The appearance of a limit cycle in the a+ phase plane does not guarantee sustained oscillations in the concentration histories. Just as stationary states may be stable or unstable, so may limit cycles. An unstable cycle does not lead to physically realizable oscillations. Instead, it acts as a “watershed” in much the same way as the separatrices of a saddle point influence the motion across the phase plane. Any point starting within an unstable limit cycle tends to the stable stationary state, which also lies inside. A system with initial conditions outside the cycle must move to a different stationary state. Thus, whenever unstable limit cycles arise, their size is important. Stable limit cycles correspond to physically observable sustained oscillations. Their sizes are representations of the oscillatory amplitude, although they give no indication of the period. For our system with no B in the inflow, Hopf bifurcations occur always at some residence time so long as 72

2 16

(25)

For the limit cycles to be stable, however, we require that r2 > ca. 28.5 The frequency wo of the oscillations at the birth of the stable limit cycle is finite and nonzero. At the point of Hopf bifurcation we have purely imaginary eigenvalues for the stationary state A’ = f i w o

3. Computed Results. The general route to evaluating the size and period of a stable limit cycle for given values of the residence time etc. is to integrate the parent equations (9) starting at a suitable initial point. A series of such “numerical experiments” establishes the dependence of the oscillatory behavior on these different parameters. In the vicinity of the critical residence time T*,,, the concentrations vary sinusoidally with relatively high frequency and small amplitude. As the residence time is increased, the pulses become more angular and are separated by longer periods. 4. Phase Portraits. The numerical computations also allow us to complete the phase portraits around the stationary states. Different combinations of (i) the number of solutions and (ii) their character and stability along with (iii) the occurrence of stable or unstable limit cycles give rise to the different patterns. In all, a total of nine rearrangements are displayed by the cubic autocatalysis with first-order decay of B. These are shown in Figure 8. Stable limit cycles surround either unstable stationary states or unstable limit cycles; an unstable limit cycle must enclose a stable stationary state. Any system will show changes between the different portraits as T ~T ,~or, Bo is varied. For instance, the appearance of sustained oscillations about the highest extent of conversion accompanies a change from pattern b to portrait e (Figure 8) as an Hopf bifurcation occurs. 5 . Extinction of Oscillations. A stable limit cycle may dis-

Figure 7. Dependence of a, on T~~ for a system with = 40 and Po = 1/15, showing variation in character and stability of stationary states and

the amplitude of sustained oscillations about the unique unstable focus. The symbol 0 represents point of Hopf bifurcation.

Figure 8. Different phase portraits in a-P plane corresponding to different combinations of stable and unstable stationary states and limit cycles: (O), sn or sf; ( X ) , sp; (O), un or uf. Stable limit cycles are shown by solid curves, and unstable ones by dashed curves.

appear in a number of different ways. If the residence time is decreased back through the point of Hopf bifurcations (P,), the limit cycle shrinks back to zero size. The corresponding oscillations decrease smoothly to infinitesimally small amplitudes as the unstable focus becomes stable. There are three possible fates for the cycle as T~~ is increased. First, the oscillations may attain a maximum amplitude and then decrease in size again as the residence time approaches a second point of Hopf bifurcation. The two other routes lead to the sudden quenching or extinction of fully developed oscillatory behavior. Either the stable limit cycle “collides” with an unstable cycle growing from inside and surrounding the same stationary state (see e.g. Figure 8) or it may grow sufficiently large to touch the separatrix of a neighboring saddle point. In the latter case a homoclinic orbit is created exactly at the point of touching.34a Both of these mechanisms lead to a ”bursting” of the stable cycle and a dramatic change in the behavior of the system. mas elk^^^^ has provided an interesting commentary on some of the bifurcations of isothermal systems.

VI. Aperiodic Behavior in Autocatalytic Systems: The Springs of Chaos The simple, two-variable scheme studied so far leads to regular oscillatory behavior. The oscillations have constant amplitude and frequency for given values of 72, T,,,, and Bo. The majority of experimental observations of periodic behavior fall into this class. However, under some circumstances, more complex patterns can be found. Orban and E p ~ t e i nhave ~ ~ discovered transitions (or bifurcations) from apparently regular, large oscillatory excursions to ones with one large peak followed by one small peak and so on, in the chlorite-thiosulfate reaction. Subsequent transitions to trains with yet smaller pulses interspersed between (34) (a) Merkin, J. H.; Needham, D. J. Proc. R. SOC.London, Ser. A, in press. (b) Maselko, Chem. Phys. 1982, 67, 17. (35) Orban, M.; Epstein, I. R.J. Phys. Chem. 1982, 86, 3907.

30 The Journal of Physical Chemistry, Vol, 89, No. I , 1985

Gray and Scott

the large excursions occur as the residence time is increased. Furthermore, certain experimental conditions give rise to nonsteady behavior which has no apparent periodicity. Such aperiodic or chaotic states have also been observed for the Belousov-Zhabotinskii reaction by various authors.36 Chaos requires that a system has at least three independent variables. A purely homogeneous, isothermal chemical scheme would require three independent concentrations. A very readable introductory study of a chaotic response for a nonisothermal system with thermal feedback having two independent reactions coupled by temperature variations has recently been presented by Jorgensen and A r k 3 ’ However, is is still an open question whether chaos arising from an homogeneous chemical mechanism has been observed experimentally or whether it comes from the imperfect control of external features, such as small variations in the ambient temperature of the flow rate. Let us consider that the use of a peristaltic pump might superimpose a small periodic perturbation on the flow rate, so that the residence time is given by

The determinant vanishes at the turning points of the S-shaped dependence of the stationary-state concentrations on T:,. Thus, eq 30 can only be satisfied if Po < I / * . In the limit of zero inlet concentration of the autocatalyst, eq 29-31 arc satisfied by

= rEs/(l

E COS

ires

2xw7)

where t is a small number and 7ES is the mean value of T,,,. Substituting this into eq 9 has the effect of introducing the time explicitly into the right-hand side of the governing equations

dB - cup P +Po - P _ --

+

dT

7:

4 P o - P ) cos 27rwt

(26)

K

=2

These equations are called nonautonomous. It is possible to show that for certain values of the various parameters 7 2 , 7Es, Po, t, and w this system will describe aperiodic behavior. We first examine the autonomous equations obtained in the limit of vanishingly small e, t 0. These are the same as eq 9; however, it will be convenient here to study the case where 72 is not constant but varies with the residence time in such a way that

-

-

is a constant. Multiplying through (26) by 7Es and introducing 0 the new time variable t ‘ = 1/72, = t/tres, we find as E da = -7:,ap* (1 - a)

~

S

= S 1

Bo - KPs,

we require Tr(J) = ( 1

+ 2K)pS,2- 2(1 + Po)Pss+ ( I + .)/.E,

=0

(29)

where p,, is given by the solution of

(36) See e.g.: Vidal, C. In “Chaos and Order in Nature”; Haken, H., Ed.; Springer-Verlag: New York, 1981. Swinney, H . F. Phys. Rea. Leu. 1982, 49, 245. (37) Jorgensen, D. V.; Aris, R. Chem. Eng. Sci. 1983, 38, 45 (38) Keener, J. P. SIAM J . Appl. Math. 1982, 67. 25.

P,, = y4 7Es

and r2 are equal,

Again we find that a characteristic half-life of the autocatalyst such that 16k2 = klaoZ is of great significance. In section I11 we have shown that this is the minimum value of 72for which isolas arise in the 7ss-7res plane. It is also the threshold value of 7 2 for which Hopf bifurcation can occur (section V), and now we see that it is the minimum value for chaotic behavior.39 This novel result is clearly of great interest.

VII. Quadratic Autocatalysis (Simple Branched-Chain Reactions) with Complex Termination We have so far discussed exclusively the cubic autocatalysis (1b). The simpler rate law (la) with a quadratic dependence on concentration may be more widespread: it is certainly also of interest and of importance. The step A+B-2B has already found a p p l i ~ a t i o n l ~inq ~elementary ~ studies of isothermal branched-chain reactions and also in many biochemical or biological fields. When this rate law is coupled with a first-order decay (4) of the autocatalyst B for systems with no catalyst inflow (Po = 0), there may be two stationary-state solution^.'^ One of these corresponds to no conversion, ass= 1 and p,, = 0, and exists for all residence times. The other solution corresponds to nonzero extents of conversion and exists for residence times such that 7res

= klaOtres

> 72/(72

- 1)

(32)

where

+

The dependence of yss = 1 - asson 72, governed by (28) is an S-shaped curve provided Po is less than I/*. Mushrooms and isolas are not found when 7 Z s and r2 are linked by eq 27. The conditions under which the perturbed, nonautonomous equation (26) with t # 0 will have chaotic solutions are those for which the trace and the determinant of the Jacobian matrix of the autonomous system (28) become zero simultane~usly.~~ Thus, noting that

= 16 and

From eq 27 we see that K = 2 implies that hence 7 2 = 16

7:s

7:,

7Es

72

= k,ao/k2

is the appropriate form of the characteristic catalyst lifetime. When both these stationary states are realistic, i.e. when ( 3 2 ) is satisfied, the upper nonzero state is stable (either a stable node or a stable focus (see section IV)) and the lower state with zero conversion is an unstable saddle point. For short residence times, such that inequality 32 is reversed, the zero solution becomes stable and the physically unrealistic nonzero state (which now corresponds to negative catalyst concentrations, p,, < 0) becomes a saddle point. The effect of including some catalyst in the inflow (Po > 0) is to remove the bistability found above. A unique nonzero stationary-state extent of conversion exists for all residence times. This solution is always stable, although it may change character from node to focus as the residence time is varied. Irrespective of the inlet concentration of the autocatalyst (Po L 0), there is no Hopf bifurcation with this rate law and first-order removal of B; hence, this system cannot display sustained oscillations or stable limit cycles. Similarly, there is no possibility of a transition to chaotic behavior even if the flow rate is perturbed periodically. We shall now consider the effects of coupling a quadratic autocatalytic step with a removal or deactivation of the catalyst B whose rate increases less quickly than linearly with the catalyst concentration B C rate = k2b/(l + rb) (33)

-

(39) Gray, P.; Knapp, D.; Scott, S. K., paper in preparation.

The Journal of Physical Chemistry. Vol. 89, No. 1, 1985 31

Feature Article Such a rate law may arise for instance if the conversion of B to C occurs on a surface where the number of free sites available for adsorption may become significantly reduced (LangmuirHinshelwood kineticsa) or if the conversion is catalyzed by an enzyme which might be immobilized within the reactor (Michaelis-Menten kinetics41). This modification introduces sufficient nonlinearity to the system to allow a stable focus to become ~ n s t a b l eand ~ ~ for ,~~ sustained oscillatory behavior to grow from such a point of Hopf bifurcation. In dimensionless form, the mass-balance equations for this scheme are da _ - -4 -(11 a) d7 71,

+

+

where p = rag is the dimensionless measure of the saturation term in (33). For simplicity we shall concentrate on systems with no catalyst in the inflow, Po = 0. Equation 34 then has three stationary-state solutions. One of these is the state of no conversion, am = 1 and P,, = 0. The other two solutions are given by the roots of a quadratic equation, one of which corresponds to positive concentrations of autocatalyst and the other to negative values of &. Only the positive root is physically acceptable. Again we find that this solution does not exist for all residence times, only those for which inequality 32 is satisfied. This requirement on T , is~ determined only by the catalyst lifetime 72 and is independent of the size of p . The dependence of ym = 1 - asson residence time is shown in Figure 9. The stabilities of the stationary states along the two branches are also indicated. The state of zero conversion is a stable node (sn) at the shortest residence times but changes to a saddle point as the second, nonzero branch emerges. The upper solution is, at first, a stable node. As the residence time is increased, there are changes in character between stable node and stable focus (sf). Ultimately, there is a Hopf bifurcation at a critical residence ~ a stable limit cycle emerges to surround an unstable time T * and focus. Sustained oscillations are then observed in the concentrations of A and B computed by integrating eq 34. Hopf bifurcation is only possible for systems in which the dimensionless value of the saturation term is greater than the characteristic lifetime of the catalyst, Le. for which P

> 72

The bifurcation is always supercritical; Le., it gives rise to stable limit cycles and sustained oscillations whose amplitudes grow smoothly from zero as the residence time is increased beyond T * , ~ The dependence of the critical residence time T*= on the saturation term p is shown in Figure 10 for various values of T ~ In . general, we find a peninsula within which oscillatory behavior is found. These peninsulas move to longer residence times as the catalyst becomes more stable, Le. as T~ lengthens. The point of Hopf bifurcation is given a p p r o ~ i m a t e l yby~ ~

VIII. Discussion and Conclusions The results presented in the previous sections demonstrate clearly that complex kinetics and nonisothermal effects are not necessary requirements for a system to display a vast range of (40) Langmuir, I. J . Am. Chem. SOC.1918, 40, 1361. (41)Michaelis, L.; Menten, M. Biochem. Z . 1913,49, 333. (42) Merkin, J. H.; Needham, D. J.; Scott, S.K., submitted for publication in Proc. R . SOC.London, Ser. A .

$* ‘J

Tres

Figure 9. Quadratic autocatalysis with surface termination of catalyst: dependence of stationary-stateextent of conversion and its character on residence time. Hopf bifurcation occurs at leading to sustained oscillations and stable limit cycles. The symbol 0 represents point of

Hopf bifurcation.

~

(35)

Although this expression is strictly valid only in the limit of large T~ (and hence, from (35), of large p), it is good to better than 5% accuracy over much of the range of p for r 2 as low as 4.

1

I 0

1

I

‘Ogl,&P 3

Figure 10. Dependence of residence time for Hopf bifurcation on saturation term p for various values of catalyst lifetime r2. The dashed line is the approximate bifurcation locus in the text.

complex and interesting (perhaps even bewildering) behavior. The analyses require a certain amount of algebraic manipulation, but the use of flow diagrams in particular helps illustrate many of the stages much more than any purely mathematical approach. We conclude this article with a list of discussion points which will draw out the major conclusions and indicate areas of future research. (1) Kinetically elaborate models are not necessary to model extremely varied chemical behavior. Even the simplest situation, a single-step reaction (Id), may lead to two different patterns for the dependence of the stationary-state extent of conversion on residence time. One finds either unique solutions or multiplicity, with associated ignition, extinction (“washout”), and hysteresis. (2) When the autocatalyst B is not infinitely stable but decays in a first-order manner (k2 # 0) there may be three patterns of stationary states vs. T,,: unique solutions, isolas, and mushrooms. Isolas have two points of extinction but no ignition; mushrooms have two extinctions and two ignitions, with two regions of hysteresis. (3) Considerations of the local, dynamic stability further subdivide these patterns. There are two forms of unique solutions, differing as to whether the stationary state may become unstable or not, three types of isola, and four varieties of mushroom. These different patterns divide the P0-q plane into nine regions. (4) The points at which a stable focus becomes unstable are the same as those at which limit cycles emerge in the phase plane. Stable limit cycles correspond to conditions under which sustained oscillations are observed in the species concentrations. Unstable limit cycles may also arise. These play a similar role to the separatrices of a saddle point acting as a frontier around a stable stationary state. Initial conditions starting within the cycle tend

32 The Journal of Physical Chemistry, Vol. 89, No. 1, I985 toward that state, while those starting outside move to a different singularity. ( 5 ) The onset and disappearance of oscillations may arise in a number of ways as the conditions (e.g. the residence time) are changed. The most straightforward of these has a smooth growth of the limit cycle and, hence, of the oscillatory amplitude from infinitesimally small values. If the limit cycle becomes too large, it may touch the separatrix of a saddle point. The oscillations then cease abruptly and the system moves to the lower (stable) stationary state. If the residence time is then decreased again, the system may jump back into large oscillations as the stable solution merges with the saddle point. An alternative method by which oscillations are extinguished is the merging of a stable and unstable limit cycle which both surround the same stable stationary state. Decreasing the residence time reveals a range of conditions over which there is hysteresis between a steady reaction and an oscillatory mode. The system eventually jumps back into fully fledged oscillations. (6) All the above is achieved with the simple, irreversible autocatalytic step of Z e l ’ d ~ v i c hand ~ ~ a decay reaction. These schemes are not elementary steps but are adequate representations of some real systems over certain ranges of concentrations. To match the behavior over wider ranges, we may add additional reactions, but the principal interesting features still survive. The simplest embellishment2” of the mechanism is to invoke the lower order, quadratic autocatalysis (IC) and the uncatalyzed process ( 2 ) simultaneously with the cubic step (Id). We may also recognize that the reverse reactions may have nonzero rates. The conditions under which multiple stationary states may exist over some range of the residence time for this augmented scheme are readily established by considering the tangency of the R and L curves on a flow diagram. If we at first ignore the reverse reactions, R and L become tangential twice, corresponding to extinction and ignition, provided where k,, k,, and k, are the rate constants for the cubic, quadratic, and uncatalyzed steps, respectively. If the quadratic step is neglected (kq 0), multiplicity requires k,ao2 > 27ku (37)

-

while if the uncatalyzed step can be ignored, the condition becomes kcao > kq (38) The inclusion of these extra steps reduces the range of k, and a. over which multiple stationary states can exist, but it does not remove multiplicity completely. The reverse reactions also help to restrict the conditions for more than one solution, but again multiplicity can still be a feature of these autocatalytic systems. We may note that the principle of detailed balance requires that each of the steps IC, Id, and 2 have the same equilibrium constant K,, because they have the same stoichiometry. If kq is very small, then inequality 37 becomes for the reversible system kcaOZ> 27(1 K;’)2ku

+

( 7 ) The rigorous analyses employed in the present paper may lead to results which differ quantitatively and even qualitatively from the “pool chemical” approaches. In the latter, the concentrations of certain species, in particular the reactants, are assumed to vary so slowly compared with the time scale of the experiment or oscillations that they may be treated as constants. In this way, Karmann and Hinze have concluded incorrectly%that multiplicity is not a real feature of the system considered in ( 6 ) . B. F. Gray has shown4 that requirements on the magnitude of different rate (43) Zel’dovich, Ya. B. Zh. Tekh. Fiz. 1941, 11,493. Zel’dovich, Ya. B.; Zysin, Y. A. Zh. Tekh. Fiz. 1941, 11, 501. (44) (a) Gray, B. F.; Aarons, L. J. Symp. Faraday SOC.1974, No. 9, 129; J . Chem. Phys. 1976, 65, 2040. (b) Gray, B. F.; Morley Buchanan, T. J . Chem. SOC.,Faraday Trans. I , in press.

Gray and Scott constants for the pool chemical approximation to be valid for the “Bru~selator”’~ are inconsistent with the conditions for the onset of oscillatory behavior (Hopf bifurcation) but that the model of Field and NoyesI6 is not inconsistent. The use of flow systems and the machinery of sections 11-VI1 provide a much safer route to the analysis of chemical systems. (8) The kinetic schemes studied have displayed sustained oscillations only in open systems. In the limit of infinitely long m, it is a “common-sense” expectation that residence times, t , the model should go over smoothly to the required behavior of a closed system, i.e. that there should be a unique stationary state of stable nodal character which tends toward the chemical equilibrium composition. For the simplest models of sections 111-VI1 this is not, however, the case. The cubic rate law has a unique stable node at this limit, but this corresponds to a state of no conversion rather than the complete conversion of equilibrium. With the quadratic catalysis of section VII, the felony is further compounded: not only does the system tend to a state different from chemical equilibrium in this limit, it also retains its focal character. For the simple first-order decay of B, or with values of the saturation term such m: for larger that p < 72rthe solution is a stable focus as T~ values of p, the system is an unstable focus surrounded by a stable limit cycle. The inclusion of the lower order and uncatalyzed steps, even with the smallest nonzero rates, puts this right, however, and all the systems approach equilibrium as t, becomes infinite. (9) The complexity of the behavior exhibited by these schemes arises from the balance between the reaction term R and the flow term L. The cubic rate law leads to a highly nonlinear form for R which is typical of many real systems. Although this complexity in R may arise from a termolecular step, it might equally be the result of a subsequence of simpler processes. It cannot be emphasized too strongly that our approach does not require belief in a single elementary reaction. The cubic form for R requires only a linear dependence of L on the concentration to lead to the array of exotic patterns outlined in the previous section. If L can become a more complex function itself, then simpler dependences of R would suffice. (10) An interested reader will find many similar analyzes applied to a simple exothermic reaction in the chemical engineering l i t e r a t ~ r e . ~ ~Recent , * ~ studies include the discoveries of both new patterns in the dependence of stationary-state extents of conversion on residence time5 and new phase-plane portraits.46 Isothermal oscillations in biological systems have been demonstrated by Chance and co-authors4’ and discussed by various a ~ t h o r s . ~ ~ - ~ ~ Further investigationsof the cubic autocatalytic rate law of sections 111-V are being made by Lignola and D’Anr~a’~ using the catastrophe theory approach of Golubitsky and Keyfitz” as adapted by Balakotaiah and LUSS.~

-

-

Acknowledgment. We thank Professors B. F. Gray of Sydney, I. Epstein of Brandeis, C. Escher of Aachen, and particularly K. Showalter of West Virgina for stimulating criticism and comments. We are also grateful to the SERC and the EEC,to British Gas for research grant support, and to the ARGC for the award of a Queen Elizabeth I1 Fellowship to S.K.S. (45) Vaganov, D. A.; Samoilenko, N. G.;Abramov, V. G. Chem. Eng.Sci. 1978.33, 1282.

(46) (a) Williams, D. C.; Calo, J. M. AIChEJ. 1981,27, 514. (b) Kwong, V. K.; Tsotsis, T. T. AIChEJ. 1983, 29, 343.

(47) Chance, B.; Eastbrook, R. W.; Ghosh, A. Proc. Narl. Acad. Sci. U.S.A. 1964, 51, 1244. (48) Higgins, J. Ind. Eng. Chem. 1967, 59, 19. (49) Edelstein, B. B. J . Theor. Biol. 1970, 26, 227. (50) Seelig, F. F. J. Theor. Biol. 1971, 30, 485. (51) Sel’kov, E. E. Eur. J . Biochem. 1968, 4,19. (52) Lignola, P.-G.; D’Anna, A. Proc. R. SOC.London, Ser. A , in press. ( 5 3 ) Golubitsky, M.; Keyfitz, B. L. SIAM J . Marh. Anal. 1980, 11, 316.