Network disorder and nonaffine deformations in marginal solids

0.93 MB·English

by  Alessio Zaccone

Network disorder and non-affine deformations in marginal solids
Alessio Zaccone, Jamie R. Blundell and Eugene M. Terentjev
Cavendish Laboratory, University of Cambridge, JJ Thomson Avenue, Cambridge CB3 0HE, U.K.
(Dated: January 12, 2012)
The most profound effect of disorder on the elastic response of solids is the non-affinity of local
displacementswherebytheatoms(particles,networkjunctions)donotsimplyfollowthemacroscopic
strain,astheydoinperfectcrystals,butundergoadditionaldisplacementswhichresultinasoftening
of response. Whether disorder can produce further effects has been an open and difficult question
duetoourpoorunderstandingofnon-affinity. Herewepresentasystematicanalysisofthisproblem
by allowing both network disorder and lattice coordination to vary continuously under account of
2
non-affinity. In one of its limits, our theory, supported by numerical simulations, shows that in
1
lattices close to the limit of mechanical stability the elastic response stiffens proportionally to the
0
degreeofdisorder. Thisresulthasimportantimplicationsinavarietyofareas: fromunderstanding
2
theglass transition problem tothe mechanics of biological networkssuch as thecytoskeleton.

n
a PACSnumbers: 62.20.de,83.80.Fg,81.05Kf
J

1
1 I. INTRODUCTION ical study of the elastic rigidity (described through the
shear modulus) as a function of the degree of structural

] disorder, which varies from the perfect crystal all the
t Inadisorderedlatticesubjecttoasmallexternaldefor-
f way up to the completely disordered lattice. We focus
o mation,everyatomornetworkjunctionmovesinitiallyin
our analysis to central-force systems10 where the phys-
s an affine way, that is, displaced in geometric proportion
. ical bonds between particles have zero bending rigidity,
t to the macroscopic strain. At the same time, the atom
a andtosystemswhichareclosetothevanishingofrigidity
experiencesforcesfromitsnearest-neighbors(NN)which
m (marginal solids). A wide class of physical systems are
areinitially displacedaffinely too. Inany disorderedlat-
- captured by these central force models, such as flexibly-
d tice atoms will therefore experience a non-zero force due
linked biological networks and metallic glasses. In some
n to the fact that their NN are randomly oriented. By
cases,suchascovalentlybondedmaterials,junctions can
o contrast, in a perfectly ordered lattice the forces trans-
possess bending rigidity which makes the transition to
c mitted by the NN may balance each other by symmetry
[ a rigid phase occur at lower connectivity11–13. Without
and the resultant local force in that case can be zero.
loss of generality, our results provide the first direct evi-
1 The finite forces acting on each atom in the disordered
dencefortheexistenceofabroadrangeofstructuraldis-
v networkdriveadditionallocaldisplacements,callednon-
5 affine displacements1–3. As these are accompanied by order where increasing the disorder causes lattices near
2 theisostaticlimit(i.e. thelimitwherethenumberofgeo-
a release of energy, needed in order to reach mechani-
3 metricconstraintsexactlymatchesthenumberofdegrees
cal equilibrium on each atom, the effect of non-affinity
2 of freedom) to stiffen. This result has deep implications
induced by the disorder on the macroscopic rigidity is
.
1 inviewofthefactthatallliquidsontheirwaytowardthe
always to soften the solid with respect to a purely affine
0 solid state do pass through this isostatic limit, whereby
deformation. This softening is a well established phe-
2 even a small difference of mechanical stability with re-
nomenonobservedmanytimesinsimulationsandexperi-
1 spect to internal stresses generated by vitrification can
: mentalstudies4,5. However,sinceaconsistenttheoretical
v decide their ultimate fate, i.e. whether to become per-
descriptionofnon-affineelasticresponsehasbeenelusive
i manently ordered (crystal) or fully disordered (glass).
X sofar,ithasbeendifficulttopredicttherigidityofdisor-
deredsolids(includingglass,granularmedia,cytoskeletal
r
a networks, etc.) as a function of the degree of structural
order/disorder. Apart from its obvious importance for II. THE GAUSSIAN BOND-ORIENTATIONAL
materials science, rigidity is the defining property of the DISORDER MODEL
solid state of matter6 and furthermore the concept of
(generalized) rigidity as the epitome of cooperative be- Let us start by introducing the bond-orientation vec-
havior is essential for various other properties, including tors between nearest-neighbor atoms i and j defined by
how information is transmitted in many-body systems7. n = (cosφsinθ,sinφsinθ,cosθ). The affine part of
ij
In spite of this, the relation between rigidity and disor- the shear modulus of harmonic lattices can be derived
der has remained largely unknown8. This may well be from the expansion of the free energy in the harmonic
the reason for our present lack of understanding of the approximation and is given in this limit by the well-
glasstransitionphenomenon,whereinrigidityemergesin known Born-Huang formula14, for instance, CA =
xyxy
a (disordered) supercooled liquid without any apparent (1/V)κR2 c nxny nxny,wherethelatticesumgoes
change in the symmetry of disordered structure7,9. over all 0nearijesti-jneiijghibjoirj piajirs ij. The atomic force
P
Herewereportonasystematictheoreticalandnumer- constant κ is defined by the harmonic pair-potential,
2
U(r ) = 1κ(r R )2, between nearest-neighbors. R
ij 2 ij − 0 0
is the interatomic distance at rest in the reference frame
and V is the volume of the system before deforma-
tion. Here and below the Roman indices are used to
label atoms while Greek indices are used to label Carte-
sian components of vectors. c is the connectivity ma-
ij
trix with c = 1 if i and j are nearest neighbors and
ij
c = 0 otherwise in our simple case. Let us consider
ij
the shear modulus G of a simple cubic lattice in d = 3
in the affine approximation. Every atom has z = 6
nearest neighbors and shares a bond with each of them.
With N atoms in the lattice, on averageeach atom con-
tributes z nxnynxny to the lattice sum, where the
h ij ij ij iji
bond-orientation averaged factor is evaluated with an
appropriate,normalizedorientationdistributionfunction
(ODF) f(θ,φ):
nxnynxny = f(θ,φ)sin4θcos2φsin2φsinθdθdφ
h ij ij ij iji
Z
(1)
Using this averaging the Born-Huang formula can be
written:
CA =(1/V)κR2z nxnynxny . (2)
xyxy 0 h ij ij ij iji
FIG. 1: (a) Primitive cell of the (perfectly ordered) simple
Evaluating f(θ,φ) for different physical situations is one
cubic lattice. The bonds lying on the horizontal plane are
of the main exercises of our theory, leading to different called azimuthal bonds (az); those on the vertical axis are
forms of the resulting shear modulus. For the simple zenithal (ze). (b) Schematic model of cubic lattice with ran-
cubic system,out ofthe six orientationsof the bonds at- domly added extra bonds (dashed lines). (c) Marginal solids
tachedtoeachatom,onlytwoareindependent,whilethe with varying degree of disorder parameterized by the degree
other 4 can be obtained from these by symmetry opera- of disorder σ definedin thetext.
tions that leave products of nα terms (α = x,y,z ) un-
ij
changed. Thesesymmetryoperationsareeitherrotations
(belonging to the octahedral rotation symmetry group
O ) about the zenithal (z-)axis, or reflections through
h
any plane of symmetry of the cubic cell. One of the two orientations occur with the same probability. This re-
independent orientations lies along the zenithal axis and sults in a completely de-correlated (bond-orientational)
is parameterized by θze = 0 and φze, while the other disorderdescribedbytheuniformODFfUD(θ,φ)=1/4π
one lies in the azimuthal plane and∀ is parameterized by whichleadstohnxijnyijnxijnyijiUD =1/15,whereweusethe
θ = π/2 and φ = 0 (see Fig.1a). The former occurs subscript UD for “uniform disorder”. Hence using this
az az
with probability 2/6 and the latter one with probabil- result in Eq.(2) gives G = CA = (1/30)(N/V)κR2z,
xyxy 0
ity 4/6. Hence the normalized ODF of the simple cubic which is evidently finite for any coordination number
(SC) lattice canbe rewritten, for our purposes,in a sim- z. Thus the affine approximation here dramatically fails
pler normalized form which is going to be used through- to comply with the vanishing of rigidity at isostaticity
out (see the Appendix A for full details): f (θ,φ) = point (i.e. at z = 6), in contrast with both Maxwell’s
SC
(3/2)[(2/6)δ(θ)+(4/6)δ(θ π/2)δ(φ)]. Using this ODF constraint-counting and the empirical evidence8. This
in Eq.(1) immediately lea−ds to nxnynxny = 0 represents the well-documented failure of the affine ap-
and so the affine shear modulushfoijr tihjeijSCijilSaCttice is proximation to describe the elasticity of disordered lat-
GA = CA = (N/V)zκR2 nxnynxny = 0. This tices8,19,20.
xyxy 0h ij ij ij ijiSC
resultgivenbythe affinetheoryforthe perfectlyordered
simple cubic lattice agrees with Maxwell’s counting of Let us now introduce a way to quantify arbitrary de-
degrees of freedom15 which indeed predicts G = 0 for grees of structural disorder. We again consider the sim-
isostatic lattices with z = 2d = 6, and also with evi- ple cubic lattice as our starting point and introduce dis-
dence from polonium-210 in its alpha form, and the su- ordered realizations, always keeping z = 6 for the time
perconductingcalcium-III,whichbothfeaturesimplecu- being. The orientations of the bonds now deviate from
bic structure and a vanishing shear stability at ambient thecubiclatticeaxes,andweassumethattheirdistribu-
pressure16–18. tion is Gaussian about these axes. A particularly conve-
In the opposite limit of a completely disordered net- nientchoicetoparameterizeGaussianbond-orientational
workofharmonicsprings,one may assumethatallbond disorder in the network of this topology is the following
3
ODF: f =(2/6)fze +(4/6)faz , with:
GD GD GD
fze =N (σ2)e−θ2/2σ2
GD ze
(3)
fGazD =Naz(σ2)e−(θ−π2)2/2σ2e−φ2/2σ2
where the subscript GD stands for “Gaussian disorder”
and σ2 is the variance, which acts as our quantitative
measure of lattice disorder. N (σ2) and N (σ2) are
ze az
normalizationfactorswhichalso depend onσ2. We have
verified that using this compact form for the ODF, in
termsofthetwoindependentdirectionsonly,yieldsquan-
titativelythe sameresultsthatonewouldobtainbycon-
sideringexplicitlyallthe 6bondorientations. ThisODF
allows one to span all intermediate degrees of disorder
by varying the parameter σ2. In the limit σ2 this FIG. 2: The theoretical structure-dependent factor in the
→ ∞ shearmodulus,plottedasafunctionofthedegreeofGaussian
ODF flattens and one recoversthe uniform orientational
disorderσ calculatedusingEqs. (1)and(3). Theasymptotic
distribution f (θ,φ) f (θ,φ) = 1/4π. In the op-
GD → UD value of the strong-disorder plateau and the scaling at weak
posite limit of σ2 0, delta-functions are developed at
→ disorder are shown by dashed lines.
the angles of the cubic lattice so that one recovers the
simple-cubic ODF, f (θ,φ) f (θ,φ).
GD SC
→
given strain14, which leads to the Born-Huang formula
for the elastic moduli already used above.
III. ANALYTICAL THEORY OF NON-AFFINE
The negative non-affine contribution can be handled
ELASTIC DEFORMATIONS
analytically(furthertechnicaldetailscanbefoundinthe
Appendix A). It has been recently shown by Lemaitre
A. Simple cubic and uniformly disordered lattices
and Maloney by normal mode decomposition3 that the
non-affine contribution to the modulus takes the form
The ultimate cause of non-affinity lies in the fact that
the particles bonded to a tagged particle, upon the ex- 1 3N (Ξ v )(Ξ v )
ternallyappliedstrainareinitiallydisplacedaffinely,and CNA = ιξ· k κχ· k , (5)
ιξκχ V λ
owing to their displacement they exert a net force on k k
X
the tagged particle. If the particles were placed in an
ordered fashion around the tagged particle, such as in whereΞιξ aretheaffine(force)fieldsactingoneachparti-
crystals, the resultant sum of these forces would be zero cle, while vk andλk arethe eigenvectors(which coincide
duetosymmetry. Inadisorderedsolidsuchacancelation with the normal modes of the solid) and eigenvalues of
does not occur. Hence the resultant force is finite, and the Hessian matrix, respectively. The sum runs over all
it induces an additionaldisplacement onthe taggedpar- translational modes (degrees of freedom) of the atoms.
ticle in order to keep the mechanical equilibrium, which Both the affine force fields and the eigenvectors vk de-
adds to the affine displacement dictated by the imposed pend on the orientations nij of the bonds in the solid in
strain. In this way the mechanical equilibrium is pre- a way geometrically similar to the way the bonds react
servedonallparticles. Clearly,thereisaworkassociated elastically (affinely) to the imposed shear deformation.
with these non-affine displacements, which bears a neg- Accordingly, the non-affine term is also proportional to
ative sign since the overall energy is reduced. Formally the same hnxijnyijnxijnyiji factor that has been present in
this work represents an additional term in the Born free the affine model discussed above. As shown in the Ap-
energy expansion and can be written as pendix A, an allowed and convenient form of the eigen-
modes can be found which leads to the exact evaluation
uNA of the non-affine terms in the sum in Eq.(5). Taking the
1
−W =Xi Z0 fiαdxαi =−2Xij Hiαjβxαixβj(cid:12)(cid:12)x=uNA(4.) atoveGraNgAe o≡veCrxNtyAhxey:disorder of the terms in the sum leads
(cid:12)
The sum runs over all particles and the path integral is GNA =(3N)κR2 nxnynxny /V. (6)
overthenon-affinedisplacementsuNA whichtheparticle 0h ij ij ij iji
i undergoes under the action of the force field f exerted
The important result is that the non-affine part of the
by the neighbors. Hαβ = ∂2U/∂rα∂rβ is the Hessian shear modulus is in fact proportional to the same angu-
ij i j
matrixandxα arethecoordinatesofparticlei. Thenthe lar average as the affine terms discussed above (because
i
free energy expansion for a disordered solid accounting the affine fields which are responsible for the non-affine
for non-affinity is: δF = δF W where δF is the termhave the same dependence on disorderas the affine
A A
−
affine (Born) part of the free energy expansion under a part),andalsotothetotalnumberofdegreesoffreedom
4
in the network, 3N, which all take part in the force re- tors to the shear modulus, z 6 (the coordination) and
laxation. The non-affine effect is, therefore, independent nxny nxny (thedisorde−r). Animportantextension
h ij ij ij ijiGD
of the connectivity z because the non-affine relaxation ariseswhenbondsbetweenfurtherapartatoms(network
involves the degrees of freedom rather than the bonds. junctions) are added on a pre-existing lattice. This sit-
In the limit of uniform disorder we may recall that uation describes higher coordination networks and oc-
GA =(1/30)(N/V)κR2z. Subtractingthecorresponding curs, for example, in cytoskeletal networks where fila-
0
uniform disorder-averaged expression for the non-affine ments are superimposed on an underlying network of a
part, GNA, we obtain the shear modulus as differenttypeoffilaments. Clearly,thebondorientations
oftheaddedfilamentswillbedistributeddifferentlyfrom
1 N 1 3N N
G= κR2z κR2 = κR2(z 6). (7) the underlying network and thus have a different degree
30V 0 − 15 V 0 30V 0 −
of disorder. In order to capture this extension within
This equation is derived in the Appendix A. This mod- our theory, we consider an underlying isostatic network
ulus correctly vanishes at z = 6, that is, at the isostatic (z = 2d = 6) with a variable degree of disorder onto
point of central-force solids where the number of con- which additional bonds are superimposed and placed at
straints(bonds)equalsthenumberofdegreesoffreedom. random at the network’s junctions. This scheme fol-
For z < 6 there is a finite number of floppy modes and lows the most common model employed in the study of
the systems is no longer rigid15. marginalsolids1,22,25andconsistsoftakinganunderlying
Equation (7) is based on the requirement that the ordered isostatic lattice and then randomly place addi-
Hessian averaged over the orientational disorder be a tionalspringstypicallyconnectingnext-nearest-neighbor
3N 3N matrix where the 3 3 submatrices coincide or third-body neighbor atoms (TBN).
× ×
with the identity matrix, as shown in the Appendix A.
Thisrequirementissatisfiedforboththesimplecubiclat-
Let us identify the bonds belonging to the underlying
tice and the uniformly disordered lattice that have been
isostatic lattice as isostatic bonds, while the additional
definedabove,i.e. forboththelowerandtheupperlimit
bonds will be referred to as excess bonds (see Fig.1b).
ofthe Gaussiandisorderspectrumthat weintroducedin
On average, out of z bonds per atom there are 6 iso-
the previous section. Even though Eq.(7) may not hold
static bonds of a base cubic lattice, which we assume is
exactlyfordegreesofdisorderintermediatebetweenthese
formed from the variable (Gaussian) disorder, and z 6
two limits, we shall continue to assume that Eq.(7) still −
excess bonds per atom. Let us assume the excess bonds
applies as an approximate relation in between the two
are placed between TBN atoms along the diagonal of
limits. In the worst case this may give a reasonable in-
the cubic cell, so their average length is of the order of
terpolation over the whole disorder spectrum. As this
√3R . The four third-body neighbor excess bonds in
is a rather uncontrolled assumption, it will be tested in 0
theupperhalf-planearedefinedbyθ =π/4, φ,andφ=
comparisonwith numericalsimulations later in this arti- ∀
π/4,3π/4,5π/4,7π/4 while the four TBN bonds in the
cle.
Since the same orientation average nxnynxny ap- lowerhalf-planeareparameterizedbyθ =π/4+π/2,∀φ,
h ij ij ij iji and φ=π/4,3π/4,5π/4,7π/4. Hence, onaveragea ran-
pears as a factor in both GA and GNA, it follows that
domly placed third-body neighbour contributes a factor
it plays the same role for an arbitrary degree of disor- proportional to nxny nxny = 0.06944. This value
der, that is, G nxnynxny , where the average h ij ij ij ijiEX
∝ h ij ij ij ijiGD should be affected by the disorder and thus should vary
encodes the information about the lattice disorder via
with σ. However, it is evident that this value can vary
f (θ,φ). The dependence of the microstructure func-
GD onlyinthe verynarrowrangefrom0.06944(atσ =0)to
tion nxny nxny upon the degree of disorder as pa-
h ij ij ij ijiGD 0.06667 (at σ = ). Therefore, the variation of this
rameterized by σ2 is explicitly shown in Fig. 2. In the term with σ does∞not lead to a significant change in
limit σ2 , the completely uncorrelated uniform dis- the shear modulus as a function of σ, also because the
order cas→e is∞recovered,hnxijnyijnxijnyijiGD →1/15. In the σ dependence of the ”isostatic” terms is comparatively
opposite limit of weak disorder, the structure function much stronger. This consideration justifies our setting
goes to zero proportionally to σ2, recovering the simple- nxny nxny =const=0.06944 independent of σ in
cubic lattice at σ = 0. The deviation from this propor- thheijfoiljlowijinigjicEaXlculations.
tionalitystartstobesignificantwhenσ2 1,afterwhich
≈
the uniform-disorder plateau (=1/15) sets in.
With this decomposition, we can proportionally split
the affine partof the shear modulus as GA =GA +GA
iso ex
B. Gaussian bond-orientational disorder with where GAiso ≡ 12(N/V)κR02[6/z]zhnxijnyijnxijnyijiGD is due
third-body neighbors to the isostatic bonds which occur with frequency 6/z,
and GA 1(N/V)κR2[(z 6)/z]z nxnynxny is
ex ≡ 2 0 − h ij ij ij ijiUD
We have therefore established that in the simplest sit- due to the excess bonds that occur with frequency (z
−
uation, where all z NN atoms participate of the same 6)/z. We can similarly split the non-affine contribution
degree of disorder, the effects of coordination and disor- into two parts, one due to the 6 isostatic bonds and the
der are decoupled and they contribute two distinct fac- second due to the additional (z 6) randomly placed
−
5
TBN bonds: tionofatomsonthesurfacescaleswiththetotalnumber
of atoms as N−2/3 (i.e. sub-extensively). It is diffi-
1
GNA =GNA+GNA = (N/V)κR26 (8) cult to precisely determine the number of atoms (the
iso ex 2 0
thickness of the surface layer) which follow the strain
6 z 6
nxnynxny +3 − nxny nxny of the box in a purely affine way. To account for this
× z ij ij ij ij GD z ij ij ij ij EX uncertainty, we introduce a parameter B, of order of
(cid:26) (cid:27)
(cid:10) (cid:11) (cid:10) (cid:11) unity, and write this additional affine surface term as
We alsoneedto considera surfaceterm to accountfor GA = 1(N/V)κR2zBN−2/3 nxny nxny . With all
the fact that a finite-size system has the atoms on the S 2 0 h ij ij ij ijiGD
these additional contributions put together, the total
surface being displaced affinely as they are constrained
shear modulus takes the form:
to follow the deformationimposedon the box. The frac-
G= GA +GA GNA GNA +GA
iso ex− iso − ex S
(cid:8)1 6 (cid:9) (z 6)2
= (N/V)κR2 (z 6) nxnynxny +3 − nxnynxny
2 0((cid:20)z(cid:21) − h ij ij ij ijiGD z h ij ij ij ijiEX) (9)
1
+ (N/V)κR2zBN−2/3 nxnynxny .
2 0 h ij ij ij ijiGD
When the number of atoms N is not infinite, the most IV. NUMERICAL SIMULATIONS
interesting situation arises for marginal solids as z 6.
→
In this limit, and in the limit of a large solid (N 1)
≫
where we can assume that the surface term is negligible, To test this conclusion, we performed numerical sim-
wecanalsoneglectthequadraticterm(z 6)2 andhave, ulations where the disorder is introduced as a perturba-
approximately, G (N/V)κR2(z 6) n−xny nxny tiononareferencecubiclattice(Fig.1c). Themagnitude
≈ 0 − h ij ij ij ijiGD
for the non-affine marginal solid with a finite degree of of random perturbation for each atom follows the Gaus-
disorder measured by σ. Guided by the Fig.2, we can sian distribution with a variance σ2 along each of the
split the disorder spectrum into a weak disorder regime three Cartesian axes. The parameter σ, therefore, con-
with 0 < σ < 0.1, and a strong disorder regime with trolsthedegreeofquencheddisorderinthenetwork,with
σ > 1. In the strong disorder regime the shear modulus σ =0correspondingtothecubiclatticeandσ 1corre-
≥
ispracticallyconstant. Remarkably,inthislimitourthe- spondingtoacompletelydisorderedlattice. Thenetwork
ory recovers the result obtained independently from the topologyisformedbyintroducingz bondsplacedatran-
numerical analysis of random packings21 and depleted domonthenearestavailableatomsiteratively,untileach
regular lattices10: G (1/30)(N/V)κR2(z 6). atom has z neighbors. A potential energy quadratic in
≈ 0 −
the displacement is associated with every bond and the
Ontheotherhand,forintermediateandweakdisorder totalfreeenergyiscomputedbyaddingthecontributions
theshearmodulusisstronglydependentuponthedegree from all the bonds. One should note that the perturba-
of disorder. For weakly disordered cubic networks, with tion of the lattice leads to bond lengths that differ from
0<σ <0.1, the shear modulus is proportional to σ2: the equilibrium length R , an effect which is not consid-
0
ered by the theory. However, the contributions from all
thebondsaresummedupandthereforethecontributions
N of shorter and longer (than R ) bonds compensate each
G 0.175 κR2(z 6)σ2 (10) 0
≈ V 0 − other yielding something not dissimilarfromthe equilib-
rium bond length upon summing up the contributions.
A further consequence of assuming that all bonds are
Iftheunderlyingisostaticlatticeisperfectlyordered,i.e. at the equilibrium length is the automatic vanishing of
σ =0,thegeneralEq. (9)hastheonlycontributionaris- the first-order terms in the free-energy Born expansion.
ing from the excess bonds, giving G = G (z 6)2. These terms are related to internal stresses but also in
ex
∝ −
Remarkably,thisscalingagreeswiththerecentresultob- this case, upon averaging over the whole network, the
tained from the numerical study of square lattices with various contributions cancel to give a small net contri-
randomlyplacedextrasprings22. We havenowseenthat bution (a well-documented fact in the glass simulation
the theory predicts that the rigidity of partially disor- literature, see e.g. Tanguy et al.4). In weekly connected
dered lattices near the isostatic point increases upon in- systems (where z <6 and bond-bending interactions are
creasing the degree of structural disorder. active) the situation might be different as discussed in
6
FIG. 4: The shear modulus of simulated networks, for a
fixed value of the disorder variance σ = 1, plotted as a
function of connectivity z, taking only integer values in this
simulation. The arrow indicates the direction of increasing
number of particles: the lattice size from top to bottom is
N = 43,53,63,83,103,123,153 and 203. The solid line is the
theoretical scaling19 G∼(z−6), valid for z →6.
plateau at σ 1. Further, it is seen that non-affinity
≈
causes a substantial softening of the solid (lower shear
modulus) with respect to the purely affine case even
though the dependence upon σ is the same for both
FIG. 3: The shear modulus of simulated isostatic networks affine and non-affine simulations. Varying the connec-
(z = 6) as a function of the disorder variance σ. Red tivity at a fixed degree of disorder, the simulated net-
lines: affine, blue lines: non-affine results. In plot (a) the works display a rigidity transition at the isostatic point.
arrows show the curves changing with the increasing net- As shown in Fig.4, the fundamental scaling G (z 6)
work size, N = 43,53,63,73,83,93 and 103, respectively. due to non-affinity is asymptotically approach∝ed in−the
In plot (b) the arrows show the convergence of simulation large N limit. In Fig.5 theory and simulationresults are
(N=216 as an example) for increasing number of steps, t =
shown together as a function of the degree of disorder
10,30,100,1000,104,3104 and 105.
upon varying the connectivity. The theory is in fairly
good agreementwith the simulations for both z =6 and
z = 7. It should be noted that the surface term (GA in
earlier works (e.g.1). S
Eq. (9)) plays an essential role in this comparison be-
Upon submitting the system to small external defor- cause it ensures that the solid remains marginally rigid
mationsofsimpleshear,themechanicalequilibriumcon- with N finite also at z = 6. Further, it is seen that
dition on each atom/node is enforced such that it fol- upon increasing the connectivity z above the marginal-
lows the imposed deformation by moving along paths of solid limit (z =6) the shear modulus becomes gradually
minimum mechanical energy. This implies, in turn, that lesssensitivetodisorderandeventuallybecomesindepen-
eachatommovesnon-affinelyinresponsetothe external dentofσ. This isanotherimportantandnewresult: the
strain (see Appendix C for detail). From the total free rigidity of central-force lattices with z 7 is very little
energycomputedinthis waywemeasurethe shearmod- affectedby the degree ofstructuralorde≥r/disorderofthe
ulus of the lattices upon varying the degree of disorder lattice. Inthewindowofdisorder0.05<σ <0.5thethe-
σ and the connectivity z. The results are reported in ory appears to underestimate the simulation data which
Fig.3. First of all, these plots confirm that our simula- exhibit a large hump before reaching the plateau corre-
tionsarevalid,bothintermsofthescalingG/N N−2/3 spondingtostrongdisorder. Althoughtherecanbemany
∼
of Eq.(9) (Fig.3a) and the convergence to equilibrium reasonsfor this disagreement,the simulationbeing com-
(Fig.3b) tested. It is especially instructive to note that pletely independent, and using a subtly different way to
when only a short time is allowed for force relaxation, define disorder σ, one possibility also lies in the assump-
the non-affine modulus is very close to the affine result. tion (Appendix A) that the submatrices of the Hessian
Thesimulationresultsreproducethefundamentallaws which determine the form of the eigenvectors (required
predicted by the theory, as summarized in Eq. (9). For to calculate the non-affine term) reduce to identity ma-
the isostatic (marginal-solid) limit with z = 6, and for trices upon averaging, which is an exact procedure only
low to moderate disorder, the shear modulus increases for the two opposite SC and UD limits. We have thus
in proportion to σ2, and then reaches a strong-disorder shown that the rigidity of a cubic lattice increases upon
7
packed (hcp) lattice with z =12 we find that
HCP
1 27√3+3
nxnynxny = 0.0019. (11)
h ij ij ij ijiHCP 16π 512 !≈
As a result, the ratio between the hcp modulus and the
modulus of a uniformly disordered z-coordinated lattice
takes the form:
GHCP = hnxijnyijnxijnyijiHCP 0.0285zHCP, (12)
G nxnynxny ≈ z 6
UD h ij ij ij ijiUD −
where we used that G 0 for the ideal hcp lattice.
NA
≡
If, for instance, the disordered lattice has z = 7, then
FIG. 5: The shear modulus from theanalytical theory (solid
lines)andsimulations(opensymbols)uponvaryingbothdis- GHCP/GUD 1/3, i.e. the shear modulus of an ideal
≈
order variance, σ, and lattice connectivity, z. Red symbols hcp lattice with z = 12 (for which the non-affinity is
refer to the affine moduli, blue symbols to the moduli ac- certainly negligible) is almost three times smaller than
counting for non-affinity, for the two discrete cases z = 6 the shearmodulus ofa lattice withuniformdisorderand
and z =7. The continuous variation of z above the isostatic z = 7, in spite of the hcp lattice being more strongly
point z =6 shows how the effect of disorder diminishes, and connected. The difference becomes more pronounced
gradually disappears at higher connectivity.
upon increasing the coordination of the disordered lat-
tice. Therefore the disorder-induced stiffening of crystal
lattices is a generaleffect not limited to the cubic lattice
introducing structural bond-orientational disorder, with studied in detail above.
a law of direct proportionality to the disorder variance. Secondly,thesefindingsmayhaveconsequencesforthe
Clearly, the disorder-induced stiffening is an effect con- glass transition of supercooled liquids. Descriptive the-
trolled by the orientation average nxnynxny , which ories of the glass transition have emphasized the role of
h ij ij ij iji
hasasimplegeometricalmeaning: itmeasurestheextent local icosahedral symmetry which emerges at the glass
to which the bond-vectors of the lattice are misaligned transition26. Icosahedral clusters have also been identi-
with respect to the shear-deformation axes, x and y. In fied with the cooperatively rearranging regions of mean-
fact, the product nxnynxny vanishes for the simple field(e.g. Adam-Gibbs)theories27. Thegrowthandjam-
h ij ij ij iji
cubiclatticewhereallbondsareparalleltoeitherthexor ming of such icosahedral clusters lead to the sudden in-
theyaxis. Ontheotherhand,theproduct nxnynxny crease of the viscosity and eventually to the structural
h ij ij ij iji
is maximized and saturates to its highest value (=1/15) arrest. In the icosahedral clusters an atom at the center
whenthebondstendtoberandomlyoriented. Inthecu- is bonded to its 12 nearest neighbors, forming an icosa-
bic lattice the vanishing of nxnynxny and the shear hedron where none of the bonds can be parallel with ei-
h ij ij ij iji
modulus with it, reflects the fact that there is no restor- ther the x or y directions. Instead, clusters with hcp
ing force if all parallel columns in either the x or y di- symmetry could be formed which have the same num-
rection are tilted (maintaining the constant separation ber of bonds as the icosahedra but 6 of these bonds
of neighbors), which is an example of a floppy mode. In can be aligned with the x and y axis. The same ap-
disordered lattices, with bonds misaligned with respect plies to fcc clusters. Therefore it follows that the factor
to x and y, the longitudinal components of the bond dis- nxny nxny is much larger for the icosahedra than for
h ij ij ij iji
placements cannot be neglected (they are no longer sec- eitherhcporfccclustersandthustheicosahedraaremore
ond order with respect to the displacement magnitude) rigidthanthehcporfccclusters. Thisimpliesthat,dur-
which implies a finite restoring force and thus a finite ing the supercooling process, the icosahedra, which are
shear modulus. From a different perspective, the fac- the precursor of the rigid glass state28, are able to stand
tor nxny nxny measurestheaveragelocal(geometric) the internal stresses generated by the vitrification1 bet-
h ij ij ij iji
frustrationofthe bonds and it takesa finite value when- ter than the hcp or fcc clusters which are the precursor
everthebondsnolongercanbedisplacedwithoutpaying of the ordered crystal. This is an unprecedented obser-
a finite energy cost. vationthatshedslightonthephysicaloriginoftheonset
of rigidity at the glass transition.
In further developments, this theory can be applied
to analyze the onset of rigidity in systems which possess
V. CONCLUSIONS medium-rangeorder29,suchasthetechnologicallyimpor-
tantamorphoussemiconductorsandmetallicglasses30,31.
Firstofall,theconclusionthatbond-orientationaldis- In principle, the orientation distribution function for
order causes stiffening of an ordered lattice can be ex- partially disordered lattices (e.g., those containing de-
tended to other crystal lattices by re-evaluating the fac- fects such as dislocations causing distortion of the bond-
tor nxnynxny . For example, for the hexagonal close orientations,ortheabovementionedsolidswithmedium-
h ij ij ij iji
8
range order) can be extracted from structural investiga- sumoverthe3N eigenmodesoftheHessianordynamical
tion of the material in question (e.g., scattering or mi- matrix of the lattice. The 3N 3N Hessian matrix for
×
croscopy techniques) and used as input in our theory to harmonic lattices is given by
investigate the elastic properties as a function of struc-
ture. Hiαjβ =δij κcisnαisnβis−(1−δij)κcijnαijnβij (A2)
In summary, by combining theory and simulations we s
X
have established a fundamental law which governs the
Eq.(A2) follows from replacing the harmonic potential
rigidityofsolidsasa functionofstructuraldisorder. Ac- U(r ) = 1κ(r R )2 in the definition of the Hessian
cordingtothislaw,thestiffnessoflatticesincreaseswith ij 2 ij − 0
matrix: Hαβ ∂2U/∂rα∂rβ. R is the rest length
the degree of bond-orientational disorder. Our finding ij ≡ i j 0
of the bonds. v and λ in Eq.(A1) are eigenvectors
suggests a new way of increasing the mechanical sta- k k
and eigenvalues of the Hessian, respectively. The inner
bility of technologically important marginal-solids such
product (Ξ v ) is the projection of the affine force
as e.g. Po210 and the superconducting Ca-III, both of ιξ · k
field Ξ (i.e., the force field exerted on every atom by
which exhibit simple cubic lattices which therefore can ιξ
the affine motions of its neighbours) on the eigenvec-
be stabilized by introducing structural disorder (e.g. in
tor v . The analytical form of the affine fields is given
the form of defects). Furthermore, our findings suggest k
by3: Ξα = R κc nαnκnχ. Thus, the evalua-
a new way of looking at the unresolved problem of the i,κχ − ij ij ij ij ij
j
glass transition. In light of our results, vitrifying liquids tion of the non-aPffine term in the elastic moduli reduces
uponcrossingthemarginal-solidityorisostaticthreshold to the task of evaluating the eigenmodes of the Hessian,
find themselves in disordered configurations which are v = a w where r = 1...N and l = x,y,z. In gen-
k r l
intrinsically more rigid than the ordered ones. eral,there areno analyticalroutesto evaluatethe eigen-
N
modes. This becomes possible however if one chooses
w = e where e is the standard Cartesian basis of R3,
Acknowledgments. l l l
as we show below. This treatment is exact in the case
of the simple cubic (SC) lattice where the Hessian is ex-
This research has been supported by the Swiss actly given by Hαβ = κ δ c (1 δ )c δ
NationalFoundation (PBEZP2-131153)and the EPSRC ij 3 ij j ij − − ij ij αβ
TCM Programme Grant. We appreciate useful discus- andclearlythe vk =ar (cid:16)el arPeeigenvectorsofthe(cid:17)Hes-
sions and other input from S.D. Guest, T.C. Lubensky, sianbecausethesubmatricesoftheHessianarediagonal.
N
M. Warner, H. J. Herrmann, and R. Blumenfeld. In the case of the uniformly disordered (UD) (isotropic)
lattice as defined in the main article, it can be shown
thattheeigenvectorshavethesameform. Ageneralthe-
orem states that if v and λ are respectively the eigen-
k k
Appendix A: Non-affine deformations vectors and eigenvalues of a matrix A, the same eigen-
vectors are shared by a matrix f(A), for any function
f of the matrix. The eigenvalue equation for the lat-
In the following we give a brief outline of the analyti-
ter matrix is then given by: f(A)v = f(λ )v . Let us
cal method that we devised to calculate the elastic con- k k k
take f(...) = ... , where ... denotes the average over
stantsofdisorderedharmoniclatticeswhichaccountsfor
h i h i
non-affinity3,19. The method makes use ofthe non-affine the bond-orientational disorder, as usual. Then we have
linear response formalism for disordered solids which is f(H) = κ δ c (1 δ )c δ , where δ is
3 ij j ij − − ij ij αβ αβ
described in detail by Lemaitre and Maloney3. One of theKroneck(cid:16)er’sdelta. Furthermore,(cid:17)λ =mω2 whereω
thecentralresultsoftheformalism,whichrepresentsthe P k k k
arethenormalmodefrequenciesofthe solid. Thelattice
starting point of our derivation, is the exact lattice sum
with uniform disorder is an isotropic solid, which means
for the elastic moduli accounting for non-affinity3
that the normal mode frequencies depend only on scalar
C =CA CNA (A1) quantities. Clearly the same applies to the λk which in
ιξκχ ιξκχ− ιξκχ turnimplies: f(λ )=λ . Thisresultestablishesthatthe
k k
R2κ 1 3N (Ξ v )(Ξ v ) eigenvectorsofthe form: v =a e are alsoeigenvec-
= 0 c nι nκnξ nχ ιξ· k κχ· k k r l
2V ij ij ij ij ij − V λ tors of the Hessian for the uniformly disordered lattice.
Xij Xk k Hence, in both the SC and UD limNits one can write the
The superscripts A and NA denote the affine and non- eigenvalue equationfor the Hessian as: (H˜ I)(a e)=
⊗ ⊗
affinepartsofthemoduli,respectively. Theaffinepartis λ(a e), where H˜ κ δ c (1 δ )c and I
the standard Born-Huang expression in terms of a lat- ⊗ ≡ 3 ij j ij − − ij ij
denotes the 3 3 identit(cid:16)y matrix. Then the inne(cid:17)r prod-
tice sum over nearest-neighbors (NN) i and j, where × P
ucts of the affine fields with the eigenmodes become19:
n =(cosφsinθ,sinφsinθ,cosθ)istheunitvectordefin-
ij
ingtheorientationofbondij. cij istheadjacencymatrix (Ξ v )= (A3)
(c = 1 for two NN atoms and c = 0 otherwise) and ιξ· k
ij ij
κ the harmonic spring constant of the harmonic inter- = κ2R02 arar′crscr′s′nlrsnιrsnξrsnlr′s′nκr′s′nχr′s′
action between two NN atoms. The non-affine part is a rs,r′s′
X
9
where the sum runs over two pairs of NN atoms to evaluate the orientation average in Eq.(A5). In the
at the time, rs and r′s′. Consistent with the worstsituation,thisapproximationstillshouldprovidea
main approximation of this theory, we replace the reasonableinterpolationsince it correctly describes both
orientation-dependent terms with their orientation av- the lower (SC) and upper (UD) limits of the disorder
erage. One finds that nl nι nξ nl nκ nχ = spectrum. The validity of the approximation has been
h rs rs rs r′s′ r′s′ r′s′i
(δrr′δss′ −δrs′δsr′)Bl,ιξκχ where Bl,ιξκχ =hnιijnξijnκijnχiji cahrteicckleedanbdy caosmexppaercistoendwititphrosvimiduelsataiognooddatdaeisncrtihpetiomnaionf
and ... = ...f(θ,φ)sinθdθdφ denotes the averageover
h i the data apart for a window of disorder 0.05 < σ < 0.5
bondorientationsaccordingtosomeorientationdistribu-
R located in the middle of the spectrum.
tion function f(θ,φ). Replacing this averagein Eq.(A3),
Appendix B: Gaussian disorder
one obtains:
(Ξ v )(Ξ v )= (A4) Thedistributionfunctionfortheorientationofthelat-
ιξ· k κχ· k
tice bonds (ODF) f(θ,φ) is defined by the following re-
= κ2R2B a2c c a a c c lation
0 l,ιξκχ r rs rs− r s rs sr!
rs rs
X X
... = ...f(θ,φ)sinθdθdφ (B1)
Furthermore, we have that c2 = c c = c and the h i
rs rs sr rs Z
N N N
identities a2r crs− arascrs = aras[( crj)δrs− To evaluate the shear modulus we need to evaluate its
r s rs rs j structure-dependent part
P P N P P P
c (1 δ )] = 3 a a H˜ , where we defined H =
rs − rs κ r s rs
rs nxny nxny = f(θ,φ)sin4θcos2φsin2φsinθdθdφ
H˜ I and I is thPe 3 3 identity matrix. Recalling h ij ij ij iji
⊗ × Z
N (B2)
that H˜rsas = λar we obtain (Ξιξ·vk)(Ξκχ·vk) = All the information about the structure of the lattice is
s
contained in f(θ,φ). For isostatic lattices in d = 3, as
κR02λPkBl,ιξκχ UsingthisresultinEq.(A1),thenon-affine explained in the main body of the article, the degree of
part of the shear modulus follows as:
structural disorder can be varied continuously from the
1 3N (Ξ v )(Ξ v ) simplecubiclatticetothecompleteuncorrelateddisorder
GNA xy· k xy· k by using the following Gaussian ODF
≡ V λ
k
k
X 2 4
= 1 N 3 3κR02λkBl,ιξκχ (A5) fGD(θ,φ)= 6fGzeD(θ,φ)+ 6fGazD(θ,φ) (B3)
V λ
k
Xk=1Xl=1 with
=3NκR2 nxny nxny /V
0h ij ij ij iji fze (θ,φ)=N (σ2)e−θ2/2σ2
GD ze (B4)
whichis akeyresultusedin the mainpartofthe article. faz (θ,φ)=N (σ2)e−(θ−π/2)2/2σ2e−φ2/2σ2
This theory strictly applies to the SC and UD lim- GD az
its only because only in these specific cases the angular
The normalization factors N and N are defined by
average nαnβ = δ /3 is exact. However, we make ze az
h ij iji αβ
the additional assumption that for the intermediate de- N (σ2)fze (θ,φ)sinθdθdφ=1
ze GD (B5)
greesofdisorderthis resultstillholdsapproximately,i.e. N (σ2)faz (θ,φ)sinθdθdφ=1
nαnβ δ /3for0<σ < . Withinthisapproxima- R az GD
h ij iji≃ αβ ∞ R
tion Eq.(A5) can be used to evaluate the elastic moduli CombiningEq.(B3)-(B5)weobtainthe followinganalyt-
of lattices with disorder (parameterized by variable σ), ical expression for the microstructure-dependent factor
by using the Gaussian ODF defined in the main article in the shear modulus:
10
nxnynxny = f (θ,φ)sin4θcos2φsin2φsinθdθdφ
=h ije−ij12σij2 ij1i/GσD2√σ2(EGrDf[(π 5iσ2)/√2σ2]+Erf[(π+5iσ2)/√2σ2]
+5{e8σ2(Erf[(π 3iσR2)/√2σ2]− Erf[(π+3iσ2)/√2σ2]
2e4σ2(Eprf[(π− iσ2)/√2σ2]−Erf[(π+iσ2)/√2σ2]+2iErfi[ σ2/2])
− − −
+2iErfi[(3√σ2)/√2]) 2iErfi[5√σ2/√2])
− } p
384(Erf[((π+iσ2) 1/2σ2] 2iErfi[1/ 2/σ2]
×{ −
+iErfi[( 1/σ2(iπ+σ2))/√2]) −1
+ e−20σ2(5e8σ2(2e4σp2(Erf[(π} 2iσ2)/2√p2σ2]+Erf[(π+2iσ2)/2√2σ2]) (B6)
{ p −
+Erf[(π 6iσ2)/(2√2σ2)]+Erf[(π+6iσ2)/(2√2σ2)])
−
+Erf[(π 10iσ2)/(2√2σ2)]
+Erf[(π+−10iσ2)/(2√2σ2)])( 2+2e8σ2Erf[(√2π)/√σ2]
−
+Erfc[(√2(π 2iσ2))/√σ2]+Erfc[(√2(π+2iσ2))/√σ2]))
− }
768Erf[π/√2σ2](Erf[(π 2iσ2)/2√2σ2]
×{ −
+Erf[(π+2iσ2)/(2√2σ2)]) −1
}
Eq.(B6) is a function of the Gaussian variance σ2 only where we identified:
and has been plotted in Fig.2 in the main body of the
article. fGazD ≡fG{Dπ2,0} ≡Naz(σ2)e−(θ−π2)2/2σ2e−φ2/2σ2 (B8)
fze f{0,∀φ} N (σ2)e−θ2/2σ2
GD ≡ GD ≡ ze
Derivation of the Gaussian ODF
i.e. the first terms of the two sums in Eq.(B7), defined
by p = 0 and q = 0, respectively. Let us start with the
Here we provide a proof of the expression given by
orientations in the xy plane. The principal azimuthal
Eq.(2) in the main article (Eqs.(B3)-(B4) in the sec- orientation faz , with p = 0, contributes to the average
GD
tion above) for the Gaussian ODF used to vary the de- nxny nxny the following integral:
gree of disorder in our model. We start by consider- h ij ij ij ijiGD
ing the orientations of the bonds in the primitive cell of
the simple cubic (SC) lattice. There are 4 bonds lying sin4θcos2φsin2φe−(θ−π2)2/2σ2e−φ2/2σ2sinθdθdφ
in the xy plane that are defined by pairs (θ,φ) of an- Z
gles (π/2,0), (π/2,π/2), (π/2,π), and (π/2,3π/2). Fur- + sin4θcos2φsin2φe−(θ−π2)2/2σ2e−(φ−2π)2/2σ2sinθdθdφ
ther, there are two bonds lying along the polar z-axis
Z
defined by (0,∀φ) and (π,∀φ). The Gaussian model =2 sin4θcos2φsin2φe−(θ−π2)2/2σ2e−φ2/2σ2sinθdθdφ
realizes the situation where each of these bonds is dis-
Z
tributed around the simple-cubic angle pairs mentioned (B9)
above according to a Gaussian distribution f with whereweomittedthenormalizationfactorandwherethe
GD
variance σ2. In the limit σ 0 the model devel- equality holds because the functions in the two integral
→
ops Dirac deltas at the SC bond orientations thus re- on the l.h.s. subtend the same area within the same
covering the perfectly ordered SC lattice. We wish to intervalofintegration[0,2π]. Itshouldbenotedthatthe
demonstratethatinordertoevaluatetheangularaverage second term on the l.h.s. is strictly required in order to
nxnynxny = ...f (θ,φ)sinθdθdφ, one needs have a Gaussian centered on φ = 0. Without this term,
h ij ij ij ijiGD GD
only to consider the two principal orientations (π/2,0) since integrationgoes from 0 to 2π, one would have only
and (0, φ) (from whiRch the other ones can be obtained half of the Gaussian and would not properly count the
uponap∀plicationofthesymmetryoperationsofthegroup contributions from the angles which lie close to φ = 2π
O ), which we denote as azimuthal and zenithal respec- (or equivalently on the negative axis close to φ=0−).
h
tively,becausealltheotherorientationscontributeterms The secondtermin the firstsum in Eq.(B7), with p=
which are identical to either of these two. In other 1, contributes to the average the following term:
words, we want to show that for the purpose of evalu-
ating hnxijnyijnxijnyijiGD the following identities hold: πsin5θe−(θ−π2)2/2σ2dθ 2πcos2φsin2φe−(φ−π2)2/2σ2dφ
Z0 Z0
3 1 (B10)
f 1 f{π2,pπ2}+ f{qπ,∀φ} = 4faz +2fze In this integral,let us change to the variable Φ φ π.
GD ≡ 6 p=0 GD q=0 GD ! 6 GD 6 GD Upon applying the well-known trigonometric r≡elati−on2s:
X X
(B7) sin(φ+ π)=cosφ and cos(φ+π)= sinφ, the integral
2 2 −