# Topological persistence and dynamical heterogeneities near jamming

###### Abstract

We introduce topological methods for quantifying spatially heterogeneous dynamics, and use these tools to analyze particle-tracking data for a quasi-two-dimensional granular system of air-fluidized beads on approach to jamming. In particular we define two overlap order parameters, which quantify the correlation between particle configurations at different times, based on a Voronoi construction and the persistence in the resulting cells and nearest neighbors. Temporal fluctuations in the decay of the persistent area and bond order parameters define two alternative dynamic four-point susceptibilities, and , well-suited for characterizing spatially-heterogeneous dynamics. These are analogous to the standard four-point dynamic susceptibility , but where the space-dependence is fixed uniquely by topology rather than by discretionary choice of cutoff function. While these three susceptibilities yield characteristic time scales that are somewhat different, they give domain sizes for the dynamical heterogeneities that are in good agreement and that diverge on approach to jamming.

###### pacs:

45.70.n, 47.55.Lm, 64.70.Pf, 02.40.Pc## I Introduction

The dynamics in both thermal and driven systems become increasingly heterogeneous Ediger (2000); Glotzer (2000); Cipelletti and Ramos (2005) on approach to jamming, where a control parameter such as temperature, packing fraction, or driving rate is varied to the point that all rearrangements cease Cates et al. (1998); Liu and Nagel (2001). Such spatially-heterogeneous dynamics (SHD) may take the form of intermittent strings or swirls of neighboring particles that follow one another in correlated fashion from one packing configuration to another. Observations of SHD have been reported for coarsening Durian et al. (1991); Mayer et al. (2004) and sheared Gopal and Durian (1995) foams, for dense colloidal suspensions Marcus et al. (1999); Kegel and van Blaaderen (2000); Weeks et al. (2000), and for driven granular systems Pouliquen et al. (2003); Marty and Dauchot (2005). Quantification of SHD may be accomplished by a dynamic four-point susceptibility, , defined in terms of density correlations across intervals of both space and time Lacevic et al. (2003). Bounds on have been inferred from standard measurements of dielectric susceptibility for glass-forming liquids and of dynamic light scattering for dense colloidal suspensions Berthier et al. (2005). Direct measurements of have been reported for quasi-two-dimensional granular systems, where particle motion was excited by shear Dauchot et al. (2005) and by a fluidizing upflow of air Keys et al. (2007). In the former, was measured vs space and time at a fixed packing fraction; in the latter, was measured vs time at fixed for several packing fractions. In both Refs. Dauchot et al. (2005); Keys et al. (2007), the dynamic susceptibilities are based on temporal fluctuations of an overlap order parameter; however, the former employs a Gaussian cutoff function while the latter uses a step function.

In this paper we conduct further measurements and analyses of spatially-heterogeneous dynamics (SHD) in a system of air-fluidized beads similar to Ref. Keys et al. (2007). Our primary extension of Ref. Keys et al. (2007) concerns the cutoff function in the overlap order parameter. First we explore the length-dependence of by varying the cutoff distance . Second, we remove the arbitrariness in choice of cutoff function by using topological measures of overlap. In particular we build upon the concept of persistence, introduced to help characterize the coarsening of spin domains in magnets or of bubbles in foams Bray et al. (1994); Tam et al. (1997). As a topological measure of overlap for particulate systems, we define a normalized, dimensionless persistent area as the fraction of space that remains inside the same Voronoi cell after a given time interval. We also introduce a second topological measure of configurational overlap in terms of persistent bonds, i.e. the fraction of nearest-neighbor pairs that remain neighbors after a given time interval. While both persistent area and persistent bond order parameters vanish with time, the latter can exhibit a much slower decay if a string of particles follow one another over a long distance. After discussion of such cutoff functions, we use the resulting four-point susceptibilities both to explore the connection of SHD to local structure and to observe the evolution of SHD on approach to jamming. Finally we consider possible artifacts due to finite experiment duration.

## Ii Methods and Background

The granular system under study here is a quasi-two-dimensional monolayer of spheres that roll without slipping due to a fluidizing upflow of air, as in Refs. Abate and Durian (2005, 2006); Keys et al. (2007). The spheres are steel with an equal number of cm and cm diameter sizes. This is the same system as in Ref. Keys et al. (2007), though here most run durations are 120 minutes rather than just 20 minutes. The steel beads are confined in a circular flat sieve with a diameter of 17.7 cm and a mesh size of 150 m. Energy is injected by a uniform flow of air through the sieve at a flux which is high enough to randomly drive the balls by turbulence () without causing levitation. The flux is increased from 560 cm/s to 700 cm/s as the area fraction is decreased from to , to ensure that the low density systems are fluidized. Three feet above the sieve are six incandescent lights aligned in a ring one foot in diameter. At the center of the ring is a 120 Hz Pulnix 6710 video camera with a array of 8BIT square pixels and a zoom lens. Light reflects specularly off the tops of the beads and is imaged by the camera. Bright spots are tracked so that the positions are known for all beads and for all times during the entire experiment. Further details are available in Refs. Abate and Durian (2005, 2006).

By increasing the fraction of area occupied by spheres, at fixed gas flux, this system exhibits a jamming transition at point-J, such that the average bead kinetic energy vanishes linearly on approach to random close packing at about Abate and Durian (2006). Near this transition, the bead dynamics exhibit the usual sequence of crossovers from ballistic, to sub-diffusive, to diffusive motion as a function of delay time. Spatially-heterogeneous dynamics occur in the sub-diffusive regime, particularly at delay times such that the beads are breaking out of their cages and beginning to diffuse. This is illustrated graphically in Fig. 1(a-c), which show example average velocity vector fields, , for three different delay times ; the area fraction and delay time for each image is chosen to highlight motion in each of the three regions of phase space, as specified in Fig. 1d. For both very short and very long delay times, corresponding to ballistic and diffusive regimes, the magnitude and direction of the average velocities for neighboring beads are completely uncorrelated. However, at intermediate delay times, Fig. 1b, corresponding to subdiffusive motion, a few string-like clusters of neighboring beads move in swirls while other beads move relatively little. As time evolves, such swirls come and go in different regions of the sample so that the dynamics become completely ergodic. This is best seen by a video of the average velocity field Abate and Durian (2007).

Now we begin the central task of this paper, which is to quantify the spatially-heterogenous nature of the dynamics depicted qualitatively in Fig. 1. For this, a standard tool is the four-point dynamic susceptibility , as reviewed in Ref. Lacevic et al. (2003). This may be approximated from temporal fluctuations in the instantaneous self-overlap order parameter, defined here as

(1) |

where is the number of beads and is a cutoff function that equals 1 if the displacement of bead remains less than across the whole time interval and that equals 0 otherwise. Whereas here and in Ref. Keys et al. (2007) is a step-function cutoff, in Ref. Dauchot et al. (2005) it is a Gaussian. In either case the average self-overlap order parameter and the four-point susceptibility are given by moments of averaged over all times as follows:

(2) | |||||

(3) |

Note that is independent of system size because the variance scales as . Also note that grows in proportion to the variance, or heterogeneity, in the displacements of different beads across a given time interval.

Example results for the space and time dependence of and are shown in Fig. 2, for one particular area fraction . In particular, these two functions are plotted vs delay time for several different values of . By construction decays from one to zero as a function of . For increasing , the observed decay becomes somewhat sharper, but more importantly the location of the decay moves to longer delay times, since the beads take longer to move a larger distance. By contrast the susceptibility is a peaked function of , which vanishes at both early and late times. By comparison of Figs. 2a-b it is evident that the peak of is located at the characteristic decay time of . Furthermore, the height of the peak depends on , such that it vanishes for large and small and becomes a maximum for on the order of the bead size. Thus in Ref. Keys et al. (2007) the overlap threshold was set equal to the radius of the small beads, , and the dynamic susceptibility was considered only as a function of time, .

Results for and vs are displayed in Figs. 3a-b for a sequence of area fractions. This reproduces Fig. 3c of Ref. Keys et al. (2007), with the identical system but with a longer run duration. As seen before, the location of the peak in moves to longer times on approach to jamming. Furthermore, the peak height increases concurrently. Thus the dynamics become not only slower but more spatially heterogeneous on approach to jamming. The size of the spatial heterogeneities represents a dynamical correlation length, , which diverges at jamming; such behavior was first reported in Keys et al. (2007), and will be discussed further in Section IV. Note that for all area fractions, the value of the overlap order parameter is nearly constant, , when the susceptibility reaches its peak.

The dynamical correlation length may be deduced from the value of the self-overlap order parameter and the dynamic susceptibility at its peak, using the following physical picture. For a system of total beads, we wish to know the typical number of heterogeneities and more importantly the typical number of beads in each of these fast-moving domains. Since the contribution to the self-overlap order parameter is from the highly-mobile beads inside domains, and is from the other less-mobile beads, the average self-overlap order parameter at is . However the instantaneous value of varies with time because of counting statistics in the number of heterogeneous domains actually present at any instant. Therefore, the peak dynamic susceptibility is . By eliminating in these expressions for and , the typical number of particles in each fast-moving dynamic heterogeneity is found to be

(4) |

Since is nearly constant, as observed in Fig. 3, the number of beads per domain is directly proportional to the peak susceptibility . If the domains are compact, then the dynamical correlation length scales as , where is the bead diameter and is dimensionality Dauchot et al. (2005); Keys et al. (2007). If the domains are strings, then the dynamical correlation length is . If fluctuations in domain size are also included, then the left-hand side of Eq. (4) is slightly modified to .

To our knowledge, the counting arguments leading to Eq. (4) have not been previously published. While it seems common knowledge that the peak susceptibility depends upon domain size, we have found no arguments or calculations in prior literature to support or quantify the relation.

## Iii Topological Persistence

The dynamic susceptibility discussed in the previous section relies upon a somewhat arbitrary choice for a cutoff function, , in Eq. (1). Furthermore, since is the more interesting parameter, it also relies upon a particular choice of cutoff length, , so that spatially-heterogeneous dynamics may be studied as a function of delay time. In this section we propose two alternative dynamic susceptibilities, intended for the same purpose, in which the role of the cutoff function and cutoff length are determined uniquely by topological considerations. Both are based upon the Voronoi tessellation constructed from the particle positions at each instant of time.

### iii.1 Persistent area

An order parameter for quantifying the degree of overlap between particle configurations may be constructed by considering the fraction of the sample which remains inside the same Voronoi cell across a time interval . In particular, we define the instantaneous persistent volume as

(5) |

where is the sample volume, if the point remains inside the same Voronoi cell across the whole time interval , and if the point becomes enclosed by another Voronoi cell. Since our experiments are in two-dimensions, we shall refer to Eq. (5) as the instantaneous persistent area. By construction, it decays monotonically from one to zero as the beads move a distance comparable to their size. Thus the persistent area is similar to the previous self-overlap order parameter with a cutoff length set by bead size.

A graphic illustration of the contribution of each point in the sample to the average over space given by Eq. (5) is shown in Fig. 4 for a sequence of increasing delay times . Pixels are colored black for and are colored white otherwise; therefore, the sample is all black for and it becomes progressively white for increasing . As illustrated in Fig. 4, white regions first appear at short near the boundary of the initial Voronoi cells. At larger the white regions thicken as the persistent black regions at the cell cores shrink. Eventually each black spot vanishes. The graphics here can be compared with results for a coarsening foam in Ref. Tam et al. (1997), where domains are defined physically by actual bubbles rather than by a Voronoi construction. For dynamics that are spatially uniform, the persistent area vanishes uniformly at the same rate throughout the whole sample, with only random uncorrelated variation between neighboring regions. This is not the case in Fig. 4, where the spatially-heterogeneous nature of the dynamics is evident in the long swaths of white caused by chains of fast-moving beads. This is most pronounced for delay times between about 5 and 50 s, and corresponds closely with structure in the average velocity field seen for example in Fig. 1.

The average persistent area and a dynamical susceptibility quantifying spatial heterogeneity may now be defined, in analogy with the previous section, by moments of averaged over all times :

(6) | |||||

(7) |

The variance is multiplied by the number of beads in the sample, so that is independent of system size. As an example, the decays of vs plotted in Fig. 5a display considerable variation for different choices of starting time . The resulting susceptibility, , plotted underneath in Fig. 5b, displays a peak as expected near the characteristic decay time for the average persistent area. The height of this peak, , and the corresponding average persistent area , give the typical number of beads in the dynamic heterogeneities as

(8) |

by the same arguments giving Eq. (4). Here, is the contribution to the persistent area by the beads in the fast-moving heterogeneities and is the contribution from the remaining less mobile beads. Based on graphics like Fig. 4, we estimate and . To remove this uncertainty, an alternative approach might be to enforce and by defining an overlap order parameter such that the contribution from each bead is a step function that drops from one to zero when all its persistent area vanishes.

The behavior of and are displayed in Fig. 6 for different bead packing fractions, , using the same sequence of values and the same symbol/color codes as in Fig. 3. The decay time of and the peak position of coincide and become longer for increasing on approach to jamming. For all area fractions, the value of the persistent area is nearly constant, , when the susceptibility reaches its peak. The peak height also increases as the dynamics slow down, indicative of a diverging dynamical correlation length. This behavior strongly parallels that in Fig. 3 based on an overlap order parameter with a particular step-function cutoff. In fact, the ratio of peak heights is nearly constant for all area fractions: . The only striking qualitative difference is that exhibits a two-step decay at high , whereas does not. This is because the persistent area begins to decay at very short delay times, when the beads experience small-displacement ballistic motion, rattling within a “cage” of unchanging nearest neighbors. The shoulder in corresponds to the crossover from ballistic to subdiffusive motion observed in the mean-squared displacment. For a typical cutoff comparable to bead size, this so-called relaxation is not detected by the traditional overlap order parameter , which decays only due to the cage-breakout -relaxation process. Furthermore, no shoulder is evident in for any choice of in Fig. 2. Besides being naturally and uniquely defined by topology, the persistent area order parameter thus has the added advantage of being able to capture the two-step and relaxation processes characteristic of glass-forming systems.

### iii.2 Persistent bond

A second order parameter may be defined naturally from the set of nearest neighbors, or “bonds”, specified by the Delaunay triangulation dual to the Voronoi construction. In particular, we define the instantaneous persistent bond number as the fraction of all bonds at time that remain unbroken across the interval . The progressive breaking of bonds throughout the sample is illustrated in Fig. 7, at the same area fraction as in the persistent area fraction diagrams of Fig. 4. By comparison, the persistent bond number decays more slowly, due to neighboring beads that move a large distance while remaining next to one another. Also by comparison, the string-like swirls of the dynamic spatial heterogeneities are less evident by casual inspection of the persistent bond diagrams than for the persistent area diagrams. Nevertheless, bond lifetimes and the spatial heterogeneity of broken bonds have been used to characterize simulations of glassy systems Yamamoto and Onuki (1998).

The average persistent bond number, the related dynamic susceptibility, and the size of the dynamic heterogeneities, respectively, are given in precise analogy with the previous sections as follows:

(9) | |||||

(10) | |||||

(11) |

The trends in this order parameter and susceptibility are displayed vs in Figs. 8a-b for the same sequence of area fractions as in the analogous plots based on self-overlap and persistent area. Here, as before, the order parameter decays more slowly and the peak height increases as the area fraction is increased towards jamming. Similar to , but in contrast with , the average persistent bond number exhibits a one-step decay. By inspection, we estimate and . Note that for all area fractions, the value of the persistent bond is nearly constant, , when the susceptibility reaches its peak. Also, the ratio of peak susceptibility heights is nearly constant: .

## Iv Comparisons

Three dynamic susceptibilities, , , and , have now been discussed for characterizing spatial dynamic heterogeneities in terms of the variance of order parameters based respectively on self-overlap, persistent area, and persistent bond. The nature of these order parameters, and their possible relation with local structure, are compared in Fig. 9 by snapshots of the system at eight different instants in times. Here the time increment is chosen as , half the time delay at which reaches its peak, so that there is a reasonable balance of continuity and evolution between successive snapshots. While each row represents a different time, the first three columns represent local structure and the last three represent local dynamics. In terms of structure, Voronoi cells are shaded according to (a) their number of sides; (b) their circularity shape factor, equal to perimeter-squared divided by times area Moucka and Nezbeda (2005); Reis et al. (2006); Abate and Durian (2006); and (c) the reciprocal of their area, as a measure of the local density. In terms of dynamics, the average velocity vectors and the persistent areas are depicted at and the persistent bonds are depicted at , so that the spatial dynamic heterogeneities are maximally emphasized. By inspection, we note the following points. First, the spatial structure of the dynamics revealed by the average velocity vectors and the persistent areas is clearly related. Both measures give a clear feeling for string-like swirls of neighboring beads that is intermittently excited against a background of less mobile beads. By contrast the persistent bonds give no such feeling of motion. And while there are more broken bonds near the fast moving domain, the correspondence isn’t one-to-one because neighboring mobile beads may remain nearest neighbors. Second, there is no apparent correlation between any of the measures of structure and the location of the spatial dynamic heterogeneities. It is not evident how to predict when or where a dynamic heterogeneity will occur based on usual measures of local structure.

Next we compare the characteristic times and domain sizes for dynamic heterogeneities, as deduced from the three order parameters and susceptibilities. The growth of these scales as a function of increasing bead area fraction are displayed in Fig. 10. Note that the time scales are different for the three susceptibilities. As expected, is the largest while is only slightly larger than . However all three times appear to grow in proportion to one another, suggesting that they probe similar physics. Further reinforcing this point, the domain sizes are indistinguishable for the three susceptibilities. Thus the overlap order parameter with discretionary choice of step-function cutoff at is not as arbitrary as it may seem, since it reveals the same picture of spatially heterogeneous dynamics as the persistent areas and bonds, which are uniquely defined by topology.

Next we consider the functional form of the growth of dynamical heterogeneities on approach to jamming, exactly per Ref. Keys et al. (2007). Thus in Figs. 10a-b, data for the peak times and domain sizes are shown on log-log plots vs . With a single choice of , the data can all be well-fit to a power of where the exponent is for the domain sizes and is for the peak times. Such power-law divergences on approach to below random-close packing are consistent with a mode-coupling theory, where the mode-coupling temperature is above the glass-transition temperature. The value of the correlation length exponent is consistent with prior simulation results based on finite-size scaling analysis of random close packing fraction O’Hern et al. (2003), on the size of the disturbance away from a perturbation Drocco et al. (2005), and on scaling analysis of stress and strain rate Olsson and Teitel (2007). Slightly smaller exponents for the dynamical correlation length have been predicted for the approach from above Wyart et al. (2005); Silbert et al. (2005); Schwarz et al. (2006). The observed time scale exponent is less than that for simulation results based on stress relaxation Durian (1995, 1997), for the same model as in Ref. Olsson and Teitel (2007). However the power-law description is not unique. As for the glass transition in molecular liquids, the growth of characteristic time and length scales can also be well-fit to Vogel-Tammann-Fulcher (VTF) equation, . Thus in Figs. 10c-d, data are shown on semi-log plots vs . With a single choice of , which corresponds to random close packing, the data follow the VFT equation where for the peak times and for the domain sizes. Altogether, as reported in Ref. Keys et al. (2007), the growth of dynamical time and length scales on approach to point-J in the system of air-fluidized beads happens in close quantitative analogy to the glass transition in liquids.

Before closing we make one final comparison of the four-point dynamic susceptibility for runs of different duration. In particular, we compute peak heights for subsets of duration of the full 120 minute runs and normalize by the long-duration value . Results are plotted in Fig. 11 vs duration time, normalized by the long-duration peak time , for all area fractions considered above. Note that the normalization causes all data to collapse, even though peak heights and times grow on approach to jamming. Also note that the long-duration value of the peak height is attained only if the run is about ten times longer than the peak time . Intuitively, if a run is not very long compared to the decay time of the overlap order parameter, then too few configurations are sampled and the variance is underestimated. This effect will limit the proximity to which point-J may be approached. For example, if run duration is held fixed and the area fraction is gradually increased, then the peak height and the corresponding dynamical correlation length will initially grow but will eventually appear to decrease as grows beyond about one-tenth of the run duration. This can be seen in the inset of Fig. 11. To eliminate such erroneous finite-time artifacts, the experimental observation time must be at least a decade longer than the time it takes for the overlap order parameter to vanish. Conceivably, an alternative might be to take a smaller cutoff in the self-overlap order parameter or to alter the proportionality constant relating susceptibility and domain size according to the experimental conditions.

## V Conclusion

Spatially-heterogeneous dynamics may be characterized successfully using four-point susceptibilities based on many choices for the configurational order parameter. For a self-overlap order parameter with the seemingly ad-hoc choice of a step-function cutoff at , the characteristic dynamical correlation length for the size of intermittent mobile domains is indistinguishable from those based on persistent areas and bonds in a unique topological description. This agreement requires a physical picture for how to deduce domain size from the peak of the susceptibility, as in Eqs. (4,8,11). Nevertheless, some differences still exist in the detailed form of the three susceptibilities. For example the peak times all differ, though they always remain in constant proportion to one another. At early times, below the peak, the persistent area is the first to begin its decay since self-overlap and persistent bonds are insensitive to the small ballistic motion of grains as they rattle in a cages of fixed nearest neighbors. At late times, past the peak, the persistent bond is the last to decay since a few neighboring beads remain in close contact as they diffuse through the sample. Thus, while all three susceptibilities give the same picture of the spatially heterogeneous dynamics, they offer different possibilities for characterizing short and long time dynamics. As applied to the change in behavior with increasing bead packing fraction, all three give characteristic time and length scales for the dynamical heterogeneities that appear to diverge in accord with simulation results for supercooled liquids and for dense athermal systems of soft repulsive particles.

###### Acknowledgements.

We thank Maksym Artomov, Sharon Glotzer, Aaron Keys, Andrea Liu, and Steve Teitel for helpful discussions. Our work was supported by the National Science Foundation through grants DMR-0514705 and DMR-0704147.## References

- Ediger (2000) M. D. Ediger, Annu. Rev. Phys. Chem. 51, 99 (2000).
- Glotzer (2000) S. C. Glotzer, J. Non-Cryst. Solids 274, 342 (2000).
- Cipelletti and Ramos (2005) L. Cipelletti and L. Ramos, J.Phys.-Cond. Matt. 17, R253 (2005).
- Cates et al. (1998) M. E. Cates, J. P. Wittmer, J.-P. Bouchaud, and P. Claudin, Phys. Rev. Lett. 81, 1841 (1998).
- Liu and Nagel (2001) A. J. Liu and S. R. Nagel, eds., Jamming and Rheology: Constrained Dynamics on Microscopic and Macroscopic Scales (Taylor and Francis, NY, 2001).
- Durian et al. (1991) D. J. Durian, D. A. Weitz, and D. J. Pine, Science 252, 686 (1991).
- Mayer et al. (2004) P. Mayer, H. Bissig, L. Berthier, L. Cipelletti, J. P. Garrahan, P. Sollich, and V. Trappe, Phys. Rev. Lett. 93, 115701 (2004).
- Gopal and Durian (1995) A. D. Gopal and D. J. Durian, Phys. Rev. Lett. 75, 2610 (1995).
- Marcus et al. (1999) A. H. Marcus, J. Schofield, and S. A. Rice, Phys. Rev. E 60, 5725 (1999).
- Kegel and van Blaaderen (2000) W. K. Kegel and A. van Blaaderen, Science 287, 290 (2000).
- Weeks et al. (2000) E. R. Weeks, J. C. Crocker, A. C. Levitt, A. Schofield, and D. A. Weitz, Science 287, 627 (2000).
- Pouliquen et al. (2003) O. Pouliquen, M. Belzons, and M. Nicolas, Phys. Rev. Lett. 91, 14301 (2003).
- Marty and Dauchot (2005) G. Marty and O. Dauchot, Phys. Rev. Lett. 94, 015701 (2005).
- Lacevic et al. (2003) N. Lacevic, F. W. Starr, T. B. Schroder, and S. C. Glotzer, J. Chem. Phys. 119, 7372 (2003).
- Berthier et al. (2005) L. Berthier, G. Biroli, J. P. Bouchaud, L. Cipelletti, D. E. Masri, D. L’Hote, F. Ladieu, and M. Perino, Science 310, 1797 (2005).
- Dauchot et al. (2005) O. Dauchot, G. Marty, and G. Biroli, Phys. Rev. Lett. 95, 265701 (2005).
- Keys et al. (2007) A. S. Keys, A. R. Abate, S. C. Glotzer, and D. J. Durian, Nature Physics 3, 260 (2007).
- Bray et al. (1994) A. J. Bray, B. Derrida, and C. Godreche, Europhys. Lett. 27, 175 (1994).
- Tam et al. (1997) W. Y. Tam, R. Zeitak, K. Y. Szeto, and J. Stavans, Phys. Rev. Lett. 78, 1588 (1997).
- Abate and Durian (2006) A. R. Abate and D. J. Durian, Phys. Rev. E 74, 031308 (2006).
- Abate and Durian (2005) A. R. Abate and D. J. Durian, Phys. Rev. E 72, 031305 (2005).
- Abate and Durian (2007) A. R. Abate and D. J. Durian, Chaos ?, to appear (2007).
- Yamamoto and Onuki (1998) R. Yamamoto and A. Onuki, Phys. Rev. E 58, 3515 (1998).
- Moucka and Nezbeda (2005) F. Moucka and I. Nezbeda, Phy. Rev. Lett. 94, 040601 (2005).
- Reis et al. (2006) P. M. Reis, R. A. Ingale, and M. D. Shattuck, Phys. Rev. Lett. 96, 258001 (2006).
- O’Hern et al. (2003) C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E 68, 011306 (2003).
- Drocco et al. (2005) J. A. Drocco, M. B. Hastings, C. J. O. Reichhardt, and C. Reichhardt, Phys. Rev. Lett. 95, 088001 (2005).
- Olsson and Teitel (2007) P. Olsson and S. Teitel, arXiv:0704.1806 (2007).
- Wyart et al. (2005) M. Wyart, S. R. Nagel, and T. A. Witten, Europhys. Lett. 72, 486 (2005).
- Silbert et al. (2005) L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. Lett. 95, 098301 (2005).
- Schwarz et al. (2006) J. M. Schwarz, A. J. Liu, and L. Q. Chayes, Europhys. Lett. 73, 560 (2006).
- Durian (1995) D. J. Durian, Phys. Rev. Lett. 75, 4780 (1995).
- Durian (1997) D. J. Durian, Phys. Rev. E 55, 1739 (1997).