Critical Scaling of Shear Viscosity at the Jamming Transition
Critical Scaling of Shear Viscosity at the Jamming Transition
Peter Olsson1 and S. Teitel2
Department of Physics, Ume̊a University, 901 87 Ume̊a, Sweden
Department of Physics and Astronomy, University of Rochester, Rochester, NY 14627
(Dated: July 30, 2021)
We carry out numerical simulations to study transport behavior about the jamming transition of
a model granular material in two dimensions at zero temperature. Shear viscosity η is computed as
a function of particle volume density ρ and applied shear stress σ, for diffusively moving particles
with a soft core interaction. We find an excellent scaling collapse of our data as a function of the
scaling variable σ/|ρc − ρ|
∆, where ρc is the critical density at σ = 0 (“point J”), and ∆ is the
crossover scaling critical exponent. We define a correlation length ξ from velocity correlations in
the driven steady state, and show that it diverges at point J. Our results support the assertion that
jamming is a true second order critical phenomenon.
PACS numbers: 45.70.-n, 64.60.-i, 83.80.Fg
Keywords:
In granular materials, or other spatially disordered
systems such as colloidal glasses, gels, and foams, in
which thermal fluctuations are believed to be negligible,
a jamming transition has been proposed: upon increas-
ing the volume density (or “packing fraction”) of parti-
cles ρ above a critical ρc, the sudden appearance of a
finite shear stiffness signals a transition between flowing
liquid and rigid (but disordered) solid states [1]. It has
further been proposed by Liu and Nagel and co-workers
[2, 3] that this jamming transition is a special second or-
der critical point (“point J”) in a wider phase diagram
whose axes are volume density ρ, temperature T , and ap-
plied shear stress σ (the latter parameter taking one out
of equilibrium to non-equilibrium driven steady states).
A surface in this three dimensional parameter space then
separates jammed from flowing states, and the intersec-
tion of this surface with the equilibrium ρ − T plane at
σ = 0 is related to the structural glass transition.
Several numerical [3, 4, 5, 6, 7, 8, 9, 10], theoretical
[11, 12, 13, 14] and experimental [5, 15, 16, 17, 18] works
have investigated the jamming transition, mostly by con-
sidering behavior as the transition is approached from the
jammed side. In this work we consider the flowing state,
computing the shear viscosity η under applied uniform
shear stress. Previous works have simulated the flow-
ing response to applied shear in glassy systems at finite
temperature [19, 20, 21], and in foams [4] and granular
systems [10] at T = 0, ρ > ρc. Here we consider the
ρ − σ plane at T = 0, showing for the first time that,
near point J, η−1(ρ, σ) collapses to a universal scaling
function of the variable σ/|ρc − ρ|
∆ for both ρ < ρc and
ρ > ρc. We further define a correlation length ξ from
steady state velocity correlations, and show that it di-
verges at point J. Our results support that jamming is a
true second order critical phenomenon.
Following O’Hern et al. [3], we simulate frictionless soft
disks in two dimensions (2D) using a bidisperse mixture
with equal numbers of disks of two different radii. The
radii ratio is 1.4 and the interaction between the particles
V (rij) =
ǫ(1− rij/dij)
2/2 for rij < dij
0 for rij ≥ dij
where rij is the distance between the centers of two par-
ticles i and j, and dij is the sum of their radii. Particles
are non-interacting when they do not touch, and inter-
act with a harmonic repulsion when they overlap. We
measure length in units such that the smaller diameter
is unity, and energy in units such that ǫ = 1. A system
of N disks in an area Lx ×Ly thus has a volume density
ρ = Nπ(0.52 + 0.72)/(2LxLy) . (2)
To model an applied uniform shear stress, σ, we first
use Lees-Edwards boundary conditions [22] to introduce
a uniform shear strain, γ. Defining particle i’s position as
ri = (xi+γyi, yi), we apply periodic boundary conditions
on the coordinates xi and yi in an Lx × Ly system. In
this way, each particle upon mapping back to itself under
the periodic boundary condition in the ŷ direction, has
displaced a distance ∆x = γLy in the x̂ direction, re-
sulting in a shear strain ∆x/Ly = γ. When particles do
not touch, and hence all mutual forces vanish, xi and yi
are constant and a time dependent strain γ(t) produces a
uniform shear flow, dri/dt = yi(dγ/dt)x̂. When particles
touch, we assume a diffusive response to the inter-particle
forces, as would be appropriate if the particles were im-
mersed in a highly viscous liquid or resting upon a rough
surface with high friction. This results in the following
equation of motion, which was first proposed as a model
for sheared foams [4],
dV (rij)
x̂ . (3)
The strain γ is then treated as a dynamical variable,
http://arxiv.org/abs/0704.1806v2
obeying the equation of motion,
LxLyσ −
dV (rij)
 , (4)
where the applied stress σ acts like an external force on
γ and the interaction terms V (rij) depend on γ via the
particle separations, rij = ([xi−xj]Lx +γ[yi−yj]Ly , [yi−
yj ]Ly), where by [. . .]Lµ we mean that the difference is
to be taken, invoking periodic boundary conditions, so
that the result lies in the interval (−Lµ/2, Lµ/2]. The
constants D and Dγ are set by the dissipation of the
medium in which the particles are embedded; we take
units of time such that D = Dγ ≡ 1.
In a flowing state at finite σ > 0, the sum of the inter-
action terms is of order O(N) so that the right hand side
of Eq. (4) is O(1). The strain γ(t) increases linearly in
time on average, leading to a sheared flow of the particles
with average velocity gradient dvx/dy = 〈dγ/dt〉, where
vx(y) is the average velocity in the x̂ direction of the par-
ticles at height y. We then measure the shear viscosity,
defined by,
dvx/dy
〈dγ/dt〉
. (5)
We expect η−1 to vanish in a jammed state.
We integrate the equations of motion, Eqs. (3)-(4),
starting from an initial random configuration, using the
Heuns method. The time step ∆t is varied according to
system size to ensure our results are independent of ∆t.
We consider a fixed number of particles N , in a square
system L ≡ Lx = Ly, and vary the volume density ρ by
adjusting the length L according to Eq. (2). We simulate
for times ttot such that the total relative displacement per
unit length transverse to the direction of motion is typ-
ically γ(ttot) ∼ 10, with γ(ttot) ranging between 1 and
200 depending on the particular system parameters.
In Fig. 1 we show our results for η−1 using a fixed
small shear stress, σ = 10−5, representative of the σ → 0
limit. Our raw results are shown in Fig. 1a for several
different numbers of particles N from 64 to 1024. Com-
paring the curves for different N as ρ increases, we see
that they overlap for some range of ρ, before each drops
discontinuously into a jammed state. As N increases, the
onset value of ρ for jamming increases to a limiting value
ρc ≃ 0.84 (consistent with the value for random close
packing [3]) and η−1 vanishes continuously. For finite
N , systems jam below ρc because there is always a fi-
nite probability to find a configuration with a force chain
spanning the width of the system, thus causing it to jam;
and at T = 0, once a system jams, it remains jammed
for all further time. As the system evolves dynamically
with increasing simulation time, it explores an increasing
region of configuration space, and ultimately finds a con-
figuration that causes it to jam. The statistical weight
10-3 10-2 10-1
eta^-1
eta^-1
eta^-1
eta^-1
eta^-1
beta=1.6
volume density  #c " #
(b)$ = 10
#c = 0.8415
% = 1.65
!"1~ |# " #c|
0.78 0.79 0.80 0.81 0.82 0.83 0.84
eta^-1
eta^-1
eta^-1
eta^-1
eta^-1
volume density  #
(a)$ = 10"5
FIG. 1: (color online) a) Plot of inverse shear viscosity η−1
vs volume density ρ for several different numbers of particles
N , at constant small applied shear stress σ = 10−5. As N
increases, one see jamming at a limiting value of the density
ρc ∼ 0.84. b) Log-log replot of the data of (a) as η
−1 vs
ρc − ρ, with ρc = 0.8415. The dashed line has slope β = 1.65
indicating the continuous algebraic vanishing of η−1 at ρc with
a critical exponent β.
of such jamming configurations decreases, and hence the
average time required to jam increases, as one either de-
creases ρ, or increases N [3]. In the limit N → ∞, we
expect jamming will occur in finite time only for ρ ≥ ρc.
In Fig. 1b we show a log-log plot of η−1 vs ρc − ρ, us-
ing a value ρc = 0.8415. We see that the data in the
unjammed state is well approximated by a straight line
of slope β = 1.65, giving η−1 ∼ |ρ − ρc|
β in agreement
with the expectation that point J is a second order phase
transition.
If point J is indeed a true critical point, one expects
that its influence will be felt also at finite values of the
stress σ, with η−1 obeying a typical scaling law,
η−1(ρ, σ) = |ρ− ρc|
|ρ− ρc|∆
. (6)
Here z ≡ σ/|ρ−ρc|
∆ is the crossover scaling variable, ∆ is
the crossover scaling critical exponent, and f−(z), f+(z)
are the two branches of the crossover scaling function for
ρ < ρc and ρ > ρc respectively.
In Fig. 2 we show a log-log plot of inverse shear vis-
cosity η−1 vs applied shear stress σ, for several different
values of volume density ρ. Our results are for systems
large enough that we believe finite size effects are negli-
gible. We use N = 1024 for ρ < 0.844 and N = 2048 for
ρ ≥ 0.844. Again we see that ρc ≃ 0.8415 separates two
limits of behavior. For ρ < ρc, log η
−1 is convex in log σ,
decreasing to a finite value as σ → 0. For ρ > ρc, log η
is concave in log σ, decreasing towards zero as σ → 0.
The dashed straight line, separating the two regions of
behavior, indicates the power law dependence that is ex-
pected exactly at ρ = ρc (see below). Similar power
law behavior at ρc was recently found in simulations of a
three dimensional granular material [23].
10-5 10-4 10-3 10-2
shear stress   # 
$=0.830
$=0.834
$=0.836
$=0.838
$=0.840
$=0.841
$=0.842
$=0.844
$=0.848
$=0.852
$=0.856
$=0.860
$=0.864
$=0.868
#=0.0012
FIG. 2: (color online) Plot of inverse shear viscosity η−1 vs
applied shear stress σ for several different values of the vol-
ume density ρ. The dashed line represents the power law
dependence expected exactly at ρ = ρc and has a slope
β/∆ = 1.375. Solid lines are guides to the eye. Points la-
beled σ = 0.0012 correspond to densities ρ = 0.870, 0.872,
0.874, 0.876, and 0.878.
In Fig. 3 we replot the data of Fig. 2 in the scaled
variables η−1/|ρ−ρc|
β vs σ/|ρ−ρc|
∆. Using ρc = 0.8415,
β = 1.65 (the same values used in Fig. 1b) and ∆ = 1.2,
we find an excellent scaling collapse in agreement with
the prediction of Eq. (6). As the scaling variable z → 0,
f−(z) → constant; this gives the vanishing of η
−1 ∼ |ρ−
β at σ = 0. As z → ∞, both branches of the scaling
function approach a common curve, f±(z) ∼ z
β/∆, so
that precisely at ρ = ρc, η
−1 ∼ σβ/∆ as σ → 0 [24].
This is shown as the dashed line in both Figs. 3 and 2. A
similar scaling collapse of η has been found in simulations
[20] of a sheared Lennard-Jones glass, as a function of
temperature and applied shear strain rate γ̇, but only
above the glass transition, T > Tc. By comparing the
goodness of the scaling collapse as parameters are varied,
we estimate the accuracy of the critical exponents to be
roughly, β = 1.7± 0.2 and ∆ = 1.2± 0.2.
That the crossover scaling exponent ∆ > 0, implies
that σ is a relevant variable in the renormalization group
sense, and that critical behavior at finite σ should be
in a different universality class from the jamming tran-
sition at point J (i.e. σ = 0). The nature of jamming
at finite σ > 0 will be determined by the behavior of
the branch of the crossover scaling function f+(z), that
describes behavior for ρ > ρc. From Fig. 3 we see that
f+(z) is a decreasing function of z. If f+(z) vanishes only
when z → 0, then Eq. (6) implies that η−1 vanishes for
ρ > ρc only when σ = 0, and so there will be no jamming
at finite σ > 0. If, however, f+(z) vanishes at some fi-
nite z0, then η
−1 will vanish whenever σ/(ρ− ρc)
∆ = z0;
there will then be a line of jamming transitions emanat-
ing from point J in the ρ − σ plane given by the curve
ρ∗(σ) = ρc + (σ/z0)
1/∆. If f+(z) vanishes continuously
10-2 10-1 100
z = %/|# " #c|
#c = 0.8415
$ = 1.65
& = 1.2
f"(z)
# < #c
f+(z)
# > #c~ z$/&
FIG. 3: (color online) Plot of scaled inverse viscosity η−1/|ρ−
β vs scaled shear stress z ≡ σ/|ρ − ρc|
∆ for the data of
Fig. 2. We find an excellent collapse to the scaling form of
Eq. (6) using values ρc = 0.8415, β = 1.65 and ∆ = 1.2.
The dashed line represents the large z asymptotic dependence,
∼ zβ/∆. Data point symbols correspond to those used in
Fig. 2.
at z0, jamming at finite σ will be like a second order
transition; if f+(z) jumps discontinuously to zero at z0,
it will be like a first order transition. Such a first order
like transition has been reported in simulations [20, 21]
of sheared glasses at finite temperature below the glass
transition, T < Tc. However, recent simulations [10] of
a granular system at T = 0, ρ > ρc, showed that a sim-
ilar first order like behavior was a finite size effect that
vanished in the thermodynamic limit. With these obser-
vations, we leave the question of criticality at finite σ to
future work
The critical scaling found in Fig. 3 strongly suggests
that point J is indeed a true second order phase transi-
tion, and thus implies that there ought to be a diverg-
ing correlation length ξ at this point. Measurements
of dynamic (time dependent) susceptibilities have been
used to argue for a divergent length scale in both the
thermally driven glass transition [25], and the density
driven jamming transition [17]. Here we consider the
equal time transverse velocity correlation function in the
shear driven steady state,
g(x) = 〈vy(xi, yi)vy(xi + x, yi)〉 , (7)
where vy(xi, yi) is the instantaneous velocity in the ŷ
direction, transverse to the direction of the average shear
flow, for a particle at position (xi, yi). The average is over
particle positions and time. In the inset to Fig. 4 we plot
g(x)/g(0) vs x for three different values of ρ at fixed σ =
10−4 and number of particlesN = 1024. We see that g(x)
decreases to negative values at a well defined minimum,
before decaying to zero as x increases. We define ξ to be
the position of this minimum. That g(ξ) < 0, indicates
10-2 10-1 100 101
z = %/|# " #c|
#c = 0.8415
$ = 0.6
& = 1.2h"(z)
# < #c
h+(z)
# > #c
~ z$/&
0 5 10 15 20
rho=0.83
rho=0.834
rho=0.838
# = 0.830
# = 0.834
# = 0.838
N = 1024   % = 10-4
FIG. 4: (color online) Inset: Normalized transverse velocity
correlation function g(x)/g(0) vs longitudinal position x for
N = 1024 particles, applied shear stress σ = 10−4, and vol-
ume densities ρ = 0.830, 0.834 and 0.838. The position of
the minimum determines the correlation length ξ. Main fig-
ure: Plot of scaled inverse correlation length ξ−1/|ρ− ρc|
scaled shear stress z ≡ σ/|ρ− ρc|
∆ for the data of Fig. 2. We
find a good scaling collapse using values ρc = 0.8415, ∆ = 1.2
(the same as in Fig. 3) and ν = 0.6. Data point symbols
correspond to those used in Fig. 2.
that regions separated by a distance ξ are anti-correlated.
We can thus interpret the sheared flow in the unjammed
state as due to the rotation of correlated regions of length
ξ. Similar behavior, leading to a similar definition of
ξ, has previously been found [26] in correlations of the
nonaffine displacements of particles in a Lennard-Jones
glass, in response to small elastic distortions.
As with viscosity, we expect the correlation length
ξ(ρ, σ) to obey a scaling equation similar to Eq. (6). We
consider here the inverse correlation length ξ−1, which
like η−1 should vanish at the jamming transition, obey-
ing the scaling equation,
ξ−1(ρ, σ) = |ρ− ρc|
|ρ− ρc|∆
. (8)
The correlation length critical exponent is ν, but the
crossover exponent ∆ remains the same as for the vis-
cosity.
In Fig. 4 we plot the scaled inverse correlation length,
ξ−1/|ρ−ρc|
ν vs the scaled stress, σ/|ρ−ρc|
∆. Using ρc =
0.8415 and ∆ = 1.2, as was found for the scaling of η−1,
we now find a good scaling collapse for ξ−1 by taking the
value ν = 0.6. By comparing the goodness of the collapse
as ν is varied, we estimate ν = 0.6±0.1. From the scaling
equation Eq. (8) we expect both branches of the scaling
function to approach the power law h±(z) ∼ z
ν/∆ as
z → ∞, so that ξ−1 ∼ σν/∆ as σ → 0 at ρ = ρc [24].
This is shown as the dashed line in Fig. 4. Our result
is consistent with the conclusion “ν is between 0.6 and
0.7” of Drocco et al. [7] for the flowing phase, ρ < ρc.
It also agrees with ν = 0.71 ± 0.08 found by O’Hern et
al. [3] from a finite size scaling argument. Wyart et al.
[14] have proposed a diverging length scale with exponent
ν = 0.5 by considering the vibrational spectrum of soft
modes approaching point J from the jammed side, ρ > ρc.
While our results cannot rule out ν = 0.5, our scaling
collapse in Fig. 4 does seem somewhat better when using
the larger value 0.6.
This work was supported by Department of Energy
grant DE-FG02-06ER46298 and by the resources of the
Swedish High Performance Computing Center North
(HPC2N). We thank J. P. Sethna, L. Berthier, M. Wyart,
J. M. Schwarz, N. Xu, D. J. Durian, A. J. Liu and
S. R. Nagel for helpful discussion.
[1] Jamming and Rheology, edited by A. J. Liu and
S. R. Nagel (Taylor & Francis, New York, 2001).
[2] A. J. Liu and S. R. Nagel, Nature 396, 21 (1998).
[3] C. S. O’Hern et al. Phys. Rev. E 68, 011306 (2003).
[4] D. J. Durian, Phys. Rev. Lett. 75, 4780 (1995) and Phys.
Rev. E 55, 1739 (1997).
[5] H. A. Makse, D. L. Johnson and L. M. Schwartz, Phys.
Rev. Lett. 84, 4160 (2000).
[6] C. S. O’Hern et al., Phys. Rev. Lett. 86, 000111 (2001)
and 88, 075507 (2002).
[7] J. A. Drocco et al., Phys. Rev. Lett. 95, 088001 (2005).
[8] L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. Lett.
95, 098301 (2005) and Phys. Rev. E 73, 041304 (2006).
[9] W. G. Ellenbroek et al., Phys. Rev. Lett. 97, 258001
(2006).
[10] N. Xu and C. S. O’Hern, Phys. Rev. E 73, 061303 (2006).
[11] J. M. Schwarz, A. J. Liu and L. Q. Chayes, Europhys.
Lett. 73, 560 (2006).
[12] C. Toninelli, G. Biroli and D. S. Fisher, Phys. Rev. Lett.
96, 035702 (2006).
[13] S. Henkes and B. Chakraborty, Phys. Rev. Lett. 95,
198002 (2005).
[14] M. Wyart, S. R. Nagel, T. A. Witten, Europhys. Lett.
72, 486-492 (2005); M. Wyart et al., Phys. Rev. E 72,
051306 (2005); C. Brito and M. Wyart, Europhys. Lett.
76, 149 (2006).
[15] V. Trappe et al., Nature (London) 411, 772 (2001).
[16] T. S. Majmudar et al., Phys. Rev. Lett. 98, 058001
(2007).
[17] A. .S. Keys et al., Nature physics 3, 260 (2007).
[18] M. Schröter et al., Europhys Lett. 78, 44004 (2007).
[19] R. Yamamoto and A. Onuki, Phys. Rev. E 58, 3515
(1998).
[20] L. Berthier and J.-L. Barat, J. Chem. Phys. 116, 6228
(2002).
[21] F. Varnik, L. Bocquet and J.-L.Barrat, J. Chem. Phys.
120, 2788, (2004).
[22] D. J. Evans and G. P. Morriss, Statistical Mechanics of
Non-equilibrium Liquids (Academic, London, 1990).
[23] T. Hatano, M. Otsuki and S. Sasa, condmat/0607511.
[24] In general, one should consider nonlinear scaling vari-
ables. In our case, the most important correction would
be to replace ρ− ρc in Eq. (6) by gρ(ρ, σ) ≡ ρ− ρc + cσ
this could lead to noticeable corrections to our scaling
equation near ρ = ρc. However, since we find ∆ > 0.5,
our conclusion that η−1 ∼ σβ/∆ at ρ = ρc remains valid.
See, A. Aharony and M. E. Fisher, Phys. Rev. B 27, 4394
(1983).
[25] L. Berthier et al., Science 310, 1797 (2005).
[26] A. Tanguy et al., Phys. Rev. B 66, 174205 (2002).
