Critical Scaling of Shear Viscosity at the Jamming Transition

Critical Scaling of Shear Viscosity at the Jamming Transition

Peter Olsson Department of Physics, Umeå University, 901 87 Umeå, Sweden    S. Teitel Department of Physics and Astronomy, University of Rochester, Rochester, NY 14627
July 6, 2019

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 , where is the critical density at (“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.

45.70.-n, 64.60.-i, 83.80.Fg

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 increasing the volume density (or “packing fraction”) of particles above a critical , the sudden appearance of a finite shear stiffness signals a transition between flowing liquid and rigid (but disordered) solid states LiuNagel1 (). It has further been proposed by Liu and Nagel and co-workers LiuNagel2 (); OHern2 () that this jamming transition is a special second order critical point (“point J”) in a wider phase diagram whose axes are volume density , temperature , and applied 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 intersection of this surface with the equilibrium plane at is related to the structural glass transition.

Several numerical OHern2 (); Durian (); Makse (); OHern1 (); Drocco (); Silbert (); Ellenbroek (); Ning (), theoretical Schwarz (); Fisher (); Chak (); Wyart () and experimental Makse (); Weitz (); Majmudar (); Durian2 (); Swinney () works have investigated the jamming transition, mostly by considering 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 flowing response to applied shear in glassy systems at finite temperature Yamamoto (); Berthier (); Varnik (), and in foams Durian () and granular systems Ning () at , . Here we consider the plane at , showing for the first time that, near point J, collapses to a universal scaling function of the variable for both and . We further define a correlation length from steady state velocity correlations, and show that it diverges at point J. Our results support that jamming is a true second order critical phenomenon.

Following O’Hern et al. OHern2 (), 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 is,


where is the distance between the centers of two particles and , and is the sum of their radii. Particles are non-interacting when they do not touch, and interact 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 . A system of disks in an area thus has a volume density


To model an applied uniform shear stress, , we first use Lees-Edwards boundary conditions LeesEdwards () to introduce a uniform shear strain, . Defining particle ’s position as , we apply periodic boundary conditions on the coordinates and in an system. In this way, each particle upon mapping back to itself under the periodic boundary condition in the direction, has displaced a distance in the direction, resulting in a shear strain . When particles do not touch, and hence all mutual forces vanish, and are constant and a time dependent strain produces a uniform shear flow, . When particles touch, we assume a diffusive response to the inter-particle forces, as would be appropriate if the particles were immersed 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 Durian (),


The strain is then treated as a dynamical variable, obeying the equation of motion,


where the applied stress acts like an external force on and the interaction terms depend on via the particle separations, , where by we mean that the difference is to be taken, invoking periodic boundary conditions, so that the result lies in the interval . The constants and are set by the dissipation of the medium in which the particles are embedded; we take units of time such that .

In a flowing state at finite , the sum of the interaction terms is of order so that the right hand side of Eq. (4) is . The strain increases linearly in time on average, leading to a sheared flow of the particles with average velocity gradient , where is the average velocity in the direction of the particles at height . We then measure the shear viscosity, defined by,


We expect 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 is varied according to system size to ensure our results are independent of . We consider a fixed number of particles , in a square system , and vary the volume density by adjusting the length according to Eq. (2). We simulate for times such that the total relative displacement per unit length transverse to the direction of motion is typically , with ranging between and depending on the particular system parameters.

In Fig. 1 we show our results for using a fixed small shear stress, , representative of the limit. Our raw results are shown in Fig. 1a for several different numbers of particles from to . Comparing the curves for different as increases, we see that they overlap for some range of , before each drops discontinuously into a jammed state. As increases, the onset value of for jamming increases to a limiting value (consistent with the value for random close packing OHern2 ()) and vanishes continuously. For finite , systems jam below because there is always a finite probability to find a configuration with a force chain spanning the width of the system, thus causing it to jam; and at , 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 configuration that causes it to jam. The statistical weight of such jamming configurations decreases, and hence the average time required to jam increases, as one either decreases , or increases OHern2 (). In the limit , we expect jamming will occur in finite time only for . In Fig. 1b we show a log-log plot of vs , using a value . We see that the data in the unjammed state is well approximated by a straight line of slope , giving in agreement with the expectation that point J is a second order phase transition.

Figure 1: (color online) a) Plot of inverse shear viscosity vs volume density for several different numbers of particles , at constant small applied shear stress . As increases, one see jamming at a limiting value of the density . b) Log-log replot of the data of (a) as vs , with . The dashed line has slope indicating the continuous algebraic vanishing of at with a critical exponent .

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 obeying a typical scaling law,


Here is the crossover scaling variable, is the crossover scaling critical exponent, and , are the two branches of the crossover scaling function for and respectively.

In Fig. 2 we show a log-log plot of inverse shear viscosity 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 negligible. We use for and for . Again we see that separates two limits of behavior. For , is convex in , decreasing to a finite value as . For , is concave in , decreasing towards zero as . The dashed straight line, separating the two regions of behavior, indicates the power law dependence that is expected exactly at (see below). Similar power law behavior at was recently found in simulations of a three dimensional granular material Sasa ().

Figure 2: (color online) Plot of inverse shear viscosity vs applied shear stress for several different values of the volume density . The dashed line represents the power law dependence expected exactly at and has a slope . Solid lines are guides to the eye. Points labeled correspond to densities , , , , and .

In Fig. 3 we replot the data of Fig. 2 in the scaled variables vs . Using , (the same values used in Fig. 1b) and , we find an excellent scaling collapse in agreement with the prediction of Eq. (6). As the scaling variable , constant; this gives the vanishing of at . As , both branches of the scaling function approach a common curve, , so that precisely at , as nonlinear (). This is shown as the dashed line in both Figs. 3 and 2. A similar scaling collapse of has been found in simulations Berthier () of a sheared Lennard-Jones glass, as a function of temperature and applied shear strain rate , but only above the glass transition, . By comparing the goodness of the scaling collapse as parameters are varied, we estimate the accuracy of the critical exponents to be roughly, and .

Figure 3: (color online) Plot of scaled inverse viscosity vs scaled shear stress for the data of Fig. 2. We find an excellent collapse to the scaling form of Eq. (6) using values , and . The dashed line represents the large asymptotic dependence, . Data point symbols correspond to those used in Fig. 2.

That the crossover scaling exponent , 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 transition at point J (i.e. ). The nature of jamming at finite will be determined by the behavior of the branch of the crossover scaling function , that describes behavior for . From Fig. 3 we see that is a decreasing function of . If vanishes only when , then Eq. (6) implies that vanishes for only when , and so there will be no jamming at finite . If, however, vanishes at some finite , then will vanish whenever ; there will then be a line of jamming transitions emanating from point J in the plane given by the curve . If vanishes continuously at , jamming at finite will be like a second order transition; if jumps discontinuously to zero at , it will be like a first order transition. Such a first order like transition has been reported in simulations Berthier (); Varnik () of sheared glasses at finite temperature below the glass transition, . However, recent simulations Ning () of a granular system at , , showed that a similar first order like behavior was a finite size effect that vanished in the thermodynamic limit. With these observations, 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 transition, and thus implies that there ought to be a diverging 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 Berthier2 (), and the density driven jamming transition Durian2 (). Here we consider the equal time transverse velocity correlation function in the shear driven steady state,


where is the instantaneous velocity in the direction, transverse to the direction of the average shear flow, for a particle at position . The average is over particle positions and time. In the inset to Fig. 4 we plot vs for three different values of at fixed and number of particles . We see that decreases to negative values at a well defined minimum, before decaying to zero as increases. We define to be the position of this minimum. That , indicates 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 Tanguy () in correlations of the nonaffine displacements of particles in a Lennard-Jones glass, in response to small elastic distortions.

Figure 4: (color online) Inset: Normalized transverse velocity correlation function vs longitudinal position for particles, applied shear stress , and volume densities and . The position of the minimum determines the correlation length . Main figure: Plot of scaled inverse correlation length vs scaled shear stress for the data of Fig. 2. We find a good scaling collapse using values , (the same as in Fig. 3) and . Data point symbols correspond to those used in Fig. 2.

As with viscosity, we expect the correlation length to obey a scaling equation similar to Eq. (6). We consider here the inverse correlation length , which like should vanish at the jamming transition, obeying the scaling equation,


The correlation length critical exponent is , but the crossover exponent remains the same as for the viscosity.

In Fig. 4 we plot the scaled inverse correlation length, vs the scaled stress, . Using and , as was found for the scaling of , we now find a good scaling collapse for by taking the value . By comparing the goodness of the collapse as is varied, we estimate . From the scaling equation Eq. (8) we expect both branches of the scaling function to approach the power law as , so that as at nonlinear (). 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. Drocco () for the flowing phase, . It also agrees with found by O’Hern et al. OHern2 () from a finite size scaling argument. Wyart et al. Wyart () have proposed a diverging length scale with exponent by considering the vibrational spectrum of soft modes approaching point J from the jammed side, . While our results cannot rule out , our scaling collapse in Fig. 4 does seem somewhat better when using the larger value .

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. Schrö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 variables. In our case, the most important correction would be to replace in Eq. (6) by ; this could lead to noticeable corrections to our scaling equation near . However, since we find , our conclusion that at 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).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description