PLACE IN RETURN BOX to move thIo checkout from your record.
TO AVOID FINES return on or balm duo due.

 

DATE DUE DATE DUE DATE DUE

 

W

 

 

 

.1

l

'n:: 7 ’. fj-
» 2 6 7 ..
m msozaoz I

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU Is An Affirmdivo Action/Equal Opponunity IMMIOI'I

 

 

VIBRATIONAL DYNAMICS OF DISORDERED DIPOLAR

AND IONIC SOLIDS AND LAYERED SILICATES

BY

Joel Moritz Gales

A DISSERTATION

Submitted to

Michigan State University

in partial fulfillment of the requireuents

for the degree of

DOCTOR OF PHILOSOPHY

Departnent of Physics and Astronouy

1989

ABSTRACT

VIBRATIONAL DYNANICS OF DISORDERED DIPOLAR

AND IONIC SOLIDS AND LAYERED SILICATES

By

Joel Moritz Gales

In this thesis we study the vibrational dynamics of two
classes of materials. In the first part we look at
disordered solids in which long-range Coulomb forces are
present. We first present the basic features of these
forces in crystalline materials, in particular, the LO-TO
splitting, and the techniques used with special emphasis on
the Ewald method. The computational technique, the
Equation-of-Hotion (EON) method is described which we apply
to a study of depleted simple cubic arrays of point dipoles.
We show that a trucated approach is numerically equivalent
to the exact Ewald method but computationally simpler and
faster. The main results of this section are given in a

study of the LO-TO splitting in (III-V) IV mixed crystals

semiconductors. Using both Raman and infrared responses, we
show that this splitting vanishes when the III-V ionic
material no longer percolates. We demonstrate this for both
a one- and two-mode behavior material and compare our
results with recent experimental data on these systels.

In the second part, we present a theoretical survey of
the lattice dynamics of layered silicates or clay. We
present the zone center frequencies, the phonon dispersion
curves and the elastic behavior. We derive expression for
the low-frequency zone center modes that should be useful in
the extraction of force constant values from experimental
data. We also study the torsional mode dynamics of ternary
intercalated alloys using the EON method and develop a
simple theory for the frequency as a function of varying
concentration that incorporates the effects of layer

rigidity.

DEDICATION

To my parents,

who listened to it all.

11

ACKNOWLEDGMENTS

I would like to thank my advisors, Prof. H.F. Thorpe
and Prof. S.D. Hahanti for their guidance and. support and
for our mutual perseverance through sometimes difficult
periods.

I am grateful to Prof. S.A. Solin, Prof. H. Repko, and
Prof. H. Abolins for their service on my guidance coumittee.

For his careful presentation, both in person and in
print, of the somewhat subtle features of the Coulomb
interaction, I would like to thank Dr. 3.". deLeeuw and for
her assistance in tasks both large and small throughout my
time here, I thank Hs. Janet King.

Finally, for their friendship, humor, and refuge, I

wish to thank Roy Day, Rich Kennan, and Dr. G. Pollack.

iii

TABLE OF CONTENTS

LIST OF FIGURESOOOOOOOOOOOOOOOO ...... OOOOOOOOOOOOOOOOOIOV1
LIST OF TABLESOIOOOOIOOODOOOIOIDIIOIOO‘OIOOOICOOOOIIIIIOix
MI.

I. IntrOductionOOOCOOOOOOOOOOOOIOOOOOO‘OO0.0.0.00002
II. Equation of Motion for Dipolar Systems.........6

III. Dynamics of Periodic Dipolar Arrays
A. Ewald Method for Periodic Lattices...... ....... 9
E. Zone Center Modes - Depolarizing Factors......l7
IV. Dynamics of Disordered Dipolar Systems
A. IntrOductionllIOOOOIOOOOIOOOOI0.00.00.00.0000039
B. Equation of Motion / Linear Response Theory...43
C. Implementation fo the Coulomb Interaction.....52
DO Effects Of DisorderOOOOO...0.00.00.00.0000000062
V. Coulomb Effects in (III-V)bm IvgxiMixed
Crystals
A. IntrOductionOOOOOOOOOO0.0.00.0...0.0.0.000000073
8. Implementation of Coulomb Interaction for
Ionic syste.80000.0.0.0000...0.0.0.0600000000077
C. Equation of Motion / Simulation Techniques....80
D. Construction of Mixed Crystal - Percolation
' Of DinerBOOOOO00.00....000......0.0.00.000000086
E. One-Mode Behavior (GaAs) Ge..................92
F. Two-Mode Behavior (InSb) Ge.. ...... ..........98
VI.‘ Summary and Conclusions......................107

iv

APPENDICES

A. Fourier Sums for Three-Dimensional Periodic
Functional.OOOOOOOOOOOOOOOO00.0.0000000000000109

B. Relation between the L0 and TO Responses
for Mixed Crystals...........................112

£11111

I. Introduction.......................... ...... .116
II. Crystal Structure............................120
III. Phenomenological Force Constant Models.......126
IV. Lattice Dynamics of the Periodic Model System

A. Zone Center Modes............................l32

B. Phonon Dispersion Curves.....................l4l

C. Elastic Behavior.............................lSl
V. Torsional Mode Dynamics in Ternary

Intercalant Alloys

Introduction.................................160
Description of Model.........................162
Torsional Mode Frequency / Layer Rigidity....165
Simulation Results...........................l76

cocoa»

REFERENCES.............................................l77

2(a)
2(b)

2(c)

5(a)
5(b)

5(c)

10

11

LIST OF FIGURES

2111.1.

Model System for Dynamical Calculations ........ 25
of Ionic Systems

Energy Dispersion for SC Dipole Lattice........30
Energy Dispersion for 300 Dipole Lattice.......31
Energy Dispersion for FCC Dipole Lattice.......32
Brillouin Zones for Cubic Lattices..... ........ 33
Dispersion for Simple Cubic (100) using........34
Truncation and Ewald Methods
Density of States for SC Dipole Lattice........36
Density of States for BCC Dipole Lattice.... ..36
Density of States for FCC Dipole Lattice.......37
'003 for Disordered Simple Cubic (p=0.9)........59

Im 61(3) for Disordered Simple Cubic (p=0.9)...59

Im

II

Im

Im

:.(E) for Disordered Simple Cubic (p=0.9)...60

41(E) for Disordered Simple Cubic (p=0.9)...60

ct(E,E)

for Disordered Simple Cubic.........61

(p=0.9) f=<2n/L) k

for Disordered Simple Cubic.........61

(p=0.9) I=<2n/2L) i

vi

12(a)
12(b)
12(c)

13

14

15

16

17

18

19

20

21
22
23

24(a)

24(b)

25
26(a)

26(b)

DOS for Disordered SC Lattice.............

DOS for Disordered BCC Lattice............

DOS for Disordered FCC Lattice............

DOS for Disordered Cubic Lattices (p=0.05)

and Random Dipolar Gas

DOS for Dipolar Gas with Lorentzian Tails.

from Dipole Pairs

DOS for Quenched Dipolar Liquid (Pick)...

Conventional unit cell of zincblende.....

structure

Percolation Statistics for Dimers on Diamond.

(III-V) IV "Mixed Crystal" on Hexagonal.
Lattice

Raman Response for (GaAs) Ge............

Peak Positions of Raman Response for....
(GaAs) Ge

Infrared Response for (GaAs) Ge.........
Raman Response for (InSb) Ge............
Raman Measurement for (GaSb) Ge (Expt.)

Peak Positions of Raman Response........
(GaAs) Ge (Expt.)

Peak Positions of Raman Response........
(InSb) Ge (Theory)

Infrared Response for (InSb) Ge.........
Longitudinal Response EOM vs KR (x=0.2)

Longitudinal Response EOM vs KK (x=0.8)

vii

..63

..67

..67

.102

.103

.104

.105

.114

.114

3(a)

3(b)

6(a)

6(b)

8(a)

8(b)

10(a)
10(b)

11(a)

11(b)

12

13

2.112111
Structure of Model Layered Silicate........

Low Frequency Zone Center Modes............

w(f) for I (0,0,k).......................

“(3) for E

(0’k,0)oeoooeoooaaoomoooaoaaao

Structure for y-axis Dispersion Model......

Density of States of Model Layered Silicate.

Sound Velocity as calculated from Model....

Sound Velocity for Muscovite as measured...
by Brillioun scattering

w(i) for f = (0,k,0) for 31205 Layer........

Raman Measurements of Torsional Mode.......
Torsional Mode Frequency vs x (Expt.).....

Two-Dimensional Kagomé-Intercalated Lattice

Infinitely Rigid Clay Layer.................

Floppy Clay Layer..........................

Torsional Mode Frequency vs x (Theory)....

(Sub-linear Behavior)

Torsional Mode Frequency vs x
(Super-linear Behavior)

(Theory)....

Raman Response for Torsional Mode..........
calculated with EOM

Density of States for Kagomé-IC calculated.
via Direct Diagonalization and EOM

viii

..121

..133

..144

O .145

..147

..150

..154

..154

..159

..161

O .161

..164

..169

..169

..171

..172

..174

..175

LIST OF TABLES

Table 1 Structural parameters of the model clay.......123

Table 2 Phenomological Force Constant Values..........131

ix

WWQW
mmmm

mem

In their paper "Coulomb Effects in Disordered Solids",
Thorpe and de Leeuw (Th 86) developed a general technique
for computing the vibrational response functions of
disordered solids (such as glasses and mixed-crystal alloys)
in the presence of long-range Coulomb forces. They derived
moment sum rules for the transverse and longitudinal
dielectric responses, a generalized LST (Ly 41) relation for
the static and high-frequency dielectric constant, and most
importantly, a general relation between the transverse and
longitudinal responses allowing one to be calculated once
the other is known. They illustrated these results by
computing the response functions for a model AX: glass
roughly corresponding to vitreous silica. These simulations
showed, unambiguously, the existence of longitudinal
optic - transverse optic (LO-TO) splitting in disordered
materials for which there had been experimental evidence for
nearly a decade (Ga 76, Ga 83). It had been hoped initially
that these dielectric responses could be calculated
theoretically, albeit approximately, as opposed to being
merely computed. However, the long-range nature of the
interactions coupled with the disordered nature of the
materials makes this a difficult problem and we have been

unable to make much progress in this regard. We must

content ourselves, therefore, with results essentially
empirical and descriptive in nature.

We'll begin with a look at a prototypical Coulombic
system which displays the characteristic longitudinal
optic/transverse optic splitting of ionic systems, and
provides a good introduction to the techniques used in
studying these materials. We study crystalline cubic
lattices of Einstein oscillators coupled via long-range
dipole interactions. For these purely "optic mode” systems,
we introduce the Ewald method, which transforms the slowly
convergent dipolar interaction term (V(r) + l/r') into an
equivalent, rapidly convergent, expression. In addition, it
demonstrates explicitly the characteristic splitting near
the zone center (i=0). We'll next show how the macroscopic
shape of the sample can be incorporated through the
depolarizing factors and reveal the connection between the
longitudinal and transverse modes at small wavevector and
the response of the system, of a given shape, at the zone
center. Dispersion curves, w(E) vs. E, will be presented
for the three cubic lattices along with the vibrational
density of states.

With this foundation established we can consider the
effects of disorder. In this part, numerical simulations
play a major role. We use the equation-of-motion (EOM)

method for virtually all computations on the disordered

materials. This technique computes the response of a system

to various external probes and will be used primarily to
calculate the density of states (DOS) and dielectric
(infrared) response (Th 86). We will see how the Coulomb
forces can be implemented once again using the Ewald method
and introduce a simpler, more physically intuitive approach
that gives essentially the same results as the Ewald method
but is faster and easier to implement (Ga 89a). We use this
technique to study the effects of disorder on the DOS of
point dipole systems, by considering depleted cubic
lattices.

This section culminates with a look at a more realistic
ionic system, that of a model (III-V) IV mixed crystal
semiconductor. For x=0, we have a crystalline. partially
ionic, zincblende structure with well-defined zone center
splitting of the optic modes, for x=l, a purely covalent
material. We will compute the dielectric (IR) and Raman
response for two systems, one exhibiting one-mode behavior
based on (GaAs)‘__xGe2x and the other with two-mode

corresponding roughly to (GaSb)t«Ge and show that the

2x’
LO-TO splitting vanishes not at x=1 when the effective
charge disappears but rather when the (III-V) material no

longer percolates. Thus we see that a local,

nearest-neighbor property can have dramatic consequences for

an effect usually associated with, and indeed dependent on,

the presence of long-range forces.

mwummmm

Consider first an isolated pair of particles with equal
and opposite charges, q. They are coupled together by a
Hooke spring with force constant k. We neglect explicitly
the Coulomb interaction between the two charges, however the
"short-ranged" part of this force can be considered to be
included in the harmonic potential. Fixing the center of
mass in space to prevent translational motion, the

Lagrangian of the system is given by:
) , (1)

where 3+ and 3_ are the displacements from equilibrium of

the two charges. From the center of mass condition,
m+ap+ m_3;=0, and the expression for the dipole moment,
a = q (3 - 3_), we can re-write eqn. (1) in terms of the

f

more convenient dynamical variable 3:

_ L m 42 _ ‘1 k +2
‘- - 2 Ir)“. 2 [flux ‘2’

where m is the reduced mass given by, m+m;/(m++ m_). We now
add the dipolar interaction between this dipole, the 1th,

and all the others in the system, which we consider to be

identical. This is given by (Ja 74):

 

 

+ + + +
P.‘ u. 3 (pa 3..) (us 3..)
vaa(1)= E [ " ’ ‘ "’ ’ "] (3)

8 - 5
. . R..
Jflt. t.)

where Ru is the displacement vector between the 1th and jth
dipole. Adding this term to the Lagrangian in eqn. (2) and

using the Euler-Lagrange equation (Ma 70):

d 0 L 0 L
__ _. --—=0, (4)

 

dt 0 z 0 3
we get the equation of motion:
.. qz z! 3 (:13 3..) N”
9 z 9 J J ‘J U
u... p.._Z _. =0, (5)
t. o I. I 9 8
.. R.. R..
33“. 1.3 I.)
where ab: = k/m. Because of the linear nature of this
equation, we have harmonic solutions of the form:
fist) = Bieiwt. Substituting this into eqn. (5) we get:
(wf- mi) 3 = Q,3 where D, is a matrix containing the
.’
structural information and p = (p‘x,p‘y,uu,

,HN*,HNV,HN8). Since Q,is independent of mo, so are its

eigenvalues and thus changing w: merely shifts the origin of

the spectrum when expressed as a function of 03.

Furthermore, by introducing the plasma frequency, defined

by:
N z N
q, a d s
z _ 4 n i = 4n9 L_ ‘L 4mg
we- 02,-, 0 2[m++m_] 8 mp (6)
i. i

where Np(Nd) is the number of particles (dipoles) in the
volume 0 and p is the number density of dipoles, we can

write eqn. (5) as:

 

 

z 9 9
fi.+w’fi.+ " ——’—- ’ " " =0 (7)
t. O t. 4719 3 R5
13“. ij ij .

The term in brackets in eqns. (5) and (7) is (minus) the

electric field at the 1th dipole due to the Jth dipole
located a distance Ru from 1 along the unit vector

rn= RH/R . We can therefore write eqn (7) as:
t: L: t1

2
9 2 9 up
4' (o -

u. 0 Hi

 

E(R.) (7a)
1. I.

4np

where 3(3) is the total electric field at the 1th dipole due

to all the other dipoles.

mmuwmm

memmm

Periodic arrays of point dipoles are an important
special case of the harmonically coupled systems we have
discussed above. Rather than dealing with an interaction
matrix D of dimension 3Nd, we need only consider a matrix of
dimension 3s, where s is the number of dipoles per unit
cell, that is, the number of elements in the smallest
periodically repeated unit making up the system. In this
section we will consider only Bravias lattices, where s=1,
and will develop the Ewald method appropiate to this case.
We return to eqn. (3), assuming there is a dipole at every
lattice point defined by the set of direct lattice
vectors, R. The interaction energy of a given point dipole,

Do, with this periodic lattice of point dipoles is given by:

F___o-R 3d) 3(30- fluid!)- i!)
“(140) :32 - ] (a)

 

where the sum is over all direct lattice vectors and

R a IRI. Because of the periodicity we consider solutions

iE'R

of the form: 3(R)=3°e , therefore:

10

 

(9)

As it stands, the sum converges too slowly to be of
practical use. The Ewald technique re-expresses the single
sum as two rapidly-convergent sums, one in real-space, the
other in reciprocal space.

We start the derivation with the integral

representations:

1m
1 2 ix: “-83%
.3 = — I dt t e (10a)
R f; o
1 4 an axe -I2t.
-3 = -—— J'dt t e (10b)
R an”? 0
which follow from the substitution u = Rt‘flz. Substituting

these relations into eqn. (9), we get the two terms:

*3

as
«It

a.»

11

 

 

2“: * R m
= ——‘—’- 2 e‘k' In: t‘” e “ (11a)
a Fee 0
I
{II-R 3(“0' it)
I a —2 e
s R5
Rec
4 9 3 m 2
= - — 2 ed" (30- i!)2 In r’” e"" ‘ (11b)
{5 iso °
0
4 . 9 . _ a
= — (30° (’82 2 e‘(“?) if Int t’” e “ (11c)
{E Karo ° ?-o
We split the integrations into two parts, from O to n2 and

n2 to at where n is a freely adjustable parameter. We will

deal with the four resulting terms one at a time beginning

with 1;” where:

 

s
.. _2
o 2 e‘kij‘dt ti/se mi (12)
o

and

12

sz 9 R 1m e
1;” 5 -——° 2 ea. In t": e" ‘ (13)
' a 3:0 n2 .

With the substitution u E R t‘AL the integral in (13) can

be written as:

Q

 

2 z -u2 -s am.)2 1’1?
—; Idu u e = R nR e + -§-erfc(nR) (14)
R 0‘
where we have integrated by parts using a = u and
2
dt = u e“u du. Thus we have:
2:42 eff-i! z
1‘” -..- —-2 2 mt e"”" + E erfc(nR) (15)
° 4.? n' 2
3:10

Because both terms in the brackets go exponentially to zero

as R-90, this term is rapidly convergent. We now consider

(1)

I.

and begin by interchanging the order of the sum and

integral. This gives:

a
3 n
2p ,9. _ 2
1:95 o I“: ti/zz .ikRent
n o

l3

s 2
2 n 2 n
2p ,9. . _ 2 2p
= o I“: ti/z eik R e mi __ o Idt tax: (16)
f n o f n o

where we have added the 380 term to the sum, and compensated

with the second term. The sum in the first term is
9 9

unchanged when k -9 k + 3 where C is a reciprocal lattice

vector and can therefore be expressed as a three dimensional

Fourier sum in reciprocal space. This gives (see

Appendix A):

8/2

9

gait-R e-Izt ___ :7 [g]

ge'm’z'zl“ (17)

where the 3's are the reciprocal lattice vectors and u is
the volume per dipole in real space and thus equal to the

inverse of the number density. From eqn. (16) then:

 

2
s n 2
Zap _ 9_ 2 4p
1;» z o o I“; e Ik GI /at _ o T). (18)
o 34; .

As we shall see shortly, the first term cancels with an

identical term from T;” which we now evaluate.

14

2
n 2
T“) = :1 dt t3/2 ”9‘ _ V )2 et(:+?)-§ e-m t.
5 fi 0 0 {
i=0
s
n 9 9 2 ’0 1 -lk’-Cl2/at
n2
_ 2 1 -|E-C|'/at
21.10 ; Iodt t— e (19)

where we have applied eqn. (17) to the real space sum.

Thus:

- 2
Ico+ T“) = 42 o e (k/an) +

9 9 2
4—32 It‘é‘k‘a’“ e_|f_a|2/‘n2
o |g_a|2

-54" n g (20)

:1 lot»
m

duo

where we have separated the irregular 3-0 term from the sum.
As 3 goes to zero, we see that the value of this term
depends on how this limit is taken, as it is a function of

the angle between 30 and E. For the longitudinal mode

15

(30 | E), this term gives 4nu2/o, while for the transverse
mode (30 1.3)! it gives zero. Thus we have explicitly
demonstrated the zone center LO-TO splitting. The ambiguity
in the limit is a consequence of the long-range nature of
these forces which requires that the boundary shape be
specified to completely define the problem. This is not
surprising when we recall that problems in electrostatics
are essentially boundary-value problems. We shall see
shortly how the sample shape can be incorporated in the

Ewald sum through the use of convergence factors.

(2)

We return now to the evaluation of T5 . From
eqn. (11b):
(2) 4 9 2 tk’R O s/2 -mzt.
T5 = - —-- 2040' R) e Idt t e . (21)
f n a’o ”2

The integral can be evaluated once again through the

substitution u

R t}/2 and integrating by parts using
2

s I u' and dt 8 u efi‘du. This gives finally:

 

2 9 2

1..., 1...: _2 2 a” { [52 _ 3W. ‘5’]
D B I 8
*fi'i—o R R

16

3
x [72R e-mm + £3 erfc(nR)] } -

(i5 - fi)’ (22)

Equations (20) and (22) taken together give the Ewald
expression for the dipolar interaction term. It is valid
for non-zero E but as stated earlier, remains undefined for
E = 0. The adjustable parameter n controls the convergence
of the real- and reciprocal-space terms. For large n, the
convergence of the first term is increased while that of the
latter is decreased. For small n, the reverse is true. Its
numerical value is chosen to optimize the total convergence.
For the calculations on the cubic arrays, which we present
in the next section, we use 7; I 2.0 a“, where a is the
linear dimension of the conventional unit cell. Before
presenting results for these three lattices, the simple
(SC), body-centered (BCC), and face-centered cubic (FCC), we

first discuss the procedure for handling the zone center

(E = 0) modes.

l7

al.22ns.antst.hsdea.:.anelatizins.tastsrs.
We begin the analysis of the zone center modes with the

modified interaction:

 

2 9 s
M 3m - R)
.9 O 0
V (p 'E=O) = lim -—- — x
dd 0’ s90 :2 [R' R. ]
a: a2 3:
exp -s —;’+ -%’+ —; (23)

a b c

where the second term is an ellipsoidal convergence factor
that imposes a "shape" onto the sum through a directionally
dependent, exponentially decaying, weighting factor. The
rigorous justification for the mathematical form of this
factor is given in de Leeuw et a1. (dL 80) where it is shown
that the infinite sum in eqn. (23) is equivalent to a finite
sum over the supercells with a surface defined by s[(x/a)z +
(y/b)z + (2/c)2] = 1. We obtain the desired expression for
infinite sample size in the limit of the scale factor s,
going to zero. Proceeding as before, we use the integral
representations for inverse distances, eqns. (10), and split
the integrals into two parts. We need only concentrate on
the evaluation of the integral with limits from 0 to n2 as
the other term simply lead to eqn. (22) in the limit I 9 0,

once the limit s 9 0 is taken. Evaluation of T;” for this

18

case, leads to an expression similar to eqn. (16) but with

the sum over R replaced by:
_ s a _ s a _ s a
gexp[ (t + is)Rx]exp[ (t 4- Ss)Ry]exp[ (t + EHRZ] (24)

This sum can be re-expressed in terms of a three-dimensional
Fourier sum. From Appendix A, part II, we get for the a = 0

term:

3/2
T!

 

(25)

(t + £2)" (t + fie)" (t + gn‘”

We note that for finite C, in the limit s 9 0, we recover
eqn. (17) with the C = 0 term excluded and with E set to

zero. Equation (25) gives:

 

tn _
Ia (8s0) -
1
2n ' n ta/s
11' 7E0 I dt 2 4/2 s 1/2 a 1/2 (26)
s90 o (t + :s) (t + 52) (t + 3')

Note that if the limit is taken, this leads to a logrithmic
divergence. This cancels with an identical divergence of

opposite sign from the term T;”(C=O) which is given by:

19

(u _ _
1's (3-0) -
2
fl a/s
4n t
lim —— dt x

 

.90 0° (t+:‘)s/s(t+%.)u2(t*%.)s/2

 

2
9 . 2 -( /4(t+s)
(yo Vz) e =0 (27)
The second term of the integrand gives:
#2 u' #2
_—_xs 4» __L_s + ___2‘ (28)
(t+;¢) (ti-~52) (ti-Es)
Adding eqns. (26) and (27) gives finally:
3
9 9 4n 2
vuwozk-OJ-O) = —; find “a (29)
all
(a = x.y.2)
with
2
n 1/2
abc t
Bx B 11' . 2 I“; 2 ex: 2 V2 2 v2
s90 o (a tts) (b t+s) (c t+s)

20

a)

 

abc 4“
3 I 2 9/2 2 1/2 2 us ( 3 O )
o (a +u) (b +u) (c +u)
where we have used the substitution t a s/u. Similar
expressions exist for the other depolarizing factors,

By and Bz. The expression in eqn. (29) is the corresponding
term, for the zone center modes (E=0), of the first term in
eqn. (20) for finite E. For the case a=béc (spherical
geometry), it is easily shown that Bx = By 8 B2 = $9 whereas
in slab geometry with the slab face perpendicular to the
2-axis (a,b 9 m), we have Bx = By = 0 and B2 = 1. Thus by
comparing eqn. (29) with the first term in eqn. (20) in the
limit E 9 0, we see that the longitudinal mode corresponds
to a uniformly polarized array with 30 perpendicular to the
slab face and the transverse mode, to an array polarized
parallel to the slab. Therefore, when E = 0 and hence the
difference between longitudinal and transverse is undefined,
we can then refer to the direction defined by the normal to
the slab face in order to re-establish this distinction.

For the cubic lattices, in the case of spherical
geometry, one can show that the interaction energy,
eqn. (23), vanishes for the case of E = 0. To do so, we

first note that the lattice vectors in these structures are

given by:

21

8.1
R = a (l,m,n)
l,m,n are integers with no restrictions (31)
ESE.
R = ; (l,m,n)
l,m,n all even OR all odd (32)
ESL
R = % (l,m,n)
l+m+n even a all even OR 1 even, 2 odd (33)
where a is the conventional unit cell length. We consider

the ordered triple (l,m,n), all non-negative, and its
various permutations and sign changes, which we designate by
<1mn>. The set of lattice points so defined forms a
spherical shell a distance R = (12 + m:I +n2)1/2 Ro from the
origin, where R0 is a for the SC and a/2 for the BCC and
FCC. (We neglect the case where l=m=n=0.) If we represent

30 by (x,y,2), then the interaction energy of the central

dipole with this shell is given by:

 

 

V = 2 xz+ yzi- 22 _ 3(xl‘t’Lm‘I-zn)a
shell. (lam) R9 R5

22 22 22
= -l; 2;» (xz+73+sz) - 3‘” 1 +y ' :2 n )
R <1 (R/Ro)

 

22

.. ——6——2' 2 (xy 1m 4- xz 1n 4- ya mn) (34)
(R/R ) (l n)
o 0
Now:
(12/30)” = (1’+-’+n’) .. 3 (1' or e’ or n') (35)
(1 n) (l n) (I >

where the last expression follows from the cubic symmetry.
Thus we see that the first sum in curly brackets in
eqn. (34) vanishes. In addition, a typical term in the

second sum can be written as:

x7 2 1m =2 [1m + l(-m)] = 0 (36)
< n)

(Inn)

Thus the second sum is also zero. Since the entire lattice
can be summed up in this manner, we have demonstrated that
for a spherical cubic symmetric lattice, the dipolar
interaction energy is zero. Because this corresponds to the
case where Bx = B = B = - in eqn. (29), we have shown

y 2 a’
that:

9 4n 2
V(k=0) = 53% + Z[f=o]

sphere
3:0

2
u
+ Z[E=o] - in. = 0 (37s)
n
in

“I.

23

where the middle two terms signify the reciprocal and real
space (shape-independent) sums in eqns. (20) S (22)
evaluated at i=0. Thus from the previous discussion on the
connection between the longitudinal and transverse modes and
the slab geometry, the total interaction energy (a a?) for

these mode (with 30 along the z-axis) is given by:

V‘L’dsm = 5’3- p2 + E [#0] (37b)
0 0
slab
3:10

2 2
u u
9 4 o s __ _ 1 =3 _2
+ Z [3:0 ] " 5- —fl 7) -' [‘7' (1 .) 871/3] 0
ho
and
V"’(f-0) = 0 + 2 [2:0 ] (37c)
slob
8:10
' 2 2
u u
+ Z[f=o] - 3—3-7? = [4n(0-§)--4n/3]—3-
1'!
3:0

LMMM;WW
Up to this point, we have taken a purely microscopic
approach and our results have been fairly formal. We now

show how they can be interpreted physically in the context

24

of simple macroscopic electrostatics. In what follows, we
consider the lattice to be composed of a spherical ”near”
region where the dipoles are treated discretely imbedded in
a ”far” region where the material is treated as a
homogeneous polarizable medium with an ellipsoidal outer
surface. We show this schematically in Figure 1. If the
medium is now uniformly polarized, a surface polarisation
charge forms, both at the outer ellipsoidal one and the
fictitious, spherical, inner one. The areal charge density
is given by 3-3, where P is the polarization, in this case
equal to polo, and Q is the unit normal to the surface. For
spheres, this charge density is given by [Pl cos 6, where 6
is measured from the direction defined by the polarization
vector, while for slabs, it is IPI and 0, for P
perpendicular (longitudinal) and parallel (transverse) to
the slab face, respectively. There is no volume
polarization charge as 9-3 = 0 for constant P. The electric

field due to polarization charge at the center of the inner

sphere is given by:

 

9, . 9
it?) = ‘3" ’ “1 (r? - 2') (30)
I? - f‘l’
surface

where the vector surface element is given by, dz - n da.
The integration can be easily carried out for the special

cases of the sphere and slab, giving: E” = + (Mt/3) P

25

 

 

 

Figure 1. Model System for Dynamical Calculations of
Ionic Systems

26

where the minus (plus) sign corresponds to the field arising

from the polarization charge on the outer (inner) surface,
and gang-L) = -4n 3 / Isiah", = 0 where ‘J.’ and <|> refer
to the polarization perpendicular or parallel to the slab

face. Since the interaction energy between a point dipole

and an electric field is given by, -po°E, we see that
eqn. (29) can be physically interpreted as the energy due to
the surface polarization charge. The net electric field
from the ”far" region, which is the sum of the fields due to
polarisation charge on the (real) outer surface and the

(fictitious) inner surface is given by:

(netJar )

sphere 0 (39.)
(net..fo.r)
“aw = ~(8n/3) 3 (39b)
(netJor)
““4. = (411/3) 3 . (39c)

There is also, in general, an electric field (and hence
interaction energy) due to the discrete charges within the
sphere (the so-called ”near” field), however, in the
analysis above we have shown that this vanishes for the case
of cubic symmetry. Therefore the net interaction energy of
the central dipole with all the other dipoles of the system

in this macroscopic approach is given by:

27

vdd - “0 total - - “o . outer + inner) a
r 0 (spherical surface) (40a)
2
A 8n/3 -3 (PL slab face / longitudinal) (40b)
2
“o
L -4n/3 ;—- (P I slab face / transverse) (40c)

 

We have thus recovered the previous results, eqn. (37),
without recourse to the Ewald method, using simple
macroscopic electrostatics. As we shall see, this approach
has a greater validity than demonstrated here, and can be
used to calculate the two most important quantities, the
density of states (DOS) and the dielectric (IR) response,
for disordered systems. This results in a significant
savings in computational time when compared with the Ewald
method. In fact, for the study of the mixed crystal
semiconductors, this allowed us to use the on-site computer

(IBM 3090VF) rather than remote supercomputing facilities.

medtmmmhtmunm

Assuming a harmonic solution, uo(t) = poeudht,

eqn. (7a) gives:

28

w
+ J— 3030) = 0 (41)

[w(E)' - ail;
412p .

0

For general E, -E(3°) is obtained from eqns. (20) and (22),

the Ewald result, with the substitutions: i:}9 Do,
(50- f)'—. (60- Inf. and (Do- §)'-+ (Do- an. Thus

eqn. (41) leads to a 3 x 3 secular determinant (one equation
for each of the three components: (#0)”, (p0)y, (Mo’sL’ from
which the eigenvalues (0(3)2 can be determined numerically.
Results for the three cubic lattices, described by
eqns. (31)-(33), are presented in Figure 2(a-c). They give
the dispersion curves along various paths in the first
Brillouin zone, connecting the symmetry points of each
lattice as labeled on the lower axis. We further identify
the symmetry lines along the coordinate axes, (100), the
face diagonal, (110) and the body diagonal, (111). The
locations of these special points and lines within the
Brillouin zone are shown in Figure 3(a-c) which have been
reproduced from Callaway (Ca 76). The LO-TO splitting is

clearly apparent at the zone center (P) for all three

lattices. From eqns. (39), the net electric field is given
#3 “my
by, -8n/3 —0 [471/3 T] for the longitudinal (transverse)

branch at the zone center. Thus:

(a: - (0:) .. £00: (1.2.)

2 2 _ _ 1
(wk - w5)- . w (42b)

Note that this holds for all three lattices. To get an idea
of the utility of the Ewald method for finite I, we plot in
Figure 4 the longitudinal and transverse branches along the
(100) direction for the simple cubic lattice using the Ewald
method (solid line) and the truncated dipolar interaction,
eqn. (7), for various cut-off radii, in units of the lattice
spacing. We see that the latter all start at zero at the
zone center, for both the longitudinal and transverse
branches, and then approach the Ewald values at varying
rates. For the smallest E, even truncating at ten lattice
spacings, corresponding to over 4000 dipoles, is inadequate.
However, as we approach the zone boundary, the values for
small cut-off (Rmneu'= 2, for example), are reasonably
good. Thus we see that the Ewald method is most important
at small, but non-zero(I), wavevector, For modes at the
zone center or near the zone boundary, the truncated
interaction is potentially still useful. We will see this
empirical rule further validated when we consider the
problem of disorder.

For the disordered dipolar systems, we are most
interested in the density of states (DOS). Thus to provide
a basis for comparison, we calculate the DOS for the three

periodic cubic lattices. We do so by randomly choosing a

30

 

 

 

 

 

 

1-0 E : : E I 3
0.5 r
CL
I.“
-
I.“
0 _.
_0_5 s a s s
r X M r R M X R
(100) (110) (111)

Figure 2(a). Energy Dispersion for SC Dipole Lattice

31

 

1.0

0.5 r

E/Ep

 

-0.5

 

 

 

(100)

Figure 2(b).

(110) (111)

Energy Dispersion for BCC Dipole Lattice

 

32

 

1.0 : f 2 : g ; I

E/Ep

 

 

 

-0.5 i. E 5 i E E E i
(100) (110) (111)

Figure 2(c). Energy Dispersion for FCC Dipole Lattice

 

33

 

 

 

(a) SC

 

 

 

 

(b) BCC

 

 

 

(c) FCC

 

Figure 3. Brillouin Zones for Cubic Lattices

3).

 

1.0

0.5

E/Ep
-\

 

 

 

 

-0.5

Figure 4.

(100)

Dispersion for Simple Cubic (100) using

Truncation and Ewald Methods

35

-9 9 9 9
point in reciprocal space, k = x‘b‘ + xzb2 xsba, where
x‘,xz,x3 are random numbers between zero and one, and
3;,gz,g. are the reciprocal lattice vectors, and then

calculating the (three) corresponding eigenvalues, ”(3)2.
Note that each point in this region can be mapped to a
unique and equivalent point in the first Brillioun zone.
The values are then distributed into one of 150 bins from
-0.5 0): to 1.0 «3:. For the results shown, in the solid
line, in Figure 5(a-c), one million points were sampled for
each lattice. Although all three show a ”transverse" peak
on the left, there are quite large difference between the
simple cubic and the body-centered and face-centered
lattices. First of all, whereas the last two have band
edges at - 3 and E, corresponding to the modes at the .zone
center, the simple cubic band is slightly broader, the edges
due to the longitudinal and transverse modes at the
point (X) in the Brillioun zone along the (100) direction.
Furthermore, both the BCC and FCC lattices have a dip at
approximately 0.12 and a broad ”longitudinal” peak to the
right. And other than the former having a slightly sharper
transverse peak and less broad longitudinal peak, the two
are quite similar. In contrast, the simple cubic lattice
has a rather flat band to the right of the transverse peak
other than a smaller peak at zero which is due to the modes

at point (B) at the zone boundary along the (111) direction.

36

 

 

 

 

 

6 l T V r y
Ewald
—-- NM
4 — _
A
|.I.l
V
Q
a.
ll.)
2 - _
.J x
’ .\
o 1 l 1 J \ 1
-0.50 0 0.50 1.00
E/E
D
Figure 5(a). Density of States for SC Dipole Lattice

 

 

 

 

 

 

1.00

 

E/ep

Figure 5(b). Density of States for BCC Dipole Lattice

37

 

Ep 0(5)

 

 

 

 

1.00

 

Figure 5(c). Density of States for FCC Dipole Lattice

38

(The even smaller peak just to the right we feel is noise as
we can find no stationary points in the dispersion curves
corresponding to this frequency.)

To determine to what extent these structures are

related to the coordination (number of nearest neighbors) of

the lattices, we have calculated the DOS, having truncated
the dipolar interaction at the nearest neighbor. These are
plotted in dashed lines on the same set of figures. For the
SC, the result is quite flat as in the case for the full
interaction. It is symmetric, with a small peak at zero (as
before) but with no transverse peak. At first glance, the
result for the BCC is quite different from its full
interaction DOS, with an extremely steep maximum at the
origin. However, we can consider the hollow to the right of
this peak to be the origin of the dip seen previously and
the smaller side peaks to be related to the (narrow)
transverse and (broad) longitudinal peaks. For the FCC, the
nearest-neighbor results are surprisingly similar to the
full interaction results even to the extent of the band
edges being in almost exactly the same places. Two peaks
are seperated by a hollow of nearly identical depth and
position as before, although the right-hand one is narrower
for the truncated interaction. Thus these NI results are
suggestive of those for the full interaction, 2 point we

will return to when we discuss the model of Pick (Pi 77).

39

mummnnmmdmnmmmm;

Alum“

The presence of periodically repeated structural units
in crystalline solids allows an enormous simplification to
be made in the calculation of the dynamics of harmonically
coupled systems. This is possible because the
eigenfrequencies (and eigenmodes) can be paramaterized by
the wavevector, thus requiring only the diagonalization of a
3s x 3s matrix where s is the number of atoms in the basis.
When disorder is present, be it structural as in a glass or
substitutional as in an alloy, the wavevector loses precise
meaning and we are faced with the diagonalization of an
I x I matrix, where I is the number of atoms in the solid.
In practice, one works with a number sufficiently large to
characterize the disorder, usually a few hundred to a few
thousand units placed within a cube. To minimize the
surface effects, periodic boundary conditions are used. Our
disordered system is thus a periodic one but with a very
large unit cell called a supercell. One constraint imposed
by this periodicity on numerical simulations is that we are

.9.
always working at i=0 since M1(R) =- 121(0)“k R 8 [41(0) where
R is a lattice vector of the supercell. As we have seen at

the zone center, it is imperative that shape effects be

40

included in the calculations which we do so as before
through the use of depolarizing factors.

To diagonalize directly the resulting dynamical
matrices for these supercells is not a trivial task even
with today's supercomputers. The Equation of Motion (EOM)
method, however, provides an efficient means of calculating
many quantities of interest for random systems (Th 86 and
references therein). This is a variant of conventional
molecular dynamics tailored for harmonic systems. The
initial conditions at t=0 are determined by the particular
response function it is desired to calculate, and the
equations of motion are integrated forward in time with
appropriate weighting factors. We discuss this method more
fully in Sec. II-B.

In order for a numerical algorithm to be accurate, the
step size must be small enough to adequately track the
highest frequency, or equivalently, shortest period
oscillation of the motion. One the other hand, the
simulation must be run long enough that the longest period
motion is sufficiently sampled. Thus there is an advantage
in minimizing the maximum frequency as this allows for a
larger time step and hence a smaller total number of time
steps needed for the simulation. As we have seen
previously, the natural frequency of the dipoles, “5’ merely

shifts the spectrum, thus the argument above suggest that we

41

set it as small as possible. Note, however, that if set w:
smaller than the lower band edge (a: E (0:), then (.02 will
become negative at some point. This leads to exponentially
growing solutions (a —9 -1|m|, iwt —9 +|m|t) which cannot be

properly simulated as it quickly leads to overflow errors in

the computation. (Physically, this instability is related
to the ferroelectric transition in which a soft phonon leads
to a new equilibrium structure with not static dipole
moment. However, whereas in our case, m5 is an adjustable
parameter with little physical significance, the
corresponding frequency for a real solid is related to the
physical masses and force constants of the system which are,
of course, given by nature and not easily tampered with.)

If instead we solve a dynamical equation that is first-order

 

 

in time:
t o t. 4flp R 9 R 5
jfli ij ij

(E03 00:, RP: 00:)

with solutions, A“) = rote-mt, then this problem can be

avoided. Note that this equation has the same eigenvalues

DOW

as before since the dynamical matrix, ii

given by:

42

 

a s
an 60-3 3 R1.) Rii
Di) = n—s— " as g 1‘.) (2)
ij 11
remains unchanged by this switch from a? to E. However,
even for negative E, the solutions e‘Et are still
oscillatory and bounded. In addition, the complex

arithmetic features of FORTRAN 77 make the implementation of
eqn. (1) straightforward. The only disadvantage of this
transformation is that the 3's are now complex, essentially
doubling the number of variables, but this is more than
compensated by the decrease in the maximum frequency or
rather energy (Eamg). We set Eo=0, expressing results in
units of Ep 5 mi. Equation (1) is solved using a 4th order

Adams-Bashforth method which is of the form:

_ h - _ a
ym1 - y" + 57.- (San 59fn_‘+ 37f"_2 9fn_n) + O(h ) (3)
where Ens f(xn,yn) and y; = fn and h is the step size. The
previous starting values are generated using the

(self-starting) 4th order Runge-Rutta method (accurate to
order 0(h5) which is described along with the Adams method

in Koonin (Ko 86).

43

meummmmLume

We first consider the density of states as it is
somewhat simpler. Because eqn. (1) has the form of the
time-dependent Schrfidinger equation, we find it convenient
to use bra-ket notation; however this in no way implies that
we are treating the dipoles quantum mechanically. Thus we

re-write eqn. (1), with Eo=0, as:

Q.

1 —— |p(t)> = H |u(t)> (4)

n.

where |p(t)> is a Bid-dimension complex (Hilbert space)
vector and H, in the site basis [1,00, is given

by -(Ep/4np) 0:3 In this basis, we can write the initial

1‘
state |p(0)> as:

|u(0)> = 2.1,“ [1,00 (s)
1,a
where the a1 a 's are in general, complex. We can introduce
9
another basis, the energy eigenbasis, IE"), which we can

relate to the site basis using the completeness relation:
2 "3“"! = 1 (6)
h

where 1 stands for the unit matrix with the same dimension

as the energy eigenbasis. Thus:

44

|1,a> e 2<En|1,a> IE") = Zhv‘f'“ |tn> . (7)
k k
Since:
d
1 a? |En(t)> = u |ln(t)> = 2" Iran» (a)
we have: |En(t)> = e-iEnt IE“). Thus from eqns. (6-8):

-i.E t
|p(t)) = 2e 0 It") “J“‘m’ (9a)
n
_ (1,a) -iE t
-12G.1902b" e m 'En> (9b)
’ n

If we now take the inner product of this state |u(t)> with

its value at t-O, we get:

G(t) a <u(0)|p(t)> (10a)
= II Ni',a') (1,00 -tE t
e13“. '1,a2hn h“ e n (10b)
1,a n
10’“.

where we have used the orthogonality condition,
(1042") a 6",”. Choosing a1,“ - expuai’a) where 91’“ is
a random number between 0 and Zn, we have

< " ' > a .5 6 h th b k t i if
a1,“ .1,a “v 1,’1 a',a‘ w ere e rac e s s gn y

the average over many runs (random initial conditions).

45

Taking the Fourier transform, we get:

m
«(2» a J- 1 dt e"Et (G(t)) (112)
a CV

av
'0

Q
= 2“: .2“’°‘" béi’m] 7:71“: e“"'»”} (11b)
1
n

a
’ -a>

a)
_ 1 ((24 )t _ _
_2[§;J-dte n ]-26(r sh). (11c)

n '0 n

which is the desired density of states, where we have

used (1):

3 5 . (12)

which expresses the orthonormality of the eigenvectors of H
(and follows from the completeness of the site basis, |1,a>,
and energy basis, IE“), vectors.) and (2), the integral
representation of the Dirac 6-function. In practice, large
systems are self-averaging and it is unnecessary to average
over different sets of random numbers.
Two technical points need to be addressed. First, it

is not possible, of course, to integrate starting from -u)in

a computer simulation. Note, however, that from eqn. (9a):

46

eiEnt

|u(-t)> =2 ".3 <3n|p(0)> (13)

therefore: <u(0)|p(t)>' = <u(0)|p(-t)> and

O Q
G(E) . .51,- [I dt eat <u(0)|u(t)> + Idt eat <M(o)|u(t)>]
-1»

Q O
= 271? [I dt 6““: <p(0)|u(-t)> + Idt em <u(0)|u(t)>]
0 0

Re Idt e)“ <u(0)|p(t)> . (11.)
0

u
al"

Secondly, since it is impossible to integrate forever, the
upper limit should be replaced by T, the total length of the

simulation. Because I is finite, the delta-functions in
eqn. (11c) are replaced by the broadened resolution

functions:

e
A'(E) .. 2 Re Jdt e“:t (15)

e
which is peaked at E80 with ripples on the "wings". Thus

the precision of our calculations are limited by the finite

width of these resolution functions much as experimental

47

measurements are limited by the finite resolution of the
intruments. To dampen the small oscillations on the tails,
we actually calculate:

r
. _ 2
G(E) - %.Re J'dt e‘nt e h(t/T) (p(0)|u(t)>, (16)

0

where h is a smoothing parameter chosen, in our case, to

be 3.

The other quantity of interest is the dielectric
response. This is the response of the system to a
harmonically oscillating, peturbative, electric field. It
is a function of both the frequency of oscillation (or
equivalently for our first-order equation, the frequency
squared (3 E2)) and the wavelength of the displacement

pattern and is given by:
“3.00) = so + (audio) (17)

where sm’is the low-frequency electronic response of the
solid. The infinity refers to fact that this is
high-frequency when compared with typical phonon
frequencies. We consider our dipolar system to be a pure
insulator, therefore we set sm'= l, the dielectric constant
of the vacuum. The problem then becomes one of determining

the susceptibility, x(E,w), which relates the polarization

48

(the response) of the dipolar system at a given frequency
to the external electric field (the probe) oscillating at

the same frequency:

pas-210.33”. (1s)
(a

(Note that we denote the electric field by !, reserving the
letter E for the energy.) We use linear response theory
following the derivation presented in Thorpe & de Leeuw
(Th 66) but re-written in bra-ket notation. (See also
Maradudin et a1., Chap. 6, (Ma 71).) The equation of motion

for the perturbed system is given by:

Et

'9-

i |s>e"‘ = no |s>e"‘" - Z 33 Ij.fi>e-ut (19)

1.3

O-

t

where H0 is the original Hamiltonian with eigenvectors,
IE“), and eigenvalues, In. Multiplying on the left by (Enl
e -'tEt

and eliminating the common term: , we get, after some

simple algebra:

“..I’J”

"| 2 r9 (r -n
J.” h

Multiplying on the left once again, this time by <i,a|ln>,

and summing over the indice n, we get:

49
<i,a|[2 |2n><sn| ]|r> = <i,e|1|r> = <i,a|E>

<i,a|ln><En|j,fl>

= 2 2 8,, . (21)

n j,” (En-E)

 

The polarization is given by:
i 1
Pa - 5211051) = 52<i,a|l>
i i

<i,a|ln><En|j,fi>

(E -E)
n

 

DP‘
DV1

___ 2 ‘3 . (22)
(3

i J

9
fl
Comparing with eqn. (18), we therefore get the desired
result:

<i,a|Efl><En|j,fl>

x . .. (23)
a” “ (sh-s)

 

The systems we are considering are either cubic or
random and therefore the susceptibility, at first glance, is
isotropic and characterized by a single scalar quantity.
Recall, however, that for i=0, the response is
shape-dependent and that the transverse and longitudinal
responses, corresponding to polarizations parallel or
perpendicular to the slab face of a slab geometry, are

different. (Recall that the sample shape is incorporated in

50

the potential via the depolarizing factors.) Therefore the

susceptibility tensor is diagonal but with two distinct

elements: 1" = xyy E x, and x” a 1L. (For a spherical
sample the response at E is isotropic and

z a x z a 2..) Each of these elements has the form:

X8 yy

<i,a|En)<En|j,a>

 

xm = 52 (24.)
(s -s)
is.) n
n
z[2bd °° h” °°]+— . (24b)
-3)

However, because D?? is real, so are its eigenvectors, given

(1.6:)

by: (tn 1. Furthermore, we are interested in the

imaginary (absorbtive) part of the suceptibility, thus:

Im x” a a: [2 h“ '°°]I- —;— —2) (25a)

2
= 3.2 [2 h;"°°] sun-s), (25b)
n 1

where we have used the formal relation:

n 6(E -s) = lim (3 - r -i.r)-‘ . (26)
790

To actually calculate this quantity with the Equation

of Motion method, we give each dipole an identical initial

51

displacement along the direction of the field. In bra-ket

notation, we have:

Ga (27)

a1,” = <i,n|u(0)> - 6

Therefore, from eqn. (10), we have:

3 .
G(t) = <u(0)|u(t)> =2 [ 2 a1 fl b;1’n’] eflint
1 9
h

2 .
=2 [Zb;1,w] e-‘tlnt (28)
n 1

and
O
c(£)=}7aeIdt e"it G(t)
0
(ion 2
-2 [2b ' ]6(:-:)
n n
n 1

25(E) * 0 (aax,y)

e (29)
x£(E) * 0 (asz)

where the last line follows from eqn. (25b) assuming the

normal to the slab face is along the z-axis.

To obtain 2a(E,l) for non-zero I we set
.. . ,9 ..
..1th II ka e3 1.1 where k is a unit vector along k’: to obtain
9
AL ‘3‘? A
aa(3,l) we set a1,“ I ka e .i where k, is a unit vector

52

perpendicular to E. In both cases, 2 = (2n/L)*(n‘,n',n.),
where L is the supercell box length and the n’s are
integers, and :1 is the postion vector of the ith dipole

within the supercell. Thus, to summarize, we build an

initial state |u(0)> (see eqn. (5)), calculate (numerically)

its time development |u(t)>, and take the Fourier transform
of G(t) 5 <u(0)|u(t)>. We next consider the implementation,

within the numerical simulation, of the Coulomb interaction.

El.Ianlsmsnlalisn.2£.ths.anlsmh.Intszasiian.

Consider a cubic supercell of length L on a side
containing I identical point dipoles. We denote the lattice
vectors by K I L(n‘,nz,n.) where the n's are integers. The
displacement vector between two dipoles in the same
supercell, say i and J, is given by i?“ 2 3.; :5.

To implement the disorder, we construct perfect (cubic)
lattices ofll'sites within the supercell which we then
populate with point dipoles with a given probability p.
This allows us to interpolate between the two extremes of
periodicity (p81) and completely uncorrelated randomess
(p —0 0, I —9 (n, with 1]) finite.) It should be emphasised
that we are not decreasing the total number of dipole by

decreasing p, but rather increasing the number of sites each

dipole has available and thus the disorder. We adjust an

53

the number of sites, so that the total number of dipole in

the supercell, M remains fairly constant, typically

at
between 160 and 200. Furthermore, note that from eqn. (1),
we see the ”coupling constant”, in units of Ep(m mi), is
given by Hum)”. In this way, the results for different p
are directly comparable.

The dipolar interaction energy of the ith dipole in a

given supercell with all the other dipoles in the system is

 

given by:
)3 3. 31?: (f..+3)ltfi.-(? +3):
dd(i) = ‘ ’ ‘ ” ’ u (30)
9 '9 s 9 9 o
[rt j+n| Ir..+n|
n ‘J

where the prime over the sum over 1 indicates that the i=1
term for 330 is excluded. We can, as before, re-express the
slowly convergent sum into rapidly convergent ones. Because
the derivation is similar to the one for the periodic
lattice, we will not present the details. The interested
reader is refered to the article by de Leeuw et a1. (dL 80).

We merely quote the results.

 

-. __-_,_‘.’._j _ 313,- in.) 13,- 3“) 1 x
“m ZZd—L—H Ii 1' I3 I“

n i=1 " ‘5

54

2
[7:11.ti e "7315’ 9 £2! erfc(nRij)] -

 

2
19-.°mn"')(i5-3)(35-i) --‘3:u'
*3- R' i. ij j 15 316? 1
ii
I
d
«uh/n)z ..
+ :33 Z 1‘ (35.; I) Z (333 if) e": ”u
k’flo i=1
3
4n 2
+ ?)— 28°! (“1’01 (31)
asl
where it) a fir-K, It”. 2 lit)" 5 2 (ill) (n‘,n.,n.),

951.7, and not are the depolarization factors defined
previously. For our simulations we have used n I 3.0 and
have truncated the real-space sum at RU < 0.75 and the
reciprocal space sum at [3| < 3.5.

The results of Sec. (III-C) suggest that we can,
alternatively, truncate eqn. (30) at some radius, Re, about
each dipole, treating the other dipoles discretely (whose
net contribution will in general not vanish, unlike the case
for the cubic arrays) and augment this interaction with the
electric field from the ”far" dipoles. (See Figure 1.)
This field, as we see from eqn. (39-III), is proportional to

P, in this case, the polarization averaged over all the

55

dipoles in the supercell (= {5 {til L'). In essence,
therefore, we are handling the ”far” region contribution of
the Coulomb interaction using a mean-field approximation.
As we shall show, this approach is quite satisfactory for
calculating two of the more important response functions,
the density of states and the dielectric (IR) response in
the long wavelength limit (3-0) for these random systems.
It should be stressed that in both approaches, lwald and
truncated, the interaction of a dipole with all others in
the system have been included.

We test for the equivalence of the two implementations
of the Coulomb forces (the exact Ewald method and the
approximate truncated/auguented interaction) on a simple
cubic lattice with J'8 6' = 216 sites in the supercell
choosing p-0.9. We feel this provides a rather severe test
of the approximation as it includes a large number of low
symmetry configurations. Simulations were run, on a single
configuration, for ten periods, where one period is the time
needed for an oscillation of energy IP to be completed.
time steps equal to one-fifteenth (:5) of a period. (cycle)
were used for the DOS and the 3-0 dielectric responses. For
251(3,“ smaller steps were needed to insure numerical
stability, i th and i th of a cycle for the truncated and
Ewald runs respectively. for the runs using the truncated

method, BC was chosen to be 0.‘ L(- 2.‘ s) and 0.5 L('3.0 a)

56

where a is the nearest-neighbor distance of the simple cubic

lattice sites. These radii contain 52 and 102 dipoles,
respectively. The truncated runs ran at least 50 times
faster than those using the lwald method. This might be
overstating the time advantage somewhat as the Ewald code
was not fully optimized, although it was partially
vectorized. (The truncated method was applied successfully
to ionic systems (with discrete charges) by Boltjes & de
Leeuw (Ga 89) although the time improvement was not nearly
as dramatic. As we did not review the code in this case, we
cannot explain the reasons for this. lonetheless, a quick
comparison of eqn. (41) with eqn. (40) reveals that the
truncated method is far simpler to implement and code.)

The results are shown in Figures (6-10). There is some
discrepency in the 00$ for the Re a 0.‘L run but by
RC 3 0.5L the convergence to the Ewald result is nearly
complete. (The overall "peakiness" is due to the finite
size of the model. It is the comparisons between the three
results that is of interest to us here.) For the 2.0
dielectric responses, the finite size is of less consequence
and other than the case of the longitudinal response,
Im 5(1) where the truncated results are slightly narrower,
the agreement between the two techniques is excellent.

A more serious discrepency arises for Tm st(t,l) for

3*(2n/L) h1_and is analogous to the problem with truncation

57

that we found when calculating the dispersion for the simple
cubic lattice with the same method. (See Figure A) In the
truncated results the energy of the mode is underestimated
and the response broadened. however, as Re is increased,

the results improve and for sufficiently large radius would
eventually agree with the Ewald result. In Figure ll we
show the truncation results for §8(2n/2L) ii! A 8 2L. There
is a smaller change from the Re I 4 to the Re a 5 result,
indicating that the convergence to the true (Ewald) result
is farther along. As we see in Figure 6, this is also the
case for the simple cuhic results.

This failure of truncation suggests that correlations
between dipoles separated by more that three bond lengths
(the size of Re) are important for the small wavelength
modes. Note that for non-zero : the net polarisation nearly
vanishes and it is the ”near” region that provides most of
the interaction. As : increases, the more rapid alternation
in the phases of the dipoles leads to a faster cancelation
of the interaction from more distance ”shells”. For the
008, the random phases of the dipoles also leads to the
vanishing of 3 and in this case, the cut-off procedure is
even better. In the case of the dielectric response for
2.0, this polarisation, and the macroscopic field it gives
rise to, is quite important as it leads to the splitting.

In addition, the success of the truncation method shows that

58

the variation of this field throughout the supercell is
quite small and is accurately represented by the ”average”
given by F. Encouraged by these results we now examine more
systematically, the effects of disorder on the density of

states.

59

 

 

2.0 r I l
Ewald
/‘ ——r-05L
. I —-—r'-0.4L .
“7 /
K“1.0? I \ \ -
“I

 

 

0
-1.0 -0.5 0 0.5 1.0

Figure 6.

 

 

 

 

 

 

 

9 l I l
n Ewald
6 n. — - P . 0.5 L ‘
'— —-— P I 0.4 L
E
:2?
.5.
3 - _
O —- _J
-1.0 -0.5 O 0.5 1 0
E/Ep
Figure 7.

Im 11(E) for Disordered Simple Cubic (p30.9)

60

 

 

 

 

 

 

 

 

 

 

 

 

9 T l l
—— Ewald
__... — r - 0.5 L
__—_- FIDOu‘L
6 P —J
E?
‘3
2.
.5.
3 _ ..
o 1 ' ‘
-1.0 -0.5 0 0.5 1 0
E/Ep
Figure 8. Im s.(E) for Disordered Simple Cubic (p80.9)
9 I l I
6 _ .__-— r-IOJSL
ET.
2‘.’
.5.
3 .—
o ‘ '
-1.0 -0.5 0
E/Ep
Figure 9.

Im 53(E) for Disordered Sisple Cubic (p80.9)

 

61

 

lmlflk. 5)]

 

 

 

0
-1.0

 

.’
Figure 10. Im 8l(k,E) for Disordered Simple Cubic (p80.9)

f-(zn/L) i

 

 

 

 

 

 

 

 

\
m...“ l
'— 5p. ——-rmax-0.5 I‘-
W I
“~53 ,’ I
.§ 3 — I )-
I
l
91.0 -0l.5 40 0.5 1.0
E/Ep

Figure 11. Im clan!) for Disordered Simple Cubic (p80.9)

32(2n/2L) E

62

Lmummer.

In Figures 12(a-c) we give the results for the density
of states, calculated via the truncation method with

RC = 0.5L and averaged over ten configuration, for the three

cubic lattices. He chose p=0.75 (—————-), p=0.50 (——— -——),
p=0.25 (—— ° ——), and for conparison, give the results for a
completely random dipole gas (-— -<-) in which the
equilibrium positions were chosen at random. The disordered
DCC and FCC lattices have nearly identical DOS, the
”transverse” peak remaining identifiable even at p-0.5 and
the ”longitudinal" shoulder having filled in the dip that
existed at p=1.0. Together with the SC, all show a gradual
broadening and, not unexpectedly, all approach the random
result as the number of possible sites becomes larger and
the discreteness of the lattice less important. Eventually,
symmetric tails form, as seen in Figure 13 for the case of
p-0.05 and the DOS for the three lattices have nearly
converged to that for the dipole gas. (Also note the nearly
factor of five increase in the bandwidth with increasing
disorder.) These tails are a result of isolated,
strongly-interacting dipole pairs which become more
plentiful as the lattice becomes more disordered. We can
explain this quantitatively by diagonalising the 6 x 6

interaction matrix of a single pair with 3 8 Rs.

63

 

 

 

 

I I T I
—— p-0.75
-- - 9-0.50
—r— 9-0.25
—-— PDHOOIII
A 2 r- ‘fl
Ill
V
Cl
0.
I.l.|
d
1.5

Figure 12(a).

 

DOS for Disordered SC Lattice.

 

 

 

 

 

T I l T
-— 9-0.75
— - p-0.50
—— p-0.25
—‘— PDHOOR
A 2 '-
Ill
V
O.
0.
III
1 L—
/
P’-
o A L
-1.5

Figure 12(b).

 

DOS for Disordered BCC Lattice

64

 

3 T I I I I
--—-—— p-0.75
— — 9-0.50
—— 9-0.25
—--— random _
A 2 -
|.I.l
v
Q
a.
u:
1 .—
/./
-l'f

 

 

 

 

0
-1.5

Figure 12(c).

 

DOS for Disordered FCC Lattice

65

 

 

(31)
I

where 0 represents a 3 x 3 matrix with all the element equal
to zero and A is a 3 x 3 diagonal matrix with A“: A”: 1
and A”: -2. The eigenvalues are AL 2 -2,2 for the
in-phase,out-of-phase longitudinal modes and A, a -1,1 for
the doubly degenerate (3 along § or ;), transverse modes,
where we have given the values in units of (lp/dnpi'). Thus
the spectrum consists of a symmetric set of 6-functions and
that for a given mode, there exists a one-to-one
correspondence between the energy and the radius of

seperation, that is:

h I,
E (R) = . (32)
A Anna.

 

We can therefore write a relation between the density of

states from a single mode, px(l) and the density of dipoles

a given distance away, p(R). He have:

p(R) d3 = “Px(l) d2 (33)

where p“) I lmlzp (p is the number of dipole per unit

volume) and the minus sign reflects the fact that the energy

66

decreases as the radius increases. Therefore
px(l) = -p(R)/(d!/dR) where from eqn. (32),

dF/dR = -3x FP/(énpi‘). Thus:

(kapn') ' XI
px(2) = -——————- = -—§ (3‘)
3X 33
where the above result is for a single mode. The actual

density of states is an average of all the modes of the pair
of a given sign, thus :9 p(3) 8 RP (px(l)> I (k) / 3(1/39)’,
where (A) 2 3 x 1 + a x 2 = 3. (There are 2 transverse
modes, x, = 1 and 1 longitudinal mode, AL I 2 per sign, out

of a total of 6.) Our final result is then:

22’

:9 p(E) = ;;§ (35)
which we plot in Figure lb along with the results for the
dipole gas. The agreement is quite good even for relatively
small I but breaks down eventually as the probability of an
isolated pair with large B, corresponsing to a low I,
vansishes as larger dipolar units come to dominate the
spectrum. Equation (35) is the one analytical expression we
have succeeded in applying to the output of our computer

calculations because it involves solving the tractable

problem of a dipole pair. The three-dipole problem, and

67

 

 

 

 

 

 

 

 

 

 

 

 

0.50 I I T I I I Y
i
i
g —— random
A I — —' SC
‘19 —-— acc
Q ---— FCC
a. 0.25 - %
Ill
I
: /'
0 i L L 1 1
-10.0 -5.0 0 5.0 10.0
E/Ep
Figure 13. D08 for Disordered Cubic Lattices (p=0.05)
and Random Dipolar Gas
0.50 I .lgr I
i ‘ Pandas
93 l I
Q \
a \
Ill
-10
Figure 14. D08 for Dipolar Gas with Lorentzian Tails

from Dipole Pairs

68

those involving larger clusters, are in principle solvable
but we no longer have a simple relationship between a single
radius and the energy, which was the key for the dipole
pair.

Although we have been considering the random dipole gas
where a continuous range of radii exist, we can apply the
pair analysis to the disordered discrete lattices. For
example, the prominent peaks that exist on the tails for the
simple cubic in Figure 13 are a result of the. longitudinal

modes for dipole pairs seperated by the nearest-neighbor

distance. From eqn. (32) we have:

 

2 Sp
E (R ) = . (36)
1. an 3
(mp (Rm)
However, p = N/V =p luau/V and (RM)3 = VII-u“, thus:
1
E (R )/E = -— a: 3.18 (37)
I. rm 9 2119

for p-0.05. The peaks due to the transverse modes (at
Flip z 1.59) are somewhat obscured by the shoulders and
indeed are responsible for the deviation of the simple cubic
DOS from that of the dipole gas on the left-hand side. For
the DCC and FCC lattices, I!“" = sf3/2 and R 8 elf?

“h

respectively, where a is the linear dimension of the

69

(ICC)
l .
mites

conventional unit cell. Also, a 2V/sa and

('00)
sites

= 4V/a'. Therefore the longitudinal mode for the

nearest-neighbor dipole pair occurs at:

 

1
thump/I: = — a: 2.12 (ace) (38a)
P 3np
1
EL(Rm)/Ep = up a: 2.25 (FCC) . (38b)

These peaks are also seen in Figure 13 although masked by
the asymmetry of the DOS for these lattices that persists to
some extent even for large disorder, p=0.05.

He close with a comment on the previous results of Pick
and Yvinec (Pi 77) and its relation to our work. Although
they considered the dynamics of point dipoles, their model
started from a molecular dynamics simulation of liquid
argon, interacting via a Lennard-Jones potential, which was
then quenched to give the equilibrium positions of the
dipoles. Because of the excluded volume effects due to the
hard-core repulsion, their results, obtained using a
truncated dipole interaction with no ”extra” polarisation
term which as we have seen is quite adequate for the DOS,
differ quite profoundly from our calculation on the dipole
gas in which there is no restriction on the minimum

separation of the particles. And as we have seen, the

70

presence of the Lorentzian tails, missing in the Pick
results, shown in Figure 15 (reproduced from Pi 77), is a
directly consequence of this. However, a liquid, to a first
approximation, has a close-packed structure with twelve

nearest-neigbors, as is also the case for the FCC lattice

which indeed does exhibit a two-peak DOS, even when only
nearest-neighbor dipole interactions are considered. The
effect of disorder would be to smooth the peaks and give a
result much like Pick's or our own disordered FCC lattice.
(See Figure (12c) for p=0.75.) Although for identification
purposes it is useful, as Pick has done, to describe the
sharper low-frequency peak as belonging to the ”transverse”
mode and the other to the "longitudinal" mode, the splitting
is an attribute of the dielectric response 8(w) not the DOS.
Note, for example, that in the case of the simple cubic,

where the spliting at the zone-center is quite clear, the

DOS shows no clear two-peak structure but rather an
off-center "transverse" peak with a slowly decending
plateau.

In this section we have concentrated on the DOS as this
shows the overall effect of disorder on the system. In the
next section, we will consider the more selective dielectric
and Raman response and the effect of disorder on the Coulomb
splitting. We do so for an ionic mixed crystal system with

discrete charges, studying, in particular, the effects of

7l

..ZJO

 

0.90 0.95 1 we no Inf-F;

Figure 15. DOS for Quenched Dipolar Liquid (Pick)

72

diminishing ionicity on the LO-TO splitting as
positive-negative charged pairs are replaced by neutral,

covalently bonded atoms.

73

memmug_xuuaimm

A... W111.

Recently there has been a good deal of interest in
mixed crystal semiconductors of the type (III-V)tmIva
where a cation is taken from column III of the periodic
table, an anion from column V, and a neutral (covalent)
element from column IV. Examples of these include
(GaAs)1_xGezx and (InSb)‘_xGeu. In the mixed crystals that
have been more extensively studied, such as the alkali
halide labflxx Cl or the III-V semiconductor Ga Pvasx,
either the cation or anion of the pure crystal is
substituted with another in the same column of the periodic
table so that the basic character of the material reaains
unchanged (see the review article by Chang Ch 71). However,
in the (III-V-IV) mixed crystals, we go from a (partially)
ionic (III-V) semiconductor to a purely covalent (IV) one.
Accompanying this shift is a structural transition in which
a zincblende structure abruptly changes to a diamond
structure as the order parameter, defined as the number of
cations on one sublattice minus the number on the other,
goes to zero. There has been much controversy about the
exact structural nature of these alloys and various models

have been proposed for their formation and growth, each

with a different value for xc, the critical concentration of.

74

the transition. While not addressing these concerns
directly, the techniques we have developed for calculating
the dynamics of disordered ionic systems is well suited to
the study of an important aspect of this problem as there is
experimental data by McGlinn et a1. (Hc 89) suggesting that
the LO-TO splitting in optical response vansishes at this
transition point. We perform our calculations on a dimer
percolation model that exhibits precisely the structural
transition that occurs in the actual materials, although at
a critical concentration different from the experimentally
determined value. Nonetheless, this model system displays
many of the same optical features found experimentally
giving us confidence in the applicability of our techniques.

A crystal of pure GaAs or InSb has a zincblende
structure, in which the lattice points form a diamond
lattice, made up of two FCC sub-lattices separated by the
basis vector: 3 = (a/4) x (x + ; + ;) where a is the side
length of the conventional unit. (See Figure 16 reproduced
from Us 70.) All the cations are on one sub-lattice, call
it the A sub-lattice, and all the anions on the other (D
sub-lattice). As x increased, clusters of GaAs (or InSb)
form,. separate from the percolating infinite cluster,
isolated by a ”buffer” of germanium. These clusters can
have a different registration from the percolating cluster,
that is, cations on the B sub-lattice and anions on the A.

The surrounding Ge prohibits the occurence of any Ga-Ga or

75

 

 

 

 

 

Figure 16 Conventional unit cell of zincblende structure

76

As-As nearest neighbors which are energetically unfavorable
because of the large Coulomb repulsion between identical
species. At a critical value of x, the largest cluster of
the III-V material no longer percolates and indeed all

clusters of this type become isolated. There is no longer a

prefered registration scheme, defined previously by the
infinite cluster, and hence there is an equal probablity
that a given cluster is either phase (A-B) or anti-phase
(B-A), where the first letter gives the cation sub-lattice
and the second, the anion sub-lattice. In the zone-center
optic modes, for the perfect crystal, the two sub-lattices
”beat” against each other, and are in phase with the
macroscopic electric field that is set up during the
oscillation. With isolated clusters of mixed registration,
charges of both types would be on a given sub-lattice which
would destroy the coherence, leading to a degradation of the
electric field and hence the splitting. Thus what might be
called ”dynamical frustration” occurs as the mode pattern
favors one motion while the electric field coupling forces
the opposite. This non-rigorous arguement leads to the
ansata that the splitting disappears not with the ionicity
(xul), but rather when the (III-V) material no longer
percolates, which occurs in our model, as we shall see, at
xc = 0.69. Using the equation-of-motion technique to
compute the optical response, both IR and Raman, we verify

this for a system with one-mode behavior, roughly

77

corresponding to (GaAs) Ge, and also for a model material
based on (InSb) Ge which exhibits two-mode behavior. By
these terms, we mean that the net optical response of the
first is given by a single (broadened) peak, whereas in the

second, the InSb and Ge structures (peaks) are distinct, at

each concentration. We begin with a discussion of the
implementation of Coulomb forces for disordered ionic

systems.

anmmnmlminnunmmmnm

The ideas that were developed for the dipolar systems
carry over with only minor changes to ionic systems. Most
importantly, the long-range Coulomb interaction can once

again be separated into a near region where the discrete

ions are treated exactly, and a far region that is
incorporated through the average polarization of the
supercell using the proper depolarization factors to

introduce the macroscopic shape. Coworkers at the Univesity
of Amsterdam, 3. BoltJes and s.w. de Leeuw, have computed
the density of states and dielectric response (for both :30
and 3 finite), for a structurally disordered AX: ionic
glass, modeled on silica, and have obtained essentially the
same results as ours for dipolar systems (Ga 89a). That is,
the DOS and 0(a) can be calculated quite accurately with the

truncation scheme giving the same results as the Ewald

78

method for a sufficiently large cut-off radius, on the order
of four bond lengths. Likewise, the truncated results for
s(k,w) for the saallest possible non-zero f converge,
somewhat slowly, to the proper Ewald result. (For details
on implementing the Ewald method for ionic systems,
see (dL 80).)

We consider a cubic supercell of linear dimension L
containing N ions of charge qi (i=1,2,...,N) with

instantaneous positions ? For the mixed crystal we are

1.
considering, qi is either +1, -1, or 0 unit charge.
Electrical neutrality is always maintained by replacing a
dimer of (III)o”(V)b° by one of Gez. These supercells
form a simple cubic lattice described by translation vectors

3=(nx,ny,nz)L. The electrostatic interaction energy of the

ith ion with the other ions in the "near” region is given

 

by:
Q
- i
vnmar - ‘12.: + 3 I (1)
J 11 ii
where 315 gives the difference in the equilibrium positions

..

of the ith and jth ions and u11 = 31 - 3 the relative,

J
instantaneous deviation froa equilibrium. Fxpanding in a

power series of the small parameter (u/R), we get:

79

 

 

2 '9 2
q u 3(u - R )
1 11 _ 11 11 s
vnear = -q1 27 R, R5 ] + 0(a) (2)
J 11 11

where R11 8 |R11|. Because the equilibrium positions are

chosen so that there is no net static force on any of the
ions, there is no resultant linear term. The total
effective force ?1 on the ith ion due to the other charges.
both in the near and far regions. is equal to

F1 = -61 Vtou" and given by:

 

+ +
i? ___ q {:3 “11 _ ““11 in’i‘u]
1 1 “‘

2 s s
1 “11 R11
3
4n 1
+2—2[§'-Ba]?a} (3)
4m
asl
where
l
l +
Ia'az‘lj“: (4)
i=1

is the dynamical polarization of the supercell. Note that
the first ters in eqn. (3) is very similar to the expression
for the point dipole interaction, falling off as R-.. We
choose a cut-off radius of about three bond lengths:

eighty-six neighboring particles are included in the ”near"

region lattice sum.

80

LWflMLWW

We now discuss briefly the Equation of lotion technique
used to calculate the optical response of the system. In
these simulations we use the usual second-order (Newtonian)
dynamical equations which lead to slight modifications in
the formalism. The details can be found in Thorpe &
de Leeuw (Th 86). We merely quote the main results. The
central entity of this technique is the Green function

defined by:

 

 

1 (“um :guml
G (i 1'») = (5)
an ’ ’ a/s :2 2 _ a
('1'1) w Kn

where t and A are the eignevectors and eigenvalues of the
dynamical matrix of the disordered supercell. The m’s are
the particle masses. (In the previous section on the TO!
for the first-order equation, we were also dealing with
Green functions. We just never called them by their proper
name.) Both the density of states and the dielectric
response, or more precisely, the susceptibility x, can be
written in terms of this function. From Sec. II of Thorpe 6

de Leeuw, we have,

- .. ESL .
pm) - 3“ 1):“ m1 Im cw(1,1.m) (6)
0

81

”(w) =

DP‘

io‘qz Wu 1; co) . (7)
1a,

Note that the results are now a function of frequency rather
that energy (a up). Because we are interested primarily in
the imaginary (absorptive) part of the susceptibility, we
need only the imaginary part of the Green function for both
these quantities. It is shown in Sec. VI of (Th 86) that

this imaginary part can be computed numerically from:

r
z
_ a -X(t/T)
G(w) - "J;dt cos(wt) e 12; A1,“ “i,a(t) (8)
0

where u1 a(t) is the displacement along the nth axis of the
9

ith particle. Also mentioned earlier, the smoothing

factor A lessens the ripples in the convolution function

introduced by the finite T truncation. If we let
Ai,a ani and “i,a anil'i’ where E is a unit vector in

the direction of the electric field of the incident light,
then we get, using slab geometry implemented via the

depolarizing factors:

Im (stag) - a: m] = g— G(w), E parallel to slab (9a)
2::2 ‘
Im [6"(0) - em] = 53— G(w), E perpendicular to slab. (9b)

The Raman response of a material is a function of the

polarizability of the electronic- degrees of freedom in

62

contrast to the ionic polarizability that gives rise to the
dielectric (IR) response discussed previously. As such, a
fundamental expression for it involves the complicated
quantum mechanical calculation of the electronic structure
of a solid. For our purposes, a phenomological approach
will suffice, based on a local bond-polarizability model of
Alben et a1. (Al 75). The perturbing interaction between
the external radiation and the phonon is described by the

polarization tensor:

9 2‘1) [’15 ° ’11 " 31] “1 ’11 (1°)

i,j
where r11 is the unit vector pointing from site i to site 1,
1 is the 3 x 3 unit matrix, and 3 is the displacement of

i
the atom on site i. Note that the symbol 0 signifies the

outer product of the two unit vectors, giving, of course, a

tensor. In general, the coupling coefficient A11, can
depend on the bond, however for simplicity, we take A11 8 A
to be the same for all bonds. Thus in the perfect

tetrahedral arrangement that exists in the zincblende (and
diamond) structure, the second term in brackets vanishes.

That is:

2r“ =0 (11)
J

for the unit vectors given by:

83

3‘ = :(1/73) (+1,+1,+1) (12.)
$2 = :(1/13) (+1.-1,-1) (121)
:3 = :(1/15) (-1,+1,-1) (12.)
;‘ a 2(1/13) (-1.—1,+1) (121)

.where the two signs correspond to the two sites in the
zincblende basis. For this same reason any diagonal element
of the polarization tensor (which gives the response of the
system along the same direction of the incident radiation

field) also vanishes. In contrast, for the off-diagonal

elements, say 01“)", we have for the contribution due to the
ith site:
(xy) (x) (y) (on (on + (on
011 or Zrij r“ r:L1 “i at _6“ u1 (12)
j

where the plus (minus) sign applies if the ith site is on
the A (B) sub-lattice. Experimentally, this corresponds to

the arrangement where the scattered radiation field is

measured perpendicular to the incident radiation. In

practical terms, for the ROM, this corresponds to
8+ 0 3+

A1,“ '6a: and ui,a ‘6asf'i’ thus even in covalent

materials such as germanium with little effective charge,
the Raman response is non-zero. lote that the induced
motion of the atoms is along the direction of the incident
beam and thus the Raman technique gives (primarily) the

longitudinal response of the system. This complements the

84

IR response where the particle displacements are along the
electric field direction and hence transverse to the
wavevector. The EON method, however, can compute either the
longitudinal or transverse response for either situation (IR
or Raman) simply by setting the initial displacements either
along or perpendicular to the normal of the slab face as

defined by the depolarizing factors.

The equations of motion are numerically integrated
using the Verlet algorithm (Ve 67) that can be easily

derived using Taylor expansions. We have:

x(t+At) = 11(1) + 11:11)“ + 311””) (At)2 + out)’ (13.)

x(t-At) = 11(1) - x'(t)At + ix"(t) (At): - out)”. (131)

Adding these two equations and solving for x(t+At) gives:
nun) . 211(1) - 11(1-111.) + 11"(1) (111;)2 + 0(At)‘ (11.)

where x"(t) is the acceleration of the particle given by

the force on the particle divided by its mass. From
eqn. (14), we see that both the current and previous values
are needed, and therefore the first step, x(At), must be
determined before the Verlet algorithm can be started.

First note that since we do not have velocity dependent

85

forces, the dynamical equations are invariant under time
reversal. Secondly, we start the particles from rest.
Therefore x(-At) = x(At) and from eqns. (13a-b), we see that

all the odd order terms vanish. Thus:
11(111) = 11(0) + §x"(0) (At). + 0(At)‘ (15)

which is of the same accuracy as the Verlet expression,
eqn. (14).

In contrast to our dipolar simulations which were

purely long-ranged, we must include short-ranged
nearest-neighbor forces in the present calculations. We
model these interactions with a Born model which to

second-order is described by (Th 76):

- _ _ * . 1
V11 - 2(0 3) Eu“ r11] +2B(1j) (15)

where 31’ 3 61 - 3’, 211 is the unit vector along the
direction of the bond, and a and 3 are the bond-stretching
and angle-bending force constants, respectively. The first
term gives resistance only along the direction of the bond
and for that reason is called the central-force term. By
itself, the four-coordinated diasond lattice would be

unstable. The second, isotropic, term provides the extra

constraints needed for stability. (B.B.: There are problems

86

with this model that arise because of the non-rotational

invariance of the second term. This is mainly a problem for
acoustic modes and those optic modes involving rotations of
molecular units. For our purposes here, these are not of

concern where the ease of implementation is of greater

importance. This subject will be dealt with in much greater
detail in the second part of this thesis on the lattice

dynamics of layered silicates.)

qummzwum

Our supercells consist of a 5x5x5 a 125 cubic array of
conventional cubic unit cells each containing eight sites
arranged in a diamond lattice. This gives a total of 1000
sites, 500 on the A (FCC) sub-lattice and 500 on the B
sub-lattice. We start off with all the sites occupied by
neutral particles which represent the germanium atoms that
we replace with dipolar dimers, representing, for the
purpose of this section, the GaAs and consisting of a
positively and negatively charged particle separated by a
distance equal to the bond length of the diamond lattice.
In contrast to the case of monomer percolation, the
placement procedure is a two step process. First, a site on
the A sub-lattice, chosen at random, must be available, that
is, occupied by a Ge atom. If this is satisfied, one of the

four orientations, defined by eqn. (12), and pointing to a

87

site on the B sub-lattice, is chosen randomly. If the
corresponding B-site is free, that is, occupied by a Ge,
then the dimer can be ”parked" on the lattice with the
positive (Ga) atom always on the A sub-lattice and the
negative (As) atom on the B sub-lattice. If this B-site is
filled, a new A-site is choosen and the procedure repeated.
In this way, a concentration of up to 802 GaAs can be
achieved. Although the mean static polarization, averaged
over many samples, is zero, for any given sample there is a
significant probability that such a polarization does exist.
However, if we include only the dynamical polarization, as

given by eqn. (4), in the force calculations (see eqn. (3)),
then this spurious, finite-size effect, is of no
consequence.

With these rules established we can calculate the
percolation statistics. We do so on diamond lattices of two
different sizes, one with 63 x 8 = 1728 sites, the other
with 73 x 8 = 2744, to test for convergence. Both are
larger than the simulation lattices in order to obtain
better statistics. We scan the ”parameter” space from
x = 0.5 to 0.8, averaging over 100 configurations for each
point. Once the mixed crystal is built, we enumerate the
clusters of GaAs. To do this, we randomly pick a Ga from
the A sub-lattice, then find its As neighbors. In turn,
each of these is then checked for its Ga neighbor, making

sure not to ”retreat” to the Ga already found. This

88

procedure is repeated until the only neighbors left are Ge
atoms, which define the edges of the cluster. We then
choose another Ga, not on any previously determined cluster,
and repeat the operation until all the Ga and As have been
acounted for.

We next determine whether the largest cluster spans the
supercell, which is the finite-size criterion for
percolation. This spanning can occur along any one of the
three directions and in any combination. We choose, in
sequence, a plane perpendicular to a given axis and
determine whether the cluster has a member within that
plane. If it does not, then the cluster does not also does
not span the supercell and hence, does not percolate. If a
cluster does have a member in the plane and every other in
the supercell (there are four such planes for each
conventional unit cell), then the cluster does percolate
along that direction and is given the probability of
percolation of 3. We check the other two coordinate axes in
the same manner. Thus the largest cluster can have a
percolation probability of zero to one in increments of
one-third. This procedure is performed for a given x, on
100 configurations, and the average taken to give the
percolation probability for that value of x. For an
infinite cluster, the (second-order) percolation transition
would have a sharp step function profile. Because of the

finite-size of our samples, the transition is smoothed out.

89

We determine xc, the concentration at which percolation
occurs, using the midpoint of the probability of percolation
which we show in Figure 17. From the tabulated data we
determine the critical concentration for percolation of
dimers on a diamond lattice to be xc = 0.69. The other
curve, gives P(x), the probability of a Ga or As atom being
on the infinite (largest) cluster and is calculated by
simply dividing the size of the largest cluster by the total
number of Ga and As atoms at a given x. It should be
rigorously zero for x > xc, however because of finite-size
effects, it tails off at the percolation point. In
Figure 18 we show a two-dimensional analogue of our mixed
crystal on a hexagonal net. The open symbols represent the
percolating phase and the closed symbols, the isolated
anti-phase clusters. The bare vertices are the neutral
germanium which separate the two phases. The squares
represent the cations and the hexagons, the anions. This
whole unit cell is charge neutral and periodic boundary
conditions hold both vertically and horizontally. (For
clarity, we have not drawn the germanium bonds. It must
stressed, however, that the structures we work with are

fully connected.)

1.00

90

 

 

 

 

 

I I T I r
\.
\.
\.
\
\
0.75 - \
8
L 0 50 -
C1
—- per: (1728)
—- - lnfc (1728) ‘
—-— per: (27“) \\
0.25 _ — — in": (2744) \\
\\
l l l l 1 I L
8.50 L 0.60 0.

Figure 17 Percolation Statistics for Dimers

 

Diamond

91

3
”a
H
”\a

C :1 D: 9 II D:
c 9 0 n o
1: a 1. :1
'3 9 L‘ :0 1: :1
: ‘° '3 .0 1:
’3 *9 0 c :1

 

fl

(0
We 7 (O
‘ )3

Figure 18 (III-V) IV "Hixed Crystal" on Hexagonal Lattice

92

LWWW
Since gallium, germanium, and arsenic are consecutive
elements in the periodic table with masses, 69.7 amu,

72.6 amu, and 74.9 amu, respectively, for the model
(GaAs) Ge mixed crystal we give the same (unit) mass to all
the particles. The disorder therefore lies in the
arrangement of the charges and the resulting electrostatic
forces. In the units of our simulation we choose a s 1 and
fl = 0.4. The unit of frequency is the maximum frequency of
the uncharged diamond lattice given by
(so = [g (a + 2(3) 1"2 = 2.191, where the triply degenerage
eigenmode corresponds to the two-sublattices "beating"
against one another along any one of the three coordinate
axes. The plasma frequency of the pure GaAs crystal, which
is closely related to the splitting (see eqns. (42a-b-III)),
is given by (I): = (Lg- H qz. There are eight particles in a
conventional unit cell of side length a =4/i3 since we
choose the nearest-neighbor distance to be the unit
distance. (Recall that the displacement vector between the
two FCC lattices making up the diamond lattice is given by
d I (a/4) x (§ + ; + £).) We take |q| = 0.3, thus we get a
numerical value for w: of 0.735 = 0.153 mi. We set ‘a)= 1.
The GaAs on the isolated clusters can have the same

phase as the infinite cluster, that is, Ga on A and As on B,

or the anti-phase, Ga on B and As on A. Since there is no

93

natural preference for one or the other, we feel the most
natural arrangement is half and half, therefore we randomly
reverse isolated clusters until 50% of the GaAs, on the
isolated clusters, is of the anti-phase form.

Only one run is neccesary for each simulation as the

results vary only slightly for different configurations of
the same concentration. We used a time step of
one-twentieth of a period defined by T a ZR/wb and each
simulation ran for sixty such periods which gives a
resolution comparable to the experimental value of
approximately 6 cm.‘ (full width at half-height). The run
times ranged between thirty minutes (IBM 3090 CPU) for the
pure (III-V) system where the Coulomb force (electric field)
had to be evaluated for every particle to about three minute
for the germanium where there are no long-range

interactions.

In Figure 19 we plot the Raman response, both
transverse and longitudinal, as a function of w / as

starting at the bottom with pure GaAs and proceeding, with
decreasing ionicity, until we reach pure germanium at the
top. Note that at all values of x (percent Ge), the
response is given by a single peak. As x increases, that
is, as the GaAs is replaced by Ge, the peak position shifts
upwards for the transverse response (dotted line) and
downward for the longitudinal (solid line) until at x=0.7,

when the GaAs no longer percolates, the two essentially

94

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.8 0.9 1.0

Figure 19 Raman Response for (GaAs) Ge

. 1

95

merge. This is shown more dramatically in Figure 20 where
we plot the peak positions (0) as a function of x, with the
solid line to guide the eye. The frequency resolution in
our units is about 0.01.

We can give greater support to the premise that it is
mis-registration of the cation-anion sites that leads to the
loss of the splitting if we now change the character of the
isolated clusters by making all of them with the same phase
as the percolation cluster. We see in Figure 20 that the
response (+) is "repeled" compared to the previous (50-50)
response. That is, the transverse frequency is somewhat
lower and the longitudinal somewhat higher. Note that for a
given concentration, the isolated cluster positions are
identical for all three variations as they were constructed
using the same random number seed.

Further support can be garnered from the dielectric or
infra-red response shown in Figure 21. Unlike the Raman
response where the total weight remains constant for all x
(due to the fact that all bond types have equal
polarizability), in the IR, there is a decrease in the total
response proportional to w;(x) = (1-x) (wile, a result which
follows from a sum rule derived by Thorpe & de Leeuw (See
eqn. (62) in Th 86.) and due to the decreasing number of
charged species. In addition, there is a small transfer of
response to the lower frequencies, somewhat more pronounced

in the transverse response. Up to x=0.6, however, which is

96

 

 

 

 

T I l r I r I T l
1.04 - ~
1.02 - -
O
3
s " .
+
1. +
1.00 - -4
1
0.98 *- 4
L l i l 1 l l l l
0 0.2 0.4 0.6 0.8 1.0

Figure 20 Peak Positions

X

of Raman Response for (GaAs) Ge

97

 

 

 

 

_4

 

0.8

Figure 21

 

 

 

 

\

 

1.
(1)/(1)0

Infrared Response for (GaAs) Ge

1

98

the last point for which a clear peak exists for both
responses, the splitting is virtually the same as that of

the Raman data.

El.Ixnzflads.fisha1121.Llnfihlzfis.

If we move one row down the periodic table for both the
third and fifth column, we arrive at InSb. Indium has a
mass of 114.82 amu and antimony, 121.75 amu. We represent
these in the model with particles of equal mass given by E .

The ratio of the ”bare" (uncharged) frequency of InSb, given

by (ow. = ((201): + (op/31‘”, to the maximum frequency of Ge

is about 185 cn"/ 302 on" = 1.63. (These and other

optical data were taken from the Landolt-Bérnstein Handbook
La 82). To achieve approximately this same ratio in our
simulation we chose the In-Sb force constants to be 63% of
the strengths quoted in the previous section for the
(GaAs) Ge system. For the indium-germanium and
antimony-germanium bonds we use a common convention for
mixed crystal computations and set the strengths to be
(0.63)"2 of the Ge-Ge bonds. There are, of course, no
In-In or Sb-Sb bonds.

We plot the Raman response in Figure 22 using the same
convention as before, solid lines for longitudinal, dotted
for transverse. Each response has (at least) two peaks,

thus the term, two-mode behavior. One the right we have the

99

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

l l l l
x=1.0 L.—
x=0.9
x=0.8
x=0.7
x=0.6A
7 - x=0.54,l\
. x=0.4
m x=0.3
P N Mgo'z
/' Q x=0.l
1L 11 111'; (”a”)
0.5 0.8 0.7 0.8 0.9 1.0 1.1
(1)/u)o
Figure 22 Raman Response for (InSb) Ge

100

peak(s) due to the germanium which is essentially the same
for both the longitudinal and transverse responses. It is a
single peak from x=1.0 to about x=0.5 after which multiple
structures form. These are due to the the isolate clusters
of germanium of various sizes and thus different ”normal
mode" frequencies. For x=0.1 it is difficult to generate
the mixed crystal by replacing 90% of the germanium with
(III-V) dimers, therefore we start with a pure InSb crystal
and place at random, ten percent of the sites with germanium
dimers. Unlike the other concentrations, we have no
germanium monomers. This can be considered a slight
weakness in the model but it is not serious as it occurs far
from the transition point at x80.69. However, it is
interesting to note that the two smaller peaks on the right
for x=0.1 can be attributed to the modes of isolated dimers
of germanium, with the taller one due to the two transverse
modes (displacements perpendicular to the bond) and the
smaller, which is roughly half the size, due to the single
longitudinal mode with motion along the bond. We verify
this by calculating the local density of states of a single
Ge dimer in InSb using the 20H, which involves merely
displacing the two Ge atoms with all others in their
equilibrium positions and then following the dynamical
behavior as before. These two dimer frequencies are marked
along the x axis by (T) at w/wb = 0.760 and (L) at

w/wb = 0.824. Clearly, the two germanium peaks for x30 are

lOl

converging to these values in Figure 22.

The TO-LO splitting, on the other hand, is apparent on
the left side of the figure. The responses broaden,
decrease in height and shift to lower frequency as x

increases, while simultaneously, the splitting diminishes.

Raman measurements have been made on a related material,
(GaSb) Ge by HcGlinn et a1. (Hc 88). We reproduce their
data for the longitudinal response in Figure 23. The
similarities between the measurements and the simulations
are striking although our results show greater structure for
the germanium which is perhaps due to our particular model.
(See the last paragraph.) In the next figure, (24 a-b), we
plot our peak postions (a) for both the InSb and Ge, and
theirs (b) for GaSb and Ge. Qualitatively, they are quite
similar, the main difference being that the actually system
(GeSb) Ge undergoes the zincblende to diamond transition at
around x=0.3 whereas for our percolation model, this occurs
at about x=0.7. In addition, because of the ambiquity in
our germanium result for x < 0.5, we do not plot these
positions in (a) but do mark the germanium dimer frequencies
along the y axis. In general though, we feel confident that
we have captured the essentials of the process with our
model and have given support to the argument that it is this
loss of long-range coherence, as defined by an infinite
cluster of (III-V) semiconductor of a given registration,

that leads to the vanishing of the splitting. We plot our

102

 

 

 

 

 

x = 0.95 _/ K _.
31:0.80 A
~—:::FI73—4‘ —4'-*-“’////fi\\\_______-

Ramon Inlensily
I
’.‘
.0
D
U)
4

M

A
‘—

 

 

 

 

 

 

__ _J
M
x : q13 A“
,_ -1
GaSb 1 . g
100 150 200 250 300 350

 

 

Ramon Shifl (cm“)

Figure 23 Raman Heasurement for (GaSb) Ge (Expt.)

Raman Shift (cm'1)

 

 

1) : {1; +11}: 'GaSb"
ZZOIFO ++++ + +

0 4O 60 80
(3290 Concentration (percent)

A T l 1 T T T 7 T T ’
2500L 11:
++ .

280» "C39" + a

+ +
+11++ Jr

ZGOF+ -
. .
240 r Lo -

-

 

100

Figure 24(a) Peak Positions of Raman Response

(GaAs) Ge (Fxpt.)

 

104

 

 

 

 

Figure 24(b)

1.1 r I 1 I I r *T I I
1 "Ge" [3 I]
E] El
E] E]
0.9 r ‘
L _-
7.1
0.7 - ‘
A
A ”InSb"
A
J?- A
+ + A A
‘l' + + ¢ A a A
L l I 1 l l l l
0.50 0.2 o 4 0.5 o a 1 o
X

Peak Positions of Raman Response
(InSb) Ge (Theory)

105

 

I "2,9) I l l

 

x=0.8

 

X=0£7

 

X=0.6

 

x=0.5

x=0.4

 

/\ x=0.3

 

 

‘ x=0 2
A

1 1 x=0.1
/ +_ +

k M

 

' ‘ x=0.0

, l l J 1
0.5 0.6 0.7 0.8 0.9 1.0

 

U) (I)
/ 0

Figure 25 Infrared Response for (InSb) Ge

 

 

106

simulations of the infra-red response in Figure 25. The
results are by now, self-explanatory. These peak postions

once again are identical with those of the Raman data.

107

mmmmma.

In this first part of the thesis, we have considered
the implementation and effects of Coulomb forces in
disordered materials. Because of the subtleties involved
with this long-range interaction, we have presented,
somewhat deliberately, the basic ideas and techniques.
These include, of course, the LO-TO splitting and the Ewald
method. An important result of this part, is the
demonstration, empirically, that for the important
zone-center modes, the exact Ewald method can be replaced by
an approximate truncated interaction technique if care is
taken to include the long-range bulk and surface
contributions from the dynamic polarization field set up by
the particle motion. We have shown this using a pure
dipolar disordered system for the density of states and the
dielectric responses.

Once established, we have applied this technique to a
disordered ionic system based on the class of mixed crystal
semiconductors of the type (III-V) IV. We have shown, with
both Raman and IR response, that the LO-TO splitting
disappears at the percolation concentration for dimers on a
diamond lattice for both one-mode and two-mode behavior
materials. An account of this work will be submitted to
Physical Review Letters. Although simulation techniques

have long been applied to dynamical problems in condensed

108

matter physics, most systems under consideration have
involved only short—ranged bonding forces. It is hoped that
the techniques developed in this thesis will play a role in
encouraging more research and ultimately, greater insight

into the properties of disordered ionic materials where

long-range Coulomb forces are present.

109

Appendix A

mmmwmm

In this appendix we give a derivation of

eqn. (17-III). Consider the function:
.4. _+
H?) = ge‘k (if ’) nil—i") (Al)

where the sum is over the complete set of lattice vectors
and f(?) is an arbitrary function. It is easy to see that
F(?) is periodic, that is, F(?+R') = F(?) for any lattice
vector R'. Thus F(?) can be expressed in a

three-dimensional Fourier series:
+ 13'?
1(1) . gF-o e (12)
where the sum is over all reciprocal lattice vectors G and

, .9
= ~151M- e": ’ 11(3) (113)

where the integral is over the Wigner-Seitz cell and u is

the volume of this cell. Substituting eqn. (Al) into (A3)

110

gives:

. '9 . _
1‘3 = 13§I113 8.111%) (3 1)) Mil-f) (14)
c

43-? = eat-(1L?)

where we have used e since G-R is an

integral multiple of Zn. From the combined effect of the

sum and the integral, we have:

. .‘9 ,
r» = ‘— 11? 9-1124?) ’ “-3) . (15)
a o
All. Space
9 -tr2
We now choose f(r) = e . The integral is easily

evaluated in spherical coordinates giving:

9 z
e-'k+al /‘t . (A6)

The desired result, eqn. (17), then follows from eqns. (A1),
(A2), and (A6). after taking 3 = 0 and making the change in

the summation from G to —-G. If instead we choose:

a a a a z a
((3) = e-(t+s/a )x ‘e-(t+s/b )y e-(t+s/c )2 (A7)

then from eqn. (A5), we have for the case of 380:

111

(I)
. 2 2
h = 1_ I 1111 e-ucxx) e-(t+s/a )11 x
a 0
-O
a) O
I ., H ,. I .. H
-a1 -a>

where the integrals over y and z are of the exact same form
2

as the that over x. Because e Ax is an even function, we

can replace the complex exponential by its real (symmetric)

part, cos(Gxx). The resulting integral is evaluated in

Dwight (Dw 61) as:

It 2 2
[ ] ‘G / 4 (t+s/a )

 

(A9)
t + s/aa
This gives:
ex:
1 1'!
F» (s) = —
° ° (tn/3’)“ (tn/15‘” (tn/o”)MI
z s s z z z
e-Gx1/ 4(t+s/a ) e-Gy,/ 4(t+a/b ) 2-631/ 4(t+s/c ) (A10)

In the limit, a -+ 0, we recover eqn. (A6) for the case 3:0,

For 3 -+ 0, we get the desired result, eqn. (25).

112

Appendix B

mmmmmmmmmm

In this appendix we look at an interesting general
result obtained by Thorpe S de Leeuw (Th 86) that relates
the longitudinal response to the transverse. The derivation
is worked out in Sec. IV of their paper and is a consequence

of the fact that the shape-dependent effects, contained in

the depolarizing factors, are separable from the other terms

of the potential. We merely state the final result.

st(w) s

8(1) 61(0))

 

Taking the imaginary part of both sides we get:

- .. 2 _1_
Im [81(0),] - so Im [85(0)]

Im [51(w)]
a s

 

8
. (82)
°° (Re 113(11):)” + (II 15(11):)“

Since the real part of 11(w) can be determined from the

imaginary part, through a Kramers-Kronig (El 71) given by:

113

a) 1m [stuo' )]

 

11(11')2 , (113)

2”"

Re [stun] = em +

I s s

o (w') + m

then the longitudinal response, which can be determined
separately by simulation, can also be obtained from the
simulation results of the transverse response. We have
verified this for the two values of x = 0.2, far from the
percolation transition , and x = 0.8, when all the GaAs is
in isolated clusters. We present the results in
Figures. (26a-b). Except for the optic peak for the case of
x = 0.2, the agreement is quite good. Barrow peaks give
difficulties because the computed peak heights are sensitive
to the resolution of the simulation and hence not always

accurately reproduced.

114

 

 

 

 

 

 

 

 

o I I l A
0 0.2 0.4 0.6 0.8 1.0 1.2

(1)/(1)o

Figure 26(a) Longitudinal Response EOH vs KR (x=0.2)

 

 

 

 

 

 

 

0 0.2 0.4 0.6 0.0 1.0 1.2
(1)/u)o

Figure 26(b) Longitudinal Response EOH vs KR (x=0.8)

115

mwumm

ll6

mehsim
Clays or layered silicates are a ubiquitous component
of many soils and sediment deposits throughout the world and

hence have long been of interest to soil scientists and

geologists. In addition, the ability of a class of these
materials, the smectite clays, to absorb a wide range of
intercalants, both inorganic and organic, and act as
catalytic agents, has made them of great interest to
chemists (Pi 83). Recently, physicists have begun to study
clays as they are layered materials with highly anisotropic
properties, and provide an additional system, somewhat
complementary to graphite, for the study of intercalation
and two-dimensional phenomena (So 88). Both materials
provide hexagonal ”cages” for the intercalants; however,
whereas each graphite layer is a simple honeycomb lattice
with only two atoms per unit cell, a single clay layer is a
complicated structure with a large unit cell, consisting of
both tetrahedral and octahedral molecular units. This makes
the study of the dynamics of layered silicates, in
particular the intercalant dynamics, more difficult as we
must contend with a much larger number of force constants
and isolate the low-frequency modes of greatest interest
from the much more numerous and less important ”internal”
modes of the clay. There have been previous studies of the

zone center (3.0) modes, relevant to infrared and Raman

117

spectra, by Ishii et al. (Is 67) and Gupta et a1. (Gu 88).
In this paper, we extend these studies to finite wavevector,
but more importantly, we explain the results in terms of
simple dynamical models and thus gain insight into the
intercalant-clay interaction. Specifically, we derive
approximate semi-empirical relations for the relevant
low-frequency zone center modes (IR and Raman) in terms of
the intercalant force constants and show how these models
can describe the essential features of the phonon spectra
both within and perpendicular to the clay layer. Although
this conceptual "simplification" of the clay is the primary
result of our study, we hope that the more detailed results,
in particular the phonon dispersion curves we provide, will
both encourage and guide those making measurements on these
materials.

The model clay used in this study is typical of many
naturally occuring minerals. After presenting the crystal
structure of these layered silicates, we describe the
phenomenological force constant approach used to model the
atomic interactions. We base our calculations on both
bond-stretching and bond-bending potentials, but in contrast
to the study by Gupta et al. (Gu 88), we avoid the use of
tangential two-body bond-bending potentials as these give
spurious and inconsistent results for both the phonon
spectra and the elastic (long-wavelength acoustic) response.

Using the eigenmode information from 3=0 calculations, we'

118

formulate a simple four-unit linear model for the clay (one
for the intercalant, two for the tetrahedral layer and one
for the octahedral layer) and derive the frequency relations
described in the previous paragraph. 'We next present phonon

dispersion curves along the symmetry axis of rotation

(y-axis), within the plane of the clay layer, and along the
reciprocal lattice vector (z-axis) perpendicular to the
xy-plane. As a summation to the section on the crystalline
materials, we present the results for the elastic behavior,
illustrating the highly anisotropic nature of the crystal
within the plane perpendicular to the clay layer and the
nearly isotropic hexagonal nature within this layer. These
results compare quite well with experimental Brillouin
scattering measurements by Vaughan and Guggenheim (Va 86)
and a pure continuum (macroscopic) model giving us
confidence in the soundness of our microscopic approach.

In the last part of this section we present a simple,
heuristic model for the torsional mode frequency of a
ternary alloy system. Using a probabilistic approach we
derive expressions for the frequency in terms of the
effective force constants. We introduce the effects of
layer rigidity using simple rules that account for
variations in local environment. To test these ideas, we
simulate the disorder using the equation-of—motion method,
described previously, on a two-dimensional intercalated

Ragomé lattice, to calculate the Raman response of the

119

alloy. A good agreement with the above phenomenological

model is found.

120

LL-lealatmmu.

In Figure 1 we show the structure of a typical layered
silicate. We define the x-axis to be within the plane in
the direction of the bias (slant). The s-axis is

perpendicular to the plane. Each layer consists of three
parts, the upper and lower of tetrahedral units and the
middle of octahedral units. Such materials are designated
as 2:1 phyllosilicates. (There are also layered silicates
with one tetrahedral and one octahedral section, the
so-called 1:1 silicates.) The tetrahedral sections consist
of AO‘ tetrahedra, connected at the vertices so as to form a
planar Kagomé lattice, with the apex oxygen atoms forming a
coplanar hexagonal lattice above. A hydroxyl group at the
center of each hexagon gives rise to a triangular lattice
providing the bottom faces of AO‘(OH)z octahedra. An

inverted tetrahedral layer is then attached to the top faces

of these octahedrals to complete the clay layer. The
intercalants are positioned between two clay layers
"encaged" by the hexagons of the Kagomé lattice. These

crystals are monoclinic (space group C:h(C2/m)) with the
upper tetrahedral layer shifted 2a/3 relative to the lower
one where a is the 0-0 distance. Although there is freedom
in the stacking of the layers, for the calculations
presented in this paper the layers are placed such that the

Kagomé lattice of the lower tetrahedral layer lies over the

122

Kagomé lattice of the upper tetrahedral layer below. The
intercalant therefore interacts with the clay layers through
twelve equivalent oxygen atoms. The structural properties,
including the bond lengths, lattice vectors and density are
summarized in Table 1. The (t) and (0) refer to the- 0-0
bond distances in the tetrahedral and octahedral layers.
The 0-0 (t) bond distance is computed from the Si-O bond
distance (1.62 A) taken from Liebau (Li 85).

We now discuss the chemical species ' within the
tetrahedral and octahedral units. These clays are
classified in part by their octahedral components. In the
dioctahedral mineral group, one out of every three octahedra
is vacant while in trioctahedral clays, all such sites are
occupied. Typically the tetrahedra contain silicon with an
occasional partial substitution of aluminium. The
octahedral sites contain either magnesium (tri), aluminium
(di), or iron (di), exclusively or in mixture. For a
listing of structural formuli, the reader is refered to
Pinnavaia (Pi 83). In naturally occuring clays, the
intercalant is typically an alkaline earth or alkali ion,
often hydrated. The charge of the intercalant balances an

equal but opposite charge within the clay layer giving rise

123

 

 

 

 

 

 

Table 1 Structural parameters of the model clay.
Bond Distance (in A)
0-0 (t) 2.64 (E a)
0-0 (0) 2/153 = 3.05
Si-O 13/8a = 1.62
Cs-O V3/2a = 3.23
3 2a x
6 273a y
-’ A A
c 2/3a x + 3.99s 2
s
p 3.34 gm/cm

 

 

 

 

124

to a fairly strong ionic component to the binding which
results in stiff, brittle clays such as muscovite (mica).
In the pyrophyllite-talc group, however, there is no
intercalant due to charge neutrality of the clay layer

itself. Clays of this type are very soft and powdery. For

the calculations in this paper we use the structural unit
cell formula Cs((Siz)Os}2[Hg3(OH)z], where the curly
brackets represent the tetrahdral (silicate) layer and the
square brackes, the octahedral (brucite) layer. This is not
a charge neutral compound and does not occur in nature as
such, but by neglecting the Coulomb force and the species
variation within the clay, we can essentially reduce the
number of atoms needed for a lattice dynamics calculation by
half. As we shall see later, for the "soft” clays such as
talc and the smectites of ”intermediate” stiffness, the
correction to the low frequency zone center modes, where the
”long-range” Coulomb effects are greatest, is of a few
percent. However, for the ”hard” mica-like clays, with
large effective charges, the effect can be ocloser to ten
percent. The reason for not including them at the beginning
is that it is difficult to assign accurate effective charges
to the individual atoms. Furthermore, an exact lattice
dynamical calculation with Coulomb interactions would
require the use of the Ewald method which would greatly
complicate the computation for such a large unit cell

material. (Note however that the ”short-range" part of the

125

Coulomb interaction is already included in the
phenomological forces.) Finally we note that a mismatch
between the natural lengths of the tetrahedral and
octahedral layers causes the tetrahedra to rotate a few
degrees, distorting the Kagomé lattice. (The planar layer
goes from a C to a C symmetry.) We neglect this

av 3v

distortion in our present calculations.

126

ULWMMM

In these lattice dynamics calculations we have used a
combination of phenomenological force constant models to
represent the bond-stretching and bond-bending components of
the interatomic interactions. Bond-stretching is

incorporated through a central-force term of the form:

C
II
NI”
9
|_'
'34
l
11

" I — 1151 — 1111]: (1)

CF

_ 1 z a
- E a (ri uij) + 0(u )
+ e
where a is the radial force constant, r1 = R1 + 111 is the
instantaneous position of the ith atom, rij is the unit
vector between the two equilibrium particle sites, and 311

is the relative displacement of the two atoms from their
equilibrium positions. Note that only motion along the
equilibrium bond contributes to ch.

We model the angular bonds with a three-body
bond-bending potential first suggested by Resting (Re 66).
This potential is particularly well suited for describing

covalent bonds with strong directional character. It is

given by:

127

z
[(3 ° 3 ) - Rz cos6 ] (2)

where r is the force constant, R is the equilibrium bond

distance, and 6 is the equilibrium bond angle. Here, the
ith particle is the vertex atom, the other two being the
”leg” atoms. Expanding to second-order in the displacements

we get:

+ 2(1..- '13.) (r. - 3.)]1» 0(u”) (3)
13 1. 1.): 1.1

In our initial calculations, following Gupta et a1., we

used a two-body bond-bending (de Launey) term (dL 56) given

by:

4 + s C + s
u ° u

- s 2 - m
) - ’23 [111-1 ( r11 “11) ] (4)

where 3 is the angular force constant and 61L is the
displacement perpendicular to the equilibrium bond. (The
central-force and de Launey potentials generate the

so-called radial and tangential forces, respectively, and

128

when combined are sometimes described as the Born model.)
Although easier to implement than the Keating potential, it
suffers from a number of problems. First of all, the de
Launey term is not rotationally invariant and the pivoting

of a rigid tetrahedron, such as occurs with the torsional

mode, coats energy not only due to the distortion of the
basal plane but in addition, and incorrectly, due to the
fixed rotation of the internal O-Si-O bonds. To compensate
for this, Gupta et a1. used an unphysical, negative value
for the 0-0 de Launey tern. With these simple two-body
force terms, they were able to reproduce reasonably well the
values of the 3:0 phonon frequencies calculated by Ishii et
al (Is 67) who used the more complicated Urey-Bradley
potentials. when we used this same force model to calculate
phonon dispersion curves (qln) for the full clay, it led to
spurious softening in finite wavevector modes involving
motion of the octahedral layer. This is due to the
instability introduced through the negative force constants.
The lack of rotational invariance of the de Launey potential
has an additional detrimental effect on the sound velocity
results which we will discuss in greater detail in the
section on the elastic behavior of the clay. The Keating
potential described previously, however, is rotational
invariant, and suffers from none of these disadvantages.
For the reasons mentioned in this section, virtually all

angular forces have been modeled using this potential.

129

These include the O-Si-O and the Si-O-Si bonds within the
silicate layer, the O-Hg-O bonds within the brucite
(octahedral) layer, the 0-0-0 bonds in the Kagomé basal
plane and the O-Cs-O intercalant-silicate layer bond. The
one exception is the O-H angular bond responsible for the
librational motion of the hydroxyl unit for which we use the
de Launey potential.

The values for the silicon-oxygen and magnesium-oxygen
force constants are taken from Ishii et :1. One shortcoming
of our model is the absence of dihedral force terms that are
of importance in the accurate calculations of isolated
tetrahedron vibrational frequencies and thus of some
relavance here. However, we are more interested in the
intercalant-oxygen interaction and the (analytic) results we
give are relatively insensitive to the precise values of the
tetrahedral and octahedral bonds. It is only necessary that
these be significantly stronger than the intercalant bonds.
This requirement is well satisfied. (The Si-O
bond-stretching force constant is nearly two orders of
magnitude larger.) The values for the Cs-O and O-Cs-O
intercalant bonds and the 0-0-0 basal plane bond are chosen
to give ”reasonable" results for the low-frequency modes but
we emphasize that it is not the values of the force
constants in themselves that are important but rather the
relationship between them and the frequencies. The hydroxyl

group bonding parameters were determined by fitting the

130

typical librational and vibrational experimental values of
340 cm-‘.and 3675 cm-‘ respectively (Is 67). The force
constants used in the calculations are given in Table 2.
Once the geometry of the crystal has been described and
the values of the force constants given, it is neccesary to
construct and then diagonalize the dynamical matrix, which
relates the force on each atol to the displacements of the
other particles. Both these steps are carried out on the
computer, the latter using a standard IHSL routine, the
former using a program written by one of us that can be used
on any structure whose dynamics are described by 2- and
3-body potentials. The output of this program, in the form
of w.(3) (eigenvalues) and 3.(i) (eigenvectors), gives the
phonon dispersion curves, elastic properties, and density of
states of the model clay. In the following section, we
describe the low frequency zone center phonons and analyze

then with a simple model.

131

 

 

 

Table 2 Phenomenological Force Constant Values
Bond Strength (in 10‘ dyne/cn)
Si-O 35.08
Hg-O 3.20
Cs-O 0.40

H-O 75.50
O-Si-O 3.89
Si-O-Si 1.50
O-Hg-O 0.40
O-Cs-O 0.60
0-0-0 0.02
H-0 0.60 (2-body)

 

 

 

 

132

mmwummmm

menfln

In general the eigenmodes (particle displacements) of
our full clay are quite complicated, given the large number
(22) of atoss involved. However for the low frequency
aodes, our lattice dynamics calculations show that to a
first approximation, we can group the atoms into four
structural units, indeed, the same four units described in
the crystal structure section. These are the intercalant,
the two silicate layers, and the middle brucite (octahedral)
layer. Coupled together with effective force constants,
this simplified model can be used to derive analytic
expressions for the lowest infrared and Raman mode
frequencies. The (approximate) displacement patterns are
shown in Figure 2(a-d). The masses of the intercalant,
silicate layer, and brucite layer are m (133 amu for Cs),

I, (136 amu), and m. (107 amu), respectively. The effective

force constants are designated by k and k

I 3, the first

giving the coupling between the intercalant and the silicate

layer and determined directly from the structural geometry

and the actual intercalant force constants, “Esq: and
7M“), while the second gives the coupling between the

silicate and brucite layers and is determined empirically

(although it is related to au94Y) Because the effective

 

 

 

133

 

 

 

 

 

 

 

 

 

 

 

Figure 2

Low Frequency Zone Center Hodes

/ \ u / \
e e e o
my) /\ s *— m\ /\ / ‘
I \ / u kI \ /
I
/\ /\ /\ /\ l
[I \v/ \ . / \v/ \
(xy) (2)
. / \ “g k. / \
a
/e\/e\ ‘ ‘— /e\/e\ +
(xy) (3)
x \ / ux kx \ /
(fly) 0 I (s) O I
I / \ k! / \
(a) ir (planar) (b) ir (z axis)
/ \
us
\\\f///\\\f///
.———.
ks
\ /
/\ /\
/ V
k. / \ u
s
/e\/e\ '
kx \ : /
kI / \
(c) Raman (planar) (d) Torsional (Ra)

(S)

(B)

(S)

(I)

 

134

couplings are different for the planar and z-axis motion, we
further distinguish these by the superscripts.

We now show how the effective force constants are
related to the bond strengths (Table 2), considering first

the central-force term. Let 3 8 'J.‘ - 35 be the relative

A

displacement of the two particles defining the bond and r
the unit vector along the bond direction. The force on the
1th particle is given by: F = -a (; ' 3) E. In the
approximate modes we consider, the net restoring force is
along the displacement vector, therefore we need only
consider the force component along this vector. Because in
all cases the displacement is along one of the coordinate
axes, we have: Ft = -ar:u, where i = (x.y.z). Therefore the
contribution per bond to the effective central-force
coupling is given by kcIF = a rt. For the Resting potential,
for the case where 3‘ = 3. = 3, the force on the vertex atom
is given by: FK= -y I; - 3] E, where E a ;s+ ;.. For the

angular bonds we consider, (r‘)x = -(r.)x, (rA)y = -(r.)y,

and (rs): = (r.)'. Therefore, (Fx)x = (FR)y = 0 and

(FR): = -47 r: u, hence, h“ = 47 r: 6“. where 61: is the
Kronecker delta. The intercalant is bonded to a given
silicate layer through six central-force and three angular
bonds. The brucite layer, although distributed over seven
atoms is also considered as a single unit and also coupled

to the silicate layer, through the six Hg-O (apex) bonds.

In both cases the six unit vectors defining the bonds are

135

given by:

3‘ = [:E, o, E] (5.)
32 = [:E, :12,

Thus we get for the total effective force constants:

II“

1 (5...

km" = 2 a (5.)
I “-0

k‘“ - 2 a 4 (6b)
1 " as-o ro-cs-o

(xy)_

km - 2 (am-8 (6°)
(2)

k. = 2 a._. . (6d)

We see that the coupling in the x- and y-directions are
equal and therefore the corresponding modes should be nearly
degenerate. In addition, we find that the z-axis coupling
is further stiffened by the action of the intercalant
angular bond and hence this mode should be greater than the

planar ones. From our force constant values we get:
(z)

¢ xy)

k1

- 0.8 x 10‘ dyne/cm and k = 3.2 x 10‘ dyne/cm.
we now consider these modes in greater detail.
Figure 2a shows schematically the planar mode in which the
intercalant beats transversely against the clay layer. If

we consider the intercalant to have a positive charge and

136

the clay layer to have a compensating negative charge, then
we see that this motion sets up an oscillating dipole moment
within the xy-plane and therefore couples to electromagnetic
radiation polarized with wavevector along the s-sxis.

Figure 2b shows the corresponding beating mode with the
dipole noment along the z-axis. For both of these cases,

the equations of motion for the intercalant and the brucite

layer give:

< xy.2)

z
'xwxur - 2kI (u!+us) (78)
_ (xy.z) _
m.wbu. - 2k. (u. us) (7b)
where “1’ “3’ u. > 0 and «3 = w; 5 m. We define: us E

n3 “1 and u. a n. “1' With these substitutions, dividing

eqn. (7a) by eqn. (7b) gives:

( xy.m)

 

7, [11].; ("I-"'5’ m
m m (ny)
r k! (1 + n8)
Because there is no center-of-mass motion for this a = 0
mode, we have:
2 m n + n n = m (9)

Solving for n. and substituting into eqn. (8) gives the

following quadratic equation for n.:

137

26-3-1 n: + {[263] - 1]- (ti-:1 + 11“} .,
. {(2.33 - 1} = o (10)

c xy,z)

/ k!

( xy,z)

where K E (k
m

). He deternine h. from the
planar mode (Figure 2c) in which the silicate layers beat,
in oppostion, against the stationary intercalant and brucite
layer. Because no dipole moment is generated this mode is
not IR active, but because the polarizability of the
electrons within the silicate-brucite-silicate layer are
different at the 90' and 270' phases of the motion, this
mode is Raman active. From the equation of motion of the

silicate layer for this mode we get:

3 = m who I . (11)
The (average) frequency for this mode is 88.0 cm",
therefore we get k:‘y’ = 5.40 x 10‘ dyne/cm, implying that
ch. 8 2.70 x 10‘ dyne/cm. This value is reasonable given

that dug-o =- 3.2 x 10‘ dyne/cm. lots that this is also the

(2)

value of k-

Returning to eqn. (10), we have:

KW = 6.75 and Kw = 1.69, which gives (from the
quadratic formula using the negative branch): 77:”) = 0.33

and n;”'= 0.25. From eqn. (7a) we have:

138

1/2
2k!
(DIR = ;— (1+Y)8) (12)
I
thus (0:1)" = 52.1 cm_1 and not; = 101.0 cm”. These should

be compared with the actual values from the full lattice

dynamics calculation: 0);?" = 51.2 cm.1 and

(2)

I. = 99.6 cm-1. ‘The relative errors are 1.76% and 1.41%,

respectively. Considering the simplicity of the model this
agreement is very good. More importantly, eqns. (6) and
(10)-(12) can be used to extract the values for the
phenomenological intercalant force constants, 0110-0 and
ro-rc-o (and as-s)’ from experimental measurements of the
low frequency zone center modes. If we simplify the model
further by considering the silicate-brucite-silicate layer

to love as a single unit, we have:

(xy,z>

z
mxcoxuI - 2kI (1+7?) uI (13)
where n = 'r/(H-'1) = 1.35 (M = 512 amu). This gives
0);” s 52.5 cm“ and not; = 105.0 cm“. The first agrees

quite well with the previous value while the second is quite
a bit different. We infer from this that modes in which the
particle motion is perpendicular to the plane probe to a
greater degree, the internal structure of the clay layer,

which is certainly reasonable.

139

We now present a derivation of the frequency of the
torsional mode. This node is unique in that its frequency
can be calculated exactly due to the simplicity of the
eigenmode, shown in Figure 2d. We see that only the oxygen
atoms in the silicate basal plane love, and only within the
plane. Concentrating on the coupling of the oxygen atoms to
the (fixed) intercalant, we find there is no contribution
from the intercalant-basal plane angular force and a
contribution to the total effective force constant from the
intercalant-oxygen stretching bonds of §~a . For the

Cm-O

planar 0-0-0 angular bonds of the Kagoné lattice we have:

:.=[:-.::.01

was»)
[we] - ransom“ as]
35 [[w] - ra-sonu- es)

where u is the magnitude of the displacement of a given
oxygen atom. To find the contribution from these bonds we

first calculate the potential energy of this distortion from

- 1 3
eqn. (3). We have. vo-o—o - 2 [92' ] u . We have two

of these bonds per basal plane oxygen, thus the effective

140

torsional mode force constant is given by:

4
k - 5- ace-o + 18 y -o (15)
and
2 (TOP)
cor" - k / “Oxy . (16)

Note that from eqn. (15), we see that the basal plane
angular force is much more effective than the stretching
force in stiffening this mode. In addition, if “c.43 has
been previously determined, then eqns. (15)-(16) will give
the value of ro-o-o‘

Finally we address the problem of the the Coulomb
forces. The "long-range" portion of this interaction gives
rise to an additional force on the particles due to the
macroscopic electric field generated by the particle
displacements. We can estimate this effect quantitatively
by noting that the shift in frequency (squared) is of the
order of the plasma frequency (squared), wi, given by
(4n/0) Z (qt/mi), where 0 is the unit cell volume (equal to
2.54 x 10.22 cm“) and the sum is over the atom in the cell.
For our purposes, we consider the clay to consist of two
dynamical units, the intercalant with a mass of 133 amu
(Cs)and the remainder (silicate-brucite-silicate layers)
with mass 379 amu. Assuming a relatively small charge

transfer‘ between the two units of 0.4 e (e = 4.3 x 10'“’

statcoulombs), we get Aw z:w: / 2w equal to approximately

141

1.5 cm-‘ for a) = 75.0 cm”, a typical frequency of interest.
For muscovite, where the effective charge is as large as
1.0 e, we get a shift of about 9.5 on" for the same base
frequency. For such materials it is probably advisable to

include the Coulomb interactions. Note however that not all

modes will be affected. For example, the Raman active
frequencies do not create a dipole moment and thus do not
produce an electric field to couple to the charges. The
results of this section could still be applicable to these
materials if it were possible to subtract the contribution
of the long-range forces so that we could work with only the
”bare” frequencies. This is relatively simple for cubic (or

isotropic) structures (once the effective charges are known)

where, for the zone center modes, we have shown,
a): = w: - 3 w: and a): I to: + go): where T and I. refer to the
transverse and longitudinal polarizations and ab is the
frequency without the charges. For our monoclinic
structures, however, we are unfamiliar with the

corresponding expressions. We stress though, that overall,
many of the results of this paper are qualitative and quite
robust, with relevance to materials of both small and large

ionicity.

Lzhamnimnimmms.

We now discuss the finite wavevector response of the

layered silicate, shown in Figure 3(a,b). Although our clay

142

has over sixty branches, we will consider only a few,
concentrating entirely in the region below 250 cm”. For a
solid such as this, consisting of many weakly bound
structural units, most of the modes will be relatively

dispersionless and more characteristic of the seperate units

themselves. Our interest is with modes that involve
inter-unit motion, in particular, the motion of the
intercalant against the clay layer. Through the use of

simple models, including the ones introduced in the previous
section, we can explain many of the prominent features that
occur in the phonon dispersion curves.

We first consider the phonons with propagation vector
perpendicular to the clay layer, along the z-axis. One
striking feature, that is characteristic of phonons in
layered materials, with wavevectors perpendicular to the
planes (Be 88), is the large number of virtually
dispersionless modes. But in addition, there are two groups
of curves, one in two sections, the other in four, that are
characteristic of the zone-folded acoustic mode of a linear
chain. Specifically, these features are well described by
the two-unit (intercalant - clay) and four-unit (intercalant
- silicate - brucite - silicate) chain models used
previously. These two models relate to the transverse and
longitudinal modes of the full clay respectively. The
dispersion curves from these models, using the same

effective force constants determined previously from the

143

zone-center analysis, are plotted on Figure Be, as dashed
lines for the first model and as dash-dotted lines for the
second. We see that the acoustic modes are very well
reproduced, the lower optic modes moderately reproduced, and
the higher optic modes only qualitatively reproduced.
However, in general, the approximation of grouping the
individual atoms into larger dynamical units remains useful

even when extended to finite wavevector.

We next consider phonons with wavevector within the
clay layer along the two-fold rotation axis. The dispersion
curves as shown in Figure 3b demonstrate a quite different
behavior from those along the z-axis, as one might expect
for an anisotropic material. In the previous case, the
phonons probe the interaction between the layers, ignoring
for the most part, the structure within a given layer. In
this case, however, the phonons are greatly influenced by
the interaction of the individual intercalants with the
hexagonal "cage” of the clay, and therefore the detailed
structure of the intercalant - silicate basal plane boundary
must be adequately modeled. Indeed, to first approximation,
we can concentrate exclusively on the dynamics of a single
interface. The atoms of the clay layer are strongly
connected to each other as well as being coupled indirectly
through their weaker bonds to the intercalant which are
placed in the ”interstitial” regions. The model structure

having these features is shown in Figure 4. The values of

(1) (cm“)

250

200

150

100

144

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 3(a)

.10 0.15 0.20

qa/ 1t

w(k) for k = (0,0,k)

0.25

(0 (cm")

145

 

250

 

 

 

 

 

 

 

200

 

 

150 - ’7/

 

 

 

 

 

qa/ K

Figure 3(b) w(k) for k = (O,k,0)

146

the stretching and angular clay-intercalant bonds are chosen
to give the infrared zone-center frequencies calculated
using the two-unit model (eqn. 13). The ”atoms” within the
clay layer are coupled via central-force and de Launey terms

'(tbe latter is neccesary for transverse rigidity) with

strengths equal to the 81-0 bonds. Note that this model
allows for both transverse and longitudinal modes. The
resultant dispersion curves are overlayed on Figure 3b with
the transverse modes displayed in dashed lines and the
longitudinal in the dash-dotted pattern. The characteristic
features of the acoustic and optic branches of this model
are clearly present in the full clay. The quantitative
agreement with both the transverse and longitudinal optic
modes and the longitudinal acoustic mode is good. The
characteristic optic and acoustic branches can be explained
as a hybridization between the acoustic mode of a monoatomic

linear chain made up of the larger ”clay” atoms and the

dispersionless optic mode characteristic of isolated
intercalants vibrating between them. The interaction
between these two modes that gives rise to the

aforementioned hybridization leads to the greatly flattened
acoustic mode and the strongly dispersive, sigmoidal shaped,
optic mode. We can conclude that one of the characteristic
structural features of the clay, namely, the loose ”traping”
of the intercalant in the hexagonal holes of the basal plane

has a definite effect on the underlying structure of the

147

\ /

lntercalants

Figure 4 Structure for y-axis Dispersion Hodel

148

phonon dispersion curves for in-plane propagation vectors.
The one shortcoming of this model is its failure to properly
describe the transverse acoustic modes of the full clay.
Obviously, it is neccesary to treat the transverse rigidity
of the clay layer with more subtlety than this simple ”ball
and spring” model allows. We consider these dynamics, in
particular the softening of the transverse acoustic mode
with particle displacement perpendicular to the layer (TAL),
in greater detail in the second part of the next section.

To give a general overview of the lattice dynamics we
plot the density of states for all modes less than 1000 cm'4
in Figure 5. This was calculated by the usual technique of
calculating the frequencies for a given k chosen at random
in the Brillioun zone and creating a histogram. We took a
sample of 5000 points in reciprocal space which gives
approximately 100,000 total frequencies. We neglect in this
graph the single peak at 3676 cm-‘ due to the O-H stretching
mode of the hydroxyl. It shows a number of sharp peaks as a
consequence of the relatively dispersionless intra-molecular
modes. Just below 1000 cm-‘:Ls the peak due to the Si-O

bond stretching of the tetrahedron. At about 330 cm-‘

is
the O-H librational mode which has a strong neutron signal
due to the presence of the hydrogen. Heasurements by
Reumann et a1. (Re 89b) show this structure is sensitive to

the presence (or lack of) vacancies in the octahedral layer

giving a two-peak structure for dioctahedral phylosilicates

149

such as montmorillonite and a single-peak for the
trioctahedral type materials (this model). At present, the
simple de-Launey force model used for the O-H bond is too
crude to duplicate this behavior. It appears possible that

the double peak is a "mixed-crystal" effect as the number of

vacancies might vary from cell to cell and give rise to
different local environments for the hydroxyl group and
hence slightly different frequencies. Finally note that
there is a strong peak at about 50 cad due to the flat

longitudinal acoustic mode, essentially due to the

intercalant "rattling" within the hexagonal cage.

150

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1 1 1 .J_ 1 L 1 l
o 200 400 500 800 1000
(0 (cm")

Figure 5 Density of States of Model Layered Silicate

151

LWW

Elasticity theory is concerned with the long wavelength
acoustic dynamics of a solid. In this regime, A >>a where a
is the typical atomic separation and A the wavelength of the
phonon, and all atoms are moving essentially in phase.
Therefore the relative distortion of the atoms within a
given unit cell is small and as a result, we can ignore the
discrete nature of the material, describing the distortion
from equilibrium by the "average” or center-of-mass motion
of a macroscopically small though microscopically large
region. If we assume a linear relation between the
distortion (strain) and the resulting restoring force
(stress), the harmonic normal modes of the system are given

by the solutions of the secular determinant (Lo 34):

A A 2
' cijkl. kjkk " 9" 65,1, I = 0 (17)

(ijkl = mm)

A

where ci is the elasticity tensor, k5 is a component of

jkl
the unit vector along the propagation vector, p is the mass
density, and v is the velocity of propagation. Because each
indice can take on three different values, the "dynamical"
matrix in this case is 3 x 3 with three eigenvectors (and

values) corresponding to the three acoustic branches.

Although the tensor is of fourth-rank and has 81 components,

152

the number of independent elements is twenty-one for the
most anisotropic crystal and only thirteen for the
monoclinic clay structure. Furthermore it is possible to
group pairs of indices and identify these by a single
number. We use the convention: xx +4 1, yy +4 2, zz +4 3,
yz +4 4, zx +4 5, xy +4 6. In this form, the components of
the elasticity tensor are called the elastic moduli. In
principle, these moduli can be related analytically to the
force constants (see Born and Huang (Bo 54)), although this
procedure is difficult to carry out for all but the simplest
lattices. Alternately, they can be determined
experimentally through ultrasonic or Brillouin measurements
(Va 86). We, however, can determine these velocities from
the slopes of the acoustic modes, generated by our full
lattice dynamical calculations.

This data, for our clay, is presented in Figure 6a, for
propogation vectors both within and perpendicular to the
layer. Note the nearly isotropic behavior within the layer
(x to y) and the highly anisotropic response between layers
(z to -z). As a result of monoclinic (hexagonal) symmetry,
one can show that the two lowest velocities along the z-axis
must (nearly) correspond with those from the lowest curve
along the x- and y-axes. It was the failure of our previous
models with de Launey angular forces to satisfy these
symmetry requirements, that alerted us to the problems with

these potentials. In fact, the Keating potential was first

153

introduced precisely to correct this very problem with the
de Launey (Born) dynamical models, although it was the
simpler diamond structure that was initially considered
(Re 66).

In Figure 6b, we reproduce the acoustic velocity data
from Vaughn and Guggenheim for muscovite, a common mica.
(N.B.: The y-axis in Figure 6b is incorrectly labeled. The
velocities are actually in km/s.) The symbols are the
measurements (from Brillion scattering) and the solid curves
are the calculated velocities from an elastic constant
model. That is, the elastic moduli are determined from a
subset of the experimental measurements and these moduli
then inserted in eqn. (17) to determine the velocities
(eigenvalues) for any given k. The fact that our curves,
computed from the full dynamical matrix, are quite similar,
gives us confidence that our sicroscopic model is capable of
correctly describing the dynamics of these layered
silicates. For our structure with no torsional distortion,
the simpler hexagonal symmetry model (with only 5
independent constants) would have been sufficient, but in
general this will not be the case. Indeed, the experimental
data shows greater anisotropy within the plane for a real
material where this distortion is present.

We now briefly sketch how the elastic moduli (for. the

hexagonal case) are determined from the propagation

velocities. The value for cum: cu) (in gigapascals) is

 

154

 

 

 

 

 

 

 

 

 

 

 

 

 

 

7
L
a s —
E
§.
> i
3 _ '
i
I . I
: 1 '
1 1 l l
x y z x -2
Figure 6(a) Sound Velocity as calculated from Model
.0 1 ‘ 1 D 1 1 1 C. 1 1 L I “C.
0?: o m o Mnntumu M0000 " no |
=1 ; \ é ' ‘~ *
I \\.. j IO// U
. 2' ~ 2 {l \
3:1] \‘1 i .. 2/ a
= 7 \ "
g . “%
> .
3 . 3
.f E Tunsverie Modes I
-0 '0 1° 00 IO ’0 3 =0 )0 0| CO ’9 0 '0 80 4O 00 ’0 00 'Ol '10 '30 10° '00 'IO
Acausuc Propagauon Direcuon (Dogma from Ann
Figure 6(b) Sound Velocity for Muscovite as measured

by Brillioun scattering

155

where p is the density in kg/ma and

iven b 2
8 Y 9" Lazy) V Lucy)

is the velocity of longitudinal waves within the xy-plane in

2
meters per second. Likewise, c” - pvu” where vb“) is
the velocity of longitudinal waves along the z-axis. From

the transverse waves with wavevector within the plane, we

2 2

get c“ - pvflxy) and c“<z 0.55) - pvTm where vnxy) and
vTm are the velocities polarized within and perpendicular
to the xy-plane respectively. Finally, c1: can be

determined from either one of the acoustic modes propagation

along the k = V1/2(l,0,l) direction and polarized within the

xz-plane. Specifically, if we define x a pvz,
A '=" (cn+ c“)/2, B E (c”+ o“)/2, and C 2 (our c“)/2,
then C2 = h: - (A+B)x +AB, from which we get cla. The
values of the elastic constants, obtained from our
calculations, are culc22 a: 161.1, c” a: 50.1,
c«/c55 a: 14.4, c“ a: 55.0, and sale” a: 8.4, where the
figures quoted are in gigapascals. Note that these
materials are much stiffer within the plane than

perpendicular to it, as expected.

In the preceeding paragraphs of this section, we
analyzed the long wavelength region of the acoustic modes in
terms of an elasticity theory that assumed a linear relation
between stress and strain. This leads to a second-order
equation of motion and hence a linear relation between
frequency and wavevector. To explain the softening in the

lowest transverse acoustic (TAL) branch (see Figure 3b);

156

requires that this assumption be augmented. For small qa
(x >> a), this mode exhibits an anomalous quartic term, that
is, a? is proportional to A.(qa)z+B (qa)‘, where A and B are
both positive. This is in contrast to the usual case where

B < O, which follows from a small as expansion of the
typical acoustic mode dispersion w a:sin(qa). This
behavior, which occurs to an even greater extent in graphite
(So 88) and is present in the neutron scattering
measurements on muscovite by Wada et a1. (Wa 86), reflects a
distinctive elastic property of layered solids. That is, in
the continuum limit where elasticity theory is applicable,
these materials should be considered as laminar plates
embedded in a three-dimensional matrix. The relevant
equation of motion for these plates is that for transverse
waves on a beam, which is fourth-order in the derivative of
the displacement, leading to a quadratic dispersion relation
and hence quartic term for “F. It is the compression (and
extention) of the layer (within the plane) that gives rise
to the restoring force, thus this term in the dispersion is
related to the Young's modulus (c“,czz). The linear part
of the dispersion is due to the restoring forces
perpendicular to the plane and gives a measure of (c“,c’5).

We can determine the relative importance of each of the
microscopic force constants on the two kinds of elasticity

(linear or quadratic) by calculating the change in the ratio

B/A as each force constant is varied (t 202), one-at-a-time,

157

from its "normal" value found in Table 2. To determine A
and B, we performed a quadratic least-squares fit (with zero
constant term) of a} vs (qa)z along the y-axis from the .zone
center to qa/n = 0.1 at intervals of 0.05. We find that the
variation of the O—Si-O bond-bending term is strongly
correlated with the ratio B/A while the Cs-O stretching term
is strongly anti-correlated, that is B/A decreases as this
term is increased. Furthermore the Si-O, Hg-O, and O-Hg-O
terms are rather weakly correlated and the Si-O-Si, O-Cs-O,
and O-O-O terms virtually uncorrelated with this ratio. It
is not surprising that the Cs—basal layer coupling plays a
greater role in the usual elastic behavior as it acts as the
matrix "glue" between the clay layers. Likewise, since
bending the clay layer has a greater effect on the
tetrahedra, it is the silicate bonds that most strongly
determine the quadratic behavior of the material.

For a final look at this system we briefly discuss the

dynamics of a single Sizo5 tetrahedral layer. In Figure 7

we can see quite clearly the strongly dispersive
longitudinal and transverse optic modes without the
complications due to the coupling octahedral layer. Also

note that the transverse perpendicular acoutic mode shows
very plainly quadratic dispersion. Because we have only a
single layer, the two-dimensional nature is greatly enhanced
and without the addition octahedral and silicate layers and

the z-axis coupling between layers, the linear aspect is

158

absent. To test the idea that the lack of coupling between
layers of finite thickness interacting internally through
harmonic forces leads to exact quadratic dispersion, we
calculate the dispersion of two monatomic chains staggered
and coupled so as to form a triangular truss. .This system,

with two atoms per unit cell and with displacement allowed
only within the plane of the structure, has a 4 x 4
dynamical matrix that can be solved analytically in the
elastic limit. Although the algebra is somewhat tedious and
will not be shown here, if we look for solutions of the
form: w = A (qa), that is, acoustic modes, we find one with
non-zero A corresponding to the standard linear elastic mode
with longitudinal displacements and the other with A=0 for
the (quadratic) mode with transverse displacements. This
supports the premise that we can attribute all the linear
dispersion to inter-layer coupling and the quadratic

dispersion to the infra-layer (bending) dynamics.

0) (cm“)

 

250 1 l x

 

 

 

 

 

 

 

 

l T F ' I ' I r

200 *- .4
150 L 4
100 F -1
50 -

0 1 l l 1 l L 1

0 0.10 0.20 0.30 0.40 0.50
qa/ 7t

Figure 7 (M?) for k. = (0,k,0) for $1205 Layer

160

mewmmwm

memm

In the previous sections we have discussed the dynamics

of crystalline layered silicates in which the unit cell is
periodically repeated. We now consider the case where the
intercalant atom can be one of two different species, placed
at random within the hexagonal cages. The chemical formula
for such a material is given by Arme [Clay], where A and B
are the two intercalant species and the [Clay] designates
the layered silicate, for example, vermiculite.
Experimental measurements of the torsional mode of these
materials show an interesting variation in the frequency as
the concentration x of the intercalant B is varied. In
Figure 8a we show the Raman measurements by York et a1.
(Yo 85) on Cs“x be Vermiculite; Figure 8b gives the
torsional frequency as a function of x. We see quite
clearly the non-linear behavior although also note that this
is a very small effect in this case as the end member
frequencies are within a few percent of each other. In this
section, we present a simple, probabilistic model that can
explain the qualitative features of the (w vs x) curve in
terms of local force constant variation in connection with
the effects of layer rigidity. To test this model we

compute exactly, using the Equation-of-Hotion technique

161

 

 

RAMAN INTENSITY (ARBITRARY UNITS)

 

 

 

l l l

75 87.5 100 1 12.5 125
RAMAN SHIFT (cm")

 

Figure 8(a) Raman Measurements of Torsional Hode

 

107
CsM‘beV

106

Raman Shift (cm‘I)

 

 

 

105 1 1 1 1
(10 02 (14 03 (18 L0

 

Figure 8(b) Torsional Hode Frequency vs x (Expt.)

162

described in the first part of this thesis, the Raman
response of the torsional mode for a mixed crystal Kagome
lattice where two species of intercalants are present rather
than the signle one in the crystalline case. We feel that

this approach gives more physical insight into the situation

than simply fitting the phonon frequencies (Gupta et a1.
Gu 88) by adjusting the intercalant force constants of an
averaged "virtual” crystal defined by taking the mass of the
intercalant to be the average of the two end members. This
study (Gu 88) was further weakened by the use of a deLauney
force model (dL 56) with unphysical negative force

constants.

81.qu

As discussed previously, in the torsional mode, the
tetrahedral units pivot around an axis perpendicular to the
basal plane such that only the oxygen atoms within this
plane move. All other atoms, including those in the
octahedral brucite layer and the intercalant, remain
stationary. This allows us to study the mode and the
effects of alloying, with a simplified system consisting of
a single basal plane layer with one intercalant in the
center of each of the hexagonal opening and within the
plane. We show a representative sample of the model in

Figure 9. Although the oxygen atoms are weakly bound to

163

each other, the effective coupling via the silicate bonds
within the tetrahedrons is quite strong. To maintain the
triangular shape in the basal plane, we make the O-O
(stretching) bonds relatively large with unit strength. The
oxygens are also given unit mass. In contrast, the
intercalant-oxygen bonding is assumed to be weak with
strength 0.040 for the O-A bond and 0.024 for the 0-8 bond.
Since we are modeling the system roughly for Cs-Rb
Vermiculite, we give masses 8.0 to A (Cs) and 5.0 to B (Rb).
We consider only bond-stretching; the O-IC-O and O-O-O bond
bending forces, although important when modeling the
dynamics of the entire clay, are not essential for capturing

the basic physics of the alloying effects.

164

 

 

>

 

Figure 9 Two-Dimensional Kagomé-Intercalated Lattice

 

165

1;. WMWLLMW
From Figure 8a, we see that even with alloying, the
torsional frequency can be defined, to a first

approximation, by a single number giving the peak position.

(It is also characterized by the width of the response,
which we ignore, however.) From the point of view of the
oxygen atoms, in the crystalline limit (x=0 or 1), the
intercalants act as fixed posts between which they oscillate
back and forth along the line defined by bonds. This

frequency is given by:

‘0

 

w = V'2k / m (1)
Tor oxy
where k is the O-IC (stretching) force constant and moxy is

the oxygen mass. Note that the intercalant masses do not
appear as they are stationary, nor does the O-O coupling
since all oxygen atoms move in phase so as not to distort
this bond. We now randomly substitute different intercalant
species. If we consider the oxygen to be independent of one
another, then each would oscillate at a frequency given by
eqn. (1) but with k‘ + k2 substituted for 2k, where k and

I.

k2 are the two intercalant-oxygen force constants, not
neccesarily equal. (We have implicitly assumed that the
intercalants remain essentially stationary even after the

alloying. This is strictly true for the pure cases, x=0

166

and 1, however it remains a good approximation for
intermediate x.) We now assume that the major effect of the
mutual oxygen bonding, related ultimately to the internal
silicate bonding, is simply to average the various force

constants so that the effective frequency of the torsional

mode in the alloy is given by:

 

1.1”” = V2 11(111/ m (2)
Tor o oxy
where k.(x) is the effective force constant at
concentration x, given by the average force constant k. We

can make this approximation because the 0-0 force constant
is much larger that the fluctuations in the
intercalant—oxygen coupling due to the alloying. In the
case where the O-IC force constants retain their pure case

values k and k , we have:
A a

k.(x) = k = (1-x) k‘ + x k. . (3)

This gives rise to a nearly linear relation between a)?“
and x.

We can now introduce the important effects of layer
rigidity via the force constants. In Figure 10a, we show

schematically, an infinitely rigid clay layer stop to large

intercalants, labeled A, with a smaller intercalant B

167

between them. In this case, the equilibrium bond distance
for the O-A bond remains unchanged from its value in the
pure system. However, for the O-B bond, this distance is
increased due to the pillaring effect of the surrounding A
intercalants. The effective coupling is weaker than in the
pure case and thus he should be replaced by ‘1: km where I.
is a reduction factor less than or equal to one. If the
clay layer were completely "floppy", as shown in Figure 10b,
then it (the layer) would conform to the various heights of
the different intercalants in such a way that the couplings
would always be given by their pure case, the so-called
”Vegard's Law" regime. We incorporate these ideas into our
model by adopting the following rules. If a 8 atom has an A
atom as its nearest-neighbor on the (triangular) intercalant
sub-lattice, then the O-B force constant along that bond is
reduced by the factor f.- Since the intercalants are
arranged randomly, the probability of a B atom having an A
atom neighbor is given by (l-x). The corresponding
probability of having a B atom neighbor is of course given
by x. The O-A force constant, at this stage, is unaffected

by its local environment. Thus we have:

k.(x) = k (1-x) k‘ + x [(l-x)#. + x] k.

1‘ {(l-x) + [#.x + (1-£_)ler1 (a)

168

where r a (k./kA) < 1. We plot the corresponding frequency

 

w (x) = 72k (x)/m in Figure 11a for r=0.6 and #?80.5
e e oxy I

(

 

) and r=0.8 and fin=0.25 (- - -). Both show sub-linear
behavior with the second exhibiting a minimum for x < 1.
(We will explain the origin of the symbols in the next
section.) It is easy to see that the value of k.(x) given
by eqn. (4) is always less than that given by eqn. (3). If
a larger A atom is surrounded, totally or in part, by
smaller B atoms, this can lead to a pillaring effect which
produces an enhancement in the O-A force constant due to the

pushing down of the layer on the intercalant, then we get:
k.(x) = k‘ [ (l-x) {(1-x)°1 + xf‘} +

k
f!- x{x-1 + (l-x)fi.}] (5)

A

where we have assumed this effect is the same for all six
oxygen-intercalant bonds about the A atom I and is
proportional to the number of B nearest-neighbors on the
intercalant sub-lattice. The enhancement factor 4‘ is
greater than or equal to one. In order for a polynomial to
show both sub- and super-linear behavior, it must have a
point of inflection, that is, the second derivative must
vanish at some value of x. Since eqn. (5) is only second
order (quadratic) in x this is impossible and thus can give

only pure super- or sub-linear behavior. In Figure llb we

169

”Clay"

0

Intercalants

 

Figure 10(a) Infinitely Rigid Clay Layer

"Clay"

   

lntercalants

Figure 10(b) Floppy Clay Layer

170

choose the parameters r=0.6, £‘=1.6, £.=0.5, to give a
super-linear variation with x. Our basis approach can
continued, considering, for example, the effects of the
next-nearest-neighbors, which gives a cubic equation for

k.(x), with each step leading to a higher-order polynomial

approximation and richer behavior for the w vs. x curve.
The cost is that new parameters are introduced at each
stage. Nonetheless, we can test for the validity of this
parameterization procedure and the effective force constant
approximation, eqn. (2), by comparing our predictions
against the actual behavior of the model system, calculated

using the equation-of-motion (EOH) method.

Lilnlamfisulm

We use a 32 x 32 intercalated Kagomé lattice with 1024
sites consisting of 256 intercalants and 768 oxygens. The
intercalant were placed randomly within the hexagonal sites
for x from zero to one in increments of 0.1 and the force
constants determined using the local environment rules
presented in the previous section. We assume that the
torsional mode induced in the alloy is the same periodic
torsional mode as in the pure system. This is equivalent to
assuming the polarizability of the intercalant bonds is
unchanged by the alloying and the same for both O-A and O-B
bonds. The torsional frequency was determined using the

peak height of the response, an example of which is shown in

171

 

 

 

 

 

 

 

I l T I I I I l I
1'0 b-O.6 r-o.so ‘
—— r-O.8 f-0.25
0.9 -
o
3
\
3
0.8 r
0.7 1 1 n 1 1 L 1 1 1
0 0 2 0 4 0.6 O 8 1.0
X
Torsional Mode Frequency vs x (Theory)

Figure 11(a)

(Sub-linear Behavior)

172

 

1 .

0.

w/w0

0 .

 

 

 

o.711h111111

 

Figure 11(b) Torsional Hode Frequency vs x (Theory)
(Super-linear Behavior)

173

Figure 12. These results are represented as the individual
symbols on Figure 11(a-b) for the same choices of
parameters. The fits are rather good although both the sub-
and super-linear model overestimate the actual frequencies

for x < 0.5. Note that the minimum for the dotted curve in

Figure 11a is rather nicely reproduced. The resolution in
the figures is given by $0.03 frequency units. The simple
assumptions we made in deriving our probabilistic model
appear to be reasonable ones although admittedly this model
provides more of a convenient (but physically motivated)
parameterization of the problem rather than a fundamental
explanation, which could be provided, for example, by
expressions for the adjustment factors I‘ and I. in terms of
the local spacing between layers. It should be emphasized
that these numerical calculations are exact (within the
limits of the model) and make no assumptions about fixed
intercalants or oxygen-oxygen bonds. To demonstrate the
accuracy of this method we plot in Figure 13, the density of
states for a pure, periodic Kagomé-IC lattice using both the
direct diagonalization of the 8 x 8 dynamical matrix sampled
throughout the Brillioun zone and the EOH using a 56 x 60
supercell. All features are quite well reproduced by the

simulation.

174

 

 

 

 

 

P
—
"' X
I
d
O
>
I—
P
h.

 

O4 06 0.8 10 12 14
(1)/(Do
Figure 12 Raman Response for Torsional Hode

calculated with EOH

175

 

 

5 I I I I I I I r I I I
4 - -
Diag
—°-— EOM 3360
3 .-

 

 

 

 

 

 

 

 

 

Figure 13 Density of States for Kagome-IC calculated
via Direct Diagonalization and EOH

176

Him

In this part of the thesis we have surveyed the lattice
dynamics of a typical layered silicate, with particular
emphasis on the intercalant-clay layer dynamics. With the
help of simple models we have developed approximate formulas
for the low frequency zone-center modes in terms of the
intercalant-basal plane force constants providing a
procedure for these quantities to be determined from
experimental data. We have extended these calculations to
finite wavevector and studied the elastic properties of
these materials for the first time from a microscopic
viewpoint with a particular emphasis on the dispersion of
the transverse acoustic mode. In the second part we have
introdued a simple model for the torsional mode frequency
for ternary alloy systems that we have checked with exact
simulation results. In general, much work needs to be done
in the experimental realm and it is hoped that this study
will encourage experimentalists, expecially those involved
in neutron spectroscopy, to make more detailed measurements

of this class of layered solids.

BIBLIOGRAPHY

Ab

Be

80

Ch

dL

dL

Dw

E1

Ga

Ga

Ga

Ga

Gu

Is

75

88

54

71

56

80

61

71

76

83

89a

89b

88

67

177

REFERENCES

R. Alben, D. Weaire, J.E. Smith, Jr., and M.H.
Brodsky, Phys. Rev. B 11, 2271 (1975)

M. Bernasconi et al., Phys. Rev. B 354 12089
(1988).

M. Born and X. Huang, Dynamical Theory of Crystal
Lattices (Oxford University Press, London, 1954).

I.F. Chang and 8.8. Mitra, Adv. Phys. 20” 359
(1971)

J. de Launey, Solid State Physics, Vbl. 2, edited
by F. Seitz and D. Turnbull (Academic Press, New
York, 1956), pp. 219-303.

S.W. de Leeuw, J.W. Perram, and E.R. Smith, Proc,
R. Soc. London, Ser. A 111, 27 (1980).

H.B. Dwight, Tables of Integrals and other
Mathematical Data, 4th ed. (Macmillan, New York,
1961).

.R.J. Elliot, J. A. Rrumhansl, and P.L. Leath, Rev.

Mod. Phys. 464 465 (1971).

F.L. Galeener and G. Lucovsky, Phys. Rev. Lett.
324 1477 (1976).

F.L. Galeener, A.J. Leadbetter, and M.W.
Stringfellow, Phys. Rev. B 214 1082 (1983).

J.M. Gales, B. Boltjes, M.F. Thorpe, and S.W. de
Leeuw, (submitted to Phys. Rev. B) (1989).

J.M. Gales (unpublished).

H.C. Gupta, S.D. Mahanti, and S.A. Solin, Phys.
Chem. Minerals 16, 291 (1988).

M. Ishii, T. Simanouchi, and M. Hakshira,
Inorganica Chemica Acta 11;, 38 (1967).

La

Li

Lo

Ly

Ja

Ke

Ma

Mc

Pi

Pi
So

Th

Th

82

85

34

41

74

66

70

71

88

89

89b

77

83

88

76

86

178

Landolt-Bérnstein, Numerical Data and Functional
Relationships in Science and Technology New
Series Group III Vol. 17a, (Springer-Verlag,
Berlin, 1982).

F. Liebau, Structural Chemistry of Silicates,
(Springer-Verlag, Berlin, 1985).

A.E.H. Love, A Treatise on the Mathematical Theory
of Elasticity, 4th ed. (Cambridge University
Press, Cambridge, England, 1934).

R.M. Lyddane, R.G. Sachs, and E. Teller, Phys.
Rev. 524 673 (1941).

J.D. Jackson, Classical Electrodynamics, 2nd ed.
(Wiley, New York, 1974), Chap. 7.

P.H. Keating, Phys. Rev. 151, 637 (1966).

J.B. Marion, Classical Dynamics of Particles and
Systems, 2nd ed. (Academic Press, New York,
1970).

A.A. Maradudin, E.W. Montroll, G.H. Weiss and J.P.
Ipatova, Solid State Physics Supplement 3
(Academic, New York, 1971), Chap. 6.

T.C. McGlnn, M.V. Klein, L.T. Romano and J.E.
Greene, Phys. Rev. B 11, 3362 (1988).

Newman et al., Phys. Rev. B 124 657 (1989).

Dan Neumann (private communication).

R.M. Pick and M. Yvinec, in Proceedings of the
International Conference on Lattice Dynamics,
Paris, 1977, edited by M. Balkanski (Flammarion,
Paris, 1977).

T.J. Pinnavaia, Science 2294 365 (1983).

S.A. Solin and H. Zabel, Adv. Phys. 11, 87 (1988).
M.F. Thorpe, in Physics of Structurally Disordered
Solids, 1976, edited by S. Mitra (Plenum, Mew
York, 1976).

M.F. Thorpe and S.W. de Leeuw, Phys. Rev. B 11,
8490 (1986).

Va

Ve

Wa

We

Yo

86

67

86

70

85

179

M.T. Vaughan and S. Guggenheim, J. Geophysical
Res. 21, 465 (1986).

L. Verlet, Phys. Rev. 112g 98 (1967)

R. Wada, W.A. Kamitakahara, M. Suzuki, in
Proceedings of Symposium K, 1986 Fall Meeting of
the MRS, edited by M.S. Dresselhaus, G.
Dresselhaus, S.A. Solin (Materials Reasearch
Society, Pittsburgh, 1986), p. 174.

Wert and Thompson, Physics of Solids, 2nd ed.
(McGraw-Hill, New York, 1970).

B.R. York, S.A. Solin, N. Wada, R.H. Rayathatha,
I.D. Johnson, T. Pinnavaia, Solid State Commun.
11‘ 475 (1985).