39"?" .
' “a,

u
x

to -'$u
.. an.

m “3::

’1 3%?

mm“,
439.?

3

m

 

”'52::

' “‘47

7

~

.v-fi.

.’

a i
4‘ 12%.]
.‘ .flfig ,

L. ’ ‘ ‘

A

.4

w...
’0'.“
.- w m.
77......
MM... mum;
33:15:”! 2:53,.
team
.. ‘

5

 

......
‘ 4::

-.

TH ESlS

                                                                                                                             

3 1293 01555 025

This is to certify that the

dissertation entitled

"Propert ies of Disordered

Continuous and Discrete Soli

presented by

Bran islav Djordj evic

has been accepted towards fulfillment
of the requirements for

Ph . D . degree in PhESiCS

him

Major professor

 

 

Date _May_l._129_6__

MS U i: an Affirmative Action/£4 ual Opportunity Institution 0- 12771

 

 

 

LIIRARY
Michigan State
University

 

 

 

PLACE IN RETURN BOX to romovo this checkout from your rocord.
TO AVOID FINES roturn on or before date duo.

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

l

 

 

 

 

 

 

 

 

MSU It An Affirmative AdloNEquol Opportunity Initiation
Want

 

 

PROPERTIES OF DISORDERED
CONTINUOUS AND DISCRETE SOLIDS

By

Branislav R. Djordjevié

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
1996

ABSTRACT

PROPERTIES OF DISORDERED
CONTINUOUS AND DISCRETE SOLIDS

By

Branislav R. Djordjevié

In this thesis I have studied both continuous and discrete disordered solids.
The focus is on the calculations of the effective conductivity and the corresponding
spectral function in random binary composites, and on the computer modeling of
amorphous materials. The spectral function is fully determined by the geometry
of the two—phase medium. I examine the case of circular inclusions randomly
distributed in a conducting sheet. I calculate the effective conductivity to second
order in the concentration of inclusions f. The spectral function for this problem
is found to be a truncated Lorentzian with weight f and with the width that is
proportional to f, and the final result is written in a very simple form. Vari-
ous two dimensional problems are discussed and it is shown how one can use the
conformal mapping technique to calculate the effective conductivity, when the in-
clusions are polygonal holes. I have also studied the properties of a topologically
disordered networks — amorphous solids. I computer generated a very large (4096
atom) model of amorphous diamond using the Wooten—Weaire method. Starting
from the perfect diamond crystalline super—cell, I first randomize the structure by

introducing topological disorder through random bond switches at a sufficiently

high simulation temperature. Subsequently, the simulation temperature is re-
duced in stages, and the topological and geometrical relaxation takes place using
the Keating potential. After a long annealing process, a random network of com-
paratively low energy is obtained. I calculate the pair distribution function, and
all other relevant quantities for the final relaxed structures. I minimize the total
strain energy by adjusting the density of the sample. I compare the results for
amorphous diamond with the results I obtained for amorphous silicon, and also
with the recent experimental measurement of the structure factor for (predomi-
nantly tetrahedral) amorphous carbon. Finally, I examine the elastic properties
of amorphous networks. The amorphous diamond model is bond—diluted, and a
sequence of amorphous structure with different mean coordination (r) is obtained.
I calculate the bulk modulus as a function of the mean coordination. I compare
my results with those that were obtained using the crystalline super—cell as an
underlying structure on which the bond dilution is performed. My simulations
show that the bulk modulus as a function of the mean coordination obeys the
same power low, no matter what the underlying lattice is - amorphous or a crys-
talline. All results of my research are presented on the comprehensive historical

background of the modeling of amorphous networks.

I dedicate this Thesis to my late father Rista Djordjevié

who suddenly died in July 1995.

I equally dedicate this Thesis to my wife Dragana Vasiljevié.

ACKNOWLEDGEMENTS

I should like to thank my thesis advisor Professor Michael F. Thorpe for
his guidance.

I am grateful to Professors M. Berz, S.D. Mahanti, J. Hetherington and J.
Cowen for serving on my guidance committee. I am indebted to Mrs. Janet King
for her readiness to help me whenever it was necessary. Her knowledge, experience
and kindness, made my life in this department much simpler. My gratitude goes
also to Dr. George Perkins who was always available and willing to help me with
computers.

I wish to thank Dr. Normand Mousseau for numerous discussions we had
at all stages of this work, for his always valuable comments and suggestions, espe-
cially regarding Chapter 8. I am grateful to Professor F. Wooten for introducing
me in the exciting field of computer modeling of amorphous networks, to Dr. A.R.
Day for his helpful discussions and comments, especially regarding Chapter 3, and
to Dr. D. Jacobs for the fruitful discussions we had and for his great and infectious
enthusiasm for hard work.

I had a very good time working in my office thanks to my wonderful
colleagues and office mates: H. Xiao Yuqing, S. Hyun, B. Chen and T. Nagy.

I want to acknowledge the moral support given to me by my friends N ;
Mousseau, G. Overney, F. Liu, J .-H. Pinto, S. Bezborodnikov, J. F algas, M. Roura,
E. Pekker, F. Leganes Nieto, F. Fernandez Pirias, and J. Bredeck. Their friendship
marked my stay in East Lansing, and made my life here much more enjoyable and
interesting.

I feel happy to acknowledge that I was able to complete my research and

iv

this thesis thanks to the generous financial and moral support of the following
friends of mine, who lent me money during the last difficult year of my studies,
and thus made this thesis possible. Their names in alphabetical order are: G.
Adam, A. and V. Bakié, S. Bezborodnikov, D. Bogdanovié, E. BoZin, Z. Bucalo,
V. Bukovinski, D. Cvetkovié, G. Debelnogié, L. and E. Gawarecki, V. Grigorescu,
B. Hare, I. Jasiuk, Z. and K. Jovanovic’, M. Ledvij, M. Markovié, T. Nagy, B.
Nakanishi, M. Ostoja-Starzewski, M. Palaékovic', J.-H. Pinto, T. Pokoudina, D.
Pokudin, Lj. Popovié, M. Popovié, A. Radichevich, B. and M. Schall, V.B. Vasil-
jevié, O. Vulié and J. Welch.

Finally, I want to express my deepest gratitude and love to my wife Dra-
gana Vasiljevié, who stoically endured not only my continuous personal absence,
but also the absence of my income during the last twelve months of my stud-
ies. This thesis has largely been done thanks to her unselfish and unconditional
saCrifice, support and love.

Financial support from Michigan State University, the National Science
Foundation and the Center for Fundamental Material Research are gratefully ac-

knowledged.

Contents

Abstract ..................................

Acknowledgments ............................
1 Introduction

2 Continuous Disordered Media

2.1 Introduction ..............................
2.2 Effective Conductivity ........................
2.3 Spectral Functions ..........................
2.4 Two-dimensional Case ........................

2.4.1 Reciprocity Theorem .....................

2.4.2 Effective Conductivity ....................
2.5 Conclusion ...............................

3 Review of the Results for various two-dimensional systems

3.1 Introduction ..............................
3.2 Elliptical Inclusions ..........................
3.3 Conformal Mapping Technique ....................
3.3.1 Example ............................
3.3.2 Schwarz—Christoffel Transformation .............

vi

iv

10
11
13
14

15
15
16
21
22
25

3.4 Circular Hole ............................. 29

3.5 Slit ................................... 34
3.6 Effective Conductivity for Polygonal Holes ............. 36
3.6.1 Equilateral Triangle ...................... 37
3.6.2 Square ............................. 38
3.6.3 Regular Polygon ....................... 40
3.6.4 Diamond Shape ........................ 42
3.6.5 Star—shaped Holes ...................... 45

3.7 General Case of Inclusions with Sharp Corners ........... 47
3.8 Random Resistor Network ...................... 50
3.8.1 Dilute Limit .......................... 53
3.8.2 Effective Medium Theory .................. 55

3.9 Conclusion ............................... 58
4 Pair Terms for Circular Inclusions 61
4.1 Introduction .............................. 61

4.2 General Form for the Two—dimensional Conductivity to the Second

Order in f. .............................. 63

4.3 Dipole Moments ............................ 66
4.4 Perfectly Conducting Inclusions ................... 71
4.5 Effective Conductivity ........................ 75
4.6 Analytic Continuation of the Pair Term ............... 80
4.7 Spectral Representation ....................... 83
4.8 Conclusions ............................... 90

5 Discrete Disordered Media 93

vii

5.1 Introduction .............................. 93

5.2 Topological Disorder ......................... 93
5.3 Distribution Functions ........................ 97
5 4 Structure Factor vs. Pair Distribution Function .......... 106
Historical Background of Modeling Random Networks 111
6.0.1 Hand—built CRN Models ................... 112
6.0.2 Computer—built CRN Models ................ 122
6.1 Wooten—Weaire Method ....................... 125
6.1.1 Keating Potential ....................... 128
6.1.2 Wooten—Weaire Bond Switch ................ 130
6.1.3 Initial Randomization of the Crystalline Lattice ...... 132

6.1.4 Annealing: Geomet" ‘al and Topological Relaxation . . . . 135

6.2 Further Developments ........................ 140
6.3 Conclusions .............................. 146
Modeling Amorphous Diamond 147
7.1 Introduction .............................. 147
7.2 Wooten-Weaire Method ....................... 148
7.3 Discussion of Amorphous Carbon .................. 150
7.4 Results and Comparisons ....................... 152
7.5 Discussion ............................... 157
7.6 Conclusion ............................... 160
Elastic Properties of Random Networks 162
8.1 Introduction .............................. 162
8.2 Rigid and Floppy Networks ..................... 163

8.3 Elastic Constants of Random Networks ............... 166

8.4 Bond—dilution of Amorphous Diamond Network .......... 170
8.5 Strain Energy and Bulk Modulus .................. 172
8.6 Conclusion ............................... 186
Conclusion 188
Bibliography ............................... 193

ix

List of Figures

3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
3.10
3.11
3.12

4.1
4.2
4.3
4.4
4.5

2d example of a random composite. .................

Example of conformal mapping. ...................
Schwarz—Christoffel transformation of polygonal boundaries.

Mapping of the half—strip onto the half—plane. ...........
A sketch of the flow around the circular hole .............
Slit perpendicular to the current flow .................
Equilateral triangle ...........................
Square ..................................

Diamond. ...............................

Spectral function for triangular inclusion of general 7 ........
Spectral function from effective medium theory. ..........

Spectral function for random resistor network. ...........

Dipole images for two circles in parallel filed .............
Dipole images for two circles in perpendicular electric field .....
Total dipole moments. ........................
Voltage between two circles. .....................

Interpolation formula for F(7). ...................

25
26
30
35
37
39
43

49
57
59 l

67
70
73
76
81

4.6
4.7
4.8
4.9

5.1
5.2

5.3
5.4
5.5
5.6
5.7
5.8
5.9
5.10
5.11
5.12

6.1
6.2
6.3
6.4
6.5
6.6
6.7

Imaginary and real part of F(;r). ..................
Spectral function for several values of f ...............
Loretzian vs. exact spectral function. ................

The function ImiM(.r) in the dilute limit .............. i.

A piece of a crystalline diamond structure. .............
Schematic sketches of the atomic arrangements in three states of
matter ..................................
A piece of an amorphous diamond structure .............
Crystalline and amoprhous form of a binary compound .......
Diffraction patterns of crystalline and amorphous structure.

Typical Pair Distribution Function for a crystal and gas. .....
Typical Pair Distribution Function for amorphous materials. . . .
Typical shapes of Radial Distribution Functions ...........
The origin of structural features in RDF . .............
Typical Pair Correlation Function for an amorphous solid .....
Typical shape of Reduced Radial Distribution Function . .....

Structure Factor for a crystal and amorphous structures. .....

Zachariasen’s CRN model for amorphous silica. ..........
The first model of amorphous structure. ..............
RDF for the model of Bell and Dean for amorphous silica.

The first model of tetrahedraly bonded amorphous solid. .....
RDF for Polk’s model of amorphous silicon. ............
RDF for the amorphous germanium model of Steinhardt et al.. . .

Two examples of Guttman’s bond switches . ............

xi

87
89

91

94

95
96
97
98
100
101
102
104
105
105
108

113
115
116
118
119
121
126

6.8 Crystalline diamond cubic cell .....................
6.9 Keating potential . ..........................
6.10 Wooten—Weaire bond switch ......................
6.11 2d schematic representation of a topological bond switch ......
6.12 Pair Correlation Function for different amounts of randomization.

6.13 The Pair Correlation Function at three points in the annealing pro-

6.14 Wooten and Weaire PCF for a — Ge. ................
6.15 A computer—generated picture of the crystalline/ amorphous silicon

interface model . ...........................

7.1 A piece of computer generated amorphous a — C structure .....
7.2 Pair Distribution Functions for a — C' and a — Si models ......
7.3 Structure Factors for a - C and a — Si models. ..........

7.4 Comparison of the Structure Factor with experimental results. . .

8.1 Elastic moduli for random bond-depleted crystallyne lattice. . . .
8.2 Total strain energy as a function of mean coordination. ......
8.3 Jump from one to another metastable state. ............

8.4 Two different energy minima. ....................

127
130
131
133
136

168
174
176
177

8.5 Configuration of the two different states of the same diluted network. 179

8.6 Polymeric chains which give rise to the two states of the same
random network. ...........................
8.7 Bulk modulus B vs. mean coordination (r). ............

8.8 The optimal supercell size vs. mean coordination. .........

182
184

8.9 Evidence of the network collapse from the pair distribution function.185

xii

List of Tables

3.1

4.1

6.1

7.1
7.2

Summary of results for regular polygonal holes. .......... 42
Taylor series for circular inclusions. ................. 79
List of linear sizes of the super—cell box for Si and C ........ 128
The values of the parameters used in the models of C and Si. . . 149
Theistructural parameters for all C and Si models. ........ 153

xiii

Chapter 1

Introduction

In the last three decades or so, disordered properties of matter attracted a lot of
attention in the condensed—matter physics community. The development of new
mathematical and computational tools resulted in the major conceptual turning
point in our understanding of the properties of matter. Even though most phys-
ical systems that exist in nature are disordered, traditional physics was mainly
interested in simple, crystalline, regular or, in general, well ordered systems. This
was not only because it was easier to treat such ordered systems with the existing
mathematical tools (an enormous work had needed to be done to acquire satisfac-
tory understanding even of such simple systemsl), but also because the concept of
nature held by most physicists was still that of the older times, when the idea of
order, or rational ordering of nature, was the main driving force for physicists to
try to discover that hidden order of nature behind the seemingly chaotic appear-
ance of things. Generally accepted aesthetics, the measure of beauty of physical

laws which governs the Universe, was the aesthetics of order. Here is one example

[\D

of this acceptance in words of John Tyndall 1 in one of his public lectures:

To many here present a block of ice may seem of no more interest and
beauty than a block of glass: but in reality it bears the same relation to
glass that orchestral harmony does to the cries of the marketplace. The
ice is music, the glass is noise; the ice is order, the glass is confusion.
In the glass, molecular forces constitute an inextricably entangled skein;

in ice they are woven to a symmetric texture...

This attitude has been changed in the last three decades when many
difficult problems in disordered systems have become tractable.

In this thesis, I will present some of the new concepts which are used when
dealing with both continuous and discrete disordered systems. The emphasis will

be on two classes of media: random binary composites, and amorphous materials.

 

lJ. Tyndall, ”Heat = A Mode of Motion.” Longmans, Green, London, 1863. The reference
is taken from: Wooten and Weaire (1987).

Chapter 2

Continuous Disordered Media

2. 1 Introduction

In this chapter I will give a brief account of the research that has been done in
calculating the effective values of various physical quantities in disordered con-
tinuous media. The system in consideration is a continuous media host which
contains inclusions of various shapes, randomly distributed over the sample. The
inclusions are of macroscopic size, but small compared to the size of the sample.
They can have sharp corners (e.g. polygons in 2d), or they can be of round shape
(e.g. elliptical, or circular in 2d, or spherical in 3d). One wants to calculate the
effective quantities in different kind of problems: electrical conduction, magnetic
permeability, heat conduction, etc. All these problems are mathematically equiv-
alent, and I will speak only about the electrical conductivity calculation in such
systems. Instead of the language of conductivity adopted here, one could use
the language of dielectric constants. All results are straightforwardly translatable

from one language to another.

 

O~r1

 

     

 

 

 

Figure 2.1: 2d example of a random composite. The system has effective conduc-
tivity a, and contains regions with conductivities 00 (white) and 01 (black).

2.2 Effective Conductivity

The conductivity of random composites has long been a subject of interest, dating
back to the time of Maxwell (1873) and perhaps even further. The number of
exact results is limited, but these few exact results have a special importance in
testing and evaluating approximate theories and judging the accuracy of computer
simulations.

The system in consideration is schematically represented in Figure 2.1.
The black inclusions can in principle be of an arbitrary shape, but most of the
earlier work had been done for the case of spherical particles immersed in a matrix
of uniform conductivity thus creating a statistically homogeneous and isotropic

suspension. The problems of electrical conduction, electric permittivity, magnetic

permeability and heat conductivity are mathematically equivalent, and knowing
the solution of one of those problems, one can get the solution of any other one
by appropriately renaming the symbols appearing in it. Maxwell (1873) solved
the problem to first order of the concentration f, taking into account only the
induced dipole moments of the spherical particles. The effective conductivity of
the suspension he found, can be written in the following form

l=1+37f, (2.1)
0'

0

where f is the concentration of spheres and 7 is dimensionless parameter defined
as
01 — 00

= —. 2.2
7 01+200 ( )

The numerical factors 3 and 2 which respectively appear in Equations (2.1) and
(2.2), are actually d and d — 1, where d is the dimension of the system. We shall
come back to this when we discuss 2d case in details.

Lord Rayleigh (1892) took into account the contribution of induced oc-
tupole moments to the effective conductivity and thus showing that Maxwell’s
formula is accurate only for small volume fractions.

Peterson and Hermans (1969) calculated the effective dielectric constant
of a suspension of parallel cylinders and spheres to the order of f2, using the
method of reflections to solve the two—sphere problem. Jeffrey (1973) calculated
the average heat flux of the random suspension of spheres also to the order of
f2, using the general method deviced by Batchelor (1972) for the averaging of
interactions between spheres to obtain bulk properties of suspensions. To solve

the two—sphere problem Jeffrey used twin spherical expansions first formalized by

Ross (1968).

2.3 Spectral Functions

The use of the spectral function probably dates back to Foster (1924) for dis-
crete resistor networks. This approach was applied to single cubical inclusions
by Fuchs (1975). The general theory for continuous media was further developed
by Bergman (1978) and by Milton (1981). It was subsequently put on a more
rigorous basis by Golden and Papanicolaou (1983). Extensive numerical evalua-
tions of the spectral function in the low-concentration limit have been performed
by Felderhof and coworkers (1982, 1989, 1992). A discussion of the discrete re-
sistor case, especially near percolation, has been given by Straley (1979). Day
and Thorpe (1996) further developed the calculation of the spectral function for
discrete resistor networks, for arbitrary concentrations.

Despite this activity, spectral functions are not mu- 11 used in this area,
and the subject has become a bit of a backwater. This is unfortunate as this is
potentially a very powerful formalism, that could lead to important insights into
system behaviour. The attraction is that there is a single real function that char-
acterizes the behavior of the system, for a fixed geometry. Even if the (complex)
conductivities of the two components change, by varying the frequency or the
temperature, the same spectral function describes the response of the compos-
ite. Thus the spectral function is determined by the geometry of the two—phase
medium, and so it should be possible to find numerical algorithms to calculate
it directly and efficiently. This has not yet been possible, although significant
progress has been made recently by Day and Thorpe (1996).

It has been shown first by Bergman (1978, 1986) that for any two—phase

composite medium, provided Maxwell’s equations of electrostatics hold, the effec-

tive conductivity (or dielectric constant) can be written as

 

a 1 u clu
__ z 1 _ M, (2.3)
00 o t — u
where parameter t is defined as
00
t = . 2.4
00 __ 01 ( )

Equation (2.3) holds in any dimension d, and serves to define the spectral func-
tion g(u) as a transform on the conductivity 0, which is the function of the
material parameter t.

The spectral function 9(a) is a real, positive semidefinite function for
0 < u < 1, and zero outside this range. The spectral function is fully determined
by the random geometry of the two-phase medium. Quite generally, for a system

which on average is uniform the total weight of the spectral function is:

[019mm = f, (2.5)

where f is the concentration as before, or in other words, the volume fraction of
the phase 01 in the host 00.
If in addition the system on average is isotropic, then the first moment of

the spectral function is given by:

[01 ug(u)du = f—(l—dlfl, (2.6)

where d is the dimension.

It is convenient sometimes to define the normalized moments < u" > of

the spectral function by

fol < u" > g(u)du
1019006“

 

<u" >=

so that from (2.6) we have for the first moment

1-f

<u>=——-—. (.
d

[v
00
V

In the case when the inclusions are spheres (3d) or circles (2d), if we ap-
proximate the spectral function by a 6 function with weight and location such
that the sum rules (2.5) and (2.6) are satisfied, then (2.3) leads to the well known
Clausius—Mossotti or Maxwell Garnett formula [see the Equation (2.12) in Felder-

hof and Jones (1989) for 3d case]. Explicitly, in d dimensions,

 

QCMW) = f6(u+£(1- n) (2.9)
leads to
UCM _ df .
00 —1+—_1—f—dt° (2.10)

Instead of using t, we can define another material parametar 7, which will be later

used throughout the rest of the text. Again, in d dimensions,

 

0’1 — 00
= , 2.11
7 01 + (d —' l.)0'0 ( )
which in 2d reduces to (2.2)
From (2.11) and (2.4) follows
dt = 7—_l-. (2.12)

7
Substituting (2.12) in (2.10) we have another useful expression for the effective

conductivity of the mixture

3;: 1+(d-1)f7. (2.13)

0'0 1—f7

 

Assuming that the volume fraction f is small, we can expand (2.13) in Taylor

series and get the effective conductivity for spherical (circular) inclusions to the

first order in f

0
—=1+df7+~--. (2.14)

O

I will show later how one can go further in this expansion of the effective conduc-
tivity, to the second order in f in the case of circular inclusions.

It would be desirable to know the spectral function directly from the given
random geometry, and then calculate the effective conductivity from (2.3), but at
this stage of the theory, this is not possible. Instead, in order to collect information
about spectral functions, one must first calculate the effective conductivity, and
then extract the form of the spectral function from that, by taking the imaginary

part of the effective conductivity 0/0‘0. Indeed, if we change the sign in (2.3)

—=1+/ gm” (2.15)

u-t

 

and recall the well known formula from the theory of analytic functions (see e. g.

Mathews and Walker 1970)

SW .
u — (t) — ie —p/u——)—— 7rg(u), (2.16)

where so is the Cauchy principal value, or in the somewhat symbolic but useful

form

 

1

—— = ' — t . .
u—t-ie pu— +17r6(u ) (217)
Substituting (2.17) into (2.3), and separating the real and imaginary parts it

follows that imaginary part of the effective conductivity can be written as

Im (i) = 7rg(u), (2.18)

00

or solving for spectral function

g(u) = l1772(1) . (2.19)

71' 0'0

10

To calculate the spectral function one thus has to know the effective con-
ductivity. The effective conductivity can always be found in principle from the
induced dipole moment on the whole sample, that is produced by the induced
charge on the boundaries between the regions of 00 and 01, shown in Figure 2.1
(Hetherington and Thorpe 1992, Thorpe 1992). We shall see later how one can
calculate the induced dipole moment in various cases.

There is another useful form in which we can write Clausius-Mossotti

formula, in terms of the mean polarizability of a single inclusion. In d dimension

 

o - o n"_
0 = —*3, (2.20)

0' + (d -‘ l)0’0 d
where n is the number of the inclusions per unit volume V0 in d—dimesional space,

and 3 is the mean polarizability of a single inclusion. Solving (2.20) for 0/00, and

using the following relations

1'3 2 dV07 (2.21)

and

f = nVo, (2.22)

and one arrives back at the formula (2.13).

2.4 Two—dimensional Case

Since my work which deals with the continuous random media was done on two—
dimensional systems, I will now present some characteristics of the 2d systems,
as well as derive some formulas valid only in 2d which will be later extensively
used. By two-dimensional system we do not mean that the physical space of
the problem is two—dimensional, but only that there is no dependence on one

Cartesian coordinate, as for example in the case of aligned fibers.

11

2.4.1 Reciprocity Theorem

The reciprocity theorem states that if the conductivity of the two—phase medium
in 2d is 0(00, 01), and if the geometry is kept fixed but 0'0 and 01 are interchanged,

to give a conductivity 0(01. 00), then these two conductivities are related by
0(00701)U,(01900) = 0001. (32-23)

We are interested in a special case of (2.23) when one component, say 01 is present
in small amounts, so that the other one with 00 can be considered as a host. In
other words, this is to say that the volume (area) fraction of the component 01,
f, is small. Then, we can write an expression for the effective conductivity which

is, to first order in f, equivalent to (2.23)

 

a
— =1+7G(7)f+---. (2.24)
00
where
0'1 — 00
= , 2.25
7 01 + 0'0 ( )

which is the general expression (2.11) with d substituted by 2, and C(7) must
be an even function of 7. To see that, rewrite (2.24) for the two interchanged

geometries, and with the explicit a dependence on 00 and 0‘1

0(ao,al) = on + 007G(7)f, (2.26)
where 01 is present in small quantities, and

0(al,oo) = 01 — 017G(—7)f, (2.27)

where the geometry was inversed, so that 00 now is present in small quantities.

Multiplying (2.26) by (2.27), we see that, to first order in f, the reciprocity theorem

12

is restored, ifthe function G(7) is even, so that linear terms in f cancel each other
on the right—hand side, leaving the expression 0001 as in (2.23).

From the other side, writing the effective conductivity ad hoc in the form
(2.24), may seem to be too hard to guess a priori. To get a hint about the correct
form of the conductivity in the low concentration limit of one of the components,
we shall solve the problem first in the weak scattering limit, and then, generalize
it. The weak scattering limit occurs when the two conductivities 0‘1 and 00 do
not differ much, i.e. when the difference 01 — 00 is small quantity. Then, using
perturbation theory (Landau and Lifshitz 1960) to second order in 01 — 00, one

can solve the problem exactly to obtain (Landau and Lifshitz 1960)
a = '6 — —, (2.28)

where the bar represents an average. This formula is valid for 3d case, and it is
easy to see that in 2d case which we have here, we have to replace 3 by 2 in (2.28).

Explicitly,

(l-f)00+f01

02 = (1—f)a§+faf . (2.29)

QI
||

02—0' .
=‘— . 2.30
o a 25 ( )

 

Substituting the expressions (2.29) into (2.30) we get

f(1— f)(01 - 00V (2.31)

U=UO+f(01-UO)_ 2[0'0+f(0'1—0'0)].

Now, we expand the expression (2.31) in terms of small quantity A defined as

A z ”I " ”0 (2.32)

00

 

13

to get

 

f—2 f__2)\2 f2 A3
— =1 A — . 2.‘ :
00 + f 2 +— 2 + 2 ( 33)
Expressing /\ through 7 as defined 1n (2. 25), and expanding the new expression

assuming 7 is small quantity, and keeping only terms to first order in f we get
a 9 .9
E; = 1+ f7[2 — 0(7) + - - -] + 0(f‘). (2.34)

Comparing (2.34) to (2 .24), we see, first that in the weak scattering limit, the

function 0(7) is generally given as
G(‘7) = 2 + 0(72), (235)

no matter what particular shape inclusions have.

We shall come back to this when we present our results for the case of
circular inclusions where we shall include f2 terms in (2.24) using again the re-
strictions imposed by the reciprocity theorem.

The quantity 7 defined in (2.25), for a given finite conductivity of the
host 00 is always within the range [-1,+1]._ When 7 = -1 we have a hole as
inclusion, i.e. 01 = 0. In the other limiting case, when 7 = +1 we have a perfectly
conducting inclusion, i.e. 01 —> 00. Of course, when 7 = 0, that means that there

are no inclusions at all.

2.4.2 Effective Conductivity

In 2d it is convenient to change variables in the basic equation for the effective
conductivity (2.3). Using the variable 7 defined in (2.25) instead the variable t

defined in (2.4), and with r’ E u — % and h(:c’) E g(:1:’ + %), we may rewrite (2.3)

as

— =1+ 2 (2.36)

—- i1+27:”

14

where
A2 ('(I'WI' = f. (2 37)
and
i f2
j; r h(.r )dx = ——2-. (2.38)
The first moment (2"), defined generaly in (2.8) is now given by
1 - f 1 f .
(fl— d —§——? am)

for d = 2.
The spectral function is now a real positive function in the interval —% <

:r < l and zero outside this range. The sum rules (2.37) and (2.38) are, of course,

2’

special case of (2.5) and (2.6) in 2d.

2.5 Conclusion

In the next chapters I restrict myself on two—dimensional systems, and show how
one can calculate the spectral function in various cases, in the low concentration
limit. For the case of circular inclusions I will show in details how one can go one
step further from the low concentration limit to include pair terms in the induced
dipole moment, i.e. f2 terms in conductivity. I will explain the importance of the
conformal mapping techniques for the case of polygonal and star-shaped holes.
Finally, I will briefly explain the use of Fredholm equation for the general case
of polygonal inclusions, give an example of the discrete resistor network in the

low—concentration limit, and the results of the effective medium theory.

Chapter 3

Review of the Results for various

two—dimensional systems

3.1 Introduction

The general case of the inclusions of arbitrary 7 and of arbitrary shape is, of
course, very difficult. Practically the only solvable problems are the circular and
elliptical inclusions, and we were able to include the pair terms in the condictivity
only for circular inclusions. Even simple shapes of the inclusion like an equilateral
triangle make the problem too difficult if 7 is arbitrary. Only for the particular
values 7 = 21:1 can one relatively simply solve the problem for polygonal inclusions
using conformal mapping techniques. These limiting values of 7 correspond to a
hole (7 = —1) or superconducting inclusion (7 = 1). If 7 is arbitrary then one
has to solve numerically the Fredholm equation for the induced charge around the
polygonal inclusion (Hetherington & Thorpe 1992). Nevertheless, it is very useful
to have solution of the problems with the holes of various shapes. This has been

extensively studied by Thorpe (1992) who calculated the effective conductivity for

15

16

various shapes of inclusions with straight edges. Of course, when we fix the value
of 7, as for the holes or superconducting inclusions, our goal is only to calculate
the effective conductivity. The calculation of the spectral function, requires a
variable 7, which in most cases is an impossible task.

We shall first present an example where the exact solution exists, which
is for elliptical inclusions. Then we shall follow the work of Thorpe (1992) in pre-
senting the various results for different shapes of the inclusions when the straight
lines make their edges. For that purpose we shall introduce the conformal map-
ping technique, and in particular, Schwarz—Christoffel transformation which is

very useful for the mapping of the regions with straight edges.

3.2 Elliptical Inclusions

A useful example that can be solved exactly is provided by the case of elliptical
inclusions. The three—dimensional case of a dielectric ellipsoid in a uniform ex—
ternal electric field is solved in Landau & Lifshitz (1960. p. 44). If a, b and c
are the semi—axes of the ellipsoid, and if the external field is along a-axis (i.e.
:c—axis), then the field inside the ellipsoid is given by (Landau & Lifshitz 1960,

p.44, equation (8.8)):
1 + (2.} — 1) n.’

where 61 and co are dielectric permeabilities of the inclusion and of the host,

E;

(3.1)

 

respectively, and n3, is the depolarization coefi‘icient given by (Landau & Lifshitz

1960, p.26, equation (4.25))

n” _ 20 c o (s + a2)3/2(3 + (fill/209 + GUI/2,

 

17

where a, b and c are semi—axes of the ellipsoid. The integral in (3.2) is an elliptic
integral.

The 2d limit of the ellipse we get letting c ——+ co in (3.2). The depolar-
ization coefficient is then:

n —-1-ab‘/oO d8
I — 2 o (s+a2)3/2(s+b2)1/2'

 

(3.3)

The integral in (3.3) can be solved in terms of elementary functions to give

b
a+b’

 

(3.4)

n3:

where a and b are semi—major and semi—minor axes of .the ellipse respectively. If
the external field is along the b semi—axis, the depolarization coefficient is

a

ny—a+b.

 

(3.5)

Now, we can substitute (3.4) into the expression for the field inside the elliptical

inclusion (3.1) to get
_ c b '
1+ (2% ‘ 1) m

We see that the field inside the elliptical (as well as ellipsoidal) inclusion is uniform,

 

E;

(3.6)

which makes the exact solution possible in this case.

Imagine now that we have a small concentration of elliptical inclusions,
randomly distributed in the host, but all oriented along the external field E0. _
We can now consider the electric field averaged over volumes which are large
compared with the scale of the inhomogeneities. The mixture is a homogeneous
medium with respect to such an average field, and it can be characterised by an
effective dielectric permeability e. If E and D are the field and induction averaged

in this way, then, by definition of the effective dielectric permeability e:

D,B = (SE-1., (3.7)

18

where we assumed that external field is along the :r—axis.

To find 6 we write:
1 — — _
‘7- ](Dr — foEr) E DI — COEr- . (3-8)

The integrand in (3.8) is zero except within inclusions, where E, 2 E2, and

DJr = 61E; Taking the uniform inside field E;, outside of the integral, we get
f(€1 —' £0113; = 51‘ — €0Ex» (3.9)

where f is the volume fraction of the inclusions AV/V.

We have already found the field inside a single elliptical inclusion, E;,
which is placed in an external field E0 (3.6). Here, there are many elliptical inclu-
sions randomly dispersed in the host medium, but oriented in the same direction.
This mixture gives rise to a mean field which we can consider as being the true
external field with respect to each elliptical inclusion in the emulsion. Thus, using
(3.6), with the mean field Er as the external field. and substituting this expression

into (3.9), we get the proportionality coefficient between D: and Ex

6 61—60
_=1+___
co bel+aeo

a + b) f, (3.10)

or translated in terms of the conductivites,

2'1: .222. b
00 1+b01+aao(a+ )f. (3.11)

Equation (3.11) gives the effective conductivity of the mixture when the external
field is parallel to the semi—major axis a of the elliptical inclusions.

If the external field is aligned with the minor semi-axis of the ellipse b,
the effective conductivity is given by

:1. =1+ L10" (a + m. (3.12)

19

Formula (3.11) is given in the paper by Hetherington & Thorpe (1992), but with
a and b permuted. This becomes obvious if we consider the special case when the
inclusion is a hole, i.e. when 01 = 0 (7 = —1). Then from our formula (3.11) we
get

$711 a + b

:1—
00

 

f. (3.13)

where f is the area fraction of the ellipses and can be expressed as
f = ns, (3.14)

where n is the number of ellipses per unit area, each with area 3 E 7rab. Substitut-
ing this into (3.13) and finding the same hole—limit of (3.12) we get the following

two expressions

fl

= 1— n7rb(a + b),
00
ii = 1— n7ra(a + 5), (3-15)
0’0

which are correct [see equation (40) in the paper by Thorpe (1992)].
Averaging over the two orientations of the elliptical inclusions, which is
equivalent to an isotropic average, we obtain the conductivity in the form (2.24)

with
2

1<—)

Thus, for the elliptical inclusions, the function G(7) is obviously even, as it should

 

0(7) (3.16)

be.

Thus, for circular inclusions, to first order in volume fraction f, we have

the exact result
GU) = 2- (3.17)
This result is correct for all 7 in the case of circular inclusion, while (2.35) is

correct only in the weak scattering limit, but for any shape of inclusions.

20

With the use of (2.36) and (2.24) we can easily calculate the spectral
function h(.r) for this problem (we use the dummy variable .1: here instead of 23’).
Thus, for the randomly oriented elliptical inclusions, to the first order in f, the

spectral function is

f 1 a—b 1 a—b .
h(:1:)= 216(x_2(a+b)) +6($+2(a+b))j , (3.18)

where a and b are the major and minor semi—axes of the ellipse.

 

 

It is easy to check that the result (3.18) is compatible with the first sum
rule (2.37) but NOT with the second sum rule (2.38). This is because we did
not include pair terms representing interactions between elliptical inclusions. The
inclusion of the pair terms is not straightforward and the only case when this can
be done exactly is the case of circular inclusions, which I will discuss in detail
later.

The function (1(7) in (3.16) has two simple poles. In the limit of a slit
(b = 0), the two poles have equal weight and are at a: = :1:1 which are the outer
limits of the allowed spectral region. In the limit when the ellipses become circles
(a = b), the two poles come together to form a single pole at the origin with
weight f.

In the limit of a slit, note that perpendicular conductivity (3.15) becomes

U—L =1 — mraz, (3.19)
00

where a is the half length of the slit. We shall get this result again later, using the
conformal mapping technique.

Using (2.24) together with (3.16) we can calculate the effective conduc-
tivity in the limit of 11 randomly oriented slits, each of length 2a. The slit limit

(stick—like holes) is achieved as before with 0'1 = 0 and b —+ 0. And because

21

f E nrab for ellipses, the effective conductivity becomes

1 In” 2 3 90
a I Q .u
00 2 ( )

3.3 Conformal Mapping Technique

The conformal mapping technique is very powerful for solving many problems in
two—dimensions which involve Laplace’s equation (Churchill et al. 1974; Nehari
1952; Smythe 1939).

As we shall see it is very convenient to use the conformal mapping tech-
nique for the calculation of the effective conductivity when the inclusions are either
holes (7 = —1), or perfectly conducting (7 = 1). For arbitrary 7, and for inclu-
sions with straight edges, conformal mapping is not useful and one has to solve
Fredholm integral equation for the induced charge density around the inclusion.

The properties of the conformal mappings
w = w(z), where z = :r + iy w = u + iv, (3.21)

make them particularly useful for the solution of many problems in electrostatics,
magnetostatics, heat flow, hydrodynamics, elasticity, and other fields.
First, the real and imaginary coefficients of every analytic function of a

complex variable 2 are harmonic functions of :1: and y, i.e. they satisfy Laplace’s

equation
6% 6’21 62v 0222
Second, The Cauchy—Riemann conditions imply
a“ a” a“ a” = (Vu) . (Vv) = o, (3.23)

aa+aa

22

which means that the vectors Vu and V0 are perpendicular to each other. These
vectors are perpendicular to the curves u = constant and v = constant, respec-
tively, so that the curves themselves are mutually perpendicular. Usually, one
takes one of the functions, u or v, for the potential function, while the other one
is then the stream function. The choice for the potential is arbitrary and depends
on the given problem. The function w = u + iv is then a complex potential.
Third, the transformation w = w(z) is conformal, which means that the
angle of intersection of two curves in the z—plane stays the same in the w~plane

after the mapping. This follows from the analytic property
w(z) — w(zo) z w (zo)(z — :0), (3.24)

when 2 is close to 2:0. The mapping is not conformal at a singular point of the
transformation, that is, a point where the function w(z) or its inverse w(z) is not
analytic.

The power of conformal mapping technique is that it is usually possible
to find such a transformation which maps a complicated problem in the z—plane,
into a very simple problem in w-plane which can be solved by inspection. One
looks for a transformation which produces much a simpler boundary cOnfiguration

in the w—plane.

3.3.1 Example

As a simple, but instructive example we can use conformal mapping to find the
electrostatic potential of the charged rectangular conducting boundary as in Fig-
ure 3.1.

Let assume such boundary conditions that potential 4) equals to zero along

the axes a: = 0 and y = 0. The transformation

23

 

Figure 3.1: Application of the conformal mapping w(z) = .22. Inside corner of

conducting sheet (a) maps into the upper w half—plane (b).

24

n‘

)= 2" (3.25)

w (

transforms the inside corner in Figure 3.1(a) into the upper half of the w—plane,
in Figure 3.1(b). In other words, the rectangular boundary is transformed into
the real 11 axis in the w—plane. As before, 2 = .r + iy, and w = u + iv. The
new problem in the w—plane is now trivial. The electric field in the w—plane is
uniform, and it is obvious that the solution for the complex potential Q in the

w—plane is

Q=—Cw=—Cu—in=w+i¢, (3.26)

from which we see that
w = -Cu <2 = —Cv, (3.27)

where C is a real constant We choose the imaginary part of the complex potential
Q as the potential function, and the real part is now the stream function which
describes the electric field lines. It satisfies transformed boundary conditions that
it is zero along the real u—axis, i.e. when v = 0.

Using (3.25) the complex potential becomes
Q = —sz,
= —C(;z:2 — y2 + i2xy). (3.28)
(3.29)
Thus, we now have the solution of our original potential problem in the z—plane

a = 4ny (3.30)

with the constant C = E = T/eo, where E is the homogeneous electric field in the
w—plane, and r is the charge density on the surface of the conducting lower half

of the w—plane, and e is the dielectric permeability.

 

An

 

AN

 

 

 

 

 

 

z - PLANE =5 w- Fume

Figure 3.2: Schwarz—Christoffel transformation maps polygonal vertices A), onto
pornts ak on the real axis in the w—plane.

3.3.2 Schwarz—Christoffel Transformation

The Schwarz—Christofell transformation is very powerful when the region to be
mapped has the boundaries which consist entirely of line segments (see Silvermann
1974). The obvious example is a polygon with n vertices. The Schwarz—Christoffel
transformation maps the interior or exterior of a polygon in z-plane, onto the
upper half of the w—plane, as shown in Figure 3.2.

For our purposes it is enough to start from the general form of the trans-
formation and to understand its applications in different cases. Generally the

Schwarz—Christroffel transformation is given by the following integral

2 = c Cw(w — a1)2«L"1(w — aft-1 - . . (w — aft-law + (:2, (3.31)
1

where the real numbers on are the images of the polygonal vertices An, and an

 

 

 

 

.8
A
v8
V;

 

Z - PLANE w- prune

Figure 3.3: Mapping of the half—strip —7r/2 < 36(2) < 7r/2 onto the upper half-
plane Im(w) > 0.

are the internal angles of the polygon, as shown in Figure 3.2. With no loss of
generality we can fix the constant C1 to zero. Namely, changing C1 amounts
only to changing C2. There is a theorem (see Silvermann 1974) which proves
that specification of three points (11, a2, a3 on the real u—axis corresponding to
three vertices A1, A2, A3 of the polygon automatically determines the remain-
ing points a4, - - - ,an and the constants C and 02. The actual determination of
a4, - - - ,an, C, 02 is the main difficulty when using the Schwarz—Christoffel trans-
formation.

It can be proved that the main formula (3.31) still holds even if one of
the real points ak is infinite, and also if one of the vertices of the polygon A), is

r
infinite. In those cases, one has simply to include only remaining (finite) points

27

ak, or vertices Ak. If one of the vertices of the polygon is at infinity, we have the
case of an open region bounded by the straight segments in Figure 3.3.

The internal angle on, is defined as the angle which belongs to the region
mapped, which is, by convention, always the region to the left as we pass along
the polygonal boundary in the given direction. In Figure 3.3 we pass along the
rectangular boundary of the half—strip in the anti—clockwise direction, so that the
desired mapping region is to the left. We can regard the half—strip as a ” degenerate
triangle” with a vertex at infinity. In practice wejust have to include the two finite
vertices A1 and A2 into transformation (3.31). We choose points a1 and oz to be
—1 and +1 respectively. The internal angles 01 and 012 are both equal to 7r/2.

Substituting these values into (3.31) we get

2 = C'/Ow(w+1)_l/2(w—1)—1/2dw+C1, (3.32)
01'

w d
z=C/()¢—1:—w;w_7+C1=Carcsinw+Cl. (3.33)

To determine the constants C and C1 we note that the real points —-7r / 2 and +7r / 2

in the z-plane, map into the real points —1 and +1 in the w—plane respectively.

In other words from (3.33)

: _Cg '1' Cl,

wlamla

= 0% + cl, (3.34)
which gives C = 1, CI = 0. Thus the mapping finally becomes
2 = arcsin w, (3.35)

with inverse

w = sin 2, (3.36)

28

01‘

u + iv 2: sin(r + iy),
= sin(sc) cos(iy) + cos(r) sin(iy), (3.37)

= sin(r) cosh(y) + cos(;r) sinh(y).
Separating the real and imaginary part in the last formula we get

u = sin(r) cosh(y),

v = cos(a:) sinh(y). (3.38)

Dividing the first equation in (3.38) by cosh(y), and the second by sinh(y), squar-
ing and then adding them, we get

u2 v2

cosh2(y) + sinh2(y)

 

= 1. (3.39)

Similarly, dividing by sin(x) and cos(:c), squaring and substracting, we get

u2 v2

 

sin2(:r) cosz(:r)

 

= 1. (3.40)

We see that according to (3.39) the curves in the w—plane on which y = constant
are confocal ellipses, while looking at (3.40) we see that the curves on which
a: = constant are confocal hyperbolas in w—plane. In this way, if we consider a

uniform complex potential Q in the z—plane
Q = —|E|z, (3.41)

depending on the choice of the potential and stream function in the z—plane, we
obtain the solution for either the potential around a charged conducting strip of
width 2 which is perpendicular to the w—plane, or the potential due to the slot in

a conducting sheet.

29

Going back to the Figure 3.1 we can see that the transformation w(z) = 2:2
is just a special case of Schwarz-Christoffel transformation. Indeed, considering
the inside corner in Figure 3.1 as a degenerate polygon, and applying the equation

(3.31) with a1 = 0, and 011 = 7r/2, we get

w dw ‘
z _ CA W + 0,. (3.42)

Because 0 in the z—plane maps onto 0 in the w—plane, we get C1 = 0. If we choose

that 1 in the z—plane maps onto 1 in the w-plane, we get
1 dw
C / —,
o (A?
1 2 2C, (3.43)

from which we have C = 1/2. Substituting constants we found into (3.42), solving

the integral. and inverting the equation, we get
w = 22, (3.44)

which is the transformation we used before.

3.4 Circular Hole

The example of the circular hole is very useful as a starting point in treating
general case of a polygonal inclusion. It has been shown (Thorpe 1992) that the.
effective conductivity of a sheet containing a few randomly distributed holes can

be written to the first order in f as
‘7 2
—=1—af+0(f). (3.45)
00

In all cases the total induced charge around the hole is zero. In the case of

circular-hole inclusion the induced charge only has a dipole term, that leads to a

30

 

 

 

 

 

 

 

 

 

 

 

 

41
/\ -
g x
- o. o a. -7 6 i 71.
Z- PLAN ‘
E W- PLANE

Figure 3.4: A sketch of the flow around the circular hole in the upper half of the
z—plane. The upper half of the z—plane outside the circle maps onto the upper
half of the w—plane. z = in maps into w = :l:1, and z = ia maps into w = 0.

potential that decays as r"1 where r is the distance from the hole. For holes of
other shapes there are also an infinite number of higher-order multipole moments,
but only the dipole moment, i.e. the dominant term in the far field, contributes
to the coefficient a in (3.45) (Thorpe 1992).

Thus, our goal is to find the reduction in conductivity caused by the hole,
i.e. the coefficient a. Imagine that there is a circular hole present in a conducting .
sheet, as shown in Figure 3.4. Due to the symmetry of the hole, it is enough to
solve for the flow in the upper half of the z—plane.

The conformal mapping between the upper halves of the z— and w-planes

(see Figure 3.4) is given by ,
w = 1;.(3 + 2), (3.46)

a Z

31

where a is the radius of the hole. In the w—plane the complex potential Q is given
as

(b = -E0w '2 —E0u — iEOU = CD + iL‘, (3.47)

where E0 is the electric field far from the hole. Substituting the conformal mapping
(3.46) in (3.47) we get

E 2
<1> = -—2—°: (1+ 9-3). (3.48)

The actual potential o is taken here as the real part of Q, and noting that as

before 2 = .r + iy = rcos 0 + ir sin 0, the potential (b is given as
E 2
d) = —-2—0rcosd (1+ :3),
I (12
at = —E0rcosd (1+ 7) , (3.49)
r

where E; is the new scaling factor for the electric field at infinity. Using E. =
-8¢/8r we can calculate the radial component of the electric field in the z—plane
which is

2
E, = El, c050 (1 - L) . (3.50)

There is obviously no normal current flow into the hole, and because of symme-
try there is also no normal current flow into the .r—axis of the z—plane. These
simple Neumann boundary conditions (Churchill et al. 1974) make this problem

straightforward to solve. We see from (3.50) that
Ema = 0, (3.51)

as it should be at the circular boundary.
We shall find the effective conductivity of the system by calculating the
change in the Joule heating upon the introduction of the hole (see Thorpe 1992).

Imagine a large concentric circle of radius R drawn around the circular hole. First,

by

.
W E

32

if there is no hole in the system the potential at any point on that circle is given
by
cbno hole = —E0RCOS 9, (3.52)

while in the presence of the hole, using (3.49) with r = R, the potential is

. I a2
(”with hole = -EoRcos¢9 (1 + 7,3) - (23-53)

Assuming that far from the hole potential is fixed, no matter whether the hole is

present or not, we can equate (3.52) and (3.53) to get

13.. = E30 + f), (3.54)
where
no2
f -— fi—R; (3.55)

is the area fraction of the hole.

Consider the volume V inside the large circle with radius R, but excluding
the volume of the circular inclusion itself. The Joule heating ,7 within that volume
is given by

J = /V(j-E)dV, (3.56)

where j is the current density, and E is the electric field, integrated over the whole

volume. Using the known relations E = -V¢, V - j = 0, and
WM) = ¢Vj +jV¢, (3.57)

we get for Joule heating
j = _ [V v . (j¢)dV. (3.58)

Transforming the volume integral into a surface one we finally get

J = - [5 3(5 - n)dS, (3.59)

.
u"
n L!

Be

33

where the integration goes over the outer surface only because there is no current
flow into the inner surface, and n is the unit normal on the surface element on
the outer surface.

Now, the change in the Joule heating upon the introduction of the hole

can be written as

M = —/56(6j . n)d5. (3.60)

Usingj : 00E, and noticing that 6j - n = aoAEr, with the help of (3.50) we get
(Sj-n = —-‘20E(',fcos 0, (3.61)

where 00 is the conductivity of the host. Substituting (3.61) into (3.60), using
(3.52) and performing the integration, we get the final expression for the change

in the Joule heating upon the introduction of the hole
6,7 = —2E0Eg,ao7ra2. (3.62)
Balancing the energies we get
NR20E3 = szaoEg — 2E0Eéaoa2, (3.63)
from which using (3.54) we get

i =1-— 2mm2 =1— 2f, (3.64)

00

where n is the number of inclusions per unit area.
The factor a2 in (3.64) comes from the coefficient in the dipole term 1/2
of the conformal mapping in (3.48).

Now, if we assume that the conformal mapping is

z=c[w+g+o(i—,)], (3-65)

the:

if; ll
9': ‘

34

then inverting (3.65)

 

~
~
(j ~ 2
" 4'

.4672 1
w = — [1— + 00 (7)], (3.66)

and comparing (3.66) with (3.48), we see that, in general, instead of having a'2
in the expression for conductivity (3.64) we have the product —AC'2. Thus the

general result for the conductivity is

E. = 1+ 265/162 + 0(62). (3.67)

00

This simple result is possible only because the higher order multipole terms in-
tegrate out to zero, and only the dipole term contributes to the conductivity
(Thorpe 1992). The formulas (3.65) and (3.67) are very useful for calculating the

effective conductivity when polygonal holes are present.

3.5 Slit

We shall apply the Schwarz-Christoffel transformation (3.31 to the slit that lies

in the z-plane as shown in Figure 3.5, and use the result (3.67) of the previous

section to find the effective conductivity when the current flow is perpendicular

to the slit. The slit is of the length 2a and only the upper half of it is shown in

Figure 3.5. Due to the symmetry it is enough to solve for the flow in the upper

half of the z-plane. The current flow is perpendicular to the axis of the slit.
Applying the Schawrz—Christoffel transformation (3.31) we get

2 = C C(w+1)zir3"1(w)27r£'l(w—1)£vil’ldw+C'2,
l

= C —-'5’——dw+cz.
01 1122-].

= C\/w2 - 1+ Cg. (3.68)

35

216'

 

   

“/2.

v3

4V

 

vQ

 

 

o - 4
z- PLANE

o 4
va— Pmue

Figure 3.5: Mapping the upper half of the z—plane with a slit onto the upper half

of the w-plane.

36

If the point w = 1 in the w—plane corresponds to the point z = 0 in the :—plane,
we see from (3.68) that C2 = 0. From the other side, if the point ia in the :—plane

maps into w = 0, it follows that C = a, and the mapping becomes
2 = a w2 — 1. (3.69)

Expanding the mapping (3.68) for large w we get

A 1
z=C[w+—+O(—2)]. (3.70)
w w
Comparing (3.68) with (3.67), and using C = a, we get for the conductivity

2; =1— mraz, (3.71)

00

which is the same result as (3.19) we got before as a limit from the solution for

the elliptical inclusion.

3.6 Effective Conductivity for Polygonal Holes

The effective conductivity for the systems with polygonal holes was studied by
Thorpe (1992), and in this section we will present the main results of this work.
For regular n—gons it is not necessary to do any rotational averaging because the
symmetry of regular n-gon is high enough that the conductivity tensor reduces to _
a constant 0 times the unit tensor. In all cases there is a rotation axis of order 3
or higher, which is the necessary condition for isotropy. In the case of anisotropic
inclusions like diamond-shaped or rectangular, one has to calculate separately the
parallel and the perpendicular effective conductivity, and then to do the rotational
averaging. In all cases there is a reflection plane and it is centered on the art-axis,

so that it is enough to solve for the flow in the upper half—plane.

37

 

 

 

 

a .\\\

 

 

 

\\

W plane

Figure 3.6: The region outside the upper half of an equilateral triangle in the 2—
plane maps into the upper half of the w—plane. The three vertices of the triangle,

0, ax/3/2 + ia/2, a\/3/2 map into the points 0, 1, and k (k > 1) in the w-plane,

respectively.

3.6.1 Equilateral Triangle

The equilateral triangle of side a in the z-plane is positioned as shown in Figure
3.6.
Using the Schwarz—Christofell transformation (3.31) for this triangle, the

conformal mapping becomes

10 __ 2/3
2:0 (w 1)

c. tux/6(6) _ szdw + 02. (3.72)

 

Equating the distances between corresponding points in z— and w—planes we can

'P‘
W}
:0

l

4 V0

i)
\L

38

obtain the relevant constants k and C in (3.72).

k _1) 2/3 (1_ ,)2/3
a = QC]1 (‘0 w=C/l( U“) dw (3.73)

w1/6(k _ w) )1/2d .(wl/G k _ w) )1/2

 

 

Equation (3.73) is an equation for k. Calculating the integrals on both sides (in
terms of Beta and Hypergeometric functions) gives 1: : 4/3, and the right—hand
side of (3.73) gives

8772

=CW»

where F(%) is a gamma function (Abramowitz & Stegun 1970). Expanding the

(3.71)

conformal mapping (3.72) for large w we get

 

k2 k-—--
Z=C[w +(%k—§)lnw—8— w 9+O(lw 2)+const.]. (3.75)

Substituting k = 4/3 in (3.75) the term in lnw vanishes as it should, because
the net induced charge around the inclusion is zero. Comparing the expansion
(3.75) with the general form of expansion (3.65) gives A = —%. Using this result
together with (3.74) and substituting the obtained expressions for A and C into

the general form for conductivity (3.67), we get

0' 736(6) ..
—=1—a3f=1————f=1—(2.5811...)f, (3.16)

0'0 8713

where the area fraction f is here f = inx/3-a2 (Thorpe 1992).

3.6.2 Square

The upper half of the square is shown in Figure 3.7. The reflection plane is again
positioned along the x—axis. The three vertices of the square in the z—plane,
—7“5, 0 and 1% map into the points —1, 0 and 1 in the w—plane, respectively.

Schwarz—Christofell transformation (3.31) gives the following conformal mapping
1
2

z=C Cw( w )idw+02. (3.77)

 

w2—1

39

 

 

 

 

3 l)“
1' a/vz
3‘ u
-°/v: 0 an; -4 o 4
2 - PLANE w- PLANE

Figure 3.7: The region outside the upper half of a square in the z—plane maps
into the upper half of the w—plane. The three vertices of the square, —a/\/2, 0,
and a/\/2, map into the points —1, 0, and 1 in the w—plane, respectively.

40
As before, equating the distances between the corresponding points in :-

and w—plane we get
1

 

1 102 f
a = 67/0 (1_w2) dw, (3.73)
and solving the integral in (3.78) we get the constant C
27731
a = —. (3.79)
1‘20)
Expanding the conformal mapping (3.77) for large 10 leads to
~-C[ 1+0(1))+ t 380
.-. — w 4w w? cons ., ( . )

which gives for the dipole term coefficient A = —i. Using this and (3.79), and
then comparing (3.80) with (3.67) we get (Thorpe 1992)

r4 i-
53:1_a4f=1——8—7E§2f=1—(2.1844...)f, (3.81)

with f = naz.

Generally, for a regular N -gon, the conductivity can be written as

i =1—aNf+0(f2), (3.32)

00

where we want to calculate an.

3.6.3 Regular Polygon

In the case of a regular N —gon it is more convenient to use a different mapping
than in the previous paragraphs if we want to find the general expression for 0m
in (3.82). We use the transformation (Thorpe 1992) which maps the exterior of

a regular n—gon in the z—plane into the exterior of a unit circle in the w—plane

41

(Muskhelishvili 1963). The N vertices of the N—gon map into N points on the

unit circle at the N roots of unity e.2:p(27rik/N). The mapping is

w 1 73,7
z = C (1— —) dw + 02. (3.83)

C1 102
Expanding the mapping (3.83) for large w gives

C 1 + 2 + N 7 2 + + t 3 ‘4
Z = 1U) . . . . . . ,
1V(1V — 1)w‘v 1V2(1V _ 1).w2N C0118 ’ ( b )

 

 

where C = 1 if potential is to be the same far from the N —gon and the circle. There
is no dipole term in (3.84), so one cannot use (3.67) here. Instead, one should
use the mapping (3.83) itself. Equating the distance between the corresponding

points (3.83) gives for the side a of the N-gon

,, 2 I"? l + .1-
a =/ (2515 0)T’~7d6 = 277M. (3.85)
. o 1“ (1 + 7%,)
The ratio of the area of the N —-gon to the area of the circle is
a2 77
4—77 cot (N) . (3.86)

Using the expression for the conductivity for the circle (3.45) and comparing with
(3.82) (f = ns) we find that the coefficient cm is just twice of the inverse of (3.86)

tan(§) I“(fi)
ON: 27rN 172(737).

 

(3.87)

The general result for my (3.87), of course, reproduces the previous results for the
equilateral triangle (N = 3), and for the square (N = 4), which can be seen by
calculating the right-hand side of (3.87) using the relations between the gamma
functions with different arguments (Abramowitz & Stegun 1970). It can be also
checked that in the limit as N —+ co, the equation (3.87) gives doc = 2, as it must

be for the circle.

la

(1

Table 3.1: Showing the exact results for the coefficient ow in the conductivity
o/oo- —— 1 — aNf + 0(f2) for various regular polygonal holes from (3. 87).

shape n ow

 

triangle 3 2.5811...

square 4 2.1884. . .
pentagon 5 2.0878. . .
hexagon 6 2.0486. . .
septagon 7 2.0298. . .
octagon 8 2.0197. . .

9 2.0137. . .
10 2.0099. . .
11 2.0074...
12 2.0057...
circle 00 2

 

 

 

Summary of the results for the effective conductivity in the case of regular
polygonal holes given generally by formula (3.82) with an given by (3.87), is shown
in the Table 3.1 (Thorpe 1992).

3.6.4 Diamond Shape

As an example of the anisotropic holes we shall present here results for the inclu-
sion of diamond shape. In the similar manner, finding the appropriate conformal
mapping, one can solve the problem of the rectangular hole (Thorpe 1992).

In the case of a diamond shaped inclusion, one has to treat separately the

43

 

 

 

A
f

-a O a. -1 o 4
Z - PLANE w - PLANE.

Figure 3.8: The region outside the upper half of a diamond in the z—plane maps
into the upper hal of the w—plane. The three vertices of the diamond, —a, ib,
and a, map into the points —1, 0, and 1 in the w—plane, respectively.

two limiting orientations of the diamond. If the longer diagonal of the diamond
is 2a, and the shorter 2b, then we first treat the case when the longer diagonal is
aligned with x-axis, i.e. when the current flow is parallel to the longer dimension
of diamond. When the current flow is perpendicular to the longer dimension, we
can get the result by interchanging a and b in the result for parallel orientation.
The parallel orientation of diamond is shown in Figure 3.8.

Writing the Schwarz—Christoffel transformation for this case we get the

conformal mapping as

 

w w2 i
Z = C ( ) dw + Cg, (3.88)

C; 1122—1

where the angle 0 is given by tan 0 = b/a. The only relevant constant C is obtained

by 61

11.6

1th

44

by equating the distances between corresponding points to get

 

1 ,2 %
Val-Fl)2 = C/ ( u ) dw,
0

1_ .wz
F(1-%arctan(%)),1‘(1—%arctan(§))

= C fi . (3.89)

The dipole term coefficient can be found as before by expanding the mapping

 

(3.88) for large w
0 1
z = C (w — — + O (7) + const.), (3.90)
w

which gives A = —6/77. Combining these-results and using (3.67) leads to the
expression for the effective conductivity

91 :1_ fF”, (3.91)

00
where

7r(a/b+ b/a)arctan(b/a)
F = ..
ll P2(1—(1/7r)a.rctan(a/b))[‘2(1_(1/7r)arctan(b/a))a (3 92)

 

and the area fraction is f = 2nab where n is the number of inclusion per unit area.
When the diamond is oriented perpendicular to the :r—axis, the conductivity 01
can be obtained by interchanging a and b in (3.92).

The results for the case of rectangular inclusions are also obtained by
Thorpe (1992). In all cases the results obtained for holes can be also applied
to superconducting inclusions using the reciprocity theorem (Keller 1964). The
reciprocal problems for all examples in this section have the hole replaced by a su-
perconducting inclusion, and the equipotential lines become the current flow lines.
The results for the reciprocal problem are the same, except that the conductivity
is replaced by the resistivity and the sample is rotated by 77/2. This rotation

matters only for the anisotropic inclusions like diamond, and there one should

fig

13'.

45

 

Figure 3.9: The geometry of an (n = 3) star. The shape is determined by the
ratio of the radii of the inner circle b to outer circle a (Thorpe 1992).

interchange the subscripts for parallel II and perpendicular .1. in the expressions

for conductivity.

3.6.5 Star-shaped Holes

The problem of a N-pointed star-shaped hole was solved by Hetherington &
Thorpe (1992). An example of the 3—pointed star is shown in Figure 3.9.

The conformalimapping of the star onto the unit circle, such that N -
points around the outer circle of radius a and N points around the inner circle of
radius b, are mapped onto the 2N points uniformly spaced around the unit circle

in the w—plane, is

w w” _ 1 1-20/1r w” l-29/1r—2/N
= d C , 3.93
2 C 01 ( w” ) w” +1 w + 2 ( )

 

 

11111?

poi

all

46

where
bsin(7r/N)
a — bcos(7r/N)'

 

tanO =

(3.91)

The cute angle of the star is 20. Equating the distance between corresponding

pomts we get

 

' l _.2/n I‘(—{;+07r)
bsm (1V) —2 77F J§)F(%)i (3.95)

The area of the star is

N112 sin(7r/N) sin(7r/N + 0)

 

 

- (
sin0 (3'76)
and comparing with (3.82) we get
2 area circle 2sin(7r/N) sin0 3 9"
a. = -— = . .
‘V area star 24/NN77 sin(7r/N + 9) ( i)

In terms of the number of inclusions n : f/area star, the conductivity is given

by
2

-U——1—27rn aF(1—1/n-0/7r)
ao " 22/nr(1—1/n)r(1—6/7r) '

 

(3.98)

When 0 = 0, the star reduces to n sticks which are pointing outwards radially

from a central point. The conductivity (3.98) reduces to

o 217a2N

_=1_

0'0 271-

 

(3.99)

When n = 2 the result (3.99) reduces to the already obtained result (3.20) for a
dilute collection of randomly oriented sticks, each of length 2a (Hetherington and

Thorpe 1992; Thorpe 1992).

47
3.7 General Case of Inclusions with Sharp Cor-

ners

Results for polygonal inclusions presented in previous sections were possible to
obtain only for the limiting case when inclusions are either holes (7 = —1) or
perfectly conducting (7 = +1). In this case only it is possible to use the conformal
mapping technique to solve the problem (Thorpe 1992).

In order to get results for a general polygonal inclusion (arbitrary 7), it is
necessary to use a numerical method. The induced charge p(s) on the boundary

of the polygon is given by the Fredholm equation (Hetherington & Thorpe 1992)

p(s )=27n(s )2E0+ 17%) v 16 | r(s )— 1"(s’) | p(s')ds', (3.100)

where the integration is around the polygonal boundary. The vector r(s) goes
to a point 3 on the polygonal boundary, which has a unit normal 13(3), and .3
measures the distance around the perimeter to the point r. The external applied
field is E0 and the logarithmic term is the 2d Coulomb propagator. Equation
(3.100) determines the induced charge and hence the dipole moment from which
we can find 0(7) in the expression for conductivity (2.24). In solving (3.100),
great care must be taken at the corners of the polygons where the current density
is singular (see Hetherington & Thorpe 1992). It has been found by Hetherington
and Thorpe (1992) that the results for regular polygonal inclusions with arbitrary
7 can be well fit by the interpolation formula

2

0(7) = W, (3-101)
where the constant
c = 1 — —42- (3.102)

48

is chosen to get the correct result at 7 2 21:1. The interpolation formula (3.101)
was found to be quite accurate for all real 72 S 1, which covers the region where
both the conductivities 00 and 01 are real and positive. The form (3.101) is only
an interpolation formula and it would be misleading to conclude that there are
branch cuts at 7 = :l:1/c (which corresponds to .2: = :tc/2 = 21:0.316), which are
outside the region for which (3.101) is valid. Hetherington and Thorpe attempted
to find the spectral function directly by examining the simple poles that occur
on the real axis in the complex 7 plane when the charge density in (3.100) is
discretised. The result for the equilateral triangle is shown in Figure 3.10. It
consists of 64 poles on the real :r-axis that have been. broadened by convolution
with a Lorentzian (Figure 3.10 solid curve). It can be seen that the spectral
function h(:r) is even as expected. When the interpolation formula (3.101) is used
to derive the spectral function a similar curve, broadened with a Lorentzian, is
obtained (Figure 3.10). The similarity of these two results gives some confidence
in the general shape of the spectral function. This shape is very different from the
result for the ellipse, and is broad, rather than having one or two delta functions.

If the same procedure was followed for regular polygons, it is reasonable to
expect the spectral function to become narrower, with the spectral weight confined
to I 3: IS l/N, as the number of sides N increases, and therefore approaching a
delta function at the origin in the circle limit (N -> 00). The spectral function is‘

determined by the spectrum of the operator
f1(s) . Vln | r(s) — r’(s’) I . (3.103)

The only parameter in this operator is the boundary r(s). A method of going
directly from r(s) to the spectral function would be most attractive. Furthermore,

more work is also needed on the nature of the singularities at the edges of the

49

 

 

 

 

I I I II I I I I I r r I I I I I II I I I
I , l . l I , 1 .
2.0 — . Triangular inclusions I -—
1- r
_ 1 , ‘ I ‘
_ l, '\ '1 .
1: - '.‘ ' ’ "
V l- ‘ I 7
a D ' i ’ 3' d
a . 1' ‘- ’
8 - ‘~ ' ‘
a “ ’ ‘
n \ I
'3 ' -
3 - 1 1 7
3. t I 1 7
m 1 1 ‘
—
b . 7
D | I «I
I l
00 1 l l 11 1 L L L L L 1 1 1 L 1 l 1 1 11 1 l 1 L

Figure 3.10: Spectral function for triangular inclusions in the dilute limit. The
area under the curves is f but is shown as unity. The solid curve is from the
numerical solution of the integral equation (3.100). The dashed lines are specu-
lated spectral limits at a: = 21:1 / 3. The dot—dash curve is the spectral function
derived from the interpolation formula (3.101), which has band edges at :l:c where

c = 0.316. Both results are broadened with a Lorentzian with a full width of 0.06
at half the maximum height (Thorpe, Djordjevié and Hetherington 1994).

50

spectral function which is unknown at present.

3.8 Random Resistor Network

All our previous discussion has referred to the continuum problem. Random re—
sistor network is an example of a discrete random media for which one can apply
the general method developed in Chapter 2, to calculate the effective conductin
ity and the spectral function in the low concentration limit. Basic equations for
the effective conductivity and spectral function (2.24)—(2.3) still hold, if the area
fraction f is now the fraction of the resistors (conductances) that are present in
the network (Foster 1924; Straley 1979). This approach has been criticized by
Straley (1979) as not being much of use around percolation, which occurs for
larger concentration of conductances f. Recently, Thorpe and coworkers (1994)
calculated the spectral function for the random resistor network in the low con-
centration limit. They also obtained the result for effective medium theory, for the
square net and for the whole range of concentration f. Finally, Day and Thorpe
(1996) presented a derivation of the spectral function which identifies it as a den-
sity of states. They did direct numerical calculations of the spectral function of
two dimensional square lattice for the whole range of concentartion f, applying
the Y — A transformation as an algorithm to obtain the effective conductivity of
the random network. The analytic expressions for the first five moments of the
spectral function are also obtained.

The discrete random system in question, say the square lattice, can be
created in two different ways, i.e. with two different distribution geometries.
In the first geometry each bond is randomly assigned a conductance 0'0 with

probability 1 — f, and a conductance 01 with the probability f. In this way the

51

random system is created by random bond substitution. In the second geometry,
which is referred to as random site substitution, each site is assigned label A with
, probability 1 - f and label B with probability f. If the two neighboring sites are
labelled differently, i.e. if they create bond A — B, then half of that bond, which
is adjacent to site A, is assigned a conductance 00, while the other half, which is
adjacent to B, is assigned a conductance 01. Thus the bond A - B consists of
the two conductances 00 and 01 connected in series. Thus the conductance of the
A — B bond is or), : cool/(00 + 01).

It is convenient in this section to use the spectral function which is con-
fined in the interval [0,1] instead in the symmetric interval [—1/2,+1/2]. The
effective conductivity is then given by (2.3) which we shall rewrite here in the

more explicit form

 

t, 1 h , d -

7< f) =1_/ _(_x_f_>_~’l:. (3,104)
o t — a:

The spectral function h(a:, f) (we show here the explicit dependence of the spectral

function on concentration f) can be calculated by

h(t,f) = 1lim 1m [MI—l) , (3.105)

77 we 0'0
where t and e are real, and the variable t is given as before by (2.4). Thus, to
obtain the spectral function one has to calculate the effective conductivity for
complex t.

Unlike in the continuum problem, here for a given lattice substitution
type, the geometry is completely determined by the fraction of the second phase

f, and we have the exchange relation

”(00,011f)=0(01,00,1—f)- (3-106)

which connects the conductivities of the original problem and the one where the

52

roles of the host and inclusions are inverted. Similarly, it has been shown (Day
and Thorpe 1996) that the corresponding spectral functions are connected by the

following relation
th.(t,f) = (1 — t)h(1 — t,1 — f) + I'V(f)t6(t), (3.107)

where the weight of the delta function W(f) is

1h(.r,1— f)d.r

1—3:

 

W(f)-_-1—/0 (3.108)

The identity (3.107) applies generally to any uncorrelated random lattice in any
dimension, and also to infinitely interchangeable cellular materials, but it does
not apply to most continuum problems where the geometry of the host is very
different from the geometry of the inclusion.

From the other side, forthe electrically isotropic lattice the reciprocity

relation (Straley 1977) is
0(00,0‘1,L)0'(01,0‘0,D) = 0001, (3.109)

where L denotes the original lattice and D its dual lattice. The dual lattice is
constructed by placing a bond with conductance 00, 01 perpendicular to every
bond of the lattice L with conductance 0'1, 00. The 2d square net with random
bond substitution (dual lattice is not relevant for the site substitution when there
are three types of bonds in the system), is electrically isotropic and self dual and
(3.109) becomes

0(0‘0, 0'1)o(ol, do) = 0001. (3.110)

Combining (3.106) and (3.110) it follows

0(007017f)0(0130011_f):00017 (3-111)

53

or in terms of the variable t

a(t.f)o-(t.1-f) :2. (3.112)

00 00 t

 

At percolation when f = f6 = 1/2, (3.111) and (3.112) reduce to the

following respective equations

0(00,0'1) = #0001 (3.113)

and

o(t,1/2)_ t—l
00 — t . (3.114)

From this one can extract the exact result for the spectral function at percolation

h(t,1/2)=l 1. \ (3.115)

77 t

3.8.1 Dilute Limit
It has been shown (Thorpe and Tang 1987) that in the low concentration limit,
f < 1, the effective conductivity can be written as

(01—00)
f1(§:—1)+1’

where f, is determined from the initial slope of the effective conductivity when

 

”(00,011f)=00+f (3-116)

01 = 0 and depends on the lattice. For 2d square net bond substitution f1 = %,
while for site substitution f1 = 1 -— i. The equation (3.116) holds for both bond
and site substitution, if the site substitution is defined as described above. Note
that other definitions of the site substitution are possible (Thorpe and Tang 1987;
Thorpe, Djordjevié and Hetherington 1994) where the site defect introduces all
bonds adjacent to it to have a conductivity 01, if the host bonds have a conductiv-

ity 00. In this case equation (3.116) holds only for bond substitution, and one has

54

to write another equation for the site substitution [see equation (6.3) in Thorpe,
Djordjevic' and Hetherington 1994].

Using the t variable equation (3.116) can be rewritten as

 

0(t.f) _ f . ..
00 —1 t—fI’ (3.111)

from which using (3.105) the spectral function is obtained
h(t,f) = f6(t — f1). (3.118)

In the high concentration limit, 1 — f < 1, we can invert the roles of the

host and the inclusions, and write similarly

(00—01)

 

 

 

1, 0,1- = 1 1— -

0(0 0 f) 0 +( f)f1(§‘11—1)+1 (3119)
and the effective conductivity normalized again by 00 is given by

00.7) _ _ 710 -f) __(f-f1)

a. ‘1 (lemma—fl) (I—mt’ (3m)
which gives for the spectral function
_f1(1"f) _ _ (f-fI) ‘
h(t,f) — (1-f1) 6(1 t f1) + (1 _f1)6(t). (3.121)

The equation (3.121) is an example of the general exchange relation between the
spectral functions (3.107).

The poles of the normalized effective conductivity correspond to the reso-
nances of the network (Day and Thorpe 1996). This can be physically understood
by considering a network of pure capacitors capacitors 00 = iwC and inductors
01 = 1/(iwL). Then q = 01/00 is real and negative, and the effective conduc-
tivity will diverge at the resonance frequencies of the system. If either 00 or 01

has any resistive component (a real part), q will have an imaginary component

5 5

and the system will be damped and the effective conductivity will not diverge.
So the singularities are confined to q real and negative. In this case the variable
t = 0'0/(0'0 — 01) is real too. Nevertheless, to obtain the spectral function h.(t,f)
one has to calculate the effective conductivity for complex t, by adding a small

imaginary part c to it (Day and Thorpe 1996).

3.8.2 Effective Medium Theory

The effective medium theory for the conductivity can be constructed from the
results for the dilute limit in the usual way (Briggeman 1935; Day and Thorpe
1996). Substituting the effective conductivity 0 instead of 00 in (3.116) and sum-
ming Over the all components in the system, we get the self—consistent effective

medium theory equation

 

fn(0n"”0) __
27.:(1-'f1)<7'l'f10n _0’ (3'12?)

from which, in the case of only two components in the system, 00 and 01, we get

(01-0)
1+f1(§,’*—1)

 

 

(00—0)
1- =0. 3.123
+( f)1+f1(£afi—l) ( )

This leads to a quadratic equation for the normalized effective conductivity 0/0'0
from which one can extract the general form of the spectral function (see Day and
Thorpe 1996 for details).

For random bond substitution and the square net f1 = % the spectral

function is given by (Day and Thorpe 1996)

 

 

h(t,f) = ““113:in ,for 0 s f(1— f) — (3 —1/2)2

= 0 , otherwise (3.124)

56

In the effective medium theory the percolation threshold is fc = 1 — f1, but for
random bond substitution this coincides with the exact percolation threshold,

fC = i" The equation (3.124) becomes

1 l—t
—— (3.125)

Mt» fc) = ; t
which is the same expression as the exact result (3.114) obtained before. This
suggests that effective medium theory might be a useful guide at higher con-
centrations for deducing the overall form of the spectral function. An example
of the system which obeys formulas (3.113), (3.114) and (3.115), is an infinite
chess board with an equal number of black and white squares. We note, that
all formulas given here in the dilute limit and in the effective medium theory
are equivalent to the corresponding formulas for bond substitution in Thorpe,
Djordjevic' and Hetherington (1994), which can be easily seen if we introduce the
variable 7 == (01 — 00)/(01 + 00) instead the variable t = 0'0/(0'0 — 01), and then

express all results in term of the new variable a: defined as .7: = —%. Thus, for

example, one can rewrite formula (3.124) in the following form

129.7): M .for22£f(1—f)

«(Ll-2.1:)
= 0 , otherwise. (3.126)

and, at percolation, equation (3.125) can be expressed as

(3.127)

 

In Figure 3.11, we plot the spectral function (3.126) for several concentrations f,
including the the spectral function at percolation concentration f = 1/ 2 (3.127)
(see Thorpe, Djordjevic' and Hetherington 1994).

Day and Thorpe (1996) have performed computer simulations for both site

and bond substitution on a square lattice, and compared the obtained spectral

01
*1

 

2.5

 

 

 

 

b I 1: I l I I I I I I U I I I I I I r l I I r I l I I I I -'
: 1 Effective medium theory I
- I .,
20 _ l —-1
A b | d
3, - , {-0.50 «
n: " I ..
,5 1.5 — 1 -—
6 - . 7
I: ' j q
a - 1 r-o.zo .
'3 1— ! ._:
‘3 1.0 P 1 3
3. I ' I
m . ' {-0.05 4
0.5 — ' -
b I ’ 'l
1- | \ q
i ' N 3

0.0 1 l 1 1 1 1 1 1 1 1 1 1 I 1 1 1 1 . h.

-0.0 -0.4 -0.2 0.0 0.2 0.4 0.0

Figure 3.11: Spectral function h(x, f) for various concentrations f from the effec:
tive medium theory, as given by equations (3.126) and (3.127) (Thorpe, DJOI‘ClJeVlC

and Het herington 1994).

58

function with the predictions of effective medium theory. For illustration, we shall
give here just results for the bond substitution shown in Figure 3.12. For more
details about the simulation technique and their analytical calculation of the first

five moments of the spectral function, the reader is referred to the work of Day

and Thorpe (1996).

3.9 Conclusion

All results presented in this Chapter are valid only in the low concentration limit
f < 1, when the inclusions can be regarded as well separated and thus non—
interacting. In other words, the effective conductivity and the spectral function is
correctly obtained to the first order in concentration f. The only exception to this
is the computer simulation results of Day and Thorpe (1996), and the results of
the effective medium theory, which are valid in the entire range of concentration
f. We have shown that exact solution exists only in a few cases, namely for
elliptical (circular) inclusions. The conformal mapping technique turned out to
be a very powerful tool for calculating the effective conductivity for particular
type of inclusions (holes or perfectly conducting inclusions) which may have sharp
corners. The problem becomes much more difficult for inclusions with arbitrary 7
if they contain sharp corners. In this case one has to numerically solve a Fredholm
equation for the induced charge density around the perimeter of the inclusion. It
appears that the only direct numerical calculation of the spectral function [from
the general definition (2.19) or (3.105)] has been done in computer simulations
of Day and Thorpe (1996) on the discrete random square net. This has been
possible because of the use of a powerful Y — A transformation technique (Frank

and Lobb 1988) which calculates the effective conductivity of a random network

59

 

 

0.8 ,_ u L
* 9'0-93 b 9-0.50
* L
' F
a '1 - 5.5 )_
s 1 »
).
0.1 ,_ u L
p
i J
.3 ‘ 1 7 l 1 0.. . 1 1 m

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 3.12: The spectral function h(t, f) for a two-dimensional square net with
random bond disorder at different inclusion concentrations f. The solid lines are
data from the simulations. The dashed lines are the effective medium theory
given by (3.124). For f > 0.5 the dotted lines are the spectral function at (1 — f)
transformed by (3.107) (Day and Thorpe 1996).

60

very efficiently. Finally, the effective medium theory seems to give very good
results for intermediate concentrations around the percolation treshold fc, and
thus proves itself useful for deducing the overall shape of the spectral function.
In the next Chapter we shall present our results for another case when one
can go beyond the linear approximation in concentration f to include f2 terms,
or in other words, the pair interaction between the inclusions. So far, this was
possible to do only for circular inclusions (Djordjevié et al. 1996), and for random

resistor network (Day and Thorpe 1996).

Chapter 4

Pair Terms for Circular

Inclusions

4. 1 Introduction

The conductivity (or dielectric behavior) of a binary composite system is con-
veniently expressed in terms of a spectral function, which is determined by the
geometry of the composite. In this Chapter we examine the case of circular in-
clusions in a conducting sheet and keep terms up to second order in the inclusion
concentration f. The two inclusion problem can be solved exactly using multiple
images (Djordjevié et al. 1996), and we use this solution to construct the spec- .
tral function. We show that the spectral function is a truncated Lorentzian, that
can be calculated in a simple closed form. Both the weight and the width of the
spectral function are linear in f.

Exact results for the properties of composite materials are rare, and so the
few that exist are particularly important as aids to our understanding. Although

spectral functions are widely used in physics, their use in composite materials has

61

62

not been widespread. This is mainly because this attractive formalism is difficult
to apply to particular composite geometries, because of the numerical difficulty of
computing the spectral function. In this paper we examine a non-trivial example
of a spectral function that can be calculated exactly. This example involves the
calculation of the spectral function of a sheet containing circular inclusions with
a different conductivity from that of the host. Using the well-known technique of
multiple dipole images, the terms in the effective conductivity up to second order
in the concentration of inclusions f can be found. We have been able to express
the spectral function in a rapidly convergent series. We show that the spectral
function is a truncated Lorentzian with weight f and a width that is pr0portional
to f, and we are able to give the final result in a very simple form.

Here, we will apply the general formalism of Chapter 2 to a two dimen—
sional sheet containing circular inclusions, and obtain the spectral function to
second order in the area fraction f of inclusions. In other words, we shall add one
more term (f2) in the effective conductivity (2.24) and then calculate its imagi-
nary part to obtain the spectral function. The layout of this Chapter is as follows.
In the next section we restrict the form of the conductivity 0 by using the reci-
procity theorem (2.23) that is valid in 2d, generalizing the previous result (2.24)
from Chapter 2. In section 4.3, we use the method of multiple images to write
down the position and strengths of the dipoles in a pair of circular inclusions. In
section 4.4, we examine the results in more detail, particularly for the case of per-
fectly conducting inclusions. In section 4.5, we obtain the effective conductivity
by averaging over the dipole moments of the pairs that we have found previously.
This solution is then tortured into a form from which the spectral function can be

derived in sections 4.6 and 4.7. Finally in the conclusions, we discuss the uses

63

of this kind of approach.

4.2 General Form for the Two—dimensional Con-

ductivity to the Second Order in f.

The problem of circular and polygonal holes in two-dimensional conducting sheet
was recently investigated (Thorpe 1992) up to the first order in the area fraction
of inclusions, using the conformal mapping technique. From the other side, the
general problem for circular inclusions had been solved earlier up to the second
order (Peterson & Hermans 1968) using bipolar coordinates.

We use the method of electrical images to solve the problem of the ef-
fective conductivity in the conducting sheet containing circular inclusions of an
arbitrary conductivity. This method is not limited to holes and perfectly conduct-
ing inclusions as is the conformal mapping technique.

We first consider one isolated pair of identical circles with radius a sep-
arated by a distance R. Considering the two independent cases of parallel and
perpendicular external electric field, we construct the infinite set of dipole, or
dublet images (Binns & Lawrenson 1962). Each dipole image is found in a simple
continued-fraction-form. We find and solve the appropriate difference equation for
successive dipole images to get all dipole moment images, now in closed form, in
terms of hyperbolic functions. All further calculations are based on these dipole
images. Naturally, the total dipole moment of the two circles is given by the sum
of all image-dipole moments in the case of parallel external field, and by the same

sum, but with alternating signs in the case of perpendicular external field.

The reciprocity theorem (Keller 1964; Mendelson 1974) (2.23) for the two

64

phase medium in 2d, states that following relation must hold (Hetherington 85

Thorpe 1992) [we are rewriting here the equation (2.23) from Chapter 2].
0(ao,ol)o’(ol, 00) = 0001, (4.1)

where 0'0 and 01 are conductivities of two phases, and 0(00,01) and 0701.00)
are the effective conductivities of the original system and the one with the same
geometry but 00 and 01 interchanged.

we find the generalization of the effective conductivity (2.24) (See Hether—
ington & Thorpe 1992) to the second order in the concentration f of the inclusion
phase 1 in the host phase 0. The effective conductivity then must be of the

following form

 

010 = 1 + 776(7) + 72%) + 003). (4.2)
where
_ 0'1 — 0’0
7_ 01+00, (4.3)

and the function C(7) is an even function of 7 (Hetherington & Thorpe 1992) as
follows from the reciprocity theorem (4.1). Similarly, for interchanged phases 0

and 1, (7 2 —7) the effective conductivity 0’ is given as

= 1-f7G(7) +f2H(-7) +0(f3)- (4-4)

1

Multiplying (4.2) by (4.4) and using the reciprocity theorem (4.1), we get to the

second order in concentration f
H(7) + H(-7) - 7202(7) = 0. (4.5)
from which it is clear that function H (7) must generally be of the form

H(7) = grazuwnm. (4.6)

65

where the function L(7) is an even function of 7.

The weak scattering limit (see also Hetherington & Thorpe 1992) gives
L(7) = 7F(7), where 17(7) is an odd function of 7 (proportional to 7 for small
7), and F(’0) is zero. So, finally (4.2) becomes

£— : 1+f7G(~/)+$272037)+f272F(7)+0(f3), (4-7)

0

where C(7) and F(7) are an even and odd function of 7, respectively. The area

fraction f can be written as
f = 717702, (4.8)

where n is the number of inclusions per unit area, each with area 7ra2. Two usual
limits are obtained when the inclusion is a hole (7 = —1) and when the inclusion
is perfectly conducting (7 = 1). Equation (4.7) is valid generally but for circular

inclusions C(7) = 2 as derived before (3.17) (See Thorpe 1992). So, for circles as

inclusions (4.7) simplifies to

7:; = 1+ 27 f + 2W" + 772777) + 0(f3). (4-9)

and our goal is to calculate the function F (7) which contains the pair terms in
the effective conductivity.

We relate the function F (7) to the pair contribution to the polarizability
AB. This is defined as the mean change in the polarizability of a single inclusion,
due to the presence of the other inclusions. To leading order, AB can be obtained
from the dipole moment of pairs. Indeed we have stressed in previous papers
(Thorpe 1992, Hetherington & Thorpe 1992) that the conductivity of a sample
is directly related to the dipole moment of all the inclusions. Once we have

all dipole moment images in the aligned and perpendicular external field, it is

66

straightforward to find the polarizability change AH by averaging the parallel and
perpendicular two-circle total dipole moment and integrating over the whole sheet.

Then we use relation

72f2F(7)=nA3. (4.10)

4.3 Dipole Moments

We shall first consider the case of an aligned (parallel) external electric field with
respect to the line connecting the centers of the two circles, as shown in Figure 4.1.
Circles are not charged and have an arbitrary conductivity 01. In the external field
they become polarized, i.e. there is an induced dipole moment on each circle. So,
it is intuitively clear that we have to deal with dipoles in constructing the infinite
set of images which solves this electrostatical problem.

The images of a dipole can be derived from the combination of the separate
images of the pair of charges forming the dipole (see Binns & Lawrenson 1962,
pp.48-53). By definition, a two dimensional dipole is a combination of a positive
line charge and a negative line charge, an infinitesimal distance apart, so arranged
that the product of charge and distance remains finite.

One circle in a uniform external field is completely described by one dipole
at the center of the circle. We start from that dipole at the center of one circle 7

which we denote by p1, given by

p1 = 2777an2. (4.11)

Assuming that we started from the circle to the right, and taking in

account the presence of the other circle to the left, we find the dipole image p2 of

67

 

1 .91
m

 

 

 

 

 

 

 

 

 

Figure 4.1: Orientation of the dipole images for two circles in a parallel applied
field E0, showing the notation in the text (Djordjevié et al. 1996).

68

P1 in the circle to the left. Now, we find the next dipole image p3 of p2 in the circle
to the right. Continuing in that way we construct an infinite set of dipole images.
Odd numbered images are all inside the right circle, and all even numbered images
are inside of the left circle. But, due to the symmetry between the two circles,
there must be another infinite set of images generated by placing the first dipole
p1 at the center of the left circle. The complete solution of the problem consists of
two infinite sets of images, redistributed in such a way that each circle now has the
same infinite set of images p1,p2,p3..., which are given in the simple, continued

fraction form
a 2
p2 = 7 (E) p17

2
a
P3 = 7( )P2, (4.12)

 

 

P4=7

where 7 is given by (4.3), a is the radius of each circle, R is the separation between
their centers, and p1 is given by (4.11). All images of dipole moments have the
same direction, along the applied field, as shown in Figure 4.1.

The case of the perpendicular external field is shown in Figure 4.2. Due to
the direction of the applied field, all dipole moments are now perpendicular to the
horizontal connecting the centers of the two circles. The first dipole image p2 now
has the opposite direction with respect to the initial dipole p1. The alternating

direction of successive images is the only new thing here, so we can repeat essen-

 

69

tially the same procedure we used for the parallel field case, to get exactly the
same infinite set of dipole images whose magnitudes are given by (4.11) - (4.12),
as before. Alternating orientation of these dipoles are represented symbolically in
Figure 4.2.

Positions of the first few images in both circles, measured from their re—
spective centers, as shown in Figures 4.1 and 4.2, are found in the similar continued

fraction form.

 

$1 = 09
a2
3:2 = E1
(4.13)
02
$3 = 2 9
R — 7.-
a2
$4 = 2 a
R - —2—.—
R-a-

These positions hold for both the parallel and perpendicular case.
Combining (4.11)-(4.13), it is easy to see that the dipole moments satisfy

the following difierence equation

1 1 R 1

 

 

 

 

+ = — . (4.14)
VPH-l 7vps-1 a 7R!
Solution of (4.14) gives all dipole images in terms of hyperbolic functions,
sinha 2
Pa = 277130027, ( . ) , (4.15)
Slnh sa

\

 

 

 

 

 

 

 

 

 

Figure 4.2: Orientation of the dipole images for two circles in a perpendicular
applied field E0, showing the notation in the text (Djordjevic' et al. 1996).

71

where a is given by the relation
R
cosha = Z. (4.16)

The equation (4.15) is our main result for dipole moment images, and all further

calculations are based on it.

4.4 Perfectly Conducting Inclusions

As an interesting sideline, we look at some properties of a single pair of inclusions.
The formulae for this case have a particularly simple form. We calculate the total
dipole moment of the two perfectly conducting (7 = 1) circular inclusions in an
external field, separated by distance R. It is conceptually obvious that in the
parallel field, the total dipole moment is given by the sum of all dipole images

multiplied by factor 2 (two circles)
P“ = 22:“ (4.17)
3:1

Using (4.15) we have for perfectly conducting circles

 

0° . 2
I)” :: 477E002 z ( Slnha ) , (4.18)
5:1

sinh'sa
which is the total dipole moment of two circles in the parallel external field.
In a perpendicular external field, the only difference is the alternating
orientation of the consecutive dipole moments, while they still have the same

magnitudes as in the parallel field case. So, we can write

00 - 2
Pi = 477an22 (-1)3+1 (El-Pi) . (4.19)

,=1 sinh 307
Using (4.18) and (4.19), we numerically calculate the total dipole moment

dependence on the separation R/ a of the two circular inclusions, in both cases.

:1
[0

These graphs are shown in Figure 4.3. When the circles touch, (4.18) and (4.19)

give exact values for the two dipole moments (with a ——+ 0).

°° 1 2
P” : 477an2 Z 3—2 : 477E0a2C(2) : 5E0a2773, (4.20)

3:1

1 1
'3 = 277E002C(2) = —E0Cl27l'3, (4.21)

Pi = 477an2 Z (-1)’+l
3:, s 3

where the sums are expressed in terms of the Riemann zeta function (Gradshteyn
& Ryzhik 1965).

In this particular case of perfectly conducting circular inclusions, and in
parallel external field, each consecutive dipole image introduces some change of

the potential on each circle. The first few contributions in the decreasing potential

of the circle 01, which correspond to the first few dipole images p1,p2, p3, ..., are
V1 = 0.
1 P1
V = ———,
2 277 R
(4.22)
1 P2
V : —— ,
3 277 R — %
1
W _ "—é— p3 02 7
77 R — —;
3"“?

where p1,p2,p3, are given by (4.11) and (4.12). In general

= W321, (4.23)

8
277a2

 

'7‘" UrUUIU'UUlITUIIIIVU'UIr

 

mmamm

 

'1
fl

 

r l L l L LLL L L I l 1 l L l l l l L l l l L L l l l l L
0 2 4 0 8 10 12
' R/e

 

 

Figure 4.3: Total dipole moment, P” and P; of two perfectly conducting circles
in parallel and perpendicular external field. The two diamond symbols when the
circles touch at R/a = 2, are from equations (4.20) and (4.21) (Djordjevic’ et al.

1996)

74

In the absence of inclusions, the external electric field E0 creates a potential

difference between the two centers given by
AVE = EOR. (4.24)

The total potential difference, i.e. voltage AV, between the two perfectly con-

ducting circles in the external field is given by

AV = AVE—22ml. (4.2.5)
3:2

Combining (4.12), (4.13) and (4.22) the total voltage can be expressed in

terms of a sum over hyperbolic functions

 

(:2 (sinhoz)2 )
AV = an —-—2Z . (4.26)

3—2 sinh sa sinh (s —- 1)a
The series in (4.26) is summable to a simple form. To show this rewrite this sum
with the help of (4.16) and in terms of the new variable y = e". (4.26) then

becomes

1 °° 1
3” .IJ '2 (3" 9) §( (y _y-.)(y._. _y-(.-.))]- (4'27)

Rewriting in terms of partial fractions (4.27) becomes

1 1 °° 1 1

y 3:2

 

AV = E00

 

 

Writing explicitly the first few terms in the sum it is easy to see that only .
one term in the second partial fraction survives for s = 2, while all others cancel

each other to give
AV 2 an (y — —) (4.29)

and hence,

AV = Eo\/R2-4a2. (4.30)

75

This remarkably simple result for the voltage between two perfectly con-
ducting circular inclusions is shown graphically in Figure 4.4. Presumably (4.30)
can be derived by a more elementary method, but we have been unsuccessful in

doing 50.

4.5 Effective Conductivity

As we mentioned in the Introduction, the problem is to find the function F(7) in
(4.9). The function F(7) contains the information about the long range Coulomb
interaction between all pairs of circles. Obviously, we have to connect the function
F ('7) with our solution for dipole moment images (4.15). But the dipole moment
p is related to polarizability by the relation p = 6E0, where B is the polarizability.
So, we first relate 17(7) to the pair correction in the polarizability.

Introducing an effective conductivity 0 in the observed region containing
11 circles per unit area, and the mean polarizability of a single inclusion 23., we find

the Clausius-Mossotti type result [see the general equation (2.20)]

(4.31 )

 

Which defines 3 (see the general equation (2.20)) The mean polarizability 3 con-
tains higher order effects, so we solve (4.31) for 0/00, and expand it to second

order in n

0' _ 712—2
—=1+nfl+ 3

+ . . .. (4.32)
0’0 2

 

Separating the single site and the correction due to pairs in 23—, we put

8 = 27ra27 + M, (4.33)

 

Voltage V

 

LlllllllLLLllllllLlllAlllll

 

o I l 1 l LLL l l L 1 LL] L1 LA I l l I I l l l l L
0 2 4 O O i 10
Whammtm IV.

on.
(9

F igure 4.4: Voltage between the two perfectly conducting circles in the external
fie d is shown as a solid line. The dashed line would be the voltage in the absence
of the conducting circles, (i.e., just host material) (Djordjevié et al. 1996).

77

where A43 is the change in the polarizability due to the presence of other inclusions

nearby. Therefore we find (4.9) in the form

(7
E_=1+27f+272f2+nA6 (4.34)
0

which is consistent with (4.10).

The change to the average dipole moment can be written as
_ 21r 00
P’ = n f / (Pl; cos2 9 + P1 sin2 6) d6 RdR, (4.3.5)
0 2a

where Pll and P1 [and hence 13’] don’t contain p1 [as designated by the primed
quantities], but are summed over all the other induced dipoles as in (4.18). Terms
containing p1 are already taken in account, in the single site part of polarizability
(terms proportional to f in (4.34)).

The series (41) is only conditionally convergent as the leading term invov-
ing p2 is divergent unless the angular integration is done first. Integration over 0
gives factor 1r, and even dipole images cancel out in the integrand, leaving only odd
terms doubled. The elimination of the even terms is necessary in order for F (7)
to be an odd function of 7 as required by the reciprocity theorem and discussed
in section II. It may also be possible to show that the even terms vanish in (41)
using methods similar to Felderhof, Ford and Cohen (1982), but this is a lengthy
procedure that we have not attempted as we know that the even terms in (41)
must integrate to zero using the reciprocity theorem. In three dimensions, the
reciprocity theorem cannot be invoked and hence it is necessary to do the analysis
surrounding the analogous equation (3.35) in Felderhof, Ford and Cohen (1982).

Recalling the relation between the polarizability and the dipole moment,

we can write

2f cc

13002 20

 

P
Aflz-E-O-z (p3+p5+p7+...)RdR. (4.36)

78

The behavior at large B in the integrand in (4.35) is dominated by the n. =
3 term in the summations in (4.18) and (4.19) which leads to the assymptotic form
of the integrand in (4.35) going as R“, which leads to an absolutely convergent
integral in (4.35) and (4.36). With the help of (4.10) the function 17(7) is written

as

 

 

A, _ ‘2 ”(p3+p5+p7+...) ._
F(,)_ «Em/1.. 7, RdR. (4.31)

Using the hyperbolic function form of dipole images (4.15), we have the following

expression

 

0° °° sinh3 a cosh 01
= 16 f 28-1 d . 4.38
0 Z ’7 sinhz (2s + 1)oz oz ( )

Finally, (4.38) can be written in the form
F(7) = 16 (Kn + 14273 + 14375 + K477 + . . .) , (4.39)

where K, is given by

 

K: /0°° sinh3 a cosh a

4.4
sinh2( 2s +1)a ( 0)

Substituting (4.39) into (4.9) we have found the expression for the effective
conductivity to the second order in the area fraction of the inclusions, f, for any
value of 7 between 0 and 1:

a

a— = 1+ 27f + [2‘12 + (1614373 + 161(275 +16K377 + - - )] f2, (4-41)
0

where K, are given by (4.40). It is useful to simplify the coefficients K... First,

we transform the integral (4.40) using the identity

1 _ °° 4j
sinh2 [(23 + 1)a] _ 2 exp [2ja(2s + 1)]

J:

(4.42)

 

 

Table 4.1: Taylor series for circular inclusions of conductivity ratio parameter 7.

Taylor series

 

 

Circle F(7) = 0.66667 7 + 0.05567 73 + 0.01317 75 + 0.00464 77 + 0.00204 79
+ 0.00104 7” + 0.00058 713 + 0.00035 715 + 0.00022 7” + 0.00015 7”
+ 0.00010 721 + 0.00007 723 + 0.00005 725 + 0.00004 727 + 0.00003 729
+ 0.00002 731 + 0.00002 733 + 0.00002 735 + 0.00001 737 +

 

 

which can be easily checked (Gradshteyn Sc Ryzhik 1965). Substituting (4.42)

into (4.40) we have the following expression for the coefficients K,
K, = 4 Zj/ sinh3acoshaexp [—2joz(2s +1)] da (4.43)
i=1

which was previously found by Peterson and Hermans (1968), using bipolar coor-
dinates.

Substituting (.0 = exp[—2ja(2s +1)] in (4.43), carrying out the integration
over a, and using 8.362 (1) (Gradshteyn & Ryzhik 1965), we express the result in

terms of the digamma (\II) function

K, = mahogwwtx.)Mia-30:31)]

(4.44)

 

 

 

 

which was also found by Peterson and Hermans (1968). The tabulated values of
the digamma function were used to construct Table 4.1.

In Figure 4.5, we show a plot of F (7), using the Taylor expansion shown in
Table 4.1, with enough terms so that convergence was obtained even at 7 = 1. This

required about 100 terms to get 4 figure accuracy. We have used an approximate

80

interpolation form to fit the result for F(7).

2 7
F = —,—————. -
(7) 3(1_ WV (4 45)
The coefficients (3 and c were chosen to reproduce the result at 7 = 1, and known

behaviour of F(7) for small 7.
2
F(7) = -3-7 for small 7, and F(1) = 0.744989676... (4.46)

Here the following a number denotes that it is known exactly to that number
of digits.

Peterson and Hermans (1968) contains a numerical error. The last two
given coefficients in their formula (67) are 0.0194 and 0.00505. That should be

corrected to 0.0198... and 0.00697..., respectively.

4.6 Analytic Continuation of the Pair Term

We shall need the both the real and imaginary parts of F (7) for the calculation
of the spectral function in the next section. The imaginary part can be obtained
by analytic continuation. Starting from (4.38) we rewrite it in slightly different

form:

723-1

= °° ' 3 d . 4.47
F(7) 16/0 smh( c)o(sh a”: sinh2[( 2.3 +1) (1] a ( )

 

The sum in the integrand

oo 2.9—1
_ 7
[4(7) = 3;; [60(2‘H'1) __ e-a(2s+1)] (448)

 

can be rearranged and expanded in the series

81

 

 

 

 

 

P
A D
F .—
17 .
P
P
r
P
0.0 L1 L L L L l A l 1 l L l l l l l l l L l 1 1
0.0 0.8 0.4 0.0 0.8 1.0
7

Figure 4.5: Interpolation formula (smooth curve) for the function F (7) (crosses).
The diamond symbol when 7 = l is from equation (4.39) (Djordjevié et al. 1996).

44(7) = $723“ [6‘2“‘2’H’+2e“40(2’+”+3e‘60(25+”+---]. (4.49)

3:1

Introducing the new index I/ E s — 1 we can transform the (4.49) into

_ e—sa .7 _7_ 2ue—12a _7_ 2" —lSa
— 7: 62—; + 670 +3 :60 e +---.

 

 

 

 

 

11:0
(4.50)
The summation over 1/ is possible for each term in the brackets to give
6—60: 26—120: 36-180:
4(7) = 7 —:—2+ 2+ 2+"' (451)
1- (:2?) 1- (72.) 1- (2%)
which can be expressed in terms of the variable :1: : —1/27.
a: 00 me—6ma
A(:1: —§ 2( -..,... 2. (4.52)
m=1(’~’—e( 2 )
We write the function F(7)‘in (4.47) 1n terms of the new v; able :1:
0° . 3 oo me-Gma
F(:c) = —322:/ Slnh (a) cosh(a) Z 2 2 da, (4.53)
0 ml ,2 _ (£223.)
or
1/2( $I)2[( 213 )1-/2m __ (2$,)1/2m]3 [(231)-1/2m + (2:13,)1/2m] I
F( (x) : -2:c WEE/02 $2 _ £2 dzc,
(4.54)

where we have introduced another new variable x’ : liexp(-2ma). Finally, trans-

forming this expression into partial fraction form, we get after some manipulations

 

F(:c) : [”2 Zm=l m(|$ x’|)da:’, (4.55)

-1/2 a: — x’

83
where (writing the dummy variable 96’ as 9:)

4.0 x I) = 4132(0er 0“" — (2 I x 0””)3 ((2 I 3 0‘1”": + (2 a .2: 01”")

(4.56)
and hence, reading from (4.55) we obtain imaginary part of F(.z:) as
1mm) -—-— 2 4.0 x l). (4.57)
m:l

for —1/2 < a: < 1/2. This derivation was done by J. Hetherington. Near :1: : 0

there is a cusp
ImF(:r) = 1—2|x|+0(|:1:|%) (4.58)

which is caused by the large number of widely separated pairs of circles. We use

ImF(a:), given in (4.57), to construct the spectral function in the next section.

4.7 Spectral Representation

It is convenient to employ here the variables used for the spectral function in
Chapter 2. Its properties were summarized in equations (2.36) — (2.38), where
the effective conductivity is expressed through the spectral function which is sym-
metric with respect to the origin :1: : 0, and the first and second moment of the
spectral function are given. For simplicity for the rest of this section, we drop the

prime on 2:. We shall rewrite formulas (2.36)—(2.38) here

a %h(:1:)d:c
— = 1 2/ ,
+ 7 —%1+27:r

 

l1 ’ h(a:)da: = f, (4.59)

84

where f is concentration (volume fraction) of the phase 01 in the host of 00.

The first moment (1:), as we have seen (2.39) is given by
(.r) : ——. (4.60)

The higher moments of the spectral function h(;1:) depend on the detailed geometry
and are not general like the first moment given by (4.60). For example the second

moment for the present case is given by
(4’2) = 2K1f, (4.61)

where 1611] is the leading term in (4.38) and from (4.44) has the numerical value

K1 : 313. Hence we have

(1:0) : if? (4.62)

The spectral function h(:1:) is formed by taking the imaginary part of
0/00 given by (4.9) and rearranged so that a self energy is in the denominator.
Such rearrangement can be accomplished formally [see for example Felderhof and
Jones (1986) and Cichocki and Felderhof (1985)]. We have been content to write
the conductivity in the form (4.63), which when expanded is equivalent to (4.34)

to second order in the concentration f:

a 1

— : . (4.63)

_ 271
0° I 14410-44171]

 

Taking the imaginary part of (4.63) (substituting 7 = —1/2x as before) we have

h(.:) = 1: lilmflfll ,, (4.64)
7' [x + g + £ReF(:c)]2 + [£1mF(fL‘)]

 

where ReF(:v) is the real part of the function F. We obtained it from the known

imaginary part in (4.56) and (4.57) using Kramers-Kronig relation. Both ImF(:c)

85

and ReF(.c) are shown in Figure 4.6. We plot the spectral function given by
(4.64), for several concentrations f in Figure 4.7, using the complete expression
(4.64). Of course, the form (4.63) is not unique to second order in concentration
f. Many other forms can be written down which when expanded will give the
conductivity 0 correct to order f2 [see (4.2)].

We note that Felderhof and Jones (1989) also use a continuued fraction
in their Eq. (3.11) identical to our (4.64). We prefer the form (4.63) as it leads
to a single smooth spectral function as shown in Figure 4.7. The single delta
function in the single inclusion limit, is broadened by the interactions between
the inclusions at higher concentrations, so that the width is proportional to f as
shown in Figure 4.7.

Similar behavior is in the spectral function associated with the vibrations
of light point defects in solids (Elliott, Krumhansl and Leath 1974). We note that
there are no isolated poles in the expression (4.63) and hence no delta functions
in the spectral function, as can be seen from Figs. (4.7) and (4.8). This can easily
be checked as our spectral functions obeyed the sum rules given in (4.59). We
interpret the spectral function in the low concentration limit as being analogous
to the broadening of single impurity lines, due to interactions between the pairs
of light mass defects in solids (Elliott, Krumhansl and Leath 1974). We note that
(Felderhof and Jones 1989) found two isolated delta functions in the spectral func-
tion in three dimensions, in addition to a continuous part. These delta functions
we believe are unphysical and result from the neglect of higher order multipoles,
even at the pair level. Because of the simplifications possible in two dimensions
we have been able to sum up all the contributions at the pair level, via the infinite

set of dipole images. As a result there are no unphysical delta functions in the

86

 

 

 

 

 

 

 

 

Figure 4.6: ImF(x) and ReF(:r) for circles, with a: : —1/(27). The real part is
obtained from the known imaginary part (4.56) and (4.57), using Kramers-Kronig
relation. Note that ImF(z) : 0 for 2:2 > 1/4 (Djordjevié et al. 1996).

 

 

 

 

l I I I I I I I I I I I I I I I I I I I I I I I I d
8.0 {-0.20 . ._
, .
1 {-0.10 .
g 1.5 : r-o.oo «
_ , . ‘1
{-0.01 .
' -(
1.0 ‘ -
i 1
0.5 L l -'
,
.. - __ JIM
-o.a -o.a -o.1 0.0 0.1 0.2

Figure 4.7: Spectral function for several values of inclusion area fraction, using

(4.64), f : 0.01;0.05;0.10;0.20 (Djordjevié et al. 1996).

88

spectral function. If we have truncated the expansion (4.39) after the first few
terms, we would have obtained a spectral function with many delta functions, that
only form a single continuous function as the series (4.39) is summed to infinity.

In the dilute limit the spectral function becomes Lorentzian. For small
concentrations f, and for small :13, using (4.46) and (4.58), we can make the ap—
proximation F(.v) : i, i.e. ReF(.r) : 0 and ImF(1:) : 1, and the spectral

function becomes

L
h(:c) : i 4 .. (4.65)

" (I + f)? + (1)2

In Figure 4.8 we show the spectral function calculated from the full ex- '

 

pression (4.64) and the Lorentzian form (4.65), for two small concentrations. It
can be seen that the differences between these two forms only become apparent at
higher values of the area fraction f, where the pair approximation breaks down.

It is easy to find the real part of 0/00 in the dilute limit:

 

12.».(3) z 1- 2f(2:c+ f) (4.66)

2 .
0° (4) + (2x+f)2
Then, assuming complex 2:, it is straightforward to show that complex

effective conductivity has the following form
(4.67)

where a: is now an arbitrary point in the complex plane. Formula (4.67) reproduces
correct imaginary and real part as given by (4.65) and (4.66), respectively. Of
course, the form (4.67) is not correct as the single pole is off the real axis, which is
not allowed. Also the Lorentzian extends out to :too, which violates the fact that

the spectral function is only non-zero between 21:]. Nevertheless, the Lorentzian

89

 

IIIIIIIIT'fIIIlIIII'IIIIII’I

 

 

 

 

 

1.0 '-
L , ‘
q '1 ‘ r-om -
3’ - l ‘ . 4
‘00 _ I‘O.“ ’ ‘ -'"
. . I 1
1 ‘ i
D ‘ q
n I \
I ‘ ‘
P ’ '1
0.0 r- ’I . . _
b ’ q
I
l I a
h I, .1
l’ 4
0.0 L L l A l l I A A
-0. 1” -0.070 4.50 4.080 0.000 0.025

Figure 4.8: Lorentzian shape of spectral function h(:z:) from (4.65) compared with
the exact result (4.64). The Lorentzian is shown as the dashed line, and the exact
result as the solid line, for two concentrations f : 0.05 and f = 0.10 (Djordjevic’

et al. 1996).

90

form, and hence (4.67), form a good approximation for small f and 3:, as shown
in Figure 4.8.

We also define the function

 

 

0' -— 0'0
M E .
(7) 0+00 (Jr-68)
It can be shown that imaginary part of 111(r) is given by
tImF ~
Im1M(.r) = 4 4 2 (x) 2 (4.69)
“ (a: + £R€F($)] + [£1mF(1‘)]
which in the dilute limit becomes
[mill/[(1) : é—i—f, (4.70)
£2 + (i)

which is both Lorentzian and an even function of :r.

In Figure 4.9 we contrast plots of the full expression (4.69) and Lorentzian
form (4.70) of function ImM(:1:), for two concentrations. Comparing Figure 4.9
with Figure 4.8 we see that function I mM (2:) given by (4.69) keeps its approximate
Lorentzian form up to relatively higher concentrations than does the spectral

function h(x).

4.8 Conclusions

We have shown that the spectral function for a sheet containing holes can be
calculated correct to second order in the area fraction of inclusions f. The result
is a truncated Lorentzian that for all practical purposes is a Lorentzian in the
range of interest. A knowledge of the spectral function permits the (complex)
conductivity, 0, for any (complex) values of 0’0 and 01 to be found, by doing a

single integral. This is particularly useful if there is a parameter in the composite

91

 

 

 

 

 

 

 

I I I I I I I I I I I I I
2.0 f- l ‘ l r
. ' 1 q
1 , <
D ' d
1
1.0 1- 1 1 ‘ _]
.- l .7
. 1 “ ,
1 1
A b
11 1 1 1 ‘ ‘
V P I , | ‘ .‘
3 to l- , 1' 1 ‘ -0.00 —
. 1
1 1
h i 1 1 “ ..
1' 1 l 1 -o.1o ~
- 1 ' 1 .
0.6 1— ’I _.
I I ‘
oo («4...— L . L 1 1 A L 5*: -
-0.1 0.0 0.1

Figure 4.9: The function I mM (:13) in the dilute limit. The solid curve is the exact
result from (4.69), and the dashed curve is the Lorentzian approximate (4.70).
Results are shown for two concentrations f : 0.05 and f : 0.10 (Djordjevié et

aL 1996).

92

that can be varied at fixed geometry. Examples would be the frequency and the
temperature.

The calculation in 2d was possible because an exact solution to this prob-
lem is available in terms of multiple images. We have managed to do an analytic
continuation which has allowed the spectral function to be found which is finite
only over a narrow region at low concentration.

Similar spectral functions can also be used in discrete systems made up
of resistors (Foster 1924 and Straley 1979). Spectral functions are widely used
in physics to study electronic and vibrational excitations. They have proved es-
pecially useful in systems containing impurities, where local modes can occur
(Elliott, Krumhansl and Leath 1974). The weight in the absorption spectrum
goes like the defect concentration f and the width is determined by the interac-
tion between defects and goes also as f. The reason for this close analogy between

these two somewhat disparate situations is not entirely clear at the present time.

Chapter 5

Discrete Disordered Media

5.1 Introduction

In the preceding chapters we have studied the properties of the continuous dis-
ordered media. The only exception to that was the case of the discrete random
resistor network, which we presented at the end of Chapter 3. We shall now
turn our attention to another wide class of discrete disordered systems. These are

topologically disordered networks, or amorphous solids.

5.2 Topological Disorder

A perfect crystal is the structure in which the atoms, or groups of atoms, are
arranged in a. pattern that repeats itself periodically in the three dimensions to
an infinite extension. For instance, crystalline diamond makes a lattice which is
produced by periodical repetition of the 8—atom diamond cubic cell in all three
spatial direction. A piece of the crystalline diamond structure is shown in Fig-

ure 5.1. Amorphous materials do not have the periodicity (long range order)

93

94

 

 

 

 

 

 

Figure 5.1: A piece of a crystalline diamond structure. One of the characteristic
sixfold rings .in a crystalline diamond structure is shown by the highlighted atoms.

95

 

 

 

 

 

 

Figure 5.2: Schematic sketches of the atomic arrangements in (a) a crystalline
solid, (b) an amorphous solid, and (c) a gas. (Taken from: Zallen 1983.)

characteristic of a crystal. Starting from a crystalline structure and slightly dis-
placing every atom in a random manner, the periodicity is certainly destroyed, but
the random structure obtained is not yet amorphous. An amorphous structure
is topologically different from a crystalline one. Thus, to obtain an amorphous
structure from crystalline one, it is necessary not only to introduce randomness
in the atomic positions, but also to change the topology of the original perfect
lattice. We can schematically sketch atomic arrangements in three typical cases:
a crystal, an amorphous solid, and a gas, as in Figure 5.2.

In Figure 5.3 we show a piece of the 3d amorphous diamond structure
which was obtained in our computer _simulations, and which we shall talk about
in the next Chapter. The two networks in Figures 5.1 and 5.3 are topologi-
cally different: while the crystalline diamond network possess only 6-fold rings of
atoms, the amorphous network also has 5-fold, T-fold, and higher order rings, and

consequently the number of 6-fold rings is reduced from its crystalline value.

96

 

 

       

  

{15
J flins‘giw

   

Figure 5.3: A piece of an amorphous diamond structure (Djordjevié et al. 1995).

97

 

(I) (b)

Figure 5.4: (a) Hypothetical crystalline compound A285. ((1) Amorphous form of
the same compound. (Taken from: Elliott 1984.)

Similarly, binary compounds like SiOg, can be found both in the crys-
talline and in the amorphous state, the later being topologically different from-
the crystalline Si 02. Schematic two—dimensional representations of the crystalline
and amorphous structures for a hypothetical A283 compound are shown in F ig-
ure 5.4. In the rest of this dissertation we shall be mainly interested in elemental
amorphous structures, even though we shall briefly mention how one can make a

computer model of an amorphous binary compound.

5.3 Distribution Functions

The experimental evidence of a fundamental difference between the crystalline
and amorphous structures is given by the diffraction patterns obtained when an
incident monoenergetic collimated beam of either X—rays, neutrons or electrons,
is scattered from the atoms of the structure. In the case of a crystal, a beam of,

say. X—rays is elastically scattered from regular Bragg planes, according to the

98

 

Figure 5.5: Typical diffraction patterns for a crystal (a), and for an amorphous
material (11). (Taken from Elliott 1984.)

Bragg condition for the constructive interference of the outgoing beams
2dsin0 :nA, (5.1)

where d is the spacing of the two Bragg’s planes, 20 is the scattering angle and A
is the wavelength of the incident X —ray beam. As a consequence, a characteristic
diffraction pattern with a series of a sharp spots is obtained [Figure 5.5(a)]. In-
amorphous materials such sharp structure in the scattering pattern is absent, and
all one sees is a concentric sequence of diffuse rings, as shown in Figure 5.5(b).
A collection of a large number of small randomly oriented perfect crystals would
produce a sequence of concentric sharp rings. Thus, unlike in the case of crystals,
the variation of the intensity in the diffraction pattern of amorphous materials
cannot be used to determine a well—defined structure. The most one can determine
is the distribution in the separation of pairs of atoms, usually expressed in the
form of a distribution function.

There are at least four different (related) distribution functions that are

99

used in the literature, and there is a certain amount of confusion in how different
authors name them. We shall introduce here all four of them and explain their
mutual relation. Ditribution functions serve as the most common way of describing
amorphous structures. They can be calculated from the theoretical models, or
deduced from scattering experiments.

We shall start with defining Pair Distribution Function (PDF). If the
average number density of the structure is no : N/V, where N is the total

number of atoms, and V is the volume of the sample, then we can generally write

1

PM = n N1}: 60 —1R.~— Ran, 15.2)
0 #1

 

where bold—faced letters denote vectors, angular brackets represent the ensemble
or configuration average, and the vectors R.- represent the positions of all the
atoms in the structure. The sum in (5.2) picks up a delta function whenever any
two atoms anywhere in the sample are separated by R,- — Rj. This step is done
for each possible configuration, and a weighted average of the results is taken in.
the usual thermodynamic way. The result is made dimensionless by dividing by
N no. In the crystalline structure the result will actually depend on direction, and
in this case we should write P(r). However, for amorphous structures the result
will be isotropic, so that it is enough to write P(r) for the PDF. In any case, the
meaning of the function P(r) is the following: if there is an atom at the origin,
P(r) describes the probability that there will be an atom at r, averaged over all
atoms being at the origin. For example, if we compute P(r) for a crystal, taking
1' in the direction of one of the edges of the unit cell, we find a narrow Gaussian
distribution about each lattice site. In the extreme case of the perfect gas, where
an atom at the origin has no effect at all on the uniform distribution of other

atoms, P(r) : 1. This is shown in Figure 5.6. Typical shape of the PDF for

100

 

6 . Crystalline GaAs at 10° K

 

Pair Distribution Function
5

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.‘ 11ll11]]11,/)1

I T

0.0 ' 0.5 1.0 ' 1'5 ' 2T0 ' 2f5 ' 3.0 3.5 4.0 4.5
rldo

 

Q

Figure 5.6: Typical shape of the Pair Distribution Function for a crystal, computed
along the one of the edges of the crystalline unit cell. This is the PDF calculated
for the crystalline GaAs at 10K (Chung and Thorpe 1996). Gaussian peaks are
due to the thermal broadening of delta—functions in (5.2), and are taken in the
account by performing thermodynamic average. The horizontal line is PDF for a
perfect gas. The distance has been normalized by the nearest—neighbor distance

amorphous structures is shown in Figure 5.7.

If we ignore the thermal broadening, we can simply write for the PDF

1
P(r) = mtg 6(1‘ - 1‘5), (5.3)

where the r, are the lattice translation vectors. If we integrate both sides of (5.3)

over the volume of the-isample we get

no / P(r)d3r -_- i230 / 5(1- — r,)d3r = N— 1, (5.4)

where the (—1) comes from explicitly excluding the atom at the origin in both
(5.3) and (5.2). Thus, noP(r)d3r gives the probability of finding an atom at a

point d3r which is a distance r from another particle at the origin.

 

 

 

 

 

 

 

 

err—flu 1:— . 1 . 1
7- .
a-SI -

1: .
.2 ”i —
‘2' 4
a 5“ «
LL .
c 4
.2 44 .
H
a 1 1
'5 3‘ 1
.2 +
D
g 21 ‘
(o .
0- 1. —c_,~ L

1

o l I I I I

0 1 2 a 4 5 0

rld0

Figure 5.7: Typical shape of the Pair Distribution Function for an amorphous
sample. This is the PDF calculated for the amorphous Silicon (Djordjevié et al.

1995).

Assuming that the sample is isotropic, we can calculate the average num-
ber of atoms in a spherical shell centered on a particle at the origin, and dr thick,

thus introducing a new distribution function '
[h 11 noP(r)d3r : 41rr2noP(r)dr : noG(r). (5.5)
S e

The equation (5.5) defines G(r), which is called the Radial distribution Function
(RDF) '
C(r) : 4rr2P(r). (5.6)
We see that, on average, the RDF increases parabolically with the distance r.
Typical shapes of the RDF in a crystal, glass and a gas are shown in Figure 5.8.
The Radial Distribution Function is the most used distribution function in
the literature on amorphous solids. The appeal of the RDF is that it gives a rather

direct physical picture of the spatially averaged structure. However, one has to

102

CRYSTAL GLAS GAS
(0) , (b) (c)

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 5.8: Typical shapes of the Radial Distribution Functions for (a) a crys-
talline solid, (b) an amorphous solid, and (c) a gas. These distributions schemat-
ically correspond to the atomic arrangements sketched in Figure 5.2.(Taken from

Zallen 1983.)

keep in mind that the RDF (and other distribution functions from this section)
are only one-dimensional representations of a three—dimensional structure, and
thus contain incomplete information about it. Nevertheless, the RDF is a useful
tool for investigating the complex structure of amorphous materials.

Integrating the RDF over the distances which include just the first peak

of G'(r) we get the total number of atoms in the first coordination shell Z

la... peak G(r)dr = z. (5.7)

The position of the first peak in RDF, r1, gives a value for the nearest-neighbor
bond length, and the position of the second peak r2 gives the next—nearest-
neighbor bond length. It is easy to see that the knowledge of r1 and r2 gives

a value for the bond angle 0

0 = 2arcsin (52-) . (53)

7‘1

103

The width of the first peak is generally caused by thermal vibrations, but in the
static limit (ignoring thermal and zero-point effects) for amorphous solids, it is
due to the deviations in bond lengths. Similarly, the width of the second peak
is caused by the static deviations in the bond angles in amorphous solids. The
origin of structural features in the RDF is schematically shown in Figure 5.9.
There are two more distribution functions used in the literature. The first

one, C(r), is called Pair Correlation Function and is defined through RDF as

 

(5.9)

and thus goes linearly with r. The Pair Correlation Function (PCF) has two ad-
vantages over RDF. First, it permits structural data to be reasonably presented
over a greater range of r, because of its linear behavior. Second, more fundamen-
tally, it is in PCF where experimental broadening is symmetric and r independent
(Wooten and Weaire 1987; Etherington et al. 1982). Thus, PCF represent bet-
ter way of comparison between the theory and experiment, then RDF. A typical
shape of the PCF is shown in Figure 5.10.

The second distribution function, D(r), and the last one we shall introduce

here, is called Reduced Radial Distribution Function (RRDF) and is defined as

 

- 47rrno. (5.10) .

It called reduced because it oscillates about zero since the average density back-
ground curve has been subtracted. Because it is a difference curve and is weighted
by a factor r, in many cases it displays structural correlations more clearly than

the RDF itself (Elliott 1984). A typical shape of an RRDF is shown in Figure 5.11.

104

 

RDF

 

Figure 5.9:. Schematic illustration of the ori 'n of structural features in RDF
for an amorphous solid. Atoms are shown as ying on sharply def :1 rings, for
simplicity. (Wooten and Weaire 1987.)

 

 

 

 

 

 

1 I r 1
6» -
41- _
21- ..
OM L 1 1
O 2 4 6 8 IO
r(K)

Figure 5.10: Typical shape of PCF for an amorphous solid. (Wooten and Weaire
1987)

IT» I
2.-

1

 

 

 

i
1
1

6(1) (atoms/A”)
N v— C 1— N u & UI
l

 

11(4)

Figure 5.11: Typical shape of the RRDF for an amorphous solid. (Taken from
Elliott 1984.)

106
5.4 Structure Factor vs. Pair Distribution Func-
tion

It is well known that the static Structure Factor (SF), for, say, X—ray scattering

in the amorphous sample is given by

N
S1Q) = 4, 1 Ewan-R.) 12. (5.11)

where Q is the momentum transfer vector and R,- are atomic positions in the
structure.

The Structure Factor defined in (5.11) is an important quantity in both
experimental measurements and theoretical calculation. In crystals, 5 (Q) gives
the Bragg peaks and is at the foundation of crystallography. In amorphous ma-
terials, from the experimental measurr nts of the intensity in X —ray scattering,
one can usually determine the structure factor S (Q), from which one can find
PDF, P(r), and compare it with the predictions of theoretical models. Further-
more, S (Q) is an important quantity for establishing the criteria of randomness
in modeling amorphous structures (Wooten and Weaire 1986). Whatever method
we use to create a random structure, it is useful to calculate S (Q) for the struc-
ture, along one particular direction in Q—space. As we said, in a crystal, S (Q)
is obtained as a set of delta functions located at certain values of Q, which cor-
respond to the Bragg peaks along the chosen direction in Q—space. The sample
is well randomized if S (Q) becomes smeared over all values of the wave vector Q
equally, so that, ideally,there are no remnants of the Bragg peaks in the original
crystalline structure, which become small and comparable to values of S (Q) for
any other Q. We have used this criterion in creating our computer models of

amorphous diamond, and amorphous silicon. In Figure 5.12 we show 3 (Q) for a

107

crystlline, and also for several randomized structures (Wooten and Weaire 1986).

We now want to connect the Structure Factor, (Q) with the Pair Distri-
bution Function, P(r). For that, we have to calculate the configurational average

of S(Q) for the real sample (see Cusack 1987, pp.56-60).

Let us introduce the single—particle density function, 1/(r)

u(r) = Ebb—R.) (5.12)

1:1
It is obvious that

/u(r)dr = N. (5.13)

Introducing 1/(r) into (5.11) it can be shown that S(Q) can be alternatively

written as
5(Q) =1iV/V(r1)€$P(iQ°r1)dr1 >< /V(r2)€l‘P(iQ-r2)dr2- (5-14)
Rewriting (5.14) as
S(Q) = j—tf-/u(r1)u(r2)exp[iQ - (r1 - r2)]dr1dr2, (5.15)

changing the variables so that r2 —+ r’, r1 — r2 —+ r, and after some algebraic

manipulation we get
1 . ..
5114) = ,7 / F(r)e:1vp(1Q- r), 1416)
where, using the properties of delta function, the function F (r) is obtained as

F(r) = 260 — (R.- — Rn]. (5.17)

So, to calculate the average S (Q) from (5.16), one has to calculate the

average of the function F (r)

108

 

 

1,0 -101 _4

0.0 )-
0.0 b

222.1 1-._1_--L_1__-1-__1-_§ 1

T 1 1 1 - 1
1 l . l l I

 

 

 

 

0.“
0.04
0.“
0a
0.01

. 1...11..,1..1111,.1.11.11.11111I].ud.uu._1.£

((31 I i

 

 

q

0.010

0.010

- (11111191111111.114114:
.... - 1 ..

«35.111 1.1 1111.141111111

 

 

T'I

 

 

 

 

 

 

 

L
10 1.

«(1111 1A“)

 

 

 

 

Figure 5.12: Values of the S (Q) computed along (111) direction of the 216-atom
periodic Si unit cell. (a) The crystalline diamond cubic structure. (b) - (11) More
and more randomized structures (Wooten and Weaire 1986).

109

Going back to the definition of PDF, (5.2), we note that the atom at the
origin is excluded in the definition of PDF, while all atoms are taken in the sum
in (5.17). Then it is obvious that we can write

(F(r)) = (2 5f? - (R. - RM).

id
= (Z 6[r— (R, —R,-)]) + Nam. (5.18)
#1
The first term in the second line of (5.18) is just equal to Nn0P(r), as given by

(5.2), so that we can write
(F(r)) = noNP(r) + N6(r), (5.19)
or for the average of the Structure Factor,

1510)) = (172.11%) — 1) + n. + 6(r)]erp(iQ - r)dr (5.20)

which gives, using the properties of the delta function,

(S(Q)) = 1+ no /V[P(r) — 1]e:1:p(iQ - r)dr + (2«)35(Q). (5.21)

The term 6 (Q) in (5.21) implies strong scattering intensity in the forward direction
with I Q I: 0 : 0 [see Eqn. (5.1)]. This peak would fall in the incident beam and
not be observed in pratctice. In any case, this term should be replaced by another
term which properly takes in account the shape and the size of the sample in the .
integration (5.20) (see Guinier 1955), but we can neglect this term anyway because
it does not contain information about the structure and is negligible unless Q is
very small in which event 0 would be too near the incident beam direction for
observation.

Transforming (5.21) into spherical coordinates and separating the vari-

ables we finally get the connection between the Pair Distribution Function and

110

the averaged Structure Factor. 1(Q).

 

[(62) a (5(a)) = 1 + 11.0/(:0[P(r) — 1]8i22?r47rr2dr. (5.22)

we assumed in (5.22) that the system is isotropic. Noting that P(r) —+ 1 as
r —100, we can see from (5.22) that [(Q) —+ l as Q —+ 00.

Finally, we note that inferring [(62) from the observed intensity in scat-
tering experiments is not always straightforward. For the detailed analysis of

the problems regarding this problem, the reader is referred to Etherington et al.

(1982).

Chapter 6

Historical Background of

Modeling Random Networks

Modeling of amorphous structures can be divided in two large classes. The first
one is the Dense Random Packing of hard spheres, (DRP), in all its varieties, and
the second is based on the Continuous Random Network, (CRN), concept. The
DRP model was very successful in explaining the amorphous structure of some
liquids, amorphous metals, ionic solids and molecular organic solids. Atoms in
these systems are connected with each other predominantly by non—directional
forces, which makes DRP a suitable model for the description of these materials.
DRP modeling has been pioneered by Bernal (1959,1964), and subsequently de-
veloped by Finney (1970), Bennett (1972) and others. For a review of modeling
amorphous metals with DRP models see an excellent and comprehensive article
of Cargill (1975). In this Thesis, however, we shall focus our attention only on
the Continuous Random Network models which are suitable for describing the
structures of most non—metallic materials with predominantly covalent directed

bonds, namely, covalent glasses. Covalent glasses posses much higher degree of

111

112

short-range order than metallic glasses and it is dominated by covalent bonding
with definite number of nearest neighbors, bond lengths and bond angles. The
prototypical covalent amorphous solids are amorphous silicon, (a — S i), and amor-
phous germanium, (a — Ce) which belong to the semiconductors of a Group—IV.
Every atom in these structures has for nearest neighbors, i.e. it is tetrahedrally
coordinated.

We note here that the invention of the term continuous in connection to
random networks is somewhat unfortunate, because random networks are, in fact,
paradigmatically descrete in their structure. By ”continuous” it is meant that one
can build an amorphous network continuously, starting from the short—range order
unit, to indefinite size without including unsatisfied bonds or breaking the struc-
ture. In other words, ”continuous” means that there are no identifiable boundaries

in the structure, separating regions of distinctly different type of structure. '

6.0.1 Hand—built CRN Models

The first attempts to explain the structure of covalent glasses, around 1930, were
based on, then very natural, hypothesis that amorphous materials actually con-
sist of very large number of elemental microcrystals, randomly arranged into a
very fine polycrystalline structure, which appears as amorphous. All attempts to
fit the experimental data with such microcrystallite models failed, and the early
idea of a continuous random network of Zachariasen (1932) became soon the only
viable alternative. It was proposed by Zachariasen that ”...the atomic arrange-
ment in glass is characterized by an extended three—dimensional network which
lacks symmetry and periodicity.” In Figure 6.1 we show the historical schematic

diagram of Zachariasen which represent his CRN model of amorphous silica. But

113

 

Figure 6.1: Zachariasen’s schematic 2d diagram represents his CRN model for
amorphous silica. Open circles are two—fold coordinated oxygen atoms, and black
dots are S' i atoms which are four-fold coordinated in 3d (Zachariasen 1932).

114

the first models based on the Zachariasen’s ideas were built much later. One of
those was the one for Si04, done by Bell and Dean (1966,1972). This was a
hand—built model and is shown in Figure 6.2. This model was quite successful in
explaining the diffraction data. In building the model, Bell and Dean used basic
units consisting of a central Si atom and four rigid wires in perfect tetrahedral
arrangement. The O atoms were attached by inserting the wires into pre—drilled
holes in the 0 atoms. The model was built to have a mean Si — 0 - Si bond angle
close to 160°, but without any particular bond angle distribution. It maintained
full coordination without dangling bonds in its interior. The model consisted of
614 atoms. Atomic coordinates were determined with the help of photographic
techniques. From those coordinates the rms bond—length deviation was small
in agreement with the experiment, while the bond—angle deviation was approxi-
mately i30° around the mean value of 153°. The Radial Distribution Function
was computed and was in good agreement with the diffraction data of Mozzi and
Warren (1969), as shown in Figure 6.3.

Very soon, amorphous elemental semiconductors were intensively studied.
The amorphous Group-IV semiconductors, a — Si and a - Ge, came to be regarded
as prototypical covalent amorphous solids. It was not clear in the begining whether
tetrahedral bonds in such systems can be randomly connected ad infinitum in
a disordered manner, or not. It was believed by Phillips (1979,1981), that for '
instance, amorphous silicon, due to its overconstrained structure (a notion that
will be explained in detail later), must necessarily have a discontinuous structure,
in which regions of continuous random network are limited in size. But it is
obvious that, in principle, there is a very close connection between a — Si and

amorphous silica. A random network model for one can be easily converted into

115

 

Figure 6.2: The hand—built model of Bell and Dean for amorphous silica (Bell
and Dean 1966).

116

 

    
 

Cfl'iGBELLBDEAN

 

RDF i .

«1-
-1)-
«1)-

 

 

 

L l
O 2 , 4 6 8
. . .
r 1 A)
Figure 6.3: Comparison of the RDF for S 10.. derived from the model of Bell and

Dean (1972) with the experimental reults of Mozzi and Warren (1969). (Taken
from: Zallen 1983.)

 

117

a corresponding model for the other, by adding or substracting the oxygen atoms
which link the silicon atoms and introducing a weak bond—bending force at oxygen
sites.

The first CRN model for a tetrahedrally bonded amorphous solid was built
by Polk (1971) and is shown in Figure 6.4. There were no particular rules used in
building Polk’s model, except that there were no dangling bonds left within the
interior of the structure. Whenever a choice about the bonding of a new atom
should be decided, simply an arbitrary, random choice was made, just as in the
work of Bell and Dean. The model contained substantial number of five—, six—
and seven—membered rings which ensured generation of a non-crystalline struc-
ture. There were no dangling bonds in the interior of the structure. Bond—length
deviations werere restricted within 1%, and bond—angle deviations were within
i10° about the tetrahedral angle of 109°. The model consisted of 440 atoms. In
Figure 6.5 we show the RDF for Polk model compared to the experimental result
of Moss and Graczyk (1969).

The original Polk model was physically extended to 519 atoms and refined
by Polk and Boudreaux (1973). The refinement was in reducing the deviation of
nearest-neighbor bond lengths from 1% to 0.2% using precise laser—beam mea-
surements of the atomic coordinates and computer techniques. Steinhardt et al.
(1974) and Duffy et al. (1974) used 519 model of Polk and Boudreaux (1973) and
refined it further. This refinement consisted of introducing the Keating (1966)
potential in the structure and relaxing the system toward the minimum of the
elastic energy in terms of the bond lengths and bond angles.

The method of relaxation adopted by Steinhardt et al. (1974), which we

also used in creating our models of amorphous diamond, is to move each atom

118

 

5 original hand—built model for amorphous silicon (Polk 1971).

1

Figure 6.4: Polk

119

 

‘—
B-
K1
0:

 

 

 

WERIKH

 

 

 

 

 

 

-
—

 

 

 

Figure 6.5: Comparison of RDF for Polk model of amorphous silicon (1971) with
the RDF deduced from the experimental measurements of Moss and Graczyk

(1969). (Taken from: Zallen 1983.)

120

separately and sequentially, while keeping the coordinates of the other atoms
fixed, until there is no force on it. This process is iteratively repeated until
the convergency is achieved. In this way Steinhardt et al. found new, relaxed
coordinates of 519—atom model of Polk and Boudreaux. On the other hand.
they also built a new 201—atom model starting from a 21—atom cluster with new
atoms serially added. Whenever a small group of atoms was added to the growing
cluster their positions were relaxed by computer to minimize the forces on them.
Thus, the connectivity of this network is determined by the relaxation procedure.
They did not find any significant structural difference between the two models.
Duffy et al. (1974) used the same 519—atom model of Polk and Boudreaux but
adopted different relaxation procedure. Each atom was simultaneously moved in
the direction of force on it, a distance proportional to that force. After each move
the forces are recalculated, the atoms moved again, and the process iterated until
the forces and displacements become small. In this way they obtained results
which were pretty much the same as those of Steinhardt et al.. In Figure 6.6 we
show the RDF obtained from the relaxed 519—atom model of Steinhardt et al.
(1974), compared with the experimentally derived RDF.

Provoked by the controversy involving the presence of five—membered
rings in amorphous III — V compounds like GaAs, which would, of course im-
ply the occurence of energetically unfavorable, ”wrong” bonds between like atoms,
and to check if different ring statistics would produce different RDFs, Connell and
Temkin (1974) built a model of a —— Ge with no odd—membered rings. The model
consisted of 238 plastic and metal tetrahedra. This model did not produce any
improvement of the agreement with experiment. Our own calculations of RDF

for a - Si and a — C show that RDF is very similar for different ring statistics

 

8

N
U

RDF (nous/i)

 

 

 

 

Figure 6.6: Comparison of the RDF of a relaxed CRN model of a — Ge by Stein-
hardt et al. (1974) with the experimentally derived curve. (Taken from: Zallen

1983.)

122

(Djordjevié et al. 1995).

It is interesting to mention that there was an early effort to build a model
for amorphous diamond. This was hand—built model of Kakinoki et al. (1960)
which consisted of about 100 atoms. To construct an ”amorphous” structure,
the authors used two basic units, the first with tetrahedral symmetry of diamond
lattice, and the second with the trigonal symmetry of graphite. The C -— C
bond lengths were different in those two units, and equal to those in crystalline
diamond and graphite, respectively. To reproduce electron diffraction data the
individual regions in the structure dominated by either diamond—like or graphite-
like structure, were supposed not to be of larger size than a few Angstroms. To
our knowledge no detailed analysis of this model have been done, but we could
argue that this model is conceptually wrong for describing an amorphous solid.

All these early hand-built models suffered from the their finite sizes, or
in other words, from the lack of periodic boundary conditions. In these relatively
small models, a large fraction of the atoms were necessarily located close to the
surface. This complicates the analysis of the bulk properties, and calculation of
the RDF. The first hand—built model which had the periodic boundary conditions

implemented was made by Henderson (1974). The model consisted of 61 atoms.

6.0.2 Computer-built CRN Models

Maybe the very first computer modeling of this kind was done by Kaplow et
al. (1968). They tried to simulate the structure of amorphous selenium and
to compare it to the experimental results. The method Kaplow et al. adopted
in creating computer model is a good example of a misconception that one can

make an amorphous structure just by randomizing the atomic coordinates in a

123

crystalline structure. This is exactly what the authors did in this early work. They
computer—generated coordinates of crystalline arrangement of about of 100 atoms,
and then used Monte Carlo procedure to randomize these coordinates. Only
those which improved RDF were accepted. After about 10° moves a ”reasonable”
agreement with experiment was achieved. The problem with this procedure, as we
have already mentioned before, is that obtained structure is topologically identical
to the crystalline one, and thus can not be expected to give good description of
the amorphous selenium.

Henderson et al. experimented with computer algorithms for creating
structures with periodic boundary conditions, but they did not succeed in getting
good agreement with experiment. Nevertheless, the approach of Henderson, with
generating periodic structures and with using a simple computer algorithm, in-
spired later work on computer modeling. In contrast to this approach, Shevchik
(1973) and Polk and Boudreaux (1973) tried to imitate in their algorithm the
hand—building process and that proved to be very cumbersome procedure.

As we have already stressed, a truly amorphous structure created by a
successful randomization of the crystalline atomic positions, must be topologically
different from the structure of its crystalline parent, and it must have perfect
connectivity. One obvious technique to ensure these two requirements, for instance
in‘ creating amorphous silicon, is to start from the diamond cubic structure, remove
one atom at random, and rejoin the four dangling bonds in pairs. In this way the
number of atoms in the structure clearly decreases, but the topology is changed
due to the fact that six—membered rings of the crystalline structure are destroyed,
while five—, and seven—membered rings are introduced. Four-fold connectivity of

the network is preserved. This method was used by Alben et al. (1975) who used

124

64—atoms crystalline supercell as a starting point, and ended up with a 58—atoms
model of amorphous silicon. In an unpublished work of Wooten and Weaire (see
Wooten and Weaire 1987) similar studies were done on models that contained
700 — 1000 atoms. Results were less then impressive. Calculated RDFs are in
poor agreement with experiment. The structures obtained in this way have large
bond—angle distortions, and it was concluded that they are not satisfactory as
models of amorphous silicon.

The other technique which illustrates the topological connection between
the structure of vitreous silica and amorphous silicon, was employed by Evans et al.
(1974, 1975). They used a decoration transformation (see Wright et al. 1980) of
the vitreous silica random network model of Evans and King (1966) which involved
removing the oxygen atoms, rejoining the two dangling bonds and adjusting the
new network to have the meand bond length appropriate for amorphous silicon.
The model achieved relatively good agreement with experiment in the first two
peaks, but the model pair correlation function exhibits too much structure at
higher separations.

The first approach to computer generation of periodic random networks
was extensively developed by Guttman (1974, 1976). His aim was to construct
an ideal random network with four—fold coordination which would describe the
structure of amorphous germanium. Guttman does not use a single computer
algorithm, but otherwise his approach is quite similar to the Wooten and Weaire
method (1987) which we used in our modeling of amorphous diamond. Guttman
starts from the crystalline super—cell which contains many unit cells of the crys-
tal], periodic boundary conditions are applied. There are no bonds in the model

initially. Bonds are introduced in the stochastic way: four neighbors are selected

125

at random from the given set of atoms which is assigned to each atom in the lat-
tice. Thus four bonds are constructed at a time, and the process is repeated until,
hopefully, the desired network is created. Guttman noticed that the probability
of constructing a fully four-fold coordinated network in this way was decreasing
with the size of the model. The size of the model varied from 64 atoms to 512
atoms. The structure was relaxed using a simple Hooke’s—law bond—stretching,
and bond—bending term. Guttman also varied the super—cell size to find the op-
timal density which minimizes the total strain energy. He found that the optimal
model density was 15% greater than the observed crystalline density, while the
density of a —- Ge experimentally was found to be within the few percent of the
crystalline density. Furthermore, a large compressive stress was introduced by
random bonding assignments, and the rms angular deviation from the perfect
tetrahedral angle (00 : 109.47) was found to be about 17°; much greater than the
experimental value of abut 10° (see Etherington et al. 1982). Guttman then used
a specific topological rearrangement to generate random network. In Figure 6.7
we show two examples of Guttman’s bond rearrangements, or bond switches. The
switch is accepted if it lowers the energy of the network, and rejected if it does not.
In this way RDF was slightly improved, but the models still did not agree very
well with experiments, because all Guttman’s models contained many dangling

dangling bonds.

6.1 Wooten—Weaire Method

Before we present the Wooten—Weaire method which we adopted in our computer
modeling studies of amorphous diamond, we shall make a small digression here to

describe the characteristics of the different sizes of models made for tetrahedrally

126

 

Figure 6.7: Two examples of bond switches used by Guttman (1976). Six-
membered rings of the crystalline lattice are shown, and four atoms are selected.
A path 1 - 2 - 3 — 4 is defined and bonds 1 - 2 and 3 - 4 are exchanged with
non-bonds 2 - 3 and 4 - 1 respectively.

 

 

 

 

 

 

 

 

Figure 6.8: Crystalline diamond cubic cell.

bonded amorphous semiconductors which in the crystalline phase have diamond
lattice structure. In figure 6.8 we show the unit cubic cell of the crystalline dia-
mond. A super-cell of the crystalline diamond can be created by adding additional
unit cellsialong the three coordinate axis. Thus, if there are 2 unit cells per axis,
that creates a super-cell of 8 unit cells, each containing 8 atoms, which gives
64-atom model. Similarly, if there are 3 unit cells along each axes, that gives
3 x 3 x 3 x 8 : 216-atom model. In this way, models consisting of 64, 216, 512,
1000, . . . , 4096 [appear in the literature. Thus it is easy to tabulate the linear sizes
of the various super-cells. For instance, in Table 6.1 we show the list of linear

sizes, L, of the super—cell cube for Si and C, given in Angstroms. The numbers

128

Table 6.1: Linear super-cell sizes, (L),in Angstroms, for Si and C used in modeling
of a - Si and a — C (Djordjevié et al. 1995). N is the number of atoms in the
model.

[V L 5, LC

 

 

216 16.2813 10.7004
512 21.7084 14.2672
1000 27.1354 17.8340
4096 43.4167 28.5344

 

 

in Table 6.1 are used to calculate the number density N / L3 of the model before
one allows the size of the super—cell to vary untill the total strain energy is at

minimum.

6.1 .1 Keating Potential

In the Wooten and Weaire model of a — Si and a — Ge, the Keating potential is
adopted to describe the interaction between the atoms. Keating potential was in-
troduced eariler to fit the elastic and vibrational properties of Group-I V elements
(Keating 1966). It represents a semiempirical description of bond—stretching and
bond—bending forces and involve only few parameters. It consists of two terms,
the first being the bond-stretching one, and the second the bond-bending.
2

V = T36% a: (1‘11 ' 1'11 - T1202 + 2% 1021;} (1'11 ° 1‘11! + 2T3) , (6.1)
where a and ,6 are bond-stretching and bond—bending force constants, respec-
tively, and d is the strain—free equilibrium bond length in the diamond structure.

For Si its value is d : 2.35A, and for Ge d : 2.46A. Note, however, that the

129

second term in (6.1) depends not only on the angle between the bonds, but also
on the bond lengths, and thus contains cross terms. The first sum in (6.1) is on
all atoms 1 and their four neighbors i. The second sum is on all atoms and pairs
of distinct neighbors and thus represents three—body part of the potential.
Central bond-streching forces are not enough to stabilize the structure,
because the number of constraints imposed by maintaining fixed bond lengths
is smaller than the number of degrees of freedom. This will be explained in
more detail in Chapter 8, where we shall discuss the elastic properties of diluted
amorphous diamond. At the moment it is enough to stress the importance of the
bond—bending forces for creating a stable network. Bond—bending three-body
force could be, in principle, simulated by introducing an extra fictitious bond
for each 5 independent angular constraints. This is a sometimes useful strategy
for investigating rigidity problems in networks, but for simulating the amorphous
structures it is more realistic to use angular term in Keating potential. In practice,
bond-stretching forces are much larger than bond—bending forces. Wooten and
Weaire followed Martin (1970) by taking fl/a : 0.285 in building their model of
a — Si and a — Ge. Due to the dominance of the a term, when the model is
relaxed to minimize energy, all bonds relax to within a few percent to their ideal
unstretched length, but because this is not enough to stabilize the structure, the 3
term defines the stable structure which minimizes the bond—bending energy. Thus,
the total energy is roughly proportional to 3, and relaxed structure is independent
of fl/a (Wooten and Weaire 1985; 1987). The Keating potential is schematically

show in Figure 6.9.

130

 

Figure 6.9: The Keating potential bond—stretching (a), and bond—bending ([3)
terms. (Taken from: Wooten and Weaire 1987.)

6.1.2 Wooten-Weaire Bond Switch

Wooten et al. (1985) created a new model of amorphous silicon and germanium
which significantly improved the agreement with experiment. It consisted of 216
atoms. A new be 4-, 'itch was invented which does not produce large bond—-
length and bond-angle dostortions, and is simple and the only one used in the
two—stage process which creates the model. This new bond switch is shown in Fig-
ure 6.10. The bond switch presented. in Figure 6.10 ischaracterized by switching
two second—neighbor bonds that are parallel to each other in the perfect diamond
structure. This is a local rearrangement of bonds that creates minimal strain
into originally strain-free perfect diamond structure, which was not the case for
Guttman’s switch (Figure 6.7). Atoms are numbered, their coordinates are listed,
and neighbor-table is created. Bonds 1 — 2 and 5 - 6 in Figure 6.10 are parallel,

and can be a candidates for switch. The switch consists of breaking those two

131

 

Figure 6.10: Wooten—Weaire bond switch. (0) Configuration of bonds in the
diamond cubic structure. (b) Relaxed configuration after switching two bonds.

132

bonds, and creating two other bonds, 1 — 5 and 2 —6. The switch is not acceptable
if new bonds are longer than 1.7 times the unstretched nearest—neighbor bond,
which is the length slightly larger than the second—neighbor distance, as it should
be if a successful bond switch is to be made in the original crystalline structure.

One can define a ring as any closed non-returning path of bonds. The
number of n—membered rings defined in this way increases with n. An irreducible
ring is defined as one that has no shortcuts across it. That is, given any two
atoms (vertices) on the ring, there is no shorter path between the two atoms
(as measured by the number of bonds along the path) than a path on the ring
itself. One advantage of such rings for topological purposes is that the number of
n—membered rings goes to zero for large n, but no topologically important rings
are omitted, and the complete table of ring statistics remains finite (Wooten and
Weaire 1996).

The perfect diamond structure has only 6— membered irreducible rings.
The introduction of a bond switch changes the ring statistics by destroying even—
membered rings and creating odd—membered rings, 5-, 7—, 9—, etc. membered
rings. This process is schematically shown in 2d case in Figure 6.11.

When the structure is already randomized, bonds that are switched are
chosen to be approximately parallel by ensuring that they do not belong to the

same 5—, 6—, or 7—membered ring.

6.1.3 Initial Randomization of the Crystalline Lattice

Taking the crystalline super—cell of 216 atoms, the initial randomization is achieved
by introducing a large number of bond switches in a random manner, at the tem-

perature just above the melting point. Note that the temperature used in the

133

 

Figure 6.11: 2d schematic representation of a Wooten-Weaire topological bond
switch. There are four 6—membered rings initially (diamond structure), and after
a bond switch is made two 5—, and two 7—membered rings are created.

134

present context is the effective simulation temperature and it does not correspond
to a real temperature. It is meaningful only within the simulated annealing en-
vironment which will be described in the next Section. Similarly, in the context
of the simulated annealing, the melting point is defined as that temperature for
which the mean square displacement of the atoms from their positions in the
diamond structure increases linearly with time and the structure factor [(1.2), in-
troduced in Chapter 5. approaches a value of order l/N. This differs from the
true melting behavior of a - Si and a — Ge in that tetrahedral bonding is pre-
served. For a - Si and 1’3/01 : 0.285 the melting point is kT : 1.0eV (Wooten and
Weaire 1987). The temperature appears here through the inclusion of Boltzmann
factor, erp(—E/kT), which describes the probability of the acceptance of a given
bond switch. The inclusion of a Boltzmann factor in the initial randomization
process is inspired by its use in the second stage of the model building which
will be described in the next Subsection. The melted structure has large bond—
angle deviations, and the subsequent process of annealing is supposed to decrease
those deviations to their experimentally observed values. That means that the
melted structure has a large amount of strain which is partially released in the
annealing process. Bond-length deviations are not so big in the melted structure
because local cluster of atoms around every new bond switch is relaxed under
Keating potential. We stress again that the initial randomized structure is not
really melted; it still has tetrahedral bonding and it is still solid. We use the term
melting conditionally, only within the simulated annealing environment.

One bond switch introduced into the crystalline diamond super-cell af-
fects the topology of the structure by destroing 12 six—membered rings, and si-

multaneously creating 4 five—membered, and 16 seven-membered rings. During

135

the process of initial randomization, RDF becomes gradually featureless, reflect-
ing, from the very beginning, the amorphous nature of the created structure. In
Figure 6.12, Pair Correlation function is shown for various stages of the initial
randomization (Wooten and Weaire 1987). The initial randomization must be
accomplished at the temperature which is high enough to prevent the returning

to the crystalline phase after the subsequent annealing process.

6.1.4 Annealing: Geometrical and Topological Relaxation

Geometrical relaxation can be achieved in different ways. The conjugant gradient
method is an obvious choice, but for these tetrahedrally bonded structures the
simplest relaxation method, used first by Steinhardt et al. (1974), is sufficient. In
this method, as we have already briefly mentioned, after a bond switch is made,
the force (vector) on each atom involved is calculated, and also the force gradient
(tensor), and each atom is sequentially moved toward its equilibrium position
while other atoms are kept fixed. The vector displacement, Ar, of each atom is
found from a Taylor series expansion for the force about the minimum energy
atomic position

W, = —(V2Vr)Ar (6.2)

which requires a matrix inversion of the gradient tensor. This process is repeated
many times until the convergence to equilibrium is achieved.

After the initial randomization produced a disordered, but highly strained
model, it is necessary to continue with creating new bond switches searching for
those which will lower the strain energy. Of course, the lowest energy possible is
zero, i.e. that of crystal. So one has to be careful to avoid returning to the crys-

talline state during the process of annealing. Wooten and Weaire (1987) showed

136

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

111 ' ' V
12 _ (bl
l 01-
9‘7
8 1r
: ‘1- -
1 .=
0 r— "
2 "’ ._.
ol 1 1 _ o l l 1
0 2 4 I I 10 0 2 4 0 I 10
. I r .
(01 l l (d) V r I 1
'9-
2 P —t 2 -- __‘
0 L l L o. 1 L L
0 2 4 0 I ‘0 0 2 4 0 I 10
. r I r T . I I I I
101 10
1" -1 I r- _.
$7
‘ "" " E ‘ l- -I
3
o L 1 1 o L 1 1
0 2 4 I I 10 0 2 4 0 I 10
11M 11M

Figure 6.12: Pair Correlation Function for different amounts of randomization. (a)
The ideal crystalline diamond structure. (b)-( f), PCF for different concentrations
of bond switches per atom, c : 0.06, 0.12, 0.18, 0.24, and 0.30, respectively. (Taken
from: Wooten and Weaire 1987.)

137

that the system indeed returns to the original crystalline diamond structure if in
the initial process of randomization the number of introduced bond switches was
not large enough. They found that if the number of bond switches per atom that
randomizes the lattice in the beginning is c > 0.6, annealing does not return the
system back to the crystalline structure, and the energy of the final relaxed struc-
ture is comparatively low. Here c is the number of bond switches per atom. For
c z 0.3 the system did return to the crystalline state. So, to make sure that the
long process of annealing will lead to an acceptable final amorphous structure, one
has to deal with two initial parameters: theomelting temperature, and the number
of bond switches that is introduced in the crystalline structure at the beginning.
If the initial temperature is too high the subsequent annealing is not able to bring
the system down to an acceptably low energy, and the obtained model has too
much strain in it. In other words, with too high initial temperature the system is
trapped in one of the numerous metastable states of relatively high energy. From
the other side, the number of bond switches has also to be chosen optimally: if
there are too many of them in the beginning, than again, the system will not be
able to release a sufficient amount of energy during the annealing process, and the
RDF of the final structure will show that the system has too large bond—length
and bond-angle distortions. Wooten and Weaire found that for creating a - 52'
model the optimal melting temperature for the initial process of randomization
has to be kT = 16V. As we explained in the previous Section, this 10, 000K high
temperature is meaningful only within the simulated annealing environment.
The annealing process itself is based on Metropolis Monte Carlo algo-
rithm (1953) already applied in the literature to general problems of optimization

by simulated annealing (Kirkpatrick et al. 1983; Vanderbilt and Louie 1984). The

138

Metropolis algorithm is generally used as an efficient simulation of a collection of
atoms in equilibrium at a given temperature. In the process of creating a model
of an amorphous structure, the introduction of topological rearrangements (bond
switches) change the total strain energy of the lattice at given temperature. If
that change, AE = E,.,,,- — E, (where E is the strain energy after a bond switch
is introduced) is negative, i.e. if the introduction of a bond switch lowers the
total energy, this bond switch is accepted. On the other hand, if the AE is posi-
tive (increases the total energy), that bond switch is still accepted, but with the
probability P(AE) = exp(—AE/kT). In practice, a random number uniformly
distributed in the interval (0,1) is chosen each time when a bond switch is intro-
duced, and AE is calculated. If the picked random number is less than P(AE)
than the new configuration of atoms is accepted. If the random number is greater
than P(AE) the new configuration of atoms is rejected and the original one (be—
fore the bond switch was introduced) is used to start the next step. Of course, the
next step always consists of introducing the new bond switch, i.e., in repeating
the described basic step. When this is done many times, the procedure simulates
the behaviour of the system at given temperature T, as it reaches the thermal
equilibrium. The choice of P(AE) in the exponential form has the consequence
that the system evolves into a Boltzmann distribution. It should be noted here,
that the temperature T in the Metropolis algorithm is simply a control parame-
ter, not the real temperature of the system. In making a model of an amorphous
structure one essentially deals with the static structures where all dynamics the
non—zero physical temperatures introduce in the system are neglected. The only
place when the finite temperature effects come into a consideration is the broad-

ening of the calculated raw RDF by convoluting it with Gaussian to simulate the

139

existing thermal broadening of the measured RDF in experiments (Etherington
et al. 1982). In all other instances when we speak about temperature, we mean
by it just a parameter in the probability function in Metropolis algorithm.

The introduction of the Boltzmann factor is essential in the process of
annealing because it is able to lift the structure out of a metastable state of
higher energy, thus enabling it to find another path in the configurational space
toward the lower-energy states. Without it the system would not be permitted to
ever increase its energy, and thus it would be easily trapped in a metastable state
with no way out. It has been noted (Wooten and Weaire 1987) that allowing
the presence of the 4—membered rings in the annealing process facilitated the
convergence toward the energy minimum, and helped to avoid trapping into high-
energy metastable states.

Thus, after the structure has been initially randomized by the introduction
of large number of bond switches, the temperature is then reduced in small steps,
and at each temperature the system is kept long enough by introducing new bond
switches and using Metropolis algorithm until the system is in the equilibrium.
As we already mentioned this is followed by the Keating relaxation of the local
regions affected by a bond switch. A sequence of the temperatures is chosen, and
at the end, when the temperature is zero, one hopes to obtain the structure with all
characteristics of amorphous materials. Wooten and Weaire (1985; 1987) obtained
in this way good model for a — Si and a — Ge. Their final 216—atom structure
had the rms bond-angle deviation about 12.6°, which is a significant reduction
from about 22° which the initial randomized structure had. Consequently, the
energy of their amorphous models was comparatively low. We shall see in the

next Chapter that our 4096-atom model of a — Si exhibits even lower rms bond—

.140

angle deviation (10.5° — 11°) which corresponds to the lower strain energy per
atom than in Wooten and Weaire models. Nevertheless, some physicists object
that there is still too much strain present in the amorphous structures obtained
by using Wooten—Weaire method, and that a modification of it is needed if one
wants to produce more realistic, less strained amorphous structures (Mousseau
1996 (a)).

In Figure 6.13 we show the Pair Correlation Function as it evolves in the
annealing process (Wooten and Weaire 1987). The Wooten and Weaire model
of a - Ge and a — Si (1985) has a slight remnant of the original crystal visible
in the structure factor, but nevertheless, there was no serious discrepancy with
experiment (Wooten and Weaire 1987). The slight memory of the diamond cubic
structure that remained in the amorphous structure is the consequence of the
insufficient initial randomization and not from the annealing process. It should
be noted here that the only essential difference between the a - Si and a — Ge
.is a simple scaling in proportion to the bond lengths (d = 2.35Afor Si, and
d = 2.46Afor Ge). Otherwise, the experimental data are practically identical for
a - Si and a — Ge. In Figure 6.14 we show the pair correlation function obtained in
Wooten and Weaire model for a — Ge compared with experiment. The discrepancy
in PCF at high r are related to the presence of inhomogeneities in the real sample
which are not taken in account in fully tetrahedrally connected models, which

always have slightly higher density than samples produced experimentally.

6.2 Further Developments

Wooten and Weaire (1989) used their method to generate a model of the interface

between the crystalline and amorphous silicon. Their model consisted of 512 Si

141

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

r(A)

Figure 6.13: The pair correlation function at three points in the annealing process
(Wooten and Weaire 1987). The rms bond-angle deviations, from tOp to bottom,
are 15.6°, l3.2°, and 12.6°, respectively.

142

 

 

 

_ Calculated

_ Experiment
0
I ' "‘ 1Ii 1 1 fl

0 2 4 6 8 10
r (A)

 

 

 

 

_

Figure 6.14: Comparison of PCF for a — Ge. Thin line is the Wooten and Weaire
model (1985), and thick line is experimental.

H3

 

Figure 6.15: A computer—generated picture of the crystalline/amorphous silicon
interface model (Wooten and Weaire 1989).

atoms and was partitioned into two regions. One region of 191i atoms bounded by
( 100) interfaces was held at zero temperature. The remaining region of 320 atoms
was randomized at [CT = 16V and then annealed.‘ In Figure 6.15 we show the
generated interface model. The structure factor drops to a value characteristic
of the bulk amorphous region over a distance of about 3A. The rms bond—angle
deviation in the amorphous region is‘ 12°. After the annealing process was done.
the whole structure is relaxed under the Keating potential. Atoms in the first
crystalline layer had the rms displacement of about 0.21.3“ It was noted that if
the cell was held at half the melting point for a time much longer than the time

required to make good interface model. the amorphous region crystalizes.

144

Another application of the Wooten—Weaire method was done by Mousseau
and Lewis (1990) in their computer simulations of amorphous silicon hydrides.
The starting point is a model of amorphous silicon obtained by Wooten—Weaire
method. Instead of the Keating potential the authors use Stillinger—Weber poten-
tial (Stillinger and Weber 1985) to describe the interaction between atoms. They
introduce pairs of H atoms in the a — Si structure in such a way that the the
fourfold coordination is preserved. The nearest—neighbor Si atoms were allowed
to have to each bond to either one, or two H atoms. After each introduction
of hydrogen pairs, the structure is Monte Carlo relaxed using a finite number of
Wooten—Weaire bond switches, and relaxed using the Stillinger-Weber potential.
The process is repeated until the desired concentration of H atoms is achieved.
Mousseau and Lewis calculate Si — Si, Si - H, and H — H partial pair distribu-
tion functions, as well as corresponding static factors, and find that the agreement
between the model and experiment improves with H concentration.

As we have seen, the Wooten—Weaire annealing method for generating
amorphous networks does not have a firm prescription of how one can obtain
a good model at the end of a long simulation process. The initial randomizing
(melting) temperature, the number of the bond switches introduced during initial
randomization, the time system is to be held at a given temperature, the sequence
of temperatures, these are some of the parameters one has to play with during
the creation of a model. This requires certain amount of intuition and experience
that one has to develop to create a successful model. The authors colloquially
call this method shake and bake (Wooten and Weaire 1987). There have been
attempts to improve the efficiency of Wooten—Weaire method by designing an

optimal temperature schedule for the annealing, as well as to increase speed of the

145

whole process (O’Mard 1993) by allowing only 5—,6—, or 7—membered rings in the
structure. The annealing algorithm was performed on parallel computer system.
The optimal temperature schedule was found to be an inverse logarithmic function
(O’Mard 1993). The author reports that the parallel algorithm was run using ten
processors, and the 512-atom model took 140 hours of real time. The author
argues that his parallel algorithm with the temperature schedule and the ring—
type restriction, works better than the original Wooten—Weaire method. Indeed,
O’Mard obtained rms bond—angle deviation of 11.65° in his model, which is lower
than 12.34° - the value he obtained using the original Wooten-Weaire method.
Our objection to the claims that O’Mard’s algorithm is qualitatively bet-
ter than the original Wooten—Weaire method is based on our own results obtained
in modeling amorphous silicon and amorphous diamond. We have used the origi-
nal Wooten-Weaire method in constructing 216-, 512—, and 4096—model of a —- Si.
In the next Chapter it is shown that our 512 model of a - Si has even smaller
rms bond-angle deviation than O’Mard’s model. If 4—membered rings are ex-
cluded our model has rms bond—angle deviation of 11.06°, and with 4—membered
rings allowed we got even smaller value of 10.54°. Furthermore, the real time that
was needed to create our 512 models was certainly less that 140 hours reported
in O’Mard’s simulations, although due to the shake and bake nature of Wooten~
Weaire method we do not have a precise number for the elapsed time in creation of
our models which were all produced on ALPHA-DEC workstation. So, it seems
that the small numerical improvement O’Mard achieved in his simulations can
not be attributed exclusively to the temperature schedule used, and to the re-
striction of ring types. It still remains to be determined if the ’improvements’ in

the O’Mard algorithm really improve the Wooten-Weaire algorithm.

146
6.3 Conclusions

Early attempts to build a model of an amorphous structure by hand, even though
clumsy, laid the path for the implementation of computer algorithms which led
to the creation of very successful models. The idea of a microcrystallite nature
of amorphous materials was abandoned in favor of continuous random network
concept. Periodic boundary conditions were included, and increasing awareness
of the importance of topology in amorphous materials led to the invention of the
Wooten—Weaire topological bond switch, which, implemented into a Metropolis
simulated annealing algorithm, proved to be a powerful way to create a topo-
logically realistic amorphous structure using a computer. Designed originally
for tetrahedrally bonded amorphous elemental compounds, the Wooten—Weaire
method is adaptable to treat other systems as well, like vitreous silica, or amor-
phous mixture of diamond and graphite—like carbon atoms. The continuous ran-
dom network model which well describes covalent glasses should be understood
as an idealisation of realamorphous structures, in the same sense as a perfect
crystal is an idealisation of real crystals one can find in nature. CRN models can
serve as reference ideal amorphous structures to describe more realistic amorphous

structures which usually contain voids and various kind of defects.

Chapter 7

Modeling Amorphous Diamond

7 .1 Introduction

As we have seen in Chapter 6, a great deal of work has been done on modelling
amorphous silicon and germanium. The structure is well understood and widely
accepted as being a continuous random tetrahedral network (Polk 1971, Steinhardt
et al. 1974). In the early seventies it was still not completely clear if continuous
random networks could be continued outwards indefinitely without introducing
dangling bonds in the structure, or large vacancies. Subsequently, it' was under-
stood that this was possible by allowing slight variations in bond lengths and
bond angles which accommodate all the connectivity strains. Ziman, however, ,
wrote in 1979, that ”carbon bonds would not allow these strains so ideal ’glassy
diamond’ does not seem to exist” (Ziman 1979). In this Chapter we shall show
the opposite, presenting a good CRN model of amorphous tetrahedrally bonded
carbon. There has been growing interest recently in tetrahedral amorphous dia-
mond (Drabold et al. 1994), largely because McKenzie and coworkers (McKenzie
et al. 1991), and others (Lossy et al. 1992, Gaskell et al. 1992) have successfully

147

148

grown an amorphous form of carbon, which is predominantly (85%-90%) fourfold
coordinated.

In this Chapter we present our computer-generated continuous network
model of tetrahedrally bonded amorphous diamond (Djordjevic' et al. 1995‘), con-
structed by a method introduced by Wooten, Weaire and Winer (Wooten et al.
19885, Wooten and Weaire 1987) as described in the previous Chapter. We hope
this large model will serve as a useful reference point for discussions of the ob-
served structures of various forms of amorphous carbon, which in reality (to date)
always contain varying amounts of graphitic bonding. Our periodic supercell cell
consists of 4096 = 8 x 83 atoms. For the sake of comparison, we have constructed
similar models for amorphous silicon. Furthermore, we make two different struc-
tures for both cases; one where four-membered rings are allowed, and, the other
where four-membered rings are forbidden.

The layout of this Chapter is as follows. First, we shall give a short
overview of relevant properties of amorphous diamond. We than present and
compare results for all four models, and with the experimental measurements.

Discussion and conclusions follow.

7 .2 Wooten-Weaire Method

The method is described in detail in the previous Chapter and we give here just a
brief summary in which we refer to our amorphous diamond model (Djordjevic' et
al. 1995). The starting structure is a perfect diamond crystal. Periodic boundary
conditions are imposed so that the structure is fully tetrahedrally coordinated, i.e.,
there are no dangling bonds anywhere. The bond-bending and bond-stretching

forces are described by a Keating potential (Martin 1970, Keating 1966) as given

149

Table 7.1: The values of the parameters used in the models of C and Si (Djord-
jevié et al. 1995). These were obtained from measurements of phonon dispersion
relations and elastic constants in crystalline materials (Martin 1970).

 

 

 

 

Silicon Carbon
0 (dyn/cm) 0.485x103 1.293x103
fl/a 0.278 0.655
r0 (A) 2.35 1.535
L (A) 43.42 23.53

 

by the Equation (6.1) and rewritten here

3 a 2 2 36 1 2 2
“mymwol +§¥,§§,(r~'m+§%l ”-1)

where a and [i are diamond bond-stretching and bond-bending force constants
respectively, and r0 is the strain-free equilibrium bond length in the diamond
structure. There are 4096 atoms placed in a box of size L. The numerical values
of all the parameters used in our work are given in Table 7.1 (see Martin 1970).
The first step after generating the initial crystal structure is to randomize
the network sufficiently, so that subsequent annealing will not lead the system
back to the crystalline state. The randomization process is realized by the in-
troduction of many bond rearrangements (bond switches) (Wooten and Weaire
1987) at a temperature just above the melting point for the model. The melted
structure had A9 z 22° for the rms bond angle deviation. A useful criterion as to
whether the system is melted or not, is the structure factor I (k) associated with
the (111) direction for the diamond cubic structure. (We use here k to denote
wave vector, instead of Q which was used in Chapter 5.) If this quantity is of

roughly comparable intensity for all k values, then the structure is well randomized

150

(Wooten and Weaire 1987), and all memory of the original crystalline structure
is gone. The bond rearrangements maintain four-fold coordination at each atom
thus introducing large strains in the structure.

In the second stage of the procedure, the temperature is reduced in small
steps, and at eachnew temperature thermal equilibrium is established. The sys-
tem is relaxed geometrically (the release of the strain energy allowing the stretch-
ing and bending of bonds, according to Keating potential (Keating 1966)), and
topologically (creating more bond switches in the system). As" a consequence, the
total strain energy, rms bond-angle deviation and rms bond-length deviation are
decreased until an optimized amorphous structure is created as the temperature
approaches zero. A piece of the obtained structure is shown in Figure 7.1.

The presence of four-membered rings in the structure makes a success-
ful annealing process easier, while the absence of four-membered rings increases
the probability that the system will be trapped in an unsatisfactory metastable
state. When four-membered rings are allowed, the final structure contains only
a tiny fraction of four-membered rings, even though there was a large number of
such rings at the beginning of the process (in the melted structure after the first

randomization).

7 .3 Discussion of Amorphous Carbon

Since McKenzie and coworkers (McKenzie et al. 1991) have obtained 85% — 90%
four-fold coordinated amorphous carbon by plasma-arc deposition, there has been
an increasing interest in amorphous diamond-like carbon. These diamond-like
carbon films are found to be hard, optically transparent, and chemically inert,

which makes them potentially important for applications in coating technology

151

 

   

   

- ‘3 '6' "‘9 ‘51
' w A s
'l ' ' toil—5 e.‘

      
  
  
 

  
   

 

‘1— . " 1. ‘fiaéi ' .
‘ \ ’ ‘ " ‘E ‘
”9 “gigggflf’w‘gfii

 

Figure 7.1: A piece of computer-generated tetrahedral amorphous diamond struc-
ture (Djordjevié et al. 1995).

152

and also for use as wide band gap semiconductors. It is known that amorphous
silicon and amorphous germanium exhibit a striking similarity in the shape of
their radial distribution functions (see Etherington et al. 1982 and Wooten et al.
1985). In both of these materials the bond streching forces are dominant compared
to the weaker bond-bending forces. Amorphous carbon is different with respect
to both a-Si and a-Ge because of very strong chemical bonding, and because the
bond-bending forces are much larger relative to the bond-streching forces. The
ratio of bond-bending to bond-stretching forces fl/a is 2.4 times larger for carbon
than for silicon (see Table 7.1). One might have expected a different amorphous
structure for carbon (i.e. different ring statistics etc.) compared with silicon and
germanium, because of the greater bond stiffness of carbon. This was the initial
motivation for this work, but our results show that the Pair Distribution Function
(PDF) of tetrahedral amorphous diamond is remarkably similar to the PDF’s of

amorphous silicon and amorphous germanium.

7 .4 Results and Comparisons

We constructed four different models: two for amorphous carbon, and two for
amorphous silicon. We were interested in the ring statistics, and whether the
presence of the four-membered rings makes any significant difference in the fi-
nal structure. So, one amorphous carbon (silicon) model was constructed with
four-membered rings present and the other without four membered rings. We
calculated the total strain energy, rms-bond length, rms bond-length deviation,
rms angle deviation, and the mean angle between the bonds. Also, we calculated
the number of 4—,5-,6-, and 8-membered irreducible rings present. The definition

of a irreducible ring is given in Chapter 6, together with the discussion of their

153

Table 7.2: The structural parameters for the models discussed in this Chapter
(Djordjevié et al. 1995), at two densities - the crysalline density (p/po = 1.000)
and the optimal density. The number of irreducible rings per site n, is independent
of the density as the topology is not changed when the density is varied, and it
drops to zero for large r..

 

 

 

 

Si(4) Si C(4) C
p/po 1.000 1.025 1.000 1.033 1.000 1.050 1.000 1.054
<r>/r0 0.996 0.988 0.997 0.987 0.996 0.980 0.997 0.980
<9>° 109.24 109.23 109.25 109.23 109.31 109.30 109.31 109.30
Ar/ro(%) 2.52 2.46 2.65 2.57 4.39 4.31 4.40 4.31
219° 10.51 10.54 11.02 11.06 9.29 9.31 9.06 9.08
E(eV)/atom 0.336 0.334 0.367 0.362 0.757 0.741 0.731 0.715
n4 0.015 0 0.009 0
725 0.491 0.523 0.472 0.457
n, 0.698 0.676 0.760 0.800
127 0.484 0.462 0.472 0.498
n8 0.156 0.164 0.156 0.143

 

 

 

 

 

importance for topological purposes All results are presented in Table 7.2.

We tested the influence of the sample density on the strain energy, and -
found that the crystalline density is not the optimal one for the amorphous struc-
ture. In all cases, a slightly larger density gave a smaller strain energy when using
the Keating potential. The Stillinger-Weber potential (Stillinger and Weber 1985)
goes in the opposite direction for silicon: slightly smaller density gave the small-
est strain energy (Wooten and Weaire 1987). The small density differences should

therefore not be taken seriously, although it is proper to optimise the density to

154

give the smallest strain energy within a given model. The optimal density was
determined by fitting a parabola to a plot of strain energy as a function of density,
and finding the minimum, given in Table 7.2.

From the coordinates of the final amorphous structure obtained for the

optimal density, we calculated the PDF, which is given as:

P(r) — G(r) (7.2)

— 47rr2p0

 

where G(r) is the radial distribution function defined in Chapter 5, and p0 is the

number density of the system:

N
P0 = E (7-3)

where N = 4096, and L is the actual size of the supercell (either the original
supercell size as given in Table 7.1, or the optimized one). For comparison, we
rescale the PDF of silicon to that of carbon, with the distance scaled by r0, the
nearest neighbor distance in the crystalline form. The four PDF’s are shown in
Figure 7.2.

Using a Fourier transformation, the structure factor can be written (see

equation (5.22) as

1(k)=1+4‘7”’0

 

f0” r[P(r) — 1] sin(kr)dr (7.4)

where r0 is the strain free nearest neighbor distance for the crystalline structure
as given in Table 7.1. We calculate the structure factor for the four models which
are shown in Figure 7.3. Because there is no contribution to the integral (5) when
P(r) approaches unity, it can be seen from Figure 7.2 that this integral should

(and does) converge well, with only a few termination ripples at small 11:.

 

 

 

 

 

 

 

).I I r I I I I I I I I I I I I I I I I I I L I I I I I I I I I [I I r I I I I I I I I1.

P d d

I- u)- -(

- dP d

h d. d

D d- II

- q- .1

P I- u:

‘ — —_ —

V -II- —

b -h -l

- -- q

)- -1- -1

8 — __ —

_ 0 _ 1 _.

.. W1
1- b i

q

o 1 l l 1 1 11 l 1 1 l

,_ , I I I I I I I I I I I L

 

 

 

 

U)
H‘

 

 

Pair distribution function

0
'IIIIITIITIIIIIIIII

 

111111111111111

IIIjIITIIIIIIIIII

  

if?

  

 

 

 

 

 

 

 

111111111111'11

I I 4 0
Distance r/ro

 

l
M1IIILLIIIIIL 11

8 3 4

O
in

Figure 7.2: Pair distribution functions P(r) for amorphous carbon and amorphous
silicon models, with, and without four-fold rings (Djordjevié et al. 1995). All
models are at the optimized density.

156

 

  
   

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

IIIIIIIIITIII IIII IIII f III IIII IIIr IIII IIII FF
2.0 l l I 1 n1 1 l l
’ C(4) fl
1.5 ’ I
1.5 h i V ’H
H . - . ..
3 E U ’ Iii ’ \/ 2
O 0.5' ’ fl U -:
.9. E : f :
1 - , -
2 0.0 5““ “be -’
I 1 ..
3 11.5 . g , _‘
o : ‘1 . I ‘1 . I
a , ,. S1(4): .1 SI -
...a 1.5 ' l ” J
m I :1 I 1’ q
l 1' I ." i
1.0' ‘ " - A ' ’
| , -
0.5I . ‘-
| -
g i
0.0! lllllllllllllllllllllllll 11L1111L11111|11L1l1111111.
o 5 1o 15 so 0 5 1o 15 so 25

Wavevector kro

Figure 7.3: Structure Factors I (k) caculated from the corresponding pair distri-
bution functions shown in Figure 7.2 (Djordjevié et al. 1995).

157
7 .5 Discussion

The most striking feature of the results in Table 7.2 is that there are no major
differences between amorphous carbon and amorphous silicon. We considered two
densities, the crystalline density, and the density which minimizes the total strain
energy. The rms bond lengths, < r > /ro, of carbon and silicon are about the
same, and slightly smaller than in the crystalline diamond structure. The rms
bond-length deviation, Ar/ro, is larger for carbon by 44%. This is the one of
the most significant differences between carbon and silicon. The rms bond-angle
deviation, A6°, is smaller for carbon by 18%, as would be expected because of
the larger angular force reflected through the larger value of fl/a. Ring statistics
are also very similar, except that there are more six-membered irreducible rings
(Wooten and Weaire 1987) in carbon, and fewer five-membered rings than in
silicon.

The optimal density results given in the second sub-columns in Table 7.2,
show the following features. The energy is minimized, in all cases, at a slightly
higher density than the crystalline one. The rms bond length, the mean bond
angle, and the rms bond-length, are slightly decreased, while the rms bond-angle
deviation is slightly increased.

Looking at Figure 7.2, one can see a surprisingly good matching between
the PDF’s of carbon and silicon. Their general shape is the same. The first peak
is shifted towards to a value slightly smaller than unity; the crystalline value.
This is true even at the crystalline density, and corresponds to the rms relative
bond-length, which is smaller than one in the crystalline state, i.e. the mean
bond length in the amorphous structure is slightly smaller than in the crystal

(1.544A). The carbon first peak is broader than the silicon first peak, which is

158

correlated with the larger rms bond-length deviation in carbon. This represents
a different aportioning of the strain energy between bond length strain and bond
angle strain in carbon and silicon. The ratio of the strain energy associated with
the bond angle (second term in the Keating potential) and the bond length (first
term in the Keating potential) for carbon is found to be 1.57, while for silicon it
is 2.73. This means that relatively less energy is stored as bond bending strain in
carbon, when compared with silicon. Nevertheless the overall features of PDF’s
are essentially unchanged between carbon and silicon.

Figure 7.3 shows again great similarity in the structure factor for carbon
and silicon. Furthermore, it fits relatively well the measured structure factor
for amorphous, diamond-like carbon, as shown in Figure 7.4 (see Gilkes et al.
1995, Figure 1). It should be remembered that these experiments are on samples
containing 10% - 15% graphitic bonding so that a detailed comparison between
our model and experiment is not appropriate at this time. Graphitic regions
can be incorporated into computer generated models using molecular dynamics
and other techniques (Wang and Ho 1993). However these models are much
smaller than ours [the model in (Wang and Ho 1993) contains 216 atoms] so that
a detailed comparison of these models with experimental structure factors is also
premature. It must also be realised that carbon films (Gilkes et al. 1995) contain
some hydrogen which will have some effect on the PDF and structure factor.

Our very large (4096 atom) and realistic structural model of amorphous
diamond is potentially very appealing for addressing some of the other important
issues of the physics of amorphous solids. One of those is the nature of the
band tails in the electronic density of states in amorphous materials (Mott and

Davis 1979). Recently, Dong and Drabold (1996) used our structural model of

159

 

 

 

 

 

 

 

 

 

2.5[r
f
I
.l.
1
l
2.01-
i
f... 4
O 1
+4 I
g '1'
.H 1.5:-
Q) 4
i-I 1
5 1
"63
1.0
5 .
I-o
.1.) '1
U) .1
A
0.51-
.1
f f
. I1 .-
0.01111111111111-111111

50 100 150 200

k (nm'l)

Figure 7.4: A comparison of the structure factor from our computer-generated
model of amorphous tetrahedral diamond with four-fold rings (dashed line) (Djord-
jevic’ et al. 1995), and the structure factor measured by Gilkes ct al. (1995) (solid
line), on samples containing 10% - 15% grahitic regions.

160

amorphous diamond to study the nature of the electronic localization in three
dimensions in the presence of topological disorder. The crystalline diamond cell
has a defect—free gap between the valence and the conduction bands, with sharp
band egde. Dong and Drabold (1996) found that our amorphous diamond model
has extended band tailing at both valence and conduction bands, as well as few
defect states in the middle of the band gap. The defect states which are highly
localized, are caused by the atoms with the few largest bond angle distortions.
The band tail states change from highly localized states that surround some defect
cores, to less localized states that are within larger clusters, finally to the extended
states that are distributed throughout the cell, when the energy moves from mid—
gap into the valence states. The authors present a very simple picture that suggest
that any mechanism causing approximately Gaussian bond angle disorder must

lead to the exponential tails (Dong and Drabold 1996).

7 .6 Conclusion

Although carbon has much stiffer bond-bending forces than silicon (i.e. a much
larger value of B/a, the ratio of bond-bending to bond-stretching force constants
in the Keating model), this doesn’t influence the pair distribution function signif-
icantly (Djordjevié et al. 1995). The bond-angle deviations are slightly smaller .
in amorphous tetrahedral diamond compared with amorphous silicon and amor-
phous germanium, but they cannot be too small, as the atoms are constrained to
be in a fully tetrahedrally coordinated random network. Indeed our experience in
constructing random networks leads us to postulate that a rms bond-angle devi-
ation of at least 7° is required if the network is not to contain crystalline regions

and so be a two phase mixture. Our very large and realistic structural model

161

of amorphous diamond also provides a good tool for studying other important
issues of the physics of amorphous materials, for instance, the nature of electronic

localization in three dimensions in the presence of topological disorder.

Chapter 8

Elastic Properties of Random

Networks

8.1 Introduction

As we have seen in Chapters 6 and 7, continuous random network (CRN) models,
first introduced by Zachariasen (1932), are very successful in describing networks
of covalent glasses. Covalent glasses consist of atoms that make from two to four
covalent bonds, generally with specified bond lengths and bond angles. A network
glass can be regarded as a mechanical structure and one can analyse its mechanical
and elastic properties in terms of the network stability. This approach dates back
to Maxwell’s work in the mid—nineteenth century (Maxwell 1864). It has been

applied to glasses by Phillips (1979, 1980, 1981).

162

163

8.2 Rigid and Floppy Networks

hilaxwell discussed the stability of a hinged connected frameworks in terms of the
relation between the number of constraints and the number of degrees of freedom.
If the simple network consists of Nb bars connected at N hinged joints in 2d,
the number of constraints NC is equal to the number of bars Nb. The number of
internal degrees of freedom is 2N—3, where 2 rigid rotations and l rigid translation
are subtracted as being not internal degrees of freedom. Maxwell concluded that

the stability of 2d network is described by the following inequalities

Nb > 2N—3, stable,

Nb < 2N — 3, unstable. (8.1)
In 3d stability condition (8.1) becomes

Nb > 3N—6, stable,

Nb < 3N—6, unstable. (8.2)

These two conditions describe so called overconstrained (rigid), and undercon-
strained (floppy) networks respectively. For large networks, one can neglect the
rigid translations and rotations, and in 3d, the conditions which separates rigid

from floppy networks become

N. > 3N, rigid,

Nb < 3N, floppy. (8.3)

The number of constraints is equal to the number of bars only approximately, and
Maxwell’s condition of stability (8.2) should be understood, in the spirit of mean
field theory, as an average global condition - which is, however, very good one

(Thorpe 1983, 1985).

164

Another way to look at Maxwell’s result is to replace rigid bars by Hooke
springs and to calculate the vibrational frequencies of the network. T his is the
standard problem in classical mechanics. The eigenfrequencies are calculated by
diagonalising the dynamical matrix. If the network is floppy, there will be certain
number F of eigenfrequencies that are zero. These modes are called floppy modes

(Thorpe 1995). The number of floppy modes F is given in 3d by
F = 3N — NC, (8.4)

where 3N is the number of degrees of freedom, and NC is the number of linearly

independent constraints. Recalling equation (8.3) we can conclude that
if F > 0, network is floppy. (8.5)

If the network is _overconstrained, i.e. if constraints become more and more depen-
dent, equation (8.4) would give negative F - the quantity which is by definition

non—negative. In practice, good results can be obtained assuming that
F = max[3N — N00]. (8.6)

In covalent network glass, (represented by a CRN model), the interatomic
forces are well described by a Keating type of potential (7.1), which contains bond—
stretching part (proportional to a force constant), and bond—bending part (,8 force
constant). Thus, constraints imposed on the network are those associated with
the fixed bond lengths, and fixed bond angles. There is one a—force Eonstraint per
bond and 27' — 3 independent fl-force angular constraints around an atom with
r 2 2 bonds. If there are n, atoms in the network with r bonds (r = 2, 3, 4), then

the total number of constraints is given by

N. = in, [g + (2r — 3)] . (8.7)

r=2

165

Noting that

4
IV = Z n,, (88)

r=2

the mean coordination of the network is defined as

_ 23:2 Tn?

_ 23:2 "’7' c

The fraction f E F/(3N) is given by

(7‘)

f = [323274. - 2n. [g + (2r — 31]] /(3N) (8.10)

r=2 ""

which can be rewritten in terms of mean coordination (r) (8.9) as

fzz—gg) win

Thus, from (8.11), the number of floppy modes f is zero, at the mean coordination

(r) = 2.4. The Maxwell conditions of stability (8.3) can be thus reformulated as

if (r) > 2.4 network is rigid,

if (r) < 2.4 network is floppy. (8.12)

This concept that network glasses can be divided into two different types was
first introduced by Phillips (1979,1981). It was further discussed by D6hler et al.
(1983) and subsequently made more rigorous by Thorpe (1983). For the covalent
glasses which do not have dangling bonds (r = 1), the mean coordination is in the -
range from 2 to 4. When (r) = 2, f = 1 / 3, which means that one third of all the
eigen—modes are floppy. This corresponds to the floppy network which consists of
long polymeric chains. As new bonds, i.e. cross links are added, f drops toward
zero when the network becomes rigid. In the region 2 < (r) < 2.4, there are rigid
and floppy regions in the structure, but rigid regions do not percolate. For (r) >

2.4 the network becomes amorphous solid in which rigid regions percolate. This

166

result can be corrected for dangling bonds (Thorpe 1995). Thus, at the critical
mean coordination, a rigidity phase transition occurs. An extensive research has
been recently done by Thorpe at al. (1996), Jacobs and Thorpe (1995, 1996) and
Jacobs and Hendrickson (1996), to determine the universality class of this rigidity

percolation phase transition by calculating the critical exponents.

8.3 Elastic Constants of Random Networks

The ideas from preceeding section have been tested numerically by He and Thorpe
(1985). They calculated both the number of floppy modes f, and three elastic
moduli, (711,044 and bulk modulus B, for random networks with varying mean
coordination. These physical quantities couple most directly to the rigidity phase
transition.

To build large number of random network He and Thorpe used bond—
dilution of the crystalline diamond lattice. They start from 512—atoms crystalline
diamond supercell that has (r) = 4. The atoms interact via the Keating potential
(7.1) with ,B/a = 0.2. Then, bonds are removed in random manner, maintaining
the structure with periodic boundary conditions. When a bond is removed all the
a and 6 terms associated with it are removed from the Keating potential. No
dangling bonds were permitted, so an arbitrary atom can have 2, 3 or 4 bonds.
After certain number of bonds is removed, the structure is relaxed using the
Keating potential. The initial crystalline structure, which serves here as a kind
of a reservoir of bonds, is a stress—free network. After some bonds are removed, a
small stress is introduced, but subsequent relaxation releases this stress and the
system is again at zero strain energy. In this way a sequence of random networks

was obtained in the whole range of allowed mean coordination. For each of those

167

mean coordinations there were three samples with different statistics, generated by
using three different random number generators. Also, one of these three samples
had an enhanced number of 2~fold coordinated atoms at the expense of 3—fold
sites. This served to test the general prediction that the properties of the network
depend mainly on (r) and not on other parameters. All results were averaged over
these three statistical configurations.

The computed f drops linearly as predicted by the mean field theory (8.11)
but it has a small tail around the mean field transition point. Elastic constants
C11 and C44 (averaged over the 2:, y and 2 directions), and the bulk modulus
B = %(Cu + 2012), are computed with use of standard techniques (Feng et al.
1985; Feng and Sen 1984). The 512—atom supercell is redefined by an external
strain 6 and the elastic modulus C is obtained from the elastic energy %Cc2 after
the network has been fully relaxed. The result is shown in Figure 8.1. He and
Thorpe made a least—squares fit to the computed data and found that, within the
range 2.4 < (r) < 3.2 all three elastic constants are very well approximated by

the power low of the form (see Figure 8.1)
C = const.((r) — 2.4)1'5. (8.13)

The equation (8.13) does not fit the data well at larger values of (r). In earlier
studies, when the simple central—force models were used, elastic constants were
linear in (r) in the whole range, with the exception of a small tail near critical
region (Feng et al. 1985).

Recently, Franzblau and Tersoff (1992) performed computer simulations
similar to those of He and Thorpe (1985) to calculate the elastic constants of
random networks. They employed a direct algebraic approach in calculating elastic

constants (Franzblau and Tersoff 1992). Their model has 216 atoms, but all results

168

 

 

tsF r F T I q
H311 °
o—oc“ a
H.
1.01)- o
i .
Q
.-
gut

 

 

 

 

Figure 8.1: Elastic moduli averaged over three different samples as a function of
the mean coordination (r) for fi/a = 0.2. Solid lines are power low fit curves. (He
and Thorpe 1985).

169

are obtained in the wide range of the force—constants ratio I‘C E 1’3/0 over 4 orders
of magnitude. The authors also use the same method of bond—depletion like He
and Thorpe (1985), and the underlying structure is again crystalline diamond
lattice. The main graphical results are given for FC = 0.655 which is the value for
diamond lattice. This is significantly larger value than one used in work of He and
Thorpe which is I} = 0.2. Franzblau and Tersof found that for this particular,
diamond value of I}, all three elastic constants they computed (Cu, C44 and shear
constant C, E (C11 - C12)/2) follow very well the following power low in the whole
range of (r)

C = const.((r) — 2.394)1‘40. (8.14)
The elastic constant C44 behaves very accurately according to power low in the

whole range of parameter I‘c
IOFC < P < I‘C/IOOO. (8.15)

The exponent in the power low (8.14) vary from 1.35 to 1.89 within the range of
Fe.

The other elastic constants exhibit more complex behavior. For either
very small or very large Fe, C11 deviates considerably from the power low. Also,
C 3 deviates for large I}, but it obeys the power low closely for small I}.

Note that results of Franzblau and Tersoff (1992) are partially consistent
with earlier results of He and Thorpe (1985) which were done for small value of
PC = 0.2. According to Franzblau and Tersoff, Cu deviates significantly from the
power low for small I}, and so it is not surprising that He and Thorpe (1985)
were not able to fit C11 and B constants in the whole range of (r) with the power
low. But, it remains unclear why it was not possible to fit the elastic constant

C44 with the power low in simulations of He and Thorpe (1985), when according

170

to Franzblau and Tersoff (1992) that was possible in their simulations within the

wide range of Fe.

8.4 Bond—dilution of Amorphous Diamond Net-
work

All existing computer simulations of elastic properties of random networks are
performed by using a crystalline lattice as the underlying structure. In other
words, the crystalline lattice is transformed into a random network by removing
certain number of bonds at random, and subsequently relaxing the diluted struc-
ture. As we have stressed in Chapter 6, one has to be very careful in creating
models of amorphous materials to avoid any crystalline memory in the final ran-
dom network. Thus, one can be suspicious about the true amorphousness of the
random network obtained by the random removal of bonds from an ideal crys-
talline structure. Conceptually, it is more appropriate to start from an amorphous
fully connected network, and then perform bond-dilution on it, thus creating a
sequence of random networks with different mean coordinations. We use our
512-atoms CRN model of amorphous diamond as this more realistic underlying
structure. Repeating the basic guidelines presented in the work of He and Thorpe
(1985), we calculate the bulk modulus B E (Cu + 2012)/3 as a function of the
mean coordination (r) (Djordjevié and Thorpe 1996).

Our amorphous diamond network is fully 4—fold coordinated and has a
large stress stored in bond—stretching and bond—bending parts of Keating po-
tential (7.1). This is the main difference with respect to earlier work where the

starting structure was a stress—free crystalline lattice. In the present case, removal

171

of certain number of bonds releases significant amount of elastic energy, and after
subsequent relaxation the structure still has a significant amount of stress. Thus,
while the strain energy is always zero when the crystalline lattice is diluted. here
the strain energy drops from the value it has in our amorphous diamond model,
to zero near the transition point (r) = 2.4. In fact, we do not expect that the
strain energy in particular, must become exactly zero at the mean field critical
mean coordination, because at the critical point there are still many isolated rigid
regions which do not percolate, but contain a certain amount of strain energy.
It is a different question how those isolated rigid regions effect the behavior of
the bulk modulus B. Being isolated, it may happen that they do not make any
impact on the macroscopic stiffness of the structure which is described by the bulk
modulus. However, to our knowledge, this question remains inconclusive at this
stage of understanding of rigidity percolation phase transition.

We designed programs that remove bonds from the amorphous diamond
network in several different ways. The obvious one does not impose any restrictii :11
except that it does not allow creation of dangling bonds, thus fixing the possible
atomic coordination to 2, 3, or 4. As a consequence, when a large number of bonds
is removed, i.e. when (r) -i 2, there are many long linear polymeric chains in the
structure. Of course, if one wants to decrease the mean coordination to 2, linear
chains must form. But if we are mainly interested in the region 2.3 < (r) < 4.0, '
than we can construct various random networks with the same (r). Thus, at
(r) = 2.59 we have made four different structures. The first is already mentioned
— it has long linear polymeric chains. The second has only short linear chains
that contain not more than three 2-fold coordinated atoms. The third one has

only two—atom linear chains, and in the fourth one all 2—fold coordinated atoms

172

are isolated from each other. It was possible to construct random networks with
different type of chains even for lower values of (r), but we did not explore this
in more details. For our calculation of bulk modulus we used the first, non—
restrictive style of bond removal, so the polymeric chains are present in all our
random structures at low (r).

The initial structure is fully connected amorphous diamond network with
the mean coordination (r) = 4.0. To avoid any possible correaltions between
our random networks in a sequence of different (r), we adopted the following
strategy. First, we remove small number of. bonds, say 20, relax the lattiice using
the Keating potential, and calculate the mean coordination (r), and the total
strain energy E. We repeat that process 10 times using different random number
generators, and so produce 10 samples for a given (r), over which we perform
the average. To further decrease the mean coordination, we do not continue from
the structure that has been previously obtained, but rather we go back to our
fully connected amorphous diamond network, and remove now 40 bonds, relax
the diluted system an calculate new, smaller, value of (r) on 10 different samples.
In the next step, we go back to the initial structure and remove 20 more bonds
(i.e. 60 bonds), etc. . . . Each time we increase the number of bonds we remove,
the greater stress is released from the initial amorphous diamond structure and

the total elastic energy drops down.

8.5 Strain Energy and Bulk Modulus

It has been already reported (He and Thorpe 1985) that the relaxation process
becomes very difficult for approximately (r) < 2.6, even when the crystalline

stress—free lattice is used for bond—dilution. In that range of mean coordinations,

173

He and Thorpe used an extrapolation technique to obtain the behaviour at longer
times (i.e. larger number of steps in the relaxation procedure) than could be
probed. In our calculation it is even more difficult to perform the relaxation at
small values of (r), but we did not need to use any extrapolation techniques for
mean coordination as low as 2.20, which was small enough for our purposes. In
Figure 8.2 we show the total strain energy as a function of the mean coordination
(r). For large mean coordinations the relaxation is very fast, and only 500 relax-
ation steps are needed. For small mean coordinations the number of relaxation
steps has been considerably increased to achieve the convergency. The required
convergency has been kept constant, and the number of relaxation steps was vari-
able so that, at small mean coordinations, it increased to 50, 000 and more.

The bulk modulus B has been calculated in a standard way by redefin-
ing the size of the supercell uniformly along all three coordinate axes. For each
mean coordination, and for each of the ten configurations, we have included the
supercell size as a variable in the relaxation subroutine. In this way the structure
is relaxed not only with respect to the 01— and ,B—forces in the Keating potential,
but simultaneously also with respect to the size of the supercell. In this way we
obtained the reference point for bulk modulus calculations. Starting from that
optimal supercell size, we redefined the supercell four times — two times increasing
the supercell size sequentially, and two times decreasing the supercell size sequen-
tially with respect to the original optimal size which minimizes the total elastic
energy. By sequentially, we mean that we take coordinates of the atoms when
the supercell is of optimal size, and apply a small uniform strain which slightly
compresses the lattice. Then we relax the structure with Keating potential, cal-

culate the new coordinates of all atoms, and apply the same strain again, using

174

 

 

 

 

 

08 r I l l l 1 I 1 I
0-7 - E/512 = 0.4093*(<r>-2.4)1°15
N 0.6-
,_
l0
\
“'1 0.5-
E
.9
‘3 0.4..
0)
G.
>
g’ 0.3-
C
a) ..
c .
'5 0,2.
.1:
U) a
0.1.
0.0- W
2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0

<r>

Figure 8.2: Dots represent the total strain energy averaged over 10 different con-
figurations, as a function of the mean coordination (r). The solid line is the best
fit power-law curve. The energy is given in eV per atom (Djordjevié and Thorpe
1996).

175

the updated coordinates from the previous configuration. In this way we obtain
the values of the total strain energy at five different supercell sizes. These points
should, in principle, form a parabola, which, of course, is the reflection of the
harmonic nature of the Keating potential. The bulk modulus B is, up to the
constant prefactor, equal to the curvature of that parabola. The problem that
arises when diluting our amorphous underlying structure is common for all amor-
phous materials - one can spontaneously and unpredictably jump from one to
another metastable state with different energy. If that happens in the program,
one has to reject that situation, and to repeat the search for parabolic energy
dependance again. This occurs very often and makes the computation painfully
slow, if one does not first find the minimum energy with respect to the supercell
size (Mousseau 1996 (b)). This is why we implemented simultaneous minimiza-
tion of energy with respect to the supercell size. When that minimum is used
to redefine the supercell, it is very rare that system jumps from one to another
metastable state. In Figure 8.3 we show a typical jump in energy which occurs in
amorphou.l networks. With careful examination we were able to reconstruct the
’missing’ parts of the two parabolas, thus obtaining two different minima which
correspond to two metastable states of the system, as shown in Figure 8.4.

We have investigated the nature of the two states that appear in our bond—
depleted amorphous structure in more details. It became clear that the value of
the energy minimum does not depend only on the supercell size, but also on the
particular history of changes that were imposed upon the system. We were able
to identify the atoms whose motion is responsible for the existence of these two
energy states. Comparing the atomic coordinates of the structure in these states

we made a picture which shows 512 atoms in one of the states while black bars

176

 

 

 

0.0546,,V,,,.,,,v,.,.,,,T
0.0544 4 «”259 q
C
T 1
0.0542. s
d 4
0.05404 , a
(TV- 1 \ . 1
l0 .\
\ 0.0538- .\ / / -
Lu
. O\ /.
0.0536. ‘\.\ /o q
4 Ko‘, /0 ‘
0.0534. \ ./' _
”yr ‘
0.0532

 

13.80113f85'13:90'13f95'fii00'14:05'14i10'14l15'14i20'14r25 ' 14.30
L (angstroms)

Figure 8.3: The system is systematically redefined by increasing the size of the
supercell (arrows). The energy values all follow a parabolic curve, until a jump
to a lower energy state occurs, and the second parabolic curve is created. Energy
is in eV per atom, and mean coordination is (r) = 2.59 (Djordjevic’ and Thorpe

1996).

177

 

0.0486 . <r> = 2.59 /0 q
4 . -
0.0434 . ./ .
A 00432 - , ./'/ ‘
> .
3 .\o /’/
N k.‘ 0/ /
:5 1 . WV Q ,
\ 0.0473 . \ / .
I.“ O
1 \.\ /
0.0476 4 .\ o/ _
4 '\ ./'/ ,
0.0474 4 My“ .

 

 

 

13.35 ' 13T90 ' 13f95 -14.“, ' 14'05 ' 14'10 ' 1435 ' 14320 ' 14T25 ' 14.30
L (angstroms)

Figure 8.4: The picture shows two different energy minima which are characteristic
for amorphous materials. Energy is in eV per atom, and mean coordination

(r) = 2.59 (Djordjevié and Thorpe 1996).

178

that are drown represent the displacement of the given atom to another position
which it has when the structure is in other energy state. Thus relatively small
fraction of atoms is involved in repositioning which gives rise to the existence of
a second state. This is shown in Figure 8.5.

Using the neighbor list for the structure we idendified the connectivity of
the atoms involved in that motion. It turned out that all atoms there moved the
most belong to the linear polymeric chains, as shown in Figure 8.6.

We spent some time on searching for two states of particularly small dif-
ference in energy, hoping to find so called tunneling modes that exist in amorphous
systems. A common characteristic of all amorphous materials is a quasilinear term
in the specific heat at low temperatures. This has been successfully interpreted
phenomenologically in terms of a tunneling model, in which atoms, or groups of
atoms, tunnel between the two lowest—energy states in a double potential well
(Phillips 1972; Anderson et al. 1972; Phillips 1985). The tunneling is a very
subtle process which involves energy differences of about 10'4eV (Cusack 1987).
Detailed microscopic models for the tunneling process is still missing. Amorphous
elemental compounds like a - Si, a - Ge, or a — C, are of particular interest
because of their relatively simple structure. Although it once seemed that the
rigidity of the fully coordinated tetrahedral networks would prevent the possibil-
ity for tunneling states to exist, the experimental evidence is quite strong that '
tunneling states exist in these materials (Duquesne and Bellessa 1983). There
were attempts to create a tunneling model which would explain the anomalous
properties of vitreous silica (Vukcevich 1972), and an early computer simulation
of tunneling states in random network glasses for a — Si and a — Ge (Smith 1979).

In spite of all these attempts a satisfactory microscopic model for the tunneling

179

 

 

 

 

Figure 8.5: The picture shows 512 atoms in their positions which correspond to
one of the states described by parabolas in Figure 8.4. The black bars represent
the displacements of some atoms which give rise to the existence of the second
state. Mean coordination is (r) = 2.59 (Djordjevié and Thorpe 1996).

180

 

I

i O

l ///'7§"f3 /
/ W _ \ '

z'\ “\o‘ $13750th

MO—

/ .Q

 

\
O “fl
/ V/‘/'
1;;
{9' \’
nil“,
11")?“

W
"'4
lg
/~

\7/
L
\

g;

O
W]-
([917.
i
i

 

 

 

 

‘ \gf.‘ \
2% éO/‘é‘x‘: {Iii/1 !\ \o/[o‘
. §:- 'l/é“ jt\cfl'o \‘/ /§~

 

 

 

 

A} \'1 \‘-
fli /1/')!}:r

1v 2" ‘L..‘

Figure 8. 6: The picture shows polymeric chains of atoms that are involved 1n
the motion which gives rise to the existence of two states for the structure with
(r )= 2. 59 (Djordjevié and Thorpe 1996).

181

modes still does not exist. We were not successful in our search for tunneling
states in our random networks. From Figure 8.3 one can see that the typical
total energy difference between the two states we have in our structures, is about
512 x 0.0002 = 0.1eV which is three orders of magnitude larger energy than one
typically involved in the tunneling modes. Nevertheless, detailed analysis of the
energy minima in our random networks helped us to develop a program which eas-
ily finds the regions of parabolic energy dependence on the applied strain, which
was crucial for extracting the value of the bulk modulus B.

In Figure 8.7 we show our computed bulk modulus B as a function of
the mean coordination (7'). It is clear from Figure 8.7 that the underlying initial
structure from which a random network is obtained by bond dilution, does not
make significant difference in the behavior of the bulk modulus. Our amorphous
starting network, upon the bond—depletion produced the bulk modulus which
almost exactly follows the power law of Franzblau and Tersoff (1992). The small
discrepancy for the middle range of mean coordinations could be the consequence
of the bad statistics we had in our calculation having only 10 samples for averaging.
The other possibility is the fact that F ranzblau and Tersoff (1992) do not give the
explicit value of the force constants ratio I‘c apart from stating that this is the
diamond value. We have used PC = 0.655 which is the value we had in our
anorphous diamond model.

An interesting result from our simulations is the behavior of the optimal
supercell size ngn, i.e. that one which minimizes the total strain energy at given
(7‘). In Figure 8.8 we show that as more and more bonds are removed from
the initial amorphous diamond structure, the optimal supercell size first slightly

increases until about (1') = 2.7. Further reduction of the mean coordination is

182

 

 

 

 

 

 

1.2
1.04
A 0.8
e 'l
m 1
3
'3’ 0.6..
'U
0
E
41 0.44
3
CD
0.2-
00 ' I I ' I ' I ' I ' I I ' F '
2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0

<r>

Figure 8.7: Dots represent our result for the bulk modulus B as a function of the
mean coordination (r), for bond-depleted amorphous diamond lattice (Djordjevié
and Thorpe 1996). Each dot is obtained from the curvature of a parabola like
one of those in Figure 8.4, and averaged over 10 configurations. The solid line
represent the bulk modulus obtained by bond—depletion of the crystalline diamod
lattice (Franzblau and Tersoff 1992). FG = 0.655.

183

followed by a sudden decrease of the optimal supercell size.

We interpret this behavior as a collapse of the network due to the presence
of the long linear polymeric chains. The two-fold coordinated atoms in these
chains can come arbitrarily close to each other, because of the lack of an exclusion
volume that initially existed due to the tetrahedral coordination of the amorphous
network. To test this hypothesis we calculated the pair distribution function
for four different mean coordinations, (r) = 2.44,2.75,3.06,3.38 (Figure 8.9).
At small (r) = 2.44 there are many peaks close to the origin, before the first
and second peak that correspond to the nearest—neighbor, and the next nearest—
neighbor distance in the fully connected amorphous diamond, respectively. These
additional peaks that occur at very small distances, are the consequence of the
collapse of the network when the two—fold coordinated atoms in linear chains
come very close together. This looks like a phase transition which occurs much
before the rigidity percolation phase transition occurs at about (1‘) = 2.4. If we
artificially introduced an exclusion volume around each atom in the network, the
collapse would not occur. Collapse is triggered by the large stress that is present in
the intial amorphous diamond network. As we said, at small (r), large number of
bonds were first randomly removed from the fully connected amorphous diamond
network. Because of that, the stress in the network suddenly drops and the atoms
are free to move toward the equilibrium positions which occurs in the process of
relaxation. Because the exclusion volume of an atom is much smaller for small
mean coordination, those atoms that belong to linear chains can approach each
other very closely. In earlier work of He and Thorpe (1985) and of Franzblau and
Tersoff (1992), where the underlying structure for bond removal was a stress—free

crystalline diamond lattice, such a triggering mechanism was absent because the

184

 

 

 

 

14'4 r I fl I ' I ' I ' I ' I I I
j 4
14.2- ' o o .
‘ o O o 0.0.... . .
o ' ° 0 . *
14.0“l -I
0’
13.8.4 0 .
, J
E 13.6- . ..
E 1
4 1
13.44 -
4 . d
13.2- 4
‘ o
13.0- -
‘ o
12.3 . , . , . , . , . , . , - , . , .
2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0

<r>

Figure 8.8: The optimal supercell size which minimizes the total energy, Lmin, as
a function of the mean coordination (r). The length is given in A(Djordjevié and
Thorpe 1996).

185

 

 

 

 

 

 

 

 

 

 

 

 

 

7 I V I I ' r I r r 7 1 I v , . r v I
‘ M j
64 4 o (b)1
5 l l
5" 4 54 4
g . .
IL ‘1 4 ‘1 "
3 34 3( 4
S < (
.92
O 24 -l 24
.2 ‘
a 1
CL 14 1.1
o f ' I ‘ o a Ia I ' I
00 05 10 15 20 25 30 35 40 00 05 10 15 20 25 30 35 40
7 T fir * I I f V 1 I f 7 fl 7 v , v y r v 1 W f r f
‘ (cl‘ ‘ (d)
.4 T .4 ‘
4 l l
' 5~ “ 54 t
( 1
IL 44 ‘ ‘4 "
. l l
3'4 ‘l 34
' 2. 1 24
Q.

 

 

 

 

 

 

 

 

‘4 Id '
J 1 4
0 T— 0 v ‘1 c 4

00 013 1'0 1‘5 I 2'0 275 3'0 - 35 4° 00 05 1'0 (SI 2'0 ' {sfs'o ' 35 40
Distance rido Distance r/do

 

 

Figure 8.9: Pair distribution function calculated for (r) = 244(0), 275(0), 306(0),
3.38(d) (Djordjevié and Thorpe 1996).

186

random network has zero strain energy for all mean coordinations.

The slight increase of the me in the region from 2.7 < (r) < 4.0 could
be explained by the fact that the mean angle in our amorphous diamond model
is, (as it is usual in amorphous networks), slightly smaller than the tetrahedral
angle in the crystalline diamond structure which is 90 = 109.47°. Our amorphous
diamond network has the mean bond angle (0) = 109.31°. Thus, in the beginning
of the bond removal, the remaining bonds tend to relax toward the tetrahedral

angle which means the slight expansion of the supercell.

8.6 Conclusion

Earlier calculations of the elastic properties of amorphous networks were done
assuming that the bond dilution af a crystalline lattice and the subsequent re-
laxation, gives an acceptable sequence of amorphous networks, which is suitable
for the calculation of the elastic properties. In our simulations, we have used a
true amorphous structure for bond dilution, and calculated the bulk modulus as
a function of the mean coordination. Our results show that in both cases the
bulk modulus obeys practically the same power law, which suggests that for the
elastic properties, the initial structure (before the bond dilution takes place) is
not of a great importance. We have described the difficulties one encounters in
dealing with a true amorphous structure which contains large strain, which are
connected with different energy minima, i.e. different metastable states, that exist
in amorphous materials. Finally, we showed that for small mean coordinations a
collapse of the lattice cell occurs, due to the presence of many two—fold coordi-
nated atoms, which significantly reduces the exclusion volume of an average atom,

while the large strain in the structure triggers the collapse after a large fraction

 

187

of bonds is removed.

Chapter 9

Conclusion

In this thesis I have studied disordered properties of both continuous and dis-
crete random media. I focused on the two—phase composites in two dimensions,
random resistor networks, and on computer modeling of topologically disordered
networks — amorphous solids. Results for the effective conductivity and the spec-
tral function calculations, presented in Chapters 2, 3 and 4, are valid only in the
low concentration limit f < 1, when the inclusions can be regarded as well sep-
arated and thus non-interacting. In other words, the effective conductivity and
the spectral function is correctly obtained to the first order in conCentration f.
The only exception to this is the computer simulation results of Day and Thorpe
(1996), the results of the effective medium theory, which are valid in the entire
range of concentration f, and the results I obtained for the case of circular in-
clusions (Djordjevic' et al. 1996). I have shown that exact solution exists only
in a few cases, namely for elliptical (circular) inclusions. The conformal map-
ping technique turned out to be a very powerful tool for calculating the effective
conductivity for particular type of inclusions (holes or perfectly conducting inclu-

sions) which may have sharp corenrs.The problem becomes much more difficult

188

189

for inclusions with arbitrary 7 if they contain sharp corners. In this case one has
to solve numerically F redholm equation for the induced charge density around the
perimeter of the inclusion. It appears that the only direct numerical calculation of
the spectral function (from the general definition (2.19) or (3.105)) has been done
in computer simulations of Day and Thorpe (1996) on the discrete random square
net. This has been possible because of the use of a powerful Y — A transformation
technique (Frank and Lobb 1988) which calculates the effective conductivity of
a random network very efficiently. Finally, the effective medium theory seems
to give very good results for intermediate concentrations around the percolation
threshold fc, and thus proves itself useful for deducing the overall shape of the
spectral function.

I have shown that the spectral function for a sheet containing holes can
be calculated correct to second order in the area fraction of inclusions f. The
result is a truncated Lorentzian that for all practical purposes is a Lorentzian in
the range of interest. A knowledge of the spectral function permits the (complex)
conductivity, 0, for any (complex) values of 0'0 and 01 to be found, by doing a
single integral. This is particularly useful if there is a parameter in the composite
that can be varied at fixed geometry. Examples would be the frequency and the
temperature. The calculation in 2d was possible because an exact solution to
this problem is available in terms of multiple images. An analytic continuation
has been successfully done, which has allowed the spectral function to be found
which is finite only over a narrow region at low concentration. Similar spectral
functions can also be used in discrete systems made up of resistors (Foster 1924 and
Straley 1979). Spectral functions are widely used in physics to study electronic and

vibrational excitations. They have proved especially useful in systems containing

190

impurities, where local modes can occur (Elliott, Krumhansl and Leath 1974).
The weight in the absorption spectrum goes like the defect concentration f and
the width is determined by the interaction between defects and goes also as f.
The reason for this close analogy between these two somewhat disparate situations
is not entirely clear.

Chapters 5, 6, and 7, deal with the discrete random media, namely amor-
phous solids. I have studied computer modeling of amorphous tetrahedral semi-
conductors, particularly of amorphous diamond. Early attempts to build a model
of an amorphous structure by hand, even though clumsy, laid the path for the
implementation of computer algorithms which led to the creation of very success-
ful models. The idea of a microcrystallite nature of amorphous structure was
abandoned in favor of continuous random network concept. Periodic boundary
conditions were included, and increasing awareness of the importance of topology
in amoprhous materials led to the invention of the Wooten—Weaire topological
bond switch, which, implemented into a Metropolis simulated annealing algo-
rithm, proved to be a powerful way to create a topologically realistic amorphous
structures on computer. Designed originally for tetrahedrally bonded amorphous
elemental compounds, the Wooten—Weaire method is easily adaptable to treat
other systems as well, like vitreous silica, or amorphous mixture of diamond and
graphite—like carbon atoms. The continuous random network model which well
describes covalent glasses should be understood as an idealisation of real amor—
phous structures, in the same sense as a perfect crystal is an idealisation of real
crystals one can find in nature. CRN models can serve as a reference ideal amor-
phous structures to describe more realistic amorphous structures which usually

contain voids and various kind of defects.

191

I have computer generated a very large (4096 atom) model of amorphous
diamond (Djordjevié et al. 1995), and calculated all relevant quantities that de-
scribe the obtained amorphous structure. My results show that although carbon
has much stiffer bond-bending forces than silicon (i.e. a much larger value of
,3/0, the ratio of bond-bending to bond-stretching force constants in the Keating
model), this doesn’t influence the pair distribution function significantly (Djord-
jevic’ et al. 1995). The bond-angle deviations are slightly smaller in amorphous
tetrahedral diamond compared with amorphous silicon and amorphous germa-
nium, but they cannot be too small, as the atoms are constrained to be in a fully
tetrahedrally coordinated random network. Indeed our experience in constructing
random networks leads us to postulate that a rms bond-angle deviation of at least
7° is required if the network is not to contain crystalline regions and so be a two
phase mixture. Our very large and realistic structural model of amorphous dia~
mond also provides a good tool for studying other important issues of the physics
of amorphous materials, for instance, the nature of electronic localization in three
dimensions in the presence of topological disorder (Dong and Drabold 1996).

Finally, in Chapter 8, I have studied the elastic properties of amorphous
networks, using my large model of amorphous diamond as a starting structure.
Earlier calculations of the elastic properties of amorphous networks were done
assuming that the bond dilution af a crystalline lattice and the subsequent re-
laxation, gives an acceptable sequence of amorphous networks, which is suitable
for the calculation of the elastic properties (He and Thorpe 1985; Franzblau and
Tersoff 1992). In my simulations, I have used a true amorphous structure for bond
dilution, and I have calculated the bulk modulus as a function of the mean coordi-

nation. My results show that in both cases the bulk modulus obeys practically the

192

same power low, which suggests that for the elastic properties, the initial structure
(before the bond dilution takes place) is not of a great importance. I have de-
scribed the difficulties one encounters in dealing with a true amorphous structure
which contains large strain, which are connected with different energy minima, i.e.
different metastable states, that exist in amorphous materials. Finally, I showed
that for small mean coordinations a collapse of the lattice cell occurs, due to the
presence of many two—fold coordinated atoms, which significantly reduces the ex-
clusion volume around an average atom, while the large strain in the structure
triggers the collapse after a large fraction of bonds is removed.

I believe that this thesis could be a valuable comprehensive resource for
anybody interested in the current state and the historical background of these
particular fields of research. At the same time, I hope that the contribution I made,
especially in computer modeling of amorphous networks, will be useful for the
researchers who are interested in having large and realistic models of an amorphous
solid for calculations of electronic and other properties of these materials. I also
hope that my computer models will inspire more work on computer modeling of
amorphous materials itself, because I have no doubts that the era of computer

modeling is still at its beginnings.

Bibliography

Abramowitz, M. and Stegun, I. A. (1970), Handbook of Mathematical Functions.
New York: Dover.

Alben, R., Weaire, D., Smith, Jr. J.\E. and Brodsky, M. H. (1975), Phys. Rev.
B 11, 2271.

Anderson, P. W., Halperin, B. I. and Varma, C. M. (1972), Philos. Mag. 25, 1.
Batchelor, G. K. (1972), J. Fluid Mech. 52, 245.

Bell, R. J. and Dean, P. (1966), Nature (London) 212, 1354.

Bell, R. J. and Dean, P. (1972), Philos. Mag. 25, 1381.

Bennett, C. H. (1972), J. Appl. Phys. 43, 2727.

Bergman, D. J. (1978), Phys. Rep. 43, 377.

Bergman, D. J. (1986), in Homogeneization and Effective Moduli of Materials
and Media, edited by Ericksen, J. L., Kinderlehrer, D., Kohn, R. and Lions,
J. - L. New York: Springer, p. 27.

Bernal, J. D. (1959), Nature 183, 141.

Bernal, J. D. (1964), Proc. Roy. Soc. ondon A 280, 299.

193

194

Binns, K. J. and Lawrenson, P. J. (1973), Electric and Magnetic Field Problems.

Oxford: Pergamon Press.
Bruggeman, D. A. G. (1935), Ann. Phys. (Leipzig) 24, 636.
Cargill, G. S. (1975), in Solid State Physics 30, 227.
Chung, J. S. and Thorpe, M. F. (1996), private communication.

Churchil, R. V., Brown, J. W. and Verhey, R. F. (1974), Complex Variables and
Applications. New York: McGraw-Hill.

Cichocki, B. and Felderhof, B. U. (1988), J'. Stat. Phys 53, 499.
Connell, G. A. N. and Temkin, R. J. (1974), Phys. Rev. B 9, 5323.

Cusack, N. E. (1987), The Physics of Structurally Disordered Matter: An Intro-
duction. Bristol: IOP Publishing Ltd.

Day, A. R. and Thorpe, M. F. (1996), submitted to Phys. Rev. B.

Djordjevié, B. R., Hetherington, J. H. and Thorpe, M. F. (1996), accepted in
Phys. Rev. B.

Djordjevié, B. R., Thorpe, M. F. and Wooten, F. (1995), Phys. Rev. B 52, 5685.
Djordjevié, B. R. and Thorpe, M. F. (1996), unpublished.

Dohler, G. H., Dandoloff, R. and Bilz, H. (1983), J. Non—Cryst. Solids 42, 87.
Dong, J. and Drabold, D. A. (1996), submitted to Phys. Rev. Lett.

Drabold, D. A., Fedders, P. A. and Stumm, P. (1994), Phys. Rev. B 49, 23,
16415.

195

Duffy, M. G., Boudreaux, D. S. and polk, D. E. (1974), J. Non-Cryst. Solids 15,
435.

Duquesne, J, Y. and Bellessa, G. (1983), J. Phys. CzSolid State Phys. 16, L65.

Elliott, S. R. (1984), Physics of Amorphous Materials. Longman: London and

New York.

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

465.

Etherington, G., Wright, A. C., Wemzel, J. T., Dore, J. C., Clarke, J. H. and
Sinclair, R. N. (1982), J. Non-Crystalline Solids 48, 265.

Emns, D. L. and King, S. V. (1966), Nature 212, 1353.
Evans, D. L., Teter, M. P. and Borrelli N. F. (1974), A.I.P. Conf. Proc. 20, 218.

Evans, D. L., Teter, M. P. and Borrelli N. F. (1975), J. Non-Cryst Solids 17,
245.

Felderhof, B. U., Ford, G. W. and Cohen, E. G. D. (1982), J. Stat. Phys. 28,
649.

Felderhof, B. U. and Jones, R. B. (1989), Phys.Rev. B 39, 5669.
Felderhof, B. U. and Jones, R. B. (1986), Zeit. fiir Physik B 62, 215.
Feng, S. and Sen, P. N. (1984), Phys. Rev. Lett. 52, 216.

Feng, S., Thorpe, M. F. and Garboczi, E. (1985), Phys. Rev. B 31, 276.

Finney, J. L. (1970), Proc. Roy. Soc. A 319, 479, 495.

196
Foster, R. M. (1924), Bell System Technical J. 3, 651.
Frank, D. J. and Lobb, C. J. (1988), Phys. Rev. B 37, 302.
Franzblau, D. S. and Tersoff, J. (1992), Phys. Rev. Lett. 68, 2172.
Fuchs, R. (1975), Phys. Rev. B 11, 1732.

Gaskel, P. A., Seed, A., Chieux, P. and McKenzie, D. R. (1992), Philos. Mag.

B 66, 155.

Gilkes, K. W. R., Gaskell, P. H. and Robertson, J. (1995).. Phys. Rev. B 51,
12303.

Golden, K. and Papanicolau, G. (1983), Commun. Math. Phys. 90, 473.
Goodstein, D. L. (1985), States of Matter. Dover: New York.

Gradshteyn, I. S. and Ryzhik, I. M. (1965), Table of Integrals and Products. New

York: Academic Press.

Guinier, A. and Fournet, G. (1955), Small Angle Scattering of X—rays. London:

Bell.
Guttman, L. (1974), A.I.P. Conf. Proc. 20, 224.
Guttman, L. (1975), A.I.P. Conf. Proc. 31, 268.
He, H. and Thorpe, M. F. (1985), Phys. Rev. B 54, 2107.
Henderson, D. (1974), J. Non—Cryst. Solids 16, 317.
Hetherington, J. and Thorpe, M.F. (1992), Proc. R. Soc. Lond. A. 438, 591.

Hinsen, K. and Felderhof, B. U. (1992), Phys.Rev. B 46, 12955.

197

Hirai, H., Tabira, Y., Kondo, K., Oikawa, T. and Ishizawa, N. (1995), Phys.
Rev. B 52, 6162.

Jacobs, D. J. and Hendrickson, B. (1996), to be submitted to Siam J. Comput..
Jacobs, D. J. and Thorpe, M. F. (1995), Phys. Rev. Lett. 75, 4051.

Jacobs, G. J. and Thorpe, M. F. (1996), Phys. Rev. E 53, 3682.

Jeffrey, D. J. (1973), Proc. R. Soc. Lond. A. 335, 355.

Kakinoki, J., Katada, K, Hanawa, T. and Ino, T. (1960), Acta Cryst. 13, 171.
Kaplow, R., Rowe, T. A. and Averbach, B. L. (1968), Phys. Rev. 168, 1068.
Keating, P. N. (1966), Phys. Rev. 145, 637.

Keller, J. B. (1964), J. of Math. Phys. 5, 548.

Kirkpatrick, S., Gelatt, Jr. C. D. and Vecchi, M. P. (1983), Science 220, 671.

Landau, L. D. and Lifshitz E. M. (1960) Electrodynamics of continuous media,

§9. London: Pergamon Press.
Lee, C. H., Lambrecht, R. L. and Segall, B. (1994), Phys. Rev. B 49, 11448.

Lossy, R., Pappas, D. L., Roy, R. A., Cuomo, J. J. and Sura, V. M. (1992), Appl.
Phys. Lett. 61, 171.

Marks, N. A., McKenzie, D. R. and Pailthorpe, B. A. (1996), Phys. Rev. B 53,
4117.

Marks, N. A., McKenzie, D. R., Pailthorpe, B. A., Bernasconi, M. and Parrinello,
M. (1996), Phys. Rev. Lett. 76, 768.

198
Martin, R. M. (1970), Phys. Rev. B 1, 4005.
Mathews, J. and Walker, R. L. (1970), Mathematical Methods of Physics, p. 481.
Maxwell, J. C. (1864), Philos. Mag. 27, 294.
Maxwell, J. C. (1873), Electricity and Magnetism, Oxford: Clarendon, p. 365.

Mckenzie, D. R., Martin, P. J., White, S. 8., Liu, Z., Sainty, W. G., Cockayne,
D. J. H. and Dwarte, D. M. (1988), in Proceedings of the European Materials
Research Society Conference, Strasbourg, edited by P. Koidl and P. Oelhafen,

p. 203.

McKenzie, D. R., Muller, D. A. and Pailthorpe, B. A. (1991), Phys. Rev.
Lett. 67, 773.

Mendelson, K. S. (1974), J. of Appl. Phys. 46, 917.

Metropolis, N ., Rosenbluth, A., Rosenbluth, M., Teller, A. and Teller, E. (1953),
J. Chem. Phys. 21, 1087.

Milton, G. W. (1981), J. Appl. Phys. 52, 5286.; ibid. 5294.

Moss, S. C. and Graczyk, J. F. (1969), Phys. Rev. Lett. 23, 1167.
Mousseau, N. and Lewis, L. J. (1990), Phys. Rev. B 41, 3702.
Mousseau, N. (1996 (3.)), private communication.

Mousseau, N. (1996 (b)), private communication.

Mozzi, R. L. and Warren, B. E. (1969), J. Appl. Cryst. 2, 1905.

199

Muskhelishvili, N. I. (1963), Some Basic Problems of the Mathematical Theory
of Elasticity. Groningen: Noordhoff. pp. 376-380.

Nehari, Z. (1952), Conformal Mapping. New York: Dover.
O’Mard, L. P. (1993), Modelling Simul. Mater. Sci. Eng. 1, 485.
Peterson, J. M. and Hermans, J. J. (1969), J. Composite Materials 3, 338.
Phillips, J. C. (1979), J. Non—Cryst. Solids 34, 153.

Phillips, J. C. (1980), Phys. Status SolidilOl, 473.

Phillips, W. A. (1985), J. Non—Cryst. Solids 77 8c 78, 1329.
Phillips, W. A. (1972), J. Low Temp. Phys. 7, 351.

Phillips, J. C. (1981), J. Non—Cryst. Solids 43, 37.

Polk, D. E. (1971), J. Non-Cryst. Solids 5, 365.

Polk, D. E. and Boudreaux, D. S. (1973), Phys. Rev. Lett. 31, 92.
Rayleigh, Lord (1892), Phil. Mag. 34, 481.

Ross, D. K. (1968), Aust. J. Phys. 21, 817.

Shevchik, N. J. (1973), Phys. Status Solidi B 58, 111.

Silverman, R. A. (1974), Complex Analysis with Applications. New Jersey:
Prentice-Hall.

Smith, D. A. (1979), Phys. Rev. Lett. 42, 729.

Smythe, W. R. (1939), Static and Dynamic Electricity. New York: McGraw-Hill.

200
Steinhardt, P., Alben, R. and Weaire, D. (1974), J. Non-Cryst. Solids 15, 199.
Stillinger, F. H. and Weber, T. A. (1985), Phys. Rev. B 31, 5262.
Straley, J. P. (1979), J. Phys. C 12, 2143.
Thorpe, M. F. (1992), Proc. R. Soc. London A437, 215.
Thorpe, M. F. (1983), J. Non—Cryst. Solids 57, 355.
Thorpe, M. F. (1985), J. Non—Cryst. Solids 76, 109.
Thorpe, M. F. (1995), J. Non—Cryst. Solids 182, 135.

Thorpe, M. F., Djordjevié, B. R. and Hetherington, J. H. (1994), Physica A 207,
65.

Thorpe, M. F., Jacobs, D. J. and Djordjevié, B. R. (1996), to be published in
the Proceedings of the XIX International Workshop on Condensed Matter

Theories, held in Caracas, Venezuela, June 1995.
Thorpe, M. F. and Tang, W. J. (1987), J. Phys.: Condens. Matter 20, 3925.
Vanderbilt, D. and Louie, S. G. (1984), J. Comp. Phys. 56, 259.
Vukcevich, M. R. (1972), J. Non—Cryst. Solids 11, 25.
Wang, C. Z. and Ho, K. M. (1993), Phys. Rev. Lett., 71, 1184.
Wooten, F. and Weaire, D. (1984), J. Non—Cryst Solids 64, 325.
Wooten, F. and Weaire, D. (1986), J. Phys. C:Solid State Phys. 19, L411.

Wooten, F. and Weaire, D. (1987), Solid State Physics 40, 1.

201
Wobten, F. and Weaire, D. (1989), J. Non—Cryst. Solids 114, 681.

Wooten, F. and Weaire, D. (1996), in Adaptation of Simulated Annealing to

Chemical Problems, edited by J. Kalivas, Elsevier Science, (in press).
Wooten, F., Winer, K. and Weaire, D. (1985), Phys. Rev. Lett. 54, 1392.

W'right, A. C., Connell, G. A. N. and Allen, J. W. (1980), J. Non—Cryst. Solids 42,
69.

Zachariasen, W. H. (1932), J. Am. Chem. Soc. 54, 3841.

Zallen, R. (1983), The Physics of Amorphous Solids. New York: John Wiley &

Sons.

Ziman, FRS. J. M. (1979), Models of Disorder. Cambridge: Cambridge Univer-

sity Press.

   

"I111111111!liltillt