Entropy and Chaos in a Lattice Gas Cellular Automata

Entropy and Chaos in a Lattice Gas Cellular Automata Franco Bagnoli1 and Ra´ ul Rechtman2 1

Dipartimento di Energetica, Universit` a di Firenze, Via S. Marta, 3 I-50139 Firenze, Italy Also CSDC and INFN, sez. Firenze [email protected] 2 Centro de Investigac´ıon en Energ´ıa, UNAM, 62580 Temixco, Morelos, Mexico [email protected]

Abstract. We find a simple linear relation between the thermodynamic entropy and the largest Lyapunov exponent (LLE) of an discrete hydrodynamical system, a deterministic, two-dimensional lattice gas automaton (LGCA). This relation can be extended to irreversible processes considering the Boltzmann’s H function and the expansion factor of the LLE. The definition of LLE for cellular automata is based on the concept of Boolean derivatives and is formally equivalent to that of continuous dynamical systems.

1

Introduction

The relation between thermodynamics and the underlying chaotic properties of a system is of great relevance in the foundations of statistical mechanics [1] and has attracted much interest. For a family of models of simple liquids, including Lennard–Jones, a simple relation between the Kolmogorov–Sinai entropy and the thermodynamic entropy exists [2]. The connection between chaotic dynamical properties and transport coefficients has been established for Lorentz gases [3]. Investigations on extended or multi-body systems are difficult, and most of previous studies were restricted to non-interacting particles on some disordered configuration of scatteres. Moreover, continuous or partially continuous systems require coarse-graining that introduce another degree of arbitrariety in the analysis. These difficulties have lead to the study of rather abstract models (maps) [4]. In this paper we investigate if and to what extent similar relations apply to completely discrete systems like cellular automata, and in particular to lattice gas cellular automata (LGCA). LGCA are simple models with hydrodynamical behavior [6]. The D2Q9 LGCA is a two dimensional model with nine velocities and is one of the simplest models where temperature can be defined [7]. In this paper we find that the thermodynamic entropy density is proportional to the largest Lyapunov exponent (LLE) of a deterministic D2Q9 LGCA. Furthermore, in a simple irreversible process, H. Umeo et al. (Eds): ACRI 2008, LNCS 5191, pp. 120–127, 2008. c Springer-Verlag Berlin Heidelberg 2008 

Entropy and Chaos in a Lattice Gas Cellular Automata

121

A C I

I

C B

Fig. 1. (Top) All these states have the same number of particles, momentum and energy. An arrow represents the presence of a particle with the velocity of the arrow, the open circle a particle at rest. (bottom) Example of entries in a reversible collision table.

there is a linear relation between Boltzmann’s H function and the expansion factor of the LLE. The sensitivity to initial conditions for discrete dynamical systems like cellular automata are has to be understood with care. One cannot perform here the usual limit of vanishing initial distance between configurations. Nevertheless, Lyapunov exponents can be defined and cellular automata which exhibit complex space time patterns have a positive LLE [8]. This definition is closely related to the one used in continuous maps [9].

2

Model

The D2Q9 model is defined on a two dimensional square lattice. The evolution is in discrete time steps and unit mass particles at every site r can move with one of nine velocities c0 = (0, 0), c1 = (1, 0), c2 = (0, 1), c3 = (−1, 0), c4 = (0, −1), c5 = (1, 1), c6 = (−1, 1), c7 = (−1, −1), c8 = (1, −1). The state of the automaton is given by the set of occupation numbers s(t) = {sk (r, t)}, where sk (r, t) = 1 (0) indicates the presence (absence) of a particle with velocity ck at site r and time t. An exclusion principle forbids the presence of more than one particle in a given site, at a given time with a given velocity. The time evolution of the system is given by collision and streaming operations. In the collison operation, particles at a given site collide conserving mass and momentum. In the streaming operation particles move to neighboring sites according to their velocities. Since the number of states for a given site is finite (29 ), the local collision operator C is generally implemented as a look-up table. Given a local configuration, the conservation constraints may not define the outgoing state completely as we show in Fig. 1 (top). One may assume that the C look-up table has several columns for all the possible output states. In order to make the automaton deterministic, we assign at the beginning an integer random number η(r) to each

122

F. Bagnoli and R. Rechtman

site, to be used in choosing the column from which the output state is chosen (quenched disorder). That is sk (r, t + 1) = Ck (s0 (r, t), . . . , s8 (r, t), η(r)).

(1)

The choice of the quenched disorder is analogous to the random disposition of scatterers in a wind-tree or other similar Lorentz gases [10]. We have checked in the various conditions studied in this paper that the system is self-averaging: a simulation on a large enough system gives the “average” value over disorder. This is illustrated for instance in Figure 2 (left): the plot is a collection of singlerun simulations for different values of a parameter, using different realizations of the quenched disorder. One can notice that fluctuations are quite small. Moreover, the collisions can be made reversible, an interesting feature for the discussion of irreversibility. To do so, states s must satisfy the condition that in every column s = ICIC(s) holds where I denotes the operator that inverts the velocities. This means that a collision C followed by an inversion I and another collision and inversion leaves the state unchanged when it is taken from the same column in the collision table as we show in Fig. 1 (bottom).

3

Lyapunov Exponents

We recall here the salient points of the definition of largest Lyapunov exponent for cellular automata [8]. Let us simplify the notation by using the index n to indicate both the position r and the velocity k, with n = 1, . . . , 9L, where L is the number of sites of the automata. The evolution of the state sn (t) of the LGCA can be seen as an application of a set of Boolean functions sn (t + 1) = Fn (s(t)).

(2)

The functions Fn and differ in the velocity index k and quenched disorder η(r). However, since the distribution of the disorder is uniform, and, as shown in the following, the correlations among variables decay very fast, the system is translationally invariant at a mesoscopic level. Let s(0) and x(0) be two initially close configurations, for example all the components of x(0) may be equal to those of s(0) except for one. We define the bitwise difference between these two configurations with the term “damage”. The smallest possible damage is one and the damage vector v(0) is one in the component where s(0) and x(0) are different and zero in all the others. If this damage grows in average during time, the trajectory is unstable with respect to the smallest perturbation. However, due to the discrete nature of LGCA, defects may annihilate during time evolution, altering the measure of instability of trajectories. The correct way of testing for instability is that of considering all possible ways of inserting the smallest damage in a configuration, using as many replicas as the number of components of the configuration, and letting them evolve for one time step counting if the number of damages has grown or diminished. The ensemble of all possible replicas with one damage each is the equivalent of the tangent space for discrete systems.

Entropy and Chaos in a Lattice Gas Cellular Automata

123

The task of computing the evolution in tangent space is clearly daunting, but by exploiting the concept of Boolean derivatives [8,11], it is possible to develop a formula very similar to the one used in continuous systems. The Boolean derivative is defined by ∂Fn (s) = Fn (. . . , sp , . . . ) ⊕ Fn (. . . , 1 − sp . . . ), ∂sp

(3)

with n, p = 1, . . . , 9L and where ⊕ represents the sum modulo two (XOR). This quantity measures the sensitivity of the function Fn with respect to a change in sp . The Jacobian matrix J(s) has components Jnp = ∂Fn (s)/∂sp . We now consider the map v(t + 1) = J(s(t))v(t),

(4)

 with v(0) as mentioned above. It is easy to check that |v(t)| = i vi (t) is the number of different paths along which a damage may grow in tangent space during time evolution, i.e., with the prescription of just one defect per replica [8]. If there is sensitivity with respect to initial conditions, one expects that |v(T )|/|v(0)| ∼ exp(λT T ) for large T where λT is the largest finite time Lyapunov exponent. It then follows that λT =

T −1 1  log u(t) = log uT , T t=1

(5)

where u(t) = |v(t)|/|v(t − 1)| is the expansion factor of the LLE. The definition should include the limit when T → ∞ but in practice we always evaluate the finite time LLE. The LLE depends in principle on the initial configuration s(0), initial damage v(0) and quenched disorder η, but in practice it assumes the same value for all trajectories corresponding to the same macroscopic observables when T is sufficiently large. The LLE of CA as defined above has been used to classify elementary and totalistic Boolean cellular automata [8,12].

4

Entropy

We define the entropy S as S=−



p(s) log(p(s)),

(6)

s

where s denotes a possible configuration of the system and p(s) is the corresponding probability, computed averaging over a set of replicas (statistical ensemble), analogous to the configuration used for computing the LLE. In equilibrium, S is just the the logarithm of the number of possible configurations satisfying the constraints of constant mass and energy. It is possible to extend the ensemble definition of entropy also to non-equilibrium conditions, S = S(t).

124

F. Bagnoli and R. Rechtman

The single particle velocity distribution functions fk (r, t) are defined as the average number of particles at site r with velocity ck at time t over R samples that share the same quenched disorder η and macroscopic constraints with different microscopic initial configurations. Since the model satisfies Fermi statistics and correlations are unimportant we approximate the probability distribution by   1−sk (r) fk (r)sk (r) 1 − fk (r) , (7) p(s)  r,k

and the entropy can be approximated by the Boltzmann function H, defined as        fk (r) log fk (r) + 1 − fk (r) log 1 − fk (r) . (8) H=− r,k

Defining the connected correlation function Ckl (η) as  2 1 1 sk (r, t)sl (r + η, t) − sk (r, t) , Ckl (η) = L r L r

(9)

we checked that the single-site two-particle correlation function Ckl (0) factorizes into the product of single particle distributions before and after the collision, both during the relaxation (out of equilibrium) phase and in the stationary (equilibrium) conditions. In equilibrium, the correlation function Ckl (η) decays to zero for η = 1. Out of equilibrium, starting with very different configurations in the two halves of the system, the correlation function, still being quite small, exhibits a correlation length of some lattice spacings for a short time. This corresponds to the coherent motion of particles in a shock wave, where the local density is near to zero or nine. However, this correlation quickly disappears; although the motion is correlated at a macroscopic level, as soon as the local density of particles is different from the extremes (for which the collision table has few output configurations) the velocities quickly decorrelate. In equilibrium, the distribution functions do not depend on r or t and we may write 8 8   fk , e= k fk , (10) n= k=0

k=0

where n and e are the number and energy densities respectively, 0 = 0, 1,2,3,4 = 1/2, 5,6,7,8 = 1, 0 ≤ n ≤ 9, and 0 ≤ e ≤ 6, Also, f1 = f2 = f3 = f4 and f5 = f6 = f7 = f8 . There is a lower and an upper bound on e: the Fermi energy per site eF and the maximum energy per site eM . For example, if n = 4, eF = 3/2 that corresponds to a configuration where at every site there is a particle at rest and the other three have velocities along the axes. On the other hand, eM = 4 corresponds to a state where every site is occupied by four particles with velocities along the diagonals. The existence of these bounds has a deep impact on the entropy of the model: it cannot be an always increasing function of e since the system is

Entropy and Chaos in a Lattice Gas Cellular Automata 0.7

1

0.7

125 1

0.9

0.6

0.95

0.65

0.8 0.5

0.7

0.9

0.6

1

0.6

S

0.5 0.3

S λ

λ

0.4

0.95

H

0.85

0.3



0.8

H

0.4

0.2

0.85

0.9

0.55

0.8

0.5

0.75

0.45

0.65

0.75

0.7

0.2 0.1

0.4

0.45

0.5

0.55



0.1

0

0 1.5

2

2.5

3

3.5

4

0.6

0.65

0.4 0

500

e

1000

1500

0.7

0.7 0.65 2000

t

Fig. 2. (Left) Entropy per site S (solid line) and LLE λ (dashed line) for n = 4, simulations with R = 40 in a 40×40 lattice. The simulations are performed for different values of e from eF (n) to eM (n). The LLE λ is roughly proportional to S. (right) Boltzmann function H (thick line) and Lyapunov expansion factor  (oscillating thin line) versus time t, simulations with nL = 7.2, eL = 4.8, nR = 1.8, eR = 1.2, and R = 40 on a 40 × 40 lattice. The inset shows that, disregarding the oscillations of , there is a linear relation between these quantities. The dashed line is the best fit H = 1.15 + 0.2.

ordered both for small and large e. The three unknown distribution functions f0 , f1 and f5 satisfy the conservation of mass and energy expressed in Eqs. (10) and maximize the entropy which gives an additional constraint [13] f12 (1 − f0 − f5 ) = f0 f5 (1 − 2f1 ).

5

(11)

Results and Discussion

In Fig. 2 (left) we show the entropy density s = S/L and the LLE λ as a function of e for n = 4. The entropy density grows for eF ≤ e ≤ e∗ and then decreases for e∗ ≤ e ≤ eM with e∗ = 2n/3. The largest Lyapunov exponent λ shows the same behavior and is rougly proportional to the entropy s. There are small systematic deviations as shown in the inset of the Fig. 2 (left), where λ is plotted as a funciton of s for several values of n. This proportionality can understood with a simple stochastic model, in which the mesoscopic dynamics (coarse-grained in space or averaged over quenched disorder) can be approximated by a Markov chain. We assume that for given macroscopic observables, there are M compatible configurations, and that the probability of observing a configuration x is P (x) = 1/M (microcanonical distribution). The entropy is therefore S = log(M ). The (chaotic) dynamics is approximated by a Markov matrix W (x |x) with αM equal entries per row and column with value 1/αM , and we assume that W is irreducible (the asymptotic state is unique). The Kolmogorov-Sinai entropy K per unit time is   K=− P (x) W (x |x) log W (x |x) = log(αM ). x

x

126

F. Bagnoli and R. Rechtman

In the thermodynamic limit, the dominating term is α log M . On the other hand, K is roughly proportional to the largest exponent λ and if we assume that the form of the Lyapunov spectrum does not change dramatically when the macroscopic variables change, K  S + log(α). The value of the LLE is related to the number of ones in the Jacobian matrix J defined by Eq. (3). This Jacobian matrix contains the linearized effects of the streaming and collision operators. Streaming corresponds to a scrambling of the components of the tangent vector v(t) and therefore does not alter its norm. This is left to collisions, when more than one output configurations are possible. The number of equivalent output configurations in the collision table is small for the extreme values of number and energy densities, and larger for intermediate values. Similar considerations apply to the number of equivalent configurations for a given macroscopic distribution of density and velocities, and constitute the microscopic origin of the proportionality between statistical and dynamical quantities. We now discuss an irreversible process where a square lattice is initially in an equilibrium state with the left and right sides having different number and energy densities nL , nR , el , and eR . The system evolves toward equilibrium by means of dumped travelling waves. Boltzmann’s H function is defined in the same way as the entropy, Eq. (8), with a minus sign and distribution functions fk (r, t) that are found as averages over R replicas. In the same R numerical simulation, the average Lyapunov expansion factor  = (1/R) i=1 log u(i) is calculated. As shown in Fig. 2 (right), the two quantities exhibit similar behavior. The Lyapunov expansion factor exhibits more marked oscillations, indicating that this quantity is more sensible to the local variations in density. The inset of Fig. 2 (right) shows that, disregarding oscillations,  is linearly related to H for all the relaxation phase. The computation of  is performed using a set of tangent vectors, Eq. (4), and these vectors constitute a sort of local memory of the past state. In systems with local variations of density, as in our system in the presence of travelling waves, statistical quantities like H depend on the instantaneous state of the system, while dynamical ones like  depend also on the variations of this state. This factor may be the origin of the different relation between statistical and dynamical quantities in equilibrium and during the relaxation phase. The D2Q9 reversible LGCA model we have discussed exhibits hydrodynamical and thermodynamical behavior, and is therefore one of the more “realistic” models for which a simple relation can be established between its chaotic dynamical properties and its macroscopic behavior in irreversible processes and in equilibrium.

Acknowledgments The authors thank Stefano Ruffo for helpful discussions. Partial economic support from CONACyT projects U41347 and 2516, and from the Coordinaci´ on de la Investigaci´on Cient´ıfica de la UNAM is acknowledged.

Entropy and Chaos in a Lattice Gas Cellular Automata

127

References 1. Lebowitz, J.: Phys.Today, p. 32 (September 1973); Zaslavsky, G. M.: Phys. Today, p. 39 (August 2004) 2. Dzugutov, M., Aurell, E., Vulpiani, A.: Phys. Rev. Lett. 81, 1762 (1998) 3. Gaspard, P., Nicolis, G.: Phys. Rev. Lett. 65, 1693 (1990); Gaspard, P., Dorfman, J. R.: Phys. Rev. E 52, 3525 (1995); Evans, D.J., Morris, G.P.: Statistical Mechanics of Nonequilibrium Liquids. Academic Press, London (1990); Hoover, W.G.: Computational Statistical Mechanics. Elsevier, Amsterdam (1991); Dorfman, J.R., van Beijeren, H.: Physica A 240, 12 (1997) 4. Falcioni, M., Palatella, L., Vulpiani, A.: Phys. Rev. E 71, 016118 (2005); Latora, V., Baranger, M.: Phys. Rev. Lett. 82, 520 (1999) 5. d’Humi´eres, D., Lallemand, P., Frisch, U.: Europhys. Lett. 2, 291 (1986); Chopard, B., Droz, M.: Phys. Lett. 126A, 476 (1988); Chen, S., Lee, M., Zhao, K., Doolen, G.D.: Physica B 37, 42 (1989) 6. Hardy, J., Pomeau, Y., de Pazzis, O.: J. Math. Phys. 14, 1746, (1973); Frisch, U., Hasslacher, B., Pomeau, Y.: Phys. Rev. Lett. 56, 1505 (1986); Frisch, U., d’Humi´eres, D., Hasslacher, B., Lallemand, P., Pomeau, Y., Rivet, J.P.: Complex Systems 1, 649 (1987). Reprinted in Lecture Notes on Turbulence, Herring, J.R., McWilliams, J.C.(eds.): p. 649. World Scientific, Singapore (1989) . 7. Bagnoli, F., Rechtman, R., Zanette, D.: Rev. Mex. F´ıs. 39, 763 (1993) 8. Bagnoli, F., Rechtman, R., Ruffo, S.: Phys. Lett. A, 172, 34 (1992) 9. Ott, E.: Chaos in Dynamical Systems. Cambridge University Press, Cambridge (1993) 10. Ehrenfest, P., Ehrenfest, T.: Enzyclop¨ adie der Mathematischen Wissenschaften, vol. 4 (Teubner, Leipzig, 1911) English translation: Moravcik, M. J.: Conceptual Foundations of Statistical Approach in Mechanics. Cornell University Press, Ithaca (1959). Ernst, M. H., van Velzen, G. A.: J. Phys. A: Math. Gen. 22, 4611 (1989) 11. Bagnoli, F.: Int. J. Mod. Phys. 3, 307 (1992) 12. Bagnoli, F., Rechtman, R., Ruffo, S.: Lyapunov Exponents for Cellular Automata. In: L´ opez de Haro, M., Varea, C. (eds.) Lectures on Thermodynamics and Statistical Mechanics, p. 200. World Scientific, Singapore (1994) 13. Salcido, A., Rechtman, R.: Equilibrium Properties of a Cellular Automaton for Thermofluid Dynamics. In: Cordero, P., Nachtergaele, B. (eds.) Nonlinear Phenomena in Fluids, Solids and Other Complex Systems, p. 217. Elsevier, Amsterdam (1991)