Role of local response in manipulating the elastic properties of disordered solids by bond removal
RRole of local response in manipulating the elastic properties of disordered solids bybond removal
Daniel Hexner ∗ The James Franck Institute and Department of Physics, The University of Chicago,Chicago, IL 60637, USA and Department of Physics and Astronomy,The University of Pennsylvania, Philadelphia, PA, 19104, USA
Andrea J. Liu
Department of Physics and Astronomy, The University of Pennsylvania, Philadelphia, PA, 19104, USA
Sidney R. Nagel
The James Franck and Enrico Fermi Institutes and The Department of Physics,The University of Chicago, Chicago, IL 60637, USA
We explore the range over which the elasticity of disordered spring networks can be manipulated bythe removal of selected bonds. By taking into account the local response of a bond, we demonstratethat the effectiveness of pruning can be improved so that auxetic (i.e., negative Poisson’s ratio)materials can be designed without the formation of cracks even while maintaining the global isotropyof the network. The bulk, B , and shear, G , moduli scale with the number of bonds removed andwe estimate the exponents characterizing these power laws. We also find that there are spatialcorrelation lengths in the change of B and G upon removing different bonds that diverge as thenetwork approaches the isostatic limit where the excess coordination number ∆ Z → . I. INTRODUCTION
Manipulating the elastic properties of solids is an im-portant problem with broad applications [1, 2]. Themost common approach in designing mechanical meta-materials is based on a periodically repeating unit cellthat is carefully constructed to yield a given elastic prop-erty or function. Recently, a novel design principle hasbeen introduced based on “pruning” disordered springnetworks [3]. This exploits the broad distribution of howdifferent bonds contribute to the elastic moduli in suchsystems – by selectively removing a bond that contributesmore to one modulus than to another, one can prune asystem to achieve desired elastic properties. Disorderprovides two clear advantages over periodic lattices: 1)disordered systems are isotropic on large scales; 2) disor-der allows materials to be designed with inhomogeneousand even local responses [4–7].To demonstrate the potential flexibility that pruningprovides, consider the effect on an elastic modulus of theremoval of a single bond in a spring network with N b bonds. The modulus could characterize the cost of aglobal deformation, such as compression, or a local defor-mation, such as the pinching together of two nodes. If thesystem is periodic, removing a bond in one unit cell willresult in the same change in the elastic moduli as the re-moving the equivalent bond in any other unit cell. How-ever, for a disordered network, removing different bondsleads to different responses. Naively, the removal of thefirst bond results in N b different possible responses while ∗ [email protected] removing N r bonds leads to N b ! /N r ! ( N b − N r )! whicheven for small N r can be an enormous number.To realize this large range of possible designs it is im-portant for the system to obey two properties. First,removing a bond, i , affects different moduli differently inan uncorrelated way. Indeed this was shown to be thecase for bulk deformations: removal of bond i changesthe bulk modulus by ∆ B i and the shear modulus by ∆ G i , where ∆ B i and ∆ G i have nearly vanishing cor-relations [8]. Second, the change in moduli upon theremoval of a bond must have a broad range. This is alsothe case for bulk deformations, where ∆ B i and ∆ G i atsmall values scale as a power-law [8].The systems we design are based on disordered net-works derived from jammed packings [9, 10]. Soft repul-sive spherical particles are placed randomly in space andthe energy is minimized to attain force balance. The cen-ters of the spheres are then connected by springs to forma network and the equilibrium spring length is set to theequilibrium distance between nodes, thus removing allstresses. For simplicity, all the spring constants, k , arechosen to be equal. We characterize ensembles of suchnetworks by the coordination number per node, Z , andthe excess coordination number, ∆ Z = Z − Z c , where Z c is the critical value of Z at which rigidity is lost: in aninfinite system, Z c = 2 d [11–13] .If a bond length, r i , between two nodes is differentfrom the equilibrium length, r i , there is an energy costof kδ r i and a tension τ i = kδr i where δr i ≡ r i − r i .Since the networks are initially unstressed, compressingthe system results in an energy B(cid:15) B , where B is thebulk modulus and (cid:15) B is the compression strain; similarlythe energy cost of a shear is G(cid:15) G where G is the shearmodulus and the (cid:15) G is the shear strain. Since the energy a r X i v : . [ c ond - m a t . s o f t ] A ug is additive, U = (cid:80) i kδ r i this allows us to decomposethe bulk and shear modulus into a sum over their singlebond contributions: B = (cid:80) i B i and G = (cid:80) i G i where B i and G i are the contributions to B and G respectivelyof bond i [3].In an isotropic system, the Poisson’s ratio, ν , isa monotonically decreasing function of G/B . In themethod introduced by Ref. [3], bonds are pruned thattarget either G or B to yield the desired value of ν . Toattain, for example, a large G/B , two strategies may beenvisioned. Bonds that contribute a large amount to thebulk modulus (large values of B i ) can be successively re-moved; due to the relatively weak correlations between B i and G i , this results in only a moderate decrease in G but a large drop in B . Alternatively, bonds that con-tribute little to G (small values of G i ) can be removed sothat G does not change appreciably but again due to theweak correlations, B decreases more steeply. The firstapproach has the undesired effect of creating cracks [14]since the removed bonds carry a lot of stress under thedeformation, which is distributed to its neighbors whenremoved.We recently pointed out [8] that there is a differencebetween the contribution of a given bond i to the mod-ulus, M i , and the change of the modulus if bond i isremoved, ∆ M i . Clearly, the evolution of the M underpruning depends on the latter quantity. Ref. [3] uses B i and G i as proxies for predicting ∆ B i and ∆ G i in orderto tune the values of B and G . This approach was quitesuccessful–it was found that a large G/B can be attainedby removing the bonds with the largest B i [3]. Similarly,a very small G/B can be attained by either removingthe bonds with maximal G i or minimal B i . However, re-moving the minimal G i strategy fails and both B and G change in a correlated manner as the bonds are removed.Here, we show that this failure results from the use of G i as a proxy for ∆ G i , and that when minimal ∆ G i bonds are removed, one can obtain very large values of G/B . Moreover, consideration of ∆ M i instead of M i al-lows us to estimate theoretically the scaling exponents of G/B as a function of the number of bonds removed. Wealso report correlations in ∆ B i and ∆ G i as a function ofthe distance between bonds with correlation lengths ξ ∆ B i and ξ ∆ G i respectively. These correlation lengths divergeas power laws as the network connectivity decreases to-wards the rigidity threshold ∆ Z = 0 . II. TUNING
G/B
Here we tune
G/B by removing bonds based on ∆ B i and ∆ G i . We explore different pruning strategies thattarget bonds that either have a maximal or minimal valueof ∆ B i or ∆ G i .The quantities M i and ∆ M i are related via a linear-response relation [8]: ∆ M i = M i /S i , (1) where kS i is local modulus characterizing the cost of achange in the equilibrium length of bond i .To employ Eq. 1 to tune G/B , one needs to com-pute B i , G i and S i for each bond. Evaluation of S i foreach bond requires N b calculations, however this is doneonly once for the initial unpruned network. Thereafter,the evolution of the spring network is efficiently com-puted using methods described in Appendix A. We notethat there are ( d ( d +1)2 − independent shear moduli in d dimensions, which are denoted as G ( j ) . Pruning basedon their average, denoted by G = d ( d +1)2 − (cid:80) j G ( j ) , al-lows tuning of G/B in a manner that leaves the systemisotropic. We note that ∆ G i is a linear average over all ∆ G ( j ) i , which will be important for determining its dis-tribution.In Fig. 1 the evolution of G/B is shown for differ-ent pruning strategies in which we remove bonds with: max ∆ B i , min ∆ B i , max ∆ G i or min ∆ G i . Results in di-mensions d = 2 and d = 3 are shown. The number ofremoved bonds is characterized by ∆ Z = ∆ Z − N r /N where ∆ Z is the initial excess coordination number and N r are the number of bonds removed. The pruningprocedures based on min ∆ B i and max ∆ G i result in avery small G/B ratio while those based on max ∆ B i and min ∆ G i result in a very large G/B . We emphasize thatonly a few percent of bonds are removed, yet the changein
G/B can be almost ten orders of magnitude. -3 -2 -1 ∆ Z -8 -6 -4 -2 G / B ( a ) 2 D max ∆ B i min ∆ B i max ∆ G i min ∆ G i -3 -2 -1 ∆ Z ( b ) 3 D Figure 1. The evolution of
G/B as a function of ∆ Z whichsignifies the number of bonds removed (a) in two dimensions.(b) in three dimensions To characterize the variation of
G/B versus ∆ Z , we fitthe values of G/B for each pruning strategy to a power-law ∆ Z µ . As illustrated by the larger exponents shownin Table I, pruning strategies based on ∆ B i and ∆ G i aremore effective in changing G/B than those based on B i and G i . Moreover, pruning the min ∆ G i bonds yieldsauxetic ( ν < ) behavior, which had not been achievedusing minG i . pruning method µ DB + µ DB − µ DG + µ DG − µ DB + µ DB − µ DG + µ DG − B i , G i − .
36 1 .
27 3 . - − .
96 1 .
01 1 . - ∆ B i , ∆ G i − . . . − . − . . . − . Table I. A comparison of the exponents, defined by GB ∝ ∆ Z µ ,when bonds are selected based on ∆ B i and ∆ G i values versuswhen they are selected by their B i and G i values, taken fromRef. [3]. The subscript of µ designates the targeted moduluswhile + ( − ) marks the maximal (minimal) values targeted. Equation 1 explains why the procedure based on re-moving the bond with minG i was unsuccessful. Becauseof the appearance of S i in the denominator, a bond witha small G i does not in general have a small ∆ G i . Ref-erence [8] shows that G i ∝ S i at small S i . Thus bondsthat seem unimportant and carry little stress may in factbe important; their removal can vary G significantly.To highlight the difference between ∆ G i and G i wenow consider their distributions. Ref. [3] shows that thedistribution of G i at small values scales as a power-law,which to a good approximation in three dimensions isgiven by, G ≈− . i . This suggests that there are manybonds which can be removed with little change to theshear modulus. We now argue that this is not the case,as will be inferred from the ∆ G i distribution.To compute the distribution of ∆ G i we employ theanalysis of Ref. [8] which studied the distribution of ∆ G ( j ) i for any shear direction. It was shown numerically,and supported by theoretical arguments that to a goodapproximation: P ∆ G ( j ) i (cid:68) ∆ G ( j ) i (cid:69) = y = 1 √ π y − e − y . (2)The distribution of ∆ G i is then given by the sum over n = d ( d +1)2 − different shear directions, ∆ G ( j ) i . Assum-ing that these are independent, the distribution of ∆ G i iscomputed in VI B and found to be a Gamma distribution: P (cid:18) ∆ G i (cid:104) ∆ G i (cid:105) = y (cid:19) = (cid:16) n (cid:17) n/ (cid:0) n (cid:1) y n − e − ny , (3)where Γ ( n ) is the Gamma function. The important ob-servation is that at small values P (∆ G i ) ∝ ∆ G d ( d +1)4 − i ,so that in three dimensions P (∆ G i ) ∝ ∆ G + i in contrastto the P ( G i ) ∝ G ≈− . i which has the opposite sign inthe exponent. Thus, most bonds thought to be unim-portant based on their G i actually lead to a substantialdecrease in G . -5 -4 -3 -2 -1 -8 -7 -6 -5 -4 -3 -2 -1 P ( ∆ B i › ∆ B i fi ) 2 DP ( ∆ B i › ∆ B i fi ) 3 DP ( ∆ G i › ∆ G i fi ) 2 D P ( ∆ G i › ∆ G i fi ) 3 DP ( ∆ G ( j ) i D ∆ G ( j ) i E ) 2 DP ( ∆ G ( j ) i D ∆ G ( j ) i E ) 3 D Figure 2. The distribution of ∆ B i , ∆ G i and ∆ G ( j ) i in bothtwo and three dimensions. The dashed curves overlaying P (cid:16) ∆ G ( j ) i (cid:17) is the prediction from Ref. [8] : √ π y − e − y . Thedashed curves overlaying P (∆ G i ) in two and three dimensionare the prediction based on Eq. 3. III. ESTIMATING THE EXPONENTS
To understand why the curves of
G/B versus ∆ Z looklike approximate power laws for the four different pruningstrategies, we approximate Eq, 1 as a differential equa-tion: dMd ∆ Z ≈ − α M ∆ Z , (4)where d (∆ Z ) = N is the change in the coordinationnumber when a bond is pruned and α ≡ ∆ Z N τ i S i . (5)Here τ i = M i /M is proportional to the energy on bond i ,but normalized so that (cid:80) i τ i = 1 . Typically, the stressesdue to a global deformation are not localized and there-fore τ i ∼ N . We also note that [8, 15] S i ∝ ∆ Z . Forthese two reasons, α should not depend on the systemsize or ∆ Z . If α is constant then the solution to thisequation is M ∝ ∆ Z − α . If bonds are chosen with somespecific rule and the distribution of α is stationary, thenits average remains constant as the system is pruned. Inthis case, M ∝ ∆ Z − α , where α denotes the average of α and depends on the pruning procedure.In Ref. [8] the distributions of ∆ B i and ∆ G ( j ) i weremeasured for different pruning procedures and it wasshown that for different pruning strategies the distri-bution of ∆ G ( j ) i is universal and is given by Eq. 2.While the starting distribution of ∆ B i is initially differ-ent, it evolves to this universal distribution for all prun-ing strategies discussed here except for the case when thebonds with min ∆ B i are targeted. Based on the defini-tion of α in Eq. 5 α ∝ ∆ M i , and therefore α will havethe same distribution P ( α ) = 1 √ πα α − / exp (cid:18) − α α (cid:19) . (6)The only free parameter is α and it can be evaluated bynoting that if bonds are removed randomly, α = 1 [3].By requiring that ´ dαP ( α ) α = 1 we find that α = 1 .We begin by considering the exponents associated withpruning the min ∆ M i . Numerically, it is found that inall cases the change M is very small. Due to the weakcorrelation between ∆ B i and ∆ G i , pruning the min ∆ B i ( min ∆ G i ) results in the decrease in G ( B ) as if a ran-dom bond is removed. Therefore, G ∝ ∆ Z ≈ . when the min ∆ B i bonds are removed and B ∝ ∆ Z ≈ . when the min ∆ G i bonds are removed.The case of pruning max ∆ M i is very different. In gen-eral, these bonds carry a lot of stress which is then redis-tributed upon their removal. Ref. [14] shows that thesestresses are redistributed on the length scale of ζ whichdiverges in the limit of ∆ Z → . Therefore, if the systemsize is smaller than ζ , then bonds are removed approxi-mately homogeneously throughout the system, while forsystems larger than ζ a system spanning crack forms. Allour data and analysis concerns the first, homogeneouscase.We begin by estimating the exponent associated withthe change of the bulk modulus B ∝ ∆ Z α max ∆ Bi for the max ∆ B i procedure. This requires a calculation of the av-erage maximal value of α , which depends on the numberof independent bonds. Assuming N independent randomvariables, the distribution of α max is given by: Q ( α max ) = N (cid:18) ˆ α max dyP ( y ) (cid:19) N − P ( α max ) . (7)Its average can be estimated numerically, yielding α max ∆ B i ≈ . for N = 4000 for comparison to sim-ulations. This is somewhat greater than the measuredvalue of . in two dimensions and . in three dimen-sions.A similar estimate can be found for the case of pruning max ∆ G i . In two dimensions the distribution of α , is thengiven by P ( α ) = exp ( − α ) , where the average is chosento be unity. The distribution Q ( α max ) can be computedanalytically Q ( α max ) = N (cid:0) − e − α max (cid:1) N − e − α max (8) ≈ N e Ne − αmax e − α max (9)where the asymptotic form is the Gumbel distribu-tion, with α max ≈ logN + γ , and γ ≈ . is theEuler–Mascheroni constant. This an instance of ex-treme value theory which predicts that when N is large, Q ( α max ) will always be given by a Gumbel distributionas long as P ( α ) decays fast enough [16]. In d the shear modulus should therefore scale as G ∝ ∆ Z α max ∆ Gi withthe exponent α max ∆ G i ≈ . (for N = 4000 as used insimulations). This is the same order but larger thanthe value ≈ . measured numerically. A similar analysisin three dimensions yields α max ∆ G i ≈ . which is also greater than the value of . found numerically.Since in all cases the estimated exponents are largerthan those found numerically we test the validity of ourassumptions. First, we consider the distribution of α andcheck if it is indeed stationary, focusing on the large α values crucial for the max ∆ M i pruning. To probe thetail of P ( α ) we measure α max as a function of bondsremoved. This is found to be a constant, after severalbonds are removed.Interestingly, the α max is found to be smaller than thatpredicted from Eq. 7. We note that α max has a logarith-mic dependence on the system size, N , since it samplesthe exponential tail of P ( α ) . Indeed, we find numeri-cally that the α max ∆ B i grows slightly when the systemsize is increased. However, contrary to our assumptionnot all bonds are independent, and as demonstrated inSection IV there are significant spatial correlations in thesystem. We believe these correlation reduce α max , sincethey re-normalize N to a smaller value. Naively, if thespatial correlations are the dominant contribution thenthe independent number of bonds should scale as N/ξ d f ξ where ξ is the correlation length associated with ∆ B i or ∆ G i and f ξ are the number of independent bonds within ξ d .We conclude by noting that despite the system sizedependence there are features that are universal, stem-ming from the α distribution being independent of thepruning protocol. This is demonstrated by comparingpruning max ∆ B i to max ∆ G (1) i , where G (1) is the sim-ple shear modulus. Since these share the same P ( α ) , weexpect that B (∆ Z ) should have the same behavior as G (1) (∆ Z ) after several bonds are removed. Indeed, inFig. 3 these two curves shown to be almost parallel inthe small ∆ Z regime. -2 -1 ∆ Z -5 -4 -3 -2 -1 BB G (1) G (1)0 Figure 3. A comparison between B (∆ Z ) for the max ∆ B i pruning to the G (1) (∆ Z ) for the max ∆ G (1) i pruning. Afterseveral bonds are removed these are almost parallel, suggest-ing universality. Here G (1) denotes the simple shear modulusand B and G (1)0 denote the bulk and shear modulus beforeany pruning. IV. SPATIAL CORRELATIONS
In this section we discuss spatial correlations of ∆ B i and ∆ G i . This is most easily done in Fourier space, forwhich we define: (cid:68) | ∆ M ( q ) | (cid:69) = 1 N σ M (cid:42)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) i ∆ M i e − iqr (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:43) (10)where σ M is the variance of ∆ M i which is a conve-nient normalization. The average is performed over alldirections in q space and over about 10 realizations of dis-ordered unpruned networks. Aside from normalization,this is the Fourier transform of the correlation functionof (cid:104) ∆ M ( r ) ∆ M (0) (cid:105) − (cid:104) ∆ M (cid:105) . Since identifying growinglength scales requires systems with a large linear dimen-sion, we focus on two-dimensional systems.Figure 4 shows (cid:68) | ∆ B ( q ) | (cid:69) and (cid:68) | ∆ G ( q ) | (cid:69) in twodimensions for different values of ∆ Z . There are tworegimes. At large values of q both quantities vary as q − γ with γ ∆ B ≈ . for the bulk modulus and γ ∆ G ≈ . for the shear modulus. This indicates spatial power-lawcorrelations which scale as, r − d + γ for r smaller than acorrelation length, ξ . At small q there is a crossover toa constant which signals a transition to an uncorrelatedstate. The crossover value of q is identified as the inversecorrelation length and its ∆ Z dependence is found by adata collapse where the two axes are scaled with a powerof ∆ Z . We find that the ∆ B i correlation length is givenby ξ ∆ B ∝ ∆ Z − . ± . while the ∆ G i correlation lengthscales as ξ ∆ G ∝ ∆ Z − . ± . . This differs from the usualpicture [17], in which the bulk modulus is associated with the cutting length (cid:96) ∗ ∝ ∆ Z − and the shear moduluswith the transverse length scale (cid:96) † ∝ ∆ Z − . -2 -1 q › | ∆ B ( q ) | fi (a) -2 -1 q › | ∆ G ( q ) | fi (b) -2 -1 q ∆ Z − . -1 › | ∆ B ( q ) | fi ∆ Z . q − . (c) ∆ Z = ∆ Z = ∆ Z = ∆ Z = ∆ Z = -2 -1 q ∆ Z − . -2 -1 › | ∆ G ( q ) | fi ∆ Z . q − . (d) Figure 4. The spatial correlations in two dimensions measuredby the angle average of ∆ M ( q ) defined in Eq. 10. In panel ( a ) and ( b ) we show the unscaled distributions while in panel ( c ) and ( d ) we collapse by rescaling the axis with powers of ∆ Z . The ∆ B i correlation length is consistent with ξ ∆ B ∝ ∆ Z − . while the ∆ G i correlation length is consistent with ξ ∆ G ∝ ∆ Z − . . Since ∆ M i depends on S i we also measure the cor-relations of S i . Using the same definition in Eq. 10we compute (cid:68)(cid:12)(cid:12) ∆ S i ( q ) (cid:12)(cid:12) (cid:69) by replacing ∆ M i with S i .These correlations are shown in Fig. 5 and the diverg-ing length scale is consistent with ∆ Z − . ± . . Interest-ingly, at large q , (cid:68)(cid:12)(cid:12) ∆ S i ( q ) (cid:12)(cid:12) (cid:69) ∝ q ≈− . as in the case of (cid:68) | ∆ B ( q ) | (cid:69) . This suggests that these correlations havethe same source. We also note the length scale ∆ Z − . has been observed in the tension profile resulting fromsqueezing a bond [18, 19]. This should be the same lengthscale. -2 -1 q ∆ Z − . -1 › | S i ( q ) | fi ∆ Z . q − . ∆ Z = ∆ Z = ∆ Z = ∆ Z = Figure 5. The spatial correlations in two dimensions of S i measured using (cid:68)(cid:12)(cid:12) ∆ S i ( q ) (cid:12)(cid:12) (cid:69) . The correlation length agreeswith ∆ Z − . V. CONCLUSIONS
In this paper we have shown that pruning bonds withmaximum or minimum values of ∆ B i and ∆ G i pro-vides an effective method for designing disordered meta-materials with targeted properties. Previous protocols [3]relied on targeting the contributions to the moduli, B i and G i . The current approach targets the change in themoduli when a bond i is removed (i.e., targeting ∆ B i and ∆ G i ). This procedure allows us to tune the network toa nearly maximally negative Poisson’s ratio by targetingbonds with very small ∆ G i . This was not possible withthe previous approach and is a significant improvementbecause it allows the network to be pruned to the auxeticlimit without developing cracks.We have also shown that if the system is kept isotropic,the possibility of removing a bond which contributes verylittle to the shear modulus is greatly reduced. This prob-lem is more severe in higher dimension, yet even in threedimensions there are enough small ∆ G i bonds to yieldan auxetic material for min ∆ G i pruning.We have provided a rationale for why moduli tend toscale as power laws with ∆ Z when bonds are removed.We estimated the exponents for different pruning strate-gies, and found agreement to within 30% of the measuredvalues. Our analysis suggests that the exponents havesome universal features, however, they depend weakly onthe system size. We also argue that spatial correlationsmay reduce the exponents.Finally we have examined the spatial correlations of ∆ B i and ∆ G i in two dimensions and identified diverginglength scales. The correlation length for ∆ B i is givenby ξ ∆ B ∝ ∆ Z ≈− . while the correlation length for ∆ G i is given by ξ ∆ G ∝ ∆ Z ≈− . . This contrasts with the intuition [17] that associates the (cid:96) ∗ = ∆ Z − with thebulk modulus and (cid:96) † = ∆ Z − with the shear modulus.We also find that S i has a correlation length of ∆ Z − . .These spatial correlation reduce the range of possible de-signs in the system. ACKNOWLEDGMENTS
We thank C. P. Goodrich for important discussions.We acknowledge support from the Simons Foundation forthe collaboration “Cracking the Glass Problem” award
VI. APPENDIXA. An efficient algorithm for recomputing theelastic properties resulting from bond removal
In this section we provide an efficient algorithm forrecomputing the elastic response due to the removal ofa bond. We employ the notation of Ref. [8, 20]. Theenergy cost of a deformation is given by: U = 12 k (cid:88) i,j e i [ S i ] j e j and the tension in a bond is given by: t i = k N (cid:88) j =1 e j [ S i ] j . (11)Here e i is the affine extension of a bond and [ S i ] j is thetension in bond j resulting from a unit change in theequilibrium length of bond i . All information regardingthe elastic behavior depends only on [ S i ] j and thereforewe would like to compute how it changes with the re-moval of a bond. In particular S i ≡ [ S i ] i and M i = t i V (cid:15) where V is the volume and (cid:15) is the strain. Eq. 5 of Ref.[8] allows to compute the change in the tension when abond is removed for any deformation. The correspondingaffine extension for this deformation is given by e j = δ ij .Assuming bond k is removed the element [ S i ] (cid:48) j is thengiven by: [ S i ] (cid:48) j = [ S i ] j − [ S i ] k [ S k ] j [ S k ] k . Thus the evolution of the whole elastic response is easilycomputed when a bond is removed, including B i , G i and S i . B. Distribution of the sum of random variableswith the universal form
We compute the sum of n independent ran-dom variables, Z = (cid:80) ni =1 y i , where P ( y i ) = √ πy y − / i exp (cid:16) − y i y (cid:17) . We exploit, the fact that x i = √ y i is a Gaussian random variable. P ( x i ) = 1 √ πy exp (cid:18) − x i y (cid:19) . (12)To compute the distribution of Z , we first compute thedistribution of R = √ Z = (cid:112)(cid:80) x i and then use a trans-formation of variables to compute P ( Z ) . The distribu-tion of R is the sum of Gaussian variables and thereforestraightforward: P ( R ) = S ( n )(2 πy ) n/ R n − exp (cid:18) − R y (cid:19) (13) where S ( n ) = π n Γ ( n ) is the surface of a n-dimensionalhypersphere and Γ is the Gamma function. A transfor-mation of variables results in: P ( Z ) = P ( R ) dRdZ = 12 S ( n )(2 πy ) n/ Z n/ − exp − Z y . (14) = 1(2 y ) n/ Γ (cid:0) n (cid:1) Z n/ − exp − Z y . (15)The average of Z can be computed by noting that it isa sum of n identical random variables, and thus is givenby n (cid:104) y i (cid:105) = ny . [1] S. R. Nagel, Rev. Mod. Phys. , 025002 (2017).[2] P. M. Reis, H. M. Jaeger, and M. van Hecke, ExtremeMechanics Letters , 25 (2015).[3] C. P. Goodrich, A. J. Liu, and S. R. Nagel, Phys. Rev.Lett. , 225501 (2015).[4] J. W. Rocks, N. Pashine, I. Bischofberger, C. P.Goodrich, A. J. Liu, and S. R. Nagel, Proceedings ofthe National Academy of Sciences , 2520 (2017).[5] L. Yan, R. Ravasio, C. Brito, and M. Wyart, ArXiv e-prints (2017), arXiv:1708.01820 [physics.bio-ph].[6] H. Flechsig, Biophysical Journal , 558 (2017).[7] L. Yan, R. Ravasio, C. Brito, and M. Wyart, Proceedingsof the National Academy of Sciences , 2526 (2017).[8] D. Hexner, A. J. Liu, and S. R. Nagel, ArXiv e-prints(2017), arXiv:1706.06153 [cond-mat.soft].[9] C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel,Phys. Rev. E , 011306 (2003).[10] A. J. Liu and S. R. Nagel, Annu. Rev. Condens. MatterPhys. , 347 (2010).[11] D. J. Durian, Phys. Rev. Lett. , 4780 (1995). [12] C. P. Goodrich, A. J. Liu, and S. R. Nagel, Phys. Rev.Lett. , 095704 (2012).[13] C. P. Goodrich, S. Dagois-Bohy, B. P. Tighe, M. vanHecke, A. J. Liu, and S. R. Nagel, Phys. Rev. E ,022138 (2014).[14] M. M. Driscoll, B. Gin-ge Chen, T. H. Beuman, S. Ul-rich, S. R. Nagel, and V. Vitelli, ArXiv e-prints (2015),arXiv:1501.04227.[15] M. Wyart, Ann. Phys. Fr. , 1 (2006).[16] S. I. Resnick, Extreme values, regular variation and pointprocesses (Springer, New York, 2008).[17] L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev.Lett. , 098301 (2005).[18] E. Lerner, E. DeGiuli, G. Düring, and M. Wyart, SoftMatter , 5085 (2014).[19] D. M. Sussman, C. P. Goodrich, and A. J. Liu, SoftMatter (2016), 10.1039/C6SM00094K.[20] T. C. Lubensky, C. L. Kane, X. Mao, A. Souslov, andK. Sun, Reports on Progress in Physics78