'13. ‘3.

‘w
.n ‘u‘
' “ya A. ”hawfi
affix-
i 533- .
S: v “3' j
’ ‘ ‘3.

       

:14
v-L't

r23:
4‘.

    
   

 
 

‘ .
.. v
v:
‘ ‘
.

   

WIII iv!“ ,~
\.I .. .

y"
h"-
‘ .

 
  
       

    

  

. . , , .- I v.‘.- -..,-»;:_’-f.-.-,..',-:.;-';1np' .‘w'.‘:-'}..Iv
‘ ‘, I 9. ‘ .f. “1'”- “h l :2 ' ”19'?" H '-. "ill. ., ‘ .g__ .(I , \ ‘
v ~ .. . ~ ~ ~ l "r" ~.., '1 W»--' -, v . . . '
. ‘ .‘ w". ‘ ! w W ,. ,1. - —
¢ “~ WM“ ‘., ‘. .. .:n;.. 1... m.- III .|
"~ I 4 ‘ 1‘ 's . 1 H. . v .n“ ;I.ul,“"_ "I“ nil ' "mu m A". H
9.24- ‘ mmpf. ‘. * .w . ..."'.'.al' .-"'~‘Ma "W511.” _..' -

 

 

   
 
 

 

TH 5515 L1 3 may l
Michigan state
University l

 

f w

This is to certify that the

dissertation entitled

THEORY OF THE ELASTIC AND VIBRATIONAL

PROPERTIES OF CENTRAL FORCE RANDOM NETWORKS
presented by

Edward Joseph Garboczi
has been accepted towards fulfillment

ofthe requirements for

Ph. D. degree in Physics

WW4

Major profefl
M. F. Thorpe

 

 

DmeFebruary 20, 1985

 

MSU is an Affirmative Action/Equal Opportunity Institution 0 12771

 

MSU ‘

LIBRARIES
.—_,—-

 

RETURNING MATERIALS:
Place in book drop to
remove this checkout from
your record. FINES will
be charged if book is
returned after the date
stamped below.

 

 

 

 

 

 

THEORY OF THE ELASTIC AND
VIBRATIONAL PROPERTIES OF

CENTRAL FORCE RANDOM NETWORKS

BY

Edward Joseph Garboczi

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Physics and Astronomy

1985

ABSTRACT

THEORY OF THE ELASTIC AND
VIBRATIONAL PROPERTIES OF
CENTRAL FORCE RANDOM NETWORKS
BY

Edward Joseph Garboczi

Central force random networks which are generated by
randomly removing bonds on a lattice are simple models with
which to study the rigidity transition where the elastic
moduli of a system go to zero as a function of decreasing
mean connectivity. A simple constraint counting argument
gives a prediction for the critical threshold at called the
rigidity threshold, and f(p), the fraction of zero frequency
modes, where p is the fraction of bonds present. Chapter 1
presents evidence from three different networks that shows
these predictions to be in excellent agreement with numer-
ical simulations. Also it is shown that the elastic moduli
for nearest neighbor central forces scale linearly with p,
in excellent agreement with a simple effective medium theory.
Chapter 2 demonstrates how this coherent potential approx-
imation based effective medium theory also gives a quite good

description of the entire vibrational density of states for

these same networks. Chapter 3 extends the ideas of Chap-
ter 1 to systems with first and second neighbor central
force bonds. Evidence is presented from both numerical
simulation and effective medium theory that a fixed point
exists--the ratio of Cu/C..always converges to a constant
along a given track in phase space independent of its start-
ing value before removing bonds. The dependence of this
critical value along the rigid-floppy boundary in phase
space is derived under effective medium theory. It is
demonstrated that f for this system depends only on the
total number of first and second neighbor bonds present

and not on the details of the phase space track. It is shown
how these results can be extended to systems with central

forces of arbitrary range.

ACKNOWLEDGEMENTS

First of all I would like to thank my advisor Professor
M.F. Thorpe for all of his careful patient guidance over the
last five years. His knowledge of solid state physics and
his approach to research have been a real inspiration to me.
I wbuld also like to express my appreciation to Professor
J. Bass for a year of carefully taught solid state lectures,
to Professor S.D. Mahanti for a year of invaluable statis-
tical mechanics lectures and many helpful conversations, to
Professor W.W. Repko for not a few helpful conversations
about computational and mathematical topics, and to all three
of the above for serving on my guidance committee. I
gratefully acknowledge Dr. J. Parkinson's many contributions
of ideas and advice to the work described in Chapter 1.

H. He, Dr. D. Heys, Dr. D. Sahu, and Dr. J.M. Thomsen
have been_good friends and the source of much help through-
out my graduate career.

Many thanks are due to my son Joshua, whose birth
inspired me to work very hard, and who kindly allowed me to
convert his changing table back into a desk at night in
order to write this thesis. I would like also to express

my gratitude to my dear wife Alice for all her patience,

ii

encouragement, and assistance over these last four and one

half years.
And finally, my gratitude and praise go to the Lord
God, the infinite-personal Creator of the universe, Who
with His incredible skill has made it such a delight to study

His marvellously intricate creation.

iii

TABLE OF CONTENTS

List Of TaleSOIIOOOOOOOOOOOOOOIOOOOOOOOOOOOOOOOOOOOOOOOVi
List Of FigureSOOOOOOOOOO0.00.00.00.00.0.0.0.... ....... Vii

Chapter 1. Rigidity Percolation on Central Force
Elastic Networks................... ....... ...1_

Section 1.1. Introduction and Definition
Of PrOblem...OOOOOOOOOIOOOOOOOO0.0.0.01

Section 1.2. Mean Field Constraint Counting and
Fraction of Zero Frequency Modes......8

Section 1.3. Lattice Potentials and Elasticity
Theory for Central Force Crystals....14

Section 1.4. Triangular Net--Numerical
Simulation ResultSOIOOOO00......00000017

Section 1.5. FCC--Numerical Simulation Results....39
Section 1.6. BCC--Numerical Simulation Results....50

Section 1.7. Effective Medium Theory-~Static
MethOdOOOOOIOOIOOOOOOO0.0.0.0...00.0.57

Section 1.8. Effective Medium Theory--CPA
MethOdOOOOOOOOOIOOOOOII...0.0.00.00.0063

Section 1.9. Conclusions..........................67

Chapter 2. Effective Medium Theory Vibrational
Density of States for-Central Forces. ...... .69

Section 2.1. Introduction.........................69
Section 2.2. Negative Eigenvalue Method...........71
Section 2.3. CPA Equation.........................76
Section 2.4. Pure System..........................79

Section 2.5. Results..............................84

iv

Section 2.6. Conclusions..........................97

Chapter 3. Rigidity Percolation on Elastic Networks
with lst and 2nd Neighbor Central Forces....98

Section
Section

Section

Section
Section
Section

Appendix A.
Appendix B.

Appendix C.

3.1. IntIOductionOOOIOOOOOOOO0.00.00.00.0098
3.2. Details of the Model................100

3.3. Constraint Counting and Fraction of
Zero Frequency Modes................104

3.4. Effective Medium Theory.............107
3.5. Numerical Simulation Results........112
3.6. conCIuSionSOOOOOOOOOOOOOOO00.000.00.123

Trimming and Supertrimming Rules for
Central Force Networks............... ..... 126

Extrapolation Technique for Elastic
Energy RelaxationOOOOIOOOOOOO0.0.00.000000129

Effective Medium Theory of Percolation
on Central Force Elastic Networks.........131

Table
Table

Table.

Table

bUNI-J

LIST OF TABLES

Comparison of pc and fit.......................13
. 5

Critical values of p1 for tracks l-3.........105

YVL( ratios used for numerical simulations.113

Extrapolation vs. relaxation...... ........ ...130

vi

Figure

Figure
Figure
Figure

Figure
Figure

Figure
Figure

Figure
Figure.
Figure

Figure

Figure

Figure

Figure
Figure

Figure

9.

10.
ll.
12.

13.
14.

LIST OF FIGURES

A schematic view of the random substitution
of selenium atoms into amorphous
germaniumOOOOOCOOOOOOOOOOOOOOOOO0.0.0000......2

Network unit with no restoring force
against external strain.......................7

Constraint counting estimate for f(p) = the
fraction of zero frequency modes.............11

Rectangular unit cell for triangular
net (20x22)0......OOOIOOIOOOOOOOOOOOOOOO0.018

20 x 22 triangular cell with p = O.65........25

Network from Figure 5 after trimming.........26

.Examples of three types of supertrimmable

unitSIOOOOOOOOOOOOOOOOOOOOOCOOOOOIO0.0.0.0...27

Network from Figure 5 after trimming
and supertrmingIOOIOOOOOOOOOOOO0.0.0.00000029

P' vs. p for 20 x 22 triangular network......30
Illustration of the diode effect............32
Two bond elastic "network"..................33

C“ and C..vs. p for the triangular
networROOOOOOOOOOOOOOOOOOOOOOOOOOOOO0.00.0.036

f(p) vs. p for the triangular network.......37

P' vs. p for 500 and 2048 site FCC
networks...........OIOOOOOOOOOIOOOO0.00.0...42

C“ and C.,vs. p for FCC network.............44
C“ and K vs. p for FCC network..............45
f(p) vs. p for perturbed and unperturbed

108 site FCC networks.......................47

vii

Figure

Figure

Figure
Figure

Figure

Figure

Figure
Figure

Figure

Figure
Figure

Figure

Figure
Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

18.

19.

20.
21.
22.

23.

24.
25.
26.

27.
28.
29.

30.
31.
32.

33.

34.
35.

36.

37.
38.

f and Af vs. p for 108 site
FCC networks.................. .............. 49

P' vs. p for 432 and 2000 site BCC
networks.OOOOOOOOOOOOOOOOOO0.0.0....O. ..... .53

Cu and ngvs. p for BCC network.............54
Cn,and-K vs. p for BCC network....... ..... ..55

"Wrong bond" substitution in
triangular netOOOOOOOOOOOOIOOOOOOO. ...... .0058

Sparse storage indexing scheme for the
negative eigenvalue method..................75

Density of states for the triangular net....80
Dispersion relation for the triangular net..81

First Brillouin zone for the
triangular netOOOOOOOOOOOOOO0.00.00.00.0000082

0.85.00.00.00000.00.0...86

CPA vs. NEM for p
CPA vs. NEM forp=0.70.;00000000000000.0.088

Effect of supertrimming ondensity of

states for p = 0.70..... ...... ...... ..... ...89
CPA vs. NEM for p = 0.50. ................. ..91
CPA vs. NEM for p = 0.20 ..... ....... ....... .93

Real and imaginary parts of effective
force constant vs. energy...................95

Diagram of square net with first and
second neighbor springs....................101

Phase space for square net model...........102

Fraction of zero frequency modes for
traCRSIand BOOOIIOOIOOOOOOOOOOOOOOO00.0.0106

Elastic modulus ratio along critical
line (EMT)OOOOOOOOOOOOOOOOOO0.0.00.00000000110

C...andC..vs.pfortrack1................114

Cwand Cnvs. p for track 2................115

viii

Figure
Figure

Figure

Figure

Figure

Figure

39.
40.
41.

42.
43.

44.

Elastic moduli for tracks 1 and 2 ......... .116
C" / C.,flow diagram for track 1...........118

C. / C... flow diagram for track 1
(small system)....... ........ . ............. 119

C. / C., flow diagram for track 2. ......... 121

C. / C.,f1ow diagram for track 2
(small system)...COCOOCOOOOOOOOOOOO ...... .0122

Two types of connections between floppy
units and the rest of the network .......... 127

ix

Chapter 1. Ri idity‘Percolation on Central Force
E astic Networks

Section 1.1. Introduction and Definition of Problem

It is widely believed that the structure of amorphous
solids whose bonding is primarily covalent is well described
by some form of covalent random network (Zallen 1983). A
good example is amorphous germanium (a-Ge). Crystalline
germanium (c-Ge) has the diamond structure, where each atom
is tetrahedrally connected to four nearest neighbors, with
the nearest neighbor bond lengths and angles between pairs
of nearest neighbors bonds having specified values. A-Ce
has short range order similar to that of c-Ge, in that all
atoms have four nearest neighbors in a roughly tetrahedral
arrangement. However, the nearest neighbor bond lengths
and angles in a-Ge take on values in a narrow distribution
centered on their crystalline values. This occurs in such
a way that there is no long range order present. Since
it takes energy to stretch bonds and change bond angles,
a-Ge has a higher energy than c-Ge. A-Ge cannot be'
continuously deformed into c-Ge, however, since the two
structures are topologically inequivalent. C-Ge only has
rings with an even number of members, while a-Ge has rings

with odd and even numbers of members. This implies that

an energy barrier exists between the two states, and there-
fore a-Ge is a metastable state of solid germanium.

The structure of two component amorphous solids may also
be modelled with a covalent random network. An example of
this type would be GexSel_x. A selenium atom naturally
forms covalent bonds with two other selenium atoms, which
results in pure Se being made up of isolated polymer chains
where each atom has two nearest neighbor bonds. We may"
then think of producing a compound like GexSe1_x by starting
with pure a-Ge and randomly replacing Ge-Ge bonds with
Ge-Se-Ge bonds until the desired atomic fractions of x and
l-x are achieved. A schematic picture of this process is

shown in Figure 1.

Q
d

 

 

 

Figure 1. A schematic view of the random substitution of
selenium atoms into amorphous germanium.

It is then entirely natural to ask what effect adding
Se atoms has on the elastic properties of the material.

It is known experimentally that nearest neighbor bond stretch-
ing forces (which tend to maintain constant bond lengths)

are roughly four to five times stronger than the bond bending
forces (which tend to preserve angles between pairs of near-
est neighbor bonds) in these kinds of covalent random net-
works (for a review of the experimental evidence see

Zallen 1983). In the substitution shown in Figure 1, we are
replacing a strong bond stretching force with a relatively
weaker bond bending force. Thus we would expect the elastic
moduli of GexSe1_x to soften as x decreases from 1, or
equivalently, as (r) , the mean coordination, decreases
from four. The theory of this is worked out in Thorpe (1983)
hereafter designated as MFT.

The key result from MFT is a prediction for rp, the
critical value of the mean coordination.<r) where the elastic
moduli of the network go to zero. The elastic moduli can
actually soften all the way to zero even when the network
is still well connected because with only nearest neighbor
bond bending and stretching forces specified there can be
units in the network which do not contribute to the elastic
strength of the material. An example would be a selenium
chain of length 9 6 (Thorpe 1983). This transition we call
rigidity percolation as opposed to ordinary connectivity
percolation, since the network will still be well connected

but have zero elastic moduli. This is a result of only

specifying nearest neighbor bond stretching and bending
forces. If we included a dihedral angle force (very tiny
in real glasses) all connections become rigid and the
problem would become a geometrical percolation type problem,
where the elastic moduli become zero only when the system
is actually geometrically disconnected.

To experimentally achieve the desired range of x so that
one can study solids with values of‘(r) both above and below
rP is possible but rather difficult. One must prepare
enough samples to be able to cover the relevant range of x
properly while maintaining tight chemical control of the
value of x. The Ge-Se system is perhaps the best choice for
this work. One could put sulfur into germanium but the
sulfur tends to form eight-fold rings which would not de-
crease (r) in a random enough way, since we want spatial
homogeneity. Experiments are currently being planned to
try to test the rigidity percolation ideas, but there is no
good experimental evidence as of yet. These are not easy
experiments to perform, mostly due to the difficulty of
sample preparation.

Lacking any experimental data, we wished to computer-
generate covalent random networks with various values of
(r) and compute their elastic moduli via numerical simulation.
In this way we could attempt to verify the mean field
prediction for the value of rp given in MFT.

In doing numerical simulations, one must choose one's

model carefully not only for its physical content but more

importantly for its numerical tractability. It is possible
to computer-generate covalent random networks of the
GexSe1_x type (He 1983) which have reasonable bond length and
angle distortions. One randomly inserts two-coordinated
sites into a previously generated four-coordinated covalent
random network, and then relaxes with bond stretching and
bending forces until a local minimum of the potential is
found producing a metastable state. However, since elastic
energies are small deviations from the local minimum of

the potential, in order to do accurate computations of
elastic moduli one must know the energy of this minimum
with great accuracy. It turns out to be very difficult

to impossible to obtain that sort of accuracy on these
networks (He 1983).

We are therefore led to try to find a simpler model which
exhibits the basic physics we wish to explore. What is
needed? An adequate model must have the following prop-
erties: l) unstrained potential energy should be known
exactly, 2) mean coordination must be able to be decreased
without interfering with l), 3) rigidity percolation must
be able to take place, and 4) the cost of numerical
simulations must not exceed the limits of our computer
budget!

The class of models whose elastic properties we chose
to explore are central force random networks whose randomness
is of a special type. We start with a Bravais lattice, the

triangular net in two dimensions, the face-centered cubic

and body-centered cubic in three dimensions, where the only
interatomic forces are Hooke's law springs between nearest
neighbors. We then simultaneously reduce the mean coord-
ination and introduce randomness by randomly removing bonds.
The initial elastic energy is always zero, as the lattices
have no initial strains in them. The elastic moduli are then
found numerically as a function of p= fraction of bonds
present, and 5; the bond fraction where the elastic moduli
appear to go to zero, is compared with the mean field
prediction obtained from constraint counting arguments.

To satisfy 3), there are units in the network which,
while fully connected, make no contribution to the restoring
force under an external strain. A simple example of one of
these units can be seen in Figure 2. More complicated
examples will be shown and discussed in Section 1.4.
Therefore we claim that results for these central force
models will have meaning for the covalent random networks
which after all are more physically reasonable models.

Condition 4) was also met--with a few cents to spare!

Figure 2. Network unit with no restoring force against
external strain

Section 1.2. Mean Field Constraint Counting and Fraction
of Zero Frequency Modes

As was stated in Section 1.1 one can construct an esti-
mate for rp, the value of the mean coordination where the
elastic moduli vanish for a covalent random network, using
a mean field constraint counting argument. This can also
be done for the central force models, which are essentially
”covalent" though without any bond bending forces. A change
of notation should be noted here. While (r) , the mean
coordination, is the natural variable to use for the covalent
random networks, for our central force models where bonds are
being randomly removed it is more natural to track the

changing coordination through p = the fraction of bonds

present. The relationship between p and (r) is simple:

(1:) = zp

r = z .
p P

z = n.n. coordination
in the crystal

The constraint counting argument for a central force
network in d dimensions goes the following way, which I have
adapted from MFT. We define the following quantities for a

network with N particles:

Nc - the number of independent constraints
in the system Which must be satisfied
in order to have a zero frequency mode

Mo = the number of zero frequency vibrational
modes of the system

Nd = the total number of degrees of freedom
or modes

f a Mo / Nd = the fraction of vibrational modes
with zero frequency

To have a zero frequency mode, one must be able to make
some sort of collective motion of the N particles such that
all the bonds remain at their equilibrium length 10' because

the lattice potential is

V= 72 “ 2 (1X34); . (1)

(all Lands)

We then have that

M =Nd-Nc

O
(2)

f =,1 - Nc / Nd

One way of seeing the validity of equation (2) is to note
that Nc - the rank of the dynamical matrix D(R), so that the
nullspace of D(R) has just Nd - Nc elements. Perhaps the
easiest way to think of equation (2) is to consider the Nd
equations of motion for the network involving the variables

Q

ri . If we have Nc equations of constraint (each bond length

10

that must independently remain constant gives an equation
like ('?i - 33 )2 = 102 ), then we can define Nd - Nc

new variables?j , which are linearly independent and which
automatically satisfy the equations of constraint. We
therefore must have M0 = Nd I Nc zero energy modes involving
these variables.

For the central force case, with a lattice potential
like (1), Nc is just equal to the number of bonds which
independently must have their equilibrium length so that
all bonds will be unstretched and thus a zero frequency mode
made possible. Our mean field estimate will be that Nc is
just equal to the number of bonds actually present, szp,

so that f becomes
:=f<p)=1-zp/2d

.A graph of the mean field prediction for f(p) vs. p
is shown in Figure 3. When p = 0, meaning no bonds are
present, the system consists of N free particles so that all
modes have zero energy. As p increases past fit the mean field
prediction for f(p) would become negative. A negative value
for f is meaningless, however, since f counts numbers of
modes. This is evidence of the fact that our mean field
estimate of Nc is too high, that each bond does not give an
independent constraint. An obvious example in the triangular
net is the last bond needed to complete a hexagon. The

length of this bond, if all the other bond lengths are

11

 

1 _!Guimaiomenmmhwamnnne—

 

 

 

the

Figure 3. Constraint counting estimate for f(p)
fraction of zero frequency modes

12

specified, is not an independent quantity. It is fixed by
the other bonds' lengths. So under our mean field theory,
f(p) remains zero for p.) p' . We call p. the rigidity thresh-
old, and use its value for our mean field estimate of where
the elastic moduli of a central force network will go to
zero.

It should be noted that in reality, the elastic moduli
will be non-zero when f is non-zero as well. In other words,
the rigidity transition will take place when f is actually
finite. A network can be macroscopically glqlg but still
have localized floppy islands. This is the idea of rigidity
percolation . In MFT a simple model is solved exactly
showing that at the critical threshold f is finite and the
third derivative of f is discontinous, implying a third order
phase transition. The dotted line in Figure 3 shows an
estimate of a more realistic f(p).

Table 1 gives the values of 6' for the three lattices
studied. It also gives the known connectivity percolation
threshold pc for each lattice. In all cases p. is much
larger than pc. In fact, a similar mean field estimate
for pc (Kirkpatrick 1973) gives pc = 2 / z = p./ d . So the
constraint counting theory predicts that the central force
networks will have zero elastic moduli while being still

very much connected.

13

Table 1. Comparison of pc and p’

g a
Lattice as E
triangular 0.3473 0.6666
fcc 0.119 0.5000'
bcc 0.179 ' 0.7500

@ Values for pc taken from Zallen (1983)

14

Section 1.3. lattice Potentials and Elasticity Theory
for Central Force Crystals

Before we can understand the elastic moduli of a central
force random network, we must understand and calculate the
elastic moduli of the crystal from which it is generated.
The most general way to do this is to set up the dynamical
matrix D(R), calculate the energies of the long wavelength
vibrational modes and pick off the values of the elastic
moduli. In the remainder of this section D(R) will be
derived for a central force potential, and then in the
following sections the triangular net and body-centered
cubic (BCC) and face-centered cubic (FCC) lattices will be
looked at in detail.

We start with the potential V from equation (1) and

rewrite it focussing on the sites rather than the bonds:

V: V‘I‘Z (la'i’l ”(9)1 (4)
(4’!)

where (13') means that J! and {limust be nearest neighbors,
and 1i is the total displacement of the atom whose '
equilibrium position is )3 . The factor of 1‘ includes a
factor of 8 for a spring potential and a factor of k to take
care of double counting. .

db "” ’ A ,
Writing 7'“ 4‘3; and '1')!" “1 , equation (4) is

15

l
expanded to second order in the “1 (harmonic approximation).

We then have

v= is}: [(41.9 - (MW

(1!!) ’37!”

Making the double sum unrestricted and using a well known
transformation (Ashcroft and Mermin 1975), the potential

takes on the form

(5)

ad
where SAN is the alpha component of the unit vector from

I
J! to J! , and Rt is the vector from site 4 to one
of its 2 nearest neighbors. Fourier transforming, we then

obtain

d .... ”3;.“
may: :57, k.'(/-— e ")

(6)

m a mass of atoms on lattice

A
K. = unit vector from a site to one of its
nearest neighbors.

16

1,, = nearest neighbor bond length = 1,
making k also dimensionless

D(R) is the dynamical matrix, whose eigenvalues and
eigenvectors are respectively the normal mode energies and
atomic displacement patterns for the phonon spectrum. Since
we are considering only Bravais lattices, D(R) is a d x d
matrix. We can now go on to diagonalize D(R) for k -+ 0
and extract the values of the elastic moduli for various
lattices of interest, thus giving the starting values of

the elastic moduli for our.bond removing computations.

17

Section 1.4. Triangular Net--Numerical Simulation Results
Details of the Lattice

The triangular net was the first lattice on which we did
numerical simulations. It is a two dimensional Bravais
lattice with six nearest neighbors per site (2 = 6). One
can construct a rectangular unit cell for the triangular net,
a picture of which is shown in Figure 4. For the simulations,
one would like to pick a cell that is as square as possible,
so as not to have a preferred direction after cutting bonds.
All the elastic moduli simulations were performed on a 20 x 22
supercell, giving 440 sites. Taking the nearest neighbor
bond length to be one, this makes the linear dimensions of
the supercell 20 x 19, which deviates from a square by only

. \

5 %.

Elastic Moduli for the Crystal
The dynamical matrix D(R) for the triangular net when

k -* 0 is given by

an x’u ’/4 K: 3/9 4”! ‘9-

1
’/4 K): I} 3/t' It” J/‘t 4')

Solving the eigenvalue equation

- a."
[)(FDo a. - as an (7)

18

AAAAAAAAAAAAAAAAAAA
VAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVA
AVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAV
VAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVA
AVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAV
VAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVA
AVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAV
VAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVA
AVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAV
VAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVA
AVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAYAVAVAVAV
VAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVA
AVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAV
VAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVA‘VA
AWAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAV
VAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVA
AVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAV
VAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVA
AVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAV
VAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVA
AVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAV
VAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVAVA

VVVVVVVVVVVVVVVVVVV

Figure 4. Rectangular unit cell for triangular net (20 x 22)

19

we find the two solutions are

u": qd/IM kt
to,” 3‘/'P\ 4"

To identify the character of the sound waves associated with
these two frequencies, we plug 0.‘ and u: back into equation
(7). Inserting 00." the result is Tc x ii = 0 which implies
that vi = 9-!/ 8m (LONGITUDINAL). Inserting 0%: into
equation (7) we get 'R'UE = 0 which implies that

v2 = 3-‘/ 8m (TRANSVERSE). Since the direction of'i was

t

arbitrary, we therefore have proved that the triangular
net is isotropic; i.e., v1 and vt are the same in all
directions. This condition implies that there will only be
two independent elastic moduli, C» and 4., (see for example

Nye 1957). The value of these is obtained from the relations

vi a: G./€ V: = (”V/C I

(a = density of triangular net

which give the results

 

3“ at - 3 4
CL: 3 C9” "' _—
4 9
Numerical Simulation Techniques
There are two general methods by which one can numer-

ically compute the elastic moduli for a system of mass points

20

connected by various forces. One method emphasizes computing
the stress for a known strain, and thereby extracting the
values of the moduli. This method probably requires rigid
boundary conditions (Feng and Sen 1984). The other method

' depends on computing the elastic energy for a known strain,
which can be used with periodic boundary conditions. Periodic
boundary conditions almost always give better results for

a given system size compared to rigid boundary conditions.
The energy method with periodic conditions was used for all
computations of elastic moduli in this work. There is also
more than one way to handle periodic boundary conditions
computationally. The method used in this work was that

of having an extra boundary layer of "imaginary" sites surr-
ounding the 440 ”real" sites, where the imaginary.sites were
periodic images of real sites, thus maintaining the proper
environment for each real site. There is only a small cost
in extra computer memory involved in this technique, as
opposed to a somewhat larger savings in computational ease
and efficiency.

The general method of computing the elastic moduli via
the strain energy for a given bond fraction p is as follows:
one first initializes the coordinates of all the sites, then
bonds are removed randomly until a fraction p are left. Then,
depending on whether (a or Gm is being computed, the
coordinates of the sites are transformed corresponding to a
certain external strain.

To compute C; , one makes the transformation

21
X" (If Ear)“
Y" (I- 6,»!

where 6; , 65'are not both non-zero at the same time.
For the triangular net, Cu was computed as an average over
the x and y directions by taking the following combinations

Of 6; and Er :

(x) 5x =- 10‘5 57 = o

_ I (a)
(y) €.=o €y=105

Averaging over the x and y directions helps to get rid
of finite size noise, because for an infinite triangular
network which starts out isotropic, randomly removing bonds
should not break that isotropy. Averaging over x and y
helps to restore this symmetry which can be broken for a
cell of finite size.

To compute (an, , the necessary coordinate transformation

3
x: X+€

x” (9)
ylzy+ €y.x

5

where 6x a 10- and 57 = 0. In two dimensions one need

not average over direction, as the requirement that the

system have a net torque of zero guarantees the same result

5

as when 6, = o, 57 - 10" (Kittel 1967).

22

Whatever the applied strain is, one then proceeds to
allow each site to move to its equilibrium position under
the condition that its nearest neighbors' positions are fixed.
Of course, as a site's neighbors are themselves relaxed,
the original site will no longer be at the equilibrium
position, but as this process is repeated many times the
system as a whole will converge to its equilibrium config-
uration. It is important to note here that periodic
boundary conditions are maintained throughout this process
which reflect the new shape of the unit cell under the
given external strain. If the new shape of the unit cell
were not maintained via the periodic boundary conditions
the sites would just relax back out to their unstrained
positions. The relaxation process finds the equilibrium
configuration of the cell for the strained cell shape.

In other words we find the microcopic atomic configuration

which gives rise to the macroscopic elastic modulus.

 

By summing over all bonds we then obtain the final
relaxed elastic energy per unit volume U, from which the
values of G and Cw are extracted according to the

following relations:

1,6..6‘

C‘.
II

or

U = 8 (Walgz

depending on which kind of external strain was used.

23

By varying p one proceeds to map out (1 and ‘29 as
functions of p. When this process is done one then repeats
it by choosing the bonds to be removed with a different
set of random numbers. We used three sets of random numbers,
and averaged the elastic moduli over these three configur-
ations. This damps out finite size noise, giving a better

approximation to the thermodynamic limit.

There is one aspect of the numerical simulation technique
which must be mentioned in connection with rigidity percol-
ation. As was mentioned before, when the sites in the random
network are connected by only Hooke's law springs, certain
configurations of sites and bonds arise upon bond removal
which, while still quite thoroughly connected to the main
network, have no contribution to the final relaxed elastic
energy. The simplest example is the dangling bond which can
always stretch to its equilibrium length no matter where
the rest of the sites relax to. All dangling bonds can be
successively removed without affecting the final
relaxed strain energy.

A non-trivial unit of the type mentioned above is the
site with two bonds. By simple geometry one can see that
no matter where its two neighbors move to (for small strains),
such a site can always adjust its position so that its two
bonds have their equilibrium length. So all such sites,
along with their bonds, may be successively removed until none

remain without affecting the final relaxed elastic energy.

24

This has a major effect, since such sites, as opposed to
dangling bond sites, play a real role in the connectivity

of the network. This process of iteratively removing all

one and two coordinated sites we call trimming. Figure 5
shows a network with 65% of its bonds present after random
cutting. Figure 6 shows the same network after trimming.

One can clearly see that many bondshave been removed, even
bonds belonging to sites which originally were 3,4,5 or 6
coordinated. This is an effect of the nonlocal aspect of
trimming, that one can only tell whether a certain bond will
be trimmed by looking at other bonds more than one or two bond
lengths away. One might wonder whether pairs of colinear
bonds which exist in the triangular net can always be trimmed
as well, because they can buckle uner compression but are
rigid under tension. This is an important point, and will

be discussed fully under Diode Effect.

There are higher order units which are functionally just
generalizations of the one and two coordinated sites. An
example of each is shown in Figure 7 a) and b). These also
will not contribute to the final relaxed elastic energy.

The unit shown in Figure 7 c) is somewhat different. Any
cluster of sites which is only connected to the rest of the
network by three bonds can also relax fully, since such a
cluster always has three degrees of freedom, two translational
and one rotational, so the three constraints represented

by maintaining the attached bonds'equilibrium lengths can

always be achieved. .All of the above configurations are

'25

 

V VAV‘ VAVAVA A‘A AVAA
If; ‘A'AVA AVA'AVA VV
A‘VA VIAVAVIIA
A AVA VA'AVA‘V’ VAV VI A A
A‘ VAV AVAV V‘ A
. AH'AVAV'I VVAV VA“
VV VAVAVAVIIAV'A VA A‘V'A‘

92‘ 19' ' “at“
V9 9..

 

 

 

 

AVAVA AA
:VAV-V AVA A V
VA‘VAVXA‘+AVAV- AVA
éI-V. ‘V I AAV
VAVA I‘VAVA‘V AVA
AV VAVVA AA
)VA"A AVA A AA‘VI AV

1I9\ as: w‘ m ‘9»: .

“VAVV V'IAI

 

 

 

Figure 5. 20 x 22 triangular cell with p=0.65

26

A

   
 
 
 
 
 

   
  

A A ‘ A A A
VAVAVA A‘V VA A
IAVA‘A AVA'AVA
A A‘VA VIAVAVII

  
   
 

‘V AVAV V V - A
'AVAVII VAV VA V
AVAVAVIIAVI VA A‘VIA
AVAV A VIA'
A AV V A‘ VAV VA
A‘V

A AVA
VIA AV-V AVA ‘

1% II“"" ”I
am a s f”
V AVA AVA‘V A

VAV A“VAVAV'VA
AVAVA'AV A A VI. A
A I-V ‘ V ‘ A ‘

A AV VA‘VA ‘ A‘V A I

, V AVAV VII A

VV \V\

 

 

 

I
‘

 

Figure 6. Network from Figure 5 after trimming

27

VA
le

 

Figure 7. Examples of three types of supertrimmable units

28

removed by hand before relaxing the network by a process
we call supertrimming. The removal of such units sometimes
creates new one and two coordinated sites, which leads to
further trimming, which can lead to further supertrimming,
_ etc. These two processes are iterated sequentially until
no new units can be identified by these criteria. There are
certainly other units which could be removed as well, but we
have not been able to identify them. A complete set of rules
for all supertrimmable units which we have found is listed
in Appendix A. Figure 8 shows the same network as in Figure 5
but now supertrimmed as well as trimmed.

The effects of trimming and supertrimming are to greatly
accelerate the relaxation process and save computer time
by not wasting relaxation steps on sites which will not
contribute to the elastic energy anyway. One can also find.
a lower bOund for p. by finding the point at which the net-
works breaks up under the combined effects of trimming and
supertrimming. Figure 9 shows p. as a function of p, where
p. = the fraction of bonds left after trimming or trimming
plus supertrimming. This gives a lower bound for p'of
about 0.61. A more careful computation has been done, using
ten configurations of 60 x 72 triangular supercell. One
trims and supertrims until the network breaks up in all
directions. The value of p for which this occurs is the
lower bound for pt After averaging, one finds that

9.20.62.

29

VAVAVA.

AVA A AVA AVA
V VIAVAVII

A A A‘V VAV V

  

A‘ ‘V. A‘V
AVA
VIA VAV-V AVA
‘AVAV
'A-AVA

AV;A\\VAVAV V v /
AVAVA'AV A
VA I-V VA V \ VAA\
AA) VA‘VA ‘ A‘V
V AVAV VII

VV \V\

 

Figure 8. Network from Figure 5 after trimming and
supertrimming

30

 

 

 

 

8.8 “
as . ‘/
, /
P /
a4 - f
I jUPERTRIMMED
6.2 4 TRIMMED/‘I
I
x,
0 I I I . 7" ‘ r I I I
a

Figure 9. P’vs.

p for 20 x 22 triangular network

 

31

Diode Effect

As was mentioned previously, when one cuts bonds on the
triangular lattice, there can exist sites with two bonds
which are colinear. These kinds of sites introduce a non-
linearity into the system, since under compressive strain
this connection will buckle, providing no restoring force,
while under tension this connection does have a restoring
force and so is rigid. This point is illustrated in Figure
10. We call this problem the "diode effect." The difficulty
is summarized in the expression "You can't push things
with a rope!" A system with these kind of connections will
then have different elastic moduli depending on the sign of.
the strain, which is a violation of classical elasticity
theory. What must be done to remedy this situation is
perhaps best illustrated by considering a simple example.

Consider a two bond "system" as shown in Figure 11 with

a potential energy

v-.- $3 ~ (4.4.? + $80-$31

It!

We can now put a horizontal strain 5 on this system,
where 6) 0 means compression and 6 < 0 means tension.
The elastic "modulus" C which couples to this strain can be

calculated with a little algebra to give

Z-rp slate.
c= —~ , ,1
1.6! Co; {0. + 1%5115199)

32

_____O__. ummmeo

m 4— —-+ TENSION

_/O\_ conrlfmoN

Figure 10. Illustration of the diode effect

33

 

Figure 11. Two bond elastic "network"

34

Now, to make this example illustrate the situation in
the triangular net, we need to take the two limits/5 —b 0
and 9. —v 180°. There are two orders in which to take these

limits of course.

lim lim emu!) = o (10)
9401' [#0
lim lim emu?) = «I, (11)
‘fao BEN?

Equation (11) illustrates the diode effect. C will
be finite under tension and zero under compression. Equation
(10) is then the proper order in which to take the two
limits. It is clear how we are to do the calculation on the
triangular net. Taking ,3 = 0 gives us our central force
potential. Then we perturb the triangular lattice so that
the angle between all pairs of straight bonds is less than
180°. This perturbation must be of the order 67", so as
to allow all straight pairs of bonds to fully relax under
both tension and compression. As the strain-+ 0, so will
the perturbation, and we will get back the proper result
.for the triangular net.

‘Computationally, we use 6 = 10"5 4‘ 1, and do the
simulation by simply trimming all two coordinated sites and
leaving the rest of the lattice unperturbed. For small

strains this is equivalent to the ideal limiting process.

35

Numerical Simulation Results

The actual results for Cu and (w as a function of p
for the triangular net are presented in Figure 12. These
results are averaged over three configurations, where each
network was trimmed by a computer program then supertrimmed
by hand interactively. The straight line in the graphs is
the result of the effective medium theories described in
Sections 1.7 and 1.8. The agreement between theory and
"experiment" is clearly very good, much better than we usually
get for simple effective medium theories. The data is almost
exactly linear nearly all the way down to $1 The constraints
counting prediction for p‘ = 2 / 3 is very accurate. The small
tail past p‘ has contributions from finite size effects and
incomplete relaxation in this region. Near fifi each connected~
site in the network had to be relaxed many thousands of times
to even approach convergence. A confirmation of our result
has recently come from Tremblay (1984). He and his co-
workers used finite-size scaling to obtain a value of
0.649 for 83 with a critical exponent of 1.4 i 0.2, which
is consistent with our results.

Figure 13 presents the numerical result for f(p). f(p)
was computed for three configurations of a 12 x 14 triangular
network by setting up the full 336 x 336 dynamical matrix
and numerically diagonalizing it, then counting how many
modes had zero frequency. These networks were small because
of central memory limitations in the MSU CYBER 750, and

were untrimmed and unperturbed. The difference that a lattice

HHO

550

Figure 12.

1.4

1.2

8.6

8.4

a.2

8.6

8.5

8.4

8.3

8.2

8.1

36

(~46 SITE TRIM PET-1HREE (ION?IC~$UI~2fiTIO|"1-_L

d

 

 

 

I I I I I

  

 

       

l I
8 8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8
-_—440 SITE TRIM PET—THREE COI‘FIGURRTIO‘E—

I

8.9

 

 

 

 

a 8.1 9.2 9.3 c.4'a.s 8.6 8.7 6.3 (3.9

Cu and C.” vs. p for the triangular network

 

37

1aasn31mnmcufinrenamesnnazcasuansnawn

 

 

 

 

Figure 13. f(p) vs. p for the triangular network

 

38

perturbation like that described earlier makes in f(p) will
be shown explicitly in the FCC results.. The difference

is very small. The numerical result for f(p) matches the
mean field prediction very well indeed, as can be seen by
comparing Figures 3 and 14. The slight curvature around
3:; 0.65 is due partly to finite size effects and partly

to real critical effects, as we do not expect mean field

theory to be exact (see Thorpe 1983).

39

Section 1.5. FCC Numerical Simulation Results
Details of the lattice

One of the two lattices in three dimensions upon which
we have tested the constraint counting theory is the face
centered cubic (FCC) Bravais lattice. Every site in this
lattice has twelve nearest neighbors (z = 12). One.can
think of the FCC lattice as being made up of four inter-
penetrating simple cubic sublattices, where each site's
twelve nearest neighbors are on a different sublattice from
itself. The nearest neighbor bond lengths are all a /423
where a is the edge length of the standard cubic unit cell.
For computational purposes it was convenient to take a =12:

so that the n.n. bond lengths all became equal to unity.

Elastic Moduli for Crystal
Using the dynamical matrix formalism for central forces
set up in Section 1.3, and letting k-¢ 0, we can write out

the elements of D(k) for the FCC lattice:

‘ a. z 1
[xur:'EE? {<k1f‘}v> Lay: if} 18:“?
D”: 3‘(k‘fxyv Dre 3 E—ng’sz

t
D22: 357656;) Dr" '55- 9’3

40

I will now take the mass m = l as this was what was done in
the numerical simulations. The FCC density then becomes

6’ = 4 / a3, since there are four sites per standard cubic
unit cell. D(k) is particularly easy to diagonalize in the
(100), (110), and (111) directions. We can find vi and vi,
the longitudinal and transverse sound velocities in these
three directions and then find C, ’, (n- , and (09 using

a result from Kittel (1967) relating the moduli and sound
velocities. For cubic symmetry there are only three inde-_

pendent elastic moduli which it is standard practice to take

to be C. , C. , and Cu, . For the FCC lattice we find that

C... ‘ Qd/A Cu: “/4 (lie: “/4

Letting Va =V-2- and“=II 1 to match up with the units used in the
numerical simulations, the final form of the elastic moduli

for the FCC crystal is

(IRIS.- (n: (99" Q's/2-

Q

The equality of a»: and (it. comes from a Cauchy relation.
For a Cauchy relation to hold, all sites must have inversion
symmetry and the forces between sites must be purely central

(Love 1944).

Numerical:Simulation'Results

All the elastic moduli computations for the FCC networks

 

41

were done on a cubic cell which was made up of 125 standard
cubic unit cells (5 x 5 x 5), giving 500 sites altogether,

with 3000 bonds present initially. Periodic boundary
conditions were maintained throughout the computation

in a manner similar to that used on the triangular net.

The values of the bond fraction p ranged from 1.0 to 0.46.
C). and as were averaged over the x,y, and 2 directions
for each configuration, and these results plus the bulk
modulus K were averaged over three configurations. (n. was
obtained by using the computed values of G, and K in the

relation

k:'3L{C.n+24a.) .

The diode effect is present in FCC as well as in the
triangular net. Here, however, one cannot escape dealing
with the problem directly as up to four bonds attached to one
site may all lie in the same plane, thus giving different
elastic response for in-plane tension or compression. This
difficulty is resolved by perturbing the FCC lattice through
moving each sublattice so that all bonds become-non-planar

5).

by an amount proportional to €Y‘( 5 = applied strain = 10-
In effect, each site is "puckered" slightly out of each of
the three cubic faces of which it was formerly a member.
This is a small perturbation, and only appreciably affects the
elastic moduli very near the critical region.

Trimming and supertrimming are both possible for FCC.

fiznmm v

nsz‘O '0

8.6

8.4

8.6

8.4

8.2

Figure 14.

42

 

 

588 SITE FCC CaL-F'IVE CONFIGURRTIOP&

 

 

I I V r j I I

62 sea me 3:: as ad' as

 

 

2848 SITE FCC CELL—FIVE COWIGLRRTIQ‘S

 

6.9

 

P’vs. p for

500 and 2048 site FCC networks

 

 

43

In three dimensions any site with three bonds or less will

not contribute to the final relaxed elastic energy so can

be removed. There are certainly higher order units which

do not contribute to the strain energy, as in two dimensions,

but these we have not identified because we cannot draw three

dimensional pictures. It fortunately turns out that super-

trimming appears to be of lesser importance in three dim-

ensions, as a 500 site network is broken up by trimming alone

on the average at about p = 0.43, which is already quite

close to {2' = 0.50. It would seem that ordinary trimming

is more effective in three dimensions than in two, if we

compare Figure 14 with Figure 9. Figure 14 shows the small

dependence of trimming on the size of the network. For

2048 sites, the trimming lower bound for p.is about 0.44.
Figures 15 and 16 present the results of the numerical

computations of (a. , G... , Co: , and K. For all cases, the

p dependence of the elastic moduli is very close to linear

all the way down to p‘ =0.50, with a small tail which is

partly a finite size'rounding of the transition and partly

a real critical effect. It is very difficult to completely

relax the network for values of p close to 65 so that

incomplete relaxation could also have made a contribution

to the tail. In fact, for pé 0.64, an extrapolation

scheme was used to extend the relaxation procedure. The

details of this scheme can be found in Appendix B. The

straight lines drawn in Figures 15 and 16 are the result

for Co" (p) obtained from the effective medium theories

44

 

588 SITE FCC CELL-THREE CONFIGURATIONSee

8 8 -
O

tar-O

8.4

8.2 ‘

 

 

a I r I I :‘I

 

I I I j
8 8.1 8.2 8.3 8.4 8.5 8.6 8.? 8.8
8.8

 

S88 SITE FCC CELLr-THREE CONFIGURRTIONS

8.7 d

8.5

550

8.4 ‘

6.3 .1

8.1 .

 

 

 

 

 

Figure 15. C“ and C44 vs. p for FCC network

45

 

5% SITE FCC CELL—“REE CG‘FIGURRTIOI‘S—

8.5

NPO

8.4 "

8.3 -

8.2 "

901 "

 

 

 

I I I T I I I

 

I I

8 8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 8.9 1

1.2 M SITE FCC CHI—THEE CWIGRRTIM

 

8.8 .

8.4

802 "'

 

 

 

8 8.1 8.2 8.3 8.4 8.5 8.6 8.? 8.8 8.9 1

Figure 16. Cr: and K vs. p for FCC network

 

46

described in Sections 1.7 and 1.8. The agreement of num-
erical simulation with theory is seen to be very good indeed,
and the agreement of the mean field estimate for p'with the
computational results is also quite good. It is interesting
to note that the near perfect linearity of GW'and Chvall the
way into the critical region implies that the Cauchy relation
(It’c‘fll is also obeyed with high accuracy over a wide range
of p. It is not known why this should be so, as inversion
symmetry at each site is lost as bonds are cut.

Figure 17 shows the numerical result for the fraction
of zero frequency modes f for both perturbed and unperturbed
networks. These curves were obtained by setting up the
full 3N x 3N real space dynamical matrix for the FCC network
at_various values of p, then diagonalizing it and counting
the actual number of zero frequency modes. The three Gold-
stone modes are subtracted out at the beginning. Because
of central memory limitations on the MSU CYBER 750 computer,
this computation was done for a 108 site cubic cell (3x3x3),
and averaged over three configurations. The constraint
counting prediction for f(p) agrees quite well with the
numerical simulations. There are only small differences at
the 2% level around p' due in part to finite size effects
and actual critical effects, since we do not expect mean
field theory to be perfectly accurate.

Figure 18 at the top shows a blow-up of the region
in which f(p) for the perturbed and unperturbed lattices
differs at all. The bottom graph shows fun(p) - fper(p).

47

 

 

 

 

 

Figure 17. f(p) vs. p for perturbed and unperturbed 108
site FCC networks‘

48

The maximum difference comes at about of The unperturbed
f is always larger than the perturbed f since for central
forces, out of the plane motion of a site with coplanar
bonds has zero frequency. This symmetry is broken for
the perturbed lattice giving this kind of motion a small

finite frequency.

MOZI'IflI'I‘I'IHU

Figure 18.

8.2

8. 175

8.15

8.125

8.1

.875

0%

. 815

fand

49

W—FRRCT ION 0F ZERO FREQUENCY NODES—FCC 188 SITES

    
  

. -- . . U
PERTURBED

 
 

o
o
‘
.Q

 

.--..O.--.§“
.-\‘ .....
“‘ ‘.--“--

 

 

I I V I I

8.4 8.45 8.5

8.55 8.6 8.65

W—FRQCTION OF ZERO FREQUEI‘CY ”DOES—FCC 188 SITE}.

 

 

8.7

 

 

A f vs. p for 108 site FCC networks

 

50

Section 1.6. BCC Numerical Simulation Results
Details of the Lattice

Numerical simulations to compute Cc, (p) were also
done on the body centered cubic (BCC) lattice. The BCC
lattice can be thought of as consisting of two inter-
penetrating simple cubic sublattices. A site on sublattice
A has all its eight nearest neighbors (z = 8) on sublattice '
B and vice versa. The nearest neighbor bond lengths are all
‘43'a / 2, where a = edge length of the standard cubic
unit cell. For computational purposes it was convenient to
take a = 2 /45'so that the nearest neighbor bond lengths

were normalized to one.

Elastic Moduli for the Crystal
Following the same procedure as was done for FCC, we

write out in the k - 0 limit 062):

z z.
Dxxgbflebaz: fk a

l

Dtv’bxr ’§' “'7‘
DI} ,DIY: g-eriat
' 2

on ' D3! ’5" K7,!“

We diagonalize D(k) in the (100), (110), and (111) directions

and calculate the elastic moduli similarly as for FCC. This

gives the somewhat surprising result:

51

CI,=C1= (”qr-3:1;-

When we let a = 2 /43 to rescale the nearest neighbor bond

lengths we obtain:

(“II2 ('1: (”Ve?‘

The equality of Ca. and GM is due to the existence
of a Cauchy relation (Love 1944). C. andCtbeing equal is
a consequence of the BCC lattice symmetry. A curious result
of this second degeneracy is that vt2’ the transverse sound
velocity in the (110) direction where the atomic displacements
are in the plane, is identically zero. In fact, one can
easily show by considering D(k) for.a general-k that one
of the three branches of 405(k) is identically zero in the
(110) direction. Thus even with all bonds present the BCC
lattice still has non-trivial zero frequency modes. The
number of these is insignificant in the thermodynamic limit,
but for the small systems used in the simulations their
fraction was large enough to make the determination of f(p)
meaningless. It was therefore pointless to compute f(p),

but the behavior of the elastic moduli was still of interest.

Numerical Simulation Results

The cell used for the BCC simulations contained 216

standard cubic cells (6 x 6 x 6) for a total of 432 sites.

52

There were 1728 bonds present when p = 1.0. Periodic bound-
ary conditions were maintained, and the values of p invest-
igated ranged from 1.0 to 0.67. Regular trimming was done
in all cases, with effectiveness comparable to that in the
FCC networks. Regular trimming caused the breakup of the
elastic backbone on the average at about p = 0.64. The
effect cell size had on this number can be seen in Figure 19
where p. vs. p is graphed for both a 432 site and a 2000 site
cell. P' a the fraction of bonds left after trimming.

All other numerical techniques were the same as those used
for the FCC networks.

The diode effect is present in the BCC networks. There
are planes of four bonds attached to one site which will
have a different elastic response depending on whether the
in-plane local strain is tension or compression. Such bond
configurations can persist after trimming. We therefore
perturb the lattice by shifting the two sublattices with
respect to each other by an amount proportional to 3%,

where 6 - the applied strain - 10‘s. This breaks up all

the planes of bonds, assuring the same result for the elastic
moduli for i 6 . As in the case of FCC, this lattice
perturbation has its only significant effect near the critical
region.

Figures 20 and 21 present the results for CR , 62,,
Ch, , and K. The results are comparable in quality to those
for the triangular net. The straight line on the graphs

come from the effective medium theories described in Sections

53

 

 

 

1 T_———-432 SITE BCC CELL-18 CGTIGLRQTIO‘IS
as ' ,,’
P x”
006 CI ’l”
P x”
R 3”
I ."O
N ”I.
E ,x
604 ‘I '3"
8.2 " '0’..
a ' I r I r I T7 r I I

 

8 8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 8.9

 

 

 

 

 

1 2000 SITE BCC CELL—SEVEN CON'IGURRTIONS
8.8 ‘
p
8.6 ..
P
R
I
H
E
8.4 "'
902 "'
8 I I I I I

 

 

8 8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 9.9

Figure 19. P’vs. p for 432 and 2000 site BCC networks

 

54

 

8.8 432 SITE BCC CELL-THREE CONFIGURRTIONS

8.7
8.6 -

8.5

PJ‘F)

8.3 .1

8.1 -

 

 

 

a r r I I I I "i; I I
8 8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.3 8.9
8.8

 

432 SITE BCC CELLr-THREE CONFIGURHTIONS__

8.7
8.6 q

805

Jb&(3

8.4 4

08.3 -

 

 

 

 

 

8 8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 8.9 1

Figure 20. C” and CW vs. p for BCC network

55

8.8

 

432 SITE BCC CELL-THREE CONFIGURRTIONSe

8.7
8.6 -

8.5

NHO

8.3 -

 

 

 

O
I I I r I I I I I

8 8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 8.9
8.8

 

432 SITE BCC CELL-THREE CONFIGURRTIONS

a? f
as -

9.5

8.4 ‘

8.3

8.2 7

8.1 -

 

£AA°

 

 

 

8 8.1 8.2 8.3 8.4 8.5 8.6 8.? 8.8

Figure 21. CI; and K vs. p for BCC network

56

1.7 and 1.8. The agreement is quite good between theory

and experiment, though not quite as good as for FCC. The
constraint counting prediction for p'is 0.75. This prediction
agrees quite well with the data, although the size of the
critical tail for BCC is larger than that for FCC. As

noted before, the tail is no doubt due to finite size rounding
of the transition and relaxation difficulties in the critical
region.- The Cauchy relation Cm‘cn. and the 4.3 (n. relation
are maintained with a high degree of accuracy all the way

into the critical region. The reason for the persistence of
this behavior over a large range of p when random bond

cutting should have long ago destroyed the summetries on

which these degeneracies depend is unclear at the present

' time.

57

Section 1.7. Effective Medium Theory--Static Method

There are two ways that we know of to develop an effect-
ive medium theory that gives a functional form for the
elastic moduli as a function of p, and that gives a predic-
tion for of The method described in this section we call the
static method. It is based on an effective medium theory
(EMT) used by Kirkpatrick (1973) to describe the dependence
of the conductivity of a resistor network on the fraction of
resistors randomly present. This static method was developed
by S. Feng, and is described in Feng (1984), which has been

included in this thesis as Appendix C.

Let us start with a triangular net, with all interatomic
force constants equal to 9‘ . Now apply a uniform stress

to the system, which can be hydrostatic, uniaxial, etc.

For sake of definiteness let us consider a hydrostatic
stress, so that all bonds have the same relative displacement
uh. Now let us focus on two sites, labelled l and 2, where

the spring constant between 1 and 2 has been changed to Y'
with Y“< . (This last assumption is not needed. Y can
have any value. Note that in the numerical simulations,

Y' a 0.) After this substitution, the spring between 1 andiz
will compress by an additional amount <84 , since it is not
strong enough to resist the rest of the lattice and maintain
a relative displacement of only uh. This scenario is depicted
in Figure 22. Now if we were to apply an external force F

to move sites 1 and 2 back to their original positions, F

58

 

 

 

 

 

.5

Figure 22. "Wrong bond” substitution in triangular net

59

would have to be such that

F+ Yuh:q’uk

or F=/¢’-T)“A

The external force F "stiffens" the weak spring Y. . Now
suppose we were to apply the same F to sites 1 and 2 without
the uniform stress being applied to the system first. This
force would induce a relative displacement dflv between 1

and 2. By the principle of superposition,
Jul: (in.

The relation between F and J“ can be obtained in the
following manner. When all the spring constants are fit and
we apply F to sites 1 and 2, we can define an effective spring

constant between 1 and 2 by

“effiF/J“

We can write “eff as "eff 3 °(/ ace , where O l. acen 4 l,

n
since at eff will be greater than " because all the other

indirect connections between 1 and 2 will give a positive
contribution to 0‘ e ff. If we then remove the 4 spring
between 1 and 2 and replace it with a r spring, :4

eff
becomes

oi == 0‘ - at
eff / acen + r

60

Equating the two expressions for F,

F = uh (‘“ - Y )

(12)
F= «eff Ju

we have that

Juan

 

(d/acen-x+ Y)

J“ can be thought of as the ”fluctuation" in the relative
displacement of sites 1 and 2 that resulted from changing the
spring constant from 0‘ to Y . The key idea of the eff-
ective medium theory enters the argument at this point. Let
us renormalize a< to ~&., and try to choose GA in such a way
so that a network where all bonds are d; has the same elastic
properties as the original network with 1' bonds substituted
in for some of the x bonds. The obvious way to choose flu.
would be so that Jh is zero for every bond. This is clearly
not possible. However, we can try to choose fl; so that
the average value (In) of {a is zero. The average is

defined by the following formula:

61

(Mn) =5 {Pmmwx

where P(Y)= probability
distribution for
values of r‘.

This gives the usual kind of effective medium or mean field

result, with fluctuations treated in only an average sort of

way. Setting ‘(43) = 0 gives the following equation:
d"/‘tea.' ‘4‘ " Y‘ (13)

where we note that 0‘ has been'replaced by 4':- . The

 

probability distribution of interest for percolative pro-

cesses is

P(H= “Hr-a) + (I—,)J‘(r)

Evaluating (l3) and simplifying gives the result

4... = d (1’: “"1-
I’ ‘6!“

Since we know that all the elastic moduli scale with

 

at , this shows that they all will be linear in p, and go to

zero at p - acen‘ The value of acen can be calculated

exactly (see Appendix C) and turns out to agree with the

62

constraints counting result:

2d / z

0)
II
ll

'6

cen

So the final result for the elastic moduli is

CLvW6£)='C:9 (1) ”-JP‘
I- r'

63

Section 1.8. Effective Medium Theory--CPA Method

The coherent potential approximation (CPA) effective
medium theory gives the same result for C;i(p) and p’as does
the static method. The CPA gives more than this, however.
As will be seen below, the CPA method gives a self-consistent
equation for the effective medium density of states at all .
frequencies. This effective medium theory was worked out
by M.F. Thorpe. The presentation which follows is taken
from Feng (1984). This paper may be consulted in Appendix C
of this thesis ( see Elliot, Krumhansl, and Leath for a
general reference on the CPA).

One starts with the Hamiltonian

HBHO+V

where
H= 211: * {‘5’ [IE-:3.)-fi-,.]z
.: "" 03)
and V a fifth!) [(E-ZlaJI

O

Hois the Hamiltonian of the system with all bonds present, and
V contains the effect of changing the force constant on the
bond between sites 1 and 2 from V to Y . We define the

Green's function for HD by:

‘55.. 2: gamma» __ @FIWIE'Q

‘1 1‘ U‘Un‘tw. u’W’Wo

64

H
where {“03 is a complete set of energy eigenstates. G

13'
is defined similarly for the full Hamiltonian. It should
fill.
be noted here that i, j run over sites, and‘aij and Pij
are d x d matrices for each pair i, j. V can be rewritten

in matrix form as

H 4 4
Vij a (r“) r.‘f.‘ "xi

where

my ’ (6.!" J5 * a}; 6;,- - J}; {‘7‘ - J”- J,»

.0
It can be shown that .5, .13, and V satisfy the Dyson equation
which is of the form

..g on.”

(5 -3 + pvc (14)

We are only considering a single defect, so that summing over
all repeated scatterings from this same bond we can rewrite

(14) as

a
G=P+PTP

where '71" is called the "t-matrix" (for obvious reasons) and

is given by

"' Y-x 4 ..
Tij .1: runs ”3‘
a a 7

1'2(F‘)fit(&‘n);n

65

G-.

fin.
It should be noted that Tij has the same matrix form as Vij
I

but with a different overall multiplicative factor. The
heart of the coherent potential approximation now is
setting the average t-matrix equal to zero, allowing a< to

be renormalized-to 0!. , the effective medium value:

<: Y-4Mm( :>"- C)
1- 2(r—.z...)1:..(‘r:.--‘r;)an (15)

 

The physical idea here is that we have somehow eliminated the
leading order corrections to our effective medium, defined

by (15). The average is again defined as in Section 1.7:

(Mn) -- [f(r) Amalr

where P") is the probability distribution for the values
of Y . Using the percolative distribution we obtain the
result

f(d-G’n) + (I'f) (“d-Q = 0

i-z/«-e.>-=..(‘r:-‘i.>v.. 1 . u. :..zz-m. as)

We can use the equation of motion for P11, given by

a".
nut?“ = i + 2—1‘ gr; (80"0)?l$
J (17)

66

to simplify (16) . It should be noted that d-* w... in the
definition of‘ng as well as in the average t-matrix equa-

tion. Combining (16) and (17) we obtain finally

a [P 2'5 ("“‘fi"’)]"’~ [1‘ g/"J’V'fl ' 0 (18)

If we let tfit -¢ 0, one can readily see that (18) becomes

(m: d P-" 9‘ 21
l‘f‘ ) f ‘3

 

the~same result as obtained in Section 1.7. Equation (18)

serves as the starting point for the analysis of Chapter 2.

67

Section 1.9. Conclusions

 

The original purpose of the work that has been described
in Chapter 1 was to attempt to verify the predictions of a
simple constraints counting argument for fit the point at
which the elastic moduli of a central force random network
vanish. The data presented in Sections 1.4, 1.5, and 1.6
show clearly that this simple argument works remarkably
well. An extra bonus in this work has been the excellent
agreement between the numerical results for the elastic
moduli and effective medium theory.
I should note at this point that numerical simulations
for the bulk and shear modulus of the triangular net and
FCC networks were first done by Feng and Sen (1984). Using
a potential that was strictly harmonic, they obtained the
result 5' a 0.58. Their value for p. is less than ours because
in their model pairs of straight bonds could contribute
to the elastic energy and so the network could hold together
elastically for lower values of p. They essentially did a
different problem than we did because their potential was
the harmonic approximation to ours. Our work was carried out
independently from Feng and Sen.
The application of the constraints counting argument
‘to more realistic models of glasses now seems hopeful. The
(zentral force models, while very simple and unrealistic.
have the essential physics of the glasses: decreasing
aVerage connectivity produces units in a random network that

‘3<> not contribute to the elastic energy under an applied

68.
external strain. Therefore rigidity percolation, of which
these are an example, is a different class of problem from
ordinary percolation. Just how different remains to be
seen as work progresses. A more realistic model of a
glass has been analyzed using numerical simulations, and
the results agree quite well with the constraints counting

ideas (He 1984). The central force models provide a base for

this and further work in this rapidly developing area.

Chapter 2. Effective Medium Theory Vibrational Density of
States for Central Forces

Section 2.1. Introduction

 

The success of effective medium theory in describing the
behavior of the elastic moduli in central force random networks
led us to see if other applications would be as successful.

The values of the elastic moduli essentially give us the small
I01 behavior of the vibrational density of states (see last
section of Appendix C). The CPA-based effective medium theory
described in Section 1.8 can give the density of states for

all out under this approximation. The so-called static method
of Section 1.7 can also be modified to give this information
(Feng 1984 and Appendix C) but it is simpler to use the already
existing CPA result. This method will be described in

Section 2.3.

The numerical "experimental data" for the vibrational
density of states 3(6):) with which the CPA theory is to be‘
compared has been obtained using the negative eigenvalue

method (Bell 1972, 1976', 1982: Dean 1972), modified for use
With large randomly sparse matrices. This method gives 3(03 )
V3. «1‘ in histogram form. The accuracy of the histogram

is llimited only by the size of the system and the amount of

comPuter time one wishes to invest. This method will be

69

70

described in Section 2.2.

All the work in this chapter was done on the triangular
net rather than on the BCC or FCC networks. The reasons for
this choice were the following: 1) we can get effectively
larger unit cells in two dimensions thus reducing finite
size effects, 2) the size of the dynamical matrix is only
2N x 2N instead of 3N x 3N, where N= the number of sites,

3) the number of non-zero elements in the dynamical matrix
which actually have to be stored is proportional to the near-
est neighbor coordination, and 4) the integrals over the first
Brillouin zone which must be done in solving the CPA equation
are much cheaper to compute in two as opposed to three

dimensions.

71

Section 2.2. Negative Eigenvalue Method

The following brief description of the standard negative
eigenvalue method is drawn from Bell (1982).

We start from the usual eigenvalue equation for LU; ,

the phonon frequencies:

2".

D.szw u,

where‘Q is the real space dynamical matrix (sometimes called

the force constant matrix):

~33; "’ Nil“?
4? 4 g A, . . .
Dfii' 5J5). .5). ”(.7
‘1
0 o‘HICrwfse, (1)

$3.: unflvcchr ("n A: h 7.

All masses are equal to unity. We can rewrite (l) as

(2)

72

The eigenvalues of A will be equal to those of kaut shifted
downward by an amount -uf’. Since the eigenvalues of Q are
all non-negative, this implies that the only eigenvalues of
A which are negative will be those corresponding to eigen-
values of 2 less than wt . So the number of vibrational
modes with frequency less than :03 ‘will just be equal to the
number of negative eigenvalues M(A) of A. To compute this
number without actually having to diagonalize A_is made poss-
ible by the negative eigenvalue theorem.

This theorem in its simplest form relates M(A) to M(AQ
and M(§;, where A,and A,are two smaller matrices related to a

partition of A, We partition A_(symmetric because 2 is symm-

etric) in the following way:

A. Q.
' (3)

9* 9.

where A’is n x n, AJis m x m, §,is (n-m) x (n-m), and 9,18
n x (n-m). The result of the negative eigenvalue theorem

is that

M(A)= MAJMMz)

(4)

where

73

We can continue this process with 5,,eventually giving the
result

1
N(A)=S WA.)

1"

where we have iterated the partitioning process 1 times.

It is most convenient to choose the next A;to be 1 x 1, so that

H(Ax)31 "f ,4120
f” (Ab): C7 )2; ’41‘2‘9

Processing the matrix A_in this way, after n steps (if A
is n x n) we will have computed M(A). -

In realizing this process on the MSU CYBER 750 computer,
one quickly runs into the problem of limited central memory.
In symmetric storage mode, which requires N(N+l)/2 words of
storage for an N x N matrix, N must be at most about 400 or so.
For the triangular net this means one can only compute the
density of states for a system of 200 sites or less. By way
of contrast, the results for the elastic moduli were obtained
on a 440 site network, more than twice as large. A fairly
large number of sites is crucial for these density of states
computations in order to see the details of the spectrum.

It happily turns out that the matrix Q’is very sparse
because we have nearest neighbor forces only and therefore
so is 5, Less happily, this sparseness is not predictable in

any useful way. One must then store the non-zero elements in

74

a one dimensional array, with an associated array giving
their position in A according to some numbering scheme. The
numbering arrangement used in this work is illustrated in
Figure 23 for a system with 3 sites (N=3) which gives a 6 x 6
dynamical matrix. There are two entries in A per neighbor
plus two self terms for each degree of freedom, which gives

14 entries per row and there are 2N rows. We only need to
store roughly half of the entries because A,is a symmetric
matrix. Therefore the initial storage requirement is approx-
imately 14N for an N site triangular net. This figure of
course decreases as we cut bonds. This storage requirement
only increases by about 50% at most during a typical computation
so that values of N up to 440 can easily be handled. The
storage requirement does not go up nearly as much as one might
expect, because although each Aghas some new non-zero elements
in it, the size of the Age is decreasing as well. These
computations are quite long, however, and require a lot of
central processor time. But they are possible on the MSU
CYBER 750 system. In this method, one essentially trades

the "hard" limitation of central memory availability for

the ”soft" limitation of computing time (and cost).

75

3 4 5 6

 

8 ‘1 10 11

lZ 13 14 15'

 

 

 

E] 16. 17 15

 

 

 

19 2.0

(:1 I 21

ID

I I U I \7 W
I
E]

 

 

Figure 23. Sparse storage indexing scheme for the
negative eigenvalue method

 

76

Section 2.3. -CPA Equation
The starting point for the calculation of the CPA
density of states, 9‘(u}), is equation (18) of Section 1.8,

which is reproduced here as equation (5):

P+ P' («Nu-I) " 4’5. [1 * RNA-1)] : 0

where p- fraction of bonds present
m.= effective force constant, which is
complex in general
to a phonon frequency squared
P, =- complex Green's function for the
effective medium (all bonds present
with force constant 4. )

Equation (5) really consists of two equations, one for the
real and one for the imaginary part, since R. is complex
and IV... will acquire an imaginary part in general. P" is

complex according to the prescription

 

R./w‘="‘?: "

N‘ I")! 63- «4:13)

P“R s ,6“, Kg 8.0.3596?!)

‘pea

P“: 3 A,“ In Pu /““x7}

‘40

(6)
I
9:60"): 17!; go [(01)

'chr the triangular net d=2 and z=6, so that equation (5) becomes

77

9+ 31./m...» - «..[1 +321wa = a

What does equation (7) mean? Equation (7) is a pair of coupled
equations in the variables 9/...“ and V: . They are self-
consistent in the sense that we cannot explicitly solve for
or,“ and v.1 in terms of the parameters of the problem.
This is because R: and ..I depend on the real and imaginary
parts of V». in a very complicated way, involving sums over
the first Brillouin zone in k-space. Therefore we must solve
them by iteration at every p and. 69'. In the usual way with
the CPA method, on. becomes energy dependent.

The procedure for solving equation (7) is the following.
Starting with a trial solution for d: and V..." , we calculate
Pu‘ and PHI . This is done by changing the k-space sum
of equation (6) into an integral then performing the integration
by Gaussian quadratures. The small imaginary part 7' of the
energy is taken to be zero, since for p # 1.0, at: ‘will be
non-zero so that no extra convergence factor is needed. The
imaginary part‘of 47.. will broaden the delta functions
in the integral in equation (6) into Lorentzians which can be
evaluated numerically. One then feeds the values of n‘ and
R}; into equations (7a) and (7b), and solves this pair of
linear equations simultaneously for the new values of 0’: and
(9/...I . Equations (7) are linear since Pu“ and 85 are
taken to be parameters of the equations. One then begins

again with these new values, iterating the entire process

78

until of. and Pg. are seen to converge reasonably well.
This typically took 5-12 iterations depending on the energy
and on the initial guess. The final calculation of 3.1/00.)

is saved, with the CPA density of states given by
it’d): .i'L PHI‘W‘) '

79

Section 2.4. Pure System

A few comments about the properties of the triangular
net with all bonds present are appropriate as background for
the p(l.0 results presented in the next section.

For the harmonic lattice potential given in Section 1.3,

we can solve for 67?) for the triangular net, with the result:

(07(3): ng i IVE—J

Fc3- (askxa- Zeos U44 Con’gk’a

(8)

c= («w we» may. MW we“

One should recall that m-l and a.=1 for the pure net, while
a(.. replaces 0( when “(a appears in the effective medium
Green's functions used in the CPA method. The maximum frequency
squared is equal to 6‘95 .
Given this dispersion relation (8), one can easily gen-
erate the density of states 360’) by randomly picking points
in k-space and counting how many fall into a given range of
a}. . If enough points are selected (typically about 100,000)
the result converges to the correct value. Figure 24 shows
the density of states for the triangular net which was gener-
ated by this procedure. The four discontinuities in slope
labelled on the graph refer to Figure 25, which shows the
dispersion relation (8) graphed along the first Brillouin
zone path FkA-B-r'defined in Figure 26. Features (1), (3),
z

and (4) clearly arise from local extrema in lat/F). At U =4.5

the lower branch (obtained by always taking the (-) sign

80

 

a (w’)

 

 

 

 

 

Figure 24. Density of states for the triangular net

81

 

 

 

 

 

 

 

r A' I <3 . F
WAVEVECTOR

Figure 25. Dispersion relation for the triangular net

82

 

 

..__-____._.\ A

 

Figure 26. First Brillouin zone for the triangular net

83

in (8)) has its maximum, not explicitly shown in Figure 25,
which gives rise to feature (2).

The fact that 3(6) -0 constant as wt --)0 is perhaps
a little disconcerting at first if one is used to looking at
a three dimensional density of states, but one can easily

prove that in d dimensions.

5/60") M cod-” M (ML-’0

so that 300‘) ~ constant as «9" -v0 for the triangular net.
Figure 24 then is our starting point. As we remove bonds
we will be able to see how g/w‘) evolves from its pure

system form.

84

Section 2.5. Results
Choice of values ofjp

Four values of p, the fraction of nearest neighbor bonds
present on the triangular network, have been selected at which
to compare the CPA density of states 5,. (00‘) to the negative
eigenvalue method "exact experimental" results. These four
values of p are: l) p=0.85, 2) p=0.70, 3) p=0.50, 4) p=0.20.
The first value was chosen to be well above the critical region
and close enough to the crystal so that the CPA should do well
there if it does well anywhere, and yet far enough from the
pure system so that some real differences from the crystal
might be seen. Point 2) was chosen to be just above p‘to give
a stringent test of the accuracy of the CPA in the critical
region.’ Also, trimming and supertrimming have a fairly large
effect at p=0.70, so we desired to see how the density of
states divided itself between the rigid "backbone" and super-
trimmable pieces which together make up the network. The
third value of p investigated is interesting because when 50%
of the bonds are present, the triangular network is still
connected geometrically but is disconnected elastically.
It should be noted that p=0.50 is less than even Feng and Sen's
(1984) result for 550.58, so that there is no question about
the elastic moduli being zero or not. And finally, point 4)
is below p‘for regular connectivity percolation. The network

consists of isolated clusters only at p=0.20.

85
1) p=0.85

Figure 27 shows the two results for the density of states
of the triangular network. The dashed line is the CPA result
while the histogram with area normalized to l is the negative
eigenvalue method (NEM) "experimental " result. It should be
noted here that all the NEM results are an average over three
configurations of a triangular network with 440 sites, in fact
the same three used in Section 1.4 for the elastic moduli
computations. The agreement between NEM and CPA in Figure 27
at the level of accuracy of the histogram is quite good. One
can see the memory of the crystal in the peaks which remain

at ae‘

- 2 and 5, while the effect of bond cutting shows in the
increase in the small frequency density of states, due to the

fact that

gas) ~ (4'. 6.2% (p-r'f' a4 ..2. o

All these features including the decrease of the right band
edge are reproduced well by the CPA. The integrated weight

under the CPA curve is equal to 1, as it should be.

2)¥p-0.70

The density of states results for p=0.70 are shown in
Figure 28. The agreement between CPA and NEM (untrimmed) is
similar to that in the previous figure. The CPA has a 3251
narrow peak that goes up off the graph to about 1.9 and then
comes back down to the correct €01 =0 EMT result. The bin
size used in the NEM was not small enough to check if this

behavior was real or an artifact of the CPA. To get higher

DENSITY OF STATES

Figure 27.

86

 

(L3

9
a:

P
d

 

p = 0.85

 

CPA vs. NEM for p = 0.85

 

 

 

 

87

accuracy with the NEM, one should go to smaller bin sizes
algng’with larger networks so as to keep the noise level tol-
erable in each bin. The average number of modes pg; bin

is really what determines the noise in the histogram. In
Figure 28 one can clearly see the beginning of the divergence
of g(0) as the elastic moduli go to zero. The right band
edge continues to come in from the pure system value of 6.0
with the CPA tracking this behavior quite well.

Figure 29 gives the NEM results for the p=0.70 networks
as well as these same networks trimmed and supertrimmed. The
supertrimmed histogram is normalized to the fraction of gitgg
present in the remaining backbone (about 0.85) with all the
.zero frequency modes from the missing sites subtracted out.
The untrimmed histogram's integrated weight is normalized
to l. The solid lines are the contribution from the "rigid"
backbone while the dashed lines show the contribution from the
supertrimmed parts. The effect of supertrimming seems to be
to deplete 3(6)” by about the same amount across all frequencies
except for a larger amount at low frequency. This is perhaps
not so surprising when we recall that supertrimming can remove
sites with all possible connectivities (one to six) so that
the complicated pieces of the network removed could support
the whole range of allowable frequencies. As a final comment,
one should recall the statement made in Section 1.4 that the
supertrimming rules we formulated were not complete, but were
‘probably almost so. This was seen clearly to be true in compar-

ing the NEM data for the p=0.70 networks. After supertrimming,

88

 

 

 

 

 

 

1.2 l I l I r l
I
p=OJO
1r
{3
:2 (L8
'5
u'os
O
Em:
m
E
0 Q2
L l
(h 1 2 3 4 5 6
w2

Figure 28. CPA vs. NEM for p = 0.70

89

 

common 0" TRIMED MD WRIH‘ED DOS
0.3

a? -5
as #1

6.5

UIOO

0.4 7

 

as J

 

 

 

 

Figure 29. Effect of supertrimming on density of states
for p I 0.70

90

the networks still had a few (one or two) non-trivial zero
frequency modes left. Therefore supertrimming, as formulated,
does not identify all the possible floppy regions in the

network.

3)4p-0.50

At p=0.50, the network is connected geometrically but
disconnected elastically. Figure 30 presents the CPA and NEM
results for this value of p. Again, the overall shape of the
NEM density of states is reproduced by the CPA to the level
of accuracy of the histogram.' The one systematic deviation
is at the right band edge where the CPA cuts off before the
actual band edge. The CPA shows a gap opening up at small
frequencies, with the new band edge at about to" ==0.05.
The histogram does not have fine enough resolution to see this.
It is very difficult to solve the CPA equation for small ‘0‘ ,
since in this region, to get the band gap, the solution has
to jump from the Riemann sheet that was valid above the trans-
ition to a new solution sheet in the gap. It is easier to see
the small gap in the p=0.20 results. At p=0.50, the value of

f, the fraction of zero frequency modes, should be
{(P):1- P/P' '3 0.25—

TPherefore the weight in the band should be 0.75. The weight
llnder the CPA density of states comes out to this value, and

the NEM results agree. The contribution of the zero frequency

91_

 

 

 

 

 

0.5
0.4
a:
' 1'3
.5 0.3
:0
II-
o
::(12
a
z
m 0.1
o
0
0
wz

Figure 30. CPA vs. NEM for p = 0.50

 

 

92

modes to the density of states will be a delta function at the
origin with weight 0.25. This is not shown in Figure 30.
However the CPA does predict this delta function accurately,

as will be shown in the p = 0.20 results.

4) p = 0.20

At p = 0.20, the network consists of isolated clusters of
sites and bonds, with the most probable clusters containing
one or two bonds. This is reflected in the NEM results shown
in Figure 31. The large bins at u" = 2.0, 1.5, and 2.5
come from one and two bond clusters respectively. These are
isolated cluster modes, and obviously a simple effective med-
ium theory cannot be expected to get these right. However the
CPA does not do too badly in reproducing the other features
of the density of states. There is a band gap at about

h;' = 0.3 which the CPA puts at ‘0‘ = 0.5. The low density

of states in the first bin probably shows a Lifshitz(l964) tail,
which the real infinite system must have. This is an expo-
nential tail which comes from the band edge into zero. The
right band edge has moved in to about hf’ = 4.3 which the CPA
puts at 4.0. The overall shape of the density of states is
reproduced fairly well by the CPA. At p = 0.20 the weight at
zero frequency is 0.7, so that the weight in the band is 0.3.
The weight under the CPA and NEM curves is equal to this value
within computational error.

The CPA below 13‘ predicts a delta function at the origin
with weight 1 - p /‘ p' (see Section 1.2) . This comes about in

the following way. Assume "w"; —to for p¢p’. Then eq. (5) is

0.25

.0
:0

0.15

DENSITY OF STATES
'o 9
a. .s

Figure 31.

93

 

 

 

 

 

 

paoao

l

 

 

 

CPA vs. NEM for p = 0.20

or P“: (I- P/F:)-;!:

When we take the imaginary part of P" to get a, (w‘) , we

3.0%“ 3%: Im W (.1— ’7’“) “w"

get

which agrees with constraint counting. In Figure 32 the real
and imaginary parts of 'TL are graphed as functions of 09‘
for p-O.85 and 0.20. For p=0.20 it is clearly seen that bgth
01,: and '5 go to zero with w‘ , so that the CPA does predict
the right delta function at the origin for p¢pf The graph of

«n for p-O.85 shows the real part going to

0.55 : Ll. /
,:¢Ltr

which is the correct result as was seen in Chapter 1. One
should also note on both the p=0.85 and p=0.20 graphs that

or: is zero outside the band. This is because the only

way 3J0.” can be non-zero is if "I is non-zero, and that
will occur only when «g, has a finite imaginary part.

Other features of interest in Figure 32 are the cusp-like
behavior of d: (p=0.20) for small frequencies. This is an
indication of the change of solution sheets going from in the
band to below the band, as was mentioned previously. The

solution

95

 

(18
(LB
(L4

 

 

 

 

 

 

Figure 32. Real and imaginary parts of effective force
constant vs. energy

96

or,“ /w‘=0) = Cr.
l‘f‘

remains a solution below p.as well as above, but below p7
it is not the correct solution as it gives elastic behavior
at small frequencies. For p¢pt at small frequencies one must
carefully search numerically for the other solution which is
the correct one in this region. The less prominent cusp-like
feature in s.“ and «.3 (p=0.20) at 00‘ =4.o is probably at least
partially a numerical artifact as it was difficult to solve

the CPA equation at the right band edge for plpt

Section 2.6. Conclusions

It has been demonstrated that not only does a simple
effective medium theory give a good description of the zero
frequency elastic properties of nearest neighbor central
force random networks, but it also gives a quite good descrip-
tion of the finite frequency behavior as well. This is true
over a range of values of p where the network varies from
' being geometrically and elastically connected to being geo-
metrically connected and elastically disconnected and finally
becoming completely disconnected. Band edges and the zero
frequency delta functions are also reproduced by the CPA
as well as the overall shape of the density of states.

In the next chapter the results of Chapter 1 are extended
to longer range central forces with again a spectacular
success of effective medium theory. The finite frequency
behavior of these systems has not yet been investigated, but
it is interesting to speculate that the CPA equation, suitably

generalized, might be quite successful in that case as well.

Chapter 3. R1 idity Percolation on Elastic Networks with
Ist and 2nd Neighbor Central Forces

Section 3.1. Introduction

In Chapter 1 rigidity percolation was studied on various
random networks which had nearest neighbor central forces only.
Effective medium theory and constraint counting arguments
described the elastic moduli and rigidity threshold very well
indeed for these systems. It was of interest to us to see
if the same kind of results would be obtained for systems with
longer range central forces. Perhaps the success of our theory
was due only to the nearest neighbor character of the atomic
interactions in our model. Also, since more realistic models
of covalent glasses have angle forces, higher neighbor central
forces are interesting too because in the right context, a
second neighbor central force is not so much different from
a nearest neighbor angle force.

The second thrust of this work was to explore some ideas
that have appeared recently concerning the universal critical
behavior of the ratio of C"/(qq as p—wpf Bergman (1984) and
Schwartz (1984) have shown that in certain models, C"/Cw
approaches a limiting value as the critical threshold is
approached that does Egg depend on the initial value of the
ratio.

Section 3.2 presents the details of the model we studied,

98

99
Section 3.3 describes the constraint counting arguments plus
numerical results for the fraction of zero frequency modes
as a function of p for this model, Section 3.4 derives the
effective medium theory we used, Section 3.5 presents the results
of numerical simulations for the elastic moduli and
ratio and compares these to the effective medium theory,

and Section 3.6 draws conclusions from this work.

100

Section 3.2. Details of the Model

The model we studied was a square net with first and second
neighbor Hooke's law springs. The first neighbor force con-
stant was 01 and the second neighbor force constant was Y'.
Figure 33 shows a piece of this lattice with all bonds present.
The solid lines are first neighbor bonds and the dashed lines
are second neighbor bonds. The full potential for this

model is

Vargaz)‘ + i rfaj—ir

(1)

where the sums over i and j respectively are over all first
and second neighbor bonds present. First neighbor bonds are
present with probability p.and second neighbor bonds with

probability pf It should be noted that p'and p‘are independent

 

probabilities. The phase space for this system is then the
rp‘plane, with 05p,p51.0. Figure 34 is a diagram of this
phase space. The lines labelled 1,2, and 3 are the three
phase space tracks actually used in the various simulations.

The three tracks are defined as follows:

Track 1: gfp'
Track 2: pf2.37prl.37 (2)
Track 3: pf2.37p.

The reasons for the choice of these three particular tracks

are given in Section 3.5.

101

 

Figure 33. Diagram of square net with first and second
neighbor springs

102

 

(L4

(L2

 

 

 

 

Figure 34. Phase space for square net model

103

Using the methods described in Chapter 1 we can easily

show that the elastic moduli for the pure system are:

C” a: N 4' Y
Cg”: (1;: r (3)
Since we still have only central forces, and each point is

a center of inversion, (“I and (a. are related by a Cauchy

relation (Love 1944) equating them.

104

Section 3.3. Constraint Counting and Fraction of Zero
Frequency_Modes

As in Chapter 1 the fraction of zero frequency modes f

1c: Mo/dN .- dAl- xv,
4A!

is given by

 

 

j. _ W+M3)

div
F 1-(fa+F‘) (4)

where z, z‘are the first and second neighbor coordinations

in the pure system and are both equal to four in this case,
and N,and N,are the number of first and second neighbor bonds
actually present. The critical line in the p.-p‘ phase plane

is given by f=0 or

P."+P.’=1

This result is shown as the solid line in Figure 34. Note that
equation (4) predicts that f depends only on the sum of p and
ptand not on their individual values. In fact, if we use
P'8(g+g), the actual fraction of 511 bonds present, we can

rewrite (3) as

FOO): I‘ZP

We have done numerical simulations to test the validity of (3).
Using the negative eigenvalue method as described in Section
2.2 we have selected values of g and p’along tracks 1 and 3

and computed f by counting the number of modes with frequencies

105

less than some very small number x (x:0.00001 ). These
results were averaged over three configurations and graphed
vs. p. The result for f along tracks 1 and 3 is shown in
Figure 35. The straight line is the prediction from equation
(5). Note that the agreement between tracks 1 and 3 is almost
perfect except for some very small deviations in the critical
region around fi;0.50, and both separately agree quite well
with equation (5). As at least some of the curvature in the
critical region is due to finite size effects, we can conclude
that the fraction of zero frequency modes is dependent only
on the sum 303! with only very small fluctuations in the
critical region from individual variations in p and p‘.

One final note: the critical line being defined as in (4)
implies that fi'for each of the three tracks will be as given
in Table 2.

Table 2. Critical values of p? for tracks l-3

' Track 2':
l 0.50
2 0.703

3 0.703

106

 

 

 

 

1 | I l
a max #1
°'° ' o TRACK #3
0.0 -
-
0.4 ..
OJIP
0

 

 

0 0.25 0.5 0.75

Figure 35. Fraction of zero frequency modes for
tracks 1 and 3

 

107

Section 3.4. Effective Medium Theory

We use the static method described in Section 1.7 to
develop an effective medium theory for our model. Since we
have two independent effective force constants, d;,and Y;.,
we end up with two equations like equation (13) in 1.7:

Vn‘ d {ff—4')

'. a| (63)

r..= r (w <6!»

1"at

 

where a'and alare defined similarly as in Appendix C:

L... 5 am :“W/fi) 6770]

N2, in? We)
a): 2r. 2’ 7,171-; ‘25 "U04” 0 07)] (7b)

”28. “1”

A

J; - lst neighbor unit vectors

3 - 2nd neighbor unit vectors

a - lst neighbor distance = 1

b a 2nd neighbor distance =V2'

D(k) is the inverse of the dynamical matrix of our model but
with at and Y replaced by d. and {a . The dynamical

:matrix D(k) is easily shown to be

D" (if? 2“. (I‘lu’xa) J- 2 Y,“ (l-(uha 011;.)
017(t)‘: Dy:[[)‘ 2 C. {”1614 {0.164,

108

.Dy7 (F) ‘-‘ )4. (1' (0513.)). 2 Y,“ {/‘toer‘ c.4174)

The elastic moduli go to zero at fiéaIand éEaz, and since the

following sum rule is easily proved from (7)

that implies that

fi‘i- ff: 1

which agrees with the constraint counting result of Section 3.3.
The existence of the sum rule (8) means that only a,need be
independently computed. The definition of a,in (7) involves
a lst Brillouin zone summation which can be changed to an
integral and evaluated numerically. This double integral
cannot be evaluated analytically but instead was computed with
a 22 x 22 point Gaussian quadrature. The explicit form of
the integral for a,is not important but it should be noted
that afa'( 5%..) only. ,

To solve equation (6) along a given track, we divide
(6b) by (6a) and solve for p, eliminating ptby the relation
(2) of interest and azvia the sum rule (8). This gives an

equation of the form

[9,: A “Va/J

109

This is a ngnrself consistent equation for pr We can substitute
in values of 'T/n., starting with the initial value of {Ad
and continuing until the known critical value of ’3/2. is
reached. Each of these points give a value of p which can
then be substituted back into equations (6) along the
calculated value of a.to calculate “L and r; thus the
elastic moduli via (3). The value of YLAZ. at the critical
point is known, since on the critical line figaf Thus for a
given track p? is known and so
X" l = a" 0”)
06-.

m1“: [nu
Figure 36 shows how the value of 60/6" varies along the
critical line. The elastic modulus ratio is simply related

to YF/k. by

C"/(.,., =' l * (53')

All of the above theory can be easily generalized to
account for an arbitrary number of Hooke's law forces of arbi-
trary range. If there are zibonds from a given site to its
ith neighbors, each randomly present with probability g,
and if we have forces out to the nth neighbor, then equations

(6) will be generalized to

d1: dfl.(,"aa') 24:: I’K

 

110

 

 

 

 

 

 

Figure 36. Elastic modulus ratio along critical line
(EMT)

111

ay -' '35 Z Tr[(1-¢""‘?'33')/?.Z°)'157W]
Nix A53.

Sum rule (8) can be generalized to

I!

2' 2,0,: 20!

4°31
so that the critical surface becomes defined by
A e
2: {cf} : 24’
in
This also the result which would be obtained from constraint

counting in this general case.

112

Section 3.5. Numerical Simulation Results

We present results for the elastic moduli computed along
tracks 1 and 2 in the p1 - p2 phase plane. The results are
displayed as a function of p = 15(p1 + p2), the total fraction
of bonds present. The rigidity threshold is then p”3 0.50
for both tracks. All the simulations were done on 40 x 40
square nets with periodic boundary conditions. Five inde-
pendent configurations were averaged over for each track.
Regular trimming of one and two bond sites was performed,
but supertrimming, though possible, was not attempted. The
regular trimming took care of the "diode effect" as it did
for the triangular net in Section 1.4. All other details
of the simulations were as described in Chapter 1.

Three different values of ‘Db were used for each track.
These are given in Table 3 along with the related initial
value and critical value of (”/(‘W .

Figures 37 and 38 present the results for 5N and Cu .
Within the scatter of the simulations they are equal over
the whole range of p. 'As in Chapter 1 this is a surprising
result, since inversion symmetry at each site is lost as bonds
are cut, so that the Cauchy relation equating. (w and Cu.
should become invalid. It is possible that Cauchy's thoerem
may be true under more general conditions than have been
previously assumed.

Figure 39 shows the simulation results for Cu and Cw
for each track and each value of ”Al. The lines shown are

from the effective medium theory (EMT) of Section 3.4.

Table 3.

Y/x

Track

113

ratios used for numerical simulations

7/4

0.25

1.0
2.0

0.1

0.25

2.5

:‘m‘i‘fl

CL/th

 

11.0
5.0
1.4

0 fl. 0.",

(h/QVV

114

 

2 TRACK#1 ' V 9
‘cu ‘3 C12 via-2.0 a
‘ 5)- °°44 0°12 yin-1.0 a
' °°44 V012 via-0.25 '
a

Elastlc Moduli
:‘

 

 

 

 

Figure 37. C." and C“ vs. p for track 1

115

 

 

 

 

2.5 . ° fl
TRACK #2 ' '
“344 '3 C12 7/0: - 2.5
2 - °°44 061': war-0.25 '
= °°44 V312 via-0.10
3 a
3 1.5 -
g
3 1 - a
. a
0.5 -
o 0
0 J O V V
0 0.25 0.5 0.75

Figure 38. Cu. and C... vs. p for track 2

 

116

Elastic Moduli

 

 

 

 

 

 

 

 

9° 0.». c 0 9w.
9. a . a

u . 3.8.. e» .
”.3 v \I I ”.0

a .M M: J
..n . 2 .

g . .
0.. u 1

O P p

O 0.». 0.. Pd.

wwosnm um.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

:50... 3 . J . ..u . 3.5.: 3 . .
., .
a 105.0 . zeroes
0 o: . 0.. .0 n:
. oArs Palace. 4
Pa . 03.
P b O n b
o 0.». ob 0.3 . o 0.3 o... 9.3 .
.b . .:SMme» . . .
.
isle;
oh ab 0.:. .
oh 9. oAWb .
¢.D 0.. L
oh oh .
o p .
o 0.3 o... 0.: .
MHmmneo somoww won «Hmorm H mam N

117

In all six graphs the agreement between EMT and simul-
ation is excellent, and perhaps even more impressively so
than in Chapter 1, as the shape of the curves is no longer
simply linear. As the rigidity threshold.pr= 0.5 is closely
approached, however, there are what appear to be very small
systematic deviations of the simulations from EMT, with the
experimental points being below the EMT curves. These
deviations mean that EMT is probably'ngt exact for this
model, but since they Egg very small, EMT describes this
sytem quite well indeed.

The EMT predicts that the ratio 6. law should go to a
fixed value at the critical point of any given track which
is independent of its initial value. The simulation results
for this ratio along tracks 1 and 2.are presented in Figures
40-41 and 42-43 respectively.

For track 1 the EMT predicts the critical value of (n/(Ly
to be 2.0 ( 'Aflt. - 1). Figure 40 shows the numerical
simulation result. Within the numerical scatter, the sim-
ulations appear to agree with the EMT result. Figure 41
shows the same result for the ratio but computed on three
configurations of a 20 x 20 square net. The finite size
effects encountered in using too small systems are clearly
seen. For track 1, the difference between Figures 40 and 41
are not that great but the larger system with more config-
urations clearly is an improvement. One should note that
even with the excellent agreement of the moduli with EMT,

small finite size errors in the moduli are amplified in the

118

 

' men? #1 ' ' l
o via-0.25

. war-1.0

4 .- o Via-2.0

 

0.5 0.0 0.7 0.8 0.9

Figure 40. 6.. la", flow diagram for track 1

 

119

 

 

 

 

 

 

5 ‘ FLOW TO ‘IUDPOIUT RR TRACK “’1
4.5 ‘ O V/d'OJur
n ”('10 .
fi [[4 ' 1. 0
4 a
.I
c 3.5 . '
1
1
I
C 3 '
4
4 I
2.5 . '
I
O
z I! a -P—-
EE 0 B'B“‘*’--.:._,,: - 1 ""
1.3 . ° 0 u-~u--.,.--r,v_.fi__°_
1 I l I I I J
0.5 0.5 0.7 9.9 3.5 1
P

Figure 41. C.. ICW flow diagram for traCk 1 (small system)

120

ratios because we are dividing two very small numbers.
Finite size noise has a much greater effect for the
.track 2 results. Figure 42 shows that within numerical
scatter, the ratio Ch/hm seems to become asymptotically
equal to 5.0, in agreement with the EMT result for this track.
Figure 43 shows the smae ratios computed using the smaller
systems described above. Based on Figure 43 aigng, we
would have said that there is large numerical disagreement
with the EMT result. Because we are computing critical
ratios, our results are quite sensitive to finite size
noise in the critical region, when the appropriate correl-

ation length becomes of the order of the size of the system.

121

 

12

' moi}: ' r '
aria-0.1
10 P aY/a-O.25
oyla-2.5
o
8 " 0 °
0: 0°° 0
\P o b O O
0' ‘0 .‘A. ‘
H a ..
2 ..
o l I | | I

 

 

 

0.5 0.0 0.7 0.8 0.9

Figure 42. cu /(m, flow diagram for track 2

122

 

 

 

 

 

 

 

 

 

12 ~Hou
TOFDCDPOINTFm'mAcx .3.
I 0 .
10 ~ /
O
6 ' a,“
a . . ‘0..'
c . o . o . . .'...'
1 ' 0 O . . ,0.O".' 0 VII . 0'1
1 ,»"° 770304:
/ A ,. "
c s . ' v 4 n. . 2.5
2 -v‘""."6 '7, v
{if ’ J.— I.
4 _' \s‘ . ..... L : t
‘\\‘A ‘
\
s“‘
‘m
2 .. «‘““
\fl—fi—n
a . I I l
i
3.5 3.5 9.7 3.8 a 9
' 1
p

Figure 43.

("/00 flow diagram for track 2 (small

system)

123

Section 3.6. Conclusions

We can conclude that EMT does as well or better for
first 32g second neighbor central forces as it did for
nearest neighbor bonds only. This agreement between
EMT and simulation is far better than in most other problems
where EMT has been applied. The reason for this may be
that since the EMT always gets the initial slope right for
gny_problem. and in this case EMT also predicts the
critical threshold quite accurately, that with both ends
of the curve fixed there is little room for error in
between. This of course does not address the question of
why the EMT or constraint counting theory prediction for
p'is so accurate. 4

To properly understand the implications of our result for
the critical behavior of the “ICU-I ratio, one should recall
that there are two different classes of problem whose
elastic properties are being studied currently. Class 1
includes those problems where the elastic moduli go to zero
only when the system becomes physically disconnected. This
includes all percolation problems where enough microscopic
forces have been specified so that geometrical connection
implies elastic connection. Class 2 is the class of all
rigidity percolation problems: that is, all problems in
which geometrical connection does not necessarily imply
elastic connection. As was stated in Chapter 1, all problems

studied in this thesis are Class 2.

124

Our EMT predicts that CS/QmIgoes to a fixed point value
at the rigidity threshold which does not depend on the initial
value of the ratio but does depend on geometry, i.e., the
track followed in phase space and thus the critical value of
p1 (see Figure 36). The numerical simulations confirmed
this prediction quite well. This is the first and only such
finding so far for a class 2 problem, analytical or numerical.
Bergman and Kantor (1984) predicted, based on mean field
theory and an exactly solvable fractal model for percolation
backbones, that the critical value of 6" kw would go
to a universal value dependent only on dimension. Our result
clearly contradicts this; however, their prediction was
really only for class 1 problems. Bergman (1984) showed
for a honeycomb lattice class 1 problem that Ct/a~ did go
to a critical value independent of the starting value.
However he had no way of varying the geometry of his problem.
Schwartz (1984) found the same type of behavior for the
Born model, which gives a class 1 problem. Thorpe (1984)
showed, using an EMT, for the class 1 problem of randomly
punching elliptical holes in an elastic sheet, that the
critical value of (flfiw depended on the ratio of the
lengths of the semimajor and semiminor axes of the hole.

This is the only indication thus far that Bergman and
Kantor's prediction might be incorrect for a class 1
problem. There are no numerical results as of yet for
different geometries of class 1 problems that would support

or contradict their prediction. Feng et a1. (1984) confirmed

125

Bergman's (1984) result for the critical exponent of Cu by
doing simulations on a square net, but gave no results for
“/4... since they apparently did not compute (w . Finally,
the node-link picture developed by Kantor (1984) for class 1
problems assumes that the elastic properties of a lattice
will become isotropic (i.e.,C'Cu’ 2"“! ) as the critical
threshold is approached. There is no good numerical

or analytical work for class 1 problems to support or deny
this assertion thus far. For the class 2 problem studied in
this chapter, the square net does 323 become isotropic at

at §‘- 0.5 except perhaps at one isolated point along the
critical line. The EMT would predict this point to be

p1 - 0.6, p2 = 0.4, an isolated point of no apparent signif-
icance.

These questions of the universality of the ratios of
critical amplitudes and of isotropy may turn out to be
another area in which clear distinctions between ordinary
connectivity percolation (class 1) and rigidity percolation

(class 2) can be made.

APPENDICES

APPENDIX A

APPENDIX A

Trimming and Supertrimming_Rules for Central Force Networks

The distinction between trimming and supertrimming is
made clear by noting that trimming involves single sites
only. Every other kind of removable unit (denoted a floppy
unit) is classified under supertrimming.

The complete set of rules for determining whether a site
is trimmable or not are summarized in the following state-
ment :

A single site in a central force
network is trimmable if, in d dim-

ensions, it is connected to the rest
of the network by d or fewer bonds.

The justification is clear: any mass point has d degrees
of freedom, so it can always satisfy d or fewer constraints,
a bond being a single constraint in this case.
Units which are supertrimmable can be classified accord-
ing to the type and number of their connections to the rest
of the network. A connection by a single bond we call type
1, with the number of these given by r1. A connection by
a pin joint we call type 2, with the number of these given
by r2. Figure 44 shows an example of both types of connection

for the triangular net.

126

127

 

Figure 44. Two types of connections between floppy
units and the rest of the network

128

It is important to realize that type 1 connections
quench one degree of freedom of the attached floppy unit
while type 2 connections quench two floppy degrees of free-
dom. Therefore in two dimensions, which is really the only
case of interest as far as supertrimming goes, the rule to

determine whether a given unit can be supertrimmed or not is:

if r1 + 2r2 5 3 supertrim
if rl + 2r2 > 3 cannot
supertrim

Three is the number of degrees of freedom of any extended body
in two dimensions--two translational and one rotational.

Note that rl + 2r2 >23 for a given unit'gggg_ngt necessarily
mean that the unit is rigid. There can be internal floppy
modes which still allow the unit to be removed. That is

why trimming and supertrimming as formulated are incomplete.
As these internal floppy modes can be extremely complicated,

we have been unable to rigorously identify higher order rules.

APPENDIX B

APPENDIX B

Extrapolation Technique for Elastic Energy Relaxation
Empirically we found that as a network relaxed the elastic

(a)

modulus C near the final stages n of relaxation obeyed the

asymptotic form

C(n): C(‘fl't a b’h

where a > 0, b > 1, and b itself was slowly decreasing with n
but at‘a rate mggh slower than that of C“: Each step in n
really represents 50-100 actual relaxations per site. Call
this a iéggg relaxation step. If N of these large steps
were used at a given value of p, then by fitting the

'I’ -
exponential form above to C” , C“ u, and CM) , we can solve

for C”) to get

u) ((014) - (Cu-u) L

(

(3°) C:
(1‘ -¢ ("
CH‘2C(~)+C”U

C :-

 

Ca') will only be an upper bound, however, since the

extrapolation technique assumes b constant while in fact the

slow decrease in the value of b will result in a final

129

130

correctly relaxed value of the elastic modulus that is
somewhat smaller than C”’. In effect the extrapolation
technique is just a cheap way of getting several thousand more
small relaxation steps. Table 4 gives an example where the
extrapolated modulus can be compared with actual further
relaxation steps. This data is taken from a BCC computation

for the shear modulus, p=0.77.

Table 4. Extrapolation vs. relaxation

 

 

 

No. of small Value of modulus
relaxation steps before extrapolation Value after
800 0.0745 ' 0.0732

4800 0.0736 0.0702

APPENDIX C

APPENDIX C

EFFECTIVE MEDIUM THEORY OF PERCOLATION
0N CENTRAL FORCE ELASTIC NETWORKS

SHECHAO FENG' and M. F. THOR? “
Schlumherger-Doll Research
PO. Box 307
Ridgefleld, CT 06877
E. GARBOCZI
Department of Physics and Astronomy

Michigan State University
East Lansing, MI 48824

Abstract

We show that effective medium theory gives an excellent description of regular
lattices when nearest neighbor central force springs are present with probability p.
Efl’ective medium theory shows that all the elastic constants go to zero at p“; - 0 pc
where d is the dimension and 1),, is the effective medium estimate or the ordinary
percolation threshold.

Present address: Department 0! Physics. Harvard University, Cambridge. MA 02138

" Permanent addren: Department 0! Physics and Astronomy. Michigan State University. East
Lansing. MI 48824

131

132

1. INTRODUCTION

There has been considerable interest in the elastic properties of random systems
recently‘“. While much of this attention has focused on the critical properties of the
random elastic network as it breaks up, it is of considerable interest to describe the
overall behavior. We focus our attension in this paper on the central force elastic
percolation rnodel.I We show that the simplest efl'ective medium theories give
excellent agreement with numerical simulations for most values of p, the bond
occupation probability, except for a very small range of values of p near pm, on both
the triangular network and on the f. c. c. lattice.

Although it is unlikely that the present model will apply directly to any real physical
system, it is an important one as it is the simMest member of a general class of
models. These models are made up of units that can be connected but still transmit no
. elastic restoring forces. In the present case the units are bonds and a pair of two and
just two connected bonds form a 'free hinge“ as there are no angular forces. The
whole network has finite elasticity for high values of p because the high connectivity
'locks" all the free hinges. As bonds are removed from the pure system, the local
My regions are crested which eventually prevent the rigid regions (i.e. the parts of
the network that have finite elasticity) from percolating and the system loses its
elasticity. Thus the breakup of the system is determined by rigidity flotation and not

ordinary anectivity ngcolatign.

The most important physical manifestation of this phenomenon is in the c_o__valent
rand__9_m networks’ where its local movement is not a free hinge but the dihedral
angle _—associated with three. connected bonds. The energy required to change the
dihedral angle is small compared to that needed to change the bond lengths and bond
angles, so it is reasonable to neglect it as a first approximation. This leads to a division
of covalent random networks into low co-ordinated Elvmeric g_la_s_s_es_ and high co.-

ordinated amgmhous glids.S

The other class of models for the elastic percolation phenomenonm'4 involve
specifying a suflcient number of microscopic forces so that all connected parts of the
lattice are rigid by themselves and the elastic percolation transition occurs at the
ordinary percolation threshold pc, albeit with difierent exponents.2~3"

The system under consideration is made up of Hooke springs connecting nearest
neighbor sites i,j to give a potential .

V-a — 22 [(E- uj) rulzpi, (1)

2<li>

where the angular brackets denote a sum over nearest neighbor pairs which are
connected by springs with spring constant a and pi, is a random variable that is
associated with each bond and is 0. l with probability l-p, p. The Ti. are the
displacements from equilibrium and ii,- is a unit vector connecting nearest neighbor

133

pairs in equilibrium.

In the next section we show how a constraints counting argument can be used to give
an estimate of pa. at which the elastic constants of the syStem vanish. In section 3
and 4 we develop effective medium theories for the. elastic constants for p > pm
from two difierent viewpoints, both of which lead to the same result. The first one is
a direct generalization of the method developed in the study of electrical condution
near percolation threshold, and the second one applies the coherent potential
approximation to the present problem. ,

2. CONSTRAINTS METHOD

The simplest way to estimate where the transition takes place is to use a constraints
argument’. When p is small, the system consists of disconnected pieces and hence
has many zerg frgguencv modes whose number is given by the number of degrees of

freedom (dN) minus the number of constraints (é-zbip). Thus, the efi'ective medium
estimate of the fraction f of zero frequency modes is given by

r- (dbl - é-szl/dN

_ -22
_1_2d (2)

so that f goes to zero at p,- - 2d/z.

Next consider comparison of this result in the numerical simulation. We have
computed f numerically for a‘ 168 atom triangular network (see inset in figure 1) and a
108 atom f.c.c lattice (see inset in figure 2). Both lattices had periodic boundary
conditions and f was obtained by directly diagonalizing the dynamical matrix formed
from (1) and counting the number of modes with eigenvalues of zero. It can be seen
that equation (2) describes the results well except near pm where the very small
deviations from (2) are due to a combination of finite size and critical efi‘ects.

A similar constraint counting argument for one, rather than d, degrees of freedom per
site, which is the case for the electrical conduction problem, leads to the effective
medium estimate for the ordinary percolation (i.e. the connectivity problem)‘ of
pc - 2/2 and hence

pa. - d pg (3)

The transition described here takes place well before ordinary percolation occurs
because many connections in the network produce no elastic restoring force’. A
simple example is shown in figure 3. Configurations like this correspond to one

example of the Mg: rggigns.’

1'34

3. STATIC METHOD

For p >pm, it is of interest to develop a mean field theory for the elastic constants.
We have done this in two ways; both of which lead to the gang result. The one that
is simpler conceptually is to adapt an argument used by Kirkpatrickz for the electrical
resistor network. Suppose a uniform stress (uniaxial, hydrostatic, etc.) is applied to a
lattice where all the springs are a, and that the atoms labeled 1, 2 in figure 4 have a
relative displacement, along in of Sun. Now substitute a single 'wrong bond' on as
shown in figure 4 and imagine an e_x_tgg external [4355 f applied to 1 and 2 as shown to
restore l and 2 to their positions before a was substituted for am. It is then easily
seen that f should be:

f- Bunk, - a) (4)

From the superposition principle, the relative displacement au between 1 and 2
induced by f when the system is unstressed is the same as the e_x_tg displacement
between 1 and 2 when there is an applied uniform strain an, but no f. Then an is the
'lluctuation' in relative displacement of 1,2 'due to the introduction of the wrong bond
0:. The relation between f and 80 can be obtained in the following way. If the force f is
applied to l, 2 when gil_ the springs are am, there will be an efiective spring mutant
0“ - afila‘ where 0 < a' < 1 takes account of all the connections between 1', 2
including the direct one in this uniform system. We will calculate a‘ later, but for now
treat it as a known constant. If the a, spring between 1 and 2 is removed, then the
new efi'ective spring constant a" between 1, 2 is

a", - anla‘ - an ' (5)

\

as shown in figure 4. If a is added in parallel to a’. and the force f applied, then the
relative displacement Bu of l and 2 is given by

f- 50(0',‘ + a) (6)
or
Eu - 3%(0, - e)/(au/a' - an 4' a) (7)

The effective medium result is obtained by choosing a, so that the average value
- <50) of an is zero to give '

< 0“-" >-o (3)

 

135

For a distribution N0) of spring constants, this leads to

P(a)da _ ' .
I l-ae(1_a/an) 1 (9)

 

For the present case of interest P(a') - p5(a - 0') + (1 - 0) 8(0') gives

a 1-3‘

which goes to zero when p“. - a'.

The quantity a‘ is obtained for the perfect lattice as follows. The force P, on the atom
i is given by

- ‘ 8V - ..
F3 " ‘ '7. " ’ D--'ll (11)
an, 12 u l

where

l . . .
.5 Tamriirij jail
" “a 2818: 5"
{ jii

(12)

These equations can be inverted by Fourier transformation to give
u, - - '5“(TE)~'I-‘, (13)

where D(k) is the dxd dynamical matrix for a Bravais lattice, and 15k and Ti,‘ are the
Fourier transforms of the F, and Ti, respectively.

566) - 25,.xpn‘io', -r,>1
ii

- c.2[l -- exp(iaE-3)l§§ (14)
a

where 5 is a unit vector in the direction o_f_ one of the z nearest neighbors and a is the
nearest neighbor separation. Putting Fj - f ium, - 532), the negative of the

136

externally applied force as shown in figure 4, we find that

a, -r, - if; ;[2 - exp(iaI-iu) - .xp(-iaT-r,,)1'n“(T)-r,2 - (15)
k .
which defines
a‘ - 55 212 - exp(iak-fu) - exp(-ia-k'~i'u)liu'D"(kl-in (16)
E

As all bonds are equivalent, we can replace in by any of the nearest neighbor bonds
8, sum over all nearest neighbor bonds and divide by z to give

a' - 3‘1?- z Tr[(l - .xptiaT-ii)(iS)-B“(Tn
N2 u

2- 2d
Nz g Tr(1) z p... (17)

. Thus, the result (11) can be written’*8 as

32-222.— (13)
a l-pm

Since all the elastic constants C,,- of the network are linear functions of a“ only, we
conclude that they all go to zero linearly at pm, and any ratio of two elastic constants
is independent of p. In particular the ratio of the bulk modulus to the shear modulus
of the triangular network is alwag 2.

Next consider comparison of this result with numerical simulation. The straight lines
from eqation (18) are plotted in figures 1 and 2 for the triangular and the f.c.c. lattices
and show excellent agreement with simulations except for a very small critical region
near pg... Such close agreement between simulations and efi'ective medium theories
is rare although it does occur for the conductivity problem‘ in 2d. The deviations are
larger in 3d for the conductivity problem‘ than in the results for springs in figures 1
and 2. Indeed 2d/z appears to be a superior estimate for pm than 2/2 is for pc. This is
particularly apparent in comparing the results for the f.c.c. lattice in this paper and
those in ref. 6. In general we would expect efi'ective medium theories to do better in
higher dimensons.

The numerical simulations were done by removing bonds at random in a lattice with

137

periodic boundary conditions, imwsing an appropriate external small strain 5 (which
redefines the vectors that define the large periodic cell) and relaxing as fully as
possible by moving atoms towards positions where there is no force on them. The
elastic constants were obtained by computing the potential energy of the system via

equation (1) and using V - i—Cez. This procedure has some advantages over that

used previously’ in that all atoms are treated equivalently and there are no “surface
atoms“. Typical strains used were s - 10".

There are some subtle and important difi'erences in the numerical simulations done in
this work and those done previously‘. In figure 3, the two bonds can be removed as
they do not produce a restoring force. (Dangling bonds are also removed as in the
conductivity problem‘) However, when 0 - 180° in figure 3, there are two ways to
proceed. In ref. 1, it was assumed that there was a restoring force as given by the
potential (1), whereas in the present work we assume there is none and remove such
configurations. This has the effect of increasing P“. from 20.58 nearer to 2/3. We
remove such configurations because 'buckling' can occur in compression rendering
the connection in figure 3 inefi'ective even if 9 - 180°. This model can be visualized
as replacing all 180‘ bonds by 180‘ - A, evaluating the elastic properties for a small
strain ¢< <A, and letting A-O. Both models (i.e., ref. 1 and the present one) are
equally valid and of course the efi'ective medium theories cannot distinguish between
them. In the numerical simulations to calculate f, as described previously, it was found
that f was slightly larger for the A -' .0 case than for the small A One. However this
was only appreciable around pm and too small an efi‘ect to show on the scale of the
inserts in figures 1 and 2. ‘ ‘

4. COHERENT POTENTIAL APPROXIMATION

We now consider a multiple scattering approach in which all repeated scatterings from
the same bond are summed to give a coherent potential approximation (C.P.A.)°. We
find that same result (18) is obtained from this method. Writing the Hamiltonian for
the system as '

H - Ho + V (19)
where
3i: “m
H -‘ -— + —- Ira, -n'.)-r--12 (20)
° ,2 2M 2 (1:5, I '1

and

in,

v - 3:75am, - spent! (21)

represents a single “defect“ bond in an efi‘ective medium a, as shown in the left hand
part of figure 4. It is convenient to define Green‘s functions for the system described
by Ho as

<o|n,|n><nra,uo> <oln,|n><n|ri,|o>l

o-u.+ao o+u,-oo

 

 

'1'5ij - 2 (22)
I

where T5,, is a d x d matrix. A similar quantity 23,,- is defined for the system described
by H. We rewrite (21) in matrix form with

vii - (a - C‘)f‘2f1zmfi (23)

where mi,- - (8‘58" + 82,823 - 8,,82; - 82,8") arises from the translational invariance of
each bond. It is easy to show that the Dyson equation is satisfied, i.e.,

§-?+$VE ' (3)

where ‘17, ‘6 and V are considered as matrices in'both site indices OJ) and in a d
dimensional ‘vector space. Remembering that we are considering a single defect bond,
Eq.(25) can be rewritten as

G-3+53$ am
where

Ti] , a .- a,ll

 

l " 2(0 "’ Gn)f12(Pu ‘ Puff” U 12 J

has the same form as Vi,- but with a renormalized coemcient. Setting <T> - 0, we
find

a " “m
1 " 2(a " amfiu’fin "' $12)?”

 

> - 0 (28)

which is analogous to (8) except that it applies at all frequencies a and cum). Note
that all elements of the T matrix are set equal to zero by (28) as all the T matrices

139

have the form of a scalar multiplying m,,-. The method of section 3 can also be easily
generalized to finite frequencies.

The pure syst_em Green functions Pu obey the equation of motion (noting that
Pu ' P22 and P12 "' P21) '

- 2 1““ . - - ‘
M0 P" - 1 4’ Tfu'(Pn "’ P12)?” (29)

where P" is the magnitude of the isotropic site diagonal Green function. As (92"0,
equation (28) reduces to the previous result (8) with again a“ - 2d/z.

5. CONCLUSIONS

There is a slight difi'erence between these two versions of effective medium theory for
p> pm. In the first, the elastic constants Cit are calculated, and the various masses are
irrelevant. In the second, the sound velocities are calculated in the long wavelength
limit. This leads to the elastic constants if it is assumed that the efi‘ective mass is
independent of p. In reality, the effective mass will depend on p, but this dependence
is expected to be much weaker than the dependence on p of the elastic constants near
pm, analogous to the conductivity case. An account for such efi’ects is obviously
beyond the capabilities of efiective medium theories.

.To summarize, we have shown that simple efi'ective medium theories give a
remarkably good overall description of the dilute elastic systems with Hooke springs.
The detailed behavior around pa. is a subject still under study.

We should like to thank P. Sen for many useful conversations on this problem, and
also R. 1. Elliott, B. Halperin, B. Nickel and L. Schwartz for clarification on particular
points. Two of us (S. F. and M51.) would like to thank Schlumberger-Doll
Research Center for their hospitality in the summer of 1984. The support of the
O.N.R. and NSF. by MIT. and 5.6. is gratefully acknowledged.

f"

140

REFERENCES

. S. Feng and P. N. Sen, Phys. Rev. Lett. 52, 216 (1984)
D. J. Bergman and Y. Kantor, Phys. Rev. Lett. 53, $11 (1984)
. Y. Kantor and I. Webman, Phys. Rev. Lett. 52, 1891 (1984)

. 8. Alexander, preprint; S. Feng, P. N. Sen, B. I. Halperin, and C. J. Lobb,
preprint; 1... Benguigui, preprint;

M. F. Thorpe, J. Non Cryst. Solids 57, 355 (1983)

6. S. Kirkpatrick, Rev. Mod. Phys. 45, 574 (1973)

This result was first given, without proof, in G. R. lerauld, L. E. Scriven and H.
T. Davis, preprint and also by M. Sahimi, private communication.

. M. Sahimi, B. D. Hughes, 1.. E. Scriven, and H. T. Davis, J. Chem. Phys. 78,
6849 (1983)

For a review, see R. 1. Elliott, J. A. Krumhansl and P. 1... Leath, Rev. Mod.
Phys. 46, 465 (1974) '

141

Figure Captions

Figure 1

The elastic constants C" and C“ averaged over three configurations for a 440
atom triangular network. For the pure system (p-l), C“ - Cu/3 - 0‘5/4.
The inset shows the fraction of zero frequency modes f for a 168 atom triangular
network averaged over three configurations. ‘The straight lines are from the
effective medium theories described in the text.

figure 2

The elastic constants C", C“ and B- (Cu + 2Cu)/3 averaged over three
configurations for a $00 atom f.c.c. lattice. For the pure system (p-l).
C" - 2C1; - 2C“ - cuff/b where b is the nearest neighbor separation. The
inset shows the fraction of zero frequency modes f for a 108 atom f.c.c. lattice
averaged over three configurations. The straight lines are from the effective
medium theories described in the text.

Figure 3

Showing a two coordinated bridge connecting two regions. This bridge is
inefi'ective in transmitting any elastic restoring force and can be trimmed. . '

Figure 4

Showing the notation for constructing the efi’ective medium theory. On the left a
single bond between sites 1 and 2 is modified by having a spring constant a rather
than a“. the strain existing before the modification can be reestablished by
applying a force f across the bond as shown. On the right, we show an equivalent
circuit for the bond joining 1 and 2 as described in the text. The springs a" and
a are connected in parallel.

Elastic Moduli

142

 

 

0
0

 

0.1

 

d'fl'

" Triangular

 

 

 

 

     

. l

l l l i l
0.2 0.3 0.4. 0.5 0.6 0.7 0.8 0.9
D

Figure I

-. Elastic Moduli

1.6

143

 

 

1.4—
1.2}-
1.0..

0.8 -

 

0.6 -
0.4 -

0.2 '-
0

 

1

p .

1 L 1
0 0.1 0.2 0.3 0.4 0.5 0.6 0

D

Figure 2

 
  

.L
.7 0.8 0

      

.9

 

 

1.0

144

145

 

 

 

 

 

 

0W. x ./
-1. bna
: x I, WV i
3
N

-n

Figure A

LIST 01“ REFERENCES

LIST OF REFERENCES

R. W. Ashcroft, N. D. Hermin, S lid State Ph aics, Bolt Rinehart
and Winston, Philadelphia (1976).

R. J. Bell, Reports on Progr. in Phys. ;§, 1315 (1972).

R. J. Bell, Methods in Camp. Phys. 12, 215 (1976).

R. J. Bell, in Excitatiens in Diserdered steegs, ed. by M. F. Thorpe,
Plenum, New York (1982).

D. J. Bergman, to be published in January 1985 issue of Phys. Rev. B.
D. J. Bergman, Y. Rantor, Phys. Rev. Lett. g;, 511 (1984).
P. Dean, Rev. Mod. Phys. 44 127 (1972).

R. J. Elliott, J. A. Krumhansl, P. L. Leath, Rev. Mod. Phys. fig, 465
(1974).

S. Feng, P. R. Sen, Phys. Rev. Lett. 2;, 216 (1984).

S. Feng, P. N. Sen, B. I. Halperin, C. J. Lobb, Phys. Rev. B ;Q, 5386
(1984).

S. Feng, H. P. Thorpe, B. Garboczi, Phys. Rev. B e;, 276 (1985).
B. J. Garboczi, H. P. Thorpe, submitted to Phys. Rev. B (1985).
B. Be, H. P. Thorpe, submitted to Phys. Rev. Lett. (1984).

Y. Kantor, I. Webman, Phys. Rev. Lett. 2;, 1891 (1984).

S. Kirkpatrick, Rev. Hod. Phys. 42, S74 (1973).

C. Kittel, Ingeedeceion te Selid Stete Phyeics, 3rd ed., J. Wiley and
Sons, New York (1967).

I. H. Lifshitz, Adv. Phys. 1;, 483 (1964).

A. B. H. Love, A Treetise en the Maehemeticel Theory 0: EleseieiEI,
Dover, New York (1944).

146

147
J.F. Nye, Ph sical Properties of Crystals, Clarendon Press,
Oxford (19 T.

L. Schwartz, P. Sen, 8. Feng, M.F. Thorpe, to be published
(1984).

M.F. Thorpe, J. Non-Crys. Sol. El, 355 (1983).

R. Zallen, The Ph sics of Amorphous Solids, J. Wiley and
Sons, New YorE (1983).

 

ATE

ll (HillIiiiililiiilililliliiliiliIES
0 6 1 1 5 0 7

93 03

 

mililiilillilli
3 12