Approximate analytical description of the nonaffine response of amorphous solids

Alessio Zaccone1 and Enzo Scossa-Romano2

1
Department of Physics, Cavendish Laboratory, University of Cambridge,
JJ Thomson Avenue, Cambridge CB3 0HE, United Kingdom and
2
Department of Chemistry and Applied Biosciences, ETH Zurich, CH-8093 Zürich, Switzerland
(Dated: November 2, 2018)
An approximation scheme for model disordered solids is proposed which leads to the fully analyt-
ical evaluation of the elastic constants under explicit account of the inhomogeneity (nonaffinity) of
the atomic displacements. The theory is in quantitative agreement with simulations for central-force
systems and predicts the vanishing of the shear modulus at the isostatic point with the linear law
µ ∼ (z − 2d), where z is the coordination number. The vanishing of rigidity at the isostatic point is
arXiv:1102.0162v2 [cond-mat.soft] 22 Aug 2011

shown to be the consequence of the canceling out of positive affine and negative nonaffine terms.

PACS numbers: 46.25.-y, 64.60.aq, 45.70.-n

I. INTRODUCTION butions to the elastic constants. This leads to a fully

analytical description of the elastic constants which ac-
counts for the microscopic nonaffinity of the atomic dis-
Disordered or amorphous solids represent a great part
placements.
of ordinary matter (e.g., glass), including biological mat-
ter (e.g., cytoskeletal networks) [1–3]. Yet, the relation-
ship between rigidity and disorder has remained elusive II. FORMALISM
and no theory has hitherto proved able to correctly de-
scribe their elastic constants. The rigidity of disordered
solids is also intimately related to the fundamentally un- In the following, Roman indices are used to label atoms
solved problem of the glass transition [4]. Elastic rigid- while Greek indices are used to label Cartesian compo-
ity in supercooled liquids emerges from the fluid state nents. The summation convention over repeated indices
at the glass transition without any detectable lowering holds throughout for Greek indices. Bold characters de-
of the symmetry (apart from translational and replica note vectors in dN -dimensional space (N =total number
symmetry-breaking), as opposed to what happens in ”or- of atoms). We closely follow the notation of Lemaitre
dinary” liquid-solid phase transitions [5]. In recent years and Maloney [8] and we start our analysis from the def-
it has been recognized that the intrinsic ”softness” of dis- inition of a Bravais cell for the disordered lattice. The
ordered solids is related to the nonaffinity of the atomic cell is described by three Bravais vectors and thus by a
displacements [1, 4, 6, 7]: the atoms in a strained disor- matrix h. The potential depends on the particle position
dered solid are not displaced proportionally to the global ri and on the shape of the cell which enforces bound-
strain. The calculation of the contribution to rigidity due ary conditions. Macroscopically imposed deformations
to such nonaffine displacements poses formidable difficul- of the cell are described by changes in the Bravais vec-
−1
ties because it requires the analytical knowledge of the tors through a linear map F = h h̊ and the relation
eigenmodes of the dynamical or Hessian matrix of the h = F h̊. We denote quantities in the reference frame, as
system [8, 9], which is a sparse random matrix. The well as quantities that are measured with respect to the
simplest disordered solids where nonaffinity is supposed reference frame, with a circle. Accordingly, the unit cell
to play a significant role are those with nearest-neighbor in the reference configuration is given through the ma-
central-force interactions [9, 10]. In these systems, it is trix h̊. In the language of continua F is the deformation
well-known that the shear modulus µ vanishes at the iso- gradient tensor. The position of atom i after an affine
static point where each particle has an average number of deformation is given by
nearest neighbors (mechanical contacts) z = 2d [9, 10].
Although it seems reasonable from constraint-counting riα = Fαβ Riβ (1)
arguments that rigidity is lost when z = 2d [11], the As a unique exception to the ring notation, we denote
linear law µ ∼ (z − 2d) observed in simulation studies with Ri the position of the atoms in the reference cell.
of completely disordered solids [12] has remained un- This relation illustrates the definition of affine displace-
explained and the physics behind it represents a long- ment as atom i is displaced proportionally to the external
standing problem [9, 10] where the role of nonaffinity is deformation. The position of the atom after a strain is
yet unclear. Further, the well-documented inadequacy denoted by riα (F ) and differs from Eq.(1) if the nonaffine
of affine theories to describe the elasticity and transport displacement is not zero. It is useful to introduce the
properties of amorphous materials calls for an improved particle position r̊i for an atom which undergoes both
theory beyond the affine approximation [13]. In this Let- affine and nonaffine displacements, defined by
ter, we propose an approximation scheme which gives a
well-defined deterministic limit for the nonaffine contri- riα (F ) = Fαβ r̊iβ (F ). (2)
 2

Remark also that Ri = r̊i (0). Hence, for an affine dis- The elastic constants are defined by
placement the position vector r̊i is kept fixed while the
new position is determined by F . Any additional (non- 1 ∂2U
Cιξκχ = (6)
affine) displacement is thus parameterized for a given V̊ ∂ηιξ ∂ηκχ η=0
strain in terms of r̊i (F ). With these definitions we can
express the potential U({ri }, F ) in the coordinates of the In order to account for the nonaffine relaxation in the cal-
reference frame as Ů({r̊i }, F ) which is defined by culation of the elastic constants, the derivatives in Eq.(6)
have to be taken along a trajectory of (locally) minimum
Ů({r̊i }, F ) = U({F r̊i }, F ). (3) energy. Following Lemaitre and Maloney [8], we obtain
" !#
It is convenient to introduce the Cauchy-Green strain 1 D ∂ Ů ∂ Ů Dr̊i
tensor η = 21 (F > F − I) to describe the deformations Cιξκχ = +
V̊ Dηιξ ∂ηκχ ∂r̊i Dηκχ
η=0
since the elastic constants are defined in terms of second !
derivatives in η. The deformation is completely described 1 2
∂ U 2
∂ U Dr̊i
by this tensor since the total internal energy of the solid = +
V̊ ∂ηιξ ∂ηκχ η=0 ∂r̊i ∂ηιξ η=0 Dηκχ η=0
being deformed can be expressed as U({|rij |}), i.e. as a
functional of the set {rij } of relative distances between 1 ∂2U 1 Dr̊i
= − Ξi,ιξ
atoms rij = |rij | = |ri − rj |, and η describes affine trans- V̊ ∂ηιξ ∂ηκχ η=0 V̊ Dηκχ η=0
formations in terms of relative interatomic distances ac- = A
Cιξκχ − NA
Cιξκχ
cording to |F Rij |2 = |Rij |2 + R> ij η Rij . (7)
A homogeneous strain of the cell F will first bring each where it is evident that the true elastic constant is given
atom to its affine position riα = Fαβ Riβ . In this affine by the affine (Born-Huang) elastic constant Cιξκχ A
cor-
position the total force acting on atom i is in general NA
rected by the nonaffine term −Cιξκχ . Following Lemaitre
not zero since the neighboring atoms may exert a non- and Maloney [8], and using Eq. (4) in Eq. (7), one derives
vanishing force-field due to their affine motion. This is the following expression for the nonaffine correction,
especially true for disordered solids where the nearest
NA αβ −1 β
neighbors are placed at random around atom i so that Cιξκχ = Ξα
i,ιξ (Hij ) Ξj,κχ > 0 (8)
they transmit unbalanced forces to i (in an ordered lat-
tice the transmitted forces balance by symmetry such The last inequality in Eq. (8) is justified in view of the
that this effect is often negligible). It is in response to Hessian matrix being semi-positive definite at mechanical
these virtual forces that the atom undergoes an addi- equilibrium. Hence it follows that the correction due to
NA
tional motion after it has been displaced affinely such the nonaffine relaxation, −Cιξκχ < 0, necessarily gives a
that the energy released in the process reestablishes (lo- negative contribution to the total rigidity [8].
cal) mechanical equilibrium. Hence the system under
reversible strain evolves adiabatically along a trajectory
r(η) that minimizes the mechanical energy for a given III. APPROXIMATION SCHEME
D
strain η. If we denote by the derivative with respect
Dηκχ
A. The Cauchy bonded-network model
to adiabatic changes of the strain under the constraint of
mechanical equilibrium, one obtains an equation of mo-
tion of the nonaffine displacement by differentiating the Let us consider the disordered Cauchy solid, defined
∂U
force fiα = − ∂r α evaluated in the true position (where it
by the following properties [1]: (i ) atoms interact pair-
i
wise and only with their nearest neighbors; (ii ) the in-
vanishes). In the limit η → 0 the equation reads
teraction potential is a central-force harmonic potential;
(iii ) the reference state is unstressed, i.e all springs (in-
X Dr̊jβ
αβ
Hij = Ξα (4) teratomic bonds) are relaxed in the minimum of the
i,κχ
j
Dηκχ harmonic well; (iv ) disorder is spatially decorrelated.
η=0
The equivalence with a random network of harmonic
where we have introduced the Hessian Hijαβ
and the affine springs is evident [9, 14].
P Hence, the total free energy is
α
force field Ξi,κχ given respectively by given by U({rij }) = hiji Vij (rij ) where the sum runs
over all pairs of nearest-neighbors hiji. The pair in-
αβ ∂ 2 Ů ∂2U teraction potential is given by the harmonic potential
Hij = = V (rij ) = κ2 (rij − R0 )2 . κ is the atomic force constant
∂r̊iα ∂r̊jβ η=0 ∂riα ∂rjβ η=0 and R0 is the interatomic distance at rest in the refer-
ence frame. Under these conditions the Hessian matrix
∂ 2 Ů X ∂ 2 Ů ∂Fββ 0 β 0
Ξα
i,κχ = − =− R becomes
∂r̊iα ∂ηκχ j ∂r̊iα ∂r̊jβ η=0 ∂ηκχ j X
η=0 αβ β α β
Hij = δij κcis nαis nis − (1 − δij )κcij nij nij (9)
(5) s
 3

where we used the identity ∂/∂rij = nij ∂/∂rij , with implies sacrificing some details of the vibrational spec-
nij = rij /rij . Further, cij is the (random) occupancy trum. This problem is addressed in section IV.B and
matrix with cij = 1 if i and j are nearest neighbors and IV.C where we show what details are lost and we study
cij = 0 otherwise. cij is a matrix where each row and the validity and limitations of the approximation.
each column have on average z elements equal to 1 dis-
tributed randomly under the constraint that the matrix In d = 3 it is nij = (cos φij sin θij , sin φij sin θij , cos θij )
be symmetric. Using this form of the Hessian one obtains and the pair of angles φij and θij univocally specifies the
the affine part of the elastic constant as orientation of the bond hiji. The orientation-dependent
β
factors in the Hessian, nα ij nij in Eq. (9), for a large
R02 κ X
A
Cιξκχ = cij nιij nξij nκij nχij (10) system with uncorrelated isotropic disorder (where every
2V̊ ij bond can take any orientation in the solid angle with the
same probability 1/4π), can be replaced with its isotropic
which is the well-known Born-Huang formula [1, 7, (angular) average, i.e. nα β
ij nij ⇒ δαβ /d. Within this ap-
8]. Further, we also obtain a microscopic expres- proximation, the Hessian becomes
sion for the affine force field from Eq.(5) as Ξα i,κχ =
κ χ
− j Rij κcij nα
P
ij nij nij . We can now turn to the non-
NA
affine part of the elastic stiffness, Cιξκχ . The Hessian is  
a dN × dN symmetric semi-positive definite matrix with αβ κ X
Hij = δij cij − (1 − δij )cij  δαβ (12)
d eigenvalues equal to zero which are due to the global d j
translational invariance of the solid. Eq.(4) can be solved
by normal mode decomposition which leads to [8]

NA 1 X (Ξιξ , vk )(Ξκχ , vk ) Remark that this is still a sparse random matrix because
Cιξκχ = (11)
V̊ λk of the positional disorder in the random coefficients cij .
k
λk 6=0 According to Eq.(12), let us define H = H̃ ⊗ I where
where vk are the eigenvectors of the Hessian (which are I is the d × d identity matrix (which represents δαβ )
orthogonal since the Hessian is symmetric), λk the corre- and H̃ is the matrix which multiplies δαβ in Eq.(12).
sponding eigenvalues and (, ) denotes the normal scalar Denoting with {aq }q=1..N the set of eigenvectors of H̃,
product on RdN . In the next section, we shall evalu- which is an orthonormal basis (ONB) of RN , and with
ate the deterministic limit of Eq.(11), which is a self- {el }l=1..d the standard Cartesian basis of Rd , it follows
averaging quantity [8]. that (H̃ ⊗ I)(a ⊗ e) = λ(a ⊗ e) and thus the dN dimen-
sional set {aq el }q=1..N,l=1..d is an ONB of eigenvectors
B. The approximate Hessian and the affine field
of H as given by Eq.(12). This allows us to write (with
projection on its eigenmodes v = a el for some a ∈ {aq }q=1..N ):

A rigorous derivation requires one to first determine ! !

the eigenmodes vk of the Hessian matrix given by Eq.(9) N
X N
X



in order to calculate the projection on them of the affine Ξιξ , v Ξκχ , v = ar Ξr,ιξ el ar Ξr,κχ el
fields in Eq.(11). Thereafter the average over the dis- r r
NA X
order is taken to get the thermodynamic limit of Cιξκχ . = κ2 R02 {(ar ar0 crs cr0 s0 )
However, the Hessian is a random matrix and both vk r s r 0 s0
and λk depend on the realization of disorder. Also, be-
× (nlrs nιrs nξrs nlr0 s0 nκr0 s0 nχr0 s0 )}
ing the Hessian sparse, there are no analytical forms even
(13)
for its statistical spectral distributions. Nevertheless, the
With our isotropic approximation, we replace the
deterministic limit of Eq.(11) can be calculated analyti-
orientation-dependent terms with their isotropic angular-
cally within the following approximation that we propose
averaged values which gives nlrs nιrs nξrs nlr0 s0 nκr0 s0 nχr0 s0 =
here.
(δrr0 δss0 − δrs0 δsr0 )·Bl,ιξκχ where the Bl,ιξκχ are geomet-
Our approximation consists in performing a disorder-
ric coefficients resulting from the angular average. For
average of the orientation-dependent part of the Hessian
d = 3 and d = 2 they are as follows
first and then use the result to calculate the eigenmodes,
their inner products with the affine fields and finally the
nonaffine correction. Inverting the sequence of ”calculat-
ing” and ”averaging” is sometimes referred to as an ef- d=3 P d = 2P
fective medium approximation and is not at all unusual l x y z l x y l
1 1 1 1 5 1 3
in dealing with disordered systems [15] since it is often Bl,xxxx 7 35 35 5 16 16 8 (14)
1 1 1 1 1 1 1
the only strategy to keep the treatment analytical. Us- Bl,xyxy 35 35 105 15 16 16 8
ing some averaged form of the Hessian matrix necessarily 1 1 1 1 1 1 1
Bl,xxyy 35 35 105 15 16 16 8
 4

Substituting in Eq.(11) we obtain

0.02
Ξιξ , v Ξκχ , v = κ2 R02 Bl,ιξκχ



!
X X
× a2r crs crs − ar as crs csr
rs rs
0.01
N

µ
dX
= κ2 R02 Bl,ιξκχ ar as H̃rs theory
κ rs
simulations
(15)
where we used that c2rs = crs csr = crs and the identities 0.00
PN 2 P P PN PN 0.0 0.1 0.2 0.3 0.4 0.5 0.6
r ar s crs − rs ar as crs = rs ar as [( j crj )δrs − z − 2d
d
PN
crs (1 − δrs )] = κ rs ar as H̃rs . Recalling that
PN


s H̃ a
rs s = λar , we obtain Ξ ιξ , v k Ξκχ , vk =
FIG. 1: (color online). The theoretical predictions for the
dκR02 λk Bl,ιξκχ .
shear modulus, Eq. (17), and the simulation results of
Ref.[12]. No fitting parameter is used in the comparison.
κ = 1, R0 = 1.
IV. RESULTS AND DISCUSSION

A. Elastic moduli
The predictions of Eq.(17) for the shear modulus, with-
out fitting parameters, can be compared with the simula-
Hence, we have shown that within the isotropic ap- tions of Ref.[12] of d = 3 disordered packings of (monodis-
proximation of the Hessian, Eq.(12), the nonaffine part perse) compressible spheres interacting via harmonic re-
of the elastic stiffness, Eq.(11), has the following thermo- pulsion in Fig.(1). In the simulations the spheres, at
dynamic limit T = 0, are slightly compressed to packing fractions φ just
N d above the so called jamming point at φJ = 0.64 which is
NA 1 X X dκR02 λq Bl,ιξκχ
hCιξκχ i= also an isostatic point with zJ = 2d [9, 12]. Further, the
V̊ q=1 λq jamming point is a zero-stress point [12], and the effect
l=1
(16) of stress on the global rigidity is therefore small. How-
d
N X
ever, it seems from this comparison, and from our results,
=d κR02 Bl,ιξκχ .
V̊ l=1
that stresses, in general, are not likely to affect the qual-
itative behavior of the global rigidity as they play no role
The affine part of the elastic constants for the disordered in arriving at the fundamental scaling law Eq.(18).
Cauchy solid can be obtained by performing the disor-
A Finally, we note that the vanishing of K at the isostatic
der average of Eq.(10), hCιξκχ i, where h.i denotes the
point predicted by our theory does not agree with the
angular average, and we always use the isotropic distri-
scaling of K observed in soft-sphere packings where it
bution of the bond orientations. In d = 3 we thus obtain
1 N remains finite at zJ = 2d [12, 14]. It agrees however
µA = hCxyxy
A
i = 30 2
V κzR0 , for the shear modulus, and with the behavior of random networks where K vanishes
K A = 13 (hCxxxx
A A
i + 2hCxxyy i) = 181 N 2
V κzR0 , for the bulk linearly at zJ = 2d [14]. The reason for this might be
modulus. Therefore, using these affine moduli together tentatively identified with the fact that in our theory, just
with the coefficients of Eq.(14) and with Eq.(16) we de- like in networks, excluded volume effects are irrelevant,
rive expressions for the shear and bulk modulus of the whereas they are important in packings [14].
d = 3 disordered Cauchy solid, respectively as
1 N
µ = µA − µN A = κR02 (z − 6)
30 V (17) B. Vibrational density of states and validity of the
1 N approximation
K = KA − KNA = κR02 (z − 6)
18 V
For d = 2 we obtain The isotropic Hessian matrix introduced for the calcu-
1 N lation of the nonaffine contribution to the elastic moduli
µ = µA − µN A = κR02 (z − 4) can be used to obtain the density of vibrational states
16 V (18)
5 N (DOS) numerically. Recall that the approximate Hessian
K = KA − KNA = κR02 (z − 4) is given by the following random matrix
48 V
Generalizing this result to arbitrary space dimensions  
gives the following scaling for the moduli in d dimensions αβ κ X
Hij = δij cij − (1 − δij )cij  δαβ
µ ∼ K ∼ (z − 2d) (19) 3 j
 5

The Hessian is defined by the coefficients cij which de-

pend on the realization σ and are, therefore, random
variables. The eigenvalues and thus the eigenvalue dis-
tribution of a random matrix are also random quantities.
The explicit calculation of the eigenvalues as functions
of the matrix elements is not possible. The approach to
the eigenvalue problem in random matrix theory makes

D(ω)
use of the self-averaging assumption that the eigenvalue
distribution becomes deterministic in the limit of an in-
finite system size [16]. As analytical solutions for the
eigenvalue distribution in our case are not possible (due
to the sparseness of the Hessian), we resort to a numer-
ical analysis assuming that the self-averaging property
holds. Therefore, we can define a limiting eigenvalue dis-
tribution ρ(λ) as follows
ω
lim hρN (λ)iσ = ρ(λ)
N −→∞

Then we have for all realizations σ that FIG. 2: (color online). 1: The DOS calculated for a system of
N = 10000 particles according to Eq. (1) (solid line) and z =
lim ρN (λ)[σ] = ρ(λ) . 8. 2: Simulations for repulsive unstressed harmonic packings
N −→∞
of [17] with φ − φc = 0.1 which corresponds to z − 6 ' 2.
Remark that for a finite N the eigenvalues distribution
is discrete and given by
N quite well reproduced, this was somewhat expected be-
N 1 X cause sacrificing information about the bond orientations
ρ (λ)[σ] = δ(λ − λi )
N i=1 does
p not change the single-particle parameters (of order
κ/ω) which dominate the highly (Anderson) localized
where δ is the Dirac delta function. In the limit N −→ high-ω modes. On the other hand, it is now clear from
∞, the set of eigenvalues has infinite elements, and the this comparison what details of the vibrational spectrum
distribution becomes continuous. The factor 1/N is nec- went lost in our isotropic approximation of the Hessian,
essary to normalizeRthe density given that the normal- Eq. (12): the spectrum calculated with the isotropic
∞
ization condition is 0 ρ(λ) dλ = 1. To analyze numer- Hessian is significantly depleted of low-frequency modes
ically the eigenvalue distribution we calculate the eigen- whose density is significantly underestimated in compar-
values sets {λi }r for large systems ( N ' 10000 ). Then ison with the DOS of the packings. This is probably re-
we create an histogram of the eigenvalues set and we fit it lated to the anisotropic character of the local correlations
with a continuous curve which approximates the limiting between particles in the packings (ultimately related to
eigenvalue distribution. The eigenvalue distribution is excluded-volume and static effects) which play an impor-
usually described in terms of the vibrational DOS which tant role in the low-frequency modes and that are lost in
we denote as D(ω), where ω is the vibrational frequency. the approximation. This observation hints again at the
The latter is related to the eigenvalue distribution by the possibility that overall our model describes random net-
change of variables works (where excluded-volume are absent) better than
r sphere packings. The relationship between these obser-
λ vations and the excluded-volume effects should be more
λ→ω= and its inverse ω → λ = mω 2
m systematically investigated in future studies. This was al-
ready noted in relation to the bulk modulus prediction of
Hence, with dω = 2√1mλ dλ and dλ = 2mω dω we get our theory which agrees indeed with the random-network
that the density distributions of λ and ω are related by: scaling (K ∼ z − 6), though not with the sphere packing
q one, for arguably similar reasons.
λ
D( m ) In view of these considerations, it is natural to ask
2
D(ω) = ρ(mω ) 2ω and ρ(λ) = √ why our theory, which seems to underestimate the low-ω
2 mλ
modes, still yields a correct prediction of the shear mod-
We compare the so obtained DOS from Eq. (12) with ulus for sphere packings. This question is related to the
the DOS from simulations of unstressed harmonic pack- issue of the role played by the low-frequency modes in the
ings [17] where φ − φc = 0.1 (φc ≈ 0.64 is the jamming nonaffine response. While this issue is a very open and
packing fraction), corresponding to z − 6 ' 2. The com- unsolved one in our current understanding of amorphous
parison is shown in Fig.2. While the upper end and the solids [10], a tentative, and certainly incomplete, answer
main features of the spectrum (width and average) are to this deep question is proposed in the next section.
 6

C. On the role of low-frequency modes in the In the case of the bulk modulus, simulations [8] give
nonaffine response practically the same behavior for the correlator Γιξκχ (ω)
as for the shear modulus, with the low-frequency modes
To assess the relative importance of different regimes contributing to the nonaffine response to a minor extent.
of the vibrational spectrum in the nonaffine elastic re- In this case, the failure of the theory in predicting the
sponse, it is instructive to rewrite the elastic moduli with correct scaling for sphere packings (despite being suc-
the nonaffine correction in the continuous frequency do- cessful for networks) is more likely to be ascribed to the
main. According to the nonaffine linear response formal- geometric attenuation of the random affine fields under
ism of Lemaitre and Maloney Ref. [8], in the thermody- hydrostatic pressure due to excluded volume, as we spec-
namic limit one has ulated in section IV.A. However, this hypothesis has to
Z ∞ be tested in future work by means of ad hoc numerical
A D(ω)Γιξκχ (ω) studies as the bulk modulus scaling of packings is a prob-
hCιξκχ i = hCιξκχ i − dω (20)
0 mω 2 lem currently under debate [14].
where the correlators on the frequency shells are defined
by
V. CONCLUSION
Γιξκχ (ω) = h(Ξιξ , vk )(Ξκχ , vk )iωk [ω,ω+dω] (21)
We have developed an approximate, fully analytical
The function Γιξκχ (ω) thus represents the projection of theory of the nonaffine elastic response of amorphous
the affine fields on the frequency shells and its magni- solids which explicitly takes into account the nonaffin-
tude gives the importance of the contribution of each fre- ity of the atomic displacements. We have applied the
quency shell to the nonaffine response. From Eq.(20) it is nonaffine linear formalism in the formulation of Lemaitre
evident first of all that in the zero-frequency limit ω → 0 and Maloney [8] to the so-called Cauchy bonded-network
the moduli diverge to minus infinity, i.e. hCιξκχ i → −∞, model [1], i.e. to networks of harmonic central-force
unless either D(ω = 0) = 0 or Γιξκχ (ω) = 0. At the iso- springs. In order to evaluate the nonaffine correction
static or jamming point of sphere packings one has that to the elastic moduli analytically, an approximation of
the DOS develops soft modes with D(ω = 0) 6= 0. As the the Hessian matrix has been proposed where the bond
nonaffine linear formalism is an exact theory, it is then orientation-dependent factors in the Hessian are replaced
strictly necessary that with their isotropic average (isotropic Hessian). Even
though the isotropic Hessian has a density of states
lim Γιξκχ (ω) = 0 (22) which significantly lacks low-frequency modes in com-
ω→0
parison with sphere packings, our approximation yields
i.e. the zero-frequency modes must not contribute to the predictions of the shear modulus in excellent quantita-
nonaffine response. This is what has been observed in- tive agreement with simulations of sphere packings [12].
deed in the numerical simulations of Ref. [8] where, in The good agreement is explained with the observation,
the case of a Lennard-Jones glass, the function Γιξκχ (ω) supported by simulations in the literature [8], that the
measured in the simulations goes to zero at ω = 0. Fur- low-frequency modes, underestimated by our approxi-
thermore, in the same simulation study [8], it was found mation, play a relatively minor role in the nonaffine re-
that Γιξκχ (ω) not only goes to zero at zero frequency, sponse, which is controlled by the upper part of the vi-
but is a monotonically growing function of ω in the en- brational spectrum (that is well reproduced by our the-
tire domain, such that it has significantly lower values ory). While our approximation is not suited to accurately
at low ω than in the middle and upper part of the spec- describe transport properties of disordered solids in the
trum where it reaches its maximum value (cfr. Fig.5 in low-connectivity and low-frequency limits [18], it seems
Ref.[8]). Hence, the simulation results of [8] indicate on the other hand successful in accurately describing the
that the contribution of low-frequency modes to the non- elastic response to shear of disordered solids. Further-
affine response is small whereas the leading contribution more, our theory provides a completely new insight into
comes from the high-frequency modes. This observation the linear vanishing of shear rigidity at the isostatic point
is also in agreement with physical intuition: the source (z = 2d) of disordered solids: this happens because the
of the nonaffine response is given by the projection of the nonaffine correction at the isostatic point becomes equal
affine fields on the eigenmodes which has a higher value in absolute value, but with opposite sign, to the affine
the more energetic the modes are. part of the shear modulus.
Based on these observations, one can conclude that Aknowledgements. A.Z. acknowledges financial sup-
the low-frequency modes play a relatively minor role in port by the Swiss National Science Foundation (project
the nonaffine response as compared to the high-frequency no. P BEZP 2 − 131153). Enlightening discussions with
modes. This explains why our theory, which underes- V. Vitelli, M. Warner, E. M. Terentjev, and R. Blumen-
timates the low-frequency modes in the case of sphere feld are gratefully acknowledged. We are very thankful
packings, still yields correct predictions for the shear to M. Morbidelli for generously supporting this project
modulus in excellent agreement with simulations (Fig.1). in the initial phase.
 7

[1] S. Alexander, Phys. Rep. 296, 65 (1998). 699 (1985).

[2] P. M. Goldbart, H. E. Castillo, and A. Zippelius, Adv. [12] C. S. O’Hern et al., Phys. Rev. E 68, 011306 (2003).
Phys. 45, 393 (1996). [13] H. A. Makse et al., Phys. Rev. Lett. 83, 5070 (1999);
[3] C. Heussinger and E. Frey, Phys. Rev. Lett. 96, 017802 H.A. Makse et al., Phys. Rev. E 70, 061302 (2004).
(2006). [14] W. G. Ellenbroek et al., EPL 87, 34004 (2009).
[4] A. Widmer-Cooper et al., Nat. Phys. 4, 711 (2008); E. [15] J.C. Phillips, in Rigidity Theory and Applications p.155,
Del Gado et al., Phys. Rev. Lett. 101, 095501 (2008). Eds. M.F. Thorpe and P.M. Duxbury (Kluwer Academic,
[5] P.W. Anderson, in Ill-Condensed Matter Les Houches New York, 1999); L.V. Kantorovich, Quantum Theory
Session XXXI, Eds. R. Balian, R. Maynard, G. Toulouse of the Solid State: An Introduction, p. 257-259 (Kluwer
(North-Holland, Amsterdam, 1979). Academic, Dordrecht, 2004); G. K. Batchelor and R. W.
[6] A. Tanguy, et al. Phys. Rev. B 66, 174205 (2002). O’Brien, Proc. R. Soc. London, Ser. A 355, 313 (1977).
[7] B. A. DiDonna and T. C. Lubensky, Phys. Rev. E 72, [16] I. M. Lifshitz, S. A. Gredeskul, and L. A. Pastur, Intro-
066619 (2005). duction to the theory of disordered systems (New York,
[8] A. Lemaitre and C. Maloney, J. Stat. Phys. 123, 415 Wiley, 1988).
(2006). [17] N. Xu, V. Vitelli, A. J. Liu, and S. Nagel, EPL 90, 56001
[9] M. van Hecke, J. Phys.: Condens. Matter 22, 033101 (2010).
(2010). [18] N. Xu et al., Phys. Rev. Lett. 102, 038001 (2009); V.
[10] M. Wyart, Ann. Phys. (Paris) 30, 1 (2005); W. Ellen- Vitelli et al., Phys. Rev. E 81, 021301 (2010).
broek et al., Phys. Rev. Lett. 97, 258001 (2006).
[11] J.C. Phillips and M.F. Thorpe, Sol. State Commun. 53,