Complex dynamics in systems with many degrees of freedom - CaltechTHESIS

We now summarize briefly some of the work that has been done using Eq. (1.15).
Kaneko’s original work consisted primarily of a classification of the behavior as a
function of the coupling strength e¢ and the nonlinear parameter in the local map
function f(x), e.g., ‘a’ if f(x) is the logistic map. He found a variety of different phe-
nomena, ranging from the formation of complex patterns, to the correlated motion
of defects [32].

Several groups [33, 34, 35] have used coupled map lattices to study the phe-
nomenon of spatiotemporal intermittency [15]. An interesting conjecture by Pomeau
[36], proposing to describe intermittent transitions in terms of directed percolation,
has also led to various investigations of nonequilibrium phase transitions using cou-
pled map lattices [37].

Bunimovich and Sinai [38] have shown that for weak coupling and certain condi-
tions on the map function f(z), a mapping to a 2D statistical mechanics spin model
is possible, where the discrete time index plays the role of the second dimension.

Several other interesting and ambitious questions have been addressed using

12

coupled map models, with results that still need to be considered as rather specu-
lative. For example, Crutchfield and Kaneko [39] have constructed a coupled map
model similar to Eq. (1.15), which displayed chaotic motion over a transient time
which increased exponentially with the system size. The authors then argue that
many complex phenomena in nature, including turbulence, may really be transients
where the asymptotic state is never observed. Another interesting suggestion, based
on numerical investigations of coupled map models, is that chaos is a purely local
phenomenon [40, 41], i.e., spatially coarse-grained variables will not be chaotic.
Finally, there is an interesting neural network model consisting of infinitely many
sites [42], where every lattice site couples to all the other sites. The dynamics in
this case can be reduced to a self-consistent equation for a single site, allowing an
exact determination of the largest Lyapunov exponent. The model shows a sharp
transition from stationary to chaotic behavior as the nonlinear map parameter is
tuned, and represents one of the few model systems with many degrees of freedom

where exact results with chaotic behavior have been obtained.

1.2.3. Conservation Laws

We constructed our own coupled map model based on Eg. (1.15), but we also in-
corporate a conservation law [43]. There are a number of important reasons for
imposing conservation laws.

Firstly, conservation laws occur naturally in many situations we might try to de-
scribe. For example, in the Rayleigh-Bénard convection experiment described in the
introduction, the average fluid density is conserved. Similarly, in the surface wave
experiments we described, the average level of fluid in the container is conserved.
Indeed, in the derivation of the Navier-Stokes equation conservation laws play the
central role, since mass, momentum, and energy conservation persist through coarse
graining procedures.

Secondly, apparent conservation laws can arise through broken symmetries. An

13

example is the superfluid velocity in a Bose condensed system v, « V@ with ¢
the phase of the macroscopic wavefunction: the dynamics of v, is governed by a
conservation law, deriving in fact from the equation for ¢ that describes the broken

gauge symmetry:

so that nV, = VF{vs,..-}, (1.16)

which has the form of a conservation law. Here F cannot depend on ¢ itself by
the basic phase symmetry. An analogous case in nonequilibrium systems is the

Kuramoto-Sivashinsky equation [12]:
O,u(z,t) = —O2u — O4u — udzu. (1.17)

This equation may be derived as the dynamics of a phase variable ¢@ which here
characterizes displacements of a periodic pattern (broken translational symmetry)
with u = 0,¢.

Finally, it is known that in equilibrium systems interesting long length scale
dynamics usually derives from conservation laws: in the absence of these, long
wavelength properties often reduce to simple averages over independent, rapidly
fluctuating microscopic regions, usually yielding uninteresting results. We expect
the same to be true for dissipative systems, although there may still be interesting
questions on how to characterize the short length scale behavior in systems without
conserved quantities.

With these considerations, we modify Eq. (1.15) to the following model which

satisfies a discrete version of the continuity equation:

ye = YP + outa) — oUF YR4)- (1.18)

14

The role of the current term in the continuity equation is played by the scalar
function g. For simplicity, we impose g(z, y) = —g(y, x) to make the model invariant

under the parity transformation i ~ —i. We enforce periodic boundary conditions:

Vien = Yee (1.19)

The conserved ‘mass’ Q is evidently given by:

Q=d iy. (1.20)

We consider two different choices for g, which we will refer to as ‘model A’ and

‘model B’. Model A is defined by:

where f(x) = ar+bz(1—-z)
|x] (mod 1),

with z

where a and 6 are constants. The functional form of f is shown in Fig. 1-1.

Model B is defined by:

2f(x — y)
1+27+y?,
with f(z) = ax + bsin(cz),

g(x,y) (1.22)

where a, 6, and c are constants.

The change from Eq. (1.15) to Eq. (1.18) introduces some slight complications for
the choice of map function. For example, in the non-conserved model Eq. (1.15) it is
easy to ensure that the lattice amplitudes are always bounded: if f(x) and € range
between 0 and 1, the lattice amplitudes are also bounded between 0 and 1. In the

conserved model Eq. (1.18) this is no longer the case, and the map function has to be

15

1 T T T T T T T T T T T T 7 T T T

0.5

f(x)
ro)

—0.5

T ] T T T T | T T T T | T T T T

L i i £ | i L i 1 | 1 i it it I i 1

_ I 4 i 1 | 1 4 1 i | 1 i 13 i | i i L 3

—2 -1 0 1

Figure 1-1: Map function for model A.

chosen appropriately not only to ensure the boundedness of the lattice amplitudes,
but also to avoid the approach of a trivial fixed point or periodic orbit. Model A
and model B will remain bounded if the diffusion constant a satisfies 0 < a < 1/2.
Avoiding trivial periodic orbits is not as easily arranged. We have found empirically
that the most likely periodic orbit is a ‘zig-zag’ state where adjacent lattice sites
have amplitudes z; and z, which are exchanged every iteration, leading to a state
of period two. The linear stability of this state can easily be determined from the
equations of motion, and the map parameters can then be chosen to make this orbit
unstable. For more complicated periodic orbits a direct stability analysis is not
practical, and we rely on trial and error and the methods described in the next

section to ensure chaotic motion.

16

1.3 Characterization of Spatially Extended Chaos

1.3.1 Lyapunov Exponents

We now describe simple numerical experiments that have shed some light on the
nature of the extended chaotic state and which motivated the search for a simple
long wavelength description.

For b = 0, model A (cf. Eq. (1.21)) evidently reduces to a finite difference
approximation to the one-dimensional linear diffusion equation. Almost all initial
configurations decay to a spatially uniform state y? = const. As we increase b, we
observe oscillations of period two, and at larger b oscillations of period four. At a
critical value b, we measure a positive Lyapunov exponent, signaling the onset of
chaos. The onset of chaos is in the form of ‘local defect turbulence’, i.e., most of the
lattice sites are still oscillating with period four and only a small number of sites
shows chaotic behavior. In Figure 1-2 we show how the largest Lyapunov exponent
Amaz depends on the nonlinear control parameter b. We find that just above the
onset of chaos, Amar ~ (b — b-)'/?.

It is tempting to regard \j,az as a kind of order parameter (or rather disorder
parameter) [44], which is zero below the onset of chaos. Indeed, for the logistic
map Eq. (1.2) it has been shown that the envelope of the curve A(a) has the scaling
form A(a) ~ (a — a,)’, with a universal constant 7 = 0.4498... [45]. We should
emphasize, however, that the scaling relations seen in the logistic map are a result
of the period-doubling cascade. This is quite different from the scaling relations
observed in second order phase transitions, which occur only in the thermodynamic
limit of infinite volume.

We are mainly interested in the situation well above the onset of chaos, where
a whole set of positive Lyapunov exponents appears. We therefore not only require
the largest exponent, but the full Lyapunov spectrum. Considerable care is required

to obtain accurate numerical results and we detail our computational procedure in

17

0.3 T T T T T 7 T t | T T T T T T T T

Ir Model A 7

0.2 [-

A max

0.1

Figure 1-2: Maximum Lyapunov exponents.

the Appendix.

In Figure 1-3 we show the Lyapunov spectra for two different lattice sizes. These
spectra show several interesting features.

One obvious feature is that the sum of the Lyapunov exponents is greater than
zero. For a set of N autonomous ordinary first order differential equations of the

form:
ad .
<= fi(x1,-+-5%w), i=1,2,...,N (1.23)

it is easily shown that the sum of the Lyapunov exponents gives the average expan-

sion of a volume element:

Dri = (V-f), (1.24)
i=1
where (...) denotes time averaging. This in particular leads to the usual definition of

a dissipative system as one in which the phase space volume contracts on average,

i.e., the sum in Eq. (1.24) is negative. For discrete maps there does not exist a

18

1 TOT PT ] TT FT | TOT TT | TOT tT | TT T FT
|. 4
Pes. Model A |
0.5 oe Hay. “|
[ he 4
a, - :
—0.5 a=0.4, b=1.3 to
[ +++ L=50 1
be 4
~1- - L=100 _
—-1.5 L 1 1 i 1 | J it i 2 | L 1 i 1 | on one | 1 | 1 1 L i |
0 0.2 0.4 0.6 0.8 1

i/L
Figure 1-3: Lyapunov spectra.

continuous trajectory in phase space, and Eq. (1.24) does not apply. Nevertheless,
it remains an interesting question how to characterize the effective dimension of
the chaotic attractor. In general, fractal objects must be characterized by a whole
family of generalized dimensions [46]. We mention only the information dimension
Dy, since a relation between D; and the Lyapunov exponents has been conjectured
[47]. To obtain D;, we divide phase space into N small boxes with linear dimension

€, and determine the probability p; that the dynamical system is found in box 2.

Then
Ly Pilog(p:) _,. Se)
= ] . = _~
Dy lim Neves log(e) lim log(e)’

(1.25)

where S(¢) can be interpreted as the ‘information content’ of the partition [48]. The
Lyapunov dimension D, is defined as:
J r;

Dr, = J + isl 5)
lAj+1l

(1.26)

19

where j is the largest integer such that:

3 A; > 0. (1.27)

i=1
It can be rigorously shown that D; < Dy [49], and it was suggested that in fact
an equality holds [47]. Counterexamples to this conjecture are known [50, 26], and
we must evidently consider our coupled map model as yet another one. This is
somewhat unfortunate, since numerical dimension estimates are only feasible for
low-dimensional systems (d ~ 5), because the numerical complexity increases expo-
nentially with d [51]. As described in the Appendix, the numerical complexity of
Lyapunov spectra is ~ d?, allowing calculations with much larger d. However, it
is not yet clear how useful the knowledge of one or more generalized dimensions is
when the number of degrees of freedom is as large as in the coupled map models,
and we do not address this issue further in this work.

Another interesting feature of the Lyapunov spectrum (cf. Figure 1-3) is an
approximate point of inflection near 4 = 0. Some mathematical arguments for
singularities in the density of Lyapunov exponents exist for certain models of in-
termittent turbulence [52], but there is no clear connection of those arguments to
the situation here. Obviously, one exponent is always trivially zero because of the
conservation law (i.e., a vector for which all components are equal does not change
under the dynamics and is thus always an eigenvector with eigenvalue one), but
perhaps a connection can be established between small but non-zero exponents and
the slow modes of the system. A point of inflection also appears in the spectrum of
model B, indicating that it might be quite generic. An obvious possibility is that
the negative exponents are simply due to the diffusive modes. This is most easily
seen by solving the linearized equations, for example model A with 6 = 0, for which

we can find the Lyapunov exponents analytically:

Ap = log(|1 — 4asin?(k)]), (1.28)

20

where k ranges from 0 to the zone boundary k = 0.5 in steps of the reciprocal lattice

vector 1/L. For k < 1 we then obtain:
Ap & —4an?k?, (1.29)

clearly showing parabolic behavior.

In the strongly nonlinear chaotic regime, the notion of Fourier analysis is at first
sight not useful, because it is essentially based on the principle of linear superpo-
sition. For example, as pointed out by Eckmann and Ruelle [26], the experimental
observation of a continuous power spectrum can either indicate that the system ex-
plores an infinite number of dimensions in phase space, or that the system evolves
nonlinearly on a finite-dimensional attractor. Nevertheless, it may happen that
some of the degrees of freedom can be separated out and subjected to a traditional
Fourier analysis. In our coupled map models, a trivial example is the Fourier mode
with zero wavevector, which is not affected by the nonlinear dynamics due to the
conservation law. We can attempt a similar analysis by numerically measuring the
temporal decay of the spatial Fourier transform of the space-time correlation func-
tion for finite wavevectors, as discussed in section 1.4 below. We indeed find that
Fourier modes with wave vector k decay exponentially in time, with time constants
T, ~ k®. However, the magnitudes of the 7, do not agree with the inverse of the
Lyapunov exponents of Figure 1-3, suggesting that further work using a different
approach is necessary.

Perhaps the most interesting feature of Figure 1-3 is that the exponents for two
different system sizes lie on top of one another, if their coordinates are rescaled by
the system size. This implies that the fraction of exponents which lie in a given range
does not change with the system size. For example, the mean rate of information
production, called the Kolmogorov-Sinai entropy [26], is believed to be given by
the sum of the positive Lyapunov exponents. In Figure 1-3 this quantity is then

clearly additive. The existence of additive quantities hints at the possibility of a

21

thermodynamic formalism; we will explore this question further in section 1.6 below.
Another implication of additive quantities is simply that the system size is much
larger than a certain correlation length, such that averaging over the system amounts
to averaging over many uncorrelated regions. We can then apply the central limit
theorem to these averages and predict Gaussian distributions, as discussed in the

next section.

1.3.2 Correlation Lengths

To investigate whether our coupled map models indeed have short range correlations,

we can measure the equal time correlation function S(r) [3, 53]:

S(r) = ((u? - (u) (urs. — (Y)))» (1.30)

where (...) denotes averaging over both 7 and n. Before presenting our results, we
list a few caveats that need to be observed in interpreting correlation functions.

In many cases corresponding to short range correlations, we expect S(r) ~
exp(—r/€) for large r, defining the correlation length €. However, the actual form of
S(r) may be considerably more complicated. For example, unless one length scale
completely dominates all the others (as happens near second order phase transi-
tions), there is no reason to even expect S(r) to take on the form f(r/) with a
uniquely defined €. In addition, S(r) may display complex oscillations about zero as
a function of r, not simply describable as ‘antiferromagnetic’. We therefore believe
that a reasonable definition of ‘short range correlations’ needs to be quite flexible,
perhaps in terms of the distance over which correlations can be distinguished from
experimental noise. In this context, we also emphasize that algebraic decays in cor-
relation functions are not necessarily an indication of long range correlations, as is
often supposed. The power laws observed near second order phase transitions are

the result of, and not the cause of, a diverging correlation length. For example, a

22

10° Fe T TTT TR TT I T FTT ] TT T 7 | TTT T TT T TH

L Model A J

107! ° 3

F a=0.4, b=1.28 ;

_ L L=512 4

S 107° E- o 8 i) er

a E :

107° ° ° ; =
1074 i bout i { i 1 i i | oe oe | i iJ ft | i Lt L | | aoe ee ee

0 2 4 6 8 10 12

Figure 1-4: Spatial correlation function |S(r)].

correlation function which decays as (r/€)~%, where € is a microscopic length, can
very well reflect short range correlations. One might object that in this case higher
order moments will diverge, but it must be remembered that correlation functions
are statistically averaged quantities, and that higher order moments may need to be
averaged over a different distribution, because experimental noise sources are not
necessarily independent.

We now return to our coupled map models, and show |S(r)| for model A in
Figure 1-4. While |S(r)| does not have a simple exponential form, it decays rapidly,
indicating a finite correlation length of approximately two lattice spacings.

Appealing to the central limit theorem, we now expect that fluctuations of spa-
tially averaged quantities will have a Gaussian distribution, since the average will be
over many uncorrelated regions. This will be true even if the decay of correlations

is algebraic rather than exponential, provided the relevant power is large enough.

23

SPT TTY TTT ttt
e k = 10/1024 a=0.4 _
a 4 Gaussian Fit b=1.3 —
< b— —
gob =
ay,

Sok _

x Bb —

a) — ~|
oLtt tiprlii dy LI
-0.4 -0.2 0 0.2 0.4

y;,/k (Real Part)
SET TTT TTT ty TTT Tt TT
E- Model A —-
1.6 |— a=0.4 7
L_ b=1.3 —

$b a

a -- —
0.5 —_
og LLL bitrltlitit Lil
-1 -0.5 0 0.5 1

Yi

Figure 1-5: Distribution of (a) Fourier amplitudes, (b) Lattice amplitudes.

24

As a simple example, we consider the spatial Fourier transform of the lattice

amplitudes:

LS onise
Ye = se Deny? (1.31)
VL 7

where k is a multiple of the reciprocal lattice vector 1/L. From the central limit
theorem, we now expect the fluctuations of yg to have a Gaussian distribution, even
for short wavevectors. The only requirement is of course that the system size L is
large enough. Since yf is complex, we have computed both the real and imaginary
parts of the fluctuations of the Fourier amplitudes Ay? = y? — (yg), averaged over
n. In Figure 1-5a, we show the distribution D(Ay,) of the real part of y?, together
with a Gaussian fit, which evidently describes the distribution well. The imaginary
part of Ay? also has a Gaussian distribution, with the same variance. We have also
computed the distribution of AS? = |Ayf|?, finding an exponential distribution
as expected from the sum of the squares of two independent Gaussian distributed
variables. D(Ay,) should be contrasted with the distribution D(Ay;) of fluctuations
of a single lattice amplitude, shown in Figure 1-5b. The latter distribution is of
course strongly affected by the choice of map function, and cannot be expected to
have a simple form. It is not even symmetric about y = 0, since the equations of
motion Eq. (1.21) do not have that symmetry.

Also note that the width of D(Ay,) scales with wave vector. This is related to
the vanishing of the structure factor S(k) as k — 0 in this model, which we will

discuss further in section 1.4.

1.4 Structure Factor and Langevin Equation

In section 1.3.2 we saw that Fourier modes display Gaussian fluctuations. This
suggests that it may be possible to construct stochastic model equations which
contain simple Gaussian distributed noise terms, and which can reproduce some of

the behavior of the maps. A general ansatz for such stochastic equations is the

25

Langevin equation [55, 56}:
Oy;

a F;(y) + n(t), (1.32)

where y is a random field, and F is in general a nonlinear functional of y. The
label 1 may be discrete or continuous. 7 is a stochastic noise term, which is delta
correlated in time and independent of y. Without loss of generality, we can assume
(n) = 0. If we furthermore assume a Gaussian distribution for 7, all the moments

of 7 are completely specified by the second moment:

(ni(t)nj(t)) = 2Di,6(t — t'), (1.33)

where Dj; is a symmetric positive definite diffusion matrix [56]. We emphasize that
Eq. (1.32) is a stochastic differential equation for the stochastic process y, and the
solution will consist of the moments (y(t)), (y?()), etc., averaged over a suitable
ensemble. Instead of Eq. (1.32) it is also possible to derive a master equation for
the probability distribution P(y,t) [57]. This results in a Fokker-Planck equation
for P from which the moments of y can be computed, giving equivalent results to
the Langevin approach.

To apply these general results to our coupled map models, we need to derive
the form of F in Eq. (1.32). For the moment, we will assume that all nonlinear
terms can be neglected. Since we are interested in the long wavelength behavior, we
assume in addition that we can coarse grain the discrete lattice and time indices to
obtain partial differential equations, rather than finite difference equations. Finally,
we assume that we can expand the spatial dependence in a power series of the
wavevector k. The lowest power of k that can occur in the expansion is k?, because
the constant term is absent due to the conservation law, and the term linear in k
violates the parity symmetry k ~ —k. In addition, we also need to expand the noise
term 7. In principle, the lowest power of k allowed is now k, but as we shall see

more directly below, that term is absent. This can also be deduced from Figure 1-

26

5a, where it is seen that the width of the Gaussian distribution of the fluctuations
scales with the wavevector, introducing an additional factor of k.
We can now write down the resulting Langevin equation:

Oy(k, t)

Sen = -H(Dy(h,t) + 01h 0), (1.34)

where D is a diffusion constant and the noise term 7 has correlations of the form:
(n(k, t)n(k’,t)) = Ad(k + h’)6(t — 2), (1.35)

where A is an effective noise amplitude.
To give a more transparent but less general derivation, we start directly with the

definition of the spatial Fourier transform Eq. (1.31):

yet ae ertiiky an (1.36)

We now substitute for yet? using the equations of motion for model A, Eq. (1.18)

and Eq. (1.21):

1 £ .
yet = VEX ek fy + f(y?) — 2F (yr) + F(yeea)}

4 4sin?(tk)

= Ub — TE ys) (1.37)

where we have made use of the periodic boundary conditions. In the long wavelength
limit, we can approximate yz! — y? by Oy, and sin?(xk) by 1?k2. If we set f(yp) ~
cy? +, we directly arrive at Eq. (1.34). The Langevin equation is therefore basically
the original equation, with all nonlinear terms replaced by random noise, which also

leads to a renormalized diffusion constant. In real space, it is a simple linear diffusion

27

equation:

dy = Dy + Fn, (1.38)

where the noise correlation takes the form:
(n(a,t)n(2!,t!)) = Asa — 2’)é(t—#). (1.39)

We now investigate whether nonlinear terms are important in our Langevin equa-
tion. Naively, one would expect their effect to be small at long wavelengths, since
the noise averages to small values there. In recent years, Langevin models similar
to Eq. (1.38) have been studied intensively because of the possibility that nonlinear
terms can in fact lead to non-integral powers of k in the correlation function and
long time tails in the dynamics through the buildup of the noise at long wavelengths.
A famous example is the random-force driven Burgers equation [58], which has been
studied in many contexts, in particular that of the behavior of stochastically driven
interfaces [59]. In the present case, simple arguments show that no relevant nonlin-
ear terms are expected to appear in Eq. (1.38). As an example, we add to Eq. (1.38)
the nonlinear term 02y? which is the simplest term satisfying both the conserva-
tion law and parity symmetry. To study the effect at long wavelengths, we use the
method of ‘power counting’ [3] as applied in the context of renormalization group
theory.

We perform a scaling transformation + — bz, t — b*t, and y — b*y, where z and
x are scaling exponents to be determined. Including the nonlinear term, Eq. (1.38)

then becomes:
Ory = Db? A2y + ABKt-2A2y? 4 BF/2-X-5/27G 2, (1.40)

where we have used Eq. (1.39) to determine the scaling of the noise term. For \ = 0,
scale invariance is evidently achieved for z9 = 2 and x9 = —3/2. When \ # 0, a sim-

ple scaling solution is no longer possible. Power counting consists of simply looking

28

at the effective scaling of the nonlinear term, which is b*°+#°-? = -3/? in our case.
As 6 is increased, i.e., the length scale considered gets larger, the effective strength of
the nonlinear term decreases, therefore becoming irrelevant in the renormalization
group sense. The power counting argument can be made rigorous by constructing a
perturbation expansion for the nonlinear term, and showing that the perturbation
becomes irrelevant in the limit of small wavevectors [3]. It is easy to check that
all other analytic nonlinear terms consistent with the symmetries one might add
to Eq. (1.38) are similarly irrelevant. Relevant nonlinearities will therefore only be
expected to appear if the conservation law or the parity symmetry is broken. This

is the case with the random-driven Burgers equation, which has the form:
Oy + ydsy = DAZy + Aen, (1.41)

where the noise correlation is as in Eq. (1.39). Power counting now shows that
the nonlinear term scales as b'/?, and is therefore relevant at large length scales.
Eq. (1.41) has been analyzed by Forster, Nelson, and Stephen [58] using dynamic
renormalization group theory.

Returning to the linear Langevin equation Eq. (1.38), we can now make sev-
eral predictions. The most important quantity we can compute from the Langevin

equation is the space-time correlation function:

S(r,t) = ((y(r’ +7, +t) — (y))(y", #7) — (y)) (1.42)

where averaging is over r’ and ¢’. The Fourier transform of S(r,t) is the dynamic

structure factor S(k,w):
S(k,w) = / ” dk / © dus e728) Sp, 4), (1.43)

The Wiener-Khinchin theorm [55] shows that S(k,w) = (|y(k,w)|?). Using this

29

10° Errit Li T T TTTTTy] T T TP TTTE T T TS

E 0.3 crrpreeeprreperpery :

107 F 4 =

F 0.2 |~ a 5

eeloo£ ; =

10 E 0.1 am E

1073 L of CO hbbid a

= E 0 0.1 0.2 0.3 0.4 0.5 E

wh 1 j

1074 Ee ae

F Model A E

10° a=0.4, b=2.0 4

E L=8192 7

10°° 3

107° rial 1 i it rroitil ik L 1 croicil 1 J |
1078 107 107!

Figure 1-6: Static structure factor.
relation in Eq. (1.38), we find directly:

ADK‘

(1.44)
As shown in Figure 1-6 and Figure 1-7, we find numerical results consistent with:

S(k,t) = S(k)e7P#*, (1.45)
S(k) = Ak?, (1.46)

Performing Fourier transforms in time, these results agree with Eq. (1.44). We also
find that the effective diffusion constant D is somewhat larger than the diffusion
constant expected from the linearized equation.

The vanishing of the static structure factor S(k) as k — 0 for model A is related

to the absence of a term linear in k in front of the noise term. This is due to

30

TT TT TOT ET | TT TF | TT TT TET T

EF Model A 4

a=0.4, b=2.0
L=1024

S(k,t)/S(k)

107!

k=14/L

PeTT TTT
roti

k=40/L k=20/L

1072 iui | | | en on a | |e ee eae | { oj) i 4 i Joo of 2

0 200 400 600 800 1000
Time

Figure 1-7 Dynamic structure factor.

the particular form of model A, and is not a result of parity symmetry. Indeed,
performing the same analysis leading to Eq. (1.37) on model B, we find that the

expected Langevin equation now takes the form:
Oy = Dd?y + Oqn, (1.47)

with the correlation of 7 still given by Eq. (1.39). The structure factor for model
B will therefore have a factor of k? missing compared to model A. A numerical
computation of S(k) for model B (cf. Figure 1-8) indeed shows that S(k) — const.
ask — 0.

To summarize, we have derived a Langevin equation for our coupled map models,
which can correctly predict the long wavelength form of the dynamic structure factor.
We have shown that the conservation law and parity symmetry of the system are

sufficient to exclude nonlinear terms from the long wavelength description. As far as

31

T T T T T rT li T T T ret T T ly T | T T T T

” Model B
0.5 - —
[- a=0.2, b=-0.6 4
L c=3.1416 |
L=256
O Jt i | | oj of 1 i Es Dene SN | J} | 1 | a oe
0 0.1 0.2 0.3 0.4 0.5

Figure 1-8: Static structure factor for model B.

two point correlation functions are concerned, we have been able to replace a very
complicated nonlinear deterministic system by a simple stochastic one. A similar
analysis has been carried out by Zaleski [60], who showed numerically that at long
wavelengths the deterministic Kuramoto-Sivashinski equation Eq. (1.17) has two

point correlation functions like the random-driven Burgers equation Eq. (1.41).

1.5 Linear Response and Fluctuation-Dissipation
Theorem

In the previous section, we derived a Langevin equation which contained a Gaussian
noise term. In thermodynamic applications, this noise typically has a thermal origin,
and this provides us with additional information. A famous example is the case of

Brownian motion of a single particle in one dimension in a fluid at temperature T

32

[55]. The Langevin equation for this system takes the form:

dv

me + Vet) = n{t), (1.48)

where the mass of the particle has been scaled to unity, 7 is a damping constant,
and 7 is a Gaussian random force due to the fluid, which has zero mean and a second

moment given by:

(n(t)n(t')) = Pé(é - ¢’), (1.49)

where I is at this stage undetermined. Performing ensemble averaging we immedi-

ately obtain:
(u(t)?) = vee? + 5 —e~ 77"), (1.50)

where vp is the initial velocity. At large times, the particle will be in equilibrium
with the fluid, and (v(co)?) = kgT, where kg is Boltzmann’s constant. We thus

obtain for the noise amplitude:

T = 2QykpT. (1.51)

This equation is a simple example of the fluctuation-dissipation theorem [61].
The general fluctuation-dissipation theorem (FDT) is a relation between equilib-
rium correlation functions and linear response functions. In Fourier space, we can

define the linear response function x(k,w) by the relation:

(y(k,w)) = x(k, w)h(k, w), (1.52)

where fh is a small field applied to y. The FDT for near equilibrium thermodynamic
systems provides a relation between the temperature T,, the dynamic structure factor

S(k,w), and the susceptibility y(k, w):

wS'(k, w)

= Smith) (1.53)

33

In near equilibrium systems this relation is true for all k and w. Using the dispersion

relation [2]

x(k,w = 0) = Pf deb) (1.54)

where P denotes the principal value, we can integrate Eq. (1.53) over w and relate

T to the static structure factor and static response functions:

S(k,t = 0)

T= (kw = 0)"

(1.55)

The FDT relation Eq. (1.53) was originally derived for Hamiltonian equilibrium
systems, but several authors have suggested generalizations to more complicated
situations. A particularly interesting idea is due to Leith [62], who proposed to apply
the FDT to climate prediction. For example, knowledge of the natural fluctuations
in the atmosphere would be sufficient to predict the climatic response to a small
increase in carbon dioxide concentration, with the help of the FDT. However, the
relations derived by Leith are only exact for a system of harmonic oscillators. In
addition, determining the covariance matrix of climate fluctuations to sufficient
accuracy to usefully apply the FDT may be a difficult problem [63].

In another attempt at generalizing the FDT, Hohenberg and Shraiman [53] have
suggested that one might be able to use the FDT to define a useful ‘temperature’
even for nonequilibrium systems, by simply using the right hand side of Eq. (1.53),
perhaps in the limit k — 0 and w — 0. If this limit exists and indeed defines an
effective ‘temperature’ T, even for nonequilibrium systems, there may be several
useful applications. For example, T, might provide constraints on the effective noise
term in the Langevin equation. Another useful application for T, would be the
analog of a thermometer: If two coupled systems have the same value of T., they
should be in equilibrium, i.e., there should be no net flow of the conserved quantity
between them.

Accordingly, we study the response of a coupled map lattice to a small applied

34

4 TOY TT | TOE TT I TORT CT TTT FUT TT FT FT
2 3b
ms L.
56 L
rr bee 5 + } i +
a 2 Model A —|
oO L. 4
wn
3 L.
a) | a=0.4, b=2.0 |
Oo
5 L 4
S tp _
eB) i 4
@) | 1 i i 1 | 1 1 is i | 1 i 1 1 | 1 | i 1 | i L 1 1 ]
0) 0.1 0.2 0.3 0.4 0.5

Figure 1-9: Static susceptibility.

field, chosen so as to preserve the conservation law and reflection symmetry of the

map. The equation of motion Eq. (1.18) is thus modified to:
yeh = yf + outa oP) — YP YiLa) + AE — QAP + Ry, (1.56)

where we concentrate on static fields h? = hcos(2rjk) with fixed wavenumber k.

The static susceptibility y, is then:
. 2
Xe = lim 2ue) (1.57)

where averaging is over n. Numerical calculations of x, for model A and model B
are shown in Figure 1-9 and Figure 1-10 respectively. At long wavelengths, we find
xz to be essentially independent of k.

This result can easily be understood in terms of the Langevin equation by includ-

35

rn ee FUT TT | tT TT TTT TT T T TT

10 }- >

r Model B 7

a=0.4, b=—0.6
c=3.1416

Static Susceptibility y(k)

6) |) i 4 | ee oe | | | ae ee oe | | i) tf J | it ot
QO 0.1 0.2 0.3 0.4 0.5

Figure 1-10: Static susceptibility for model B.

ing an additional term 02h on the right hand side of Eq. (1.38). Since the applied
field h should be uncorrelated with the effective internal noise 7, one immediately
finds that y(k) = 1/D. We have checked the consistency of this result by comparing
the values of D obtained from x(k) and S(k,w) for model A (cf. Figures 1-9 and
1-7). We find numerically that D = 0.45 + 0.02 from y(k) and D = 0.470 + 0.008
from S(k,w), showing good agreement.

Returning to the fluctuation-dissipation ratio Eq. (1.55), we now deduce for
model A that T, ~ k?. In the long wavelength limit, this effective temperature
always vanishes, indicating that T, does not have the desired properties. At the
simplest level, the k? dependence of the effective temperature in model A comes from
the extra gradient in the noise term of Eq. (1.38) as compared with the conventional
ansatz for diffusional Langevin equations. However, the breakdown of the constant

ratio between response and correlations can be traced to a more general flaw in the

36

formulation, namely our lack of knowledge of how the applied field couples to the
probability functional of the system. For Hamiltonian systems, the applied fields
are known to appear in the thermodynamic equations together with their conjugate
fields; e.g., a magnetic field 7 appears in the internal energy through a relation of
the form HdM, where M is the magnetization. We therefore also know the form
of the probability distribution e~*". For nonequilibrium systems, we do not have a
corresponding formalism.

A simple example may be useful here to illustrate the problem. Consider the
long wavelength description of the dynamics of a single conserved variable in an
equilibrium system. Following Hohenberg and Halperin [3], the thermodynamics

may be defined through the entropy S and its density s:
S= Sot [i s'6a(2) + 5/5 [6a(2))’ — h(x)6q(x)}d4x, (1.58)

where s’ = (08/0q)|g=q, and K = (878/0q?)|qg=q, (which plays the role of an elastic
constant) and we have expanded in small fluctuations 6g = q — qo about a uniform
equilibrium value go. For conserving dynamics, one obtains from Eq. (1.58) the

equation of motion:

6S
éq = WG) + Vn(z,t)
= 7KV6q + Vn(z,t) — yV7A(z, t), (1.59)

where y is a kinetic coefficient, and 7 is a Gaussian stochastic noise term with
strength F:
(n(x, t)n(a',t’)) = F(x — 2')6(t — ¢’). (1.60)

For an applied field h(x,t) = hy e**-* the average response 6q = dqy,e1*?

37

calculated from Eq. (1.59) yields a susceptibility:

6 dw
k =
x( w) his
vk?
= ———_,. 1.61
—iw + yKk? (1.61)
Similarly, one obtains the dynamic structure factor:
S(k,w) = ((5dbw)”)
keF
= ——— a 1.
w? + (Kk)? (1.62)
and consequently the static structure factor:
S(k) = ((6qx)2) = =. (1.63)
vi
However, averaging over the probability distribution e~*/" shows that
((6qx)°) = 2T/K (1.64)
so that the noise strength is related to the thermodynamic temperature via
F = 2yT. (1.65)

Using Eqs. (1.61), (1.62), and (1.65), the dynamic fluctuation-dissipation theorem
Eq. (1.53) is immediately verified in the long wavelength limit.

Let us now repeat this analysis in the way we are forced to for a nonequilibrium
system. Instead of defining the coupling to an external field through Eq. (1.58), we
need to introduce a field h into the dynamics; and we do not know the division of

the diffusion constant D into stiffness and kinetic parameters:

6q = DV*5q + Vn(2,t) — V?A(2, t). (1.66)

38

From this we obtain:

- _ 6 kw _ k?
X(F,w) 7 Pes ~ —tw + Dk?
k°F

so that for the proposed definition of an effective temperature

wS'(k, w)

= 1.68
Dim(x(k, 0) nes)
we find T, = F/2.
Comparing with Eq. (1.39) for the equilibrium system we obtain:
T. = 7T (1.69)

depending on the kinetic coefficient y. Without further knowledge about y and its
k and w dependence, T, cannot be expected to measure T. These are exactly the
results found for our coupled map systems.

Although in the case of model B, S(k) and x(k) both approach a constant as
k + 0, so that S(k)/y(k) — constant, we would not expect this value to have the
same importance as in equilibrium systems. It is also clear comparing Figure 1-8

and Figure 1-10 that S(k)/y(k) is k dependent for large k.

1.6 Thermodynamics

1.6.1 Postulates of Thermodynamics

In section 1.5 we have seen that the fluctuation-dissipation theorem is not directly
applicable to nonequilibrium systems. We therefore turn to an alternative point of

view, which only depends on the concept of thermodynamic equilibrium, rather than

39

on small fluctuations about equilibrium.

Thermodynamics can be formulated as a series of basic postulates which are
independent of the presumably underlying concepts of statistical mechanics. Indeed,
the foundation of statistical mechanics is still a subject of intense investigation
[65], while thermodynamics is logically complete. There are alternative ways of
formulating the postulates of thermodynamics, and we will follow very closely the
presentation by Callen [64], with comments on how the postulates may relate to our

coupled map models.

Postulate I There exist equilibrium states of simple systems that, macroscopically,
are characterized completely by the internal energy U, the volume V, and the mole

numbers N,, No,---,N, of the chemical components.

We can hypothesize that our coupled maps are ‘simple systems’, whose macro-
scopic states are completely characterized by the conserved quantity Q, and the
number of lattice sites N. A ‘simple system’ in this context means that the model
system is isotropic and translationally invariant, as we have ensured with the basic
equations of motion of the coupled map lattices. To experimentally determine and
characterize the equilibrium state is of course a major challenge, and below we will

describe several computer experiments with this goal in mind.

Postulate II There exists a function (called the entropy S) of the extensive pa-
rameters of any composite system, defined for all equilibrium states and having the
following property. The values assumed by the extensive parameters in the absence
of an internal constraint are those that maximize the entropy over the manifold of

constrained equilibrium states.

The question of the existence of an entropy function that is maximized is the
main motivation for our investigations. Consider, for example, two different coupled

map lattices whose conserved quantities are Q, and Qz2 respectively. Suppose that

40

the two maps are allowed to interact in such a way that Q can flow from one
map to the other, with the constraint that Q,; + Q2 = remains constant. Then
maximizing the entropy will select a particular Q{ and Q5. One difficulty of course
is how to implement the interaction. The interaction has to impose no constraints
of its own, so that the dynamics can sample all possible states consistent with the
conservation law in order to find the particular state with maximum entropy. As
we shall see, ensuring this in a computational model can be quite difficult. Another
more obvious requirement is that the entropy of the ‘interaction’ be small compared
to the entropies of the two systems, to ensure the additivity required in the following

postulate:

Postulate III The entropy of a composite system is additive over the constituent

subsystems. The entropy is a continuous and differentiable function of the energy.

Using these postulates, it is now possible to introduce the concept of tempera-
ture. Consider again two systems whose equilibrium states are characterized by the
conserved quantities Q; and Qs. The total entropy S of the two systems is, by addi-
tivity, simply the the sum of the entropies of the two systems: S = S1(Q1)+S2(Qz).
Maximizing S with the constraint that Q; + Qo = Q is constant immediately gives

the result that at equilibrium:

dS; dS

dQ: dQo

= 1.7

defining the temperature T.

1.6.2 Heat Bath Algorithms

In this section we describe some numerical experiments which test the postulates
described in section 1.6.1. One of the consequences of the existence of a temperature
is the transitivity of thermal equilibrium: If X is in equilibrium with Y, and Y is

in equilibrium with Z, then X is in equilibrium with Z. We test this transitivity by

41

using three different coupled maps X, Y, and Z, obtained from model A by different
choices of the nonlinear control parameter 6 (cf. Eq. (1.21)). The lattice sizes Lx,
Ly, and Lz are chosen such that Lz = Ly < Lx. We implement coupling between
two lattices by modifying the equation of motion Eq. (1.18) at the endpoints of the
lattices to contain a current term linear in the difference of the lattice amplitudes
in the two lattices. For example, if X and Y are coupled, the equation of motion

for the endpoint of X becomes:

yet (Lx) = yx(Lx) — ox(y¥X(Lx), vx(Lx -— 1)) + d(yy(1) — yx (Zx)), (1.71)

where different lattices are distinguished by subscripts, and the lattice positions
are given in brackets. We arrange the two lattices on a ring, so that we obtain
analogous equations to Eq. (1.71) for yx(1), yy(1), and yy(Ly). The coupling
strength is parametrized by d.

Starting with a random initial configuration with total Q = 0, we in turn couple
Y and Z to X and iterate the new equations of motion until steady state is achieved.
The steady state is characterized by new values of the conserved quantities Qy and
Qz in Y and Z, while the change in the mean value of Qx should be small because
of the large lattice size Ly.

In Figures 1-11 a, b, and c we illustrate the sequence of events described above
by plotting the time evolution of Q for three situations: The thermalization of two
small lattices Y and Z with a large lattice X is shown in Figures 1-1la and 1-11b;
in Figure 1-11c we test for equilibrium of Y and Z when initialized with the values
of Q generated from the equilibration with X. For convenience, we normalize Q by
the system size and denote the result by q, e.g., gx = Qx/Lx.

The transitivity of equilibrium appears to hold very well, with the small shifts
in gy and qz actually seen being qualitatively consistent with those expected due
to the finite size of X. In Figure 1-11d we show the time average of the local lattice

amplitudes for maps X and Z, accumulated after steady state had been reached.

42

CONS TTTTT] TITTY TT TTY TITTY 08 PITT TT TT TTY TTT TTT

L_ (a) = L (b) Map X =

— Map Y 7 — nn

0.01 }— 0—— =

a a ay=az=0.4 _

L- 4 - by=1.3 4

% 0.005 |-— 70.05 F— bz=2.0 —_

4 Ly=4096 4

0 — ~0.1 aa

-o.oo5 -LLitititititritiiri iy) -~o.15 -LLLLLI I PT rr =|
O 500000 1x10® 1.5x10® 2x10 0 500000 1x10® 1.5x10® 2x10

Time Time

0.08 —7T777 771 TT ro ied

L (ce) Map Y 4 Cy) Map X Map Z =

i = = _

0}+— aa 0.1 }-— —

L 4 ¢ LC a

L ay=0.4 4 é a =

o —0.05 -— by=1.25 — - 0 —

5 Ly=Lz=128 4 9 x 4

Lo a 8 z a

-o.1 fe _ -0.1 fe a

Po Map Z - - |

9.15 Ree bp -ogMtdiriti tii tii his

) 10000 20000 3600 3800 4000 4200 4400

Time Lattice Position

Figure 1-11: Time evolution of g (a) Y coupled to X (b) Z coupled to X (c) Y
coupled to Z. (d) Time average of local means.

43

T FT TUM T T UTTER T TO ETO T T ETT
Map Y
| a é & $ 1
0 = 4
ion
-O.1- +
F Map Z 4
i; g 3 Fy ® 7
L i toil 1 i trol rl i too i J} | iyi
1074 1073 107°? 107! 10°

Figure 1-12: q as a function of coupling strength d.

Except very close to the interface, the local means are constant throughout each
lattice, supporting the idea that each lattice has reached local equilibrium.

We also test whether the equilibrium values of g reached are independent of the
coupling strength d. In Figure 1-12 we show the equilibrium values of gy and qz as
a function of d; no significant dependence is found for d < 0.1.

In addition to the coupling strength d, it is also important to check different
functional forms of the coupling. We therefore replace Eq. (1.71) by:

yxt (Lx) = yk(Lx) — 9x(yk(Lx), y(Lx — 1))
+ dtanh(a(yy(1) — yx (Lx))), (1.72)

where the parameter a parametrizes the nonlinearity of the coupling function. If
there indeed exists a thermodynamic limit, then for large system sizes L and small

coupling strengths d, the equilibrium values of g should be independent of a. In our

44

TT TF TT VT rn a a | rn | T rT TT

EE Rg
= 0017 i 7

Hey

O a=1.0 ay=ay=0.4
X a=3.3 by=1.28 |
t e a=10 by=1.5

0.016 } Z

L t 4

L a oe oe | i |e ne loner | | ji | 3 | Vrrbesened, 1 I 1 i |

0 0.02 0.04 0.06 0.08 0.1

Figure 1-13: q as a function of coupling function.

numerical experiments, we choose two coupled maps X and Y with equal length L,
and with an initial g = 0. For a fixed a and d we then determine the equilibrium
value of q as a function of system size. We extrapolate the results to infinite system
size (the convergence goes roughly as 1/L) and plot the extrapolated values of q as
a function of d and a in Figure 1-13. The error bars are rather large due to the
uncertainties in the extrapolations, but the results seem to converge to a unique
limit.

An alternative test for thermodynamic equilibrium consists of coupling three
maps arranged on a ring, and monitoring the current flowing around the loop. If
the three maps reach equilibrium, there should be no net current flow. In our
numerical simulations, we do find persistent currents, but they are small compared
to the shifts of the means between two lattices, as is evident in Figure 1-14, where we

show a typical time average of the local means. The persistent currents are barely

45

F t T T t T T ij | T 7

0.1 7

Map Z Map X Map Y 4

© L |

3s OF 7

Q L 4
fo)

- ay=ay=az=0.4 4

r by=1.4, by=1.3, bz=2.0 4

|- Ly=800, Ly=Ly=200 4

-O.1b d=0.02 Sn

L 1 1 1 i | i i 1 | L 1 4

0 500 1000

Lattice Position

Figure 1-14: Persistent currents in three maps on a ring.

discernible as a non-zero slope of the profile of the means across each lattice. We
find numerically that the magnitude of the persistent currents scales roughly as d/L,
where L is the system size and d is the coupling strength. A current scaling as 1/L
leads to a constant offset in g, and would question the thermodynamic equilibrium
picture, but our numerical resolution does not allow any firm conclusions.
Assuming that the thermodynamic picture does hold, we can test it further by
checking whether the distribution of g between the three maps is consistent with that
found in the temperature bath experiments where two maps are coupled at a time.
To do this, we increase the size Ly of one of the three lattices, and determine the
corresponding rearrangement of g. The infinite system limit is approached approx-
imately as (1/Lx)!/? (cf. Figure 1-15). Within error bars, the extrapolated means
agree with those obtained from the heat bath experiments described previously.

While the numerical experiments described strongly suggest that thermodynamic

46

0.05 ee ee ee ee a TT T Tord rT

sa)
ne}
can
Yo

om —0.05

~O.1

ree tp fe Ne

—0.15

Figure 1-15: Dependence of gy and gz on Lx.

considerations may be applicable to the coupled map models considered here, the
straightforward application of ‘standard’ thermodynamic ideas is difficult. As an
example of these difficulties, we attempt to apply the Einstein fluctuation theorem
[66]: If there exists an entropy function S which depends only on the conserved
quantity Q, then small fluctuations of @ about its equilibrium value should have a
Gaussian distribution of width 079/0Q?. Measuring the fluctuations as a function
of Q should then provide an indirect way to measure the entropy. Now, if S and Q
are éxtensive quantities, it follows that the mean square fluctuations should scale
as 1/L. However, for model A we find that the fluctuations instead scale as 1/L?.
This is shown in Figure 1-16, where we plot the mean square fluctuations of @ in
map X, while coupled to map Y, as a function of the system size Lx = Ly = L.
This result can be traced to the k? behavior of the structure factor S(k) at long

wavelengths (cf. Figure 1-6). To show this, consider the mean square value ay

47

10°? T T TTTMTf£ T Li T FT TTTT T F T TTTT”
103% - +
a — 4
A C 5
oe a 4
Vv
NN r |
A 1074 & 4
(eg - E
Vv C _
| L 1
iz L |
10°b 4
1078 l |. rotiil i i Liiiil i i |e ce
10! 10° 10° 104

Figure 1-16: Mean square fluctuations as a function of system size.

of the coarse-grained lattice amplitude fluctuations, averaged over N sites, where

1l12% °
ow = ( woe ), (1.73)

where, for simplicity, we assume the fluctuations are about zero mean. By simple
manipulations, oy can be related to the structure factor S(k):

112 sin?(1Nk)

= TET

For model A, S(k) ~ sin?(zk) (cf. Figure 1-6), and therefore oy ~ 1/N?. If an
entropy function indeed exists, the simple Taylor series approach underlying the
Einstein fluctuation theorem appears therefore not to be directly applicable. Model

B does not suffer this difficulty because S(k) ~ const (cf. Figure 1-8).

48
1.6.3. Discussion

In section 1.6.2 we have presented numerical results both in favor and against the
existence of a thermodynamic limit of our coupled map models. One of the main
problems is how to distinguish partitions of Q due to bulk effects from partitions due
to surface effects. As an example how surface effects can occur, we briefly discuss a
system that has been suggested as a counterexample to the idea that the division
of Q is determined by maximizing a bulk entropy [67].

We again start with two maps X and Y of equal length L, but now we place
them alongside each other such that every site in X couples to a corresponding site

in Y. The intra-lattice dynamics is the same as in model A:

yet) = yf) + AGG + 1)) - 2A(y?@) + Aly? - 1), (1.75)
and = yg (i) = yB(i) + folyS(i + 1) — 2fo(yh (0) + fo(y3 (i — 1)). (1.76)

The inter-lattice coupling is introduced by periodically adding the following terms

to the equation of motion:

git (i) yt (4) + af fo(ys(2)) + f(y (é + 1)) — 2fi(yP(2))}, (1.77)
and = ys (i) = yb (é) + of fi(yt (a) + folys(i + 1)) — 2fe(y3(d)}, (1.78)

where a is a constant. It is clear that the inter-lattice coupling violates the parity
symmetry 7 < —2. More important, however, is that a division of Q is likely to occur
without any bulk effects being responsible. An easy way to see this is to consider

the case where the nonlinear functions f are the sum of a linear and a random term:

fi(z) = art+bhié, (1.79)
fo(z) = aga + bef, (1.80)

where the €; on different lattice sites are uncorrelated. Taking the time average of

49

Eq. (1.77), setting (y?(2)) = qi and (y$(7)) = qe, we find that equilibrium is attained

when:

2g2 — 4191 + (bz — b1)(€) = 0. (1.81)

We will therefore quite generally find gq; 4 go, and it is obvious that no bulk effects
are responsible.

Returning to couplings Eq. (1.71) and Eq. (1.72) investigated in section 1.6.2,
we see that they do not have the flaw of the system described above. However, there
is another potential problem with these couplings, and we therefore introduce yet
another candidate for the flux term between the two systems.

We would expect that in the limit of uncoupled systems (i.e., the coupling
strength goes to zero), the coupling operator itself should average to zero. One
way to ensure this is to require the coupling term to depend on the local currents of
the uncoupled system, rather than on the local amplitudes. We can implement this

by replacing Eq. (1.72) by:

yx (Lx) = yk(Lx) — ox(ye(Lx), ye(Lx - 1) (1.82)

+ C (9x (yx (Lx), yx (Lx — 1)), gv (y?-(2), 9f-(1))),

where C(z,y) is an arbitrary function satisfying C(z,y) = —C(y,x) to maintain
parity symmetry 7 —2.

One difficulty with implementing Eq. (1.82) is that only the nonlinear com-
ponents of the coupling function C(z,y) can lead to macroscopic currents. In our
numerical simulations, we are hampered by large fluctuations, which often cause the
coupled maps to enter a periodic orbit. We have therefore not been able to obtain
conclusive results with this coupling, although there are indications that equilibrium
is established for a different division of Q than for the coupling Eq. (1.71), casting
some doubts on our heat bath algorithms.

We have been able to obtain some results for model B, with C(z,y) = dzy,

50

where d is the coupling strength. This coupling does not satisfy parity symmetry,
but we believe that in the limit d — 0 the same results as for a parity symmetric
coupling should be obtained. We have used these results to investigate a rather
speculative question, namely whether the rate of information production should
be at an extremum when two maps are in equilibrium. The rate of information
production is called the Kolmogorov entropy h [26]. It can be shown that the sum
of the positive Lyapunov exponents provides an upper bound to h [68], and it is
a strong possibility that an equality holds for systems having a physical measure
[26]. For a given map, we therefore assume that we can compute / as a function
of Q from the Lyapunov spectrum. We then compute h, and he for two maps with
different map parameters, for conserved quantities Q; and Qo = Q. — Qi, where Q,
is fixed. We then plot hj +2 versus g; = Qi/L, and check whether the plot shows an
extremum at that value of gq, which arises when the two maps are coupled together
with an initial value of Q, for the conserved quantity. We show one particular result
in Figure 1-17, where the initial value of q was 0.5.

The plot shows indeed an extremum for q; + 0.68 + 0.01, compared with a value
of q, * 0.65 + 0.03 observed when the maps were coupled. The large error bars
on the latter value are due to the large fluctuations during the evolution of the
coupled system. This is qualitatively consistent with Figure 1-17, which shows that
hy + he depends only very weakly on q;. We emphasize that these results have to
be considered as rather speculative, since there is no clear understanding whether
the Kolmogorov entropy should take on an extremum. We also have not been able

to consistently reproduce the behavior of Figure 1-17 in model A.

1.7 Analysis of an Experiment

We now present a more detailed discussion of the convection experiment by Ciliberto

and Bigazzi [9] already mentioned in section 1.1.1 above. Recall that the experiment

ol

27.8 T 7 T T ] t T T T T T T T U T T t
| Model B :
27.6 Fb a
oo 274 =
& [" 4
+ [ a,=0.4, b,=—-0.8 1
< 2a72b ap=0.4, be=-0.4
i C,=Co=3.1416 |
27 L,=L,=32 ~
26.8 i L i f | | i 1 1 1 | 1 i 13 ! | i 1 it i :

0.4 0.5 0.6 0.7 0.8
qi

Figure 1-17 Kolmogorov entropy as function of q.

consists of convection in an annulus, with convection rolls along the radial direction
(cf. Fig. 1-18).

In the experiment, the vertically averaged temperature gradient perpendicular
to the radial direction is measured at N equally spaced points located at a fixed
radial distance from the center of the cell. N is typically between 128 and 256, and
the sampling rate is approximately 1/10 of the oscillation time of one convection
roll.

If we denote the experimental measurements as y"(i), where 7 labels the position
of the probe, and n labels the sampling time, we can interpret yz) as forming a
one-dimensional coupled map lattice, with a discrete time evolution given by the
sampling rate. We wish to explore the similarity between this ‘lattice’ constructed
from experimental observations, and the coupled map models we have studied in

this chapter. An important difference between the experimental situation and our

Tt AT

Figure 1-18: Convection experiment in an annulus.

models is that the experimental data points are strongly driven externally by the
motion in the whole convection cell, while our models are entirely self-interacting.
We assume that the strongly nonlinear map functions in our model can emulate this
external driving process.

The experimental y”(z) also satisfy an approximate conservation law, since the
temperature gradient integrated around the cell vanishes. Because of the discrete
sampling, this conservation law is only approximate.

In the experiment, the temporal fluctuations of y"(7) at a fixed position 7 show
a complicated distribution. By contrast, the fluctuations of the Fourier transforms
y"(k) at a fixed wavevector k show approximately a Gaussian distribution in the
chaotic regime. This result agrees with our discussion in section 1.3.2 of the applica-
bility of the central limit theorem due to finite correlation lengths. Below the onset
of spatiotemporal chaos, the distribution of the y"(k) is not Gaussian, as we would
expect from the spatial periodicity of convection roll patterns at small driving.

Finally, the authors measure quantities they term ‘energy’ and ‘entropy’. The
y; 8. PY

o3

energy E” is defined as:
E* = (y"(i)) (1.83)
i=1

Using Parseval’s theorem, we interpret the time average E of E” as the static struc-
ture factor S(k) summed over k. From our Langevin description of the structure
factor in section 1.4, we would therefore expect E to provide a measure of the ‘noise’
in the system, and to increase with the nonlinearity. In the experiment it is indeed
found that in the chaotic regime F increases approximately linearly with the applied
vertical temperature gradient.

The entropy o” is defined as:

N/2

3 "(k) loa "(K)), (1.84)

c= —

do
where ¢"(k) = |y"(k)|?/E", and oo = log(N/2). A rough interpretation of o = (o”)
is as an indicator how equipartitioned the static structure factor S(k) is as a function
of k. When the structure factor is a constant independent of k, o takes on its
maximum value of 1, while a structure factor peaked at a single value of & results
in the minimum value of 0.

In the chaotic regime, @ is found to remain close to its maximum value, indicating
a static structure factor similar to model B. (cf. Figure 1-8). Unfortunately, the
authors do not report the dynamic structure factor, which would provide a more
stringent test of the Langevin description of section 1.4.

While these observations do not provide rigorous quantitative tests, they do
indicate that some of the results obtained with our coupled maps models may apply
to physical experiments. Whether our heat bath algorithms of section 1.6.2 can
be tested in an actual experiment remains to be elucidated and is certainly worth

investigating.

54

1.8 Conclusion

Our study of coupled maps with conservation laws has led to some interesting re-
sults, while at the same time raising several questions which require further investi-
gations. We have derived simple Langevin equations which can adequately describe
many long wavelength features of strongly chaotic systems. We have shown why the
fluctuation-dissipation theorem will in general not allow the definition of a useful
effective temperature in a nonequilibrium system. Finally, we have presented em-
pirical evidence that some ideas from equilibrium thermodynamics may carry over
to chaotic systems without an underlying Hamiltonian.

The most important question that our work has not been able to resolve com-
pletely is how to conclusively distinguish surface effects from bulk effects in the
distribution of a conserved quantity between two chaotic systems. Even if this ques-
tion can be resolved, and something analogous to the maximum entropy principle
of thermodynamics can be established for chaotic systems, one needs to identify
experimentally measurable quantities to test these ideas. We have only been able to
present rather preliminary results using the Kolmogorov entropy (cf section 1.6.3),

which can in principle be measured experimentally.

Appendix

In this appendix we briefly outline the numerical procedures we used to compute
Lyapunov exponents. Since in our coupled map models the equations of motion and
therefore the Jacobian matrices are known, we could, in principle, use Eq. (1.13) to
determine the Lyapunov spectrum. However, this would quickly lead to numerical
difficulties due to the exponential divergence of trajectories, which results in singular
values of the products of the Jacobians. A better approach consists of following
the evolution of an orthonormal set of vectors, which is periodically renormalized

to avoid collapse onto the most rapidly growing directions [69, 70]. This can be

59

accomplished with two mathematically equivalent, but numerically slightly different
methods.

In the first method, one starts with a random orthonormal set of vectors, and
operates on that set successively with the instantaneous Jacobians. The vectors
will therefore start to align themselves along the most rapidly growing directions.
One then periodically applies the Gram-Schmidt renormalization to these vectors.
During Gram-Schmidt renormalization, the. space spanned by the first k vectors is
left unchanged, where k ranges from one to the number of vectors available. The
Gram-Schmidt procedure therefore replaces an almost degenerate complete set by an
equivalent orthonormal one, solving the computational difficulties. The first vector is
not rotated at all during renormalization, and can therefore align itself onto the most
rapidly growing direction. The second vector is orthogonalized to the first vector,
and can therefore only seek the second-most rapidly growing direction. Likewise,
each succeeding vector will be orthogonal to all the previous ones. The logarithms of
the scale factors needed to normalize the vectors to unit length, averaged over time,
converge to the Lyapunov exponents. One potential difficulty with this method is
that the Gram-Schmidt procedure gets more and more inaccurate numerically as
the number of vectors to be orthonormalized is increased. _

The second method is based on the QR algorithm [26] and is numerically more
robust. The QR algorithm factors a matrix A into an orthonormal matrix Q, and an
upper triangular matrix R. The matrix Q forms an orthonormal basis for A, and is
identical to that obtained by the Gram-Schmidt procedure. As in the first method,
one starts with a random set of orthonormal vectors, and applies the Jacobian to that
set. The resulting matrix is factored according to the QR algorithm, and at each
succeeding iteration, the Jacobian is applied to the orthonormal matrix resulting
from the previous QR factorization. This generates a sequence of upper triangular
matrices, R,, R2,.... The product of these matrices is also upper triangular, and

the logarithms of the diagonal entries, averaged over a large number of R;, converges

56

to the Lyapunov exponents. Evidently only the diagonal entries of the R; need to
be kept during simulations. The numerical complexity of both algorithms increases
as the cube of the number of exponents computed.

The QR algorithm is computationally slower than the Gram-Schmidt procedure
by approximately a factor of two, and for fewer than ~ 50 vectors we have ob-
tained satisfactory results with a modified Gram-Schmidt algorithm [71], without
appreciable sacrifice in accuracy.

It is also possible to compute the Lyapunov spectrum from a time series of
experimental data [72], although it appears that only the positive exponents can be
reliably obtained. The complexity of these algorithms increases exponentially with

the dimension, and they are therefore only applicable to low-dimensional attractors.

ov

References

{1] D. J. Tritton, Physical Fluid Dynamics, Oxford University Press (1988).

[2] D. Forster, Hydrodynamic Fluctuations, Broken Symmetry, and Correlation
Functions, W. A. Benjamin Inc. (1975).

[3] P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 435 (1977), and

references therein.
[4] J. L. Lebowitz, Physica 140A, 232 (1986).

[5] For a review of the early experimental and theoretical work on Rayleigh-Bénard
convection, see S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability,

Oxford University Press (1961).
[6] A. C. Newell and J. A. Whitehead, J. Fluid. Mech. 38, 279 (1969).

[7] For an extensive review, see M. C. Cross and P. C. Hohenberg, Pattern Formation

outside of Equilibrium, to be published in Reviews of Modern Physics.

[8] G. Ahlers and R. P. Behringer, Prog. Theor. Phys. (Kyoto) Suppl. 64, 186 (1978);
J. P. Gollub, A. R. McCarrior, and J. F. Steinmann, J. Fluid Mech. 125, 259
(1982); S. Ciliberto, J. Phys. Cond. Matter 2, SA483 (1990).

[9] S. Ciliberto and P. Bigazzi, Phys. Rev. Lett. 60, 286 (1988).

[10] E. N. Lorenz, J. Atmos. Sci. 20, 130 (1963).

08

[11] J. H. Curry, J. R. Herring, J. Loncaric, and S.A. Orszag, J. Fluid Mech. 147,
1 (1984).

[12] Y. Kuramoto and T. Tsuzuki, Prog. Theor. Phys. 55, 356 (1976); G. I. Sivashin-
sky, Acta Astronautica 6, 569 (1979).

[13] P. Manneville, “Lyapunov exponents for the Kuramoto-Sivashinsky model,”
in Macroscopic modeling of turbulent flows, ed. O. Pironneau, Lect. Notes in

Physics 230, Springer Verlag (1985).
[14] S. Ciliberto and M. Caperoni, Phys. Rev. Lett. 64, 2775 (1990).

[15] Spatiotemporal intermittency in hydrodynamics refers to the random switching

between laminar and turbulent flow. See, e.g., Ref. [1] above, p. 285.
[16] M. Faraday, Phil. Trans. R. Soc. London 121, 299 (1831).

[17] For a review of parametrically excited surface waves, see J. Miles and D. Hen-

derson, Ann. Rev. Fluid Mech. 22, 143 (1990).

[18] N. B. Tufillaro, R. Ramshankar, and J. P. Gollub, Phys. Rev. Lett. 62, 422
(1989). See also J. P. Gollub and R. Ramshankar, “Spatiotemporal chaos in
interfacial waves,” in New Perspectives in Turbulence, ed. L. Sirovich, Springer

Verlag (1991).

[19] For a review of maps on the unit interval, see P. Collet and J. P. Eckmann,

Iterated Maps of the Interval as Dynamical Systems, Birkhauser, Boston (1980).

[20] For an extensive discussion of the logistic map and its history in modeling

population dynamics, see R. M. May, Nature 261, 459 (1976).

[21] M. J. Feigenbaum, J. Stat. Phys. 19, 25 (1978), and J. Stat. Phys. 21, 669
(1979).

59

[22] Periodic-doubling chaos has been observed in several convection experiments,
e.g., J. P. Gollub and S. V. Benson, J. Fluid Mech. 125, 15 (1980); J. Maurer
and A. Libchaber, J. de Phys. Lett. 41, L515 (1980).

[23] For a review, see N. B. Abraham, J. P. Gollub, and H. L. Swinney, Physica D
11, 252 (1984).

[24] W. Lauterborn and E. Cramer, Phys. Rev. Lett. 47 1445 (1981); M. Kitano,
T. Yabuzaki, and T. Ogawa, Phys. Rev. Lett. 50, 713 (1983).

[25] See, e.g., L. Glass, M. R. Guevara, A. Shrier, and R. Perez, Physica D 7, 89
(1983).

[26] For an introduction to the mathematical foundations of the theory of chaos,

see J.-P. Eckmann and D. Ruelle, Rev. Mod. Phys. 57, 617 (1985).

[27] S. M. Ulam and J. von Neumann, J. Bull. Am. Math. Soc. 53, 1120 (1947); J.
von Neumann, Collected Works 5, 768 (1961).

[28] Physical measures are discussed in Ref. [26], p. 626.

[29] See, e.g., M. Hénon, “Numerical exploration of Hamiltonian systems,” in

Chaotic Behavior of Deterministic Systems, ed. G. Iooss, North-Holland (1983).
[30] V. I. Oseledec, Moscow Math. Soc. 19, 197 (1968).
[31] B. Mandelbrot, The Fractal Geometry of Nature, Freeman, New York (1983).
[32] K. Kaneko, Physica D 37, 60 (1989) and references therein.
[33] K. Kaneko, Prog. Theor. Phys. 7, 1033 (1985).
[34] J. D. Keeler and J. D. Farmer, Physica D 23, 413 (1986).

[35] H. Chaté and P. Manneville, Physica D 32, 402 (1988).

60

[36] Y. Pomeau, “Front motion, metastability, and subcritical bifurcations in hy-
drodynamics,” in Spatio-temporal coherence and chaos in physical systems, eds.
A. R. Bishop, G. Griiner, and B. Nicolaenko, North-Holland (1986); See also
Ref. [54] below.

[37] See, e.g., P. Grassberger and T. Schreiber, Physica D 50, 177 (1991).

[38] L. A. Bunimovich and Y. G. Sinai, Nonlinearity 1, 491 (1988)

[39] J. P. Crutchfield and K. Kaneko, Phys. Rev. Lett. 60, 2715 (1988).

[40] T. Bohr, G. Grinstein, Y. He, C. Jayaprakash, Phys. Rev. Lett. 58, 2155 (1987).

[41] G. Grinstein, Y. He, C. Jayaprakash, and B. Bolker, Phys. Rev. A44, 4923
(1991).

[42] H. Sompolinsky, A. Crisanti, and H. J. Sommers, Phys. Rev. Lett. 61, 259
(1988).

[43] M.S. Bourzutschky and M. C. Cross, Chaos 2, 178 (1992).

[44] See, e.g., the review by J. P. Crutchfield and J. D. Farmer, Physics Reports 92,
45 (1982).

[45] B. A. Huberman and J. Rudnick, Phys. Rev. Lett. 45, 154 (1981).
[46] A. Renyi, Probability Theory, North Holland (1970).

[47] J. L. Kaplan and J. A. Yorke, “Chaotic behavior of multi-dimensional difference
equations,” in Functional Differential Equations and Approximation of Fixed
Points, eds. H. O. Peitgen and H. O. Walther, Lect. Notes in Math. No. 730,
Springer Verlag (1979).

[48] D. Farmer, E. Ott, and J. A. Yorke, J. Diff. Eqns. 49, 185 (1983).

[49] F. Ledrappier, Comm. Math. Phys. 81, 229 (1981).

61

[50] V. N. Shtern, Phys. Lett. 99A, 268 (1983).

[51] For a detailed discussion of dimension algorithms, see J. Theiler, Quantify-
ing Chaos: Practical Estimation of the Correlation Dimension, Ph.D. Thesis,

California Institute of Technology 1988 (unpublished).
[52] D. Ruelle, Commun. Math. Phys. 87, 287 (1982).
[53] P. C. Hohenberg and B. I. Shraiman, Physica D 37, 109 (1989).

[54] P. Manneville, Dissipative Structures and Weak Turbulence, Academic Press

(1990).

[55] N. G. Van Kampen, Stochastic Processes in Physics and Chemistry, North-
Holland (1981).

[56] U. Deker and F. Haake, Phys. Rev. A11, 2043 (1975).

[57] R. Graham, “Statistical Theory of Instabilities in Stationary Nonequilibrium
Systems with Applications to Lasers and Nonlinear Optics,” Springer Tracts in

Modern Physics, Springer Verlag (1973).
[58] D. Forster, D. R. Nelson, and M. J. Stephen, Phys. Rev. A16, 732 (1977).
[59] M. Kardar, G. Parisi, and Y-C. Zhang, Phys. Rev. Lett. 56, 889 (1986).
(60] S. Zaleski, Physica D 34, 427 (1989).
(61] H. B. Callen and T. A. Welton, Phys. Rev. 83, 34 (1951).
[62] C. E. Leith, J. Atmos. Sci. 32, 2022 (1975).
[63] T. L. Bell, J. Atmos. Sci. 37, 1700 (1980).

[64] H. B. Callen, Thermodynamics, John Wiley & Sons (1961).

62

[65] For a recent review, see The dynamic origin of increasing entropy, M. C.

Mackey, Rev. Mod. Phys. 61, (1989).
[66] See, e.g., L. Landau and E. M. Lifshitz, Statistical Physics, Pergamon (1980).
[67] J. Miller, private communication.
[68] Ya. B. Pesin, Russian Math. Survey 32, 55 (1977).

[69] I. Shimada and T. Nagashima, Prog. Theor. Phys., 61,1605 (1979); G. Benettin,
L. Galgani, and J. M. Streclyn, Phys. Rev. A14, 2338 (1976).

[70] A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano, Physica D 16, 285
(1985).

[71] G. H. Golub and C. F. van Loan, Matric Computations, Johns Hopkins Uni-
versity Press (1989).

[72] J.-P. Eckmann, S. O. Kamphorst, D. Ruelle, and S. Ciliberto, Phys. Rev. A34,
4971 (1986); R. Brown, P. Bryant, and H. D. I. Abarbanel, Phys. Rev. A43,
2787 (1991); and Ref. [70] above.

63

Chapter 2

Self-Organized Criticality

An einer Theorie ist es wahrhaftig
nicht ihr geringster Reiz,

daf sie widerlegbar ist:

gerade damit zieht sie feinere Kopfe
an.

F. Nietzsche

64

2.1 Introduction and Outline

Some of the most interesting and difficult problems in condensed matter physics
involve a large range of spatial and/or temporal scales. Quite often in these problems
physical variables display scaling behavior which can be expressed mathematically
by an homogeneity relation. An example is the Gibbs free energy G(T, H) of a
magnetic system, when the temperature T approaches a critical temperature T, at
which a second order phase transition occurs. Introducing the reduced temperature
e=(T-T,)/T., and letting G(T, H) — G(e, H), it is found [1] that G(e, H) obeys

a generalized homogeneity relation:
G(\%e, ®H) = AG(e, H), (2.1)

where a and @ are critical erponents and X is an arbitrary parameter. Choosing

\ = H-V/8 | we find:
He/8’

Gle, H) = HV/G( 1), (2.2)

so that the Gibbs free energy only depends on the ratio «/H®/®, apart from a mul-
tiplicative factor. The parameter space is therefore reduced from a function of two
variables to a function of one variable, and two critical exponents. Indeed, the
critical exponents depend only on the dimensionality and symmetry of the under-
lying system, not on the details of the microscopic interactions. The fundamental
reason for this is that near the critical point the correlation length diverges, which
effectively makes other microscopic lengths irrelevant. The renormalization group
approach developed by Wilson [2] provides a systematic procedure for computing
the critical exponents, with predictions that generally agree quite well with exper-
iment [3]. It is generally believed that the scaling behavior observed can only be
achieved by tuning one or more parameters, e.g., the critical temperature T,.

The situation in nonequilibrium systems is considerably less clear, for example

the ubiquitous occurrence of ‘1/f’ noise in situations ranging from resistance fluc-

65

tuations in semi-conductors to movements in the stock market remains one of the
biggest puzzles in science [4]. Motivated by the problem of 1/f noise, Bak, Tang, and
Wiesenfeld have proposed a class of cellular automata models, referred to as ‘sand-
piles’, because they mimic sand flowing down an incline [5]. Bak et al. argue that
these models display what they call self-organized criticality, which is quite different
from the criticality in equilibrium statistical mechanics: No tuning of parameters is
necessary, but the scaling behavior arises instead through the propagation of noise
through ‘minimally stable’ states. Consider, for example, a hypothetical incline on
which sand is continuously dropped. Because of friction, the sand can only slide
down the incline when the local slope exceeds a certain threshold. When sand ac-
tually does slide, it rearranges the distribution of sand on the incline. According
to Bak et al., the distribution of sand eventually reaches a minimally stable state
where disturbances can propagate arbitrary distances, leading to avalanches of all
sizes. Numerical simulations of sandpile models indeed show nontrivial power law
distributions, but there is ongoing debate as to whether these power laws really
reflect fundamentally new dynamic principles, or whether they can be understood
within established frameworks.

Central to the problem of self-organized criticality is the question of under what
circumstances scaling behavior can occur. In this chapter we will investigate this
question from several different points of view. In Section 2.2 we review the basics
of dimensional analysis, which often leads to the prediction of scaling relations with
very little theoretical effort. We then describe a simple one-dimensional sandpile
model which motivated us to study a piecewise linear diffusion equation, similar to
one solved by Barenblatt [6]. We use this equation to show how non-trivial scaling
behavior can arise in partial differential equations, with power laws that cannot be
obtained by dimensional analysis alone.

In Section 2.3 we focus directly on sandpile models. We start by providing a

summary of the extensive research that has been done on self-organized criticality.

66

We then describe detailed numerical investigations of two particular sandpile models,
one with and one without a conservation law. We derive a mean field theory for these
models, and comment on the importance of conservation laws and other invariance
principles for the self-organized critical state.

In Section 2.4 we investigate the notorious cellular automaton model ‘Game of
Life’ introduced by Conway [7], which has been suggested to show self-organized
criticality without an underlying invariance principle [8].

We summarize our results in Section 2.5. In the Appendix we provide an elemen-
tary introduction to Boolean algebra and describe efficient computational algorithms

for the simulation of cellular automata.

2.2 Nonlinear Diffusion: A Simple Example

2.2.1 Dimensional Analysis and Scaling

Dimensional analysis plays an important part in both experimental and theoretical
science. Measurements of physical quantities are always expressed as ratios with
certain fundamental units. The complete set of fundamental units used to describe
a given class of phenomena is the system of units of measurement. For example, the
Systéme International (SI) standard provides a set of units to measure all physical
phenomena deemed worth measuring. Evidently the definition of the fundamental
units is quite arbitrary and is usually dictated by convenience [9]. However, once
a system of units has been chosen, the units of any physical quantity are products
of powers of the fundamental units, independent of the magnitude of those units.
For example, the unit of force in the SI system is MLT~?, where M, L, and T are
the fundamental units of mass, length, and time respectively. This is true whether
one measures mass in kilograms (as in the MKS system) or in grams (as in the cgs
system). The MKS and the cgs system really belong to the same class of units,
namely the class that takes M, L, and T as the fundamental units.

67

The subject of dimensional analysis derives mainly from the observation that
the unit of a physical quantity is always expressible as a product of powers of the
fundamental units. Any physical measurement consists of establishing functional

relations of the form:

a= f (dj, Q2,..., Gk, Ak41,---5@n), (2.3)
where the quantities a;,...,@, are called the governing parameters of a [6]. Of these
n parameters, in general only a subset consisting of k members a;,..., a, will have

independent dimensions. The dimensions of the quantity a, and the dimensions
of the quantities a,41,...,@,, can then be expressed as products of powers of the
dimensions of the members of that subset. Using this information, one can then
introduce dimensionless ratios by dividing a, @,41,...,;@n by appropriate products of
a,,...,a,. Let these dimensionless ratios be I, Ih,...,1,—,. Then the Buckingham

Pi theorem [10] states that there exists a function F' such that:

II = F(th,...,Un-x). (2.4)

As a famous example to illustrate the Pi theorem, we briefly describe how G.I.
Taylor deduced the energy released in a nuclear explosion from a set of pictures
showing the radius of the fireball at several instants in time [11]. He assumed that
the radius r of the fireball depends mainly on the energy of the explosion EF, time

t, and the ambient air density p:

r=r(E,t,p). (2.5)

The dimensions of E, t, and p are ML?T~?, T, and ML~ respectively. Obviously
these dimensions are all independent, and from the Pi theorem we deduce that the

right hand side of Eq. (2.4) is a simple constant in this case, since n — k = 0. By

68

forming the product of EF, t, and p that has the dimension of length we obtain:
r/(E®/p) =C, (2.6)

where C is a constant. Plotting log(r) vs. log(t) one should thus obtain a slope of
2/5, as was indeed observed [11]. In order to obtain the actual energy E, however,
one also needs to know the proportionality constant C. The big advantage of a
scaling relation like Eq. (2.6) is that no reference is made to the actual size of the
explosion. We then expect the constant C' to be same in a nuclear explosion as in a
conventional chemical explosion, and that is indeed how Taylor deduced the energy.

There are of course several caveats in applying dimensional analysis, since the
choice of governing parameters usually cannot be deduced from first principles. For
example, in the analysis of the nuclear explosion above, there is no a priori reason
not to include the linear dimension / of the bomb itself as one of the governing
parameters. Applying the Pi theorem, we then find that the right hand side of
Eq. (2.6) is no longer a simple constant, but an arbitrary function of the ratio
l/(Et?/p)'/>. Our intuition is that this ratio is ‘small’ and that we can let it go to
zero, making the arbitrary function a constant. In general this cannot be rigorously
justified, since the limit z — 0 of an arbitrary function f(x) may not exist. Indeed,
non-trivial scaling relations, which cannot be deduced from dimensional analysis
alone, can arise precisely because of subtleties in taking such limits [6].

Another simple example, which will be the starting point of the more complicated
problem discussed in section 2.2.3, is provided by the one-dimensional linear diffusion

equation

Ay(e,t) = Dédzy(z,t), (2.7)
y(z,0) = Qé(x). (2.8)

Here D is the diffusion constant, and the domain z is unbounded. The parameters

69

z, t, D, and Q have dimensions L, T, L?/T, and EL respectively, where we have
assigned y an arbitrary dimension of &. Only three of the four parameters have
independent dimensions, and applying the Pi theorem we find that y should be of

the form

ul2,t) = of(a/VDi), (2.9)

where f is an arbitrary function. Substituting Eq. (2.9) into Eq. (2.7) we obtain an
ordinary differential equation for f(n) in the single variable n = 2/V Dt:

1 1
f" +5nf' +5f =0. (2.10)
2 2
The solution which vanishes at infinity is evidently

f(n) = Aexp(-F1?), (2.11)

where A is at this stage an arbitrary constant. The complete solution to Eq. (2.7)

is now easily found to be

y(z,t) = so exp aD (2.12)

which is of course the well known solution to the linear diffusion equation with
an impulsive initial condition. The main advantage of dimensional analysis in this
simple example was the reduction of a partial differential equation to an ordinary
differential equation. As in the example of the nuclear explosion above, we expect
the general scaling form Eq. (2.9) to remain valid in the long time limit, even if the
initial condition is not a delta function, but spread over some distance /. Again this
corresponds to assuming that the additional dimensionless ratio 1/ VDt becomes
irrelevant for small / and large t.

The time dependent prefactor of Eq. (2.9) will in general depend on the initial

conditions. In section 2.2.3 betow we will attempt to analyze the behavior of a

70

diffusive system in a semi-infinite domain with vanishing slope at the origin and a
vanishing spatial integral. The expected scaling form is then not t~/?f(x//t), but
rather t~3/? f(2/ Vt).

We now describe the ‘local limited’ sandpile model, our motivation for studying

a simple nonlinear version of the linear diffusion equation considered above.

2.2.2. The Local Limited Model

The local limited model is a one-dimensional cellular automaton, with rules which
are very similar to many of the sandpile models studied for self-organized criticality
[13, 12]. The model consists of integral variables h(7), arranged on a one-dimensional
lattice indexed by 2, where 7 ranges from 1 to the system size L. The integer h(?)
can be interpreted as the number of ‘sand grains’ at position 7, where each grain has
a height of unity. The local slope o(7) is defined by the height difference between

two nearest neighbors

a(t) = h(i) — h(@ +1). (2.13)
The boundary conditions are

o(0) = 0 (2.14)

h(i) = 0 fori>L. (2.15)

The dynamics consists of randomly selecting a site ¢ and adding a sand grain by
incrementing h(z). If the local slope o(7) exceeds a threshold value o,, an ‘avalanche’

is initiated by the process:

h(i) = h(i)—n (2.16)

Ati +1) = A(Git1)4+n, (2.17)

71

where n is an integer (if i = L, n grains drop off the pile and h(L + 1) is kept
at zero). Eq. (2.16) is then applied repeatedly to all lattice sites that exceed the
threshold until all the local slopes are below or equal to o,.

Several interesting quantities can then be measured for this model, for example
the distribution of sand leaving the system at 1 = LD, and the distribution of the
number of sites participating in an avalanche. Extensive numerical simulations
indicate that these quantities may have multi-fractal distributions [13, 14].

It is convenient to recast the entire dynamics in terms of the slope variables o(7)
alone. Dropping a sand grain corresponds to a distortion in the slope, but the sum
of the local slopes over the system remains unchanged. The dynamics when the

threshold o, is exceeded becomes

o(i) = o(i)—2n, (2.18)
o(it1) = ofitl)+n, (2.19)

which is very much like a discretized diffusion equation.

As our model system we will now take the equation:

Ay(x,t) = Dé?y(z,t), (2.20)
D = 1 fory>0

l—e fory<0.

We interpret y as being similar to the slope variable in the local limited model,
and we have implemented the nonlinear threshold condition by the discontinuous
diffusion coefficient D. We will be interested in the decay of initial conditions
corresponding to a local distortion of y, e.g., taking the form of a single period of
a sine wave. It is clear that Eq. (2.20) is nonlinear for « # 0, and one may even
question the existence of a solution. A theorem by Kamenomostskaya [15], however,

guarantees the existence of solutions which are both continuous and have continuous

72

first derivatives.

We should emphasize that Eq. (2.20) cannot be considered as a good description
of the local limited model, but is rather an attempt to capture features which may
be relevant in similar models. For example, it can be shown that avalanches in the:
local limited model always terminate at ‘troughs’, which are sites where the local
slope o(i) satisfies o(i) < 0. — n [12]. Numerically, it is found that the density of
such troughs scales with the system size as L~'/3, The trough dynamics are clearly
not captured in Eq. (2.20). Some of the consequences of the trough dynamics are

analyzed in [16].

2.2.3. The Barenblatt Equation

Barenblatt [6] has studied an equation very similar to Eq. (2.20), and we will there-
fore refer to Eq. (2.20) and similar piecewise linear diffusion equations as Barenblatt
equations. Barenblatt considered the case where the diffusion coefficient D is a dis-
continuous function of 0,7y, rather than y itself, and where the initial condition is
a simple pulse rather than the sine pulse considered here. The Barenblatt equation
has also been used to illustrate singular perturbation theory, if one tries to expand
about € =0 [17].

Let us first consider the trivial case € = 0. As indicated in the previous section,
we expect scaling solutions of the form t~7/?f(x/./t). This is most easily seen from

the exact solution for arbitrary initial conditions and zero slope at the origin:

west) = 7 [* deule){exp(—(€ + 2/41) + exp(—(€ ~ 2)°/48)} (2.21)

If u(x) falls off fast enough as + — oo, we can obtain the asymptotic solution for
large t by expanding the exponentials. Requiring the spatial integral of u to vanish,

we directly obtain the desired scaling form.

73

Consider now the case e # 0. Applying dimensional analysis naively, we would

expect a scaling solution of the form:
y(a,t) = 1! f(x/ Vi, €). (2.23)

Substituting this ansatz into Eq. (2.20) we obtain an ordinary differential equation

for f in terms of the single variable n = 2/ Vt
it 1 Ul 3
Df + gif + af = 0. (2.24)

We now need to separately solve the cases f < 0 and f > 0 and match the solutions
at f = 0. The solution to Eq. (2.24) which is positive and has zero slope at the

origin can be expressed in terms of a parabolic cylinder function U [18]:
f(n) = —Aexp(—7?/8)U(-5/2, n/V2), (2.25)

where A is a constant. Similarly, the negative solution to Eq. (2.24) which goes to

ZeTO as 7 — 00 Is:

F(a) = Bexo(-s"™ = U(-5/2, nf 209), (2.26)

where B is a constant. To match the solutions we need to locate a common zero of
Eq. (2.25) and Eq. (2.26). However, since U(—5/2, x) has only one zero [18], there
is no common solution unless € = 0. We conclude that there are no scaling solutions
of the form assumed in Eq. (2.23).

Surprisingly, a scaling solution is still possible. Indeed, we need to reexamine
our assumption that we can neglect the spatial extent | of the initial condition. In
our previous analysis, we have assumed that the scaling function f(x/Vt, 1/4, €)
will become independent of / in the limit of large t and small J. Another possibility,

however, is that f is a singular function of I, taking on the form f ~ (I//t)~*,

74

with a being an arbitrary exponent which cannot be determined from dimensional

analysis alone. We then write our new scaling ansatz as:
yl, t) = 720+ F(a/Vt). (2.27)
The differential equation for f then becomes:
i 1 7 3
f+ gif + 51 +a)f =0. (2.28)

We now follow the same procedure of matching solutions for f > 0 and f < 0. The
solution for f > 0 having zero slope at the origin can be expressed in terms of a

confluent hypergeometric function M [18]:

f(n) = M(3(1 + @)/2, 1/2, —n?/4), (2.29)

where we have fixed an arbitrary overall normalization by setting f(0) = 1. The

solution for f < 0 which goes to zero as n — oo is now:

Fn) = Aexp(— gr )U(-80 — 5/2, 0/201 6). (2.30)

Solving f(7) = 0 simultaneously for both Eq. (2.29) and Eq. (2.30) then leads to

the following nonlinear eigenvalue problem:

M(3(1 + a)/2, 1/2, —n§/4) = 0 (2.31)
U(—3a — 5/2, m/Y2(1—-e)) = 0. (2.32)

Once a and 7 are determined for a given €, the constant A in Eq. (2.30) is fixed by
requiring continuity of the first derivatives of f.
In Figure 2-1 we plot the numerical solutions to Eq. (2.31) and Eq. (2.32).

The solutions for 0 < € < 1 are unique, and a increases monotonically with e. To

75

1 TOTP T TOT TT TT Tt rs ee i TOF TT

0.8

0.6

0.4

0.2

rere tp ye Py pt tt

(0) J I I. 4 i i | i 1 1 1 | i 1 1 i | J Leek ra
0 0.2 0.4 0.6 0.8

Figure 2-1: Solution to nonlinear eigenvalue problem.

rigorously prove that the scaling solution so obtained is an attractor for arbitrary
initial conditions appears very difficult. We limit ourselves to a simple numerical
simulation, where we integrate the nonlinear equations for « = 0.5 and the initial
condition y(x,0) = 1 for 0 < x <1, and y(z,0) = —-1 for 1 < x < 2. In Figure
2-2 we plot y(0,t) vs. t, which clearly shows a power law decay from which we
deduce a@ = 0.1320.002, in excellent agreement with a = 0.1324 determined from
Eq. (2.31).

We have also verified that the solution y(x,t) rapidly approaches the scaling
form Eq. (2.27), for several initial conditions.

To summarize, we have found scaling solutions to a simple nonlinear diffusion
equation, with a non-trivial exponent determined by solving a nonlinear eigenvalue
problem. As in the linear diffusion case, the scaling form corresponds to the solution

for a particular initial condition, which acts as an attractor for a large class of other

76

on?
1 t ToT erty T T TTT 7 TOT TPTit T TOP TTtt

th

y(0,t)

4 peu 1 coil L Lieu! it ctu i oul Pivilt

eS
Oo

T TUTTI T TTT T TTT T TUTTI TPT 7 TUT T TTT

10-9 f i rootul i 1 rit it i roonl 1 tT
107! 10° 10! 10° 10°

Figure 2-2: Anomalous exponent in nonlinear diffusion equation.

initial conditions. In the linear case, this initial condition is simply a delta function,
for the nonlinear case, this initial condition is more complicated. It should be
emphasized that the conservation law was not important in arriving at the scaling

form.

2.3 Sandpile Models

2.3.1 Background

Sandpile models were introduced by Bak, Tang, and Wiesenfeld [5] who believe
them to be simple model systems that display self-organized criticality. The original
sandpile models are defined as d-dimensional hypercubic lattices of size L¢, where
each lattice site z;,i = 1,...,L%, takes on a discrete set of positive integral values.

The integer z; is referred to as the number of ‘sand grains’ at site 2. When z; < 2d,

77

the site 2 is called stable. The dynamics is defined by the rules:

1. When the configuration is stable, i.e., all z; < 2d, a site 7 is selected at random

and z; is incremented, z; — z;+1.

2. If the configuration is unstable, then synchronously for each unstable site z;

we apply the rule:

7 on a) 2d,

zy 2341

where 7 ranges over the nearest neighbors of 7. We will refer to the application of
Rule 2 to a given site 2 as a ‘toppling’ event. The boundary conditions are typically
chosen such that if 7 falls outside the lattice, the corresponding sand grain is ‘lost’.
Rule 2 is applied until all sites are stable, after which a sand grain is added again
randomly, using Rule 1. Because sand grains can ‘fall off’ at the boundary, the
number of grains in the lattice will remain bounded.

The event of adding a sand grain and relaxing the lattice to equilibrium will
be referred to as an ‘avalanche’. For each avalanche, the following quantities can
be defined. The ‘flip number’ F is the total number of toppling events during
an avalanche. The ‘size’ S is the number of distinct sites that topple during an
avalanche. Because a given site can topple more than once during a single avalanche,
we evidently have S < F. Finally we define the ‘duration’ T of an avalanche as the
number of synchronous updates required by Rule 2 to relax the lattice to a stable
equilibrium configuration.

Extensive computer simulations by numerous groups [19] have clearly established
that the distribution functions D(F’), D(S), and D(T) have power law forms for large
system sizes L. In Figure 2-3 we show D(F) for two-dimensional lattices of several
system sizes. Defining the ‘toppling exponent’ 7 by the relation D(F) ~ F!~7, we
find 7 = 2.124 0.05, in agreement with other investigations [20].

78

~-1
10 TPT TOT TTn TOT TTTnn TOT TETiT TOT UTINT TPT EIT

2D Conserved Sanpile

D(F)

1071! L Leuul 1 cituul I fwannit wt Lil i Ltt! Loeb LE
10' 10% 10° 10% 10° 108 ~ § 10%

Figure 2-3: Distribution of flip number for 2D conserved sandpile model.

Various modifications and generalizations of the sandpile models have been con-
sidered in the literature and we give only a brief summary of those which, we believe,
lead to the most important quantitative differences. Among the sandpile models,
there seems to be clear differentiation between those operating in one spatial di-
mension, and those operating in higher dimensions. In fact, the sandpile model
described above leads, in one dimension, to a trivial critical state where all sites
assume value one. Nontrivial behavior in one-dimensional models is only observed
for more complicated rules, such as the local limited model described in section
2.2.2. These one-dimensional models generally do not show simple power laws, but
more complicated distributions, which may be multi-fractal [13, 14]. For a particu-
lar one-dimensional model it has been shown that the dynamics can be mapped to
a nonlinear diffusion equation, where the diffusion coefficient diverges as a power of

the system size [21].

79

As has been pointed out by Dhar [22] the sandpile model defined above has an
Abelian property, in that the final equilibrium configuration reached after a set of
perturbations does not depend on the order in which the perturbations are applied.
Using the Abelian property, it can be shown [22] that after an initial transient, the
configuration space of the lattice consists of a set of recurrent configurations, which
all occur with equal probability. The number of recurrent configurations diverges
exponentially with the number of lattice sites N, but is somewhat less than the total
number of stable configurations (2d)". Dhar has also shown that the expectation
value of the flip number F obeys the relation (F) ~ L? in arbitrary dimensions [22].
However, even with the Abelian property no explicit expressions for the exponents
have been obtained.

The Abelian property is easily destroyed, for example, by toppling rules that
depend on the local slope, rather than the local height. Examples of such models
have been considered in [13] and [23]; such models usually display power-law behav-
ior as well, but with different exponents. In sections 2.3.2 and 2.3.3 we will study a
sandpile model which violates the local conservation law and is not Abelian.

Experimental work on real sandpiles [24, 25] usually shows dynamic behavior
which is too complex to be described by simple power laws. Large fluctuations in
avalanche sizes have been observed under certain conditions [26] but the precise
nature of those conditions and how they may relate to the sandpile models has not
been elucidated.

There have been a number of theoretical approaches to the sandpile models,
but apart from Dhar’s work on Abelian sandpiles [22], none have led to conclusive
results. One of the difficulties is that there seems to be little universality in the
exponents observed for different models. Another problem is that the concept of
an avalanche is difficult to define within the context of a continuum field theory.
To illustrate these difficulties, we briefly review the work by Hwa and Kardar [27],

which has motivated extensive research using similar approaches [28, 29, 30].

80

Hwa and Kardar [27] construct a continuum model of a sandpile that incorporates
several symmetry considerations. Sand flows in a preferred ‘downhill’ direction,
while the system is translationally invariant in the transverse direction. Sand is
locally conserved, except for the external driving, which is assumed to take the form
of an uncorrelated Gaussian noise source. Using these considerations, Hwa and
Kardar [27] argue that the simplest equation for the fluctuations h of the slope of

the sandpile about is critical slope is:
2 2 r- 19

where the first terms describe relaxation of the height through surface tension, split
into two parts perpendicular and parallel to the transport direction. 7 is a Gaussian

noise term with zero mean and second moment:
(n(x, t)n(x’, t')) = 2D6%(x — x')6(¢ — 1), (2.34)

where D is the noise amplitude. In the long wavelength limit, the two-point corre-

lation function C(x — x’,t — t’) = (h(x, t)h(x’,t’)) is assumed to take the form:

t 2
C(x, t) = eaers =) (2.35)

where €, z, and ¢ are the roughness, dynamic, and anisotropy exponents respectively.
Using the dynamic renormalization group [31], the exponents for Eq. (2.33) can in
fact be determined exactly. Above a critical dimension of d = 4, the nonlinear
term in Eq. (2.33) becomes irrelevant and ordinary diffusive behavior with z = 2 is
obtained.

Note that the symmetries incorporated into Eq. (2.33) are quite different from
those of the original sandpile model. In particular, the original model does not

have a preferred direction. Nevertheless, an exact solution for a particular directed

81

sandpile model is known [32], with predictions that deviate considerably from those
obtained from Eq. (2.33).

There have been several attempts to modify Eq. (2.33) to incorporate the correct
symmetries of the various sandpile models [28, 29] but none have led to reliable
qualitative predictions. Another problem with trying to apply Eq. (2.33) is that the
effective noise term in the original sandpile models may have complicated temporal
correlations. The effect of these correlations has indeed been considered [33] but
again no definite results have been established.

We conclude that while noisy diffusion equations of the form Eq. (2.33) are
interesting models in their own right, they may not provide reliable information
about the sandpile models. In particular, the question of which invariance principles,
if any, are at work in the sandpile models probably needs to be approached from a
different point of view. We will present such a view in section 2.3.2 and discuss our

numerical work in section 2.3.3.

2.3.2 Mean Field Theory

For complicated many-body problems, useful insights can often be obtained from
mean field theories, in which correlations are ignored. In the context of second order
phase transitions mean field theories usually give correct results above a critical
dimension d,. The mean field limit of noisy diffusion equations like Eq. (2.33)
corresponds to neglecting the nonlinear terms.

Tang and Bak [34] have obtained a mean field theory for sandpiles by assuming
that neighboring sites in the lattice are independent. This leads to an evolution
equation of the average lattice amplitude, including a random term due to the
driving. The validity of the approach by Tang and Bak has been questioned in
[35], where the importance of edge effects is emphasized, since sand can only be
lost through the edges. Furthermore, the approach in [34] does not lead to a direct

prediction of the exponents.

82

In our mean field theory we will assume that within a given avalanche, each top-
pling event leads to its own independent avalanche. This is equivalent to assuming
that during any toppling, sand grains are not dropped to the nearest neighbors, but
to sites arbitrarily far away. Mathematically, this leads to the well-known subject
of branching processes [36], originally developed to study the fate of family trees in
the British aristocracy [37].

In the following discussion, we will refer to the grains that are being dropped
due to topplings as children. Let Zo, Z,,Z2,... denote the number of children at
successive generations, where each generation corresponds to a synchronous update
of the lattice. We have Zp = 1, since the first child is simply the initial grain added to
the lattice. Let each toppling event generate Z children and let p be the probability
that a child causes another toppling event. By our assumption of independence, p
is the same for each child at every generation. The average number of children m
that a single child will have is then obviously given by m = pZ. If m < 1, then each
child will have less than one child on average, and it can be shown rigorously [36]
that the avalanche will die out in a finite number of steps with probability one. On
the other hand, if m > 1, there is a finite probability that the avalanche will last
forever, i.e., none of the Z, is zero for arbitrarily large n.

We now suggest that the borderline situation of m = 1 is precisely the ‘critical
point’ of the avalanche dynamics. This critical point should be quite independent
of how the avalanche actually takes place, in particular, we do not expect local
conservation laws to play a significant role. The only assumption is that the system
is homogeneous, so that the probability p is indeed the same for each child.

For the special case m = 1, several results have been obtained for branching pro-
cesses, which we can use to determine the flip number exponent 7, and the duration
exponent y, defined by the relation D(T) ~ T~¥. The flip number distribution D(F)
is evidently just the distribution of F = Z) + Z; +--+. It was shown by Otter [38]

83

that:
D(F=r)~r 3, ro, (2.36)

from which we immediately deduce 7 = 3/2+ 1 = 5/2. Finally, it was shown by
Kolmogorov [39] that the extinction probability P(Z, = 0) ~ 1/n for large n. Since
the extinction probability is the cumulative distribution function of the relaxation
times, we deduce that y = 2. These predictions agree with exact solutions of the
sandpile model on a Bethe lattice [40]. In section 2.3.3 below we present numeri-
cal evidence that a sandpile without conservation law or other obvious invariance
principles shows criticality, and we compare its behavior with the corresponding

conserved model.

2.3.3 Numerical Simulations

The non-conserved model we will investigate has almost the same equations of mo-

tion as the conserved model described in section 2.3.1. The rules are [41]:

1. When the configuration is stable, i.e., all z; < 2d, a site 7 is selected at random

and z; is incremented, z; — z; +1.

2. If the configuration is unstable

Zz 0

25 _— zj +1.

Evidently only Rule 2 is different from the conserved model. When a toppling site
has more than 2d grains, z; — 2d grains are lost from the system.

In Figure 2-4 we show the flip number distribution for two-dimensional lattices.
We clearly observe a power law, from which we deduce 7 ¥ 2.44 + 0.02. While this
result is very close to the value 7 = 2.5 predicted by mean field theory, we find

numerically that the exponent y © 1.68 + 0.05 (cf. Figure 2-5), quite different from

84

10 TT TTT Trinny TTT TT ETITIT T TTI

10° .
2D Non—Conserved Sandpile
107°
10~*

io?

D(F)

1078
1077

1078

107°

L=1000

ii TOT T Try T TTT T TUTTE TTT

1 Lossutl 1 toil L Loutul i rid tL LEV Ii
10° 10° 104 10° 108

107°
10

ay

Figure 2-4: Distribution of flip number for 2D non-conserved sandpile model.

the mean field value y = 2.

The exponents of the non-conserved model are thus quite different from those of
the conserved model, which has r + 2.12 and y = 1.18 [19].

Comparing Figures 2-3 and 2-4 it is also clear that considerably larger system
sizes are needed in the non-conserved model to observe criticality. We can trace some
of this behavior to the differences in space-time dynamics between the conserved and
non-conserved models. In Figure 2-6 we show the duration T of an avalanche as a
function of the linear extent D of the avalanche. We find that T ~ DS with ¢ + 1.
This suggests that avalanches in the non-conserved model are very anisotropic. For
computer simulations, this result indicates that finite size effects will be much more
pronounced than in a system whose behavior is closer to ordinary diffusion with
¢ = 2. For the conserved model, one finds ¢ = 1.5 [42].

As with the theory of critical phenomena, we expect sandpile models to approach

the mean field limit as the spatial dimension is increased, because avalanches will

85

107} TT TRACT

T ToT TTT TT TT rTtt TUT TTTty.

T TPT
Pip tititn

2D Non-Conserved

10°* Sandpile

T UTTTTATT
1 roti

io3

T TP TOT
ii bail

D{T)

1074

T UPerrayy
4 resriul

107 L=3000

1 ciasut

L=1000

t |

1078 i I. retirl J J. Loitiul i i Littul i Joi ft tt
10° 10! 10” 103 10+

Figure 2-5: Distribution of relaxation times for 2D non-conserved sandpile model.

3000 T T F T T T T T ¥ T T F
I 2L Sandpile |
2000 f- _
7 i: 4
00 |-
10 : L=3000 |
0 t i i 1 | i i rt i | 1 i i i
0 500 1000 1500

Figure 2-6: Duration of avalanches as a function of linear extent for 2D non-
conserved sandpile model.

86

3D Non—Conserved Sandpile

D(F)

L=250

19719 4 rreiul l rool 4 weet 4 col 1 boul ioe eeenl
10° 107 10% 10% 10% 10° 108

Figure 2-7. Distribution of flip number for 3D non-conserved sandpile model.

be ‘less interacting’ due to the number of different directions available to them. In
critical phenomena, the mean field limit is usually reached at finite upper dimension
d,, which is typically four or six. Whether there exists an upper critical dimension
for sandpile models is not clear at present. For the conserved model, numerical
results for d = 4 [19] are quite close to the mean field exponents derived here. It
has been claimed that the exponents are different for d = 5 [42], but only very
small systems could be simulated casting some doubt on these predictions. Our
own investigations are hampered by finite size effects as well, and even for d = 3
we cannot obtain very accurate results. In Figure 2-7 we show the flip number
distribution for three-dimensional lattices. It is clearly very difficult to obtain a
precise exponent from these data: we estimate 7 + 1.50+0.20. We therefore do not
have enough numerical resolution to determine whether in large enough dimensions

the mean field exponents are obtained.

87

To summarize, we have shown that critical fluctuations can occur in models
that do not have an obvious invariance principle in their equations of motion. We
have presented a very simple mean field theory, which predicts exponents roughly
consistent with those seen in computer simulations of conserved models of sandpiles
in four dimensions. For the non-conserved sandpile model, the exponents obtained
in two dimensions are very close, but not quite equal to the mean field predictions.
We do not have sufficient computer power to obtain conclusive results for our non-

conserved model in higher dimensions.

2.4 The Game of Life

2.4.1 Cellular Automata

Cellular automata were introduced by John von Neumann [43] as simplified compu-
tational models for biological systems. The main defining characteristic of a cellular
automaton model is that it consists of a discrete lattice of cells, where each cell
can only take on a discrete set of values [44]. An equation of motion is defined
in discrete time, where all lattice sites are updated synchronously according to a
specified rule. This rule is typically deterministic and described by Boolean algebra,
but probabilistic cellular automata, for which the rules are stochastic, can also be
considered. The rule is generally local in space and time, i.e., for any given time
step, each site is only affected by a small number of neighboring sites, and only a
small number of previous time steps.

It is clear from their definition that cellular automata are mainly designed for
efficient simulation on digital computers, and at first sight it is not obvious why they
should provide useful models for physical systems. However, many ‘traditional’
models are already very similar to cellular automata. A simple example is the
Ising model, which, since its introduction in 1925 [45], has been one of the most-

studied models for many body phenomena, with over 500 research papers involving

88

the Ising model published in 1991 alone. The Ising model and its generalization,
the Potts model [46], already satisfy most of the criteria for cellular automata,
except that the original models were purely static. Dynamic generalizations of
the Ising model typically involve continuous time, rather than discrete time [47,
48]. However, there are many direct relations between Ising models and cellular
automata. For example, Creutz [49] has shown that it is possible to construct purely
deterministic cellular automata whose dynamics generate the same distribution as
equivalent static Ising models. These deterministic cellular automata allow much
faster computer simulations of Ising models than more conventional Monte Carlo
algorithms [50]. It is also well known that cellular automata models in d spatial
dimensions can be mapped to Ising models in (d + 1) dimensions [51].

The close connection between Ising models and cellular automata is perhaps
not too surprising. More interesting is the observation by Frisch, Hasslacher, and
Pomeau [52] that it is possible to construct cellular automata which, in the limit of
large system size, can simulate the Navier-Stokes equation. This result is not trivial,
since it is not immediately obvious how to incorporate the correct conservation laws
and symmetries, like Galilean invariance, into a cellular automaton. The massively
parallel algorithms available for cellular automata might therefore provide the com-
putational capability necessary to simulate fluid flow at high Reynolds number. A
potential disadvantage of cellular automata is that coarse graining over many lattice
sites may be necessary to obtain a physically relevant picture [53].

Another advantage of cellular automaton models is that they allow simulations
of collective behavior in large spatial dimensions. This allows numerical verification
of predictions that above a certain critical dimension many body systems display
trivial ‘mean-field’ behavior. Evidence for nontrivial collective behavior of cellular
automata in five spatial dimensions has been presented in [54], where it is argued
that this behavior may require new theoretical tools for its description. However,

only numerical evidence is available at this time, and the argument must be regarded

89

as somewhat speculative.

Finally, the study of cellular automata for their own sake can be very instruc-
tive. The behavior of any given cellular automaton can be classified roughly into
four categories [44]. The first category corresponds to evolution to a spatially homo-
geneous state, analogous to a fixed point in dynamical systems. The second category
evolves to spatially periodic patterns, analogous to limit cycles in dynamical sys-
tems. The third class shows chaotic evolution, characterized by positive Lyapunov
exponents [55]. The fourth class is the most complicated — cellular automata in this
class are capable of universal computation. A computationally universal system is
one that can be programmed, through its initial conditions, to simulate any digital
computation. Computational universality has also been demonstrated for classical
mechanical systems consisting of a set of billiard balls [56]. Computationally univer-
sal systems are necessarily nonergodic [57], since the indefinite future depends in an
arbitrarily programmable way on the initial conditions. Computational universality
is of course only possible in an infinite system: spatially finite systems can only
evaluate ‘spatially-bounded’ functions [58]. We therefore expect computationally
universal systems not to show scale invariance in general, since arbitrarily complex
functions may become computable as the system size increases.

Although all computationally universal systems are identical in their capabilities,
there may still be considerable qualitative differences between them if regarded as
dynamical systems [59]. For example, suppose we have a two-dimensional computa-
tionally universal cellular automaton with two possible states per site, one or zero.
It may happen that nontrivial computations only result from initial configurations
where every other site has value zero, in checkerboard fashion. For random initial
conditions, these configurations occur with probability zero, so that this particular
automaton cannot serve as a model for ‘self-organization’ in nature.

The Game of Life [7| is a two-dimensional cellular automaton with two states

per site that has been shown to be computationally universal [60]. The proof was

90

accomplished by explicitly constructing ‘wires’ out of structures which propagate
without dissipation, and ‘logical gates’ that can perform the Boolean NAND opera-
tion [60] (see also the discussion of Boolean algebra in the Appendix). The Game of
Life is therefore an extremely complex system, and computer simulations on finite
lattices may not provide results that have statistical significance for large systems.
Nevertheless, using arguments based on computer simulations on square lattices up
to size 100 x 100, Bak et al. [8] maintain that the Game of Life naturally evolves
to a self-organized critical state, characterized by power laws in various distribution
functions. We will discuss the Game of Life in more detail in sections 2.4.2 and

2.4.3, where we present numerical evidence that the Game of Life is not critical.

2.4.2 Rules of the Game

The rules of the Game of Life are deceptively simple. The dynamics takes place
on a square lattice, where each site can be either ‘dead’ or ‘alive’. At each time
step, every site is only affected by its eight nearest and next nearest neighbors. In

particular
e A live site will stay alive if it has two or three live neighbors.
e A dead site will become alive if it has exactly three live neighbors.

The Game of Life can thus be considered as a simple model for a biological system,
where excessive births are counterbalanced by deaths through overcrowding, and
very small populations cannot survive.

The Game of Life has been one of the most popular computational pastimes,
and a considerable amount of empirical data have been obtained from computer
simulations [61]. When the initial density of live sites on an N x N lattice is
between 0.2 and 0.7, the lattice evolves in ~ N? steps to a configuration containing
a fraction of © 0.03 live sites. This configuration generally consists of stable patterns,

or patterns oscillating with a short temporal period. In Figure 2-8 we show a typical

91

section of a lattice that has reached equilibrium.

The dark circles show live sites, the hollow circles show dead sites that will
become alive at the next iteration. The largest temporal period in Figure 2-8 is
two. The most striking feature of Figure 2-8 is the clustering of live sites into stable
configurations containing only a small number of sites, separated by relatively large
distances. This is already an indication that a simple mean field theory may not
be capable of predicting the asymptotic density. Possibly the simplest mean field
theory consists of assigning a live probability p" to each site at time n, and then
computing the probability p”*! that a site will be alive at the next time step [62].
We then obtain:

p"** = 28(p")9(1 — p")°(3 — p"). (2.37)

This equation correctly predicts a stable fixed point for p = 0, but it also predicts
another stable fixed point for p = 0.37, which is never observed in simulations.
We can improve the simple mean field theory by computing p”*? from p”. This
requires the complete solution of the Game of Life in a 5 x 5 neighborhood, since
n+2

the interaction propagates one lattice spacing in each time step. We can write p

in the form:
25
n+2

pr? = 2 ae(P") (1 — pre, (2.38)
i=
where a; is the number of configurations containing 7 live sites for which the central
site is alive after two time steps. We show the resulting return map in Figure 2-9;
the only fixed point remaining in this approximation is p = 0.

Higher order return maps clearly become impractical to evaluate exactly, and
one needs to resort to Monte Carlo simulations to arrive at the asymptotic density
of = 0.03 [61].

So far we have ignored the presence of propagating structures, which are the

main reasons the Game of Life can be computationally universal. By far the most

common propagating structure is the glider (cf. Figure 2-10), which moves by one

92

dd
se
eee
fe)
eee
fe)
3 e
ee
ee
oe0 ee
e ee
oe
ee

ese
e ° @
o80 0e0
e e
ese
°, 8
Se
ee
ee
ee
ee
ee
ee
° °
o@0

ee

ee

oe0

Figure 2-8: Typical Game of Life equilibrium configuration.

93

1 ToT TT TOT tT TT TT TT TF T T TTT

0.8

0.6

n+2
Pp

0.4

0.2

PeTeT OT TTT TT TT PT TT

[an ee | a oe oe 2 ! | | i 4 | | ee oe | Joo}

0 1 | ot | Ee ne OO | | i o| i i | ee ae
0.2 0.4 0.6 0.8

oO
eR

Figure 2-9: Return map for Game of Life.

diagonal square in four time steps.

Gliders are generated at a rate of approximately one per 2000 sites. Their mean
free path is = 25 lattice spacings [63] in a lattice that has otherwise reached equi-
librium. There are particular initial configurations called glider guns which can
continuously generate gliders [61], so that the population in an infinite system can
increase indefinitely. When we speak of equilibrium in a finite lattice, we therefore

mean configurations where propagating structures have either escaped to infinity,

if |
TI itt TT rit

Figure 2-10: The glider.

94

or have been halted at the boundary.

Bak et al. [8] have presented numerical evidence that if an equilibrium Life
configuration is repeatedly perturbed, the distribution of the relaxation times D(T)
before equilibrium is reestablished follows a power law D(T) ~ T-1%. This would
suggest that the equilibrium configuration rearranges itself to a minimally stable
state, where perturbations can propagate an arbitrary distance. In the next section,
we will present numerical evidence that the distribution D(T) is in fact exponential
of the form D(T) ~ exp(—T'/To), where To is a characteristic time scale of an

avalanche [64].

2.4.3 Numerical Simulations

In any numerical investigation of scaling behavior it is important to be aware of the
relevant scales involved. The most obvious scales in a finite lattice are the lattice
spacing on the one hand, and the lattice size L on the other hand. Identifying
intermediate scales will depend on the particular system in question. From Figure 2-
8 it appears that stable clusters in the Game of Life are separated by 10 or more
lattice spacings, so that we may expect intermediate length scales of that order.
The obvious time scale for the Game of Life is the time it takes for a random
configuration to relax to equilibrium; this time is always finite in a finite system.
If the driven system indeed displays self-organized criticality, it must do so over a
time scale larger than the relaxation time of the undriven system.

In Figure 2-11 we show the distribution of relaxation times for two lattices size
256 x 256 and 512 x 512, starting from a random initial configuration, without
external driving.

Asymptotically we find the distribution function to take the form D(T) ~
exp(—T/To), where To * 1850 independent of the lattice size. We conclude that
To is an important time scale for the undriven system, and to convincingly demon-

strate scaling behavior for the driven system we need to accumulate data over several

95

10°73
TT T T T T T T T T T T T T T T T T T T

T TT TTT

in oe a i

1074

T ToT TTITTy

L boil

107°

D(T)
TTT TTT
1 Jot ri

107°

T TTT
i reo

tov? Leet

0 5000 10000 15000 20000 25000

Figure 2-11: Life relaxation times, undriven system.

factors of Tp.

In Figure 2-12 we present our results for D(T) for the driven system, where each
relaxation time T is the result of a single perturbation to a system which has relaxed
to equilibrium from a previous perturbation.

We again find an exponential distribution for large T’, with the same time con-
stant Ty as the undriven system. The data for lattice sizes 256 x 256 and 512 x 512
are virtually identical, indicating that the finite relaxation time is indeed the result
of a finite correlation length, rather than a finite size effect.

In addition to the relaxation time To, it is also useful to measure the correlation
length Lo to compare with the system size L. We can obtain a rough estimate of
Lo by measuring the net number A of sites that are changed during an avalanche,
which should be approximately proportional to the area covered by the avalanche.

In Figure 2-13 we show that A is roughly proportional to the duration T indicating

96

107?

1073

1074

107°

10-8

1077

1078

1 cout 1 Liu L Lisl 4 iui i Lieut | uel Lot titi:

T TUT T PUTT TTT

107? { L L i | 1 i i | 1 1 1 i | | i L 1

5000 10000 15000 20000

con)

Figure 2-12: Life relaxation times, driven system.

that avalanches are propagating similar to an ordinary diffusive random walk. Note
that the y-intercept is not zero because even at equilibrium there is a background
activity on the lattice due to local oscillators, which we have not subtracted out.

To obtain an estimate of the actual spatial size of an avalanche, we use our
knowledge that the equilibrium density is ~ 0.03. We find that an avalanche with
duration To + 1850 covers an area of = 8000 sites, leading to an estimate for the
correlation length Lo = 50.

To determine Lp more precisely, we perturb the lattice not at a single point at
a time, but along a whole edge, which we call a ‘hot’ boundary. We then measure
the active and passive densities as a function of the distance away from the hot
boundary (the active density is defined by live sites not alive six steps ago). We
show our results in Figure 2-14, for a lattice size 256 x 488 with the hot boundary

on the left, a closed boundary on the right, and periodic boundary conditions in the

97

7000 T T F T T T T T T T ¥ T T 7 r

6000

5000

4000

3000

COT TTT TT PTT TTT TTT

poe tee ee Pa fe |

2000 i i 1 i | i 1 1 1 | i i 1 ik | i i t ik
0 5000 10000 15000 20000
Time

Figure 2-13: Number of sites changed during avalanche.

vertical direction.

The passive density approaches the characteristic equilibrium value ~ 0.03 away
from the hot boundary, while the active density decreases exponentially with a char-
acteristic length Lp = 42 +3. We obtained equivalent results for lattices of vertical
size 128, thereby verifying that our results are not affected by vertical correlations.

To summarize, we have obtained strong numerical evidence that the Game of
Life is subcritical with a large but finite relaxation time Ty ~ 1850, independent
of whether the system relaxes to equilibrium from random initial conditions, or
from repeated perturbations at isolated points. We estimate a correlation length
Io © 42 +3. This suggests that the numerical results presented by Bak et al. [8]

for system size 100 x 100 are strongly influenced by finite size effects.

98

ttt

Passive

POUT
veya

TT THT

) 4 vin

Density

i) uuu

1074

TTT TT |

ij sist

1075

T TTT)

pope Eee tt

0 100 200 300 400
Distance

Figure 2-14: Active and passive densities away from hot boundary.

2.4.4 A Universe Model

We emphasize that the numerical results presented in section 2.4.3 do not exclude
the possibility that the Game of Life may be metastable with respect to nucleation
events too rare to have occurred in any computer simulation to date. However, our
simulations strongly suggest that these metastable states, if they indeed exist, do
not occur ‘naturally’ as a result of repeated driving, as required by self-organized
criticality.

In this context we briefly describe another cellular automaton, called the universe
model [65], which we believe illustrates another way one might erroneously interpret
numerical simulations as showing self-organized criticality. In contrast with the
Game of Life, the universe model does not display a definite correlation length
independent of the system size. Instead, the distribution of relaxation times decays

exponentially with a time constant that increases with system size. At first sight

99

oT

Figure 2-15: Glider in universe model.

this suggests that the model is indeed critical, since the system size is the only
relevant length scale. However, we will argue that in this case the relaxation times
result from relaxations of states which are only stabilized by finite size effects. Such
states do not occur in the infinite size limit, and it is therefore meaningless to talk
about critical fluctuations around these states.

The universe model is defined on any d-dimensional hypercubic lattice, with

three possible states per site, ‘active’, ‘passive’, and ‘empty’. The rules are:
1. Active sites at time ¢t ‘burn out’ and become passive at time t + 1.

2. Passive sites are annihilated and become empty when they have exactly one

active neighbor.

3. Empty sites become active when the have exactly one active neighbor which

has a passive site at the opposite position.

Rule 3 allows the existence of ‘gliders’ as in the Game of Life. The configuration
shown in Figure 2-15 propagates to the right with unit velocity, where the filled
circle represents an active site, while the empty circle represents a passive site (in
the figure, all other neighbors are empty sites).

In their simulations, Chen et al. [65] start with a random initial distribution
of active, passive, and empty sites and let the system come to rest with no active
sites remaining. As in the Game of Life, the system is then repeatedly perturbed
by adding a single propagating particle, and the relaxation times to return to rest
are measured. As this process is repeated, the system is believed to evolve to a
statistically stationary state, with power law distributions in the relaxation times

[65].

100

In our simulations, we continually start with a random initial configuration,
rather than repeatedly perturb one that has come to rest. We find that the dis-
tribution D(T) of relaxation times has the form D(T) ~ T exp(—T'/To), where Ty
increases rapidly with system size. In Figures 2-16a and 2-16b we show D(T) for
system sizes 60 x 60 and 92 x 92. In Figure 2-16c we plot the relaxation time To for
system sizes ranging from 28 x 28 to 156 x 156.

Figure 2-16c shows that To increases approximately exponentially with system
size. This result suggests that configurations with no active sites will have negligible
probability in the large system limit. By extensive simulations on a system of size
384 x 384 we estimate that the steady state density of active sites is + 0.0088.

To summarize, we have shown that finite size effects are probably responsible
for the apparent scale invariance observed in two cellular automaton models. In
this context we mention that numerical evidence has been presented [66] that yet
another cellular automaton, called the ‘forest fire’ model, does not seem to show

nontrivial scaling behavior as had previously been suggested [67].

2.5 Conclusion

Our investigations of self-organized criticality have included two broad classes of
models. The first class contains the sandpile models, for which we have presented
a simple mean field theory which predicts power law behavior under very general
conditions. The second class contains cellular automata which can be considered as
crude models of biological systems. We have presented numerical evidence that the
scale invariance reported for two such models was due to finite size effects and does
not appear to hold for larger systems.

We have also presented a third scenario under which nontrivial scaling behavior
can occur, by studying a variant of the Barenblatt equation. The scaling behavior

is particularly complicated because the scaling exponent depends continuously on a

101

- ~6
ed 10 rTP TTT ryt

“5
10 L=60 L=92

1078

D(T)/7
D(T)/T

1077

| TTT I TTTTea TT

107°
(a) (b)
to? LL tt bret yutirey to@ LLL It re
0 500 1000 1500 0 2000 4000 6000
T T
10° E774 PTT TT] Trt gt tt 1S
F (ce) =
104 _
E é q
ee 10° —
107 : =
[lg 4
io! LLL LI ee
0 50 100 150 200

Figure 2-16: Universe model (a) D(T) for L = 60 (b) D(T) for L = 92 (c) Tp versus
L.

102

parameter in the model. This scenario does not seem to apply to any of the models

studied under self-organized criticality.

Appendix

In this Appendix we briefly review the basics of Boolean algebra as it applies to the
numerical simulation of cellular automata.

Boolean algebra consists of binary operations carried out on binary values which
can be either 0 or 1. For example, the result of the operation A AND B gives 1 if
and only if both A and B are 1. Similarly, the operation A OR B gives 1 if either
A or B is 1. Another important operation is the XOR operation. A XOR B gives 1
if either A or B is 1, but not both. Finally, the unary operation NOT reverses the
value of its operand, e.g., NOT 0 gives 1.

Any dynamical rule that can be stated as an operation involving Boolean algebra
can be efficiently implemented on digital computers since they typically carry out
several Boolean operations in parallel. This has been utilized in simulations of the
Ising model [50], where large lattices sizes can be simulated with the help of this
parallelization.

Here we describe how the Game of Life can be formulated in terms of Boolean
algebra. We first need a way to represent the number of live neighbors n of a given
lattice site c. Since n can range from zero to eight, and a Boolean variable can only

take on two values, we write n in binary notation as
n=bg xX 8+box4+b, x2+b x1, (2.39)

where bo,...,63 are all Boolean variables. If we treat c as a Boolean variable taking

on the value 0 if the cell is dead and the value 1 if the cell is alive, the rules of the

103

Game of Life (cf. section 2.4.2) can be implemented as:

c= (c OR by) AND b, AND (NOT b) AND (NOT b3), (2.40)

where c will have the correct value after the logical operations. We also need to
implement the process of addition as a set of Boolean operations. This can be
accomplished by observing that if we are given two Boolean variables t; and to,
we can add them in binary notation by the operations b; = t; AND to, and by =
t; XOR t». Using this procedure, we can add the values of the eight neighbors to
arrive at Eq. (2.39).

In principle, we could use a similar procedure to simulate the sandpile models
of section 2.3.3. However, since in this case the actual number of sites toppling at
any given time is usually much smaller than the number of lattice sites, it is more

efficient to simply update a list of the sites that need to be toppled.

104

References

[1] See, e.g., H.E. Stanley, Introduction to Phase Transitions and Critical Phenom-
ena (Oxford University Press 1971).

[2] See, e.g., K.G. Wilson, Rev. Mod. Phys. 47, 773 (1975).

[3] See, e.g., S.K. Ma, Modern Theory of Critical Phenomena (Benjamin Publishing
Company 1976).

[4] For reviews of 1/f noise, see P. Dutta and P.M Horn, Rev. Mod. Phys. 53,
497(1981); M.B. Weissman, Rev. Mod. Phys. 60, 537(1988).

[5] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 59, 381 (1987) and Phys.
Rev. A38, 364 (1988).

[6] G.I. Barenblatt, Similarity, Self-Similarity, and Intermediate Asymptotics, Con-
sultants Bureau, New York, (1979).

[7] J.H. Conway, unpublished (1970).

[8] P. Bak, K. Chen, and M. Creutz, Nature 342, 780 (1989).

[9] P.W. Bridgeman, Dimensional Analysis, Yale University Press (1931).
[10] E. Buckingham, Phys. Rev. 14, 345 (1914).

[11] G.I. Taylor, Proc. Roy. Soc. 201, 175 (1950).

105

[12] J.M. Carlson, J.T. Chayes, E.R. Grannan, and G.H. Swindle, Phys. Rev. A42,
2467 (1990).

[13] L.P. Kadanoff, S.R. Nagel, L. Wu, and S. Zhou, Phys. Rev. A39, 6524 (1988).

[14] T.C. Halsey, M.H. Jensen, L.P. Kadanoff, I. Procaccia, and B.I. Shraiman,
Phys. Rev. A33, 1141 (1986).

[15] S.L. Kamenomostskaya, Dokl. Akad. Nauk SSSR 116, 18 (1957).

[16] L.P. Kadanoff, A.B. Chhabra, A.J. Kolan, M.J. Feigenbaum, and I. Procaccia,
Phys. Rev. A45, 6095 (1992).

[17] N. Goldenfeld, O. Martin, Y. Oono, and F. Liu, Phys. Rev. Lett. 64, 1361
(1990).

[18] Handbook of Mathematical Functions, M. Abramowitz and I. Stegun eds., Dover
Publications 1970.

[19] B. Grossmann, H. Guo, and M.Grant, Phys. Rev. A41, 4195 (1990); P. Grass-
berger and S. S. Manna, J. Phys. 51, 1077 (1990).

[20] S. S. Manna, J. Stat. Phys. 59, 509 (1990).

[21] J.M. Carlson, J.T. Chayes, E.R. Grannan, and G.H. Swindle, Phys. Rev. Lett.
65, 2547 (1990).

[22] D. Dhar, Phys. Rev. Lett. 64, 1613 (1990).

[23] S. S. Manna, Physica A 179, 249 (1991).

[24] P. Evesque, Phys. Rev. A43, 2720 (1991).

[25] H.M. Jaeger, C.H. Liu, and S.R. Nagel, Phys. Rev. Lett. 62, 40 (1989).

[26] G.A. Held, D.H. Solina, D.T. Keane, W.J. Haag, P.M. Horn, and G. Grinstein,
Phys. Rev. Lett. 65, 1120 (1990);

106

[27] T. Hwa and M. Kardar, Phys. Rev. Lett. 62, 1813 (1989).

[28] A. Diaz-Guilera, Phys. Rev. A45, 8551 (1992).

[29] G. Grinstein and D.-H. Lee, Phys. Rev. Lett. 66, 177 (1991).

[30] G. Grinstein, D.-H. Lee, and S. Sachdev, Phys. Rev. Lett. 64, 1927 (1990).
[31] D. Forster, D.R. Nelson, and M.J. Stephen, Phys. Rev. A16, 732 (1977).
[32] D. Dhar and R. Ramaswamy, Phys. Rev. Lett. 63, 2659 (1989).

[33] T. Hwa and M. Kardar, Phys. Rev. A45, 7002 (1992).

[34] C. Tang and P. Bak, J. Stat. Phys. 51, 797 (1988).

[35] J. Theiler, unpublished.

[36] T. E. Harris, The Theory of Branching Processes, Springer-Verlag (1963).

[37] H. W. Watson and F. Galton, J. Anthropol. Inst. Great Britain and Ireland 4,
138 (1874).

[38] R. Otter, Annals of Math. Stat. 20, 206 (1949).
[39] As quoted in ref [36], p.21.
[40] D. Dhar and S. N. Mayundar, J. Phys. A 23, 4333 (1990).

[41] C. H. Bennett and M. S. Bourzutschky, unpublished; H. J. S. Feder and J.
Feder, Phys. Rev. Lett. 66, 2669 (1991).

[42] K. Christensen, H. C. Fogedby, and H. J. Jensen, J. Stat. Phys. 63, 718 (1991).
[43] J. von Neumann, Collected Works 5, 288 (1961).

[44] For a general review of cellular automata, see §. Wolfram, Rev. Mod. Phys. 55,
601 (1983).

107

[45] E. Ising, Z. Phys. 31, 253 (1925).

[46] For a review of the Potts model, see F. Y. Wu, Rev. Mod. Phys. 54, 235 (1982).
[47] R. J. Glauber, J. Math. Phys. 4, 294 (1963).

[48] K. Kawasaki, Phys. Rev. 150, 285 (1966).

[49] M. Creutz, Phys. Rev. Lett. 50, 1411 (1983), J. Stat. Phys. 42, 823 (1986).
[50] B. Hede and H. J. Herrmann, J. Phys. A 24, L691 (1991).

[51] A. M. A. Verhagen, J. Stat. Phys. 15, 213 (1976); E. Domany and W. Kinzel,
Phys. Rev. Lett 53, 311 (1984).

[52] U. Frisch, B. Hasslacher, and Y. Pomeau, Phys. Rev. Lett. 56, 1505 (1986).
[53] S. Wolfram, J. Stat. Phys. 45, 471 (1986).

[54] H. Chaté and P. Manneville, Europhys. Lett. 14, 409 (1991).

[55] N. Packard and S. Wolfram, J. Stat. Phys. 38, 901 (1985).

[56] E. Fredkin and T. Toffoli, Int. J. Theor. Phys. 21, 219 (1982).

[57] C. H. Bennett and G. Grinstein, Phys. Rev. Lett. 55, 657 (1985).

[58] S. Wolfram, Physica D 10, 1 (1984).

[59] S. Wolfram, Physica Scripta T 9, 170 (1985).

[60] E. R. Berlekamp, J. H. Conway, and R.K. Guy, Winning Ways for Your Math-
ematical Games, Vol. 2, Academic Press (1982).

[61] M. Gardner, Mathematical Games, Sci. Amer. 224, Feb. 112 (1971), and
Wheels, Life and Other Mathematical Amusements, W. H Freeman (1983).

[62] L. S. Schulman and P. E. Seiden, J. Stat. Phys. 19, 293 (1978).

108

[63] C. H. Bennett, private communication.

(64] C. H. Bennett and M. S. Bourzutschky, Nature 350, 468 (1991).
[65] K. Chen and P. Bak, Phys. Lett. A 140, 299 (1989).

[66] P. Grassberger and H. Kantz, J. Stat. Phys. 63, 685 (1991).

[67] P. Bak, K. Chen, and C. Tang, Phys. Lett. A147, 297 (1990).