Shocks near Jamming

Shocks near Jamming

Leopoldo R. Gómez    Ari M. Turner    Martin van Hecke    Vincenzo Vitelli Instituut-Lorentz for Theoretical Physics, Leiden University, Leiden NL 2333 CA, The Netherlands.
Department of Physics, University of California, Berkeley, California 94720, USA.
Kamerlingh Onnes Lab, Universiteit Leiden, Postbus 9504, 2300 RA Leiden, The Netherlands.
July 20, 2019

Non-linear sound is an extreme phenomenon typically observed in solids after violent explosions. But granular media are different. Right when they jam, these fragile and disordered solids exhibit a vanishing rigidity and sound speed, so that even tiny mechanical perturbations form supersonic shocks. Here, we perform simulations in which two-dimensional jammed granular packings are dynamically compressed, and demonstrate that the elementary excitations are strongly non-linear shocks, rather than ordinary phonons. We capture the full dependence of the shock speed on pressure and impact intensity by a surprisingly simple analytical model.

Granular materials exhibit a wide range of complex collective behaviors, making them an important testing ground for the physics of amorphous materials Epitome ()-Xu2010 (). The confining pressure is perhaps the most important parameter controlling their properties. Strongly compressed granular media are, in many aspects, simple solids in which perturbations travel as ordinary phonons. However, when the confining pressure is lowered to zero, or the amplitude of the disturbance is much higher than the initial compression, the mechanical response of granular media becomes increasingly anomalous.

Several insights have been obtained by studying a simple model of granular media comprised of soft frictionless spheres just above the jamming point Epitome ()-Xu2010 (). The jamming point corresponds to the critical density at which the grains barely touch and vanishes Epitome (). The first insight is that the vibrational modes of jammed packings resemble ordinary phonons only below a characteristic frequency scale that vanishes as goes to zero Wyart ()-Souslov (). Above , the modes are extended but strongly scattered by disorder Xu ()-mattransport (). Second, as a direct consequence of the nonlinear dependence of the local contact force on the grain deformations, the sound speed vanishes as goes to zero Nagel ()-mattransport (): linear sound cannot propagate when the particles barely touch. Third, the range of validity of linear response vanishes when goes to zero. This is intuitive since the material is about to fall apart when the pressure vanishes Xu2010 ().

As the pressure (or density) is lowered towards the jamming point, there are thus three anomalies that forbid the propagation of ordinary phonons: disorder disrupt phononic transport for all frequency scales, the sound speed vanishes and linear response is no longer valid. The vanishing of the sound speed and absence of a linear range clearly suggest that the excitations near jamming will be strongly nonlinear. Nevertheless, most numerical and analytical studies of energy transport have been carried out in solids just above the jamming point, within a vanishingly small window of linear response. By explicit design, these studies cannot probe non-linear energy transport because the dynamics of the system is solved through a normal mode expansion Xu (); mattransport (); Vitelli (); Somfai (). Therefore, with the exception of theoretical and experimental studies on solitons in one dimensional granular chains, started with the seminal work of Nesterenko NesterenkoBook (); Nesterenko1984 (); Sen (); Daraio (); DaraioII (), non-linear energy transport in granular packings remains largely unexplored.

Numerical model.   To probe how elastic energy is transported close to the jamming point, we performed molecular dynamics simulations of a piston-compression experiment carried out in two dimensional polydisperse amorphous packings of soft frictionless spheres, whose radii, , are uniformly distributed between 0.8 and 1.2 times their average . Particles and at positions and interact via a non-linear repulsive contact potential Somfai ():


only for positive overlap , otherwise when . Here, the interaction parameter is expressed in terms of the effective Young’s modulus of the two particles, , see Ref. Somfai () for more details. The case corresponds to Hertz’s law. Lengths are measured in units of average particle diameter. The unit of mass is set by fixing the grain density to unity. The effective particle Young modulus is set to one, which becomes the pressure unit. These choices ensure that the speed of sound inside the grain, is one Somfai ().

We prepare Hertzian packings at a fixed pressure , or equivalently, an average particle overlap . They are then continuously compressed by a piston which moves with a constant velocity in the direction throughout the simulation, see Fig.  1. The subsequent motion of the particles is obtained by numerical integration of Newton’s equations of motion subject to periodic boundary conditions in the direction and a fixed boundary on the right edge of the system. We use two dimensional packings in the range of to particles with various width to length ratios.

Figure 1: Snapshots of the piston-compression simulation. A massive piston moves to the right at a constant velocity , resulting in the formation of a compression front traveling at a speed . Color indicates the local pressure at each grain. The average particle overlap is to the right of the front and to the left.

Phenomenology.   The piston compression leads to the formation of a front that separates two states. Ahead of the front, we find a region where the particles are at rest having the initial overlap , whereas behind it there is a compressed region with particles moving on average with the piston speed and an overlap . Figure  2a shows typical profiles for the longitudinal particle velocity (in the direction) as a function of , obtained upon averaging velocity fluctuations in the direction.

Two qualitative features of the shocks stand out for all the amorphous packings probed in this study: the fronts are smooth and stable. The smoothness can be contrasted with the typical shock profile that arise in ordered lattices of grains. Figure  2b, obtained for a triangular lattice of grains with zero initial overlap, shows large coherent pressure oscillations caused by the in-phase motion of the crystalline planes. These peaks are washed out by disorder in the amorphous packings.

Second, we have systematically tested the stability of the front against sinusoidal perturbations (in the -direction) of varying amplitudes and wavelengths in disordered packings under various pressures. This was done through directs simulations Gomez () as well as by performing a Dyakov’s stability analysis Gomez (); Dyakov (); Swan (). A typical result from our simulations, illustrated in Fig. 2c, shows how the front remains stable due to a classic stress focusing process, where particles “left-behind” experience a large compression, pushing them to catch up with the rest of the front. In the light of these observations, the shocks can be treated as one dimensional front propagation phenomena.

Figure 2: (a) Profiles of a shock wave at different times obtained by averaging the particle velocity in the direction (symbols). The red lines shows the fits of the fronts to an empirical fit formula. (b) Oscillatory velocity profile of a shock-like wave propagating through an hexagonal array. The inset shows the in-phase motion of the crystalline planes that gives rise to the pressure oscillations. (c) Representative snapshots of the focusing and flattening of an initially curved front generated by a sinusoidal piston (time progresses from left to right).

Front speed.   Once transients have died out, the front propagates with constant speed in the amorphous packings. Upon using conservation of mass across the shock front, we derive a one dimensional relation between the characteristic velocities and , through the average radius of the particles, , and the average compression in the shock , and ahead of it, :


Since the particle compression is typically much less than its diameter , Eq. (2) implies that . This is consistent with our numerical findings summarized in Fig. 3a where the dependence of on is explored systematically for different compressions.

Inspection of Fig. 3a reveals two distinct regimes. For low , the front speed is nearly independent of - in this (quasi)linear regime, is simply controlled by the initial pressure . The strongly non-linear shock waves regime is reached for high compression speed , where depends on , but not on .

The data for can be collapsed onto a single master curve, as shown in Fig. 3b. We achieve this upon rescaling the axis by , the numerically determined value that the front speed attains in the limit of vanishing (see Fig. 3a). The axis is rescaled by a pressure-dependent velocity scale , obtained by matching the low and high asymptotes in Fig. 3a (see arrow): marks the crossover between linear acoustic waves and shocks.

Scaling analysis.   The pressure dependence of can be rationalized using scaling arguments. We expect that reduces to , the speed of linear longitudinal sound waves. To determine the scaling of with pressure, note that , where the bulk modulus and . The change in volume scales linearly with , the average overlap between particles, while the energy scales as , see Eq. 1. Upon setting , we obtain the pressure dependence of the longitudinal speed of sound valid for Hertzian interactions Somfai (). Figure 3c shows that the numerical data for , represented by red symbols, is consistent with the scaling, which is shown as a continuous red line.

We now turn to the regime of high piston speeds, , when the front speed becomes nearly independent of . Since , and are all known, we need one additional relation which combined with Eq. (2) will make a definite prediction for the shock speed. We note that for strong shocks, the propagating front generates a characteristic compression and a corresponding increase in the kinetic energy. By assuming that the kinetic and potential energies are of the same order, we obtain . We have tested numerically that this non-trivial proportionality relation exists for strong deformations, see Fig. 3d. Upon combining the balance between kinetic and potential energy with Eq. (2), one readily obtains the power law , plotted as a dashed line in Fig. 3b. This scaling relation is clearly consistent with our numerical data for the speed of strongly non-linear shock waves.

Figure 3: a) Speed of the front versus particle velocity measured in units of , the sound speed within the grain, for decreasing particle overlap . b) Same plot as in (a) but with normalized by and normalized by the crossover particle speed : and are indicated in panel (a). The dashed line indicates the power law characteristic of a sonic vacuum. The black line indicates the theory developed here to describe the universal transition from weakly to strongly non-linear waves in systems close to jamming. c) Variation of and with distance to the jamming transition parameterized by the initial average overlap . The dashed lines indicate the power laws , . d) Variation of the kinetic energy with potential energy in dimensionless units - same color code as in (a-b). The dashed line indicates the linear relationship observed for strong shocks.

We deduce the dependence on compression of the crossover speed by smoothly matching the two asymptotic relations for the front speed and . This leads to the power law relation (continuous blue line in Fig. 3c) that is consistent with our numerical values (blue symbols). Note that the data collapse in Fig. 3b depends only on the scaling and it is not sensitive on the precise definition of the crossover speed. Upon using the conversion relation , the intuitive expectation that the crossover takes place when is confirmed.

We conclude that by controlling or P, which parameterize the distance to the jamming point (at and ), we can tune and the onset of the strongly non-linear response of the packings. Our key numerical findings on the shock velocity summarized in Fig. 3 can be grasped from scaling near the jamming point.

Analytical model.   In order to account for the dependence of on and the smoothness of the shock profiles, we construct the simplest possible 1D model that quantitatively accounts for the trends observed in Fig. 3 and sheds light on the role of disorder.

In the continuum limit, we obtain the following equation governing the dynamics of the system in terms of the strain field comment ():


To gain some intuition for the physics behind Eq. (3), note that by setting , one recovers a linear dispersive wave equation, with speed proportional to in the long wavelength limit. By contrast, when a non-linear wave equation is obtained. Nonlinearities and dispersive effects gives rise to finite amplitude waves: either solitary waves or shocks are possible depending on the drive NesterenkoBook ().

Shock propagation is modeled by the combined strain , where gives the shape of the shock and . Upon inserting this ansatz into Eq. (3), we obtain the conservation law , where is given by


This conservation law can be interpreted as describing the total energy of an effective particle at position rolling down a potential well , shown as a red line in Fig. 4a (here maps to time so that is the kinetic term of the particle) Herbold ().

One of the key ideas of our work is that disorder can act as an effective viscosity for the shock: the energy imparted unidirectionally by the piston is redistributed among other degrees of freedom, reducing the energy propagating with the shock front. In our mapping, this implies that the effective particle, initially located at the maximum of the potential , moves to the minimum of the potential well (see Fig. 4a). Thus, upon setting , we can obtain a relation between propagation velocity and induced compression in the front


that is independent of viscosity, even if an infinitesimal amount of dissipation is necessary to obtain a steady state solution of Eq. (3).

Figure 4: (a) Simulations of an ordered chain with small viscosity display large coherent oscillations in the front profile (black line) Herbold (). If the viscosity is large enough, one obtains an homogeneous shock profile, shown as a green line, similar to the profile in Fig. 2a. (b) The presence of an effective viscosity will induce the oscillation of the particle (black trajectory) towards the bottom of the potential , shown as a red line. If the viscosity is large enough, the particle can move directly to the minimum without performing any oscillations, see the green trajectory corresponding to the homogeneous shock profile of panel (a).

Together, Eqs. (2) and (5) can be seen as a parametric relation between front and particles velocities, where the overlap produced by the passage of the front is the parameter. Such a parametric plot of versus is drawn as a continuous curve on the numerical data in Fig. 3b. This comparison shows that Eqs. (2) and (5) are in excellent agreement (without any fitting parameter) with the results of our numerical experiments on shock propagation.

Discussion.   The shock formation explored in the present study is a generic phenomenon independent of the dimensionality of the sample that relies purely on the presence of a nonlinear law between grains (for any ) and not on the presence of friction. Experimentally, this can be tested by impacting a box of (frictional) glass beads with a heavy mass, for a range of impact speeds and pressures - preliminary experimental results for the front speed compare favorably to our theoretical predictions in Fig. 3 Martin ().

We note, however, that in frictional granular media, a second type of densification front can be observed, which is often referred to as plowing Jaeger (); Mandi (). Whereas our shock waves always propagate with speeds above the linear sound speed, and continue to propagate even after the driving stops, plowing fronts are generally much slower (in Jaeger (), of order 1 m/s), and stop almost immediately when the driving stops. We believe that the underlying difference is that our shocks are dynamical phenomena, set by a balance of potential and kinetic energies, whereas plowing is in essence a quasistatic phenomenon, dominated by dissipation. In the dynamic case, the change in packing fraction induced by the shock is associated with grain deformations, whereas in the quasistatic case, densification is dominated by grain rearrangements and compaction.

Outlook.   The shocks that arise in grains near jamming are just one representative of a broader class of strongly nonlinear excitations that emerge near the marginal state of suspensions, emulsions, wet foams and weakly connected fiber networks MvHecke (), mattPRL (), chase (). Close to losing their rigidity, all these materials exhibit a vanishing range of linear response, so that almost any amount of finite driving will elicit an extreme mechanical response in the form of rearrangements, yielding and flow Xu2010 (), brianPRL (), Gollub (), Zaccone (). It remains an open question whether all these phenomena can be successfully described in terms of simple scaling near jamming.

Acknowledgments We acknowledge inspiring discussion with X. Campman, V. Nesterenko, W. van Saarloos, A. Tichler and N. Upadhyaya. LRG and AMT acknowledge financial support from FOM and Shell.


  • (1) C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E 68, 011306 (2003).
  • (2) L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. Lett. 95 098301 (2005).
  • (3) M. Wyart, S. R. Nagel, and T. A. Witten, Europhys. Lett. 72, 486 (2005).
  • (4) W. G. Ellenbroek, E. Somfai, M. van Hecke, and W. van Saarloos, Phys. Rev. Lett. 97, 258001 (2006).
  • (5) A. Souslov, A. J. Liu, and T. C. Lubensky, Phys. Rev. Lett. 103, 205503 (2009).
  • (6) M. van Hecke, J. Phys.: Condens. Matter 22, 033101 (2010).
  • (7) C. H. Liu, and S. R. Nagel, Phys. Rev. Lett. 68, 2301 (1992).
  • (8) X. Jia, C. Caroli, B. Velicky, Phys. Rev. Lett. 82, 1863 (1999).
  • (9) J. R. Royer, B. Conyers, E. I. Corwin, P. J. Eng, and H. M. Jaeger, Europhys. Lett. 93, 28008 (2011).
  • (10) M. M. Bandi, T. Tallinen, L. Mahadevan, Europhys. Lett. 96, 36008 (2011).
  • (11) H. A. Makse, N. Gland, D. L. Johnson,, and L. Schwartz, Phys. Rev. E 70, 061302 (2004).
  • (12) E. Somfai, J. N. Roux, J. H. Snoeijer, M. van Hecke, and W. van Saarloos, Phys. Rev. E 72, 021301 (2005).
  • (13) N. Xu, V. Vitelli, M. Wyart, A. J. Liu, and S. R. Nagel, Phys. Rev. Lett. 102, 038001 (2009).
  • (14) V. Vitelli, N. Xu, M. Wyart, A. J. Liu, and S. R. Nagel, Phys. Rev. E 81, 021301 (2010).
  • (15) M. Wyart, M., Europhys. Lett. 89, 64001 (2010).
  • (16) N. Xu,, V. Vitelli, A. J. Liu, S. R. Nagel, Europhys. Lett. 90, 56001 (2010).
  • (17) L. Gomez et al., in preparation.
  • (18) V. F. Nesterenko, J. Appl. Mech. Tech. Phys. 5, 733-743 (1984).
  • (19) V. F. Nesterenko, Dynamics of Heterogeneous Materials (Springer-Verlag, New York, 2001).
  • (20) S. Sen, J. Hong, J. Bang, E. Avalos, and R. Doney, Physics Reports 462, 21-66 (2008).
  • (21) A. Spadoni and C. Daraio, Proceedings of the National Academy of Sciences 107, 7230 (2010).
  • (22) N. Boechler, G. Theocharis, and C. Daraio, Nat. Mat. 10, 665 (2011).
  • (23) S. P. Dyakov, Sh. Eksp. Teor. Phys. 27, 288 (1954).
  • (24) G. H. Swan and G. R. Fowles, Phys. Fluids 18, 28 (1975).
  • (25) Equation (3) is simpler than the one originally introduced by Nesterenko NesterenkoBook ().
  • (26) E. B. Herbold and V. F. Nesterenko, Phys. Rev E 75, 021304 (2007).
  • (27) S. van den Wildenberg, R. van Loo and M. van Hecke, in preparation.
  • (28) M. Wyart, H. Liang, A. Kabla, and L. Mahadevan, Phys. Rev. Lett. 101, 215501 (2008).
  • (29) C. P. Broedersz, X. Mao, T. C. Lubensky, and F. C. MacKintosh, arXiv:1011.6535 (2010). To be publised by Nature Physics.
  • (30) B. P. Tighe, E. Woldhuis, J. J. C. Remmers, W. van Saarloos, and M. van Hecke, Phys. Rev. Lett. 105, 088303 (2010).
  • (31) K. N. Nordstrom, E. Verneuil, P. E. Arratia, A. Basu, Z. Zhang, A. G. Yodh, J. P. Gollub, and D. J. Durian, Phys. Rev. Lett. 105, 175701 (2010).
  • (32) A. Zaccone and E. Scossa-Romano, Phys. Rev B 83, 184205 (2011).
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