Anomalous diffusion and stretched exponentials in heterogeneous glass-forming
liquids: Low-temperature behavior
J. S. Langer and Swagatam Mukhopadhyay
Dept. of Physics, University of California, Santa Barbara, CA 93106-9530
(Dated: November 5, 2018)
We propose a model of a heterogeneous glass forming liquid and compute the low-temperature
behavior of a tagged molecule moving within it. This model exhibits stretched-exponential decay
of the wavenumber-dependent, self intermediate scattering function in the limit of long times. At
temperatures close to the glass transition, where the heterogeneities are much larger in extent than
the molecular spacing, the time dependence of the scattering function crosses over from stretched-
exponential decay with an index b = 1/2 at large wave numbers to normal, diffusive behavior with
b = 1 at small wavenumbers. There is a clear separation between early-stage, cage-breaking β
relaxation and late-stage α relaxation. The spatial representation of the scattering function exhibits
an anomalously broad exponential (non-Gaussian) tail for sufficiently large values of the molecular
displacement at all finite times.
I. INTRODUCTION
A growing body of evidence, both experimental [1, 2, 3]
and numerical [4, 5], points to the existence of intrinsic
spatial heterogeneities in glass forming liquids at tem-
peratures slightly above the glass transition, and the rel-
evance of these heterogeneities to stretched exponential
decay of correlations and anomalous, non-Gaussian diffu-
sion. Here we propose a simple model of a heterogeneous
glass former and use it to compute observable proper-
ties of a tagged molecule moving within it. In particular,
we show that this model naturally predicts stretched-
exponential decay when the scale of the heterogeneity is
much larger than the molecular spacing; and we compute
deviations from Gaussian displacement distributions dur-
ing intermediate stages of the relaxation process.
Our model is loosely motivated by the excitation-chain
(XC) theory of the glass transition proposed recently by
one of us (JSL) [7, 8] This theory suggests that a glass-
forming liquid at a temperature T not too far above the
glass transition temperature T0 consists of fluctuating do-
mains of linear size R∗(T ), within which the molecules
are frozen in a glassy state where they have little or no
mobility. R∗(T ) is a length scale that characterizes dy-
namic fluctuations in the equilibrium states of these sys-
tems. The theory predicts that R∗(T ) diverges near T0
like (T − T0)−1 and decreases to zero at the upper limit
of the super-Arrhenius region, T = TA.
Characterizing the more mobile material that lies out-
side the hypothetical frozen domains remains one of the
deeper challenges in glass physics. In the language intro-
duced by Widmer-Cooper and Harrowell [9], the mobile
molecules reside in regions of high “propensity.” Those
authors perform molecular dynamics simulations of a
glasslike, two-dimensional, binary mixture. They average
the mean-square displacement of each molecule over an
isoconfigurational ensemble in which the initial positions
of all the molecules are fixed, but the initial velocities are
selected at random from a Boltzmann distribution. The
propensity of a molecule is its mean-square displacement
during a structural relaxation time. Contour maps of
propensity as a function of initial molecular positions do
indeed exhibit well defined domains of low and high mo-
bility.
Widmer-Cooper and Harrowell compare their propen-
sity maps with corresponding plots of the potential en-
ergy, which is taken to be indicative of the local struc-
ture. Interestingly, the propensity appears to be almost
entirely uncorrelated with the potential energy. In retro-
spect, this lack of correlation between short-range struc-
ture and mobility may not be surprising. Such a corre-
lation would be expected for a slowly coarsening system
consisting of domains of two structurally distinct phases,
for example, a slightly supersaturated (mobile) fluid in
which (frozen) crystallites are growing. In that case, the
domains are very nearly but not quite in thermodynamic
coexistence with each other. A true glass-forming liquid,
however, is in a state of thermodynamic equilibrium. The
temperatures and chemical potentials are uniform every-
where, independent of propensity; and the local pair cor-
relations, which ought to be determined primarily by the
temperature in an equilibrium state, appear to be statis-
tically uniform as well.
A local quantity that does correlate with the propen-
sity is the mean-square vibration amplitude of a molecule
averaged over times much longer than its period of oscil-
lation but much shorter than the time required for irre-
versible rearrangements to occur. A large mean-square
vibration amplitude – or, equivalently, a small Debye-
Waller factor – implies that a molecule is participat-
ing in a soft, low frequency, elastic mode. It is well
known that soft modes are abundant in amorphous mate-
rials near jamming transitions.[10] Moreover, it appears
(for example, in Fig.4 of [10]) that the largest displace-
ments in such modes lie predominantly along one dimen-
sional, chainlike paths. Thus, regions of high propensity
may be elastically soft, and may somehow be correlated
with excitation-chain activity. Of course, local variations
of elastic stiffness are ultimately structural in nature;
but soft modes are collective phenomena involving many
molecules, and are not easily detected by measuring near-
http://arxiv.org/abs/0704.1508v2
neighbor pair correlations.
In the model proposed here, glassy domains of low
propensity are separated from each other by regions of
higher propensity in which many of the molecules are
mobile. This domain structure undergoes persistent fluc-
tuations on time scales of the order of τ∗α = R
∗2/Dα,
where Dα = ℓ
2/τα is the diffusion constant associ-
ated with α relaxation, ℓ is approximately the average
molecular spacing, and τα is the strongly temperature-
dependent, super-Arrhenius, α relaxation time. In con-
trast to the slow fluctuations of domain boundaries, indi-
vidual molecules within high-propensity regions are ran-
dom walkers with a diffusion constant DM = ℓ
2/τM ,
where τM is a temperature-dependent time scale – pre-
sumably not super-Arrhenius – that characterizes mobile
molecular displacements.
To a first approximation, a tagged molecule at any
given time is either frozen in a glassy domain or is mobile.
If the former, then it waits a time of the order of τ∗α before
it is encountered by mobile neighbors and itself starts
to undergo displacements. If the latter, then it diffuses
for some distance until it finds itself frozen again in a
glassy domain. That distance is at least a few molecular
spacings ℓ; but it seems more likely to be proportional to
R∗ if the domains are large and if diffusive hopping takes
place more readily within the mobile regions than across
their boundaries and into the glassy domains. Therefore,
we assume that the time elapsed during mobile motion
is of the order of R∗2/DM = (R
∗/ℓ)2 τM . The ratio of
these two time scales is τα/τM ≡ ∆. At temperatures
well below TA, we expect that ∆ ≫ 1.
A similar picture of diffusion in fluctuating regions of
varying mobility emerges in kinetically constrained mod-
els [11, 12] and in a recent model of a gel [13]. Our
analysis is most closely related to the diffusion model of
Chaudhuri et al [14], which has been a valuable start-
ing point for the present investigation. We go beyond
[14] by including a mechanism for producing late-stage
stretched-exponential relaxation, and by distinguishing
that mechanism from the diffusion of mobile molecules.
In contrast to [14], our molecules switch back and forth
between frozen and mobile behavior and thus explore the
geometry of the domain structure. As a result, our model
exhibits a clear separation between the slow α relaxation
and the faster mobile motions that we interpret as β re-
laxation. It also respects time-translational symmetry,
which seems to be violated in [14].
The preceding discussion pertains just to the case of
large-scale heterogeneities at temperatures only a little
above T0 and well below TA. Some of the most interesting
and experimentally accessible physics, however, occurs
near the crossover to small-scale heterogeneity near TA.
Constructing a theory of diffusion in the crossover region
requires a temperature-dependent analysis, which will be
described in a following paper.
II. BASIC INGREDIENTS OF THE MODEL IN
THE LIMIT OF LARGE-SCALE
HETEROGENEITY
As outlined above, there are two different mechanisms
to be described probabilistically in this sytem, each op-
erating on its own characteristic time scale. The separa-
tion between time scales, characterized by a large value
of ∆, suggests that the natural mathematical language in
which to discuss this system is that of a continuous time
random walk (CTRW) in which a tagged molecule al-
ternates between long waits in glassy domains and faster
displacements in mobile regions. In this Section, we com-
pute the waiting time distributions for molecules in glassy
and mobile domains, and the corresponding probability
distributions for the mobile displacements.
We immediately encounter a problem, however, be-
cause there are more than just two characteristic time
scales in this system. Analog experiments [3] and numer-
ical simulations [4] typically follow the motions of parti-
cles starting from their initial positions, and resolve these
trajectories on length scales smaller than the interparti-
cle spacing, and on time scales of the order of the shortest
vibrational periods. This is a continuous range of scales.
At its long-time limit, it includes the slow modes that
are involved in molecular rearrangements, whose ampli-
tudes may exceed the Lindemann melting criterion in the
mobile regions. Thus, the cage-breaking, β-relaxation
mechanism may be part of the same continuous range of
time scales that includes the intra-cage vibrational mo-
tions. Indeed, no sharp distinction between vibrational
and cage breaking time scales seems to appear in exper-
imental data. If these are not clearly distinct kinds of
events, then the CTRW approximation is not valid in
the short-time limit.
In view of this difficulty, we choose here to consider
only the β and α time scales, τM and τα respectively. By
neglecting the shorter time scales, i.e. those of the order
of τ0 ∼ femtoseconds - picoseconds, we imply that we are
not resolving length scales much smaller than the molec-
ular spacing ℓ. We then assume that the probability that
a molecule has left its cage at times of order τM ≫ τ0 is
determined by the probability that it has mobile neigh-
bors and participates in local molecular rearrangements.
At the longest time scales τα ∼ seconds (by definition,
the viscous relaxation time at the glass temperature),
the probability that a molecule has escaped from a large
glassy domain is determined by the probability that it
has been encountered by the boundary of the domain.
Consider first the slowest motions, i.e. those associated
with fluctuations of the domain boundaries on length
scales R∗ and time scales τ∗α. Define t
∗ ≡ t/τ∗α, where
t is the physical time in seconds; and let ψG(t
∗) be the
normalized probability distribution for the time that a
molecule spends in a glassy domain before entering a mo-
bile region. Write this distribution in the form
dρW (ρ)
, (2.1)
where ρ is the linear size of a domain in units R∗. W (ρ)
is a normalized distribution over these sizes, and the re-
maining factor inside the integrand is a normalized dis-
tribution over t∗. The quantity ρ−2 appearing in the
exponential in Eq.(2.1) is the lowest eigenvalue of the
diffusion kernel for a molecule moving in a domain of
size ρ. Note the similarity to the trapping models dis-
cussed, for example, in [15, 16, 17]; but also note that
the “trap” here is a two dimensional subspace bounding
a three dimensional domain. The quantity that is fluc-
tuating in a normal diffusive manner is the mean square
of the distance between the tagged molecule and the do-
main boundary. Because the boundary surrounds the
molecule, the direction of the diffusive drift is irrelevant;
any point of contact on the boundary is equivalent to any
other. Therefore, it makes little difference whether the
molecule is the diffuser and the boundary is the target, or
– as in this case – the boundary diffuses and the molecule
is the target.
Strictly speaking, Eq.(2.1) is a long-time approxima-
tion; at shorter times, the higher eigenmodes of the diffu-
sion kernel make non-negligible contributions. However,
this approximation is qualitatively adequate for our pur-
poses at shorter times as well because ψG(t
∗) is well be-
haved at small t∗ and, as mentioned previously, we are
not trying to include the short-time behavior in this func-
tion. If we assume that the distribution over values of ρ
is Gaussian, W (ρ) ∝ ρ2 exp(−ρ2) in three dimensions,
∗) = 2 e−2
t∗ . (2.2)
That is, we find a rudimentary but nontrivial stretched
exponential of the form exp (−const.× t∗b) with b = 1/2.
There is little reason to expect that W (ρ) remains
Gaussian out to the large values of ρ that determine
the long-time behavior of ψG(t
∗). Deviations from the
Gaussian would produce different indices. For example,
as the temperature increases and the glassy domains be-
come smaller, their size distributions might cut off more
sharply than a Gaussian. Thus, if W (ρ) ∝ exp(−ρm),
with m ≥ 2, then b = m/(m + 2) → 1 as m → ∞. The
problem of computing this distribution or, equivalently,
values of m from first principles may eventually become
solvable as we learn more about the statistical mechan-
ics of glass forming liquids; but that problem is beyond
the scope of this investigation. For present purposes, the
important points are that a plausible distribution W (ρ)
produces stretched-exponential behavior, and that the
mechanism by which this happens could produce a tem-
perature dependent index b.
Next, consider the waiting-time distribution ψM (t
for a molecule in a mobile region. This situation is qual-
itatively different from that of a molecule in a glassy do-
main because the distance traveled by the molecule is
related to the time during which it remains mobile. As
stated earlier, we assume that a mobile molecule diffuses
a distance of order R∗ before reentering a glassy region,
so that its residence time in the mobile region is of order
τM (R
∗/ℓ)2 which, in t∗ units, is simply ∆−1. Then, for
simplicity, assume an exponential waiting-time distribu-
tion:
ψM (t
∗) = ∆ e−∆ t
. (2.3)
Compared to ψG(t
∗) in Eq.(2.2), ψM (t
∗) is sharply
peaked near t∗ = 0 if, as expected, ∆ is large. To com-
plete the model of mobile motion, we need the conditional
probability pM (r
∗, t∗) for diffusion over a scaled distance
r∗ = r/R∗ in time t∗. The natural choice is
pM (r
∗, t∗) =
(2 π∆ t∗)3/2
2∆ t∗
, (2.4)
which is a normalized, three dimensional distribution
over r∗.
III. CONTINUOUS-TIME RANDOM WALKS
The next step is to translate these physical ingredients
of the model into the language of continuous-time ran-
dom walks.[19, 20] Define two different probability distri-
bution functions, nG(r
∗, t∗) and nM (r
∗, t∗), for molecules
starting, respectively, in glassy domains or mobile re-
gions, and moving distances r∗ in times t∗. Each molec-
ular trajectory consists of a sequence of transitions be-
tween glassy domains and mobile regions. A single tran-
sition in which a molecule starts in a glassy domain at
time t∗1 and ends in a mobile region at time t
2 (without
having changed its actual position) occurs with probabil-
ity ψG(t
2−t∗1). Similarly, a transition in which a molecule
starts in a mobile region at time t∗1 and position r
1, and
ends in a glassy domain at time t∗2 and position r
2, has
probability ψM (t
2 − t∗1) pM (|r∗2 − r∗1|, t∗2 − t∗1). The prob-
ability that a molecule starts or arrives in a glassy do-
main at time t∗1 and is still there at the final time t
∗ − t∗1), where
∗′) dt∗
= (1 + 2
t∗) e−2
t∗ . (3.1)
Finally, the probability that a molecule starts or arrives
in a mobile region at time t∗1 and position r
1 and is still
in that region at the final time t∗ and at position r∗, is
φM (t
∗ − t∗1) pM (|r∗ − r∗1|, t∗ − t∗1), where, in analogy to
Eq.(3.1),
φM (t
ψM (t
1) dt
1 = e
−∆ t∗ . (3.2)
The multiple convolutions that describe each trajec-
tory are, as usual, converted into products by computing
Fourier-Laplace transforms of each function:
ñj(k, u) =
e−ik·r
∗, t∗) dt∗, (3.3)
where j = G,M . For jumps starting in glassy regions,
we need:
ψ̃G(u) =
∗) dt∗, (3.4)
φ̃G(u) =
1− ψ̃G(u)
. (3.5)
For the mobile regions, the spatial Fourier transform of
the conditional probability distribution pM (r
∗, t∗) is
p̂M (k, t
∗) = exp
; (3.6)
therefore, define
fM (k, u) ≡
p̂M (k, t
∗)ψM (t
∗) dt∗
u+∆(1 + k2/2)
, (3.7)
gM (k, u) ≡
p̂M (k, t
∗)φM (t
∗) dt
u+∆(1 + k2/2)
. (3.8)
Putting these pieces together, and summing over in-
definitely many individual jumps in each trajectory, we
ñG(k, u) =
φ̃G(u) + gM (k, u) ψ̃G(u)
1− ψ̃G(u) fM (k, u)
NG(k, u)
W (k, u)
; (3.9)
ñM (k, u) =
φ̃G(u) fM (k, u) + gM (k, u)
1− ψ̃G(u) fM (k, u)
NM (k, u)
W (k, u)
; (3.10)
where
NG(k, u) =
1− ψ̃G(u)
(1 + k2/2) + u/∆; (3.11)
NM (k, u) = 1− ψ̃G(u) + u/∆; (3.12)
W (k, u) = 1− ψ̃G(u) + k2/2 + u/∆. (3.13)
In both Eqs.(3.9) and (3.10), in the numerators of the
first expressions, the first term corresponds to trajecto-
ries that end in glassy domains, and the second to those
that end in mobile regions.
-6 -4 -2 0 2
FIG. 1: Intermediate scattering function Fs(k, t
∗) with PG =
0.5 and ∆ = 100, for k = 0.8, 1, 1.2, 1.5, 2, 3, 5, 10, reading
from top to bottom.
IV. INTERMEDIATE SCATTERING
FUNCTION
The conventional way to study the diffusion mecha-
nisms discussed here is to measure the self intermediate
scattering function Fs(k, t), which, in the present no-
tation, is the mixed Fourier-time representation of the
weighted average of nG and nM . That is,
Fs(k, t
∗) = PG n̂G(k, t∗) + (1− PG) n̂M (k, t∗), (4.1)
where PG is the probability that a molecule starts its
motion in a frozen, glassy environment; and
n̂j(k, t
∫ +i∞
ñj(k, u), j = G, M. (4.2)
The inversions of Laplace-transforms in Eq.(4.2) are
nontrivial because the stretched-exponential function
ψ̃G(u), defined in Eq.(3.4), has an essential singularity
at u = 0. For t∗ > 0, the integrations over u in Eq.(4.2)
must be performed by closing the contour in the negative
half u-plane, which requires that ñj(k, u) be analytically
continued to points u = |u| exp(iθ). To evaluate ψ̃G(u)
at such points, using Eq.(3.4), let t∗ = y2 and then ro-
tate the y contour away from the positive real axis to
a line y = ξ exp(−iθ/2), 0 < ξ < ∞ so that the inte-
grand always decreases rapidly at infinity. For example,
for u→ −w ± iǫ, θ = ±π,
ψ̃G(−w ± iǫ) = −4
ξ dξ e−wξ
2±2iξ
e−1/w
≡ A(w) ∓ i B(w), (4.3)
where Erfi is the imaginary error function. The imagi-
nary part of ψ̃G(−w ± iǫ) is the discontinuity across a
cut along the negative u axis. With this formula, it is
straightforward to compute the discontinuity across the
cut for the ñj(k, u) in Eqs.(3.9) and (3.10), and in this
way to compute the scattering functions by closing the
contour around this cut. That is:
n̂j(k, t
e−w t
Nj(k,−w + i0)
W (k,−w + i0)
, j = G, M. (4.4)
The integrands are well behaved as w → 0+, therefore it
is possible to let the lower limit of integration be w = 0.
The graphs in Fig.1 show Fs(k, t
∗), computed nu-
merically from Eq.(4.4), as a function of log10(t
∗), for
∆ = 100, PG = 0.5, and for a sequence of wavenumbers
k. These are low-temperature scattering functions, where
∆ ≫ 1, which means that molecules spend much longer
times frozen in glassy domains than they do when moving
in mobile regions. The choice PG = 0.5 is made primar-
ily for clarity; but it seems likely that this value of PG
describes a system that is well into its non-Arrhenius,
anomalous-diffusion regime. The set of curves in Fig.1
closely corresponds to those shown, for example, in [6].
Howewver, as explained in Sec.II, we do not claim to
resolve the small-time behavior accurately.
For the larger values of k shown in Fig.1, the scat-
tering functions exhibit two-stage relaxation. The first
drop, which we identify as the β relaxation, occurs be-
cause the initially mobile molecules are diffusing beyond
their cages. Then, after much longer times indicated by
the plateaus, the scattering curves cross over to typical α
behavior in which all molecules, independent of whether
or not they were initially mobile, make slow transitions
back and forth between mobile and frozen states. This
two-step relaxation disappears at small k, where the scat-
tering functions are averages over distances larger than
R∗ and exhibit normal, diffusive behavior.
The most interesting feature of these scattering func-
tions is that they exhibit a continuous range of stretched
exponential relaxation modes. To see this, we next look
at the long-time behavior of n̂G(k, t
∗) for relatively large
values of t∗ in the range 1 < t∗ < 1000, and for a sequence
of values of k. In Fig.2, we plot − log10 [− log10 n̂G(k, t∗)]
as a function of log10(t
∗), so that the slopes of the curves
are equal to (minus) the stretched-exponential index b.
The results are again for ∆ = 100.
Each of the curves in Fig.2 has a constant slope over
about two decades in t∗, indicating apparently well de-
fined values of b in n̂G ∝ exp (−t∗b). In the limit of large
k, n̂G(k, t
∗) is indistinguishable from the glassy waiting
time distribution φG(t
∗) defined in Eq.(3.1) and shown
by the red curve in the figure; and b is accurately equal to
1/2 for large t∗. b increases toward unity as k decreases
toward values of order unity or smaller. For any nonzero
k, these curves also exhibit a crossover from b > 1/2 at
small t∗ to b = 1/2 at large t∗.
-1 0 1 2 3
FIG. 2: (Color online) Intermediate scattering function
n̂G(k, t
∗) for k = 0.7, 1, 1.5, 5, 10, reading from top to bot-
tom. The red curve is the glassy waiting time distribution
∗), which in this plot is indistinguishable from the large-
k limit of n̂G(k, t
To see this behavior analytically, we can deduce from
Eq.(4.3) that A(w) ≈ 1+ 3w/2 for the very small values
of w that are relevant at very large t∗. Therefore, for
large ∆,
W (k,−w + i0) ≈ k
+ i B(w); (4.5)
and, because B(w) ≪ 1 for w ≪ 1, the integrand in
Eq.(4.4) is sharply peaked at w = k2/3 so long as k2 ≪ 1.
For t∗ ≫ 1, the integrand in Eq.(4.4) has another sharp
peak at the maximum of the function exp (−w t∗)B(w),
i.e. at w = w∗ = 1/
t∗. If k2/3 < 1/
t∗ ≪ 1, then the
diffusive peak is dominant, and n̂G ∼ exp (−k2 t∗/3). On
the other hand, if 1/
t∗ < k2/3 ≪ 1, then the anoma-
lous peak at w∗ is dominant, and a saddle-point esti-
mate yields n̂G ∼ 2
t∗ exp (−2
t∗) as expected from
Eq.(3.1). Thus, at any fixed large time t∗, n̂G becomes
purely diffusive in the limit of small k. Conversely, at
fixed small k, the diffusion becomes anomalous at large
enough t∗.
A complementary set of behaviors is illustrated
in Fig.3, which is the analog of Fig.2 for the
initially mobile molecules. That is, Fig.3 shows
− log10 [− log10 n̂M (k, t∗)] as a function of log10(t∗), but
for a broader range of times extending in two cases
down to t∗ ∼ 10−3. At large times, the preceding
analysis remains valid; for large k, n̂M ∼ φG(t∗) ∼
t∗ exp (−2
t∗). As may be expected, however, and
consistent with the behavior seen in Fig.1, n̂M decays
diffusively at small times (with slope −1 in this graph),
-3 -2 -1 0 1 2 3
FIG. 3: (Color online) Intermediate scattering function
n̂M (k, t
∗) for k = 0.7, 1, 2, 5, reading from top to bottom.
The red curve is the glassy waiting time distribution φG(t
and then enters a flat plateau before crossing over to the
long-time anomalous behavior characteristic of both of
these scattering functions. For large enough k and/or ∆,
the early-time diffusive behavior can be deduced analyti-
cally. In this case, the relevant values of w are large, and
A(w) ≈ −2/w becomes negligibly small. Therefore,
W (k,−w + i0) ≈ 1 +
+ i B(w), (4.6)
where, for large w, B(w) ≈ 2
π/w3/2 again becomes
small. Now the integrand has a sharp peak at w = ∆(1+
k2/2). This peak dominates the integrand at large w for
n̂M , but its amplitude is smaller by a factor ∆
−1 for
n̂G. Thus, for small t
∗, n̂M ∼ exp [−(1 + k2/2)∆ t∗] is
normally diffusive, while no such behavior occurs in n̂G.
V. SPATIAL DISTRIBUTION FUNCTIONS
Yet another view of the normal and anomalous dif-
fusive behaviors is obtained by looking at the spatial
distribution functions themselves, that is, by comput-
ing nG(r
∗, t∗) and nM (r
∗, t∗). The Fourier transforms in
Eqs.(3.9) and (3.10) can be inverted analytically, yielding
∗, t∗) = φG(t
∗) δ(r∗)− 1
4 π r∗
∗, t∗), (5.1)
where
∗, t∗) =
∫ +i∞
2π iu
ψ̃G(u)κ(u) e
−κ(u) r∗ ,
(5.2)
κ(u) =
1− ψ̃G(u) + u/∆
, Reκ ≥ 0. (5.3)
0 2 4 6 8 10
Displacement x*
FIG. 4: Weighted spatial distribution function F̄s(x
∗, t∗) as a
function of the scaled displacement x∗ for scaled times t∗ =
0.03, 0.1, 0.3, 1.0, 3.0, 10.0, reading from left to right along
the bottom of the graph.
Similarly,
nM (r
∗, t∗) = −
4 π r∗
ΓM (r
∗, t∗), (5.4)
where
ΓM (r
∗, t∗) =
∫ +i∞
2π iu
κ(u) e−κ(u) r
. (5.5)
It is convenient to project these distributions onto, say,
the x∗ axis, that is, to integrate out the perpendicular
directions. The formulas analogous to Eqs.(5.1) and (5.4)
for the projected distributions, denoted n̄G(x
∗, t∗) and
n̄M (x
∗, t∗), are
n̄G(x
∗, t∗) = φG(t
∗) δ(x∗) +
∗, t∗), (5.6)
n̄M (x
∗, t∗) =
ΓM (x
∗, t∗), (5.7)
where the functions ΓG and ΓM are the same as those
given in Eqs.(5.2) and (5.5) with r∗ replaced by x∗.
The δ functions in nG(r
∗, t∗) and n̄G(x
∗, t∗) are the
results of our having neglected the intra-cage motion dis-
cussed earlier. In this approximation, with spatial res-
olution only of the order of ℓ, molecules in frozen do-
mains remain exactly at their initial positions with de-
caying probability φG(t
∗). The distributions nM (r
∗, t∗)
and n̄M (x
∗, t∗) have no δ-function contributions because
diffusion starts immediately in mobile regions. In graphs
of experimental or computational data, these δ functions
are visible as narrow Gaussian displacement distributions
centered at x∗ = 0, whose root-mean-square widths, say
ℓ∗C = ℓC/R
∗, are the average intra-cage displacements
that we are neglecting here.
As a first look at the spatial distribution functions,
in Fig.4, we show log10[F̄s(x
∗, t∗)] as a function of the
displacement x∗ for a sequence of times t∗. In analogy to
Eq.(4.1), we write
F̄s(x
∗, t∗) = PG n̄G(x∗, t∗) + (1−PG) n̄M (x∗, t∗), (5.8)
For plotting these graphs, and for this purpose only, we
replace the δ function in Eq.(5.6) by a normalized, Gaus-
sian, intra-cage distribution
δ(x∗) → p̄C(x∗) =
(2 π ℓ∗C
2)1/2
2 ℓ∗C
. (5.9)
We make this replacement primarily because the result-
ing graphs look – and indeed are – more realistic this
way; but we emphasize that this is not a systematic cor-
rection of the small-time behavior because the intra-cage
fluctuations are not otherwise included in n̄G(x
∗, t∗) or
n̄M (x
∗, t∗). Again, we choose ∆ = 100 and PG = 0.5;
and, so that the small-x∗, early-time behavior be visible
in the figure, a relatively large but possibly realistic value
of ℓ∗C = 0.1. Because of the factor φG(t
∗) in Eq.(5.6), this
central peak in n̄G(x
∗, t∗) disappears after times of the
order of τ∗α.
At small times t∗, and for sufficiently large displace-
ments x∗, the graphs in Fig.4 exhibit the broad exponen-
tial tails reported elsewhere in the literature (see [3, 14]).
The exponential tail is a robust mathematical feature of
this class of models, closely associated with the decou-
pling of the mobile and glassy behaviors, and not depen-
dent on any details of the glassy waiting-time distribution
∗). The degree of decoupling is reflected here by the
magnitude of the parameter ∆, which must be large for
strong decoupling. The only requirement on ψG(t
∗) is
that it decays rapidly enough that ψ̃G(u) becomes van-
ishingly small at large u.
To see how the exponential function emerges, evaluate
the u-integrations in Eqs.(5.2) and (5.5) by first integrat-
ing around a circle of radius u0 centered at the origin,
and then closing the contour around the cut on the neg-
ative u-axis from u = −u0 to u → −∞. (The curves
in Fig.4 were computed numerically with u0 = 2.) If
∆ ≫ 1, we can choose 1 ≪ u0 ≪ ∆, so that ψ̃G(u) ≈ 0
and κ(u) ≈
2 everywhere around the circle. Then,
for small t∗, the integration around the circle yields
∗, t∗) ≈ ΓM (x∗, t∗) ≈ exp (−
2 x∗). This limiting
behavior is most apparent in Fig.4 for the two earliest
times, t∗ = 0.03 and 0.1, where the slope has the pre-
dicted value of −
At larger times, and in the limit of indefinitely large ∆,
the preceding argument implies that all of these curves
approach the same slope at large x∗. On the other hand,
at very early times t∗ ≤ ∆−1, the integration is domi-
nated by the discontinuity across the cut near u = −∆,
and the displacement distribution is dominated by the
early diffusive motion of the initially mobile molecules.
Another interesting case is the limit of large t∗ and
x∗ ≪ t∗. The integrands in Eqs.(5.2) and (5.5) vanish
like exp (−1/w) in the limit−u = w → 0+. Therefore, we
can use the long-time, small-w approximation, A(w) ≈
1 + 3w/2, to estimate
∗, t∗) ≈ ΓM (x∗, t∗) ∝ exp
; (5.10)
thus the spatial distribution reverts to a Gaussian.
Eq.(5.10) fits the curve in Fig.4 for t∗ = 10 reasonably
well for x∗ < 5.
VI. MOMENTS OF THE SPATIAL
DISTRIBUTIONS
Another way of looking at these spatial distributions
is to compute their time-dependent moments 〈r∗k(t∗)〉.
We use the identities
〈r∗2(t∗)〉 = −3
∂2Fs(k, t
, (6.1)
〈r∗4(t∗)〉 = 15
∂4Fs(k, t
. (6.2)
When inverting the Laplace transforms as in Eq.(4.2),
we again close the contour of integration in the negative
u plane, but in this case it is mathematically essential to
include the circle around the origin in order to account
for the singularities that occur there. The mathematical
situation becomes clear by inspection of the formula for
〈r∗2(t∗)〉M :
〈r∗2(t∗)〉M = 6
∫ +i∞
κ2(u)
∫ +i∞
1− ψ̃G(u) + u/∆
. (6.3)
Because ψG(u) ≈ 1−3 u/2 for small u, this integrand has
a first-order pole at the origin. Higher-order moments
have higher-order poles.
We can use the analogs of Eq.(6.3) for higher moments
and for both the nG and nM distributions to obtain ex-
act asymptotic results in the limits of very small and
very large t∗. For vanishingly small values of t∗, close
the contour around a circle at large values of u where
ψG(u) ≈ 2/u. In this limit, the integrand behaves like
exp (u t∗)/u2, and
〈r∗2〉M ≈ 3∆ t∗. (6.4)
Similar small-t∗ calculations yield:
〈r∗4〉M ≈ 15∆2 t∗2; (6.5)
-6 -5 -4 -3 -2 -1 0 1 2
FIG. 5: (Color online) log
(t∗)〉F ] as a function of
(t∗) for PG = 0, 0.5, 0.75, 0.9, 0.97, 0.999, and 1.0 (red
curve), reading from top to bottom.
〈r∗2〉G ≈ 3∆ t∗2; (6.6)
〈r∗4〉G ≈ 10∆2 t∗3. (6.7)
For very large t∗, close the contour in Eq.(6.3) on
a vanishingly small circle around the origin and use
ψG(u) ≈ 1−3 u/2. The result is that, after times so long
that the tagged molecule no longer remembers whether it
started in a mobile or a glassy region, both 〈r∗2〉M and
〈r∗2〉G converge to the same slowly diffusing Gaussian
distribution for which 〈r∗2〉 ≈ 2 t∗ (with corrections of
the order of 1/∆).
Note several features of these results. As expected, the
initially mobile molecules exhibit rapid (large ∆) Gaus-
sian diffusion at small times. On the other hand, the
displacements of the initially frozen molecules are non-
Gaussian and “pseudo-ballistic” with 〈r∗2〉G ∼ ∆ t∗2.
This behavior has nothing to do with early-stage bal-
listic motion of molecules within their cages but, rather,
is an intrinsic feature of intermediate-stage, anomalous
diffusion in this model.
We turn finally to the full time dependence of the mo-
ments. The results for the weighted average,
〈r∗2(t∗)〉F = PG 〈r∗2(t∗)〉G+(1−PG) 〈r∗2(t∗)〉M , (6.8)
are shown in Fig.5 for a series of different values of PG,
as functions of log10(t
∗). At very short times, these mo-
ments rise as functions of t∗ according to the estimates
in Eqs.(6.4) - (6.7). The pseudo-ballistic behavior is ap-
parent only for PG ≈ 1 because, for smaller PG, these
weighted moments are dominated by the displacements
of the initially mobile molecules. After the initial rise, at
roughly t∗ ≈ ∆−1 = 10−2, the moments cross over from
-3 -2 -1 0 1 2
FIG. 6: Non-Gaussian parameter α2 as a function of log10(t
for PG = 0.
fast β relaxation to slow α decay as the tagged molecules
are repeatedly trapped and then escape from glassy do-
mains. Finally, in Fig.6, we show the non-Gaussian pa-
rameter
3 〈r∗4(t∗)〉M
5 〈r∗2(t∗)〉2M
− 1 (6.9)
only for PG = 0 because, in this approximation where
we have neglected intra-cage vibrational motions, only
the mobile molecules exhibit Gaussian displacement dis-
tributions at early times. Here we see explicitly that the
non-Gaussian behavior occurs during the crossover from
β to α relaxation. Both of these last two figures are ar-
tificial in the sense that we are independently varying
the parameter PG, which in fact ought to be determined
uniquely by the temperature along with the parameters
R∗ and ∆.
VII. CONCLUDING REMARKS
The model proposed here contains a relatively simple
mechanism for producing stretched exponential decay of
molecular correlations, as observed via the self interme-
diate scattering function. That mechanism emerges di-
rectly from the spatial heterogeneity of a glass-forming
liquid. The model illustrates how anomalous diffusion,
as exemplified by a broad exponential (non-Gaussian)
tail of the molecular displacement distribution, is related
to heterogeneity and – indirectly – to the accompanying
stretched exponential behavior. Yet another feature of
the model is the way in which it illustrates how the tran-
sition from relatively fast β relaxation to slow α decay
correlates with the onset of anomalous diffusion. Our
analysis neglects the very short time, strongly localized,
intra-cage fluctuations. The transition between the latter
motions and cage-breaking events seems to us to be out-
side the range of validity of the continuous-time random-
walk theory used here.
There are, of course, many missing ingredients, almost
all of which pertain to the temperature dependence of
the various parameters introduced here. One of us (JSL)
plans to address those issues in a subsequent report.
Acknowledgments
This research was supported by U.S. Department of
Energy Grant No. DE-FG03-99ER45762.
[1] M.D. Ediger, Ann. Rev. Phys. Chem. 51, 99 (2000)
[2] L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, D.
El Masri, D.L’Hote, F. Ladieu, and M. Pierno, Science
310, 1797 (2005).
[3] E. R. Weeks et al, Science 287, 627 (2000)
[4] W. Kob et al, Phys. Rev. Lett. 79, 2827 (1997)
[5] R. Yamamoto and A. Onuki, Phys. Rev. E 58, 3515
(1998)
[6] W. Kob and H. C. Andersen, Phys. Rev. Lett. 73, 1376
(1994); Phys. Rev E 51, 4626 (1995); Phys. Rev E 52,
4134 (1995); W. Kob in “Slow relaxation and nonequi-
librium dynamics in condensed matter,” Les Houches ses-
sion LXXVII (2002)
[7] J. S. Langer, Phys. Rev. E 73, 041504 (2006)
[8] J. S. Langer, Phys. Rev. Lett. 97, 115704 (2006).
[9] A. Widmer-Cooper and P. Harrowell, Phys. Rev. Lett.
96, 185701 (2006); J. Chem. Phys. 126, 154503 (2007).
[10] L.E. Silbert, A.J. Liu, and S.R. Nagel, Phys. Rev. Lett.
95, 098301 (2005).
[11] Y. Jung, J. Garrahan, and D. Chandler, Phys. Rev. E
69, 061205 (2004).
[12] L. Berthier, D. Chandler, and J. Garrahan, Europhysics
Lett. 69 (3), 320 (2005).
[13] P. Hurtado, L. Berthier, and W. Kob, Phys. Rev. Lett.
98, 135503 (2007).
[14] P. Chaudhuri, L. Berthier, and W. Kob, Phys. Rev. Lett.
99, 060604 (2007).
[15] P. Grassberger and I. Procaccia, J. Chem. Phys. 77, 6281
(1982).
[16] C. Monthus and J.-P. Bouchaud, J. Phys. A: Math. Gen.
29 3847 (1996).
[17] J. Bendler, J. Fontanella and M. Shlesinger, Chem. Phys.
284, 311 (2002).
[18] D. Kivelson, G. Tarjus, X. Zhao, and S.A. Kivelson,
Phys. Rev. E 53(1) 751 (1996).
[19] E. W. Montroll and M. F. Shlesinger in Studies in Sta-
tistical Mechanics XI, editors J. L. Lebowitz and E. W.
Montroll, North-Holland Physics Publishing (1984)
[20] J.-P. Bouchaud and A. Georges, Physics Reports 195,
127 (North Holland, 1990).

  We propose a model of a heterogeneous glass forming liquid and compute the
low-temperature behavior of a tagged molecule moving within it. This model
exhibits stretched-exponential decay of the wavenumber-dependent, self
intermediate scattering function in the limit of long times. At temperatures
close to the glass transition, where the heterogeneities are much larger in
extent than the molecular spacing, the time dependence of the scattering
function crosses over from stretched-exponential decay with an index $b=1/2$ at
large wave numbers to normal, diffusive behavior with $b = 1$ at small
wavenumbers. There is a clear separation between early-stage, cage-breaking
$\beta$ relaxation and late-stage $\alpha$ relaxation. The spatial
representation of the scattering function exhibits an anomalously broad
exponential (non-Gaussian) tail for sufficiently large values of the molecular
displacement at all finite times.

Anomalous diffusion and stretched exponentials in heterogeneous glass-forming
liquids: Low-temperature behavior
J. S. Langer and Swagatam Mukhopadhyay
Dept. of Physics, University of California, Santa Barbara, CA 93106-9530
(Dated: November 5, 2018)
We propose a model of a heterogeneous glass forming liquid and compute the low-temperature
behavior of a tagged molecule moving within it. This model exhibits stretched-exponential decay
of the wavenumber-dependent, self intermediate scattering function in the limit of long times. At
temperatures close to the glass transition, where the heterogeneities are much larger in extent than
the molecular spacing, the time dependence of the scattering function crosses over from stretched-
exponential decay with an index b = 1/2 at large wave numbers to normal, diffusive behavior with
b = 1 at small wavenumbers. There is a clear separation between early-stage, cage-breaking β
relaxation and late-stage α relaxation. The spatial representation of the scattering function exhibits
an anomalously broad exponential (non-Gaussian) tail for sufficiently large values of the molecular
displacement at all finite times.
I. INTRODUCTION
A growing body of evidence, both experimental [1, 2, 3]
and numerical [4, 5], points to the existence of intrinsic
spatial heterogeneities in glass forming liquids at tem-
peratures slightly above the glass transition, and the rel-
evance of these heterogeneities to stretched exponential
decay of correlations and anomalous, non-Gaussian diffu-
sion. Here we propose a simple model of a heterogeneous
glass former and use it to compute observable proper-
ties of a tagged molecule moving within it. In particular,
we show that this model naturally predicts stretched-
exponential decay when the scale of the heterogeneity is
much larger than the molecular spacing; and we compute
deviations from Gaussian displacement distributions dur-
ing intermediate stages of the relaxation process.
Our model is loosely motivated by the excitation-chain
(XC) theory of the glass transition proposed recently by
one of us (JSL) [7, 8] This theory suggests that a glass-
forming liquid at a temperature T not too far above the
glass transition temperature T0 consists of fluctuating do-
mains of linear size R∗(T ), within which the molecules
are frozen in a glassy state where they have little or no
mobility. R∗(T ) is a length scale that characterizes dy-
namic fluctuations in the equilibrium states of these sys-
tems. The theory predicts that R∗(T ) diverges near T0
like (T − T0)−1 and decreases to zero at the upper limit
of the super-Arrhenius region, T = TA.
Characterizing the more mobile material that lies out-
side the hypothetical frozen domains remains one of the
deeper challenges in glass physics. In the language intro-
duced by Widmer-Cooper and Harrowell [9], the mobile
molecules reside in regions of high “propensity.” Those
authors perform molecular dynamics simulations of a
glasslike, two-dimensional, binary mixture. They average
the mean-square displacement of each molecule over an
isoconfigurational ensemble in which the initial positions
of all the molecules are fixed, but the initial velocities are
selected at random from a Boltzmann distribution. The
propensity of a molecule is its mean-square displacement
during a structural relaxation time. Contour maps of
propensity as a function of initial molecular positions do
indeed exhibit well defined domains of low and high mo-
bility.
Widmer-Cooper and Harrowell compare their propen-
sity maps with corresponding plots of the potential en-
ergy, which is taken to be indicative of the local struc-
ture. Interestingly, the propensity appears to be almost
entirely uncorrelated with the potential energy. In retro-
spect, this lack of correlation between short-range struc-
ture and mobility may not be surprising. Such a corre-
lation would be expected for a slowly coarsening system
consisting of domains of two structurally distinct phases,
for example, a slightly supersaturated (mobile) fluid in
which (frozen) crystallites are growing. In that case, the
domains are very nearly but not quite in thermodynamic
coexistence with each other. A true glass-forming liquid,
however, is in a state of thermodynamic equilibrium. The
temperatures and chemical potentials are uniform every-
where, independent of propensity; and the local pair cor-
relations, which ought to be determined primarily by the
temperature in an equilibrium state, appear to be statis-
tically uniform as well.
A local quantity that does correlate with the propen-
sity is the mean-square vibration amplitude of a molecule
averaged over times much longer than its period of oscil-
lation but much shorter than the time required for irre-
versible rearrangements to occur. A large mean-square
vibration amplitude – or, equivalently, a small Debye-
Waller factor – implies that a molecule is participat-
ing in a soft, low frequency, elastic mode. It is well
known that soft modes are abundant in amorphous mate-
rials near jamming transitions.[10] Moreover, it appears
(for example, in Fig.4 of [10]) that the largest displace-
ments in such modes lie predominantly along one dimen-
sional, chainlike paths. Thus, regions of high propensity
may be elastically soft, and may somehow be correlated
with excitation-chain activity. Of course, local variations
of elastic stiffness are ultimately structural in nature;
but soft modes are collective phenomena involving many
molecules, and are not easily detected by measuring near-
http://arxiv.org/abs/0704.1508v2
neighbor pair correlations.
In the model proposed here, glassy domains of low
propensity are separated from each other by regions of
higher propensity in which many of the molecules are
mobile. This domain structure undergoes persistent fluc-
tuations on time scales of the order of τ∗α = R
∗2/Dα,
where Dα = ℓ
2/τα is the diffusion constant associ-
ated with α relaxation, ℓ is approximately the average
molecular spacing, and τα is the strongly temperature-
dependent, super-Arrhenius, α relaxation time. In con-
trast to the slow fluctuations of domain boundaries, indi-
vidual molecules within high-propensity regions are ran-
dom walkers with a diffusion constant DM = ℓ
2/τM ,
where τM is a temperature-dependent time scale – pre-
sumably not super-Arrhenius – that characterizes mobile
molecular displacements.
To a first approximation, a tagged molecule at any
given time is either frozen in a glassy domain or is mobile.
If the former, then it waits a time of the order of τ∗α before
it is encountered by mobile neighbors and itself starts
to undergo displacements. If the latter, then it diffuses
for some distance until it finds itself frozen again in a
glassy domain. That distance is at least a few molecular
spacings ℓ; but it seems more likely to be proportional to
R∗ if the domains are large and if diffusive hopping takes
place more readily within the mobile regions than across
their boundaries and into the glassy domains. Therefore,
we assume that the time elapsed during mobile motion
is of the order of R∗2/DM = (R
∗/ℓ)2 τM . The ratio of
these two time scales is τα/τM ≡ ∆. At temperatures
well below TA, we expect that ∆ ≫ 1.
A similar picture of diffusion in fluctuating regions of
varying mobility emerges in kinetically constrained mod-
els [11, 12] and in a recent model of a gel [13]. Our
analysis is most closely related to the diffusion model of
Chaudhuri et al [14], which has been a valuable start-
ing point for the present investigation. We go beyond
[14] by including a mechanism for producing late-stage
stretched-exponential relaxation, and by distinguishing
that mechanism from the diffusion of mobile molecules.
In contrast to [14], our molecules switch back and forth
between frozen and mobile behavior and thus explore the
geometry of the domain structure. As a result, our model
exhibits a clear separation between the slow α relaxation
and the faster mobile motions that we interpret as β re-
laxation. It also respects time-translational symmetry,
which seems to be violated in [14].
The preceding discussion pertains just to the case of
large-scale heterogeneities at temperatures only a little
above T0 and well below TA. Some of the most interesting
and experimentally accessible physics, however, occurs
near the crossover to small-scale heterogeneity near TA.
Constructing a theory of diffusion in the crossover region
requires a temperature-dependent analysis, which will be
described in a following paper.
II. BASIC INGREDIENTS OF THE MODEL IN
THE LIMIT OF LARGE-SCALE
HETEROGENEITY
As outlined above, there are two different mechanisms
to be described probabilistically in this sytem, each op-
erating on its own characteristic time scale. The separa-
tion between time scales, characterized by a large value
of ∆, suggests that the natural mathematical language in
which to discuss this system is that of a continuous time
random walk (CTRW) in which a tagged molecule al-
ternates between long waits in glassy domains and faster
displacements in mobile regions. In this Section, we com-
pute the waiting time distributions for molecules in glassy
and mobile domains, and the corresponding probability
distributions for the mobile displacements.
We immediately encounter a problem, however, be-
cause there are more than just two characteristic time
scales in this system. Analog experiments [3] and numer-
ical simulations [4] typically follow the motions of parti-
cles starting from their initial positions, and resolve these
trajectories on length scales smaller than the interparti-
cle spacing, and on time scales of the order of the shortest
vibrational periods. This is a continuous range of scales.
At its long-time limit, it includes the slow modes that
are involved in molecular rearrangements, whose ampli-
tudes may exceed the Lindemann melting criterion in the
mobile regions. Thus, the cage-breaking, β-relaxation
mechanism may be part of the same continuous range of
time scales that includes the intra-cage vibrational mo-
tions. Indeed, no sharp distinction between vibrational
and cage breaking time scales seems to appear in exper-
imental data. If these are not clearly distinct kinds of
events, then the CTRW approximation is not valid in
the short-time limit.
In view of this difficulty, we choose here to consider
only the β and α time scales, τM and τα respectively. By
neglecting the shorter time scales, i.e. those of the order
of τ0 ∼ femtoseconds - picoseconds, we imply that we are
not resolving length scales much smaller than the molec-
ular spacing ℓ. We then assume that the probability that
a molecule has left its cage at times of order τM ≫ τ0 is
determined by the probability that it has mobile neigh-
bors and participates in local molecular rearrangements.
At the longest time scales τα ∼ seconds (by definition,
the viscous relaxation time at the glass temperature),
the probability that a molecule has escaped from a large
glassy domain is determined by the probability that it
has been encountered by the boundary of the domain.
Consider first the slowest motions, i.e. those associated
with fluctuations of the domain boundaries on length
scales R∗ and time scales τ∗α. Define t
∗ ≡ t/τ∗α, where
t is the physical time in seconds; and let ψG(t
∗) be the
normalized probability distribution for the time that a
molecule spends in a glassy domain before entering a mo-
bile region. Write this distribution in the form
dρW (ρ)
, (2.1)
where ρ is the linear size of a domain in units R∗. W (ρ)
is a normalized distribution over these sizes, and the re-
maining factor inside the integrand is a normalized dis-
tribution over t∗. The quantity ρ−2 appearing in the
exponential in Eq.(2.1) is the lowest eigenvalue of the
diffusion kernel for a molecule moving in a domain of
size ρ. Note the similarity to the trapping models dis-
cussed, for example, in [15, 16, 17]; but also note that
the “trap” here is a two dimensional subspace bounding
a three dimensional domain. The quantity that is fluc-
tuating in a normal diffusive manner is the mean square
of the distance between the tagged molecule and the do-
main boundary. Because the boundary surrounds the
molecule, the direction of the diffusive drift is irrelevant;
any point of contact on the boundary is equivalent to any
other. Therefore, it makes little difference whether the
molecule is the diffuser and the boundary is the target, or
– as in this case – the boundary diffuses and the molecule
is the target.
Strictly speaking, Eq.(2.1) is a long-time approxima-
tion; at shorter times, the higher eigenmodes of the diffu-
sion kernel make non-negligible contributions. However,
this approximation is qualitatively adequate for our pur-
poses at shorter times as well because ψG(t
∗) is well be-
haved at small t∗ and, as mentioned previously, we are
not trying to include the short-time behavior in this func-
tion. If we assume that the distribution over values of ρ
is Gaussian, W (ρ) ∝ ρ2 exp(−ρ2) in three dimensions,
∗) = 2 e−2
t∗ . (2.2)
That is, we find a rudimentary but nontrivial stretched
exponential of the form exp (−const.× t∗b) with b = 1/2.
There is little reason to expect that W (ρ) remains
Gaussian out to the large values of ρ that determine
the long-time behavior of ψG(t
∗). Deviations from the
Gaussian would produce different indices. For example,
as the temperature increases and the glassy domains be-
come smaller, their size distributions might cut off more
sharply than a Gaussian. Thus, if W (ρ) ∝ exp(−ρm),
with m ≥ 2, then b = m/(m + 2) → 1 as m → ∞. The
problem of computing this distribution or, equivalently,
values of m from first principles may eventually become
solvable as we learn more about the statistical mechan-
ics of glass forming liquids; but that problem is beyond
the scope of this investigation. For present purposes, the
important points are that a plausible distribution W (ρ)
produces stretched-exponential behavior, and that the
mechanism by which this happens could produce a tem-
perature dependent index b.
Next, consider the waiting-time distribution ψM (t
for a molecule in a mobile region. This situation is qual-
itatively different from that of a molecule in a glassy do-
main because the distance traveled by the molecule is
related to the time during which it remains mobile. As
stated earlier, we assume that a mobile molecule diffuses
a distance of order R∗ before reentering a glassy region,
so that its residence time in the mobile region is of order
τM (R
∗/ℓ)2 which, in t∗ units, is simply ∆−1. Then, for
simplicity, assume an exponential waiting-time distribu-
tion:
ψM (t
∗) = ∆ e−∆ t
. (2.3)
Compared to ψG(t
∗) in Eq.(2.2), ψM (t
∗) is sharply
peaked near t∗ = 0 if, as expected, ∆ is large. To com-
plete the model of mobile motion, we need the conditional
probability pM (r
∗, t∗) for diffusion over a scaled distance
r∗ = r/R∗ in time t∗. The natural choice is
pM (r
∗, t∗) =
(2 π∆ t∗)3/2
2∆ t∗
, (2.4)
which is a normalized, three dimensional distribution
over r∗.
III. CONTINUOUS-TIME RANDOM WALKS
The next step is to translate these physical ingredients
of the model into the language of continuous-time ran-
dom walks.[19, 20] Define two different probability distri-
bution functions, nG(r
∗, t∗) and nM (r
∗, t∗), for molecules
starting, respectively, in glassy domains or mobile re-
gions, and moving distances r∗ in times t∗. Each molec-
ular trajectory consists of a sequence of transitions be-
tween glassy domains and mobile regions. A single tran-
sition in which a molecule starts in a glassy domain at
time t∗1 and ends in a mobile region at time t
2 (without
having changed its actual position) occurs with probabil-
ity ψG(t
2−t∗1). Similarly, a transition in which a molecule
starts in a mobile region at time t∗1 and position r
1, and
ends in a glassy domain at time t∗2 and position r
2, has
probability ψM (t
2 − t∗1) pM (|r∗2 − r∗1|, t∗2 − t∗1). The prob-
ability that a molecule starts or arrives in a glassy do-
main at time t∗1 and is still there at the final time t
∗ − t∗1), where
∗′) dt∗
= (1 + 2
t∗) e−2
t∗ . (3.1)
Finally, the probability that a molecule starts or arrives
in a mobile region at time t∗1 and position r
1 and is still
in that region at the final time t∗ and at position r∗, is
φM (t
∗ − t∗1) pM (|r∗ − r∗1|, t∗ − t∗1), where, in analogy to
Eq.(3.1),
φM (t
ψM (t
1) dt
1 = e
−∆ t∗ . (3.2)
The multiple convolutions that describe each trajec-
tory are, as usual, converted into products by computing
Fourier-Laplace transforms of each function:
ñj(k, u) =
e−ik·r
∗, t∗) dt∗, (3.3)
where j = G,M . For jumps starting in glassy regions,
we need:
ψ̃G(u) =
∗) dt∗, (3.4)
φ̃G(u) =
1− ψ̃G(u)
. (3.5)
For the mobile regions, the spatial Fourier transform of
the conditional probability distribution pM (r
∗, t∗) is
p̂M (k, t
∗) = exp
; (3.6)
therefore, define
fM (k, u) ≡
p̂M (k, t
∗)ψM (t
∗) dt∗
u+∆(1 + k2/2)
, (3.7)
gM (k, u) ≡
p̂M (k, t
∗)φM (t
∗) dt
u+∆(1 + k2/2)
. (3.8)
Putting these pieces together, and summing over in-
definitely many individual jumps in each trajectory, we
ñG(k, u) =
φ̃G(u) + gM (k, u) ψ̃G(u)
1− ψ̃G(u) fM (k, u)
NG(k, u)
W (k, u)
; (3.9)
ñM (k, u) =
φ̃G(u) fM (k, u) + gM (k, u)
1− ψ̃G(u) fM (k, u)
NM (k, u)
W (k, u)
; (3.10)
where
NG(k, u) =
1− ψ̃G(u)
(1 + k2/2) + u/∆; (3.11)
NM (k, u) = 1− ψ̃G(u) + u/∆; (3.12)
W (k, u) = 1− ψ̃G(u) + k2/2 + u/∆. (3.13)
In both Eqs.(3.9) and (3.10), in the numerators of the
first expressions, the first term corresponds to trajecto-
ries that end in glassy domains, and the second to those
that end in mobile regions.
-6 -4 -2 0 2
FIG. 1: Intermediate scattering function Fs(k, t
∗) with PG =
0.5 and ∆ = 100, for k = 0.8, 1, 1.2, 1.5, 2, 3, 5, 10, reading
from top to bottom.
IV. INTERMEDIATE SCATTERING
FUNCTION
The conventional way to study the diffusion mecha-
nisms discussed here is to measure the self intermediate
scattering function Fs(k, t), which, in the present no-
tation, is the mixed Fourier-time representation of the
weighted average of nG and nM . That is,
Fs(k, t
∗) = PG n̂G(k, t∗) + (1− PG) n̂M (k, t∗), (4.1)
where PG is the probability that a molecule starts its
motion in a frozen, glassy environment; and
n̂j(k, t
∫ +i∞
ñj(k, u), j = G, M. (4.2)
The inversions of Laplace-transforms in Eq.(4.2) are
nontrivial because the stretched-exponential function
ψ̃G(u), defined in Eq.(3.4), has an essential singularity
at u = 0. For t∗ > 0, the integrations over u in Eq.(4.2)
must be performed by closing the contour in the negative
half u-plane, which requires that ñj(k, u) be analytically
continued to points u = |u| exp(iθ). To evaluate ψ̃G(u)
at such points, using Eq.(3.4), let t∗ = y2 and then ro-
tate the y contour away from the positive real axis to
a line y = ξ exp(−iθ/2), 0 < ξ < ∞ so that the inte-
grand always decreases rapidly at infinity. For example,
for u→ −w ± iǫ, θ = ±π,
ψ̃G(−w ± iǫ) = −4
ξ dξ e−wξ
2±2iξ
e−1/w
≡ A(w) ∓ i B(w), (4.3)
where Erfi is the imaginary error function. The imagi-
nary part of ψ̃G(−w ± iǫ) is the discontinuity across a
cut along the negative u axis. With this formula, it is
straightforward to compute the discontinuity across the
cut for the ñj(k, u) in Eqs.(3.9) and (3.10), and in this
way to compute the scattering functions by closing the
contour around this cut. That is:
n̂j(k, t
e−w t
Nj(k,−w + i0)
W (k,−w + i0)
, j = G, M. (4.4)
The integrands are well behaved as w → 0+, therefore it
is possible to let the lower limit of integration be w = 0.
The graphs in Fig.1 show Fs(k, t
∗), computed nu-
merically from Eq.(4.4), as a function of log10(t
∗), for
∆ = 100, PG = 0.5, and for a sequence of wavenumbers
k. These are low-temperature scattering functions, where
∆ ≫ 1, which means that molecules spend much longer
times frozen in glassy domains than they do when moving
in mobile regions. The choice PG = 0.5 is made primar-
ily for clarity; but it seems likely that this value of PG
describes a system that is well into its non-Arrhenius,
anomalous-diffusion regime. The set of curves in Fig.1
closely corresponds to those shown, for example, in [6].
Howewver, as explained in Sec.II, we do not claim to
resolve the small-time behavior accurately.
For the larger values of k shown in Fig.1, the scat-
tering functions exhibit two-stage relaxation. The first
drop, which we identify as the β relaxation, occurs be-
cause the initially mobile molecules are diffusing beyond
their cages. Then, after much longer times indicated by
the plateaus, the scattering curves cross over to typical α
behavior in which all molecules, independent of whether
or not they were initially mobile, make slow transitions
back and forth between mobile and frozen states. This
two-step relaxation disappears at small k, where the scat-
tering functions are averages over distances larger than
R∗ and exhibit normal, diffusive behavior.
The most interesting feature of these scattering func-
tions is that they exhibit a continuous range of stretched
exponential relaxation modes. To see this, we next look
at the long-time behavior of n̂G(k, t
∗) for relatively large
values of t∗ in the range 1 < t∗ < 1000, and for a sequence
of values of k. In Fig.2, we plot − log10 [− log10 n̂G(k, t∗)]
as a function of log10(t
∗), so that the slopes of the curves
are equal to (minus) the stretched-exponential index b.
The results are again for ∆ = 100.
Each of the curves in Fig.2 has a constant slope over
about two decades in t∗, indicating apparently well de-
fined values of b in n̂G ∝ exp (−t∗b). In the limit of large
k, n̂G(k, t
∗) is indistinguishable from the glassy waiting
time distribution φG(t
∗) defined in Eq.(3.1) and shown
by the red curve in the figure; and b is accurately equal to
1/2 for large t∗. b increases toward unity as k decreases
toward values of order unity or smaller. For any nonzero
k, these curves also exhibit a crossover from b > 1/2 at
small t∗ to b = 1/2 at large t∗.
-1 0 1 2 3
FIG. 2: (Color online) Intermediate scattering function
n̂G(k, t
∗) for k = 0.7, 1, 1.5, 5, 10, reading from top to bot-
tom. The red curve is the glassy waiting time distribution
∗), which in this plot is indistinguishable from the large-
k limit of n̂G(k, t
To see this behavior analytically, we can deduce from
Eq.(4.3) that A(w) ≈ 1+ 3w/2 for the very small values
of w that are relevant at very large t∗. Therefore, for
large ∆,
W (k,−w + i0) ≈ k
+ i B(w); (4.5)
and, because B(w) ≪ 1 for w ≪ 1, the integrand in
Eq.(4.4) is sharply peaked at w = k2/3 so long as k2 ≪ 1.
For t∗ ≫ 1, the integrand in Eq.(4.4) has another sharp
peak at the maximum of the function exp (−w t∗)B(w),
i.e. at w = w∗ = 1/
t∗. If k2/3 < 1/
t∗ ≪ 1, then the
diffusive peak is dominant, and n̂G ∼ exp (−k2 t∗/3). On
the other hand, if 1/
t∗ < k2/3 ≪ 1, then the anoma-
lous peak at w∗ is dominant, and a saddle-point esti-
mate yields n̂G ∼ 2
t∗ exp (−2
t∗) as expected from
Eq.(3.1). Thus, at any fixed large time t∗, n̂G becomes
purely diffusive in the limit of small k. Conversely, at
fixed small k, the diffusion becomes anomalous at large
enough t∗.
A complementary set of behaviors is illustrated
in Fig.3, which is the analog of Fig.2 for the
initially mobile molecules. That is, Fig.3 shows
− log10 [− log10 n̂M (k, t∗)] as a function of log10(t∗), but
for a broader range of times extending in two cases
down to t∗ ∼ 10−3. At large times, the preceding
analysis remains valid; for large k, n̂M ∼ φG(t∗) ∼
t∗ exp (−2
t∗). As may be expected, however, and
consistent with the behavior seen in Fig.1, n̂M decays
diffusively at small times (with slope −1 in this graph),
-3 -2 -1 0 1 2 3
FIG. 3: (Color online) Intermediate scattering function
n̂M (k, t
∗) for k = 0.7, 1, 2, 5, reading from top to bottom.
The red curve is the glassy waiting time distribution φG(t
and then enters a flat plateau before crossing over to the
long-time anomalous behavior characteristic of both of
these scattering functions. For large enough k and/or ∆,
the early-time diffusive behavior can be deduced analyti-
cally. In this case, the relevant values of w are large, and
A(w) ≈ −2/w becomes negligibly small. Therefore,
W (k,−w + i0) ≈ 1 +
+ i B(w), (4.6)
where, for large w, B(w) ≈ 2
π/w3/2 again becomes
small. Now the integrand has a sharp peak at w = ∆(1+
k2/2). This peak dominates the integrand at large w for
n̂M , but its amplitude is smaller by a factor ∆
−1 for
n̂G. Thus, for small t
∗, n̂M ∼ exp [−(1 + k2/2)∆ t∗] is
normally diffusive, while no such behavior occurs in n̂G.
V. SPATIAL DISTRIBUTION FUNCTIONS
Yet another view of the normal and anomalous dif-
fusive behaviors is obtained by looking at the spatial
distribution functions themselves, that is, by comput-
ing nG(r
∗, t∗) and nM (r
∗, t∗). The Fourier transforms in
Eqs.(3.9) and (3.10) can be inverted analytically, yielding
∗, t∗) = φG(t
∗) δ(r∗)− 1
4 π r∗
∗, t∗), (5.1)
where
∗, t∗) =
∫ +i∞
2π iu
ψ̃G(u)κ(u) e
−κ(u) r∗ ,
(5.2)
κ(u) =
1− ψ̃G(u) + u/∆
, Reκ ≥ 0. (5.3)
0 2 4 6 8 10
Displacement x*
FIG. 4: Weighted spatial distribution function F̄s(x
∗, t∗) as a
function of the scaled displacement x∗ for scaled times t∗ =
0.03, 0.1, 0.3, 1.0, 3.0, 10.0, reading from left to right along
the bottom of the graph.
Similarly,
nM (r
∗, t∗) = −
4 π r∗
ΓM (r
∗, t∗), (5.4)
where
ΓM (r
∗, t∗) =
∫ +i∞
2π iu
κ(u) e−κ(u) r
. (5.5)
It is convenient to project these distributions onto, say,
the x∗ axis, that is, to integrate out the perpendicular
directions. The formulas analogous to Eqs.(5.1) and (5.4)
for the projected distributions, denoted n̄G(x
∗, t∗) and
n̄M (x
∗, t∗), are
n̄G(x
∗, t∗) = φG(t
∗) δ(x∗) +
∗, t∗), (5.6)
n̄M (x
∗, t∗) =
ΓM (x
∗, t∗), (5.7)
where the functions ΓG and ΓM are the same as those
given in Eqs.(5.2) and (5.5) with r∗ replaced by x∗.
The δ functions in nG(r
∗, t∗) and n̄G(x
∗, t∗) are the
results of our having neglected the intra-cage motion dis-
cussed earlier. In this approximation, with spatial res-
olution only of the order of ℓ, molecules in frozen do-
mains remain exactly at their initial positions with de-
caying probability φG(t
∗). The distributions nM (r
∗, t∗)
and n̄M (x
∗, t∗) have no δ-function contributions because
diffusion starts immediately in mobile regions. In graphs
of experimental or computational data, these δ functions
are visible as narrow Gaussian displacement distributions
centered at x∗ = 0, whose root-mean-square widths, say
ℓ∗C = ℓC/R
∗, are the average intra-cage displacements
that we are neglecting here.
As a first look at the spatial distribution functions,
in Fig.4, we show log10[F̄s(x
∗, t∗)] as a function of the
displacement x∗ for a sequence of times t∗. In analogy to
Eq.(4.1), we write
F̄s(x
∗, t∗) = PG n̄G(x∗, t∗) + (1−PG) n̄M (x∗, t∗), (5.8)
For plotting these graphs, and for this purpose only, we
replace the δ function in Eq.(5.6) by a normalized, Gaus-
sian, intra-cage distribution
δ(x∗) → p̄C(x∗) =
(2 π ℓ∗C
2)1/2
2 ℓ∗C
. (5.9)
We make this replacement primarily because the result-
ing graphs look – and indeed are – more realistic this
way; but we emphasize that this is not a systematic cor-
rection of the small-time behavior because the intra-cage
fluctuations are not otherwise included in n̄G(x
∗, t∗) or
n̄M (x
∗, t∗). Again, we choose ∆ = 100 and PG = 0.5;
and, so that the small-x∗, early-time behavior be visible
in the figure, a relatively large but possibly realistic value
of ℓ∗C = 0.1. Because of the factor φG(t
∗) in Eq.(5.6), this
central peak in n̄G(x
∗, t∗) disappears after times of the
order of τ∗α.
At small times t∗, and for sufficiently large displace-
ments x∗, the graphs in Fig.4 exhibit the broad exponen-
tial tails reported elsewhere in the literature (see [3, 14]).
The exponential tail is a robust mathematical feature of
this class of models, closely associated with the decou-
pling of the mobile and glassy behaviors, and not depen-
dent on any details of the glassy waiting-time distribution
∗). The degree of decoupling is reflected here by the
magnitude of the parameter ∆, which must be large for
strong decoupling. The only requirement on ψG(t
∗) is
that it decays rapidly enough that ψ̃G(u) becomes van-
ishingly small at large u.
To see how the exponential function emerges, evaluate
the u-integrations in Eqs.(5.2) and (5.5) by first integrat-
ing around a circle of radius u0 centered at the origin,
and then closing the contour around the cut on the neg-
ative u-axis from u = −u0 to u → −∞. (The curves
in Fig.4 were computed numerically with u0 = 2.) If
∆ ≫ 1, we can choose 1 ≪ u0 ≪ ∆, so that ψ̃G(u) ≈ 0
and κ(u) ≈
2 everywhere around the circle. Then,
for small t∗, the integration around the circle yields
∗, t∗) ≈ ΓM (x∗, t∗) ≈ exp (−
2 x∗). This limiting
behavior is most apparent in Fig.4 for the two earliest
times, t∗ = 0.03 and 0.1, where the slope has the pre-
dicted value of −
At larger times, and in the limit of indefinitely large ∆,
the preceding argument implies that all of these curves
approach the same slope at large x∗. On the other hand,
at very early times t∗ ≤ ∆−1, the integration is domi-
nated by the discontinuity across the cut near u = −∆,
and the displacement distribution is dominated by the
early diffusive motion of the initially mobile molecules.
Another interesting case is the limit of large t∗ and
x∗ ≪ t∗. The integrands in Eqs.(5.2) and (5.5) vanish
like exp (−1/w) in the limit−u = w → 0+. Therefore, we
can use the long-time, small-w approximation, A(w) ≈
1 + 3w/2, to estimate
∗, t∗) ≈ ΓM (x∗, t∗) ∝ exp
; (5.10)
thus the spatial distribution reverts to a Gaussian.
Eq.(5.10) fits the curve in Fig.4 for t∗ = 10 reasonably
well for x∗ < 5.
VI. MOMENTS OF THE SPATIAL
DISTRIBUTIONS
Another way of looking at these spatial distributions
is to compute their time-dependent moments 〈r∗k(t∗)〉.
We use the identities
〈r∗2(t∗)〉 = −3
∂2Fs(k, t
, (6.1)
〈r∗4(t∗)〉 = 15
∂4Fs(k, t
. (6.2)
When inverting the Laplace transforms as in Eq.(4.2),
we again close the contour of integration in the negative
u plane, but in this case it is mathematically essential to
include the circle around the origin in order to account
for the singularities that occur there. The mathematical
situation becomes clear by inspection of the formula for
〈r∗2(t∗)〉M :
〈r∗2(t∗)〉M = 6
∫ +i∞
κ2(u)
∫ +i∞
1− ψ̃G(u) + u/∆
. (6.3)
Because ψG(u) ≈ 1−3 u/2 for small u, this integrand has
a first-order pole at the origin. Higher-order moments
have higher-order poles.
We can use the analogs of Eq.(6.3) for higher moments
and for both the nG and nM distributions to obtain ex-
act asymptotic results in the limits of very small and
very large t∗. For vanishingly small values of t∗, close
the contour around a circle at large values of u where
ψG(u) ≈ 2/u. In this limit, the integrand behaves like
exp (u t∗)/u2, and
〈r∗2〉M ≈ 3∆ t∗. (6.4)
Similar small-t∗ calculations yield:
〈r∗4〉M ≈ 15∆2 t∗2; (6.5)
-6 -5 -4 -3 -2 -1 0 1 2
FIG. 5: (Color online) log
(t∗)〉F ] as a function of
(t∗) for PG = 0, 0.5, 0.75, 0.9, 0.97, 0.999, and 1.0 (red
curve), reading from top to bottom.
〈r∗2〉G ≈ 3∆ t∗2; (6.6)
〈r∗4〉G ≈ 10∆2 t∗3. (6.7)
For very large t∗, close the contour in Eq.(6.3) on
a vanishingly small circle around the origin and use
ψG(u) ≈ 1−3 u/2. The result is that, after times so long
that the tagged molecule no longer remembers whether it
started in a mobile or a glassy region, both 〈r∗2〉M and
〈r∗2〉G converge to the same slowly diffusing Gaussian
distribution for which 〈r∗2〉 ≈ 2 t∗ (with corrections of
the order of 1/∆).
Note several features of these results. As expected, the
initially mobile molecules exhibit rapid (large ∆) Gaus-
sian diffusion at small times. On the other hand, the
displacements of the initially frozen molecules are non-
Gaussian and “pseudo-ballistic” with 〈r∗2〉G ∼ ∆ t∗2.
This behavior has nothing to do with early-stage bal-
listic motion of molecules within their cages but, rather,
is an intrinsic feature of intermediate-stage, anomalous
diffusion in this model.
We turn finally to the full time dependence of the mo-
ments. The results for the weighted average,
〈r∗2(t∗)〉F = PG 〈r∗2(t∗)〉G+(1−PG) 〈r∗2(t∗)〉M , (6.8)
are shown in Fig.5 for a series of different values of PG,
as functions of log10(t
∗). At very short times, these mo-
ments rise as functions of t∗ according to the estimates
in Eqs.(6.4) - (6.7). The pseudo-ballistic behavior is ap-
parent only for PG ≈ 1 because, for smaller PG, these
weighted moments are dominated by the displacements
of the initially mobile molecules. After the initial rise, at
roughly t∗ ≈ ∆−1 = 10−2, the moments cross over from
-3 -2 -1 0 1 2
FIG. 6: Non-Gaussian parameter α2 as a function of log10(t
for PG = 0.
fast β relaxation to slow α decay as the tagged molecules
are repeatedly trapped and then escape from glassy do-
mains. Finally, in Fig.6, we show the non-Gaussian pa-
rameter
3 〈r∗4(t∗)〉M
5 〈r∗2(t∗)〉2M
− 1 (6.9)
only for PG = 0 because, in this approximation where
we have neglected intra-cage vibrational motions, only
the mobile molecules exhibit Gaussian displacement dis-
tributions at early times. Here we see explicitly that the
non-Gaussian behavior occurs during the crossover from
β to α relaxation. Both of these last two figures are ar-
tificial in the sense that we are independently varying
the parameter PG, which in fact ought to be determined
uniquely by the temperature along with the parameters
R∗ and ∆.
VII. CONCLUDING REMARKS
The model proposed here contains a relatively simple
mechanism for producing stretched exponential decay of
molecular correlations, as observed via the self interme-
diate scattering function. That mechanism emerges di-
rectly from the spatial heterogeneity of a glass-forming
liquid. The model illustrates how anomalous diffusion,
as exemplified by a broad exponential (non-Gaussian)
tail of the molecular displacement distribution, is related
to heterogeneity and – indirectly – to the accompanying
stretched exponential behavior. Yet another feature of
the model is the way in which it illustrates how the tran-
sition from relatively fast β relaxation to slow α decay
correlates with the onset of anomalous diffusion. Our
analysis neglects the very short time, strongly localized,
intra-cage fluctuations. The transition between the latter
motions and cage-breaking events seems to us to be out-
side the range of validity of the continuous-time random-
walk theory used here.
There are, of course, many missing ingredients, almost
all of which pertain to the temperature dependence of
the various parameters introduced here. One of us (JSL)
plans to address those issues in a subsequent report.
Acknowledgments
This research was supported by U.S. Department of
Energy Grant No. DE-FG03-99ER45762.
[1] M.D. Ediger, Ann. Rev. Phys. Chem. 51, 99 (2000)
[2] L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, D.
El Masri, D.L’Hote, F. Ladieu, and M. Pierno, Science
310, 1797 (2005).
[3] E. R. Weeks et al, Science 287, 627 (2000)
[4] W. Kob et al, Phys. Rev. Lett. 79, 2827 (1997)
[5] R. Yamamoto and A. Onuki, Phys. Rev. E 58, 3515
(1998)
[6] W. Kob and H. C. Andersen, Phys. Rev. Lett. 73, 1376
(1994); Phys. Rev E 51, 4626 (1995); Phys. Rev E 52,
4134 (1995); W. Kob in “Slow relaxation and nonequi-
librium dynamics in condensed matter,” Les Houches ses-
sion LXXVII (2002)
[7] J. S. Langer, Phys. Rev. E 73, 041504 (2006)
[8] J. S. Langer, Phys. Rev. Lett. 97, 115704 (2006).
[9] A. Widmer-Cooper and P. Harrowell, Phys. Rev. Lett.
96, 185701 (2006); J. Chem. Phys. 126, 154503 (2007).
[10] L.E. Silbert, A.J. Liu, and S.R. Nagel, Phys. Rev. Lett.
95, 098301 (2005).
[11] Y. Jung, J. Garrahan, and D. Chandler, Phys. Rev. E
69, 061205 (2004).
[12] L. Berthier, D. Chandler, and J. Garrahan, Europhysics
Lett. 69 (3), 320 (2005).
[13] P. Hurtado, L. Berthier, and W. Kob, Phys. Rev. Lett.
98, 135503 (2007).
[14] P. Chaudhuri, L. Berthier, and W. Kob, Phys. Rev. Lett.
99, 060604 (2007).
[15] P. Grassberger and I. Procaccia, J. Chem. Phys. 77, 6281
(1982).
[16] C. Monthus and J.-P. Bouchaud, J. Phys. A: Math. Gen.
29 3847 (1996).
[17] J. Bendler, J. Fontanella and M. Shlesinger, Chem. Phys.
284, 311 (2002).
[18] D. Kivelson, G. Tarjus, X. Zhao, and S.A. Kivelson,
Phys. Rev. E 53(1) 751 (1996).
[19] E. W. Montroll and M. F. Shlesinger in Studies in Sta-
tistical Mechanics XI, editors J. L. Lebowitz and E. W.
Montroll, North-Holland Physics Publishing (1984)
[20] J.-P. Bouchaud and A. Georges, Physics Reports 195,
127 (North Holland, 1990).
