T HE UNIVERS IT Y OF MI CHI GAN COLLEGE OF ENGINEERING Department of Nuclear Engineering Technical Report DAMPING THEORY AND ITS APPLICATION TO THE INTERPRETATION OF SLOW NEUTRON SCATTERING EXPERIMENTS A, Ziya Akcasu ORA Project 04856 under contract with~ NATIONAL SCIENCE FOUNDATION GRANT NOo G-20037 WASHINGTON, D Co administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR July 1963

This report was also a dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in The University of Michigan, 1963

TABLE OF CONTENTS Page ABSTRACT v CHAPTER I. INTRODUCTION 1 IIo DAMPING THEORY 7 A, Basic Theory 8 B. Evaluation of the Inversion Integrals 18 Co Transition Probability Per Unit Time 19 III. GENERAL FORMULATION OF NEUTRON SCATTERING BY MACROSCOPIC SYSTEMS 26 IV. LATTICE DYNAMICS 32 V. DIFFERENTIAL CROSS SECTION FOR SCATTERING OF SLOW NEUTRONS BY AN ANHARMONIC CRYSTAL 40 A. Cross Section Formula 40 B. Computation of Matrix Elements of exp(io.x) 42 C. Width and Shift Formulas 45 D. Discussion of the Cross Section Formula 54 VI. SUMMARY AND CONCLUSIONS 64 APPENDIX: SIMPLIFICATION OF CUBIC ANHARMONIC POTENTIAL 67 REFERENCES 69 iii

ABSTRACT The damping theory, which includes the decay of quantum states, is developed in terms of the projection operators. The concept of probability that a transition from a given initial state will be towards a specified final state is introduced, and the probability is calculated by using the damping theoryo.When the decay of the final state is neglected, this probability reduces to the conditional probability, introduced by Heitler, that the system will be found in a specified final state after a transition from the given initial state has definitely taken placeo With this new probability concept, an approximate expression is obtained for the transition probability per unit time from the initial to the final stateo This expression is the starting point of the present line shape theory, which is applicable to both neutron and photon spectrao In this work, only the scattering of slow neutrons by an anharmonic crystal is consideredo However, the differential scattering cross section formula is obtained first for an arbitrary macroscopic medium, and then applied to a crystal. The cross section formula in the harmonic approximation follows from the present cross section formula when the width and shift of lines are neglected. The novel feature of this derivation is that it does not use Bloch's theorem. Explicit formulas for the width and the shift of the observed peaks in the energy spectrum of the inelastically scattered slow neutrons by crystals are obtained, and compared in the case of a single-phonon event to those by Maradudin and Fein. It'is found that the shift formulas agree exactly at all temperatures, whereas the width formulas agree only in the zero temperature limit. The discrepancy in the width formulas at finite temperatures is discussed in terms of the concept of phonon lifetime. The application of the present line shape theory to the interpretation of the optical experiments is also discussed. v

___ ______ _ ___ _____I_

CHAPTER I INTRODUCTION In recent years the scattering of slow neutrons by macroscopic systems has received increasing attention as a research tool, providing information about the dynamical conditions of the scattering medium, in addition to the determination of nuclear scattering lengths of the constituent atoms. For example, the scattering of slow neutrons by crystals provides a useful complement to other techniques for the determination of the energy levels of the crystal. The general theory of neutron scattering by macroscopic systems, in terms of the time-relaxed space correlation function, provides a convenient computational tool when the energy eigenstates of the scatterer are known. However, except for some special physical systems, eog., an ideal gas, the eigenfunctions of the Hamiltonian Hs of the scatterer are not available. In such cases the evaluation of the thermal averages appearing in the cross section formula have to be performed in a different complete set generated by a certain Hamiltonian The latter is chosen as a part of H, viz., HS = + H', (1.1) for which the eigenvalue problem ( s 2 - E) |n> = O (1.2) 1

2 can be solved. In some applications, it may be a good approximation to replace H by j-s, thereby ignoring H' completely, in order to explain the dominant aspects of a certain scattering experiment. However, interpretation of the details of the recent neutron scattering experiments requires a more refined analysis than that obtained by ignoring HI completely. As an example,. consider again the scattering of slow neutrons by crystals, which will be our main concern in the subsequent chapters. Here, H' is taken to be the cubic and the higher order terms in the Taylor expansion of the crystal's potential energy. The harmonic approximation which corresponds to ignoring H' predicts sharp lines in the energy spectrum of neutrons scattered inelastically by the crystal, and fails to yield any information about the structure of the observed peakso Since the width of these peaks can be measured by present-day experimental techniques, the need for a refined theory of neutron scattering by crystals, one which is capable of predicting the observed width and shift, is apparent. Valuable information about the nature of anharmonic forces in crystals can be gained by comparing the experimental and computed widths once a reliable theory is available. The study of neutron line shape may also prove to be a very suitable microscopic probe in the exploration of intermolecular forces in liquidso Several theoretical studies of'the neutron scattering by anharmonic 2-6 crystals have appeared in the literature in recent years. The most recent one, by Maradudin and Fein,6 contains a brief survey of the preceding works as well as a comparison between theoretical and experi

5 mental results for the width of one phonon line in lead. The agreement is reported to be only in the order of magnitude. As the authors point out, however, this discrepancy may be due either to the computational techniques used to represent 6 function, or to the simple crystal model used to describe lead, rather than to the approximation made in the derivations. In this regard, an independent formulation of the theory of neutron line shape would be very desirable, The potential of the study of the optical line shape as a noninterfering diagnostic tool for investigation of the physical properties of the surrounding of an emitting atom has long been recognized and used in astronomical observationso In the past fifteen years, the optical line width techniques have found application also in high temperature plasma experiments as a probe for measuring temperatures and ion concentrations within a plasma, where the use of material probes is ruled out because of the high temperatures involved. Consequently, the optical line shape theory has enjoyed rapid progress in recent years. Accuracies better than 20% in determining the ion density with the line width technique have been reported.* In spite of this numerical success, however, the existing line shape theory contains some intuitive arguments, such as the folding of Doppler and pressure effects, as well as some gaps. It treats the effect of ions in the quasi-static limit by assuming the ions to be at rest, whereas it treats electron. effects *A list of the original papers and several excellent review articles can be found in Ref, 7,

4 with the impact approximation. It does not offer any good treatment for the intermediate region, where neither of these two limiting approximations is valid (see p. 504 of Ref. 7). The starting point of the existing theory is the quantum mechanical Fourier integral formula which expresses the optical spectrum as the Fourier transform of the autocorrelation function of the time-dependent dipole operator. In the present work, an attempt has been made to develop a line shape theory that can be applied to both optical and neutron spectra. This theory involves approximations of a general nature in contrast to the previous optical and neutron line shape theories, which involve approximations appropriate to specific physical systems, e.g., plasma and crystals, and therefore have a limited range of applicability. The present theory provides a general line shape formula for arbitrary media, i.e., the surrounding of the emitting atom in the photon case and the macroscopic scatterer in the neutron case, which reduces the line shape calculations immediately to a computational form. The nature of the physical system under consideration first enters at this computational stageo The present formulation of line shape is based on an entirely different approach than the correlation formalism used by the existing line shape theories for both neutrons and photons. This approach uses the damping theory originally due to Heitler,8,9 and is essentially a perturbation approach which includes the decay with time of the quantum states of the system. This approach was chosen for the present study

5 because, as will be demonstrated by two specific applications, it provides a computational framework which is more systematic or more interpretable physically than any used in similar applications, and which can be applied directly to any line shape problem without approximations of the physical nature appropriate to a particular system. In the case of the study of optical line shape, the present theory illuminates the manner in which the various simultaneous contributions to a line shape, i.e., Doppler effect, natural broadening, and pressure broadening, arise and are combined; and it provides a method of computation which may improve upon existing approximations, ioeo, quasi-static and impact approximation, in the intermediate range mentioned aboveo In the case of the study of line shape in the energy spectrum of inelastically scattered neutrons by an anharmonic crystal, the present approach yields, in a simpler and more interpretable way, formulas for the width and shift of the lines in the zero temperature limit which are identical to those presented by Maradudin and Feino. However, the present line width formula at finite temperatures differs from theirs, and predicts a larger width. This discrepancy, if it does not arise from an arithmetical error, may improve the order of magnitude agreement reported by Maradudin and Fein, In Chapter II, the damping theory is developed and extended by introducing the concept of probability that the transition from an initial state will be toward a specified final state; this is valid.even if the decay of the final state is not negligible, as is the case in

6 crystals. The remaining chapters contain the application of the basic formula derived in Chapter II for the transition probability per unit time, to the scattering of slow neutrons by an anharmonic crystal. The application of the present theory to the study of optical line shape is not included in this work because it has been presented as a separate report.1 It is hoped that the present work will contribute to the theoretical understanding of line broadening and line shift phenomena in general. Such understanding is essential for the success of the line shape studies as a convenient probe to explore intermolecular forces in a macroscopic medium.

CHAPTER II DAMPING THEORY The interpretation of almost any experiment can be reduced in quantum mechanics, without introducing any serious approximation, to evluation of the quantity Wmn Umn(t) 12/t, (2o1) where the numerator is the probability of finding the system under consideration in the state [m > at time t, knowing that it was in the state In > at the initial time t = Oo The purpose of this chapter is to derive an expression for Wmn which is suitable to the study of line shape problems in general. Such an expression can be obtained by using the damping theory, which, in contrast to the conventional perturbation theory, takes into account the decay in time of the quantum states. Although the first attempt to include the finite lifetime of the states in the perturbation theory was made by Weisskopf and Wigner,ll the first systematic development of the damping theory is due to Heitler.8'9 The present derivation differs from Heitler's original derivation in two respects: First, it is developed in terms of the projection operator, and second, it employs a different iteration procedure. Although the results are essentially the same as those obtained by Heitler, the present approach has several appealing features, It is an extention of the quantum treatment of the evolution of a decaying state presented by 7

8 A. Messiah.12 A. BASIC THEORY The temporal development of a quantum mechanical system is determined by the "time-evolution" operator U(t): It > = U(t) >, where 10 > and It > are the state vectors at t = 0 and at time t. When the Hamiltonian H of the system is not an explicit function of time, U(t) is given by U(t) = e-itH (2.2) where the units are chosen so that ~ = 1. Assume that H can be split into two parts as H = H + V (2.3) in such a way that the eigenvalue problem (Ho-En) In > = 0 (2.4) can be solved. The eigenvectors In > are assumed to be complete and orthonormal, and thus provide a basis for the physical problem under consideration. The eigenvalues En, which may be discrete or continuous, are allowed to be degenerate. We want to calculate the matrix elements Unn(t) and Umn(t) for t > 0O However, we shall first consider the

9 operators PnU(t)Pn and PmU(t)Pn, where Pn and Pm are projection operators on the subspaces fn and m spanned by the eigenvectors belonging to the eigenvalues En and Em respectively. The projection on the complementary subspace of n will be denoted by Qn. The following relations are immediate consequences of the preceding definitions: Pn = nv > < nvl, Pn = 1, (2.5) v n PnHo = EnPn = oPn (2.6) Pn + Qn = 1, QnHo = HoQn (2,7) In (2.5), the symbol v indicates the multiplicity of the eigenvalue En. 12 13 Similar relations hold also for Pm. We may recall'13 that any hermitian operator satisfying the relation P2 = P (2.8) is a projection operator. It will be apparent at the end that Pn and Pm can be defined as projections on some subspaces of the subspaces en and i respectively. In particular, they can be chosen as the projectors on the initial and final state vectors. In this case, and also when the eigenvalues are non-degenerate, the present approach reduces to that of Heitler's except for the iteration procedure. We now introduce the resolvent of H, i.e,, G(z) - (2.9)

10 where z is a complex number. Since the eigenvalues of H are real and positive, the singularities of G(z) all lie on the positive real axis. The discrete and continuous portions of the spectrum of H give rise to simple poles and branch cuts on the real axis. The evolution operator U(t) can be expressed in terms of G( z) as Woo+itz U(t) = 2 +idz G(z)e, > 0 (2.10) i -00+iE or i1 +oo —,2 where G( x) lim 1 = PP 1 i rb(x-H). (2.12) -+0 x-H~iE x-H In the last expression PP denotes the principal part. It may be noted in passing that the contribution of the second term in (2011) vanishes for t > 0, since G(z) is analytic in the entire lower half plane as well as upper half plane. Our problem is now reduced to calculating -Pnn( zP) P G()Pn, n(Z) PmG(z)Pn ~ (2.14) Following Messiah,* we introduce *See p, 994, Vol, 2 (second printing in English) of Refo 12,

11 H = H' + H" where H' and H" are defined by H' I PnHPn + Qn HQn = Ho + PnVPn + QnVQn (2.14a) H" PnHQn + QnHPn = PnVQn + QnVPn (2.14b) The following relations can be verified by using the foregoing definitions: [PnH'I = 0, [Qn,H'] = 0, (2.15a) PnH" = Hn H"', QnH = H"P PnH"Pn = QnHQ = 0 (2l15b) The following operator identities will be needed: 1 1 1 1 1 A- A A A-B A-B A + (2.16) ~AP 1p i = P 1.[R,P] = 0, (2.17) z-R z-R z-PR= (27) where A and B are two arbitrary operators, P is any projection operator, i.eo, hermitian and P2 = P, and R is any operator commuting with P. The proof of these identities is straightforward. Using (2.16), one can write G(z) from (2.9) as follows: G( z) =' H" G( z) ( 2. 18a) 1 1 H" 1 118b) - z-H' z-H. z-H +..H " HIG(z) (2.l8b)

12 We can now compute n (z) by operating on both sides of (2.18b) by Pn, noting that the second term vanishes by (21.5b), and by using (2.17). The result is'frn(Z) 2= (2.19) /hn z-En-rnn(z) where the symbol rnn(z) is defined by rnn(z) PIn[V + V(z-QnHQn)V]Pn o (2.20) We now attempt to compute n(z). It proves more convenient to compute first QnG()Pn, and then to obtain mn(z) by multiplying the latter from the left by Pmo Thus, using (2,18a) for G(z), one finds: m(z) = Pm( ZQnHQn) oQnV n( z) (2.21) Noting that Qn= P cltn one rewrites (2.21) as follows: ~mn(Z) = Pm(z QnHQn) PmV nn( z) + 7 Pm( zQnHQn) "'PSV n( z) H^,"X7lm ~(2.22) It is noted that all the foregoing results are exact. However, they still contain the operator (z'QnHQn), which can not be calculated

13 exactly.* At this point, one resorts to an iteration procedure which consists of applying the foregoing analysis to (z-QnHQn)'1 by defining ^ — =QnHQn, -i + tp (2.23a) o0 = QnHoQn, Q= QnVQn, (2.2b) and noting that PmPo = PHoQn = EmPm = oPm. (2.24) Thus, Pm(z-)l'Pm in (2.22), which can be labelled as -mm(z) according to (2.13), can be readily obtained from (2.19) by replacing Ho and V by an and A a(mm Z) z-Em-r m(z)' (2.25) where rPm(z) - Pm[V + V(z-QmQnHQnQm)- V]Pm (2.26) The calculation of Pm( z-QnHQn) -1P appearing in (2.22) requires closer attention, because of the order of Pm and Po, We note that in this operator, Pm denotes the initial subspace whereas in (2.13) it denotes the final subspace. Therefore, the order of the initial and final sub*Messiah approximates this operator by (z-Ho)-1 by ignoring the QnVQn in QnHQn = Ho + QnVQn [cf. (2.14a)], and thus obtains a closed form for Smn( z) from (2.21). But this approximation excludes automatically the width and shift of the final state Im > as will be apparent presently [cf. (2.28) ].

14 spaces in m( Z) Pm( z-QnHQn)' 1 is reversed as compared with that in n( Z) To avoid any confusion, one can start with the second operator identity in (2o16) to compute t(z). The result can be easily seen to be imj(z) )= mm(z) iQm( Z-Qm Qm) -P e (2.27) The result of the first iteration can now be obtained by substituting (2.26) and (2.27) into (2.22). The resulting formula will contain operator (z-QmtPQm)1, which can be iterated once more. (We shall give the result after the second iteration.) Omitting the arguments (z) for typographical simplicity, one finds~ = V + 7 V V + o.. (2.28) ~mn.mm / nn A mm ft nn where lt, 1( — 2 (2.29) 1 L" 1 1z1-E rfT r" = P I[V + V(z-QmQnHQnQmQVL) 1V]PL. (2.30) It is observed that (2,28) represents formally an expansion in powers of Vo The first term corresponds to the direct transitions from the initial subspace to the final subspace; the second term represents the transitions via an intermediate subspace; and so on, In the present analysis of line shape, we shall consider only direct transitions, and hence only the first term in (2,28),

15 In evaluating the operators rn, rm, r1', etc., one may approximate the operators (z-QnHQn), (z-QmQnHQnQm) (z-QQmQnHQnQmQ), etc., appearing in their expressions, by replacing H by Ho. Thus, rnn and Fm become nn(z) = Pn[V + VQn(-Ho) V]Pn, (2.31) rFm(z) = Pm[V + VQmQn(z-Ho)'-V]Po m (2.32) The approximate expressions for r", etc., can be obtained in a similar way, When the eigenvalues En, Em, E~, etc., are all non-degenerate, or when the operators Pn, Pm, P., etc., are defined as the projectors on the individual eigenvectors, the operators in = PnGPn, m = PGPn, etc., can be replaced by their matrix elements with respect to the eigenvectors In >, Im >, etc.: inn(z) = [z-En-rnn( z)]-, (z) = [z-Em-r ( z) ]1 =mn( ) =' (z~)vm (z). (2.35) methodl0 of iteration except for the fact that the quantity rm in (2.34), which will be interpreted below as the width and shift function of the final state, excludes the transition back to the initial state

16 as well as the transition into itself. However, this difference does not seem to be of any significance in the applications we shall be concerned with, When the eigenvalues are degenerate, the above formulas, ie,, (2.33), (2.34), and (2,35), can still be used, as we shall do in the subsequent analysis. However, when the damping theory is formulated in terms of the projection operators and the projection operators are defined as the projectors on the entire subspaces spanned by degenerate eigenvectors, there are new possibilities in the case of degenerate eigenvalues, Although we are not going to explore or make use of these possibilities, it may still be worthwhile to discuss briefly some of the relevant aspects. Let the degenerate states be labelled as In>= IEnce> and im > = iEm > o The matrix element of PmG(z)Pn follows from (2.28) as < EmIG(z) IEna > = (2036) < EmplImm(z) IE-ml> < Ems IviEnC'> < En' lnn(z) IEn > - Suppose that the operators jm( z) and 1( z) can be diagonalized with respect to the degenerate states, Then, (2~36) reduces to < EmlG(z) IEn > = < Elt z) |E > < E|v >< E |n( z) |E > (2.37)

17 The first and the third terms on the right hand side of this equation are simply the eigenvalues of the non-hermitian operators n( z) and inn(z) respectively. It follows from (2.25) and (2.19) that one has to diagonalize the operators rFm(z) and rnn(z) to obtain the foregoing eigenvalues. To carry the discussion, consider rnn(z), which is defined by (2.20). We are interested in the behavior of this operator in the neighborhood of branch cuts [cf. (2.12) ]: nn(x~ie) - Sn(x); Yn(x), (2.58) where Sn(x) - Pn[V + VPP(x-QnHQn)-V]Pn, (2.39) yn(x) 2tPnVb(x-HQHQn)VPn (2.40) Note that both Sn(x) and yn(x) are hermitian operators. The nonhermitian operator nn(x+) can be diagonalized if it is normal, i.e., if Sn and 7n commute, in the case of degeneracies of finite order. In some applications, such as those involving photon emission, the subspace cn can be written as the tensor product of two spaces An' and In",, where,n is a finite dimensional space and in, is a continuum. The former represents the degeneracies of the discrete energy levels of the emitting atom and the latter corresponds to the continuous energy spectrum of the surrounding medium. In such cases, the non-hermitian operator rnn(x~) can be diagonalized1 in the subspace n',. The ad1See the footnote on p. 996 of Ref. 12.

18 vantage of this diagonalization procedure can be appreciated only when one actually performs the summations and averages on the final and initial states in the interpretation of a given experiment. Some of these points are discussed in Refo 10; we shall not dwell on the matter because in the present application, ioe, scattering of neutrons by crystals, the formulas (2o33) (2o34), and (2o35) prove to be sufficient. Bo EVALUATION OF THE INVERSION INTEGRALS We are now in a position to compute the matrix elements of the evolution operator U(t) viZo Y Unn(t) and Umn(t), by substituting nn (z) from (2 33) and - (z) from (2o35) into either (2o10) or (2.l11)o First consider the diagonal element Unn(t) o By combining (2.38), (2.33), and (2.11) one obtains the following integration: 1 dx- ixt (x) =nn(t) Fn(t) = j x.-ixt e.....x.~ (t F ) -[x-En-Sn(x)]2 + Yn 7(x) (2.41) We have introduced the symbol Fn(t) for a later use. We shall evaluate this integral approximately by treating Sn(x) and yn(x) as a constant~ This is justified because Sn(x) and yn(x) are slowly varying in the vicinity of the point x = En where the integrand attains its maximum, and because the dominant contribution to the integral comes from this region. Taking the constant values of Sn(x) and yn(x) as Sn(En) and yn(En), one finds from (2,41) the following: Fn(t) x exp[- 2 7nt] exp[-iEnt ] (242)

19 where - En + Sn, S7n - Sn n(En) (2.43) Now consider the off-diagonal matrix element of U(t), which can be written as Un(t) = 1 J dze) Vmn n(Z). (2.44a) 2 ti -oo0+i This inversion integral can be expressed as the convolution of Fn(t) and Fm(t), which are approximately given by (2.42): t Umn(t) Vmn J Fm(t-T) Fn(T) d (2.44b) 0 o The latter can be evaluated immediately. Since we are interested in the probabilities rather than in the matrix elements, only IUmn(t) 12 will be given: 2 Iu(t) 2 - — 2IVmnI2 -e mt e nt e t -i(En-Em)t UEmnJ( ) tn-E )2+( 1/4)(72-45n)2 (2.45) C. TRANSITION PROBABILITY PER UNIT TIME This section is devoted to the discussion of the behavior of IUmn(t) 2 in different time intervals by using (2.45). 1, ynt ~ 1, ymt ~ 1, t > l/IEn-Em| For such an interval of time to exist, lEn-Em| must be much greater

20 than both yn and ym. The behavior of Umn(t) i in this interval can be obtained by setting 7n = 7m = 0 and taking the limit t - oo in (2.45). Noting also that -iXt 2 lim I1e I = 2s( x), t-o x2t one obtains the well known first order perturbation result, i.eo, Iumn(t) 2 = 2ertIVmnl 2(k-k) (2.46) It is noted that in this time interval, the meaningful concept appears to be the transition probability per unit time, i.e., IUmn(t) 12/t, rather than the probability of finding the system in the state Im > at time t. 20 ynt > 1, 7 t < 1 (or vice versa) This time interval exists when 7n >> m or vice versao In this interval, (2,45) reduces to iUmnl = n, (2.47) (Em-En -Sn) 2+( 1/4) where Sm is defined by Snm - S - Sm' (2.48) Note that in obtaining (2,47), 7m has been neglected as compared with 7n in the denominator of (2045) so that the equation will be consistent with the condition Yn >> 7mo

21 It may be pointed out here that (2.47) rigorously follows4 from (2.44b) if rmm(z) = 0. This can be seen as follows: When rm(z) = 0, Fm(t) can be evaluated rigorously [cf. (2.41) ] as Fm(t) = iEmt t>0 Substituting this into (2.44b) gives t i2 Umn(t) 2 = / dT e Fn(T) (2.49) o To obtain the behavior of I|Un(t) 12 for large times, i.e., ynt >> 1, one has to consider the limit of (2.49) as t + oo In this limit, the integral involved in (2.49) defines the Laplace transform of Fn(T), which is n n(z) given by (2.33). Thus, replacing z by Em in (z) immediately gives (2.47). However there is a slight difference between (2.47) and the result obtained with the rigorous derivation. The Sn and 7n are evaluated at x = En in (2.47), whereas they are evaluated at x = Em in the latter case. But since the difference between En and Em is of the order of a line width, and furthermore since Sn(x) and 7n(x) are slowly varying functions of x, no distinction will be made between 7n(En) and yn(rEm) or between Sn(En) and Sn(Em). In conclusion, one may state that when the widths of both initial and final states are much smaller than the transition frequency (En-Em), the transition probability per unit time appears to be a meaningful concept. On the other hand, when yn > 7y, the meaningful concept is the conditional probability of finding the system in the state Im > after it

22 is certain that a transition from the initial state In > to any other state has taken place (ynt > 1). 3. Ym Yn When yn and ym are comparable, there iS no time interval in which IUmn(t) 12/t or IUmn(t)12 can be approximated by a time-independent quantity. The conditional probability defined above becomes time-dependent because the decay of the final state cannot be ignored when 7m is of the same order of magnitude as n'. However, there is a probability concept which can be defined and computed unambigiously in all the above three cases. This is the probability mn that a transition from the initial state In > will be into the final state Im >. To compute, one first observes that the approximation of constant 7n and 7m implies that the states decay exponentially at all times. This can be seen from (2.42), which can also be written as Unn(t) 2 = e7nt. (2.50) Using this information and the probability?mn' one can calculate IUmn(t) 2 in an alternative way: iur2 n(t) 1 e7nt' Yndt' n e-m(t' ) (2.51) Umn(t) 2 = 7n ( - emnt n 2 7n-Tm

23 In (2.51), the factor exp(-ynt')yndt' is the probability that a transition will occur in dt' about t', mn is the probability that this transition will be into Im >, and the last factor is the probability that the system will remain in Im > in the time interval (t-t')o Note that mn may be a function of t'o In obtaining (2o52), it is explicitly assumed, as an additional approximation, that m is independent of time. We shall now compare (2.52) with (2.45), which is obtained directly by the damping theory. For the sake of definiteness we assurge that 7n > 7m. When (yn-Ym)t > 1, the foregoing comparison yields the following expression for Imn~ V0 = | 12 7mn (2-53) Jmn Yn (n- n)2 n where we have introduced Ymn = l7n-7ml/2 (2.54) It is noted that mn reduces to (2o47) when 7m << n The expression (2.53) of mn suggests that n mn be interpreted as the transition probability per unit time, Wmn, from In > into Im >, because 7n is the probability of decay of the initial state per unit time [cf. (2.50)]: W' 12 =Ymn (2.55) Wimn ~/n~kn =2ImnI (_m —En)2 \2+mn It is again observed that Wmn reduces to the conventional form (2.41) when Ymn ~ (Em-En). The crux of the matter appears to be that the en

24 ergy conserving delta function in the conventional form of the transition probability per unit time is replaced in (2.55) by the last factor, which is a peaked function. Equation (2055) is the starting point in the present theory of line shape. The subsequent chapters will illustrate its application to the study of line shape in neutron spectrumo In view of its importance in the present work, we shall give an alternative derivation for (2.55) with the hope that it may reflect more clearly the nature of the approximation inherent in (2.55). We start with Wmn = |Umn(t) 12/t, which appears directly in the interpretation of many experiments in quantum mechanics. Since U(t) is unitary, one has I Umn(t) 2 = 1. Using m this relation, and multiplying the numerator and the denominator by 2 IUnn-UmmI, one can write Wmn as follows: I Umn 2 lUntn12 7 IUn'm12 (1.6 mWn = En2 IL t- L - (2L56) lunnUmml n'#n n'tm We have omitted the arguments of Umn(t) in (2o56); the equation is of course exact, We introduce the first approximation by evaluating the first term, in the time interval IYn-Ymlt > 1, using (2,45) and (2.50). This approximation implies also retaining only the first term in the expansion of n- [cf. (2o28) 1, and assuming that 7n and m are slowly varying. To approximate the second factor in (2o56), we first let yn and 7m tend to zero and then let IEn-Elt be large. In this approximation, the second factor is equal to ]n-Ym|l, as can be seen from (2.44a)

25 and (2.40). With the foregoing approximations, (2.56) reduces to (2.55). The justification and the implication of these approximations are beyond our reach. Therefore, the validity of the basic formula (2.55) will have to be demonstrated by its success in various applications.

CHAPTER III GENERAL FORMULATION OF NEUTRON SCATTER: BY MACROSCOPIC SYSTEMS The expected number of neutrons having a momentum ~kf and spin Tf in a system composed of neutrons and the scattering medium is given by X(f, Tf) = Tr[p(kf,Tf)D], (3.1) where p is the number operator for neutrons of the specified kind, and D is the usual density operator. The Hamiltonian H of the system consists of the kinetic energy of neutrons, Hn, the Hamiltonian Hs of the scatterer, and the interaction potential Vn between the two, viz., H = ES + tf + Vn (3.2) The interaction between neutrons is neglected. Let H~ and V denote the unperturbed Hamiltonian and the perturbation respectively, and let them be chosen as H~ = js + Hn, (3.3a) V = Vn + H', (353b) where S and H' are defined by (1.1). By computing the trace in (3.1) in a representation I{ n > } where H~, p, and the projection of the spin of neutron are diagonal, one can express10,l5 the rate of change in the 26

27 neutron number of the given kind as follows: X ~ Pni(lf-i)iWnfni (3.4) ninf In this equation,!ni > and |nf > are the initial and final states of the system; Pni, If, and qi are the diagonal elements of the density matrix and the number operator, viz., Pni - Dnini r'i = Pnini, Qf Pnfnf and, finally, Wnfni is defined by Wnfni - (l/t) | < nflexp(-iHt/i) Ini >12 (5) The eigenstates Ini > and Inf > are labelled in detail as follows: Ini > = 1hi >ITi > Ii >, Inf > = Ilf >jTf > If > where ITi > and JTf > denote the spin states of neutrons, and Ii > and if > are the eigenstates of R including the spin states of the scatterer, i.e., (f-E) > =, ( = i,f). (3.6) The occupation numbers for neutrons are either zero or one, Hence, (qf-1i) in (3.4) is either +1 or -1. The terms with positive sign correspond to the scattering of neutrons into the momentum state likf and

28 and spin state Tf, whereas those with negative sign correspond to the scattering out of the state (kf Tf). In a scattering experiment with a monoenergetic beam of neutrons, the initial occupation numbers are given by i(kT) = 1, for k = ki T = Ti, = 0, otherwise. Henceforth, the spatial part of neutron states will be denoted by |k >. The rate of change in the expected number of neutrons having a momentum in E3d3kf about'fkf can be computed from (3.4) as mj(/f )d3kf= (1/2) Pi Wnfni (3-7) k6f d3kf iTi f, f where Pni has been approximated by PiPTiPki. Furthermore, Pki has been replaced by unity since the neutrons are prepared in 1k. > initially, and PTi has been set equal to (1/2), assuming that the initial beam is unpolarized. The differential cross section per particle for scattering of neutrons from ki into dG about a = (kf/kf); dg about f - (2kf/2m) is obtained from (357) by dividing both sides by the initial neutron current (,Tki/m2) as well as by the number of particles N: _ _2 kf m2 n,_) - i ) k v 2 Pk Wmnni. (3.8) (i,f, rf,'Ti) In (3.8), 2 is the volume of normalization and m is the neutron mass.

29 Except for writing Pni as the product of PiPriPki and neglecting the off-diagonal matrix element of the density matrix, (3.8) does not contain any approximation. The main approximation will be made at this stage by substituting Wnfni from (2~55), which is provided by the damping theory. In this approximation, Wnfni is proportional to |< nfIVn+H' Ini >|, as can be seen from (2.55). Since H' does not depend on neutron coordinates, and since the initial and final neutron states are not identical, < nf|H' Ini > = 0. Thus, one has only to evaluate < nfiVn(r,R) Ini >, where r is the position of the neutron and R denotes the totality of the positions of the particles in the scatterer. To evaluate this matrix element, one may approximate Vn by the Fermi-pseudopotential,16 i.e., N V = (2J2/m) ~ a~8(r-R,), (359) ~=1 where a, is the scattering length of the ~th nucleus. In general, this scattering length is spin-dependent: aI = A~ + Be(s.S), where s and S are respectively the neutron and nuclear spin, and A and B are nuclear constants. Use of the Fermi-pseudopotential restricts the following analysis to the scattering of slow neutrons (ei. ev) whose wavelength is large compared with the range of nuclear forces. With the foregoing remarks, the differential cross section can be obtained from (308) as follows:

50 I fTfl ale - liTi > Ynfni ar( ) 2TrNki, 2i i Ti TTf i,f (Ei-Ef+e-Snfni) +7Ynfni (3.10) where K k. - kf 1 -- ti - Jf' i f =E (Qh2/2m) k2 i f if i19f Snfni = Re[lim (n-ni)] (5.1a) Ynfni = m[lim -rnil ( 3elb) 1< nIV..+H...ni >1|.2(3. 2) ni = < niiVn+H' ni > L < n +H' ni (1.12) n/ni E-Eni+ic Equation (5310) is the desired differential cross section formula for the scattering of slow neutrons by an arbitrary macroscopic system. The matrix elements appearing in this formula are expressed in the representation generated by. If the latter is chosen as HS, i.e., the Hamiltonian of the scattering medium rather than a part of it, then the shift and the width in (3.10) are caused solely by the neutron interaction. They can then be ignored when the intensity of the incident neutron beam is not too high. Hence, in the representation generated by HS, one may take the limit Ynfni + 0 and obtain the conventional scattering formula.

31 When the eigenfunctions of Hs are not available, (3.10) offers an approximation which is sufficiently accurate to include the broadening and the shift of the lines which may be present in the energy spectrum of the inelastically scattered neutrons. Whether such spectral lines exist or not depends, of course, on the natureof the scatterer. If the energy spectrum of the latter has discrete portions, then the neutron spectrum has a line structure. However, these lines may not be recognizable in the actual spectrum as a result of the broadening and overlapping of the adjacent lines. At any rate, (3510) will be valid. The subsequent chapter contains the application of (3o10) to the neutron scattering from crystals.

CHAPTER IV LATTICE DYNAMICS* We shall give in this chapter a brief review of crystal dynamics. The purpose of this review is to present those aspects of crystal dynamics that are relevant to the subsequent analysis of neutron scattering. To simplify the presentation, we shall consider a monatomic simple lattice with one atom in each unit cell (Bravais Lattice). The generalization to polyatomic crystals is straightforward. The Hamiltonian of a crystal can be written as follows: la laI'a' I(4.1) H = H3) + H(4) +... where H(3) = - x~axR a' x"a",U~, R,3'a", (4.2) ae'a'?,"a" HGI - )o X, a'or'xrt"Xr'"aa', Uer"a"a" H(4) = a! X, ~x~Ix~,a,x~,,~,,x~,,,,,,Ua ~, i, ~,,(,, ~,,,c,,, la,'a', R Qa", 2" o a"' (4.3) *For detailed discussion, see Refs. 18-20. 32

33 In these expressions, M is the atomic mass,;~ is the oCth Cartesian component of the momentum of the ~th atom, x~a is the oth Cartesian component of the displacement from equilibrium of the Rth atom and UaQR a', are the partial derivatives of the crystal's potential energy, viz., uar = X u u (4.4) i 0 R a:o bx xa6x' x:0=O The meaning of the remaining symbols are self explanatory. The H(3) and H(4) are the cubic and quartic anharmonic potentials. The displacements are measured from the equilibrium positions denoted by ~. The latter is given in terms of the lattice vectors ai as = al 1 + a2i2 + a3~3 (4.5) where ~1, ~2, and ~3 are integers known as the lattice indices. We shall always denote them by a single letter I for simplicity. The physical interpretation of U~i n, a, is apparent from (4.1). It is the vector force acting on the atom located at I when the atom at ~' is displaced by an amount x~,. Since a homogeneous displacement of all atoms in any direction causes no net force, the following condition should hold: X, ~a,c,' - 0 O By considering a uniform expansion which does not destroy the transla

34 tional invariance of the crystal, one finds* that C Uoa,' Q', "C" = o~ (4.6) and U~,,l, Q," aT" depend on (T'-~) and (h1"-~) because of the translational invariance of the lattice. We shall now express the crystal Hamiltonian in terms of phonon creation and destruction operators. For this purpose, one first expands X~ and el into normal modes as = E //2M Nwk e e -i (a++ah) (4.7) and 7.____ * -iqc. - = L_ /TwLM/2N e e -- (a+a_k). (4.8) In these expressions, N is the number of atoms within the periodic boundaries. The running index \ denotes the pair (q,j), where q is one of the N allowable wave vectors and j is the polarization index which takes the values j = 1,2,3, in a Bravais lattice. Note that -A = (-q, j), and indicates a plane wave traveling in the opposite direction to that corresponding to ~. The frequency of the normal mode described by % is denoted by wo, Since we allow negative as well as positive values of q, wX is always positive and satisfies wX = wk. The latter condition will be used very frequently in the following *See p. 3355 of Ref. 17 or Refo 6.

35 analysis without being pointed out. Finally, aA and aA in (4.7) and (4.8) are phonon creation and destruction operators* which satisfy the known commutation relations [a,a,] = 0, [a,a,] = (4.9). It is straightforward to verify that the commutation relation [xOa, 7,a'] = i~fb~ tQiol is satisfied owing to the following orthogonality relations: e e ja eqjt, = (4.10) and Z ei -- ) = N= ~ o (4.11) q Substitution of (4.7) and (4.8) into (4.1), (4.2), and (4.3) yields 7~s = Taw%(a ak + ) 12) 3 = CE,(^Xr~~ ff5a ) i (4.12) H(4)=3) a-il i+ai) (4 13) V11\2,X34 1=, H( 4) 4 -. Gj. sk4 (a+ki + (4o14),, +a.,.. 4.14 i=l *The symbol a was also used for the scattering length a~ in (35.9). The different subscripts will distinguish between these two usages.

36 In these expressions, the symbol G is defined by \L4\3 \2MN^, —-i ( q. hi A+q3' h2) E, e - - - Uoahl',hS1' ekl0eoal03a l ee" ~ (4 15) hlh2 (coICxT) The definition of G. is similar to that of G x. The explicit form G\1\2%3%4 will not be needed. The symbol 6(g,qcl+q2+q3) in (4.15) is a kronecker delta expressing the conservation of the wave vectors in phonon-phonon interactions. The vector g is a lattice vector in the reciprocal lattice. Enormous typographical simplification can be achieved if we agree to drop \ in the subscripts whenever no confusion arises. According to this convention, we have the following abbreviations: Wki = Wi G X3= Gl 3, ai3 = G1 3 a = a, a_ = a-i The following symmetry properties of G1,2,3, which play an important role in the subsequent analysis, can be proved directly from (4o15) and the translational invariance of U,, Lo,, "a", * G1,2,3 = G2,1,3 = Gl,3,2 0 (4.16) *See po 504 of Ref, 1l8

37 It also follows readily from (4.15) with e-,) = ethat G_,1,,- 3 = G,2,3. (4.17) The special case of G1,2,-2 will be encountered very often. According to (4.15) it is proportional to E exp[-iq2-(Lh-h2) ]Uou, hl,,h2,, = -iq2Eh Uoc",ha'h' hlh2 h h' The summation on h' vanishes as a consequence of the condition (4.6). Thus, one obtains the important result* Gl,29-2 = 0 o (4.18) We now return to H(3) as defined by (4.13). The eight terms containing the combinations of three phonon creation and destruction operators can be arranged, by virtue of the foregoing symmetry properties of G,2, 3, as follows: 3) = 7 G1 2,3[(a+la+2a+3+ala2a3) + 3(ala+2a3+a+la2a3) ] 1,2,3 (4.19) The derivation of (4.19) is given in the Appendix. The relative magnitudes of the various anharmonic potentials, in particular those of the cubic and quartic anharmonic potentials, will be *See for example p. 37, second footnote of Ref. 17, or footnote 22 of Ref. 6.

38 needed for a consistent scheme of approximations. Van Hove* has pointed out that the expectation values of H(n) decrease very rapidly with increasing n at temperatures which are low compared with melting point. In fact, he has given the order of magnitude of the expected value of the nth order anharmonic potential as w(u/ro) n-2 per unit volume for n _ 2. Here, w is some mean vibrational frequency of the crystal, u is an averaged atomic displacement at a given temperature, and ro is the nearest-neighbor separation in the lattice. As pointed out also by Maradudin and Fein,6 the anharmonic Hamiltonian can be written schematically as Ht = tH^ )+ 2 ) +.. (4.20) where [ = (u/ro) To complete these introductory remarks, it may be mentioned that the eigenstates of t, which form the basic set in the present formulation, can be labelled most conveniently in terms of the phonon occupation numbers. Thus, an eigenstate |i > can be specified as | > |n,n,.0o n,..>, (4.21) where nt is the number of phonons described by X = (q,j). The mean *See p, 11 of Ref. 21.

39 energy in this state is E = Xj iw(n+2) Since all the operators in the cross section formula (3.10) can be expressed in terms of the creation and destruction operators, it is sufficient to note the following relations to evaluate all the relevant matrix elements: a+Jnl, o.o n,. =/ Inl,..o n + 1, o.>, (4.22a) alInl,... n%,...> = I ~n l,.. n\ -1,...>. (4.22b)

CHAPTER V DIFFERENTIAL CROSS SECTION FOR SCATTERING OF SLOW NEUTRONS BY AN ANHARMONIC CRYSTAL A. CROSS SECTION FORMULA We are now in a position to write down the neutron cross section formula for a Bravais crystal. Substituting RI = _+xI in (3.10), assuming that the spin and spatial states of the scatterer are separable, and denoting the spin states by Isi >, |sf > and spatial states by |i >, If > [note that |i > and |f > in (3.10) include both spatial and spin states], one obtains = -NkiL ei2( ) E 2 AIAl Pi IlI II' SiSf i,f Ti Tf x. —--— nf (5.1) Ynfni ( X ~wmwO-Snfni) +Ynfni where AI - < TfSfia~siTi >, (5.2) I = < f|e - i > o (5.3) In this expression, mX denotes the difference between the final and initial occupation numbers, viz,, 40

41 f i m; = nk - nk. Henceforth we shall specify the final states in terms of mX rather than f i in terms of no, Then, the superscript i on n; becomes redundant and hence will be omitted. Note that the summation on the final spin states appearing in (51o) cannot be performed because the width Ynfni and the shift Snfni are also spin-dependent. However, when the width and the shift due to the neutron interaction are neglected, then the last factor in (5.1) becomes spin-independent, and the summation on spin states reduces to 1 < ^il^^ \Q~r~ > ~ (5o4) - X < SiTi|aja+, I siTi si,Ti The foregoing discussion indicates that the spin dependence of the neutron-nucleus interaction does not affect the shape of the lines. (This point will be more apparent in the subsequent section.) Therefore, we shall ignore spin effects and replace the average in (5.4) by a2. The cross section formula (5o1) then reads kfa2 ei T) Z 7nfni r( =kfa2 i( ) i I, i*(ni'' -' ~1' i,f 1' 2 if -Z wm\'-Snfni nfni (5.5) Our next task will be to compute the IQ, ynfni' and Snfni appearing in the cross section formula.

42 B. COMPUTATION OF MATRIX ELEMENTS OF exp( i.-xX) This section is devoted to the computation of IQ defined by (5.3) Using the expansion (4.7), one can write x~ also in the following form: xi = (ce aI + CP,), (5.6) where C = (i/2MNwN) 1/2 (e. K) exp(iq. ). (5 7) The matrix element to be computed is I = < f7l exp[i(Ceaa + caa%) ]li > Using the operator identity A+B eAeB e-(1/2) [A,B] e = e=e e-, (5.8) where A and B are operators satisfying [A,[A,B]] = [B,[AB]] = 0. and expressing the initial and final states in terms of occupation numbers, one get s I1 = t e -X / < mm+njlei a e n n >, (5o9) where X\ denotes x -- 2Ic 2 = (2/2MNwN) IeL%_Kj2 o (5.10) Note that X\ is independent of 2. To compute the matrix element in (5.9), one uses the following relations:

43 n e i In> = j (i>t In-s >, (3.11) ^-j s. 1 1) * s=O 00 eia*a+ in~s > 7 (n-s-t)! (ic*)t e I|n-s > 5 Y (n-s) (t 1 In-s+t>. (5.12) (n-s) t! t=O The result is 00 00 I T e-X%/2 Y (i ( iatX) nl(n%-s+t) eti n~:Fm%., n%- s+t't=O s=O t. (n%-s)! nmnor n. I2 =X )~( eXXI~2 ( r*nX] ( *-XS e- iClk) - nn (nk+mA)... ";; s=O s'(n%-s)!(mr+s) (5.13) as I e-X/2 (iCe)m AX(X) n ( 5m14 which is the desired formula. he last expression can be written in a compact3) follows from (recalling the definition of associate Laguerre polynomials, viz., (5 n ( = X) ( n+m)_ s=O (n-s)l(n+s)! S! as I2 LF X eX- %/2 ** n% (X2\), (5.lh) which is the desired formula. 11 = 7- e - - $ (nmX), (5.15) A.

44 where n(n,m,x) - xm n e [L(x) ]. (5.16) (n+m) We shall now calculate the thermal average of "(n%,mX,X;) with respect to the initial occupation numbers, i.e., < (nm hX^)> - (l-y) E Y (n (5.17) n2=O where Y - -(~ew/2kT) (5.18) Dropping the subscripts and arguments, and combining (5.16) and (5.17), one finds 00 <4> = (l-y2)xme-X y n n [I:(x) ]2 n (n+m)! n=0 To compute the summation in the last expression, one may use the following expansion, given by Magnus (po 85, Ref. 22): 00 x2S m+2s n!- M F X (2x -n [!T(x) Ln-s ( 2x) (5.19) (n+m) =0- s; (m+s)e s-= The result is 00 00 2S < (l-y2)e-x e sx 7+s) lIS (2x)y2n s=O n=0 (5.20) The summation on n can be performed by using the generating function for the associated Laguerre polynomials (p. 84, Ref, 22), viz.,

45 -xt/(l-t) (lt)e E L n(x)tn, n=0 as follows: 00 MT 2 2s -xy2/(1-y2) +S(2x)y2n = y 2m++. (5.21) (l_y2)m+2s+l n=O Inserting the last formula into (5.20), one finally obtains* < > = -m e-x(l+y2)/(l-y2) I 2y (5.22) where Im(x) is the modified Bessel function of the first kind. This formula will be used in Section D when we are discussing the cross section formula. C. WIDTH AND SHIFT FORMULAS In this section, we compute Ynfni and Snfni. They are defined by (3.11) and (3.12). It is observed from the latter that they involve the matrix element 1< nIVn+H' ni >12, which can be written as follows: n \< 4nl(K) |i >| + Q 2Re[ < ln'(o) li > < HIH' |i >] 5(k,k-) + I< tH' Iie >I2 (k,kf ), (5.23) where ( (K) is the Fourier transform of Vn, i.e., /(2-n(K) N 1 / iK.r 2 n4i2 iK. R - e — n(r) a -- Q (5.24) *An alternative derivation of (5.22) was also given by 0. Ruehr, the Radiation Laboratory, The University of Michigan (private communication).

46 and where K = k.-k. It follows from (5.24) that (n (0) is independent -"1 - of the atomic coordinates R~; thus the matrix element < j L,9(0) i > is proportional to 6io. The second term in (5o23) is then nonvanishing only for n = njo Since the summation in (3.12) excludes the case n = ni, there is no contribution from this term to the width and shift. Accordingly, one can split rni in (3.12) into two parts as = n p Fni = n + Ii, (5.25) where rn Vni n ESn-ni I, (5.26) ni - nini -E+i 1 Hi En-Eni+i r H.. (5.27) The first term, rni, represents the broadening and the shift due to transitions caused by neutrons. It corresponds to the natural broadening of the optical lines. The second term, I?, is due to transitions caused by phonon-phonon interactionso In many applications, the neutron width is expected to be small compared with rio Therefore, only the latter will be considered in the remainder of this work, As indicated by the subscript i, FP depends only on the spatial part of the crystal states provided H' does not include spin dependenceo Now consider (5o27) First note that H(') has no diagonal matrix elements since it involves terms containing an odd number of creation

47 and destruction operators [cf. (4.19)]. The first term in (5.27) starts with H(4), which is second order in., as indicated by (4.20). The remaining off-diagonal matrix elements in (5.27) start with HN3). But, they are also of the second order in 4 because of the square. Hence, r`. can be given, in the lowest order, as rP = H - )< <\l-EH~ iH >. (5.28) This is the quantity, provided by the damping theory, which gives the width and the shift according to (3.11). The shift is related to the real part of PI. Hence, one finds that the quartic anharmonic potential contributes only to the shift in the lowest approximation, whereas the cubic anharmonic potential contributes to both the shift and the width. The remaining task is to compute the matrix elements in (5.28). First consider H(. It is clear from (4.14) that the nonvanishing terms in the expansion of H are those that contain creation and destruction operators in pairs. The possible permutations are listed below: aa2+( aa2 + a2al) = 2N1N2 aala2+a2 + a+a2a+al = N1(2N2+1) a+ala2a2 + a+a2ala+ = 2N1(N2+l) ala+a2a+ + aaaala+ = (Nl+l)(2N2+l) ala2a+a + ala2a+a+ = 2(N+l) (N2+l) alata.a2 + a2aa2+aI = N1(N2+l) + (l+N1)N2

where Ni is the number operator, i.e., Ni = atai. We recall that ai and Ni denote aki and Nxi. Since (4.14) involves a summation on %i and %2, the terms containing a+ and a+ in the reversed sequence should not be considered separately. Using the invariance of G1,2,3,4 with respect to the interchange of the subscripts, one obtains Hii = G1,-1,2,-2(nln2+nl+n2+3) ~ (5.29) 1,2 (4) (4) As indicated by (3511), the shift depends on Hff -Hii The matrix element in the final state can be obtained easily from (5.29) by replacing the occupation numbers ni by ni+mi. Hence, (4) (4) 1 Hff - Hii = L G1,-1,2,-2 (2nlm2+mlm2+mi). (5.30) 1,2 We now consider the second term in (5.28). To compute the matrix element involved, we first calculate H(3) i >. The use of (4.19) yields: H(3)| i > G= [(l+n_l)(l+n.2)(l+n.3) ]l/2 nl,+l,n-2+,n-3+l> 30 + [nln2n3]/2 |n,-l,n2-l,n3-l > + 3[nl(l+n-2)n3] / In,-l,n,2+l,n3-l > + 3[(l+nl)n2(l+n.3) ]1/2 In_+l,n2-l,n-.3+l >. (5.31) Since the operator ( -Ei+ie) is diagonal in (li >, it does not alter the occupation numbers in the intermediate state H() i >. Therefore, the changes in the occupation numbers caused by the first H(3) have to

49 be removed by the second H(3). It is convenient to write the latter by reversing the sign of the subscripts and using G-1,-2,-3 = G1,2,3 as follows: H(3) 1 ~+ +G~ 2 [+ +\++ (H = 1 G,2,3[alta-a+a-ia,2a-3+++ata-2a)] 1,2,3 (5.32) Each term in (5.31) has to be matched by the appropriate term in (5.32). For example, the first term in (5o31), which corresponds to the creation of three phonons characterized by -Al, -\2, and -s3, will be paired by the second term in (5~32), However, since two sets of three numbers can be paired in 3! ways, this term has a factor 35. Similarly, one can see that the combinatorial factors for the second, third, and fourth terms in (5.31) are 3', 21, and 21 respectively. With the foregoing remarks, one finds: <i =3) (7S-Ei+iE) (3) > = i Gl2,3 2 t(l+nl)(l+n- 2)(1+n-3) 1,25 3 wl+w2+W3+ie nln2n3 nl(l+n-2)n3 (l+nl) n2( +n3) +3 +3 (5 33) W1+W2+W3-ie -W1+W2-W3+iE W1-W2+W3+ie In obtaining (5.33), we have not considered the special cases in which a phonon is created and absorbed at the same vertex (instantaneous phonons), or in which two or three intermediate phonons are identical. Instantaneous phonons correspond to those terms in (4.19) in which a pair of creation and destruction operators belong to the same mode, viz., \2 = -\3 or \.2 = -hlo Since these terms are multiplied by G1,2,-2,

50 they vanish by by virtue of (4.18)o Hence, instantaneous phonons do not occur at three phonon vertices in a Bravais lattice.6 Identical intermediate phonons arise when two or three creation (or destruction) operators in (4.19) are of the same kind. As an example, consider the second term in (4.19) and assume that \2 = \3 This term would give rise to nl(n2-l)n2 wl+2w2-ic which is different from the corresponding term in (5.33) obtained by setting \2 = X30 The appropriate corrections for these special cases can be included in (5033)o But their contribution can be ignored, as will be discussed presently. The matrix element in the final state can be obtained from (5o33) again by replacing nX by nk+mo The difference between the two matrix elements is then found as: < flH(3) (s-Ef+ie)-' 3) f > - < ilH(3)( Ei+i)1H(3) i > = E G,G1 3 1 F3ml( l+n2) ( +n3) +3mlm2( l+n3) +m1m2m3 3,'2 Wl+W2+W3+iC 1,2,, 3 _ 3mln2n3+3mlm2n3+m m2m3 Wl+W2+W3-ie 2ml(l+n-2) n3+m-2nln3+2mlm-2n3+mlm2(l+n2)+m)+m-2m3 + 3 -W1+W2- W3+i + 3 2mln-( (l+n3) +m-2( l+nl) (l+n3) +2mlm-2( l+n3) +mlm3n-~2+mlm-2mi3 Wl-W2+w3+ie J (5.34)

51 The shift and the width can now be obtained by inserting (5o30) and (5.34) into (3.11): OP -1 1 = 2 E G,-1,2,-2 m(2nl+l+m2) - G1 2 232 1,2 1,2,3 x PP fmi( l+n2+n) +mlm2 + 2ml(n3-n-2) -m-2( +nl+n3) -2mlm-2+mlm3 wl+W2+W3 W2-W1-W3 (5. 35) 1,2,3 + mlm-2( 2n3+l) +mm3( 2n-2+1) +2mlm-2m3 16(w-w2+w3) o (5 36) In the case of creation or annihilation of phonons of one kind, the above formulas simplify considerably. Then, the final state differs from the initial state in the occupation number of the excited (or deexcited) mode only. Let the latter mode be denoted by \o. Substituting mu = mO, for \ = o, = O, otherwise, one obtains the shift and the width in thise case as follows: =?- Go,-o0,,-l(2nl+l) - m IGO21l2 1 1,2 x p{ l+nl+n2 + n-1-n2 _ l+n n-2 Lo+wl+W2 Wo-W1+w2 Wo-W1-w2

52 i = i [I0 G0.,2 12 [(l+nl) (+n_2)+n-1n-2 ](Wo-W1-W2) 1,2 + 2[(l+n._l)n2+nl(l+n2) ]6(wo-wi+w2) (5.38) Note that we have disregarded the terms containing the product of two or three mx in obtaining (5.37) and (5038). These terms are nonvanishing only in the special cases discussed previously. Consider the terms, both in (5035) and (5o36), which contain mlm_2. They are nonvanishing only when \i = -2 = Zoo But in that case, GXO-o0.2 vanishes by (4o18). The terms containing mlm2 (or mim3), which are proportional to Gko\o\ly do not necessarily vanish. Their contribution to the shift and the width in (5.37) and (5.38) however, may be neglected as compared with the remaining terms because they involve one more constraint than the other terms, Equations (5o37) and (5o38) are the shift and width formulas of the present theory. The following physical interpretation for the width seems to be in order. Suppose that a phonon (qo,jo) has been created by neutron interaction. This phonon is of course indistinguishable from the already existing phonons in that particular normal mode. In view of this indistinguishability, the concept of lifetime of the created phonon requires careful consideration. One can define the lifetime of an additional phonon in a given normal mode as the reciprocal of the absolute value of difference between the widths associated with two crystal states which differ from each other by 1 in the occupation

55 number of the mode described by qo, o. The lifetime defined as above depends on the initial crystal state. The mean lifetime can be obtained by taking the thermal average. According to the foregoing definition, the lifetime of an additional phonon is determined by the transitions involving the normal. mode q0,jo. There are four types of such transitions if the interaction between phonons is approximated by the cubic anharmonic potential. These are shown schematically as (2oJo) = (q1,ai) + (qSJ2), (5.39a) (,jOo) + (q1, A ) - (2, j2) * (53.9b) The first and the third of these transitions correspond to the annihilation, and the remaining two to the creation, of a phonon in the normal mode under consideration. An inspection of (5.38) reveals that yP confi tains just these four phonon processes with appropriate amplitudes. The concept of lifetime as introduced above applies equally well to the annihilation of a phonon. In the zero temperature limit, the foregoing definition of phonon lifetime becomes more concrete. In this case, the width of the initial state, i.e., the ground state of the crystal, is zero. Since there is only one phonon in the final state, the lifetime of the latter is actually equal to the mean decay time of the created phonon.

54 D. DISCUSSION OF THE CROSS SECTION FORMULA We are now in a position to discuss some aspects of the cross section formula (5.1). Since we have decided to specify the final state in terms of mk, i.e., nA = m%+nk, the summation on final states in (5.1) can be replaced by summation on (ml,m2,...mX,..m3 N). The limits of these summations are -nk and +oo, But, as indicated by (5.15) and (5.16), the matrix elements of exp(i-.xg) vanish automatically when mk < -nh. Therefore, one can replace the lower limits by -o. Since the limits then become independent of the initial occupation numbers, the order of the summations on initial and final states can be interchanged. With these remarks, the cross section formula takes on the following form: 5N +4oo 0 a(w,) - ki (-y 7)ynF(mj) =1 m=-oo nk=O ~x _r,, (mn),, (5.40) w- j wkm -S(m,n) +r2(m,n) where we have introduced F(mk~ ~) -= i ( _ -m) o(~-~') eF(mi(-m )( ), (5.41) mn' w -- (i/2m) (k-kf) (5o42) r(m,n) - F(m=,m2,. o.m3N;n o,n2 ooon3N) - (1/) 7fi, (5.43) S(m,n) = S(ml,m2.. oom;nnl,n2..m n) - (1/)SP f. (5.44)

55 The expression for Sgi and fi are to be substituted from (5.35) and (5.36). All the quantities entering the cross section formula have now been defined explicitly. The remaining task is to consider the various limiting cases, such as zero and high temperature limits, and to extract some information about the shape of the observed peaks in the energy spectrum of scattered neutrons. First, we want to show that (5.40) reduces to the cross section formula for scattering of neutrons by harmonic crystals, when r and S tend to zero. Indeed, when r = S = 0 the last factor in (5.40) approaches a delta function and becomes independent of the initial occupation numbers. Then the thermal average of ~(m%,n%,Xk) can be performed independently, which follows from (5.22) as 00 j (ly_)y2n mknX) = Y e 2W I(PW), (5.45) n~=O where we have introduced for brevity 2W - X(l+y)/(l-y2), (5.46) P 2Xy/(l-y2) (5.47) The factor exp(-2W) is known as the Debye-Waller factor. With the foregoing remarks, (5.40) reduces to +~0 c(w,) = a2(kf/ki) e- 2W6 w - w j wm J F(m,^ )ym Im(Px), m=-oo X (5.48)

56 where m denotes the set (ml,m2, oo o3N). This formula is identical to the differential cross section formula obtained by Zemach and Glauber.*23 The novelty about (5~48) is that it is obtained, without resorting to Bloch's theorem to perform the thermal averages, by the mathematical identities proved in Section B. The present derivation seems to be more straightforward. It might be informative to mention that we could not use Bloch's theorem in the original cross section formula (5.1) because of the last factor in the formula which appears as a delta function in the conventional formulation of the cross section. We could not perform the summation on final states, and this is needed for the application of Bloch's theorem. Second, we investigate the cross section in the zero temperature limit, where all the occupation numbers are zero. As mentioned previously, the crystal is in the ground state in this limit. Equation (5.40) reduces to the following form as T + 0: 3N +0o a a2 kf -2W(o)(k a( w,F ) (e Lx/m j! m=0 r(m,o) x r(mo) (5.49) X 2 w -E wm-mS(m,O) +r2(m,O) where we have used C (m OX ) o= e X( m/ ) *See also p. 52 of Ref. 24.

57 In the last formula, r(m,O) and S(m,O) are to be obtained from (5.35) and (5.36) by setting n1 = n2 = n3 = 0 everywhere. It is noted that (5.49) contains only the positive values of mk as a result of m%! in the denominator. Physically, this implies that the crystal cannot give energy to neutrons when it is in the ground state. The cross section for one-phonon excitation, which is probably the most interesting case from the experimental point of view, can be obtained from (5.49) by setting mk = 1, for \ = ko = 0, otherwise. The result is,)A -a2f kf K-7_ r C ) = 2JtM ki (g, o) WO (W S) +r wo (W-W0-S)2+r (.50) where S and r are obtained from (5.57) and (5.38) by substituting mo = 1 and ni = 0: S = j Go- I + X IG OG122 L 1 1 1,2 (5.51) r = Ij |Go,,212 (Wo-Wl-W2). (5.52) 1,2 Apart from differences in notations, this formula is identical to that obtained by Maradudin and Fein.

58 We now consider the cross section at a finite temperature. Equa — tion (5.40) expresses the cross section as a weighted superposition of a sequence of Lorentzian distributions with different widths and shifts. For each initial and final state there corresponds a Lorentzian distribution. For a specified final state, the average of the appropriate Lorentzian distributions over initial states gives a peaked curve which corresponds to an observed line in the neutron spectrum. Although the individual distributions have a Lorentzian form, their thermal average will not in general be a Lorentzian distribution. However, as an approximation, one may replace the resultant curve by a Lorentzian distribution with an average width and shift. Denoting the average width and shift by T(m) and S(m), one obtains the following cross section formula: +00 A a2 kf ^ ^ r^-^ - - m= ki W w X -( +r (5-53) The average shift and width are to be obtained from (5o35) and (5.36) by replacing the occupation numbers nX by their mean values, viz., n = y /(l-y\ ) Note that the approximation made in getting (5.53) from (5.40) is equivalent to replacing the average of a function by the function of the averaged variables. The cross section for one normal mode interaction is obtained again by choosing the set m = mim2o meN in (5o53) as m = (O,O,..m,O,...) The result is

59 o(w,~) a= - a2 F(mo,y)Y 0 I-2W(o) -( Im l) r ki [w-moWo-sw(mo) ]2+r2( Imo l) x C Io(P) (5.54) where S(mo) 2 IE GO,-0,1-1 (2n1+1) - Go 1 22 L 1 1,2 x pp 71+1+nln2+__n2 l+nl+n2 x PP i + n l ]. (5-55) Wo+Wl+W2 WO-W1+W2 Wo-lW-JW r( Imol) = mO I |o | GO 1,2 fI [( l+ni) (1+2) +nn2 ]5(Wo-W1-w2) 1,2 - + 2[(l+n1i) 2+En(l+E2) ](wo-wl+w2). (5.56) Equation (5.54) can be further simplified in many applications by using the asymptotic value of Im(x) for small arguments, i.e., Im(x) ~ (x/2) Iml 1 Iml! as follows: ca(w,) =a2 kf -2W F(mo, I) 1 _O I m a(,) - We F-MOK) ki r(| ).051 1 ~ i t ki - Imol|: L2MNwo j e-PWo( Imo I-mo)/2 Im x(Ie-BhwO) Imo l C,-,om,-~cm,, ]2 ~r( 7mo I )

60 In the zero temperature limit, this cross section vanishes when mo < 0 as a result of the factor exp[-fio(w [mo|-mo)/2]. This result is in agreement with (5.49). As a matter of fact, (5o57) reduces exactly to (5049) when the temperature approaches zero, although it involves one more approximation over (5.49), i.e., replacing the width and the shift by their mean values. The reason for this is that the latter approximation becomes exact in the zero temperature limit. The dependence of the cross section on moo, i.e., the number of phonons exchanged between the crystal and the neutron, is exhibited by (5.57). Because of the exponential factor in the numerator, the cross section for the energy transfer from the neutron to the crystal is greater than the cross section for an energy transfer to the neutron by a factor exp[-Piwo moj]. At high temperatures, both processes become equally probable. The shape of the peaks in the neutron spectrum depends on Imo0| Hence, the shape is the same for both the creation and annihilation of phonons. Furthermore, since the width is proportional to Imol, the peaks are broadened and reduced in height when the energy exchange gets larger. The shift of the lines is proportional to mo, and thus depends on the direction of the energy transfer, There will be two identical peaks in the neutron spectrum on both sides of the incident neutron energy, corresponding to the creation and annihilation of equal number of phonons. The cross section formula agrees examuj^y wj~l ~L1au uuD ~41~^ oy Maradudin and Fein6 except for the value of the width (5.56). The dis

61 agreement lies in the sign of the terms hTli2 and 1n(l+n2) in (5.56). These two terms arise from the imaginary part of the last term in (5.34) when -+O. One observes that the real parts of the last two terms which enter the shift formula are subtracted, whereas the imaginary parts are added. Since our shift formula agrees exactly with that obtained by Maradudin and Fein, it is very unlikely that the foregoing discrepancy in the width formula is just an arithmetical error in the present calculations. Moreover, the terms appearing in the width formula (5.56) have simple physical interpretation. As mentioned earlier, the first and third terms correspond to transitions in which a phonon is annihilated, whereas the remaining terms whose sign is the point of disagreement correspond to those transitions where a phonon is created. In either case the number of phonons in the mode under consideration changes. It is not clear to us why the transitions corresponding to an increase in the phonon number should tend to decrease the width, and thus to prolong the lifetime. As a result of the foregoing discrepancy, the temperature dependence of the width at high temperatures is quadratic in our case, whereas it is linear according to the result obtained by Maradudin and Fein. The expression for the shift and width in the high temperature limit is obtained by replacing the mean occupation numbers in (5.55) and (5.56) by kT/iwk. Hence, one obtains

62 SM= o kT 2 Go-o0,1-i - o S(mo) _1 mW1W2 x PP WW2 + W2-W1 Wl+W2. (5.58) Wo+wl+w2 Wo-Wl+W2 Wo-Wl-W2 r( |mo|) = Irmol ( ) j i 2 —L -- [(Wo-Wl-w2) +5(Wo-Wl+w2)] W1W2 (5.59) As a final remark we note that the cross section formula at a finite temperature, ie., (5.55), was obtained from (5.40) by replacing the width and the shift by their mean values. As discussed in Ref. 10, this approximation is justified when the width of the individual Lorentzian distributions corresponding to different initial states is large compared with their shift. If this is not the case, the observed width will be dominated by the statistical spread in the locations of the narrow Lorentzian distributions, A measure of this statistical spread is the standard deviation of S(m,n) defined by (5.44), viZo, rP(m) S2 - s (5.6o) As discussed in Chapter IV-D of Ref. 10, rs(m) may be added to F2(m), as a first approximation, in order to include the statistical broadening in (5 54). The relative magnitudes of T(m) and T2(m) will depend on temperature. We shall not attempt here to write down the expression for the statistical broadening because it is rather long and

65 not very informative. However, it is straightforward to calculate r2(m) from (5.44). We only note that it is a linear, homogeneous function of the variance a2 of the occupation numbers, i.e., at - no - n = (l ew )2 (5.61) X? (\1-e PhW%) At zero temperature all o are zero, indicating that the statistical broadening, as expected, is not present. At high temperatures, a2 and thus r2(m) are proportional to T2, whereas according to (5.59), T2(m) is proportional to T4. Therefore, it may be expected that the statistical broadening will not be significant at any temperature in a crystal. It is interesting to recall that the width of the optical lines from a plasma is determined by the statistical broadening in the quasistatic limit, i.e., when the motion of the perturbers is sufficiently slow. In the impact limit, the width is determined by the finite lifetime of the emitter states. It thus appears that the width of the observed peaks in the neutron specrum is caused primarily by he finite lifetime of the crystal states rather than by the statistical spread of the lines.

CHAPTER VI SUMMARY AND CONCLUSIONS In this work, a general theory of line shape which is applicable to the study of both neutron and photon spectra has been developed and used to investigate the scattering of slow neutrons by an anharmonic crystal. The starting point of the present theory is (2.55), viz., Wmn = mn(t) 2 -IVm 12 1Y7nYmI mn I Vmn 12 [(Em+Sm) -(En+Sn) ]2+(i/2) (Yn-7m)2 This formula reduces to the conventional expression of the transition probability per unit time from the initial state In > into the final state Im >, i.e0, Wmn = IMI2 V(Em-En) Wmn Vmn ( E)' when the states are sharp. Equation (2.55) is an extension of Wmn = IVmn 12 (Em-En-Sn) 2+(/2) 2n which has been obtained by Heitler9 to investigate the natural broadening of the optical lines. The latter follows from (2o55) by setting Sm = Ym = 0, and thus applies only to the cases where the width and shift of the final state is negligibleo In the study of the optical lines, the interaction of the emitting atom with its surrounding in the final state can often be neglected since the atom is more tightly 64

65 bound in the final state, which has a lower/ energy. This is particularly true if the final state is a ground state of the atom, In such applications, the neglect of Sm and Yn is justified. However, in the study of line shape in the neutron spectrum, or in the study of optical line shape involving transitions between two excited levels, the decay-of the final state cannot generally be ignored. In such applications, the use of (2.55) is imperative. When applied to the scattering of slow neutrons by an anharmonic crystal, the present theory yields, in the lowest approximation, the following expressions for the width r and the shift S of the observed peaks [cf. (5.28) and (3.11)]: r = < Im[lim ((-f- P)] >T S = < Re[lim ( r I)] >T where i - <n (H(4 - (3) (1 3 )) l> \ if: -Ei +i~ where the final state differs from the initial state by 1 in the number of the occupation number of the normal mode under consideration, and where the symbol <.**>T denotes the thermal average. At any temperature, the shift S calculated by the foregoing formulas agrees exactly, apart from differences in notations, with the shift formula given by Maradudin and Fein.6 However, the width formulas agree exactly only in the zero temperature limit. At finite temperature, the width

66 formula ($.56) differs from (5.5b) of Ref..6 in the sign of the second and fourth terms. The present theory predicts a larger width than that predicted by Maradudin and Fein. Furthermore, the temperature dependence of the width at high temperatures is quadratic according to (5.59) of the present work, whereas it is linear in Ref. 6. The present width formula may yield a better agreement with the experimental data than the agreement reported by Maradudin and Fein, because their calculated values seem to be smaller than the experimental values obtained by Brockhouse et al.1 When applied to the study of the optical lines,10 the present theory yields a line shape formula which contains the two limiting approximations, i.e., quasi-static and impact approximations, as special cases. It treats the electrons and ions on the same basis, thereby including the motion of ions in the calculations. It may be concluded in general that the present theory of line shape is more systematic, simpler, and more interpretable than the existing line shape theories both for neutrons and photons. The theory can now be applied with reasonable confidence to the interpretation of a variety of experiments involving photons and neutrons, e.g., lasers and scattering of neutrons by liquids.

APPENDIX SIMPLIFICATION OF CUBIC ANHARMONIC POTENTIAL Write (4.13) explicitly~ _3 I + + + + + + + H(3) G1 2,3 (a-la-2a-3+ala2a3+a-la-2a3+a-la2a-3 1,2,3 + ala-~2a++ala2a+3+alat2a3+at a2a3) (A. 1) Consider the third term, and observe the following identities: G1 2,3 ala 2a3 = G1,2 3 af at3a2 1,2,3 1,2,3 = i G1,2,3 a+la2a+3 + G1,2,3 af[at+3,a2] 1,2,3 1,25,3 (A (A. 2) Note that the first term of the last line in (Ao2) is identical to the fourth term in (A.l). Repeat the same procedure for the fifth, sixth, and eighth terms in (A.l), and get H(3) G, 2 23[(at+at2at3+ala2a3) + 3(ata2a+3+alat2a3) ] 67

68 The second summation in (A.3) vanishes owing to the commutation relations (4.9), viz., [a-3,a2 = -62,-3,[a2,a+] = 6-1,2,[a3,a+2] = 53,-2 [at2,al] = -1,-2 Using these in (A.3), one finds Z [-G1,-3,2 a-i +G1,-1,2 a-2 +Gi1,22 a1 -G-_ll,2 a2] 1,2 which vanishes as a result of the invariance of G1 2 3 under the interchange of the subscripts. The same result follows also from (4.18), which indicates that the individual terms in the last expression are zero.

REFERENCES 1. B. N. Brockhouse, T. Arase, G. Caglioti, M. Sakamoto, R. No Sinclair, and A.D.B, Woods, Inelastic Scattering of Neutrons in Solids and Liquids, po 531, International Atomic Energy Agency, Vienna (1961)o 2. L. van Hove, No Mo Hugenholtz, and Lo P. Howland, Quantum Theory of Many Particle Systems, Benjamin, New York (1961). 3. J. Krieger, Phys. Rev, 121, 1388 (1961) 4. J. Kokkedee, Physica 28, 374 (1962) 5. G. Baym, Phys. Rev. 121, 741 (1961). 6. A. A, Maradudin and A, Eo Fein, Phys. Rev. 128, 2589 (1962). 7. M. Baranger, "Spectral Line Broadening in Plasmas," p. 493, Vol. 13 in the series Pure and Applied Physics, ed. by D, Ro Bates, Academic Press, New York (1962). 8. E. Arnous and W. Heitler, "Theory of Line-Breadth Phenomena," Proc. Royal Soc., Series A; Matho and Physo Sciences, 220-1142, p. 290510 (1955). 9. W. Heitler, Quantum Theory of Radiation, 3rd ed,, Oxford (1954). 10. A. Ziya Akcasu, "A Study of Line Shape with Heitler's Damping Theory," Technical Report 04836-1-T, The University of Michigan, April, 19630 11. Vo Weisskopf and E. Wigner, Zeit. f. Phys. 63, 54 and 65, 18 (1930). 12. A. Messiah, Quantum Mechanics, Vols. 1 and 2, Interscience, New York (1963). 13. P% R. Halmos, Hilbert Space, Chelsea, New York (1957). 14o R. Co O'Rourke, "Damping Theory," TNTRL-5315 (1959) o 15. Ro Ko Osborn and E. H. Klevans, "Photon Transport Theory," Annals of Phys 15, 105 (1961), 69

70 REFERENCES (Concluded) 16. E. Fermi, Ricerca Sci. 7, 13 (1936); English translation, USAEC Rept. NP-2385, 17. R. E. Peirls, Quantum Theory of Solids, Oxford (1955). 18. M. Born and K. Huang, Dynamical Theory of Crystal Lattices, Oxford (1954). 19. Jo Mo Ziman, Electrons and Phonons, Oxford (1960). 20. A. A. Maradudin, E. W. Montroll, and G. H. Weiss, Theory of Lattice Dynamics in the Harmonic Approximation, Academic Press, New York (1963) 21. L. van Hove, Phys. Rev. 95, 249 (1954). 22. W. Magnus and Fo Oberhettinger, Formulas and Theorems for the Functions of Mathematical Physics, Chelsea, New York (1949). 23. A. C. Zemach and R. J, Glauber, Phys. Rev. 101, 118 (1956). 24. S. Yip, R. K, Osborn, and C. Kikuchi, Neutron Acoustodynamics, Univ. of Mich., Industry Program of the College of Engineering, Ann Arbor, unpublished (1961).