J . Phys. Chem. 1989, 93, 2811-2816
2811
Molecular Fluctuations on Limit Cycles and Chaotic Attractors Joel Keizer* and James Tildent Institute of Theoretical Dynamics and Department of Chemistry, University of California, Davis, California 9561 6 (Received: August 19, 1988)
We continue our examination of molecular fluctuationson nonstationary trajectories using the statistical theory of nonequilibrium thermodynamics. We generalize previous work near equilibrium and steady states to orbitally asymptotically stable limit cycles and chaotic attractors. On limit cycles we show analytically that the conditional covariance matrix grows linearly with the number of periods after an initial transient. The corresponding time-average covariances on chaotic trajectories are shown by numerical calculationsto grow exponentially with a time constant determined by the largest Liapunov characteristic exponent. This latter result appears to be generic, and both results are illustrated with numerical calculations for the Rossler equations including fluctuations. Implications of these results for macroscopic measurements are discussed.
I. Introduction A great deal has been learned about the nature of molecular fluctuations in large systems over the past two decades.',2 Near equilibrium the linear theory of Onsager and Machlup3 and Fox and Uhlenbeck4 provides a thermodynamic description of fluctuations in extensive variables. Away from equilibrium that theory has been generalized in a way that permits molecular fluctuations in a great variety of macroscopic systems to be treated in a unified fashione2 This theory is valid in the thermodynamic limit and has been used successfully to analyze a variety of experiments, ranging from electrical noise near instabilitiesS to light scattering in temperature gradients6 and concentration fluctuations in bistable systems.2 The theory can also be used to treat fluctuations due to collisions for the Boltzmann equation and has been derived in various ways from both Hamiltonian' and kinetic descriptions.'J The theory allows one to calculate conditional fluctuations in the extensive variables rather simply if the deterministic average trajectories are known. The calculations are especially easy in the neighborhood of asymptotically stable steady states, where one finds a natural generalization of the Onsager-Machlup theory.9 It is also possible to use the theory to calculate conditional or unconditional averages and fluctuations around these averages for nonstationary ensembles2 These calculations are of interest because one knows that nonlinear effects increase dramatically the size of both single-time and conditional fluctuations around nonstationary trajectories.2 Experimentally, such effects have been seen in the Gunn oscillator,I0 spinodal decomposition,I' laser instabilities,I2 and bistable chemical systems.13 In this article we extend our studies of fluctuations on nonstationary trajectories to include average trajectories that are either limit cycles or chaotic attractors. We are especially interested in the question of whether or not fluctuations around such trajectories become large enough to interfere with normal measurements of the average. For trajectories that begin in the domain of attraction of an orbitally asymptotically stable limit cycle, we show that, following an initial transient, the covariance matrix can be written as the time intergral of a periodic function. As a consequence, the asymptotic behavior of the covariance, when evaluated at multiples of a period, increases linearly with the number of periods of oscillation. This implies that the fluctuations diverge, on average, with the square root of the time. Nonetheless, because fluctuations in the extensive variables are proportional to the square root of the system size, while the average of the extensive variables is proportional to the system size, itself, the fluctuations will not become comparable to the average for macroscopic systems on the laboratory time scale-except, possibly, close to critical points. Comparable arguments for systems with chaotic attractors suggest that the covariances of the fluctuations should increase exponentially with time, with a characteristic time constant equal to twice the largest Liapunov exponent. These predictions are examined numerically using the Riissler equations,I4 'Permanent address: 33/4 M U l , Patong Kathu, Phuket 83121, Thailand.
0022-365418912093-2811$01.50/0
extended to include fluctuations.2 As predicted, we find a linear divergence of the covariances for parameters in the regime of limit cycles and an exponential divergence of the covariances, suitably averaged, in the parameter regime of chaos. The implications of these results for measurements on limit cycle oscillators and chaotic systems are discussed. 11. Fluctuations on Nonstationary Trajectories
The basic formulas that describe fluctuations in the extensive variables in the thermodynamic limit are well-known.2 Although they have been derived in various ways from other assumpt i o n ~ , ' ~we~ will ~ ~ ~not* be concerned with that here, but rather will take the equations in this section as our starting point without further comment. Let the column vector n(t) represent a realization of a vector-valued stochastic processes whose elements are extensive variables. Two-time conditional averages will be represented by ii(no,to,t) '
lP2(no,t01n,t)n dn
where P2(no,t01n,t)is the two-time probability density for n at time t , conditioned on a precise value of n = no at time t o . In the thermodynamic limit the conditional average satisfies deterministic equations of the form dii/dt = Cwk(VK+- VL) S(ii,t) R(ii,t) (2)
+
K
where ii(no,to,tO)= no. The first term represents2 a sum over the forward and reverse rates of elementary molecular processes ( Vz*) evaluated at ii, and S represents systematic, nondissipative contributions to the time rate of change. Conditional fluctuations are defined by 6n = n - n and, in the thermodynamic limit, are Gaussian with a covariance matrix2
that satisfies (1) van Kampen, N. G. Stochastic Processes in Physics and Chemistry; North-Holland: Amsterdam. 1981. (2) Keizer, J. Statistical Thermodynamics of Nonequilibrium Processes; Springer-Verlag: New York, 1987. (3) Onsager, L.; Machlup, S.Phys. Reo. 1953, 91, 1505. (4) Fox, R. F.; Uhlenbeck, G. E. Phys. Fluids 1970, 13, 1893. (5) Keizer, J. J . Chem. Phys. 1981, .74, 1350. (6) Kirkpatrick, T. R.; Cohen, E. G. D.; Dorfman, J. R. Phys. Reo A 1982, 26, 950. (7) Kirkpatrick, T. R. On the theory of light scattering from fluids in nonequilibrium steady states. Doctoral Dissertation, Rockefeller University, 1981 ...
(8) Keizer, J. J . Math. Phys. 1977, 18, 1316. (9) Keizer, J. J. Chem. Phys. 1976, 65, 4431. (10) Kabashimi, S.; Yamazaki, H.; Kawakubo, T. J. Phys. Soc. Jpn. 1976, 40, 921. (11) Elder, K.; Rogers, T.; Desai, R. Phys. Rev. E 1988, 38, 4725. (1 2) Arecchi, F. T. In Pattern Formation in Dynamic Systems and Pattern Recognirion; Haken, H., Ed.; Springer-Verlag: Berlin, 1979; pp 28-42. (13) Kramer, J.; Ross, J. J . Chem. Phys. 1985, 83, 6234. (14) Rossler, 0. E. Phys. Lett. 1976, 57A. 4701.
0 1989 American Chemical Society
2812 The Journal of Physical Chemistry, Vol. 93, No. 7, 1989
du/dt = Ha
+ a# + y
Keizer and Tilden
(4)
with u,(to,to) = 0 and the superscript T representing the transpose, and where
HI, = (dRl(ii,t)/dB,)
(5)
is the Jacobian matrix based on the conditional average and y is the covariance matrixZ Y,, =
C % ( ~ #++vK-)u,K
(6)
5
x
The solution of (4) is a(t",t) = &'K(t',t)
y(to,t? f l ( t : t )
dt'
(7)
0
2
4
where the matrix K(t',t) solves
7 =
dK/dt = H K
(8)
with K(t',t') = I , the identity matrix. The probability distribution in any conditional ensemble is completely determined by (2)-(8). Furthermore, the single-time probability density, Wl(n,t), can be obtained from the conditional density using the relationship2 Wl(n,t) =
1Wl(no,to) Pz(no,toln,t) dno
(9)
where Wl(no,to)is the single-time density at time t o . Using the above equations one can obtain a useful expression for the single-time covariance2 (6n(t) 8nT(t)), 5 (n(t) nT(t))l - (n(t))l(nT(t))l
(10)
+ fWl(no,to)
(6n(t) 6nT(t)), = ~ W l ( n o , t ou(t",t) ) dno
-
s
Wl(no,t") ii(t) dno
s
8
10
At
Figure 1. Scaled conditional variance of the number density fluctuations of B molecules for the reaction A + B s 2B. The initial concentrations of B (in arbitrary units) are shown as pea. Other parameters are p = 1, k+ = 3, k- = 1, and X = 3 .
value pB' Since pA + pB = p = constant, one has (6pA(t)2)" = - ( 6 p A ( t ) 6 p B ( t ) ) " = (SpB(t)*)" 1 u(t",t). Thus to characterize the conditional fluctuations, it suffices to find u(t",t), which according to (4) and (8) solves du/dt = 2 [ k + p - 2(k+ + k - ) & ] U V-"l[k+pB(p - De) k-p~'] (15)
+
+
Equation 15 shows that u scales like VI,and it is sensible, therefore, to define a = u/u* with u* proportional to VI,which then solves an equation of the form da/dt = h(t)a g ( t ) (16)
+
where (...)l represents J...Wl(n,t) dn. One finds that2
ii(t) iiT(t) dn"
6
X
Wl(no,to) iiT(t) dn" (11)
The first term here represents the contribution from the conditional covariance averaged over the initial probability density, while the second term represents the propagation of the initial covariance due to the dependence of the conditional average on the initial condition, no. Because the rates of the elementary processes, Vx*, are proportional to the system size2 (e.g., the volume, V), the conditional fluctuations are often negligibly small with respect to the average. In fact, in stationary ensembles or ensembles that evolve to stationary ensembles, the conditional covariance is many orders of magnitude smaller than the covariance due to the experimental dispersion in initial conditions, except near critical points.2 To help illustrate these facts and to set the stage for our treatment of fluctuations on limit cycles and chaotic trajectories, we close this section with a simple example. Consider the elementary autocatalytic chemical reaction
where 8 ( t o , t o ) = 0 and h(t) and g(t) can be gleaned from (15). Equation 16 is a nonautonomous scalar differential equation which is easy to solve either analytically or numerically. Numerical solutions, using the PLODIS implementation of the Gear algorithm,16 are given in Figure 1 with u* = ue, the equilibrium variance. Notice that u remains of the order of its equilibrium value 6= ( p - pBe)pBe/pV,except transiently for initial conditions near pBo = 0. Notice, also, that ue112/pBeis of the order of (PV)-'/~. Thus unless the volume contains only a few A and B molecules, the conditional fluctuations can be neglected. Thus for a macroscopic sample the conditional fluctuations in this example will be negligible with respect to macroscopic uncertainties in the initial condition, which are governed by Wl(pB",t). This implies that the second terms in (1 1) will dominate the variance for long periods of time, as can be verified explicitly by approximating P2(no,t01n,t)by the delta function 6(n-ii(no,to,t)). Doing this, one finds that (9) can be integrated to give2 W l(n,t) = Wl(no(n,t),tO)ldnO /en1 (17)
(14)
where n"(ii,t) = no is the inverse of the mapping between no and the conditional average and ldno/dnl is the absolute value of the determinant of the Jacobian of the inverse mapping. The difference in size of the single-time and conditional variance can be seen by comparing Figure 1 with Figure 2. In Figure 2 the single-time probability density is shown at four times based on (17) using an initial Gaussian centered at pB/pBe = 0.1 and an This initial initial variance of (6pB(0)2)11/2/pBe= 3 X variance is macroscopically large, large enough, in fact, to effect profoundly the average, as can be seen from the skewing of the single-time density as time proceeds. Notice, also, that the initial spreading of the probability density is not a consequence of a diffusive term but only the different drift rates of the initial points. For very long times the single-time covariance calculated in this fashion becomes comparable to ue, and by then the effect of the first term in (1 1 ) can no longer be neglected.
where X = (k+ k-)pBe, pBe = k + p / ( K + k-), with k* the usual mass action law rate constants and p the total number density of A and B. For all 0 < pBo 5 p , pB approaches its equilibrium
(15) Kahaner, D.; Barnett, D. PLOD, version 4.9. National Bureau of Standards: Washington, DC. (16) Gear, C. W. Math. Comp. 1967, 21, 146.
A+Be2B
(12)
For this reaction the canonical thermodynamic form of the forward and reverse rates are2
v,+ =
v Q
exp[(pA + P B ) / ~ B T ~
V.- = VQ exp[2/.bB/kBq
(13)
with p l representing the chemical potentials of A and B, k B Boltzmann's constant, T the absolute temperature, and Q the intrinsic rate of the reaction. Assuming an absence of inputs, Le., S = 0 in (2), and assuming that the reaction occurs in ideal solution, the conditional average number density of B for this reaction is easily found to be2 RB/ V
E
p~(p~",= t)p
+
~1 -~ ( 1[ - p ~ ' / p ~ "exp(-Xt)]-' )
+
The Journal of Physical Chemistry, Vol. 93, No. 7 , 1989 2813
Molecular Fluctuations
A more useful characterization of this connection between orbital asymptotic stability and fluctuations is in terms of the Liapunov characteristic exponentst3
0,
3 x
s = 5
5oi
1
v i
z=
Pk
0200
0400
exp [pk(t-t 71Z(kj)(t-t’;t’) 0600
0800
1000
Figure 2. Single-time probability density for the reaction A + B s 2B, as described in the legend of Figure 1 , for the four scaled times 7 = At. The initial probability density is a Gaussian centered at pe/pee = 0.1 with a variance of 9 x lo4.
111. Molecular Fluctuations on Limit Cycles
Fluctuations for models whose deterministic trajectories contain limit cycles have been emplored by using both stochastic differential e q ~ a t i o n s ’ ~and - ’ ~ the master equation formalism.2*2z It is known that the conditional covariance of fluctuations around an orbitally asymptotically stable limit cycle initially tends to increase linearly with time and that, asymptotically, in two dimensions the probability density is shaped like a cratert2 whose rim coincides with the deterministic limit cycle. This implies that after sufficiently long times the systems loses all knowledge of its initial phase on the limit cycle. For a macroscopic system the width of the “crater” is vanishingly small compared to its diameter, so that the system stays very close to the limit cycle as time proceeds but simply loses its phase coherence. Nonetheless, in this long-time limit the fluctuations become comparable in magnitude to the average. Little is known about the time scale on which this asymptotic distribution is approached or how rapidly the fluctuations increase in size. Using the statistical theory of nonequilibrium thermodynamics described in the previous section, we have found a simple relationship between the asymptotic orbital stability of a limit cycle and the asymptotic growth of the fluctuations. To make things definite, we consider a system whose conditional average satisfies an autonomous equation like (2) which possess a orbitally asymptotically stable limit cycle,z1nlc(t),with period, T . Thus nl, solves dn,c/dt = R ( n 3
(18)
and nlc(n0,t+T) = nlc(no,t)whenever no is on the limit cycle. The limit cycle is one of an uncountable infinity of possible conditional average trajectories, ii, all of which solve (18). A necessary and sufficient condition that a limit cycle be orbitally asymptotically stable is that its variational equation have all of its Floquet multipliers, p,, save one, of magnitude less than l.23 The connection of this result with fluctuations is that the variational equation is dAn/dt = H(t)An
(19)
(20)
It is well-known that a fundamental solution to (19) is given by the solution matrix to (7). Moreover, this fundamental solution can be expressed as a linear combination of matrix-valued functions of the formz3
T = 125
0000
(1/T) In bk
(21)
where k labels the distinct eigenvalues of the matrix K(0,T); 1 Ij Iflk, with f l k the multiplicity of the elementary divisor of the kth eigenvalue and Z(kn(t-t’;t? a polynomial in its first argument with periodic coefficients of degree at most nk - 1 and periodic in both t and t’. This result permits a simple characterization of the behavior of conditional fluctuations on an orbitally asymptotically stable limit cycle. To see this, we recall that the statistical nonequilibrium thermodynamics theory in section I1 shows that conditional fluctuations are Gaussian with an average satisfying (1 8). For a system that starts from an initial point on the limit cycle, this average will move along the limit cycle, nlc(no,t). Since we are considering only autonomous equations, we can select the initial time t o = 0, so that the covariance matrix of the conditional fluctuations around the limit cycle is given by (7), i.e. a(t) = l f0K ( t ’ , t )y(t’) f l ( t ’ , t ) dt’
(22)
dK(t’,t)/dt = H ( t ) K(t’,t)
(23)
with and K(t’,t’) = 1. Since the covariance matrix of the random terms, y(t’), inherits2 its time dependence from n,,(t), it is a periodic function of period T . Thus the integrand of (22) consists of a sum of terms of the form exp[(pk + p y ) ( t - t’)]A(k,y)(t-t’;t’)c (pk(for k > 0. Thus there is a unique Liapunov exponent that vanishes, namely po, and the real parts of all the other Liapunov exponents are negative. According to (24), this implies that for t large enough all terms in the integrand (22), except the (0,O)term for which (po po) 1 0, will lead to convergent integrals. The (0,O) term, on the other hand, can be written
+
xfA(osO)(t-t’;t’)c
