Elasticity of randomly diluted honeycomb and diamond lattices with bending forces

0.83 MB·

by  Danilo B. Liarte

Elasticity of randomly diluted honeycomb and
diamond lattices with bending forces
6
Danilo B. Liarte
1
0 LaboratoryofAtomicandSolidStatePhysics,CornellUniversity,Ithaca, NY,
2 USA
InstituteofPhysics,UniversityofS˜aoPaulo,S˜aoPaulo,SP,Brazil
r
a
M
O. Stenull

1 DepartmentofPhysicsandAstronomy,UniversityofPennsylvania,
Philadelphia,PA,USA

]
n
n Xiaoming Mao
-
DepartmentofPhysics,UniversityofMichigan,AnnArbor,MI,USA
s
i
d
. T. C. Lubensky
t
a
DepartmentofPhysicsandAstronomy,UniversityofPennsylvania,
m
Philadelphia,PA,USA
-
d January2016
n
o
c Abstract. We use numerical simulations and an effective-medium theory to
[ study the rigiditypercolation transition of the honeycomb and diamond lattices
when weak bond-bending forces are included. We use a rotationally invariant

2 bond-bending potential, which, in contrast to the Keating potential, does not
v involve any stretching. As a result, the bulk modulus does not depend on the
7 bending stiffness κ. We obtain scalingfunctions forthe behavior of someelastic
2 moduli in the limits of small ∆P = 1−P, and small δP = P−Pc, where P
1 is an occupation probability of each bond, and Pc is the critical probability at
6 whichrigiditypercolationoccurs. Wefindgood quantitative agreement between
0 effective-medium theoryandsimulationsforbothlatticesforP closetoone.
.
1
0
6
1
: Submitted to: J. Phys.: Condens. Matter
v
i
X
r
a
2
1. Introduction Bending forces effectively couple next-nearest
neighbor sites, thereby increasing the effective z
Concepts associated with the rigidity percolation to values above the central-force critical value, z ,
c
transition of random elastic networks [1, 2] have providing elastic stability, and eliminating all but the
been applied in many branches of science, such as trivial zero modes of rigid translation and rotation.
amorphous solids [3, 4], granular materials [5, 6], In contrast to Keating potentials [23], bending forces
mineralogy[7],networksofsemi-flexiblepolymers[8,9] do not depend on the length of bonds and, therefore,
and the mechanics of living cells [10, 11, 12] . The do not contribute to the bulk modulus at = 1.
archetypeofthistransition[13,14,15,9,16,17]occurs The Keatingpotentialis rotationalinvariant. BPending
inperiodiclatticesinwhichbondsconsistingofcentral- forces are also, as we show in the Appendix.
force springs are populated with probability . For
P
below a threshold , the lattice loses rigidity, and (a)
c
P P
some or all of its elastic moduli vanish because floppy
regions prevent rigid ones from percolating.
The homogeneous honeycomb and diamond lat-
tices (Fig. 1) with only nearest-neighbor bonds are
strongly under-coordinated. The average coordination
number z of both the honeycomb and diamond lat-
tices is less than the Maxwell limit [18, 19] z =2d for
c
central forces, where d is the spatial dimension, below
which lattices under periodic boundary conditions de-
velop zero-frequency “floppy” modes. The honeycomb
lattice with 2 sites per unit cell and z = 3 has a defi-
ciencyofonebond,andthediamondlatticewith2sites
per unit cell and z = 4 has a deficiency of two bonds
(b)
per unit cellrelative to the Maxwelllimit. As a result,
they have an extensive number of zero modes and no
resistancetosheardistortions;butcuriouslybecauseof
the specialgeometryoftheir lattices,they bothhavea
non-zero bulk modulus B. Thus extra forces, such as
next-nearest-neighbor central forces or bending forces
favoringaparticularanglebetweenpairsofbondsshar-
ing common endpoints, are required for mechanical
stability. Here, using both effective medium theory
(EMT) and numericalsimulations, we study the prop-
erties of randomly bond-diluted central-force honey-
comb and diamond lattices with added bending forces
characterized by a local bending stiffness κ, focussing
in particular on behavior near zero dilution ( 1 )
P ≈
andneartherigiditythresholdat = . Wefindthat
c
P P
thebulkmodulusofbothlatticescanbeexpressednear
=1 as a scaling function of κ/(1 )n with n=1, Figure 1. a) Crystal structure of the honeycomb net. b)
P −P Conventional cubiccellofthediamondlattice.
muchliketheshearmodulusinthedilutedkagomelat-
tice with bending forces [16] where n = 2 rather than
1. The shear moduli of the honeycomb and diamond
2. Model
lattices, on the other hand exhibit, no simple scaling
form and approach zero even at = 1 as κ 0 as
P → The two-dimensional honeycomb lattice (Fig. 1a) is
required. We find that the bulk moduli, much like the
a triangular Bravais lattice with a two-point basis
shear moduli in the Mikado model [20, 21, 22] and the
[24]. It is defined by a set of primitive vectors, e.g.
dilutedtriangular[9]andkagomelatticeswithbending
a = √3(1,0), and a = (√3/2)(1,√3), along with
[16],exhibitcrossoverfromstretchingdominatedaffine 1 2
the positions of the atoms within each primitive cell,
responsetobendingdominatednonaffineresponsewith
c = (0,0), and c = (0,1). The length of each
decreasing , reaching a maximum nonaffinity at the 1 2
P bond in the lattice is ℓ = 1. The three-dimensional
rigidity threshold. The shear modulus, which vanishes 0
diamond lattice can be represented as a face-centered
withκ,ontheotherhand,alwaysexhibitsbendingand
cubic lattice with a two-point basis (Fig. 1b). It can
thus non-affine response.
3
be defined by the primitive vectors a =(1/2)(0,1,1),
1
a = (1/2)(1,0,1), a = (1/2)(1,1,0), along with (a)
2 3 (b)
the two-point basis vectors c = (0,0,0), and c =
1 2
(1/4)(1,1,1). Bonds connecting nearest-neighbor sites
are of length to ℓ =√3/4.
0
We consider the following lattice energy of
interaction:
E =E +E . (1)
stretch bend
The first term is a sum of central-force interactions
between nearest-neighbor pairs of sites i and j, which
in the harmonic limit are given by: Figure 2. a) Density plot in the qx×qy plane of one of the
honeycomb’s acoustic modes for k = 1, and κ = 0.01. b)
1
Eij = k[(u u ) rˆ ]2, (2) Dispersion curves along some symmetry lines for k = 1, κ = 0
stretch 2 j − i · ij (dashedcurves),andκ=0.1(solidcurves).
where rˆ is the unit vector connecting sites i and j
ij
in the undeformed lattice, and u is a displacement
i
vector. Each unit cell in the honeycomb lattice has
three independent bonds (e.g., the blue lines of Fig.
1a), and in the diamond lattice, four independent
bonds. The second term in Eq. (1) is a sum of
bendingenergyinteractionsassociatedwithtwobonds,
terminatingonacommonsitelandconnectingnearest-
neighbor pairs of sites (l,m), and (l,k):
κ
Eklm (sinβ ∆β )2, (3)
bend ≡ 2 0 klm
where β is the equilibrium angle between bonds, and
0
∆β represents the difference between the angle
klm
between the bonds (l,m) and (l,k) and β . There are Figure 3. Three-dimensional plot of the acoustic branch
0
frequenciesofthehoneycomb latticefork=1andκ=0.01.
six such terms per primitive cell for the honeycomb
lattice, and twelve for the diamond lattice. The
bending interactions effectively couple next-nearest- lattice has six phonon branches, of which two are
neighborssites,asillustratedbythereddashedlinesof floppy when κ 0. Figure 4a shows diamond-lattice
Fig.1a. Thefactorsinβ0 isamatterofconvenience. It dispersion curv→es for κ = 0.01, along symmetry lines
isequalto√3/2forthehoneycomblattice,and2√2/3 ΓX and ΓL. A sketch of the first Brillouin zone of
for the diamond lattice. Equation (3) may be written the diamond lattice with five high symmetry points is
asacombinationofdisplacementvectors,uptosecond displayed in Figure 4b. Notice that two largest and
order in u, as the two smallest eigenvalues are degenerate in both
Eklm = κ u [rˆ (rˆ rˆ )rˆ ] lines. These degeneracies can be broken along other
bend 2ℓ2 lm· lk− lk· lm lm less symmetrical lines.
0
n
2
+ulk [rˆlm (rˆlk rˆlm)rˆlk] , (4) a)
· − ·
where ulk =uk ul. o b)
−
The honeycomblattice hasfour phononbranches, L
Γ
corresponding to the four degrees of freedom within Ω K
X
each unit cell. In the limit κ 0, one of the two W
→
acoustic branches is floppy with zero frequency for
all wavenumbers q = (q ,q ) in the Brillouin Zone. L G X
x y
Figure 2a shows a density plot of one of the acoustic q
branch frequencies for homogeneous k = 1, and for
κ = κ˜ = 0.01. Hereafter k = 1 is assumed unless Figure 4. a) Dispersion curves of the diamond lattice along
symmetry lines L−Γ and Γ−X The dashed lines represent
otherwise noted. Figure 2b displays dispersion curves
κ=0andthefulllinesκ=0.01. b)SketchofthefirstBrillouin
along symmetry lines ΓM, ΓK, and KM, κ = 0
zoneofthediamondlatticeandhighsymmetrypointsΓ,L,X,
(dashed curves), and κ = 0.1 (solid curves). Figure K,andW.
3 show a 3d plot of the two acoustic branches as a
function of q and q , for κ = 0.01. The diamond
x y
4
3. Effective-medium theory and simulations analysis,asinglepairofbonds,sharingacommonsite,
contributesoneconstraintduetobendingforces. Each
We use an adaptation of the EMT developed in additionalbondsharingthatsiteisconstrainedtohave
Ref.[25,10]inwhichthespringconstantsofindividual a particular orientation,and thus adds a total of d 1
bonds and the bending constants of individual bond constraints, which is the number of angles needed−to
pairs are treated as independent random variables . specify a unit vector in d-dimensions. Therefore, the
‡
Theprobabilitythatagivenbondisoccupiedis ,and total number of bending-force constraints (for r 3)
the probability distribution for the spring constaPnt for reads, ≥
any bond is
nTh(r+1)=d 1+nTh(r) (9)
P (k′)= δ(k′ k)+(1 )δ(k′). (5) B − B ⇒
spr P − −P nTBh(r)=(d−1)r−(2d−3), (10)
Both bonds of a bond pair must be occupied in order
with
for there to be a bending energy associated with the
z z
pair. Followingthe approximationofRefs.[25, 10], we
r = in n , (11)
i i
set the probability that a given bond pair exists equal !, !
i=2 i=2
to 2. The probability distribution for the bending X X
conPstant of an individual bond pair is where ni is the number of i-coordinated sites. Note
that inapplyingEq.(10), one hasto ensurethat there
Pbend = 2δ(κ′ κ)+(1 2)δ(κ′). (6) are neither isolated nor 1-coordinated sites. Thorpe’s
P − −P
The joint probability for both k′ and κ′ is then counting, including r/2 stretching constraints, gives
P(k′,κ′)=P (k′)P (κ′). rc = 6(d 1)/(2d 1) = 2 and 12/5 = 2.4 for the
spr bend − −
honeycomb and diamond lattices respectively.
Inthis EMTtheory,eachoccupiedbondandeach
In effective-medium theory, the diluted lattice
occupied bond pair constitutes a constraint. Thus if
is modeled by a homogeneous lattice with stiffness
each site has z neighbors, there are (z /2)N bond
constraintsand(z(z 1) 2/2)N bond-paPirconstraints constants km and κm satisfying a set of self-consistent
− P equations that depend on a probability distribution
in a diluted lattices of N sites. Since each site has d
P(k′,κ′). We use an adaptation of the effective-
translationaldegrees offreedom, the Maxwellcountof
medium theories proposed in Ref. [13] to derive the
the number of zero modes per site is
set of effective-medium theory equations for k and
m
z z(z 1)
f =d P − 2,. (7) κm [25, 10, 16]:
− 2 − 2 P
k a∗
The EMT rigidity threshold is obtained by setting m = P − , (12)
k 1 a∗
f =0 to produce −
κ 2 b∗
m
= P − , (13)
1 8d(z 1) κ 1 b∗
= 1+ − 1 . (8)
c −
P 2(z 1)"r z − # where
−
aTvheursa,geincotohredihnoantieoynconmubmblaetrtiacte tPhcres≈ho0ld.6,0zan=d zthe, a∗ = z1 dv˜qTr Ds,q·Dq−1 , (14)
c P Z1BZ 0
is1.8;inthe diamondlattice, 0.56andz 2.24. (cid:0) (cid:1)
c c
P ≈ ≈ 1 dq
The EMT threshold should be compared with b∗ = Tr D D−1 , (15)
two other estimates of zc based on the Maxwell (z−1)z Z1BZ v˜0 b,q· q
zero-mode count. Phillips [27, 28], who considered where v˜ is the volume o(cid:0)f the first B(cid:1)rillouin zone. In
0
general off lattice networks, also treated bond pairs Eqs. (14) and (15), D is the translational invariant
q
as independent and obtained an equation identical to dynamical matrix in Fourier space:
Eq. (7) but without the factors and with z r
interpreted as the averagePcoordination number,→r = Dq,q′ =Nδq,q′Dq, (16)
z. His estimate leads to r = √2d or r = 2 ( =
c c c D =D +D , (17)
P P q s,q b,q
0.67)andr =2.45 ( =0.61)for the honeycomband
c
P where D , and D are contributions from the
diamond lattices, respectively. Thorpe later showed s,q b,q
stretchingandbending interactionsrespectively. They
[1] that treating each bond pair independently over-
may be written as
counts the number of bending constraints. In his
z
‡ ThisversionofCPAiswidelyusedandacceptedintreatments D =k Bs Bs , (18)
of rigidity percolation. We divide the interactions up into s,q m n,q n,−q
equivalency classes and average, so our effective Hamiltonian is nX=1
writableasasumofhomomorphicparts. ThisversionofCPAis where,
equivalent to what Yonezawa and Ogadaki call the HCPA [26],
whichyieldsananalyticalphysicalsolution. Bns,q = e−iq·fnen,−en , (19)
(cid:8) (cid:9)
5
and, 2 off-diagonal elements to be equal to γ, where γ is
κ the magnitude of the deformation which we set to
D = m Bb(i) Bb(i)
b,q ℓ02 1≤mX<n≤z(cid:16) mn,q mn,−q cγon=jug0a.0t1e. grTahdeienntwealrgeolraixthtmhethδauti pursoinvgideassutasndwaitrhd
+Bb(ii)Bb(ii) , (20) the equilibrium non-affine displacements δuna in the
mn,q mn,−q i
presence of applied deformation. Feeding these back
(cid:17)
where ℓ is the lattice spacing, and
0 into the elastic model energy density (1), we obtain
Bb(i) = e−iq·fme⊥ +e−iq·fne⊥ , the shear and bulk moduli as function of and κ. In
mn,q nm mn P
addition to the elastic moduli, we also compute the
(cid:8) −(e⊥mn+e⊥nm) , (21) so-called non-affinity parameter Γ which measures the
Bmb(nii,)q = −(e⊥mn+e⊥nm),e(cid:9)iq·fme⊥nm degree of non-affinity in the system under the applied
(cid:8) +eiq·fne⊥mn . (22) deformation,
1
The vectors en and fm (cid:9)connect nearest-neighbor sites Γ= Nγ2 (δunia)2 , (27)
and cells, respectively, , and
i
§ X
e⊥ =e (e e )e . (23) where N is the total number of sites. Our numerical
mn m− m· n n results will be displayed and discussed together with
These definitions imply thata∗ andb∗ arefunctions of
our EMT results as we move along.
the dimensionless ratio:
κ
m
κ˜ , (24)
m ≡ k ℓ2 4. Results
m 0
rather than of κm and km separately. Equations (14), In our EMT, the bulk and shear moduli of the
(15), and (17) lead to the following important relation honeycomb and diamond lattices have the same form
between a∗ and b∗: asfunctionsoftheeffectivemediumspringandbending
a∗+(z 1)b∗ =(2d)/z. (25) constants km and κm:
−
1
At the rigidity threshold, km = 0 and κm = 0. B =ABkm = (C11+2C12),
Equations (12) and (13) then require a∗ = and 3
b∗ = 2, and Eq. (25) reduces to Eq, (7) at Pf = 0 µ =A kmκm/ℓ20
and yiPelds Eq. (8) for c. µβkm+γκm/ℓ20
P
A
= µk κ˜ (1+(γ/β)κ˜ )−1 =C
3.1. Simulations β m m m 44
A κ
In the numerical portion of our work, we generate µ m, if κ˜ 1,
β ℓ2 m ≪
diluted honeycomb and diamond lattices on a = 0 (28)
 A
µ
computer. The systems sizes that we simulate range  γ km, if κ˜m ≫1,
up to 1002 unit cells for the honeycomb lattice and
203 unit cells for the diamond lattice. In all our whereCij are the standard Voigt elastic constants for
simulations, periodic boundary conditions areapplied. a cubic crystal, κ˜m = κm/(kmℓ20), AB ≡ AB,H = 3/4,
To facilitate the computations, we split up the elastic Aµ AµH = 27/2, β βH = 2, and γ γH = 9
≡ ≡ ≡
displacement ui into an affine and a non-affine part, for the honeycomb lattice, and AB ≡ AB,D = 1/12,
A A = 144, β β = 27, and γ γ = 192
u =ηx +δu , (26) µ ≡ µ,D ≡ D ≡ D
i i i for the diamondlattice. Thus B andµ aredetermined
where x is the equilibrium position of site i in as a functionofκ and once the EMTequations (12)
i
P
the absence of any applied deformation, η is the and (13) are solved.
deformation gradient tensor, and δu is the non-affine Figure 5 shows plots of numerical solutions of the
i
displacement. We fix the non-affine displacement of EMTequations(solidlines)andsimulations(symbols)
an arbitrarily chosen lattice site to be zero so that for the bulk (left) and shear (right) moduli of the
spurious zero modes associated with rigid translations honeycomb lattice, as a function of the probability ,
P
ofthe lattice aresuppressed. We apply shearandbulk for κ = 1, 10−2, 10−4, and 10−6 (in blue, red, yellow,
deformations by choosing η accordingly. For example, and green respectively). We will keep the same color
to apply shear to the honeycomb lattice, we chose definitions in all subsequent plots. Figure 6 shows
the 2 diagonal components of η to be zero and its similar plots for the diamond lattice. In both cases,
simulationsandtheEMTresultsagreewellnear =1.
§ Wehavechosene1 =−c2,e2=a2−c2,e3=a2−a1−c2, P
In the vicinity of the rigidity threshold = , the
f1=c1,f2=a2,andf3=a2−a1 forthehoneycomb lattice, P Pc
and e1 = a1−c2, e2 = a2−c2, e3 = a3 −c2, e4 = −c2, simulations display a decay that is different from that
f1=a1,f2=a2,f3=a3,f4=c1 forthediamondlattice.
6
PHLΚB, 0111.000000---.11975 ææææììòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàììòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòææàòôæà PHLΜΚ, 0111.000000---.11975 ææìæìæàìæàìæàìæàìæòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàììòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàòæà
0.5 0.6 0.7 0.8 0.9 1.0 0.5 0.6 0.7 0.8 0.9 1.0
P P
Figure 5. Bulk (on the left) and shear (on the right) moduli of the diluted honeycomb lattice as a function of P, for k =1, and
κ=1(bluecircles),10−2 (redsquares),10−4 (yellowdiamonds),and10−6 (greentriangles),frombothsimulations(symbols),and
effective-medium theory(solidlines).
1
PHLΚB, 1110000.0---1864 ææìææàìæàæìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàììòæìàòæìàòæìàòæìàòæìàòæìàòæìàòæìàòæìàòææàòôæà PHLΜΚ, 1110000.0---11864 æææææææìæàæìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàìòæàììòæìàòæìàòæìàòæìàòæìàòæìàòæìàòæìàòæìàòæàòæà
0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.4 0.5 0.6 0.7 0.8 0.9 1.0
P P
Figure6. Bulk(ontheleft)andshear(ontheright)moduliofthediluteddiamondlatticeasafunctionofP,fork=1,andκ=1
(bluecircles),10−2 (redsquares),10−4 (yellowdiamonds),and10−6 (greentriangles).
foundintheEMT.Wefind 0.5inthesimulations at small κ˜ where κ˜ is approximately equal to κ˜.
c m
P ≈
for the diamond lattice. We can use Eq. (11) and, Thus we expand a∗ and b∗ to linear order in κ˜ .
m
When κ˜ = 0 exactly, a∗ is ill defined because D
z z m s,q
z = in n , (29) has a nullspace [dimension 1 (2) for the honeycomb
i i
P !, ! (diamond) lattice)], and its inverse does not exist.
i=0 i=0
X X
When κ˜ is small, the projection of D onto the
along with a combinatorial calculation: m q
null space of D is proportional to κ˜−1, but the
z s,q m
n =N i(1 )z−i, (30) contribution of this projection to the trace is zero
i
i P −P
(cid:18) (cid:19) because by definitionDs,q is zeroinits ownnullspace.
to show that at P = Pc, r = rc ≈ 2.5, which is close Thus, a∗(κ˜m) has a well-defined limit as κ˜m → 0 and
to the value r = 2.4 that was obtained by He and a well-defined power series in κ˜ . To linear order
c m
Thorpe [14] using a different protocol. a∗ =1 ακ˜ + (κ˜2 ), (32)
The EMT equations provide analytic expressions − m O m
whereαis3.39and4.62fortheHoneycomblatticeand
for κ andk as a function ofκ˜ and inthe vicinity
m m P diamond lattices, respectively. A similar analysis can
of =1 and = , where
P P Pc be applied to b∗(κ˜m), but with Db,q rather than Ds,q
κ
κ˜ , (31) having a nullspace. It is easier, however, to use the
≡ kℓ20 relation, Eq. (25) between b∗ and a∗ to obtain
is a unitless measure of the bending stiffness. We 2d z α
b∗ = − + κ˜ + (κ˜2 ). (33)
begin with near one and seek expressions valid z(z 1) z 1 m O m
P − −
7
For small ∆ =1 and κ˜, these results imply inresponseto isotropiccompressionandthatresponse
P −P
km ∆ is nonaffine. Figure 7 plots B/B0 as a function of τ
=1 P forthehoneycomblatticeobtainedfromfullnumerical
k − 1 a∗
− solution of the EMT equations, the scaling function
∆
=1 P + (∆ κ˜m) of Eq. (40), and numerical simulations. The three
− ακ˜ O P
m curves agree extremely well near = 1. In addition,
∆ k ℓ 2 P
=1 P m 0 + (∆ κ˜ ,∆ 2), (34) the simulation and full EMT curves follow each other
m
− ακ O P P closely, and both break away from the scaling curve
whereinthe lastlinewehaveusedEqs. (13),(24)and at increasing values of ∆ as κ˜ decreases. The shear
P
(33). Thus, modulus µ cannot be expressed as a scaling function
−1 of τ. It approaches
k 1
m
k ≈ 1+ ακ˜∆P , (35) µ0 =(Aµ/β)k(1+(γ/β)κ˜)−1 (41)
(cid:18) (cid:19)
(36) as ∆ 0.
→
To obtain analytic expressions for the bulk and
and similarly,
shear moduli near the rigidity threshold at =
c
κ 2z(z 1) P P
m 1 − ∆ . (37) [given by Eq. (25)], we begin with the fact that km =
κ ≈ − z2 2d zακ˜ P κ =0 at that point. Thus
− − m
a∗(κ )= , b∗(κ )= 2, (42)
c Pc c Pc
where κ ( 0.91 for the honeycomb lattice and
c
≈ ≈
0.30 for the diamond lattice) is the value of κ˜ at
0.11 ìììììììììààààààààààìàààààààààààààààààæàææææææææææàææææææææææææææææææ eitnhqteuearretiisgotinedditfyoinrthκwrchesawhtoithlhdapoPpbcetnagsiinvteeodnlbobywyessEotlqvo.i(nr2dg5e)rt.hienWaδe∗m(κarc=e)
(cid:144)BB0 0.00.0011 òòòòòòòòììòìììììììììòìììììììì Path∗e−=cPoPeccffi−acscie(itnκ˜t)iδncPc(rκ˜ea)ansadessby∗fer=tomuPnc2zd(eeκ˜rto)e.r+mTci(onκ˜e)tdδh.PisT/e(hnzed−n,1ww)eePwusitsehet
1100--54 æòòòòòòòòòòò aEthqfesu.nli(cm1t2iiot)naδnodf δ(P13.)0.tTohBoisubttraabitneioctahmuesuersatotifaopEκpqmrso/.a(c(κ˜h4k2κm),cℓ/20bκ˜)otaihns
10-6 10-4 0.01 1 100 P →
the numerator and the denominator of this ratio are
Τ
proportional to δ and the ratio itself depends on c
P
but not δ in the limit δ 0. This limit equation
P P →
Figure 7. Scaling collapse for the bulk modulus of the determines c(κ˜):
honeycomblatticenearP =1. Thesolidlinesarethenumerical
κ 2(z 1) s κ˜
solutions to the complete EMT equations, the dashed line the c(κ˜)= c− − Pc c (43)
scaling solution of Eq. (40), and the data points are from κc+κ˜sc
simulations.
where
is a fκumnctisioanfoufntchtieonraotifo∆P and κ separately, but km sc = (z 11)−(1Pc 2) (44)
− −Pc
κ˜ The equations for k and κ then become
τ = (38) m m
∆
k 1 c(κ˜)
P m
only. We can then define a scaling function Ξ for the − δ , (45)
k ≈ 1 P
c
bulk modulus B near =1: −P
P and
B =B Ξ(τ), (39)
0 κ 2(z 1) +c(κ˜)
m c
− P δ . (46)
where, κ ≈ (z 1)(1 2) P
− −Pc
1 −1 Equations (45) and (46) predict linear behavior of
Ξ(τ)= 1+ (40)
the elastic moduli near = . However, scaling
ατ c
(cid:18) (cid:19) collapse plots of the simPulationPresults suggest that
and B0 = Aµk. When τ → ∞, [∆P ≪ κ/(kℓ20)] , B and µ scale with δ 1.6, for both honeycomb and
B approaches the undiluted, κ-independent limit of P
diamond lattices, as it is shown in Figure 8 for the
B = A k characterized by affine compression of all
0 B shear modulus. The exponent 1.6, which is near the
bonds. When τ 0 [∆ κ/(kℓ2)], B κ˜/(α∆ ),
→ P ≫ 0 → P value 1.5 found by He and Thorpe [14], is obtained
indicating that bonds in the diluted lattice are bent
as the best fit of our data near δ = 0. However,
P
8
4 ìò ìàò 1.0 ìò à
à ìò à
ìò
3 à 0.8 ìò à
ŽΚ ìàò ŽΚ ìò à
ŽHL(cid:144)Μ+Κ1.3 012 ìæàòìæàòìæàòìæàòìæàòìæàòìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò æ æ æ æ ŽHL(cid:144)Μ+Κ0.15 0000....0246 ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò ìæàò æà æ æ æ æ æ æ
0.00 0.05 0.10 0.15 0.20 0.00 0.05 0.10 0.15 0.20
∆P ∆P
Figure 8. Scaling collapse for the shear modulus for the honeycomb (left) and diamond (right) lattices near P = Pc, for k = 1,
and κ = 1 (blue circles), 10−2 (red squares), 10−4 (yellow diamonds), and 10−6 (green triangles). The red line is a best fit near
δP =0. Thesmallervalues atκ˜=1areconsistent withtheEMTexpressionforµ(Eq. (28)). Thesimulationspower-lawdecay is
characterized byanapproximateexponent of1.6forbothlattices.
we note that a log-log plot of this data is very noisy The Poisson ratio at threshold of the honeycomb
neartherigiditypercolationthreshold,andthatamore lattice,
precisedeterminationoftheexponentrequiresacareful
1 (µ/B)
analysis of finite-size scaling, which is beyond the goal σP = − = 0.233, (49)
1+(µ/B) −
of this manuscript. Finally, the simulation nonaffinity
is thus negative. The Poisson ratio of the diamond
ratio for the bulk modulus Γ , as a function of , is
B
P lattice, which continues to have macroscopic cubic
shown in Fig. 9.
symmetryintheEMTdependsondirectionofstresses,
and we do not calculate it.
1
PHLGΚ,B 01.000.000-.1114ìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàììòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàòôæàììòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôæàìòôàìòôàìòôàìòôàììòôàìòôàìòôàìòôàìòôàìòôàìòôàìòôàìòôàòôà 5Tpwainpr.oitotpehSprrsoaeuurxcnmmttieimimeaomsranaeatsrsoi.riotzfy-neTn,tehhiawegensheedbhholoaarnrtvoeteiyitncacetstoestimoruaandbrcaieteliladouynnndpsdihneiodrvnn-iaacoromtinahoonreandtndihndabalretaementldatdoi,sicntnieiingccs
æ òô the sense that their average coordination numbers are
æ
æ
10-5 æ ææ smaller than the isostatic limit zc = 2d. They present
0.5 0.6 0.7 0.8 0.9 1.0 oscillation modes which reduce to floppy modes over
P
the entire Brillouin zone when the bending energy
coupling constant vanishes. We implement disorder
byassigningaprobabilitydistributiontothe existence
Figure 9. Nonaffinity ratio of the bulk modulus ΓB for the
honeycomb latticeasafunctionoftheprobabilityP. of each bond in the lattices. When the number
of diluted bonds hits a critical value, there is a
In our effective medium theory, the ratio of µ to rigidity percolation phase transition at which both
B is a function of κ˜ : bulk and shear moduli vanish. We employ numerical
m
−1 simulations and an effective-medium theory to study
µ A β
= µ +γ . (47) scaling behavior near = 1, where the effective-
B A κ˜ P
B (cid:18) m (cid:19) medium theory reproduces the exact results for the
At the rigidity threshold at which κ˜m = κc, this ratio undiluted lattice, and near the critical probability
becomes = . The scaling behavior predicted by the
c
P P
µ 1.608 Honeycomb EMT near = 1 is fairly well satisfied by the
= (48) P
B 6.097 Diamond numerical simulations. However, the EMT predicts a
(cid:26)
different decay near = , with the elastic moduli
c
P P
9
proportional to δ , in contrast with the simulation where i and j are Cartesian indices x,y,z. The final
results, where B,µP δ 1.6. forms are the long-wavelengthcontinuum limits with
∼ P
1
u = (∂ u +∂ u +∂ u ∂ u) (A.5)
Acknowledgement ij 2 i j j i i · j
the usual rotationally invariant nonlinear Lagrangian
This work was supported by the Brazilian agencies
straintensorandu(x) the displacementfieldatspace-
Fapesp and Capes (DBL), and the NSF under Grants
point x. These limits were obtained using r
0
No. DMR-1104707 (TCL) and No. DMR-1120901 ≡
x as the reference point and u r ∂ u(x) and
a ai i
(OS). DBL thanks James Sethna for useful comments →
u r ∂ u(x). v and v are the forms that
b bi i a b
and suggestions. →
normally appear in central-force models and v is
ab
the quantity that Keaton introduces. The continuum
Appendix A. Derivation of the bending energy limits guarantee that a continuum elastic energy
in terms of displacements constructed from the bending energy of Eq. (A.4)
will be a function of u only, as it must be to be
ij
Thisappendixwillshowthatbendingpotentialscanbe
rotationally invariant. It should be noted, however,
castasamanifestlyrotationallyinvariantcombination
that the complete bending energies have terms higher
of displacement vectors that obey the Keating Rules
order in derivatives.
[23]. First some notation: Let r be the equilibrium
µ We expressthe quantities in Eq.(A.4)interms of
reference position of site µ, and let R = r +u be
µ µ µ v , v and v :
a b ab
the position of site µ after distortion, where u is the
µ R R 2 =R2R2 (R R )2 = r r 2+2V
displacement vector. Now, consider two bonds, which | a× b| a b − a· b | a× b| ab
welabelaandb,sharingacommonsite,whichwelabel Vab =(ra2vb+rb2va−2ra·rbvab)+2(vavb−va2b) (A.6)
as site 0. Bond a connects site 1 to site 0 and bond b
Theformofthisexpressiondependsonwhetherr and
a
connects site 2 to site 0. Then define
r are parallel or not: if they are parallel, r r = 0
b a b
×
ra =r1 r0; rb =r2 r0; (A.1) and β0 = 0,π. This limit applies to models for
− −
filamentous lattices [8, 20, 22]. The other limit β =0
R =R R r +u ; u =u u 0
a 0 1 a a a 1 0 6
− ≡ − applies to the honeycomb and diamond lattices. It
R =R R r +u ; u =u u .
b 2− 1 ≡ b b b 2− 0 is important to note that Eq. (A.6) and thus the
r is the equilibrium vector and R the stretched bending energy E depends on both the Keating part
a,b a,b b
vector for bond a (b). v and the “central-force” parts v and v , but by
ab a b
We seek an energy that depends on the angle β construction it depends only on the the rotation angle
ab
between bonds a and b. To define the β , we assume between bonds and not on any stretch in the bonds.
ab
that it lies between 0 and π so that its sine is positive The Keating energy depends on stretching as well as
(extensiontonegativeβ is possiblebut notrelevantto bending and contributes to both the bulk and shear
our current interest). Then moduli. The bulk modulus does not depend on a
pure bending energy of the type we discuss. There is
R R r r
a b a b
sinβab = | × |; sinβ0 = × , (A.2) no change in the angle between bonds under uniform
R R r r
a b a b compression. With the above definitions,
whereR = R andr = r ,andwesetthebending
a a a a
energy of the| bo|nd-pair ab|to| Ra Rb 1+2Vab/ra rb 2 1/2
| × | = | × | (A.7)
E = 1κ˜(β β )2, or (A.3) RaRbsinβ0 (cid:18)(1+2va/ra2)(1+2vb/rb2)(cid:19)
b 2 ab− 0 and
1 R R r r 2 V v v
= κ˜ sin−1 | a× b| sin−1 | a× b| . β β tanβ ab a b (A.8)
2 R R − r r − 0 ≈ 0 r r 2 − r2 − r2
(cid:18) a b a b (cid:19) (cid:18)| a× b| a b(cid:19)
It is clear that this energy is rotationally invariant. in the small-displacement limit.
We can now express these energies in terms of EquationsEq.(A.4)and(A.6)provideacomplete
1 1 expression for bending energies in terms of nonlinear
v = (R2 r2)= (2r u +u u ) r r u
a 2 a− a 2 a· a a· a → ai aj ij functions of the nonlinear “discrete” strain functions
1 1 va,vb,andvab. Weareofteninterestedintheharmonic
vb = 2(Rb2−rb2)= 2(2rb·ub+ub·ub)→rbirbjuij limitofthesefunctions. Webeginwiththecaseβ0 >0,
1 and we expand to lowest order in ua and ub
v = (R R r r )
ab 2 a· b− a· b 1
v r u , v (r u +r u ), (A.9)
1 a → a· a ab → 2 a· b b· a
= (r u +r u u u ) r r u , (A.4)
a b b a a b ai bj ij
2 · · · →
10
and [11] Martin Lenz. Geometrical origins of contractility in
disordered actomyosin networks. Physical Review X,
V r2r u +r2r u r r (r u +r u )
ab → a b· b b a· a− a· b a· b b· a 4(4):041002, 2014.
= r2r u⊥b+r2r u⊥a, (A.10) [12] PierreRoncerayandMartinLenz. Connectinglocalactive
b a· a a b· b forces to macroscopic stress in elastic media. Soft
where u⊥b is the projection of u onto the space Matter,11(8):1597–1605, 2015.
a a
perpendicular to r , i.e. u⊥b = Pbu , where Pb = [13] S. Feng, M. F. Thorpe, and E. Garboczi. Phys. Rev. B,
b ai ij aj ij 31:276,1985.
δ rˆ rˆ , with rˆ =r /r , is the projectionoperator
ij− bi bj b b b [14] H.HeandM.F.Thorpe. Phys. Rev. Lett.,54:2107,1985.
onto the plane perpendicular to rb. Then with the aid [15] X. Mao and T. C. Lubensky. Phys. Rev. E, 83:011111,
Eqs. (A.9) and (A.10) and the relation 2011.
[16] X. Mao, O. Stenull, and T. C. Lubensky. Phys. Rev. E,
1 1 rˆ Pbu rˆ u 87:042601, 2013.
sin2β − ai ij aj − aj aj [17] X. Mao, O. Stenull, and T. C. Lubensky. Phys. Rev. E,
(cid:18) 0 (cid:19) 87:042602, 2013.
cosβ
= 0 rˆ Pau , (A.11) [18] J.C.Maxwell. Philos. Mag.,27:294,1865.
sin2β bi ij aj [19] A. Souslov, A. J. Liu, and T. C. Lubensky. Phys. Rev.
0
Lett.,103:205503, 2009.
the energy Eb to harmonic order is [20] DavidA.Head,AlexJ.Levine,andF.C.MacKintosh. De-
κ˜ u⊥a u⊥b 2 formationofcross-linkedsemiflexiblepolymernetworks.
Ehar = ˆr a +ˆr b (A.12) Phys. Rev. Lett.,91(10):108102, 2003.
b 2sin2β0 (cid:18) b· ra a· rb (cid:19) [21] D.rAeg.iHmeeasd,oAf.eJl.asLteicvinree,spanondsFe.aCn.dMadceKfoirnmtoasthio.nDmistoidnecst
Note that this energy depends only on displacements of cross-linked cytoskeletal and semiflexible polymer
perpendicular to equilibrium bond directions and thus networks. Phys. Rev. E,68(6):061907, 2003.
it does not induce any bond compression. Setting [22] Jan Wilhelm and Erwin Frey. Elasticity of stiff polymer
networks. Phys. Rev. Lett.,91(10):108103, 2003.
κ˜ =κsinβ ), a=lk, b=lm,andr =r =l leadsto
0 a b 0 [23] P.N.Keating. Phys. Rev.,145:637, 1966.
Eq. (4). If only the vab part were kept, Ebhar reduces [24] N. W. Ashcroft and N. D. Mermin. Solid state physics.
to the form used in reference [14]. Hold,Rinehart,andWinston,NewYork,1976.
When r = r e and r = r e are anti-parallel, [25] M. Das, F. C. MacKintosh, and A. J. Levine. Phys. Rev.
β =π,sin2βa =2Va /R2R2b,and−thbelinearpartofV Lett.,99:038101, 2007.
0 ab a b ab [26] F. Yonezawa and T. Odagaki. Analytic extension of the
vanishes: coherentpotentialapproximationtoclusters. SolidState
Communications, 27(11):1199 –1202,1978.
V(1) =r2r e u r2r e u
ab b a · a− a b · b [27] J.C.Phillips. J. Non-Cryst. Solids,34:153,1979.
+r r (r e u r e u )=0. (A.13) [28] J.C.Phillips. J. Non-Cryst. Solids,43:37,1981.
a b a b b a
· − ·
The quadratic part of V is
ab
1
V(2) = (r u⊥+r u⊥)2 (A.14)
ab 2 a b b a
When u⊥ = (δ e e )u , and E becomes
ai ij − i j a,j b
κ˜V(2)/(r2r2) to harmonic order. If we had chosen
ab a b
r =r r ratherthanr =r r ,thenβ =0,and
b 0 2 b 2 0 0
− −
there would be a minus sign in Eq. (A.14).
References
[1] M.Thorpe. J. Non-Cryst. Solids,57:355,1983.
[2] S.FengandP.N.Sen. Phys. Rev. Lett.,52:216,1984.
[3] S.Alexander. Phys. Rep.,296:65,1998.
[4] M.Wyart. Ann. Phys. (Paris), 30:1,2005.
[5] S.F.EdwardsandD.V.Grinev. Phys.Rev.Lett.,82:5397,
1999.
[6] A.V.TkachenkoandT.A.Witten. Phys. Rev. E,60:687,
1999.
[7] K.D.Hammonds,M.T.Dove,A.P.Giddy,V.Heine,and
B.Winkler. Am. Mineral.,81:1057,1996.
[8] C. P. Broedersz and F. C. MacKintosh. Modeling
semiflexiblepolymernetworks. Rev.Mod.Phys.86,995,
86,2014.
[9] C. P. Broedersz, X. Mao, T. C. Lubensky, and F. C.
MacKintosh. Criticality and isostaticity in fibre
networks. Nat. Phys.,7:983–988, 2011.
[10] M. Das, D. A. Quint, and J. M. Schwarz. PLoS One,
7:e35939,2012.