Multiscaling at Point J: Jamming is a Critical Phenomenon
J. A. Drocco1, M. B. Hastings2,3, C. J. Olson Reichhardt3,
and C. Reichhardt2,3
1Department of Physics, University of Notre Dame, Notre Dame, Indiana 46656 2Center for Nonlinear Studies and 3Theoretical Division,
Los Alamos National
Laboratory, Los Alamos, New Mexico 87545
(Received 13 October 2003; revised manuscript received
5 November 2004; published 19 August 2005)
We analyze the jamming transition that occurs as a function
of increasing packing density
in a disordered two-dimensional assembly of
disks
at zero temperature for "Point J" of the recently proposed jamming
phase diagram.
We measure the total number of moving disks and
the transverse length of the moving region,
and find a power law divergence as the packing density increases toward
a critical jamming density.
This provides evidence that
the T = 0 jamming transition
as a function of packing density is a second order
phase transition.
Additionally we find evidence for multiscaling, indicating the
importance of long tails in the velocity fluctuations.
Distribution of Velocities Conclusion References
There has been a surge of activity in jamming phenomena for
T=0 systems such as granular materials, foams, and
colloids, where jamming is
defined as the onset of a nonvanishing yield stress in
a disordered state[1,2].
Liu and Nagel proposed a jamming phase diagram containing three
axes: temperature T, the inverse packing fraction 1/ϕ,
and shear σ [1]. The
system is jammed within a three-dimensional dome;
above the jamming transition, the system behaves
as a rigid solid.
Recently O'Hern et al. [3] studied the
area of the jamming phase diagram
at T = σ = 0
as a function of packing density ϕ,
and
showed that a well defined sharp jamming transition appears
at "Point J."
Elsewhere on the jamming phase diagram, the boundary between
jammed and unjammed states is not sharp since its definition is
sensitive to the experimental time scale.
It was noted in Ref. [3] that, although Point J does not exist
for liquids, the behavior near Point J
could be relevant to the physics near the glass transition.
A key question is whether Point J is a true
continuous phase
transition with a diverging spatial correlation length.
Since the physics at Point J
is non-thermal, any continuous phase transition behavior could be
dominated by rare fluctuations.
If this is the case, simple scaling close to the critical density
may not occur even if there is an
underlying phase transition.
Here, we consider the motion of a single
disk driven through an arrangement of disks
as the jamming transition is approached,
a system which,
as suggested in Ref. [3],
should provide a direct test of
the nature of the jamming transition.
The reader is invited to try the following experiment:
place a large collection of coins flat on a desk, so that they are almost
touching. Then,
push one coin and observe what other coins move.
At low packing fractions
ϕ, the driven disk or coin occasionally contacts other disks, which can in
turn contact still other disks, but the total number of disks moving
is small. At a certain density ϕ = ϕc,
however, the entire system must move simultaneously as a unit, and jams.
Force chains can be observed emanating from the driven disk [4].
In these jamming problems, the tail of the
particle velocity distribution P(v) at small velocities
controls the behavior of the system, since once the system enters
a jammed state, with everything stopped in the thermodynamic limit,
it can never exit the jammed state. The best method for characterizing
the tail of P(v) is via multiscaling, which we employ
in this Letter.
For fixed driving force,
a small v implies that a large number of other particles are dragged
with the driven particle, but we can also characterize the number of particles
moved by the driven particle by examining the force transmission. We will see
that this measurement provides direct evidence for a diverging length
scale.
Our data also suggests that the jamming transition
for very large systems coincides with the random close packed density.
Figure 1: Sample geometry for L=60 and (a) ϕ = 0.656,
(b) ϕ = 0.811, (c) ϕ = 0.837.
Large black dot: the probe disk (size exaggerated for clarity);
gray dots: moving disks; small dots: non-moving disks.
We simulate a binary mixture of two-dimensional disks
with stiff spring repulsive interactions of radius rA
and rB at T = 0
(see Fig. 1).
For all the densities
we consider here, we find << 1% overlap in the radii
as the disks interact, indicating that the spring constant is sufficiently
large to provide a good approximation to hard core disks.
The bimodal disk distribution, with a size ratio of 1.4:1,
is chosen to create a disordered arrangement and
avoid formation of a regular lattice.
We also performed
simulations with a size ratio of 1.7:1 with substantially identical results.
We employ overdamped dynamics such that the velocity of each
disk
is proportional to the force acting on it.
The equation of motion for the disks is
ηvi =
∑ i ≠ j
k(|rij| − reff)
rij
|rij|
+ fd .
(1)
Here k = 200 is the strength of the stiff restoring spring [5],
rij = ri − rj, and for all but one
of the disks, the external driving force fd=0.
The single driven disk (large dot in Fig. 1) has fd=0.1.
Disks only interact if they are separated by
a distance smaller than the sum of their radii,
|rij| < reff = r(i)+r(j), where
r(i) equals either rA or rB, depending on the given
particle.
We tested several different values
of fd and chose the value small enough to minimize
overlap as mentioned above.
We use periodic boundary conditions with a range of linear system sizes L=24
to 60. For
L=60, the
system contains N=2600 background disks
at a density of ϕ = 0.8395.
We prepare the system by randomly dropping
nonoverlapping disks up to ϕ ≈ 0.6.
To reach higher densities, we add disks at
randomly selected interstitial locations,
reduce the radii of all disks,
and then increase the radii back to the initial values while allowing the
system to evolve under a small temperature. This produces
nonoverlapping configurations at densities up to ϕc.
To reach the maximum possible density,
√3π/6 ≈ 0.907, requires phase separating the system.
At
the significantly lower ϕc ≈ 0.839 we find no phase separation.
Since the system is not jammed below ϕc, we believe that
we successfully uniformly sampled all available states,
and
that phase separation does
not occur until some ϕsep > ϕc.
To study the velocity curve,
we drive a single large disk along a 45° angle
over a distance of √2L/5.
The time required for this motion is much longer
than the time scale of any brief initial transients in the velocity.
We measure
the effective velocity of the driven disk
over the length of the simulation.
For computational reasons, it is not possible to simulate
a disk
moving at an arbitrarily slow velocity through the system. In these
cases, we declared that the system had jammed; this definition
agreed well with the critical ϕc ≈ 0.839 given by
curve fitting below.
In Fig. 1 we show snapshots of simulation results
for different densities, with
disks counted as moving and shaded in gray
if they were connected via a force
contact to the driven disk.
At the lowest
ϕ = 0.656,
in Fig. 1(a),
the driven disk moves easily and interacts
with
only a
few neighbors.
In Fig. 1(b) for
ϕ = 0.811,
close to jamming, a larger portion of the disks are moving,
while in Fig. 1(c) at ϕ = 0.837, the system is
sufficiently close to ϕc that the entire (finite-size) sample is moving.
Figure 2: Time series of velocity for L=60 and ϕ = (a) 0.656, (b) 0.747,
and (c) 0.837.
Figure 2 shows example time series
for the velocity v parallel to the applied drive
at ϕ = 0.656, 0.747, and 0.837.
Not only does the average
velocity decrease as ϕ increases,
but the time series becomes more intermittent. At ϕ = 0.837,
the velocity is usually very small, but there are
occasional bursts of much higher
velocity.
We emphasize that
for ϕ < ϕc
the threshold force
required to move the driven particle vanishes.
Instead,
v is proportional to the drive fd at any moment of time. This
follows from dimensional grounds, since for a hard-core system, there is
no physical parameter with units of force.
This is very different
from the behavior shown in
a recent numerical study of single probe particles driven through a
system with softer, long-range forces [7],
where there is a non-vanishing, but finite,
depinning force
at all densities due to the long-range nature of the force.
Distribution of Velocities-
Newton's third law, combined with Eq. (1), leads to the
result that η∑ivi=fd. Thus, if the driven
particle is not in contact with any other particles, it moves at
v=fd/η. If there are a total of n particles moving together,
including the driven particle, then they each move at
v=fd/n.
For ϕ > ϕc, all N particles in the system
move
together, and
v is vanishingly small in the thermodynamic limit.
For ϕ < ϕc, the number of
particles moving is finite in the thermodynamic
limit and diverges as ϕ approaches ϕc.
One measurement which may
show scaling as ϕ approaches ϕc is the average velocity
v of the
driven disk. However,
v is a random, time-dependent
variable, and
for any density ϕ < ϕc there is a nonvanishing probability of observing
any given, arbitrarily small v.
Since the other disks are distributed randomly, there is still some
probability that in any local region the density will
exceed ϕc [8] (this argument is inspired by Lifshitz
tails). If this region
encloses n disks,
the
probability of finding such a
region is exponentially
small in n, but it is nonvanishing. When the driven disk
hits such a region, it behaves as if it jams, and
v slows to
a value of order 1/n until the disk either escapes the region or
pushes other disks
into neighboring
regions of lower density. Thus, with exponentially small probability in
n, v ∼ O(1/n). Once the particle slows, however, it takes a long
time for it to leave the region. Thus, one expects to see
an intermittent v.
To measure the intermittency, we use multifractal[9] scaling.
For a given ϕ, let
p(v) dv be the probability of measuring a given
v at a single instant in time. Define the q-th inverse moment by
M(q)=∫dvp(v) v−q as the time average of the q-th power of the
inverse
velocity. We then define a set of multifractal exponents
τ(q) by
M(q) =
⌠ ⌡
dvp(v) v−q ∝ (ϕc−ϕ)−τ(q).
(2)
If v
were constant throughout the simulation for a given ϕ, we would
have M(q)1/q=M(1).
If instead v had relative
fluctuations of order unity about some characteristic velocity, v0(ϕ),
then M(q)1/q would
not equal M(1) in general, but
τ(q)/q would still be independent of q.
We will instead find that
τ(q)/qdepends on q. This implies that
v does
not simply fluctuate about a single
v0, but is instead much more intermittent,
with the particle
sometimes getting stuck for a long time at a slow velocity, then traveling
much more rapidly.
Figure 3: Plot of M(q,ϕ)1/q versus ϕ for q=−1, 1, 2, and 4
(from bottom to top). The curves with larger q
are noisier. Inset: M(q,ϕ)1/q versus ϕc−ϕ on
a log-log scale for q=−3,−1,1,3. A curve fit shows the power law scaling.
Figure 4: (a): Plot of τ(q)/q against q. The q dependence of
τ(q)/q shows the existence of multiscaling in this system.
(b): Plot of M(1) (crosses) and M(3)1/3 (circles) against M(−1)−1.
We compute the moments M and extract the exponents τ(q)
by fitting the moments to the form c*(ϕc−ϕ)−τ(q).
The
value of ϕc determined from these fits is
independent
of q to high accuracy for q ≥ −2, and thus we can fix
ϕc ≈ 0.839 as the onset of jamming.
In Fig. 3 we plot M(q,ϕ)1/q versus ϕ for L=60 and
various q, averaged over three realizations of the system.
In the inset we show the log-log plot of the curves
in the main panel
with solid lines indicating power law fits. The
curves not only
fail to overlap, but also
have different slopes, indicating the presence of multiscaling.
In Fig. 4(a), we plot τ(q)/q versus q. For large positive
and negative q, the plot asymptotes. The asymptote at large positive q
reflects the scaling of the typical (disregarding exponentially
rare possibilities discussed above) slowest velocity
of the system which goes to zero as (ϕ−ϕc)limq→∞τ(q)/q. Similarly, the asymptote at negative q reflects the typical
largest velocity. Since limq→ −∞τ(q)/q is very close
to zero, the typical largest velocity is largely independent of ϕ, as
seen in Fig. 2. In fact, it is likely that at sufficiently negative q the
exponent τ(q)/q becomes zero.
The error in τ(q) due to statistical fluctuations is negligible.
If we compare the τ(q) obtained by curve fitting M(q,ϕ) averaged
over all realizations to τ(q) obtained by using only a single
realization, the
difference in exponents is of order 0.001 to 0.002. Instead,
the major
source of error is finite size effects. The straight lines in the
inset to Fig. 3 are based on fitting over a certain range of ϕ. For
ϕ closer to ϕc than the endpoint of the
fitting lines, the statistical
noise in the data increases significantly and some of the realizations
jam while others do not.
Thus, for finite N the
jamming threshold is not well defined, and these ϕ are so close to
ϕc that finite size effects may be important. The error bars in
Fig. 4 are based on this consideration of finite size effects. The low end
of the error bars corresponds to the τ(q) obtained by fitting only
over the range shown in the inset to Fig. 3, while the high end corresponds
to including all ϕ < ϕc. The low end tends to underestimate the
difference in τ(q)/q for different
q. τ(q)/q is clearly q-dependent since the difference in
τ(q)/q between q=−1 and q=2, for example,
is definitely outside
the error bars.
In Fig. 4(b) we plot M(q)1/q as
a function of M(−1)−1 for q=1,3 over the
same range as in
Fig. 3 to illustrate that the data obey extended self-similarity [10].
The slopes on the log-log plot are
1.25±0.15 and 1.45 ±0.2
respectively; as expected, this is within error bars of
the ratio of the exponents −τ(1)/τ(−1) and −(τ(3)/3)/τ(−1).
It is difficult to establish that these slopes are different from
each other, but the slopes are definitely different from unity, which
already indicates multiscaling.
The fact that τ(−1) is
less
than one makes the size of the scaling regime seem small, since a wide
scaling range in ϕ leads to only a small range in τ(−1);
however, the curve fit is very good over the available range.
Figure 5: Plot of nmoving against ϕ. Solid line: L=60.
Dashed line: L=48. Dotted line: a power law fit to the L=60 curve.
Left inset: The L=60 curve and the fit against ϕc−ϕ,
showing scaling. Right inset: finite size scaling.
Figure 1 shows that as the jamming transition is approached,
the number of moving disks increases and there is a diverging
length scale as the jamming is approached. In order to quantify this,
we show in Fig. 5
the number of moving disks nmoving
vs ϕ for systems
of linear size
L=48 and L=60.
Here a clear divergence appears as the
critical density is approached. The divergence is cut off when
nmoving equals the total number of disks.
We have also considered smaller systems
for different
parameters
of drive and disk radii and again
observe a divergence; however, these smaller systems give
a much lower resolution and hence a larger error on the estimated exponent.
In Fig. 5 we fit a power law to the largest
system with (ϕc − ϕ)τ, where τ is between 1.2 and 1.46.
We also measured the sum over moving disks of the squared distance
between the driven disk and the moving disk,
a quantity nmovingl2, and find that this number
diverges as (ϕc−ϕ)−σ, with σ between 2.35 and 2.6.
The number of moving disks cannot directly be compared to the
velocity v of the driven disk, as
the driven disk can push other disks
normal to the drive so that nmoving may be much larger than 1/v.
We obtained the exponents τ and σ by a time average of the
number of moving disks, and did not perform a
multiscaling calculation, though we can extract a length
from a simple scaling analysis. If the moving disks form a cluster
with dimension d and length scale ξ diverging as
ξ ∝ (ϕc−ϕ)−ν, then τ = νd and σ = ν(d+2).
Given the values of τ, σ, and our own qualitative observations,
we find d=2, so that ν is between 0.6 and 0.7. This provides
a direct measurement of the diverging length scale, in reasonable agreement
with the finite-size scaling result [3] ν = 0.71 ±0.12.
We note that Fig. 5 in our case is completely consistent with finite-size
scaling as shown in the right inset.
Scale the y-axis by L2 (the scaling for
the total number of particles in the system), and scale ϕc−ϕ by
L−2/τ. Then, the curves for different L automatically
collapse within the scaling region, and since the divergence is
cutoff at nmoving ∝ L2 for all L, these curves also
collapse close to ϕc.
Finally, our
ϕc ≈ 0.839 is very close to the critical packing density found in [3].
Conclusion-
We have studied a system
of a T=0, zero shear 2D disordered assembly of disks
at densities below and up to
Point J in the recently proposed jamming phase diagram.
A single probe disk is pushed with a constant drive through
the other disks. Upon increasing the packing density,
we find a jamming transition associated with a power
law divergence in the number of moving disks
and a diverging spatial correlation
length, indicating that
Point J is a true
continuous phase transition.
In addition, we
show
that the tails of the disk velocity distribution
play an important role in this
transition, since due to the non-thermal nature of the
system, once the particle stops moving,
it cannot restart.
Using the multifractal moments τ(q), we chart the tails
and show that the τ(q) do not depend linearly on q;
thus, the system exhibits multiscaling.
We also find evidence for a diverging correlation length
of the force contacts
as the jamming transition is approached from below.
Our results
show that there is an underlying
second order non-thermal phase transition
at the
jamming transition.
We suggest that experimental
work should consider both the velocity time series and
the number of moving particles in the surrounding media,
and should test for
multiscaling behavior
as the jamming transition is approached.
Particular experiments would include driving a single disk on a
flat surface through a disordered assembly of other disks for increasing
density, or driving individual colloids through glassy assemblies of
other colloids.
Acknowledgments-
We thank E. Weeks for useful discussions.
This work was supported by the U.S. DOE under Contract No.
W-7405-ENG-36.
