Soliton attenuation and emergent hydrodynamics in fragile matter

Disordered packings of soft grains are fragile mechanical systems that loose rigidity upon lowering the external pressure towards zero. At zero pressure, we find that any infinitesimal strain-impulse propagates initially as a non-linear solitary wave progressively attenuated by disorder. We demonstrate that the particle fluctuations generated by the solitary-wave decay, can be viewed as a granular analogue of temperature. Their presence is manifested by two emergent macroscopic properties absent in the unperturbed granular packing: a finite pressure that scales with the injected energy (akin to a granular temperature) and an anomalous viscosity that arises even when the microscopic mechanisms of energy dissipation are negligible. Consistent with the interpretation of this state as a fluid-like thermalized state, the shear modulus remains zero. Further, we follow in detail the attenuation of the initial solitary wave identifying two distinct regimes : an initial exponential decay, followed by a longer power law decay and suggest simple models to explain these two regimes.


I. INTRODUCTION
The defining property of solid materials is their ability to resist external shear deformations, as embodied in their finite elastic shear modulus. However, it is the geometry and topology of a material's architecture that ultimately controls the strength of the elastic response. For example, if we remodel a solid block of steel into a granular aggregate of steel balls just in contact with their nearest neighbors, both the linear elastic moduli and sound speed of the aggregate drop to zero [1,2]. This extreme softness, which we term fragility, can have a local or global origin.
Locally, fragility stems from the strongly nonlinear interaction between macroscopic grains in contact. It is independent of the material the grain is made of, which could be as hard as steel or as soft as bubbles. On the other hand, weakly connected harmonic networks of cross-linked polymers or covalent glasses can also become fragile, irrespective of the individual bond strength, provided that the number of mechanical constraints is too low to maintain rigidity [3,4]. As a simple example, consider a square lattice of particles joined by Hookean springs that exhibits a vanishing shear modulus. In this case, fragility emerges as a property of a weakly connected network and it does not arise from nonlinear interparticle interactions. Granular packings exhibit both mechanisms close to the jamming threshold, where the grains are just touching: There are just enough forces between grains to hold the material rigid, and, at the same time, these forces are nonlinear. Moreover, the situation is further complicated by the presence of disorder that makes the response of a fragile material inherently nonaffine [5,6].
The critical point at which a packing of soft repulsive grains become fragile (unjams) can be accessed experimentally by progressively reducing the pressure, or alternatively, the average overlap, to zero [1]. Right at the critical point, even the tiniest mechanical strains must propagate in the form of supersonic solitons and shocks since the linear sound speed is zero [7,8]. Nesterenko coined the term sonic vacuum to designate this state in a granular chain and demonstrated the existence of solitonic excitations both experimentally and theoretically [2,7,9]. The fate of these strongly nonlinear excitations in amorphous two-or three-dimensional granular media near the (un)jamming transition or in weakly connected polymer networks is only gradually emerging [8][9][10][11][12][13][14].
The very existence of these nonlinear excitations and the associated critical points impinges on a very precise tuning of some geometrical parameters like the particle overlap that must vanish at random close packing or the coordination number of the network that must be tuned exactly to maintain rigidity. It is natural to inquire how fluctuations around these special points change their critical behavior. However, the study of fragile nanomechanical systems affected by quantum fluctuations is still in its infancy. Even the impact of thermal fluctuations on the unjamming transition is not very well studied, partly because granular media, which represent the key experimental arena for these investigations, are clearly athermal.
In this article, we study numerically and analytically the decay of solitary-wave excitations in two-dimensional mass-disordered or amorphous packings of grains that are just in contact with their nearest neighbors; see Fig. 1 and the movie in the Video 1. For hexagonal packings of grain with average mass m 0 and a small standard deviation Δm, an impulse excitation leads to the formation of a solitarylike wave whose amplitude decays exponentially at early times; see Figs. 2(a) and 2(b). We characterize this regime as the weakly disordered regime. Numerically, we find this regime to be approximately in the range of mass disorder 0 < Δm=m 0 < 0.45, where Δm is drawn from a normal random variable with mean m 0 and variance ϵ 2 that we vary to probe the different regimes. In the long-time limit, the initially well-defined solitary wave soon transitions into a triangular shocklike profile, decaying as a power law, with a universal exponent consistent with a value close to 1=2 and independent of the amount of disorder.
For hexagonal lattices with large variation in the grain masses (Δm=m 0 > 0.45), as well as in jammed amorphous packings, even at early times, a solitarylike wave is no longer identifiable and a triangular profile soon emerges, whose amplitude decays with the same exponent that is nearly 1=2. We characterize this regime as the strongly disordered regime, where the power-law attenuation is the dominant mechanism of decay.
Finally, we study the fluctuating state (similar to a thermalized state) that emerges after the complete disintegration of the initial solitary wave; see the inset of Fig. 2(b) and the trail of the soliton in Fig. 2(c). The resulting fluidlike state has two emergent macroscopic properties absent in the unperturbed granular packing: (i) a finite pressure that scales with the energy originally injected as a solitary wave and (ii) a finite viscosity that arises despite the fact that no microscopic dissipative mechanisms are present in our simulations.

II. ATTENUATION OF NONLINEAR SOLITARY WAVES
To study the dynamical response near the critical point, we perform molecular-dynamics simulations for N ¼ 4000 soft frictionless particles arranged in a hexagonal lattice or as a jammed amorphous packing. The particles interact pairwise with the Hertz potential VðδÞ ¼ Kð δ 2R Þ 5=2 , where K sets the energy scale and the particle overlap δ is measured in units of the average diameters of the disks a ¼ 2R. By discretizing the packings into bins along the x direction (defined as the longitudinal direction here), we impart an initial speed V 0 to the leftmost bin and obtain the subsequent dynamics by integrating the equations of motion using the velocity Verlet method [15]. We follow how the initial energy imparted to the packing evolves into a wave profile that propagates along the direction of initial impact (see the schematic in Fig. 2) by averaging out the transverse speeds over the bins.
In Fig. 3, we show the time evolution of the solitary-wave excitation generated in response to an impulse for three cases: (a) a hexagonal packing with grain masses, distributed as a normal random variable with a small standard deviation ϵ, (b) the same hexagonal packing with a large variation in the masses of the grains, and (c) a jammed amorphous packing prepared at a vanishingly small pressure P ∼ 10 −6 in units of K a 3 . The figures show the snapshots of the longitudinal (x direction) velocity field (solid black and gray curves) normalized by V 0 , taken at different instances of time, and reveals how the initial excitation evolves as it propagates through the medium along the x direction. In all three cases, we find that in response to an impulse, the mechanical excitation initially evolves into a well-defined solitary wave with a spatial extent around five grain diameters that propagates along the direction of impact. The solitary-wave excitation is, to a good approximation, the solitonlike solution first seen by Nesterenko in a one-dimensional granular chain of macroscopic noncohesive particles interacting via the Hertz potential [2,[16][17][18][19][20].

A. Exponential decay
As the solitary wave propagates through a weakly disordered hexagonal packing [see the schematic in Fig. 3(a)], it begins to attenuate, and the initial stages of this attenuation are well approximated as an exponential decay. The isotropic elasticity of the hexagonal packing allows us to model the solitary wave as a one-dimensional excitation-we make use of the quasiparticle approximation to understand its decay in this regime [21]. Consequently, we model the interaction of the solitary-wave quasiparticle with the mass inhomogeneity as an elastic collision process (conserving momentum and energy) that results in the disintegration of the solitary wave into two new solitary waves during each collision. After propagating through n grains (undergoing multiple collisions), we find that the mean energy of the solitary wave is given by (see Appendix A) where T 0 is the initial energy of the solitary wave, T n is the solitary-wave energy after traversing n beads, and ϵ is the The time evolution in response to an impulse in a hexagonal packing with weak mass disorder characterized by the standard deviation in mass distribution ϵ. An impulse response generates a solitary-wave excitation that decays exponentially at early times. In the long-time limit, a shocklike triangular profile emerges whose leading edge decays as a power law with an exponent x −0.5 (solid red line). (b) The time evolution in response to an impulse in a hexagonal packing with strong mass disorder ϵ. An impulse response generates a solitary-wave excitation, but the regime of exponential attenuation is barely identifiable. (c) Time evolution in response to an impulse in a jammed amorphous packing prepared at low pressure. The initial impulse still generates a solitary wave excitation but like in (b), a shock like triangular profile soon emerges whose leading edge decays as a power law with an exponent x −0.5 (solid red line).
standard deviation in the masses of the beads, assumed to be distributed as a Gaussian random variable. As shown in the inset to Fig. 3(a), we find the analytical estimate of Eq. (1), plotted as a continuous line, to be in very good agreement with our numerical data for hexagonal packings with ε ¼ f0.1; 0.25; 0.35g, plotted as black squares, red circles, and blue triangles, respectively. Our findings are consistent with previous studies conducted on weakly disordered granular chains [22].

B. Power-law attenuation
By contrast, for hexagonal packings with strong mass disorder [ Fig. 3(b)], the exponential regime is barely identifiable and the solitary wave begins to evolve into a shocklike triangular profile. For the decay in this regime, we find a striking similarity in all three cases: the long-time decay in a hexagonal packing with weak mass disorder, the dominant regime of decay in a hexagonal packing with strong mass disorder (large ϵ), as well as the dominant mechanism of decay in amorphous jammed packings [see Figs. 3(a)-3(c)]. In all cases, a triangular shocklike profile emerges, whose leading edge decays as a power law x −r with an exponent of approximately r ≈ 0.5 (solid red line), independent of the amount of disorder.
Notice that the shocklike profile is obtained as an envelope over several smaller excitations that span the system up to the leading propagating edge. Thus, unlike the initial solitary-wave excitation, whose spatial extent is around five grain diameters, the shock solution spans several hundreds of grains. In a recent study, the existence of such long-wavelength triangular shocklike solutions in a onedimensional chain of granular particles was established [23]. Here, we find that an initial excitation in two-dimensional disordered packings naturally evolves into long-wavelength shock solutions as a result of material inhomogeneities. As we outline below, the power-law decay follows approximately from the conservation of energy.
The broad spatial extent of the triangular profiles makes a long-wavelength description appropriate. In the longwavelength approximation, the average strain field δ satisfies the continuity equation δ t þ δ 1=4 δ x ¼ 0, where subscripts denote partial derivatives with respect to time t and space x [23]. This equation has a similarity solution δðx; tÞ ∼ ð x t Þ 4 . Upon integrating once with respect to x and differentiating once with respect to t, we obtain the corresponding solution for the particle velocity field ϕðx; tÞ ∼ ð x t Þ 5 . Since the energy of the beads enclosed within the shock envelope is conserved, Thus, at the location of the front, the jump in the velocity field scales as ϕðx f ; t f Þ ∼ As shown in Figs. 3(a)-3(c) (solid red lines), we do find good agreement with the numerically determined exponent r ≈ 0.5 for the decay of the shock front in the power-law regime.

C. Effects of dimensionality on wave evolution
The existence of a uniaxially propagating solitarylike wave in mass-disordered two-dimensional hexagonal packings is a novel feature (a consequence of its isotropic elasticity) and brings out an interesting similarity with solitary-wave propagation in one-dimensional granular chains with mass disorder [2,7,19,22], where our derivations for the exponential and power-law decays provide a very good match with simulations.
Like the one-dimensional case, in weakly disordered hexagonal packings, the amplitude of the solitary-wave excitation decays exponentially at early times, with a rate that depends upon the amount of disorder, while in the long-time limit, we observe a transition to a power-law regime, with a universal exponent of 1=2 that is independent of the amount of disorder. However, when we define a uniaxially propagating front in two dimensions, we inevitably ignore a part of the energy that goes into exciting transverse velocity fluctuations that effectively act as an additional source of dissipation. Nevertheless, our simulations reveal that these random excitations are small as long as a well-defined solitarylike wave and triangular shocklike profile can be identified against the background noise, and thus our results from one dimension still provide a reasonably good analytical approximation to the process of attenuation in two-dimensional amorphous packings.

D. An analogy between disorder and viscoelastic dissipation
In Fig. 3(a), we observe the formation of a two-wave structure in response to an impulse imparted to a hexagonal packing with weak mass disorder. Thus, here, an initial solitarylike wave is followed by a broader shocklike secondary pulse. As Nesterenko's leading solitarylike wave attenuates, the secondary structure evolves. Eventually, the leading wave completely disappears, leading to a new single-wave structure (which is still nonlinear) with a triangularlike profile.
This behavior has a striking similarity with the evolution of an impulse in one-dimensional granular chains with viscous dissipation, where similar two-wave structures are observed in simulations and experiments [24,25]. Although the disordered granular packings in our simulations do not include any source of real dissipation, these findings further reenforce the role of disorder as an effective viscosity that we capture analytically in the exponential and power-law attenuation regimes.
To explicitly confirm this analogy in two dimensions, we have compared the propagation of nonlinear waves on hexagonal granular packings with varying mass disorder and coefficients of viscous dissipation. Here, we model dissipation by penalizing relative velocities between particles with a force of the form In Fig. 4, we compare the evolution of an initial impulse imparted to a mass-disordered hexagonal lattice of grains in a sonic vacuum (left panel) with the same simulation performed on an ordered hexagonal lattice with viscous dissipation between contacts (right panel). The figure shows a snapshot of the velocity profiles at a given time, for different amounts of disorder and dissipation, respectively. Note how the effects of disorder and dissipation have striking similarities. Although energy is not dissipated behind the two-wave structure in the mass-disordered packings, it goes into fast uncorrelated oscillations, whose average does not produce any significant contribution to the velocity profiles.
As we increase the amount of disorder or dissipation in the system (indicated by colors in Fig. 4), we see the progressive shift from a solitarylike wave to a triangular shocklike structure. Above a "critical" value of disorder or dissipation, the two pulses merge into a single shocklikewave structure. (This characteristic value is not completely inherent to the structure and also depends on the energy of the initial impulse.) Note however, despite these similarities, the effects of disorder and dissipation are not exactly equivalent. In Fig. 4 right panel, we see that an explicit mechanism for viscous dissipation dampens out any relative motion between particles induced by the passage of the wave front and the resulting secondary wave structures are smooth. By contrast, disorder on its own merely traps the residual motion between particles in the form of uncorrelated oscillations as seen in the secondary wave structures in Fig. 4 left panel.

III. EMERGENT HYDRODYNAMIC STATE
The simple model for the exponential decay of the solitary wave suggests that at each subsequent collision with an inhomogeneity, the solitary wave splits approximately into two solitary waves-a leading pulse (the main degree of freedom that is damped exponentially as a result) and a smaller solitary wave that is either reflected or transmitted, depending upon the mass ratio. Therefore, once we allow sufficient time for the leading solitary wave to disintegrate completely such that a leading pulse is no longer distinguishable from the background, we expect to reach a state comprised of several smaller solitary waves with different energies. The solitary waves, in turn, interact with each other inelastically [19,27]. At the level of a continuum approximation, the inelastic interaction is a consequence of the Nesterenko equation of motion being nonintegrable [2], with the linear momentum and energy being the only two conserved quantities (see Appendix A). As a result, two solitary waves do not emerge unaffected after a collision event and lose a fraction of their initial energies into forming smaller solitary waves. Physically, the Nesterenko solitary waves are around five particle diameters wide, and numerical observations suggest that the solitary waves that emerge from a collision are reformed during the collision event and that the phase mismatch between the incoming solitary waves leads to newer excitations being generated [27]. The collision process may therefore be thought of as inherently dissipative. However, in addition to these processes (that seek to distribute the energy initially concentrated in a single solitary-wave excitation into multiple smaller solitary waves), the structural disorder in twodimensional amorphous packings also spills a part of the energy into transverse degrees of motion. Through a series of such intrinsic dissipative mechanisms, we eventually reach a fluctuating equilibriumlike state that spans the entire finitesized packing under investigation. In this equilibriumlike state, we find that the ratio of kinetic T to potential V energy satisfies the Virial relation, i.e., hTi ¼ 5 4 hVi, consistent with previous studies on mass-disordered linear chains [7].

A. Hydrodynamic modes
In order to rationalize the physics behind this fluctuating state that at first appears to be just noise, we assume that there is a well-defined long-wavelength, small-frequency collective hydrodynamical mode with wave number k in the x direction. Upon binning in the x direction, we define the coarse-grained particle-current density jðr; tÞ ¼ 1 ffiffiffi . For low disorder and dissipation, the pulses are weakly attenuated, and since the speed of the pulse depends upon its amplitude, they travel further into the packing. Comparing the two panels, we see the somewhat analogous effects of disorder and dissipation.
vector notation used for spatial coordinates whose Cartesian components are denoted by α.) Upon setting k ¼ ðk; 0Þ and α ¼ x or y, we determine numerically the corresponding longitudinal and transverse current-density autocorrelation functions as C l;t ðk; tÞ ¼ hj Ã l;t ðk; 0Þj l;t ðk; tÞi, where the angular brackets denote ensemble averaging over the initial time. The longitudinal and transverse power-spectral densities P l;t ðk; ωÞ are evaluated numerically using a fast Fourier transform from the respective current-density autocorrelation functions as Within the linear hydrodynamical approximation, the longitudinal power spectral density assumes the form where, η is the coefficient of shear viscosity and ω 0 ¼ ck is the linear dispersion relation indicating that the longitudinal modes propagate at the speed of linear sound c. By contrast, the transverse power spectral density has the form The absence of a linear dispersive term indicates that the transverse current density correlations simply decay exponentially in time at a rate equal to ηk 2 . Within the linear hydrodynamical approximation, the transverse modes are thus purely diffusive. In Fig. 5(a), we plot the longitudinal power-spectral density obtained numerically from molecular-dynamics simulations. By projecting the data on the two-dimensional ω-k plane, we obtain the dispersion curves, as anticipated by Eq. (3). A similar procedure is repeated to obtain the transverse dispersion curves.
In Fig. 5(b), we thus plot the longitudinal (red squares) and transverse (red circles) dispersion curves that we use to characterize the emergent fluctuating state. For comparison, we plot the corresponding longitudinal (black squares) and transverse (black circles) dispersion curves for a highly compressed jammed packing (far from the critical density), prepared at a pressure of P ∼ 10 −1 . For precompressed jammed packings with an average initial overlap between disks δ and an interaction potential of the form δ α , the pressure is related to the overlap via the relation P ∼ δ α−1 , while the shear modulus is expected to scale as G ∼ δ α−ð3=2Þ and the bulk modulus as B ∼ δ α−2 [1]. For Hertzian interaction with α ¼ 5 2 , the shear and bulk moduli therefore scale as G ∼ δ 1 and B ∼ δ 1=2 respectively. In addition, the energy of the packing due to precompression is E ∼ δ α , and therefore, G ∼ E 2=5 and B ∼ E 1=5 [1].
From the above scaling relations, the total potential energy of a jammed packing is related to its pressure via E ∼ P 5=3 . Therefore, by choosing an initial impact speed of V 0 ≈ 2.0 (in units of p K=m) for the weakly compressed packing (P ∼ 10 −6 ), we ensure that the energy imparted in the form of a propagating wave and consequently the energy of its final state E SW , is comparable to the static energy (E PC ) of highly compressed packings that are prepared at a pressure of approximately P ∼ 10 −1 .  5. (a) The power-spectral density for longitudinal modes that emerge in marginally compressed amorphous packings (P ∼ 10 −6 ) by the complete disintegration of an initial solitarywave excitation. (b) The red squares show the numerical data for the longitudinal dispersion curves. The slope of the linear regime scales with the energy of the solitary wave as c ∼ E 1=10 for a Hertzian interaction that we identify with long-wavelength longitudinal hydrodynamical modes. In contrast, the red circles show the numerically obtained data for the transverse dispersion curve, where no linear regime is seen. The shear mode is therefore nonpropagating. For comparison, shown are linear dispersion curves for longitudinal (black squares) and transverse (black circles) dispersion curves obtained for highly compressed jammed packings, prepared at a pressure P ∼ 10 −1 . In the inset, the speed of sound c determined numerically for different E is plotted as black squares, while the continuous red line is the scaling prediction E ∼ E 1=10 . (c) The half-width obtained numerically from the longitudinal modes as a function of wave number on a linear scale and compared against the analytical estimates HW ∼ k 5=3 (solid red line).
As seen in Fig. 5(b), highly compressed jammed packings behave as ordinary solids with a finite bulk and shear modulus, and this behavior translates into a finite soundspeed manifest in the linear regime of the longitudinal (black squares) and shear (black circles) dispersion curves. In contrast, exciting a jammed packing prepared at a very small pressures (that a priori has zero longitudinal and shear sound speeds) leads to a linear dispersion regime for the longitudinal modes, but no such regime is obtained for the transverse modes. The slope of the linear regime in the longitudinal dispersion curves corresponds to the speed of long-wavelength hydrodynamical sound modes.
Defining the sound speed c as the square root of the second derivative of the induced potential energy SW [28]. In order to see this scaling, we obtain the slope c for a range of initial impact energies E SW and show the numerical data (black squares) plotted against the estimate c ∼ E 1=10 SW (solid red curve) in the inset to Fig. 5(b), where we find good agreement between the two. Thus, the speed of the long-wavelength hydrodynamical sound modes scales with the energy of the initial solitary wave E SW injected into the system. This behavior is analogous to the scaling relation for precompressed jammed packings at a finite packing fraction δ, where the sound speed scales as c ∼ E 1=10 PC , with E PC being the potential energy due to the finite precompression δ. Thus, in so far as longitudinal sound modes are concerned, a rigidity induced by statically compressing a marginally compressed packing is analogous to the rigidity induced by exciting a marginally compressed packing with a finite energy wave. Note, therefore, that one can easily replace the source of energy by a heat bath and thereby obtain a thermally induced rigidity upon making the substitution E → k B T. However, unlike a state that is truly in thermal equilibrium, an external perturbation over the fluctuating state created by the disintegration of a solitary wave will further raise its energy due to the absence of a fluctuation-dissipation mechanism. The emergent state is thus at best described as a quasiequilibrium state.
In contrast to the longitudinal modes, the transverse modes obtained by energizing a marginally compressed packing do not show a well-defined linear regime; see Fig. 5(b) (red circles). This behavior in stark contrast to a statically compressed jammed packing where there is a linear transverse dispersion regime (owing to the finite shear modulus), with a slope that scales with the amount of precompression as c ∼ E 1=5 PC ; see Fig. 5(b) (black circles). Thus, the shear modes excited by injecting energy into a packing near its critical point are purely diffusive and the medium does not develop a finite shear modulus. Note that, prior to the impulse, the marginally compressed jammed packings used in our simulations have a vanishingly small but finite shear modulus of the order of G ∼ δ 0 ≈ 10 −4 , but the state of the packing is in mechanical equilibrium. However, the passage of the impulse disturbs this state and leads to rearrangement of particles via contact breaking, a mechanism analogous to melting of the system caused by thermalization [29]. There is therefore a profound difference between the resultant states obtained by (a) statically compressing a granular packing near its critical point versus (b) injecting energy either in the form of an excitation or thermally. In the former case, one obtains a solidlike medium with finite bulk and shear moduli, while in the latter, we fluidize the medium. The absence of linear response to shear means that even after the nonlinear wave excitations (characteristic of the jamming point) are strongly attenuated, one cannot describe the system purely in terms of the normal modes of a (linear elastic) solid because the system has de facto become a fluid. The emergence of a fluidized state is one of the key conclusions of our work, and it heralds the signal property of packings at the jamming threshold: They are fragile.

B. Emergent viscosity
In Fig. 5(c), we show how the half-width (inverse lifetime τ) of the longitudinal hydrodynamical modes depends upon the wave number k. We find that it scales anomalously as τ −1 ∼ k 1.6 . For purely hydrodynamical modes obeying the Navier-Stokes equation, the half-width is a direct measure of the viscosity and scales with the wave number as HW ∼ ηk 2 , where η is the shear viscosity; see Eqs.
(3) and (4). (The half-width evaluated from the power spectral density provides an alternate but equivalent definition of viscosity as the Green-Kubo integrals [30].) However, extensive numerical and analytical studies have shown that in one dimension, the time-correlation functions (whose long-time integrals by definition correspond to macroscopic transport properties such as diffusivity and shear viscosity) do not decay exponentially but display long-time tails, decaying as power laws instead [31][32][33]. This phenomenon indicates the breakdown in low dimensions of the standard continuum assumption embedded in the Navier-Stokes equation-the strength of fluctuations is too strong for simple coarse-grained theories to hold.
Analytic studies based on mode-coupling theory predict that for d ¼ 1 dimensions, the shear viscosity itself scales with the wave number (in the small wave-number limit) as ηðkÞ ∼ k −1=3 , and thus to leading order, we obtain a halfwidth that scales as HW ∼ ηðkÞk 2 ¼ k 5=3 . The estimated 5=3 power-law dependence is consistent with our findings in Fig. 5(c) (solid red line). These results have also been validated numerically in more recent studies on onedimensional nonlinear spring chains [28,34], and in Appendix B, we discuss in more detail how this exponent is consistent with the Kardar-Parisi-Zhang (KPZ) universality class.
Our packings are two dimensional, but the emergent hydrodynamic description is effectively one dimensional because we are injecting compressional solitary-wave fronts in the x directions only; see Fig. 1. Thus, despite the absence of any microscopic dissipative mechanism in our simulations, there is an emergent shear viscosity that scales with the wave number and gives the longitudinal and shear hydrodynamical modes a finite lifetime.

IV. CONCLUSION
To conclude, we find that in response to an impulse, hexagonal packings with mass disorder and amorphous packings initially generate a well-defined solitary-wave excitation that is progressively attenuated. For weakly disordered hexagonal packings, the amplitude of the solitary-wave excitation decays exponentially at early times, with a rate that depends upon the amount of disorder. In the long-time limit, we observe a transition to a power-law regime, with a universal exponent that is independent of the amount of disorder and is the dominant mechanism of decay in the strongly disordered case, where an exponential decay is no longer identifiable.
The physics of the solitary-wave attenuation in jammed amorphous packings at their critical point corresponds to the strongly mass-disordered regime in hexagonal packings (and strongly disordered granular chains), where an initially impulse excitation soon transitions into a triangular shocklike profile, whose amplitude decays with the same universal exponent close to 1=2. Here, even at early times, a well-defined solitary wave can no longer be identified. However, notwithstanding the additional sources of dissipation in two-dimensional packings, we find that a good qualitative understanding of the two stages of decay can be obtained from simple one-dimensional analytical models and scaling arguments.
We then study the resultant fluctuating state (similar to a thermalized state) that emerges after the complete disintegration of the initial solitary wave. This fluidlike state has emergent macroscopic properties absent in the unperturbed granular packing: (i) a finite pressure that scales with the energy originally injected in the form of a solitary wave and (ii) a finite viscosity that is consistent with the equilibrium description of a onedimensional fluid.

ACKNOWLEDGMENTS
We acknowledge discussions with J. M. J. Leeuwen and M. van Hecke and thank the referees for the helpful feedback. N. U. and L. R. G. acknowledge financial support from FOM, Shell, Universidad Nacional del Sur, ANPCyT, and CONICET.

APPENDIX A: THE EXPONENTIAL DECAY
Consider a lattice of beads of mass m, where nearest neighbors interact with a one-sided nonlinear potential of the form VðδÞ ¼ K α δ α , with δ being the average overlap between neighboring beads, k the force constant, and α > 2 the exponent of the nonlinear potential. If the initial overlap δ → 0, the effective spring constant defined as the second derivative of the potential vanishes, and as a result, the lattice does not support linear sound at zero temperature. Consequently, mechanical strains generated in response to an impulse evolve into long-lived nonlinear solitary waves whose continuum approximation is well described by the Nesterenko equation of motion where c is a material constant, R is the radius of the beads, and ξðx; tÞ ¼ −∂ x δ is the strain field. The left-hand side of Eq. (A1) is the standard inertia term; the second term on the right-hand side captures nonlinear dispersive effects, while the first arises from the restoring force as in the wave equation, if one considers that the force is not linear but depends on ξ 3=2 , according to Hertz law. Equation (A1) possesses only two conserved quantities, the linear momentum and energy, and is therefore not known to be integrable [35]. However, to a good approximation, the solitary-wave energy E and momentum P satisfy the classical relation E ¼ P 2 2m eff , suggesting a description of the solitary wave as a quasiparticle with a mass m eff ∼ 1.3m, where m is the mass of the beads comprising the granular chain [7,18]. Here, the momentum P ¼ m eff U, where U physically corresponds to the solitary-wave amplitude.
The above quasiparticle approximation of the solitary wave has been used to study the disintegration of a solitary wave across a mass interface [18,21]. Consider an interface between two regions of sonic vacuum with grain masses m 1 and m 2 , respectively. For mass ratios A ¼ m 2 m 1 close to 1, a solitary wave initially moving with amplitude P 0 is seen to split into two new solitary waves, with momentum P 1 , P 2 that may be obtained using an elastic collision model that conserves the quasiparticle energy and momentum from which we obtain N. UPADHYAYA, L. R. GÓMEZ, AND V. VITELLI PHYS. REV. X 4, 011045 (2014) 011045-8 Thus, the ratio of transmitted to incident energy is In order to study the propagation of the solitary wave in a medium with weak mass inhomogeneity, consider a chain of beads with the mass ratio of neighboring beads i; j related via m i ¼ A i;j m j , where A i;j ¼ 1 þ N i;j ð0; ε 2 Þ. The normal random variable N i;j ð0; ε 2 Þ has mean 0 and variance ϵ 2 . Upon appealing to the localized nature of the solitary wave (its width being around five bead diameters and also independent of the amplitude of the solitary wave), we treat each bead as an interface and invoke the quasiparticle elastic collision model.
Thus, the energy of the leading solitary wave after the first collision is where we have retained terms up to order ϵ 2 . Iterating this process n times, we find that the solitary-wave energy after it has propagated n bead diameters is approximately T n ≈ 1 − 1 4 ϵ 2 N 2 n ð0;1Þ ÁÁÁ 1 − 1 4 ϵ 2 N 2 1 ð0;1Þ T 0 : (A5) Retaining terms only to order ϵ 2 , we find where χ 2 ðnÞ is the chi function with an expectation value n.
Taking the mean, we obtain that the average solitary-wave energy after propagating n bead diameters reads This result agrees quantitatively with numerical simulations in 1 D [22] and qualitatively with our observations in 2 D packings.

APPENDIX B: THE KPZ UNIVERSALITY CLASS
In this section, we discuss how the anomalous dependence of the half-width on the wave number, i.e., HW ∼ ηðkÞk 2 ¼ k 5=3 , can be related to the KPZ universality class. As discussed in the main text, the anomalous scaling (exponent 5=3 instead of 2) indicates the breakdown in low dimensions of the standard mean-field approximation embedded in the Navier-Stokes equationthe strength of fluctuations is too strong for simple coarsegrained theories to hold.
The Navier-Stokes equation is conventionally expressed without any additional sources of fluctuations. However, in analogy with the Langevin equation, we may think of perturbing the Navier-Stokes equation with a random source of noise that is introduced phenomenologically. Previous studies on such a model equation indicate that for three-and two-dimensional fluids near equilibrium, the noise only produces a small correction to the dynamics of large-scale and long-time properties (k → 0, ω → 0). However, the same is not true in one dimension [36,37]. In addition, as we review below, it is no longer justified to ignore the nonlinear term that we have ignored while linearizing the Navier-Stokes equations of motion in order to obtain Eqs. (3) and (4).
To illustrate this point, we consider the one-dimensional Burger's equation stirred with a random force as a model for a one-dimensional fluid: where v is the fluid velocity, ν is the kinematic viscosity, and λ ∂η ∂x is the random stirring force such that in d dimensions (d ¼ 1 here), hηðx; tÞi ¼ 0; (B2) hηðx; tÞηðx 0 ; t 0 Þi ¼ 2Dδ d ðx − x 0 Þδðt − t 0 Þ: The Burger's equation (without random stirring) is obtained from the Navier-Stokes equation by ignoring the compressibility (pressure-gradient term). As discussed earlier, the finite compressibility leads to propagating modes in the dispersion relation and the width of the power-spectral density is measured in a frame that is moving at the speed of sound. Thus, the propagating modes only shift the peak and do not affect the scaling properties of the width [37]. Upon rewriting Eq. (B1)) in terms of the potential function v ¼ −λ ∂ϕ ∂x , we arrive at the KPZ equation [38]: Now, consider rescaling the variables: x → bx 0 , t → b z t 0 , and ϕ → b ξ ϕ 0 . In terms of the primed variables, Eq. (B4)) becomes [38]