Size and density avalanche scaling near jamming
aa r X i v : . [ c ond - m a t . s o f t ] N ov Size and density avalanche scaling near jamming
Roberto Ar´evalo and Massimo Pica Ciamarra ∗ CNR–SPIN, Dipartimento di Scienze Fisiche, Universit`a di Napoli Federico II, I-80126, Napoli, Italy (Dated: October 16, 2018)The current microscopic picture of plasticity in amorphous assumes local failure events to producedisplacement fields complying with linear elasticity. Indeed, the flow properties of nonaffine systemssuch as foams, emulsions and granular materials close to jamming, that produce a fluctuatingdisplacement fields when failing, are still controversial. Here we show, via a thorough numericalinvestigation of jammed materials, that nonaffinity induces a critical scaling of the flow propertiesdictated by the distance to the jamming point. We rationalize this critical behavior introducing anew universal jamming exponent and hyperscaling relations, and use these results to describe thevolume fraction dependence of the friction coefficient.
PACS numbers: 60.20.F; 64.60.Av; 61.43.Er
While in simple crystals particles experience the sameamount of deformation under a small, uniform appliedstress, the disorder characterizing amorphous materialsleads to a fluctuating local deformation. When this non-affine contribution to the displacement field is of the sameorder of the affine displacement itself, it strongly influ-ence the elastic and rheological properties of the materialand cannot be handled via perturbative approaches [1–5]. In amorphous materials of technological interest suchas foams, emulsions, polymeric suspensions and granularmaterials, this occurs close to jamming volume fraction φ J marking the onset of mechanical rigidity upon com-pression [2]. Indeed, close to jamming nonaffinity effec-tively dominates the elastic response being responsible,for instance, of an anomalous scaling of the shear to bulkmodulus ratio, µ/k ∝ δφ / . Nonaffinity also influencesthe rheological properties at a finite shear rate, whereits effects are known to depend on the energy dissipationmechanisms [6, 7]. Much less is known about the role ofnonaffinity on the rheological properties in the athermalquasistatic shear (AQS) limit, which is of particular in-terest as it allows for the identification of the microscopicplastic events, and for the study of their properties andcorrelations. In this limit plastic events result from sad-dle node bifurcations [8] that drive the irreversible re-arrangement of a elementary unit of particles, generallyknown as a “shear transformation zone”, STZ [9]; thiselementary relaxation event might then trigger furtherrearrangements giving rise to avalanches, whose spatialfeatures, such as their fractal dimension, are controver-sial [8, 18, 20]. The triggering process is mediated by theelastic displacement field produced by STZs, experimentsand simulations (e.g. [8, 10, 11]) frequently found to bealike that resulting from an Eshelby inclusion [12] in a lin-ear elastic solid. This suggests that in strongly non affinesystems, where the displacement field produced by an el-ementary plastic event no longer resemble that producedby an Eshelby inclusion, the flow features might quali-tatively and quantitatively change. Recent AQS simu-lations [15, 16] of harmonic disks in two dimensions re- vealed that nonaffinity induces a critical behavior of somequantities characterizing the plastic flow, such as the av-erage stress and energy drops, or the length of the elasticbranches, which is dictated by the distance to the jam-ming point. However, it is not known whereas this behav-ior is universal, as the role of the interaction potential andthat of the dimensionality have not been explored, andthe critical exponents have not been theoretically ratio-nalized. Similarly, the effect of the increasing nonaffinityon the size scaling of the flow properties, that reveal thegeometrical features of the avalanches, is unknown.Here we investigate the role of nonaffinity in the flowproperties of amorphous systems via AQS shear simu-lations of harmonic and Hertzian particles, in both twoand three dimensions, as a function of the distance to thejamming threshold, δφ , and of the system size, N . Weshow that the macroscopic flow properties result qualita-tively unaffected by the degree of nonaffinity, as the dy-namics exhibits the same size regardless of the distanceto the jamming point. Nonaffinity controls the criticalscaling of the dynamics with δφ ; we rationalize this scal-ing showing that the exponent characterizing the criticalbehavior of the length of the elastic branches is univer-sal, and introducing hyperscaling relations involving thecritical exponents and the interaction potential. Finally,we use these results to infer the behavior of the frictioncoefficient at the jamming threshold.We consider 50:50 binary mixtures of N particles withdiameters σ = 1 and σ = 1 / . d = 2 and in d = 3 dimensions. Particles interact viaa finite range repulsive potential, V ( r ) = ε (cid:16) σ − rσ (cid:17) α for r < σ , V ( r ) = 0 for r > σ , corresponding to an harmonic( α = 2) or to an Hertzian potential ( α = 3 / σ average diameter of the interacting particles. ε and σ are our units of energy and length, respectively. We in-vestigate different values of the volume fraction, consid-ering systems with N = 400 , , , N = 6400 and 12800.These systems have been previously throughly investi-gated in the jamming context [2], and that the pressure,the bulk and the shear modulus are known to scale withthe overcompression as p ∝ δφ α − , k ∝ δφ α − , and as µ ∝ δφ α − / . Here we deform them via athermal qua-sistatic shear (AQS) simulations in which the strain isincreased by small steps δγ = 10 − , and the energy ofthe system is minimized via the conjugate–gradient pro-tocol after each strain increment [17]. Results are robustwith respect to a factor 10 change in δγ . In the steadystate regime of the AQS dynamics ( γ > σ Y ,the stress ∆ σ and the energy ∆ E drops, the strain inter-val between successive events, ∆ γ , and the instantaneousshear modulus, µ . For each value of N and δφ , we haverecorded from 10 to 10 plastic events.Our results support the presence of scaling relationsof the form h X i ∝ N ω X δφ ν X , where X is one of theinvestigated quantities, and ω X and ν X its size and den-sity scaling exponents, respectively. Data supporting thevalidity of these relations are illustrated in Fig. 1 andin Fig. 2a. For sake of clarity, we approximate the expo-nents with simple integer fractions to which the measuredvalues agree within ≈ .
02, and summarize their valuesin Table I.The size scaling is of interest as a signature of the mi-croscopic features of the plastic events and of their cor-relations [8, 18, 19], that build up through the elasticdisplacement field induced by the STZs. For instance,if plastic avalanches have a fractal dimension D , then ω ∆ E = D/d in d spatial dimensions [8, 18]. Previ-ous works have investigated the size scaling via AQSsimulations of a variety of different systems, confirm-ing the intensive character of the shear stress and ofthe shear modulus, ω σ Y = ω µ = 0, and revealing a sizedependence of other quantities associated to the plas-tic events. The energy–stress (∆ E ∝ N ∆ σ ) and thestress–strain (∆ σ ∝ µ ∆ γ ) dependence lead to the rela-tions ω ∆ E − ω ∆ σ = 1, and ω ∆ γ = ω ∆ σ . These are alwaysverified in the literature, most works reporting valuescompatible with ω ∆ E = 1 / ω ∆ E = 1 / ω ∆ E = 1 / ω ∆ σ = − / ω ∆ γ = − /
2. These val-ues do not depend neither on the interaction potential,nor on the dimensionality. These results clarify that themacroscopic plastic flow properties are unaffected by thenonaffine local response of the system.We now focus on the scaling with respect to the dis-tance from the jamming point, δφ . For each investigatedquantity and system size, we have performed a power lawfit to extract the critical exponent and the critical volumefraction, φ j . The resulting values of φ j are compatible -2 -1 -3 -2 -1 -8 -7 -6 -5 -3 -2 -1 -4 -3 -2 -1 -3 -2 -1 -6 -5 -4 -3 -2 -9 -8 -7 -6 (a) Hertz (b)
Harmonic E • N - / (c) Y (e) Y • N / (f) , Y (d) FIG. 1: Scaling of the shear modulus (panels a and b), of theenergy drops (panels c and d), of the yield stress and of thestress drops (panels e and f), for harmonic (left column) andHertzian particles (right column). Symbols refer to N = 1600(squares), 3200 (circles), 6400 (up triangles) and 12800 (downtriangles). A superimposed cross distinguishes the 3 d data.The data collapse occurs with no vertical shifts. Lines arepower–law with exponents summarized in Table I. within 10 − , consistently with the observation of a weaksystem size dependence under shear [23], and have av-erage values φ j = 0 .
843 in 2 d , and φ j = 0 .
645 in 3 d .These values are in good agreement with those reportedfor sheared systems [16, 24], and we have used them todetermine δφ = φ − φ j . Fig. 1 summarizes our results forthe scaling of the shear modulus, of the average energydrops, of the average stress drops and of the yield stress.The data collapse and the numerical fits clarify that thecritical exponents characterizing the δφ scaling of thesequantities depend on the interaction potential, not on thedimensionality. On the contrary, neither the interactionpotential nor the dimensionality affect the scaling of thelength of the elastic branches, we find to be characterizedby a universal scaling exponent, ν ∆ γ ≈ /
8, as illustratedin Fig. 2a.The scaling relations we have illustrated referring tothe average values do actually work for the whole dis-tributions, that scale as P ( X ) = h X i g X ( X/ h X i ). Here g X is a scaling function that depends on the considered -3 -2 -1 -1 -2 -1 -5 -3 -1 Harmonic • N / Hertz (a) P ( / <> ) /< > (b) FIG. 2: (a) Universal scaling of the average length of theelastic branches. Symbols are as in Fig. 1. Data for Hertzianparticles are shifted by a factor 2 for clarity. (b) Scaling of the∆ γ distribution. Data from 100 simulations obtained chang-ing N , δφ and interaction potential collapse on the same mas-ter curve. quantity, but not on the dimensionality, the interactionpotential or the volume fraction. Fig. 2b illustrates thisscaling for the length of the elastic branches reporting thecollapse of 100 datasets referring to different system sizes,densities and dimensionalities. Similarly, Fig. 3 showsthe collapse of the distributions of the energy and ofthe stress drops. These two distributions have a roughlypower law initial decay, which is followed by an expo-nential cut-off. In the case of the stress drop, the initialpower law exponent is ≈ −
1, as found [16] for d = 2 har-monic particles. For the energy drop we find an initialpower law exponent ≈ − .
7, in agreement with earlierresults [20, 25] for d = 2 harmonic particles. As an aside,we note that repulsive systems sheared via a spring mech-anism are common experimental and numerical modelsfor earthquakes [26], the particles representing the faultgouge in between the tectonic plates. These studies havefound a different power law exponent for the energy dropdistribution, ≈ − .
7, in agreement with the well knownGutenberg–Richter law. This suggests that in purely re-pulsive systems inertia affects the scaling exponents andthe scaling distributions, as recently observed in LJ sys-tems [27].We now show the existence of hyperscaling relationsbetween the critical exponents. We begin by notic-ing that there are no correlations between the mea-sured quantities (not shown), so that h ∆ σ i = h µ ∆ γ i = h µ ih ∆ γ i . Inserting the corresponding scaling relationswe recover the known relation ω ∆ σ = ω ∆ τ for the sizescaling, and find a relation for the density scaling, ν ∆ σ = ν µ + ν ∆ γ . (1)We then consider that, due to the quadratic relation be-tween energy and strain, the energy released in a plasticevent is ∆ E = ρ Nµ (cid:2) σ Y ∆ σ − ∆ σ (cid:3) . Given the scalingof the two term in brackets, the condition ∆ E > -1 -5 -4 -3 -2 -1 -4 -2 -5 -3 -1 (b) ( /< >) -1 P ( / <> ) /< > P ( E / < E > ) E/< E> ( E/< E>) -0.7 (a)
FIG. 3: Scaling of the distribution of the energy drops (a)and of the stress drops (b). As in Fig. 2b, data from 100simulations referring to different N , δφ and α collapse on thesame master curve. h X i ∝ N ω X δφ ν X Harmonic Hertz Relations ω ∆ E / / ω ∆ σ − / − / ω ∆ σ = ω ∆ E − ω ∆ γ − / − / ω ∆ γ = ω ∆ E − ν µ / ν µ = α − / ν ∆ σ / / ν ∆ σ = ν µ + ν ∆ γ ν σ Y / / ν σ Y = ν µ + ν ∆ γ ν ∆ E / / ν ∆ E = ν µ + 2 ν ∆ γ ν ∆ γ / / ω , and of the overcompression, ν . The exponents do not depend on the dimensionality. Ournumerical fits yield values consistent with the reported integerfractions, that satisfy the indicated relations. expressed as N β δφ ν σy − ν ∆ σ >
1, and is always satisfied if ν σ Y = ν ∆ σ . (2)Conversely, the condition could be violated at small orat large δφ . Given Eq. 2, the above equation for theenergy drop leads to the already known relation for thesize scaling, ω ∆ E − ω ∆ σ = 1, and to a new relation forthe density scaling, ν ∆ E = − ν µ + 2 ν ∆ σ = ν µ + 2 ν ∆ γ . (3)The validity of these relations, that can be also derivedfrom a simple dimensional analysis, is easily verified fromthe results summarized in Table I. As concern the depen-dence on the overcompression δφ , we have investigated 5exponents, and derived three hyperscaling relations. Theonly independent exponents are that of the shear modu-lus, which is known to be fixed by the interaction poten-tial, ν µ = α − / ν ∆ γ = 3 / p Y . First, we notice that a power -3 -2 -1 -1/8 -3 -2 -3 -2 (b) p Y 7/8 (a) -3 -2 -1 -3 -2 p y FIG. 4: (a) Crossover in pressure dependence on the over-compression δφ , for 3d harmonic spheres. The inset and themain panel show the whole dependence and a zoom on thecrossover region, respectively. (b) Scaling of the friction coef-ficient for all considered system sizes and potentials. Symbolsare as in Fig. 1. law fit of our pressure data in the investigated volumefraction range gives a critical exponent that is slightlylarger than that observed at zero applied shear stress, α −
1, consistently with previous results for 2d harmonicdisks [16, 24]. However, the fit gives a critical volumefraction that is smaller and not compatible with that re-sulting from the other power laws fit. We interpret this asa signature of the fact that the pressure does not behaveas a simple power law, being affected by both the com-pression and the shear stress. Indeed, when plotted as afunction of δφ , as in Fig. 4a for 3d harmonic spheres, thepressure exhibits two regimes. At high compressions, thepressure roughly scales as δφ α − , as in the case of zeroapplied shear stress, while at small compressions it scaleswith a smaller exponent, q . Simulations with a smallervalue of δφ , which are difficult to obtain due to the highcomputational cost of AQS simulations close to jamming,are needed to estimate q with confidence. Nevertheless,dimensional analysis suggests q = ν σ Y = α − /
8, whichis a value compatible with our data. The behavior ofthe pressure leads to a crossover in the effective frictioncoefficient η ∝ σ Y /p Y . This is expected to approach aconstant value close to jamming, as in previous numer-ical [28] and experimental works [29], as q = ν σ Y , andto scale as δφ − / away from the transition, regardlessof the dimensionality and of the interaction potential.Fig. 4b shows that these predictions are verified by ournumerical data.We have shown that the size scaling of the plasticflow features of amorphous materials is surprisingly un-affected by the distance to the jamming point, and there-fore by degree of nonaffinity. This suggests that, contraryto the common belief, Eshelby like displacement fieldsmight not play a fundamental role in the plastic flow, butmore work is needed in this direction. Nonaffinity quan-titaively influences the flow leading to a critical scaling of the dynamics with the distance to the jamming thresold,with exponents not depending on the dimensionality. Wehave rationalized this critial behavior introducing hyper-scaling relations between the exponents, and a new uni-versal jamming exponent. This universality suggests acritical behavior of the energy landscape itself.We thank R. Pastore for discussions, and MIUR-FIRBRBFR081IUK for financial support. ∗ [email protected][1] F. Leonforte, A. Tanguy, J. P. Wittmer, and J. L. Barrat,Phys. Rev. Lett. , 055501 (2006). A. Tanguy, F. Leon-forte, and J. L. Barrat, Eur. Phys. J. E , 355 (2006).[2] C.S. O’Hern, L. E. Silbert, A. J. Liu and S. R. Nagel,Phys. Rev. E , 011306 (2003).[3] M. van Hecke, J. Phys.: Condens. Matter , 033101(2010).[4] A.J. Liu and S.R. Nagel, Annual Reviews of CondensedMatter Physics , 14.1 (2010).[5] W.G. Ellenbroek, E. Somfai, M. van Hecke and W. vanSaarloos, Phys.Rev.Lett. , 258001 (2006). W.G. Ellen-broek, M. van Hecke and W. van Saarloos, Phys. Rev. E , 061307 (2009).[6] P. Olsson and S. Teitel, Phys. Rev. Lett. , 178001(2007).[7] B.P. Tighe, E. Woldhuis, J.J.C. Remmers, W. van Saar-loos and M. van Hecke Phys. Rev. Lett. , 088303(2010).[8] C.E. Maloney and A. Lemaitre, Phys. Rev. E , 016118(2006).[9] M. L. Falk and J. S. Langer, Phys. Rev. E , 7192(1998).[10] V.V. Bulatov and A. S. Argon, Modell. Simul. Mater.Sci. Eng. , 167 (1994); , 185 (1994); , 203 (1994).[11] V. Chikkadi, G. Wegdam, D. Bonn, B. Nienhuis and P.Schall, Phys. Rev. Lett. , 198303 (2011).[12] J. D. Eshelby, Proc. R. Soc. London Ser. A , 376(1957).[13] V. Chikkadi and P. Schall, Phys. Rev. E , 031402(2012).[14] R. Dasgupta, H.G.E. Hentschel and I. Procaccia,arXiv:1207.3591v3 (2012).[15] C. Heussinger and J.L. Barrat, Phys. Rev. Lett. ,218303 (2009).[16] C. Heussinger, P. Chaudhuri and J.L. Barrat, Soft Matter , 1 (1995).[18] N.P. Bailey, J. Schi¨otz, A. Lemaˆıtre and K.W. Jacobsen,Phys. Rev. Lett. , 095501 (2007).[19] A. Lemaˆıtre and C. Caroli, Phys. Rev. Lett. , 065501(2009).[20] C. E. Maloney and A. Lemaˆıtre, Phys. Rev. Lett. ,016001 (2004).[21] A. Tanguy, F. L´eonforte and J.L. Barrat, Eur. Phys. J.E , 355 (2006).[22] E. Lerner and I. Procaccia, Phys. Rev. E , 066109(2009).[23] D. V˚agberg, P. Olsson, and S. Teitel, Phys. Rev. E [24] D. V˚agberg, D. Valdez-Balderas, M.A. Moore, P. Olssonand S. Teitel, Phys. Rev. E , 4385 (1999).[26] F. Dalton and D. Corcoran, Phys. Rev. E , 061312(2001); M. Bretz, R. Zaretzki, S.B. Field, N. Mitarai andF. Nori, Europhys. Lett. , 1116 (2006); K. Daniels andN.W. Hayman, J. Geophys. Res. , B11411 (2008);M. Pica Ciamarra, E. Lippiello, C. Godano and L. deArcangelis, Phys. Rev. lett. , 238001 (2010); M. Pica Ciamarra, E. Lippiello, L. de Arcangelis and C. Godano,EPL , 54002 (2011).[27] K.M. Salerno, C. E. Maloney and M.O. Robbins, Phys.Rev. Lett. , 105703 (2012).[28] P.E. Peyneau and J.N. Roux, Phys. Rev. E , 011307(2008). T. Hatano, Phys. Rev. E , 060301(R) (2007).[29] R. Lespiat, S. Cohen-Addad and R. H¨ohler, Phys. Rev.Lett.106