Discrete Simulation Methods for Chemical Kinetics
2379
Discrete Simulation Methods for Chemical Kinetics Jack S. Turner Center for Statistical Mechanics and Thermodynamics, The University of Texas at Austin, Austin, Texas 78712 (Received July 11, 1977)
Simulation of chemical reactions is usually performed at the macroscopic level, and typically involves the time-stepwise integration of first-order differential equations. For an important class of chemical systems, however, the underlying deterministic theory may itself become inadequate, although it may still provide an indication of qualitative macroscopic behavior. In this way chemical instability phenomena were first predicted on simple chemical models, and deterministic methods continue to be important in this regard. Since fluctuations are neglected in deriving deterministic equations, however, to study the chemical instabilities themselves requires a more refined description. Just as a detailed molecular theory is out of the question, so also are complete microscopic simulations not feasible. In recent years, however, methods of theory and simulation have been developed at levels intermediate between molecular and macroscopic. Following a discussion of recent advances in the stochastic theory of chemical kinetics, discrete simulation methods at two such levels will be presented: At the microscopic level, hard-sphere dynamics is assumed for individual molecules in a mixture. Chemistry is introduced by assigning reactive transition probabilities at the instant of collision between reactants and altering molecular properties according to an assumed mechanism if the conditions for reaction are satisfied. In a higher level simulation, species populations rather than individual particles are followed, and a Monte Carlo procedure determines the sequence of events. The reactive molecular dynamics simulation is a computer experiment, and includes elastic as well as reactive collisions. In the second method, only events of direct interest occur, so that computational efficiency is gained at the expense of molecular detail. In addition, the stochastic simulation algorithm is rigorously derived from fundamental theory. The special features of each method will be illustrated in application to model chemical reactions.
I. Introduction As the program of this Symposium attests, the computer simulation of chemical reactions is usually performed at the macroscopic level, and typically involves the timestepwise integration of coupled first-order differential ,equations derived from a deterministic formulation of the underlying chemical kinetics. Macroscopic methods have long been of value in studies of complex reaction kinetics, usually to assess the validity of possible reaction mechanisms or to reduce to a global analytical form a wealth of experimental data.’ Moreover, with the development in recent years of numerical methods2 applicable to the “stiff’ differential equations usually encountered in models of real laboratory systems, it is now possible to perform quantitative simulations of highly complex chemical reit is the success of these new methods a c t i o n ~ .Indeed, ~ that is most convincingly documented in the proceedings of this Symposium. This success notwithstanding, there is an increasingly important area of chemistry for which deterministic methods are unable to provide an adequate description, often even qualitatively, of macroscopic phenomena. This area concerns chemical systems which may undergo nonequilibrium instabilities and transitions between different states of macroscopic chemical behavior. A variety of chemical instability phenomena have been discussed at this Symposium, ranging from instabilities in chemical steady states, leading to homogeneous chemical oscillations, for example, to combustion and explosive instabilities associated with states which are involved in a macroscopic chemical evolution. Many of the contributions in this area are concerned with the mechanistic detail and computer modeling of specific experimental systems. In contrast, the present paper deals with more general aspects of chemical instabilities, focusing primarily on the microscopic mechanism of instability common to all such phenomena, and on the development of appro-
priate methods of theory and computer simulation essential to their study. According to the assumptions on which classical chemical kinetics is based, purely deterministic methods are adequate as long as deviations from the macroscopic average values remain negligible. At thermodynamic equilibrium the validity of the macroscopic description is established via statistical mechanic^:^ e.g., if the average number of molecules is of order N, then fluctuations about that average will be of order N”’.Recalling that results of this type hold only away from points of macroscopic instability such as phase transitions or critical points, one may identify two types of situations in chemical systems for which fluctuations should become important: (a) The number of molecules N is small. This is the obvious case in which fluctuations of order N1/’ cease to be negligiblea5 Examples of this type are exceptional, however, as most of chemistry concerns systems containing a large number of molecules. Such macroscopic chemical systems have, typically, N 1020,and fluctuations of order N112 1O1O about such average values remain entirely negligible. In most chemical systems, therefore, fluctuations are unlikely to be important unless (b) the system is near an instability of the macroscopic rate equations. In this case the above “square-root rule” does not apply and, just as at equilibrium transitions, fluctuations may be amplified near chemical instabilities to produce observable effects. As these introductory remarks imply, a basic problem for equilibrium and nonequilibrium transitions alike is to understand the range of validity of the purely macroscopic description in the neighborhood of an instability. To approach this problem in either case requires a molecular description of a macroscopic system. Fortunately, a complete microscopic description is not essential in the case of chemical instabilities, and a great deal can be learned in a higher-level (stochastic) approach to chemical
-
-
The Journal of Physical Chemistry, Vol. 8 1, No. 25, 1977
2308
kinetics.6 In addition, macroscopic methods remain valid away from the transition region, and can provide valuable information about the nature of the instability itself. Indeed, chemical instabilities were first predicted and verified in studies of deterministic kinetic models,’ and macroscopic methods continue to be important in this regard.s The organization of this paper follows closely the title of this Symposium: Following an introduction to chemical instability phenomena is a survey of typical reaction mechanisms which may exhibit instability. Exploiting the analogy between equilibrium and nonequilibrium instabilites, I present next a classification of simple chemical models in the language of equilibrium phase transitions and critical phenomena. Theoretical methods appropriate to the study of fluctuations in chemical systems are then discussed, with emphasis on the implications for theory, simulation, and experiment of quite recent results in this area. Recognizing that many outstanding questions are beyond even approximate resolution with available theoretical methods, especially in complex chemical reactions of experimental importance, I present the case for computers and computer simulation of systems in which chemical instabilities may occur. Having reviewed some of the unique and often surprising properties of typical systems in this new area of chemistry (sections I1 and 111), and having formulated the problem of theory and simulation in this context (section 111),I discuss methods of discrete event simulation, at two levels of description, which appear especially promising for further theoretical and experimental understanding of the phenomena. These two methods are developed in sections IV and V, largely in connection with applications to specific chemical models, and a general discussion of the methods and results follows in the final section VI.
11. Chemical Instabilities, Dissipative Structures, and FluctuationsgJo Consider a chemical system under open or closed conditions so that it is possible to maintain a time-independent macroscopic state. At or near thermodynamic equilibrium, the macroscopic description of such a system yields a unique stationary solution that is always stable with respect to arbitrary perturbations. Sufficiently far from equilibrium, however, the governing kinetic equations for nonlinear chemical systems may admit solutions other than the continuous extension of the equilibrium solution into the nonequilibrium domain. Beyond a critical distance from equilibrium, moreover, this thermodynamic branch may become unstable in response to fluctuations occurring in the medium. Beyond the point of instability, the response of the system to an infinitesimal disturbance leads ultimately to a new operating regime characterized by organization in space, time, or function. Whatever the outcome of such a macroscopic instability, therefore, the instability itself originated in the response of the system to a fluctuation. Consequently a purely deterministic description of the system in terms of mean values alone is no longer adequate, and it is essential to supplement the “average” description with a theory of fluctuations extended to nonlinear systems under far-from-equilibrium conditions. At thermodynamic equilibrium fluctuations constitute a negligible correction to the macroscopic description of matter except near instability points such as phase transitions or critical points. Similarly, in nonequilibrium situations as well one expects that fluctuations become important only near points of nonequilibrium instability?Il0 For nonlinear systems in which such instabilities and The Journal of Physical Chemistry, Vol. 81, No. 25, 1977
Jack S . Turner
transitions can occur, the role of fluctuations is even more dramatic than at equilibrium, due to the existence in general of many distinct “phases” compatible with the external conditions imposed on the system. This fact has been known for several years now, since the first (numerical) demonstration on simple chemical models of multiple solutions to the governing macroscopic differential equation^.^ The variety of possible solutions includes time-independent homogeneous and inhomogeneous states as well as states which are organized in time (homogeneous chemical oscillations) or in time and space (travelling chemical waves). Which of these “dissipative s t r u c t ~ r e s ” ~will J ~ arise under given conditions depends crucially on the specific type of perturbation to which the original state is initially unstable. In terms of spontaneous fluctuations, which are present in general over a range of wavelengths and frequencies, both the final state and the evolution toward that state will be determined by the space-time characteristics of the fastest-growing unstable mode. In recent years developments in theory and experiment have enhanced enormously the need for a detailed understanding of fluctuations in nonequilibrium systems. On the phenomenological side is the general recognition of nonequilibrium instability and transition phenomena as key elements of the self-organization in biology of life processes at all levels from subcellular to organismic.1° Moreover, experimental evidence for cooperative behavior at the molecular level is now well documented in a variety of nonbiological chemical systems, most notably the Belousov-Zhabotinskii reaction,ll and in several subsystems of biological origin as well (e.g., the glycolytic cycle12). On the theoretical side, the recognition of the generality of nonequilibrium instabilities and the transition to dissipative structures has spurred greatly the mathematical study of nonlinear ordinary and partial differential equations of the types frequently encountered in models of physical, chemical, and biological systems. Of particular significance have been rapid advances in bifurcation theory,ls the application of which permits explicit construction of the bifurcating solutions of nonlinear differential equations in the vicinity of points of instability. Recent applications of these methods have revealed for the first time the incredible variety of types of nonlinear behavior which may arise in even quite simple chemical systems under realistic external con~traints.’~ For example, analytical studies of the so-called Brusselator model, a prototype for chemical instabilities, have revealed a wealth of solutions (spatial, temporal, and spatio-temporal) far exceeding that observed in earlier numerical analyses of the nonlinear differential equation^.^ In subsequent numerical verification and demonstration of these analytical predictions,14 the crucial role played by the specific type of perturbation which initiates the instability in the development of the final solution has been strikingly demonstrated.16 Despite the increasingly obvious importance of fluctuations in nonlinear, nonequilibrium systems, as indicated briefly here, this subject has been relatively ignored until recently. As a natural outgrowth of the macroscopic investigation of nonequilibrium instabilities and transitions to dissipative structures, Prigogine, Nicolis, and co-workers have begun a systematic study of nonequilibrium fluctuations.16 In the years since the initial investigation of fluctuations in nonequilibrium chemical systems,16a fluctuations a t nonequilibrium instabilities have been of increasing interest in areas such as quantum optics17and hydrodynamicsls as well as chemical reactions.lg From the
Discrete Simulation Methods for Chemical Kinetics beginning this study has demanded the development of new theoretical methods. Indeed, the problems of dealing with nonlinear rate processes and open system conditions combine to make both classical (equilibrium) fluctuation theoryz0and classical (equilibrium) nucleation theory,21in their present state of development, largely inapplicable to the analysis of critical fluctuations that lead to the emergence of dissipative structures. Recently, however, considerable progress has been made in developing new methods for the study of nonequilibrium fluctuations, in open systems, in a number of related aspects:1° e.g., via the two-cell (mean-field) analysis of long-range spatial correlations in the multicell case,23the study of systems in which the effect of external noise is taken explicitly into account,24and the theory of metastability and nucleation for systems which exhibit multiple macroscopic states.z5 Relevant aspects of this work will be summarized in the next section. As one would expect, these new developments in theory have employed chemical models chosen more for analytical simplicity than for any connection to actual chemical reactions. Due to the mechanistic complexity of even the simplest laboratory systems of interest in this study, moreover, application of even approximate methods to more realistic situations is a formidable task. A t the same time a detailed microscopic approach to any of the simple chemical models, in terms of nonequilibrium statistical mechanics, for example, is also not feasible. Out of a similar situation in the study of classical fluids emerged the method of “molecular dynamics”, the computer simulation of a laboratory experiment in which the “exact” classical dynamics of a collection of atoms or molecules is followed as they interact according to the laws of mechanimZ6 Averages taken over accumulated microscopic data yield “experimental” information (e.g., distribution and correlation functions, rate constants) for use in the direction and testing of theory. Quite recently, the molecular dynamics computer model has been modified to include inelastic or reactive scattering as well as the elastic processes of interest at equilibrium phase transitions,n and several applications of this “reactive” molecular dynamics (RMD) method to simple chemical models involving chemical instabilities have been r e p ~ r t e d . ~In‘ ~section ~~ IV a variation of the RMD method is developed in an application to homogeneous chemical instabilities involving multiple steady states. In addition, some results obtained in a study of local fluctuations and nucleation in this model are presented. Reactive molecular dynamics may be viewed as a microscopic computer experiment inasmuch as all relevant particle processes can be included to some degree. In particular the method treats elastic as well as inelastic or reactive collisions between particles and, indeed, requires that the latter be relatively rare events in order to simulate actual chemical kinetics in a realistic way. When large enough systems are considered (e.g., thousands of particles) it should be possible to “measure” temporal and spatial correlation functions, for example, and to make quantitative the notion of critical size and amplitude of fluctuations necessary to nucleate a macroscopic transition. A first step in the latter direction is reported in section IV. Applications to date of reactive molecular dynamics methods demonstrate feasibility for the study of hundreds or thousands of particles involved in chemical reactions for which reactive probabilities do not vary too widely. The latter condition is essential if statistically significant numbers of all possible events are to occur within a
2301
practical computation period. Even if the requirement for a large excess of elastic collisions is relaxed, however, reaction rates typical of experimental chemical systems such as the Belousov-Zhabotinskii reaction,ll for example, demand simulation run times which approach the feasible limit except for quite small numbers of particles. Turning therefore to a higher level method for this type of system, one may treat numbers of particles in a small cell or volume element rather than individual particles, and invoke a Monte Carlo procedure for determining which reactive event will occur, how much time will elapse between events, and in which “cell” of the system the process will occur. Diffusion between adjacent cells is handled in a similar fashion. In such a method, only the events of greatest interest (reaction or transport) are treated, so that widely different rate constants need not be a problem. In the change from reactive molecular dynamics to stochastic simulation, much larger systems may be considered, and computational efficiency is gained a t the expense of molecular detail. In addition, problems related to numerical “stiffness” associated with the deterministic differential equationsz9are completely avoided in a stochastic simulation. Finally, a number of implementations of such techniques have been presented which give qualitatively reasonable result^.^^^^^ Of particular significance is a recent formulation due to G i l l e ~ p i e in , ~ ~which the transition probabilities and Monte Carlo implementation are rigorous consequences of the very assumptions from which chemical master equations of the birth-and-death type are derived. This fact is of great importance for the development and testing of theory. Specifically, computer experiments based on the latter formulation simulate directly a stochastic chemical evolution satisfying a master equation of the type which is central to the theoretical development of section 111. A report of initial investigations of chemical instabilities using stochastic simulation methods is presented in section V.
111. Theory and Simulation of Chemical Instabilities In this section I discuss standard methods of theory which have proved useful in the study of chemical instabilities.32 No attempt is made to be exhaustive in this summary. Rather, I focus on the main ideas and results, illustrated by application to specific models, which bear most directly on outstanding problems in the field and which provide a strong motivation for the use of discrete simulation technique in attacking these problems. A. Deterministic Chemical Kinetics. Consider an open, multicomponent mixture of species S, among which a number of chemical reactions R, may occur. The system is assumed for simplicity to be maintained isothermal in the absence of external forces. If X,(r, t ) denotes the concentration of species S, at position r a t time t , then the mass balance equations for the system have the form
The first term gives the change in X,(r, t ) due to chemical reactions
Fn( ( X I )
FvnpWy
(3.2)
with v, the stoichiometric coefficient for species S, in reaction R, and w, the rate of reaction R,. The second term gives the change in X,(r, t ) due to diffusion, which is the only transport mechanism considered, and which obeys a linear phenomenological law of the Fick type, with The Journal of Physical Chemistry, Vol. 81, No. 25, 1977
Jack S.Turner
2382
diagonal diffusion matrix, in a mixture which is taken to be ideal. For chemical reactions, however, linear kinetic laws are seldom applicable, so that in systems of interest here the chemical terms F, will always be nonlinear through the dependence of the rates w,,on the composition variables X,. It is these nonlinearities which will be responsible for qualitatively new types of chemical behavior arising far from equilibrium. Let us adopt the set of reaction-diffusion equations (3.1) as the starting point for a deterministic description of chemical kinetics, and focus on the solutions of these equations under time-independent boundary conditions. Because of their nonlinearity, eq 3.1 may in general admit a number of time-independent solutions XO,(r),one of which is the so-called thermodynamic branch, the continuous extension away from equilibrium of the unique equilibrium state.33 In order to assess the stability of such steady states, one considers the response of the system to small perturbations about a particular reference state. If perturbed states in the immediate vicinity of the steady state X:(r) are written in the form
Xn(r, t ) = X : ( r ) + xn(r, t )
(3.3)
then, in the limit of small perturbations (i-e.,Ix,/XO,I 0: The steady state is unstable; the mode with w, > 0 will grow. In either case (1) or (2), the behavior of individual modes will be oscillatory if the corresponding wi # 0 and monotonic otherwise. (3) At least one w, = 0, all other w, 0: The steady state is marginally stable with respect to the mode having w, = 0. As a solution of the linearized equations that mode will neither decay nor grow, but may oscillate if it has in addition wi # 0. The marginally stable steady state of the nonlinear system may be stable or unstable depending on the explicit form of the nonlinear terms. To this point nothing has been said about the effect of diffusion on the stability of steady states. In general one would expect diffusion to be a stabilizing influence due to the tendency of transport to damp inhomogeneities. This tendency notwithstanding, the explicit presence of diffusion in the analysis permits the existence of spatially inhomogeneous perturbations as well, and therefore increases greatly the variety of situations which may arise once an instability occurs. This point will be explored The Journal of Physical Chemlstty, Vol. 81, No. 25, 1977
With A and B assumed to be in large excess (or buffered) so that they act as time-independent reservoirs, the determininstic rate equations for this system become (cf. eq 3.1 and 3.2)
ax
-=
at
k , A - k2BX + k 3 X 2 Y- k k y + D x -a 2x
ar2
,-
-.
ay 2X - = k 2 B X - k 3 X 2 Y+ D y L
at
ar
where r is the spatial coordinate for a one-dimensional system, as in the preceding subsection, and Dx,D y are Fick-type diffusion coefficients for intermediates X and Y, respectively. It is easy to see that this system admits a single spatially homogeneous steady state solution (Xo, Yo)
which is the continuous extension of the thermodynamic branch into the far-from-equilibrium region. Under periodic boundary conditions, a linear stability analysis along the lines of the preceding subsection yields the following characteristic equation for the normal mode frequencies (cf. eq 3.5): + [-k2B + k3X02 -I- k4 + k 2 ( D x + D y ) ] o + [ - k 2 B k 2 D y + (h3X02 + k 2 D y ) ( k 4+ k 2 D x ) ]= 0 (3.9)
It is a simple matter now to find the points of marginal stability (Le., w, = 0) separating regions of stability and
2383
Discrete Simulation Methods for Chemical Kinetics
instability of the steady state (Xo,Yo)with respect to the various normal modes. If B is taken as the single independent variable (for fixed A , Dx, Dy), then there are two values of B at which w, may vanish:
1 B’ =-[k3X02 t k4 t k 2 ( D x + D y ) ] k2
(3.10a)
and
B”
=
’(”” + k2 k2Dy
1) (k4 + k2.Dx)
(3.l o b )
At B = B’, the coefficient of w1 in eq 3.9 vanishes, and both roots are either purely imaginary, if the constant term (coefficient of wo) is positive, or real and nonzero if the constant term is negative. A t B = B”, the constant term is zero, so that one root vanishes identically, while the other is real and nonzero. There are then two types of situations which may arise, depending, for fixed A , Dx,Dy (and rate constants), on the wavelength of the disturbance: (a) B’< B”: For B < B’, the homogeneous steady state (3.8) is stable, the decay of perturbations becoming oscillatory as B B’, The point B’is an oscillatory marginal state. For B > B’, therefore, an initial perturbation will rotate around the steady state as it grows. At some point in B’ < B < B”, both complex roots (with w, > 0) become real. Beyond that point all initial perturbations would grow monotonically. Passage through the point B = B” would have no effect on the stability of the solution (3.8), as one root would remain positive and real. (b) B’> B”: In this case, the first point of instability, B”, corresponds to a nonoscillatory marginal state. For B > B”, both roots are real, and one is always positive. Therefore all initial perturbations grow monotonically in this region. (The marginal state at B’does not exist, since the constant term in eq 3.9 is negative once B > B”.) In summary, then, as B is increased continuously, the homogeneous thermodynamic branch (3.8) will experience an oscillatory instability at B’if B’< B”or a nonoscillatory instability a t B” if B’ > B”. In order to understand better the nature of the instabilities possible at these two points, consider the “critical” wavelengths for disturbances which will produce instability at the smallest value of E. Minimizing the expressions in eq 3.10 with respect to h yields
-
(3.11a) and
Bh’ =
[(
E ) l I 2 X ot
(2)3 112
2
(3.11b)
In an actual chemical mixture there will usually be present spontaneous fluctuations having a wide range of wavelengths. As B increases, then, the instability which appears will correspond to the smaller of the two critical values B,. For case (a) above, eq 3.11a imply that very long wavelength perturbations will grow and oscillate, the system remaining homogeneous. For case (b), however, the possibility of appearance of finite inhomogeneities is seen from eq 3.11b to involve a cooperative effect between chemical reactions and diffusion. Indeed, whenever the
rates of diffusion and chemical reactions are “comparable”, short wavelength inhomogeneities may grow. As diffusion rates become dominant, however, the instability occurs for increasing wavelengths, and the system may again remain nearly homogeneous. Which type of instability will occur is seen from eq 3.11 to depend on the ratio of the diffusion coefficients. In fact, the possible emergence of a spatial inhomogeneity requires BF< B6, or from eq 3.11
(gy’2
B”(Ai9 in eq 3.101. The evolution of the system beyond this instability has been studied numerically.’ The most striking conclusion of the various computer “experiments” on the scheme (3.6) is that, for B > B$ there exists a new stable steady state which is spatially inhom o g e n e o u ~ . ~In ~ ~general, J~ of course, the length of the system will not be an integral multiple of the critical wavelength given by eq 3.11b. This leads to an indeterminacy in the nature of the spatial dissipative structure which bifurcates from the homogeneous steady state at B$ The system will select that wavelength, compatible with the boundary conditions, which lies nearest the wavelength of the particular disturbance which destabilizes the initial steady state. Once one or another dissipative structure is established, however, it is then stable with respect to small fluctuations, and hence in its wavelength possesses a primitive “memory” of the fluctuation from which it originated. See Figure 8a and 8b of section V for the results of a stochastic simulation. Other than the two simplest possibilities treated here, the Brusselator also admits more complicated situations, including time-independent, oscillatory, and, travelling wave solutions of the governing macroscopic equations. The possible emergence of many types of dissipative structures has been discussed in detail by Nicolis et aL10J4 The Journal of Physical Chemistry, Vol. 81, No. 25, 1977
Jack S. Turner
For fixed and no-flux boundary conditions, in particular, analytical construction of the bifurcating solutions has been carried out,14 predicting a surprising variety of instability phenomena. When nonuniform distribution of the reservoir variable A is permitted, moreover, spatially localized structures of several types have been realized. The implications of these predictions for self-organizing systems are discussed in ref 10, as are the numerical (deterministic) sir nu la ti on^^^ which verify the theoretical predictions. 2. Multiple Steady States. A Cooperative Isomerization Reaction. Turning now to instabilities which do not involve symmetry breaking in the sense of the preceding discussion, let us consider systems which may possess multiple homogeneous steady states. A convenient model is provided by a reversible isomerization reaction in a macroscopically homogeneous system (3.13)
having activation energies t for the forward reaction and ‘ r / Nper molecule of B present in the system for the reverse reaction, where N is the total number of molecules. The activation energy for the reverse reaction depends on the number N B of B molecules present, and therefore introduces a cooperativity into the mechanism. In physical terms, this cooperativity expresses the stabilizing influence on a B isomer of intermolecular interactions with other isomers of the same type. Aside from being a convenient way to introduce nonlinearity into the simplest chemical model, this procedure yields a model which is isomorphic to a Bragg-Williams (mean field) model for monomolecular adsorption on a uniform surface. The latter model of the equilibrium two-dimensional lattice gas is known to exhibit a first-order phase transition36and provides a convenient representation for an especially intuitive presentation of the stochastic theory of metastability and nucleation a t chemical first-order transition^.^^ Relevant conclusions of this work which are of a general nature will be reviewed later in this section. If the cooperativity is introduced into the rate constant for the reverse reaction via a “mean-field’’ dependence on NB, then the equilibrium constant for the isomerization reaction takes the form K = exp(-c v N B / N ) (3.14)
+
where the activation energies are expressed in units of kBT, with hB the Boltzmann constant and the T the absolute temperature. If N is fixed (closed system), then the kinetics of the reaction is given by a single macroscopic rate equation
m B / d t = ho(N- N B ) e - e - hoNBe-qNB” Defining a new time variable 7 fraction x = X B N B / N gives d x / d r = 1 - x - xee-qx
3
(3.15)
koe-‘t and the mole (3.16)
A linear stability analysis (e.g., section 1II.A) of the steady state xo of this system implies stability whenever the normal mode frequency w is negative, and instability otherwise, where xw = -1 vx - gx2. Points of marginal stability correspond to w = 0 and, according to the quadratic expression for W , there may be 0, 1, or 2 such points. It is easy to see that there is no instability for 7 < vc = 4. For 7 > vC,however, there is a region of E in which three steady states are possible. Typical steady state curves are shown in Figure 1. Here the curve passing through the critical point (vc, tC, x,) = (4,2,112) separates
+
The Journal of Physical Chemistry. Vol. 81, No. 25, 1977
1 .oo .80 .60 m
x
.40 .20
0.00 0.00
. 8 0 1 . 6 0 2 . 4 0 3 . 2 0 4.00
6 Flgure 1. Steady states of the cooperative isomerization model, showing subcritical (7 = l), critical (7 = 4), and supercritical (7 = 5) curves of mole fraction xBas a function of the forward activation energy t . The deterministic transitions for 7 = 5 are indicated by arrows; the dashed line denotes an “equal areas” construction which determines the unique equilibrium transition.
regions of one and three steady states. In the three state region, the upper and lower branches are stable, the middle branch unstable, according to linear stability analysis. If it is supposed for convenience that E is an externally variable parameter of the model, then for 7 > vc the transition between stable branches of the deterministic solution should occur discontinuously at the points of marginal stability (arrows, Figure 1). The similarity between this picture and the pressure-volume diagram for the equilibrium liquid-vapor transition is apparent. Indeed, according to the macroscopic van der Waals model of a nonideal gas,36the liquid-vapor phase change in a homogeneous fluid should occur at a lower pressure than the reverse transition. That a single equilibrium transition pressure (e.g., dashed line, Figure 1, and “equal-areas’’ construction) is observed regardless of the direction of change is a direct consequence of the response to spatially localized, finite-amplitude fluctuations in the initial fluid state. Moreover, when such fluctuations are of insufficient size or amplitude to provide a stable nucleus for the transition, then states beyond the equilibrium transition point may be realized. Such states are termed metastable, since they are stable only until a large enough fluctuation appears to form a “supercritical” nucleus. At which point along the metastable branch will destabilizing fluctuations appear, and how long will the metastable state remain “stable” prior to occurrence of critical fluctuations? Clearly such questions are beyond the scope of the deterministic description, and will be addressed from the viewpoints of theory and simulation following development of the appropriate methodology. C. Chemical Instabilities as Nonequilibrium Phase Transitions. In discussing the transitions associated with nonequilibrium instabilities, it is both convenient and instructive to adopt the language of equilibrium phase transitions and critical p h e n ~ m e n a .From ~ ~ this point of view a qualitative classification of chemical instabilities may be based simply on linear stability analyses as outlined in the preceding ~ u b s e c t i o n . ~Such ~ a classification is presented in Table I, the dgtails of which are discussed in the following paragraphs. As will be seen in the next subsection, moreover, the analogy is not limited to the macroscopic level, and has a far deeper meaning in terms of the microscopic dynamics of fluctuations near a threshold of instability. Invoking the formal analogy between equilibrium and nonequilibrium instabilities, therefore, and referring to both model predictions and
Discrete Simulation Methods for Chemical Kinetics
2385
TABLE I: Classification of Chemical Instabilities as Phase Transitionsa Number of Variables
Order of Kinetics
Model
Order Parameter
Final State
2
concentration
homogeneous steady state
amplitude
inhomogeneous steady state
Bifurcation and Transition Diagrams
Second-order transitions *Schlogl I1
[All
I
1
2
2
Brusselator [IIIBI]
3
3
3 2
Oregonator [D] Brusselator [IIIBll
2
5
3
Oregonator
2
6
2
Brusselator [IIIBI]
3
7
3
Oregonator
2
4
First-order transitions I . 8
I
[D] [D]
d A - A’
homogeneous o scillatio n s
3
not symmetry-breaking
SchlGgl I
[A21
3
9
I
Isomer
[IIIBZ]
**
10
2
Enzyme
[B]
2
concentration
homogeneous steady state
B“
First-order transitions 11, symmetry-brea king I1
1
*Enzyme (Modified)
[Cl [D]
2
concentration
homogeneous steady state
2
amplitude
homogeneous o scilla tio n s
12
3
Oregonator
13
2
Brusselator [IIIBI]
3
inhomogeneous steady state
14
2
Brusselator [IIIBll
3
travelling w a v e s
1
B - B’
(b) ”*,
----
d--
i
a Brackets following model name indicate location in this paper of additional details and references. An asterisk denotes Appendix. A double asterisk indicates that the isomer model has irreversible transition (or extinction): see appropriate .. . exponential, not polynomial, kinetics.
experimental observations, let us now review briefly some of the nonlinear chemical phenomena which may arise far from equilibrium: 1. Homogeneous Steady States. The occurrence of multiple homogeneous steady states, an obvious consequence of nonlinear kinetics, is predicted for a number of simple models (entries 8, 9, and 10 of Table I). The transitions between such states provide the most direct analogue of first-order (discontinuous) behavior. Theoretical and experimental examples of such situations are well-known in open chemical systems.40 Observations of transitions between multiple steady states have also been reported by Creel and R o s s , ~who ~ performed closed system studies, under laser illumination, of the photochemically controlled equilibrium between NOz and its dimer. In the best known example of chemical instabilities, the Belousov-Zhabotinskii (BZ) reaction, the existence of multiple steady states has not been established yet by either theory or experiment. In this versatile chemical system, however, most other types of instability phenomena have been observed or predicted;ll in a very real sense, the BZ reaction has become the experimental prototype for chemical instabilities. 2. Homogeneous Oscillations. a. Second-Order (Continuous) Transition. When the homogeneous steady state (hss) solution of the chemical rate equations first becomes unstable as the distance from equilibrium is increased, the simplest type of homogeneous instability which may occur is a Hopf bif~rcati0n.l~ At a critical value of a bifurcation parameter (e.g., p a measure of reaction affinity), the hss, previously a stable focus, bifurcates into an unstable focus and a stable limit cycle. Beyond the bifurcation (Le., critical) point, quasi-harmonic bulk oscillations increase continuously from zero amplitude, eventually becoming nonlinear as p increases further. In
the language of phase transitions, this situation, called a “soft self-excitation”, is analogous to a second-order transition, an appropriate order parameter being the amplitude of oscillation. Experimentally this type of instability has been observed only recently in BZ-like reactions,42a convenient theoretical model being the socalled Oregonator of Field and Noyes.lld Examples are entries 4 and 5 of Table I. b. First-Order Transition. In closed-system studies of the BZ reaction, both the onset and the end of oscillatory behavior have been found to be abrupt rather than smooth, corresponding to a discontinuous transition between the hss and a large-amplitude, highly nonlinear (i.e., nonharmonic) limit cycle oscillation. These observations suggest a different type of bifurcation, called a “hard self-excitation”, which is analogous to a first-order phase transition. In this case both the hss and the finite amplitude limit cycle are stable for a range of values of the bifurcation parameter, the actual point of transition between them depending crucially on the nature of the fluctuations occurring in the system. As at equilibrium first-order transitions, one expects that here a nucleation mechanism involving finite localized fluctuations is responsible for a subcritical transition to the oscillatory regime. In addition, there exists the possibility of a hysteresis in the reversible transition between these two states of different temporal symmetry. This effect has been predicted from model studies,llg~~~ and suggestive experimental observations have been reported.l1gSs0 c. Critical Phenomena. In the Field-Noyes (FN) models of the BZ reaction,lldaf it is possible to pass continuously between regions of second- and first-order behavior by varying an external ~ a r a r n e t e r . The ~ ~ ~point ~~?~~ at which the qualitative change occurs is analogous to a critical point of a second-order phase transition4 (entry 1 2 The Journal of Physical Chemistry, Vol. 8 1, No. 25, 1977
2386
Jack S.Turner
of Table I). Preliminary studies by the author suggest that arise from either super- or subcritical bifurcations (entries 6, 7, and 14 of Table I). such a point may be accessible experimentally, and would provide exceptional possibilities for the study of fluctuIn Table I several simple chemical models are classified ations. A simple one-variable model which exhibits such according to the nature of their primary transition(s) from a critical point is indicated in entry 11 of Table I. Critical the thermodynamic branch (hss). On the schematic points of a more familiar type occur in the other multiple transition diagrams at right the solid lines denote stable steady state models (entries 8, 9, and 10 of Table I). branches of solutions and the dashed lines unstable ones. d. Higher Bifurcations. For a region of parameter space In each case the thermodynamic branch begins at the left and a bifurcation parameter increases to the right. For in which the Hopf bifurcation occurs in the FN models, first-order transitions two types of diagrams are given: the a second transition appears as the bifurcation parameter first (a) delineates the region of first-order behavior, and P moves further from the primary (Hopf) bifurcation p ~ i n t . At ~ ~this , ~ point, ~ which may correspond mathethe second (b) characterizes the associated transition. Relevant model details and references to original literature matically to a secondary bifurcation, the initial limit cycle becomes unstable and the system undergoes an abrupt are found in this paper at the locations given in brackets following the model name in the table. transition to a new oscillatory regime. The new limit cycle is a highly nonlinear oscillation, similar to those usually D. Stochastic Formulation of Homogeneous Chemical observed in the BZ reaction, having period several times Kinetics. The decisive role of fluctuations near nonelonger and amplitude an order of magnitude greater than quilibrium instabilities has been emphasized repeatedly the previous cycle. The secondary transition, a qualitain the preceding sections. In recognition of this dominant tively new feature of the oscillating reaction, has been feature of nonequilibrium transitions, Prigogine and codiscussed and interpreted in the context of the Fieldworkers have introduced the concept of "order through Noyes (FN) models by the a ~ t h o r ,and ~ ~has , ~now ~ been fluctuations" to distinguish transitions to dissipative documented experimentally as well.42 structures from those which obey the ordering principle of thermodynamic equilibrium.1° Specifically this new 3. Inhomogeneous Steady States. The possibility of ordering principle for nonequilibrium systems refers to the a nonoscillatory instability in the BZ reaction, leading amplification of fluctuations beyond instability and to their ultimately to a spatially inhomogeneous time-independent state, has been a subject of controversy in the l i t e r a t ~ r e . ~ ~ ultimate stabilization by continuous exchange of matter and/or energy with the environment. In order to study Recently, however, the existence of such a spatial dissithis basic mechanism underlying all such macroscopic pative structure has been reported by Marek and coinstability phenomena, it is obvious that one must go worker~.~'The compatibility of such an instability with beyond the deterministic description and formulate the the Oregonator modellld has been verified by the author problem of chemical kinetics at the molecular level. (see Appendix D). The simplest transition to a spatial dissipative structure is analogous to a soft-mode instability In a system characterized by a large number ( N loz3) leading to a second-order (e.g., structural) phase transition of degrees of freedom, fluctuations must be taken into (entries 2 and 3 of Table I). Just as in the equilibrium case, account as soon as one adopts a macroscopic description. near such instabilities there is the possibility of long-range In such a description the state of a system is specified by order in the near subcritical region, and this expectation assigning values to a limited number (n
