Shear induced rigidity in athermal materials: a unified statistical framework
Abstract
Recent studies of athermal systems such as dry grains and dense, non-Brownian suspensions have shown that shear can lead to solidification through the process of shear jamming in grains and discontinuous shear thickening in suspensions. The similarities observed between these two distinct phenomena suggest that the physical processes leading to shear-induced rigidity in athermal materials are universal. We present a non-equilibrium statistical mechanics model, which exhibits the phenomenology of these shear-driven transitions: shear jamming and discontinuous shear thickening in different regions of the predicted phase diagram. Our analysis identifies the crucial physical processes underlying shear-driven rigidity transitions, and clarifies the distinct roles played by shearing forces and the density of grains.
pacs:
05.50.+q,83.80.Fg,45.70.-n,64.60.-i,83.60.RsI Introduction
Athermal materials such as dry grains and dense non-Brownian suspensions can respond to shear by organizing into structures that support the imposed load Cates et al. (1998): a process that has been termed shear-jamming (SJ) in grains Bi et al. (2011); Zhang et al. (2010); Ren et al. (2013); Sarkar et al. (2013), and discontinuous shear thickening (DST) in suspensions Brown and Jaeger (2013); Seto et al. (2013); Heussinger (2013); Wyart and Cates (2013); Fernandez et al. (2013); Fall et al. (2008, 2009, 2010). The nature of this self-organization process has been intensely investigated in recent months, and striking similarities have been observed between the two transitions. This is remarkable since the SJ transition occurs through a quasistatic process and refers to static states of particles interacting via purely repulsive contact interactions Bi et al. (2011); Ren et al. (2013); Sarkar et al. (2013), and DST occurs through a dynamical process that creates non-equilibrium steady states (NESS) of particles interacting via hydrodynamic and contact interactions Seto et al. (2013); Mari et al. (2014). The single most important trigger for these transitions has been identified as the proliferation of frictional contacts Bi et al. (2011); Ren (2013); Seto et al. (2013); Mari et al. (2014). For a range of packing fractions, , below the isotropic jamming density, , quasi static shearing causes frictional grains to come into contact leading to the SJ transition Bi et al. (2011); Ren et al. (2013); Sarkar et al. (2013). In a similar range of , athermal suspensions exhibit DST as increasing shearing rate leads to a loss of lubrication forces and increasing number of frictional contacts Seto et al. (2013); Wyart and Cates (2013).
Lattice models have a venerable history of identifying the core physical mechanisms driving phase transitions, and finding commonalities between seemingly disparate systems. We have constructed a non-equilibrium, driven, disordered model that focuses on the process of formation and rearrangement of frictional contacts under driving by a field. In contrast to studies that interrogate the microscopic mechanisms leading to the SJ and DST transitions Wyart and Cates (2013); Sarkar et al. (2013), we analyze an effective theory that is built on the premise that the driving field, either strain () or strain rate (), increases shear stress and promotes the formation of frictional contacts. We examine the consequences of the interplay between the driving field and the underlying disorder of the contact network on the development of a robust, force-bearing network. The model focuses solely on the force network: changes in the network of frictional contacts with their associated tangential forces strongly affect the viscosity of suspensions in the DST regime Seto et al. (2013); Mari et al. (2014); Wyart and Cates (2013).
The mapping between the parameters defining the model and and the physical parameters defining and controlling force networks in the SJ and DST transitions are summarized in Tabel 1. In the next section we develop the model starting from a rigorous mapping of grain-level stresses to spins.
Ii Model
ii.1 Rigorous Mapping
The tensor representing the stress state of a grain can be divided into an completely isotropic part that defines hydrostatic pressure and a deviatoric part that represents normal and shear stresses. The deviatoric part can be represented as an element of a vector space Behringer et al. (2014). Illustrating in 2D, the stress tensor of a grain, which is symmetric since the grain is torque balanced, can be written as:
(1) | |||||
(2) | |||||
(3) |
where is the hydrostatic pressure, is the normal stress, and is the shear stress. The deviatoric part of the stress, which excludes this hydrostatic part, is therefore an element of a 2D vector space spanned by two matrices, and Behringer et al. (2014). The components of the vector are the normal stress, , and the shear stress, , and the length, , provides a measure of the stress anisotropy of each grain.
The stress state of a grain is influenced by the local strain arising from the displacement of the neighboring grains. The displacement of the grains comprises of a homogeneous part, which can be characterized by a set of affine transformations and an inhomogeneous part, called non-affine displacement, which cannot be described through a series of affine transformation of the grain coordinates. The non-affine displacements are best characterized by a measure called , first introduced by Falk and Langer Falk and Langer (1998). To calculate , one measures the actual displacements of the grains, and then chooses an optimum affine strain tensor , which minimizes the mean squared deviation of the actual displacement from a homogeneous displacement due to the strain tensor. This minimized deviation is referred to . The measure has been applied to characterize the non-affine displacements in granular experiments Utter and Behringer (2008), and shows that the non-affine strains follow a Gaussian distribution with mean approximately zero. Additionally, it has been found Ren (2013) that the deviatoric stress vectors interact with these non-affine strains in a manner similar to how magnets interact with spins. Thus, the continuous vectors, , can be imagined as continuous spins. The vector sum of these spins maps to the global deviatoric stress tensor, and the magnetization measures the stress anisotropy of the global stress tensor.
ii.2 Mapping to a Spin-1 Ising Model
Though rigorous, analyzing the properties of such a continuous spin model with variable lengths, in the presence of a random field is difficult. We, therefore, use a threshold to map the grain-level stress to a spin 1 Ising model. Let us define . If, , we map the grain to , otherwise we map it to , depending on whether the grain points along or perpendicular to the compressive strain direction. The mapping is illustrated in Fig. 1 for the two-dimensional (2D) SJ system. In such a system, represent grains with two contacts, which have strong stress anisotropy( ) and represent grains with more than two contacts, which have a nearly isotropic stress tensor( ) and connect chain-like force networks. A similar mapping applies to the three-dimensional (3D) DST systems with referring to grains connecting chain-like networks Seto et al. (2013).
We envision the SJ and DST processes as ones where the contact network and grain-level stresses reach a force and torque balanced state in the presence of driving Sarkar et al. (2013); Seto et al. (2013); Mari et al. (2014). We model this by the zero-temperature, single-spin flips, energy minimizing dynamics of the energy function (Fig. 1) Vives and Planes (1994); Vives et al. (1995); Sethna et al. (2004):
(5) |
In SJ experiments, it is known that as the imposed strain, , is increased, the fraction of grains with small stress anisotropy ( in our model) increases. In DST, it is the strain rate that plays the same role. Our viewpoint is that this is the primary effect of or . So we map to for SJ and to for DST, whereas spin flips map to rearrangements of the contact network of particles. The dependent chemical potential incorporates the effect of () on the fraction of grains with small stress anisotropy in SJ (DST). In these systems with shear-induced rigidity, this fraction increases with increasing driving, and we therefore restrict our analysis to increasing functions of . However, the model is more general and can address other scenarios.
As seen in Fig. 1, shearing leads to significant non-affine displacements: displacements that are inhomogeneous, and cannot be described by any type of homogeneous deformation of the unstrained state Utter and Behringer (2008). These are represented by the random magnetic field , at every site. SJ experiments indicate that the distribution of the depends on but evolves little during the shear-jamming process Sarkar et al. (2013); Ren (2013), therefore, we treat the as a quenched random field chosen from a Gaussian distribution with zero mean and standard deviation . In the granular systems, the constraints of mechanical equilibrium introduce effective interactions between the stress tensors of grains. In a force chain where every grain has only two contacts, the anisotropy of the stresses of grains in the chain are highly correlated Sarkar et al. (2013). We model this effective interaction by a ferromagnetic interaction between spins.
We are interested in understanding the effects of shear in creating robust force networks through the introduction of frictional contacts. Our model, therefore, differs from other driven-disordered models in the class of Eq. 5 Sethna et al. (2004) in one crucial respect: the external field controls the average population of sites through a chemical potential (). For the current study, we model the dependence of as , which is the simplest that admits an increase in the concentration of with the magnitude of the driving field. In this work, we focus on the aspects of the model that are relevant to shear-induced rigidity, however, the statistics of avalanches and the yielding behavior exhibit interesting new features, which will be studied in the future.
For spins, we define two global order parameters:
(6) |
Here, corresponds to the fraction of grains with isotropic stress tensors (), and corresponds to the stress anisotropy (contact-stress anisotropy) in SJ (DST). The zero-temperature dynamics samples the metastable states of this disordered model, which we associate with the force networks sampled in the SJ and DST processes. In the SJ context, the history represents a history, while in the DST context, it represents a history. Since there are no thermal fluctuations in our model, averages ( ) are over metastable states corresponding to different realizations of the quenched disorder field, . To simplify notation, we eliminate the symbol in the following.
Starting from a metastable state, if is changed adiabatically such that all spins that can lower their energy by flipping do so, there could be a range of over which the original state is stable. At a certain , however, a threshold is crossed at some site (determined by the and the effective field and that spin changes its state. This, in turn, could lead to the threshold being crossed at other sites, creating a cascade of spin flips in an avalanche Sethna et al. (2004) until a new metastable state is reached.
In the granular context, we envision this exploration of metastable states to correspond to exploration of force networks that are in local mechanical equilibrium under driving. The SJ process is quasistatic and since the non-affine strain field is observed to evolve only weakly over the range of probed by the experiments, there is a clear correspondence between the sampling of metastable states in the model and the force networks in the granular assembly. In DST, however, one studies time averages in the NESS at a given . The correspondence between the ensemble average over and the time average is valid if the NESS dynamically samples non-affine strains with Gaussian statistics, and if the time to reach a force and torque balanced state is much shorter than the relaxation time of the non-affine strain field. These assumptions are validated in simulations Mari et al. (). The adiabatic assumption implies that is ramped up slowly compared to microscopic time scales Brown and Jaeger (2013).
A priori, it is not clear what experimental knob can be turned to tune . However, there are strong arguments presented below, based on comparing predictions of our model to existing experimental and numerical observations, linking increasing to a reduction in . A scaling description of DST has been constructed by invoking a stress-scale dependence of the packing fraction at which the viscosity diverges Wyart and Cates (2013). We relate the stress-scale to , and therefore, in our approach it is the stress scale that is controlled by , through , and by . If the dominant effect of on the force network is through the statistics of the non-affine strain field, then the two approaches should yield similar results. Below, we will establish specific mappings in the SJ and DST regimes by comparing our predictions to experiments, simulations, and the scaling theory.
Iii Results
iii.1 Meanfield Solution of the Model
To solve the spin-1 Ising model under meanfield(MF) approximation, we observe that the order parameters can be represented through the probability of finding a particular value of spin at a particular lattice point:
(7) |
where measures the probability that the spin takes the value ( or 0). Also, in the MF approximation, the energy of a spin is given by:
(8) |
Therefore,
(9) |
In our zero-temperature dynamics, a spin will be in the +1 state if , and . This condition is satisfied if
whence
(10) |
A similar calculation yields:
(11) |
Using these probabilities, and the definition of and (Eq. 7), we obtain:
(12) | |||||
(13) |
Here erf and erfc are the error function and the complementary error function, respectively. Special cases and several important aspects of the MF solution are detailed in the appendix.
iii.2 Meanfield Phase Diagram
Meanfield calculations of and along a forward trajectory, with monotonically increasing (Fig. 2(a)) suffice to illustrate that the phenomenology of both the SJ and DST transitions are realized in the model. The meanfield phase diagram in the space is shown in Fig. 2(b). Increasing corresponds to higher average concentration of zero spins at , corresponding to larger values of at zero driving field. In the SJ system Sarkar et al. (2013), , which maps on to the upper part of the phase diagram in Fig. 2(b). In DST, however, there are a vanishing number of frictional contacts at Seto et al. (2013); Mari et al. (2014), which maps these systems to the lower part of the phase diagram. The to mapping discussed earlier implies that decreases from left to right in Fig. 2(b).
Another parameter that influences the model phase diagram is , the rate of increase of with . As shown in Fig. 7 in the appendix, only the protocols lead to a monotonic increase of for all , a feature of the number of frictional contacts in both SJ and DST. We, therefore restrict our analysis to , and unless otherwise stated the results presented are for . For quantitative comparisons to experiments and simulations, one should obtain by comparing the meanfield predictions for to the increase of () in SJ (DST) systems at different .
The qualitative differences between different regions of the phase diagram are best characterized by , which maps on to the dependence of stress anisotropy on . As shown in Fig. 2 (c), for , has a peak () at , which approaches as is decreased, while at the same time . This prediction of the model is borne out by experimental SJ results, which show the same behavior with increasing . The weak dependence of on for is consistent with experiments Ren (2013), where this regime has very weak dependence on . In the limit of small , the DST regime of the model, the functional form of changes with , as shown in Fig. 2 (d)). As we discuss below, DST occurs in the range corresponding to , where is a monotonically decreasing function of , which explains the monotonic decrease of the stress anisotropy with observed in numerical simulations Seto et al. (2013); Mari et al. (2014).
iii.3 Scaling and Hysteresis in SJ Regime
The MF calculation shows that at small and , , and for any . Physically, this region corresponds to a system in which there are a large number of contacts even at zero driving. The peak anisotropy vanishes as (Fig. 3), identifying as the only characteristic disorder scale in this regime. The two order parameters, and are functions of both and . However, as shown in Fig. 4(c), upon definition of a dependent characteristic field: , they obey a scaling form: and , where, and are scaled variables: . The implication of this result is that in the regime, the behavior at different disorder strengths is controlled by the physics of the point (, ), reminiscent of critical phenomena Goldenfeld (1992). It was hypothesized by Bi et al Bi et al. (2011) that () is a critical point marking the end of a line separating fragile and SJ states. The critical point was characterized by the vanishing of an order parameter, which measures the anisotropy of the stress tensor. The current results, based on the spin model, are consistent with that picture. Numerically, meanfield predicts , and this exponent collapses experimental data for stress anisotropy and during a forward shear run Ren (2013), if we identify with (Fig. 4(d)).
The SJ experiments exhibit the phenomenon of Reynolds pressure Ren et al. (2013): pressure increasing quadratically with shear strain at small strains, with a Reynolds coefficient that depends only on and appears to diverge at (Table. 2). Very general arguments lead to the quadratic dependence of the pressure on shear strain Tighe (2014). If we make the logical assumption that the pressure increase is determined completely by , and that pressure increases as some monotonic function of , and hence , then our model provides a natural explanation for the observed dependence of pressure. The scaling form of implies that the pressure scales as: , where is a scaling function similar to defined above. The crucial feature of the scaling argument is the vanishing of as . From symmetry arguments, the pressure has to increase as some even function of the shear strain Tighe (2014). ( in the model), increases at least as fast as quadratically with for . Combined with the scaling form, this argument implies a divergence of the Reynolds coefficient as some power of , and therefore, as , where depends on the exponent , and the form of for small . From the perspective of the model, the source of the divergence observed in experiments is, therefore, directly related to the rapid rise in the number of contacts with shear strain as increases towards : a feature that is consistent with experimental observations.
In the mean-field approximation, there is no hysteresis for . As we show in Fig. 5(a), numerical simulations of the model in 2D exhibit hysteresis in this regime. It is to be noted that, in the simulation, the values of and , which define the SJ regime differ from the MF calculations. However, the overall structure of the phase diagram remains unchanged, as shown in Fig. 6. The model predictions for the scaling of the hysteresis loops (Fig. 5(b)) with are summarized in Table 2, and compared to dependence observed in SJ experiments Ren et al. (2013).
It is clear from the phase diagram, that the behavior of the model is completely smooth in the regime : all properties are continuous but sharp changes occur in the order parameters. This suggests that the SJ transition in dry grains with frictional coefficient Ren (2013) is not a phase transition but a crossover phenomenon at which the contact force network changes continuously both as a function of and . Preliminary analysis of experiments with lower friction coefficient between grains Wang and Behringer () suggests that with decreasing friction coefficient approaches from above, which leaves open the possibility of a true transition.
iii.4 Scaling and Hysteresis in the DST Regime
In the low regime, meanfield analysis predicts multiple solutions to and hysteresis under cyclic driving. We can identify three lines based on the multiplicity of solutions: For , (i) meanfield predicts two solutions for with accompanying hysteresis; (ii) for , multiple solutions appear for leading to multiple hysteresis loops, as shown in Fig. 8. As seen in Fig. 2(b), there is a critical point, , marking the end of these three transition lines. The and the lines are present in numerical simulations in 2D but the line is a meanfield feature. Simulations exhibit hysteresis over most of the region in Fig. 2(b), however, their characteristics change at , and . The line marks a discontinuous transition at which the peak anisotropy decreases dramatically, as shown in Figs. 3 and 8.
Identifying with the largest packing fraction, , at which one can have any flow Wyart and Cates (2013), and with the smallest packing fraction, , for the onset of DST, our results imply that two distinct types of force networks are stable in suspensions with : one with small stress anisotropy and large creating a highly connected network of force-bearing linear structures, and one with larger stress anisotropy and smaller . The networks with large also have large pressures since in our picture, the pressure is determined by .
Meanfield analysis shows that the hysteresis loops span a range , which grows as for , and at . These observations are in accord with scaling predictions of hysteresis loops in DST Wyart and Cates (2013) if we associate with .
Iv Discussion
We have constructed a driven, disordered, zero-temperature (non-equilibrium) statistical mechanics model, which captures all essential features of shear-induced rigidity transitions in granular materials and dense athermal suspensions. Our analysis highlights the distinct roles played by density and driving in the SJ and DST regimes: density controls the strength of disorder, whereas the driving field induces rigidity by increasing the concentration of frictional contacts. Based on analysis of experiments and simulations, we can assert that the observed phenomenology maps either to the () part (Table. 2) of the phase diagram based on whether at zero shear is small Seto et al. (2013); Mari et al. (2014) (large Ren (2013)). Controlling this parameter, for example Wang and Behringer () by tuning the friction coefficient of grains, provides an effective way of controlling where the system lies along the axis in our phase diagram, and probing the behavior near the critical point:.
Non-equilibrium critical points of random field models in the Ising class are characterized by avalanche distributions and crackling noise Sethna et al. (2001). Preliminary simulations in 2D indicate that the avalanche distribution exhibits a power law all along . Our model is distinguished from the Random Field Ising Model in a crucial way: can flip back to their original state, mediated by flips to even if the field is being increased monotonically, and the energy at a site does not approach the “flip” threshold monotonically. Recent studies Lin et al. (2014) show that this feature affects the yielding transition, suggesting that our model is relevant for understanding the yielding of athermal materials. We have focused on the shear-induced rigidity aspect of athermal, particulate systems. Yielding of the jammed states presumably occurs when the number of frictional contacts is saturated, and shearing does not lead to formation of new contacts. We are beginning to explore our model in this regime, where and independent of .
Acknowledgements.
We acknowledge extended discussions with Dapeng Bi, Eric Brown, R. P. Behringer, Jie Ren, Joshua Dijksman, and Dong Wang, and the hospitality of KITP, Santa Barbara, and are grateful to the Behringer group for sharing experimental data. This work has been supported in part by NSF-DMR 0905880 1409093, by the W. M. Keck foundation, and NSF-PHY11-25915*
Appendix A
In the following sections we discuss several properties and important aspects of the MF solution. Especially, we discuss the special disorders, which define different boundaries of the MF phase diagram. The special case of trajectories is also discussed here.
a.1 Zero disorder behavior
The energy of a spin in the zero disorder limit is:
(14) |
If , and if , it will flip to state and vice versa. A similar calculation can be done for . Hence, at zero disorder there is a discontinuous transition from to .
a.2 Asymptotic behavior of the model
The asymptotic, large field, behavior of the model is governed by the last two terms in Eqn. 5. Thus, the effective model governing the behavior at large positive , with can be written as:
(15) | |||||
(16) |
The first term on Eqn. 15 favors production of when , whereas the second term favors production of when . Since depends on , the asymptotic behavior of the model crucially depends on the functional dependence of on , which we refer to as a protocol. For a linear protocol as in Eqn. 16, which is the only kind we have analyzed, the asymptotic behavior depends on the slope, . If , dominates , . Conversely, if , dominates , and . If , there is no dependence and the asymptotic behavior depends on other terms in Eqn. 5. We discuss the trajectories in the following section.
a.3 Special disorders for trajectories
The meanfield equations for , and admit three lines of transitions which end at a critical point . These lines are defined by , , and in descending order of magnitude (Fig.1(a) in main text). The line , marks the transition from a single solution for for to multiple solutions over a range of (Fig. 8), whereas the line marks the transition from multiple solutions for with for to a single solution with for , as shown in Fig. 8. Notably, has a discontinuity at with the discontinuity increasing as , as seen in Fig. 3. The transition at is continuous. The transition lines, and , can be calculated analytically from the mean field equations and yields: , and . Here is the product log function, also known as Lambert’s function.
The transition at is a unique feature of our model and marks the onset of multiple solutions to , accompanied by system-size avalanches in which spins flip from to . is difficult to calculate analytically, and the line shown in Fig. 1(a) of the main text has been obtained numerically. Apart from very close to , . Fig. 8 illustrates the behavior of the system near these special disorders by comparing the -hysteresis. In the main text, these special disorders have been related to special packing fractions relevant to the DST transition.
a.4 , trajectories
For , increases and decreases as is increased, indicating that proliferate (Fig. 7 (a)). This trajectory is, therefore, not relevant to shear induced rigidity where grains with three or more contacts () proliferate, as the system is driven towards jamming.
Lying between and trajectories, trajectories exhibit an interesting dynamics (Fig. 7(b)). Since the chemical potential, changes at the same rate as , the applied field, the production of spins favored by competes equally with the production of spins favored by . For trajectories, the magnetization starts to decrease with increasing for , as depicted in Fig. 2 of the main text. In contrast, for , we observe that both the magnetization , and the fraction of zero spins asymptote to a disorder-dependent values and for . As increases, increases while decreases as shown in Fig. 9.
Simulations of the model (Eq. 5) in 2D, using zero temperature Monte Carlo dynamics, show that the asymptotic states for have a non-trivial spatial distribution of spins. As shown in Fig. 10, there is micro-phase separation between and spins. This spatial structure is reminiscent of shear bands observed in shear jamming experiments [3].
References
- M. E. Cates, J. P. Wittmer, J.-P. Bouchaud, and P. Claudin, Phys. Rev. Lett. 81, 1841 (1998).
- D. Bi, J. Zhang, B. Chakraborty, and R. P. Behringer, Nature 480, 355 (2011).
- J. Zhang, T. S. Majmudar, A. Tordesillas, and R. P. Behringer, Granular Matter 12, 159 (2010).
- J. Ren, J. A. Dijksman, and R. P. Behringer, Phys. Rev. Lett. 110, 018302 (2013).
- S. Sarkar, D. Bi, J. Zhang, R. P. Behringer, and B. Chakraborty, Phys. Rev. Lett. 111, 068301 (2013).
- E. Brown and H. M. Jaeger, ArXiv e-prints (2013), eprint 1307.0269.
- R. Seto, R. Mari, J. F. Morris, and M. M. Denn, Phys. Rev. Lett. 111, 218301 (2013).
- C. Heussinger, Phys. Rev. E 88, 050201 (2013).
- M. Wyart and M. Cates, arXiv:1311.4099 (2013).
- N. Fernandez, R. Mani, D. Rinaldi, D. Kadau, M. Mosquet, H. Lombois-Burger, J. Cayer-Barrioz, H. J. Herrmann, N. D. Spencer, and L. Isa, Phys. Rev. Lett. 111, 108301 (2013).
- A. Fall, N. Huang, F. Bertrand, G. Ovarlez, and D. Bonn, Phys Rev Lett 100, 018301 (2008).
- A. Fall, F. Bertrand, G. Ovarlez, and D. Bonn, Phys Rev Lett 103, 178301 (2009).
- A. Fall, A. Lemaître, F. Bertrand, D. Bonn, and G. Ovarlez, Phys Rev Lett 105, 268303 (2010).
- R. Mari, R. Seto, J. F. Morris, and M. M. Denn, Journal of Rheology 58, 1693 (2014), eprint 1403.6793.
- J. Ren, Ph.D. thesis, Duke University (2013).
- B. Utter and R. P. Behringer, Phys. Rev. Lett. 100, 208302 (2008).
- R. P. Behringer, D. Wang, J. Ren, and J. Dijksman, in Bulletin of the American Physical Society (2014).
- M. L. Falk and J. S. Langer, Phys. Rev. E 57, 7192 (1998).
- E. Vives and A. Planes, Phys. Rev. B 50, 3839 (1994).
- E. Vives, J. Goicoechea, J. Ortín, and A. Planes, Phys. Rev. E 52, R5 (1995).
- J. P. Sethna, K. A. Dahmen, and O. Perkovic, arXiv:cond-mat/0406320 (2004).
- R. Mari, R. Seto, and J. F. Morris, private communication.
- N. Goldenfeld, Lectures on Phase Transition (Westview Press, 1992).
- B. Tighe, Granular Matter to appear (2014).
- D. Wang and R. P. Behringer, private communication.
- J. P. Sethna, K. A. Dahmen, and C. R. Myers, Nature 410, 242 (2001).
- J. Lin, E. Lerner, A. Rosso, and M. Wyart, Proceedings of the National Academy of Sciences 111, 14382 (2014).