# Shocks near Jamming

###### Abstract

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 ():

(1) |

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.

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.

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, :

(2) |

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.

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 ():

(3) |

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

(4) | |||||

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

(5) |

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

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.

## References

- (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).