J . Phys. Chem. 1992, 96, 8702-871 1
8702
In the presence of MV2+,spectra as in Figure 2 were obtained, Le., spectra which did not possess either of the characteristic absorption bands at 395 and 600 nm of MV+. This means that the stationary MV+ concentration was always very low, from which one concludes that the electron transfer from MV+ to the silver particles was much faster than from the lead particles to MV2+. This latter electron transfer is rate determining at the pH of 4.5 at which the experiments were carried out. Thus, one has to conclude that the Fermi level in the lead particles must lie close to -0.445 V, Le., the potential of the redox system MV2+/MV+. Pulse radiolysis experiments were also carried out to test this assumption. When MV+ was generated in solutions of colloidal Ag (in the absence of Pb), the half-life for the reaction with the
silver particles was found to be 0.2 s. Since the equilibration reaction in the mixed leadsilver solutions containing MV2+takes place on a much longer time scale, the oxidation of Pb by MV2+ has again to be considered as the rate-determining step. The fact that the overall rate of the equilibration reaction is proportional to [MV2+]4,3possibly reflects an adsorption phenomenon, the limiting concentration of MV2+adsorbed on the Pb particles being reached at rather small MV2+concentrations.
References and Notes (1) Henglein, A.; Mulvaney, P.;Holzwarth, A.; Sosebee, T. E.; Fojtik, A. Ber. Bunsen-Ges. Phys. Chem. 1992, 96, 154. (2) Hunter, R. J. Foundations of Colloid Science; Oxford University Press:
London, 1986; Vol. 1, p 444.
FEATURE ARTICLE Simple and Complex Propagating Reaction-Diff usion Fronts Stephen K. Scott* School of Chemistry, University of Leeds, Leeds LS2 9JT. U.K.
and Kenneth Showalter* Department of Chemistry, West Virginia University, Morgantown, West Virginia 26506 (Received: June 2, 1992; In Final Form: August 13, 1992)
Chemical waves are common features of many reactions in which there is autocatalysis. These waves typically travel at a steady velocity and with a particular, constant waveform. The speed and waveform reflect the coupling of the autocatalytic reaction to diffusion. For some autocatalytic systems, waves can be initiated with ease, no matter how small the initiation stimulus. In other cases, however, depending on the exact form of the kinetics, the shape of the reaction zone, or the Occurrence of competing reactions that remove the autocatalybt, waves may develop only for certain ranges of the rate constants and only when a sufficiently large initial stimulus is applied. The analogous phenomenon for flame propagation is known as quenching and gives rise to the existence of flammability limits and minimum ignition energies. Depending on relative diffusivities, chemical waves propagating in two or three dimensions may spontaneously develop patterned fronts analogous to the appearance of “cellular” flames.
1. Introduction
Imagine a chemical reaction starting at a point and spreading through a reactant medium like a burning f i r e b u t with no heat. Just such “isothermal flames” arise in chemical systems when autocatalytic reaction couples with diffusion. A front or moving zone of reaction propagates into fresh reactants and leaves behind final products, typically with constant velocity and a constant waveform. The front itself is the spatial realization of the reaction event. Reaction-diffusion fronts may also exhibit complex behavior arising from inherent instabilities, behavior reminiscent of destabilized flames. In this article we consider fundamental types of propagating fronts, how they are related, and how they may exhibit complex behavior. Propagating reactiondiffusion fronts were first studied around the turn of the century as models for wave behavior in biological systems.’-* However, as recently as 10 years ago they fell into the category of “exotic phenomena”, as only a few experimental examples were known and their mechanisms were poorly understood. Today, many autocatalytic reactions are known to support propagating fronts? and more complex wave behavior is of widespread interest for modeling excitable media in biological system^.^ The theoretical treatment of reactiondiffusion fronts,
while first addressed over a half-century ago,$’ has also advanced in recent years. Many features are now well understood, and in addition, new theoretical challenges are apparent. All chemical waves, simple or complex, reflect the basic features of their fronts: propagation velocity and waveform are directly linked to the chemical kinetics of autocatalysis and its coupling with diffusion. We consider the two simplest forms of autocatalysis, quadratic and cubic, which we contrast and then join together to illustrate the general features of propagating reaction-diffusion fronts. These forms have long been thought tg be fundamentally different in character; we describe here their similarities and how one transforms into the other as the form of autocatalysis varies. Propagating fronts with quadratic kinetics have a long history, first formally studied in 1937 by Fisher’ and, independently, by Kolmogorov et a1.6 While this type of front has been considered the generic prototype for chemical waves (and propagating waves in biological systems), its mathematical description has remained an enigma-a simple one-variable reactiondiffusion equation without a general analytical solution. Propagating fronts with cubic kinetics were first considered by Semenov et al.’ in 1939; however, these fronts have received less attention over the years.
0022-3654/92/2096-8702%03.00/0 0 1992 American Chemical Society
The Journal of Physical Chemistry, Vol. 96, No. 22, 1992 8703
Feature Article
-da/dt = db/dt = k,ab
+ k,ab2
(1)
If the initial concentrations are a = a. and b = bo = 0, the rate is zero. Reaction must be initiated by seeding the system with some nonzero amount of the autocatalyst B in this simplest representation of autocatalysis. Once the reaction has started, the rate falls to zero again as it approaches complete consumption of the reactant A, a 0, b ao. In order to make mathematical progress with the problem of interest here, it is convenient to introduce dimensionless forms of the rate equations. For this we define dimensionless concentrations a and @ that are simply the actual concentrations divided by the initial concentration of A a = ala,, 0 = b/%. Thus both vary between 0 and 1, with a = 1 and B = 0 at the beginning of the reaction and a = 0, = 1 for the final state. We also need to introduce a time scale. For this we choose the chemical time tch = (k,ao + k,ao2)-'. Both terms in this sum are pseudo-first-order rate constants with units of s-I: their inverse has units of time. We choose this particular form so that tch remains finite at the extremes of both pure quadratic (k, 0) and pure cubic (k, 0) autocatalysis. If now we use T = t / t c h , eq 1 becomes
- -
a
Figure 1. Variation of reaction rate with concentration for autocatalytic systems: (a) pure quadratic autocatalysis with rate maximum at 50% consumption;(b) mixed autocatalysis; (c) pure cubic autocatalysis with maximum at a = I / , .
They are readily described analytically and, hence, have been thought of as distinctly different from quadratic "FisherKolmogorov" fronts. We link the quadratic and cubic reaction-diffusion equations by considering them together in a mixed-order description.* A smooth transition from the cubic to the quadratic form is found, with both contributing when cubic autocatalysis is dominant and only the quadratic form contributing when it is more heavily weighted. Each extreme exhibits a "special" minimum-velocity solution with a characteristic waveform and speed, and each exhibits an infinite number of higher velocity solutions depending on initial conditions. We describe the selection of the minimum-velocity solutions in terms of the concentration-gradient phase plane, with particular trajectories between the singularities characterizing each case as well as providing insights into the "other" solutions. The effects of autocatalyst decay and the diffusivities of autocatalyst and reactant are also considered. The simple quadratic and cubic forms allow a detailed examination of these features as well as the transformation to the final asymptotic state from the initial local seeding of autocatalytic reaction. We also describe how simple laminar flames are closely related to reaction-diffusion fronts, especially those with cubic autocatalysis. We pursue this analogy to discuss instabilities in isothermal fronts, which are similar to those seen in propagating flames. Much like instabilities leading to Turing patterns,"' instabilities in propagating fronts depnd on relative diffusivities and give rise to patterns that depend on a particular wavelength selection.I2
2. Simple Fronts a. Autocrtdysis in Well-StirredSystems. Reactions for which the rate increases with increasing extent of reaction are frequently observed empirically, and are usually termed autocatalytic. Two representative forms for the empirical rate law are those of quadratic autocatalysis and of cubic autocatalysis. As illustrated in Figure la, quadratic autocatalysis A + B 2B rate = k,ab = k,a(ao - a )
-
is the simplest approximation for an empirical rate law that is basically parabolic-with approximate symmetry about a maximum rate occurring at ca. 50% conversion. The cubic representation, A
+ 2B
-
3B
rate = k,ab2 = k,a(ao -
is useful for reactions in which the induction regime lasts to higher extents of reaction, with a maximum rate at ca. 2 / 3 conversion, Figure IC. Various, real situations can be adequately represented by mixed-order autocatalysis for which the rate law is simply a linear combination of these two extremes. Figure 1b shows the case with equal contributions from both quadratic and cubic. A well-characterized system in which mixed autocatalysis has been employed successfully is the iodate-arsenite reaction.13J4 In a well-stirred system, the rate equations for a general mixed-order reaction may be written in the form
-
-
-da/dT = d@/dT = pa@
+ (1 - p)aB2
(2)
where p = k,/(k, + k,ao). This latter quantity effectively measures the fractional contribution made by the quadratic channel to the total reaction process: p varies between 0 and 1, with p = 0 c o r r q m d i n g to the pure cubic and p = 1 to the pure quadratic autocatalysis. b. Reacth-Difhmh Equetiolra We now consider an unstirred system where diffusion plays an important role. The reactiondiffusion equations for mixed-order autocatalysis have the form
In general, the Laplacian operator V z involves all three spatial coordinates, but we will be concerned mainly with systems where the symmetry allows this diffusion term to be written as a function of a single spatial coordinate. We begin with a 1-D (slab) system; in section 4 we will consider 2-D (circle) and 3-D (sphere) problems where the radial distance is the only important coordinate. In section 6 we examine systems for which more than one coordinate is also involved. Equations 3a,b can be written in dimensionless form as
aa/aT = 6 aza/ap - (1 - p ) C Y ~ 2- pas aa/aT = a2p/ap + (1 - p)ap+
(44 (4b)
where a,8, T , and p have the same meaning as before, 6 = DA/DB is the ratio of diffusion coefficients, and the dimensionlesslength [ = X / ( D ~ ~ , , , )Many ' / ~ . ions in dilute aqueous solutions have similar diffusivities, so 6 1 in many cases. In biological systems, however, or if the reactant or autocatalyst is immobilized or can complex with an immobilized species, values of 6 significantly different from unity can also arise. The boundary conditions appropriate to this representation of chemical waves are, in the same dimensionless terms,
-
aff/a[ = ag/ag = o
-
for 5
with a = 1,B = 0 for E - + = a = as,@ = j3, for E
*m
(54
and
(5b) At this stage we leave the exact values of CY, and &, the concentrations of A and B behind the wave, unspecified (these values emerge from the solution of the problem and only in the simplest cases can they be specified in advance). To complete the formulation of our problem, we need to describe the initial conditions. We imagine that a = 1 everywhere at T = 0 and that reaction is initiated by the addition of a finite (and generally small) amount Bo of the autocatalyst over some restricted -a
Scott and Showalter
8704 The Journal of Physical Chemistry, Vol. 96, No. 22, 1992 1.o,
U
B o
0 and A(Olo)< 0 are given by
--OD
da/dz = u du/dz = -CU + a(1 - a)
+ c2 exp(X+z); + c4 exp(X+z)
=
t/2(-~
+ 4]1/2)
f [c2
(13)
One is negative, indicating exponential decay toward the origin; the other is positive, indicating exponential divergence of the second term. Ultimately, divergence wins: the origin is a saddle point. Also of importance are the eigenvectors e+(O,O)and e_(O.O)associated with these eigenvalues. Small perturbations approach the origin parallel to the inset eJ0l0)and depart parallel to the outset e+('.').
The Journal of Physical Chemistry, Vol. 96, No. 22, 1992 8705
Feature Article
line locally to (1,O). This path can be integrated backward across the phase plane, and it cuts the u-axis for some positive u. If we now return to the origin for this system, the trajectory departs along, or parallel to, the outset e+(OT0), which differs only marginally in slope from before. As this trajectory evolves, it will approach the inset e,(l*o)from below and along which it then approaches the state (1,O) as z +-. This thus forms an “allowable” connection between the two singularities. The behavior close to (1,O), which corresponds to the very front of the wave where little reaction has occurred, shows exponential decay of perturbations governed by the degenerate eigenvalue Ai(I,O). For any higher speed, c > c*, an allowable connection also exists. For these cases, the singularity at (1,O) becomes a regular node with distinct eigenvalues A(lvo)< A+(l,O) < 0. The final approach of a given solution trajectory to this singularity will be along the eigenvector e+(lSo)associated with the eigenvalue of lower magnitude: again, the path emanating from (0,O)approaches this inset from below and the connection remains within the allowable region of the phase plane, Figure 4c. Small perturbations close to (1.0) again decay exponentially, governed by the eigenvalue A+(l,O). We may note that the magnitude of A+(l*O) is a decreasing function of c (for all c 2 c*), so faster waves have slower decay of transients. Thus for pure quadratic autocatalysis, constant-form, constant-velocity wavefronts with any speed c greater than some minimum c,, are allowed. The minimum velocity here corresponds to that velocity a t which the singularity at (1,O) ceases to be a focus and becomes a node, c*. Hence, quadratic autocatalysis supports waves of velocity c 2 C,i” = c* = 2 (15)
-
U
0
1
o(
0
-
1
o(
Figure 4. Evolution of the eigenvectors of the singularity at (1,O) with wave sped for quadratic autocatalysis: (a) c < c*, the singularity at (1,O) is a focus with complex eigenvalues and is approached by a spiraling trajectory that crosses out of the allowed region; (b) c = c*, (1,O) is a degenerate node and can be approached by an allowable connection from (0,O);(c) c > c * , (1,O) is a regular node with distinct eigenvectors, an allowed trajectory approaches along e+(’.0).
We will be concerned with perturbations that leave (0,O)along e+(OV0).The value of the positive eigenvalue and the direction of the outset vary with the wave speed c from the above equation. The singularity at (1,O) is stable to small perturbations, with the approach governed by the eigenvalues given by
These waves differ in their form as z +-OD, Le., the beginning of the wavefront close to a = 1 (the initial state temporally). The minimum velocity solution c = cminhas the fastest decay here, as the eigenvalue A+(’*0) has its greatest magnitude for this speed: this gives a key to the special nature of this particular solution. Having determined that a semi-infinite number of solutions to our boundary-value problem exist, we now seek to find which particular solution will be selected from this array for a given set of initial conditions. The problem of selection has been discussed by various auFronts initiated by typical “localized” inputs of B (so-called “compact support”) will evolve to the minimum velocity solution with c = cminabove. Only if Bo(z) is a more slowly-decaying function spatially, of the form eTz for large z, with 0 < y < 1, does this initial distribution have sufficient effect far ahead of the front at large z (and so close to (1,O) in the phase plane). This initial seeding gives rise to a phase wave, which is effectively just a series of local clock reactions that appear as a propagating wave with velocity c = (1 + y2)/y. d. Pure Cubic Fronts ( p = 0). An analysis similar to that given above can be followed for the cubic autocatalytic system d2a/dz2
-
+ c da/dz - a(1 - a)’ = 0
- -
- -
(16)
(14)
with the boundary conditions da/dz = 0 as z f-,a 1 as z +-, and a 0 as z --OD. In this case, the singularity at (1,O) is a saddle node for all positive wave speeds, i.e., for c 2 c* = 0. There are, however, some additional complexities here that cause the situation to be different. The Jacobian matrix now has the form
For wave speeds c < c* = 2, the discriminant of this equation is negative, so the eigenvalues are complex conjugates. This corresponds to the singularity being a stable focus. Trajectories would approach (1,O) by spiraling into it. This requires the trajectory to repeatedly traverse regions of the phase plane for which a > 1 or u < 1, Figure 4a. Neither of these regions is accessible, however, because of the initial conditions chosen; therefore, such solutions and wave speeds are not allowable for this problem. With c = c*, the singularity at (1,O) becomes a degenerate node, Figure 4b, with equal exponents = = -1. The eigenvector e,(’*0)associated with this forms a path that is a straight
so for a = 1, the eigenvalues are A+(I*O) = 0, A_(’*0) = -c. For c = 0, both eigenvalues are zero and even for c > 0 there is one zero eigenvalue A+(I,O). This means that the system does not always follow the simple exponential growth or decay embodied in eq 12a. We return to this point below. A different approach is fruitful for this particular rate la^.^^^^,'^^^^ Again, we work with the a-u phase plane. The
Ai(l9O)
Ai(l.o)
=
l/z{-c
f
[$
- 4]1/2)
Scott and Showalter
8706 The Journal of Physical Chemistry, Vol. 96, No. 22, 1992
singularity at the origin is a saddle point as before, with an outset e+(OYo)associated with the positive eigenvalue A+(O*O). This determines the departure from (0,O). For nonzero wave speeds, the phase plane in the vicinity of (1,O) consists of an inset e_(Ivo) associated with the nonzero eigenvalue. Associated with the zero eigenvalue is a center manifold CM(Iqo)which approaches (1,O) tangent to the a-axis. The approach along e_(l.O)is exponential, while the approach along the manifold is nonexponential. If the inset is integrated backward in z, it describes a path a c r w the phase plane. This is the only trajectory that approaches ( 1,O) along this direction. All other trajectories approaching (1,O) do so tangent to the center manifold. For sufficiently low speeds, e_(lv0)cuts the a-axisbetween 0 and 1 and would then cut the u-axis for u < 0. As c is increased, this intersection moves toward the origin, and for sufficiently large c, the u-axis is cut for u > 0. We will see that this change is important below. In order to find an allowable solution to the boundary value problem in terms of a suitable connection in the acu phase plane, we might guess at a form for the dependence of u on a (Le., of spatial gradient on concentration). For this, the parabolic form u = ka(1 - a) is appropriate, where k is a constant yet to be determined. Substituting this into eq 16, we obtain
(a)
I
U
e_(ITo)
- 2a) + cka(1 - a) - a(l - a)’ 0 (18) The common factor a( 1 - a) cancels, leaving the condition k2a(l - a)(l
k(k
+ c ) - 1 = (2k2 - 1)a
(19)
For this to hold for all a,then both sides must be identically zero, requiring 2k2 = 1 and then yielding C t = 1/21/2
CM
(20)
for the unknown wave speed. The assumed form da/dz = ka( 1
- a)can be integrated with k = 1 /21/2to give the equation for the profile 1
-
U
= (1 + (f/2”2)-1
(21)
The profile is symmetric about z = 0, with maximum gradient at z = 0. umax= 1/(32)1/2for a = This analysis appears to be significantly different from that for quadratic autocatalysis in that a single solution has emerged with a single wave speed. There are, however, additional solutions which do not have the special parabolic path in the phase plane. These all have wave speeds greater than 1/2II2, so like the quadratic rate law this cubic system also has a semi-infinite spectrum of solutions with c 1 c,,
= C+ = 1/2’/2
(22)
The relationship of these additional solutions to the one determined above can be seen from the phase plane. With 0 < c < c,,, the inset cuts the a-axis between the two singularities, Figure 5a. Any solution path leaving (0,O)along or parallel to the outset e+(OV0) lies above this *separatrix”. Such trajectories cannot approach (1,O) along this inset, nor can they cross it to approach along the center manifold. For these low speeds, there is no allowable connection in the phase plane. For the special case, c = c,,,,,, the inset e_(I*O) to (1,O) evolves back across the phase plane and approaches the saddle a t the origin as z -m, Figure 5b. The approach to the saddle is along the outset of that singularity e+(OVo).There is thus now a connection from (0,O)to (1,O). This approaches (1,O) as z +a along the inset, not the center manifold, and is exactly the parabolic connection “guessed” above. For higher wave speeds, the inset e_(l*O) is integrated backward to cut the u-axis first. The outset from the saddle e+(OVo)lies below this path. The solution trajectory now leaves the saddle along its outset and approaches (1,O) along the center manifold, Figure 5c. Thus connections also exist for all c > cmin,but these do not have the exponential characteristics of the minimum velocity solution or of the solutions for quadratic autocatalysis. Again there is the problem of wave selection, and again the system seems always to evolve to the minimum velocity solution
-
-
CM Figure 5. Evolution of paths in the phase plane with wave speed for cubic autocatalysis: (a) c < c+, a trajectory starting from (O,O) lies above e>’*o) and does not approach (1,O); (b) c = ct, the eigenvector eL1*O) connects with the outset e+(Oso) leaving the origin to form a parabolic trajectory connecting the singularities; (c) c > c’, the outset from (0,O)lies below and the trajectory approaches (1,O) along CM. if the initial input of B has compact support. Billingham and NeedhamIg have conjectured the following: if B((,O) 5 1/(2’/2[) if B(4,O) if B([,O)
-
a/€
> u/(
C
c
for all 5
u
Cmin
( u > 1/2’/*)
(23)
no traveling wave develops
In other words, a moderate forcing of the velocity can be achieved by seeding over a sufficiently wide distance scale, but if the input of B becomes too large, no steady wave solution can emerge from the local kinetics. e. Mixed-order Fronts (0 5 c 5 1). For the general case with arbitrary p in the range 0 Ip I1, we can apply either of the methods discussed above: the quadratic method of looking for a change from focus to node at c = c* or the cubic method which assumes a parabolic connection and determines the speed ct from this. These yield c* = 2(p)’/2
and
(24a)
The Journal of Physical Chemistry, Vol. 96, No. 22, 1992 8707
Feature Article
31
/
/
/
2
have exponential decay at their leading edge, the exponent characterizing this decay changes discontinuously from c = cmin = ct to c > cmin.The evolution of the phase plane is similar to that presented earlier, Figure Sa-c, except that the center manifold is now replaced by a regular inset associated with the nonzero eigenvalue A+(I?O) for p > 0. The special case p = corresponds to the degenerate inset e,(IIo) for c = c* being such that it enters the origin as it is integrated backward in z, so c* = ct = (4/3)1/2. f. Observed Wave Velocities. We may summarize the results above by stating that for localized inputs of the autocatalyst the observed wave velocity will be cobs= (1 + p)/[2(1 - p ) ] I l 2 for 0 Ip I!I3 (25a)
P Figure 6. Wave speeds c* and c' as functions of I./ for mixed autothe system chooses c'; for < p < 1 the catalysis. For 0 < p < system chooses c'.
For c = ct, the profile is symmetric about z = 0 as for pure cubic autocatalysis, while for c = c* the profile is, in general, asymmetric. The variation of these two speeds with p is shown in Figure 6. The cubic curve ct(p) always lies above that for the quadratic solution c * ( p ) for all p but they touch tangentially at p = The two sections, ct(p) for 0 Ip I with ct 1/2Il2 as p 0, and c*(p) for Ip I1, with c* 2 as p 1, are shown as full lines for reasons we now discuss. With p = 0 (pure cubic autocatalysis), we have possible solutions for all wave speeds c L ct and ct > c* = 0. With p = 1 (pure quadratic autocatalysis), the range of allowable velocities is c 1 c*, and so waves with speeds less than required for a symmetric profile (c = ct +-) can be observed. For both cases, the system evolves from compact initial input of B with the minimum speed available (c' for p = 0 and c* for 1.1 = 1). For the mixed case, the behavior switches between these two forms, and this switch occurs as p = With p > the system is essentially quadratic in character. As c is increased to c*, the focus at (1,O) changes to a degenerate node (two equal, real eigenvalues). The inset e,(lP) traverses back across the phase plane and cuts the u-axis for some positive u. The outset from the origin starts underneath this inset, and approaches (1,O) along it, remaining in 0 Ia I1 throughout. A similar situation exists for c > c*, for which (1,O) is a regular node with distinct eigenvalues and eigenvectors. The solution path leaves (0,O)along the saddle outset and approaches (1,O) along the inset This associated with the less negative eigenvalue inset is approached from below for all c. This evolution of the phase plane with c is exactly analogous to that shown in Figure 4a-c. There is thus an allowable connection for all c 1 c,,, = c*, and of these, the minimum velocity solution is that governed by fastest exponential decay. The cubic form solution exists, with finite velocity ct > c* but is not selected unless special initial conditions are imposed. the focus at (1,O) becomes a degenerate node for For p < c = c*. Now the inset ei(l>O)integrates backward across the plane and cuts the a-axis between 0 and 1. Any trajectory originating from (0,O)lies always above this inset and so approaches (1,O) along the other segment of this inset-that lying in a > 1. Thus, although such a connection exists, it is not acceSSible with the initial conditions imposed (0 Ia I1 for all I ) . A similar situation exists for slightly larger speeds. Now the node at (1,O) is regular, but both insets and cut the a-axis and so any trajectory starting from (0,O)above these must enter (1,O) along the inset e+(l?O)from a > 1. Only when c has increased to ct does an allowable connection occur. For c = ct the cubic-type solution, the inset associated with the more negative eigenvalue, integrates backward so as to exactly arrive on the outset of (0,O). This gives the first allowable connection, as a I 1 for all z, and is governed by A_(l,o) in the vicinity of (1,O). This profile will be symmetric about z = 0. For all larger c, there remain allowable connections,which approach (1,O) now along the "slower manifold" associated with the smaller magnitude eigenvalue Thus although there is a continuous spectrum of speeds and all these
- --
-
-
cobs= 2 ( ~ ( ) l / for ~
y3 5 p I
1
(25b)
Thus, cobsvaries smoothly with p for all 0 Ip I1, but there is a discontinuity in the curvature of this graph at p = when both the cubic and quadratic channels contribute equally. These forms may be converted back to the actual velocities in laboratory coordinates and in dimensional terms to give dx/dt = (f/2Dk,ao2)l/2(1+ 2k,/k,ao) dx/dt = 2(Dkgao)1/2 for
for 0 I p 5 !I3 (26a) 5 p I1
(26b)
Thus if the cubic channel dominates (p < the wave speed depends on both the rate constants k, and k,. If we start with pure cubic autocatalysis (k 0) and then increase the contribution from the quadratic &innel, the speed increases linearly with k, until p = For systems with p > I/,, the wave speed would then increase as the square root of k, and becomes completely independent of the cubic channel k,. This latter point means that if we start with pure quadratic autocatalysis and subsequently increase k,, the wave speed will not vary until p
I / 3 , showing quadratic type behavior, only the diffusion of the autocatalyst is important,20a2'so dx/dt = 2(DBkqao)'/Z for
Ip I 1
(27)
For the cubic-type solutions, the situation is less simple. However, the observed (minimum) velocity tends to zero as D s / D A tends to zero. For Dg/DA >> 1, the dimensionless velocity grows as (DB/DA)'l2and then the actual velocity dx/dt will again depend only on Dg112and be independent of DA.20y21We return to the effects of unequal diffusion coefficients in section 6 to discuss how fronts may become unstable to perturbations.
3. Autocatalysis with Decay In this section we consider the passage of a wave through a chemical system with either quadratic or cubic autocatalysis, but in which these processes are in competition with a "decay" step where the autocatalyst B becomes inactivated by reaction to a stable product species C. The general chemical scheme can be written A
+ mB
-
nB
(m + l ) B nC
rate = k,abm
rate = k2bn
where we allow nth-order kinetics for the removal step. We will concentrate mainly on the cases m and n = 1, 2; Needham and Merkin22have discussed general values of m and n. The system is now described by two coupled reaction-diffusion equations aa/at = D a2a/ax2 - klab"
ab/at = D a2b/aX2
+ klabm- k2bn
(284 (28b)
8708 The Journal of Physical Chemistry, Vol. 96, No. 22, 1992
Scott and Showalter From the competition governed by the term j3 - y, clearly we require an initial input such that j3 > y. This is an important result. Now, the system ahead of the wave has gained a limited, local stability to small perturbations. With the decay process, small fluctuations or random seedings ahead of the wave will not initiate a new front-in contrast to the system with no decay. d. Cubic Autocatalysis with Quadratic Decay.23 With m = n = 2, the autocatalysis and decay processes compete effectively on equal terms (both depending on b2at the leading edge where a = 1). The condition for a constant-velocity traveling wave will again be y < 1. There is no critical initiation requirement: a wave develops for any nonzero Bo if y < 1. The wave velocity again tends to zero as y 1, although now as (1 - y)3/2rather than the square root form seen in (a) above. There is incomplete consumption of the reactant through the wave, with a = a,> 0 behind the wavefront.
-
0
21
0.0465
Figure 7. Variation of wave speed with dimensionless decay rate constant y for cubic autocatalysiswith linear decay: for y < 0.0465 there are two wave speeds, the higher correspondingto the stable wave. No reactiondiffusion waves exist for y > 0.0465.
where we assume equal diffusion coefficients for simplicity. In dimensionless form these become
aa/ar = a2a/ap aj3/ar = a2j3/ap + aj3m - yp
(29a) (29b)
The parameter y = nk2/k,c;gmC1-n reflects the strength of the decay relative to the autocatalysis. In the analysis of these equations, the most significant feature is the behavior at the very leading edge of the wavefront. There a = 1 and so the kinetic term in the equation for j3 involves the difference bmn- y. The sign of this term will be crucial to whether a constant form wave can be supported. a. Quadratic Autocatalysis with Linear J l e ~ a y . 2 ~Here 9~~ m = n = 1. The difference j3" - y becomes simply 1 - y, which is positive for y < 1 and negative for y > 1. In this case a wavefront propagates provided y < 1. Then the production of B, which is given by aj3 j3 at the leading edge, can outweigh the decay there which has rate yj3. The wave will tend to a constant velocity c = 2(1 - Y ) ' / ~for any y < 1. This clearly tends to the results derived earlier for y = 0, and the wave speed tends to zero as y ---* 1. Within the range 0 5 y < 1, such a wave develops for any nonzero input of B, no matter how small. Some other features are important. The dimensionless concentration of B does not remain high after the passage of the wavefront. Instead the decay reaction converts B to C, so j3 falls again. The @([)or B(z) profiles now have a pulselike form. Also, the concentration of the reactant A does not fall completely to zero. There is a nonzero residual concentration a, > 0 behind the wavefront. b. Quadratic Autocatalysis with Quadratic Decay.Z3 If the decay process is second-order in B, so m = 1 and n = 2, the competition at the leading edge where j3 is small and a is of order unity involves the difference 1 - yo. For sufficiently small j3, the tendency for autocatalytic growth thus always exceeds that for the decay process, no matter how large y becomes, Le., yj3 < 1 for any y as @ 0. Thus,a wave propagates in this case for all 7 and there is no requirement on the input of B. Also, there is complete conversion of the reactant through the front (Le., a, = 0). c. Cubic Autocatalysis with Linear Decay.23926!27 In this case (m = 2, n = l), the competition at the leading edge of the front has the form - yj3 with a 1, so the sign of this term is given by j3 - y. Reaction waves only develop successfully provided y < yn = 0.046, where ycris determined numerically. As the decay rate constant approaches this critical value, the wave velocity remains finite rather than tending to zero. In fact, there are two wave solutions for any y < ycr,one stable and one (of lower velocity) unstable. These two solutions merge at y = yn, as shown in Figure 7. Even if the condition on y is satisfied, a permanent wave solution will develop only if sufficient B is added to initiate the reaction.
-
-
J
4. Nonplanar Wave Propagation Next, we consider the wavefronts that may develop for two other geometries: a circle and a sphere. These represent the simplest 2-D and 3-D geometries. The natural response to circular or spherical initiations at the center of such reaction zones will be the development of fronts with the same underlying shape. We concentrate here on the quadratic and cubic autocatalysis,without decay. For a circle, the Laplacian term V 2 ahas the form d2a/dz2 (1/z)(da/dz) and this provides an extra term (1 /z) da/dz in eq 9. For the case p = 1, pure quadratic autocatalysis, the wave development is very similar to that for the 1-D geometry. A wave is initiated for any nonzero input of autocatalyst boat the center. The velocity of the wave tends to a constant value equal to that found previously, c 2. Once the radius of the wave has grown sufficiently,each local segment has the appearance of a 1-D plane wave and curvature effects become unimportant. This situation also holds for quadratic autocatalysis in spherical reaction zones. The cubic autocatalysis proves to be more interesting. In a circle, a circular wavefront is initiated by any nonzero input, no matter how small. However, there is now a long transient development stage and for small inputs (Bo 1 in eqs 4a,b. Thus, the front is destabilized when the diffusivity of the reactant A becomes sufficiently larger than that of the autocatalyst B. We consider systems with pure cubic autocatalysis here; in a detailed study, pure quadratic and mixed-order fronts are also considered and found to have different sensitivities to the destabilizing effects of diffusion.12 Consider a two-dimensional reaction domain, of infinite extent in the x-direction but of finite width in y. A front initiated by supplying a small amount of B is allowed to develop until a constant velocity and waveform are observed. Shown in Figure 9a is the evolution of a front from an initial stepfunction seeding to a final asymptotic form. As shown by the contours, this front, calculated with 6 = 1, is stable to perturbations in the sense that the spatial nonuniformity imposed in the initiation decays over time, with a planar waveform developing as the asymptotic state. This behavior, which is typical in chemical systems, is contrasted with that shown in Figure 9b, calculated with 6 = 5. Now, rather than the perturbation decaying, it gives rise to a waveform reminiscent of those seen in destabilized flames. The final state is a front with a spatially-periodic waveform in the y-direction. In some cases, a constant-velocity wave emerges, for which the shape of the front at any particular value of y does not change in time. Nonsteady fronts for which the pattern continuously varies in time are also observed. The wave in Figure 9b has this spatiotemporal instability. The same waveform appears with a variety of initial conditions for this particular value of 6 and domain width. The concentration profiles of a and j3 in the waveform are shown in Figure 10. Although both profiles reflect the same basic structure, the j3 profile shows greater deviations from the uniformity seen when 6 = 1. We may explain the instability of the planar front and subsequent evolution to the patterned state in much the same way as flame destabilization: instability occurs when the destabilizing effects of reactant diffusion outweigh the stabilizing effects of autocatalyst diffusion. The constant-velocity, patterned waveform develops when the growth of the perturbation just offsets the stabilizing effects of autocatalyst diffusion. The instabilities leading to patterned reaction-diffusion fronts have features in common with diffusion-induced instabilities leading to Turing patterns. The greatest similarity is the requirement of differing diffusivities of the reactant and autocatalyst. The analogy is especially striking with Turing instabilities of cubic autocatalysis schemes.lOJl Other common features include the onset of instability a t a critical ratio of the diffusivities, ,a = 2 in the case of pure cubic fronts, and patterns with intrinsic wavelengths. The wavelength of the cubic front shown in Figure 9 increases with an increase in grid width until there is sufficient
Feature Article 1.o
o( 0.0
1.o
r4 0.0 Figure 10. Two-dimensional waveform corresponding to Figure 9b.
room for another pair of oscillations, which grow in to restore the original wavelength. We anticipate an intrinsic wavelength that is insensitive to size in a sufficiently large system; however, spatiotemporal instabilities make the analysis of asymptotic behavior difficult. Further analysisI2suggests that the critical ratio of the diffusion coefficients increases for the mixed-order rate law, with a, (2 - p)/(1 - p ) for p < but that such instabilities are not supported by the quadratic-type waves with Ip I1. The experimental realization of diffusion-induced front instabilities may be possible by following methods utilized in the recent experimental demonstrations of Turing pattern^.^^*^*^^ The key lies in designing experimental systems in which the diffusivity of the autocatalyst B is retarded relative to that of a control in the case of cubic autocatalysis the reactant A.
-
7. Conclusion
All chemical waves reflect the features of their fronts. We have described in this article how the basic fronts are related and how they may give rise to differences in behavior-differences that are important in certain instances. The mixed-order description provides insights into the quadratic form, often employed as the prototype for propagating fronts, as well as the cubic form, which has features related to reaction-conduction flames. Diffusioninduced instabilities in 2-D and 3-D fronts offer interesting experimental and theoretical challenges for studies of dynamical complexity. Acknowledgment. We thank NATO (Grant 0124/89), the National Science Foundation (Grants CHE-8920664 and INT8822786) and WV-EPSCoR for financial support; S.K.S.is grateful for a Fulbright Visiting Scholarship Award. We also thank Dezso Horvath, Valery Petrov, Drs. Vilmos Ggspiir and David J. Needham, and Professors John H. Merkin, Peter Gray, and John J. Tyson for enlightening discussions.
References and Notes (1) Luther, R. 2.Elektrochem. 1906,12 (32),596. (2) Arnold, R.; Showalter, K.; Tyson, J. J. J. Chem. Educ. 1987,64,740. Showalter, K.; Tyson, J. J. J. Chem. Educ. 1987,64,742. (3)See, for example: (a) Oscillations and Traveling Waves in Chemical Systems; Field, R. J., Burger, M., Us.; Wiley: New York, 1985. (b) Epstein, I. R. J. Phys. Chem. 1984.88,187. (c) De Kepper, P.; Boissonade, J.; Epstein,
The Journal of Physical Chemistry, Vol. 96, No. 22, 1992 8711 I. R.J . Phys. Chem. 1990,94,6525.(d) De Kepper, P.; Epstein, I. R.;Kustin, K.; Orbln, M. J. Phys. Chem. 1982,86,170. (e) Weitz, D. M.; Epstein, I. R. J. Phys. Chem. 1984,88, 5300. ( f ) Gribshaw, T. A.; Showalter, K.; Banville, D.; Epstein, I. R. J. Phys. Chem. 1981,85,2152. (8) Razsa, G.; Epstein, I. R. J. Phys. Chem. 1985,89,3050. (h) Nagypgl, I.; Bazsa, G.; Epstein, I. R.J. Am. Chem.Soc. 1986,208,3635.(i) Szirovicza, L.;NagypB1, I.; Boga, E.J. Am. Chem. SOC.1989,121, 2842. (j) Garley, M. S.;Jones, E.; Stedman, G. Philos. Trans. R. Soc. London, A 1991,331, 237. (k) Showalter, K. J. Phys. Chem. 1981,85, 440. (I) Boga, E.; Peintler, G.; Nagypil, I. J. Am. Chem. SOC.1990,112,151. (4)See, for example: (a) Winfree, A. T. The Geometry of Biological Time; Springer: Berlin, 1980. (b) Winfree, A. T. Chaos 1991,1, 303. (c) Ross, J.; Muller, S. C.; Vidal, C. Science 1988,240,460. (d) Tam, W. Y.; Horsthemke, W.; Noszticzius, Z.; Swinney, H. L. J. Chem. Phys. 1988,88, 3395. (e) Kshirsagar, G.; Noszticzius, Z.; McCormick, W. D.; Swinney, H. L. Physica D 1991,49, 5. ( 5 ) Fisher, R. A. Ann. Eugen. 1937,7, 355. (6)Kolmogorov, A.; Petrovsky, I.; Piscounoff, N. Bull. Univ.Moscow, Set. h t . , Sect. A 1937,I, 1. (7)Voronkov, V. G.; Semenov, N. N. Zh. Fiz. Khim. 1939, 13, 1695. (8) Gray, P.; Scott, S.K.; Showalter, K. Philos. Trans. R . Soc. London, A 1991,331,249. (9)Turing, A. M. Philos. Trans. R . SOC.London, B 1952,237,37. (10)Murray, J. D. Mathematical Biology;Springer-Verlag: Berlin, 1989. (1 1) Dufiet, V.; Boissonade, J. J . Chem. Phys. 1992,96,664. (12)Horvith, D.; Petrov, V.; Scott, S.K.; Showalter, K. J. Chem. Phys., submitted for publica tion. (13) Hanna, A.; Saul, A.; Showalter, K. J. Am. Chem. Soc. 1982,104, 3838. (14)Saul, A.; Showalter, K. In Oscillations and Traveling Waves in Chemical Systems; Field, R. J., Burger, M., Eds.; Wiley: New York, 1985; Chapter 11. (15) Merkin, J. H.;Needham, D. J.; Scott, S . K.Proc. R . Soc. London, A 1989,424,187. (16)Merkin, J. H.;Needham, D. J. J. Eng. Math. 1989,23, 343. (17) Gray, P.; Showalter, K.; Scott, S . K. J. Chim. Phys. 1987,84,1329. (18) Gray, P.; Merkin, J. H.; Needham, D. J.; Scott, S . K. Proc. R. Soc. London, A 1990,430,509. (19)Billingham, J.; Needham, D. J. Dyn. Stab. Syst. 1991,6,33. (20)Billingham, J.; Needham, D. J. Philos. Trans. R . SOC.London, A 1991,334, 1. (21)Billingham, J.; Needham, D. J. Philos. Trans. R . SOC.London, A 1991,336,497. (22)Needham, D. J.; Merkin, J. H. Philos. Trans. R . SOC.London, A 1991,331,261. (23) Merkin, J. H.;Needham, D. J. Proc. R . Soc. London, A 1990,430, 315. Merkin, J. H.; Needham, D. J. Proc. R. SOC.London, A 1991,434,531. Needham, D. J. J. Appl. Math. Phys. 1991,42,455. (24)Kallen, A.; Arcuri, P.; Murray, J. D. J. Theor. Biol. 1985,116,377. (25) Merkin, J. H.; Needham, D. J. Nonlinearity 1992,5, 413. (26)Novozhilov, B. V.; Posvyanskii, V. S. Fiz. Goreniya Vzryua 1973.9, 225. (27)Novozhilov, B. V.; Posvyanskii, V. S . Fiz. Goreniyu Vzryva 1974.10, 94. (28)Merzhanov, A. G.; Khaikin, B. I. Prog. Energy Combust.Sci. 1988, 14, 1. (29)Clarke, J. F. Prog. Energy Combust. Sci. 1989, 15,241. (30)Williams, F. A. Combustion Theory, 2nd ed.; Addison-Wesley: New York, 1988. (31) Sivashinsky, G. I. Combust. Sci. Technol. 1977,15, 137. (32) Sivashinsky, G. I. Annu. Rev. Fluid Mech. 1983,15, 179. (33) Rakib, Z.;Sivashinsky, G. I. Combust. Sci. Technol. 1987,54,69. (34)Kaper, H. G.; Leaf, G. K.; Margolis, S. B.; Matkowsky, B. J. Combust. Sci. Technol. 1987,53,289. (35) Pojman, J. A.; Epstein, I. R. J. Phys. Chem. 1990,94,4966. (36)Pojman, J. A.; Epstein, I. R.; McManus, T. J.; Showalter, K. J. Phys. Chem. 1991,95,1299. (37)Pojman, J. A.; Nagy, I. P.; Epstein, I. R. J. Phys. Chem. 1991,95, 1306. (38) Castets, V.; Dulos, E.; Boissonade, J.; De Kepper, P. Phys. Rev. Lett. 1990.64, 2953. (39)De Kepper, P.;Castets, V.; Dulos, E.; Boissonade, J. Physica D 1991, 49, 161. (40)Ouyang, Q.;Swinney, H. L. Nature 1991,352,610. (41)Lee, K. J.; McCormick, W. D.; Swinney, H. L.; Noszticzius, Z.J. Chem. Phys. 1992,96,4048. (42)Lengyel, I.; Epstein, I. R.Science 1991, 251,650.
