Table Of Content

Effective Medium Theory of Filamentous Triangular Lattice Xiaoming Mao,1 Olaf Stenull,1 and T. C. Lubensky1 1Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA (Dated: November9, 2011) We present a new effective medium theory that can include bending as well as stretching forces, and we use it to calculate mechanical response of a diluted filamentous triangular lattice. In this lattice, bondsarecentral-force springs, andtherearebendingforces between neighboringbondson thesamefilament. Weinvestigatethedilutedlatticeinwhicheachbondispresentwithaprobability p. Wefindarigiditythresholdindependentoffilamentbendingrigidityatpb ≃0.56andacrossover 1 characterizing bending-, stretching-, and bend-stretch coupled elastic regimes controlled by the 1 central-force isostatic point at pCF =2/3 of thelattice when bendingforces vanish. 0 2 PACSnumbers: 87.16.Ka,61.43.-j,62.20.de,05.70.Jk v o I. INTRODUCTION N 7 Random elastic networks provide attractive and real- ] istic models for the mechanical properties of materials t f as diverse as randomly packed spheres [1–3], network o glasses [4–8], and biopolymer gels [9–20]. In their sim- s plest form, these networks consist of nodes connected . at by central-force (CF) springs to on average z neighbors. m Theybecomemorerigidaszincreases,andtheytypically exhibit a CF rigidity percolation transition [21–23] from - d floppy,disconnectedclusterstoasamplespanning-cluster n endowed with nonvanishing shear and bulk moduli at a o thresholdveryclosetotheMaxwellisostaticlimit[24,25] FIG.1: (Color online) Phase diagram of thediluted filamen- c z =2d, where d is the spatialdimension, at which the tous triangular lattice showing the central-force and bending [ CF number of constraints imposed by the springs equals the rigidity thresholds respectively at p = pCF and p = pb, the 1 number of degrees of freedom of isolated nodes. Gen- bending-dominated regime at small κ in the vicinity of pb, v the crossover bend-stretch regime near pCF, and the stretch- eralized versions of these networks, appropriate for the 1 ing dominated regime at large κ. description of network glasses [4, 5] and biopolymer gels 5 [13–15],includebendingforcesfavoringaparticularangle 7 1 between bonds (springs) incident on a givennode. For a reduces to that considered by others [29, 30] and suc- . given value of z, networks with bending forces are more 1 cessfully predicts a second-order CF rigidity threshold rigid than their CF-only counterparts, and they exhibit 1 at z = 4 < z (p = 2/3) with µ increas- CF max CF m 1 a rigidity transition at z =zb <zCF. ing linearly in p p near p and approaching the CF CF 1 − Though numerical calculations, including the pebble undiluted triangular-lattice value of µ at p = 1. When : v game [23, 26], have provided much of our knowledge bending forces are introduced, our EMT predicts a new i about the properties of random elastic networks, effec- second-order bending dominated rigidity threshold at X tivemediumtheories(EMTs)[27–33]haveprovidedsim- p = p = 0.56 for all κ > 0 near which µ κ(p p ) b m b ar pleandqualitativelycorrectdescriptionsofCFnetworks. for κ/(µa2) 1 and µm µ(p pb) for κ/∼(µa2)− 1, ≪ ∼ − ≫ EMTs[34–36]andheuristicapproaches[37]thatdescribe where a is the lattice spacing. Near p , κ is a relevant CF both bending and stretching forces have only recently variable moving the system away from the CF rigidity been developed. Here we present details of the deriva- criticalpointtoabroadcrossoverregime[35,37]inwhich tion of a bend-stretch EMT introduced in Ref. [35] and µ κ1/2µ1/2 as shown in the phase diagram of Fig. 1. m ∼ itsapplicationtoabonddilutedtriangularlattice,whose This crossover is analogous to that for the macroscopic maximum coordinationnumber is z =6. This lattice conductivity in a resistor network in which bonds are max has bending and stretching forces modeled on those of occupied with resistors with conductance σ with prob- > biopolymer networks of filamentous semi-flexible poly- abilitypandwithconductanceσ <σ withprobability < > mers, characterized by one-dimensional stretching and 1 p [38]. − bending moduli µ and κ [13–15], respectively. Our EMT Though the model we study has both stretching and calculates the effective medium moduli µ and κ as a bending forces, it differs in important ways from previ- m m function of µ and κ and the probability p = z/6 that a ously studied models for network glasses [4–8] and for bond is occupied. Both the EMT bulk and shear mod- filamentous gels [13–20]. The maximum coordination ulus are proportional to µ . When κ = 0, our EMT number for both of these systems is less than or equal m 2 2d, and thus neither has a CF rigidity transition for The outline of our paper is as follows. Section II re- p < 1 when there are no bending forces. As a result viewspropertiesofsemi-flexiblepolymersanddefinesour neither exhibits the bend-stretch crossover region near model for the harmonic elasticity of crosslinked semi- p that our model exhibits. Network glasses are well flexible polymers on a triangular lattice; Sec. III sets CF modeled by a randomlydiluted four-foldcoordinateddi- up our effective medium theory; Sec. IV presents the re- amond lattice in which there is a bending-energy cost, sults of this theory; and Sec. V compares our EMT with characterized by a bending modulus κ, if the angle be- other versions of bend-stretch EMTs summarize our re- tween any pair of bonds incident on a site deviates from sults. There are three appendices: App. B and App. C the tetrahedral angle of 109.5deg. The architecture of presents details of the derivation of crossover functions the undiluted diamond lattice (with z =4<2d=6) in the vicinity of p and p , respectively; and App. D max CF b is such that its shear modulus vanishes linearly with κ providesdetailsontheestimationofp usinggeneralized b [7] and elastic response is nonaffine. When diluted, it Maxwell counting arguments. exhibits a second-order rigidity transition from a state with bending-dominated nonaffine shear response to a state with no rigidity. As dilution decreases, rigidity is II. FILAMENTOUS POLYMERS ON A stillcontrolledbyκ,butresponsebecomeslessnonaffine. TRIANGULAR LATTICE Filamentous networks in two-dimensions are often described by the Mikado model [14–16] in which semi- Theelasticenergyofasinglepolymeriswelldescribed flexible filaments of a given length L are deposited with by the worm-like chain model [14, 44, 45] of an elastic random center-of-mass position and random orientation rod onatwo-dimensionalplaneandthepointswheretwofila- mentscrossarejoinedinfrictionlesscrosslinks. Asinour 1 L dl(s) 2 dˆt(s) 2 E = ds µ +κ , (1) model, thereis no energycostfor the relativerotationof 2 ds ds Z0 h (cid:16) (cid:17) (cid:12) (cid:12) i two rods about a crosslink, but there is an energy cost (cid:12) (cid:12) for bending the rods at crosslinks. This model is charac- where s is the arclengthcoordinate, (cid:12)L is th(cid:12)e unstretched terized by the ratio η L/l of the filament length L to contour length of the polymer, dl(s)/ds andˆt(s) are, re- c the average mesh size,≡i.e., the average crosslink separa- spectively,the extensionalstrainandunittangentto the tionl >aalongafilament,whereaisaminimumcutoff polymer at s, and µ and κ are, respectively, the one- c length. In the limit η , all filaments traverse the dimensional stretching and bending stiffnesses. Three sample, and the system→ha∞s finite, κ-independent shear lengthscales canbe identified in this elastic energy. The and bulk moduli: There is effectively a CF rigidity tran- first is the contour length of the polymers, L. The sec- sition at z = 4 when η is decreased from infinity. There ond,lbend κ/µ,characterizestherelativestrengthof ≡ is a transition at η 5.9 to a rigid state with nonaffine stretchingandbending. Foranelastic rodmade ofa ho- response [14, 39], a≈nd there is a wide crossover region mogeneous mpaterial, lbend simply corresponds to the ra- between η = η and η = in which the shear modu- dius ofthe rod. Thethirdlength, the persistencelength, c lus changes from being ben∞d dominated, nonaffine, and lp κ/(kBT),describesthe scaleofthermalundulations ≡ nearlyindependentofµatsmallη tobeingstretchdomi- of the polymer at finite temperature. A fourth length, nated,nearlyaffine,andnearlyindependentofκatlarge the mesh size lc characterizing the connectivity of the η. Our EMT applied to the kagome lattice [40], whose network, can be identified for crosslinked polymer net- maximum coordination number like that of the Mikado works. model is four, captures these crossovers. Interestingly, Alltheselengthsparticipateindeterminingtheelastic- 3d lattices composed of straight filaments with z =4 ity of the polymer network above the rigidity threshold. max exhibit similar behavior [41]. When filaments are bent, Consider the force constants for harmonic distortions of however,elasticresponseinonecaseatleast[20]ismore anindividualpolymersegmentoflengthlc. Atzerotem- like that of the diluted diamond lattice with the shear perature, the force constant for longitudinal stretch of 4m.odulus vanishingwithκ evenatatlargeL/lc orz near tflheectsioengmisenkt is kκk,/µl3∼(igµn/olcr,inagndnutmhaerticfaolrftarcatnosrvse).rsTehdues- the ratio l ⊥ ∼/l c k /k characterizes the relative External tensile stress (i.e., negative pressure) can strengthofbebnedndcin∼g and⊥strekt,cµhingofthis particularseg- cause a floppy lattice to become rigid [42]. Random in- p ment of length l at low temperature. At finite tem- ternal stresses can do so as well in a phenomenon called c perature,thermalundulations ofthe filamentsinduce an tensegrity [25]. Thus a lattice with internal stresses may effective series spring constant k k Tl2/l4 that re- havealowerrigiditythresholdthanthesamelatticewith sults from pulling out thermal ukn,tdhu∼latiBonsp[13c, 15, 17]. outinternalstresses[20]. Systemssuchasnetworkglasses At any given temperature, the physical parallel inverse canexhibittworigiditytransitions[8,43]: asecond-order transition from a floppy to a rigid but unstressed state spring constantis simply k−1 =k−,µ1 +k−,t1h so that k ,th followed closely by a first order transition to a rigid but dominatesathightemperaktureankdk ,µkatlow tempekra- stressedstate. Theseeffectsarebeyondthescopeofyour ture. The ratio lc/lp k /kth charackterizes the relative EMT and will not be treated. ease of stretching vs∼. be⊥ndinkg of this segment at high 3 temperature. Finally, L/l measures the connectivity of c the network. 0 Here we investigate the harmonic elasticity of a model filamentous network on a bond-diluted triangular lattice withlatticeconstantawhichisthemeshlengthl . Poly- c mers correspond to lines of connected, occupied colinear bonds and crosslinks to sites at which two polymers to 3 2 crosslinks. We use a model elastic energy obtained by 1 discretizing the contrinuum model of Eq. (1) on our tri- (a) angular lattice. To simplify notation, we use µ/a µ/l c ≡ to represent the full k rather than k only, and we set ,µ l l'' k = κ/a3 κ/l3. Tkhus the elastickenergy for a defor- m⊥ation that≡mapsc the position of lattice sites ℓ from r l' Θ ℓ to R = r +u can be written into the stretching part ℓ ℓ ℓ and the bending part (b) FIG.2: (coloronline)(a)Filamentoustriangularlatticewith E = Es+Eb bondsrandomlyoccupiedwithprobabilityp. Theunitvectors E = 1µ g (u ˆr )2, (2a) eˆk1,eˆk2,eˆk3aremarkedby“1,2,3”,the3stretchenergyvectors s 2a ℓ,ℓ′ ℓℓ′ · ℓℓ′ Bsn are marked by the 3 red single lines, and the 3 bending hXℓ,ℓ′i energy vectors Bbn are marked by the 3 green double lines. E = 1κ g g θ2 The purple dashed double line marks the bending vector Bs4 b 2a ℓ,ℓ′ ℓ′,ℓ′′ ℓℓ′ℓ′′ if the origin is marked by 0. (b) Illustration of the bending hℓ,Xℓ′,ℓ′′i energyterm 2b. Thegreen curvescorrespond tothefilament 1 κ 2 by requiring that it has to pass through lattice sites ℓ,ℓ′,ℓ′′. = g g (u u ) ˆr (,2b) 2a3 ℓ,ℓ′ ℓ′,ℓ′′ ℓℓ′ − ℓ′ℓ′′ × ℓℓ′ The ratio of the bending elastic modulus of the darker green hℓ,Xℓ′,ℓ′′i (cid:2) (cid:3) segment and the lighter green segment is κ1/κ2 = 3. The bendingenergytermatsiteℓ′isproportionalto2(cid:0)κ11+κ12(cid:1)−1 determined by minimizing theelastic energy. where u u u , ˆr is the unit vector pointing ℓℓ′ ℓ′ ℓ ℓℓ′ ≡ − from site ℓ to nearest-neighbor (NN) site ℓ in the ref- ′ erence state (i.e., the state with u = 0 for all ℓ), g ℓ ℓ,ℓ′ III. EFFECTIVE MEDIUM THEORY is equal to one if the NN bond ℓ,ℓ is occupied and ′ h i zero if the bond is empty, θ is the angle between ℓℓ′ℓ′′ We study the elasticity of our network using an bond ℓ,ℓ and ℓ,ℓ inthedeformedstateasshownin ′ ′ ′′ h i h i effective-medium approximation [27, 28] that maps the Fig. 2(b), with r , r , and r along a line. The stretch ℓ ℓ′ ℓ′′ random network to an effective-medium uniform lattice springconstantµ/aisdefinedoneachNN bond,andthe in which all bonds are present with respective stretching bending constant κ/a3 is defined for each next-nearest- andbending constantµ andκ , whichare determined neighbor (NNN) bond along a straight line (i.e., in the m m self-consistently. This approximation has been shown samefilament). We haveusedthefactthatthe filaments to be a powerful tool for the calculation of properties are straight (θ =0) in the reference state so ˆr =ˆr . ℓℓ′ ℓ′ℓ′′ of random systems, from the electronic structure of al- Thereisanimportantdistinction,whichwillbeofimpor- loys[27,46]totheelasticityofrandomnetworks[29,47]. tance in our development of an effective-medium theory, It is straightforward to write the elastic energy ∆E of betweenthestretchingenergyE andthebendingenergy s this uniform effective medium lattice in terms of the dy- E . The former is a sum of energies on individual bonds b namical matrix, D : ℓ,ℓ connectingNN sitesℓandℓ,whereasthelatter is q ′ ′ h i a sum of energies on the pair of connected bonds ℓ,ℓ′ 1 and hℓ′,ℓ′′i or, alternatively, on a phantom NNN hbondi ∆E = 2N2 u−q·Dq·uq, (3) connecting sites ℓ and ℓ . This phantom bond is not q ′′ X an independent connection; it is only present and resists where the Fourier transform of u is defined as bending if both bonds ℓ,ℓ and ℓ,ℓ are occupied. ′ ′ ′′ Wdisecuwsisllthhaevreelmatoiroentoofsoauhyrawboioruktttohtihshaint oSifeRc.eVfs.w[3h4e]nawnde uq = uℓe−iq·rℓ, uℓ = N1 uqeiq·rℓ, (4) ℓ q [36]. X X andthedynamicalmatrixcanbewrittenasasumofthe Upon dilution, each of the bonds is present with a outerproductofthestretchingandbendingbondvectors probability p, and the resulting lattice corresponds to a random network of semiflexible filaments of finite ran- 3 3 µ κ dom lengths L, whose average as a function of p is D (µ ,κ )= m Bs Bs + m Bb Bb (5) L =1/(1 p) [35]. q m m a n,q n,−q a3 n,q n,−q h i − nX=1 nX=1 4 where phantom bonds is κ (κ ) κ , where c s m − Bsn,q = (1−e−iq·eˆkn)eˆkn κc =2 1 + 1 −1, (8) κ κ Bbn,q = 2(1−cos(q·eˆkn))eˆ⊥n (6) (cid:16) s m(cid:17) where κ = 0 if the bond ℓ ℓ is absent and κ = κ if s 1 2 s h i witheˆkn theunitvectoralongbondndefinedinFig.2(a) it is present. In the latter case κc =2κκm/(κ+κm) is a nonlinear function of κ and κ . and eˆ⊥n the unit vector perpendicular to eˆkn. The sum With these rules,wecanwrmitethe perturbationto the is over the three bonds in one unit cell for Bs or the dynamical matrix resulting from the replacement of one three pairs of bonds for Bb. We used the periodicity of bond as the lattice to obtain this translationally invariant form acaslcwuillaltbioencobmeleowcl.ear when we discuss the perturbation Vq,q′(µs,κs) = a−1(µs−µm)Bs1,qBs1,−q′ In the long-wavelength limit, the energy of our undi- +a−3(κc(κs)−κm)Bb1,qBb1,−q′ luted harmonic triangular lattice reduces to the elastic +a 3(κ (κ ) κ )Bb Bb , (9) energy of an isotropic 2d medium, − c s − m 4,q 4,−q′ where we have chosen to replace the bond between the E = d2x λˆ Tru 2+µˆTru2 , (7) two sites r0 = 0 and r1 = aeˆk1 for convenience. The Z "2 # vector Bb4,q ≡ e−iq·eˆk1Bb1,q represents the bending of the (cid:0) (cid:1) pond pair connecting sites r0 =0, r1 =aeˆk1, and 2r1, as whereuisthelinearizedsymmetricCauchystraintensor shown in Fig. 2(a). The replaced stretching spring con- with Cartesian components uij, λˆ and µˆ are the Lamé stant µs and bending modulus κs follow the probability coefficients, λˆ = µˆ = (√3/4)(µ/a), which depend only distribution onµandnotonκ. Thebendingconstantκonlyappears P(µ ) = pδ(µ µ)+(1 p)δ(µ ) in the O(q4) terms in this elastic energy. s s− − s The effective medium is characterized by a stretch P(κs) = pδ(κs κ)+(1 p)δ(κs). (10) − − spring constant µ that has the same value on every m Inthisstudyweareprimarilyconcernedwiththestatic NN bond and a bending constantκ that has the same m propertiesoftheeffectivemedium,whichischaracterized value on every NNN phantom bond along a line. The by the zero-frequency phonon Green’s function self-consistency conditions determining µ and κ are m m treated differently. In the dilution procedure, bonds are G = D 1. (11) removed with probability 1 p. In the process, a bond q − −q − with spring constant µ is replaced by one with spring m Whenonebondoftheeffectivemediumisreplaced,D constant 0, and the scattering potential relative to the → D+V,andtheGreen’sfunctionisgivenbythefollowing effective-medium state arising from this replacement is expansion proportionaltoµ µ ifthebondisoccupiedand µ if m m − − mituisstvabceantrte.aBteednddiinffgerceonutplyle.sCNoNnsNidesrittewsoalboonngdaslainloenagnda GVq,q′ = (I−G·V)−q,1q′ ·Gq′ = Nδ G +G T G , (12) line, ℓℓ′ and ℓ′ℓ′′ , sharing a common site ℓ′ and com- q,q′ q q· q,q′ · q′ h i h i posed of rods with respective one-dimensional bending where the T matrix is moduliκ andκ definedinEq.(1). Thephantombend- 1 2 ing bond connecting NNN sites is now a composite one 1 T V + V G V [Fig. 2(b)] consisting of two NN bonds tied together at q,q′ ≡ q,q′ N q,q1· q1· q1,q′ the intermediate site ℓ′. A direct calculation of the min- Xq1 1 imum bending energy [Eq. (1) with dl/ds= 0], for fixed + V G V G V bending angle θ, in this composite bond yields an effec- N2 q,q1· q1· q1,q2· q2· q2,q′ tivebendingconstantκeff =2(κ−11+κ−21)−1 thatreduces + q.X1,q2 (13) toκ1 ifκ1 =κ2 and0ifeitherκ1 orκ2 =0. Anexample ··· ofthe energyminimizing bending configurationis shown InEMT,theeffective-mediumelasticparametersµ and m in Fig. 2(b). We derive the effective bending constant of κ aredeterminedbyrequiringthatthedisorderaverage m this composite bond in App. A. of GV be equal to G of the effective-medium uniform Consider sites ℓ , ℓ , ℓ , and ℓ along a line. The re- 0 1 2 3 latticewithstretchingandbendingconstantsµ andκ . m m movalofthe centralbond ℓ ℓ willmodify the bending 1 2 This is equivalent to requiring that the disorderaverage, h i modulus of that bond and thus the bending energy of bthoethNNNNbNonpdhaℓntℓom.bTonhduss,hℓth0ℓe2ibeannddinhℓg1ℓm3iodcuonlutsairneilnag- pT|µs=µ,κs=κ+(1−p)T|µs=0,κs=0 =0, (14) 1 2 h i tive to that of the effective medium of both these NNN of the T-matrix vanish. 5 IntraditionalEMTwithnobendingmodulus,theonly containscrosstermssuchasBb Bb G Bb Bb elastic parameter is the stretching spring constant µ , which do not reduce to V f1.,q 1,−q1· q1· 4,q1 4,−q′ m q,q′ and one can simplify the T matrix by making use of the equation Here we develop a new variation of the EMT to solve this self-consistency equation (14). In this method, we 1 N Vq,q1·Gq1·Vq1,q′ =−a−1(µs−µm)Vq,q′f(µm)(,15) ibnotrriondgubcoenadnpeaffirescatilvoengmfieldaimumentcsotuhpaltinsghabreetawebeonndn,eeig.gh.-, Xq1 BbBb, with a coupling constant λ , to account for the 1 4 m where f is a scalar function. However, in this case with crossterms discussed above. For convenience we express bending energy terms, which couple the replaced bond the self-consistency equation in the space defined by the withtwoneighboringbonds,theproductV G V basis Bs , Bb , Bs , in which the perturbation is q,q1· q1· q1,q′ {| 1i | 1i | 4i} (µ µ )/a 0 0 s m V˜(µ ,κ )= −0 (κ (κ ) κ )/a3 λ /a3 , (16) s s c s m m  0 λ−/a3 (κ (κ−) κ )/a3 m c s m − −   where the tilde marks the addition of the λ coupling. with m Cross terms such as Bs1Bb1 vanish because they violate d2q the space inversion symmetry, q q. H Bs HBs Bs H Bs (23) →− 1 ≡h 1| | 1i≡ 4π2/v 1,q· q· 1,q On the other hand, the λm coupling terms also need Z1BZ to be added to the dynamical matrix where we have defined the inner product as an integral overthefirstBrillouinzoneandthefactorvisthevolume +Dãq(µ3λm,κm3)2=cDosq(q(µmeˆ,κ)Bmb) Bb , (17) aponfliednagcHshb3uent≡witehceBenlb1lB|iHnsr|aBenab4dlis.BpabActese.rwmSeismmvialaennriltsyih,oHnine2dt≡heeahdrBlyib1en|raH,m|cBiocub1ai-l − m · kn n,q n,−q matrix D˜, because D˜ is even in momentum q. As a re- n=1 X sult,theprojections Bs HBb =0and Bs HBb =0. and the Green’s function is modified to G˜ = D˜ 1. This is consistent withh t1h|e|van1iishing croshs t1e|rm|be4tiween − By writing the T matrix as − Bs and Bb in V˜−1. Thus the self-consistency equation (20) closes, providing 3 independent equations for T˜ =(V˜ 1 G˜) 1, (18) the3unknowns(µ ,κ ,λ )whosesolutionwewilldis- − − m m m − cuss in the next section. the self-consistency equation takes the form p[V˜−1(µ,κ) G˜]−1+(1 p)[V˜−1(0,0) G˜]−1 =0(.19) IV. EMT RESULTS − − − Multiplying by (V˜−1(µ,κ) G˜) on the left and by Numerical solutions to Eq. (20) for any value of κ/µ − (V˜ 1(0,0) G˜) on the right we arrive at are easily calculated, and the results for the effective − − medium elastic parameters, and the results are plotted pV˜ 1(0,0)+(1 p)V˜ 1(µ,κ) G˜ =0. (20) in Figs. 3, 4. There are several properties of these plots − − − − that are worthy of note: This equation can be solved in the space of 1. µ vanishes at the CF Maxwell rigidity threshold Bs , Bb , Bb . To do this we define the notation m {| 1i | 1i | 4i} pCF = 2/3 when κ = 0 and at p = pb = 0.56 for µ all κ > 0. Simulations of the same model yield H(b ,l ) mG˜(µ ,κ ,λ ). (21) m m m m m that p 0.659 and p 0.445 [36]. (In App. D, ≡− a CF b ≃ ≃ we present a variation of the Maxwell floppy mode Fromthe definition ofG˜ it is straightforwardto see that counttoestimatetherigiditythresholdinpresence λH /oµnly.dWepeetnhdesnopnrotjheectraHtiotos bthmis≡baκsims/µm and lm ≡ ionf figolaomdeangtrbeeemndeinntgwsittihffnseimssualnatdioonbtraeisnulptbs.)≃0.448 m m 2. Increasingκincreasesµ forallpexceptatp=p . H 0 0 m b 1 H(b ,l )= 0 H H , (22) 3. For small κ/(µa2), there is an interesting and m m 2 3   0 H H nontrivial crossover near p whereas for large 3 2 CF   6 54.. κltbuλ(κopirmm/eosah(n=ttµvaivaavvao1anani2n)ondl)i,usr,ivswheµmanhehlmµeesueiascmeafrrhosaaoκprptpirCtypiw→Fdpmo>i.lbtyfuhppstobarthbeunae.abtdncyCdhaiFnecnrsoyistnvteheshasrvtelersiudusmehuecoontnniolcdoedteianhliru.olsytfIieetttcffdsoreeosxiclatsahtsstitiouvbstveraiicalettyesr-- Κm0111.000000---.11975 ìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòô crossoverbehaviornear p=pCF for small κ/(µa2). 0.6 0.7 0.8 0.9 1.0 In what follows, we discuss analytic expressions derived p from Eq. (20) that captures these behaviors. (a) Μm011.00000--.1175 ìçàáòôìçæàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçìàáçæòôáàôòìæçàáòôìæçàáòô1111111ìæç0000000àáòô-------ìæçàáòô0123456ìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòô Λm111100000.--0--1114185 ìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòô0ìæçàáòôìæçàáòô.6ìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòô0ìæçàáòôìæçàáòô.7ìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòô0ìæçàáòôìæçàáòô.8ìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòô0ìæçàáòôìæçàáòô.9ìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòô1ìæçàáòôæìçàáòô.0 0.6 0.7 0.8 0.9 1.0 p p (b) (a) Μm000001......024680 ìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòôìæçàáòô (oFbfI)Gκ..Pa4as:ra(cmoeloterrosnµalnimnde=)coTlaoκh2reκcEo/Md(µeTcaa22rδs)eopl+tuhtecio1sn,amfoeraκsmin(aF)iga.nd3(.2λ4m) 0.6 0.7 0.8 0.9 1.0 p where c and c are constants determined by the ar- 1 2 (b) chitecture of the lattice and are independent of p or κ/(µa2), and both κ and λ are linearly proportional m m to µ for small δp with corrections of order δp2. The FIG. 3: (color online) (a) Semi-log plot for the EMT m solution µm as a function of p for µ = 1 and κ = linear dependence on δp of these quantities is what one 1,10−1,10−2,10−3,10−4,10−5,10−6 from top to bottom, as would expect from and EMT. The prefactor, shows im- indicatedinthelegend. (b)Linearplotofµm asafunctionof portant dependence on κ. When κ/(µa2) 1, µm is ≪ p,withparametersandcolorcodethesameasin(a). Thered linearly proportional to κ and independent of µ whereas solidlineindicatesµmforthecaseofacentralforcetriangular for κ/(µa2) 1, µ is linearly proportional to µ and m lattice (κ=0). independent≫of κ. Thus the system behaves as if the it consistsofa “bending spring”with springconstantκ/a3 in series with a stretching spring with spring constant µ/a. This result explains the increase in slope of µ at m A. Behavior near p=pb p with increasing κ and the saturationof µ for p near b m p for large κ/(µa2). In Fig. 5 we show a comparison b Webeginwithbehaviornearp . AsdetailedinApp.C, between the asymptotic solution (24) and the numerical b µ approaches zero linearly in δp p p for all values solutions of the EMT equations, which agree very well. m b ≡ − 7 Μm111000---432çáòô çáòô ìçáòô ìçáòô ìçáòô ìçáòô ìçáòôìçáòôìçáòôìçáòôìçáòôìçáòô afdetohanerorsMdlcsfy|reuapircnb−eercxeloptadpsiCtorceFbenods|ypss∝tiecetodhx|eµehpilLnmarLesi,tst−aesiamcrehsimLd´tewyisicCenoooFdffet|iffi.estprhc.Tcmueihsesesuneffstod,sef,ctihhtnwLievEhiesiqccama.hrl(eei7nad)trgi.heueeTmxbshpoacuotmahsnn,eeltnibhnateess- 10-5ì ì à à à à àààààà raesgyimmpetsootifctfhoermdisluotfeµdmlacttaincec:haracterizedifferent elastic 1100--76à æà æ æ æ æ ææææææ 3 4 κ µ = µ∆p 1+ 1 A (28) .001 .002 .005 .01 m | |2 ± s − 9 a2µ∆pφ! | | ∆p √ µ1/2κ1/2 if κ 1, aA a2µ∆pφ ≫ 3µ∆p if |κ | 1 and ∆p>0,. FnuIGm.e5ri:ca(Cl soololurtoionnlisne()daAtsaypmopintotst)icosfolµumtionnea(drapsbh.edPalirnaems)etaenrds ≃ 3Aa2| ∆κp| if aa22µµ|κ∆∆pp|φφ ≪≪1 and ∆p<0. | | | | and color code are the same as in Fig. 3. These crossover regimes correspond exactly to those found in Ref. [37] using known behavior of the density B. Behavior near p=pCF ofstatesandmode structureofsystems nearthe CFiso- static limit and general scaling arguments. Thus, it is straightforward to see that above the WeshowinApp.BthatinthevicinityoftheCFrigid- crossoverregime,themacroscopicelasticmoduliarepro- itythresholdatsmallκ,theeffectivemediumspringcon- portional to µ, indicating that the elasticity of the net- stants for given ∆p p p adopt the scaling forms ≡ − CF workis dominatedby stretching energies,whereas below κ the crossover regime, the moduli are proportional to κ, µ = µ∆pt1g , (25a) m | | 1,± a2µ∆pφ indicating that the network elasticity is dominated by (cid:16) |κ | (cid:17) bending energies, in agreement with the fact that the κ = µ∆pt2g , (25b) m | | 2,± a2µ∆pφ systemisbelowthecentralforceisostaticpoint. Interest- (cid:16) |κ | (cid:17) ingly, in the crossover regime near pCF, the macroscopic λm = µ|∆p|t3g3,± a2µ∆pφ , (25c) elasticmoduliis proportionaltoµ1/2κ1/2,corresponding (cid:16) | | (cid:17) to strong coupling between both stretching and bending where modes. These observationsindicate a phasediagramasshown φ = 2, (26a) inFig.1. Thereis arigiditytransitionbetweenthe state withnorigidity(vanishingelasticmoduli)andrigidstate t = 1, (26b) 1 atp . Abovethisrigiditytransitionandatsmallfilament b t = 2, (26c) 2 bending stiffness κ, the elasticity of the diluted filamen- t3 = 3, (26d) touslatticeiscontrolledbyacrossoveratthecentralforce isostaticpointp . Forlargebendingstiffness(κ µa2) CF and thesystemdirectlyentersthestretchingdominated∼affine elastic regime as p is increased. 3 4 g (x) 1+ 1 Ax , (27a) 1, ± ≃ 2 ± r − 9 ! V. DISCUSSION 1 g (x) x, (27b) 2, ± ≃ 3 1 Two other approaches produce results similar to ours, g (x) B 1+ 1 4Ax − x2, (27c) and we briefly review compare them to ours. 3, ± ≃ 27 ± r − 9 ! In references [34] and [36], Das et al. present an alter- native EMT for systems with both stretching and bend- where 2.413 and 1.520 are constants evalu- ing forces. The first of these reference fails to identify A ≃ B ≃ ated in App. B. These scaling relations are analogous to either the new rigidity threshold at p and the wide b thatfoundinrandomresistornetworkswithtwodifferent bend-stretch crossover regime near p , but the sec- CF types of resistors [38], and central force spring networks ond does. Stretching forces are easily described by CF with strong and weak springs [37]. springs, which reside on bonds, each of which can have Plotsofthenumericalresultsinthesescalingformsare a distinct spring constant, whereas bending forces effec- shown in Fig. 6 along with the above analytical result. tively reside on phantom NNN bonds. Because remov- These scaling relations can also be expressed in terms of ingoneNN bondfromapairdefining aphantomNNN themeanfilamentlength L =1/(1 p)thus L =3 bending bond effectively removes that phantom bond, CF h i − h i 8 bending and stretching are not independent in the di- 104 luted lattice. This presents real challenges for the devel- opmentofaconsistentbend-stretchEMT.Ourapproach ΜmtHL-pp1CF 01.00011 æææææææææææææææææææææààæàæààààæàààæààààæààæààààìàæììàìàìæìììàììàìììììàììàìæìììàòòìòìòàòòòìòòìæòòòòòìòòìòàòòòìòôôòôòôìôôôôòôôàòôôôôòôôòôìôôôôòôôòôìôôôòôôò ô ô bWsfttdooutoerrrenwectedteheccoidainhscfngwoompnmimurtsoorhtoodbadmdlunloieulftocmiylsdaulleetskalhsa|κp|tedirnpsfoeeotetnworacmhlhaiseinoiancgtddhcigovuehpaecltinonebhnnadoeNlsdnobtueNdnientn.gnutdbdeaeeMnoirnnfftnlyogyedpcdi,mbntoiifihgolvoyncoedmpawdutoC.eeilloivrFyDsneµmrias,ffaoetrmnrerfeerdontetehcnndκahet---. 10-4 ifies the bend force constant k for both phantom NNN 10-4 0.01 1 100 104 bondsthatthatNN bondpart⊥iallydefinesinthemanner Κ described above. With this approach, we develop a con- Hp-p LΦ CF sistent EMT that includes the statistical correlation be- (a) tweenbendandstretch. Dasetal. ignorethiscorrelation andassumethatastretchspringonagivenNN bondcan ô be removed without affecting the bending energy on the 1000 ò ô phantom NNN bonds that include that NN bond and ΚmtHL-pp2CF0.0001.101æææææææææææææææææææææààæàæàààààæààæààààæààæàààààìæììàìàìæìììàììàììììàìììàìæìììàòòìòìòàòòòìòòìòæòòòòìòòìòàòòòìòôôòôòôìôôôôòôôàòôôôôòôôòôìôôôôòôôòôìôôôòôô TatwItbtnwhhreoibeatoghrtipetpebrrnebmhooaxeeavnrorinnyisadvddttlespeeo,idnnrmttaogchhwnbeaebsiatitapNonbhprdndioipNledneiurfitpgtNoyinsesxaenqiobeffmdtnaffoeehnaecnnettdctdtheitNnseaiovdgbfcNepeastlstheyhnNchnaerebetniuplbetswenottovpiirmdnoaterehntdetreNcs.lpdoehyrIfnNiontnsttobpNghwoareaNtibinhtbcirhgNelooeisrtnnawysbldowitnto1sborhna−orctddinaahnqssdnnet.., 0.001 0.1 10 1000 that the phantom bond does not exist unless both of Κ Hp-pCFLΦ the NN bonds defining it are present, Das et al. assign a probability q = p2 (p is the probability that a NN (b) bond is occupied) to the occupancy of a NNN bending bond, but continue to treat the NN and NNN bonds as statistically independent. Again the result is a set of ΛmtHL-pp3CF111100000--.01574 æææææææææææææææææææææààæàæààààæàààæààààæààæààààìàæììàìàìæìììàììàìììììàììàìæìììàòòìòìòàòòòìòòìæòòòòòìòòìòàòòòìòôôòôòôìôôôôòôôàòôôôôòôôòôìôôôôòôôòôìôôôòôôò ô ô cBsaaptielrooorvvontsvaaahesellcduudah[ee3psei5enspf(]lop,ryfsr-obiiwcemapol=hcdbnuhicslpe0(ahi0ssC.t4t.iyFy5eo)iin6nee=it)lslndde2[t3gqhap/5ouaC3n]oatFoiwdtnnii=shotganrewog0irsvor.ee6eifdlaaole5slrma1atbg.µbehreomnOaenvtetduemawo-rnstfeiatdhtnRrphaetκpettwmfsrc.(iohim.0t[a3hc.Bcu46rhl4soo]ai5tsyymths)iiiooeeuoavllnlddbpaesssr---. inthevicinityofp inqualitativeagreementwithsimu- CF 10-4 0.01 1 100 104 lations. Our approach gives values of the shear modulus inthecrossoverregionofthesameorderofmagnitudeas Κ Hp-pCFLΦ seen in simulations whereas that of Ref. [36] finds values 1.5 to 2 orders of magnitude higher at small values of κ. (c) In Ref. [37], Wyart et al. consider random off-lattice FIG. 6: (color online) (a), (b), and (c) show rescaled plots elastic networks derived from two-dimensional packings of the EMT solutions µm, κm and λm using the scaling of spheres [48] with a cordination number above the forms (25a) for µ = 1 and κ = 10−2,10−3,10−4,10−5,10−6, Maxwell CF isostatic limit of z =4 in which CF springs with color code the same as in Fig. 3, and exponents taking are assigned to each sphere-sphere contact. They use the value as in Eq. (26a). The thin black lines represent the numericalsimulationstostudythenonlinearrelationbe- asymptotic forms (27a) for small κ as solved in App.B. The tween shear stress σ and shear strain γ as springs are brown dash-dotted lines, the thick blue solid lines, and the purple dashed lines plot the functional form of µm obtained cut, thereby reducing z, and they find a scaling relation in the crossover, the stretching-dominated, and the bending- σ = γ δz f(γ/δz ), where δz = z 4, f(x) const. for dominated regimes of Eq. (28), respectively. x 0|+,|f(x) | 0| for x 0−, and−f(x) x→for x . → → → ∼ →∞ This scaling form predicts σ γ for δz γ > 0, σ = 0 for δz γ <0, and σ γ2∼for γ δz≫. Reference [37] ≪− ∼ ≫ then provides a theoreticaljustification for this behavior based upon the existence of a plateau in the density of 9 states [49] above ω δz and reasonable assumptions Therefore, we can assume the shape of the two bonds ∗ ∼ aboutstatisticalindependence ofeigenvectorsassociated 1 1 withdifferentnormalmodes inthe isostaticnetwork[50] h(x)=Ax+ B x2+ C x3 for x<0 1 1 andabout the nature ofnonaffine responseof nearlyiso- 2 6 1 1 staticsystems. Finally,theyextendthislineofreasoning h(x)=Ax+ B x2+ C x3 for x>0 (A4) 2 2 to nearly isostatic systems with extra weak bonds and 2 6 find three regimes of elastic response that are identical where the firstline correspondsto the bond between ℓℓ, ′ to those we identify in Eqs. (25a) to (27a) if the weak andthesecondlinecorrespondtothebondbetweenℓℓ . ′ ′′ bonds are of a bending type. This form already satisfies the condition that the two To summarize we developed a new effective-medium bonds are connected and smooth at their conjunction theory that can include bending energy of filaments, x=0. We have two boundary conditions at x= 1 and and we used it to study the development of rigidity x=1 − of a randomly diluted triangular lattice with central 1 1 force springs on occupied bonds and bending forces be- h = A+ B C 0 1 1 tween occupied bond pairs along a straight line. We ob- − 2 − 6 1 1 tained a rigidity threshold for positive bending stiffness h =A+ B + C . (A5) 0 2 2 and a crossover, controlled by the isostatic point of the 2 6 central force triangular lattice, characterizing bending- Under these boundary conditions we minimize the total dominated, stretching-dominated,and stretch-bend cou- elastic energy of the two bonds pled elastic regimes. Acknowledgments—We are grateful to C. P. Broedersz κ 0 d2h(x) 2 κ 1 d2h(x) 2 1 2 E = dx + dx and F.C. MacKintosh for many stimulating and helpful 2 dx2 2 dx2 discussions. This work was supported in part by NSF- Z−1 (cid:20) (cid:21) Z0 (cid:20) (cid:21) κ C2 DMR-0804900. = 1 B2 B C + 1 2 1 − 1 1 3 (cid:18) (cid:19) κ C2 + 2 B2+B C + 2 . (A6) 2 2 2 2 3 Appendix A: Bending of a composite rod (cid:18) (cid:19) This calculation leads to the minimum elastic energy In this Appendix we solve for the effective bending 6κ κ h2 rigidity of a composite rod. We solve this problem using E = 1 2 0. (A7) m a simple setup shown in Fig. 7. We assume that there is κ1+κ2 abondofbendingrigidityκ connectingpointsℓℓ,anda 1 ′ In the above calculation, if we assumed that both of bond of bending rigidity κ connecting points ℓℓ , with 2 ′ ′′ the two bonds are of bending stiffness κ , we get points ℓ,ℓ,ℓ hold in positions eff ′ ′′ E =3κ h2. (A8) m eff 0 1,h , 0,0 , 1,h , (A1) 0 0 {− } { } { } Thereforewearriveattheeffectivebendingstiffnessof the composite rod respectively in space, without losing generality. We then calculatetheshapeofthebondswhichminimizetheelas- 2κ κ 1 2 ticenergygiventhesefixedendpoints,treatingthebonds κeff = . (A9) κ +κ 1 2 aselasticrods. Inadditionweassumethatthestretching stiffness is muchgreaterthan the bending stiffness so we which reduces to κ if κ =κ as we expected. 1 1 2 can ignore the length change of the bonds. The elastic energy of a bond is then simply the bending energy of a rod with bending stiffness κ Appendix B: Asymptotic solution of the CPA self consistency equation near pCF κ d2h(x) 2 E = dx (A2) In this Appendix we solve the EMT self-consistency 2 dx2 Z (cid:20) (cid:21) equation [Eq. (20)] asymptotically near p at small κ. CF In the zeroth order, κ = 0, and the problem reduces to wherewehaveassumedsmalldeformationssotheelastic that of a central force rigidity percolation. We expand energytakestheformofanexpansioninsmallh(x). The theequationaroundthispointandsolvefortheeffective- variationalcondition that a configurationminimizes this medium spring constants asymptotically in small κ. elastic energy is Toobtaintheasymptoticsolution,werewritetheEMT self-consistency equation [Eq. (20)] as d4h(x) =0. (A3) dx4 pV˜(µ,κ)+(1 p)V˜(0,0) V˜(0,0)G˜V˜(µ,κ)=0,(B1) − − 10 κ=0. It is straightforwardto see fromthe symmetry ofV˜ 1 − andG˜ as in Eqs.(16) and(22) that the matrix equation [Eq. (20)] contains 3 independent equations. The 11- element gives FIG. 7: (Color online) Bending of the composite rod. The 3 red disks mark the fixed end points of the two bonds. The p H1(bm,lm) µ =µ − . (B2) blue solid line represent the shape of the two bonds solved m 1 H (b ,l ) 1 m m for the case of κ1 = 10κ2, and the red dashed line represent − the shape of the two bonds for the case of κ1 = κ2 which is symmetric underx→−x. To zeroth order, κ = 0, we obtain H1(0,0) = pCF = 2d/z = 2/3 from the fact that all bonds are equivalent, and Eq. (B2) reduces to the EMT equation for central- to avoid the apparent singularity in V˜ 1 at the point force rigidity percolation. − The 22- and 23-elements of the matrix equation (20) read 1 2 1 + 1 − p 1 1+ bm b H l H +(b2 +l2 )H +2b l H =0, (B3) b bm − 2 b − m 2− m 3 m m 2 m m 3 (cid:16) (cid:17) (cid:16) (cid:0) 1 (cid:1) (cid:17) l 2 1 + 1 − l H +b H +2b l H +(b2 +l2 )H =0, (B4) − m− b bm m 2 m 3 m m 2 m m 3 (cid:16) (cid:17) (cid:0) (cid:1) where b κ/(µ a2), H (b ,l ) and H (b ,l ) are the 22 and 23 elements of H(b ,l ) respectively. Thus we m 2 m m 3 m m m m ≡ have 3 unknowns µ ,b ,l (or equivalently, µ ,κ ,λ ) and 3 equations (B2), (B3), and (B4). As we already m m m m m m { } { } discussed, at κ = 0, the zeroth order solution is {µ(m0) = µp1−ppCCFF,b(m0) = 0,lm(0) = 0}, where pCF = H1(0,0). Notice that the equations up to this point are all exact within the E−MT. As κ becomes positive, µ increases,b and l become nonzero,andthe rigiditythreshold jumps to a lowervalue m m m p as shown in Fig. 3(b). For small κ, we have κ/(µa2) 1, and can assume that b ,l 1 (which we will verify b m m ≪ ≪ later), and we find that to the leading order the 3 Eqs. (B2), (B3), and (B4) become p p H (0,0)κ /(µ a2) CF 1,1 m m µ µ − − , (B5a) m ≃ 1 p CF − κ κ(2p 1), (B5b) m ≃ − κ 1 p λ κH (0,0) − (2p 1)2, (B5c) m ≃ 3 µ a2 p − m whereH (0,0)=∂H /∂b 2.413andH (0,0)=1.520. Fromtheserelations,wefindthatatp=p 1,1 1 m|bm=0,lm=0 ≃− 3 CF µ κ1/2, (B6a) m ∼ κ κ, (B6b) m ∼ λ κ3/2, (B6c) m ∼ indicating that µ κ λ and thus b ,l 1 as we assumed. Using these relations, together with the fact m m m m m ≫ ≫ ≪ that as κ 0, µ p p , we arrive at the scaling relations m CF → → − κ µ = µ∆pg , (B7a) m | | 1,± a2µ∆p2 (cid:16) | κ| (cid:17) κ = µa2 ∆p2g , (B7b) m | | 2,± a2µ∆p2 (cid:16) |κ | (cid:17) λ = µa2 ∆p3g , (B7c) m | | 3,± a2µ∆p2 (cid:16) | | (cid:17) where ∆p p p . These scaling forms agree well with our numerical solutions (See Fig. 6). CF ≡ −