J. Phys. Chem. 1993,97, 10358-10363
10358
Molecular Dynamics Simulations of a Thermochemical System in an Excitable Regime A. L. Kawczyhkift**and J. Goreekit Institute of Physical Chemistry, Polish Academy of Sciences, 01224 Warsaw, Poland, and Department of Chemistry, University of Montana, Missoula, Montana 59812 Received: April 28, 1993’
Results of molecular dynamics (MD) simulations for a model thermochemical system in an excitable regime are reported. The simulated system is very small (600 molecules), and thus internal fluctuations have an important influence on its behavior. If the excitation threshold is large, then the average concentration and temperature obtained in simulations show good agreement with phenomenology. However, as the excitation threshold decreases, MD trajectories exhibit irregular oscillations around the stationary state, rather than fluctuate close to it.
I. Introduction The growing interest in phenomena such as bistability (multistability), excitability, periodic oscillations, and deterministic chaos in chemical systems maintained far from equilibrium and governed by nonlinear dynamics raises up the question how to understand these phenomena on a microscopic level.’ Such understanding is important because it allows investigation of the role of fluctuations, which are completely neglected in the phenomenological description. For example, if the phenomenological dynamics permits coexistence of various attractors, which is typical for a nonlinear system, then fluctuations can play an important role as they may switch a system from one basin of attraction to another. Such strong coupling between fluctuationsand phenomenological dynamics also may be expected in the case of an excitablesystem with a single, globally attracting, stationary state. Small fluctuations around this state are damped by the dynamics, but sufficiently large fluctuations may be amplified to a macroscopic level before returning to the stationary state. In this paper we are concerned with the role of fluctuations in a model of a thermochemical system in an excitable regime. We use a molecular dynamics (MD)technique which allows for the most complete microscopic description of a chemical system. This technique has been applied to the study of bistable chemical systems293 as well as transient and sustained oscillation in thermochemical system^.^.^ In this paper we present a continuation of M D simulations for a thermochemical model, applying them to a system in an excitable regime. The parameters characterizing chemical reactions have been fixed in such a way that the stable stationary state is close to a region in which fluctuations may be amplified to a macroscopic level. A sufficiently large fluctuation forces a single “spike”of temperature before the stable state is approached. The appearance of such spikesis controlled by a threshold level for excitation. Simulations for different values of the threshold level are presented. The dependence of the character of evolution on the threshold level is discussed.
the walls of which are kept at a constant temperature, Tb. The system is composed of the reactant A, the product B, and the catalyst E. The following reactions are assumed to occur: kA
A
+ E-
B
+ E + energy
(la)
kw
B-A
(1b) The first reaction is exothermic with reaction heat UO> 0 and activation energy EA > 0. In M D simulations all reagents are represented by identical hard spheres, and they are distinguished by an identity parameter which has no influence on the mechanical motion. We adopt the standard line-of-center model7 for a thermally activated binary reaction of hard spheres to simulate the process (la). According to this model, those collisions for which the kinetic energy of the relative motion along the line of centers, calculated in the center of masses reference frame, exceeds EA may be reactive. If a reactive collision occurs, then the translational kinetic energy of B and E particles is increased by
uo. We assume in MD simulations that the unimolecular reaction (1b) occurs in collisions between B and the walls of the system. Both the activation energy and the reaction heat areequal to zero for this process. Reaction 1b is supposed to imitate unspecified mechanisms, in which the reactant A is supplied into the system and the product B is removed. Let us notice that the energetic balance of the cycle A B A is positive because the system is open. In the phenomenological description of the system (l), assuming ideal stirring, the rate constants for (la) and (1 b), kA and kw, depend on the temperature T in the following way:
--
( &)
kA = kAoTI/*exp --
11. Phenomenological Model of Excitability in a
Thermochemical System In order to have an excitable system, we have adjusted the model describing sustained thermochemical oscillations in a nonadiabatic system.5 This model is a modification of that of Sal’nikov,6 and it can be directly used in MD simulations. Let us consider a well-mixed thermochemical system which exchanges energy with its surroundings. It is enclosed in a cube, t Polish Academy of Sciences. 8 University of Montana.
Abstract published in Aduance ACS Absrracrs, September 1, 1993.
where kg is the Boltzmann constant whereas kAo and kwo are constants depending on the collision frequency. The square-root dependence on T i n eqs 2a and 2b describes the change in the binary collision frequency with temperature, whereas the exponential factor in eq 2a comes from the Arrhenius equation. The concentration of the catalyst E does not change in time in the system considered. Thus [E](t) = [E](t = 0) = [Elo, and the sum of concentrations of A and B remains constant, [A](t) + [B](t) = [A](O) + [B](O) = [NIo. The composition of the system is then determined by the concentration of A only. The
0022-365419312097-10358$04.00/0 0 1993 American Chemical Society
The Journal of Physical Chemistry, Vo1. 97, No. 40, 1993 10359
Model Thermochemical System in an Excitable Regime kinetic equation for [A] reads
46
0'120
( ):
= kwoT'/Z([N], - [A]) - kT'/* exp - k T[A]
Q
(3b)
where k = kAo[E]o. The system exchanges its energy with the surroundingsdue to a heat flow through the walls, which are kept at the constant temperature Tb. This process may be described by Newton's law, and thus the time evolution of temperature is given by the equation
where cvis the heat capacity of the system at a constant volume (it is assumed that the volume does not change during the evolution). It is convenient to introduce the parameters p = UO/CV and j3 = x0/Uo (see eq 5 ) , and then eq 4a becomes dt
k[A] T'/2 exp( -
$) -
@T'l2(T - Tb))
(4b)
The parameter p gives the increase in the average temperature of the system caused by a single reactive collision between A and E. The definition of the parameter j3 comes from the fact that in a hard sphere model the exchange of energy with the surroundingsis connected with collisions between molecules and the walls and therefore the coefficient x should include the factor T V Z (compare eqs 2a and 2b):
x = x0T'l2 (5) Let us define dimensionless variables: concentration a = [A]/[N]btemperature[= T/Tb,andtime~= kTb'l2t. Equations 3 and 4 can be transformed in these variables to the form da d7 = 5
-(1 k
-a) -exp
where t = EA/(kBTb). Equations 6a and 6b describe an open system (CSTR) with inflow-outflowterms kwo/k( 1 -a)and j3/k([- 1). These terms are the same as those in the original Sal'nikov model, except for the term [1/2, which has no effect on the nullcline positions but of course does affect the rate of reaction and energy exchange. The nullclines of equations6a and 6b are helpful for qualitative understanding of the behavior of the system. For eqs 6a and 6b they are
(7)
It may easily be shown that two extremes of CYTexist for B > 4 and that the nullcline CYThas an N-like shape (see Figure 1). This is the necessary condition for the appearance of both oscillations and excitability, because otherwise the vector field for 5 is contracting for all values of a and 5 on the phase plane. The right and the left branches of the N-shaped O ~ Tnullcline are attracting, while the middle branch is repelling. If the point of intersection of the nullclines lies on the repelling branch of (YT, then the system exhibits sustained oscillations of concentration and temperature.5 Otherwise the stationary state is stable;
Lc
0.080
c
.0 _ +,
F
c
K a, 0
5 0.040 0
0.000 0 0
20
4.0
60
F
tempera t u re
Figure 1, Nullclines a ~ ( 5and ) (IT(€) (see eqs 7 and 8) on the phase plane I, a. The four nullclines aA shown correspond to k ~ " T b = ~ 8.25 / ~ X 10' s-I (the bottom one), 1.24 X 106 s-l, 1.65 X 106 s-I, and 2.48 X 106 s-I (the top one), respectively.
however, the system still may exhibit excitable behavior. In this paper we consider the case in which the stable state lies on the left attracting branch of CYT(compare Figure 1). The behavior of the system is controlled by the deterministic dynamics and by thestrengthof fluctuations. If fluctuationsgenerate a state which is below the repelling branch of aT,then the deterministicdynamics directs its evolution toward the stable state. However if a state appearing as the result of fluctuations is above the repelling branch, then the phenomenological dynamics forces a single spike of temperature before the stable state is regained. The strength of the vector field in the [direction on the phase is controlled by the parameter M. For p a,temperature [ is a fast variable as compared to a. Then each trajectory starting at a point located above the repelling branch of (YT,or above the left attracting branch, but with an initial a! larger than the local maximum value, jumps to the right attracting branch of the nullcline LYT. It follows CYTto its minimum, where it jumps again to the left attracting branch and finally approaches the stable state. Such a trajectory would be composed of regions characterized by high and low temperature with instantaneous "jumps" between them. With decreasing p, the influence of the vertical vector field associated with the concentration a on the trajectory becomes more important. The parameters of our simulationsare chosen such that p = z/lsTb, which means that [ is not a fast variable. To study the effect of excitability in system 1, it is convenient to use kw0 as a parameter controlling the location of the stable state. It follows from eqs 7 and 8 that changes in kwo have no influence on CYT,whereas CYAis shifted upward when kwo increases (see Figure 1). To perform numerical simulations, four values of kwo have been selected in such a way that for the smallest value the excitation threshold is large, whereas for the largest one the threshold may be easily crossed by a fluctuation.
-
III. Results of Molecular Dynamics Simulations The reactants are represented on the microscopic level by hard sphereslabeled as different chemical substances.* The diameters and masses of all spheres are assumed to be the same for simplicity The calculations are performed for the system of 600 spheres with diameter 5 %Landmass 32 mass units. Thus the total heal capacity Cvis 900ke. The spheres move inside an impenetrablt cube of edge length 63 A, which correspondsto a packing fractior of 0.157. An additional parameter assigns chemical identity tc a sphere. Twenty spheres represent the catalyst E.
Kawczyfiski and Gorecki
10360 The Journal of Physical Chemistry, Vol. 97, No. 40, 1993 0.13
0.13
0.1 1
0.1 1
a
6 v-
c
0
0
0.09
C
L
I 0
.
0.09
.-0
]
Y
0
0
c i
4 i
5 0.07
5 0.07
0
0 C 0 0
c
0 0
0.05
0.05
0.03 0
E 3.0 a,
:
E3.0 a, Y
4
m x m
m x m
-
1
100
200
300
400
time /(ns)
I
a,
r
2.0
-
c
0
a,
L
3
z 1.0
4
a,
a S
i
0.0 0
100
200
300
400
time / ( n s )
0.12
1 6
6
L
0.08 -
C
.-0
0 .-c
2
Y
Y
2 c
Y
C a,
a, 0
0
0.04
-
0
E 0.04
-
0 -I
1, 0.5
0.00
temperature Figure 2. Molecular dynamics simulations of system 1 for kwo T&/Z = 8.25 X IO5s-': (A, top) a as a function of time; (B, middle) [as a function of time; (C, bottom) MD trajectory on the phase plane [,a (the dashed line). The solid lines show CYTand CYA.
Of course these simulations are not intended to describe any realistic system, but they do allow us to introduce a microscopic model correspondingto the phenomenologicalequationsdescribed
I
I
\
, I
,
,
~
1.0
I
,
I
1
1
1.5
I
,
I
I
1
2.0
I
,
,
I
I
2.5
I
I
,
I
,
3.0
I
I
I
,
]
3.5
temperature Figure 3. Molecular dynamics simulations of system 1 for kw0Tb1/2 1.65 X 106 S-I: (A, top; B, middle; and C, bottom) as in Figure 2
.
above. The mechanical definition of a reactive collision (lineof-center model') and the parameters describing the rates of reactions and heat transfer may be chosen in such a way that direct comparison with phenomenology is possible. The micro-
The Journal of Physical Chemistry, Vol. 97, No. 40, 1993 10361
Model Thermochemical System in an Excitable Regime
c 0
0.09
c
0 ._
z
Y Y
$ 0.07 0
c
0 0
0.05
0.03
0
200
100
300
0
400
200
100
0.12
4
a
a c O 0.08
I
0.08
c
C
.-0
0 ._
Y
Y
Y
z
4
F
c a, 0
400
i
0.12
c O
300
time /(ns)
time /(ns)
c
a, 0
8 0.04
0.04 -
0
0
0.00
LL-khA
0.5
1.0
1.5
2.0
2.5
3.0
3.5
temperature
I
-1
\
I
I
0.00
I 1
I 1
0.5
1
'
'
I
'
1.0
)
I
l
l
'
1.5
)
(
'
I
I
2.0
I
I
'
I
I
2.5
'
'
1
I " , '
3.0
3.5
temperature
Figure 4. Molecular dynamics simulations of system 1 for &w0Tb1/*= 1.24 X 106 s-l: (A, top left; B, top right; and C,bottom left) as in Figure 2; (D,bottom right) phenomenological trajectory of evolution toward the stable state. The initial values are (2.0, 0.06), (2.5, 0.06), and (3.0, 0.06).
scopic model allows for fluctuations, which are not considered in the phenomenological picture. The reaction parameters for which the system exhibits excitability must allow us to carry out molecular dynamics calculations within reasonable computer time. These criteria lead to the values of parameters we have chosen. We use e = 5 as being significantly larger than the critical E = 4 required for the existence of excitability. On theother hand, for E = 5 reactions are still frequent enough to be modeled by the molecular dynamics technique in reasonable computer time. The value of &c,O is fixed by assuming that the probability of reaction is equal to 0.5, provided that the energetic condition for a reactivecollision is met. We selected UO= 120k~Tb (P = 2/ IS Tb) because good agreement between phenomenologicalbehavior and the simulationswas obtained for this value in a regime exhibiting sustained oscillations.5 Collisions with the wall are responsible for reaction l b as well as for the exchange of energy between the system and its surroundings. Both elastic and inelastic collisions of all molecules with the walls are allowed. In an elastic collision, the velocity perpendicular to the wall is simply reversed. The parallel component of velocity has the precollision modulus, but its direction is chosen at random from the uniform distribution of spherical angles. Such a definition of an elastic collision helps to avoid problems with velocity correlations caused by the shape
of the box. In an inelastic collision, the perpendicularcomponent of velocity is randomly chosen from the distribution of velocities of molecules at temperature Tb. The parallel velocity component is chosen in the same way as in an elastic collision. The fraction of inelastic collisions is controlled by the parameter U I (0 Iul I l), which of course determines the value of x in the Newton's law (eq 4). Reaction 1b occurs in inelastic collisions of B with the walls only. The molecules of A appear in the system with the average energy correspondingto the temperatureof the walls. This reflects the nonisolated nature of the system. The parameter u2 (0 Iut I1) controlsthe fraction of inelastic collisionsof B with the walls which lead to reaction lb, and it controls kw0 if U I is fixed. Only u2 is varied here. The values of xo, kwo, and k are determined by the geometry of the system (the size of box, the radius of spheres, and their density) and by the parameters for reactive collisions and thermalization. Half of the collisions between A and E,for which the energetic condition is satisfied, are reactive, so k~OTbl/~[El = 1.4 x 109 s-1. For u1 = 0.05 the value of /3 is (107 s-l)Tb-3/2. To study excitability, we have chosen the following four values of u2: 0.0015,0.00225,0.003, and 0.0045, which correspondsto k ~ ' T b l / ~ 8.25 = X 105S-l, 1.24X 1O6Ss-', 1.65 x 106S-',and2.48 X 106 s-1. The corresponding changes in the position of the nullcline CYAare shown in Figure 1 and also in Figures 2-5C.
Kawczyfiski and Gorecki
10362 The Journal of Physical Chemistry, Vol. 97, No. 40, 1993
i
0.1 1
3.0
e,
Y
4
m x m
Q L
0
e, -c
0.09
C
4-
.-0
2.0
L
Y
0
F
?!
Y
0.07
3
u
0
?
C
a, 1.0
0
0
fi
0.05
0.03
0.0 0
time / ( n s )
100
time200 /(ns)
300
4
I
-
0.120
0.12 -
Q
Q
' i
L
O C
0.08
I
E 0.04
-
0.080 1 ]
.-0
C
.Y_ 0
Y
F
O i u C
4
c
e,
a, 0
0
0.040 0
0 1
0.000
i
0.5
I
I
I
I
I
1.0
1 I
I
I
,
I
1.5
1
,
,
,
1
2.0
1
I
I
I
1
2.5
temperature
,
I
,
I
I
3.0
,
I
, , 3.5
0.00
~
0.5
1.0
1.5
2.0
2.5
3.0
2
temperature
Figure 5. Molecular dynamics simulations of system 1 for = 2.48 X 106 s-l: (A, top left; B, top right; and C, bottom left) as in Figure 2; (D, bottom right) phenomenological trajectory of evolution toward the stable state. The initial values are (2.0, 0.08), (2.5, 0.08), and (3.0, 0.08). Quantitative differences between the MD simulations and the TABLE I: Values of Stationary Temperatures t,,, Stationary Concentrations Q Average Temperatures J;, and Average phenomenologicalmodel become more important with increasing Concentrations a , Together with Their Dispersions for kwoTb1/*through the sequence 1.24 X lo6 s-l, 1.65 X lo6 s-l, Different Values of uz 2.48 X 106 s-l (see Figures 3-5). Fluctuations around the no. uz Lt aIt CaV aav ((Afl2)"* ( ( A U ) ~ ) ~ / ~stationary states grow and drive the system away from the left 1 0.0015 1.075 0.0562 1.150 0.0568 0.123 0.0049 attracting branch of CYT. The MD trajectories exhibit irregular 2 0.0025 1.111 0.0713 1.167 0.0591 0.186 0.0117 oscillationsaround the stationary state rather than being uniformly 3 0.0030 1.146 0.0919 1.184 0.0563 0.192 0.0095 attracted toward it. For 2.48 X 106 s-l (Figure 5C) fluctuations 4 0.0045 1.216 0.0942 1.227 0.0707 0.232 0.0222 are large enough to switch the system across the repelling branch of CYTto the region where the component of the vector field Each simulation starts with CY = 0.057 (which corresponds to directs the trajectory toward the right attracting branch of CYT. the initial number of A and B molecules equal to 33 and 547, The differences between MD and phenomenology cannot be respectively)and theinitial temperature t = 1.2. The temperature for hard spheres is calculated from the total kinetic energy of all explained ifone assumes that all fluctuationsaround thestationary spheres. state are uniformly damped. In the phenomenologicaldescription it is assumed that reactive collisions are sufficiently rare that the The results of MD simulations are shown in Figures 2-5A, B, system has enough time to reach local thermal equilibrium and and C. Table I contains time average values of [ and CY,their homogeneity. One may expect that such conditions are satisfied dispersions, and their phenomenological stationary-state values. for systemswith a very large activation energy and a small reaction If kW0Tb'/' = 8.25 X los s-1, (see Figure 2) the MD trajectory heat. Moreover a mechanism for ideal stirring should operate. fluctuates around the stationary state and the average values of This is not the case here, as the activation energy is small (e = [ as well as CY are close to the stationary values predicted by 5), the reaction heat is relatively large (UB= IZOkBTb), and any phenomenology. The system remains near to the stationary state, forced stirring is absent. Any reactive collision between A and especially if one takes into account that only 600 particles are involved. Fluctuations are not able to drive the system away E releases 60k~Tbof kinetic energy per colliding molecule, which from the left attracting branch of CYT(see Figure 2C). is about 1 order of magnitude more than the average energy per
Model Thermochemical System in an Excitable Regime sphere. Many collisions therefore are necessary to dissipate the energy released throughout the whole system.First let us focus our attention on the average values of a observed in simulations,which are lower than those predicted by phenomenology, and the difference increases with increasing threshold height. In previous simulations of sustained thermochemical oscillations,5we observed that the exchange of energy with the surroundings is more efficient for an unequilibrated velocitydistribution, so thevalueof @islarger than theequilibrium 6 = (8 X lo* s-1)Tb-3/2. The increment in 6 shifts the nullcline (YT upwards, whereas the position of the nullcline CUAremains unchanged. Thus the stationary state should be characterized by an average a larger than the one predicted by phenomenology, which is just opposite to what is observed in MD simulations. The effects observed may be explained if we consider spatial inhomogeneities in our system. Fluctuations may create a region - a “hot spot”, which at some moment is richer in A. Therefore, in this part of the system, a local (Y may reach a value above the repelling branch of (YT,although the average value of a for the system as a whole is below it. Then the phenomenological dynamics forces a spike of temperature, and in consequencecauses a drop in a,before the left branch of O ~ Tis reapproached. This assumption of strong inhomogeneities explains why spikes of temperature may start at a below the local maximum of CYT. Further evolution of the system agrees quite well with the predictions of deterministic dynamics, if one takes into account “random initial conditions” generated by fluctuations. Figures 4D and 5D illustrate the relaxation from an initial state, determined by and a,which may be reached by fluctuations. Notice that the minimumvaluesof 4 and a seen in MD simulations agree well with the minimum values coming from the phenomenological equations. It is very difficult to estimate the peak temperature because the size of a subsystem undergoing yexplosion” (a hot spot) is related to the strength of fluctuations. Similarbehavior may be expected in any excitable system, because fluctuationsmay always create regions in which the deterministic dynamics forces evolution toward another attracting branch of nullcline. However if a system is well mixed and if it is regarded as an infinite one (thermodynamical limit), then the size of such a region is negligible if compared with the system as a whole and its contribution to the average values may be neglected. Here the size of a region undergoing fluctuation is of the same order of magnitude as the whole system, and therefore the coupling between fluctuations and excitable dynamics is clearly visible. Let us notice that the average values of [, and especially a, deviate significantly from their phenomenological stationarystate values for all higher values of kw0Tb1/2(see Table I, rows 2-4). Their dispersions also increase when the threshold barrier decreases. This is a consequence of the fact that MD trajectories spend significant periods of time away from the stationary state (see Figures 3-5).
The Journal of Physical Chemistry, Vol. 97, No. 40, 1993 10363
IV. Conclusions This paper is concerned with the influence of activation threshold on thedynamicsof an excitable, thermochemical system. The main reason for our study is a comparison of MD simulations with phenomenological description, especially the question of whether fluctuations can cause the system to enter its excitable region if the activation threshold is small. Being limited by the computer facilities available, we have been forced to perform molecular dynamics simulations for a system composed of 600 molecules only. For such a small system fluctuations have important influence on its dynamics. In the case of large activation threshold, the MD simulations are in quite good agreement with the phenomenological model. The behavior of the system is entirely different if the activation threshold is lowered. Local fluctuations may be enhanced by the dynamics, leading to “hot spots”, which belong to the excitable region. In this case, because of nonequilibrium effects related to spatial inhomogeneities,the usefulness of the phenomenological model based on equilibrium thermodynamics is limited. A better agreement between MD simulations and the phenomenology may be achieved for a larger system in which the relaxation processes are more effective. It requires more molecules, larger values of activation energy EA,and lower heat of reaction UO, but these requirements increase computer time demand. Our simulations suggest that for accurate description of an excitable system one has to take into account the relaxation processes related to the motion of particles and collisions. Therefore it is necessary to develop a statistical descriptionof far from equilibrium systems, which should be much richer than the one used in contemporary phenomenological theories. Acknowledgment. The authors would like to thank the University of Montana Computer Center for technical support and computer time. A.L.K. is very grateful to Prof. R. J. Field, University of Montana, for hospitality which made this work possible. References and Notes (1) Nicolis, G.; Rigogine, I. Self Organizafion in Non-Equilibrium Systems; Wiley: New York, 1977. (2) Ortoleva, P.; Yip, S. J . Chem. Phys. 1976,65, 2045. ( 3 ) Boissonade, J. Physica l982,113A, 607. (4) Gorecki, J.; Kawczyhki, A. L. J . Chem. Phys. 1990,92, 7546. ( 5 ) Kawczyfiski, A. L.; Gorecki, J. J . Phys. Chem. 1992,96, 1060. (6) Vol’ter, B. V.; Sal’nikov, I. Ye. Modelling and Oprimization of Curalytic Processes (in Russian); Nauka: Moscow, 1965; p 128. (7) Xystris, N.; Dahler, J. S. J . Chem. Phys. 1978, 68, 387. (8) Beme, B. J. Statistical Mechanics;Plenum Press: New York, 1977. Allen, M. P.; Tildesley, D. J. Computer Simulution of Liquids; Oxford University Press: Oxford, U.K. 1987. (9) Dong-Pa0 Chou; Yip, S. Combust. Flume 1984,58, 239.
