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.
pacs:45.70.-n, 61.43.Fs, 65.60.+a, 83.80.Fg
The defining property of solid materials is their ability to resist external mechanical perturbations, as embodied in their finite elastic moduli. However, it is the geometry and topology of a materials’ 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 neighbours, both the linear elastic moduli and sound speed of the aggregate drop to zero. These extreme softness, that we term fragility, stems directly from the anharmonic 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. 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.
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. 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. 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 Nesterenko_Book (); Nesterenko_1984 (). The fate of these strongly non-linear excitations in amorphous two or three dimensional granular media near the (un)jamming transition or in weakly connected polymer networks is only gradually emerging Gomez_2012 ().
The very existence of these non-linear 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 enquire how fluctuations (that we expect to be violent given the absence of linear restoring forces) around these special points change their critical behavior. However, the study of fragile nano-mechanical systems affected by quantum fluctuations is still in its infancy. Even the impact of thermal fluctuations on the un-jamming transition is not very well studied, partly because granular media, that 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 amorphous or mass-disordered packings of grains that are just in contact with their nearest neighbors, see Fig. (1 a-b). For weak mass disorder, the solitary wave excitation generated in response to an impulse, decays exponentially at early times, with a rate that depends upon the amount of disorder. In the long time limit, the initially well defined solitary wave soon transitions into a triangular shock like profile, decaying as a power law, with a universal exponent consistent with and independent of the amount of disorder. This power law decay is the dominant mechanism of attenuation in hexagonal lattices with strong mass disorder as well as in jammed packings. 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. (1b) and the trail of the soliton in Fig. (1c). The resulting fluid-like state has two emergent macroscopic properties absent in the unperturbed granular packing: (i) a finite pressure that scales with the energy originally injected in the form of solitary wave and (ii) a finite viscosity, that arises despite no microscopic dissipative mechanism are present in our simulations.
Ii Attenuation of non-linear solitary waves
To study the dynamical response near the critical point, we perform molecular dynamics simulations for soft frictionless particles arranged in a hexagonal lattice or as a jammed amorphous packing. The particles interact pair-wise with the Hertz potential , where sets the energy scale and the particle overlap is measured in units of the average disk radii . By discretizing the packings into bins along the direction (defined as the longitudinal direction here), we impart an initial speed to the left most bin and obtain the subsequent dynamics by integrating the equations of motion using the velocity Verlet method Allen_Book (). We follow how the initial energy imparted to the packing evolves into a wave profile that propagates along the direction of initial impact (see schematic Fig. (1)) by averaging out the transverse speeds over the bins.
In Fig. ( 2), 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 variance (b) the same hexagonal packing with large variance in the masses of the grains and (c) a jammed amorphous packing prepared at a vanishingly small pressure P. 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 5 grain diameters, that propagates along the direction of impact. Remarkably, the solitary wave excitation is to a good approximation, the soliton-like solution first seen by Nesterenko in a one-dimensional granular chain of macroscopic non-cohesive particles interacting via the Hertz-potential Nesterenko_Book (); Coste2 (); Daraio2005 (); Job (); Sen (); Nesterenko_1994 ().
ii.1 Exponential decay
As the solitary wave propagates through a weakly disordered hexagonal packing (see inset of Fig. (2a), it begins to attenuate and the initial stages of this attenuation is 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 quasi-particle approximation to understand its decay in this regime Tichler_2013 (). Consequently, we model the interaction of the solitary wave quasi-particle 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 grains (undergoing multiple collisions) we find that the mean energy of the solitary wave is given by (see App. A)
where, is the initial energy of the solitary wave, is the solitary wave energy after traversing beads and is the variance in the masses of the beads, assumed to be distributed as a Gaussian random variable. As shown in the inset to Fig. (2 a), we find this estimate to be in very good agreement with our numerical observations on an hexagonal packing and weakly-disordered granular chains Manju_2012 ().
ii.2 Power law attenuation.
By contrast, for the hexagonal packing with strong mass disorder Fig. (2 b), the exponential regime is barely identifiable and the solitary wave begins to evolve into a shock like 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 Fig. (2 a-c). In all cases, a triangular shock like profile emerges, whose leading edge decays as a power law with an exponent approximately (solid red line). This exponent is thus independent of the amount of disorder.
Notice, the shock like 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 5 grain diameters, the shock solution spans several hundreds grains. In a recent study, the existence of such long wavelength triangular shock like solutions in a one dimensional chain of granular particles was established Calvo_2012 (). Here, we find that due to material inhomogeneities, an initial excitation in two dimensional disordered packings, naturally evolve into long wavelength shock solutions. As we outline below, the power law decay follows as a simple consequence of energy conservation.
In the long wavelength approximation, the average strain field satisfies the continuity equation , where subscripts denote partial derivatives with respect to time and space Calvo_2012 (). This equation has a similarity solution . Upon integrating once with respect to and differentiating once with respect to , we obtain the corresponding solution for the particle velocity field . Since the energy of the beads enclosed within the shock envelope is conserved, it follows that contant, where is the position of the shock front and consequently, . Thus, at the location of the front, the jump in the velocity field scales as , yielding the exponent . As shown in Fig. (2 a-c), solid red lines, we do find good agreement with the numerically determined exponent for the decay of the shock front in the power law regime.
Iii Emergent hydrodynamic state
The simple model for the exponential decay of the solitary wave suggests that at each subsequent collision with an inhomogenity, 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 in-elastically, (as a consequence of the Nesterenko equation of motion being non-integrable Nesterenko_Book ()) and thus their interaction may 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 two dimensional 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 equilibrium-like state that spans the entire finite-sized packings under-investigation.
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 in the -direction. Upon binning in the direction, we define the coarse grained particle current density and its Fourier transform . (The bold face is a vector notation used for spatial coordinates whose cartesian components are denoted by ). Upon setting and or , we determine numerically the corresponding longitudinal and transverse current density auto-correlation functions as , where the angular brackets denote ensemble averaging over the initial time. The longitudinal and transverse power spectral densities are evaluated numerically using fast Fourier transform from the respective current density auto-correlation functions as
In Fig. (3 b), we plot the longitudinal (red squares) and transverse (red circles) dispersion curves obtained numerically from the power spectral density Fig. (3 a) in 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. Since the total potential energy of a jammed packing is related to its pressure via E, we choose an initial impact speed to generate a solitary wave in the weakly compressed packings (P). This leads to a total initial solitary wave energy that is comparable to the energy of the highly compressed packing and thus allows for a more meaningful comparison between the two states.
As seen in Fig. (3 b), highly compressed jammed packings behave as ordinary solids with a finite bulk and shear modulus and this translates into a finite sound speed 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 as the second derivative of the induced potential energy; ; leads to the relation, Zhirov_2011 (), closely matching the numerical data in Inset to Fig. (3 b), red squares. Thus, the speed of the long wavelength hydrodynamical sound modes scales with the energy of the initial solitary wave injected into the system.
This is analogous to the scaling relation for pre-compressed jammed packings at a finite packing fraction , where the sound speed scales as , with being the potential energy due to the finite pre-compression . 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, one can easily replace the source of energy by a heat bath and thereby obtain a thermally induced rigidity upon making the substitution E. 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 quasi-equilibrium 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. (3 b), red circles. This is in stark contrast from a statically compressed jammed packing where a linear transverse dispersion regime (owing to the finite shear modulus) with a slope that scales with the amount of pre-compression (and does not depend upon the solitary wave energy injected) as is expected (Fig. (3 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. 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 solid-like medium with finite bulk and shear moduli, while in the latter, we fluidize the medium. This means that, even after the non-linear 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. This is one of the key conclusions of our work and it heralds the signal property of packings at the jamming threshold: they are fragile.
In Fig. (3 c), we show how the half width (inverse lifetime ) of the longitudinal hydrodynamical modes depend upon the wavenumber . We find that it scales anomalously as . For purely hydrodynamical modes obeying the Navier-Stokes equation, the half width scales with the wave number as hw , where is the shear viscosity. 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 Zwanzig (); Ernst (); Ernst_2 (). 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 dimensions, the shear viscosity itself scales with the wave number (in the small wave number limit) as and thus to leading order, we obtain a half width that scales as hw. These results have also been validated numerically in more recent studies on one dimensional non-linear spring chains Pikovsky (); Zhirov_2011 (). The predicted power law dependence is consistent with our findings in Fig. (3 c) solid red line. Our packings are two-dimensional but the emergent hydrodynamic description is effectively one-dimensional because we are injecting compressional solitary wave fronts in the 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.
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. We find that the physics of the solitary wave attenuation in jammed amorphous packings at their critical point corresponds to the strongly disordered regime in hexagonal packings where an initially well defined solitary wave soon transitions into a triangular shock like profile, decaying as a power law. We understand the two stages of decay using simple 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 fluid-like 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, despite no microscopic mechanism of dissipation being present.
We acknowledge discussions with J. M. J. Leeuwen and M. van Hecke. NU and LRG acknowledge financial support from FOM, Shell, Universidad Nacional del Sur and CONICET.
Appendix A The exponential decay
Consider a lattice of beads of mass , where nearest neighbours interact with a one sided non-linear potential of the form , with being the average overlap between neighbouring beads, the force constant and the exponent of the non-linear potential. If the initial overlap; , the effective spring constant defined as the second derivative of the potential vanishes, and as result, the lattice does not support linear sound at zero temperature. Consequently, mechanical strains generated in response to an impulse evolve into long-lived non-linear solitary waves and the solitary wave energy and momentum satisfy the classical relation , where is the mass of the solitary wave quasiparticle as a function of the mass of the beads Nesterenko_1984 (); Job (). Here, the momentum where, physically corresponds to the solitary wave amplitude.
The quasi-particle approximation of the solitary wave has been used to study the disintegration of a solitary wave across a mass interface Job (); Tichler_2013 (). Consider an interface between two regions of sonic vacuum with grain masses respectively. For mass ratios close to 1, a solitary wave initially moving with amplitude is seen to split into two new solitary waves, whose amplitudes may be obtained using an elastic collision model that conserves the quasi-particle energy and momentum -
from which one can easily derive
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 neighbouring beads related via , where , The normal random variable has mean 0 and variance . Upon appealing to the localized nature of the solitary wave (its width being around 5 bead diameters and also independent of the amplitude of the solitary wave), we treat each bead as an interface and invoke the quasi-particle elastic collision model.
Thus, the energy of the leading solitary wave after the first collision is
where, we have retained terms up to order . Iterating this process times, we find that the solitary wave energy after it has propagated n beads diameters is
Retaining terms only to order , we find
where, is the chi function with an expectation value . Taking the mean, we obtain that the average solitary wave energy after propagating beads diameters reads
This result agrees quantitatively with numerical simulations in 1D Manju_2012 () and qualitatively with our observations in 2D packings.