Dynamics of Forced Excitable and Oscillatory Chemical Reaction

J . Phys. Chem. 1989, 93, 2164-2714

Downloaded by NEW YORK UNIV on September 12, 2015 | http://pubs.acs.org Publication Date: April 1, 1989 | doi: 10.1021/j100344a015

2764

the autocatalytic reaction is prolonged which results in broad wave fronts with somewhat higher amplitude. At high [MA] and [BrMA] (Figure lg) the fast production of Br- stops the autocatalytic step quickly, resulting in slightly smaller but very sharp wave fronts. The organic processes are also responsible for differences in the wave profiles of series A and B in Figure 1. Series A has high initial [BrMA] (0.09 M), whereas the initial solutions of series B contain only MA ([MAIB = [MA], + [BrMA],), a part of which is brominated to BrMA during an introductory period necessary for the solution layer to become excitable. This [BrMA] is far less than [%MA], = 0.09 M. Therefore in series A, vRIO > vR9, Ce(IV) produced in the autocatalytic reaction is reduced mainly in (RIO), producing immediately enough Br- to stop the autocatalytic reaction. On the other hand, the Br- production in series B is slower. Obviously this leads to a clearly distinguishable separation of the increasing Ce(1V)-concentration segment into two parts: The autocatalytic Ce(IV) production dominates in the segment of fast Ce(IV) increase. The Ce(1V) consumption in the organic processes causing a negative feedback through Br- production becomes observable in the segment of slow Ce(IV) increase. At extremely high MA concentrations there is almost no difference between series A and B because the pool of MA is very large as compared to the [BrMA] either initially present or produced during the introductory period (Figure lg). The [Ce(IV)] of the excitable solution layer in front of the waves (eleventh column in Table 11) depends significantly on the initial concentrations of Br03-, Ce(IV), H’, and MA. However, this minimum [Ce(IV)] expressed in percent of the total catalyst concentration predominantly depends on the initial concentrations of Br03- and MA + BrMA-the main oxidizing and reducing agents of the reactant solution-and H+, which affects their redox properties. In fact, by increasing or decreasing the [BrO 0 according to eq l a defining a jump in the phase space. We suppose that the jump operator A is continuously dependent on the parameter A I0 so that A is an identity operator when A = 0. Thus A can be interpreted as an amplitude of the pulse. The time t can be scaled with a characteristic time of the autonomous system (1b) (e.g., a period of oscillations). Given a value of A and T, asymptotic solutions of eq l a and 1b represent characteristic dynamical regimes. Each regime can be characterized by the firing number v. In general v is dependent on an initial condition which must be specified to obtain a solution of eq l a and lb. Thus multiple values of v can exist for a given choice of A and T. The firing number may change as A and T are varied. Hence the parameter plane contains regions where a given firing number exists and these regions may overlap. Periodic regimes (resonances) are characterized by a rational value of v ; Le., v = p / q (p and q are integers) and the corresponding region in the parameter space is called a p/q-resonance region (or a p/q-Arnold tongue); cf., for example, Figure 7. Provided that different resonance regions do not overlap, the dynamics is periodic if v is rational or quasiperiodic if v is irrational. The dependence of v on T or A is a piecewise constant continuous function called a devil’s staircase;I2 cf., for example, Figure 6. Resonances correspond to plateaus of the devil’s staircase. The ordering of resonances in the staircase is given by Farey sequence~;~’ i.e., given two rational numbers, p / q and r / s , their Farey sum (p r ) / ( q s ) is located between them. The parameter dependence of v is a multivalued function if resonance regions overlap, and one may expect a chaotic dynamics to be present in such situations.I4 However, this cchaos may be unstable (Le., the chaotic set in the phase space is nonattracting), and a coexisting stable resonance regime will be then observed instead of the chaotic one.

+

+

3. Forcing of Oscillatory BZ Reaction in CSTR 3.1. Single Pulse Experiments and PTC Curves. The periodically oscillating BZ reaction in a CSTR exhibits positive or negative phase shifts of oscillations when a single pulse perturbation by bromide ions is used.lS When the Br- pulse is delivered soon after a sharp increase of the redox potential (the refcrence event), a positive phase shift of oscillations is usually observed (cf. point A in Figure la). Negative phase shifts occur for perturbation applied close to the end of the cycle, Le., just before the reference event takes place (cf. point B in Figure la). By plotting “old” (p) and “new” (p’) phases, we obtain the phasetransition curve (PTC). The PTC’s were reported in ref 8 and 15. For low amplitudes of perturbation PTC’s of the type “1” exist, and for higher amplitudes PTC’s of the type “0” occur.I6 (11) Schreiber, I.; Dolnik, M.; Choc, P.; Marek, M. Phys. Lett. A . 1988, 128, 66.

(12) Herman, M. R. In Geometry and Topology; Lecture Notes in Mathematics 597; Springer: New York, 1977; p 271. (13) Hardy, G. H.; Wright, E. M., eds. An Introduction to the Theory of Numbers; Claredon: Oxford, 1979. (14) MacKay, R. S.; Tresser, C. Physica 1986, 190, 206. (15) Dolnik, M.; Schreiber, I.; Marek, M. Phys. Left.A 1984, 100, 316. (16) Winfree, A. T. The Geometry of Biological Time; Springer: Berlin, 1980.

2766

Dolnik et al.

The Journal of Physical Chemistry, Vol. 93, No. 7, 1989 -3 log P -I

-5

I

2

1

T

3

Figure 2. Parameter plane log A - T experimental measurements. Full circles denote periodic regimes; regions of periodic regimes with the same firing number Y = p/q form Arnold tongues. Empty circles denote aperiodic regimes. The region of experimentally observed aperiodic regimes is hatched. a


u [VI 10 50b

09 --t

F d T -

1.1

65

d

t

-k I

Figure 1. Phase-transition curves: evaluation of single pulse perturbation. (a) Point A: recording of redox potential for “old” phase 9 = 0.1, positive phase shift of oscillations. Point B: 9 = 0.95, negative phase shift of oscillations. (b) Curve AC, PTC type “0”; curve BD, PTC type “1”. Arrows identify the time of pulse perturbation. Dashed line, value of redox potential chosen as reference event.

The difference between both types of PTC’s is schematically shown in Figure lb. In our experiments with the BZ reaction we defined the reference event as a sharp redox potential jump of the height at least 70 mV exceeding at the same time the level of 1050 mV. Thus random fluctuations of the recorded signal were eliminated. Different choices of the reference event cause a translation of the PTC along the diagonal.” 3.2. Periodic Forcing Experiments. In the continuous forcing experiments both amplitude A of the pulse forcing (amount of added Br- ions expressed as a concentration in the reactor) and the forcing period TFwere varied. The description of the experimental appparatus and of the methods of data acquisition and evaluation were described in ref 8. Using some previously published results8 and new measurements, we have constructed a diagram of forcing period-forcing amplitude, shown in Figure 2. The dimensionless period of perturbation T = TF/ TBis expressed as the ratio of the period of pulse perturbations TFto the average period TBof free oscillations. Resonance regions (corresponding to periodic regimes with the firing number p / q ) and aperiodic regions (hatched) are intermingled in the diagram. Each marked point in the parameter plane represents one measurement corresponding to a record of the redox potential of the length of at least 250T,. Firing numbers p / q (p, number of oscillations; q, number of pulse perturbations), power spectra, and sequences of “periods” were evaluated from the experiments” The time evolution of the redox potential for three chosen regimes (first two periodic, with firing numbers 1/1 and 3/2, and (17) Ruoff, P.; Noyes, R. M. Phase response behaviors of different oscillatory states in the Belousov-Zhabotinsky reaction. To appear in J . Chem. Phys.

u [VI 1.0

0.9

~

-t

-03

150 !

100

50

G

0 i.2 6.4 k-’ Figure 3. Examples of the time evolution of the redox potential, sequences of the “period” 6Tk, and power spectra. (a, b) A = 5 X M, T = l.l7,periodicregimel/l. (c,d)A=5X10-5M,T=1.44,periodic regime 3/2. (e, f) A = 3.5 X lo4 M, T = 2.30, aperiodic regime. (8, h) Power spectra evaluated from sequences of the ‘period”; (9) corresponds to periodic regime with firing number 3/2 (d), and (h) corresponds to chaotic regime (f). I

the third one chaotic) are shown in Figure 3 (the amplitude of perturbations A = 5 X M corresponds to the PTC of the type “1” and the amplitude A = 3.5 X lo4 corresponds to the PTC of the type “0”). The arrows in the upper part of Figure 3a,c,e denote the time of the pulse application. The sequence of “pe.riods” of oscillations (6Tk]illustrates the regularity of the observed regimes. The “period” 6Tkhere represents the time interval between two successive sharp increases of the redox potential (the above discussed characteristic events). We observe that the 1/ l-resonance regime (forcing period TF = 65 s) is characterized by nearly constant values of 6Tk (cf. Figure 3b), which on average are equal


Excitable and Oscillatory Chemical Reaction Systems

The Journal of Physical Chemistry, Vol. 93, No. 7, 1989 2767

to TF Each pulse perturbation causes a negative phase shift (an increase of the “period”), and the resulting oscillatory regime is synchronized on the external forcing frequency. In the case of the 3/2-resonance regime (forcing period TF= 81 s) the same sequence of “periods” was repeated after three oscillations in the course of which two perturbations occurred. For example, in the recording of the redox potential and the sequence of “periods” 6Tk shown in Figure 3c,d, we can observe that the first pulse shown causes a negative phase shift (lengthening of the “period”), the second pulse causes a positive phase shift (shortening of the “period”), and this sequence repeats itself periodically. An example of a chaotic regime is shown in Figure 3e,f. The power spectra of the periodic regime shown in Figure 3d and of the chaotic regime are compared in Figure 3g,h. The presence of the broad-band noise with high power in Figure 3h confirms the existence of chaotic behavior. We note that the PTC’s (cf. Figure 1) can be conveniently used for the interpretation of resulting regimes. All points in the forcing period-forcing amplitude parameter plane corresponding to periodic regimes with the same value of the ratio p/q represent one resonance region (Arnold tongue) in Figure 2. The boundaries between the regions are approximate; their locations could be made more precise by adding further experiments. The comments on resonance regions can be summarized in the following way: (i) The occurrence of resonance regions follows the Farey sequence; only boundaries of regions with periodicities q = 1, 2, and 3 (111,413, 312, 5/3,2/1,7/3,5/2,3/1) could be determined at the current level of experimental noise. (ii) The widths of the Arnold tongues increase with the increasing perturbation amplitude for low values of amplitudes. (iii) The Arnold tongues with the periodicity q = 3 cease to exist for the value of the amplitude A approximately equal to 7 X M, for which the character of the corresponding PTC changes from the topological degree 1 to the degree 0. (iv) The Arnold tongues with the periodicity q = 2 are becoming narrower for higher amplitudes of perturbation A and finally cease to exist for A > 5 x IO-3. (v) Only a one-periodic regime (synchronized with the forcing period) exists for amplitudes A > 5 X The empty circles in Figure 3 denote observed individual aperiodic regimes (characterized by a broad-band noise in the power spectra). Generally, in the absence of noise we expect that the hatched region can contain periodic regimes of higher periodicities together with chaotic and quasiperiodic regimes. The presence of noise (inherent in any experiment) causes both quasiperiodic and higher periodic regimes also to look irregular. The difference between “noisy” periodic of quasiperiodic regimes and chaotic regimes is reflected in different levels of the broad-band noise present in the power spectra. 3.3. Modeling of the Pulse Perturbation of the B Z Reaction. 3.3.1. Simple “Oregonator”. The “Oregonator” model4 in the form appropriate for an isothermal CSTR can be written in the form

+ k3BX - 2 k P - k& = -kiBY - k2XY + fksZ - koY

dX/dt = klBY - kzXY dY/dt

dZ/dt = k3BX - k5Z - koZ

(2)

Here X = [HBrOJ, Y = [Br-1, Z = 2[Ce4+],and B = [BrO,-]; fdenotes a stoichiometric parameter, kl, ..., k5 are kinetic constants, and ko is a reciprocal value of the residence time in the reactor. Field and Noyes4 determined the range of parameters for which oscillations of a limit cycle type exist. In an earlier paper* we reported that for the studied sets of parameters we had not been able to predict positive phase shifts of oscillations. Later we found that by proper choice of values offand k5 the positive phase shifts can be modeled. That is, the kinetic constant k5 depends on the malonic acid concentration,I8 and if its concen(18) Jwo, J.-J.; Noyes,

R. M.J . Am. Chem. SOC.1975, 97, 5422.

n

0

0.5

’p

1

A

0

. 0.5

y

1

0 0.5 y ‘’ 1 0 0.5 CP 1 Figure 4. Phase-transition curves computed for simple “Oregonator” model; (eq 2), amplitudes of perturbation. (a) A = 1 X lo-’ M; ( b ) A = 1 X 10” M; (c) A = 1 X M; (d) A = 1 X lo4 M.

tration is low (in our case 0.05 M), then the rate of formation of bromide ions in the course of reduction of Ce4+ ions is initially slow. The pulse addition of Br- ions into the reactor at this moment accelerates the process of their formation resulting in the positive phase shift of oscillations. In the numerical modeling of the pulse addition of Br- ions we proceed in such a way that at the time of the pulse the course of the integration of the system (2) is interrupted and the value of Y is increased by the amplitude A; see eq la,b. This corresponds to the shift of the phase point in the phase space (X, Y, Z ) by the value A in the direction of the Y axis. The time of the pulse perturbation is determined by the time elapsed from the reference event corresponding here to an intersection of a trajectory with a reference plane in the phase space. The plane Z = 1 X lo4 M ( Z increasing) was chosen as the reference plane, in accordance with the method of evaluation of experimental data. Examples of PTC’s computed for the “Oregonator” model (2) and the values of kinetic constants kl = 18.9 M-’ s-l, k2 = 6 X lo9 M-I s-l, k3 = 3 X lo4 M-’ s-l , k4 = 4 X lo7 M-’ s-I, k5 = 0.53[MA] = 0.0265 M-I s-l, ko = 4.5 X f = 0.6, B = 0.002 M

s-l,

are depicted in Figure 4 for four different amplitudes of the perturbation. Both positive and negative phase shifts occur in the PTC’s. The negative phase shift increases with increasing amplitude A ; the positive phase shifts in the range of phases 0 < cp < 0.1 are approximately the same for A > 1 X IO”. We have observed in the experiments a positive phase shift also for values of phases cp > 0.1. As we have not been able to find such sets of parameters for the model (2) where the agreement between the experimental and modeled PTC’s would be better, we have used an expanded four-variable modeL7 3.3.2. Expanded “Oregomtor”. Ruoff and Noyes7 formulated a four-variable model based on a more detailed description of the formation of Br- ions and the reduction of Ce4+ ions. The corresponding set of equations can be written for CSTR in the form dX/dt = klBY - k2XY

+ k3BX - 2 k 8 - k&

- k2XY + k6Z’I2P - koY dZ/dt = k3BX - kTZ - k0.Z

dY/dt = -klBY dP/dt = klBY + 2kzXY

+ k#

- kSPY - k6Z’J2P - k$

(3)

where X, Y, Z, and B are the same variables as in the model (2) and P denotes concentration of HOBr; kl-k7 are kinetic constants. Ruoff and NoyesI7 tested the validity of the model for the description of pulse perturbations by Br- ions and found that our experimental results reported earlie# can be qualitatively described

2768

The Journal of Physical Chemistry, Vol, 93, No. 7, 1989 a

Dolnik et al.

1

a

-

711

1

0.5

;i

I

A-:

i

-3t

2


I

-5.

0

2.5

251

2 d

I

-

OO

0.5

’p

1

Figure 5. Comparison of experimental PTC’s and PTC’s computed for expanded 4-D “Oregonator”model (eq 3); amplitudes of perturbation. Results of experiments: (a) A = 1.43 X IO-’ M; (b) A = 5 X IO-’ M; (c) A = 1.4 X IO4 M; (d) A = 5 X IO4 M. Results of simulations: (a) A = 5 X lod M; (b) A = 1 X IO-’ M; (c) A = 2 X IOTs M; (d) A = 5 X M.

by the model’s based on the Field-Kiiriis-Noyes mechani~m.~ Still better agreement can be reached if we choose partly different values of parameters from those used by Ruoff and Noyes,” i.e., if we choose a reference cross section 2 = 1 X M (2 increasing) and the following set of kinetic parameters:

k l = 1.34 M-’ S I , k2 = 1.6 X lo9 M-I s-l, k3 = 8 X lo3 M-I s-l, k4 = 4 X lo7 M-’ s-l, k5 = 1 X IO5 M-’ s-l, k6 = 1.0 M-l12 SKI, k7 = 1.0 s-I, ko = 4.5 X s-I, B = 0.1 M

251

2

‘

2.5

3r

3

3

-

-

L 2.5 3

Figure 6. Rotation numbers p and Liapunov exponents X from 1-D model based on PTC‘s expanded ‘Oregonator” model. (a) A = 5 X lo4 M; (b) A = 1 X IO-’ M; (c) A = 2 X lo-’ M; (d) A = 5 X M.

(4) can give chaotic trajectories.” As the single pulse in our experiments affects only the actual ”period”, we can model the periodic forcing by repeated iteration of the corresponding PTC.8 The computed PTC‘s were approximated by fifth-degree polynomials, and Liapunov exponents A and rotation numbers p were computed similarly as in ref 8. The following definitions of h and p were used:

,,--

n

Here f denotes the nth iteration of the function

f ( d = P ~ + =I T + d ( d mod

These values of the parameters do not correspond exactly to modeled experimental conditions. They were chosen to obtain best agreement with the experiments. The resulting PTC’s are for four values of the forcing amplitudes depicted in Figure 5. The circles in the figure denote experimental data. We can observe that the agreement between experimental points and computed PTC’s is particularly good in the cases shown in Figure 5a,b (PTC’s of the type “1”). The results for higher perturbation amplitudes (PTC’s of the type “0”) shown in Figure 5c,d are less satisfactory. It is also important that both the experimental and modeled PTC’s contain regions where

and F ( 4 ) is a lift offlcp). By iterating (7), we obtain a sequence of phases at which the perturbations by Br- ions occur. The value of n equal to lo4 was chosen in numerical computations. The dependence of rotation numbers p and Liapunov exponents h on the period of perturbation for four different values of the amplitude A is shown in Figure 6 . The negative values of h correspond to periodic regimes and the positive values to chaotic ones.19 With an increasing amplitude we can observe the destruction of periodic regimes with high values of p and q and the appearance of chaotic

Idd/d44 ‘1 (4) as it is well-known that an iteration of maps with the property

(19) Collet, P.; Eckman, P.S. Iterated Maps on the Interval as Dynamical Systems; Birkhauser: Boston, 1980.

’

1

(7)



log Y

-4

I -6 -8

-10

2

1

3

T

Figure 7. Parameter plane log A-T results of simulations of 1-D model (eq 7). Periodic regimes form Arnold tongues. Hatched regions, quasiperiodic and higher periodic regimes; cross-hatched regions, aperiodic (chaotic) regimes.

-3,

-5

-51

-7

-7 -

r

0

.. .:

c

-3

200

LOO

t

[SI

.

. ‘. . * ’ e.

-10

-8

-6

log Y


Figure 9. Simulation of periodic pulse perturbations for expanded “Oregonator”model. (a and b) Time evolution of variable Y = [Br-] (a) and Z = [Ce4+](b). (c) Phase plane log Z-log Y. (d) Stroboscopic map. T = 1.35, A = 1.3 X M, aperiodic regime. log Y -6

-3

-3 r

-8

-5

-5

t

-1 0

I

-7

I 0

200

43

t[sl

-10

-8

-6

log Y

Figure 8. Simulation of periodic pulse perturbations for expanded “Oregonator”model. (a and b) Time evolution of variable Y = [Br-] (a) and Z = [Ce“] (b). (c) Phase plane log Z-log Y. (d) Stroboscopic map. T = 1.33, A = 1.3 X M, periodic regime 1/1. regimes; cf. Figure 6c,d. The existence of aperiodic sequences (qk]follow from high values of the slope of the PTC in the region 0 < cp < 0.15 for perturbation amplitudes 2 X M 1 X periodic regime with q = 1 exists. In the hatched region we can observe, contrary to the experimental parameter plane, only periodic regimes with higher values of p and q together with quasiperiodic regimes. Chaotic regimes were found for intermediate values of the perturbation amplitude between the resonance regions with q = 1 and 2 (cross-hatched region in Figure 7). The extent of the region of aperiodic solutions is relatively small compared to the regions of periodic solutions. The boundaries of the resonance regions can be determined accurately by means of numerical techniques based on continuation methods.Iivm However, this can require a significant amount of computer time, particularly if we are dealing with stiff systems such as the four-dimensional “Oregonator”. On the other hand, the approximate method of the location of these boundaries, based on the iteration of the one-dimensional map (7), is easy and very fast; the method can be also used to obtain estimates of starting

-3

1

-3b

L

-5

- 5t

I

I I

-7

Figure 10. Simulation of periodic pulse perturbations for expanded uOregonator”model. (a and b) Time evolution of variable Y = [Br-] (a) and Z = [Ce4+](b). (c) Phase plane log Z-log Y. (d) Stroboscopic map. T = 1.39, A = 1.3 X M, periodic regime 3/2. points for the above-mentioned continuation techniques. The predicted dynamic regimes can be confirmed by the dynamic simulation of the expanded “Oregonator”. Such regimes M in are shown for the perturbation amplitude A = 1.3 X Figures 8-10. The periodic regime 1/1 existing for T = 1.33 is depicted in Figure 8. The stroboscopic map containing one point in Figure 8d is constructed in such a way that in the phase plane log Z-log Y the state of the system immediately before the pulse perturbation is depicted. When the forcing period is increased to T = 1.35 (cf. Figure 9), the character of oscillations changes. The trajectory shown in the phase plane log 2-log Y is in the course of the winding around the limit cycle always perturbed at different positions (cf. Figure 9c). The stroboscopic map then contains large number of points; the irregular course of the trajectory is also evident from Figure 9b. The dynamic regime is chaotic. When T = 1.39 (Figure lo), a periodic regime 3/2 results from the simulations. The stroboscopic map (cf. Figure 10d) contains a pair of points. Hence the dynamic simulations confirmed predictions obtained from the iteration of the map (7), Le., the existence of a chaotic regime between the resonances 1/ 1 and 312. 4. Periodic Forcing of Excitable BZ Reaction

(20) K u b h k , M.; Marek, M. Computational Methods in Bifurcation Theory and Dissipative Structures; Springer: New York, 1983.

4.1. Experiments in a Tubular Reactor. Effects of periodic forcing on wave generation were studied in an essentially one-

2770 The Journal of Physical Chemistry, Vol. 93, No. 7, 1989

5 L

Dolnlk et al.

6 -

3

b

.

-

look

~ " " ' ~ " " ' , T"= &' s" ' , ' 500

1000

v

I

1

1500

t

Figure 12. (a) Periodic regime 1/3; T = 45 s. (b) Aperiodic regime 0.58, 0.05: T = 85 s.

""'5-EuL _.

... 'i

I

1/i

1 1


31

0.5

I I Figure 11. (a) Tubular reactor for wave generation: 1, glass tube (inner diameter 0.9 mm); 2, silver electrode; 3, platinum electrode; 4, electrode vessel; 5, microporous Teflon membrane; 6 , buffer; L, the point where the successive waves are recorded. (b) Wave-train formation. Regime with u = 1/2. An arrow represents a voltage perturbation 6Tk time L

interval between two waves.

dimensional experimental system. A modified BZ reaction mixture (composition 0.225 M HBr03, 0.05 M CH2(COOH)2,0.007 M KBr, 0.00375 M ferroin) in a thermostated (20 "C) tubular reactor (diameter 0.9 mm and length 12 cm) was used; cf. Figure 1 la. A similar reactor was used in experiments on the effects of an electrical field on wave formation and propagation; cf. ref 21a,b. The waves were initiated repetitively at the left end of the reactor by periodically switching the voltage applied between silver and platinum electrodes immersed in the reaction mixture. The electrode vessel 4 is separated from the reactor 1 by a microporous Teflon membrane 5 . The vessel can be operated as a flow-through system, but we have found that batchwise operation did not affect the measurements. N o generation of waves occurred when the silver electrode was negatively biased (-1 V) while a traveling pulse wave was formed when the silver electrode was positively biased (+2 V). The length of the positive voltage pulse was 3 s. Initiation of waves by an immersed silver wire in a Petri dish was studied by Showalter, Noyes, and Turner.22 Waves were observed as thin blue strips moving to the right from the silver electrode. Velocities of waves were approximately 1.5 mm-min-'. The wave velocity for a periodic wave pattern is in principle dependent on the period of initiation and can be ~~,~ times of characterized by a dispersion r e l a t i ~ n . ~Successive the occurrence of waves close to the location L (3 mm from the end of the silver electrode) and the time intervals 6Tk, k = 1,2, between passages of two waves were measured; see Figure 1 lb. The number of generated waves rand the number of applied pulses s within an experimental run were counted. The firing number v was determined as the ratio r/s. A periodic pattern was indicated by the periodicity of the sequence (6Tk). In this case v is equal

.--

(21) (a) SevCikovl, H.; Marek, M. Physica 1983,9D, 140. (b) Sevdkovl, H.: Marek, M. Physica 1984, IJD, 319. ( 2 2 ) Showalter, K.; Noyes, R. M.; Turner, H. J . A m . Chem. SOC.1979, 101, 7463.

(23) (a) P,agols, A.; Ross, J.; Vidal, J. J . Phys. Chem. 1988, 92, 163. (b) Marek, M.; SevEikovl, H. In From Chemical to Biological Organization; Malhus, M., Muller, S . C., Nicolis, G . . Eds.; Springer-Verlag: Berlin, Heidelberg, 1988.

'

0 '

.... 112

.I3

$ * €

Eo I

50

100

Ti51

150

Figure 13. Dependence of firing number v on forcing period T; experiments.

to p / q , where p is the period of {6Tk]and q is the number of pulse perturbations necessary to create p successive waves. Even if local velocities of waves are dependent on the location of the observation point L, the value of the firing number is not affected by an actual choice of L, provided no wave is annihilated. The annihilation of waves can occur close to the point of the wave generation; hence L has to be chosen sufficiently far. Our choice of L satisfied this requirement as it was verified experimentally. An example of two sequences { 6 T k ) ,one with the period p = 1 and the other aperiodic, is presented in Figure 12. The dependence of the firing'number Y on the forcing period T is shown in Figure 13. The full circles denote periodic regimes of the wave trains and the empty ones irregular regimes and regimes with higher order resonances (q > 3). We can observe a typical stepwise dependence similar to the devil's staircase with evident plateaus of 1/3, 1/2, and 1/1 resonances. 4.2. Modeling. A relevant model of a spatially one-dimensional experimental system described in section 4.1 is a periodically forced reaction-diffusion system of partial differential equations. However, the stiffness caused by the kinetics of the BZ reaction implies excessive computer time requirements for more thorough numerical investigations. A non-stiff kinetic model in an excitable mode was used in earlier to study the effects of periodic forcing on an excitable reaction-diffusion-convection system. We have found that the dynamics of generation of wave trains is qualitatively well described in a limit of very strong diffusion26 which leads to a model consisting of a low-dimensional system of differential equations (equivalent to a model of a single well-stirred reaction cell). A physical interpretation of this statement is that every generated pulse wave in the distributed system is represented as an "excited" state in the single reaction cell. We expect that this correspondence will be even better in (24) SevEikovl, H.; Marek, M. Physica 1986, 2ZD, 61. (25) Marek, M.; Schreiber, I.; Vroblovl, L. In Structure Coherence and Chaos in Dynamical Systems; Christiansen; Parmentiers, Eds.; Manchester University Press: Manchester, 1988. (26) Marek, M.; Schreiber, I. In Bifurcation, Analysis, Algorirhms, Applications; Kupper, T., Seydel, R., Troger, H., Eds.; Birkhauser: Basel, 1987.



for t # kT. The limiting value of u = k6(S - AgY) as k6 ~0 is obtained upon differentiating (1 1) with respect to time and comparing with (lob)

dZ/dt = h(Y,Z) = k3BX(Y) - k5Z

Each pulse perturbation applied at times t = kT, k = 0, 1,... can be decomposed into two steps. In the first one, a certain amount of silver ions is added, and in the second one the precipitation reaction immediately restores equilibrium values of Ag and Y. Thus the equations for the forced system can be written as

(8b)

where X is given by

X ( Y ) = {k3B- k2Y

+ [(k3B - k2Y)' + 8klk4BY]1/2)/4k4 (8c)


-

the case of the BZ reaction because of its strongly relaxational character. We intend to carry out the analysis of the reactiondiffu$ion model with the BZ kinetics in the future. Here we present an analysis of a periodically perturbed excitable BZ system in the strong diffusion limit. In spite of its simplicity, the presented model is in good qualitative agreement with experiments. This suggests a kinetics-driven mechanism of the wave-train generation. In our simulations we used the reduced "Oregonator" model of Bar-Eli and Noyes6 with two time-dependent variables dY/dt = g(Y,Z) =jksZ-klBY- k2X(Y)Y (8a)

Here B = [Br03-], X = [HBr02], Y = [Br-1, and Z = [Ce4+]. We have chosen the following values of constant parameters: klB = 0.2, k2 = 2 X lo9, k3B = 1000, k4 = 5 X lo7, k5 = 1, and f = 2. Under these conditions eq 8a,b yield one stable but excitable steady state;6 Le., a small perturbation causes a large-amplitude oscillation before approaching the steady state. Generation of pulse waves is controlled mainly by the concentration of Br- ions.2*22When the concentration of Br- ions is decreased below a certain critical level (approximately 10" M), an excitation of the wave occurs. In the course of the positive voltage pulse a number of electrode reactions (Ag/Ag+, Fe2+/Fe3+, O2evolution) may occur. However, as it was tested experimentally, no wave was generated if the Ag electrode was replaced by the Pt electrode. Taking into account the values of kinetic constants of relevant electrode reaction^^^.^^ and actual concentrations of reaction components in the solution, we can consider only the dissolution of Ag and the subsequent reaction Ag+

+ Br-

AgBr

as a cause of the wave generation. Denoting S the solubility product, the production of silver ions can be described29 by the equation dAg/dt = u(Ag,Y) = k6(S - AgY) (9) where Ag denotes the concentration of silver ions and S = 7.7 X M2 (20 "C). The value of the rate constant k6 varies considerably according to different a ~ t h o r s ;hence ~ ~ , ~we~ have chosen k6 as a free parameter ranging from lo4 M-' s-l t o infinity. Considering the precipitation reaction, the general model for the pulse periodic forcing has now the form (see eq la,b): the pulse perturbation: t = kT, k = 0, 1, ... Y(t+) = Y(t-)

z(t+)= z(t-) Ag(t+) = Ag(t-)

+A

(1Oa)

the autonomous system: t # k T dY/dt = g(Y,Z)

+ u(Ag,Y)

d Z / d t = h(Y,Z) dAg/dt = u(Ag,Y)

(1Ob)

where the forcing amplitude A corresponds to the concentration of added silver ions. Equation 10b exhibits excitability of the same type as eq 8a,b. A special case occurs in the limit k6 m. The precipitation reaction is always in equilibrium and Ag becomes dependent on Y according to the relation

-

Ag = S / Y

+ t4.J

P

(11)

(27) Tanaka, N.; Tamamushi, R. Eleclrochim. Acta 1964, 9, 963. (28) Dvoiik, J.; Koryta, J. Elektrochemie; Academia: Praha, 1983. (29) Ruoff, P.; Schwitters, B. J . P h p . Chem. 1984, 88, 6434. (30) Noszticzius, Z.; Mc Cormick, W. J . Phys. Chem. 1988, 92, 374.

(12)

whence

the pulse perturbation: t = kT, k = 0, 1, ... Y(t+) = Y ( r )

-a

z ( t + ) = z(t-) Ag(t+) = Ag(t-)

+A -a

(14a)

the autonomous system: t # k T dY/dt = g(Y,Z)

+ u.,,(Y,Z)

dZ/dt = h(Y,Z) dAg/dt = um(Y,Z)

( 14b)

Here A represents the concentration of added silver ions and a is a concentration of silver and bromide ions removed by precipitation. The value of a is determined by the equilibrium conditions immediately after the pulse.

s = (Ag(t-) + A - a)(Y(t-)

k k-6

S

u, = --(g(Y,Z)

- a)

(15)

This is a quadratic equation with a unique physically acceptable root u = ((Ag

+ A + Y) - [(Ag + A - Y)' + 4S]'/2)/2

(16)

Equations 14a,b have only two dependent variables, Y and 2. The time evolution of Ag is determined by the course of Y and Z. The forced BZ system described by eq lOa,b (or by eq 14a,b in the infinite precipitation rate limit) has three variable parameters: A , T, and k6. We have chosen two different values of k6, and made k6 = lo4 and k6 = lo6, discussed in the literature29930 an additional choice k6 m to be able to appreciate the sensitivity of dynamics to this parameter. For a fixed value of k6 we constructed the forcing amplitudeforcing period diagrams by calculating firing number for a set of randomly chosen points in the A-T plane. The sharp decrease of Y = [Br-] followed by the sharp increase was taken as the reference event needed for computation of the firing number. Using a Poincare representation of computed trajectories (given by a sequence of states immediately before the pulse), we determined periodic, quasiperiodic, and chaotic regimes. We have found no multiplicities. They certainly exist but in some small-scale regions of the parameter plane. The results of these computations were used to draw approximate boundaries of regions where periodic regimes with a specified firing number exist (Le., p/q-resonance regions); see Figure 14a,b. Regimes other than periodic ones are scattered among the resonance regions. Figure 14a corresponds to k6 = lo4 and is characterized by very narrow resonance regions. Low-order resonances (Y = 1/2, 1/3, 2/3) occupy small portions of the parameter plane. When k6 is increased, the size of the low-order-resonance regions increases at the expense of higher order resonances. The limiting case of infinitely fast precipitation is shown in Figure 14b. Note that the cusps of resonance regions do not extend from the T axis as in the case of forced oscillatory systems but come from infinity at finite values of A . Moreover, values of Y range from 0 to 1 only,

-

Dolnrk et al.

\

a

p0

0.51

I

O C

0 0

0

0 0

0

10

-


I

.

--

20

T

-.

C

.3/&

b

* 213

I

i

112

05r

,

113

0 5

0 1 01 1 i_____L---

1

1

30

20

10

Figure 14. Parameter plane forcing amplitude A-forcing period ulations. (a) k6 = lo4 M-I s-l . (b) kg-+ m. log y - 1 6.16

-

: i

.,*

c

91

T T; sim-

1

0.51

-

I i-

. 0

>

,/’

113

,

0

i

I

0

*

L_____pI_____pL

~

0 5 10 T 15 Figure 16. Dependence of firing number u on forcing period T; A =3 _ _ _ L 1 X M (a) k6 = IO4 M-I s-I (upper dashed line in Figure 14a) (b) 5 10 15 T (c) k6 m (dashed line in Figure 14b) k6 = lo6 M-‘ S-’ Figure 15. Dependence of one-periodic solution on forcing period T. k6 = lo4 M-’ s-l; A = 1.3 X 10” M (this dependence corresponds to the To make a qualitative comparison with experiments, we have lower dashed line in Figure 14a).

-

in contrast to forced oscillators. An important problem is to explain the nature of the boundaries of resonance regions and to locate them accurately. This can be done numerically by using continuation methods.l’~zOWe are currently working on this subject. A preliminary result is shown in Figure 15 as the dependence of period one ( q = 1) resonance on T for k6 = lo4 and A = 1.3 X lod (see the lower dashed line in Figure 14a). The change of stability occurs as two complex conjugate Floquet multipliersz0 go through the unit circle, indicating a torus bifurcation. This bifurcation corresponds to the boundary of O/ 1-resonance region. Thus we have a strong indication of the occurrence of quasiperiodicity at some points in Figure 14a. Boundaries of other resonance regions may correspond to different types of bifurcations.”

also constructed firing number-forcing period plots for all three see Figure 16a-c. values of k6 at a constant value of A = 3 X Clearly, there is no direct correspondence between the parameter A and the amount of silver ions released from the silver electrode during the voltage pulse in experiments. However, we have found that any choice of A yields qualitatively the same results. The dependences shown in Figure 16a-c suggest devil’s-staircase-like structure and Farey ordering. Plateaus of small size occur in Figure 16a while larger values of k, imply larger sizes of plateaus: cf. Figure 16b,c. The last two figures correspond well to the experimental dependence (Figure 13), which suggests that the precipitation reaction is probably very fast. Note that there is a small difference between Figure 16b and Figure 16c; i.e., the dynamics is insensitive to k6 when its values are sufficiently large. The situation in Figure 16a seems to be close to a classical devil’s



.

3


Figure 17. Dynamical simulation of periodic pulse perturbation. (a and b) Time dependence of Y([Br-1) (a) and Ag (b). Phase portraits: (c) log Y-log Z (d) log Z-log Ag, (e) log Y-log Ag. (0 Stroboscopic M, k 6 + -, regime 1/1. (Poincare) map. T = 13, A = 3 X

8 1OC

m

_

8

. 5

- -L 109 AQ

Figure 18. Simulation of periodic pulse perturbations. (a and b) Time dependence of Y(a) and Ag (b). (c) Phase portrait log Y-log Ag. (d) Stroboscopic (Poincare) map. T = 6, A = 3 X M, k6 = lo6 M-’ s-’,

regime 1/3. staircaseI2 with alternating periodic and quasiperiodic dynamics. Aperiodic oscillations indicated in Figure 16b do not seem to be quasiperiodic while no oscillations other than periodic were observed in the limiting case; see Figure 16c. W e anticipate that Figure 16b and Figure 16c correspond to a multivalued devil’s staircase with a mostly hidden chaotic dynamics occurring on small scales in the parameter space. We present some typical examples of dynamical behavior in Figure 17 corresponding to the point with T = 13 and v = 1/1 and in Figure 18-20 which correspond to chosen points from Figure 16a,b. Time evolution of bromide and silver ions in the is shown in Figure 17a,b, three planar projections case of k6 of the trajectory are found in Figure 17c-e, and the Poincare map is in Figure 17f. The temporal evolution is indicated by numbers in Figure 17e; the point 1 corresponds to the state immediately before the pulse addition of silver ions, point 2 is the state after the addition but before the removal of Ag+ ions by precipitation, and point 3 indicates the end of the pulse after which the state of the system changes continuously through points 4 and 5 back to point 1. Except for the pulse perturbation described by the path 123, the points move along the straight line in the log-log plot because the dynamics obeys eq 11. Sequences of dynamical regimes for finite precipitation rates are shown in Figures 18-20. Each figure contains the temporal evolution of log Y and log Ag and the projection of the trajectory and its Poincare map to the log Y-log Ag plane. Pulse perturbations and reference events are clearly indicated as horizontal and vertical pieces of the trajectory in the log Y-log Ag projection plane while the periodicity is best seen in Poincare maps. Thus one can find p/q-resonance periodic regimes directly from these graphs. Approximations of firing numbers of aperiodic regimes must be evaluated from long computer runs. Figures 18 and 19

-

IS0

I

200

.

,

. -8

. . - 5 109 .g-+

Figure 19. Simulation of periodic pulse perturbations. (a and b) Time dependence of Y (a) and Ag (b). (c) Phase portrait log Y-log Ag. (d) M,k6 = lo6 M-l Stroboscopic (Poincare) map. T = 13, A = 3 X s-l, regime 3/4.

Figure 20. Simulation of periodic pulse perturbations. (a and b) Time dependence of Y (a) and Ag (b). (c) Phase portrait log Y-log Ag. (d) M, k6 = lo4 M-’ s-I, Stroboscopic (Poincare)map. T = 9, A = 3 X quasiperiodic regime Y = 0.454.

correspond to two particular resonances with firing numbers u = 1/3 and u = 3/4 from Figure 16b (k6 = lo6). The 1/3 resonance (Figure 18) has a characteristic temporal evolution of the silver and bromide ion concentration. Two successive pulse additions of silver ions are almost immediately precipitated without abrupt changes of bromide ion concentration, but there is a gradual increase of Ag (see Figure 18a,b); the third pulse elicits an abrupt decrease of Y (the reference event) and the overall process repeats. The same behavior was observed for all l / q resonances. The dynamics of the 3/4 resonance (Figure 19) is quite different. Silver ion pulses cause a pause in the rhythmicity of the temporal evolution of Y. Figure 20 corresponds to an aperiodic regime with firing number u = 0.454 ...; cf. Figure 16a (k6 = lo4). The Poincare map of the aperiodic regime (Figure 20d) indicates quasiperiodicity rather than chaos with a highly nonuniform distribution on a closed curve which is, however, regularly covered by successive points. Aperiodicity does not involve a strongly irregular occurrence of reference events.

5. Discussion and Conclusions One forced CSTR and two coupled C S T R ’ S ~ Ican , ~ ~be viewed as the most simple models of forced complex biochemical and biological oscillators. A single periodically forced CSTR with the oscillating BZ reaction represents a one-directional coupling of two oscillators. Generally, phase-transition curves obtained from single pulse experiments can be used to predict even complex behavior of systems with periodic pulse forcing. Experimentally determined firing numbers and computed rotation numbers are then used to characterize the dynamics. Results of simulations based on the (31) Crowley, M. F.; Field, R. J. J . Phys. Chem. 1986, 90, 1907. (32) Dolnlk, M.; PaduSgkovl, E.; Marek, M. J. Phys. Chem. 1987, 91, 4407.

2114

J . Phys. Chem. 1989, 93, 2114-2180

FKN mechanism of the BZ reaction (expanded “Oregonator”) agree surprisingly well with the experimental data. Similar conclusions hold in the case when the system of two coupled CSTR’s with forcing of one of them is i n v e ~ t i g a t e d . ~Both ~ permanent and temporal extinction of oscillations are observed and m ~ d e l e d . ~ ~ . ~ ~ Early experimental observations of complex chemical wave trains generated by periodic perturbations are due to Showalter, Noyes, and Turner.22 They used the BZ reaction on a Petri dish periodically perturbed by voltage pulses on a Pt-Ag electrode couple. Farey ordering of periodic wave trains and a devil’sstaircase-like structure in a similar experimental system were reported in ref 25. We consider the narrow tube used here to be a more controllable experimental system than the Petri dish. It is satisfactory that the results of the experiments reported here

agree well with those published p r e v i o ~ s l y .When ~ ~ periodically forced, both oscillatory and excitable BZ systems reveal striking similarities in the arrangement of resonance regions and in parameter dependences of the firing number. Despite the lack of any apparent periodicity in the unperturbed excitable system, a “latent” periodicity is excited by the forcing. The difference between such a “latent” periodicity and a true one is reflected in the fact that the firing number ranges only from 0 to 1 in the forced excitable system while it may take on an arbitrary value in the forced oscillatory system. There seems to exist already a solid basis of experimental studies on periodic perturbations of chemical oscillations; cf. papers by Rehmus and Ross in ref 2 and the paper by S ~ h n e i d e for r ~ ~a review. We believe that a similar basis should be built for chemical excitable systems. Registry No. BrO,-, 15541-45-4; malonic acid, 141-82-2; ferroin, 14708-99-7.


(33) Dolnik, M.; Marek, M. J . Phys. Chem. 1988, 92, 2452. (34) Crowley, M. F.; Epstein, I. R. J . Phys. Chem., in press.

(35) Schneider, F. V. Annu. Rev. Phys. Chem. 1985, 36, 347.

Regular and Irregular Spatial Patterns in an Immobilized-Catalyst Belousov-Zhabotinsky Reaction Jerzy Maselko,+ John S. Reckley, and Kenneth Showalter* Department of Chemistry, West Virginia University, Morgantown, West Virginia 26506-6045 (Received: July 29, 1988)

Unusual spatial patterns are exhibited in a Belousov-Zhabotinsky reaction in which the ferroin catalyst is immobilized on cation-exchange resin. A thin layer of ferroin-loaded resin beads covered with solution containing bromate, malonic acid, and sulfuric acid exhibits propagating chemical waves for periods in excess of 100 h. The number of spontaneouswave initiation sites increases with increasing concentration of sulfuric acid or bromate, and above a critical concentration only counterrotating spirals are initiated. An overcrowding of these sites at high sulfuric acid or bromate concentrationsresults in irregular patterns with features suggestive of phase turbulence.

Introduction L ~ t h e r l was - ~ apparently the first to recognize that the diffusive spread of an autocatalyst in a chemically reacting mixture may give rise to propagating waves of accelerated reaction. Most studies of chemical waves since this early investigation have focused on initially homogeneous reaction N o investigations of condensed-phase reactions with immobilized catalysts have been reported, although this configuration is cOmmon in many chemical processes. A related configuration is found in surface-catalyzed gas-phase reactions, and pattern formation in oscillatory reactions such as the oxidation of H2 on Pt foil and supported Pt has been investigated.6 We report here studies of spatial behavior in the Belousov-Zhabotinsky (BZ) reaction in which the ferroin catalyst is immobilized on cation-exchange resin. Many factors make an immobilized-catalyst, oscillatory reaction attractive for study. A system in which reaction occurs only at the immobilized catalyst provides an experimental realization of the “pool-chemical” approximation,’ with reactants at almost constant concentrations supplied from a large reservoir. Large differences in the effective diffusion coefficients of the reacting species occur naturally in immobilized-catalyst systems, enhancing the possibility of stationary patterns8-” and spatial chaos.12J3 In addition, an open system is readily configured, offering better control of experimental constraints. Several studies of BZ chemical waves in open reactors have been recently reported. Waves propagating unidirectionally in a ring reactor made of acrylamide gel have been investigated, ‘On leave from the Institute of Inorganic Chemistry and Metallurgy of Rare Elements, Technical University, Wroclaw, Poland.

0022-365418912093-2774$01SO10

TABLE I: System Compositiona [KBr03] = 0.054.50 M [H2S04] = 0.25 M [CHz(COOH)z] = 0.025 M [ferroin] = 1.0 x mol/g resin “200.0 mL of solution mixed with 5.0 g of loaded resin.

TABLE II: Cation-Exchange Resin“ size mesh 1

2 3 4 5

>400 200-400 100-200 50-100 20-50

bead diameter/pm 38-7 5 75-1 SO 106-250 180-425 300-1 180

Bio-Rad Analytical Grade 50W-X4.

where reactants are diffusively replenished from interior and exterior reservoir^.'^ Progress has also been made in the sysLuther, R. Elektrochem. 1906, J2(32), 596. Arnold, R.; Showalter, K.; Tyson, J. J. J . Chem. Educ. 1987.64, 740. Showalter, K.; Tyson, J. J. J . Chem. Educ. 1987, 64, 742. Field, R. J., Burger, M., Eds. Oscillations and Traveling Waves in Chemical Systems; Wiley: New York, 1984. (5) Winfree, A. T. The Geometry of Biological Time; Springer-Verlag: New York, 1980. (6) Schmitz, R. A.; D’Netto, G. A.; Razon, L. F.; Brown, J. R. In Chemical Instabilities; Nicolis, G., Baras, F., Eds.; Reidel: Boston, 1984; pp 33-57. (7) Gray, P.; Scott, S. K. Ber. Bunsen-Ges. Phys. Chem. 1986, 90,985. (8) Prigogine, I.; Nicolis, G. J . Chem. Phys. 1967, 46, 3542. (9) Becker, P. K.; Field, R. J. J . Phys. Chem. 1985, 89, 118. (1) (2) (3) (4)

0 1989 American Chemical Society

Dynamics of Forced Excitable and Oscillatory Chemical Reaction

Recommend Documents