Energy transport in the three-dimensional, harmonic, isotopically disordered crystal

Material Information

Energy transport in the three-dimensional, harmonic, isotopically disordered crystal
Powell, John David, 1945-
Publication Date:
Copyright Date:
Physical Description:
vii, 101 leaves. : illus. ; 28 cm.


Subjects / Keywords:
Amplitude ( jstor )
Approximation ( jstor )
Atoms ( jstor )
Crystals ( jstor )
Geometric planes ( jstor )
Impurities ( jstor )
Low temperature ( jstor )
Phonons ( jstor )
Plane waves ( jstor )
Thermal conductivity ( jstor )
Crystals ( lcsh )
Dissertations, Academic -- Physics -- UF
Energy transfer ( lcsh )
Physics thesis Ph. D
bibliography ( marcgt )
non-fiction ( marcgt )


Thesis -- University of Florida.
Bibliography: leaves 98-100.
Additional Physical Form:
Also available on World Wide Web
General Note:
General Note:

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Resource Identifier:
000577511 ( AlephBibNum )
13986228 ( OCLC )
ADA5206 ( NOTIS )


This item has the following downloads:

Full Text







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.



ACKNOWLEDGMENTS ....................................................... ii

LIST OF FIGURES............................................... v

ABSTRACT................................. ........ .......... vi

CHAPTER I: INTRODUCTION...................................... 1

1. Phenomenological and Related Theory................... 1

2. The One-Dimensional 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
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





UNIFORM DISTRIBUTION.......................... 84


RANDOM DISTRIBUTION........................... 91

LIST OF REFERENCES........................................... 98

BIOGRAPHICAL SKETCH.......................................... 101



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



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 three-dimension-

al dielectric. The model used is the analogue of a one-dimensional

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.


The results are compared with the results obtained in one

dimension. Some discussion of problems remaining to be solved in

three dimensions is given.



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 non-equilibrium 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


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
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
early investigations of Eucken. Eucken found that for most solid

dielectrics the thermal conductivity varied as 1/T at high tempera-
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-


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 so-called phonon-phonon 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 Biermasz5-7 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
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-
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 co-workers6 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-
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 One-Dimensional Model

During the 1960's, several investigators independently

examined the possibility of using one-dimensional 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
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


It was found that the model produced a non-divergent 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 non-divergent. This plausible argument had previously been

stressed by Erdos22 who had performed calculations using an indepen-

dent method and found a non-divergent result for the isotopically

pure three-dimensional crystal.

Nonetheless, the non-divergent 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 one-dimensional problem and

concluded that these were insufficient to produce a length-indepen-

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 one-dimensional 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 three-dimensional model.

3. The Present Calculation

For reasons indicated in the last section, Allen and Ford

have felt obligated to attempt an extension of the one-dimensional

model to three dimensions. In 1969, they succeeded in calculating

the thermal conductivity for the three-dimensional, isotopically pure
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 three-dimensional 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 cross-sectional 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 one-dimensional 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-



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
cross-sectional 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 one-dimensional 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.



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 non-equilibrium 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] II-1
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)] II-2


p(t) = exp (-iHt/Ti) p(0) exp (iHt/T) 11-3

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 11-4
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

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) 11-5

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

<>= -T- f ds f dX Tr [J(0) Jz (t+iXh)p0] AT. 11-6
0 0 Z

In this equation, A represents the cross-sectional 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

PO = e-H 11-7
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). 11-8
z z

If one substitutes Eqs. 11-7 and 11-8 in Eq. 11-6 and

evaluates the trace using Dirac bra and ket notation, one finds

= -A ds f dX <9ej n> exp[-X(E -E )] exp[is(En-E )/M]
0 0 ,n n

exp(-8 E )
x __EY AT 11-9
Z exp(-B Ek z

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/



The integrals over s-and X can easily be performed to yield

= -TrA

Z <|J In> 6 p AT II-11
k,n z

exp(- E)
P = exp(-BEk)


Eq. II-11 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 II-13

* ->
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

1(x) = (ih) [H(x),H]



< In> = 6 n

For a general type Hamiltonian of the form

H [ + V. 1 II-15
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( x-x )(xk- )(ih) [ k ,V(x.)] 11-16
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


V(x) (x-x) Ua( ) U-(x ) II-17
j 2 ,x

where B(xj) is the so-called 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(x-xk) is a coarse-grained


delta function used by Hardy which he has chosen to specifically

represent by the formula

1 2
A(x-xk) = 3/2 3 exp[-x-xk /I2] II-18

In Eq. 11-18, k is distance large microscopically but small macro-

scopically. The function A(x-xk) is normalized such that

dx A(x-x) = 1 11-19

By using the commutation relation for the position and

momentum operators and Eq. 11-17, one can show that for i # j

2 4..
[p (x.),V(x.)] =-ih P(x -xj) pa(xi) U(x.) 11-20

Hence Eq. II-16 becomes

-1 + -+ + + + x P+X
J2(x) = -- A(x-xk)(x -x) ( U (x). 11-21
2 k,j k- (Xk-Xj) 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
representation. For the case of a single atom per unit cell, such

a transformation is accomplished by the well-known transformation

U(x) = Z [ B (x ) a(q,s) + B (x) a (q,s)]
q .2m W(q,s)
q,s Z



Shr(q,s) -+
4. 4 S ) + ** 4 t +
p(x) = -i 2 [ B (x) a(q,s) B (x) a (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 II-23

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 11-24

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) 1-25
p',k m/m

and are orthonormal in the sense that

--* 's' Q
*B (x,) B (x ) = 6- 6 11-26
S-qq ss

If Eqs. 11-22 are used in Eq. 11-21, one obtains

J = I I ) -x)
2 4 L A(x-xk)(xk-x.)
k,j a s s,s' /,T mj
-, -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. 11-27, those terms which do not conserve energy

have been neglected. As evidenced by Eq. II-11, they make no con-

tribution to .

The average value of the heat flux operator is given by

J2 V J dx J2(x) 11-28

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 coarse-grained delta function causes the integrand to

become negligibly small. Furthermore, if one assumes that
-r -. +-
B (xk-x ) is negligibly small when x k-x 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') 11-29
S -+, s,s


(x -x.) q's' *qs
G ,(q,q') = C (x-x) B (x ) B (x.). 11-30
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. II-30 reveals that

--** 3 1
G ,(q,q') = -G, (q',q) II-31

Second, if one performs the sum over j in Eq. 11-30, Eq. II-25 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


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
Eq. 11-30. 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) 11-32

and substituting in Eq. 11-25, one finds

- -
2 q qs c qs
m (q,s) Ca (x ) = a(x-k) CB (x) II-33

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. II-33 in the form

-> ->4
L qs qs
L B(x ,xk) Cg (xk) = 6 LB(xV,xk) C (xk) 11-34
k,B k,B

or in the more abbreviated matrix notation as

4 -
L Cqs = 6 L Cqs II-35


f- 2 4 2
La (x xk) = M 6 6 k B(x -k) II-36


6L (xxk) = EM 2 6 6kn 6 n II-37
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 right-hand side of Eq. 11-35

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 non-zero only when k and k

denote positions of impurities.

Maradudin has shown that the solution to Eq. 11-34 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
C(x) = C (x ) + Wq(x) II-38

qs s
where C is the incident plane wave and W the scattered wave.
The incident wave CO is a solution to Eq. II-34 when no impurities

are present and thus

L CO = 0 11-39

If Eq. II-38 is substituted in Eq. 11-35, one finds

s qs -
L Wqs = 6 L C + 6 L W II-40

Finally if the matrix G = L is defined, Eq. II-40 can be written

in the form

+- qs -
Wqs = G6LC + G6LWqs. II-41

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 2-v 2 +i6 w -v q +i56

1 1 44
+ 6 2 -2 2i exp[iq-x(Z,k)] II-42
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 II-43

In writing Eq. II-42, it has been assumed that the acoustic approxi-

mation is valid, i.e.

v(q,s) = v q 11-44

The integral in Eq. 11-42 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. 11-42

becomes the rather simple expression

4 4 -a exp[iwx(Z,j)/v] 1-45
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
Allen and Powell for the case where vL # VT, though for a simpler


The ultimate aim is to solve Eq. II-40 for Wqs A formal

solution can be obtained by writing

Wqs = [1-G6L]- G6LCO 1I-46

For purposes of calculation we will expand the matrix [1-G6L]-

which yields

Sqs qs
qs = G6LC0 + G6LG6LC + ... 11-47

The first term in Eq. 11-47 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. 11-32 and 11-38. Quite formally, one may write

+ 1/2 qs q
B ) = m. A(q,s) [C0 (x) + Wq 11-48
S 3

where Wqs(x) is given by Eq. 11-46 or 11-47. The normalization

constant A(q,s) can be determined by use of Eq. 11-26.

This section completes the formalism required for a

complete calculation of given by Eq. II-1i.

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. 11-29 into Eq. II-ii

2 -+ G +
=TAh2 G ,(q,q') G 's (q'
4kT2V2 s,s' l,n
4kT V q,q' ,Z
s s ,s

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>
x pO


Eqs. II-10 and 11-24 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

q = q', q" = q"', s = s' s" = s 'i
q~q~q=q ,s=s~s=

Mss ",,(q,q',q'",q ') = 0

, otherwise


c- Tiw(q,s)


In deriving Eqs. 11-51, use has been made of the fact that since

phonons obey Bose statistics

SN (q,s) P = =

e 1


2 ->
N (q,s) pOZ

S =

Use of Eqs. 11-50 and 11-51 in Eq. 11-49 yields

qq" ss"


(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 '


+ 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. 11-31.

If one wishes, the explicit expression for G ,(q,q') given
> ss
by Eq. 11-30 and for B (x ) given by Eqs. II-48 and II-47 may be sub-

stituted into Eq. 11-54. No further simplification results, however,

in the most general case and this will not be done here. Eq. II-54

will be used in Chapter IV where explicit calculations of the current

are undertaken.



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. 11-30. 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

CO (x ) = exp(iqx ) X(q,s) III-1

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


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. II-46

by Wl,then

,qs 23 q e ) *
Wqs (X a 3 exp[iqx(Zn )] exp(iq-xnr) X(q,s) 111-2
W1 (x) = a2 r -r
4wv2 r x(n )

In Eq. III-2 explicit use has been made of Eqs. 11-37, II-44 and

11-45. 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. III-2 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). III-3
1 rv2 r x(2n )

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
of stationary phase :

ff g(x',y') exp[iqf(x',y')]dx' dy'

2 2i g(x ,y ) exp[iqf(xj,y.)]
q cjc j-Y j

o Z


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.

In Eq. 111-4, x. and y. are points in the domain D at which f(x',y')
is stationary and

xI 1
x. ,y
J 3J

, 2
p y

x. ,y.


Yj 3x' y' '
.J .J


The number a. is given by

, for aB >

, for .B. >
3 J

. c .

Y2 a.
3 J

> 0

< 0


,for a".. < .
Y j j

If Eq. III- is applied to Eq. 111-3, one finds after a tedious but

straightforward calculation

-- -(x
W1 ("

i Eoqa -
=2cos a exp(iqx) (qs)
2 cos aj ep1c1 ci


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

Ja 3x'2

(. = 1

a. = -1

a. =-i

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. III-7. 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. III-7 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

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. II-47 can be written

S = -E4 a 3 exp[iqx(Znr)] l("r) -9
7Trv r x(en )

and if one assumes that the first Born term given by Eq. III-7 repre-

sents an adequate solution to the amplitude inside the scattering

region, a calculation identical to the calculation of Eq. III-7 can

be performed with the result

'c -|2 -
+ -i_*aa )+
W2 (x) 21cosa exp(iql-x ) X(q,s) III-10

Clearly this process can be continued to produce

Wq(x ) = Wi (x) = exp(iql-x ) X(q,s)
x [-ix(q) + (-ix(q))2 +(-ix(q)) + ...]


xt Eaqa
x(q) a= 21c III-12
2 cos a|

If we now assume that this expression can be resumed to produce a

result in the form of Eq. II-46, we have

i -ix() -
)= l+ix(q) exp(iql.x) (q,s) III-13

and from Eqs. III-1 and 1I-38,

+qsix(+) +(ixs) III-14

(x) = exp(iqx,) t(q,s) exp(iql ) (q,s) III-14

Eq. III-13 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. III-7 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. III-13 and III-14 do possess the following meritorious proper-


(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 well-known

optical case of electromagnetic waves re-

flected and transmitted by a dielectric

(iv) Finally,they predict values of.the Bs (

which satisfy the orthogonality condition

given by Eq. 11-26.

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.
We therefore conclude that use of Eq. III-14 for the C (x) cannot

cause significant error.

2. Scattering by Many Planes

In Eq. III-14, we have essentially found the reflection

and transmission coefficients for a single plane of impurities. De-

noting these by t and r, we have


r= -.

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
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 = (N-l)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 k-lth plane as well as that part of R which

is reflected from the left-hand side of the k-lth 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


Tk = T ktei + R re2i1
k k-1 k
for k = 2,3...,I-1

= Tkr + Rk+1tei

T = T tei + R re2i1
N N-1 +

RN = rTN

k N



Z=O Z=(k-I)a Z=(N-I)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
transmitted-down to the kth plane. That is, they refer to values
attained just to the left of the kth plane. Only a cross-sectional
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. III-16 can be solved for Rk and Tk in terms of Rk-l

and T and written in the form

T= 1

Rk_ T_-1 r


i- r2 e rp + re ,
T= te r2et T + R
k t Tk-I t - 1

for k = 2,3...N

RN = rTN

or in matrix notation as

T = 1

= M for k = 1,2...N-1
Tk+1 Tk

RN = rTN

In Eq. 111-18, M is the two-by-two matrix


(1-ix)e i



In obtaining Eq. III-19, Eq. III-15 has been used. Finally Eq.

III-18 can be applied successively to obtain the following equations:

T = 1


= MN- 1


RN = rTN



We denote the elements of the matrix MN-1 by

MN-1 a b-21
S= III-21
c d

and have from Eq. III-20,


RN = rTN

RN = aR1 + bT1

TN = cR1 + dT1 .III-22

These equations can be easily solved for R1 and TN yielding

dr b
1 a cr

ad bc
T III-23
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. III-18 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 UN-1-(0)
M 1 U N- (0)
M21 UN-l()

M U (0)
M12 UN-1() )
M22UN-l( )N-2

where the Mi are the elements of the matrix M and 0 is given by

0 = 1/2 (M11 + M22)


The UN-l(0) are Chebyshev polynomials of the second kind and are

given explicitly by the formula

UN-() =

sin [(N+1) cos-1 0]


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=





Applying these results to find the elements of the matrix



given by 111-21, we have

a = (1+ix) e-i UN2(0) UN-3(0)

b = ix e U (0)

c = -ix e UN-2( )

d = (1-ix) ei' UN2(O) UN3(0)
N-2 N-3


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) =


Similarly the total wave R(q) = RI reflected by the N planes is given





-ix UN1 (0)
R(q) = () UN
(l+ix)U (0) el U (0)
N-1 N-2

(l+ix) U (0) ei U N(0)
N-1 N-2

In obtaining Eqs. III-31 and 111-32, Eqs. 111-23, III-27 and III-29

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. III-31 and III-32 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. 111-31, and a wave reflected from the Nth plane whose amplitude

and phase are given by Eq. 111-32.

We now wish to apply the above results to find the ampli-
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

qs +
C (x

=T(q)exp[-iq(N-l)aJcos aj]exp(iq*x ) X(q,s)

for z cos a > 0

C-S(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(N-l)a cos aj]exp(iq*x

x A(q,s) for z > 0 and cos a < 0 III-


In obtaining Eqs. 111-33, we have taken the phase of the incident

wave to be unity at z = 0.

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. 11-30. If

Eqs. III-33 are substituted in Eq. 11-32, one finds that Eq. 11-26

is satisfied provided the normalization constant A(q,s) is given by

A(q,s) III-34

In evaluating Eq. 111-34, we have noted that the sum in Eq. II-26

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. II-32 and III-34 that for x

in the pure regions of the crystal

q -_
B( = Cs (x) III-35

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) III-36

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) III-37

Use of Eqs. III-37, III-35 and III-33 yields the following rather

complicated expressions for Bs ):

For z cos a > 0,

Bs() + = T(q) exp[-iq(N-l)alcos aj]

exp(iq-x) X(q,s)


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)]

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(iq-x ) X(q,3) R(q)cos 2a

- 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(N-l)ajcos all

B (Xi)

1 -. .... 4- .
- [exp(iq-x) X(q,2)

- R(q)cos 2a exp[-2iq(N-l)alcos all

x exp(iql x ) 1(ql,2) + R(q)sin 2a exp[-2iq(N-l)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(N-l)a cos al]


qB (x
B (xi)

exp(iql^'x) k (ql,3)


x exp(iql.x ) 1(ql,3) R(q)sin 2a exp[-2iq(N-l)alcos all

x exp(iqlx ) X(ql,2)] III-38c

4. Discussion

Eqs. III-38 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-
call that we noted that the quantities B (x ) were simply related 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
quantities B5(xj) given by Eqs. III-38.

We begin our discussion of Eqs. 111-38 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. III-38a.

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. III-38b 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
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 right-hand 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. III-2 by the integral

of Eq. III-3. To do so we noted that it was necessary to assume

that the phonon wavelength was much longer than the distance between









Figure 4

Directions of scattered waves. Waves in the lower left-
hand quadrant represent incident waves, those in the upper left-hand
quadrant reflected waves and those in the upper right-hand 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. 111-2 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

X >> 27naa -111-39

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.



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. 1-54 in which we found

2 G [ss(qq)G .ss(q'q') )
z 4kT2V2 > s,s' s(e-)e'-)
s q,q'
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. 11-30. The calculation of
G ,(q,q') is quite tedious and adds little to the physical under-

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. II-25 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 cross-sectional 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. 11-25

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) IV-2

Returning to the evaluation of Eq. IV-i, we note that it

is shown in Appendix B that G (q,q) is an antisymmetric function
of q. Thus the first term in brackets of Eq. IV-1 does not contri-

bute to the current. In addition, it is shown that G ,(q,q') is
non-zero only for q' = q or q' = ql where

ql = + yJ zk

From Eq. B-19, we find

S* 2 12 0 2 2 4
G,(qq')1 = 1 2() I qq)
SI ss 2 z
s,s L


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. IV-3

v() = v(q,s) = V(q)


is the group velocity of phonons having wave vector q, L is the
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. III-31 and 111-32.

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. IV-1

must be converted to integrals using the prescription

SV f dq
+ (21r)


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


6 2 6(q q )
q xq L x
xx x


Using Eqs. IV-3, IV-5, and IV-6, we find that Eq. IV-1


-37Th 2 2 4
= (2-T 2 dq f dq' 2(q) Iv(q) I [ T(q)I \ 6(q-q')


+ 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. IV-7, it has been noted that from the properties of

delta functions


for q' = q



6 (q-q) 6(qy-qy) 6(q+qz)

for q' = q1

The integral over q in Eq. IV-7 is to be done over the first

Brillouin zone.

From Eqs. III-31 and III-32, one finds

IT(q) 2 2 1 2
1 + x2() U _(0)


IR(q) =

x (q) U_1()

1 + x2(q) U 2 ()


6(qx- q) 6(qy-qy) 6(w(q)-(q'))

6(qx- q) 6(qy-qy) 6 ((q)-w(q',))

where x(q) and 0 may be obtained from Eqs. III-12 and III-30, re-

spectively. If one substitutes Eqs. IV-9 into Eq. IV-7 and performs
the integral over q', one finds

2 e
=2 dq 2 (q) IvZ( -I
Z (27) 4kT2 (e-l)2


x 2222 AT
1+ E ga U2 ()
4cos2 N-1
4cos a

where explicit use has been made of Eq. III-12. In the acoustic

approximation, we have

2 ) 2 2
W (q) = v q


Ivz(q)l = vicos a| IV-11

Making these substitutions in Eq. IV-7 and writing the integral in

terms of the dimensionless variable


one finds after considerable algebraic manipulation that

< -6k > 3 _4 e d_
= d2
z a3 0 (e-1)

x f AT
0 202(T/TD)22 2 -
01 + Do2(T/T D) sin2 (Ncos-1 0)

y2 1 02

TD is the Debye temperature for this model which has been defined as

2 tv
T IV-13
D = ka

and y = cos a. Also in obtaining Eq. IV-12, we have made explicit

use of Eq. 111-26. We may note that in terms of the parameter E ,

0 = cos[2(T/TD )y] + sin[2(T/TD)Ey] IV-14

In the limit of low temperatures, we may replace the upper bound on

the C integral in Eq. IV-12 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


1 2 sin2[2(T/TD)y] = 4(T/TD)2 2 IV-15

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. IV-12

-6kv T 3 4 e
= 2a3 T 2 d
rz a D 0 (e-i)

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-


-72kv TT
z 0 2a3 T (4) Tz I-17
Ta D

where r is the Riemann zeta function and we have used the fact that

f n e 2 d = n! (n) IV-18
0 (e-1)2

For large n, C(n) = 1. Eq. IV-18 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(Ncos-1 0) = 4N2(T/TD)2 2y2 .IV-19

If we use Eq. IV-19 in Eq. IV-16, we find

-6kv T '3 0 e_
= 3 d
S 2a3 D 0 (e -1)2

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


x [17 + 60(log(Neo T/TD- y)

where 0 is given by Eq. IV-17 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. IV-16 employing integration in the

complex plane. This procedure will not be discussed but the results

of its application to Eq. IV-16 reveal that for T/TD > 1,

= < E>[2 IV-22
S< z0 + 2 IV-22

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(Ncos-1 0) by its average value of 1/2. If this is done the

integrals in Eq. IV-16 become elementary and the result is

= [l IV-23
z z I /2-

which is in close agreement with IV-22 for small Eg.

Discussion of Eqs. IV-17, IV-21 and IV-22 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 IV-24


where a represented the fraction of impurities in any given plane.

In terms of the variable ( this condition becomes

2(T/T D)
<< 1 IV-25

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. IV-25

may be replaced again approximately by

8 T/TD
<< 1


T 1/2
T << 1 IV-26
T 8

In a typical dielectric, TD is of the order of several hundred

degrees. Hence at low temperatures, Eq. IV-26 can be satisfied

quite easily even for relatively small values of a and we conclude

that use of Eqs. III-38 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. IV-16 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 inter-impurity separation, i.e.

X 2Ta IV-27

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


q q G' ( ,q)G (q ,q')
S, s,s' z s (e l)(e -1)

x -
x S(w(q) W(q'))

which appears in Eq. II-54 does not contribute to in the lowest
non-vanishing order of . In addition, it is shown that


1 = 12L
IG 9,+ 2 2 +
s,(qq')j = 2 0 ~ (q)
zss, -)

x [v(q) -

2 a 3L, (q)Iv (q)I



to lowest non-vanishing order. If we make this substitution in

Eq. 11-54, take the limit as the size of the crystal becomes infinite,

and proceed along steps identical to the development of Eq. IV-16, we


S -3Trh2 + 2-) e
--- 2 dq w (q)
z (2T)4kT2 (eS-1)2


[ sivi- E2a3L0 G'
x IV (q)i- (q) .
z 4Tv 3

If use is made of Eq. IV-13, Eq. IV-30 may be written

= -6kv T 3
z 2 3 T


0 e
0 (e -1)2


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. IV-32 to be valid. We

initially assumed in the calculation that

qa > a'1/3 IV-33

which in terms of the variable E becomes

2(T/TD)E ) '1/3 IV-34

Again the function 8 e /(e -1)2 is sharply peaked attaining its

maximum about 5 = 8. Thus Eq. IV-34 becomes

16(T/T) > 11/3

which implies

T >



Examination of Eq. IV-32 reveals that the second term in brackets may

become quite large at temperatures satisfying Eq. IV-35 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.



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 one-dimensional

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 made-to 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

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 three-dimensional 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


calculation are given by Eqs. IV-17, IV-21 and IV-22. 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. IV-32.

If one wishes to define a thermal conductivity to be as-

sociated with these currents, one may do so by writing

= -K -- V-l
z 3z

where K is the thermal conductivity and 3T/3z the temperature gradient

which is assumed to be adequately given by

= V-2
Dz L0

Using this definition, the thermal conductivity of the isotopically

pure three-dimensional crystal (see Eq. IV-17) becomes

72kvL T3
K a (4 ) V-3
Ta D

For the case in which impurities are uniformly distributed through-

out the sample, one finds from Eqs. IV-21 and IV-22,

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 V-5

for N T/T >> 1. Finally for the case in which the impurities are

randomly distributed throughout the sample, we have from Eq. IV-32,

2 U
K = 13440 NE o'(T/T) (8)
K = K0 1 V-6

It must be strongly emphasized, however, that the values of

the thermal conductivity obtained in Eqs. V-3 through V-6 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. I-1 where

K is independent of the thickness, some explanation of these results

is necessary. The variance of the results with Eq. I-1 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. II-11 represents the first term in this expansion and

higher order terms are assumed to be negligible. Furthermore, the

gradient which appears in Eq. V-7 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. IV-17, IV-21, IV-22 and IV-32 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. IV-17. If one notes that the

specific heat of the phonon gas is given by

C=96k T
c 2 TD i () V-8
Tr a 3 D

Eq. IV-17 may be written in the form

J > CvAT V-9

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. V-9 is, of course, non-divergent

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. IV-21, 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) V-10

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
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
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


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 well-known 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 non-isotropically 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. IV-16. For that case one may also show that the cor-
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. IV-21

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 one-dimensional model where it is found that

the current obeys an equation similar to IV-21 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. IV-22


= <> [V1 + V-11
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

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. V-ll 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
qualitatively Eq. V-ll with a result obtained by Little who has

explicitly treated the problem of Kapitza resistance in solids.

Little considers two dissimilar, semi-infinite 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 one-half. 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. V-ll may then be written in the


144 TkT3 TE
= > 3 C(4) +-- AT' V-13
z 23 h3
vh h

If one further notes that the maximum value of the quantity in

brackets in Eq. V-13 is unity, whereas the maximum value of the quanti-

ty F obtained by Little is one-half, Eqs. V-12 and V-13 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
from the fact that as e - the scattering cross-section becomes

infinite and no current is able to traverse the sample.

Finally it is interesting to discuss Eq. IV-32 with regard

to the one-dimensional 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

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 )










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. IV-32 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. IV-21,

IV-22, and IV-32 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
et al.2 have just recently developed a scheme for treating this

problem in one-dimension. 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 one-dimensional or the three-dimensional model

is the important effect of anharmonic terms in the crystalline poten-

tial which give rise to phonon-phonon 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, phonon-phonon interactions would have to be included.



In this appendix, we seek an approximation to integrals of

the form

I = / f g(x',y') exp[iqf(x',y')] dx' dy' A-1
-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
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 A-2
9x' Dy'

and solving for the required values of x' and y'.
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 (x-x (Y'-Y0

1 +
+ 2 Y0(XI'-X)(Y'-y~0) + ,



a2 f(x' y') 2 f(
0 o0
0x' x0,Y0 0

2 f(x',y')
3x' y' x0,YO

Substitution of Eq. A-3 in Eq. A-i, yields

X1,I I
y 1 I ty,

I = exp[iqf(xQ,y )] g(x01Yg) f

+ ( -y '2 +
80 0

Sep /[ '-
f exp iq/2[a (x'-x02

S ,O'
Y0(x'-x0 (y'-y0)] dx' dy'

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. A-5 can be integrated
by use of integrals which have been tabulated. Performing the



integral, one finds

27i 0 '
I = g(x0 ,y) exp[iqf(xo,yo)] A-6


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)] A-7


where the (x.,y.) represent points at which f(x',y') is stationary.

For the problem at hand, we have from Eq. III-2

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. A-7 generally. This

point has been further discussed in Chapter III.


The purpose of this appendix will be to calculate the

G ,(q,q') which must be used in conjunction with Eq. IV-1 in the
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. 11-26 we found that

Ct i^Z (x' -x) qs )
G zss q,q') = (z.-z.) Bq's'x )B~ sx.) B-1
Si,j a, J/m.m.

and have previously noted in Chapter II that for w(q,s) = w(q',s')

if one performs the sum over j in Eq. B-l the resulting expression

is independent of z.. Therefore, if we further assume that

a (xi-xj) 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. B-l.

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
III-38a, we have

1 -I ()i= -1i-x (T
G ('q') (z.-) Q- ) T(q') T(q)
is nM 1,

x exp[iq(N-l)a( cos al -lcos a' l)]X (q',s')


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(N-l)a( cos al cos a' l)] X (',s') XA(q,s)

x exp[i(q'-q)*x.)] exp(iq*x ) B-3


The sum over the x and y components of the vector x. can be done
immediately by using the fact that

exp[i(q-q')*x.] = n n 6' B-4
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(q-q')x.] = n 6 6 6 B-5
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. IV-2 in the form

= AX (,s) 4 (x ) exp(iq*x ) X (',s') B-6
ca, m

where we have used the fact that

A (q,s) A (q,s') = 6S
cx cxS

By taking the

show that for


s = st

of expression B-6 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)


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

Use of Eqs. B-5, B-8 and B-9 in Eq. B-3 yields

G ss(q,') =

Sw-(q ,s) vz(q,s) IT(q)12 6 6qss




W2 'qs) 6 sS,

where L is the total length of the crystal along the z axis (to be
made infinite in the final formulae). Eq. B-10 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
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
I ->
for B '(x ) and B (x.) far to the left of the scattering region,
i.e. z., z < 0. We can write from Eq. III-38b the general expres-
1 j

B (xj) = n [exp(-iq-xj) (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] B-12

where the f's are to be obtained from III-38b. Similarly, we have

from Eq. III-38a for this case

q's' 1
Ba (x.) [T(q')exp[-iq(N-l)alcos a' ]

x exp(iq'-x.) a (q',s')


By proceeding in a manner identical to the evaluation of B-1, one


G (sq,q') = w (,s)v (q,s)T(q) exp[-iq(N-l)alcos al]
zss L 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()


and that

T(q) = T(ql)

as can be seen from Eq. 111-31.

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

G z,(sq') -L (q,s)v (q,s)T(q)exp[-iq(N-l)alcos al]
zss L 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 .


We see immediately from Eqs. B-ll, B-14, and B-15 that

G (q,q') is non-zero only for q' = q or q' = q Furthermore,
zss 1
-9- ->-
we see from Eq. B-1, that G (q,q) is an antisymmetric function of

q since

v () = v (-q


Consequently, the first term in brackets in Eq. IV-1 does not con-

tribute and we need only evaluate

IG s (q,q)1


zss 12


From Eq. B-11, we find

IG12 12L 2 q 14
I 2 t 2-+ 2 L
S(q,) -2 (q) v (q) IT(q)|
ss' L


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. B-14 and B-15 from Eq. III-39b, one finds

Gzss 5 q


12L0 2
= 2 (q)

2 )

IT(q) 12 R(q) 12

or in general

12Lo 2


x 6) + IT(q)I2 jR()|2




Vz( )

[ T(q)I

q 'q




We shall obtain an approximation for the quantity

G ,(q,q') in this appendix for the random distribution. The
calculation will involve the use of the Born approximation for the

solutions of the qs (x,) and a random phase approximation. From

Eq. 11-26,

x (x.-x.) q's' *qs
G 9s q,q') (z.-z.) B B (x ) B (.) C-1
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) C-2

If Eq. C-2 is substituted in Eq. C-l, one finds

-t -( 1 xB(i-Xj 4
G s (q,q') (zi-zj) (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


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. C-3 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. C-3, becomes

G ,(,') 1 (z.-z.) i ( sx (q',s') exp(iq'x)
zss n 9 m i M i %


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
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 C-5
zss L ,q ss

Denoting the contribution from the second term in brackets in Eq. C-h

by G ,(q,q'), we have

I -+ -* 1
G ,(q,') =
zss n


S(z.-z )
c,,, J

Q (X.-X.

i j

x I exp[iqx(inr)] exp[-iqx(jnr)]
r x(inr) x(jn )

2 4 6
E w .a
1672 V



For the case in which cos a > 0, we choose the solutions far to the

left of the sample to evaluate Eq. C-6. 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.)


4. -4*
q' =q
Ix i- nr

-+ -
Then, upon setting x =

zss (qn

xj -xr
j r

*- ^-4-
x.-x, Eq. C-6 becomes
a1 (

z aB(xm) (qs) h s(q,s')

xt (q,s)X ( ',s')