Main

At low temperature, high density and low driving, the constituent particles in supercooled liquids1, dense colloids2,3 and granular packings4,5,6, respectively, are nearly locked into a single disordered configuration in which motion is heterogeneous in space and time. Dynamics in these systems may be governed by proximity to a generic ‘jamming transition11’, beyond which rearrangements cease and the viscosity diverges. Key features of spatially heterogeneous dynamics (SHD) on approach to the transition include unusual correlations12 in which particles move in one-dimensional paths (‘strings’13) that aggregate into clusters14, and dynamical correlations as measured by a dynamic four-point susceptibility χ4 (refs 15, 16, 17). Clusters of strings arise naturally in dynamic facilitation18 theory and the random first-order theory of glasses19; their shape reflects the fractal nature of dynamical motion in these systems. Strings are also a crucial ingredient in a recent theory of liquid dynamics near the glass transition20.

Recent studies demonstrate that close-packed granular systems under slow shear exhibit SHD as well4,5,6, bolstering speculation that liquids and granular matter share dynamical similarities on approach to the jammed state. However, the universality of the jamming hypothesis has not yet been tested in terms of variation in the hallmark dynamical heterogeneities as a function of the control parameter. Here, we present the first simultaneous measurements in any experimental system of the growth of the cluster correlation length, string length, four-point correlation length and their characteristic timescales by varying the control parameter. We show that the SHD observed in a far-from-equilibrium, athermal system of air-fluidized granular beads is essentially indistinguishable from that observed in thermal systems such as supercooled liquids and dense colloidal suspensions. Moreover, we show that theoretical models developed for the glass transition can be used to describe our granular system, and predict a mode-coupling-like transition and, more importantly, the jamming transition packing fraction, known as point J (refs 9, 10), from quantities characterizing SHD.

We characterize the spatiotemporally heterogeneous nature of dynamics in an athermal, far-from-equilibrium system of air-driven steel spheres on approach to jamming. Compared with sheared or shaken granular systems, in which energy is injected at the boundaries, in air-driven systems the energy input is uniform in space and time. Our granular system consists of a 1:1 bidisperse mixture of steel beads of diameters ds=0.318 cm and dl=0.397 cm, with respective masses of 0.130 and 0.266 g, confined to a circular region of diameter 17.7 cm. The packing density is varied from an area fraction of φ=0.597 to φ=0.773 by changing the total number of beads from 1,470 to 1,904. Bead motion is restricted to rolling within a horizontal plane, and is excited by an upward flow of air at a fixed superficial flow speed of 545 cm s−1. Bead positions are identified by reflecting light from their chrome surface to a camera three feet above. The duration of the experimental runs is 20 min. In contrast with the molecules in a supercooled liquid, here the particles are macroscopic objects driven at random by a continuous input of energy. Consequently, the speed distributions are non-maxwellian, and the average kinetic energies of the two bead species are unequal. Nevertheless, as reported previously21, the system mimics a simple liquid for low φ and exhibits tell-tale changes in the average structure and dynamics at increasing packing densities.

Dynamical characteristics for an example case, φ=0.773, are shown in Fig. 1a–d. The mean-squared displacement versus time interval (delay time), Fig. 1a, averaged over all start times and all beads, is perhaps the most familiar quantity. It shows ballistic motion at short times, diffusive motion at late times and a plateau of sub-diffusive motion at intermediate times that is indicative of caging. Although informative, the mean-squared displacement cannot distinguish uniform from heterogeneous dynamics. For this, we carry out three measurements developed for supercooled liquids. The first involves clusters of beads that are ‘mobile’, that is, which have displacements ranking among the top 10% of all bead displacements in a given delay time22 (in this system, 10% gives the largest distinction between mobile and immobile beads at all packing densities). The configuration shown in Fig. 1e, where beads are colour-coded according to mobility, demonstrates that the mobile beads are not distributed at random; rather, mobile and immobile beads are clustered and spatially separated, indicating spatially heterogeneous dynamics. One measure of SHD is thus the average size of mobile clusters, Sc(t), defined as the average number of neighbouring mobile beads for a given time interval t. The motion within a mobile cluster, shown by the displacement vectors in Fig. 1e, tends to be correlated into quasi-one-dimensional paths called ‘strings13’. Thus, a second measure of SHD is the average string length, Ls(t), defined in terms of the number of beads that, within a time interval t, replace the initial positions of neighbouring beads to within a tolerance of 0.3ds (ref. 13). A third measure of SHD may be constructed from a four-point susceptibility χ4(t), which measures the extent to which the dynamics at any two points in space are correlated within a time interval15. The self-contribution χ4SS(t) dominates the general result16 and is computed from the variance of the self-overlap order parameter qS(t), which decays from one to zero: χ4SS(t)=N[〈qS(t)2〉−〈qS(t)〉2], where , w=1(0) if |ri(t)−ri(0)|<(>)0.5ds, N is the total number of beads and ri(t) is the position of bead i at time t; averages are taken over all beads and over all start times.

Figure 1: Bead dynamics at area fraction φ=0.773 as a function of delay time.
figure 1

a, Mean-square displacement 〈r2(t)〉. b, Number-average mobile bead cluster size Sc(t), c, Number-average string length Ls(t). d, Self-contribution to the four-point susceptibility χ4SS(t). e, An instantaneous bead configuration where the colour of the beads indicates the mobility over a time interval of 12 s (the timescale for both maximum cluster size and string length). The 10% most mobile beads are red; note that they form clusters. Beads moving in strings have vectors superimposed to indicate their directional motion. Note that the dynamics are spatially heterogeneous.

Full size image

The example results in Fig. 1b–d for the cluster size Sc(t), the string length Ls(t) and the four-point susceptibility χ4SS(t), all exhibit well-defined peaks as a function of time interval, as found in glass-forming liquids. The locations of the peaks indicate the time interval over which the dynamics are most heterogeneous, and the heights of the peaks indicate the spatial extent or ‘strength’ of the heterogeneities. As with glass-forming liquids14,23,24,25 and colloids3, the cluster size and string length are largest at the crossover between caged and diffusive motion, whereas χ4SS(t) (and χ4(t)) peaks later, in the so-called alpha or structural relaxation regime16,26. The athermal air-fluidized beads therefore exhibit spatially heterogeneous dynamics that is identical to thermal glass-forming systems with respect to these three measures.

Now that spatially heterogeneous dynamics is established for gas-fluidized beads, we turn to its variation as a function of control parameter. The distribution of cluster sizes at the peak time interval, shown in Fig. 2a for three different packing densities φ, approaches a power law P(n)∼n−τ,τ=1.5±0.3 as φ is increased. This is consistent with the percolation of mobile bead clusters; similar power laws have been observed in colloids3 and in simulations of supercooled liquids14,23,24 near the mode-coupling temperature. Furthermore, the distribution of string lengths at the peak time interval, shown in Fig. 2b, is exponential, P(l)∼exp(−l/l0), at all values of φ, where l0 is set by the average string length. This is consistent with behaviour reported in simulations of several supercooled liquids13,24,27. We note that the average cluster size is not much larger than the average string length, although the largest clusters observed (∼100 particles) are substantially larger than the largest string observed (∼30 particles) (not shown).

Figure 2: Size distributions of clusters and strings.
figure 2

a,b, Distribution of cluster sizes (a) and string lengths (b) for three values of control variable φ. The cluster size distribution approaches a power law with increasing φ (the solid lines indicate fits by the power law multiplied by an exponential cutoff), whereas the string length distribution is exponential.

Full size image

Results for Sc(t),Ls(t) and χ4SS(t) versus t are shown in Fig. 3a–c for a sequence of different packing densities φ. When beads are added to the system, the average effective temperature also decreases, resulting in a trajectory in the (φ,Teff) phase diagram that heads towards point J, the zero-temperature jamming transition previously found for this system at φ=0.83, which is coincident with the packing density at which the system is random close-packed. As the motion becomes more restricted, the peaks in all three measures of SHD grow and move to later times. Therefore, the dynamics not only slow down but also become more heterogeneous on approach to point J. As the SHD functions have approximately the same shape when viewed on a log–log plot (see data collapse in insets of Fig. 3a,b), this behaviour is fully characterized by the φ-dependence of the characteristic or peak timescales EquationSource math mrow mo{ msubsup mit mrow msub miS mtextc mo* mo, msubsup mit mrow msub miL mtexts mo* mspace mtextand mspace msubsup mit mrow msub miχ mn4 mo* mo} and length scales EquationSource math mrow mo{ msub miξ mrow msub miS mtextc mo( msubsup mit mrow msub miS mtexts mo* mo) mspace mo, mspace msub miχ mrow msub miL mis mo( msubsup mit mrow msub miL mtexts mo* mo) mspace mtextand mspace msub miξ mrow miχ mtext4 mo( msubsup mit mrow miχ mtext4 moχ mo) mo} . The length EquationSource math mrow msub miξ mrow msub miL mtexts mo( msubsup mit mrow msub miL mtexts mo* mo) mi∝ msub mit mrow msub miL mtexts is a correlation length for string-like motion, EquationSource math mrow msub miξ mrow msub miS mtextc mo( msubsup mit mrow msub miS mtextc mo* mo) mi∝ msub miS mtextc is a correlation length of mobile particle clusters and EquationSource math mrow msub miξ mrow msub miχ mtext4 mo= msup mrow mo( msubsup miχ mn4 mrow mtextSS mo( msup mit mo* mo) mo/ miπ mo) mrow mn1 mo/ mn2 of clusters of caged particles. On approach is a correlation length6 of clusters of caged particles. On approach to point J, both the characteristic times and the correlation lengths seem from Fig. 3 to grow without bound. This is reminiscent of behaviour for supercooled liquids as the temperature is lowered. Though very different, both types of systems seem to approach an unusual critical point where the growing lengthscale is purely dynamical, such that there is no macroscopic change in instantaneous structure1,28,29,30.

Figure 3: Variation of SHD as a function of control variable φ.
figure 3

a–c, Cluster size (a), string length (b) and self-contribution to the four-point susceptibility (c), all as a function of time interval t, for a sequence of area fractions φ. In order of increasing peak height: φ=0.597, 0.633, 0.647, 0.669, 0.693, 0.722, 0.742, 0.762, 0.773. The insets in a and b show collapse of all data sets on scaling by the peak heights and peak times.

Full size image

To further quantify this analogy, the growth of the characteristic timescales and dynamical length scales is shown in Fig. 4a–d as a function of packing density. Motivated by recent studies31,32 predicting a power-law divergence of dynamical length scales from mode-coupling theory7 (MCT), as well as earlier applications of MCT to liquids and colloids, we fit all data to a power law of the form 1/|φ/φc−1|x, where both φc and x are adjustable parameters. As seen in Fig. 4a,b, excellent fits are obtained to all data for a single value φc=0.79±0.02. This value of φ lies well above the onset of caging and is less than the jamming packing fraction, in analogy with well-established findings that the mode-coupling temperature is below the caging transition but above the glass-transition temperature7, and demonstrates for the first time a mode-coupling-like transition in a granular system. In addition to MCT, the glass transition can also be well described by the Vogel–Tammann–Fulcher (VTF) equation8; therefore, we also fit the characteristic time and length scales to the form exp(E/|φ/φ0−1|), where E and φ0 are adjustable parameters. As seen in Fig. 4c,d, excellent fits are obtained to all data for a single value φ0=0.84±0.02, consistent with the value of random close-packing for the bead system and the value of point J. As random close-packing is the densest random packing possible and the point at which all motion ceases, a VTF packing fraction of φRCP is analogous to an effectively zero-temperature ideal glass transition, consistent with the definition of point J. This is the first prediction of point J in a granular system from analysis of spatially heterogeneous dynamics.

Figure 4: Dependence of dynamic timescales and length scales on packing density φ.
figure 4

a, Log–log plot of dynamical correlation lengths versus 1/|φ/φc−1| fitted with x=0.40 for EquationSource math mrow msub miξ mrow miS mtextc , x=0.15 for EquationSource math mrow msub miξ mrow miL mtexts and x=1.70 for EquationSource math mrow msub miξ mrow miχ mtextc . b, Log–log plot of characteristic times t* versus 1/|φ/φc−1| with x=1.3 for EquationSource math mrow msubsup mit mrow msub miS mtextc mo* ,x=1.1 for EquationSource math mrow msubsup mit mrow msub miL mtexts mo* and x=1.0 for EquationSource math mrow msubsup mit mrow msub miχ mtext4 mo* . c, Semilog plot of dynamical correlation lengths versus 1/|φ/φ0−1| fitted with E=0.10 for EquationSource math mrow msub miξ mrow miS mtextc , E=0.03 for EquationSource math mrow msub miξ mrow miL mtexts and E=0.42 for EquationSource math mrow moξ mrow msub miχ mtext4 . d, Semilog plot of all t* values versus 1/|φ/φ0−1| fitted with E=0.27. See the text for fit expressions. The error bars in a and c represent the standard error of the measurement. The error bars in b and d represent the uncertainty in determining the time at which the function is maximum.

Full size image

Our study implies that the behaviour of jammed systems, both thermal and athermal alike, may be understood using the theoretical tools developed for liquids. This, in turn, highlights the importance of packing in the underlying physics of the glass transition and jamming. Our results open the door to future theoretical insight into the relationship between granular materials and supercooled liquids, which might be united by a unified theory of jamming.