Discontinuous change of shear modulus for frictional jammed granular materials
aa r X i v : . [ c ond - m a t . s o f t ] D ec Discontinuous change of shear modulus for frictional jammed granular materials
Michio Otsuki ∗ and Hisao Hayakawa Department of Physics and Materials Science, Shimane University,1060 Nishikawatsu-cho, Matsue 690-8504, Japan Yukawa Institute for Theoretical Physics, Kyoto University,Kitashirakawaoiwake-cho, Sakyo-ku, Kyoto 606-8502, Japan (Dated: August 10, 2018)The shear modulus of jammed frictional granular materials with the harmonic repulsive interactionunder an oscillatory shear is numerically investigated. It is confirmed that the storage modulus,the real part of the shear modulus, for frictional grains with sufficiently small strain amplitude γ discontinuously emerges at the jamming transition point. The storage modulus for small γ differsfrom that of frictionless grains even in the zero friction limit, while they are almost identical witheach other for sufficiently large γ , where the transition becomes continuous. The stress-strain curveexhibits a hysteresis loop even for a small strain, which connects a linear region for sufficiently smallstrain to another linear region for larger strain. We propose a new scaling law to interpolate betweenthe states of small and large γ . PACS numbers: 45.70.-n,05.70.Jk,81.40.Jj
Introduction.–
When the packing fraction φ exceedsa critical value φ c , amorphous materials consisting of re-pulsive particles such as granular materials, colloidal sus-pensions, foams, and emulsions turn into jammed solidswhich have rigidity. Such a transition, known as the jam-ming transition, has been the subject of extensive studiesover the last two decades [1, 2]. For frictionless grains,the pressure increases continuously from φ c , while the co-ordination number Z exhibits a discontinuous transitionat φ c in the hard core limit [3, 4].An assembly of frictionless grains under a simple shearexhibits a rheological continuous transition: the viscos-ity diverges as φ approaches φ c below, while the yieldstress increases continuously above φ c [5–27]. The jam-ming transition is also characterized by the appearanceof rigidity under an oscillatory strain above φ c . For suf-ficiently small strain, the elastic modulus of frictionlessgrains is independent of the strain and the critical ex-ponents for the jamming transition depend on the typeof the local interaction [3, 4, 28]. We call this regimethe linear response regime. For large strain, recent stud-ies [29–34] have revealed that the real part of the shearmodulus, the storage modulus G , of frictionless particlesdecreases with increasing strain as a result of nonlinearresponse because of slip avalanches [35, 36]. The presentauthors have proposed a scaling law of G to interpolatebetween the linear and the nonlinear responses [30]. Notethat the storage modulus as the ratio to the stress to thestrain can be used even in the nonlinear response regime.Most of previous studies, however, assume that grainsare frictionless, though it is impossible to remove contactfriction in experiments of dry granular particles and thefriction causes drastic changes in rheology such as a dis-continuous shear-thickening with a hysteresis loop near ∗ [email protected] φ c [37–52]. Somfai et al. [53] have numerically investi-gated the elastic moduli of a frictional system with theaid of the density of state, and found that G/B with thebulk modulus B in the linear response regime is propor-tional to ∆ Z , the excess coordination number relative tothe isostatic value. They have also found that ∆ Z dis-continuously appears at φ c for frictional grains [53]. Thissuggests that frictional grains with the harmonic repul-sive interaction exhibit a discontinuous change of G at φ c even in the frictionless limit contrast to a continuouschange for frictionless grains [3, 4, 28].The difference of the linear elasticity between the fric-tionless limit and the frictionless case casts doubt on theexpectation that the scaling laws for the elasticity of fric-tionless grains can be confirmed in experiments of grainswith sufficiently small friction coefficient. The accessibleshear strain in the experiment [29], however, is too large,and does not correspond to the previous theoretical stud-ies [3, 4, 30]. Also, little is known on the nonlinear elas-ticity of jammed frictional grains. To clarify the effectsof the contact friction on the linear and the nonlinearelasticities, we numerically investigate the shear modu-lus of two-dimensional frictional granular materials nearthe jamming point under an oscillatory shear. Setup of Simulation.–
Let us consider a two-dimensional assembly of N frictional granular particles.They interact according to the Cundall-Strack modelwith an identical mass density ρ in a square periodic boxof linear size L [54]. The normal repulsive interactionforce F (n) is given by F (n) = k (n) r with the compres-sion length r and the normal spring constant k (n) , whilethe tangential contact force F (t) in quasi-static motionis constrained by the Coulomb criterion | F (t) | ≤ µF (n) : F (t) = k (t) δ (t) in the “stick region” for | δ (t) | < µF (n) /k (t) with the tangential spring constant k (t) and the tangen-tial displacement δ (t) , while | F (t) | remains µF (n) in the“slip region” for | δ (t) | ≥ µF (n) /k (t) . To avoid crystalliza-tion, we use a bi-disperse system which includes equalnumber of grains of the diameters d and d / .
4, respec-tively. The system is subjected to an oscillatory shearwith the shear strain γ ( t ) = γ { − cos( ωt ) } , where γ and ω are the strain amplitude and the angular frequency,respectively. Then, we measure the storage modulus de-fined by [55] G ( γ , µ, φ ) = − ωπ Z π/ω dtσ ( t ) cos( ωt ) /γ , (1)where σ ( t ) is the shear stress. Note that the loss mod-ulus exhibits only a linear dependence on ω and its φ -dependence is relatively small [56]. We also point outthat G is almost independent of ω when the time pe-riod of the oscillatory shear is sufficiently larger than therelaxation time of the configuration of the grains. We,thus, focus on the dependence of G only on γ , µ , and φ for sufficiently small ω . Further details of our model andthe ω -dependence of the storage and the loss moduli areshown in Ref. [56]. Storage modulus for a given packing fraction.–
To be-gin with, we study the dependence of the storage moduluson the friction coefficient µ . In Fig. 1, we plot G against γ for various µ at φ = 0 .
87. For each µ , G is almostindependent of γ in the linear response regime, while G decreases with γ in the nonlinear response regime. Inthe linear response regime, G for µ > µ , but differs from that for µ = 0 [53]. In thenonlinear response regime, the storage moduli for µ = 0and µ > µ decreases, and G for µ = 10 − and 10 − have second plateaus. We confirm that storage modulifor various φ depend on the order of the limits:lim µ → +0 lim γ → +0 G ( γ , µ, φ ) = lim γ → +0 G ( γ , µ = 0 , φ ) . (2)The dependence of G on µ and γ in Fig. 1 can beexplained from the stress-strain curves shown in Fig.2, where we plot the intrinsic shear stress σ ( γ ) − σ (0)against the shear strain γ for various µ at γ = 10 − and φ = 0 .
87. For each µ , σ ( γ ) − σ (0) is proportionalto γ for sufficiently small γ . The proportionality con-stant for µ > µ , but is larger thanthat for µ = 0, which means Eq. (2). The range ofthe linear response becomes narrower as µ decreases, andthe stress-strain curve at µ = 10 − exhibits a hysteresisloop which connects the first linear region for sufficientlysmall γ with the second linear region. The gradient forthe second linear region is identical to that of the fric-tionless grains, which results in the second plateau shownin Fig. 1. This behavior has its origin in the change ofthe tangential friction F (t) from the “stick region” to the“slip region”. See Ref. [56] for the γ -dependence of thestress-strain curve. Scaling law.–
Next, we examine how the storage mod-ulus depends on φ . In Fig. 3, we plot G against φ forvarious γ at µ = 0 .
01. For the smallest strain ampli-tude ( γ = 10 − ), G exhibits a discontinuous transition ✶✵(cid:0)✷✶✵(cid:0)✁✶✵✂ ✶✵(cid:0)✻ ✶✵(cid:0)✹ ✶✵(cid:0)✷●❂❦✭♥✮ ✌✂✖ ✄ ✶✵(cid:0)✷✖ ✄ ✶✵(cid:0)✸✖ ✄ ✶✵(cid:0)✹✖ ✄ ✶✵(cid:0)✺✖ ✄ ✵ FIG. 1: (Color online) The storage modulus G against γ for µ = 10 − , 10 − , 10 − , 10 − , and 0 at φ = 0 . ✵✷ ✂ ✶✵(cid:0)✼✹ ✂ ✶✵(cid:0)✼✻ ✂ ✶✵(cid:0)✼ ✵ ✶ ✂ ✶✵(cid:0)✁ ✷ ✂ ✶✵(cid:0)✁✭✛✭✌✮✄✛✭☎✮✮❂❦✆♥✝ ✞✖ ✟ ✵✖ ✟ ✶✵(cid:0)✁✖ ✟ ✶✵(cid:0)✺✖ ✟ ✶✵(cid:0)✠ FIG. 2: (Color online) The intrinsic shear stress σ ( γ ) − σ (0)against γ with µ = 10 − , 10 − , 10 − , and 0 for γ = 10 − and φ = 0 . at φ c ≃ .
84. As γ increases, the discontinuity at φ c de-creases and the transition becomes asymptotically con-tinuous, where G is approximately proportional to φ − φ c .Here, we propose a new scaling law for the storagemodulus near the transition point: G ( γ , µ, φ ) = G (lin) ( µ, φ ) F (cid:18) γ µ b { φ − φ c ( µ ) } b (cid:19) , (3)where F ( x ) is the scaling function satisfyinglim x → F ( x ) = 1 , lim x →∞ F ( x ) ∼ x − c (4) ✵✵✿✶✵✿✷✵✿✸ ✵✿✽✹ ✵✿✽✻ ✵✿✽✽●❂❦✭♥✮ ✣ ✌(cid:0) ✁ ✶✵✂✄✌(cid:0) ✁ ✶✵✂☎✌(cid:0) ✁ ✶✵✂✺✌(cid:0) ✁ ✶✵✂✆ FIG. 3: (Color online) The storage modulus G against φ for γ = 10 − , 10 − , 10 − , and 10 − at µ = 0 . with exponents b , b , and c , and we have introduced G (lin) ( µ, φ ) ≡ lim γ → +0 G ( γ , µ, φ ) . (5)In Eq. (3), we have used the jamming point φ c ( µ ) de-pending on µ [56]. Figure 4 confirms the validity of thescaling plot characterized by Eq. (3), where we havedetermined G (lin) ( µ, φ ) by the extrapolation of the limit γ → +0 using the data for γ ≥ − . The criticalexponents used in Eq. (3), are given by [56] b = 1 . ± . , b = 0 . ± . , c = 1 . ± . . (6)It should be noted that the scaling law cannot be appliedto the region of the second plateau in Fig. 1In Eq. (3), we have assume that the critical strain γ c characterizing the crossover from the linear to the non-linear response regimes is proportional to µ b ( φ − φ c ) b .The exponents b and b in Eq. (6) may be understoodas follows. First, γ c is expected to satisfy γ c ∼ δ (t) c with the critical tangential displacement δ (t) c characteriz-ing the change of the tangential friction F (t) to the “slipregion”. Then, we deduce δ (t) c ∼ µF (n) with the aver-age contact force F (n) ∼ ( φ − φ c ) for grains with theharmonic repulsive interaction [10, 11]. Thanks to theabove relations, we obtain γ c ∼ µ b ( φ − φ c ) b with b = 1and b = 1, which are not far away from the estimatedvalues in Eq. (6). We, however, do not identify the rea-son why the evaluation of b deviates a little from thenumerical value.It should be noted that the scaling form of Eq. (3) isanalogous to that for frictionless case proposed in Ref.[30], though the µ -dependence is not included and c is1 / ✶✵(cid:0)✷✶✵(cid:0)✁✶✵✂✶✵(cid:0)✹ ✶✵(cid:0)✷ ✶✵✂ ✶✵✷●❂●✭❧✐♥✮ ✌✂✄❢✖❜☎ ✆✣ ✝ ✣❝✞❜✟ ❣✆✶✵(cid:0)✹❀ ✵✿✵✶✞✆✶✵(cid:0)✺❀ ✵✿✵✶✞✆✶✵(cid:0)✻❀ ✵✿✵✶✞✆✶✵(cid:0)✹❀ ✵✿✵✠✞✆✶✵(cid:0)✺❀ ✵✿✵✠✞✆✶✵(cid:0)✹❀ ✵✿✶✞✆✶✵(cid:0)✺❀ ✵✿✶✞✆✶✵(cid:0)✹❀ ✶✿✵✞ FIG. 4: (Color online) Scaling plot of G characterizedby Eq. (3). Each symbol in the legend is character-ized by ( γ , µ ), but each one has the data for φ − φ c =0 . , . , . , . .
01. The solid line repre-sents the scaling function (S16) in Ref. [56]. represents the crossover from the stick to the slip branch,while the previous scaling deals with the crossover fromthe slip to the avalanche branch.We have also confirmed that the storage modulus G (lin) ( µ, φ ) in the linear response regime exhibits a dis-continuous transition at φ c . To give further evidence, weplot the storage modulus at φ c defined as G ( µ ) ≡ lim φ → + φ c G (lin) ( µ, φ ) (7)against µ in Fig. 5, where G ( µ ) has a maximum value inthe frictionless limit. We have also found that G (lin) ( µ, φ )satisfies G (lin) ( µ, φ ) − G ( µ ) ∝ { φ − φ c ( µ ) } a (8)with an exponent a = 0 . ± .
01. It should be noted thatthe previous numerical results suggest a = 1 / γ → +0 G ( γ , µ, φ ) = G ( µ ) + A ( µ ) { φ − φ c ( µ ) } a (9)for the discontinuous transition in the linear responseregime andlim φ → + φ c G ( γ , µ, φ ) ∝ G ( µ ) (cid:18) µ b ( φ − φ c ( µ )) b γ (cid:19) c (10)for the continuous transition in the nonlinear responseregime, where A ( µ ) depends only on µ . We note thatthe scaling law, Eq. (10), with b = 0 .
90 and c ≃ . ✵✵✿✶✵✿✷ ✵ ✵✿✺ ✶● (cid:0)❂❦✭♥✮ ✖ FIG. 5: (Color online) The storage modulus G ( µ ) for µ > φ c in the linear response region. The smallest value of µ inthis plot is 0 .
01. The open circle at the origin represents theresult for frictionless grains. indicates G ∼ ( φ − φ c ) b c ≃ ( φ − φ c ) in the nonlinearresponse regime. This might suggest that the scaling of G is identical to that of the pressure P [30]. Discussion.–
Let us discuss our results. Tighe reportedthat the complex shear modulus G ∗ for a model of emul-sions satisfies G ∗ ∼ ω / near φ c in the linear responseregime [61]. On the contrary, we have confirmed that thestorage modulus is independent of ω , and the loss modu-lus is proportional to ω in the linear response regime forall of the packing fractions [56]. It should be noted thatthe static contact friction and the inertia of the particlesexist in our system, but they are not involved in Ref.[61], which might lead to the different dependence of G ∗ .In this Letter, we have proposed the scaling law forthe shear modulus of grains with the harmonic repulsiveinteraction. From the analogy of the frictionless case [30],we expect the exponent b = 3 / γ c ∼ δ (t) c ∼ µF (n) with F (n) ∼ ( φ − φ c ) / . The confirmation of this conjecturewill be the subject of further study.There are some studies to focus on the role of perco- lation to the jamming [62–65]. In the system exhibitinga conventional percolation, the critical exponents for thestorage modulus depends on the spatial dimension [66],but the exponents of the jamming transition are indepen-dent of the dimensionality, at least, for frictionless grains[11]. Further careful study on the relation between per-colation and jamming should be necessary.An important finding of this Letter is the scaling law(3) interpolating between the linear and the nonlinear re-sponse regimes for frictional grains. Somfai et al. foundthat the elastic moduli of frictional grains in the linear re-sponse regime is proportional to the excess coordinationnumber ∆ Z [53]. Based on this result, a recent reviewpaper suggested that there is no scaling of the mechani-cal properties as a function of φ − φ c for frictional grainsbecause ∆ Z of the frictional grains remains finite evenat φ c [2]. Our scaling plot based on Eq. (3), however,gives a counter evidence of the existence of the scalinglaw. We also note that Ref. [53] have ignored the changeof the tangential contact force to the “slip region”, butthis effect becomes significant in the vicinity of φ c for fi-nite γ as can be seen in Eq. (3). To our knowledge, thiseffect has not been indicated in previous studies. Concluding remarks.–
We have numerically investi-gated the frictional granular particles. The storage mod-ulus in the linear response regime for frictional grainsdiffers from that for frictionless grains even in the zerofriction limit, whereas they are almost identical in thenonlinear response regime. This dependence on the tan-gential friction has been explained from the stress-straincurve shown in Fig. 2. We have also proposed thenew scaling law that interpolates between the discontinu-ous transition for infinitesimal strain and the continuousone for finite strain. This scaling law has been verifiedthrough our simulation.The authors thank O. Dauchot, H. Yoshino, T. Ya-maguchi, F. van Wijland, and K. Miyazaki for fruitfuldiscussions. This work is partially supported by theGrant-in-Aid of MEXT for Scientific Research (GrantNo. 16H04025 and No. 25800220). One of the authors(M.O.) appreciates the warm hospitality of Yukawa Insti-tute for Theoretical Physics at Kyoto University and thediscussions during the YITP workshops “New Frontiersin Non-equilibrium Physics 2015” and “Avalanches, plas-ticity, and nonlinear response in nonequilibrium”, whichhelped to complete this work. [1] A. J. Liu and S. R. Nagel, Nature , 21 (1998).[2] M. van Hecke, J. Phys.: Condens. Matter , 033101(2009)[3] C. S. O’Hern, S. A. Langer, A. J. Liu, and S. R. Nagel,Phys. Rev. Lett. , 075507 (2002).[4] C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel,Phys. Rev. E , 011306 (2003). [5] P. Olsson and S. Teitel, Phys. Rev. Lett. , 178001(2007).[6] T. Hatano, M. Otsuki, and S. Sasa, J. Phys. Soc. Jpn. ,023001 (2007).[7] T. Hatano, J. Phys. Soc. Jpn. , 123002 (2008).[8] B. P. Tighe, E. Woldhuis, J. J. C. Remmers, W. van Saar-loos, and M. van Hecke, Phys. Rev. Lett. , 088303 (2010).[9] T. Hatano, Prog. Theor. Phys. Suppl. , 143 (2010).[10] M. Otsuki and H. Hayakawa, Prog. Theor. Phys. ,647 (2009).[11] M. Otsuki and H. Hayakawa, Phys. Rev. E , 011308(2009).[12] M. Otsuki, H. Hayakawa, and S. Luding, Prog. Theor.Phys. Suppl. , 110 (2010).[13] 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. , 175701 (2010).[14] P. Olsson and S. Teitel, Phys. Rev. E , 030302(R)(2011).[15] D. V˚agberg, P. Olsson, and S. Teitel, Phys. Rev. E ,031307 (2011).[16] M. Otsuki and H. Hayakawa, Prog. Theor. Phys. Suppl. , 129 (2012).[17] A. Ikeda, L. Berthier, and P. Sollich, Phys. Rev Lett. , 108001(2012).[19] E. DeGiuli, G. D¨uring, E. Lerner, and M. Wyart, Phys.Rev. E , 062206 (2015).[20] D. V˚agberg, P. Olsson, and S. Teitel Phys. Rev. E ,052902 (2016).[21] F. Boyer, E. Guazzelli, and O. Pouliquen, Phys. Rev.Lett. , 188301 (2011).[22] M. Trulsson, B. Andreotti, and P. Claudin, Phys. Rev.Lett. , 118305 (2012).[23] B. Andreotti, J.-L. Barrat, and C. Heussinger, Phys. Rev.Lett. , 105901 (2012).[24] E. Lerner, G. D¨uring, and M. Wyart, Proc. Natl. Acad.Sci. U.S.A , 4798 (2012).[25] D. V˚agberg, P. Olsson, and S. Teitel Phys. Rev. Lett. , 148002 (2014).[26] T. Kawasaki, D. Coslovich, A. Ikeda, and L. Berthier,Phys. Rev. E , 012203 (2015).[27] K. Suzuki and H. Hayakawa, Phys. Rev. Lett. ,098001 (2015).[28] M. Wyart, Annales de Physique , 3 (2005).[29] C. Coulais, A. Seguin, and O. Dauchot, Phys. Rev. Lett. , 198001 (2014).[30] M. Otsuki and H. Hayakawa, Phys. Rev. E , 042202(2014).[31] M. S. van Deen, J. Simon, Z. Zeravcic, S. Dagois-Bohy, B.P. Tighe, and M. van Hecke, Phys. Rev. E , 020202(R)(2014).[32] C. P. Goodrich, A. J. Liu, and J. P. Sethna, Proc. Natl.Acad. Sci. USA. , 5450 (2016).[35] K. Dahmen, D. Erta¸s, and Y. Ben-Zion, Phys. Rev. E , 1494 (1998).[36] K. Dahmen, Y. Ben-Zion, and J. T. Uhl, Nat. Phys. ,554 (2011).[37] M. Otsuki and H. Hayakawa, Phys. Rev. E 83, 051301(2011).[38] D. Bi, J. Zhang, B. Chakraborty and R. Behringer, Na-ture , 355 (2011).[39] S. Chialvo, J. Sun, and S. Sundaresan, Phys. Rev. E ,021305 (2012).[40] E. Brown and H. M. Jaeger, Phys. Rev. Lett. , 086001 (2009).[41] R. Seto, R. Mari, J. F. Morris, and M. M. Denn, Phys.Rev. Lett. , 218301 (2013).[42] N. Fernandez, R. Mani, D. Rinaldi, D. Kadau, M. Mos-quet, H. Lombois-Burger, J. Cayer-Barrioz, H. J. Her-rmann, N. D. Spencer, and L. Isa, Phys. Rev. Lett. ,108301 (2013).[43] C. Heussinger, Phys. Rev. E , 050201 (2013).[44] M. M. Bandi, M. K. Rivera, F. Krzakala and R.E. Ecke,Phys. Rev. E , 042205 (2013).[45] M. P. Ciamarra, R. Pastore, M. Nicodemi, and A.Coniglio, Phys. Rev. E , 041308 (2011).[46] R. Mari, R. Seto, J. F. Morris, and M. M. Denn, J. Rheol. , 1693 (2014).[47] M. Grob, C. Heussinger, and A. Zippelius, Phys. Rev. E , 050201 (2014).[48] T. Kawasaki, A. Ikeda, and L. Berthier, EPL , 28009(2014).[49] M. Wyart and M. E. Cates, Phys. Rev. Lett. , 098302(2014).[50] M. Grob, A. Zippelius, and C. Heussinger, Phys. Rev. E , 030901 (2016).[51] V. Magnanimo, L. La Ragione, J. T. Jenkins, P. Wang,and H. A. Makse, EPL , 34006 (2008).[52] H. Hayakawa and S. Takada, arXiv:1611.07925.[53] E. Somfai, M. van Hecke, W. G. Ellenbroek, K.Shundyak, and W. van Saarloos, Phys. Rev. E ,020301(R) (2007).[54] P. Cundall and O. D. L. Strack, Geotechnique , 47(1979).[55] M. Doi and S. F. Edwards, The Theory of Polymer Dy-namics (Oxford University Press, Oxford, 1990).[56] See Supplemental Material, which includes Ref. [57–60],for the details of the analysis and supplementary simula-tion results.[57] D. J. Evans and G. P. Morriss,
Statistical Mechanicsof Nonequilibrium Liquids
Numerical Recipes , 3rd ed. (Cambridge Uni-versity Press, Cambridge, 2007).[59] L. E. Silbert, Soft Matter , 2918 (2010)[60] L. E. Silbert, D. Erta¸s, G. S. Grest, T. C. Halsey, and D.Levine, Phys. Rev. E , 031304 (2002).[61] B. P. Tighe, Phys. Rev. Lett. , 158303 (2011).[62] G. Lois, J. Blawzdziewicz, and C. S. O’Hern, Phys. Rev.Lett. , 028001 (2008).[63] T. Shen, C. S. O’Hern, M. D. Shattuck, Phys. Rev. E ,011308 (2012).[64] L. Kovalcinova, A. Goullet, and L. Kondic, Phys. Rev. E , 032204 (2015).[65] S. Henkes, D. A. Quint, Y. Fily, and J. M. Schwarz, Phys.Rev. Lett. , 028301 (2016).[66] S. Torquato, Random Heterogeneous Materials: Mi-crostructure and Macroscopic Properties , 2nd ed.(Springer, New York, 2005).
Supplemental Materials: **
I. INTRODUCTION
In this Supplemental Materials, we present details andadditional results of our model. We explain the setup ofour simulation in Sec. II. In Sec. III, we demonstratethe dependence of the complex shear modulus on theangular frequency ω . Section IV deals with the stress-strain curves for various γ . In Sec. V, we explain themethod to determine the exponents for the scaling lawproposed in the main text. We present the method todetermine the transition point depending on the frictioncoefficient µ in Sec. VI. In Sec. VII, we show numericalresults for the storage modulus in the linear responseregime. II. DETAILS OF OUR SIMULATION
In this section, we explain our setup and model. Theposition, the velocity, and the angular velocity of thegrain i are, respectively, denoted by r i , v i , and ω i ˆ z ,where we have introduced the unit vector ˆ z parallel to z axis (perpendicular to the considering plane). The con-tact force F ij between the grain i and the grain j consistsof the normal part F (n) ij and the tangential part F (t) ij as F ij = F (n) ij + F (t) ij . The normal contact force F (n) ij isgiven by F (n) ij = F (n , el) ij + F (n , vis) ij , (S1)where F (n , el) ij = k (n) ( d ij − | r ij | )Θ( d ij − | r ij | ) n ij , (S2) F (n , vis) ij = − η ( n ) v (n) ij Θ( d ij − | r ij | ) n ij (S3)with the normal spring constant k ( n ) and the normal vis-cous constant η ( n ) . Here, r ij , n ij , v (n) ij , and d ij are,respectively, given by r ij = r i − r j , n ij = r ij / | r ij | , v (n) ij = ( v i − v j ) · n ij , and d ij = ( d i + d j ) / d i of the grain i . In Eqs. (S2) and (S3), we haveintroduced the the Heaviside step function Θ( x ) definedby Θ( x ) = 1 for x ≥ x ) = 0 otherwise.Similarly, the tangential friction force F (t) ij is given by F (t) ij = min (cid:16) ˜ F (t) ij , µ | F (n , el) ij | (cid:17) sign (cid:16) ˜ F (t) ij (cid:17) t ij , (S4)where min( a, b ) selects the smaller one between a and b ,sign( x ) is 1 for x > − x <
0, and ˜ F (t) ij is givenby ˜ F (t) ij = k (t) δ (t) ij − η (t) v (t) ij (S5) with t ij = ( − r ij,y / | r ij | , r ij,x / | r ij | ). Here, k (t) , η (t) , and µ are the tangential spring constant, the tangential vis-cous constant, and the friction coefficient, respectively.The tangential velocity v (t) ij and the tangential displace-ment δ (t) ij are, respectively, given by v (t) ij = ( v i − v j ) · t ij − ( d i ω i + d j ω j ) / , (S6) δ (t) ij = Z stick dt v (t) ij , (S7)where “stick” on the integral implies that the integral isonly performed when the condition | ˜ F (t) ij | < µ | F (n , el) ij | issatisfied. Additionally, we have introduced the torque T i on the grain i as T i = − X j d i F (t) ij · t ij . (S8)In this model, we apply an oscillatory shear along the y direction under the Lees-Edwards boundary condition[S1]. The SLLOD equations of motion are used to stabi-lize the uniform sheared state as d r i dt = p i m i + ˙ γ ( t ) r i,y ˆ x , (S9) d p i dt = X j = i F ij − ˙ γ ( t ) p i,y ˆ x , (S10) I i dω i dt = T i , (S11)where we have introduced the time dependent shear rate˙ γ ( t ), the peculiar momentum p i defined by Eq. (S9),the unit vector parallel to the x -direction ˆ x , the mass m i = πρd i /
4, and the moment of inertia I i = m i d i / φ I = 0 .
75. Afterrelaxing the system to a mechanical equilibrium state, wecompress the system without shear in small steps untilthe packing fraction reaches a given value φ . In eachstep, we change the linear system size and the positionof the grain i as L ( n s +1) = L ( n s ) s φ ( n s ) φ ( n s +1) , (S12) r ( n s +1) i = r ( n s ) i s φ ( n s ) φ ( n s +1) , (S13)and relax the grains to the mechanical equilibrium state.Here, L ( n s ) , r ( n s ) i , and φ ( n s ) denote the system size, theposition of the grains i , and the packing fraction at the n s -th step, respectively. The increment of the packingfraction is defined by ∆ φ ≡ φ ( n s +1) − φ ( n s ) . Here, we re-gard the state for T < T th as the mechanical equilibriumstate, where we have introduced the kinetic temperature T ≡ P i | p i | / (2 m i N ) and a threshold T th . It should benoted that the origin of the coordinate axes is located atthe center of the system.Let us summarize a set of parameters used in oursimulation. We use k (t) = k (n) and η (n) = η (t) = p m k (n) , where m is the mass of a grain with the di-ameter d . This set of the parameters corresponds tothe constant restitution coefficient e = 0 . t = 0 . τ ,where τ is the characteristic time of the stiffness, i.e., τ = p m /k (n) . The number N of the grains is 4000. Wehave checked that the shear modulus is independent of N for N ≥ T th = 10 − ( k (n) d ),∆ φ = 10 − , and ω = 10 − τ − in the main text, but wehave examined the ω -dependence of the complex mod-ulus in the next section. We have confirmed that ourresults in the main text are insensitive to T th , ∆ φ and ω if they are sufficiently small. Here, the storage modulus G is calculated from Eq. (1) in the main text with theaid of the shear stress σ ( t ) given by σ ( t ) = − L N X i X j>i r ij,x ( t ) F ij,y ( t ) − L N X i p i,x ( t ) p i,y ( t ) m i . (S14) III. FREQUENCY DEPENDENCE OFSTORAGE AND LOSS MODULI
This section deals with the dependence of the complexshear modulus G ∗ = G + iG ′′ on ω in the linear responseregime. In Fig. S1, we show the storage modulus G against ω for various φ at µ = 1 . γ = 10 − . Asshown in Fig. S1, G for each φ is almost independentof ω for ωτ < − . Hence, we have not discussed the ω -dependence of G in the main text.Figure S2 exhibits the loss modulus G ′′ against ω forvarious φ at µ = 1 . γ = 10 − . Here, G ′′ is givenby G ′′ = ωπ Z π/ω dtσ ( t ) sin( ωt ) /γ . (S15)As shown in Fig. S2, G ′′ is almost proportional to ω .These ω dependences of G and G ′′ indicate that the rhe-ological properties in our model are essentially describedby the Kelvin-Voigt viscoelasticity. This result is rea-sonable because the Cundall-Strack model relies on theKelvin-Voigt model [S2]. Moreover, the φ -dependence of G ′′ is not clearly visible. Hence, we have investigatedonly G in the main text. ✶✵(cid:0)✷✶✵(cid:0)✁✶✵✂ ✶✵(cid:0)✹ ✶✵(cid:0)✷● ✦✜✣ ✄ ✣❝ ❂ ✵✿✵✸✣ ✄ ✣❝ ❂ ✵✿✵✵✸✣ ✄ ✣❝ ❂ ✵✿✵✵✵✸ FIG. S1: The storage modulus G against ω for φ − φ c =0 . , .
003 and 0 . µ = 1 . γ = 10 − . ✶✵(cid:0)✻✶✵(cid:0)✺✶✵(cid:0)✹✶✵(cid:0)✸✶✵(cid:0)✷✶✵(cid:0)✁✶✵✂ ✶✵(cid:0)✹ ✶✵(cid:0)✷●✄✄ ✦✜✣ ☎ ✣❝ ❂ ✵✿✵✆✣ ☎ ✣❝ ❂ ✵✿✵✵✆✣ ☎ ✣❝ ❂ ✵✿✵✵✵✆ FIG. S2: The loss modulus G ′′ against ω for φ − φ c =0 . , .
003 and 0 . µ = 1 . γ = 10 − . The solidline represents G ′′ ∼ ω . IV. THE STRESS-STRAIN CURVES FORVARIOUS γ In this section, we demonstrate the dependence of thestress-strain curve on γ . Figure S3 exhibits the intrinsicshear stress σ ( γ ) − σ (0) against the shear strain γ forvarious γ at µ = 10 − and φ = 0 .
87. For γ = 10 − , σ ( γ ) − σ (0) is proportional to γ , but the constant of pro-portionality is larger than that of frictionless grains in-dicated by the dotted line. The reason for the differencein the constant of proportionality is that the increase ofthe tangential friction force | F (t) ij | in the “stick region”enlarges σ ( γ ) − σ (0) for µ > γ .As γ increases, the stress-strain curve exhibits a nonlin-ear behavior: a hysteresis loop to connect the first linearregion for sufficiently small γ with the second linear re-gion for γ > − , where the gradient of the second linearregion is identical to that of the frictionless grains. Thesecond linear region originates from the fact that | F (t) ij | remains constant in the “slip region” and does not con-tribute to enlarge σ ( γ ) − σ (0). ✵✷ ✂ ✶✵(cid:0)✼✹ ✂ ✶✵(cid:0)✼✻ ✂ ✶✵(cid:0)✼ ✵ ✸ ✂ ✶✵(cid:0)✁ ✻ ✂ ✶✵(cid:0)✁✭✛✭✌✮✄✛✭☎✮✮❂❦✆♥✝ ✞✞✟ ✠ ✸ ✂ ✶✵(cid:0)✁✞✟ ✠ ✶ ✂ ✶✵(cid:0)✁✞✟ ✠ ✶ ✂ ✶✵(cid:0)✼ FIG. S3: The intrinsic shear stress σ ( γ ) − σ (0) against γ for γ = 3 . × − , 1 . × − , and 1 . × − at µ = 10 − and φ = 0 .
87. The dotted line indicates σ ( γ ) − σ (0) for γ =3 . × − at µ = 0 and φ = 0 . V. METHOD TO DETERMINE THEEXPONENTS
In this section, we explain the method to determinethe exponents in the scaling plot proposed in the maintext. Here, the exponents are determined by Levenberg-Marquardt algorithm [S3], where we have assumed thefunctional form of the scaling function as F ( x ) = 11 + e P Nnn =0 A n (ln x ) n (S16)with the fitting parameters N n = 1, A = − . ± . A = 1 . ± .
02. We have checked the numericalestimation of b and b with N n ≥
2, but A n for n ≥ c in Eq. (4) ofthe main text is approximately equal to c ≈ A = 1 . ✵✵✿✶✵✿✷✵✿✸✵✿✽✹ ✵✿✽✹✺ ✵✿✽✺●❂❦✭♥✮ ✣❙❛♠♣❧(cid:0) ✶❙❛♠♣❧(cid:0) ✷ FIG. S4: The storage modulus G against φ at µ = 0 .
01 and γ = 10 − . “Sample 1” and “Sample 2” indicate the dataobtained from different initial configurations, respectively. ✵✿✽✷✵✿✽✸✵✿✽✹✵✿✽✺ ✵ ✵✿✷ ✵✿✹ ✵✿✻ ✵✿✽ ✶✣ ❝ ✖ FIG. S5: The plots of transition points φ c against µ . VI. DETERMINATION OF TRANSITIONPOINT
In this section, we explain the method to determinethe transition point φ c ( µ ). To estimate φ c ( µ ), we firstprepare assemblies with different packing fractions ob-tained from the same initial configuration of grains withthe friction coefficient µ . Then, we apply the oscillatoryshear to measure G for sufficiently small shear amplitude( γ = 10 − ), and plot G against φ as shown in Fig. S4,which plots the data at µ = 0 .
01 and γ = 10 − ob-tained from different initial configurations. G changesdiscontinuously around φ = 0 .
84, but the critical frac-tion depends on the initial configuration. Then, we de-fine φ c ( µ ) for each initial configuration as the minimumvalue of φ where G for a given initial configuration ex-ceeds G th = 10 − k (n) . In Fig. S5, we plot the averageof φ c ( µ ) against µ . φ c ( µ ) decreases with increasing µ ,which is consistent with the previous numerical results[S4, S5]. VII. THE STORAGE MODULUS IN THELINEAR RESPONSE REGIME
In this section, we investigate the storage modulus inthe linear response regime. Figure S6 plots G ( µ ) againstthe the excess coordination number Z ( µ ) − Z iso for var-ious µ , where Z ( µ ) and Z iso = 3 are the coordinationnumber at φ c and the isostatic value of the coordina-tion number for two dimensional frictional grains, respec-tively. As shown in Ref. [S6], G ( µ ) satisfies G ( µ ) ∝ Z ( µ ) − Z iso . (S17)It is known that Z ( µ ) continuously decreases with in-creasing µ from the isostatic value, 4, for frictionlessgrains [S4–S7]. Z ( µ ) in our model exhibits the iden-tical behavior as shown in Fig. S7. It should be notedthat Z ( µ ) in our system under oscillatory shear is quali-tatively consistent with that of the previous result understeady shear [S5]. These results explain the decrease of G ( µ ) with increasing µ shown in Fig. 5 of the main text. ✵✵✿✶✵✿✷ ✵ ✵✿✺ ✶● (cid:0)❂❦✭♥✮ ❩✁ ✂ ❩✐s♦ FIG. S6: G ( µ ) against Z ( µ ) − Z iso for µ = 0 .
01, 0 .
05, 0 . .
2, 0 .
3, 0 .
5, 0 .
7, and 1 .
0. The solid line represents Eq. (S17).
We have also found that G (lin) ( µ, φ ) satisfies Eq. (8)in the main text, which is verified in Fig. S8. In the inset ✸✹ ✵ ✵✿✺ ✶ ✶✿✺ ✷❩ (cid:0) ✖❖s❝✐✁✁✂t♦r② s❤❡✂r❙t❡✂❞② s❤❡✂r FIG. S7: The coordination number Z at the transitionpoint φ c against µ for the present system under oscillatoryshear and the system in Ref. [S5] under steady shear. of Fig. S8, we plot the exponent a evaluated by the leastsquares method against µ , which is almost independentof µ and estimated as a = 0 . ± . ✶✵(cid:0)✸✶✵(cid:0)✷✶✵(cid:0)✁ ✶✵(cid:0)✹ ✶✵(cid:0)✸✵✶ ✵ ✵✿✺ ✶✏ ●✭❧✐♥✮ ✂● ✄✑ ❂❦✭♥✮ ✣ ☎ ✣❝✖ ✆ ✵✿✵✶✖ ✆ ✵✿✶✖ ✆ ✶✿✵ ❛ ✖ FIG. S8: G (lin) ( µ, φ ) − G ( µ ) against φ − φ c ( µ ) for µ = 0 . .
05, 0 .
1, 0 .
5, and 1 .
0. The solid line represents Eq. (8) in themain text with a = 1 /
2. (Inset) The exponent a evaluated bythe least squares method against µ . The solid line represents a = 1 / Statistical Mechanics ofNonequilibrium Liquids , 47(1979).[S3] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B.P. Flannery, Numerical Recipes , 3rd ed. (Cambridge Uni-versity Press, Cambridge, 2007).[S4] L. E. Silbert, Soft Matter , 2918 (2010) [S5] M. Otsuki and H. Hayakawa, Phys. Rev. E 83, 051301(2011).[S6] E. Somfai, M. van Hecke, W. G. Ellenbroek, K.Shundyak, and W. van Saarloos, Phys. Rev. E ,020301(R) (2007).[S7] L. E. Silbert, D. Erta¸s, G. S. Grest, T. C. Halsey, and D.Levine, Phys. Rev. E65