ETIJEGY TRANSPORT IN THE THREEDIMENSIONAL, HARMONIC,
ISOTOPICALLY DISORDERED CRYSTAL
By
JOHN DAVID POWELL
A DISSERTATION PRESENTED TO THE GRADUATE
COUNCIL OF THE UNIVERSITY OF FLORIDA IN PARTIAL
FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1972
ACKNOWLEDGMENTS
The author wishes to thank the members of his supervisory
committee, Dr. D. B. Dove, Dr. J. W. Dufty, and Dr. T. A. Scott for
their kind consideration. A special word of appreciation must be
extended to Dr. K. R. Allen who suggested this problem to the author,
who patiently supervised the research, and who offered countless sug
gestions and much encouragement. Thanks are also due Mrs. Carolyn
Ott for a number of helpful discussions. Finally, the author wishes
to thank the National Science Foundation, the Graduate School of the
University of Florida, and the Physics Department who awarded him
fellowships for four years, one year, and one trimester, respectively,
to pursue his graduate education.
TABLE OF CONTENTS
Page
ACKNOWLEDGMENTS ....................................................... ii
LIST OF FIGURES............................................... v
ABSTRACT................................. ........ .......... vi
CHAPTER I: INTRODUCTION...................................... 1
1. Phenomenological and Related Theory................... 1
2. The OneDimensional Model............................. 6
3. The Present Calculation...............................
CHAPTER II: FORMALISM ........................................ 12
1. Kubo Formula .......................................... 12
2. The Average Heat Flux Operator........................ 16
3. Equations of Motion................................... 21
4. Synthesis of Results.................................. 26
CHAPTER III: THE SCATTERING PROBLEM.......................... 29
1. Scattering by a Single Plane of Impurities............ 29
2. Scattering by Many Planes............................. 36
9
3. Calculation of the s( ............................ 45
4. Discussion............................................ 48
CHAPTER IV: HEAT CURRENT.................................... 53
1. The Uniform Distribution.............................. 53
2. Approximation for the Random Distribution............. 64
CHAPTER V: SUMMARY AND DISCUSSION............................ 68
iii
TABLE OF CONTENTS (continued)
Page
APPENDIX A: THE METHOD OF STATIONARY PHASE................... 81
APPENDIX B: CALCULATION OF G ,(q,q') FOR THE
zss
UNIFORM DISTRIBUTION.......................... 84
APPENDIX C: APPROXIMATION OF G z(q,q') FOR THE
RANDOM DISTRIBUTION........................... 91
LIST OF REFERENCES........................................... 98
BIOGRAPHICAL SKETCH.......................................... 101
LIST OF FIGURES
Page
1. The model................................................... 10
2. Orientation of vectors............... ... ................. 31
3. Scattering by many planes................................. 38
4. Directions of scattered waves.............................. 51
5. / versus 0 for the sample containing
z z
uniformly distributed impurities.......................... 77
Abstract of Dissertation Presented to the
Graduate Council of the University of Florida in Partial
Fulfillment of the Requirements for the Degree of Doctor of Philosophy
ENERGY TRANSPORT IN THE THREEDIMENSIONAL, APMONIC,
ISOTOPICALLY DISORDERED CRYSTAL
By
John David Powell
August, 1972
Chairman: Kenneth R. Allen
Major Department: Physics
A formalism is developed for studying the effects of iso
topic impurities upon the conduction of energy in a threedimension
al dielectric. The model used is the analogue of a onedimensional
model introduced by Allen and Ford in 1968. The Kubo method of cor
relation functions is the particular technique employed to determine
the current conducted through the solid. Emphasis is placed upon
obtaining relevant mathematical techniques for treating the problem
in three dimensions.
Two special distributions of impurities within the solid
are considered and solutions are found for the current in each case.
A detailed treatment is given for the case of isotopes uniformly
distributed among the host atoms, taking into account the effects
of multiple scattering of phonons by impurities. A calculation in
volving the first Born approximation is undertaken for the case in
which isotopes are distributed at random throughout the crystal.
vi
The results are compared with the results obtained in one
dimension. Some discussion of problems remaining to be solved in
three dimensions is given.
CHAPTER I
INTRODUCTION
Attempts to place the theory of energy transport in solid
insulators on a firm foundation and to obtain a reliable model from
which calculations can be made have been underway since the early
1900's. Although research on this problem has been extensive, the
inherent difficulty in treating nonequilibrium systems has posed a
formidable barrier to success. Indeed, the very early investigators
were forced to content themselves with qualitative predictions or, at
best, calculations involving rather drastic approximations. Before
discussing the specific aims of this dissertation, it therefore seems
advantageous to present a historical survey of the research done on
this problem up until the present time. The intent will be to indi
cate how a more quantitative understanding of this problem has
emerged. It should be emphasized, however, that this is not intended
as a detailed review, but rather a brief and completely qualitative
synopsis of what the author considers some of the more significant
advances. Particular emphasis is given to the work of the last few
years and, where applicable, pertinent experimental evidence has been
indicated.
1. Phenomenological and Related Theory
Energy transport in solids is most conveniently discussed
by defining a quantity known as the thermal conductivity. From
experimental evidence, one knows that if a temperature gradient
3T/3z is introduced across an isotropic solid, the energy which
flows through the solid obeys the Fourier heat law given by1
J = K .1
z 3z
J is the current or energy per unit area per unit time which flows
z
in the direction of the gradient and K is the thermal conductivity.
The negative sign indicates that the current flows from higher to
lower temperatures.
The first evidence of any experimental effort to deter
mine the factors upon which the thermal conductivity depends is the
2
early investigations of Eucken. Eucken found that for most solid
dielectrics the thermal conductivity varied as 1/T at high tempera
3
tures. Debye, later in 1914, succeeded in explaining the 1/T de
pendence using kinetic theory and thermodynamic considerations. The
argument was entirely classical, not intended to be rigorous, and
provided no adequate means for determining the proportionality con
stant.
The first really significant advance was suggested in 1929
when Peierls, using the newly developed idea of quantized lattice
waves, succeeded in formulating a theory valid at low temperatures.
Peierls' formulation was based on the Boltzmann transport equation and
formally included the effects of the socalled phononphonon inter
actions. These interactions arose from the retention of terms of
higher order than quadratic in the crystalline potential; they con
sisted of the normal processes which conserve the crystal momentum
and the umklapp processes which do not. In effect, Peierls was able
to obtain an integral equation which represented a solution to the
collision term in the Boltzmann equation, though he was unable to ob
tain a solution. However, his theory did predict that the thermal
conductivity should increase exponentially with decreasing tempera
ture in the low temperature limit.
De Haas and Biermasz57 attempted experimental verifica
tion of Peierls' results. They found that the thermal conductivity
did increase with decreasing temperature though not so rapidly as
Peierls' theory had predicted. Though it was not understood at the
time, the discrepancy is now known to have arisen from the resistance
to the flow of energy caused by point defects which had not been in
cluded in Peierls' analysis.
During the next several years, a number of noteworthy
attempts were made to extend Peierls' original theory. Specifically,
investigators sought to account for the various mechanisms which
cause resistance to the flow of energy in a crystal lattice. Most
calculations involved a relaxation time approximation to the
8
Boltzmann equation. In 1938, Casimir suggested that the mean free
paths of phonons should be limited by the finite size of the sample.
This effect was known as boundary scattering and was shown to be the
dominant form of resistance to the flow of energy at very low tempera
tures. Casimir's theory predicted that at low temperatures the ther
3
mal conductivity should vary as T3 and depend on the size of the spec
imen, an effect which had been observed in the earlier experiments
of de Haas and Biermasz. Later in 1942, Pomeranchuk pointed out
that point defects, or variations in the mass distribution through
out the lattice, could cause thermal resistance. As indicated pre
viously, this accounted for the fact that the thermal conductivity
did not rise exponentially with decreasing temperature at low temp
eratures. Some years later, Klemensl0 showed that the normal pro
cesses were unimportant except at low frequencies and that the um
klapp processes were the only anharmonic processes which resisted
the flow of heat. Finally, in 1959, Callaway,1 using frequency
dependent relaxation times, included the effects of the normal pro
cesses at low frequencies.
This technique of using the Boltzmann equation to inves
tigate energy transport which was begun by Peierls and extended by
the workers noted above has come to be known as the phenomenologi
cal theory. Excellent reviews of this theory are available 2'13
and a detailed discussion will not be given here. It is interest
ing, however, to point out two important consequences of the theory.
First, calculations generally involve making a priori assumptions
concerning the manner in which the system under investigation can
be expected to return to equilibrium when disturbed from equilib
rium. Second, the theory predicts an infinite thermal conductivity
for those systems which contain no scattering mechanism for causing
their return to equilibrium. This second fact was explained by
arguing that in any realistic crystal the divergence is eliminated
by boundary scattering at low temperatures and by anharmonicities
at high temperatures.
Some excellent experiments conducted by Berman and Brockl4
in 1965, attempted to verify the results of the phenomenological
theory. Fairly good agreement with theoretical predictions was ob
tained for most of the temperature ranges considered.
We may summarize the phenomenological theory by noting
that the theory did produce results that could be fitted to the
experimental data reasonably well. However, from a theoretical
viewpoint, the theory was rather unsatisfactory. In particular,
the results obtained for systems which contained no scattering
mechanism were meaningless. Clearly there was not contained within
the theory a systematic method of examining how the effects of
phonon scattering become important and might be expected to lead to
a theoretical prediction of the Fourier heat law. Some more specif
ic faults of the theory will be noted in the following section.
In the 1950's, Kubo5 developed a technique for treating
transport problems which was not based on the Boltzmann equation
but rather which depended on the solution of the dynamical equations
of motion of the crystal. Considerable credence was afforded the
Kubo method by the subsequent work of Hardy and his coworkers6 who
managed to obtain from the Kubo formalism a perturbation expansion
for the lattice thermal conductivity. Though the analysis was rigor
ous and did succeed in placing the phenomenological theory on a more
satisfactory basis, it was quite formal and not readily amenable to
calculations. In practice, it was necessary to resort to approxi
17,18
nations which just reproduced the transport equation718 and which
led to the same results as the phenomenological theory. Furthermore,
no attempt was made to perform calculations using a specific model.
2. The OneDimensional Model
1921
During the 1960's, several investigators independently
examined the possibility of using onedimensional chains of atoms to
study energy transport. Since the calculations of this dissertation
are based on an extension of the specific model used by Allen and
19
Ford,9 a brief discussion of the model will be given. The model
consisted of an infinite chain of atoms connected by Hookeian springs.
A finite section of the chain of length L0 contained both host atoms
as well as isotopic impurities randomly distributed among the host
atoms. The infinite regions on either side of this finite section
consisted only of host atoms and were used to simulate high and low
temperature reservoirs which caused energy to flow through the length
L0. A first principles calculation using the Kubo theory was under
taken to calculate the current and consequently the thermal conductiv
ity. The model was amenable to a formally exact solution in the har
monic approximation, but approximations were necessary to permit ex
plicit evaluation of the thermal conductivity. In essence, Alien
and Ford resorted to a Born approximation solution to the scattering
problem.
It was found that the model produced a nondivergent cur
rent even in the harmonic approximation and even for the special case
of an isotopically pure chain for which there was no scattering mech
anism. Allen and Ford emphasized that the current should ultimately
be limited by the amount of energy available for transport and hence
always nondivergent. This plausible argument had previously been
stressed by Erdos22 who had performed calculations using an indepen
dent method and found a nondivergent result for the isotopically
pure threedimensional crystal.
Nonetheless, the nondivergent result was at variance with
the phenomenological theory and consequently Allen and Ford carefully
reviewed that theory. They found that some of the derivations in
the phenomenological theory were not internally consistent. Specif
ically, it was noted that some of the approximations led to the pre
diction of negative occupation numbers of the phonon states at low
frequencies, a prediction which clearly made no physical sense. It
is now believed, therefore, that inconsistencies in the phenomeno
logical theory produced the low frequency divergence in those systems
which had no tendency to equilibrate.
Though a divergent result was not obtained, the approxima
tions did fail to produce a thermal conductivity in the usual sense.
In particular, it was found that the heat current was proportional
to the temperature difference between the reservoirs rather than to
the temperature gradient across the length LO. This is equivalent
to obtaining a thermal conductivity dependent upon the length of the
sample which is experimentally known not to be the case. Allen and
Ford speculated that either the multiple scattering effects (omitted
due to the Born approximation) or the anharmonic processes might
cause sufficient resistance to the flow of energy to produce a length
independent thermal cond'ictivity.
Early in 1972, Greer and Ruben23 succeeded in calculating
the multiple scattering effects in the onedimensional problem and
concluded that these were insufficient to produce a lengthindepen
24
dent result. Some recent work done by Allen, Ott and Powell,2
however, indicates that multiple scattering is sufficient if one
assumes an appropriate distribution of impurities within the impu
rity bearing segment of the chain.
The primary conclusion that can be drawn from these cal
culations made in one dimension is that, unlike the phenomenological
theory, a meaningful thermal conductivity can be obtained even for
the harmonic system with no mechanism for causing thermal resistance.
In addition, the onedimensional model has shown considerable prom
ise of providing a systematic method whereby one can study the ef
fect on the thermal conductivity as various scattering processes are
included. Ultimately one hopes to understand the manner in which
these scattering processes produce a current proportional to the
temperature gradient rather than to the temperature difference. In
the discussion which follows, it will be seen that similar conclu
sions are possible for a threedimensional model.
3. The Present Calculation
For reasons indicated in the last section, Allen and Ford
have felt obligated to attempt an extension of the onedimensional
model to three dimensions. In 1969, they succeeded in calculating
the thermal conductivity for the threedimensional, isotopically pure
25
crystal,5 and concluded that it was possible to obtain a meaningful
result in three dimensions for the case in which no scattering mech
anism was present. It will be the aim of this dissertation to pro
vide a formalism whereby one can study how scattering becomes impor
tant in three dimensions, to introduce some techniques applicable to
calculations, and finally to perform some calculations for some rather
special cases. It should be emphasized, however, that the aim is not
a complete solution of the threedimensional problem; unfortunately,
the problem in three dimensions has been sufficiently complicated that
much work remains to be done.
The model to be used is the analogue of that used by Allen
and Ford in one dimension. The crystal is taken to be infinite and
contain a slab of infinite crosssectional area but finite thickness
L0 along the z axis. (For orientation of the coordinate system, see
fig. 1.) Isotopic impurities of mass M' as well as host atoms of
mass M are distributed throughout the slab of thickness L0 which will
henceforth be referred to as the sample. The infinite isotopically
pure regions on either side of the sample are used to simulate high
and low temperature reservoirs. They are maintained at different
temperatures T1 and T2 where T1 > T2. The resulting temperature dif
ference causes a current to flow through the sample parallel to the
z axis. This current, and finally the thermal conductivity, will be
calculated using the Kubo formalism.
The restrictions noted in the onedimensional model are
also present in this calculation. Namely, only the harmonic term is
retained in the crystalline potential and therefore both umklapp and
normal processes which are anharmonic have been neglected. However,
these should make a negligible contribution to the heat current at
the low temperatures which will be considered. In specific calcula
tions it is further assumed that the crystal has simple cubic sym
M
Z=LO
Figure 1
The model. Isotopic impurities are distributed through
out the sample which is labeled II. Regions I and III simulate
high and low temperature reservoirs maintained at temperatures T1
and T2. The sample and the reservoirs are assumed to have infinite
crosssectional area. For purposes of calculation the y axis is
directed vertically upward and the x axis into the page.
metry, that it is isotropic, and that the acoustic approximation is
valid. Finally, it is assumed that all three polarization modes in
the crystal propagate at the same velocity. Though anharmonicities
have not yet been included in either the onedimensional or the three
dimensional model, the additional assumptions made here are merely
for the purpose of simplifying the calculations. The model does not
restrict one to make these assumptions, but the author feels that
they do not sacrifice any essential understanding of the problem.
In Chapter II of this dissertation, the formalism necessary
for solving the problem is introduced and the formalism is cast into
a form amenable to calculations. In Chapter III, a solution to the
scattering problem is obtained for the case of a uniform distribution
of impurities. Multiple scattering effects are included. Attention
is given to the details of the calculation since the author believes
that the same techniques may very well be applicable to more general
distributions once the mathematical difficulties have been overcome.
In Chapter IV, the heat current is calculated for the case discussed
in Chapter III. In addition, a Born approximation calculation of the
current is made for a random distribution of impurities. Finally,
in Chapter V, the results are discussed and analogies with the one
dimensional problem are noted. The relevance of the solutions ob
tained to the experimentally observed phenomenon of Kapitza resist
ance is briefly mentioned. Finally, some discussion of work remain
ing to be done is given.
CHAPTER II
FORMALISM
The purpose of this chapter will be to present the formal
ism necessary for calculating the current for the model discussed.
In particular, the Kubo formula for the heat current will be given
and a method for obtaining the current operator which must be used
in conjunction with that formula is outlined. Third, the Green's
function technique for obtaining the equations of motion of the crys
tal is discussed. Finally a formal expression for the current in
terms of the quantities to be calculated will be obtained. Second
arily, this chapter will aim to set the notation used in the remain
der of the discussion and to present some rudimentary ideas relating
to the phonon representation. In cases where certain aspects of
the formalism have been developed elsewhere, only results are given,
though specific references are always indicated.
1. Kubo's Formula
The Kubo theory of correlation functions is a technique
which is capable of producing expressions for the density matrix for
certain nonequilibrium systems. The mathematical details of the
theory have been discussed elsewhere26,27 and will not be discussed
here. The general assumptions of the theory as discussed by Mori
et al.26 will be briefly enumerated, however.
et al. will be briefly enumerated, however.
In the Kubo theory it is assumed that the system under in
vestigation may be divided into macroscopically small segments each
characterized by a local temperature T(r). Initially, these seg
ments are not allowed to interact with one another and each is as
sumed to be in thermal equilibrium. At time t = 0, the segments are
allowed to interact and bring about a change in the density matrix
of the system. It is further assumed that the system will approach
a steady state and that the temperature as a function of position
does not change as the steady state is approached. The implications
of this last assumption will be further discussed in Chapter V.
At time t = 0, the density matrix of the system is given
by the local equilibrium distribution
p(0) = exp[fI(r) H(r) dr] II1
p(0) dr]
where Z is the partition function, H(r) the Hamiltonian density, and
B(r) = k T(r) where k is Boltzmann's constant. The time evolvement
of the density matrix is determined by the quantum mechanical equi
valent to the Liouville equation given by
ib p(t) = [H,p(t)] II2
3t
where
p(t) = exp (iHt/Ti) p(0) exp (iHt/T) 113
In these equations, H is the Hamiltonian of the system. It is de
sired to evaluate the density matrix at some time T after the steady
state has been reached but before the system approaches complete
thermal equilibrium. Thus for the model at hand, we wish to choose
T large enough so that ample time has been allowed for energy to
propagate through the sample, but small enough so that energy has
not propagated all the way out to the ends of the reservoirs. The
first condition allows the steady state to be approached. The sec
ond condition prevents energy from being reflected at the ends of
the reservoirs, travelling in the opposite direction, and causing
the entire system to approach thermal equilibrium. As a rough esti
mate, we may in general choose T such that
L0 << T << L 114
v v
where v is the phonon speed and L is the total length of the crys
tal along the z axis. Since L is to be made infinite, T may be
z
chosen as infinite.
Once p(t) has been found the statistical average of the
_>
steady state value of the current may be found by evaluating
= Tr J p(t) 115
where J is the quantum mechanical current operator. The actual deri
vation of p(t) andis lengthy and quite difficult. Since this
derivation is adequately discussed in references noted previously,
only the result will be given here.
For this model, the average value of the heat current is
given in the Kubo formalism by26
A
<>= T f ds f dX Tr [J(0) Jz (t+iXh)p0] AT. 116
0 0 Z
In this equation, A represents the crosssectional area of the sample,
AT the temperature difference between the reservoirs, T the average
temperature, and X is a parameter having units (energy)I. p is the
equilibrium density matrix given explicitly by
BH
e
PO = eH 117
Tr e
and J is the local heat flux operator averaged over the sample of
thickness L0 (to be discussed in the next section). In the Heisenberg
picture, the time dependence of the operator J (t+iXh) is given by
J (t+iXt) = exp(iHt/h) exp(XH) J (0) exp(?H) exp(iHt/h). 118
z z
If one substitutes Eqs. 117 and 118 in Eq. 116 and
evaluates the trace using Dirac bra and ket notation, one finds
= A ds f dX <9ej n> exp[X(E E )] exp[is(EnE )/M]
0 0 ,n n
exp(8 E )
x __EY AT 119
Z exp(B Ek z
k
where J E J(0). I> and In> represent eigenvectors of the Hamiltonian
having energy eigenvalues E and E and satisfy the orthonormality
n x/
condition
II10
The integrals over sand X can easily be performed to yield
= TrA
=
kT2
Z <J In> 6 p AT II11
k,n z
exp( E)
P = exp(BEk)
k
1112
Eq. II11 is the specific equation which will be used in the evalua
tion of the heat current and consequently the thermal conductivity.
2. The Average Heat Flux Operator
Hardy28 has shown that the local heat flux operator can
be obtained by requiring that it satisfy the continuity equation
H(x) + VJ(x) = 0 II13
* >
H(x) is the"time derivative of the Hamiltonian density operator
which is related to the Hamiltonian H of the system through the
Heisenberg equation of motion
Sx),H]
1(x) = (ih) [H(x),H]
1114
where
< In> = 6 n
For a general type Hamiltonian of the form
2
Pi
H [ + V. 1 II15
2mi i
1 1+VJ
Hardy has shown that the harmonic contribution to the local heat
flux operator is given by
2 +
1 p (xk)
) j A( xx )(xk )(ih) [ k ,V(x.)] 1116
Sk,j k k mk j
where the sum over k and j is to be performed over the entire crys
tal. In these expressions, xk is the position vector to the equi
librium position or lattice site of the kth particle, p(xk) its
momentum operator, and mk its mass. V(x.) is the potential associ
ated with the j atom and in the harmonic approximation is given
by
V(x) (xx) Ua( ) U(x ) II17
j 2 ,x
where B(xj) is the socalled force constant tensor obtained by
expanding the potential in a Taylor's series about its equilibrium
position and retaining only the quadratic term. Quite generally,
> 
O (x.) is a symmetric function of x.and is also symmetric with
+ *
respect to interchange of a and B. The vector U(x.) represents the
displacement of particle j from its equilibrium position. Greek
subscripts are used to denote Cartesian components of vector and
tensor quantities throughout. Finally, A(xxk) is a coarsegrained
18
delta function used by Hardy which he has chosen to specifically
represent by the formula
1 2
A(xxk) = 3/2 3 exp[xxk /I2] II18
In Eq. 1118, k is distance large microscopically but small macro
scopically. The function A(xxk) is normalized such that
dx A(xx) = 1 1119
By using the commutation relation for the position and
momentum operators and Eq. 1117, one can show that for i # j
2 4..
[p (x.),V(x.)] =ih P(x xj) pa(xi) U(x.) 1120
Hence Eq. II16 becomes
4.
1 + + + + + x P+X
J2(x) =  A(xxk)(x x) ( U (x). 1121
2 k,j k (XkXj) k
For further development of this equation, it is necessary to trans
form the operators p(xk) and U(x.) to the normal mode or phonon
kJ
representation. For the case of a single atom per unit cell, such
a transformation is accomplished by the wellknown transformation
29
equations
U(x) = Z [ B (x ) a(q,s) + B (x) a (q,s)]
q .2m W(q,s)
q,s Z
19
and
1122
Shr(q,s) +
4. 4 S ) + ** 4 t +
p(x) = i 2 [ B (x) a(q,s) B (x) a (q,s)] .
q,s
Here w(q,s) is the frequency of phonons having wave vector q and
t >*
polarization s. a (q,s) is the creation operator for these phonons
and a(q,s) the corresponding destruction operator. These operators
satisfy the commutation relation
[a(q,s),a (q',s')] = 6, 6 II23
When operating on the state having N phonons of wave vector q and
polarization s, they have the effect
a (q,s) N(q,s)> = N(q,s)+l I N(q,s)+l>
and 1124
a(q,s) N(q,s)> = N(q,s) N(q,s)l>
where IN(q,s)+l> symbolically means the state in which N(q,s) has
been increased by unity and so on. Additionally, a (q,s) is the
**qs s
Hermitian conjugate of a(q,s) and B (x ) is the complex conjugate
 9.
of Bqs( ). Finally the Bqs(x ) satisfy the equation
 > >
2 qs (x xk)B qs 
W (q,s) Ba (x) = B (Xk) 125
p',k m/m
and are orthonormal in the sense that
* 's' Q
*B (x,) B (x ) = 6 6 1126
Sqq ss
If Eqs. 1122 are used in Eq. 1121, one obtains
J = I I ) x)
2 4 L A(xxk)(xkx.)
k,j a s s,s' /,T mj
1127
, 9
q s' *qs *q's' qs t
x [B (xk)BS (xj)a(q',s')a (q,s)Ba (Xk)B (x)a (q,s')a(q,s)] .
In obtaining Eq. 1127, those terms which do not conserve energy
have been neglected. As evidenced by Eq. II11, they make no con
tribution to .
The average value of the heat flux operator is given by
J2 V J dx J2(x) 1128
where the integral is to be done only over the sample and V repre
sents the volume of the sample. For those atoms lying outside the
sample, the coarsegrained delta function causes the integrand to
become negligibly small. Furthermore, if one assumes that
r . +
B (xkx ) is negligibly small when x kx is larger than a few
lattice spaces, the average heat flux operator can be written with
negligible error as
it + > +
= 2V Gss,(q,q') a (q,s) a(q',s') 1129
S +, s,s
q,q
where
(x x.) q's' *qs
G ,(q,q') = C (xx) B (x ) B (x.). 1130
ss .,j a,B / J
The notation is used to indicate that the sum is to be done only
over the sample.
Before concluding this section, it is necessary to note
some properties of the quantity Gss,(q,q') which will be used in
future calculations. First, inspection of Eq. II30 reveals that
** 3 1
G ,(q,q') = G, (q',q) II31
Second, if one performs the sum over j in Eq. 1130, Eq. II25 may
be used to prove that the z component of the resulting expression
is independent of z, provided w(q,s) = w(q',s), which is the only
case of interest. This rather intuitively obvious result just means
that there exist no sources or sinks for the current along the z
axis.
3. Equations of Motion
The final bit of formalism to be discussed is the difficult
problem of determining the quantities BqS((x) which appear in
30
Eq. 1130. The method to be used is due to Maradudin30 and employs
the Green's function technique.
By defining
+ 1/2
qs qs
(.) m. C (x A(q,s) 1132
J J J
and substituting in Eq. 1125, one finds
 
2 q qs c qs
m (q,s) Ca (x ) = a(xk) CB (x) II33
k,B
which is just the equation of motion of the th atom of the crystal.
Physically, CS(x ) is the displacement of the th atom from its equi
librium position provided only the qth mode is excited.
If we recall that the crystal under investigation has host
atoms of mass M as well as impurities of mass M' within the sample,
we may rewrite Eq. II33 in the form
> >4
L qs qs
L B(x ,xk) Cg (xk) = 6 LB(xV,xk) C (xk) 1134
k,B k,B
or in the more abbreviated matrix notation as
4 
L Cqs = 6 L Cqs II35
where
f 2 4 2
La (x xk) = M 6 6 k B(x k) II36
and
6L (xxk) = EM 2 6 6kn 6 n II37
r r r
In these equations, is the mass defect parameter given by
S= 1 M'/M and xn is used to denote the position of an impurity.
It may be noted that 6L and thus the righthand side of Eq. 1135
are zero if no impurities are present, and one just obtains the
equation of motion of the isotopically pure crystal. For the case
in which impurities are present, 6L is nonzero only when k and k
denote positions of impurities.
Maradudin has shown that the solution to Eq. 1134 can be
written in the form of a plane wave incident on the sample plus a
wave scattered by the sample, i.e.
qs
*qs qs qs
C(x) = C (x ) + Wq(x) II38
qs s
where C is the incident plane wave and W the scattered wave.
qs
The incident wave CO is a solution to Eq. II34 when no impurities
are present and thus
qs
L CO = 0 1139
If Eq. II38 is substituted in Eq. 1135, one finds
s qs 
L Wqs = 6 L C + 6 L W II40
l
Finally if the matrix G = L is defined, Eq. II40 can be written
in the form
+ qs 
Wqs = G6LC + G6LWqs. II41
0
The matrix G is known as the Green's function and has been shown by
Maradudin to be given by30
+ + a3 qr 1 1
G (x V V) f d
G 4 =a' k 3 2 2 2 2 2 2 2
k M(2r)3 q q 2v 2 +i6 w v q +i56
L T
1 1 44
+ 6 2 2 2i exp[iqx(Z,k)] II42
t v q +i6
where "a" is the lattice constant, vL and vT the group velocities of
phonons in the longitudinal and transverse modes, 6 a small real
number, and
x(,k) = xX xk II43
In writing Eq. II42, it has been assumed that the acoustic approxi
mation is valid, i.e.
v(q,s) = v q 1144
The integral in Eq. 1142 can be performed by integration in the
complex plane. The resulting expression, however, is quite complicat
ed and will not be given here. Since our ultimate concern will be
for the case in which vT = vL, we note that for this case Eq. 1142
becomes the rather simple expression
3
4 4 a exp[iwx(Z,j)/v] 145
G a ( x 4v2 x(S,j) 6
4 7TMv
where v = vL = vT. It should be emphasized, however, that these
velocities have been set equal merely to simplify calculations and
that there is no fundamental difficulty in using the more general
Green's function. In fact, calculations have been performed by
31
Allen and Powell for the case where vL # VT, though for a simpler
problem.
The ultimate aim is to solve Eq. II40 for Wqs A formal
solution can be obtained by writing
>*
qs
Wqs = [1G6L] G6LCO 1I46
For purposes of calculation we will expand the matrix [1G6L]
which yields
Sqs qs
qs = G6LC0 + G6LG6LC + ... 1147
The first term in Eq. 1147 represents the first Born approximation
and subsequent terms are of higher order in the Born series.
Once the value of the scattered wave W (x) has been
found, the quantities Bs(x ) can be immediately obtained as can be
seen from Eqs. 1132 and 1138. Quite formally, one may write
+ 1/2 qs q
B ) = m. A(q,s) [C0 (x) + Wq 1148
S 3
where Wqs(x) is given by Eq. 1146 or 1147. The normalization
constant A(q,s) can be determined by use of Eq. 1126.
This section completes the formalism required for a
complete calculation of given by Eq. II1i.
4. Synthesis of Results
The results of the preceding three sections may be sum
marized by noting that we have presented formal methods for calcu
lating all quantities necessary to evaluate the steady state aver
age heat current . Specifically, in section two the average heat
flux operator was obtained and in section three a method was given
for finding the quantities Bq(x ) upon which that operator depends.
In this section, the results will be combined as much as possible,
without loss of generality, in order that an expression for which
is perhaps more amenable to calculations may be obtained.
Recalling that our interest is only in the z component of
, we have upon substitution of Eq. 1129 into Eq. IIii
2 + G +
=TAh2 G ,(q,q') G 's (q'
4kT2V2 s,s' l,n
4kT V q,q' ,Z
s s ,s
q,q
S1149
6(a(q,s)(qs',s'))< la(q',s')a (q,s) n>
x P
If we define
M ssss"(qq' ,q,q" ) = < a(q',s')a (q,s) In>
l,n
II50
x pO
27
Eqs. II10 and 1124 may be used to show that
(qq', q' ) = 6 s' q '6,  )2
ss s s qq q q (e 1)2
Mss 'I"sI II
q q', q''" q"'', s s', s" # s'
Mss ,,s (q'q''q ) = ( (1)
(e l)(e 1)
+ e
(el)2
q = q', q" = q"', s = s' s" = s 'i
q~q~q=q ,s=s~s=
Mss ",,(q,q',q'",q ') = 0
, otherwise
where
4
c Tiw(q,s)
kT
1152
In deriving Eqs. 1151, use has been made of the fact that since
phonons obey Bose statistics
44>
SN (q,s) P = =
9
1
e 1
1153
2 >
N (q,s) pOZ
S =
=
Use of Eqs. 1150 and 1151 in Eq. 1149 yields
qq" ss"
1151
e+l
(e+ 1)
(ec_ 1)
TrAt2 1
= 2 Gss(qq)Gzss(q'') (e (e 'l)
jz> 4kT 2 zss ) zsIs' I e 1)(e 1)
s I S, S '
q,q
II54
+ 2 e ]
IG ,( q,q')I ) (ec5)2 6(w(q,s) w(q',s')) AT
zss (e 2 z
where use has been made of Eq. 1131.
If one wishes, the explicit expression for G ,(q,q') given
> ss
by Eq. 1130 and for B (x ) given by Eqs. II48 and II47 may be sub
stituted into Eq. 1154. No further simplification results, however,
in the most general case and this will not be done here. Eq. II54
will be used in Chapter IV where explicit calculations of the current
are undertaken.
CHAPTER III
THE SCATTERING PROBLEM
The purpose of this chapter will be to obtain a solution
for the scattered wave Wqs(x) discussed in the previous chapter and
>a S
the B (xZ) appearing in Eq. 1130. It will be assumed that the
impurities are uniformly distributed throughout the sample. It is
further assumed that the wavelengths of the phonons are much longer
than the distance between successive impurities. The significance
of the long wavelength assumption in the calculation of the current
will be discussed in the following chapter. In addition, the solu
tions for the Bs~(xz) are found only outside the scattering region,
but it will later be argued that these solutions are sufficient for
calculating the current.
1. Scattering by a Single Plane of Impurities
We begin by considering a single plane of atoms located
in the infinite plane z = 0. A certain fraction U of these atoms
are impurities and these are uniformly distributed among the host
atoms. A plane wave of the form
,q
CO (x ) = exp(iqx ) X(q,s) III1
where q is the wave vector and X(q,s) the polarization vector, is
incident from the left on the plane. Since the impurity plane will
30
be treated in the continuum limit, there is perfect azimuthal sym
metry about an axis normal to that plane. Thus the propagation vector
of the incident wave can be chosen to be in the yz plane and make an
angle c with the z axis. (See fig. 2.)
If we denote the first term in the expansion of Eq. II46
by Wl,then
,qs 23 q e ) *
Wqs (X a 3 exp[iqx(Zn )] exp(iqxnr) X(q,s) 1112
W1 (x) = a2 r r
4wv2 r x(n )
r
In Eq. III2 explicit use has been made of Eqs. 1137, II44 and
1145. If it is further assumed that the distance between succes
sive impurities is much smaller than the phonon wavelength, the
sum over r in Eq. III2 can be replaced by an integral over the im
purity plane with the result
s 2 exp[iq(x(nr )+y, sin a)]
W1 (x) a J dx s X(q,s). III3
1 rv2 r x(2n )
r
In Appendix A, it is shown that integrals of this form can be approxi
mated by the following formula which is a consequence of the method
32
of stationary phase :
ff g(x',y') exp[iqf(x',y')]dx' dy'
D III4
2 2i g(x ,y ) exp[iqf(xj,y.)]
q cjc jY j
o Z
Co
Figure 2
Orientation of vectors. The incident plane wave denoted
by CO is chosen to lie in the yz plane and make an angle a with the
z axis. The impurity plane is the z = 0 plane. x, and xnr repre
sent vectors from the origin to the laxtice sites of a host atom
and an impurity, respectively.
I I
In Eq. 1114, x. and y. are points in the domain D at which f(x',y')
J
is stationary and
xI 1
x. ,y
J 3J
32f
, 2
p y
I I
x. ,y.
J J
2
Yj 3x' y' '
.J .J
III5
The number a. is given by
J
, for aB >
J J
, for .B. >
3 J
. c .
2
Y2 a.
3 J
> 0
< 0
III6
2
,for a".. < .
Y j j
If Eq. III is applied to Eq. 1113, one finds after a tedious but
straightforward calculation
9
as
 (x
W1 ("
i Eoqa 
=2cos a exp(iqx) (qs)
2 cos aj ep1c1 ci
III7
The upper and lower signs on the vector ql are to be taken when
z cos a < 0 and z
cos a > 0, respectively, and ql and ql are given
q q = q(sin aj cos ak)
+ q (sin a + cos k)
q = q = q(sin aj + cos ok).
Throughout the remainder of the discussions ql will be used to denote
32f
Ja 3x'2
(. = 1
a. = 1
J
a. =i
3
the vector qxi + qj qzk where qx ay, and qz are the components
of the vector q.
A thorough physical discussion of the scattered wave will
be given later. For the moment, however, it is sufficient to note
that a plane wave incident on the impurity plane will produce a re
flected plane wave at a reflection angle equal to the incidence angle
and a plane wave scattered in the forward direction without deviation.
This can be seen by examination of Eq. III7. Clearly the case
z cos a < 0 represents the reflected wave and z cos a > 0 the for
ward scattered wave.
The first Born amplitude, represented by the coefficient
of the exponential in Eq. III7 should be a good representation of
the actual amplitude for small values of q, , 0 and for a not too
close to . It may be noted, however, that this amplitude is singu
lar at a = The source of this singularity can be traced to the
2'
fact that as a the scattering becomes very strong and the Born
approximation is no longer valid. Since ultimately the intensity
of this scattered wave will have to be integrated over all values q,
the singularity must be removed. The removal of this singularity can
be most effectively accomplished by consideration of the higher order
terms in the Born expansion. For example, the second Born term in
Eq. II47 can be written
23
S = E4 a 3 exp[iqx(Znr)] l("r) 9
7Trv r x(en )
r
and if one assumes that the first Born term given by Eq. III7 repre
sents an adequate solution to the amplitude inside the scattering
region, a calculation identical to the calculation of Eq. III7 can
be performed with the result
'c 2 
+ i_*aa )+
W2 (x) 21cosa exp(iqlx ) X(q,s) III10
Clearly this process can be continued to produce
Wq(x ) = Wi (x) = exp(iqlx ) X(q,s)
i
III11
x [ix(q) + (ix(q))2 +(ix(q)) + ...]
where
xt Eaqa
x(q) a= 21c III12
2 cos a
If we now assume that this expression can be resumed to produce a
result in the form of Eq. II46, we have
i ix() 
)= l+ix(q) exp(iql.x) (q,s) III13
and from Eqs. III1 and 1I38,
+qsix(+) +(ixs) III14
(x) = exp(iqx,) t(q,s) exp(iql ) (q,s) III14
l+ix(q)
Eq. III13 will be used to represent the value of the wave
scattered by the single plane of impurities. Admittedly, one has no
rigorous assurance that Eq. III7 represents a valid solution for the
scattered wave at points lying in the impurity plane as is discussed
in Appendix A and thus the above method of calculating the higher
order terms in the Born series is perhaps questionable. However,
Eqs. III13 and III14 do possess the following meritorious proper
ties:
(i) They contain no singularities.
(ii) They reduce to the first Born approximation
for small values of x(q).
(iii) They predict that as a the intensity
of the transmitted wave vanishes and the
intensity of the reflected wave approaches
unity. This is analogous to the wellknown
optical case of electromagnetic waves re
flected and transmitted by a dielectric
interface.
*Qs
(iv) Finally,they predict values of.the Bs (
which satisfy the orthogonality condition
given by Eq. 1126.
In short, though the technique used to calculate these higher order
terms lacks mathematical rigor, it does produce results which appar
ently satisfy all mathematical and physical requirements. In addi
tion, for values of the temperature which will be considered, x(q)
will be very small except for values of a exceedingly close to .
For these values it will be seen that the current is approaching zero.
>s
We therefore conclude that use of Eq. III14 for the C (x) cannot
cause significant error.
2. Scattering by Many Planes
In Eq. III14, we have essentially found the reflection
and transmission coefficients for a single plane of impurities. De
noting these by t and r, we have
1
t=
l+ix(q)
III15
ix(g)
r= .
l+ix(q)
>tq
In this section we find the values of C(x ) and consequently
>q 9
Bs (x.) for the case in which scattering can occur from N planes.
The planes are uniformly spaced a distance "a" apart along the z
axis in the length L0 of the sample. Multiple scattering is account
ed for. Though several modifications have been introduced, our ana
33 34 35
lysis will be similar to that of James,33 Prins,34 and Darwin35 who
obtained formal equations for the total electromagnetic wave reflect
ed from and transmitted by any given plane of atoms in a lattice for
the case of Bragg reflection. The solution of these equations is
36
similar to that obtained by Mooney for the case of optical filters.
We begin by numbering the planes consecutively along the
z axis beginning with the first plane at z = 0 and ending with the
N plane at z = (Nl)a. (See fig. 3.) We assume that a plane wave
having unit amplitude and unit phase at z = 0 is incident from the
left (i.e. cos a > 0). The amplitude and phase of the total wave re
flected from and of the wave transmitted down to the kth plane are
denoted by Rk and Tk, respectively. These values refer to the ampli
tude and phase just to the left of the kth plane. Except for the
first plane and the Nth plane, Tk consists of that part of Tkl
transmitted through the klth plane as well as that part of R which
is reflected from the lefthand side of the klth plane. Similarly,
R consists of that part of Tk reflected from the kth plane as well
as that part of Rk+1 which is transmitted through the kth plane.
The contributions to the waves for k = 1 and k = N can be determined
by similar arguments. Thus we have
TI =1
RI = Tlr + R2te
III16
Tk = T ktei + R re2i1
k k1 k
for k = 2,3...,I1
= Tkr + Rk+1tei
T = T tei + R re2i1
N N1 +
RN = rTN
k N
Rk
Tk
Co
Z=O Z=(kI)a Z=(NI)a
Figure 3
Scattering by many planes. Planes are numbered from left
to right beginning at the origin. CO represents a plane wave inci
dent from the left making an angle a with the z axis. R and Tk
denote amplitudes and phases of the total wave reflected from and
transmitteddown to the kth plane. That is, they refer to values
attained just to the left of the kth plane. Only a crosssectional
view of the planes has been depicted.
where 4 = aqlcos a is a phase factor which accounts for the dif
ference in path travelled by the various contributions to the waves.
Eqs. III16 can be solved for Rk and Tk in terms of Rkl
and T and written in the form
ki
T= 1
Rk_ T_1 r
tei
III17
i r2 e rp + re ,
T= te r2et T + R
k t TkI t  1
for k = 2,3...N
RN = rTN
or in matrix notation as
T = 1
1
= M for k = 1,2...N1
Tk+1 Tk
RN = rTN
In Eq. 11118, M is the twobytwo matrix
ixe
(1ix)e i
III18
III19
In obtaining Eq. III19, Eq. III15 has been used. Finally Eq.
III18 can be applied successively to obtain the following equations:
T = 1
STN
TN
= MN 1
III20
RN = rTN
N N
(l+ix).e
ixe
We denote the elements of the matrix MN1 by
MN1 a b21
S= III21
c d
and have from Eq. III20,
T,=
RN = rTN
RN = aR1 + bT1
TN = cR1 + dT1 .III22
These equations can be easily solved for R1 and TN yielding
dr b
1 a cr
ad bc
T III23
N a cr
R1 and TN represent the values of the amplitude and phase of the
total wave (including all possible multiple reflections) reflected
from the first plane and transmitted down to the Nth plane.
For further simplification, we note that the matrix M of
Eq. III18 has a unit determinant. For such matrices it has been
shown 37that the elements of the matrix MN can be written in terms
of the elements of the matrix M and the Chebyshev polynomials of the
second kind. Specifically,
SM UN_ (0)U 2(0)
MN M 121 UN1(0)
M 1 U N (0)
M21 UNl()
M U (0)
M12 UN1() )
M22UNl( )N2
where the Mi are the elements of the matrix M and 0 is given by
ii
0 = 1/2 (M11 + M22)
III25
The UNl(0) are Chebyshev polynomials of the second kind and are
given explicitly by the formula
UN() =
Nl
sin [(N+1) cos1 0]
VI2
These polynomials further satisfy the recursion relation
u.(o) = 2 Uj_ (0) U_2 (0)
and the orthonormality condition
uM(0) U (o) V 0 O dO
M N 2=
III26
III27
III28
6MN
Applying these results to find the elements of the matrix
III24
43
given by 11121, we have
a = (1+ix) ei UN2(0) UN3(0)
i
b = ix e U (0)
N2
c = ix e UN2( )
d = (1ix) ei' UN2(O) UN3(0)
N2 N3
where
0 = cos c + x sin
Finally, if we denote the values of the amplitude and phase of the
total wave transmitted by the N planes by T(q) = tTN, we have
total wave transmitted by the N planes by T(q) = tT we have
T(q) =
III31
>
Similarly the total wave R(q) = RI reflected by the N planes is given
by
III32
III29
III30
ix UN1 (0)
R(q) = () UN
(l+ix)U (0) el U (0)
N1 N2
(l+ix) U (0) ei U N(0)
N1 N2
In obtaining Eqs. III31 and 11132, Eqs. 11123, III27 and III29
have been used.
The analysis above has been for the case in which the in
coming wave was incident from the left or for cos a > 0. However,
since T(q) and R(q) depend only on the absolute value of cos a ,
Eqs. III31 and III32 are also valid for cos a < 0. We therefore
conclude that a plane wave having unit amplitude and phase incident
from the right onto the Nth plane produces a wave transmitted
through the first plane whose amplitude and phase are given by
Eq. 11131, and a wave reflected from the Nth plane whose amplitude
and phase are given by Eq. 11132.
We now wish to apply the above results to find the ampli
4.
tude and phase at any arbitrary point x Using the notation of the
previous section and recalling that part of the phase of the total
wave is accounted for by the factor exp(iq x), we have
have
4.
qs +
C (x
=T(q)exp[iq(Nl)aJcos aj]exp(iq*x ) X(q,s)
for z cos a > 0

CS(x ) = exp(iq*x ) X(q,s)+ R(q) exp(iq*x ) X(q,s)
for z < 0 and cos a > 0

Cqs(x) = exp(iq.x ) X(q,s)+R(q)exp[2iq(Nl)a cos aj]exp(iq*x
x A(q,s) for z > 0 and cos a < 0 III
33
In obtaining Eqs. 11133, we have taken the phase of the incident
wave to be unity at z = 0.
>q
3. Calculation of the B(x)
We now turn our attention to the calculation of the
Biq(x) which must be used in the evaluation of Eq. 1130. If
Eqs. III33 are substituted in Eq. 1132, one finds that Eq. 1126
is satisfied provided the normalization constant A(q,s) is given by
A(q,s) III34
In evaluating Eq. 11134, we have noted that the sum in Eq. II26
may be performed only over the pure regions of the crystal. This
is allowable since this sum produces a term of the order of the
volume of the crystal. In the limit then as the crystal becomes
infinite, the contribution to the sum for points x lying within the
sample becomes negligible.
We therefore have from Eq. II32 and III34 that for x
in the pure regions of the crystal
q _
B( = Cs (x) III35
where n is the total number of atoms in the crystal. For future
calculations as well as physical interpretation of the results, it
is most convenient to express the polarization vectors of the re
flected waves as a function of the propagation vector ql. To do so
we adopt the following convention: One transverse polarization vector,
denoted by A(q,l), will be directed along the positive x axis. (See
fig. 1.) Longitudinal polarization vectors, denoted by X(q,3) will be
directed along the direction of propagation of the wave. The second
transverse polarization vector, denoted by X(q,2), will lie in the
plane of incidence and its direction defined such that
(q,2) x X(q,3) = X(q,l) III36
Using this convention, it is easy to show that
(q,l) = X(qll)
q,2) = cos 2ca (l,2) + sin 2a X(q ,3)
(q,3) = cos 2a X(qi,3) sin 2a (q ,2) III37
Use of Eqs. III37, III35 and III33 yields the following rather
complicated expressions for Bs ):
For z cos a > 0,
Bs() + = T(q) exp[iq(Nl)alcos aj]
/n5
exp(iqx) X(q,s)
III38a
For z < 0, cos a > 0,
ql1 (
B (x ) =
1 e + e i l
 [exp(iqex ) X(q,l) + R(q)exp(iq.x ) q(q ,1l)]
.q2
B (x) [exp(iq *x ) X(q,2) R(q)cos 2a exp(iql*x ) X(ql,2)
~ /n
+B (
B (xi) =
+ R(q) sin 2a exp(iql.x) ( ,3)]
S[exp(iqx ) X(q,3) R(q)cos 2a
n(
 R(q) sin 2a exp(i l.x.) ('"1,2)]
For z > 0, cos a < 0,
1 +" "4+.
S [exp(iq"x ) X(q,l)
x exp(iql*x ) X(ql,l)]
+ R(q)exp[2iq(Nl)ajcos all
B (Xi)
1 . .... 4 .
 [exp(iqx) X(q,2)
vSn
 R(q)cos 2a exp[2iq(Nl)alcos all
x exp(iql x ) 1(ql,2) + R(q)sin 2a exp[2iq(Nl)a cos all
x exp(iql.x) (ql,3)]
3B (x
B (xi) 
S[exp(iqx.) X(q,3) R(q)cos 2a exp[2iq(Nl)a cos al]
/n
III38b
qB (x
B (xi)
exp(iql^'x) k (ql,3)
48
x exp(iql.x ) 1(ql,3) R(q)sin 2a exp[2iq(Nl)alcos all
x exp(iqlx ) X(ql,2)] III38c
4. Discussion
Eqs. III38 represent the central result of this chapter.
In view of the vast amount of mathematics necessary in their evalua
tion, as well as the fact that they are of interest quite apart from
their relevance to the problem at hand, we interrupt our calculation
of the heat current to summarize the calculations of this chapter and
to discuss the physical significance of the results obtained. We re
>s
call that we noted that the quantities B (x ) were simply related to
to
the solutions of the equations of motion of the atoms of the crystal
which could be written as the sum of a plane wave incident upon the
sample and a wave scattered by the sample. We began by considering
a single plane of atoms containing isotopic impurities within the
crystal and found that such a plane produced a reflected plane wave
and a transmitted plane wave. We then considered all possible multi
ple reflections of these waves as they progressed through the sample
of thickness L0 and determined the total wave transmitted by and
reflected from the sample. This wave was then used to determine the
>s
quantities B5(xj) given by Eqs. III38.
We begin our discussion of Eqs. 11138 for the case in
which zz cos a > 0. We recall that a is the angle made by the prop
agation vector q of the incident wave with the z axis and that
+>
BqS(x ) was initially written as the sum of this incident wave and
the wave scattered by the sample. Thus if the incident wave is
incident from the left (cos a > 0) and if one observes the effect at
a point x for to right of the sample (z2 > 0), one finds that the
scattered wave has the same direction as the incident wave and com
bines with it to produce a transmitted wave. The polarization of
the transmitted wave is the same as that of the incident wave and
consequently produces a disturbance in the same direction at the
point x as would the incident wave in the absence of any scattering.
The amplitude of the transmitted wave, however, has been reduced
from that of the incident wave by the factor given in Eq. III38a.
This reduction in amplitude results from the fact that part of the
incident wave has been reflected from the sample. For the case in
which both cos a and z2 are less than zero, exactly analogous re
marks apply.
For the case in which z cos a < 0, the situation is
slightly more complicated. Eqs. III38b predict that if one examines
a point x on the same side from which the incoming wave is incident,
one observes reflected plane waves. These waves are reflected such
that the angle of incidence equals the angle of reflection. If the
incident wave is polarized transversely perpendicular to the plane
of incidence (s = 1), there is only a single reflected wave polariz
ed perpendicular to the plane of incidence. If, however, the in
cident wave is polarized longitudinally (s = 3) or transversely in
the plane of incidence (s = 2), the reflected wave contains compc
nents which are both longitudinal and transverse in the plane of
incidence. These remarks are illustrated in fig. 4 where the vari
ous incident, reflected and transmitted waves are drawn for each
type of incident polarization.
Clearly the problem is very similar to the reflection and
transmission of electromagnetic waves through a dielectric interface.
The only significant difference is that lattice waves, unlike electro
magnetic waves, are capable of propagating longitudinal modes. It
is for this reason that one finds longitudinal incident waves pro
ducing reflected transverse waves as well as reflected longitudinal
waves and vice versa. One may note that the two reflected waves, as
does the transmitted wave, always produce a net disturbance at the
3
point xZ which is in the same direction as a disturbance caused by
the incident wave. Thus if the incident wave is transverse perpen
dicular to the plane of incidence, there is only a single reflected
wave also polarized in that direction. Similarly incident waves
polarized in the plane of incidence do not produce reflected waves
having a component of their polarization vector perpendicular to the
plane of incidence. The fact that a wave reflected from the left
hand side of the sample differs by a phase factor from a wave re
flected from the righthand side results from the fact that we chose
the phase of the incident wave to be unity at z = 0. This differ
ence is of no physical importance.
Before concluding this chapter, it is necessary to discuss
the approximation of replacing the sum in Eq. III2 by the integral
of Eq. III3. To do so we noted that it was necessary to assume
that the phonon wavelength was much longer than the distance between
4,
q51,
q,l
(a)
q,3
al
(b)
4q,3
(c)
Figure 4
Directions of scattered waves. Waves in the lower left
hand quadrant represent incident waves, those in the upper lefthand
quadrant reflected waves and those in the upper righthand quadrant
transmitted waves. The wave vector and polarization index of each
wave have been indicated. Each of the three possible types of
polarization of the incident wave is represented by parts (a), (b)
and (c). Though both reflected waves are reflected at angle a,
they have been displaced slightly for clarity in parts (b) and (c).
impurities. This assumption implies that the phase changes negligibly
in Eq. 1112 as the vector xr explores the impurity plane and replace
ment of the sum by an integral should be a reasonable approximation.
If a represents the fraction of atoms in the plane which are impurities,
we must require roughly that
1/2
X >> 27naa 11139
In the next chapter it is shown that at low temperatures this condi
tion can be easily satisfied by those phonon wavelengths that contri
bute significantly to the heat current.
CHAPTER IV
HEAT CURRENT
In this chapter the heat current will be calculated from
the solutions to the scattering problem obtained in the previous chap
ter. We recall that those solutions were valid only for a uniform
distribution of impurities. Therefore, we have in section two of this
chapter obtained an approximate expression for the current for a ran
dom distribution believed to be valid at somewhat higher temperatures.
The restrictions placed on the temperature due to the necessary re
strictions on the phonon wavelengths in the solution of the scattering
problem are evaluated.
1. The Uniform Distribution
Having solved the scattering problem for the uniform dis
tribution in the previous chapter, the heat current may be evaluated
by Eq. 154 in which we found
2 G [ss(qq)G .ss(q'q') )
z 4kT2V2 > s,s' s(e)e')
s q,q'
IV1
s , 1 2 e
G ,(qq)] 6(w(q,s)w(q',s')) AT
Iss (e 2 q z
where G ,(q,q') is given by Eq. 1130. The calculation of
zss
G ,(q,q') is quite tedious and adds little to the physical under
zss
standing of the problem. We have therefore performed this calcula
tion in Appendix B. Before indicating the results, however, it is
necessary to digress briefly in order to indicate one important
property of the model which is used in the calculation performed in
Appendix B.
Ordinarily one expects to have to use Eq. II25 and solve
a secular equation to determine the normal mode frequency spectrum.
However, this necessity can be avoided in the present case by use
of a theorem due to Lord Rayleigh.38 This theorem may be stated
as follows: If a single mass M contained within a system of equal
masses connected to one another by springs obeying Hooke's Law is
reduced (increased) by an amount 6M, all frequencies are unchanged
or increased (reduced) but never by an amount greater than the dis
tance to the next unperturbed frequency. The model at hand consists
of a crystal having an infinite volume and thus a continuous fre
quency spectrum. Isotopic impurities are contained within a slab
having an infinite crosssectional area, but finite thickness L0.
Consequently, if the impurity mass M' > M, the introduction of any
number of impurities within this slab cannot cause a finite shift
in the normal mode frequencies. The same is not true for M' < M,
however, as some of the normal mode frequencies may shift to values
greater than those for the pure crystal. These may be shown to cor
respond to complex values of the wave vector and are referred to as
localized modes. At low temperatures they do not contribute to the
heat current,19 and may be neglected. We can therefore use Eq. 1125
with M' = M to obtain the normal mode frequencies. We obtain
2 + + + + +
Mw (q,s) a (q,s)= S(x) exp(iq'xk X(q,s) IV2
B,k
Returning to the evaluation of Eq. IVi, we note that it
is shown in Appendix B that G (q,q) is an antisymmetric function
zss
>*
of q. Thus the first term in brackets of Eq. IV1 does not contri
bute to the current. In addition, it is shown that G ,(q,q') is
ss
nonzero only for q' = q or q' = ql where
A A A
ql = + yJ zk
From Eq. B19, we find
2
S* 2 12 0 2 2 4
G,(qq')1 = 1 2() I qq)
SI ss 2 z
s,s L
z
IV3
x [ () 6+ + IT(q)12 IR()12 6 ]
where it has been noted that
w(q,s) = w()
for the case in which all three modes of polarization propagate at
the same group velocity. In Eq. IV3
v() = v(q,s) = V(q)
Iv4
is the group velocity of phonons having wave vector q, L is the
z
total length of the crystal along the z axis (to be made infinite
in the final formulae), and T(q) and R(q) are to be obtained from
Eqs. III31 and 11132.
As noted previously, we must now take the limit as the
crystal becomes infinite and the frequency spectrum becomes con
tinuous. Consequently, sums over q and q' appearing in Eq. IV1
must be converted to integrals using the prescription
SV f dq
+ (21r)
IV5
3 
where V/(2r)) represents the density of states in q space and
V = L LL is the total (infinite) volume of the crystal. In addi
x y z
tion, Kronecker deltas must be replaced by Dirac delta functions
using the relationship
IV6
6 2 6(q q )
q xq L x
xx x
etc.
Using Eqs. IV3, IV5, and IV6, we find that Eq. IV1
becomes
37Th 2 2 4
= (2T 2 dq f dq' 2(q) Iv(q) I [ T(q)I \ 6(qq')
IV7
+ T(q)21 R(q) 2 6(qxx) (qy) 6(qz+) )2 AT
x x z z I (e_1 2 z
In obtaining Eq. IV7, it has been noted that from the properties of
delta functions
1
for q' = q
IV8
and
6 (qq) 6(qyqy) 6(q+qz)
for q' = q1
4
The integral over q in Eq. IV7 is to be done over the first
Brillouin zone.
From Eqs. III31 and III32, one finds
IT(q) 2 2 1 2
1 + x2() U _(0)
N1
IV9
IR(q) =
x (q) U_1()
1 + x2(q) U 2 ()
N1
and
6(qx q) 6(qyqy) 6(w(q)(q'))
6(qx q) 6(qyqy) 6 ((q)w(q',))
where x(q) and 0 may be obtained from Eqs. III12 and III30, re
spectively. If one substitutes Eqs. IV9 into Eq. IV7 and performs
4>
the integral over q', one finds
2 e
=2 dq 2 (q) IvZ( I
Z (27) 4kT2 (el)2
IV10
x 2222 AT
1+ E ga U2 ()
4cos2 N1
4cos a
where explicit use has been made of Eq. III12. In the acoustic
approximation, we have
2 ) 2 2
W (q) = v q
and
Ivz(q)l = vicos a IV11
Making these substitutions in Eq. IV7 and writing the integral in
terms of the dimensionless variable
kT
one finds after considerable algebraic manipulation that
TD/T
< 6k > 3 _4 e d_
= d2
z a3 0 (e1)
IV12
x f AT
0 202(T/TD)22 2 
01 + Do2(T/T D) sin2 (Ncos1 0)
y2 1 02
TD is the Debye temperature for this model which has been defined as
2 tv
T IV13
D = ka
and y = cos a. Also in obtaining Eq. IV12, we have made explicit
use of Eq. 11126. We may note that in terms of the parameter E ,
eo(T/TD)
0 = cos[2(T/TD )y] + sin[2(T/TD)Ey] IV14
In the limit of low temperatures, we may replace the upper bound on
the C integral in Eq. IV12 by infinity since the factor e /(e l)2
causes the integral to become negligibly small for large values of
b. For the same reason, we may approximate sin[2(T/TD)Ey] by
2(T/TD)y and consequently 0 by cos[2(T/TD)y] in the term 1 02
yielding
1 2 sin2[2(T/TD)y] = 4(T/TD)2 2 IV15
We may not, however, make this approximation in the term sin(Ncos 0)
without further consideration since this term is oscillatory and con
tains in its argument the factor N which can be quite large.
These approximations yield finally for Eq. IV12
6kv T 3 4 e
= 2a3 T 2 d
rz a D 0 (ei)
IV16
1 5
x y22 zA
0 y + sin (Ncos E)
Unfortunately, this integral cannot be evaluated analytically. How
ever, some special cases may be considered.
First, if C = 0 the integral can be easily performed yield
ing
72kv TT
z 0 2a3 T (4) Tz I17
Ta D
where r is the Riemann zeta function and we have used the fact that
f n e 2 d = n! (n) IV18
0 (e1)2
For large n, C(n) = 1. Eq. IV18 represents the heat current of the
isotopically pure crystal and agrees with expressions obtained pre
viously by Allen and Ford25 and by Erdos.22
Second, in the case of extremely low temperatures, one may
show that if N T/TD < 1, the term
sin2(Ncos1 0) = 4N2(T/TD)2 2y2 .IV19
If we use Eq. IV19 in Eq. IV16, we find
6kv T '3 0 e_
= 3 d
S 2a3 D 0 (e 1)2
IV20
1 3
S 2 y dy AT
2 2 2 2 2 2 AT
0 y + N E (T/TD ) C
This integral may be approximated for small E by relatively straight
forward though extremely tedious means to produce
= <1 + N2E22(T/TD)2
IV21
x [17 + 60(log(Neo T/TD y)
where 0 is given by Eq. IV17 and y is Euler's constant
y = 0.5772....
Finally, we consider the more likely case where N T/TD >> 1.
For that case Allen has developed a procedure for approximating
integrals of the form of Eq. IV16 employing integration in the
complex plane. This procedure will not be discussed but the results
of its application to Eq. IV16 reveal that for T/TD > 1,
= < E>[2 IV22
S< z0 + 2 IV22
A far simpler though less rigorous approximation than that develop
ed by Allen may be employed for the case where EU << 1. For this
case, the term sin2(Ncos 0) will oscillate rapidly if N T/TD >> 1
and, for small cO, it should not be a bad approximation to replace
sin2(Ncos1 0) by its average value of 1/2. If this is done the
integrals in Eq. IV16 become elementary and the result is
= [l IV23
z z I /2
which is in close agreement with IV22 for small Eg.
Discussion of Eqs. IV17, IV21 and IV22 will be deferred
until the next chapter. Before concluding this section, however, it
is perhaps worthwhile to discuss some of the approximations necessary
in the solution of the scattering problem. It may be recalled that
for the solutions of the scattering problem for the uniform distri
bution to be valid, we had to require roughly that
qa << 1 IV24
1l/2
63
where a represented the fraction of impurities in any given plane.
In terms of the variable ( this condition becomes
2(T/T D)
<< 1 IV25
1/2
It may be noted that the function 4 e /(e 1)2 is sharply peaked
and attains its maximum value at about 4 = 4, falling off rapidly
for larger and smaller values of E Thus the condition of Eq. IV25
may be replaced again approximately by
8 T/TD
<< 1
G1/2
or
T 1/2
T << 1 IV26
T 8
In a typical dielectric, TD is of the order of several hundred
degrees. Hence at low temperatures, Eq. IV26 can be satisfied
quite easily even for relatively small values of a and we conclude
that use of Eqs. III38 has been appropriate for the calculation of
the heat current. As a final point, it was noted in Chapter III
that the current tended to zero as cos a 0. This may be seen by
noting from Eq. IV16 that the integrand becomes proportional to
y as y = cos a tends to zero.
2. Approximation for the Random Distribution
An approximate value for the heat current for the case in
which the impurities are randomly distributed throughout the sample
is obtained in this section. It is assumed that the calculation of
. 
the quantities G ,(q,q') may be performed using the Born approxi
mation for the BS(x9). It is further assumed that the phonon wave
lengths are not long compared to the interimpurity separation, i.e.
X 2Ta IV27
where a' represents the fraction of impurities in the entire sample.
Thus for the results to be valid at relatively low temperatures, it
must be assumed that the fraction of impurities is quite small.
Using these approximations we show in Appendix C that the
term
q q G' ( ,q)G (q ,q')
S, s,s' z s (e l)(e 1)
IV28
x 
x S(w(q) W(q'))
which appears in Eq. II54 does not contribute to in the lowest
z
nonvanishing order of . In addition, it is shown that
I
sls
2
1 = 12L
IG 9,+ 2 2 +
s,(qq')j = 2 0 ~ (q)
zss, )
L
z
x [v(q) 
z
2 a 3L, (q)Iv (q)I
4wrv3
IV29
to lowest nonvanishing order. If we make this substitution in
Eq. 1154, take the limit as the size of the crystal becomes infinite,
and proceed along steps identical to the development of Eq. IV16, we
find
2
S 3Trh2 + 2) e
 2 dq w (q)
z (2T)4kT2 (eS1)2
IV30
[ sivi E2a3L0 G'
x IV (q)i (q) .
z 4Tv 3
If use is made of Eq. IV13, Eq. IV30 may be written
= 6kv T 3
2
z 2 3 T
TTa D
TD/T
0 e
0 (e 1)2
IV31
1 NE20'(T/TD) "4 4
x y r dy
0 1
As T/TD > 0 this integral may be easily performed yielding
= 0
z z0
 13440 NE 2'(T/TD4 (8)
STc 4
Before concluding this chapter and turning to a general
discussion of the results, we shall again determine the approximate
temperatures for which one might expect Eq. IV32 to be valid. We
initially assumed in the calculation that
qa > a'1/3 IV33
which in terms of the variable E becomes
2(T/TD)E ) '1/3 IV34
Again the function 8 e /(e 1)2 is sharply peaked attaining its
maximum about 5 = 8. Thus Eq. IV34 becomes
16(T/T) > 11/3
D
which implies
,1/3T
T >
16
IV35
IV32
Examination of Eq. IV32 reveals that the second term in brackets may
become quite large at temperatures satisfying Eq. IV35 except for
very small values of N, E, and a'. As the second term in brackets
approaches unity, one may suspect that the Born approximation is no
longer valid. Thus the restriction on the values of q coupled with
restrictions brought about by the Born approximation leads us to con
clude that this treatment for the random distribution may be expected
to be valid only for rather limited values of the parameters involved.
This point will be further discussed in the next chapter.
CHAPTER V
SUMMARY AND DISCUSSION
The purpose of this concluding chapter will be to summarize
the calculations made and to attempt to discuss the motivation which
has led to these calculations. The plausibility of the results will
be discussed with reference to the results of the onedimensional
model for which calculations can be made in a more rigorous manner.
For reasons to be discussed, the results here preclude detailed com
parison with experiment or even with other theoretical calculations.
Nevertheless, where applicable, attempts will be madeto explain why
the results are reasonable from a physical standpoint.
The primary emphasis of this investigation has been the
formal exposition of a problem in three dimensions equivalent to that
19
in one dimension initially set forth by Allen and Ford. Though we
ultimately hope to be able to explain experimentally observed phe
nomena, the intent of this dissertation has not been to obtain re
sults which could be fitted to the experimental data. Indeed, the
calculations have not been restricted to those systems for which
experimental results are available. Rather we have sought to intro
duce some general techniques for treating the problem in three dimen
sions and hopefully to provide some insight into the nature of energy
transport in threedimensional solids. Thus a rather detailed cal
culation of the current is given for a crystal containing isotopic
impurities uniformly distributed throughout. The results of this
68
calculation are given by Eqs. IV17, IV21 and IV22. The mathemati
cal techniques used in this solution will be applicable to more gener
al distributions of impurities as well. This point will be further
discussed later. In addition, an approximate solution is given for
the case in which the impurities are randomly distributed throughout
the sample. The result is represented by Eq. IV32.
If one wishes to define a thermal conductivity to be as
sociated with these currents, one may do so by writing
aT
= K  Vl
z 3z
where K is the thermal conductivity and 3T/3z the temperature gradient
which is assumed to be adequately given by
9T IAT
= V2
Dz L0
Using this definition, the thermal conductivity of the isotopically
pure threedimensional crystal (see Eq. IV17) becomes
72kvL T3
K a (4 ) V3
Ta D
For the case in which impurities are uniformly distributed through
out the sample, one finds from Eqs. IV21 and IV22,
K = K 1 + N2 202(T/T2 [17 + 60(log(NEO T/TD) y) (4
for N T/TD << 1, and
2 2 5
K = KO 1 + F a V5
for N T/T >> 1. Finally for the case in which the impurities are
randomly distributed throughout the sample, we have from Eq. IV32,
2 U
K = 13440 NE o'(T/T) (8)
K = K0 1 V6
It must be strongly emphasized, however, that the values of
the thermal conductivity obtained in Eqs. V3 through V6 all depend
on the thickness L0 of the sample. Consequently, they do not repre
sent thermal conductivity in the usual sense and are analogous to the
results obtained in one dimension by Allen and Ford19 mentioned pre
viously. Since experimentally solids are known to obey Eq. I1 where
K is independent of the thickness, some explanation of these results
is necessary. The variance of the results with Eq. I1 may be ex
plained by noting that the Kubo theory obtains an expression for the
current which may be written as
J = K + ... V
z z
where the dots represent higher order terms in powers of the gradi
ent. Eq. II11 represents the first term in this expansion and
higher order terms are assumed to be negligible. Furthermore, the
gradient which appears in Eq. V7 is that which is initially imposed
on the sample. It is assumed in the Kubo theory that this gradient
is not significantly altered as the steady state is approached. It
has been pointed out that this assumption is probably reasonable in
systems for which there is a strong scattering mechanism. Otherwise,
however, as has been pointed out by Allen and Ford,19 the Kubo theory
will not produce a thermal conductivity in the usual sense. Never
theless, one may still use Eqs. IV17, IV21, IV22 and IV32 to
discuss the steady state current in terms of the temperature differ
ence between the isotopically pure regions of the crystal. We may
conclude from the results above that multiple scattering effects
from uniformly distributed impurities do not limit the flow of energy
sufficiently to produce a thermal conductivity independent of the
thickness of the sample. Some additional discussion of this point is
given later. We now turn to a discussion of specific results.
For the case E = 0, one merely has the current of the iso
topically pure crystal given by Eq. IV17. If one notes that the
specific heat of the phonon gas is given by
C=96k T
c 2 TD i () V8
Tr a 3 D
Eq. IV17 may be written in the form
3
J > CvAT V9
This equation may be obtained from the results of Allen and Ford25
and of Erdos22 by setting the velocities of propagation of the vari
and of Erdos by setting the velocities of propagation of the vari
ous polarization modes equal. Eq. V9 is, of course, nondivergent
and its meaning is rather obvious. It simply states that in the ab
sence of any scattering, the current is ultimately limited by the
amount of energy available for transport.
In considering Eq. IV21, we may note that for extremely
small values of N T/TD, E, and a the dominant term in that expres
sion is the one involving the logarithm. Thus for purposes of dis
cussion, we write
z 10 + 60N222 2(T/TD)2 log(Neo T/TD) V10
Since NEa T/T << 1, the second term in brackets is negative and the
current is reduced from that of the pure crystal as expected. It
might at first appear surprising to find that the correction term
varies as (T/T D)5 log(T/TD) since it may be observed in Chapter III
that the intensity of a wave scattered by a single impurity varies
4
as w This is just Rayleigh scattering and one would ordinarily
expect it to produce a term proportional to T as has been found in
the phenomenological theory. The reason for the weaker dependence
on T is that although the amplitude of a wave scattered by a single
2
impurity does indeed vary as w the total amplitude of the plane
wave scattered by an entire plane of impurities is a more compli
cated function of w. For small frequencies and angles not too close
to 2 the amplitude is proportional to w. This difference results
from the fact that the scattering has been treated coherently with
very definite phase relationships between individually scattered
73
waves accounted for. Clearly, this is the proper treatment of a
scattering problem for which impurities are uniformly distributed
and for which the phonon wavelengths are long compared with the
distance between successive impurities.
The dependence of the correction term on N2 log N rather
than N may be explained similarly. One may show that at temperatures
satisfying the condition N T/TD << 1, only those phonons having wave
lengths greater than the thickness L0 of the sample contribute signi
ficantly to the current. Thus the waves scattered by the planes of
impurities within the sample are very nearly in phase. For sources
emitting sound waves in phase, it is a wellknown fact in the theory
of sound that the resultant intensity is just N2 times greater than
that from any one source. The result here is perfectly analogous.
The weaker dependence on log N results from the fact that the scat
tering has been treated nonisotropically and is rather strongly
dependent on angle.
It is well to point out at this point that Allen et al.24
have obtained a result in one dimension for the case of impurities
uniformly distributed throughout the sample that is very much ana
logous to Eq. IV16. For that case one may also show that the cor
2
reaction term to the current is proportional to N where N is the
number of impurities in the sample and to T3. The difference in the
temperature dependence results from the difference in the density
of states at low frequencies in one and three dimensions.
Finally, it should be mentioned that although Eq. IV21
was derived with the assumption that the impurities were uniformly
distributed throughout the sample, the result for this special case
(N T/TD < 1) may well be approximately valid for more general dis
tributions of impurities as well. No rigorous proof of this asser
tion will be given but it can be made plausible by a physical argu
ment. For phonon wavelengths of the order of magnitude of the
sample, one would not expect that the deviations of the impurities
from their uniformly distributed positions would be strongly felt.
This results from the fact that the waves scattered by individual
impurities will be very nearly in phase regardless of how the im
purities are distributed throughout the sample. An analogous result
has been noted in the onedimensional model where it is found that
the current obeys an equation similar to IV21 even for a perfectly
random distribution at temperatures such that N T/TD < 1 where N
is the number of impurities.
For the case in which N T/TD << 1, it was found in Eq. IV22
that
= <> [V1 + V11
z 2 2j
This expression accounts for all the effects of multiple scattering
and no assumption has been made regarding the magnitudes of E and a.
We note immediately that the result is independent of the number of
planes N and thus the thermal conductivity depends linearly on the
thickness of the sample. An analogous result has been found in one
24 23
dimension by Allen et al. and by Greer and Ruben. It is also
interesting to note that the current depends on T AT, a dependence
39
which is indicative of Kapitza behavior in three dimensions. Thus
if one employs a sample in which all atoms have mass M'(O = 1),
Eq. Vll could possibly be used to predict the experimentally observed
drop in temperature as the current proceeds from one dielectric to
another at very low temperatures. At present, the model has not been
sufficiently refined to account for all the effects which contribute
to Kapitza resistance. It is nevertheless interesting to compare
h0
qualitatively Eq. Vll with a result obtained by Little who has
explicitly treated the problem of Kapitza resistance in solids.
Little considers two dissimilar, semiinfinite solids joined at the
plane z = 0 and obtains a result for the current across the inter
face which may be written in the present notation as
z 32v
h v
AT is the temperature difference between the two media and P is a
complicated integral which involves the total transmission coeffici
ent of the interface. The integral is not explicitly evaluated by
Little but he shows that its maximum value is onehalf. The maximum
occurs whenever there is a perfect match of densities at the inter
face. If one assumes in the present case that the temperature with
in the sample is given by T =(T1 + T2)/2, then the temperature dif
ference AT' between the isotopically pure region on the left and the
sample is given by AT' = AT/2. Eq. Vll may then be written in the
form
144 TkT3 TE
= > 3 C(4) + AT' V13
z 23 h3
vh h
If one further notes that the maximum value of the quantity in
brackets in Eq. V13 is unity, whereas the maximum value of the quanti
ty F obtained by Little is onehalf, Eqs. V12 and V13 are identical
in form. The examination of the above model's ability to explain
Kapitza resistance is an obvious problem for future investigation.
In fig. 5, / has been plotted as a function of
the product ca. It is seen that as e  0. This results
z
from the fact that as e  the scattering crosssection becomes
infinite and no current is able to traverse the sample.
Finally it is interesting to discuss Eq. IV32 with regard
to the onedimensional model. It is seen that the current is reduced
from that of the pure crystal by a term proportional to N and
(T/TD)7. This T dependence is the same temperature dependence for
the correction term that has been obtained in the phenomenological
22
theory and by Erdos.2 The linear dependence on the number of
planes results from the fact that the scattering has been treated
completely incoherently and is again analogous to the familiar result
of incoherent scattering of sound waves in the theory of sound.
There one finds that if the phases of the waves emitted by the vari
ous sources are completely random the resultant intensity is N times
the intensity from any one source.
24 23
Allen et al.2 and Greer and Ruben23 have recently succeed
ed in obtaining a much better approximation to the current for a
randomly disordered chain in one dimension. Their approximation ac
counts for multiple scattering and is valid at lower temperatures.
In the Born approximation of their expression, one also finds that
< z >o
z 0 )
1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.I
I 2 3 4 5 6 EC
Figure 5
/ versus EG for the sample containing uniformly
z zO
distributed impurities.
the current is reduced from that of the pure chain by a term which
is linearly proportional to the number of impurities and which de
pends on the temperature. The temperature dependence is different,
of course, due to the difference in the density of states. The
author suspects that Eq. IV32 may well be valid at temperatures
lower than those assumed for the purposes of making the calculation.
It does not seem feasible at this point to attempt any
comparison of the results obtained here with experimental data for
a number of reasons. First, the case of uniformly distributed iso
topes, for which a rather complete calculation has been given, is
except for the case of Kapitza resistance clearly not relevant to
experiment. Second, the approximation obtained in Chapter IV for
the random distirubtion of impurities can be expected to be valid
only for very limited values of the parameters involved. Multiple
scattering effects have been completely neglected. Finally, it
should be noted that so far as the author is aware, the only experi
mental data available are for those solids for which one would expect
boundary scattering to diminate at low temperatures, that is for
long rods having small diameters. In fact, most experiments con
ducted at low temperatures have been undertaken expressly to verify
the effects of boundary scattering. These effects have not been
considered in this model. Though boundary scattering is unquestion
ably important in solids of certain geometrical shapes, the author
does not believe that it is of fundamental importance in understand
ing energy transport or even of practical importance in a wide variety
of solids. It is hoped that experimental evidence on solids for
which boundary scattering does not dominate at low temperatures will
be available in the future.
It is perhaps well to summarize the results of Eqs. IV21,
IV22, and IV32 by noting that in each case the current has been
reduced from that of the pure crystal. It is evident that a formal
ism has been developed from which one can study the effects of includ
ing scattering processes in three dimensions which has been the in
tent of this dissertation. Ultimately, one hopes to obtain a com
plete understanding of how the Fourier heat law comes about in three
dimensions. The formalism initially set forth by Allen and Ford in
1968 for treating the problem in one dimension has recently shown
considerable promise for yielding such an understanding for the one
dimensional case. The author believes that a similar understanding
is possible in three dimensions. However, the mathematical techni
aues in three dimensions are sufficiently more complicated that
working out the details may require some time. We conclude by in
dicating some specific problems needed to be dealt with in the three
dimensional problem.
Of foremost importance is the necessity for being able to
treat the effects of multiple scattering for more general distri
butions of impurities. In this regard, it may be noted that Allen
24
et al.2 have just recently developed a scheme for treating this
problem in onedimension. In essence, they have found that the Born
series may be reordered in such a way that the deviations of the im
purity positions from the perfectly uniform distribution may be aver
aged separately from the sums over the uniformly distributed posi
tions. Their results clearly indicate that the manner in which the
impurities are distributed throughout the sample is of utmost impor
tance in the determination of the thermal conductivity even at low
temperatures. Unfortunately, such a procedure has not yet been pos
sible in three dimensions due to the complicated form of the three
dimensional Green's function and the fact that the waves scattered
by individual impurities are spherical rather than plane waves. In
fact, it is not clear that one can always obtain plane wave solutions
to the scattering problem for an arbitrary distribution of impurities
in three dimensions.
The second important feature that has thus far not been
included in either the onedimensional or the threedimensional model
is the important effect of anharmonic terms in the crystalline poten
tial which give rise to phononphonon scattering. It has been assum
ed that at very low temperatures, the effects of these terms are neg
ligible. However, at high temperatures, the lattice motion is an
harmonic and in any extension of this model to the high temperature
limit, phononphonon interactions would have to be included.
APPENDIX A
THE METHOD OF STATIONARY PHASE
In this appendix, we seek an approximation to integrals of
the form
I = / f g(x',y') exp[iqf(x',y')] dx' dy' A1
00 0O
The approach will be to use the method of stationary phase. This
method has been discussed for the case of single integrals by vari
ous authors and comprehensively treated for the case of double inte
grals when f(x',y') is real by Jones and Kline4l and by Kline and
42
Kay.2 Since the rigorous theory is quite complicated, our discus
sion will be confined to some rather intuitive arguments.
If the function qf(x',y') is much greater than unity, one
would expect the exponential term to oscillate rapidly and tend to
average to zero. The only place where this would not be true is in
the neighborhood of the points where f(x',y') is stationary, i.e.
points where f(x',y') obtains a maximum or minimum value. Consequent
ly, a good approximation to the integral should be obtained by ex
panding f(x',y') in a Taylor series, retaining only the quadratic
term, and integrating the resulting expression. The points at which
f(x',y') is stationary can be found by requiring that
Df(x',y') 0 and f(x') = 0 A2
9x' Dy'
and solving for the required values of x' and y'.
I I
Denoting such a point by (x0 ,0) and expanding f(x',y')
about this point, we have
' 1 2 +1 '2
f(x',y') = f(x0 00 + a (xx (Y'Y0
1 +
+ 2 Y0(XI'X)(Y'y~0) + ,
A3
where
a2 f(x' y') 2 f(
0 o0
0x' x0,Y0 0
2 f(x',y')
YO
3x' y' x0,YO
Substitution of Eq. A3 in Eq. Ai, yields
X1,I I
y 1 I ty,
00)
I = exp[iqf(xQ,y )] g(x01Yg) f
00
+ ( y '2 +
80 0
Sep /[ '
f exp iq/2[a (x'x02
I0
S ,O'
Y0(x'x0 (y'y0)] dx' dy'
I I
where we have assumed that g(x',y') can be replaced by g(x0,Y0) in
the integral. This should be a valid approximation provided g(x',y')
varies slowly with respect to the exponential as it does for the
problem at hand. The integral appearing in Eq. A5 can be integrated
43
by use of integrals which have been tabulated. Performing the
A5
83
integral, one finds
27i 0 '
I = g(x0 ,y) exp[iqf(xo,yo)] A6
where
0 = + 1 for a00> Y2 O > 0
0 = for aO > Y "O < 0
00 O 0
00 = i for a 00< Y0
In the most general case, there will be more than one point at which
f(x',y') is stationary and thus I becomes
I 2 i Y g(xj,y ) exp[iqf(x ,yj)] A7
I I
where the (x.,y.) represent points at which f(x',y') is stationary.
For the problem at hand, we have from Eq. III2
2 2 2 1/2
f(x',y') = [(x x') + y_y') + z ] + y' sin a
where we have set x = x'. Even for values of the wavelength
S= 2r/q which are many lattice spaces long, the function qf(x',y')
is much greater than unity except as both x' and x approach the
origin. For such points the above approximation may not be valid.
However, we shall assume that such points cannot substantially affect
the answer and use the approximation of Eq. A7 generally. This
point has been further discussed in Chapter III.
APPENDIX B
4.4.
CALCULATION OF G ,(q,q') FOR THE UNIFORM DISTRIBUTION
zss'
The purpose of this appendix will be to calculate the
G ,(q,q') which must be used in conjunction with Eq. IV1 in the
ZSS
calculation of the heat current. From the factor 6(w(q,s)w(q',s'))
which appears in that equation it is clear that our only interest
is in the case where W(q,s) = w(q',s'). In Eq. 1126 we found that
Ct i^Z (x' x) qs )
G zss q,q') = (z.z.) Bq's'x )B~ sx.) B1
Si,j a, J/m.m.
SJ
and have previously noted in Chapter II that for w(q,s) = w(q',s')
if one performs the sum over j in Eq. Bl the resulting expression
is independent of z.. Therefore, if we further assume that
a (xixj) is negligibly small when x.x. is larger than a few
lattice spaces, we may use the values for >'s' (x.) and qs(x.) in
1 J
either of the two pure regions of the crystal to evaluate Eq. Bl.
We begin by considering the case in which both cos a and
cos a' are greater than zero, where a and a' are the angles made by
q and q' with the z axis. For this case, it is most convenient to
choose the solutions for the Bs(x.) and '(x far to the right
J .
of the scattering region, i.e. for z. and z. > 0. Thus from Eq.
J 1
III38a, we have
1 I ()i= 1ix (T
G ('q') (z.) Q ) T(q') T(q)
is nM 1,
x exp[iq(Nl)a( cos al lcos a' l)]X (q',s')
B2
x A (q,s) exp(iq'*x.) exp(iqx.)
> *
By making a change of variable and letting x = x.x., we have
m 1 j'
S S(qq z x ) T(q') T(q)
zss nN im m m
i,m a,0
x exp[iq(Nl)a( cos al cos a' l)] X (',s') XA(q,s)
x exp[i(q'q)*x.)] exp(iq*x ) B3
.9
The sum over the x and y components of the vector x. can be done
1
immediately by using the fact that
exp[i(qq')*x.] = n n 6' B4
x.,y 1 x y qx y
where n and n represent the total number of particles along the x
x y
and y axes. Furthermore, since q = q' and both cos a and cos a' are
greater than zero, we must have qz q. Thus we can write
exp[i(qq')x.] = n 6 6 6 B5
i 0qxx yy z z
where nO is the total number of atoms in the sample of length L .
For further simplification, we write Eq. IV2 in the form
= AX (,s) 4 (x ) exp(iq*x ) X (',s') B6
ca, m
where we have used the fact that
A (q,s) A (q,s') = 6S
cx cxS
By taking the
show that for
gradient
s = st
of expression B6 with respect to q one can
+ W(q,s) = v(q,s) = (qs)
q 214M(q,s) xa,
x (x m) X exp(iq*x ) (q',s)
m
B8
where v(q,s) is the group velocity of phonons having wave vector q
and polarization s. In addition, Hardy has shown for s # s', but
W(q,s) = w(q,s') which is true in the present case that
a (x) m exp(i x ) A (q,s') = 0
m m
m
Use of Eqs. B5, B8 and B9 in Eq. B3 yields
G ss(q,') =
zss
2iL
Sw(q ,s) vz(q,s) IT(q)12 6 6qss
L
z
a,S
B9
B10
W2 'qs) 6 sS,
where L is the total length of the crystal along the z axis (to be
z
made infinite in the final formulae). Eq. B10 was derived with the
assumption that cos a and cos a' were greater than zero. However,
if both cos a and cos a' are less than zero, the same results are
obtained. We conclude therefore that
G + 2iL
G (q,q') = W(q,s) v (q,s)IT( )2 6+6ss
zss L qq' s
z
Bll
for cos a cos a' > 0 .
The case in which cos a > 0, cos a'< 0 is slightly more
complicated. For the first of these cases, we choose the solutions
4.
I >
for B '(x ) and B (x.) far to the left of the scattering region,
3
i.e. z., z < 0. We can write from Eq. III38b the general expres
1 j
sion
B (xj) = n [exp(iqxj) (q,s) + fl(q,s) exp(iqlx )
x X,(ql,s) + f2 q,s) exp(iql'x ) X(q ,3) 6s2
+ f3(q,s) exp(iq1 xj) X (ql,2) 6s3] B12
where the f's are to be obtained from III38b. Similarly, we have
from Eq. III38a for this case
q's' 1
Ba (x.) [T(q')exp[iq(Nl)alcos a' ]
a/n
x exp(iq'x.) a (q',s')
B13
By proceeding in a manner identical to the evaluation of B1, one
finds
2iL
G (sq,q') = w (,s)v (q,s)T(q) exp[iq(Nl)alcos al]
zss L z
z
x flq,s)6ss' + f2( ,S)62 s3' + f3 q,s)6 6s2] 6,
for cos a > 0, cos a' < 0
where we have used the fact that
v (q) = v()
B14
and that
T(q) = T(ql)
as can be seen from Eq. 11131.
Finally, for the case in which cos a < 0, cos a' > 0, we
again choose the solutions for to the left of the scattering region
and proceed as before to obtain
2iL
G z,(sq') L (q,s)v (q,s)T(q)exp[iq(Nl)alcos al]
zss L z
z
x [f (q',s')6s, + f2 (q',s')6 6 + f q',s')6362] 6+,
1 ss 2 s'2 s3 3 s'3 s2 q q
for cos a < 0, cos C' > 0 .
B15
We see immediately from Eqs. Bll, B14, and B15 that
G (q,q') is nonzero only for q' = q or q' = q Furthermore,
zss 1
9 >
we see from Eq. B1, that G (q,q) is an antisymmetric function of
q since
v () = v (q
Zz
B16
Consequently, the first term in brackets in Eq. IV1 does not con
tribute and we need only evaluate
IG s (q,q)1
s
and
zss 12
I
s,s
From Eq. B11, we find
2
IG12 12L 2 q 14
I 2 t 2+ 2 L
S(q,) 2 (q) v (q) IT(q)
ss' L
B17
where we have noted W(q,s) = w(q) and v (q,s) = v (q) for the case
z z
where vT = v. In a similar fashion, by obtaining the values of
the f's in Eqs. B14 and B15 from Eq. III39b, one finds
Gzss 5 q
sIs,
2
12L0 2
= 2 (q)
L
z
2 )
z
IT(q) 12 R(q) 12
or in general
2
12Lo 2
L2
z
x 6) + IT(q)I2 jR()2
qq
B18
s,s'
Gzss
2
Vz( )
[ T(q)I
q 'q
B19
APPENDIX C
APPROXIMATION OF G ,(q,q') FOR THE RANDOM DISTRIBUTION
zssl
We shall obtain an approximation for the quantity
G ,(q,q') in this appendix for the random distribution. The
zss
calculation will involve the use of the Born approximation for the
solutions of the qs (x,) and a random phase approximation. From
Eq. 1126,
x (x.x.) q's' *qs
G 9s q,q') (z.z.) B B (x ) B (.) C1
i,j a, 1 /m.m.
i j
where the qs (x.) are given in the Born approximation by
qs 3 exp[iqx(inr)]
B (xi) 1 exp(iq.xi)X(q,s) a
1 /n 4 T2v r x(inr
exp(iq' xnr) (q, s) C2
If Eq. C2 is substituted in Eq. Cl, one finds
t ( 1 xB(iXj 4
G s (q,q') (zizj) (q',s')s (q, s)
zss n i,j a, /m.m
1 J
2 >3 exp[iqx(jnr)]
Sexp(iq' x. )exp(iqx.) exp(iq'*x) 
1 J 42 1 x(jnr)
S2exp a 2a3 + exp[iqx(inr)] .
Sexp(iqXnr) 42 exp(iqx ) x(in r) exp(iq'xn)
UTnV r
C3
2 4 6 exp[iqx(inr)]exp[iqx(jnr,)]
+ 62v r, x(in ) x(jnr) exp(iq'xn)exp(iqxn, )J
For a random distribution of impurities and for cases in which
qAxnr 1, one would expect the second and third terms in brackets
in Eq. C3 to contribute negligibly due to the presence of random
phase factors. Similar statements would apply to the fourth term
except for cases in which r = r' and = Q'. Thus Eq. C3, becomes
G ,(,') 1 (z.z.) i ( sx (q',s') exp(iq'x)
zss n 9 m i M i %
Ch
2 4 6 exp[iqx(inr)]exp[iqx(jnr)]
S 162v r x(inr) x(jnr) q
The first term in brackets can be shown to produce just the value of
0
G for the isotopically pure crystal. Denoting this by Gzs, and
obtaining its value from the results of Appendix B, one finds
0 2iL
G (q,') wL tW(,s) v(qs) qq 6ss C5
zss L ,q ss
z
Denoting the contribution from the second term in brackets in Eq. Ch
by G ,(q,q'), we have
zss
I + * 1
G ,(q,') =
zss n
i,j
S(z.z )
c,,, J
Q (X.X.
i j
x I exp[iqx(inr)] exp[iqx(jnr)]
r x(inr) x(jn )
246
2 4 6
E w .a
1672 V
c6
qq'
For the case in which cos a > 0, we choose the solutions far to the
left of the sample to evaluate Eq. C6. Then, noting that if
a j ny sl *j
S(x. x ) is negligibly small unless x.x is
Ot i 1 j
less that a few lattice
spaces, we can write
exp[iqx(in )] exp[iqx(jn )] exp(iq'x x.)
where
4. 4*
xixXnr
q' =q
Ix i nr
+ 
Then, upon setting x =
zss (qn
i,m
xj xr
j r
Ixjxnr
* ^4
x.x, Eq. C6 becomes
a1 (
z aB(xm) (qs) h s(q,s')
a,g
xt (q,s)X ( ',s')
