puma-2,":

L

LllCl-i.1l:|;.l. [A

 

E

I S-

 

 

 

LAi.

._|-;..",e.= 313413; in

f

 

 

 

 

 

  

 

4 .,, . .
‘ . ,4‘,,‘.,‘,
¥ ». ‘. r.”

 

 

 

 

(llllllllllllllllHllllllllllllllllllllllllllllllllllllllll

3 1293 00908

This is to certify that the

dissertation entitled

ToroLOJiCaL mjiotwrj amt Struatuml.
Mud/j of ”DisaroCMCA gJVfP/Wl

presented by

Yong C Okl.

has been accepted towards fulfillment
of the requirements for

EA . D . degreeinwzg‘

 

    

Ma or professor

Datel{ it“: quI

MSU is an Affirmative Action/Equal Opportunity Institution 0-12771

 

 

 

l LmAFv
lMéetaigan “3m
l 0 '23.

l Unwersmi

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.

 

 

 

 

 

 

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU Is An Affirmative Action/Equal Opportunity Institution
c:\c|rc\datm.nm3-p. 1

 

' TOPOLOGICAL RIGIDITY
AND STRUCTURAL STUDY
OF DISORDERED SYSTEM

By

Yong Cai

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

1991

é54“ @906

ABSTRACT

TOPOLOGICAL RIGIDITY
AND STRUCTURAL STUDY
OF DISORDERED SYSTEM

By

Yong Cai

In this thesis we have studied the elastic and structural properties of disordered
systems. Our research provides a good understanding of the size mismatch problem
frequently encountered in intercalated multilayer systems and semiconductor alloys.
By studying the short range and long range properties of these disordered systems in
various dimensions, it is possible to draw a clear physical picture by introducing a set
of topological rigidity parameters which can be understood from an effective medium
point of view.

The observation of soft phonon and elastic properties of network glasses can be
understood in term of the rigidity percolation concept. This study adds a new example
to the application of rigidity percolation theory. It also provides an explanation for
the soft phonon modes observed in covalent glasses.

The size mismatch problem is interesting because of its wide ranging application
to semiconducting materials. Using analytic solutions of the simple model and the
effective medium treatment of realistic models, we studied the intercalated graphite
multilayer systems, 2d triangular lattice, 3d FCC lattice and 3d diamond lattice. Our
theoretical results suggest simple structural expressions for average nearest neigh-
bor distances and corresponding distributions in term of dimensionless topological

parameters. From the point of View of composition, we studied binary, ternary and

 

 

quaternary system. Indeed, it is a theoretical success that we have been able to obtain
the analytic solutions for quaternary system which contains all compounds we studied
before as special cases. This leads to a general formalism for the length mismatch
problem in semiconducting alloy systems.

The investigation of the long range properties of 3d disordered systems leads to the
concept of “effective temperature for disorder”. Huang scattering is studied and we
have a better understanding of the effect of static disorder on diffraction patterns.

The mismatch problem in 2d is especially interesting because this(2d) happens to be
a marginal dimension. We observed various instabilities in computer simulations. The
most interesting one is the “static melting”, which is a zero temperature equivalence of

thermal melting but driven by disorder. Research in this direction is still in progress.

 

ACKNOWLEDGEMENTS

With great respect and deep gratitude I acknowledge my adviser, Professor M.
F. Thorpe, for his professional guidance and endless encouragement in this work.
His scientific research attitude, his physical insight of complicated phenomena, and
his constant searching for simple truth has been and will always be the greatest
inspiration to me. I thank him for his understanding, supporting and patient advising
during the last four years.

I am grateful to Professor S. D. Mahanti, Professor W. Bauer, Professor N. Birge
and Professor M. Dubson for their service on my Ph. D. guidance committee, and
I would like to express my special thanks to Professor J. Kovacs for his patient
assistance and Mrs. Janet King for her help.

I would like to thank Professor S. D. Mahanti, Professor P. Duxbury, Professor D.
Tomanek, Professor W. Bauer, and Professor T. Kaplan for their valuable lectures,
advice and many helpful conversations. Their help and influence on my graduate
study are greatly appreciated.

My thanks also go to my family for their support and encouragement, and my dear
friends: W. Q. Zhong, J. D. Chen, N. J. Jiu, N. Mousseau, B. Pi, x. Hqu, Y. Wang,
Q. F. Zhu and others in this department for their many interesting discussions during
my studying in MSU and help in preparation of the thesis. I acknowledge Dr. Y. S.
Li and Dr. G. Overney for their valuable friendship.

Financial support from Michigan State University, the National Science Foundation

and the Center for Fundamental Material Research are gratefully acknowledged.

iii

it

Contents
1 INTRODUCTION 1
2 FLOPPY MODES IN NETWORK GLASSES ‘ s
2.1 INTRODUCTION ............................ 8
2.2 CONSTRAINT COUNTING ...................... 10
2.3 ELASTIC MODULI ........................... 12
2.4 DENSITY OF STATES ......................... 18
2.5 DEBYE WALLER FACTOR ...................... 24
2.6 CONCLUSIONS ............................. 3O
3 RIGIDITY OF INTERCALATED LAYERED SOLIDS 34
3.1 INTRODUCTION ....................... - ..... 34
3.2 THE MODEL .............................. 40
3.3 ANALYTIC SOLUTION ........................ 42
3.4 CORRELATED CASE ......................... 52
3.5 EFFECTIVE MEDIUM THEORY ................... 55
3.5.1 Random Intercalation ...................... 55
3.5.2 Correlated System ........................ 58
3.6 PROBABILITY DISTRIBUTION FOR GALLERY HEIGHTS . . . . 63
3.6.1 Distribution Function PA(d) and P3(d) ............. 65
3.6.2 Low and High Concentrations .................. 68

iv

3.7 CONCLUSIONS ............................. 70

STRUCTURAL CHARACTERISATION OF SEMICONDUCTOR

ALLOYS ' 73
4.1 INTRODUCTION ............................ 73
4.2 VALENCE FORCE MODELS ..................... 74
4.3 THEORY ................................. 78
4.3.1 Solution to Pure Length Mismatch ............... 79
4.3.2 Watson Integrals as Function of Force Ratio .......... 81
4.3.3 Effective Medium Theory .................... 81
4.4 APPLICATION ............................. 84
4.5 CONCLUSIONS ............................. 93
DIFFRACTION FROM RANDOM ALLOYS 96
5.1 INTRODUCTION ............................ 96
5.2 COMPUTER SIMULATIONS ..................... 97
5.2.1 Bragg Scattering ......................... 99
5.2.2 Huang Scattering ......................... 101
5.3 THEORY ................................. 101
5.3.1 Bragg Scattering .................... ‘ ..... 101
5.3.2 Huang Scattering ......................... 105
5.4 RESULTS ................................ 108
5.4.1 Bragg Scattering ......................... 108
5.4.2 Huang Scattering ......................... 111
5.5 DISCUSSION .............................. 111
5.6 CONCLUSION .............................. 113
TWO DIIVIENSIONAL MIXED CRYSTAL 117
6.1 INTRODUCTION ............................ 117
6.2 SHORT RANGE PROPERTIES .................... 121

 

6.3 DIFFRACTION PATTERN ......................
6.4 ANGULAR CORRELATION ......................
6.5 FINITE SIZE SCALING ........................
6.6 CONCLUSION ..............................

OVERVIEW AND PERSPECTIVE

7.1 LENGTH MISMATCH PROBLEM ...................
7.1.1 Quaternary Solution and Sublattice Decoupling ........
7.1.2 Topological Rigidity Constants .................

7.2 MODEL AND MATHEMATICS ....................

7.3 TWO DIMENSIONAL INSTABILITY .................

The Equation of Motion
Distribution Function of Multilayer System
Symmetric Distribution Case

Average NN-distance in Semiconductor Alloys

vi

130

138
138
138
148
150
155

160

163

166

168

List of Tables

4.1 force constants ..............................

4.2 Formula constants .............................

vii

 

81

 

 

List of Figures

2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9

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

Modulus versus mean coordination ................... 13
Cu with weak force on .......................... 16
C12 with weak force on .......................... 17
Density of state .............................. 21
Density of state .............................. 22
p(w)/w ................................... 23
p(w) = awzlwzo .............................. 25
The —1 moment of DOS ......................... 27
The ——2 moment of DOS ......................... 29
Intercalated Graphite ........................... 36
Intercalated Graphite with impurity‘ .............. - ..... 37
Intercalation of two type atoms ..................... 38
Stronger layer with intercalation ..................... 39
Z-curves .................................. 47
Strain energy curves ........................... 48
Width curves ............................... 49
Anticlustering intercalation ....................... 50
Clustering intercaltion .......................... 51
Z-curves .................................. 59
Strain energy ............................... 60
width curves ................................ 61

viii

3.13 probability of distribution ........................ 64

4.1
4.2
4.3
4.4
4.5
4.6
4.7
4.8

5.1
5.2
5.3
5.4
5.5
5.6

6.1
6.2
6.3
6.4
6.5
6.6
6.7

7.1
7.2
7.3
7.4

Watson integral .............................. 82
“Z” curves for III — V compounds ................... 85
“Z” curves for III — V compounds ................... 86
“Z” curves for II — VI compounds ................... 81
“Z” curves for II — VI compounds ................... 85
“Z” curves for I n1_,Ga,As ....................... 82
“Z” curve for Cd1-xanTe ........................ 9(
Length distribution of Cdo,5Zno,5Te. .................. 9‘2
Triangular lattice ............................. 98
Bragg scattering .............................. 10C
Huang scattering ............................. 102
Huang scatter(larger mismatch) ..................... 102
Mean squared displacement ....................... 10S
Wilson plot ................................ 11(
Triangular lattice ............................. 11S
Triangular lattice with mismatch .................... 12(
The density-density correlation function ........... M ..... 12.I
Diffraction pattern ............................ 12'
Triangular lattice with pleating ..................... 12l
Pleating transition ............................. 12!
Angular correlation function ....................... 13f
Quaternary Distribution ......................... 14‘
Si1_zGex ................................. 15
Ga1-,In,As ................................ 15
Next nearest neighbor distance ..................... 15

ix

7.5 Melting of 2d LJ system ......................... 157
7.6 Topological Defects in 2d ......................... 158

 

Chapter 1

INTRODUCTION

Solid state theory first was a theory for crystals, for the simple reason that one learns
walking before running. Systems with periodicity have been well studied, and we have
beautiful and powerful mathematical means (Bloch theory, Group theory) in hand.
The study on amorphous systems where long—range-order no longer exists became
active thereafter. There are strong demands from application because of the various
optical and electronic properties of the amorphous materials, and also because of the
richness of such materials with their imperfection. Examples of applications include
the use of ultratransparent optical fibers in telecommunications, the use of amorphous
semiconductors in xerography and solar cells, and the ubiquitous everyday uses of or-
ganic glasses as structural materials. Scientifically, it has been a long struggle to
understand the real beauty behind the seemingly disordered and constantly changing
world. It is necessary to first have a good understanding of crystals. Based on that,
perturbation theories of various kind could be performed. Yet new concepts and new
methods are needed to deal with disorder comprehensively. Localization and percola-
tion theory are the main ones that have met the new challenges. Some old approaches,
such as the chemical bonding vieWpoint, remain useful. The most important thing in
studying amorphous material is the scientific insight. Our approaches are based on

studying simple models which carry the main features of the corresponding physical

 

systems and also have soluble limits. The effort is to sort out the essence of the
physical phenomena that have the dominant effect on the observations. In this thesis
we will discuss the application of chemical bonding, which is related to short range
order, topological rigidity and rigidity percolation.

This dissertation is basically a collection of the works the author has been in-
volved in during the Ph. D. program. Chapter 2 is published in Phys. Rev. B40,
10535(1989), it is finished under the guidence of Prof. M. F. Thorpe. The author also
makes major contributions both in theory and computation to the works described
in Chapter 3 (Phys. Rev. B42, 8827(1990)), Chapter 4 and 7 (to be sent to Phys.
Rev. B). Prof. M. F. Thorpe worked out most of the theory described in Chapter 5
(Phys. Rev. B43, 8282(1991)) and 6 (Phys. Rev. B43, 11019(1991)) and the author
did the major work in computer simulations and part of the work on the theory.

Because each chapter can be treated as an independent paper, and as a matter
of fact many of the papers have been published, there are some notation differences
between these chapters. For example, we use “W” to denote the Watson Integrals in

it”

Chapter 3, and in other chapters we use “a . In some chapters we call the Watson
Integrals the “Topological Rigidity Parameters” because of better understanding of
its physical meaning. I hope this explanation could ease the inconvenience created
to the readers because of these changes in notation. Bibliographies also follow each
chapter as in an independent paper.

The next chapter presents the concept of rigidity percolation. The percolation con-
cept is very important in disordered systems. It is a powerful unifying construction
and an outstanding vehicle for exhibiting many of the modes of thought characteristic
of theoretical approaches to strongly disordered systems: emphasis on statistical dis-
tributions, transition between localization and delocalization, critical points, scaling
behavior, and dimensionality dependence etc.

Rigidity percolation is about the elastic response of a system. In a lattice net-

work, rigidity percolation always occurs after conductivity percolation. The reason is

because conductivity percolation concerns the connectivity while the rigidity perco-
lation concerns the constraints. Constraint on a lattice site implies the existance of
connectivity of this site to others. There has been questions whether these two types
of percolation belong to the same universal class. They are not necessarily the same
because of different symmetry. And there have been theoretical work that prove it
by presenting different critical exponents. Beyond the universality problem, rigidity
percolation is a useful concept in investigating the mechanical properties of covalent
glasses.

In Chapter 2 we discuss how the concepts from rigidity percolation can be used to
understand the low—frequency excitations and elastic properties of network glasses like
Get/131,3 61-3-1” When the mean coordination (r) = 2+ 22: + y is low, these materials
are soft and their properties are strongly influenced by low-frequency phonons. We
use a bond-depleted diamond lattice to mimic the coordination properties of the glass.
We show that a model with only covalent forces is unstable for (r) < 2.4 but can be
stabilized by small additional forces. We calculate the elastic constants, the density
of states, and the Debye—Waller factor as a function for (7'). Despite the Simplicity of
the model, the phase transition at (r) = 2.4 is consistent with recent experimental
results involving ultrasonics, inelastic neutron scattering, and Méssbauer experiments
on GexSe1-x glasses. The floppy modes, observed in inelastic neutrbn scattering,
disappear as (1') increases from 2 up to 2.4.

” There are

The later chapters mainly deal with the “length mismatch problem.
many approaches to disordered materials. People have studied the mass and com-
pressibility disorder. The calculation of electronic band structure of semiconducting
alloys needs the configuration of the compound system. Theoretical calculations usu-
ally adopt the virtual crystal model. In this model the positional disorderis averaged
out and the lattice constant is taken to be an averaged constant. One reason comes

from the simplicity after such averaging and the other is that the detailed structural

information needed is unavailable for technical reasons. The predictions are many

 

times not very good. And the need for structural information is essential for any
improvement. Our research focuses on the disorder caused by size difference. The
reason lies in the fact that the differential size effect has very dominant effect on
the structural properties of a large group of materials, like II — VI and I I I — V
ternary compounds such as I nGaAs, and quaternary semiconductor alloys such as
I nPGaAs. On the experimental side, EXAFS experiments are now capable of detect-
ing the nearest and next nearest neighbor distances. This parallel development is very
encouraging and provide checks to the ideas used in the theory. And the feedback is a
better and clearer understanding of the physical phenomena from a microscopic point
of view. The models we have used are very simple, which can be seen through all the
chapters. But they grasp the main features of the problems. For the size mismatch
problem, the simplest model is the lattice network of “bond” disorder. It is not the
best to represent an alloy system, but is pretty good in describing the intercalated
graphite multilayer system. For binary and ternary compounds, a better model is a
lattice network of “site” disorder. We studied both models. The useful concepts are
the chemical bonding and atomic radii. Pauling first studied these concepts. The
stability of chemical bonding leads to the introduction of the natural length of that
bond. It is similar to the natural length of a Hook spring. The additivity of atomic
radii to give natural bond length is very useful because our model Can be solved
exactly under such conditions. In reality, the approximation we made are good for
many systems and provide a good starting point for others where the conditions are
not satisfied. For example, there are disorder in the local compressibility, and the
natural bond lengths may not be additive from the atomic radii. It is even possible
that natural length is not a very stable concept when charge transfer is important
and varies from one environment to another. From the perturbation analysis we can
understand the condition for the validity of Vegard’s law and the trend of the devi—
ation. For instance, effective medium theory can be introduced to describe systems

when there is also disorder in the force constants.

In Chapter 3, we investigated the gallery structure of intercalated multilayer sys-
tems A]_,BIL, assuming that the intercalants A and B are randomly distributed
in each gallery flanked by host layers L. We allow for the correlation between the
intercalants in different galleries denoted as interlayer correlation (ILC). The gallery
structure is characterized by an average gallery height, by fluctuations in the gallery
height from site to site, and by a complete height distribution function. These quan-
tities have been studied by constructing a harmonic spring model which incorporates
the ILC through a one—dimensional Ising spin model. We have obtained an exact
solution for this model for arbitrary ILC when the two intercalants have the same
compressibility (kA 2 k3). In this limit, we Show that Vegard’s law is obeyed for all
values of the transverse rigidity of the host layers and arbitrary ILC. In the general
case (kA # k3), we have developed an effective—medium theory whose results are in
excellent agreement with computer simulations.

While the above gallery model is very useful in its own right in studying the in-
tercalated graphite multilayer system, it is the first step toward the length mismatch
problem we have studied in the alloy systems.

In Chapter 4 we have made an exhaustive study of all pseudobinary semiconductor
alloys (A1_,Bx)C. This includes both the III — V and the II - VI alloys. We use
a Kirkwood model with parameters derived from the elastic constants of the pure
materials. It is shown that the mean lengths are linear in the composition a: if the
spring constants for the two pure materials AC and BC are the same. This is because
the displacement fields due to the various inclusions can be superimposed. Curvature
is caused by the many body effects in the concentrated alloy that occur when the force
constants for AC and BC differ. We obtain analytic results for the mean lengths for
AC and BC bonds and for their widths using effective medium theory, which is in
excellent agreement with the results of computer simulations of the same models. The
length distribution function for nearest neighbor bonds of a particular type, is shown

to be well approximated by a Gaussian.

What we studied In Chapter 4 are only the short range order properties. It should
be pointed out that the mismatch problem can be solved in the quaternary model
when the structural information are concerned. We discuss this further in the last
chapter.

Chapters 5 and 6 deal with not only the short range properties but also with
the long range properties in a length mismatched lattice network. The studies are
performed for 3d systems and 2d systems. In Chapter'5, we calculate the diffraction
from random alloys, consisting of bonds of different natural lengths. The resulting
structure is a distorted crystal with strong correlations between displacements at dif-
ferent sites. The diffraction pattern is dominated by a set of Bragg peaks whose
intensity is modulated by a Debye-Waller factor. At the base of each of these Bragg
peaks, we find divergent Huang scattering due to the second order correlations be-
tween displacements. Higher order correlations produce a negligible non-divergent
background. Results are obtained in closed form and compared with computer sim-
ulations.

In Chapter 6 we show that the long range positional order in a two dimen—
sional crystal A148,, is destroyed by size-mismatch disorder as well as by thermal
disorder. The size-mismatch disorder is characterized by an effective temperature
kBTD = K a:(1 — 2:)(L% — L902 associated with the strain energy, where K is a force
constant, LOB -— L9, is the length mismatch and :1: is the concentration. We study a
model which has a fixed triangular net topology. By comparing with computer sim-
ulations, we show that a linearized small displacement theory is adequate for small
size-mismatches. The long range orientational ordcr remains. The positional corre-
lation function decays algebraically, which leads to power law peaks that replace the
Bragg peaks in the diffraction pattern. We argue that this model should provide a
reasonable qualitative description of real two dimensional mixed crystals, in the limit
of small size-mismatch.

When we relax the topological constraint by replacing the harmonic spring potential

by, for example, a Lennard—Jones potential, the study on the 2d alloy system leads
us to the instability problem, which we referred to as 2d length-mismatch-driven-
melting. We call it melting in the sense that the shear modulus drops dramatically
when mismatch reaches certain value and the system becomes very soft. We come
back to this discussion in the last chapter, too.

We summarize our main theoretical results in the last chapter from a mathematical
point of View and discuss the models we used. We discuss the role played by topo-
logical rigidity parameters that have frequently appared in the above chapters. After
giving a summary of the “length mismatch problem”, we discuss the 2d instability

which is very interesting and not completely solved yet.

Chapter 2

FLOPPY MODES IN
NETWORK GLASSES

This chapter is published in Physical Review B. Reference: Y. Cai and M. F. Thorpe,
Phys. Rev. B40, 10535(1989).

2.1 INTRODUCTION

Many network glasses are well described by the Continuous Random Network (CRN)
model first discussed by Zachariasen[1] more than 50 years ago. Examples are Si and
Ge that like to make 4 bonds in a tetrahedral arrangement; As that makes 3 bonds;
and 5, Se and O that make 2 bonds. A CRN can be described by a chemical formula
like G'exAsyS 61-3-3, (0 < x, y < 1). If the total number of atoms is N and there are

n, atoms with coordination r (r = 2, 3 or 4) then

N=Znn (2.1)

and we can define the mean coordination (r) by

(r) =Zmr/Zn. =2+2x+y (2.2)

We note that (r)(2 < (r) < 4) gives a partial, but very important, description of
the network. We will also see that the elastic properties depend sensitively on (r).

We regard a perfect CRN with no dangling bonds as being analogous to the perfect
crystal. In practice one can approach a perfect CRN by careful preparation techniques
in which the number of dangling bonds is kept to a minimum. In some circumstances
the CRN model is inappropriate and the network has effectively lower dimensionality
and may be better thought of as a spaghetti of twisted ribbons.[2] The ideas that we
are discussing in this chapter also apply to such networks if suitable adjustments in
the theory are made.

The standard procedure for constructing a CRN is to hand build a model using
hundreds to a few thousand appropriate chemical units.[3] More recently, such net-
works have also been made by restructuring diamond cubic networks.[4],[5] This has
led to satisfactory networks for Si and for $20 . Molecular dynamics techniques
have also been used to construct such networks, starting from a random configu-
ration of ions interacting with prescribed potentials.[6] Both these computer based
techniques are superior to hand building because they are unprejudiced by the builder
and have periodic boundary conditions. For a given size of model, periodic boundary
conditions are preferable as they maintain the coordination. While this is of some
minor help in calculating static quantities, like the radial distribution function, it is
of much greater importance for dynamic quantities, like the density of states, that
we will be mainly concerned with in this chapter. Unfortunately to date, neither
of the computer generated techniques mentioned above has been used to build non-
stoichiometric networks, although molecular dynamics could probably be adapted to
do so. Hand building a series of networks for the studies in this chapter would be
prohibitively time consuming and would yield structures with free boundaries.

Our belief and experience is that the coordination of the constituent ions is the
main determining feature of these networks for vibrational properties.[7],[8] We have

therefor used a simple bond depleted diamond cubic crystal lattice. This model

10

was first introduced by He and Thorpe[8] in a calculation of the elastic moduli as a
function of (r). Bonds are removed randomly to produce 3 and 2 fold coordinated
sites. If the removal of a bond would result in a dangling bond, that move is not
allowed. He and Thorpe showed that the values of :r,y in Eq.(2.2) were important
only through (1‘) by also removing bonds in a correlated way. That is if all the sites
were 3 coordinated (representing As), the elastic properties were essentially the same
as those for a 50 : 50 mixture of 2 and 4 coordinated sites. This assumes of course
that the forces between the atoms are independent of the chemical species. We will
comment later on this more.

The layout of this chapter is as follows. In the next section we review the constraint
counting arguments that lead to a transition at (r) = 2.4. In Sec. 2.3, we compute the
elastic moduli with the covalent forces only, and with the covalent plus weak forces.
In Sec.2.4, we calculate the density of states using the equation of motion method.
The low frequency floppy modes are seen to agree well with recent experimental
measurements in GexSe1-3. In Sec.2.5, we use the results from the density of states
to calculate the Debye Waller factor that has recently been measured in a Mdssbauer
experiment. In all cases, quantities are calculated as functions of (r) and comparisons
are made with experimental measurements on GexS 61-3,. The equation of motion
technique is described in Appendix A in two different versions; the conventional
one involving “plucking” and a variation involving a “kick start” that is necessary to

get good results at very low frequencies.

2.2 CONSTRAINT COUNTING

In covalent networks like GexAsySel_3_y , the bond lengths and angles are well
defined and small displacements from a (supposed) equilibrium structure can be de-
scribed by the potential[8]

a . .6
= 5 gm; — 113') ' r.,-]2 + 5 E; lijljk(A9ijL-)2- (2-3)
‘1] U

11

We will refer to this as the covalent potential as it has bond stretching and bond
bending terms. Here u,- is the displacement of the atom i and 1",,- is a unit vector
connecting nearest neighbor sites 2', j; 1;,- is the length of the bond i j and 0,3,, is the
angle between the bonds 2' j and j k. The bond stretching force constant a and the
bond bending force constant ,3, are the largest forces in covalent glasses. Other forces
are smaller as would be expected from the nature of the covalent bond. The potential
in Eq. (2.3) has been written as harmonic but this is not important in what follows
in this section as we shall be making small virtual displacements. The detailed form
of the potential will be important when we come to calculate spectral responses later.

We will regard the solution of the eigenmodes of the potential in Eq.(2.3) as a
problem in classical mechanics .[9] The dynamical matrix has a dimensionality 3N
which corresponds to the 3N degrees of freedom. In a stable rigid network we would
expect all the squared eigenfrequencies (.02 > 0 with 6 modes at zero frequency. These
6 modes are just the 3 rigid translations and the 3 rigid rotations and have no weight
in the thermodynamic limit. However we will find that in some covalent networks,
described by the potential in Eq. (2.3), there are a significant number of additional
zero frequency modes which correspond to internal deformations of the network. The
number of such modes can be determined by counting the constraints in the system.
There is one constraint associated with each bond which we can assign as r/2 con-
straints associated with each 7' coordinated atom. In addition there are constraints
associated with the angular forces in Eq. (2.3). For a 2 coordinated atom there is
one angular constraint, for a 3 coordinated atom there are three angular constraints
and for an r coordinated atom there are a total of 27‘ — 3 angular constraints.[7] The
total number of constraints is equal to the number of finite frequency modes which

are therefore
2 n,[r/2 + (27‘ — 3)] (2.4)

The fraction f of zero frequency modes is given by

f = {3N — an[r/2 + (27‘ — 3)]}/3N (2.5)

 

12

which can be rewritten as

5

f = 2 — 6(7‘) (2.6)

where (r) is defined in Eq. (2.2). When (1') = 2(e.g. Se chains), then f = 1 / 3; that is
one third of all the modes are at zero frequency. As atoms with higher coordination
than 2 are added to the network, f drops and goes to zero at (r) = 2.4 and the

network becomes rigid.

2.3 ELASTIC MODULI

In order to simulate a CRN , bonds were removed at random from a series of diamond
lattices in 'such a way that sites with 2, 3 and 4 fold coordination were produced. The
networks had periodic boundary conditions and the supercell contained 216 atoms.
This was the largest network that could easily be handled when the weak forces are
included. We found that the results were not very size dependent by looking at a few
networks with 512 atoms. If removing a bond would result in a dangling bond that
removal was not allowed. Using this technique, networks with (r) as low as 2.082
can be produced. The covalent potential in Eq. (2.3) is rewritten to apply more

conveniently to this situation as,

V = g: IDs-(u..- ' 1‘23)? + g X (”Pig'PMuz'j ° file + 11.1. ° fa)?’ (2-7)

(m) (iik)
where the probability that the bond i j is present or missing is included by letting ng
be either one or zero. Here all the bond lengths l are the same and Uij = u.- — Uj .
The angular force involving 6 is not quite the same in Eq. (2.7) as in Eq. (2.3) but
we have found the potential (2.7) a little easier to use in numerical calculations. It
can be seen from the insert in Fig.(2.1) that all the elastic constants 011,012 and C44

go to zero around (7‘) = 2.4.

In order to stabilize the lattice and move the floppy modes away from zero frequency,

 

13

 

 

 

 

 

 

 

- - v 1 w i g V
A i > I on
g , i n o cu
0‘ ’ a .. D 1
a » in
v 1000-- o. 4 -fi 0 ..
g 8 4 a
i -‘ n o ”
a . o f
g . o o * * + _..
0— ° 4»
m D 2 .g 1
> D ‘
’ .
0 fl - - - i r -
2 3 ,‘
(1')

Figure 2.1: The elastic moduli 011,012 and B with fl/a = 0.3 as a function of the
mean coordination (r). The insert is without the weak forces and the main graph has

the weak forces included.

14

we add a weak stabilizing potential V, to the potential V given in Eq. (2.7),
V = :27“; PiijkPkl(A¢£j,kl)2 + 92520 - P.,-)(u.-,- ' 9.5)? (2.8)
:3 t:

The first term in Eq. (2.8), proportional to 7, is the dihedral angle force involving
the dihedral angle ¢ij’k[ which by itself is sufficient to stabilize all zero frequency
modes. However we found it convenient to add a small “interchain” term proportional
to a’ . This was so that we could position the two major frequency bands in roughly
the correct places to agree with the experimentally determined density of states.
This led to or : fl : 'y : a’ = 1 : 0.3 : 0.1 : 0.03 and the overall magnitude of the force
constants was chosen to get the maximum frequency correct. As we will see later, the
maximum phonon frequency can,” is rather insensitive to composition and we will take
Lam” = 40 meV = 320 cm'1 . This frequency is independent of 7 and a’ and given by
tag,” = 8(a+ ,B)/3M so that a = 1.28 x 105 dynes/cm and fl = 0.38 x 105 dynes/cm,
where we have taken M = 72.6 amu as is appropriate to Ge; the masses of As and
Se are not significantly different as they are next to Ge in the periodic table. For
simplicity we therefor put all the masses equal. While these masses are irrelevant
for the elastic constants, they do influence other quantities that we calculate. Using
the ratios given above, the two weak force constants are 7 = 0.128 x 105 dynes/ cm
and 01’ = 0.038 x 105 dynes/ cm. We stress again that we are not seeking detailed
agreement with experiment, but are trying to get a reasonably accurate overview.
These parameters are now taken as fixed for the rest of this chapter. No separate
attempt was made to fit the elastic constants. The elastic moduli for the perfect

diamond cubic crystal in the absence of the weak forces are given by

011 = (0’ + 3fl)/4a

C44 = (am/[(0 + ma]
B = (011 + 20.2)/3 = (30 + m/m,

where I = a/x/3 is the nearest neighbor distance, which for Ge is 2.45A. Using the

15

values of a and 3 given above, we find that Cu = 1433 kbars, C44 = 696 kbars and
B = 830 kbars. These are in close enough agreement to the experimental values
C11 = 1289 kbars, C44 = 671 kbars and B = 752 kbars for our purposes.

From Fig.(2.1), we see that the weak forces stabilize the lattice as expected and
lead to elastic moduli that are non-zero for all (r). The most dramatic effects are seen
at low mean coordination where the singularity at (r) = 2.4 is completely washed out.
Note that we have used fl/a = 0.3 rather than 0.2 as in He and Thorpe. [8] One would
have hoped, with the addition of the weak forces, that some traces of the singular
behavior around (1') = 2.4 might have been retained, but this is not the case. Of
course, if we made the weak forces even smaller, then the singularity would start to
become apparent. However the weak forces are as nature has given and cannot be
altered. This Situation is analogous to any other phase transition in a small field that
couples to the order parameter; an example being the magnetization of a ferromagnet
in a small external magnetic field.

For (7') > 2.4, the network is mechanically stable and the elastic properties are
largely determined by the covalent forces a, 6. Other forces are small and can be
neglected. For (7') < 2.4, the network is mechanically unstable and we must include
other weaker forces to stabilize the struCture. In Figs. (2.2) and (2.3), we show a
comparison between these results and the experimental results of Yun et al.[11] taken
using ultrasonic data for the system GexS 61-x. Similar results have been obtained by
a number of other authors,[12] and initial discrepancies with the data of Halfpap and
Lindsay[13] have been resolved as being due to oxygen contamination. Also included
in these plots are the results for a-Ge.lt is not possible to make samples over the
range 2.6 < (r) < 4, as phase separation occurs. These results do show the softening
of both the longitudinal and transverse elastic constants of Se over Ge by more than
an order of magnitude. Also the experimental results agree with the theory in not
showing a singularity around 2.4. However the experimental variation of both C L and

CT over the range 2 < (r) < 2.6 is significantly less in the experiment[11] than in the

 

 

 

 

 

 

W l I ' l 7 l 1 1
- + Theory
- I Experiment.
’5 - +
g 1000 — "
v
.1 ,. *
D C
_ +
_ +
+
o M l 1 j l J J l l
2 3 4
<r>

Figure 2.2: The longitudinal elastic modulus C obtained by averaging the results

with the weak forces included, from Fig.2.1, is compared with experiment.

 

l7

 

 

 

1000 . . a . 1 a j
"' r + Theory
. . - I Expednunt
A ‘2‘:
E t b * .
:3 500 — *
V +
.- I- q
o .1 *
I- * d
P W I
0 IL!!! "'1" '-

 

 

J #1 J J _l J J
3

<r>

Figure 2.3: The transverse elastic modulus C obtained by averaging the results with

the weak forces included, from Fig.2.1, is compared with experiment.

18

theory. We do not have any good explanation for this except that our model is very
simple. It is hard to pinpoint the exact cause.

The model always maintains its cubic symmetry, on average, as the depletion takes
place, whereas the real systems are isotropic. We have therefore done a very simple
average over the principle directions to obtain the longitudinal CL and transverse CT

elastic moduli from the theory. [8] We set

CL = (6011,00 + 80,1,“ + 120,110) /26 (2.9)
and
CT = (60}00 + 80%“ + 60%" + 60532") /26 (2.10)

where C71 and Cm are the two different transverse moduli in the {110} direction. Ex-
pressing these quantities in terms of the cubic elastic moduli Cg; as given in Kittel,[14]
and using the relation for the bulk modulus BB = 2C11 + C12, we find the longitudinal
elastic constant,

0;, = (2701. + 518 + 68044l/78 (2.11)

and the transverse elastic constant

The results are not very sensitive to the precise way in which this averaging has been
done. For an isotropic system, in which C11 —— C12 = 2044 , Eqs. (2.11) and (2.12)

reduce to the result CL = CH and CT = C44 as expected.

2.4 DENSITY OF STATES

More complete information about these systems is contained in the phonon density
of states as a function of frequency. Such experiments have recently been completed
for 06,561-; glasses. Previously a density of states for a—Se had been obtained by

Gompf,[15] which showed a sharp low frequency peak at around w/wma, = 0.1 that

19

contains about 1/3 of the total weight. Here cam” is the top of the phonon band
which is quite insensitive to 2:, as we have remarked previously. Similar densities of
states are obtained for the crystalline forms of S e where these low frequency modes
can be associated with the torsional vibrations of the chains.[15] In the amorphous
material we prefer to think of these modes as the zero frequency modes that are
shifted to higher frequencies when additional terms are added to the potential given
in Eq. (2.7). We shall refer to these as floppy modes. Inelastic neutron scattering
measurements of the density of states in 061,581.;- glasses follow this low frequency
peak as a function of composition. The frequency of the peak position does not
shift much with 2:, but the weight in this peak rapidly decreases as 2: increases so
the the anomalous modes have vanished by about a: = 0.2 which corresponds to
(r) = 2.4. This is the behavior we would expect. The zero frequency modes produced
by the potential (2.7) are moved from zero frequency to a low frequency by the
additional weak terms like (2.8). However the weight in these floppy modes is still
given approximately by Eq. (2.6) and decreases from 1/3 for Se to 0 for (96035603 .
In order to better understand these experimental results, we have also employed the
bond depleted diamond lattice discussed in the previous section.

Using the equation of motion technique,[16] as described in Appendix A, for an
8 X 125 = 1000 atom sample with periodic boundary conditions, we have calculated
the density of states for various values of (r). The oz component of the displacement
for the ith atom u;a(t) is computed using the equation of motion that results from
the potential (2.7) or (2.8). The density of states p(w) can be written as a cosine

transform ,[16]

T
p(w) = 2] dt cos(z.ut)e“3(t/T)2 ZAgauiafl), (2.13)
71' 0

3,0:

where choosing the truncation time T determines the frequency resolution of the
result. The exponential factor in (2.13) helps with convergence. Here the appropriate

starting conditions for the density of states are,

“for = Aia = V5008 9m, (2.14)

 

20

where 0,0, is a random angle uniformly distributed in the range 0 < 0,0, < 271'. We also
set the velocities 0,0(0) = 0. We refer to this as “plucking”. Results are shown in
F ig.(2.4) for various values of (r). At small values of (r) the zero frequency modes are
broadened by the resolution function. It can even be seen that there are still a few
states with zero frequency beyond (1‘) = 2.4. This is because there is a very small,
but finite, probability of finding undercoordinated regions even when (r) > 2.4. This
is a Lifshitz-type[17] argument. It can be seen in Fig.(2.4) that the floppy modes are
moved into the low frequency peak and the weight is shuffled about within the rest of
the spectrum. The maximum frequency does not change significantly as (r) is varied
or when the weak forces in V are included.

In Fig.(2.5), we recalculate the spectra of Fig. (2.4), but using the “kick start”
method described in Appendix A, that leads directly and naturally to the density
of states divided by the frequency which accentuates the low frequency part of the
spectrum. The quantity p(w)/w can be written as a sine transform,

p(w)/w = i]: dt sin(wt)e—3(t/T)2 Z Agauga. (2.15)
The theory is compared to experiment in Fig.(2.5). The experimental results were
obtained using inelastic neutron scattering at a large wavevector transfer.[18] The
flOppy modes can be clearly seen in both the theory and the experiment-although they
are rather more pronounced in the theory. These experiments are quite difficult to do
at low frequencies. As the mean coordination (r) is increased, the floppy modes are
quenched out and vanish around (7') = 2.4. Because of the rounding of the transition
and the background phonons underneath the floppy modes, it is hard to be much more
quantitative than this. Nevertheless the agreement between the experiment and our
simple theory is impressive and convinces us that we have captured the essence of the
behavior.

In Fig.(2.6), we replot the data of Fig.(2.5), but divide by an additional factor of
the frequency. This is to accentuate the elastic behavior as well as the floppy modes.

Both the experimental and theoretical spectra should head towards a constant at very

 

21

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

<r>-2.18 I 1. (DP-2.24
z - ' 2 Q
N\ A
0' ‘ 0
0 0.5 1 0 0.5 1
l <r>-2.40 (r>-2,.4B
3 2 ' 2 -
E
0 0 L
0 0.5 1 0 0.5 1
<r>-2.60 <r>-z.66 fl
2 h P: z '- o.
o J \L 0 J L
0 0 5 1 0 0.5 1

Frequency 0

Figure 2.4: The density of states p(w) for various values of the mean coordination
(r). The dotted curve is without the weak forces, and the solid curve is with the weak
forces included. The full width of the frequency resolution function is Aw = Lama; / 20.

The frequency is in units where can,” = 1.

22

 

 

Theory Experiment
‘ r>=2.082 r>=2.08
r>=2.16 r>=2.16
r>=2.40 _ r>=2.40
r>=2.60 r>=2.80
’ r>=2.80 r>=2.80
a I
\ 20 20 - ’
A
3..
o.

 

 

 

//

 

 

 

 

 

 

0 ‘ i ‘ -
0 0.5 1 o 0.5 1

Frequency co

Figure 2.5: The density of states divided by the frequency, p(w)/w, for various values
of the mean coordination (r). The left hand set of graphs are theoretical and the
right hand set of graphs are from inelastic neutron scattering experiments.(Ref. [18])
The full width of the frequency resolution function for the theory is Aw = w...“ / 20.
The frequency is in units where 02m” = 1. cam” is the maximum 02 before the weak

forces are turned on.

p(w)/wz

 

23

 

 

 

 

 

 

 

 

 

 

Theory Experiment.
1000 1000 \g
<r>=z.0ea <r>=2.08
<r>=2.16 <r>=2.18
<r>=2.40 <r>=2.40
<r>=2.60 <r>=2.60
500 - <r>=2.80 500 - <r>=2.80
\ J“
‘ J's‘\
0 ' ' o P fl
0 0.2 0 0.2

Frequency on

Figure 2.6: The density of states divided by the frequency squared for various values
of the mean coordination (r). The data is derived from Fig.2.5. The left hand set
of graphs are theoretical and the right hand set of graphs are from inelastic neutron
scattering experiments. The full width of the frequency resolution function for the
theory is Aw = wmaz/20. The frequency is in units where wmar = 1. wma, is the

maximum of w before the weak forces are turned on.

24

low frequencies, although there are difficulties in both the theory and experiment in
getting accurate results in this part of the spectrum. The theoretical limit can be

found from the elastic constants via the integral,
(tn/.22 _. if: / 1' / 2" C-3(0 ¢)s' 0d6d¢ (216)
p Sr2 {:1 0 o i ’ 1n ’ '

where the three velocities of sound C;(0, (15) in the direction (0, qt) are found by diag-
onalizing the appropriate 3 x 3 matrix. We have evaluated this integral numerically
and the results are Shown in Fig. (2.7). A similar analysis has been performed by Yun
et al,[11] using their measured values of CL and CT and their results are also indicated
in Fig. (2.7). These provide a consistency check within both the experiments and
the theoretical model. We have had some difficulty getting good agreement between
the two theoretical approaches at low frequencies, despite strenuous efforts. Probably
larger models would be needed to take care of this point satisfactorily. However the
theoretical elastic constants are much larger than the experimental ones at small (7').
The agreement between the experiments and our simple model shows most of the
qualitative understanding of the floppy modes is contained within the general argu-
ments involving constraint counting etc. given in this chapter. We have taken care
in plotting the results this way as it is important to have the low frequency behavior

as accurately as possible in order to calculate the Debye Waller factors.

2.5 DEBYE WALLER FACTOR

The Debye-Waller factor, exp[—((k - u)2)], modifies the intensity in many scattering
type experiments. Here k is the momentum transferred to the sample. A measure-
ment of this quantity allows the mean square displacement to be extracted. The mean
square displacement in a particular direction for a single frequency is given by (23”),

where

(2:2) = __(,-, + é). (2.17)

25

 

 

 

 

' o 0 (A) Theory
1000 ~ + + (B) Theory
- I (C) Experiment
_ 00
6 n
+
500 -
. o
0 +
" o
’ I I I +.
_ 3 I4. I I I +
1 l l L 1 Y + m 0 l J
0
2 2.5 3
<r>

Figure 2.7: The coefficient a defined by p(w) = aw2 at low frequencies. The results
are (A) theory using the elastic moduli determined statically (B) theory using the
equation of motion method and (C) the elastic moduli from experiment. The units

are the same as those used in Figs. 2.4-2.6 with wmu = 1.

26

and ft is the Planck function.[19] At low temperatures, this becomes

h

2 — ——
and at the high temperatures classical limit,
kBT
2 = __
(x ) M002, (2.19)

where k3 is Boltzmann’s constant. For the many harmonic oscillators in the systems
under study here, we must average over all the eigenmodes of the system. We define

the moments
<w") = (”Meow (”Mm (2.20)

At low temperatures, we may extract the n = —1 moment and at high temperatures
we may extract the n = —2 moment. At both these two extremes the Debye-Waller
factor is very sensitive to the low frequency part of the spectrum. Indeed this is true
at all temperatures.

It is very dangerous to use the computations reported in Fig. (2.4) to calculate
these moments as the density of states in Fig.(2.4), with the weak forces included,
have a small finite value at w = 0 due to the resolution function and the cosine
factor in Eq. (2.13). We have computed the area under the curves in Figs.(2.5)
and (2.6) directly. However, it is better to compute this quantity directly using Eq.
(A.l5) in Appendix A. This is equivalent to integrating over Eq. (A.15). The result
can be obtained more directly using the single integral form in Eq. (A.l5), and we
have checked that a direct integration of the area under Fig.(2.5) leads to the same
result. The results for the r = —1 moments are shown in Fig.(2.8) for various values
of (1'). It can be seen that the vibrational amplitudes increase as (7') decreases, and
appears to Show a small discontinuity of Slope around (1') = 2.4. We should emphasize
that these results are calculated with the weak forces included. A similar procedure
can be applied directly to the experimental results shown in Fig.(2.5), to yield the

r = —1 moment and these results are shown in Fig. (2.8) also. Both these results

 

27

 

 

 

 

10
~ 0 Mossbauer (expt.)
- D + Neutron (expt.)
- 0 ' Theory
/\
3 ' D
o
E 5 - + , n
a + r V
b 0 D D D
\3/ g 3 +
.. . .
o
o l l l I I J J l
2 2.5
<r>
Figure 2.8: The r = -—1 moment (wmn / w) for various values of the mean coordination

(r). The figure shows a comparison of results obtained from Massbauer and inelastic

neutron scattering experiments with our theoretical results.

 

28

are compared with a direct measurement of this moment from the low temperature
Méssbauer data.[20]

We can obtain the r = —2 moment by directly integrating the area under the
graphs in Fig. (2.6). It is possible to write this moment as a single integral similar to
Eq. (A.15). However, this is a more difficult quantity to calculate than the r = —-1
moment as the single integral form is only conditionally convergent at large t, and
we found it unreliable to work with. We have therefore integrated directly the areas
under the curves in Fig. (2.6). This yielded fairly reliable results, but the results for
r = —2, shown in Fig. (2.9) for various values of (r), are somewhat less accurate
than those for r = —1 in Fig. (2.8). This same quantity was also obtained from the
neutron data[l8] shown in Fig. (2.6). These results are compared to the results from
the Méssbauer experiments[20] at high temperatures. The vibrational amplitudes
increase as (1') decreases, and possibly show a small discontinuity of slope around
(1') = 2.4. There is really no discontinuity in a thermodynamic sense in the results
for the 7‘ = —1 and r = —2 moments as the transition is washed out by the weak
forces. However the results of Figs. (2.8) and (2.9) do show more sharply defined
changes around (7') = 2.4 than the density of states themselves.

These results show that most of the contribution to the r = —1 moment comes from
the floppy modes in the low frequency peak making the low temperature measurement
of the Debye—Waller factor a sensitive probe of the floppy modes. The r = —2 moment
probes even lower frequencies and is sensitive to both the floppy modes and the
acoustic modes as can be seen from Fig. (2.6).

The preliminary analysis of Méssbauer results[20] seem to show behavior in qual-
itative agreement with these results in Figs. (2.8) and (2.9) and will be reported
on in more detail in other publications.[21] These experiments measure the Debye
Waller factors on trace amounts of Méssbauer active Sn atoms that are believed to
be tetrahedrally bonded into the network. Note that the masses and force constants

of the Sn are expected to be sufficiently similar to Ge and Se that (1”) would not be

29

 

 

 

 

_ . D Mossbeuer (expt.)
+ Neutron (expt.)
/\ ' . 0 Theory
N
A +
50 r-
3 D g :
E t D + +
a +
3 D ’ n
v '- 9 D D
V o
I o
O
O
O l L l J J l J— L
2 2.5
<r>
Figure 2.9: The r = —2 moment (win/w”) for various values of the mean coordi—

nation (r). The figure shows a comparison of results obtained from Mfissbauer and

inelastic neutron scattering experiments with our theoretical results.

 

30

expected to vary significantly from site to site. Our calculations involve an average
over all sites. It may be that regions of the sample that are denser in floppy modes
have larger (:32), but that remains to be seen. A more accurate analysis would also

have to take into account the heavier Sn mass.

2.6 CONCLUSIONS

In this chapter, we have shown that network glasses like GexS 61-, undergo a washed
out phase transition around a mean coordination of (r) = 2.4. The elastic moduli
Show no anomalous behavior around this point. The low frequency “floppy” modes
are quenched rapidly and disappear around this point. This behavior is reflected in
the mean square displacements as measured in Méssbauer experiments. The mean
square displacements increase as (r) is lowered and the rate of increase is larger for
(r) < 2.4. The results of experiment and a very simple theory are shown to be in

reasonable agreement, although there are quantitative differences.

 

Bibliography

[1] W. H. Zachariasen, J. Am. Chem. Soc. 54, 3841 (1932).

[2] J. C. Phillips, Phys. Stat. Sol.(b) 101, 473 (1980); J. C. Phillips, J. Non-Cryst.
Solids 34, 153 (1979); 43, 37 (1981).

[3] See for example, J. A. Erwin, A. C. Wright, J. Wong and R. N. Sinclair, J. Non
Cryst. Solids 51, 57 (1982) and references therein.

[4] F. Wooten, K. Winer and D. Weaire, Phys. Rev. Lett. 54, 1392 (1985) and
Wooten

[5] S. W. DeLeeuw, H. He and M. F. Thorpe, Sol. St. Comm. 56, 343 (1985); H. He,
Ph. D. thesis, Michigan State University (1985).

[6] P. Vashista, R. J. Kalia, G. A. Antonio and I. Ebbsj;, Phys. Rev. ”Lett. 62, 1651
(1989)

[7] M. F. Thorpe, J. Non-Cryst. Solids 57, 355 (1983).
[8] H. He and M. F. Thorpe, Phys. Rev. Lett. 54, 2107 (1985).

[9] See for example H. Goldstein, Classical Mechanics, 1950 (Addison-Wesley, Cam-
bridge).

[10] P. N. Keating, Phys. Rev. 145, 637 (1966)

31

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

32

S. S. Yun, H. Li, R. L. Cappelletti, P. Boolchand, and R. N. Enzweiler, Phys.
Rev. B39, 8702 (1989)

U. Tille, G. H. Frishat and K. J. Leers, in the Procedings of the fourth Inter-
national Conference on the Physics of Non-Crystalline Solids, Clausthal, West
Germany (Aedermannsdorf, Switzerland: Trans. Tech. SA 1977) page 631; R.
Ota, T. Yamate, N. Soga and M. Kunugi, J. Non Cryst. Solids 29, 67 (1978);
J. Y. Duquesne and G. Bellessa, J. de Phys. Paris C10, Supp 2, 46, 445 (1985):
ibid, 46, 449 (1985): J. Non Cryst. Solids 81, 319 (1986); K. Tanaka, Sol. St.
Comm. 60, 295 (1986)

B. L. Halfpap and S. M. Lindsay, Phys. Rev. Lett. 57, 847 (1986)

C. Kittel, Introduction to Solid State Physics, 1971 (J. Wiley and Sons, New
York) page 133

F. Gompf, J. Phys. Chem. 501., 42, 539 (1981)

See for example, M. F. Thorpe and S. W. DeLeeuw, Phys. Rev. B33, 8490 (1986)

and references contained therein.
1. M. Lifshitz, Advances in Physics 13, 483 (1964)

W. A. Kamitakahara, R. L. Cappelletti, P. Boolchand and F. Gompf, Bull.
Am. Phys. Soc. 32, 812 (1987); W. A. Kamitakahara et al ., Procedings of
the International Conference on Neutron Scattering, Grenoble, 1988 [Physica
B (to be published)]; W. A. Kamitakahara, P. Boolchand and R. L. Cappelletti
(unpublished).

See for example, N. W. Ashcroft and N. D. Mermin, Solid State Physics, 1976
(Saunders College, Philadelphia).

P. Boolchand and R. N. Enzweiler, Bull. Am. Phys. Soc. 32, 812 (1987); P.
Boolchand (unpublished)

 

33

[21] P. Boolchand and M. F. Thorpe, to be published. -

[22] M. F. Thorpe and Y. Cai in Disordered Systems and New Materials, Ed. by M.
Borissov, N. Kirov and A. Vavrek (World Scientific, Singapore, 1989) p. 81

 

Chapter 3

RIGIDITY OF INTERCALATED
LAYERED SOLIDS

This chapter is pulished in Physical Review B. Reference: Y. Cai, J. S. Chung, M. F.
Thorpe and S. D. Mahanti, Phys. Rev. B42, 8827(1990)

3.1 INTRODUCTION

All crystalline solid solutions show a composition dependence of the average unit cell
volume (V) which increases with the concentration of the largest constituent.[1] Inter—
calated multilayer systems such as A1-xB,L are an important class of such systems.
The linear variation of (V) with a: is the well known Vegard’s law,[2] although solid
solutions often exhibit a complex nonlinear (superlinear, sublinear, sigmoidal) behav—
ior. The value of (V) depends, at the microscopic level, on competition between the
local and global energies associated with forming the solid solution. These energies
depend upon the relative sizes and compressibilities of the intercalated atomic species
and on the overall rigidity of the system.

There is a large variety of ternary layered intercalated compounds of the form

A1_xBxL [A (B) denote small (large) intercalants and L denotes layer] which can be

34

35

treated theoretically in a relatively simple way, because the major expansion takes
place in the direction perpendicular the layers, denoted as the c-axis. [3]-[8] However
correlations between the positions of the intercalants in the multilayer case play an
important role in determining the gallery structure. These correlations arise from the
interactions between the intercalants mediated by the host layers.[9] For example,
similar-sized intercalants prefer to stay together in the same layer to minimize the
layer distortion energy. In contrast, similar—sized intercalants in adjacent galleries,
tend to repel one another.[9] To a large extent, how important a role this correlation
plays in the multilayer system, depends on the method of sample preparation. Rapid
quenching leads to a more random distribution of intercalants, whereas annealing
allows for interlayer correlations (ILC) between the intercalants to develop which
may ultimately lead to segregation and/or phase separation. Even in the rapidly
quenched samples a finite degree of short range correlation may be present.

In the past, several attempts[4],[9]-[11] have been made to study the nonlinear
:c-dependence of the c—axis expansion for a bilayer (two layers ‘sandwiching’ a single
layer of intercalants). In this case, the most recent study,[12] referred to henceforth as
I, rectified shortcomings in previous work by considering both the transverse rigidity
of the layers, the differential size, and the compressibilities of the intercalants. The
multilayer system we will discuss in this paper is a generalization of such a bilayer
system and is more appropriate for a real materials. It contains all the ingredients of
the previously studied bilayer model in I and in addition takes into account the effect
of interlayer correlations between the intercalants in different layers. It is assumed
that correlations between intercalants in the same layer are much less important. Such
a multilayer system, with a single kind of intercalant is shown in Fig.(3.1), where the
circles could represent K atoms intercalated between the C atoms in graphite. If a
few larger, say Cs, atoms are present, then the layers begin to buckle as shown in
Fig. (3.2). This buckling of the C layers increases with the concentration, and is a

maximum for the 50-50 random alloy, as shown in Fig. (3.3). A similar configuration

36

 

oxexoxexoxoxoxe

”nu/11011111111111

01.101.101.101.

IIIIIIIIIIIIIIIII”III

01.101.101.101.

IIIIIIIIIIIIIIIIIIIIIII

01.101.101.101.

nun/10017001111

01.101.101.101.

”unlnlllnnllllll

01.101.101.101.

[Inn/1111100157011

exoxexexoxexexe

 

 

Figure 3.1: The circles represent a single species of atom, intercalated between host

layers of graphite.

37

 

QQmmmth
‘ IOOOOOC‘)

' ////////Z/////J//fffJ/I

©0000...

Illllllnl

ammflmmmm
QOOOOOOO
0000000

//]/////////]///7/////

00000000

////////////////////// I

 

 

 

Figure 3.2: A few larger atoms are present in addition to the small atoms shown in
Fig.(3.1). This drawing is produced from an actual computation for a 2d multilayer

system as described in the text. Here I: = kA = k3 = k1.

38

 

 

01.101 O 01.

I’lllll’ f/7 / // // /

01.1.10 00

[Illa/”zuglllnull

OIOWIO
(I'll m010m1fi'310

”Ipllllllllllllllllll

01.101.11.101. 0

1127”” l’llllfllllllll

010101010101. O

.inll'tlllnllllinnl

0,1010fl0101010

 

 

Figure 3.3: Equal number of large and small atoms are intercalated between the layers.
The buckling of the layers is a maximum in this case. This drawing is produced from

an actual computation for a 2d multilayer system as described in the text. Here

k=kA=kB=kT.

39

 

O

///////

“W

W

.W. 0.

////////:

OIIIOIOIOIIIOEO

Figure 3.4: The same as Fig.(3.3), but with thicker layers, appropriate to silicates,
rather than the graphite layers shown in Fig.(3.1)—(3.3). This drawing is produced

from an actual computation for a 2d multilayer system as described in the text. Here

k=kA=kB=kT.

40

is shown in Fig.(3.4), but with thicker layers more appropriate to silicates (clays).
Figs.(3.2)-(3.4) are drawn from the results of actual simulations and are not just
sketches. A detailed discussion of the nature of this buckling and how it depends
upon various physical parameters such as the compressibilities of the intercalants and
the rigidity of the host galleries is given later in the chapter. Mori et.al[13] were first
to point out the existence of static displacement and charge density modulation of
the adjacent graphite layers when the disordered unregistered potassium layers order
at low temperatures in stage-II potassium-graphite.

The arrangement of the chapter is as follows. In sec. 3.2, we introduce the har-
monic spring model describing these systems. In sec. 3.3, we present the results of
an exact calculation for the case k4 = k3, where kA and k3 are the local compress-
ibilities of the intercalants. In sec. 3.4, we extend these results to take into account
correlations perpendicular to the layers. In sec. 3.5, we develop an effective medium
theory for the more general situation when 19,; # k3 . In sec. 3.6, we discuss the
complete probability distributions for the various gallery heights (or equivalently the
interlayer distances) and the effect of ILC on these distribution functions. Compari-
son with computer simulation results are presented where appropriate throughout this
chapter. Various exact mathematical aspects of the complete distribution functions

are discussed in the appendices.

3.2 THE MODEL

Consider a layered ternary system with composition A1_szL where L represents
the host layer such as graphite, dichalcogenide or vermiculite.[10] A and B are two
different types of intercalants which are assumed to occupy a set of well defined lat-
tice sites. This assumption is reasonable for most of these intercalation compounds
because of the strong interaction between the alloy (intercalants) and the substrate

(host layers). For example it is known that in ternary alkali intercalation compounds

41

M1-3M;:Cg (M, M’ = K, Rb, Cs), the alkali atoms are located in registry over car-
bon hexagons forming a (2 X 2) RO" in-plane structure.[7] For a fixed (i.e. quenched)
distribution of intercalants, the total energy of the system which controls the struc-
tural distortions, can be approximated by a sum of two major contributions; one
associated with the interaction between the intercalants and the host and the other

between host atoms themselves.

1 2 1
E = .2— : k, (z, .. 2,-” - h9_,,) + 51% 2 (.2i — z,+,)2 (3.1)
i (L6)

The first term is the host-intercalant interaction approximated by a harmonic spring
of strength k,(= [CA or kg) and equilibrium height h?(= he, or kg). The second term
is the interaction between host atoms. We only include the transverse part kT of
the layer rigidity which is sufficient for the present purpose (in general one has to
also include the flexural rigidity[12],[14] and other terms). Our notation follows I. A
diagram of the spring system associated with Eq. (3.1) is given in I for the bilayer
case. The generalization to the multilayer case is obvious.

In Eq. (3.1), i is a vector labelling the position of the ith spring and zi is the
z-component of the coordinate of the ith vertex. Here a vertex refers to a point on
the host layer which is directly above the ith intercalant and hi = 25 — 2i-” is the
local gallery height at the site of the ith intercalant. Only the z-component is needed
because the motion of the atoms is constrained to be perpendicular to the planes
in this model. The unit vector 1/ is along the c-axis and 6 is a vector in the plane

perpendicular to the c-axis connecting adjacent lattice points.

We define the average gallery heights (h), (hA), and (its) as
I
h = — h° . .

Here N is the number of galleries. Site specific averages (h A) and (113) are given by
similar equations except that i only runs over sites occupied by either A or B . Notice

that the total height of the sample is proportional to (h) as is the area or volume

42

(V) of the sample, as there are no length changes parallel to the planes. Thus all
definitions of length, both microscopic and macroscopic, are equivalent.

One can define a dimensionless height (11 through the relation:
hi = ha + (“his - hi) , (33)
where di = 0 or 1, for A or B respectively when there is no distortion. The average
gallery height ((1) is given by:
(d) = "11;; di , (3.4)

with similar definitions for (dA) and ((13). Note that by definition:

(d) = (1 - xdi) + $(d3)- (3-5)

3.3 ANALYTIC SOLUTION

When the local compressibilities are equal, we have kA = k3 = k, and the energy in

Eq. (3.1) becomes

1 1
= 5k: (zi—z.-.— hi..)”+ rim—2m)? W»
i (L6)

Unlike the treatment for the bilayer system in I, it has not been possible to express
the energy functional (6) as proportional to (hf?1 — h%)2 at the outset,-although this
result is obtained later on. It is quite difficult to deal with the above form for E,
because the coordinate z is cumulative. We therefore introduce another variable ui
which is defined as:

ui = Zi — h" (i - u), (3.7)

where the parameter h“ is an average length for each spring and is determined later
by minimizing the distortion energy; ui can be interpreted as the displacement of the
ith vertex from the average regular lattice. In terms of this new variable, the energy

is given by

1
E = 4:: (ui - Ui—u + h' - h?_,,)2 + 516T Z (“i - Ui+5)2-
i (ms)

43

= E}; + EkT . (3.8)

where E; and Eh, represent the first and second terms above.
For a fixed 12‘, the equilibrium conditions are 8E/8zi = 0 or aE/aui = 0, which
lead to:

k (221; — ui_,, — uiw) + k7 Z (ui — u5+5) = k (111,, — h?) (3.9)
6
Eq. (3.9) can be rewritten in matrix form as:

MU = (I), (3.10)

 

where U is the displacement vector of the system and the matrix M is formed from the
coefficients. The random vector (I) contains all the disorder. Since M is a nonrandom
matrix, we can solve the problem analytically. (Note that this is possible only for
(CA = 133.) A similar form was given in I [see Eq. (16)].

Denoting M'1 as the inverse of M, we obtain:
u, = kZ(M“)U (1131,, — hf) . (3.11)
.i

In the expression for the energy E in Eq. (3.8), only the first term contains h“.
Substituting for ui from Eq. (3.11) into Ek, we obtain

k .. ..
= 5 :{w >2 - 2h ha)... + mm
1

+2Z[<M ”(M—1h a] (123.. — h°)(h — hi.)
+k2 Z [(M‘x. —<M-‘). m] [(M’lm — (M—lli—ul] (hgt. - hs’xhi’... — hm .
(3.12)
We introduce a new variable 0i replacing h? such that
ho , if 0° = 1
(1?: (101+ b = A l (3.13)
h%, if 0i = -I.

44

It is very convenient to use this language, borrowed from magnetism [1-5] to describe

the locations of the intercalants. The two parameters a and b are:

_(h° - ho ____)_B
2 (3.14)
and 0 0
b = (Jig—hi) . (3.15)

Replacing h? in terms of 01 in the expression for E1, and performing the ensemble

average, we find
E1. = :21 (If)? — 2h*<h> + (12)
4k); [(M“)-(; — (M-1)I—Vj] 'az((0j-v0i-v) - (Um—u»
+12 é: [(M-‘m — (Mm-..) [(M-‘m — (M-1)1-.1]
«Quorum-» + (0101) — (0.-.,00— man-» >)} (.3 (16)

Minimization of the energy with respect to h", i e

8E _6Ek

am‘ ah, - —0 - (3.17)

 

 

gives Vegard’s law:
h'“ = (h?) = h0A(1—x)+h%x= a(1—2x)+b, (3.18)

where we have made use of the relation (0 -)- — 1 -— 2:13.

The expression for the energy contains various two spin correlation functions. For
the case of a random distribution of intercalants, the two spins are uncorrelated
(corresponding to two intercalants either in the same gallery, or in adjacent galleries

being uncorrelated) and we have

(Uioj) = (I — 2.1:)2 + 4:1:(1 — x)6ij . (3.19)

 

45

After some mathematical manipulation, we obtain simple expressions for the dimen-
sionless energy per site 6 = E / (4N azk), and also for the dimensionless lengths (d A),

(d3) and the fluctuations:

e = gen — :c) (1 — W(k)) , (3.20)
MA) = x (1 — W00) 9 (321)
(d3) = 1 - (1 " x) (1 — W00) 1 (3-22)

(d2) — (d1)? = x0 — 22) (W102) — WW)

= ((123) — (d3)2 . (3.23)

These equations are identical to those for the bilayer case with functions W(k) and
W1(k) being similar to quantities used in I, but now appropriately modified for the
multilayer case. As in I, we refer to all quantities like W(k) and W1(lc), as gener-
alised Watson integrals.[12] These Watson integrals reflect the influence of the lattice

structure on the various quantities calculated. They are given respectively by

WU?) = 2(M"1)11 - (M—1)i+ui - (M_1)i-Vi

 

 

= 111?; 2H1 — (fish. 11)} , (3.24)
and
Wl(k) = Z[2(M-1)U — (M-l)i+uj — (M_l)i-Vj]2
J
= if; 4k [1 - :2s(q ~ u)] , (3.25)

where the dispersion relation for the pure system is given by

Aq = 2k[1 — cos(q - u)] + sz(1 — 7(1) (3.26)

 

46

and -
1 .
7q = Z E e“. (3.27)
6

In Eq. (3.27), z is the coordination number in the plane. Note that we have made

use of the expanded form of (M’1)fi:

e—iq- r-U

 

1
)n _ K,- g6)“ (3.28)
From Eqs. (3.5), (3.21) and (3.22), we find that the average normalized interlayer
spacing (d) = a: and therefore we see that Vegard’s law is obeyed. It is remarkable
that Vegard’s law is independent of the number of layers and the only way one can
see the effect of multilayers is in the site-selective gallery heights (dA) and (d3) and
gallery height fluctuations through the generalized Watson integrals. In the bilayer
case the factor [1 —- cos(q - 11)] in Eqs. (3.24) - (3.26) is replaced by one and the q,
summation also gives unity. We have found that for a given strength of layer rigidity
(k7), W(k) for the multilayer is smaller than for the bilayer indicating increased
stiffness. The exact results given in Eqs. (3.20)—(3.23) are shown in Fig. (3.5)-(3.7) in
the center panels marked R (for Random). A discussion of these results is given in the
next section in conjunction with the results for the correlated case. The diagrams for
Figs. (3.2)-(3.4) (and also for Figs. (3.8) and (3.9), were constructed from computer
simulations and show a small piece of a much larger sample. The intercalant atoms are
centered midway between the lattice positions i and i + V. The thickness of the host
layers, shown as pairs of parallel lines, is chosen to roughly correspond to graphite
in Figs. (3.1)-(3.3) (and Figs. (3.8) and (3.9)), while the thicker layers in Fig. (3.4)
correspond to the silicate layers in clays. The intercalant atoms have fixed radii, one
for the A species and another for the B species. The variation in gallery heights in
exaggerated somewhat in these figures for easier visualization. This is not a serious

deficiency as the heights all scale with (11% — ha).

 

 

 

 

47

 

 

 

 

 

 

 

 

 

 

 

 

' I ‘ ' l T

1 A 7 R 7
a: I I
s - - - r-
z - «- - d
E}
m .. .1 1. d
5 a

. 1 - J

F: I /
<1 0‘4 l 1 —- 4 1 l . -—4

O 05 1 0 06 1

Figure 3.5: Showing the average heights of the A, B galleries (dA), (dB) and the
overall average height (d), in dimensionless units. The symbols are from computer
simulations in a 2d multilayer system and the solid lines are exact results. The vertical
lines represent the full width at half the height for the partial length distribution
functions, and are not error bars.The symbols A, R, and C stand for Anticlustering, -

Random and Clustering respectively. Here k = kA = k3 = kT.

 

48

 

 

 

 

 

 

 

 

 

 

 

 

I I I I l I 1 I I
0.1—A J —R — -C -—o.1
- .l .. .. .. ..
~ .1 1.. .- .- cl
1- .l .. .. .- a
r- J .- d .— .l
00 l l 1 1 l l 1 l 1 00
0 05 1 O 0.5 1 O 05 1

Figure 3.6: The strain energy 6, is shown for the 2d multilayer case. The symbols
are from computer simulations and the solid lines are exact results. The symbols A,
R, C stand for Anticlustering, Random and Clustering respectively. Here k = kA =

k3 = 1:1.

49

 

 

 

 

 

 

 

 

 

 

 

 

 

1 i T .1 l l 1 .1 3 I I .-

3
H 0.2 L— A - -- _ _. _
5* R C

/\ . -

4

"C

Y - ..

/\
51‘: ~ . ’ -

'U

V -
V

0.0 ' J ‘
o 0.5 1 _ .
X X X

 

Figure 3.7: The width, \/ (d1) — (dA)2, of the distribution function for the A gallery
height for the 2d multilayer case. The symbols are from computer simulations and the

solid lines are exact results. The symbols A, R, C stand for Anticlustering, Random

and Clustering respectively. Here k = kA = k3 = kT.

50

 

W0
01010191010101.
O O O '

Illnllllllllflnlllll'

01.1.IOIOI.I.IO

IIIII

C )C ) £5“de

I'IIIHIIIIIII

01016161010101.

”IIII'IIIIII”

0101010101.

 

 

 

 

 

Figure 3.8: Equal number of large and small atoms are intercalated between the
layers with Anticlustering perpendicular to the layers. The drawing is produced

from a actual computation for a 2d multilayer system as described in the text. Here

k=kA=kB=kT.

51

 

r fl///////

0 o 1.319}.

jJ//// [III

0 O I‘fg’fo'IO

///J//J :IIIIInlllIII

)0 O ICE-IO

Cj/j/ij/ IIIIIHIIIIII

o (M 193.10-

0 OK IQECI.

IlnlInllllIl

MCIOEOIO

IIIIIIIIIIIIII

o o % 10mm.

JJJ/J/jf

 

4

 

Figure 3.9: Equal number of large and small atoms are intercalated between the .

layers with Clustering perpendicular to the layers. The drawing is produced from

a actual computation for a 2d multilayer system as described in the text.

k=kA=kB=kT.

52

3.4 CORRELATED CASE ’

When there is finite correlation between intercalants occupying different galleries we
have to know precisely how these are correlated before we can obtain the intergallery
structure by minimizing the energy for a given configuration of the intercalants. For
the sake of illustration we show in Fig.(3.8) a situation where the larger intercalants
in adjacent galleries tend to avoid each other in order to reduce the build up of
local strain energy. In our formulation of the problem the correlation between the
intercalants appears through the Ising model spin-spin correlation functions which
must be known in order to describe the positions of the intercalants in the galleries..
For simplicity we use a 1d Ising model [15] to allow for correlation along the z—axis
only, i.e., include only interlayer correlation and neglect intralayer correlation. It is
believed that this corresponds to the actual situation in many cases. This is fortunate
as this case can be handled in a rather straightforward way.

The interlayer correlation function can be obtained from a 1d Ising spin model [15]

with the Hamiltonian given by
H = —J 2030,44 — h: 0'; . (3.29)

where the magnetic field h controls the value of (0;). Using the standard transfer

matrix method,[15] we have
(Uin) — (U{)<0'j> = 4:1:(1 — 3:)Ali—jl , (3.30)

where (0;) = 1 — 2:1: and is independent of i. In Eq. (3.30), we have used

 

_ \/4:1:(1 — x) + (1 — 2:021” — P

A — ,
\/4$(1— x) + (1 — 21:)2102 + P

 

(3.31)
where
P = e-ZJ/kBT (3.32)

and J is the nearest neighbor interaction, k3 the Boltzmann constant, and T the

temperature. The degree of interlayer correlation is controlled by the parameter

 

53

P. When J = O, i.e. in the random distribution limit, P = 1 and‘A = 0. For
ferromagnetic interactions J > 0 and P < 1, which implies that two similar ions
(large-large and small-small) prefer to stay together. Usually in intercalated systems,
similar ions in adjacent layers would be expected to repel each other, because of steric
effects, which implies an antiferromagnetic interaction, i.e. J < 0 and P > 1. Note
that in the extreme clustering (ferromagnetic) limit, P -+ 0 and A —) 1, while in the
extreme anticlustering (antiferromagnetic) limit, P —+ co and A —+ —:r/ (1 — :c).

For the fluctuations in the heights for the galleries containing the A and B inter-
calants, one needs three-spin correlation functions, which can also be calculated using

transfer matrix techniques.[15] The result is

(mag-(7k) - (0.)(aj)(ak)

= 4:2:(1 _ 3:)(1 _ 2x) (Ali-11+ Ali-kl + Alk-il _ 2A’§(li—jl+li-kl+lk-i|)) . (333)
which can be rewritten compactly as
(Uiajak) = (0) ((02%) + (Ujak) — (01:00) a (3-34)

for i _<_ j S k. These results, Eqs. (3.30) and (3.33), can be utilized in a multi-
dimensional network with correlations only in one dimension as follows; the two and
three spin correlation functions are given by Eqs. (3.30) and (3.33) if all spins are in
the same column, and are zero otherwise.

In the correlated case we also have arrived at the same simple forms for (dA), (dB),
and 5, as in Eqs. (3.20)-(3.22), except that the Watson integral W(k) is replaced by
Wc(k) where

 

 

Wows) = fig: ”“1 ' if” ' ””ch, q) ' (3.35)
and
Fennel) = 1" A (336)

1+ A2 — 2A cos(q. V)

 

54

There are two points worth mentioning here. First, just as we obtained Vegard’s
law behavior in the random case, we continue to get a linear relationship between
(d) and :1: even when there is finite ILC. This shows that the Vegard’s law is
independent of the strength of the ILC as long as kA = 1:3. This is not
surprising as Vegard’s law is obeyed in the most extreme case of ILC, i.e. complete
phase separation. Second, the simple forms of the strain energy and other physical
quantities imply that the effects of correlation can be incorporated into the generalized
Watson integral. This is significant because it provides a way to construct an Effective
Medium Theory (EMT) in the presence of the ILC. This point will be clarified after
we discuss the important role played by the Watson integrals in EMT in sec. 3.5.
The fluctuations in the A and B heights are more difficult to obtain because 3-spin
correlation functions are involved. This requires the introduction of another integral

besides modifying W1 (see Eq. (3.25)) to W01. For example,

 

 

<d31> — <d1>2 =x<1—x>(wc. —wcz)+x'~’(wcz—Wé) , (3.37)
where
1 41:2 1 — cos 1/
W01(k = )‘NZ [ ”(q )1 F a(A,q) (3.38)
q
and we need in addition
2 - . — I. . _ l i 0.. o- -'
W020.) z _1__ k [2 2cos<q 13112 2cos<q m 26.3%... mph/(1,11,11,11
N 3 q’q, Aqu’ ij,1

1 k2[2 -— 2 cos(q - u)][2 — 2 cos(q’ - V)]

 

 

 

 

 

 

— F02(q, q', A) . (3.39)
where
uv + A 1 + no uv + A 1 l
F02(q,q,A)= uv—A+A(1—uv uv—A)i(u—A+v—A) (3.40)
in which n = 6"” and v = e'q .In the random distribution limit A = 0, and all

these results reduce to those of sec. 3.3.
We have carried out a series of computer simulation for 2d multilayer systems where

z = 2. The results are shown in Figs. (3.5)-(3.7). Fig. (3.5) shows the dimensionless

 

55

heights (d), (d A), and (d3) versus the composition :3, for different amounts of correla-
tion, ranging from clustering through random to anticlustering. In Figs. (3.5)-(3.7),
the symbols AR and C, stand for Anticlustering (or antiferromagnetic), Random
and Clustering (or ferromagnetic) respectively . The corresponding parameters P
[see Eq. (3.32)] are P = 500, P = 1 and P = 0.1 respectively. In contrast to the
average interlayer separation the (d), ((1,1) and ((13) show nonlinear x-dependence
particularly for the clustering situation when the concentration of one of the inter-
calants is very small. Our theoretical study of the multilayer system clearly points

out the need for probing the x-dependence of site-selective local gallery heights in

 

layered intercalation compounds, which could be done using EXAFS.[16] Fig. (3.6)
shows the strain or deformation energy as function of composition at. This energy is
greatest for the case of clustering as the packing is most difficult to achieve in that case
because the height mismatch between adjacent columns would tend to accumulate .
This would lead to a very great strain energy and so in response to this situation,
the layers become quite flat as shown in Fig. (3.9). Indeed as P becomes infinite, the
layers are absolutely flat and all the strain energy is in the first term, involving E,
in Eq. (3.8). In the other limit of anticlustering, the strain energy is reduced as the
packing becomes more comfortable as shown in Fig. (3.8). Here the strain energy is
shared between the two terms in Eq. (3.8). Note that the strain energy is particularly
favorable when x = 1/ 2, which is the case illustrated in Fig. (3.8), when there is a
narrowing cusp in the energy as seen in Fig. (3.6). The width of the distribution also

shows a cusp in this case as shown in Fig. (3.7).

3.5 EFFECTIVE MEDIUM THEORY

3.5.1 Random Intercalation

When the two intercalants have different compressibilities kA = 1:3, we cannot obtain

an exact solution. In this case one has to develop approximate methods such as

56

effective medium theory. The effective medium theory for the random multilayer case
is a straightforward generalization of I which follows the basic ideas of Feng, Thorpe
and Garboczi.[17] The effective medium has uniform spring constants k, kT and spring
length he. The effect of the environment on a single intercalant is contained within an
effective spring constant k... When a force F is applied locally to expand the gallery
height, the resulting displacement uo equals the force divided by the effective spring
constant Ice 5 k/W(k). F

F
uo = I): = WOC)? . (3.41)

The function W(k) is the same Watson integral we have discussed in the previous
section and in I. It has a clear physical meaning here, i.e. the rest of the lattice is
characterized by another spring with compressibility k; = 13C — k and natural length
he, and the effect of this environment is contained in the Watson integral.

Let’s get back to the construction of the EMT. Replace one of the springs (k, he}
in the effective medium by {k2,hg} (a = A or B). The deformation energy then

takes the form

5 = £13.01 — 1.3,)2 + $140. — he)2 , (3.42)

where the first term is from the deformation of the a-spring when its length is changed
from hg to h, and the second term is from the the effective medium representing the
rest of the material. By minimizing e we obtain the local equilibrium length h = ha,
where
(13’ he + ho, ho)
ha, 2 3 °‘ . .
(k; + ka) (3 43)

Substituting this back in Eq. (3.42) gives the energy for a single a-impurity in the

 

effective medium,
(he - hgl2kéka
2(k; + k0,)

There are two different types of impurity springs and the total deformation energy

 

ea = (3.44)

is the sum of two types of deformation energy weighted with probabilities pa. Here

57

pA = 1 — :c and p3 = 3:. Therefore, we have

_ 1 (he - hg)2kéka
6 — 2 211,. (k; + kc) . (3.45)

a

 

Minimizing c with respect to he determines he, i.e.

h = 2.p.hgk./<k;+ka)
c 2a Parka/(k; + kc) .

We need another equation to obtain the two unknown quantities, k and he and use

 

(3.46)

the following argument. Apply a force F at the single impurity site. This produces a
displacement F/(k; + k0,) which when averaged over all sites is set equal to F/k, to

give the self consistency condition,

 

pa _ i
)3 k; + k... _ k. . (3.47)

a

Eqs. (3.43) and (3.47) solve the problem. The results are:

 

 

 

(d) = :1: + :z:(1 — 1:)Fd (3.48)
Cckek’kB
d = c 3.49
l A) 1414+ k4)(ké + k3) ( l
(l - x)kek;kA
d = 1 — 3.50
< B) We; + mm. + k3) ( l
l
E = -2'\/kAkB 23(1— :L‘)(h% — (£30217; , (3.51)
where
kck’UcB — kA)
F = c 3.52
d WC; + (“AWL + k3) ( )
and
FM}: k
F. = —“'—A—‘?- . (3.53)
kg — (CA

The fluctuations of h .4 and hB can be found by applying the Feynman-Hellman
theorem.[12],[18] The results have the same form as those obtained in I:

«(114 - (dA))2) — 3(1- xlilkAkBk3]/[k(k; + 19020,; + kB)]}2

“ Wz/(w. - W) + ((1... 433- (can/[em knew 43)] (3'54)

 

58

«d — W) = “"1 " “{[kAkBkeI/Wk; + k1)(k;+ legume/(W. — W?)
W"/ W - W?) + ((1. — was — ken/[(14 + 1:..ka 1.3)]

 

(3.55

The EMT is extremely good in describing the random case. Comparison witl
computer simulations are shown in Figs. (3.10) - (3.12). The EMT results are show:
as dashed lines, to contrast with the exact results shown in Figs. (3.5) - (3.7) where
the exact results were shown as solid lines. We see that the average height (d) m
longer obeys Vegard’s law, but is superlinear when 1:33 > kA. When k3 < kA, th!
average height (d) is sublinear. A sigmoidal shape is never found in this model. Th:
very stable configuration in the anticlustering case is also obtained when kA 31$ k3 am
2: = 1 / 2. The EMT contains a cusp in the solution at this point, as does the exac
solution obtained when kA = k3. We have not calculated the widths in Fig. (3.12
for the correlated cases as this would involve very messy summations over the three
spin correlation functions. It should be emphasized that the EMT is not exact am
although the overall agreement with simulation is excellent, there may be situation:

where this is not so.

3.5.2 Correlated System

In the random case, the Watson integral is naturally introduced into the EMT t1
account for the effect of the environment.[17] As we have discussed before, the eas:
results in the limiting case kA = k3 show that we can incorporate the ILC into the
effective rigidity of the multilayer system by using a suitably modified Watson integral
Whenever the problem is analytically soluble, i.e. when the matrix M is non-random
we have arrived at expressions for the energy and various average heights. The fac
that all the background effects, including correlations, are contained in the Watsm
integral suggests that Wg(k) is an intrinsic function of the lattice network and it ca]
be used to describe physical quantities which are related to the first moments of th

height distributions, such as the energy and the average height, because these depent

59

 

 

 

 

 

 

 

 

 

 

AVERAGE LENGTH
\\
”4:“
+—
. \.
3‘

 

 

 

Figure 3.10: Showing the average heights of the A, B galleries (dA), (dB) and t]
overall average height (d) in dimensionless units. The symbols are from comput
simulations in a 2d multilayer system and the dashed lines are from EMT. The vertic
lines represent the full width at half the height for the partial length distributic

functions, and are not error bars. The symbols A, R, and C stand for Anticlusterir.

Random and Clustering respectively. Here hp = 19,, and k3 = 5kA.

60

 

 

 

 

 

 

 

 

 

 

 

002 I 5 j I T I o
. I l 1 ° 2
r- “ I. i ‘X T
.- A J .. R d h- p . C -
.. .l .. - ‘ ‘\ -l
j I \
_ .. .1 l. ' ‘ 1
I \
0.1 ‘— —‘ '— .“ -— ”—1 \‘ "" 0.1
.- J L Q, .\ .1 l— T ‘ J
, x ; a
x"..’-'~ ‘ ‘ f x ' ~. ~ ‘
1.. I \ .( r. ’ ‘ - -" § -*
f ‘0 I. h ' ‘
h’ \J r-I \“ l, ‘c
I l \ r l \ I 1
0.0 ' ' ' ‘ 4 ‘ 0 0
0 0 5 1 0 0.5 1 0 0.5 1

Figure 3.11: The strain energy 6, is shown for the 2d multilayer case. The symbols
are from computer simulations and the dashed lines are from EMT. The symbols A,
R, C stand for Anticlustering, Random and Clustering respectively. Here kT 2' Is),
and k3 = 51%-

 

 

61

 

 

 

 

 

 

 

 

 

 

 

 

 

_- l j d ’- 1 I I q k l I I _l‘
N
> 0.2 h. A -"' —' "T.\ R -—I — O C .u-J
p ’. 3‘ . Q .
A< r- . . cl .— ” “ .1 l- . u
'5 o ’ f ‘, o 0
V .. . O . _I .\ _( Q
. y \ " d
I o . O , \ 00.
NA "’ ' ‘\
6‘ ’ ' 7' ‘5 ’ "
o ' 1
\V/ ’ ' f' j - .
0.0 i l l 1 l l l I 1
0 0.6 1 0 0.5 1 0 0 5 1

 

Figure 3.12: The width, (fldi) — (dA)2, of the distribution function for the A gallery
height for the 2d multilayer case. The symbols are from computer simulations and the

dashed lines are from EMT. The symbols A, R, C stand for Anticlustering, Random

and Clustering respectively. hr = Is), and k3 = 5kA.

62

only on two-spin correlations. It then becomes straightforward to construct an EMT
for the correlated system and arguments used for the random case above can be
extended. We replace W(k) everywhere by Wc(k), which is assumed to characterize
the effective medium when correlation is present. The new EMT is found to be
surprisingly good when we compare its results with computer simulation as shown
in Figs. (3.10) -(3.12). with P gé 1. Indeed the comparisons with simulations when
correlations are present, are about as good as for the random case. Note again that
we have not carried through the complicated summations over three spin correlation
functions that would be necessary to generate results for the correlated cases in Fig.

(3.12). In the random case, where A = O, we have the limits,

Wc(k) = W(k)

WC](k) = W1(k)
z _k26(W(k)/k)
8k

WCzUC) = W(k)2 . (3.56)

In both the random case and the correlated cases, there is no linear relation between
(d) and a: in general as demonstrated in Fig. (3.10)..

The EMT for a correlated system can be derived also like that for a random one.
The interpretation of Eq. (3.41) must be modified. Instead of applying a pair of forces
:le to the ends of a single bond, it is necessary to apply these forces to every A
gallary. To avoid a net pressure on the body, forces :FFx/ (1 — 2:) must also be applied
to the ends of all the B galleries. That is, identical forces are applied to the complete
set of the A bonds and another set to the B bonds. A generalised force constant
Wc(k) is obtained by this procedure, where the conjugate generalised displacement
is the mean extension of the A type intercalants. If the bonds are uncorrelated, the
effects of the other bonds of the same type being stretched cancels out, leading to

Eq. (3.41). If the bonds are correlated, this is no longer the case and W(k) must

 

63

be replaced by Wc(k). Note that when a single gallery is stretched with forces F at
either end, the effective spring constant is the same as for a perfect lattice order of
the static distortion and this extra distortion can be interchanged because the system

is linear. These are subtle points that require some thought.

3.6 PROBABILITY DISTRIBUTION FOR
GALLERY HEIGHTS

Fig. (3.13) shows the probability distributions PA(d) and P3(d) obtained from com-
puter simulations, for various representative cases. It is striking that these distribu-
tions are not bounded between 0 and 1, as some of the long bonds are stretched and
some of the short bonds are compressed. This did not occur in the bilayers studied in
I. The reason for this perhaps surprising effect can be seen most easily in the dilute
limit, where the distribution of the gallery heights of each intercalant is a main peak
surrounded by two satellites. If a single large intercalant is inserted into the system,
as in Fig. (3.2), then the adjacent short bonds in the same gallery are expanded,
whereas the two short bonds in the galleries above and below are compressed. More
interesting is that when kA = 1:5 , the two peaks for the A and B distributions in
the random case in Fig. (3.13) have the same shape. This unexpected result was first
suggested to us by the computer simulations and then proved rigorously. The result
is true for all compositions a: so long as both, the distribution remains random and
13A = 193. If either 19,; # kg or correlations are present this symmetry is destroyed.
This symmetry is especially interesting in the dilute limit, where it relates the host

and inclusion distributions.

 

64

 

 

 

1.0

 

P(da)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

_j I—fi 7 _I— 1 _ I 1*1 I 1 1 __ j I I Ifi I I q
'- R ’ " 1_A : — ,_ R c d
I— :: J “ \
' " ~ I ' " kA=1 -. _1
0.6 — II A h a 1_k8=5 ..
. 52 I . I: I
\ ~ . -t
0.0 1 l 1J1" 1k 0 1 I 1 Jnk 0 1 I 1 1 J ll "
0 l O 1 0 l
__ I j fi fit Ltj i ffi l J 6—fl I I I I I-J
h- 3. u )- i If 8‘: .J
'U i- :: -I i- d " :2 a
. . T [\l - .- \ILJ _
or JVLIA 1 0+ :I l l L O J. I l l l
O 1 O 1 0 1
d d (I

Figure 3.13: The probability distributions P (d A) and P ((13) are shown from computer
simulations for the 2d multilayer system, as a function of the dimensionless distance
(I. Here a: = 0.1 in the first two pairs of panels which have 1:: = kA = 193 = k1. Here
R, A stand for Random and Anticlustering respectively. The third pair of panels is

also Random with by = 13A and k3 = 515A.

65

3.6.1 Distribution Function PA(d) and PB(d) ‘

First let’s define the distribution function for CIA and d3,

 

a(d) a ((1 +26“) 6(d — do) , (3.57)

where (1 + 600/2 is the projection on A (or B) when 6 = +1 (or —1). Replacing the

5-function by a q integral, we get

. +00 _
Pc(d) = ((1 +2501) 51;;1-00 e:q(d—di)dq>

 

 

1 +00 .
_ _ 19d
— 211.00 e F((q)dq , (3.58)
where
1 i _- .
a(q) = (( +2“ ) e ml). (3.59)

and the dimensionless length di in the above equations can be obtained analytically

when kA = k3 by combining Eq. (3.3) , Eq. (3.7), Eq. (3.11), and Eq. (3.13),

di = :c + 2141in , (3.60)
.i

where

Aij = “$- [2(M—1)U — (M—l)i+uj - (M'lh—uj] . ' (3-51)

Substituting di into Eq. (3.59), we get ,

F¢(CI) = <(1 +2501) e—iqze-iq Zj Aijaj>

 

 

. 1 .
= e-qu << +2501) H(cos inj — iaj sin inj)> , (3.62)

.i
where we have made use of the identity 032 = 1. So,

1 —£z ' -
F.(q) = 56 " {<H(coquu-20581nqAn)>
J

 

66

‘

+€(H(COS inj — ioj sin qA-fi)) - ((0; cos ini — isin ini)) }3.63)
565

Notice that the product is over distinct galleries, and so expanding, we will have to
average quantities such as (010m...on), which is just (0'1)(om)...(an) in the random
case. This means that we can factorize all such products. Using (Uj) = 1 — 22:, we

have

F¢(q) = ge’iq” Ln(cos inj — i(1 —- 2x) sin un)
.fl

{ cos ini — i(1 — 2x) sin ini + e[(1 — 2:6) cos ini — isin QAii]}(3.64)

Making use of the fact that 62 = 1, we can rewrite F¢(q) as

 

 

’1 1 — 2 .
F¢(q) = + 6(2 3)) e""“’(cos in; — z'e sin ini) H[cos inj — i(1 — 22:) sin quj]
. . 535i
' _ 2 ‘ . .
= 1 + 6(12 2:) e""$"‘°‘4ii H[cos inj - i(1 — 2:13) sin in5] . (3.65)
. . j¢i

 

 

Comparing Eq. (3.61) with Eq. (3.24), we know that An = —1W(k). Therefore,

flafl) _ F+(<I)

F301) — F—(Q)

1 ... .
= ”a” . (3.66)

 

This tells us that the two distributions have a weight ratio of (1 — :c)/a: and their
positions are shifted by W(k). This is most easily seen by doing the back Fourier
Transform to obtain the appropriate probability distribution functions. Their shapes,
however, are exactly the same.

All information about the system is contained in the distribution function or its
Fourier transformation Fc(q). We have obtained the first and the second moments in

Sec. 3.3. We can also calculate the moments of the distribution from F€(q) as shown

 

67
below. Introduce an auxiliary distribution function F ‘(q) by the relation,

FA(q) = (1 — :r)e'iq”(l—W)F‘(q) , (3.67)

where

F‘(q) = e"(1’2“)q’lii H[cos inj — 7(1 — 22:) sin qAfi] . (3.68)
#i
FB(q) can also be expressed in terms of F'(q) as

F361) = xe-‘q“““-x)“'WHF‘(q) . (3.69)

Note that [1:(1 — W)] and [1 — (1 — x)(1 — W)] are the average values of (IA and
d3, respectively. The quantity F "(q) is the Fourier transformation of the reduced

distribution function P’(d) which has an average value of 0. From the equation

PM) = 2i / eiqu'(q)dq, (3-70)

71'

we obtain the moment expansion

M) = / P*<d)e-‘qdd(d)
= Z [% / P*(d>d”d(d)] q"

E Z c;q" . (3.71)

Note that (d') = 0, so that:

 

(4)77 . (3.72)

The moments Mn can be calculated from the expansion coefficients c; (see Appendix

B):

M. s ((d‘ — <d*>)"> = <(d*)"> = / P*(d>d“d<d> =

M1 = o (3.73)
M; = :c(1—:c)(W1—I/l/'2). (3.74)

 

68

These results agree with our previous results. Higher moments can also be obtained,

although they rapidly become very involved.

M3 = 83(1— $)(2x —1)ZA% , (3.75)
535
and
M4 -— 3M3 = a:(1 — :r)[1 —— 6.7:(1 — :c)] EA}? (3.76)
j¢i

We have checked the x dependence of M3 and M4 against our simulation results and
the agreement is excellent. We have not attempted to calculate the multidimensional
integrals necessary to obtain the constants involving the Ag,- in the expressions for
M3 and M4 above. This special symmetry is lost when either kA 54$ kg or correlations
are present as indicated in Fig. (3.13). The results we have obtained in this section
are general and apply to all types of lattice networks as long as 16,; = 163 and the two
species are distributed randomly.

There is an accidental symmetry for the 2d square net when 16,; = k3 = Is. It
can be seen from Fig. (3.13), when P = 1, the peak for A (or B) is symmetric
around its center. This symmetry can also be demonstrated rigorously by proving

the disappearance of all the odd moments of P’(d). (see Appendix C).

3.6.2 Low and High Concentrations

When a: is small, the B intercalants can be treated as defects. The single-defect and
double-defect problems can be solved exactly. For a single B defect at the origin
i = O, we can obtain all the A heights from Eq. (3.10), where the vector Q has only

a single non zero element:

 

 

(j) )2["11\7é:= [1 — co:( q UH cos(q - rj) , (3.77)
and
0733(0) = 71,-}: 2“ ” °::(q' ”)1 . (3.78)

 

69

Here, the superscript ‘3’ implies a single defect, which has a unique (1 associated with
it. Summing over j in Eq. (3.77) gives the distribution of dA in this limit. Eq. (3.78)
does not give the distribution of dB although it gives the position of the center of the
distribution. The double-defect configuration, say one B at i = 0 and the other at j,
gives us

2[1 - cos(q - 11)]

dB(i)='11V; ,q (1+cos(q.r,.)) (3.79)

 

Subtracting the position of the center dj;(0) from d3(j), we get the relative position

0f (13 ,i.e.,

d’BO) = da(i)- 33(0)

= dA(.l) . (3-80)

This explains why the distributions of A and B are identical in this limit. It is
interesting that different configurations contribute in similar ways to give this equiv-
alence. The above argument in the dilute limit is informative and in particular is
helpful in understanding the distributions when the ILC is turned on. In the random
case, the configuration {...ABBA...} has the same probability as the configuration
{...ABA...ABA...}. The angular brackets {....} denote a configuration along the
c direction. The {...ABBA...} configuration is highly suppressed in the antiferro-
magnetic (anticlustering) case. So we expect that the distribution of dB is affected
by this suppression. In {...ABBA...}, B is compressed to a smaller size than in
{...ABA...ABA...} [from Eq. (3.79)]. This leads to a decrease in the weight of the left
satellite of the dg distribution, while the distribution of d A is barely affected because
this effect is of order $2. This can be seen by comparing the first two panels in Fig.
(3.13). In the random case, P = 1, and we have satellites on both sides of the main
peak of the dg distribution. In the strong anticlustering situation, P = 500, and the
simulation results show that the left satellite of (13 vanishes whereas the distribution

of dA barely changes.

 

70

3.7 CONCLUSIONS *

We have investigated randomly intercalated multilayer systems described with a
harmonic spring model and have obtained an analytic solution for certain cases. We
also have reached an understanding of the probability distribution function for the
heights of the galleries containing the intercalants and the modification caused by
interlayer correlations. Our effective medium theory for both random and correlated
intercalation works extremely well in the representative situations tested. It is hoped
that our results will be tested experimentally using EXAFS in the near future.
The formalism used in this work is quite generally applicable and has been used to
describe random alloys where the displacements are not uniaxial.[19] Comparing our
theoretical results with experiment, the system that comes closest to the situation
RA = R3 is the stage-1 ternary K1_,,beCg. According to our model, this system
should show a Vegard’s law behavior for the average gallery height(average c-axis
periodity) ((1). There is however a small deviation from Vegard’s law; the origin of
which is electronic and cannot be accounted for within the present elastic model. In
intercalated layered silicates such as 121714563,r vermiculite, one observes a sigmoidel
ctr—dependence of (d). The initial sublinear x-dependence has been understood[10] in
terms of the preferential occupation of non-gallery sites by the larger Cs intercalants.
These intercalants do not contribute to the gallery expansion. If one corrects for these
sites and considers (d) as a function of X9, the concentration of Cs ion occupying
the gallery sites, then one observes a superlinear behavior. This can be understood
within the present model as the Cs ion is larger than the Rb ion. Gallery expansion
in LizCs and LixTz’Sz can also be included within the present model by treating
these systems as Vl_$LirC'6 and I/1-xLi$TiSz respectively, where V is a vacancy and
RLi/RV > 1.[20] More detailed comparison with high quality experimental results

remains to be done.

Bibliography

[1] F. S. Galasso, Structure and Properties of Inorganic Solids (Pergammon, New

York, 1970).

 

[2] L. Vegard, Z. Phys. 5, 17 (1921).

[3] J. B. Thompson, in Researches in Geochemistry, Vol. II, ed. by P. H. Abelson,
(Wiley, New York, 1967) p. 340.

[4] J. R. Dahn, D. C. Dahn, and R. R. Haering, Solid State Comm. 42, 179 (1982).
[5] J. E. Fisher and H. J. Kim, Phys. Rev. B 35, 3295 (1987).

[6] S. A. Solin and H. Zabel, Adv. Phys. 37, 87 (1988).

[7] P. C. Chow and H. Zabel, Phys. Rev. B 38, 12837 (1988).

[8] D. Medjahed, R. Merlin, and R. Clark, Phys. Rev. B 36, 9345 (1987).

[9] w. Jin and s. D. Mahanti, Phys. Rev. B 37, 8647 (1988).

[10] H. Kim, W. Jin, S. Lee, P. Zhou, T. J. Pinnavaia, S. D. Mahanti, and S. A.
Solin, Phys. Rev. Lett. 60, 2168 (1988); S. Lee, S. A. Solin, W. Jin and S. D.
Mahanti, in Graphite Intercalation Compounds: Science and Technology, Ed. by
M. Endoh, M. S. Dressellhaus and G. Dressellhaus, p. 41 (Materials Research
Society, 1988).

[11] M. F. Thorpe, Phys. Rev. B 39, 10370 (1989).

71

72

[12] M. F. Thorpe, W. Jin, and S. D. Mahanti, Phys. Rev. B 40, 10294 (1989).
[13] M. Mori, S. C. Moss, Y. M. Yan and H. Zabel, Phys. Rev. B 25, 1287 (1982).

[14] K. Komatsu, J. Phys. Soc. of Japan, 6, 438 (1951); see also S. Lee, H. Miyazaki,
S. D. Mahanti, and S. A. Solin, Phys. Rev. Lett. 62, 3066 (1989).

[15] D. C. Mattis, The Theory of Magnetism, (Harper & Row, New York, 1965),
p. 258.

[16] See for example, J. C. Mikkelson, Jr. and J. B. Boyce, Phys. Rev. B 28, 7130
(1983).

[17] S. Feng, M. F. Thorpe, and E. Garboczi, Phys. Rev. B 31, 276 (1985).
[18] R. P. Feynman, Phys. Rev. 56, 340 (1939).
[19] M. F. Thorpe and E. J. Garboczi, Phys. Rev. B42,8405(1990)

[20] S. Lee, H. Miyazaki, S. D. Mahanti, and S. A. Solin.(See Ref.14)

 

Chapter 4

STRUCTURAL
CHARACTERISATION OF
SEMICONDUCTOR ALLOYS

This chapter is to be submitted to Physical Review B. This preliminary version is

written by the author and the paper will be published under the names of Y. Cai and
M. F. Thorpe.

4.1 INTRODUCTION

Electronic properties of semiconductors can be changed by varying the concentration
of components. Structural information on semiconducting materials is of fundamental
importance in calculating and predicting these properties. Extended X-ray absorp-
tion fine structure (EXAFS) experiments have found that pseudobinary semicon-
ductor alloy system (A1_xBx)C shows a bimodal structure.[2]-[7] The first neighbor
cation-anion distance remains closer to that in the pure binary compound than to
that of the average or virtual crystal distance. The discovery inspired considerable

theoretical interest.[8]-[11] All these theories are quite successful though they focus

73

 

74

on the dilute limit. In this chapter we study all the pseudobinary semiconductor
alloys, with arbitrary concentrations within the valence force model. They are shown
to be isomorphous to the general length mismatch problem. [13],[14] The model
is also applicable to binary alloy system like GezSi1_3 and the quaternary alloys
A1_,B,,Cl_,,Dy. Our investigation provides a better understanding of the Vegard’s
law[1]. This chapter will be presented in the following way. As an introduction to
later work, we first discuss the valence force models in sect. 4.2. The force constants
of the models are fit from elastic measurements for known pure binary crystals. In
sect. 4.3 we apply the appropriate valence force model to the alloy system and give
exact solution to pure length mismatch problem when there is no force constant disor-
der. The simple result allows us appreciate the importance of the Watson Integral,
which characterizes the topological rigidity of the lattice system. We will also discuss
its importance in the construction of an effective medium theory, for use when there
is force constant disorder. In sect. 4.4 we study in considerable detail all ternary
I I I - V and II — VI semiconductor alloys using both effective medium theory and
computer simulation. It is found the theory very well reproduce the structural fea-

tures observed in experiment.

4.2 VALENCE FORCE MODELS

In this section we study the Kirkwood Modeland Keating Model[15] together in a
combined form. Valence force models have been used to study semiconductors quite
successfully. In diamond or zincblende structure, nearest neighbor interaction along
can not solidify the system. There are one third vibrational modes of zero frequency
because of the understraining. Farther range interactions are needed to stabilize the
system. Kirkwood and Keating model contain second neighbor interactions and are
all able to stabilize the zincblende structures that we are going to discuss. There are

some subtle differences between these two models. And Keating model is tested to be

 

75

better. This may be due to the fact that Keating derived his model from symmetric
point of view, while Kirkwood model does not satisfy the rotational invariance. We
will see the difference is minor. In Kirkwood model, the angular restoring part of the

deformation energy is written as the square of angular difference,

a I
Ekgrkwood = -2— 2(hgj — ha)” + -8— 2(cos 9,31 - cos 0%,)2 (4.1)
(ii) (iii)
a and ,6’ are nearest neighbor and angular force constant respectively. hij is the
bond length between atoms 2' and j. hf), is the natural length. 9.3; is the angle made

9.
1.11

by nearest neighbor bonds i j and ii. The natural angle 9 is usually taken to be
about 109° as in perfect diamond lattice. The sign () under summation excludes the
double counting.

The energy can be written in expansion form for small length mismatch. Denote
he the nearest neighbor distance of the underline crystal structure, and u.- the dis-

placement vector of atom i from its crystalline position, and keep to linear term,
hgj 2' he + fij - u,-, (4.2)

where in,- is the unit vector in crystalline structure pointing from atom j to i.

Redefine the angular force interaction parameter,

_Ei'

_. h, (4.3)

we have the expansion form(only keep to linear terms for simplicity and it is found

to be a good approximation):

(1 , . A 1 A ,.
Ekirkwoad = 52(hc—h?j+reyus)2+§ Z)[rij°uil+ril°uij+§(rij°uii+ri1'uil)l2 (4-4)
(‘1') (53"
0

In the perfect system, there is no length mismatch and he = hij, we have,

a . . . 1 . .
Ekirkwood = 5 23(18): - 11.5)2 + g- le‘ii - 11.7 + ru 418' + §(r.-,- ' “ij + Pu ° 11.7)]2 (4-5)
(ii) (if!)

 

76

To compare with Keating Model, -

a A A A
Ekeating = a 2(1'1'3' ' 11.3)2 + g- 2(1'53' ' “a + m ° uij)2 (4.6)
(‘1') ('31)

It can be seen that these two models are very similar except there are two more
terms in Kirkwood model. In fact we can combine them in a single model and study
them together,

01 . 2 fl . . A . ~ 2
EA = 5 2035 ' uij) + g le‘ij ' 11.7 + ril ° uij + 3033 ° uz'j + rel ' 11.7)] (4.7)
(if) ('3’)
where /\ = 0 stands for Keating model and z\ = 1 stands for Kirkwood model.
The elastic modulus of the combined A-model can be solved in the similar way as

in Keating’s original chapter.[15] The results are,

 

__ 1 216 612
Cu — 1(a+3,6’—3_‘+T)
_ 1 216 612
012 - 4(a-fi-—3 +—3)
_ 1 216 m2 [awe-W912
0‘“ ‘ 4{“+fl" 3 + 9 ' a+fl(1+A/3)2} (4'8)

and we can also work out the frequency of the optical phonon at k = 0 for the

combined A-model,

 

.50 = (Em + 5(1+§)2] (4.9)

where m is the mass and the lattice constant is taken to be one.

The force constants in the models are so chosen to fit Cu and bulk modulus. We
list in Table.(I) the force constants for different compounds. It can be seen that [3 is
the same for both the Keating model and the Kirkwood model. The a constants do

not change too much from Kirkwood to Keating description. The main reason that

 

77

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Kirkwood K eating
0 fl 0 fl

AIP 65.03386 14.19144 60.30338 14.19144
Al As 44.18573 8.94280 41.20480 8.94280
AIS b 35.69020 6.79200 33.42620 6.79200
GaP 48.06614 10.69252 44.50197 10.69252
GaAs 44.34032 9.25720 41.25459 9.25720
GaSb 34.29977 7.33188 31.85581 7.33188
I nP 41.72646 6.60229 39.52569 6.60229
I nAs 35.09977 5.75993 33.17979 5.75993
1713!) 31.30400 5.07011 29.61396 5.07011
Z nS 40.30599 4.78803 38.70998 4.78803
ZnSe 33.74111 4.56242 32.22031 4.56242
ZnTe 31.06783 4.66933 29.51139 4.66933
C dS 35.66434 4.75621 34.07893 4.75621
CdSe 33.18513 4.37257 31.72761 4.37257
CdTe 27.30606 2.72412 26.39802 2.72412
H gTe 30.72841 2.93264 29.75086 2.93264

 

 

Table 4.1:

 

 

 

78

we adopt the Kirkwood model is its simplicity when there is size difference in the
system. As can be seen in Eq.(4.1), the length mismatch enters only in the a part.
This makes it easy to be solved in a familiar way.[13],[14] Considering that the 3
part contributes less than the a part, we expect that the two models should generate
qualitalively similar results.

For most of the binary alloys we studied, as in table (I), the fitted force constants a
and 6 give C44 and we within 20% of the experimental values. Of course it is possible
to choose different angular force constants for ABA and BAB type bonding and this
leads to better fitting of the elastic constants and optical mode frequency. One can
further consider the charge transferring and deploy shell model. However our current
effective medium theory does not tolerate too many disorders in the force constants
and on the other hand the local structure is not very sensitive to the difference. And
more sophisticated potential like the Embedded Atom Potential can only be studied
in numerical simulations. We will adopt the one—set—parameter fitting to get simple

theoretical predictions.

4.3 THEORY

In this section we first solve the pure length mismatch problem of (A1_,B,)C sys-
tem. Then we study the Watson Integral. Finally we discuss the effective medium
method.

The structure is zincblende for pure AC and BC binary system. For example,
in the AC binary system, A and C each occupy one of the FCC sublattices in the
zincblende structure, and the same for BC binary alloy. To make (A1_1BI)C, AC
and BC alloys are mixed with appropriate concentration 1 — x and 2:. In the final
structure, C-type atoms take one of the FCC sublattices while A and B occupy the

other FCC sublattice randomly. We assume no correlations between A and B atoms.

79

4.3.1 Solution to Pure Length Mismatch ~

Our theoretical approach to the problem assumes harmonic approximation. This is
valid in small mismatch situation. Start from the energy form of Eq.(4.4). It can be

written in compact form. We discard the subscript “Kirkwood” from now on,
1
E = §U+MU — U+d> + E0; (4.10)

where U is the displacement vector and M is the structure or connectivity matrix.

The disorder vector qfl can be expressed as,

 

d’ = Zhfiulij >
(£7)

= Eqs.-,6} > (4.11)
(*1)
where
¢ij = ‘01“: — hfj) (4'12)
and E0 takes the form,
0
E0 = .5 20., — 59,.)2 (4.13)
(ii)
Define the reduced length d as,
h - (‘91

Apply the same idea as in multilayer system[13], we get the total average length,
average length of AC type bond, BC type bond and distortion energy(see Appendix
D).

The results turn out to be very similar in form to the results we reached in the

Intercalated Graphite Multilayer system[13] and the 2d triangular network.[14]

(d) :1:

((1.40) = $(1”a")

80

(1136') = 1— (1 — x)(1— a")

1 ..
e = Ex(1—:c)(1—a ) (4.15)

e is the reduced energy per site. The Watson integral a" here is for sites case
instead of bond case. It can be evaluated from the static Green’s function directly(see
Appendix D).

The simple results above should be properly appreciated because of the correlation
involved in the site mixing case. In the AB type random alloy the solution is lengthy
and more lattice integrals are involved.[17] The compact form is attributed to the
regular occupying of one sublattice by C type atoms.

To write down the length fluctuation it needs the introduction of another “Watson”
integral a;‘,[13]

26(a“*/a)

a1“ = —a T (4.16)

where a the central force constant. For example the fluctuation of AC type bond
length,
was) — (4...)? = 241— our — a“) ' (4.17)
The derivation for fluctuations is more difficult than that given in Appendix
D but shares the same old idea.[13],[14] In triangular and FCC lattice[14] a" is
independent of the force constant Q. So a‘f" = a". It is not the case in Zincblende
structure because of the angular force term in the energy expression. In general these
quantities depend on ratio of the force constants as well as the concentration. We

will address this aspects of the Watson Integral specifically next.

81

 

K eating K irkwood
C1 C2 03 C1 C2 03
a‘ 2.198 3.200 0.796 2.0058 3.200 0.880
a" 2.2786 4.600 2.629 1.2489 3.600 1.171

 

 

 

 

 

 

 

 

 

 

 

 

Table 4.2:
4.3.2 Watson Integrals as Function of Force Ratio

We need to study the Watson integrals because they are the only parameters involved
in all local properties. These integrals characterize the rigidity of the system and
contain the information about the correlation. The ratio of the force constants,
y E (ll/,6, varies from 0.1 to 0.6 under valence force model description. For the
purpose of effective medium study, we needs a simple analytic form to describe the
change of these integrals as function of the force ratio. And it is found the following

form well describes the Watson integral in the ratio range from 0.0 to 0.7,

 

1 + 013/
' = ——-——————- 4.18
a (y) 1+02y+03y2 ( )
The integral used in fluctuations associated with a“ is defined as:
. 0(99‘)
“1(9) = By (4-19)

The values of the parameters are listed in Table. (II). And Fig.4.1 demonstrates

the good representation of the Watson Integral by the above simple analytic form.

4.3.3 Effective Medium Theory

When there are disorder not only in size, but also in force constants, there is no
analytic solution. And in practice this is always the case. We need to develop methods
with more approximation. The Effective Medium Approximation we used here has

been applied by Thorpe and et. al. to many other problems. The essential idea is

82

Ke ating Kirkwood

 

 

 

 

 

 

0.0 - . ' - - 4 0.0
o 0.2 0.4 0.6 o 0.2 0.4 0.6

49/9 fi/a

 

Figure 4.1: The Watson Integral as function of the force ratio 01/6. The top panels
for a' and bottom for a". Diamond denotes a and square denotes al which relates to
the derivative of a and can be calculated from the fluctuation of distribution. These
two symbols are from computer simulation. Straight lines are the three parameter

fitting.

83

to replace the effect of the matrix on an embedded chemical bond bnnother bond
of appropriate strength and natural length so to reach self consistency. The rigidity
and correlation information are contained in the Watson Integrals which appear in
the expressions of the quantities concerned.

We will list the results only and leave the derivation to reference paper[13] for
interested readers. The results are written in reduced variables (1 instead of h, and

subscript A means AC bond and B means BC type.

 

 

 

 

 

(d) = a: + a:(1 — 3:)Fd (4.20)
$16,,ka
(d4) = ’9ka + kA)(kf+ kg) (4.21)
_ (1— $)kek;kA
(d3) ‘ 1 " We: + was; + Ice) “'22)
5 = -;'\/ICAICB 33(1— x)(h% -— (13021;: , (4.23)
where
_ 1661:2033 — 16);)
F" ' k(k;+ mac; + k3) (4‘24)
and
F. -.= iii—Vi“? . (4.25)
And the fluctuations,
2 _ 3(1- xlflkAkBkEl/[kfké + k4)2(k; + kBll}2
“d" ‘ ”A” l “ W2/(W1— W) + [(k- 616— kB)]/[(k;+ k4)(ké+ks)l (4'26)
«(1 _ ((0)2) _ $(1 _ x){[kAkBkel/[k(kiz + kA)(k:: + kB)]} Wl/(Wl — W ) (427)

_ W2/(W1 " W2) + [(k — k4)(k - [W)]/[(192 + 19.4)ka + ’63)]

In pure length mismatch problem, the EMT goes back to exact solution, and we
observe the generalized Vegard’s law, i.e., linear changing of average length with
concentration, not only for the total average length, but also for specific average length
AC and BC, even when there is correlation. However, Vegard’s law is destroyed when

the force constants a and B are different for AC and BC binary systems.

 

84

In real compounds, the difference in the ratio of force constants, a/fl, is not large.
And this leads to the non-remarkable deviation from Vegard’s law in the experimental

data. But in general, Vegard’s law is not obeyed.

4.4 APPLICATION

Out of the 36 possible ternary alloys that could be formed from Zn, Cd, H g, P,
As, Sb, and Al, Ca, In, Si, Ce, Sn, we study the 29, for which force constants
can be extracted from experimental data. We use the Kirkwood model so to get the
comparison between simulation and effective medium theory.

It is easy to solve the effective medium approximation using the function form of
the Watson Integral given in Eq.(4.18). Fig.(4.2) to Fig.(4.5) show the computer sim-
ulations and effective medium results. The excellence of EMT is well demonstrated.
Compared with experiments data availible, we plot the “Z” curves for I nxGa1_,,-As
in Fig.(4.6) and anCd1_xTe in Fig.(4.7).

We see that both simulation and EMT give fairly good account of the experimental
measurements. All the curves are essentially straight plus a miner bowing, as the

following description,
h = p1(1 — :22) + p21: + p3a:(1 — 2:) (4.28)

p1 is the length at the zero concentration and p2 on the 100% side. p3 denotes the
bowing of the curve.
The three parameters can be solved from EMT numerically. We can also derive an

approximate values of these three parameters from EMT equation.
Denote £1 = (01.4 + aB)/2, and ,6 = (6,; + BB)/2. For AC type bond,
his - hit
1+ CM/{Otza[1/a"""(/5'19/0113) -1l}

 

P1 =h0A+

 

85

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Asl—xPx Sb,_,P, Sbl—xAsx
2.450 ‘ . .
, 2.0 ~ .
2.425 5 . 2.0 -
2.400 » . 2.5 1 1 Al
l
2.5 J
2.375 1 . 2.4 r
8.450 47 4 1
‘ 2.0
2.420 » 1 2.0 .
7
2.400 r . 2.6 f . Ga
» . 2.5 .
2.375 2.4 _ 1
2.825 ’ I 2.8 i 2.80
2.000 » .7 ._
2.7 - 2 5 ,
2.575 ~ 1 l 2.70
1
2.550 _ 1 2.6 4 12.85
x
Figure 4.2: .

The “Z” curve for III — V. Solid lines are from Effective Medium Theory. The

Kirkwood model is used in the computer simulation

86

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

A11_,Gax In1_,,Alx In1_,,Gax
2.305 » . 2.5 r
P
. . 2.4 4
2.000 » 24
2.450 -
2.440 ~
2.05 ~ .
Sb
2'7 _ 2.7 -
2.04 -
0 0 5 1 0 0.5 1 o 0 5 1
K
Figure 4.3: .

The “Z” curve for III — V. Solid lines are from Effective Medium Theory. The

Kirkwood model is used in the computer simulation

87

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Cdl-xznx Hgl-xznx Cdl-xflgx

2.50 ' '

2.0 -

2.0 - 1

2.45 -

2.5 ~ - S
2.40 -

2.5 ~

2.35 - . 2.4 ' .
2.00 ~ - .

2.0 2.0 >
2.50 -

2.6 ' ' Se
2.50 - ~ 2.5 ~

2.4 - -
2:06 - i 2.010 ~
2.75 - - 2-75 ’ '2.005 -
2.70 - . 2.70 - '2000
2.05 - . 2.05 - _ 2.795

0 05 1 0 05 1
x x
Figure 4.4: .

The “Z” curve for II — VI. Solid lines are from Effective Medium Theory. The

Kirkwood model is used in the computer simulation

2.45

2.40

2.35 7

2.60 7

2.55 7

2.825 7

2.600 7

2.575 7

2.550 7

88

 

    

 

 

 

 

 

 

 

 

 

 

2.0 -

2.00 -
2.5 7 - 255 .
2.4 . 2.50
2-8 7 2:33
2.7 _ . 2.75 "

2.70 -
2.0 - -

2.05
2-8 7 7 2.00 ‘-

2.75 -
2.7 - 4

2.70 -
2'” ’ ' 2.05

05 1 0 05

Figure 4.5: .

 

 

 

 

Te1_xSex
' Zn
- Cd
.Hg
5.6 1
x

The “Z” curve for II — VI. Solid lines are from Effective Medium Theory. The

Kirkwood model is used in the computer simulation

89

 

 

 

 

 

I I T T 1 I I T I I ‘I I j I I I I I I I I I I I d
d
2.50 —- g i -—
2 ' -l
V 2.55 —- .J
8 7 2 Simulation 7
g C I Experiment '7
*3 4
.Z’. 2.50 -—-—- EM ....
a .- -
z i "
z . "
2.45 —

l- 4

J J l J l L l l A l l l l l l l I J J l l J l 1
o 0.2 0.4 0.8 0.5 1

Figure 4.6: The “Z” curve for I n1_xGaxAs. Solid lines are from Effective Medium
Theory. The Kirkwood model is used in the computer simulation. Experimental data

goes with error bar.

90

 

2.80

    

 
  
 
 

 

 

 

A
4. 2.75
V
a) ' — EM
0 . - fl
5 7 ° Slmulation ‘
_,_, 2.70 — —
E i Experiment
2
z
z 2.65
l l l l I l I l I I I l J I l l J l I I l I l J .1
o 0 2 0.4 a a 0 a 1
x

Figure 4.7: The “Z” curve for Cd1_IZnITe. Solid lines are from Effective Medium
Theory. The Kirkwood model is used in the computer simulation,experimental data

goes with error bar.

91

uux"— a -0 null-— flB—fl
p3 = {[1 — a (fl/a)]2(—B—a—-A7) + a (fl/a)(—7—‘)1<h°8 — ha) (429)
The prime means derivative over the argument inside. For the BC type bond, the
parameter p3 takes the same form. 101 takes the form of p2 for the AC bond, but with
A and B switched. The same change occurs for parameter p2.

For the total average,

P1 = hia'
P2 = hi4
p3 = {[1—a7*(B/a)1(fl9—§—“fi)(h%-h9.> (4.30)

It is interesting to analyze the approximation expression above. We can see that
the bowing of the curve for AC and BC are the same. They depend on the differences
of a and fl of the two binary components. However, the bowing of the total length
average only depends on the difference of the central force a. For most alloys the
bowing of total average and specific average go in opposite directions, i.e., the angular
force difference is dominating in the bowing in specific length average.

The length distribution is approximately Gaussian. Fig.4.8 shows the distribution
of nearest neighbor in Cdo.5Zno,5Te system. There is barely any overlap between
two peaks, in contrary to the situation in metal alloys, where distributions have very
remarkable overlap. This is mainly due to the soft angular force in semiconducting
materials, in which deformation in atomic distance is more energetic than that in the
angle of two neighboring chemical bonds. The widths of AC and BC distribution
is not exactly the same because of the difference in force constants and also because
of the bond-correlation in site disordered system. Like in correlated Intercalated
Graphite system, EMT does not give good prediction for the width of the distribution.

Of course we can get these informations from computer simulations.

 

92

Cd0.5Zno.5Te

 

 

 

 

 

 

 

 

,_ 1 I I I j 1 j j q

0.01 — _.

g .- -
"-0
a.)

g " '1
02-!
I...

Tn, "' -1
"III
2

0.00 7-7 ....

. l i . . . l L
0 1

Figure 4.8: The length distribution of Cdo,5Zno.5Te. The shapes are approximately

gaussian. Correlations and the difference in the a force constants make the two peak

different.

93

4.5 CONCLUSIONS _

We studied the properties of the Kirkwood and the Keating model by combining
them into a unified model. Then we studied the structure of the (A1_,B,)C system.
An exact solution is found for the pure length mismatch problem. And effective
medium theory is applied to multinary compounds and get good agreement. While
preliminary information is available from theory, the full configuration of multinary
system can be solved by computer simulation. The X -ray scattering calculation can
also be performed.[16] It reveals the long range properties of structure of multinary
compounds. The accomplishment allows us to appreciate more of the simplicity by
introducing the concepts of chemical bond and natural length in solid state materials.
And we point out that the bimodel properties is characteristic of all “length mismatch
problems”. In these systems the nature of the chemical bond is assumed not to change
significantly in different circumstances and a natural length for the chemical bonding
is adopted to simplify the interaction picture. The relaxation of the natural length in
mismatched system inevitably leads to the bimodel observations that the deformed

chemical bond will stay as close as possible to its natural length.

 

Bibliography

[1] L. Vegard, Z. Physik 5 (1921) 17.

 

[2] J. C. Mikkelsen and J. B. Boyce, Phys. Rev. Lett. 40, 9001978.
[3] J. C. Mikkelsen, Jr. and J. B. Boyce, Phys. Rev. B28(1983).

[4] J. C. Mikkelsen and J. B. Boyce, Phys. Rev. Lett. 49, 14121982.
[5] J. B. Boyce and J. C. Mikkelsen, Jr., Phys. Rev. B31 (1883) 6803.

[6] A. Balzarotti, in: Ternary and Multinary compounds, Eds. S. K. Deb and
Zunger(Material Research Society, Pittsburgh, PA, 1987)p.333.

[7] T. Sasaki, T. Onda, R. lta and N. Ogasawara, Japan, J. Appl. Phys. 25 (1986)
231.

[8] J. L. Martins and A. Zunger, Phys. Rev. B30 (1984) 6217
[9] A. B. Chen, A. Sher, Phys. Rev. B32 (1885) 3695

\10] C. K. Shih, W. E. Spicer, W. A. Harrison and A. Sher, Phys. Rev. B31 (1985)
1139

[11] M. Podgorny, M. T. Czyzyk, A. Balzarotti, P. Letardi, N. Motta, A. Kisiel and
M. Zimnal-Starnawska, Solid State Commun. 55 (1985) 413

[12] Y. Cai and M. F. Thorpe, Phys. Rev. B40 (1989) 10535

94

95

[13] Y. Cai, J. S. Chung, M. F. Thorpe and S. D. Mahanti, Phys. Rev. B42, 8827
(1990)

[14] M. F. Thorpe and E. J. Garboczi, Phys. Rev. B42, 8405 (1990)

[15] P. N. Keating, Phys. Rev. 145, 637(1966)

[16] M. F. Thorpe, J. S. Chung and Y. Cai, to be published in Phys. Rev. B.
[17] J. Chen and M. F. Thorpe, unpublished

[18] “Group 111: Crystal and Solid State Physics”, V22, semiconductors, Landolt-

Bornstein.

 

 

 

 

Chapter 5

DIFFRACTION FROM
RANDOM ALLOYS

This paper is published in Physical Review B. Reference: M. F. Thorpe, J. S. Chung
and Y. Cai, Phys. Rev. B43, 8282(1991)

5.1 INTRODUCTION

The elastic scattering (diffraction) from a random alloy is a complex problem that
has received considerable attention over the years. In this chapter, we describe the
diffraction pattern caused by the correlated static displacements from an underlying
crystalline structure. We find a set of Bragg peaks (determined by the average lat-
tice), modified by an appropriate Debye-Waller factor. Each Bragg peak has some
diffuse (Huang) scattering associated with it. To quote from Huang,[1] “Considering
a crystal lattice formed of randomly distributed atoms of two kinds mixed in compa-
rable proportions, it must obviously be extremely difficult to describe qualitatively
the distorted configuration”, and he goes on to consider the case of isolated impuri-
ties. Modern analytic techniques and computer simulations now allow us to treat the

concentrated alloy. Although the example worked through here is simple, we have.

96

 

97

found it instructive and we hope the reader will too. _

We use the model alloy of Thorpe and Garboczi[2] (henceforth referred to as I) for
an A1_,Bx alloy. This consists of a crystalline lattice containing two kinds of bonds,
A and B which are randomly positioned throughout the lattice with probabilities
1 — :c and :5 respectively. The static structure is obtained by minimizing the energy
associated with the potential

v: €205.- —R‘.-I—L9,-)2 (5.1)

:1
Here L2,. can take on the two values L94 and LOB with probability 1 — 2: and a: re-
spectively. The summation goes over all nearest neighbor bonds 2' j and the angular
brackets denote that nearest neighbor bonds are only counted once. The vector R,-
goes to the site i at the end of the bond 2' j from some (arbitrary) origin. The potential

(5.1) can be minimized with respect to the R. to give,
0 = 21(12- — ii.) — Lia-Ra] (5.2)
J

Here Rh- is a unit vector along the relaxed bond direction. These equations determine
the equilibrium positions of all the sites, as described fully in I. An example is shown
in Fig.5.1, which can be visualized as static concentration waves of all wavelengths
freezing out from a perfect crystalline lattice. We will take all the spring constants
to be equal in this note as it simplifies the analytic treatment as discussed in I, and

does not qualitatively effect the result.

5.2 COMPUTER SIMULATIONS

We have used an FCC lattice rather than the triangular net as in I. This is because of
the well known instability of two-dimensional lattice ordering to fluctuations.[3] While
this does not effect quantities like the mean bond length, it does have a profound effect

on the diffraction pattern, which is very sensitive to long range order. The general

 

98

 

 

 

 

 

 

 

 

 

 

 

---- on--- D...
l '\‘
0
I ’ ‘
' I \

----’ ----

\

Figure 5.1: A piece of a relaxed triangular network. This figure is reproduced from I.
The short bonds are shown as dashes and the long bonds by solid lines. The sample
shown has equal numbers of short and long bonds (:1: = 0.5) and the natural length
of the long bonds is 30% greater than the natural length of the short bonds. The two

spring constants are equal.

 

99

formalism of I, which is geared towards the calculation of local distances, applies
equally well in any dimension. We have used M x M x M FCC samples with periodic
boundary conditions containing N = 4]” = 32,000 atoms and 241” = 192, 000 bonds,
which were randomly assigned to be A or B with probability 1 — a: or :1: respectively.
The simulation program, which uses a variant of the conjugate gradient method,[4]
adjusted the positions of all the atoms and the size of the supercell to minimize
the energy. The simulation was terminated when the energy agreed with the exact
energy 3N2:(1 — x)(L% — Lay/2 as given in I, to better than 1%. The supercell was
kept strictly cubic to facilitate the analysis. The mean bond length (L), which is

proportional to the sample size, is given by Vegard’s law as discussed in I
(L) = (1 — ”Lg + :rLoB (5.3)

The diffraction pattern was computed from the relaxed coordinates using

101) = 73—1 :expaq- a)? (5.4)

where 1(q) is the scattered intensity per site.

5.2.1 Bragg Scattering

The diffraction pattern is dominated by the Bragg peaks shown in Fig. 5.2, for
different directions in reciprocal space, and for different values of the length mismatch
parameter (Lg; — L9,) / L. The Bragg peaks are Kronecker delta functions, rather than
Dirac delta functions, because of the finite size of the supercell. We have taken
out a factor N from the computed intensity, to make the connection between the
computed Kronecker delta functions and the calculated Dirac delta functions used in
the next section. These Bragg peaks have no width and are associated with a single
superlattice point in reciprocal space. It can be seen that the diffracted intensity falls
off as a Gaussian in all three cubic directions. This is certainly what would be naively

expected, and is similar to the situation with dynamic disorder caused by phonons

100

 

1.00 .:

0.75

0.50

13(9)

0.25

 

 

0.00 '

 

Figure 5.2: The Bragg scattering for a random FCC alloy with two kinds of bond
lengths, similar to that shown in Fig. 5.1. The length mismatch parameters are
4% (upper curve), 8% (middle curve) and 16% (lower curve) and the concentration
.7: = 0.5. The symbols indicate the different cubic directions and the solid line is
the Debye-Waller factor calculated in the text. The intensity is averaged over the
different equivalent directions of a single 20 X 20 x 20 cubic sample containing 32,000

atoms. The solid line is the theory (5.8) using the Debye—Waller factor (5.12).

101

that produces a Debye-Waller factor,[5] although the static case requir_es handling in
a slightly different way.

5.2.2 Huang Scattering

The computed background is artificially reduced by a factor N with respect to the
Bragg peaks because of the Kronecker delta functions caused by the supercell, and
therefore cannot be seen in Fig. 5.2. The background scattering is shown in Figs.
5.3 and 5.4 for different values of the length mismatch parameter. The noise in the
results was greater than for the Bragg peaks and so we have averaged over 9 samples.
The scattered intensity is shown at the 2M — 1 = 39 superlattice points between pairs
of Bragg peaks in the 100 direction. This Huang scattering[1] is seen to be largest at
the Bragg peaks in the 100 direction. Similar results are found around all the Bragg
peaks.

5.3 THEORY

It is convenient to separate the scattering into two parts; the Bragg scattering IB(q),

and the remainder that we will refer to as the Huang scattering IH(q), so that

I(Ci) = 13(q) + 111(01) (5-5)

5.3.1 Bragg Scattering

The cross section (5.4) can be conveniently rewritten as
1 . .
1(q) = N Z eXPlzq- R31<6Xplzq- 11.31) W»
',J

where R,- = R? + u,- , and the R? define the underlying mean lattice and the u;
represent the (small) displacements from it. The angular brackets (...) in (5.6) denote

an ensemble average over different static configurations. The size of the mean lattice

 

1.0:)

102

 

 

 

 

 

 

 

 

 

 

 

0 10 12

Qi/Zvr

Figure 5.3: The Huang scattering in the 100 direction for a random FCC alloy

with two kinds of bond lengths. The length mismatch parameter is 4% and the

concentration x = 0.5. This figure is a magnified version of the Bragg peaks shown in

Fig. 5.2 and is averaged over the different 100 directions for nine 20 x 20 x 20 cubic

samples, each containing 32,000 atoms. The solid line is the theory (5.16) given in

the text and the vertical bars mark the positions of the Bragg peaks shown in Fig.

5.2.

 

103

25

 

A‘7

207 :>

15- . .

In(q)

 

 

 

 

 

 

v
v

 

 

 

 

 

Figure 5.4: The same as Fig. 5.3, except that the lattice mismatch parameter has

been increased to 16%.

 

104

is given by Vegard’s law (5.3). Inasmuch as the diffraction pattern.is dominated
by the Bragg peaks as shown in Fig. 5.2, it is suggestive that we should make the

conventional approximation,[5]

(eXP(iq- 115)) = (eXP(iq - ui))<€XP(-iq ° um = exp(-2W) (5-7)

which assumes that the displacements at each site are independent, as with ”frozen
Einstein oscillators”, where the displacement at each site is chosen independently.
This is manifestly not the case here, as we will discuss later. Nevertheless the as-

sumption (5.7) leads to excellent answers. Inserting (5.7) into (5.6), we find that,

13(01) = eXP(-2W)5(q - g) (53)

where g = (\/2/L)(n1, n2, n3) with the integers n1 ,n2 ,n3 either all even or all odd,
are the reciprocal lattice vectors for the FCC lattice with a mean lattice spacing L
given by Vegard’s law (5.3).

In I, it was shown that the displacements u,- , could be expanded in terms of the

Green functions ng of the perfect lattice, the various natural bond lengths Lg- , and
the nearest neighbor vector directions R4,- of the perfect lattice
u.- = K E G“ -1‘1,.,,(L?,,, — (L)). (5.9)

1,711
After some manipulation, of the kind discussed in I and involving doing the statistical
averaging[6] over the (random) bonds, we may write the full Debye— Waller factor as
1 .
W = —§ Zln{1— 4:1:(1 — 3:) sinzéKq . Ga - R1m(L% — L31)”. (5.10)
l,m
This would be complex to evaluate numerically as all the lattice Green functions
Gil must be used in the summation. However for reasonable values of the length
mismatch parameter (LOB — L94) / L, the argument of the sine term in (5.10) is small

and so to leading order

1 A
W = 51(230 — $)(Li3 — 13022 :(q'Giz 'le)2
l,m

 

= %((q-u)2) = %(u2). (5.11)

which is the conventional Debye-Waller factor. We have dropped the site label from

(112) as it is the same for all sites. After some manipulation, we find the result,
‘12 2 K 2 0 o 2 2
W = W > = —,-q 2(1— $ng - L.) <1/w >. (5.12)

where (l /w2) is the mean inverse squared frequency of the FCC lattice with nearest
neighbor springs K. The “phonon” band[7] is defined by the frequencies 0 < 022 < 8K,
and we find by numerical integration that (1 /w2) = 2.38K. Thus the Debye-Waller

factor can be regarded as known through (5.12).

5.3.2 Huang Scattering

By subtracting the Bragg scattering (5.8), from the total scattering as (5.4), we find
that the Huang scattering[l] [8][9] is given by

11.01) = exp<—2W)% Zequ- R9,)[(.XP(.-q - m.» exp<2W) — 11. (5.13)

This is as hard as the original expression to evaluate, and so we approximate (5.13) by
expanding the terms in the final square bracket. The approximation (5.7) is replaced

by the improved approximation

(expfiq 7 11:5) = exp -2wl1 + ((q 7 u,)((q 7 111) + 0014)] for i# j;

:1 for i=j. (5.14)

Note that because of inversion symmetry, the correction terms to (5.14) are 0(q4).

Doing the summations, we are led to the result[6]

1H(Q) z 1-(1+2W)exp(-2W)+eXP(—2W)%Zexp(iq-Rfj)((q-ui)(q7ug'))

106

= 1 — (1 + 2W)exp(—2W) + eXPf—2W)((q ' uq)(q ' u—q»

[1 -— (1 + 2W) exp(—2W)] — exp(—2W)K:r(1 — m)(L°B — LoA)q . G(q) - q.

(5.15)

Here uq is the Fourier transform of the displacements u; . The Green function
G(q) is the Fourier transform of Gij , and diverges at small q as q2 . The term
[1 - (1 + 2W) exp(—2W)] in (5.15) is small and accounts for most of the background
away from the Bragg peaks. For this reason, it is important to include the higher
order diagonal terms, even though the higher order off-diagonal terms are neglected
in (5.14). The last term in the expression (5.15) does not lead to large scattering
near the origin as the q'2 divergence is cancelled by the q2 factor in (5.15). However
this divergence remains at all Bragg peaks. Note that this divergence is integrable
because of the phase space factor in three dimensions. We have lumped together
all the non-Bragg terms into (5.15) for convenience, although they could have been
separated into a diffuse background term (the first term in square brackets) and the
divergent scattering (last term). The result (5.15) along the principal cubic directions

becomes
IH(q) = 1 — (l + 2W) exp(—2W) + exp(—2W):r(1 — m)[q(L% — L9,)]2/ sin2 4). (5.16)

where ¢ = qL/x/S in the 100 direction and 45 = qL/\/6 in the 111 direction. The
divergences in (5.16) occur at the Bragg peaks (<15 = n7r, where n is any non-zero
integer). We show the results along the 100 direction in Figs. 5.3 and 5.4 for different
values of the length mismatch parameter (L93 —L94)/L. Although the Huang scattering
falls away from the Bragg peak as (g — q)2 in all directions, it is not isotropic around
the Bragg peak as can seen from (5.15). In some directions, as those shown in Figs.

5.3 and 5.4, the coupling is entirely to the longitudinal phonon modes, while in others

 

107

it is entirely to the tranverse modes. In an arbitary direction, the coupling is to both
kinds of mode. This is quite different from the scattering from single isolated site
impurities,[8] where the spherically symmetric strain field leads to only a longitudinal
coupling. The difference arises because we are substituting bonds and not sites in our
model system. The analysis can be extended to the more complex of concentrated
site substitution.[17]

To obtain the total integrated strength, we make the rough assumption that the
scattering is isotropic around each Bragg peak, and integrate the last term in (5.16)
inside a sphere of radius equal to half the spacing to the nearest Bragg peak in the

same principal direction. This leads to the ratio of intensities
Jam/130;) = A141 - $)[9(L% - Lil)? (5-17)

where A is a constant that is ~ 0.1. This is to be expected from sum rule considera-
tions. If the length mismatch parameter (L93 — L94) / L is increased slowly from zero,
then the intensity of a Bragg peak is decreased by a factor exp(—2W) = 1 — 2W,
and the intensity 2W = .7:(1— 1:)[g(L%— L0A)]K(1/w2) = 0.42.7:(1 -— 2:)[g(L% — L°A)]2
appears at the base of the Bragg peak as the Huang scattering. These estimates are in
agreement apart from a numerical factor of about 4, due to the rough integration that
was done to obtain (5.17). Note that the background term 1 — (1 +2W) exp(—2W) in
(5.16) is 0(q4). The q dependence of the Huang scattering 92 exp(-—(u2)g2/3) means
that it is absolute maximum at g2 = 3/(u2), when the Debye-Waller factor is down
by a factor e_1 = 0.368. However, relative to the Bragg peak, the Huang scattering
keeps on increasing in intensity as q increases as given in (5.15).

The (g — q)”2 divergence of the Huang scattering at the Bragg peaks is due entirely
to the correlations between the displacements u; . When the displacements at each
site are independent, as with “frozen Einstein oscillators”, correcting the expression
(5.7) for the self terms, leads to an intensity 1 —— exp(-2W) for the Huang scattering

which is just a non-divergent background.

108
5.4 RESULTS

To make a direct check of the result (5.12), we have obtained (“2) from the simulation
results. To do this is was necessary to superimpose and position a perfect FCC grid,
of the appropriate size as given by (5.3). The positioning was done by minimizing
Ziu? with respect to rigid motions of the grid to obtain the best registry, and is
accomplished by using the conjugate gradient method. The results are shown in Fig.
5.5, where the agreement with the parabolic composition dependence 3(1 — a2) and
overall magnitude of (uz) is verified. Note that there are no adjustable parameters in
(5.12), as we have calculated the magnitude via the inverse second frequency moment

of the FCC lattice.

5.4.1 Bragg Scattering

We show the result (5.8) for the Bragg scattering 13(q) in Figs. 5.2 and 5.6, where
the Debye-Waller factor is obtained from (5.12). The agreement, involving a total
more than 50 separate Bragg peaks, is very good. Some simulation points lie slightly
above the theoretical Debye—Waller factor, especially at intermediate q values for a 4%
mismatch, as can also be seen in Fig. 5.5. We have tested the size dependence of our
results by varying the size of the supercell and find no significant changes. This small
discrepancy at intermediate q may be real and due to higher order terms in q in the
full Debye-Waller factor as written in (5.10). We note that we have no higher order
anharmonic terms in our potential (5.1) which would also modify the Debye-Waller
factor.[11] We have computed the displacement- displacement correlation function for
nearest neighbor sites, with the result[12] that (Ui’uj‘)/(U2) = 0.08 and 3((u;-R,-,-)(uj-
R4j))/(u2) = 0.40. These two ratios are obtained for all concentrations x and for all
values of the length mismatch parameter (L2a — LED/L. While these correlations do
fall off rapidly with distance, in no sense are the local displacements associated with

the static concentration wave uncorrelated and insignificant. Nevertheless, our results

 

 

109

 

j 1 j ‘ I I 1 WI j
0.00050 — 7-

N
a D d
V 0.00025 - , -

I- I
I- d
I- I

 

 

 

0.00000 775717747
0

Figure 5.5: Illustrating the parabolic dependence of the mean squared displacement
(u2) for a random FCC alloy as a function of the composition. The squares are from
the computer simulation and the solid line is the theoretical result (5.12). The natural

bond lengths are 0.98 and 1.02, so that the length mismatch parameter is 4%.

110

 

ln(IB)

 

l

 

 

 

o 200 400 500 _ 500

0114/2702

Figure 5.6: The same as Fig. 5.2, but shown as the logarithm of the Bragg intensities

plotted against the square of the wavevector. This is often referred to as a Wilson

plot when used to analyze experimental data.

111

show that these correlations have no effect on the Bragg Scattering. __

5.4.2 Huang Scattering

We show the result (5.16) for the Huang scattering Iy(q) in Figs. 5.3 and 5.4 for
two different values of the lattice mismatch parameter. It can be seen that the
agreement between theory and the simulation results is very good, considering that
the approximation (5.14) has been made in the theory. This fit is equally good at all
q values, except that at large q the scattering is weak due to the Debye-Waller factor
(see Fig. 5.2) and the simulations become noisy. This result pleasantly surprised
us as the approximation (5.14) is only good to 0(q2) in the off—diagonal terms. It
appears that the higher order off-diagonal terms are insignificant and non- divergent.

There is no evidence for any rounding of the peak, as would occur if the (g — q)’2
divergence were associated with the wings of a Lorentzian. Equally, there is also
no evidence for additional divergent terms as might have occurred from the higher
order terms ignored in the approximation (5.14). It might have been expected that
differences might have appeared when the lattice misnmatch parameter was as large
as 16%, but the comparison between theory and simulation is about as good in Fig.

5.4 as in Fig. 5.3.

5.5 DISCUSSION

As well as the static effects considered here, the thermal vibrations about the distorted
static lattice, such as shown in Fig. 5.1, also contribute to the observed scattering.

The cross section in (5.6) is modified to give
1 . . .
1(9) = N 265qu - R?.~)<exp(zq 7 uij)><6Xp(zq - “EC-Dr (5-18)
U

where R,- = R? + ui + u? and the u? are the thermal displacements. The thermal

average[13] (...)T is independent of the static average, when both the static and dy-

112

namic displacements are controlled by the potential (5.1). The approximation (5.7),
when applied to the expression (5.18), means that the total Debye-Waller factor is
a product of the static part discussed previously, and a thermal part which has the

form[13] exp(—2W) so that W —> W + WT, where
WT = «W» = hi<1n<w>+ $1M (5.19)

and n(w) = [exp(hw/kBT) + 1]‘1 is the Planck function. All the expressions in the
previous sections must be modified using the replacement (5.19). At low tempera-
tures the thermal part of (5.19) involves only the zero point vibrations, which are
small except for very light atoms. At high temperatures (greater than the Debye

temperature), we obtain the classical result
W. = Z—BMT-u/w?) (5.20)
This is the same (l/wz) that occurs in the static Debye—Waller factor (5.12), except
that the mass M is relevant. For the static distortions, the mass was irrelevant and set
equal to 1. For the FCC lattice the phonon band is now defined by 0 < (.02 < 8K/M
so that the thermal Debye-Waller factor (5.20) is larger than the static Debye-Waller
factor (5.12) if
kBT > 5(1— x)K(L‘}a - L2,)2 7 (5.21)

which compares the thermal energy kBT to the potential energy associated with a
displacement LOB — LOA , weighted with the usual alloy factor :r(1 — 2:). We note from I
that the strain energy per bond 6 = Kzz:(1— x)(L?3— LEO/4 so that (5.21) is equivalent
to kBT > 46. Note that the condition (5.21) is quite general and not lattice specific
as all the lattice effects occur through the common factor (l/wz).

In thermal neutron scattering, the neutron has a low energy so that the “tran—
sit time” for the neutron is sufficiently long that it sees the time averaged lattice.
The phonon motion leads to Bragg scattering modified by a Debye—Waller factor
exp[—2(W + WT)], as discussed in the previous paragraph. The term invoving

113

((q - u?)(q - u$))T that arises when (5.18) is subject to the approximation (5.14),
leads to one phonon inelastic scattering. Thus there is no thermal equivalent equiva-
lent of Huang scattering and all the (g — q)2 scattering around the base of the Bragg
peaks can be attributed to Huang scattering from the static distortions.

In X-ray scattering, the high energy of the X-ray leads to a fast “transit time”
so that the structure is effectively frozen. The X-ray scattering is an average over
all such frozen structures, which means that the diffraction spectrum is an integral
over all frequencies. The quantity W is replaced by W + WT + Wp, where W7 is the
thermal part discussed in the previous paragraph and W1: is the contribution from
the atomic form factor. The Bragg scattering is still modified by the Debye-Waller
factor. The term involving ((q-u?)(q-u;f))7~ from applying the approximation (5.14)
to the thermal part of (5.18), leads to thermal diffuse scattering; the divergent piece
of which has the form[13]

—2w _ -2 fl °° Imlq7G(q,w)-ql
e <<q~u3><q~ui.>>r—-e W./_..——1_exp(—h../k.r>d‘”

At low temperatures (5.22) diverges like [g—ql around the Bragg peaks and is generally

(5.22)

expected to be weak and masked by the Huang scattering. At high temperatures
(compared to the Debye temperature), the thermal diffuse scattering (5.22) becomes

— e'szBT q - G(q,w = 0) ~q/M _ (5.23)
This has the same form as (5.15) with the two Green functions being identical except
for a factor M, so that the thermal diffuse scattering is like the Huang scattering
shown in Fig. 5.3, including the background terms. The condition that the divergent

(g — q)‘2 thermal diffuse scattering is greater than the Huang scattering is again

(5.21).

5.6 CONCLUSION

We have shown that the elastic scattering from a model alloy consists of two parts; the

Bragg scattering and the Huang scattering. The Bragg scattering is modulated by a

114

Debye-Waller factor and the divergent Huang scattering comes entirely from the sec-
ond order displacement-displacement correlations. Our central approximation (5.14)
has been shown to essentially numerically exact, and can be used with confidence in

more complex alloys and geometries.

Bibliography

[1] K. Huang, Proc. Roy. Soc. A190, 102 (1947).
[2] M. F. Thorpe and E. J. Garboczi, Phys. Rev. B42,8405(l990)

[3] See for example, R. M. White and T. H. Geballe in Solid State Physics, (Aca-
demic Press, New York, 1979) Vol. 15, page 45.

[4] W. H. Press, B. P. Flannery, S. A. Teukolsky and W. T. Vetterling, Numerical
Recipes (Cambridge University Press, Cambridge, England, 1986) page 301.

[5] See for example, N. W. Mermin and N. D. Ashcroft, Solid State Physics, (Saun-
ders College, Philadelphia, 1976), page 790.

[6] The most important result from the statistical averaging that we use in this note
is (ugauJ-g) = —K;r(1 — x)(L°B — Lara-3}“. All the lattice Green functions G' used
in this chapter are evaluated in the static limit. The Green function for the FCC
latice is defined as the inverse of the dynamical matrix and is given in E. J.

Garboczi and M. F. Thorpe, Phys. Rev. B32, 4513 (1985).

[7] L. von Heimendahl and M. F. Thorpe, J. Phys. F L87 (1975). Note that the

masses of the atoms are irrelevant for statics and we set them equal to unity.

[8] M. A. Krivoglaz, The Theory of X—ray and Thermal Neutron scattering from
Real Crystals (Plenum, New York, 1969).

115

116

[9] H. E. Cook, J. Phys. Chem. Solids, 30, 1097 (1969); W. Fenzl an_d S. C. Moss,
J. Phys. F17, 1285 (1986) and references contained therein.

[10] A. J. C. Wilson, Nature 150, 152 (1942).
[11] A. A. Maradudin and A. E. Fine, Phys. Rev. 128, 2589 (1962)

[12] Using the result in l for the total strain energy, it can be shown that the ratio
3((u; - 1‘1.,)(u,— - 11.,» /(u2) = 1 — §(1 — a‘)/(K/w2) = 0.40 for the FCC lattice
where a' = 1/2. These ratios are independent of the concentration x and the

length mismatch parameter (Lg; — LBJ/L.

[13] The required thermal averages are given by (125,213,, = —% 3°00 I mGS-B (w) / [1 —
exp(—hw/kBT)]dw, as discussed for example by M. F. Thorpe, Phys. Rev. B8,
5352 (1973).

Chapter 6

TWO DIMENSIONAL MIXED
CRYSTAL

This paper is published in Physical Review B. Reference: M. F. Thorpe and Y. Cai,
Phys. Rev. B43, 11019(1991)

6.1 INTRODUCTION

The melting transition has been extensively studied in single component two dimen-
sional systems.[1]-[5] It is known that thermal fluctuations destroy the long range
positional order at any finite temperature. This produces a phase with quasi-long-
range-order that is characterized by a power law decay in the position-position cor-
relation function at all temperatures up to some temperature T1. This phase has
power-law peaks in the diffraction pattern, where the Bragg peaks would have been.
However the long range orientational order remains all the way up to T1 above which
the angular correlations decay algebraically, in the hexatic phase.[5] Finally at some
higher temperature T2, melting is completed, and both the positional and angular
correlations fall off exponentially.

In this paper, we investigate the idea that the size—mismatch between the two

117

118

components in a two dimensional random alloy Aqur also destroys the long range
positional order in an analogous way to the thermal disorder in a single component
system, as previously discussed by Nelson.[5] We set up a triangular network with
nearest neighbor forces with spring constant K and natural spring lengths L2, and
L93. The topology of this network cannot be changed and every atom always has
the same six nearest neighbors so that the disorder is quenched in. Nevertheless, the
long-range positional order is destroyed, even though the fixed topology prohibits the
formation of either disclinations or dislocations.[5] Our approach is analogous to the
spin wave description of the thermal disordering of the classical X Y model.[6],[7]

We show that it is possible to define an effective temperature TD associated with
the disorder via,

kBTD = Kx(1— 5mg — L3,)2 (6.1)

where L93 — L2, is the length mismatch. The triangular lattice is shown in Figs.6.1
and 6.2, with two kinds of bonds A and B that are randomly located in the lattice
with probabilities l — :1: and m respectively. This mathematical simplification to bond
disorder rather than site disorder does not qualitatively effect any of our results.
For example, if the lattice is made up of randomly positioned A and B sites with
probabilities 1 — m and .2: respectively, then the disorder temperature in Eq. (6.1) is
replaced by

kBTD = $1150 — 3:)(LOBB — L3“)? (6.2)

where L933 and LEM are the lengths of the AA and BB bonds, and the AB bonds
have natural lengths LOAB = %(L°BB + LOAA).[8] Other models for the length disorder
would lead to more complex expressions for T D.

In the next section, we recall that the short range properties are entirely conven-
tional and similar to those in higher dimensions. In section 6.3, we show that the
decay of the positional correlations is controlled by the parameter 77, which leads to
power law peaks in the diffraction pattern if 7] < 2, with no peaks if r] > 2. In section

6.4 it is shown that the angular correlations do display long range order which is cal-

. —

119

A A A A A
VAVAVAVAVAVA
AVAVAVAVAVAV
VAVAVAVAVAVA
AVAVAVAVAVAV

VAVAVAVAVAVA

AVAVAVAVAVAV

VAVAVAVAVAVA
V V V V V

 

Figure 6.1:

The regular triangular lattice.

120

l i I -- no r---
I \

s
l 2“ \| .
....- .--?“ \ \‘ ‘
\. I 2“.
s l \ I \ I \
‘0' ‘ l 0 , \
on \‘ , \ v \
.09. ’
'H s f
\ I
\ a
I
.../A '
g 1
a
--- .. s _
Y'w \‘ a
‘9‘ I.“ \ I,
I

‘Q
\ "~ ...--
€ \

b
0

vs
,a

 
  

\‘

 

Figure 6.2: The relaxed triangular lattice with two different kinds of bonds with

different natural lengths. The concentration :1: = 0.5; the shorter bonds are shown

as dashed lines, and the longer bonds as solid lines. The length mismatch parameter

(L9; — L9,)/(L) is 30%.

121

culated. Throughout this paper we compare our results with computer simulations.
In sections 6.3 and 6.4, we show that the positional order is finally destroyed when
the lattice pleats.[9] The pleating is not described by our linearized approach, and is
brought about by the fixed topology that prevents a more reasonable disordering of
the lattice at a smaller length mismatch. Nevertheless, the pleating is of some interest
in its own right as it relates to the crumpling transition;[10] the difference being that
no motion out of the plane is permitted in pleating whereas it is permitted in crum-
pling. We stress that the pleating transition is unphysical, because of the absence
of any repulsive forces at short distances in our model, as well as the imposition of
a fixed topology. The pleating transition takes place at a length mismatch of about
50% which is way beyond the region of physical interest as no random solid solutions
can tolerate length mismatches of more than about 15%.[11] In section 6.5 we discuss
finite size scaling and in the conclusions, we explore ways in which the model needs
to be modified to give a more realistic description of size-mismatch effects in two

dimensions.

6.2 SHORT RANGE PROPERTIES

The triangular network is described by a potential
V = €23ng — ij)2 (6.3)
(u)
where K is the force constant between nearest neighbor atoms, denoted in the sum-
mation by the angular brackets. The length of the bond ij in the network is L5 and
the natural (unstrained) length of this bond is Lg. The bonds are chosen randomly
on a triangular network to be either type A with probability 1 — x or type B with
probability 3:.
The short range properties of this model have been well studied and are quali-

tatively the same in two and three dimensions,[9],[12] as summarized below. The

122

network is relaxed so as to minimize the strain energy. This gives fairly broad dis-
tributions of the A and B bond lengths, from which the mean lengths are given
by,[9l
:c
(L..) = L1+ 3w}. — L9.)

1 $03. — L3.) (6.4)

 

(L8) = Lia —

where the angular brackets (...) denote an ensemble average. The mean overall length,

which determines the size of the sample obeys Vegard’s law[13]

(L) E (1 - $)<LA) + x(LB)

= (1 — 2:)L21 + ng. (6.5)
The widths of the A and B distributions are equal and given by

(Li) — (L2)? = (L13) — (LB)?

2
= 52:0 — rm. — L1)? (6.6)
The total strain energy E, can also be found, and is given by
E, 1
N = 519(1— $)(Lfa — L302
1 .
= ZkBTD (6.7)

where N is the number of sites and the disorder temperature TD has been previously

defined in Eq. (6.1).

6.3 DIFFRACTION PATTERN

In three dimensions, the model (6.3) on say a face centered cubic lattice, leads to con-

ventional Bragg peaks and associated Huang scattering in the diffraction pattern.[14],[12].

123

The situation is very different in two dimensions. Assuming the scattering length of
all atoms to be unity, the elastic scattering cross section for neutrons can be written
as
1(0) = i723 e‘Q'Rg(t7:“2""*‘>- (6.8)
u
where Q is the momentum transfer. In this equation, R? is the equilibrium position of
atom i, and the average displacement (ui) is required to be zero to define a reference

lattice with nearest neighbor distance a = (L) as given by Eq. (5). The approximation
for the average needed to evaluate (6.8), up to terms 0(Q“) is,

(69%) = Cq(R) 2 (“(8%”) (6.9)

This is defined as the density-density correlation function Cq(R), where R is distance

between sites i and j. From our previous work[9],[15]
(uin) = —K:z:(1 — z)(LOB — L9,)2Gu (6.10)

which is similar to the relationship for long wavelength phonons in a single component
lattice at a temperature T, if we were to replace Km(1—x)(L9, — LOB)2 by 1337". Using
(6.10) we have,

1 K
509-115?) = 7x0 — x)<L9. — L93)? [Q - (0.. + ij — LG...) - Q] (6.11)
The Green’s function can be written in k—space using the Fourier transform,
515110

Gii-l-ij—ZGij = %ZG(k)(l —C 1J). (6.12)

For small k, we may write,

G(k) k2 (6.13)
where, for the triangular lattice,[16]
8
AT = 3AL = (6.14)

3a2K

124

and a is the lattice constant. Evaluating the integrals, we have _

 

2‘“? ' "02> = K10 — 201909. — 2.1252,} AT + W 7.2/{Tm ' k) (1 — 27‘7“?»

fa

 

                      

= Kx(1— :r)[Q Q(L° — L0 )128 (6.15)
where we have used the integral,
112119.
— U
/(1—:2—)d2k = 27rln R2]- + const. (6.16)

which is valid for large R9]. The logarithmic distance dependence in (16) comes from
the integrand around the origin, where the approximation (13) is valid. This leads us

to the result,
{Hm-nu?) = 302)

 

 

( R91)" (6.17)
and hence
CQ(R) = 8],?) (6.18)

where to leading order we can ignore the difference between R9]. and R. The exponent

1] given by,

n = 1\$(1— $)[Q(L° — L0 ,0ng

         

)

=3——} 2(19231—2)[Q( _ng (6.19)

The correlation function Cq(R) has the same form as that derived by Nelson,[5] with
Lamé constants /\ = ,u = x/3K/4 and the real temperature replaced by TD from
Eq. (6.1), plotted in Fig.6.3. The results are averaged over the six equivalent nearest
neighbor bond directions. Both the longitudinal and transverse distortions contribute
to 17 in (19) through AL and AT. A similar result to Eq. (6.19) for site rather than
bond disorder leads to a result for TD as given in Eq. (6.2). In the site case, the

distortions around each site have radial symmetry, which means that only the long

 

125

 

0.95

   

17(Q2)=0.0348

0.85 7

 

 

Figure 6.3: The density-density correlation function Cq(R) plotted against the dis-
tance R in units where the mean bond length is unity. The results are averaged over
the six equivalent nearest neighbor bond directions. The length mismatch parameter
(Lg — L9,) / (L) is 4% and the results are averaged over 4 samples, each of which is
100 X 100. The solid line is using the theory in Eqs. (18) and (19). The upper line

is for Qla = 3717 and the lower curve for Q20 = 4r; neither of which are reciprocal

lattice points.

 

126

wavelength longitudinal distortions contribute to 7] through A L.[8] The algebraic decay
of the correlation function (6.18) is checked against computer simulations as shown
in Fig.6.3, where the amplitude B (Q) is treated as an adjustable constant. We can
see that the theoretical approximation (6.18) is adequate.

The disordering of the perfect lattice is controlled by the parameter 1; which is
determined by the size-mismatch rather than the temperature.[5] From the result
(6.18), we are able to calculate the asymptotic form of the I (q) around the reciprocal
lattice vector at g=Q-q, where q is measured from the nearest reciprocal lattice vector
g. Divergent behaviour is found for peaks with n < 2. However, contrary to previous
claims[5] there is no cusp when 7} > 2. Combining Eqs. (6.8), (6.9) and (6.18), we

have,

00 21r ,
I(q) = I. [0 e'qR°°79(%)"Rdeo
= /°° Jo(qR)(%)"RdR

= a”q”—2 /°° J0(:c)xl"’d:c (6.20)
go

When n < 2, the kernel of the integral is well behaved at the low cut-off which can be
replaced by zero. The integral (6.20) contributes a multiplier and I (q) is divergent at
small q as q"(2"’) where q is measured from the reciprocal lattice vector. When 17 > 2,
we can not replace the low cut-off in (6.20) by zero. Using the fact that Jo(:r) 2 1
for small x, we see that the lower cut-off gives a constant contribution a2/ (2 — 7]),
and so we obtain an analytic result for I (q) This suggests no Q-dependence in I (Q)
and hence no cusp when 17 > 2, which is verified by our simulation results shown in
Fig.6.4.

The pleating transition occurs at a length mismatch of about 50%, so that the
pleating is very apparent at a length mismatch of 60% as shown in Fig.6.5. Just
before a 50% length mismatch, a few sites begin to pleat and the number increases

rapidly as the size-mismatch exceeds 50% as shown in Fig.6.6.

 

 

I ' II 7 - 10

i <—- —> J
& _qy qy 8

 

MI)

10

 

 

 

 

Figure 6.4: The lower panel shows the diffracted intensity I(Q) plotted against the
momentum variable Q3, where the a: direction contains a nearest neighbor bond, and
the results are averaged over the six equivalent directions. The length mismatch
parameter (L93 — L9,)/ (L) is 4% and an average is taken over 87 samples of 40 x
40 triangular lattices. The upper panel shows the diffracted intensity around the

reciprocal lattice points in the lower panel but in the perpendicular direction.

128

 

 

 

 

Figure 6.5: The relaxed triangular lattice after pleating has occurred. The length
mismatch parameter (L93 — LED/(L) is about 60%.

129

 

 

 

 

 

 

1.0 . . .
. O
0.8 - ?&
0.6 7 . ‘8’
0.4 ‘ .1
0.2 7 o
o
0 0 1 ° ° . v 4 4 0 10
o o - 0.08
J 0.08
LL ° - 0.04
° 7 0.02
4‘: ¢ :7, e " 0 O i : 0.00
1.00 ——-l——§—-i——g—¢—u-——n-—x-—-x——-x——-
o o o
0.75 7 o . /\
. o o . r7
0-50 7 '- — - Mlero v
0.25 ’ 0 Macro .
0.00 4 7 - 1
0 0.2 0.4 0.6 0.8 . 1

LENGTH MISMATCH

Figure 6.6: Showing the long range angular correlation C9(oo), the pleating P, and
the microscopic and macroscopic lengths (L) plotted against the length mismatch
parameter (L2; —— Li)/ (L). The symbols are from the computer simulations and the

dashed and solid lines are the theoretical results from Eq.(6.31) and ref.(6.5).

 

130

A pleated site is defined to have all six bonds to its nearest neighbors lying to one
side within the 180° arc of a semicircle. The network chooses this form of disordering
for large size-mismatches because it is prevented from forming dislocations and discli- ’
nations as occurs in the thermal case.[5] In the thermal case, the solid begins to melt
as first dislocations and then disclinations begin to unbind, destroying the topological
order. However, the topology is preserved in our model even when pleating occurs

because the six nearest neighbors of a site remain fixed.

6.4 ANGULAR CORRELATION

Long range angular correlations are present and can be calculated within our model.

The angular correlation function is defined as,
C9(R) E [(6690)] (6.21)

where 9;,- = 0,- — 0,- is the angle between two bonds emanating from lattice points i and
j, and R is their separation. It is convenient to measure the 0,- from the underlying
reference lattice, but it could also be measured from the line joining the sites i and j;
the result would be the same to leading order. Because long range angular correlations

exist, for large R,
(6‘60”) __) (61.69; > (e-i60j> (6.22)
and,
e-Ws s |(e-7'69-7)| : 67%«6900 (6.23)
and only the deviation 66 = 0,- — nir / 3, where n is integer, is needed

(“12 X Riz)

(L)

and (L) is the mean bond length. This gives us,

60:

 

(6.24)

_ 180112 - (11127312?)

W6 (1)2

 

 

 

131

= —36Ka:(1 — :r) [EC—L051] Tr [(1 — R‘I’ZR‘I’ZXGH — 0.12)] . (6.25)

(L)
where we have used Eq. (6.10). We know that R92 - (G11 — G12) - R92 = —d/z,

where d is the dimension and z the coordination number.[17] For the triangular lattice,

d/z =1/3, so that,

W6 = —36K:c(1 — ..) [WY [Tr(G11 — G12) — g] (6.26)

where

l 1 l
= —Nzk: [W + m] (1— 7k) (6.27)

and (.02 = D1 :h x/Dz. The expressions needed are,

7(k) = %(cos 2.2: + 2 cos 2 cos g) (6.28)
D1 = 3 — 2 cos 2x — 2 cos 2 cos y (6.29)
D2 = (cos :6 cos y — cos 22:)2 + 3 sin2 :1: sin2 3; (6.30)

where a: E kra/2 and y E \/3kya/2. The k-integration in (6.27) is done numerically

to give,

W6 = 17.22:(1 — x) [WY (6.31)

and for large R gives,
09(00) = e-2W6. (6.32)
The result (6.32) shows that the long range angular order always exists in our

linearized calculation. When pleating occurs, the long range angular correlation is

132

 

 

 

 

1.00” o o o o——o—:§—+—+——+—+——+—o—-¢—-+-—
D a 1:1 HW—a—e—a—
10%
0.75“ 207. ..
* + + r—r—LFTT-v—T-w—r—T-r—
A .
n: 0.50) x 307: -
v
a n x n x x x x x: it"? n x n
D
0.25’
50%
5: 5 5c 57 e-
57 e
0.00 7 ‘3' f 5' 6' 7 f f “'

 

 

 

0 5 10 15
R

Figure 6.7: The angular correlation function C9(R) plotted against R in units where
the mean bond length is unity. The results are averaged over the six equivalent direc-
tions defined by the nearest neighbor bonds. The results were obtained by averaging
over 9 samples each of which was 40 x 40. The solid line is the asymptotic theory and

the various values of the the length mismatch parameter (L93 — L94)/ (L) are indicated.

133

destroyed as shown in Fig.6.6. We have performed simulations for different length
mismatches and the results are shown in Figs.6.6 and 6.7. When there is no pleating,
our analytic results agree well with the computer simulations. When pleating occurs,
the computer result begins to drop below the prediction. The top panel in Fig.6.6,
shows a comparison between the calculated long range angular correlation Cg(oo) 7.:
exp(—2W6) and the simulation results. The agreement is very good up to about a 30%
length mismatch, when effects due to pleating first start to have an effect. The theory
predicts that the long range angular correlations are always present as no account is
taken of pleating. In the second panel of Fig.6.6, we show the fraction of sites P
that are pleated. This gives the clearest signature of the phase transition, although
there is rounding as a result of finite size effects. The fraction of pleated sites may be
thought of as a disordering parameter. In the bottom panel of Fig.6.6, we show the
microscopic length and the macroscopic length from simulations. The microscopic
length, defined as the average bond length, is constant as predicted by the theory
(6.5), while the macroscopic length, deduced from the sample size, decreases above
the pleating transition as would be expected. In pleating, the network folds over-on
itself, in the same manner as a pleat put in a piece of material, and hence reduces

the area.

6.5 FINITE SIZE SCALING

Because of finite size effects, the Bragg peaks are never completely eliminated in
our simulations. This is because of the artificial long range order imposed by the
periodic boundary conditions. Nevertheless it is convenient to keep these boundary
conditions as then every atom has exactly six nearest neighbors. We need to estimate
the magnitude of this effect and try to extract the infinite size behaviour from the
finite size simulations. This is easy for some quantities, like C9(R), and very difficult

for others, like I (Q) The scattering at the Bragg momentum g is, from Eqs. (6.8)

134

 

and (6.17), --
1 8(6)
I(s) = -
N )6: (Rt-)7
L (PR
o< / — (6.33)
For n > 2 and L >> a, this leads to,
I(g) o< L2"7 (6.34)
and because L2 o< N, we may write,
E oc L'” (6.35)

N

This remnant Bragg peak contains an extra factor N because it is a Kronecker delta.
rather than a Dirac delta function.[12] Otherwise the Bragg peaks vanish as the
system gets larger according to Eq. (6.35). This finite size scaling law can be verified
by computer simulation and we have roughly verified that it does work using rather
small samples with L _<_ 200. The Bragg peak of the perfect lattice (n = 0) has a
height proportional to N. As disordering is introduced by increasing TD, and hence
17, this weight is redistributed locally. In three dimensions, it becomes the Huang
scattering,[12] while in two dimensions it either becomes a power law peak if 17 < 2 or
background if 77 > 2 as shown in Fig.6.4. For large Q, the scattering tends to unity as
the phases in Eq. (6.8) become random. Note that there is always a Bragg peak at the
origin with weight N that represents a Dirac delta function with weight unity. When
L’” << 1, the scattering profile approaches that of the infinite sample. In Fig.4.4
the Bragg residue contains 30% of the weight of the first peak, while the second peak

has only a 0.75% Bragg residue and is therefore representative of an infinite sample.

6.6 CONCLUSION

We have shown that the long-ranged properties of a two dimensional mixed crystal

can be described by a linearized small displacement theory, as long as the topology

 

135

is kept fixed and no disclinations and/or dislocations are allowed. This is achieved
naturally in a spring model with a fixed triangular network topology. While there is
always long range orientational order present, the positional correlation function falls
off algebraically, leading to power law peaks in the diffraction pattern.

This situation is similar to the thermal disordering of two dimensional lattices, and
indeed many of our results can be understood if the temperature T in the thermal
case is replaced by the disorder temperature TD as defined by the strain energy in
Eq. (6.1) or (6.2). This model can be regarded as the analogue of the quadratic
02 model of Zittartz[6] and José et al.[7] for spin wave excitations in the classical
X Y model in two dimensions. Like our model, the 02 model does not appear to
contain any disclinations or dislocation type excitations, while successfully accounting
for the instability against long range positional disorder and for the initial thermal
disordering of the lattice. Our model can be regarded in a similar way. A more
complete treatment would require abandoning the constraint of having six defined
nearest neighbors, and using pair potentials for the AA, BB and AB interactions
between the two kinds of atoms A and B in the random alloy A1_3B,,. The repulsive
part of the pair potential will completely inhibit the pleating, but the possibility of
switching nearest neighbors leads to the flexibility of lowering the energy by creating
disclinations, dislocations etc. Taking the thermal analogy seriously, such excitations
would be expected to be present in the ground state as TD is increased.

Thermal effects can be incorporated into our model at low temperatures. If the
thermal and static displacements are controlled by the same force constant K, then

the thermal and configuration averages can be done separately to give,

2 4q2kB(T + TD)
3x/37rK

In the limit that TD goes to zero, the usual thermal result is recovered.

 

(6.36)

Bibliography

[1] J. M. Kosterlitz and D. J. Thouless, J. Phys. C: Solid State Phys. 6, 1181 (1973)
[2] D. R. Nelson, Phys. Rev. B18, 2318 (1978); 327, 2902 (1983)

[3] D. R. Nelson and B. I. Halperin, Phys. Rev. B19, 2457 (1979)

[4] K. Standburg, Rev. Mod. Phys., 60, 1 (1988)

[5] D. R. Nelson, in Phase fiansitions and Critical Phenomena, Ed. by C.
Domb and J. L. Lebowitz, 7, 1 (1983)

[6] J. Zittartz, Z. Physik B31, 89 (1978)

[7] J. José, L. P. Kadanoff, S. Kirkpatrick and D. R. Nelson, Phys. Rev. B16, 1217
(1977)

[8] J. Chen and M. F. Thorpe, unpublished
[9] M. F. Thorpe and E. J. Garboczi, Phys. Rev. B42, 8405 (1990)

[10] M. Paczuski, M. Kardar and D. R. Nelson, Phys. Rev. Letts., 60, 2638 (1988)

and references contained therein.

[11] W. Hume-Rothery, The Structure of Metals and Alloys, Institute of Metals,
London, 1936

[12] M. F. Thorpe, J. S. Chung and Y. Cai, to be published in Phys. Rev. B.

136

137
[13] L. Vegard, Z. Phys. 5, 17 (1921)
[14] K. Huang, Proc. Roy. Soc. A190, 102 (1947)

[15] Y. Cai, J. S. Chung, M. F. Thorpe and S. D. Mahanti, Phys. Rev. B42, 8827
(1990)

[16] E. J. Garboczi and M. F. Thorpe, Phys. Rev. B32, 4513 (1985)

[17] S. Feng, M. F. Thorpe and E. Garboczi, Phys. Rev. B31, 276 (1985)

Chapter 7

OVERVIEW AND
PERSPECTIVE

This chapter is to be submitted to Physical Review B. This preliminary version is
written by the author and the paper will be published under the names of Y. Cai and
M. F. Thorpe.

7 .1 LENGTH MISMATCH PROBLEM

We summarize the length mismatch problem[1] in this section and discuss the problem

from a mathematical point of view in the next section.

7 .1.1 Quaternary Solution and Sublattice Decoupling

Given small displacement, it is possible to expand it in a quadratic form,

1
E = §U+MU + U+¢ + E0 (7.1)

U is the displacement field vector and M the connectivity matrix. For a pure crys—

tal, d) is usually the external force. For an alloy system, (6 originates from different

138

139

types kind of disorder, such as disorder in the interaction strengths and lattice con-
stants(natural lengths). It is then interpreted as internal strain. When the mismatch
between natural lengths is small, the harmonic expansion is always valid.
The displacement field when relaxed minimizes the distortion energy, so we have
the field equation:
MU = 45 (7.2)

The pure “Length Mismatch Problem” involves only disorder in the natural lengths
so that the matrix M is free of disorder and only <13 contains randomness. Technically
this field equation can be solved by inverting M to obtain the Green’s Function G

and then evaluating G in momentum space,
U = G65 (7.3)

We will demonstrate this by solving the pure “Length Mismatch Problem” for a
quaternary system A1_zB,C1-yDy using the Kirkwood model where the disorder
vector (15 is linearly proportional to the length mismatch. The assumption needed for
obtaining an analytic solution is that the natural lengths of the chemical bonds are
additive.

The quaternary compounds A1_$B,,C1_,,Dy are interesting in the semiconductor
industry because of the two tunable parameters that can satisfy the requirement of
matching the substrate lattice constant and designing the band gap of coating ma-
terial. One example is I n1_xCaxP1_yAsy. The compound takes diamond structure.
There are two FCC sublattices in the diamond lattice. Ga and In occupy one FCC
while P and As take the other. Therefore in terms of A, B, C and D, there is no
nearest neighbor connection between A, B and between C, D, and they take ap-
propriate concentrations :1: and y in these different sublattices. This system is very
general in the sense that we can reduce other alloy system into this form. For example,
521-30% can be written as a quaternary as Si1-,Si3Gel_xGex. And ternary system

such as I n1-xCa3As can be represented by Asl_yAsyI 121-30% and so on. The idea

 

 

 

 

140 i

is applicable to two dimensional system also, e.g., by decoupling a triangular lattice
into three equivalent sublattices.

The sublattice decoupling actually transforms a system of less components to a
system of more components. The trade-off to the increase of components is the decou-
pling of randomness between nearest neighbors so that the statistics can be performed
independently in two sublattices. It is very useful when short range properties like
the distribution of nearest neighbor distance is concerned. It may not be applicable
to the problem of long range properties such as diffraction patterns because ‘we need
to decouple the system to such a degree that within the “concerned range” no lattice
points belong to the same sublattice. The importance of decoupling of randomness
will be seen in the following derivations.

We can calculate the distribution function of the nearest neighbor distance, from
which the specific average lengths(first moment) and fluctuation(second moment) can
be derived. We will proceed to give the full derivation of the distribution function.

t “0” denotes natural

We first make rules for the notation. Length with superscrip
length. Roman letter i, j, denote one of the sublattices and Greek letter a, [3,
denote the other.

Denote the atomic radii of atom A, B, C and D as rA, r3, re and U). The additive

of atomic radii means,

he“: = TA + 7‘0
ho — r +
AD — A 1"D
ho — r
BC —7 B + TC
11.ng = r3 + 17;) (7.4)

Vegard’s law is valid in the pure length disorder case. It means that the NN-

141

distance of the underline virtual crystal, or the mean NN-distance of the disorder

system, is,

h. = (1 - a=)(1- M. + (1 - 7"?)3/(1310 + 950- y)h°3c + whim

= (1 — 2:)rA + mm; + (1 -— y)r0 + yrD (7.5)

Now we introduce the projection operator on the site 2' as 0.7, and the same for the
other sublattice by replacing i by a. a is an 1d variable, i.e., when a = +1, the site
is A(or C), and when a = —1, the site is B(or D). Given the concentration in two
sublattices, (0;) = 1 — 2a: or (070, = 1 — 2y. The ensemble average is independent of i
or a.

The natural length of any bond in can be expressed in terms of the projection

operators and the atomic radii,

 

 

 

 

 

 

 

 

1+0; 1+aa 1+0.- 1—00,
ha. 2 ° 2 (77910 + 2 ' 2 him
1—0, 1+0a 1—0; 1—0‘0,
1 3 1" i 1 a 1- or
= 2072+ 2073+ +20 Tc+ 20 TD
= hm + AABUi + A0006. (7.6)

where hm = (rA + r3 + rc + rD)/2 and AAB = rA — r3, ACD = r0 — rD.
The length distribution function of certain kind of bond is defined as,

(h) -— (1 +251“ . 14.2.2006“ — h..,)) (7.7)

the 6 can be +1 for A(or C), and —1 for B (or D). So there are four combinations of
6162 for AC ,AD,BC and BD respectively.

 

P

C1C2

Replacing the 6-function by a q integral, we have,

1+ea; 1+eaa1 +°°,~ _.
P(1¢2(h) = < 21 . 22 _2_7_r_‘/;oo Cq(h haa)dq>

 

142

 

1 +00
= eiqh
—/:2.F7172(q)dq7 (7.8)
where
l+eag 1+ec7or _,-.
F....(q)=( ,‘ - ,2 e 7’77) (7.9)

Now we need the expression for hga which is supposed to be solved from Eq.(7.3).
We will utilize the Kirkwood model as an example to solve the length mismatch
problem here.

Kirkwood model with only natural length disorder takes the form,

E = EZUIU- 2+ )g-(hc)2 2((308 0.3); + 7:1?)2 (7.10)

2(0) (51*)
In this equation i, j, k can be in either sublattices. The force constant B has the same
unit as a in the above expression. The expansion form of Kirkwood model fits the
field equation of length mismatch problem in the following first order approximation,
namely,
hia = he + 120.7111... (7.11)

This linear approximation is checked to be good by computer simulation for small
mismatch. The second order terms of expansion contributes to the disorder in force
constants, which is not important in the length distribution when the disorder is
not large. The M matrix is free of disorder, thus is the same as for the crystalline
structure. Disorder only appears in the d) vector. It can be written as,

= —a 20100,)[2'0 > (7.12)

(for)

where [ia >E [i > —[a >.

So the he, is,

hi0: = he + ria ' uioz

= he —QZBGia,j,B(he 42273)

143

= he + a Z Gi3,j:Bh-15 (7.13)
.06

any element of the green’s function in the shorthand notation above is defined as,

G6,“! 5 2,737 7 (< i[— < j[)G([k > —[l >) 7 I‘M (7.14)
and
z GiagB= (7.15)
(J16) -

because an uniform displacement keeps the system invariant.
Substitute Eq.(7.6) into Eq.(7.13), and using the properties given by the sum-rule
Eq.(7.15), we obtain,

hi0: = he +777 3:9,;jfi AABUj‘l'ACDUB)
2019)

= he + 2 9.13.43 2 G13 17ij + 72-Acp Z Gia 1:73:73 (7.16)
(J3) (J3)
Note that we have decoupled the random variables a; and 00, so that the ensemble

average can be perform in the each sublattice separately. Denote

A7,- : 5131132975376
30)

A018 = 2ACD%G‘OW (7.17)
.7

where j(B) means B is the nearest neighbor of j.

We define the lattice integral a" as

a" E z Gga’jg, (7.18)
(3(J')
so,
A1; are:
—-— = —-—A 5 la“ (7°19)

144

a" is one of the lattice integral we encounter most frequently. _
With the above notations, and making use the fact that the random variables in
two sublattices can be averaged independently, from Eq. (7.9) and Eq. (7.16) we

have,

F¢1¢2(9) = 6‘4.th <1 +2610}. 6“? xi A”) <1 +2610}! C—iq 2,314.”) (720)

We only have to evaluate one of the average in the above expression and the other

is the same.

1 i --i . "a-
Ic(9) E <——-+2€ae 921A" '>

 

2 .

J

1 g . .
= < + 60' H(cos quj — :0,- sm qA;,-)> , (7.21)

where we have made use of the identity 0? = 1. So,

.7

1 . .
1.01) = '2' {(H(COS inj — 20:“ SID quj»
+6(H(COS qA,-,- — i0,- sin qA,,-)) - ((0,- cos (1A,,- — isin qA;;))}(7.22)
1?“
Notice that the product is over distinct sites and so expanding, we will have to average

quantities such as (ammucn), which is just (0‘1)(am)...(a,,) in the random case. This

means that we can factorize all such products. Using (0,) == 1 -— 22:, we have

1 . .
1.01) = 5 [11(008 924.1 — 2(1 - 233) Sin (124.7)
#z‘

{ cos qug — i(1 — 2x) sin qu'i + 6[(1 — 2:13) cos qui — isin ini]} (7.23)

Making use of the fact that 62 = 1, we can rewrite 1.;(9) as

 

1 1 — 2
[((q) = [ + 6(2 3)] (cos qA,-,- — ie sin qug) H[cos quj -— i(1 — 2x) sin QAij]
j¢i

145

 

= [1 + 602— 27)] e'i‘qA“ H[cos quj — i(1 —— 22:) sin inj] . (7.24)
1'?“

Substituting this result in the distribution function and regrouping them properly,

we obtain the solution,

_1+€1(1—2:L‘) 1+€2(1—2x)
_ 2 2

 

Fc1C2(‘1) f (q)

e{_.-q[h.+(.1-1+2:)A,,Ba”/2+(c2—1+2y)Acoa“/2)l}, (7.25)

The function f (q) is independent of the species and has the form,

“9) = gal-27”“ H [C05 inj " ‘11 — 23) Sin injl
17?“)
oe'i(1—2y)q’4°° H [cos qAag — i(1 — 2x) sin qAag], (7.26)
5050)

Note that f (q) is independent of 61 and 62 which specify the type of chemical bonds.
From Eq.(7.25) for F¢,¢,(q) we conclude that the distribution functions for all four
specific bonds are exactly the same for all given concentrations. The weights of
four peaks are (1 — 2:)(1 — y), (1 — :c)y, x(1 — y) and my for AC,AD,BC and BD
respectively. Computer simulation result checks with this as shown in 'Fig.7.1.

From the distribution function we derive the average of specific nearest neighbor

distance as,
(hm?) = k8 + a"[AAB(€1 — 1 + 22:) + ACD(62 — 1 + 2y)]/2 (7.27)
The width of the distribution are the same for all,
(hf... - (hwy) = [95(1 - $)Ai43 + 31(1 - y)A2cpl(ai‘ - a”2) (7-28)
where constant a? is related to a",
(NW/CY)

a? = —a2_—5Z!— (7.29)

146

 

I I l I

Figure 7.1: Four N N -length distribution of a Quaternary compounds, A1-,Bzcl-yDy.

l

l
l

V

L

l
l

l

l

1

fl

 

 

 

0.00

l
A
4

:1: = 0.6 and y = 0.45. The four peaks are exactly the same in shape with proper
weights 0.27 : 0.33 : 0.18 : 0.22. The symbols are from computer simulation using
Kirkwood model with fl/a = 0.2 and the solid line is a gaussian fit use the centers

and widths given in the text.

147

We can also evaluate the averaged next nearest neighbor distance using the projec-
tion method. We give below the full derivation.
12,-,- is the NNN-distance and i, j belong to the same sublattice. Upto linear ap-

proximation( same as in Eq.(7.11) ), we have,

Ru = Tim. + 333 ' u,-,-
c 3 . .
— rnnn + go.” — rja) . (uia — uja)

3 . . . ..
Tin.” + “£10.10 ° uiar + rja ' uja) _ (rial ° uja + rja ' uia)] (730)

7‘5"", = (/8/ /3h is the mean next nearest neighbor distance.
The ensemble average of R;,- can be performed by using the projection operator

P: = (1 + 600/2.

 

. . 1
R :ng = _ J. 3113;; 193212,,- 7.31
< > (Pci1><Pcc;)(P‘-’3><P > ( )

iaj means site a in the middle with i and j connected to it as its nearest neighbors.
There are 8 combinations from 616263 representing 8 different type NNN-connection.
Switching two sublattice notation we have another 8 combinations which are sym-
metric with the above expression.

Bring Eq.(7.30) and Eq.(7.3) into above equation and perform the ensemble aver-
age. Remember that the formula Eq.(7.15) effectively eliminates the constant terms
from the averaging. The statistics is straight forward. Regrouping the results prop-

erly, we end up with following expression,

. . 8 62
R 1a) : C _ 1 _—
< >CIC263 rnnn + \/:;a ACDy( y)“ (PS2)

 

 

3 ...... 61 63
+\/%b AABx(1—$)[<P.‘1>+<PF3)]

148

1 8
= Tim, ‘1" §\/;077ACD(52 + 2y — 1)

€1+€3

3
+([§b"AAB(2x+1+ 2 ), (7.32)

 

where a, i, j denote sublattices as mentioned before. I)" is another constant similar

to a" and is evaluated from the Green’s function,

b" = (E 2;, (2,3,2, (7.33)

In the expression ,8 and a(implicit) are first neighbors of i and fl is summed over.
j is the nearest neighbor of a thus the next nearest neighbor of i.

The simple model is extremely good in describing a large class of semiconductor
material. For materials where difference in the force constants is appreciable, we
resort to the effective medium theory, which appears to be good in describing these
kinds of problems. But here we will be satisfied with the qualitative conclusions that
can be drawn from the above simple expressions, obtained for the case when force

constants are the same.

7 .1.2 Topological Rigidity Constants

The constant a", for example, has been defined in term of Green’s function in Eq.
(7.18). It can be understood however from another angle. Suppose we have a pure
diamond lattice described by the Keating model or the Kirkwood model with force
constant a and ,6. Consider an atom and apply uniform radial forces f on its four
neighbors. This Opens a cage and the neighbors gain a displacement u. If there is no
matrix effect, the displacement is just f / a. But the existence of the matrix makes it

harder to open a cage. The force needed for the same displacement is,

f = a u (7.34)

afi‘

 

 

149

Note that we use the same parameter a" as above and it is straight forward to
show it to be true. Here a" acquires its physical meaning. It represents the rigidity
of the matrix which depends on the force ratio fl/a and the topology of the lattice
construction. We call it the topological constant in opening a cage. A similar ap-
proach has been used by Thorpe and Garboczi to describe a triangular net of “bond”
disorder, and a different topological rigidity constant a" is used which corresponds to
stretching of a bond by a pair of forces along that bond. We have another constant
I)" introduced in the last section. Compared to a“, which is related to the displace-
ment of first neighbor when opening a cage, b" is related to the displacement of next
nearest neighbor when opening the same cage.

In general, there are infinite number of lattice integrals in a rigid network. Using a
simple model we have shown how the local quantities are related to each other by a
certain number of characteristic topological rigidity constants. Our approach to the
semiconductor alloys has the advantage that it produces simple analytic formulae that
make the underlying physics clear in terms of these topological rigidity parameters.
The non-additivity of the atomic radii, and the variation of the force constants a and
6, can be incorporated into our approach using perturbation theory. Non-random
solid solutions can also be handled by defining concentration dependent topological
rigidity parameters.

For semiconductor alloy material, a" is around 0.75, while it is about 0.24 for
FCC metal alloys. This reflects the large number of nearest neighbors(l2 rather
than 4). However effects like charge transfer, variation in the force constants, non-
additivity of the atomic radii, and even the failure of rigid ion potentials or embeded
atom potentials, make the situation much more complex. And a similar approach,

involving nearest neighbor central forces is much less successful for some systems like

AU1_3N2:$.[5]

 

150

7 .2 MODEL AND MATHEMATICS _

Analytic solutions are useful because it can serve as a starting point. That is why
people like to construct models which are mathematically simple and physically re—
alistic. While we define the “Length Mismatch Problem” in such a simple field Eq.
(7.2), the structure of semiconductor materials are also studied by more sophisti-
cated method like ab initio[2] calculations. Fig.7.2 compares the nearest neighbor
distance of SiGe alloy. In the plot there are result from the ab initio calculation
and experiment, and the simple application of our quaternary result from Eq.(7.27)
and Eq.(7.28). The only parameter that we fit is the topological rigidity constant
a" which is about 0.707 corresponding to a value of fl/a about 0.2 in the Kirkwood
model. As we stated before, a" is one of the characteristic of the system. Once the
parameter is known, we can express local properties, such as averaged NN-distance
and N N N-distance. Fig.7.3 is the nearest neighbor information for Ga1_zIn,,As. The
characteristic constant, a“, from this diagram, is about 0.8. And given this a"
we compare the Z—curve and distribution with computer simulation and experiment.
The experiment data agrees very well with on theoretical predictions. And from our
model, we know that b“ is about one half of a". From this we can compare the next
nearest neighbor distance of the same system with experimental results. The theory
is from Eq.(7.32). In Fig.7.4 we can see the agreement is very good. We should point
it out that while both ab initio calculations and our theory give good agreement with
the experiment, ab initio calculations show 4% upward bowing in contrast to the 7%
downward bowing observed in experiment for the Si1_xGex binary alloy. [6]

Most of our models are simple from mathematical point of View. Our goal has
been to understand the physics of certain observation with minimum mathematical
complexity. Let’s take the model for intercalated graphite multilayer system as an
example. The model is a harmonic lattice model of “bond” disorder. First, the intro-
duction of natural lengths simplifies the interaction and makes it possible to write the

energy in a compact quadratic form. The next step is the harmonic approximation.

 

151

 

 

 

 

vavfltvvfiluttwjvvyulvvvvwwfi'vvvv'vvvv'vvv
.2 2.45 r Si,_,Gex x =-- 0.55:
P- . o r ‘
.c: . . 2 he 1
'3’. - I” 1
g 2.40 r- g -’ . ‘ a
«:1 y” ‘ .
d / .
O r. ’9.) . ' b
m 2.35 ' _
l4 1 l 1 l 1 1'1 1 l l l l 1 1414.14 A i‘JAAlJJAJlALL‘lAJJA

 

 

0.0 0.2 0.4 0.6 0.8 1.0
Si 1: . Ge

Figure 7.2: The left panel shows the results of our model using the Kirkwood potential
with fl/a = 0.2 for Si1_,,Gex alloys plotted as straight lines, the simulation results
for the same model(solid symbols) and the pseudopotential results of ref. [2] (open
symbol). The lengths for pure S 2' and Ge are chosen to agree with ref. [3]. The right
panel shows our computer simulation results on a 13, 824 atom sample, with periodic
boundary conditions, for the length probability distribution of nearest neighbor S i-S i,
G's-Ge and 55-06 bonds for .7: = 0.55. The solid curves are Gaussians with centers,

widths and weights determined from the equations in the text.

152

 

 

 

 

 

 

 

P. I V T Ii V U U l U Iii l V V U V I I’ V V V W V '1‘ I 1 V V j 1 V Y 1
A r Ga,_,ln,As x = 0.20 ‘
“:1 z a i r i
‘ b '1
.c: . / in- ‘
‘63 ~ (T ‘
a » /./ t
o-I‘ P (I ' 1
'U 2.5 - -
g - / .
m (P I 7 - _ ' _ 2: ‘
-1_l l A l 1 A l 1 l 1 l l l l l l A L I 1 1 1 L 44 A A l A 4 A AJ _A A
0.0 0.2 . 0.4 0.6 0.8 1.0
GaAs x InAs

Figure 7.3: The left panel compares the experimental results of Mikkelson and
Boyce[3](open symbols) for Ga1-xInxAs with simulation results using the Kirkwood
potential with ,8/01 = 0.12 (Solid symbols) and straight lines from the theory given
in the Eq.(7.27). The right panel shows our computer simulation results on a 13, 824
atom sample, with periodic boundary condition, for the length probability distri-
bution of nearest neighbor Ga-As,In-As bonds for a: = 0.2. The solid curves are

Gaussians with centers, widths and weights determined from the equation in the text

N NN distance

153

 

. ‘. I, ..
4.2 " -—
- 4

 

 

 

4.0 _-
l' 1 L 1 1 1 J 1 L 1 d
0 0.5 1

Figure 7.4: The next nearest neighbor distance from experiment (with error bars) and
our theory prediction(solid lines). The a" is about 0.8 corresponding to 0.12 of ,B/a

in the Kirkwood model. b" is about 0.4 for the same force constant.

 

154

This step can be skipped in this particular system because the field equation is ready
to be written down as MU = ((1. Generally, within the lattice model there is disorder
in the force constant, mass and natural length. The first two are important to the
phonon structure and related electronic properties. But they are less important in de-
termining structural configurations comparing to size difference. The most important
reason to stress only the size difference, from a mathematical point of view, is that M
matrix would be the same as that for a crystalline structure. This makes it possible
to invert the M matrix to the other side of the equation and give explicit solutions
for the displacement field.We can utilize all the methods we know from Bloch theory.
The disorder, with or without correlation, can be treated statistically. From this, we
conclude that as long as disorder in force constant is negligible, we have a similar path
to an analytic solution whether there is correlation or not in the intercalating process.
The same idea applies to alloy systems when natural length concept is meaningful.
In this case, more steps are however needed in the simplification because we need the
lattice model of “site” disorder instead of “bond”. Further harmonic approximation
is necessary to reach the same field equation as for the “bond” case which we define
as the pure “length mismatch problem”. We have checked with computer simulation
that all the approximation are good up to a moderate size mismatch. (see Chapter
4)

We see that “length mismatch problem” is simple from a mathematical point of
view while it carries the feature of semiconducting materials when only structural
information is concerned. It makes use of the old concept like virtual crystal, con-
sequently the old method of Fourior Transformation or K space expansion resulting
from the translational invariance.

One however has to bear in mind that analytic solution is not everything. Simplifi-
cation many time lose subtle features of the original system in the process, especially
in dynamic systems. So care has to be taken in any simplification. It is important to

have physical insight. Our another model for covalent glasses(Chapter 2) serves as

 

155

a very good example. It is not soluble analytically, but we have seen- it works very
well. In the covalent glass systems chemical bonding is quite stable and serves as the
major bonding force to solidify the system. This picture has lead to the introduction
of rigidity percolation which states that when the number of constraints provided
by these bonds drop below the requirement to fix all the degree of freedom of the
lattice network, the system becomes floppy and there will be zero frequency modes
present. This basic idea leads to the birth of the well known concept of threshold
point (r) = 2.4. [4] It is very important to have this insight that the soft phonon ob-
served in experiment actually are another form of the zero frequency modes resulting
from the underconstrainted degrees of freedom in the rigidity percolation problem.
The transformation from zero frequency modes to low frequency modes is done by
turning on the weak forces which are neglected when doing the percolation analysis.
These weak forces are not important when rigidity percolation prevails. But they
play important roles in stablizing the system when rigidity fails to percolate. So our
model is built to generate the most important ingredient of the system, i.e., the mean
coordination number (1'). This very simple model performs very well and predictions

of the model agrees with experiments.

7.3 TWO DIMENSIONAL INSTABILITY

In Chapter 6 we studied the 2d triangular network. It is found that when the mis-
match is very large, the lattice begins to pleat. The pleating is not very physical
because there are no repulsive forces between lattice vertices in the simple model.
This transition, however, shows the unstablility generally exists in 2d system. The
unphysical phenomenon originates from the topological constraint, i.e., the neighbor
being fixed. To release this restriction, we use Lennard-Jones potential which is in
principle a long range interaction. The fast decrease of interaction strength with dis-

tance allows the cutoff in computer simulation. We first study the mismatch problem

 

156

of a binary system. It is found that the system becomes unstable around 18% mis—
match. It is possible that our simulation study is related to this phenomenon on a
very basic level. The instability of 2d binary is characterized by a constant drop of the
shear modulus and the appearance of topological defects, see Fig.7.5 and Fig.7.6. The
appearance of topological defects with the softening of shear modulus is similar to the
well known KTHN Y phase transition; but the driving parameter here is not temper-
ature but the size difference. So we call this instability “2d length-mismatch-driven
melting”. Using Molecular Dynamics simulation on binary system, it is possible to
get a complete phase diagram with two axis as temperature and length mismatch.
We hope that in the near future we can define this problem clearly and carry out a

survey of this interesting phenomenon.

 

157

 

 

 

 

 

 

 

. . . .5
. 0 00’1 0.2 033 014 0 0.6
0.3 " 0° 0 00°
0° 0 " 0.4
o
0.2 " o J
o
D“ ° ) ~o.z 3‘
0.1 1 ° °°
° °° 0.0
0.0 W: 4 . n : J. :
‘ 1.00
1.04 ' ‘ o A
”/1 o 0.75 g
A v
1...] 1.02 1 00° ' , b
v 000 l“ 0 0.50 0
0°
. ' 0.25
1.00 . 00° °°
,0" 4‘ *4 A #1 9-4—4—0———J 0,00
0 0.1 0.2 0.3 0.4 0.5

 

LENGTH MISMATCH

Figure 7.5: The x-axis of all four panels is length mismatch. right top panel shows
the percentage of defect points. Defect points have non—six N N . Right low panel is the
macroscopic size of the sample after relaxation. Left top panel is the shear modulus
of the sample and lower panel the angular correlation function. The solid line in the
left lower panel is from theory. We see that the simulation result drops below theory
prediction when topological defect is about to present. All four diagram suggest the
singularity point which characterize the instability to be about 18%.

158

>
Avflb V V V v (FA'AVANAVA 4g VEBAVAVAVAvAVAVAVAwAVAV‘y v {
uvgvgg‘figfitégge .xfiahg‘i'egggxe'sfige‘efiflfiéfia
A“ h“ A‘ ‘ V V I A.“ > a
v v k A
4» 45“v uv v‘ a ‘EYAV
‘ m“ m ‘u“ ‘h’ v VA h““uv 1
19 ‘3 f V ”m" .BF‘A'AVQEVE " v53:

'4‘:
\uuv vr‘ v ’
A 4§L " AVA‘NAV‘V‘VE‘ ‘1‘ u

b

A
figfiAV‘AVA‘

 

A
anyways!“ 4&7;
VAVAVAVAuVAuvAE ggsvauuv‘uuv; AV
<uuuuuvmv4
4
Aunuufl'flgi

 

V
AVAVAVE 451""4‘
VAVAVAVA“vegexeggflahvgigzfi‘vxggvgvéxv‘v‘VAVAWAuV , 4uuuuuu
exev‘uv‘u v v“ v‘u'AVAvszh'A'A'A'AVA‘ an!” 45
..ssfissvseuwgfiu91.23% assessemseezezaae
Ehfl'fie'flvwhuifififixfifiauv Av‘v‘auufimuuuu‘iflfifiw‘
séiifivs'sfissseu slams smears. 1%.... vs: .19..
"A'A'Afihuflh 4 " '5' «V V ufl w: <hv u‘ 0
h€€¢ AV“ "6? 5.4% vgfimgfimlg'

AV» “VA 0' 1
'7 WVVVVVV gpfi
“WW“MMV‘ ‘ ‘ ‘bQ nfifi‘h namtumaug

 

 

aflggflk gm A “kl“
1 fl v10. vmfitmgmv
Ivfifivgvgfiéggpana‘uukuiuFaun4AAsunny“

Figure 7.6: The 22% size mismatch diagram. The nearest neighbor connection follow
the Voronoi algorithm. The Topological defects are points having more or less than

six neighbors. Quadrupoles and dipoles(dislocation) are seen frequently. As mismatch

increases, the defects begin to connect and percolate.

Bibliography

[1] Y. Cai and M. F. Thorpe, to be sent to Phys. Rev. B.

[2] S. De Gironcoli, P. Giannozzi and S. Baroni, Phys. Rev. Lett. 66 2116(1991)

 

[3] J. B. Mikkelson and J. B. Boyce, Phys. Rev. Lett. 49, 1412(1982)
[4] H. He and M. F. Thorpe, Phys. Rev.Lett. 54 355(1983)
[5] N. Mousseau and M. F. Thorpe, submitted to PRB.

[6] J. P. Dismukes, L. Ekstrom and R. J. Paff, J. of Phys. Chem. 68,3021(1964).

159

Appendix A

The Equation of Motion

The equation of motion techniqueis a variant of conventional molecular dynamicsthat
is useful in disordered harmonic systems. The initial conditions at t = 0 are deter-
mined by the particular response function it is desired to calculate, and the equations
are integrated forward in time with appropriate weighting factors. The equations of
motion for any harmonic potential V are given by

Miuia(t) = — Z Wigwam, (A-l)

7.3

where Ugo, is the displacement of the ith atom in the a direction. It is convenient to
think of the age, as classical displacements. The motion in this system can be written

as a superposition of normal mode vibrations. Thus,
11,0, 2 (l/M,‘)1/2 Z qn§n(i, a) cos(wnt + 5n) (A.2)

where 1.2,, is the frequency of the nth normal mode and 6,, is the phase shift that
is determined by the initial conditions at t = 0. The quantity qn is the amplitude
of the nth normal mode and €n(i,a) is the polarization eigenvector which gives the
magnitude of the nth eigenvector on particle i in the a direction. These obey the

usual orthogonality relations

2 = {710.7 a)€m(i1 0) = 67m; (A.3)

160

161

and ..

Z = 6.03 am, (3) = 6.16... (AA)

71

with the eigenvectors chosen to be real because the dynamical matrix is real and
symmetric.

In previous work, the initial conditions have corresponded to “plucking” in which
the initial velocities are set equal to zero, but the displacements are finite. Here
we use the same notation(see Chapter 2 and references there) to descibe the “kick
start”.

The initial conditions for a “kick start”, where an impulse produces an initial

velocity but no displacement, are given by putting 6,, = —7r / 2 in Eq. (A.2):

u;a(t) = (1/M,‘)1/2 Z qn{n(i, a) sin(wnt), (A.5)

so that the velocities are given by

12,0,(25) = (l/M,)1/22wnqn{n(i,a)cos(wnt), (A.6)

and the initial velocities are

1.140(0) = (1/Mi)1/2:wnqn€n(iva)a (A7)

which can be inverted using the orthogonality relation (A3) to give

wnqn = Z Mil/zaj.fi(0){n(j1 :3) (A'S)
it?
We now compute
T
G(w) = if dt sin(wt)e_3(‘/T)2Zu,a(t). (A9)
0 {a

We see that

ZAiauia(t) = ZZZ(Mj/M{)1/2€n(i, a)§n(j,fl)A,-atljg(0) sin(wnt)/wn, (A.10)

5.0: 1.3 ’1

 

162
and hence v

G(w) = Z 2; 2(Mj/Mamaa. am, fl)A.-aizm(0)[A(w — w.) — New + «val/w...
' J' (A.11)
and choosing

...-am) = A... = \/2cos 9..., (A.12)

where 0,0, is a random angle uniformly distributed in the interval 0 < 0,0, < 27r. We

have

(GM)

2 2511(1) a)[A(w - wn) - A(w + wn)]/wn,

:2 [A(w—wn) — A(w+wn)]/w, (A.13)

which is the density of states divided by the frequency p(w)/w. Note that this form

leads to an odd function of the frequency because of the sine transform.
2 T 2
p(w)/w = ;/ dt sin(wt)e_3(t/T) Z Ammo“). (A.14)
0 1,01

This form artificially forces the density of states to zero as the frequency goes to
zero. By happy chance p(w) in Eq. (A.14) goes as (.112 at low frequencies, which is the
correct form in three dimensions. Thus the form imposed by the numerical method
leads to the correct form and will aid convergence somewhat. If we integrate both
sides of Eq. (A.14) over frequency and integrate the right hand side by parts, we find

the first inverse moment
Wmas‘ 2 T —3(t/T)2
(may... = (1 /w) = ; / dte 2A,au,a(t)/t. (A.15)
0 1,01
Because of the “kick start”, the displacements u,a(t) . t at small times so that the
integral (A.15) converges at small times. It also converges at large times because the

Iu;a(t)| are bounded and so the 1/ t factor ensures that the integrand decreases.

 

 

 

 

Appendix B

Distribution Function of

Multilayer System

In this appendix, we show how the first few moments can be extracted from the
distribution function for the random case when kA = k3. Recalling Eqs. (3.68) and

(3.71), we can write down a generating function for the moments,

Z czq" = e’iu’hhAii (cos qAfi — i(1 — 2:12) sin qAfi) (B.1)
” 579i

Taking the logarithms and expanding up to fourth order, the left hand side gives

1110 + 6'17 + 63112 + 0&1" + 629‘ + ...)
= q[C'{l
2 It 1 11:2
+9 [62 - 5C1 l
3 In In :1: 1 1113
+9 [c3—c1c2+§cl]

1 1
+q41c; — 5c? — czc; + 6‘ch — Zea“ . (13.2)

163

 

 

164

From the right hand side of Eq. (B.1), we have

ln[ 6 _iu'hnAii H(cos qAfi — i(1 — 2x) sin qAU)]
#5

= q2[—2x(l — as) 2A3]

j¢i

—q3[i§:c(l — 2:)(1 — 22:) 2 A3,]
m

2
+q4[-§at(l — x)(1 — 62: + 61:?) 2A3] . (B.3)
#1
We confirm that c; = 0 and also find that

 

c; = —2:1:(l ‘3)2Ai2j
513i
It 4 3
c3: —2§$(l—$)(1—2$):AU
j¢i
c; = §x(l — x)(l — 6:2: + 62:2) 2: Aii + é-cgz. (BA)
.iaéi
Combining with Eqs. (3.24), (3.25) and (3.72) we have
2k;
M2 ‘ (4)2
= 4:1:(1 — xx: A125 — A3,)
5
= 23(1— :13) [W1(k) — W(k)2] , (13.5)

which agrees with the result obtained previously in Eq. (3.23) for the fluctuation in
random case.

We also obtain results for the higher moments,

310;
M =
3 (-i)3

 

 

165

= 8x(1;x)(2z—1)2A3, (B.6)
#i
and
M4 — 3M,2 = :z:(1— x)[1- 6x(1— 1)] 2A3,- (B.7)

#1

 

 

Appendix C

Symmetric Distribution Case

It is apparent from Appendix B that all the odd moments have the general form
(for n 2 1):

M2,,“ = Z f,,,(:z:) 2A3“ . (0.1)
(:1 j¢i
When 13,; = k3 = 1:, from Eqs. (3.28) and (3.61) and denoting kT/k as a, for the 2d

square net only, we can rewrite Aid as

-1. [1-..s(q..,)]aq-ru
2N :q: [1 — cos(q - V)] + all — cos(q - 6)] .

 

Aw) = (0.2)

Noticing that the 93(6) and 2(V) directions of q are all summed over, we can switch 6

to z/ and vice versa without changing the integration. Then we have

 

,q.,..
Add.) : -21? g [1 —— cc:[(1q-T :fisg-qdfln—e cos(q. 6)] (0.3)
From above two equations and Eqs. (3.24) and (3.61) we have
W(a) + mi) = 1 ((3.4)
Now making use of the fact that
6(q) = i: a” (0.5)
1

166

 

 

 

 

167

and using Eq. (C.2) , we obtain the following expression for n > 2: a.

 

n _ 1 1 [1 — cos(q’ - V)] . ..
Z: A“ (—2)nN"-1qqnmq(n) [1 — cos(q’ - 11)] + a[l — cos(q’ - 6)]
[1 - COS(Q"° 11)] ...
[1 — cos(q” - V)] + a[1 — cos(q” - 6)]

[1 - cos(Q‘“) - V)l
[1 — cos(qW - V)] + a[1 — cos(qW . 6)]

 

5(q/+ q" + ‘ - - + 01“”) - (C-G)

We make the replacement:

 

 

[1 — cos(q") 1 11)] = 1 _ a[1 — cos(q“) - 6)]
[1 — cos(q“) . 11)) + an — cos(q“) . 6)) [1 — cos(q“) - v)) + an — cos(q“ '(63317)

in Eq. (Q6) and use the relation between Aii and W, to get the following result.

2J3A5(a) = (Avg/«1min (gm — mg»: — (arm-13“) . (0.8)

From Eq. ((3.4) and using Aii = — W/2, we have,

>J3A5<a) = <—1)” ;A5<§) + «gnu/(a): — <-—1)“W(;1,-):1

— avg/15%) + [A5(a) — (Ar/15%)]. (0.9)

When I: = kT, which means a = 1, we have

;A5(1)—(— IVE/1:51) =A5< (1)— (-1)"A:5(1) , (0.10)

which leads to the conclusion that when n is odd, Zj A53 = A3} or 23,5 A3 is identically
zero. So the right hand side of Eq. (C.1) is zero when k = k7, and all the odd moments

of the distribution vanish for the 2d square net. For other lattices, this symmetry is

not present.

Appendix D

Average NN-distance in

Semiconductor Alloys

In this appendix we will derive (dAc) as an example. The idea is the same as that in
the multilayer system,(see Chapter 3)but the method is modified because it is site
problem here instead of bond problem.

We can neglect notation C in A— C and B —C' bonds. Because C forms a sublattice
which bridges the A and B, every A and B should be connected to four C while they
never connect to each other. We thus can make C as dummy variable in summations.
In many cases we only write those site variables that are A or B.

Define projecting operator Pf which is 1 if site i is A and 0 if i is B. The other
projector for BC bond is P}? and is complement of A projector. It is obvious that
(PAC) = 1 — a: and (Pi?) =.:c.

Minimization of energy in Eq.(4.10) with respect to the displacement vector U

gives the solution.

:- qu (D.1)

168

169

Bring Eq.(D.l) into Eq.(4.2) and perform the ensemble average for AC type bonds,

we have,

Um) = 1.1. ii—(§)<P1"¢mc1)G.-a,ma (11.2)

“0” after the summation is the ensemble average. “0” under the summation is
understood to exclude double counting.

The Green’s function GiC,mC’ is a short hand for < iCIGImC’ > Note that summa-

tion over sub-indices of GiC,mC’ is zero because of symmetry. So any constant term

from the ensemble average is effectively zero in the average length expression.

Bring Eq.(4.12) into Eq.(D.2) and simplify it,

 

(ll/1C) =he+ 1 a

Z<Pfh$nm>cwmg (D3)
7710'

The natural length can be expressed in the projectors as:

112,0, = 119,193)” + 1135,1130 (DA)

We need to work out ensemble averages (PACPX‘C') and (PACPEC').

Note that C and 0’ all belong to C type.

(PfPX‘C') = (1 — :c)2(l — 6m) + (1 — n5...

2 (1 — 2:)2 + x(1— x)6.-m (D5)

The first term is effectively zero because it does not contribute to the final result
of (hAc) after the summation. We use a notation SE to replace 2 and discard these

terms on the right side of the equations,
<P10PEC’)3£—f—x<1 — 1):... (11.6)
The same for (PACPE‘C'),

<P1CP5”C’)3{1 — a(1 — as. (11.7)

170
Bring Eq.(D.6) and (D.7) into Eq.(D.3), _.

((1.40) = he + $013 — hos) Z 0511116]:me
mC’

5 he + 3:01:34 — h%)a" (D.8)
Note that h6 = (hij) = (1 — :c)ha +xhoB, Eq.(D.8) giVes (dAC') in Eq.(4.15).And the

Waston integral is related to Green’s function in following way,

a" = Z “Gated" (D9)
CI
When apply this method to quaternary system or pure bond mismatch system, a"

shows up and has the expression:
a" = 06};an (D.10)

Same method applies for energy calculation and the introduction of projectors

simplify the mathematics.

,11111111[[1111[1111111))I"