Table Of ContentDETERMINISTIC WALKS IN QUENCHED RANDOM
ENVIRONMENTS OF CHAOTIC MAPS
TAPIO SIMULA AND MIKKO STENLUND
9
0 Abstract. Thispaperconcernsthepropagationofparticlesthrough
0
a quenched random medium. In the one- and two-dimensional
2
models considered, the local dynamics is given by expanding cir-
n
cle maps and hyperbolic toral automorphisms, respectively. The
a
J particle motion in both models is chaotic and found to fluctuate
4 aboutalineardrift. Intheproperscalinglimit,thecumulativedis-
1 tribution function of the fluctuations converges to a Gaussian one
with system dependent variance while the density function shows
]
S no convergence to any function. We have verified our analytical
D results using extreme precision numerical computations.
.
h
t
a
m
1. Introduction
[
Variants of a mechanical model now widely known as the Lorentz gas
1
v have occupied the minds of scientists for more than a century. Initially
5 proposed by Lorentz [Lo] in 1905 to describe the motion of an electron
7
in a metallic crystal, the model consists of fixed, dispersing, scatterers
8
1 in Rd and a free point particle that bounces elastically off the scatterers
.
1 upon collisions.
0
If the lattice of scatterers is periodic, the model is also referred to as
9
0 Sinai Billiards after Sinai, who proved [Si] that the system (with d = 2)
: is ergodic if the free path of the particle is bounded. In the latter case
v
i it was also proved that, in a suitable scaling limit, the motion of the
X
particle is Brownian [BuSi]. Sinai’s work can be considered the first
r
a rigorous proof of Boltzmann’s Ergodic Hypothesis in a system that
resembles a real-world physical system.
Despite its seeming innocence, the Lorentz gas exhibits a great deal
of complexity. One example is the lack of smoothness of the dynamics
caused by tangential collisions of the particle with the scatterers. An-
other one, the presence of recollisions, is a source of serious statistical
Date: January 14, 2009.
2000 Mathematics Subject Classification. 60F05; 37D20, 82C41, 82D30.
Key words and phrases. Random environments, hyperbolic dynamical systems,
Gaussian fluctuations, Lorentz gas.
1
2 TAPIO SIMULA AND MIKKO STENLUND
Figure 1. Schematic diagram illustrating the qualita-
tive features of our two-dimensional model. The medium
is composed of two different types of square cells sepa-
rated by walls which allow particles to pass only in the
direction of the arrowheads. The zig-zag line shows a
path of a particle through the medium.
difficulties that have not been overcome in the study of the aperiodic
Lorentz gas. For more background, see [Ta, Sz, ChMa, ChDo] and the
references therein.
Our study concerns an idealization of the aperiodic Lorentz gas with
semipermeable walls illustrated in Figure 1. In each cell, there is a con-
figuration of scatterers drawn independently from the same probability
distribution. In our case, the distribution is Bernoulli, so that there
are two possible configurations to choose from inside each cell. As an
important aspect, the environment thus obtained is quenched; once the
scatterer configurations have been randomly chosen, they are frozen for
good, and the only randomness that remains is in the initial data of
the particle. Between the cells are semipermeable walls that allow the
particle to pass through from left to right and from bottom to top,
as shown by the arrowheads, but not in the opposite directions. This
model may be thought of as describing the propagation of particles in
an exotic, anisotropic, medium.
Notice that, as a significant simplification, there is no recurrence;
once a particle leaves a cell, it never returns to the same cell again.
Yet a particle can occupy a single cell for an arbitrarily long time
before moving on to a neighboring one—albeit a long occupation time
has a small probability. Moreover, where, when, and in which direction
the particle exits a cell depends heavily on the scatterer configuration
DETERMINISTIC WALKS IN QUENCHED RANDOM ENVIRONMENTS... 3
inside the cell, in addition to the position and direction of the particle
at entry. Inside each cell, the dynamics is chaotic and hyperbolic.
In our one- and two-dimensional idealizations, the billiard dynamics
is replaced by discrete dynamical systems acting in each cell. In other
words, acting on the particle’s current position by a map associated
to the current cell gives its position one time unit later. In dimension
one the maps associated to the cells are expanding circle maps while in
dimension two they are hyperbolic toral automorphisms. Such choice
of maps retains the chaotic and hyperbolic nature of the problem. A
closely related model has been studied in [AySt, AyLiSt].
Our objective is to understand certain statistical properties of the
motion. More precisely, we are interested in how the particle dis-
tribution evolves with time when the initial distribution is uniform
and supported on one initial cell. We make several analytical proposi-
tions, which we verify numerically. We show that, on the average, the
particles follow a linear drift and that, after taking a suitable scaling
limit, the cumulative distribution of the fluctuations about the mean is
Gaussian. Moreover,thedriftandvariancearethesamefor(almost)all
environments drawn from the same distribution. Nevertheless, the fine
structure of the particle distribution is peculiar due to the quenched
environment. In particular, the density function does not converge to
that of a normal distribution. In fact, it does not converge at all in the
aforementioned scaling limit.
Acknowledgements. We are indebted to Arvind Ayyer and Joel
Lebowitz for stimulating discussions. Tapio Simula is supported by
the Japan Society for the Promotion of Science Postdoctoral Fellow-
ship for Foreign Researchers. Mikko Stenlund is partially supported by
a fellowship from the Academy of Finland.
2. One-dimensional model: expanding maps on the circle
2.1. Introduction. Imagine tiling the nonnegative half line [0, ) so
that each interval—or tile—I = [k,k + 1) with k N carries∞a la-
k
∈
bel ω(k) that equals either 0 or 1. Such a tiling can be realized by
flipping a coin for each k and encoding the outcomes in a sequence
ω = (ω(0),ω(1),...) 0,1 N called the environment. The coin could
∈ { }
be balanced but the tosses are independent, with Prob(ω(k) = i) = p
i
for all k. Below, P will stand for the corresponding probability mea-
p0
sure on the space of Bernoulli sequences ω. In each experiment we
freeze the environment—meaning that we work with onefixed sequence
ω at a time.
4 TAPIO SIMULA AND MIKKO STENLUND
The dynamics in our model is generated by the following definitions.
Let A ,A 2,3,... and define the circle maps T (x) = A x mod 1.
0 1 i i
An experim∈en{t compr}ises iterating many times the map R S1
R S1 : (v,x) (v +A x x,T (x)), where [v] is the×integ→er
ω([v]) ω([v])
× (cid:55)→ −
part of v and ω([v]) is the corresponding component of ω. The actual
initial condition is (x,x), with x [0,1), although in the treatment be-
∈
low we will also consider initial conditions of the form (V +x,x) with
V N arbitrary. Let P denote the Lebesgue measure (i.e., the uni-
form∈ probability distribution) on the circle S1 and E the corresponding
expectation.
The difference A x x above can be thought of as the jump
ω([v])
−
experienced by a particle. In other words, the evolution of the v-
coordinate on the line keeps track of the journey of a particle, from the
point of view of the particle. The x-coordinate is just the projection
onto the circle (and is needed for technical purposes): if (v ,x ) is the
n n
image of (v ,x ) = (V +x,x) after n steps, then x = v mod 1.
0 0 n n
The model thus describes the deterministic motion of a particle in
a randomly chosen, but fixed, environment. In probability jargon, the
particle performs a deterministic walk in a quenched random environ-
ment. The challenge is that the map determining the position of the
particle in the nearest future depends on the tile I the particle is in
[v]
through the label ω([v]) of the tile. That is, the motion of the particle
is guided by the a priori chosen environment.
Example 1. A concrete example is obtained by choosing A = 2, A =
0 1
3, and p = p = 1.
0 1 2
In general,
n−1
(cid:88)
v = [v ]+x+(A 1)x+ (A 1)T T (x).
n 0 ω([v0])− ω([vj])− ω([vj−1])◦···◦ ω([v0])
j=1
Notice that already in the expression of v , the initial V and x and the
2
environment ω are present in a complicated fashion:
A = A .
ω([v1]) ω([x+(Aω(V)−1)x])
The complication referred to earlier is obvious; for each x (and ω) we
have to compute which map the symbol A stands for. Each v
ω([v1]) n
(n 2) is generically a piecewise affine function of x and the number
≥
of discontinuities grows exponentially with n.
We anticipate that, for large values of n, v behaves statistically (in
n
the weak sense) like
v (nD,nσ2), (1)
n
≈ N
DETERMINISTIC WALKS IN QUENCHED RANDOM ENVIRONMENTS... 5
where D is a deterministic number called the drift and (nD,nσ2)
N
stands for a real-valued, normally distributed, random variable with
meannD andvariancenσ2. Inprinciple, D andσ2 coulddependonthe
environment ω, but remarkably it turns out that they do not, as long
as the environment is typical. By typicality we mean that ω belongs
to a set whose P -probability is one and whose elements enjoy good
p0
statistical properties such as the convergence of 1# k < n ω(k) = 0
n { | }
to the limit p .
0
It is reasonable to expect that the limit
v (x)
n
lim
n→∞ n
exists and has the same value for almost all x 1. Thus, we are led to
conclude that
(cid:18) (cid:19)
v (x) 1
n
D = E lim = lim E(v (x)).
n
n→∞ n n→∞ n
The final equality follows from the bounded convergence theorem.
Let us consider the (asymptotically) centered random variable
X = v nD,
n n
−
which measures the fluctuations of v relative to the linear drift. Pro-
n
vided (1) is true, X is approximately Gaussian with variance nσ2.
n
More precisely, we would like to know if √1 X converges in distribu-
n n
tion to (0,σ2). By definition, this means that, for any fixed y R,
N ∈
(cid:18) 1 (cid:19) 1 (cid:90) y
lim P X y = e−s2/2σ2ds.
n
n→∞ √n ≤ √2πσ
−∞
2.2. Markov partition. We next reduce the deterministic walk in a
random environment to a random walk in a random (still quenched)
environment which is easier to treat. This can be done using a Markov
property of the tiling that allows us, in the statistical sense, to ignore
the exact position of the particle and only keep track of the tile it is
occupying.
Let [ ] denote the integer part of a number. If we define
·
V = [v ] and x = v [v ], (2)
n n n n n
−
1Thiscannotholdforallx. Forinstance,ifx= 1 ,then1=v =v =...
kAω(0) k k+1
and the process stops.
6 TAPIO SIMULA AND MIKKO STENLUND
then the earlier dynamics with the initial condition (v ,x ) = (V +
0 0 0
x ,x ) is equivalent to
0 0
(cid:2) (cid:3)
V = V + A x
n+1 n ω(Vn) n (3)
(cid:2) (cid:3)
x = A x A x .
n+1 ω(Vn) n − ω(Vn) n
Recall our convention [0, ) = (cid:83)∞ I , where I = [k,k +1) is
∞ k=0 k k
called a tile. Suppose now that v I . This is equivalent to V = k.
n k n
∈
As before, we are interested in the probability distribution of v , but
n
this time only at the level of tiles. More precisely, we wish to know
the probability distribution of V . This is the probability vector ρ(n) =
n
(ρ(n),ρ(n),...) where the numbers
0 1
ρ(n) = P(V = k)
k n
are such that (cid:80)∞ ρ(n) = 1.
k=0 k
We now consider the dynamical system being initialized with the
condition (v ,x ) = (V + x ,x ), where V N and x [0,1) are
0 0 0 0 0 0 0
∈ ∈
random variables; x is uniformly distributed and independent of V ,
0 0
but for the moment we do not specify the distribution ρ(0) = P(V = k)
k 0
which is the initial distribution of the chain V .
n
Since each A I is exactly the union of a few of the intervals I 2, the
i k k(cid:48)
collection I is a simultaneous Markov partition for the two maps.
k
{ }
We then obtain the Markov property
P(V = k V = k ,...,V = k ) P(V = k V = k )
n n n−1 n−1 0 0 n n n−1 n−1
| ≡ |
for arbitrary histories. Thus, V is a time-homogeneous Markov chain
n
on the countably infinite state space N with the transition probabilities
(cid:40)
P(v I v I ) = 1 if l 0,1,...,A 1 ,
γ = n+1 ∈ k+l| n ∈ k Aω(k) ∈ { ω(k) − }
k→k+l
0 otherwise.
Notice that the above holds for any environment, ω, but the resulting
Markov chain does depend on the choice of ω.
Defining the transition matrix Γ = (γ ) , ρ(n) = ρ(n−1)Γ. Thus,
k→k(cid:48) k,k(cid:48)
ρ(n) = ρ(0)Γn (4)
for an arbitrary initial distribution. In particular,
ρ(0) = (1,0,0,...)
corresponds to initializing the dynamical system at (v ,x ) = (x,x)
0 0
with 0 x < 1 arbitrary.
≤
(cid:104) (cid:17)
2A maps k+ l ,k+ l+1 affinely onto [k+l,k+l+1).
i Ai Ai
DETERMINISTIC WALKS IN QUENCHED RANDOM ENVIRONMENTS... 7
In principle, (4) provides us with complete statistical understanding
of the dynamics. For instance, the drift can be expressed as
∞
1 1 1 (cid:88)
D = lim E(v ) = lim E(V ) = lim kρ(n).
n→∞ n n n→∞ n n n→∞ n k
k=0
In practice, calculating Γn for large values of n is difficult.
2.3. Drift and variance. Foreach(i,j) 0,1 2 thetransitionprob-
∈ { }
ability at time n from a tile labeled i to a tile labeled j is
α (n) = P(ω(V ) = j ω(V ) = i).
ij n+1 n
|
The analysis of this quantity is subtle, because it depends on the tiling.
For instance, if ω = (0,0,...), then P(ω(V ) = 0 ω(V ) = 0) = 1.
n+1 n
The conditional probability P P (ω(V ) = j ω|(V ) = i) equals
× p0 1 | 0
(cid:18) (cid:18) (cid:19) (cid:19) (cid:18) (cid:19)
1 1 1
α∗ = δ + 1 p +(1 δ ) 1 p ,
ij ij A − A i − ij − A j
i i i
because the elements ω(k) of the tiling are independent. Here p is the
i
Bernoulli probability of getting an i in the tiling. We think of α∗ as
ij
an effective transition probability which only depends on the statistical
properties of the tiling.
We expect the actual transition probability α (n) to converge to the
ij
effective value α∗ with increasing time,
ij
lim α (n) = α∗.
ij ij
n→∞
This is so because, as n increases, the position of the particle at time n
depends on the tiling on an increasing subinterval of [0, ) and should
∞
therefore reflect increasingly the statistics of the tiling instead of its
local details.
Moreover,
(cid:18) (cid:19)
p 1 p
lim(α∗)n = − (5)
p 1 p
n→∞
−
for a p (0,1) that can be found by diagonalizing α∗ or by solving the
∈
equilibrium equation (p,1 p)α∗ = (p,1 p):
− −
(cid:16) (cid:17)
p 1 1
p = 0 − A1 = p0A0(A1 −1) .
1 p 1 p 1 A A p A p A
− 1A0 − 0A1 0 1 − 1 1 − 0 0
For instance, in the case of Example 1 we obtain p = 4.
7
Notice that, for any probability vector (q,1 q),
−
(q,1 q) lim(α∗)n = (p,1 p).
− n→∞ −
8 TAPIO SIMULA AND MIKKO STENLUND
The probability vector
k
(cid:89)
(q,1 q) α(n) = (q(k),1 q(k))
− −
n=0
will converge to some (q∗,1 q∗), because α(n) α∗. In fact,
− →
2N 2N
(cid:89) (cid:89)
lim (q,1 q) α(n) = lim (q(N),1 q(N)) α(n)
N→∞ − N→∞ −
n=0 n=N+1
= (q∗,1 q∗) lim (α∗)N = (p,1 p),
− N→∞ −
as N .
→ ∞
We interpret the result above so that P(ω(V ) = 0) p and
n
→
P(ω(V ) = 1) 1 p as n . That is, along a given (typi-
n
→ − → ∞
cal) trajectory, the fraction of time the particle spends in a tile labeled
0 is p:
# k < n ω(V ) = 0
k
lim { | } = p. (6)
n→∞ n
Notice that p does not depend on the (typical) tiling.
2.3.1. Drift. Define the jumps ξ = V V (i 1). Then V =
i i i−1 n
(cid:80)n ξ . We also denote ξ(j) a random−variable th≥at takes values in
i=1 i
0,...,A 1 with uniform distribution. For a (typical) tiling,
j
{ − }
E(V ) (cid:80)n E(ξ )
D = lim n = lim i=1 i = pE(ξ(0))+(1 p)E(ξ(1))
n→∞ n n→∞ n −
A 1 A 1 pA +(1 p)A 1
0 1 0 1
= p − +(1 p) − = − − .
2 − 2 2
We also claim that this equals lim Vn for almost all x. In the case
n→∞ n
of Example 1, D = 5.
7
2.3.2. Variance. The variance is
(cid:18) (cid:19)
X Var(V )
σ2 = lim Var n = lim n .
n→∞ √n n→∞ n
Let us assume A A and study the process W = (cid:80)n ζ having
0 ≤ 1 n i=1 i
the i.i.d. increments ζ whose distribution is Prob(ζ = k) = p + 1−p
i 1 A0 A1
if 0 k < A and Prob(ζ = k) = 1−p if A k < A . The increments
≤ 0 1 A1 0 ≤ 1
have been chosen so that W mimics V as closely as possible and has
n n
the same variance. For instance, staying in the same tile (ζ = 0) has
1
probability p +1−p, where p is the probability of being in a tile labeled
A0 A1
0 and 1 is the probability of staying in that tile, while the second
A0
term accounts similarly for the case of label 1. Then Mean(W ) =
n
DETERMINISTIC WALKS IN QUENCHED RANDOM ENVIRONMENTS... 9
nMean(ξ ) = nD. Setting K(m) = (cid:80)m−1k2 = 1m3 + 1m2 + 1m, the
1 k=0 3 2 6
variance of W is
n
(cid:18) (cid:19)
p 1 p
Var(W ) = nVar(ζ ) = n K(A )+ − K(A ) D .
n 1 0 1
A A −
0 1
Using the values p = 4 and D = 5 obtained from Example 1, the
7 7
formula above gives 1 Var(W ) = 24.
n n 49
2.4. Sensitivity on the initial condition. Let us next consider the
Lyapunov exponent
n
1 dv 1 (cid:88) dv
n k
λ = lim ln = lim ln
n→∞ n dx n→∞ n dvk−1
k=1
which measures the exponential rate at which two nearby initial points
drift apart under the dynamics. Above, the chain rule has been used.
Recall the notation introduced in (2) and that v = x = x. As v =
0 0 k
(cid:0) (cid:1)
v + A 1 x ,
k−1 ω(V ) k−1
k−1 −
dv dA dx
k = 1+ ω(Vk−1)x +(cid:0)A 1(cid:1) k−1.
k−1 ω(V )
dv dv k−1 − dv
k−1 k−1 k−1
With probability zero v is an integer, in which case v = V ,
k−1 k−1 k−1
x = 0, and the process stops. We assume that v is not an integer.
k−1 k−1
Then dAω(Vk−1) = 0, dxk−1 = 1, and dvk = A , such that, by (6),
dvk−1 dvk−1 dvk−1 ω(Vk−1)
n
1 (cid:88)
λ = lim lnA = plnA +(1 p)lnA
ω(V ) 0 1
n→∞ n k−1 −
k=1
and is positive. Roughly speaking, the distance between two very
nearby trajectories thus grows like eλn = (cid:0)ApA1−p(cid:1)n, which is tan-
0 1
tamount to chaos.
2.5. Numerical study. In orderto study the model introducedabove
numerically we first create the random tiling (or environment) ω of
length2n+1wherenisthenumberofjumpstobeperformedinasingle
trajectory. Thisguaranteesthateverypossiblepathfitsinsidethetiling
although some computational effort could be saved by choosing the
number of tiles closer to [nD]. Notice also that while in principle new
tiles could be added dynamically to the end of the tiling as required, it
is computationally far more efficient to construct the tiling as a static
entity in the beginning of the computation.
In practice, the tiling is generated by producing a vector of pseudo-
random numbers distributed uniformly on the interval (0,1) using the
MersenneTwisteralgorithm. Thelabelofthetileω(k)isthenobtained
10 TAPIO SIMULA AND MIKKO STENLUND
3
0.717
2.5
0.716
2
n
/ 0.715
n
V
1.5
0.714
5 6
10 10
1
0 1 2 3 4 5 6
10 10 10 10 10 10 10
n
Figure 2. Position of a particle V divided by the num-
n
ber of jumps n taken for a single typical trajectory.
The horizontal line denotes the exact value for the drift
D = 5/7. The inset shows a blow-up of a part of the
main figure.
by rounding the number on each tile k to the nearest integer. The
computational tiling ωˆ is finalized by the operation ωˆ(k) = ω(k)(A
1
−
A )+A yielding a vector whose each element is either 2 or 3.
0 0
Each ensemble member (particle trajectory) is initialized by generat-
ingapseudo-randomnumbertodeterminethestartingpointx (0,1)
0
∈
of the trajectory. The subsequent particle positions are determined by
the underlying tiling. The jumping process could be performed deter-
ministically by keeping track of the exact position v of the particle.
n
However,theMarkovpropertyoftheprocessprovidesusasuperiorway
of obtaining the desired statistics stochastically. In this algorithm, be-
fore every jump, we sample a new pseudo-random number d from the
interval (0,ωˆ(k)) depending on the current tile k. Then a jump to the
tile k + [d] is made and the whole procedure is repeated n times to
produce a single trajectory.
Figure 2 shows a typical trajectory of n = 106 jumps obtained using
the above prescription. The position of the particle V divided by the
n
number of jumps n taken is clearly seen to saturate to the analytical
value for the drift D = 5/7, plotted as a straight line in the figure. The
inset shows the late-time evolution of the drift of the particle.
Figure 3 displays the computed probability density function for the
randomvariable(V nD)/√nobtainedusing107 trajectoriesoflength
n
−
