\documentclass[12pt]{article}
\usepackage{wrapfig}
\usepackage{amsmath}
\usepackage{graphicx}
\message{
*************************************************************************
 	Copyright, 1995-2017, all rights reserved, Mark Jarrell (Dept.of 
 Physics and Atronomy, Louisiana State University, LA 70803).  This 
 material may not be reproduced for profit, modified or published in any 
 form (this includes electronic redistribution) without the prior written 
 permission of the author listed above.	
*************************************************************************
}

\title{Chapter 4: Crystal Lattice Dynamics}
\author{Debye}
\input{../defs}
\def\baselinestretch{1.4}
\bdo
\maketitle
\large
\tableofcontents

%\Large
\pagebreak
A crystal lattice is special due to its long range order.  As you
explored in the homework, this yields a sharp diffraction
pattern, especially in 3-d.  However, lattice vibrations are
important.  Among other things, they contribute to 
\begin{itemize}
\item the thermal conductivity of insulators is due to dispersive 
lattice vibrations, and it can be quite large (in fact, diamond has 
a thermal conductivity which is about 6 times that of metallic copper).

\item in scattering they reduce of the spot intensities, and 
also allow for inelastic scattering where the energy of the scatterer
(i.e.\ a neutron) changes due to the absorption or creation of a phonon
in the target.  

\item electron-phonon interactions renormalize the properties of
electrons (make them heavier).

\item superconductivity (conventional) comes from multiple electron-phonon
scattering between time-reversed electrons.

\end{itemize}

\section{An Adiabatic Theory of Lattice Vibrations}

At first glance, a theory of lattice vibrations would appear impossibly 
daunting.  We have $N\approx 10^{23}$  ions interacting strongly (with 
energies of about ($e^2/\A$)) with $N$ electrons.  However, there is a 
natural expansion parameter for this problem, which is the ratio of the
electronic to the ionic mass:
\beq
\frac{m}{M} \ll 1
\eeq
which allows us to derive an accurate theory.

     Due to Newton's third law, the forces on the ions and electrons are comparable $F\sim e^2/a^2$, where $a$ is the lattice constant.  If we
imagine that, at least for small excursions, the forces binding the electrons 
and the ions to the lattice may be modeled as harmonic oscillators,
then
\beq
F \sim e^2/a^2 \sim m\omega_{electron}^2a \sim M\omega_{ion}^2a
\eeq
This means that 
\beq
\frac{\omega_{ion}}{\omega_{electron}}\sim\lep\frac{m}{M}\rip^{1/2}
\sim 10^{-3} \;\;\mbox{to}\;\;10^{-2}
\eeq
Which means that the ion is essentially stationary during the
period of the electronic motion.  
\begin{figure}[htb]
\centerline{\includegraphics[width=0.49\textwidth,clip=true]{setup.pdf}}
\caption[]{\em{Nomenclature for the lattice vibration problem.  $\s_{n,\alpha}$
is the displacement of the atom $\alpha$ within the n-th unit cell from
its equilibrium position, given by $\r_{n,\alpha}=\r_{n}+\r_{\alpha}$,
where as usual, $\r_n=n_1\aone+n_2\atwo+n_3\athree$. }}
\label{fig:setup}
\end{figure}
For this reason we may make an adiabatic approximation:
\begin{itemize}
\item we treat the ions as stationary at locations $\R_1,\cdots\R_N$ 
and determine the electronic ground state energy, $E(\R_1,\cdots\R_N)$.
This may be done using standard ab-initio band structure techniques 
(DFT, GGA, etc.).
\item we then use this as a potential for the ions;  i.e.. we     
      recalculate $E$ as a function of the ionic locations, always
      assuming that the electrons remain in their ground state.
\end{itemize}
Thus the potential energy for the ions
\beq
\phi(\R_1,\cdots\R_N)= E(\R_1,\cdots\R_N)+\mbox{the ion-ion interaction}
\eeq
We will define the zero potential such that when all $\R_n$ are at 
their equilibrium positions, $\phi=0$.  Then
\beq
H=\sum_n\frac{P_n^2}{2M} +\phi(\R_1,\cdots\R_N)
\eeq

Typical lattice vibrations involve small atomic excursions of the order
$0.1\A$ or smaller, thus we may expand about the equilibrium position of 
the ions.
\beq
\phi(\{r_{n\alpha i}+ s_{n\alpha i}\})=
\phi(\{r_{n\alpha i}\}) +
\pde{\phi}{r_{n\alpha i}}s_{n\alpha i} +
\frac{1}{2}\spdm{\phi}{r_{n\alpha i}}{r_{m\beta j}}s_{n\alpha i}s_{m\beta j}
\eeq
The first two terms in the sum are zero; the first by definition, and
the second is zero since it is the first derivative of a potential 
being evaluated at the equilibrium position.  We will define the matrix
\beq
\Phi_{n\alpha i}^{m\beta j}=\spdm{\phi}{r_{n\alpha i}}{r_{m\beta j}}
\eeq
From the different conservation laws (related to symmetries) of the system 
one may derive some simple relationships for $\Phi$.  We will discuss these 
in detail later.  
\begin{figure}[htb]
\centerline{\includegraphics[width=0.85\textwidth,clip=true]{periodicpot.pdf}}
\caption[]{\em{Since the coefficients of potential between the atoms linked by 
the blue lines (or the red lines) must be identical, 
$\Phi_{n\alpha i}^{m\beta j}=\Phi_{0\alpha i}^{(m-n)\beta j}$.}}
\label{fig:periodicpot}
\end{figure}
However, one must be introduced now, that
is, due to translational invariance.
\beq
\Phi_{n\alpha i}^{m\beta j}=\Phi_{0\alpha i}^{(m-n)\beta j}
=\spdm{\phi}{r_{0\alpha i}}{r_{(n-m)\beta j}}
\eeq
ie, it can only depend upon the distance.  This is important for
the next subsection.


\subsection{The Equation of Motion}

	From the derivative of the potential, we can calculate the
force on each site
\beq
F_{n\alpha i}= -\pde{\phi(\{r_{m\beta j}+s_{m\beta j}\})}{s_{n\alpha i}}
\eeq
so that the equation of motion is
\beq
-\Phi_{n\alpha i}^{m\beta j} s_{m\beta j} = M_\alpha \ddot{s}_{n\alpha i}
\eeq
If there are $N$ unit cells, each with $r$ atoms, then this gives
$3Nr$ equations of motion.  We will take advantage of the periodicity of 
the lattice by using Fourier transforms to achieve a significant 
decoupling of these equations.  Imagine that the coordinate $s$ of each
site is decomposed into its Fourier components.  Since the equations
are linear, we may just consider one of these components to derive our
equations of motion in Fourier space
\beq
s_{n\alpha i}=\frac{1}{\sqrt{M_\alpha}} u_{\alpha i}(\q)
e^{i(\q\cdot\r_n-\omega t)}
\label{eq:displacement}
\eeq
where the first two terms on the rhs serve as the polarization vector
for the oscillation, $u_{\alpha i}(\q)$ is independent of $n$ due
to the translational invariance of the system. 
\begin{figure}[htb]
\centerline{\includegraphics[width=0.7\textwidth,clip=true]{propagate.pdf}}
\caption[]{\em{$u_{\alpha i}(\q)$ is independent of $n$ so that a lattice 
vibration can propagate and respect the translational invariance of the 
lattice.}}
\label{fig:propagate}
\end{figure}
In a real system the 
real $s$ would be composed of a sum over all $\q$ and polarizations.
With this substitution, the equations of motion become
\beq
\omega^2 u_{\alpha i}(\q) 
= 
\frac{1}{\sqrt{M_\alpha M_\beta}}
\Phi_{n\alpha i}^{m\beta j}
e^{i\q\cdot(\r_m-\r_n)}
u_{\beta j}(\q)
\;\;\;\mbox{sum repeated indices}\,.
\eeq
Recall that $\Phi_{n\alpha i}^{m\beta j}=\Phi_{0\alpha i}^{(m-n)\beta j}$
so that if we identify
\beq
D_{\alpha i}^{\beta j} 
=
\frac{1}{\sqrt{M_\alpha M_\beta}}
\Phi_{n\alpha i}^{m\beta j}
e^{i\q\cdot(\r_m-\r_n)}
=
\frac{1}{\sqrt{M_\alpha M_\beta}}
\Phi_{0\alpha i}^{p\beta j}
e^{i\q\cdot(\r_p)}
\eeq
where $\r_p=\r_m-\r_n$, then the equation of motion becomes
\beq
\omega^2 u_{\alpha i}(\q) 
= 
D_{\alpha i}^{\beta j}
u_{\beta j}(\q)
\eeq
or
\beq
\lep D_{\alpha i}^{\beta j}- \omega^2\delta_{\alpha i}^{\beta j}\rip
u_{\beta j}(\q) 
=
0
\label{eq:secular}
\eeq
which only has nontrivial ($u\neq 0$) solutions if 
$ \det \lep \D(\q)-\omega^2\I \rip =0$.  For each $\q$ there are 
$3r$ different solutions (branches) with eigenvalues $\omega^{(n)}(\q)$
(or rather $\omega^{(n)}(\q)$ are the root of the eigenvalues).
The dependence of these eigenvalues $\omega^{(n)}(\q)$ on $\q$ is known
as the dispersion relation.

\subsection{Example, a Linear Chain}


\begin{figure}[htb]
\centerline{\includegraphics[width=0.8\textwidth,clip=true]{M1M2.pdf}}
\caption[]{\em{A linear chain of oscillators composed of a two-element
basis with different masses, $M_1$ and $M_2$ and equal strength springs
with spring constant $f$.}}
\label{fig:M1M2}
\end{figure}

Consider a linear chain of oscillators composed of a two-element
basis with different masses, $M_1$ and $M_2$ and equal strength springs
with spring constant $f$.  It has the potential energy
\beq
\phi=\frac12 f\sum_n 
\lep s_{n,1}-s_{n,2}\rip^2 + \lep s_{n,2}-s_{n+1,1}\rip^2\,.
\eeq
We may suppress the indices $i$ and $j$, and search for a solution
\beq
s_{n\alpha}=\frac{1}{\sqrt{M_\alpha}} u_{\alpha}(\q)
e^{i(\q\cdot\r_n-\omega t)}
\eeq
to the equation of motion
\beq
\omega^2 u_{\alpha}(\q) = D_{\alpha}^{\beta} u_{\beta}(\q)
\;\;
\mbox{where}
\;\;
D_{\alpha}^{\beta} 
=
\frac{1}{\sqrt{M_\alpha M_\beta}}
\Phi_{0\alpha}^{p,\beta}
e^{i\q\cdot(\r_p)}
\eeq
and, 
\beq
\Phi_{n,\alpha}^{m,\beta}
=\spdm{\phi}{r_{0,\alpha}}{r_{(n-m),\beta}}
\eeq
where nontrivial solutions are found by solving $ \det \lep \D(\q)-\omega^2\I \rip =0$.

	The potential matrix has the form
\beq
\Phi_{n,1}^{n,1}=\Phi_{n,2}^{n,2}=2f
\eeq
\beq
\Phi_{n,1}^{n,2}=\Phi_{n,2}^{n,1}=\Phi_{n,1}^{n-1,2}=\Phi_{n,2}^{n+1,1}=-f\,.
\eeq
This may be Fourier transformed on the space index $n$ by inspection, so that
\beqa
D_{\alpha}^{\beta} 
&=&
\frac{1}{\sqrt{M_\alpha M_\beta}}
\Phi_{0\alpha}^{p\beta}
e^{i\q\cdot(\r_p)}\nonumber \\
&=&
\lep\begin{array}{cc}
\frac{2f}{M_1}	& -\frac{f}{\sqrt{M_1M_2}}\lep 1+e^{-iqa}\rip \\
-\frac{f}{\sqrt{M_1M_2}}\lep 1+e^{+iqa}\rip	& \frac{2f}{M_2} 
\end{array}\rip
\eeqa
Note that the matrix $\D$ is hermitian, as it must be to yield real, 
physical, eigenvalues $\omega^2$ (however, $\omega$ can still be imaginary
if $\omega^2$ is negative, indicating an unstable mode).  The secular
equation $ \det \lep \D(\q)-\omega^2\I \rip =0$ becomes
\beq
\omega^4 -\omega^2 2f\lep\frac{1}{M_1}+\frac{1}{M_2}\rip
+\frac{4f^2}{M_1M_2}\sin^2(qa/2)=0\,,
\eeq
with solutions
\beq
\omega^2= f\lep\frac{1}{M_1}+\frac{1}{M_2}\rip
\pm 
f\sqrt{
\lep\frac{1}{M_1}+\frac{1}{M_2}\rip^2 -\frac{4}{M_1M_2}\sin^2(qa/2)
}
\eeq
\begin{figure}[htb]
\centerline{\includegraphics[width=0.7\textwidth,clip=true]{dispersion.pdf}}
\caption[]{\em{Dispersion of the linear chain of oscillators shown in 
Fig.~\ref{fig:M1M2} when $M_1=1$, $M_2=2$ and $f=1$.  The upper branch $\omega_+$
is called the optical and the lower branch is the acoustic mode.}}
\label{fig:dispersion}
\end{figure}
This equation simplifies significantly in the $q\to 0$ and $q/a\to\pi$ limits.
In units where $a=1$, and where the reduced mass 
$1/\mu=\lep\frac{1}{M_1}+\frac{1}{M_2}\rip$,
\beq
\lim_{q\to 0} \omega_{-}(q)=qa\sqrt{\frac{f\mu}{2M_1M_2}}
\;\;\;\;\;\;
\lim_{q\to 0} \omega_{+}(q)=\sqrt{\frac{2f}{\mu}}
\eeq
and 
\beq
\omega_-(q=\pi/a)= \sqrt{2f/M_2}\,.
\;\;\;\;\;\;
\omega_+(q=\pi/a)= \sqrt{2f/M_1}
\eeq
As a result, the $+$ mode is quite flat; whereas the $-$ mode varies from
zero at the Brillouin zone center $q=0$ to a flat value at the edge
of the zone.  This behavior is plotted in Fig.~\ref{fig:dispersion}.

	It is also instructive to look at the eigenvectors, since they
will tell us how the atoms vibrate.  Let's look at the optical mode
at $q=0$, $\omega_+(0)=\sqrt{2f/\mu}$.  Here,
\beq
\D=\left(\begin{array}{cc}
2f/M_1	& -2f/\sqrt{M_1M_2}\\
-2f/\sqrt{M_1M_2} & 2f/M_2
\end{array}\right)\,.
\eeq
Eigenvectors are non-trivial solutions to $(\omega^2\I-\D)\u=0$,
or
\beq
0=\lep\begin{array}{cc}
2f/\mu-2f/M_1	& 2f/\sqrt{M_1M_2}\\
2f/\sqrt{M_1M_2} & 2f/\mu-2f/M_2
\end{array}\rip
\lep\begin{array}{c}
u_1\\u_2\end{array}\rip
\,.
\eeq
with the solution $u_1=-\sqrt{M_2/M_1}u_2$.  In terms of the actual
displacements Eqs.\ref{eq:displacement}
\beq
\frac{s_{n1}}{s_{n2}}=\sqrt{\frac{M_2}{M_1}}\frac{u_1}{u_2}
\eeq
or $s_{n1}/s_{n2}=-M_2/M_1$
so that the two atoms in the basis are moving out of phase with 
amplitudes of motion inversely proportional to their masses.  
\begin{figure}[htb]
\centerline{\includegraphics[width=0.75\textwidth,clip=true]{M1M2_opt.pdf}}
\caption[]{\em{Optical Mode (bottom) of the linear chain (top).}}
\label{fig:optical}
\end{figure}
These modes are described as optical modes since these atoms, if oppositely 
charged, would form an oscillating dipole which would couple to optical 
fields with $\lambda\sim a$.  Not all optical modes are optically
active.

\subsection{The Constraints of Symmetry}

	We know a great deal about the dispersion of the lattice
vibrations without solving explicitly for them.  For example, we know that
for each $q$, there will be $dr$ modes (where $d$ is the lattice dimension,
and $r$ is the number of atoms in the basis).  We also expect (and implicitly 
assumed above) that the allowed frequencies are real and positive.  However,
from simple mathematical identities, the point-group and translational 
symmetries of the lattice, and its time-reversal invariance, we can learn
more about the dispersion without solving any particular problem.

The basic symmetries that we will employ are
\begin{itemize}
\item The translational invariance of the lattice and reciprocal     
      lattice.
\item The point group symmetries of the lattice and reciprocal       
      lattice.
\item Time-reversal invariance.
\end{itemize}

\subsubsection{Symmetry of the Dispersion}

\paragraph{Complex Properties of the dispersion and Eigenmodes}
First, from the symmetry of the second derivative, one may show that
$\omega^2$ is real.  Recall that the dispersion is determined by the
secular equation $ \det \lep \D(\q)-\omega^2\I \rip =0$, so if
$\D$ is hermitian, then its eigenvalues, $\omega^2$, must be real.
\beqa
D_{\alpha i}^{*\beta j} 
&=&
\frac{1}{\sqrt{M_\alpha M_\beta}}
\Phi_{0\alpha i}^{p\beta j}
e^{-i\q\cdot(\r_p)}\\
&=& 
\frac{1}{\sqrt{M_\alpha M_\beta}}
\Phi_{0,\alpha,i}^{-p,\beta, j}
e^{i\q\cdot(\r_p)}
\eeqa
Then, due to the symmetric properties of the second derivative 
\beq
D_{\alpha i}^{*\beta j} 
= 
\frac{1}{\sqrt{M_\alpha M_\beta}}
\Phi_{-p,\beta, j}^{0,\alpha,i}
e^{i\q\cdot(\r_p)}=
\frac{1}{\sqrt{M_\alpha M_\beta}}
\Phi_{0,\beta, j}^{p,\alpha,i}
e^{i\q\cdot(\r_p)}=
D_{\beta j}^{\alpha i}
\eeq
Thus, $\D^{T*}=\D^{\dag}=\D$ so $\D$ is hermitian and its eigenvalues 
$\omega^2$ are real.  This 
means that either $\omega$ are real or they are pure imaginary.  We will
assume the former.  The latter yields pure exponential growth of our 
Fourier solution, indicating an instability of the lattice to a second-order
structural phase transition.

Time-reversal invariance allows us to show related results.
We assume a solution of the form
\beq
s_{n\alpha i}=\frac{1}{\sqrt{M_\alpha}} u_{\alpha i}(\q)
e^{i(\q\cdot\r_n-\omega t)}
\eeq
which is a plane wave.  Suppose that the plane wave is moving to the
right so that $\q=\xh q_x$, then the plane of stationary
phase travels to the right with
\beq
x=\frac{\omega}{q_x}t\,.
\eeq
Clearly then changing the sign of $q_x$ is equivalent to taking 
$t\to -t$.  If the system is to display proper time-reversal invariance,
so that the plane wave retraces its path under time-reversal, it must have 
the same frequency when time, and hence $\q$, is reversed, so
\beq
\omega(-\q)=\omega(\q)\,.
\eeq
Note that this is fully equivalent to the statement that 
$D_{\beta j}^{\alpha i}(\q) = D_{\beta j}^{*\alpha i}(-\q)$
which is clear from the definition of $\D$.

Now, return to the secular equation, Eq.~\ref{eq:secular}.  
\beq
\lep D_{\alpha i}^{\beta j}(\q)- \omega^2(\q)\delta_{\alpha i}^{\beta j}\rip
\epsilon_{\beta j}(\q) 
=
0
\eeq
Lets call the (normalized) eigenvectors of this equation $\epsilon$.  They 
are the elements of a unitary matrix which diagonalizes $\D$.  As a
result, they have orthogonality and completeness relations
\beq
\sum_{\alpha,i}
\epsilon_{\alpha,i}^{*(n)}(\q)\epsilon_{\alpha,i}^{(m)}(\q)=\delta_{m,n}
\;\;\;\;\mbox{orthogonality}
\eeq
\beq
\sum_{n}
\epsilon_{\alpha,i}^{*(n)}(\q)\epsilon_{\beta,j}^{*(n)}(\q)=
\delta_{\alpha,\beta}\delta_{i,j}
\eeq
If we now take the complex conjugate of the secular equation
\beq
\lep D_{\alpha i}^{\beta j}(-\q)- \omega^2(-\q)\delta_{\alpha i}^{\beta j}\rip
\epsilon_{\beta j}^*(\q) 
=
0
\eeq
Then it must be that
\beq
\epsilon_{\beta j}^*(\q) \propto \epsilon_{\beta j}(-\q)\,.
\eeq
Since the $\{\epsilon\}$ are normalized the constant of proportionality may 
be chosen as one
\beq
\epsilon_{\beta j}^*(\q) = \epsilon_{\beta j}(-\q)\,.
\eeq

\paragraph{Point-Group Symmetry and the Dispersion}

     A point group operation takes a crystal back to an identical
configuration.  Both the original and final lattice must have the 
same dispersion.  Thus, since the reciprocal lattice has the same point group 
as the real lattice, the dispersion relations have the same point group 
symmetry as the lattice.

For example, the dispersion must share the periodicity of the Brillouin zone.
From the definition of $\D$
\beq
D_{\alpha i}^{\beta j}(\q) 
=
\frac{1}{\sqrt{M_\alpha M_\beta}}
\Phi_{0\alpha i}^{p\beta j}
e^{i\q\cdot(\r_p)}
\eeq
it is easy to see that $D_{\alpha i}^{\beta j}(\q+\G)=
D_{\alpha i}^{\beta j}(\q)$ (since $\G\cdot\r_p=2\pi n$, where $n$ is 
an integer).  I.e., $\D$ is periodic in k-space, and so its eigenvalues 
(and eigenvectors)  must also be periodic.
\beq
\omega^{(n)}(\k+\G)=\omega^{(n)}(\k)
\eeq
\beq
\epsilon_{\beta j}(\k+\G) = \epsilon_{\beta j}(\k)\,.
\eeq

\subsubsection{Symmetry and the Need for Acoustic modes}  

Applying basic symmetries, we can show that an elemental lattice 
(that with $r=1$) must have an acoustic model.  
\begin{figure}[htb]
\centerline{\includegraphics[width=0.49\textwidth,clip=true]{shift.pdf}}
\caption[]{\em{If each ion is shifted by $s_{1,1, i}$, then the lattice
energy is unchanged.}}
\label{fig:shift}
\end{figure}
First, look at the translational invariance of $\Phi$ .  Suppose we make 
an overall shift of the lattice by an arbitrary displacement 
$s_{n,\alpha, i}$ for all sites $n$ and elements of the basis $\alpha$
(i.e.\ $s_{n,\alpha, i}=s_{1,1, i}$).  
Then, since the interaction is
only {\em{between}} ions, the energy of the system should remain unchanged.
\beqa
\delta E&=& \frac12\sum_{m,n,\alpha,\beta,i,j} \Phi_{0,\alpha,i}^{m-n,\beta,j} 
s_{n,\alpha, i}s_{m,\beta, j} =0 \\
&=& \frac12 \sum_{mn\alpha,\beta,i,j}\Phi_{0,\alpha,i}^{m-n,\beta,j} 
s_{1,1, i}s_{1,1, j}\\
&=& \frac12 \sum_{i,j}s_{1,1, i}s_{1,1, j}\sum_{mn\alpha,\beta}
\Phi_{0,\alpha,i}^{m-n,\beta,j} 
\eeqa
Since we know that $s_{1,1, i}$ is finite, it must be that
\beq
\sum_{m,n,\alpha,\beta} \Phi_{0,\alpha,i}^{m-n,\beta,j}=
\sum_{p,\alpha,\beta}\Phi_{0,\alpha,i}^{p,\beta,j}=0
\eeq

Now consider a strain on the system $V_{m,\beta,j}$, described
by the strain matrix $m_{\beta,j}^{\alpha,i}$
\beq
V_{m,\beta,j}=\sum_{\alpha,i}m_{\beta,j}^{\alpha,i}s_{m,\alpha, i}
\eeq
\begin{figure}[htb]
\centerline{\includegraphics[width=0.75\textwidth,clip=true]{squish.pdf}}
\caption[]{\em{After a stress is applied to a lattice, the movement
of each ion (strain) is not only in the direction of the applied
stress.  The response of the lattice to an applied stress is
described by the strain matrix.}}
\label{fig:squish}
\end{figure}
After the stress has been applied, the atoms in the bulk of the sample are 
again in equilibrium (those on the surface are maintained in equilibrium
by the stress), and so the net force must be zero.  Looking at the
central ($n=0$) atom this means that
\beq
0=F_{0,\alpha,i}=-\sum_{m,\beta,j,\gamma,k}\Phi_{0,\alpha,i}^{m,\beta,j}m_{\beta,j}^{\gamma,k}s_{m,\gamma, k}
\eeq
Since this applies for an arbitrary strain matrix $m_{\beta,j}^{\gamma,k}$, 
the coefficients for each $m_{\beta,j}^{\gamma,k}$ must be zero
\beq
\sum_{m}\Phi_{0,\alpha,i}^{m,\beta,j}s_{m,\gamma, k}=0
\eeq
An alternative way (cf.\ Callaway) to show this is to recall that the 
reflection symmetry of the lattice requires that $\Phi_{0,\alpha,i}^{m,\beta,j}$
be even in $m$; whereas, $s_{m,\gamma, k}$ is odd in $m$.  Thus the
sum over all $m$ yields zero.

Now let's apply these constraints to $\D$ for an elemental lattice where
$r=1$, and we may suppress the basis indices $\alpha$.
\beq
D_{i}^{j}(\q) 
=
\frac{1}{M}\sum_p
\Phi_{0, i}^{p, j}
e^{i\q\cdot(\r_p)}
\eeq
For small $q$ we may expand $D$
\beq
D_{i}^{j}(\q) 
=
\frac{1}{M}\sum_p
\Phi_{0, i}^{p, j}
\lep 1+ i\q\cdot(\r_p) -\frac12 \lep \q\cdot(\r_p)\rip^2 +\cdots \rip
\eeq
We have shown above that the first two terms in this series are zero.
Thus,
\beq
D_{i}^{j}(\q) 
\approx
-\frac{1}{2M}\sum_p
\Phi_{0, i}^{p, j}
\lep i\q\cdot(\r_p)\rip^2
\eeq
Thus, the leading order (small $\q$) eigenvalues $\omega^2(\q)\sim q^2$.  
I.e. they are acoustic modes.  We have shown that all elemental lattices
must have acoustic modes for small $\q$.

	In fact, one may show that all harmonic lattices in which the
energy is invariant under a rigid translation of the entire lattice must 
have at least one acoustic mode.  We will not prove this, but rather
make a simple argument.  The rigid translation of the lattice corresponds
to a $\q=0$ translational mode, since no energy is gained by this translation,
it must be that $\om_s(\q=0)=0$ for the branch $s$ which contains this 
mode.  The acoustic mode may be obtained by perturbing (in $\q$) around
this point.  Physically this mode corresponds to all of the elements of 
the basis moving together so as to emulate the motion in the elemental 
basis.

\section{The Counting of Modes}

In the sections to follow, we need to perform sums (integrals) of functions
of the dispersion over the crystal momentum states $\k$ within the reciprocal 
lattice. However, the translational and point group symmetries of the crysal, 
often greatly reduce the set of points we must sum.  In addition, we often
approximate very large systems with hypertoroidal models with periodic
boundary conditions.  This latter approximation becomes valid as the 
system size diverges so that the surface becomes of zero measure.

\subsection{Periodicity and the Quantization of States}

A consequence of approximating our system as a finite-sized periodic system
is that we now have a discrete sum  rather than an integral over $\k$.  
Consider a one-dimensional finite system with $N$ atoms and periodic boundary 
conditions. We seek solutions to the phonon problem of the type
\beq
s_n=\epsilon(q)e^{i(qr_n-\omega t)}\;\;\;\mbox{where }r_n=na
\eeq
and we require that
\beq
s_{n+N}=s_n
\eeq
or 
\beq
q(n+N)a=qna+2\pi m\;\;\;\mbox{where $m$ is an integer}
\eeq
Then, the allowed values of $q=2\pi m/Na$.  This will allow us to convert
the integrals over the Brillouin zone to discrete sums, at least for cubic
systems;  however, the method is easily generalized for other Bravais
lattices.

\subsection{Translational Invariance: First Brillouin Zone}

We can use the translational invariance of the crystal to reduce
the complexity of sums or integrals of functions of the dispersion over 
the crystal momentum states. As shown above, translationally invariant
systems have states which are not independent. It is useful then to define 
a region of k-space which contains only independent states.  Sums over k 
may then be confined to this region. This region is defined as the smallest 
polyhedron centered at the origin of the reciprocal lattice and enclosed by 
perpendicular bisectors of the G's is called the Brillouin zone 
(cf.~Fig.~\ref{fig:firstbz}).
\begin{figure}
\centerline{\includegraphics[width=0.7\textwidth,clip=true]{firstbz.pdf}}
\caption[]{\em{
The First Brillouin Zone.  The end points of all vector pairs that 
satisfy the Bragg condition $\k-\k_0=\G_{hkl}$ lie on the perpendicular 
bisector of $\G_{hkl}$. The smallest polyhedron centered at the origin 
of the reciprocal lattice and enclosed by perpendicular bisectors of 
the G's is called the first Brillouin zone.
}}
\label{fig:firstbz}
\end{figure}
Typically, we choose to include only half of the bounding surface
within the first Brillouin zone, so that it can also be defined as the
set of points which contains only independent states.

From the discussion in chapter 3 and in this chapter, it is also clear that 
the reciprocal lattice vectors have some interpretation as momentum.  For 
example, the Laue condition requires that the change in momentum of the 
scatterer be equal to a reciprocal lattice translation vector.  The end 
points of all vector pairs that satisfy the Bragg condition $\k-\k_0=\G_{hkl}$ 
lie on the perpendicular bisector of $\G_{hkl}$.  Thus, the FBZ is also the 
set of points which cannot satisfy the Bragg condition.

\subsection{Point Group Symmetry and Density of States}

	Two other tricks to reduce the complexity of these sums are
worth mentioning here although they are discussed in detail elswhere.

	The first is the use of the point group symmetry of the system.
It is clear from their definition in chapter 3, the reciprocal lattice
vectors have the same point group symmetry as the lattice.  As we discussed
in chapter 2, the knowledge of the group elements and corresponding
degeneracies may be used to reduce the sums over $\k$ to the irreducible
wedge within the the First Brillioun zone.  For example, for a cubic
system, this wedge is only $1/2^33!$ or $1/48$th of the the FBZ!

	The second is to introduce a phonon density of states to reduce the
multidimensional sum over $k$ to a one-dimensional integral over energy.
This will be discussed in chapter 5.

\section{Normal Modes and Quantization}

In this section we will derive the equations of motion for the lattice,
determine the canonically conjugate variables (the the sense of Lagrangian
mechanics), and use this information to both first and second quantize
the system.

Any lattice displacement may be expressed as a sum over the eigenvectors 
of the dynamical matrix $\D$.
\beq
s_{n,\alpha,i} =\frac{1}{\sqrt{M_\alpha N}}\sum_{\q,s}
Q_s(\q,t) \ep_{\alpha,i}^s(\q)e^{i\q\cdot\r_n}
\eeq
Recall that $\ep_{\alpha,i}^s(\q)$ are distinguished from 
$\u_{\alpha,i}^s(\q)$ only in that they are normalized.  Also since 
$\q+\G$ is equivalent to $\q$, we need sum only over the first Brillouin 
zone.  Finally we will assume that $Q_s(\q,t)$ contains the harmonic time 
dependence and since $s_{n,\alpha,i}$ is real $Q_s^*(\q)=Q_s(-\q)$.

We may rewrite both the kinetic and potential energy of the system as 
sums over $Q$.  For example, the kinetic energy of the lattice
\beqa
 T & = & \frac12 \sum_{n,\alpha,i}M_\alpha\lep\dot{s}_{n\alpha,i}\rip^2\\
   & = & \frac{1}{2N}\sum_{n,\alpha,i}\sum_{\q,\k,r,s}
\dot{Q}_r(\q)\ep_{\alpha,i}^r(\q)e^{i\q\cdot\r_n}
\dot{Q}_s(\k)\ep_{\alpha,i}^s(\k)e^{i\k\cdot\r_n}
\eeqa
Then as
\beq
\frac{1}{N}\sum_n e^{i(\k+\q)\cdot\r_n}=\delta_{\k,-\q}
\;\;\;\mbox{and}\;\;\;
\sum_{\alpha,i} \ep_{\alpha,i}^r\ep_{\alpha,i}^{*s}=\delta_{rs}
\eeq
the kinetic energy may be reduced to
\beq
T=\frac12 \sum_{\q,r}\left|\dot{Q}_r(\q)\right|^2
\eeq
The potential energy may be rewritten in a similar fashion
\beqa
V&=&\frac12 \sum_{n,m,\alpha,\beta,i,j}
\Phi_{n,\alpha,i}^{m,\beta,j} s_{n,\alpha,i}s_{m,\beta,j} \nonumber\\
 &=& \frac12 \sum_{n,m,\alpha,\beta,i,j} 
\frac{\Phi_{0,\alpha,i}^{m-n,\beta,j}}{N\sqrt{M_\alpha M_\beta}}\nonumber\\
&&\sum_{\q,\k,s,r}
Q_s(\q,t) \ep_{\alpha,i}^s(\q)e^{i\q\cdot\r_n}
Q_r(\k,t) \ep_{\beta,j}^r(\k)e^{i\k\cdot\r_m}
\eeqa
Let $\r_l=\r_m-\r_n$
\beqa
V&=&\frac12 \sum_{n,l,\alpha,\beta,i,j} 
\frac{\Phi_{0,\alpha,i}^{l,\beta,j}}{N\sqrt{M_\alpha M_\beta}}\nonumber\\
&&\sum_{\q,\k,s,r}
Q_s(\q,t) \ep_{\alpha,i}^s(\q)e^{i\q\cdot\r_n}
Q_r(\k,t) \ep_{\beta,j}^r(\k)e^{i\k\cdot(\r_l+\r_n)}
\eeqa
and sum over $n$ to obtain the delta function $\delta_{\k,-\q}$
so that
\beq
V=\frac12 \sum_{l,\alpha,\beta,i,j,s,r}
Q_s(-\k)\ep_{\alpha,i}^s(-\k)Q_r(\k)\ep_{\beta,j}^r(\k)
\frac{1}{\sqrt{M_\alpha M_\beta}}
\Phi_{0,\alpha,i}^{l,\beta,j}e^{i\k\cdot\r_l}\,.
\eeq
Note that the sum over $l$ on the last three terms yields $\D$,
so that
\beq
V=\frac12 \sum_{l,\alpha,\beta,i,j,s,r} D_{\alpha,i}^{\beta,j}(\k)
Q_s(-\k)\ep_{\alpha,i}^s(-\k)Q_r(\k)\ep_{\beta,j}^r(\k)\,.
\eeq
Then, since $\sum_{\beta,j}D_{\alpha i}^{\beta j}(\k)\ep_{\beta j}^r(\k)= 
\omega^2_r(\k)\ep_{\alpha,i}^r(\k)$ and $\ep_{\alpha,i}^s(-\k)=
\ep_{\alpha,i}^{*s}(\k)$,
\beq
V=\frac12 \sum_{\alpha,i,\k,r,s}\ep_{\alpha,i}^r(\k)\ep_{\alpha,i}^{*s}(\k)
\omega_r^2(\k)Q_s^*(\k)Q_r(\k)
\eeq
Finally, since 
$\sum_{\alpha,i}\ep_{\alpha,i}^r(\k)\ep_{\alpha,i}^{*s}(\k)=\delta_{r,s} $
\beq
V=\frac12 \sum_{\k,s} \omega_s^2(\k) \left|Q_s(\k)\right|^2
\eeq
Thus we may write the Lagrangian of the ionic system as
\beq
L=T-V=\frac12 \sum_{\k,s}
\lep
\left|\dot{Q}_s(\k)\right|^2
- 
\omega_s^2(\k) \left|Q_s(\k)\right|^2
\rip\,,
\eeq
where the $Q_s(\k)$ may be regarded as canonical coordinates, and
\beq
P_r^*(\k)=\pde{L}{Q_r(\k)}=\dot{Q}_s^*(\k)
\eeq
(no factor of 1/2 since $Q_s^*(\k)=Q_s(-\k)$) are the canonically
conjugate momenta.

The equations of motion are
\beq
\frac{d}{dt}\lep\pde{L}{\dot{Q}^*_s(\k)}\rip-\pde{L}{{Q}^*_s(\k)}
\;\;\;\mbox{or}\;\;\; 
\ddot{Q}_s(\k)+\om_s^2(\k)Q_s(\k)=0
\eeq
for each $\k,s$.  These are the equations of motion for $3rN$ independent 
harmonic oscillators.  Since going to the Q-coordinates accomplishes the
decoupling of these equations, the $\{Q_s(\k)\}$ are referred to as normal coordinates.

\subsection{Quantization and Second Quantization}

P.A.M. Dirac laid down the rules of quantization, from Classical 
Hamilton-Jacobi classical mechanics to Hamiltonian-based quantum mechanics
following the path (Dirac p.84-89):
\begin{enumerate}
\item First, identify the classical canonically conjugate set of variables 
$\{q_i,p_i\}$
\item These have Poisson Brackets
\beq
\{\{u,v\}\}=\sum_i\lep\pde{u}{q_i}\pde{v}{p_i}-\pde{u}{p_i}\pde{v}{q_i}\rip
\eeq
\beq
\{\{q_i,p_j\}\}=\delta_{i,j}\;\;\; \{\{p_i,p_j\}\}=\{\{q_i,q_j\}\}=0
\eeq
\item Then define the quantum Poisson Bracket (the commutator)
\beq
[u,v]=uv-vu=i\hbar\{\{u,v\}\}
\eeq
\item In particular, $[q_i,p_j]=i\hbar\delta_{i,j}$, and
$[q_i,q_j]=[p_i,p_j]=0$.
\end{enumerate}

Thus, following Dirac, we may now quantize the normal coordinates
\beq
[Q_r^*(\k),P_s(\q)]=i\hbar\delta_{\k,\q}\delta_{r,s}
\;\;\;\mbox{where the other commutators vanish}\,.
\label{eq:qcom}
\eeq
Furthermore, since we have a system of $3rN$ uncoupled harmonic
oscillators we may immediately second quantize by introducing
\beqa
a_s(\k)&=&\frac{1}{\sqrt{2\hbar}}\lep\sqrt{\om_s(\k)}Q_s(\k)+
\frac{i}{\sqrt{\om_s(\k)}}P_s(\k)\rip\\
a_s^{\dag}(\k)&=&\frac{1}{\sqrt{2\hbar}}\lep\sqrt{\om_s(\k)}Q_s^*(\k)-
\frac{i}{\sqrt{\om_s(\k)}}P_s^*(\k)\rip\,,
\eeqa
or
\beqa
Q_s(\k)&=&\sqrt{\frac{\hbar}{2\om_s(\k)}}\lep a_s(\k)+a_s^{\dag}(-\k)\rip\\
P_s(\k)&=&-i\sqrt{\frac{\hbar\om_s(\k)}{2}}\lep a_s(\k)-a_s^{\dag}(-\k)\rip
\eeqa
Where
\beq
[a_s(\k),a_r^{\dag}(\q)]=\delta_{r,s}\delta_{\q,\k}
\;\;\;[a_s(\k),a_r(\q)]=[a_s^{\dag}(\k),a_r^{\dag}(\q)]=0
\label{eq:acom}
\eeq
This transformation $\{Q,P\}\to\{a,a^{\dag}\}$ is canonical, since is
preserves the commutator algebra Eq.~\ref{eq:qcom}, and the Hamiltonian
becomes
\beq
H=\sum_{k,s}\hbar\om_s(\k)\lep a_s^{\dag}(\k)a_s(\k) + \frac12 \rip
\eeq
which is a sum over $3rN$ independent quantum oscillators, each one
referred to as a phonon mode!

The number of phonons in state $(\k,s)$ is given by the operator
\beq
n_s(\k)=a_s^{\dag}(\k)a_s(\k)
\eeq
and $a_s^{\dag}(\k)$ and $a_s(\k)$  create and destroy phonons respectively, 
in the state $(\k,s)$
\beqa
a_s^{\dag}(\k)
\left| n_s(\k)\right> & = & \sqrt{n_s(\k)+1}\left|n_s(\k)+1\right>\\
a_s(\k)\left|n_s(\k)\right> & = & \sqrt{n_s(\k)}\left|n_s(\k)-1\right>
\eeqa
If $\left|0\right>$ is the normalized state with no phonons present,
then the state with $\{n_s(\k)\}$ phonons in each state $(\k,s)$ is
\beq
\left|\{n_s(\k)\}\right>
=
\leb
\lep\prod_{\k,s}\frac{1}{n_s(\k)!}\rip^{\frac12}
\rib
\prod_{\k,s}\lep a_{s}^{\dag}(\k)\rip^{{n_s(\k)}}\left|0\right>
\eeq
Finally the lattice point displacement
\beq
s_{n,\alpha,i}=\frac{1}{\sqrt{M_\alpha N}}\sum_{\q,s}\sqrt{\frac{\hbar}{2\om_s(\q)}}
\lep a_s(\q)+ a_s^{\dag}(-\q)\rip \ep^s_{\alpha,i}(\q)e^{i\q\cdot\r_n}
\eeq
will be important in the next section, especially with respect to
zero-point motion (i.e.\ $\left< s^2 \right>_{T=0}\neq 0$).

\section{Theory of Neutron Scattering}

\begin{figure}[htb]
\centerline{\includegraphics[width=0.75\textwidth,clip=true]{neutron.pdf}}
\caption[]{\em{Neutron Scattering.  The De Broglie wavelength of the 
neutrons must be roughly the same size as the lattice constants in 
order to learn about the lattice structure and its vibrational modes
from the experiment.  This dictates the use of thermal neutrons.}}
\label{fig:neutron}
\end{figure}

To ``see'' the lattice with neutrons, we want their De Broglie
wavelength $\lambda=h/p$ 
\beq
\lambda_{\rm{neutron}} =\frac{0.29\A}{\sqrt{E}}\;\;\;\;
\mbox{E measured in eV}
\eeq
to be of the same length as the intersite distance on the lattice.
This means that their kinetic energy $E\approx \frac12 Mv^2\approx 0.1$eV,
or $E/k_b\approx 1000K$; i.e.\ thermal neutrons.

     Since the neutron is chargeless, it only interacts with the
atomic nucleus through a short-ranged nuclear interaction
(ignoring any spin-spin interaction).  The range of this interaction 
is ~1 Fermi ($10^{-13}$cm.) or about the radius of the atomic nucleus.  
Thus
\beq
\lambda\sim\A\gg \mbox{range of the interaction}\sim 10^{-13}{\rm{cm.}}
\eeq
     Thus the neutron cannot "see" the detailed structure of the
nucleus, and so we may approximate the neutron-ion interaction
potential as a contact interaction
\beq
V(\r)=\sum_{\r_n} V_n \delta(\r-\r_n)
\eeq
i.e., we may ignore the angular dependence of the scattering
factor f.


\subsection{Classical Theory of Neutron Scattering}


Due to the importance of lattice vibrations, which are inherently
quantum in nature, there is a limit to what we can learn from a
classical theory of diffraction.  Nevertheless it is useful to
compare the classical result to what we will develop for the
quantum problems.

     For the classical problem we will assume that the lattice is elemental
($r=1$) and start with a generalization of the formalism developed in the last 
chapter
\beq
I\propto \left|\rho(\K,\Omega)\right|^2
\eeq
where $\K=\k_0-\k$ and $\Omega=\om_0-\om$.  Furthermore, we take 
\beq
\rho(\r(t))\propto\sum_n\delta(\r-\r_n(t))
\eeq
where
\beq
\r_n(t)=\r_n+\s_n(t)\;\;\mbox{and}
\;\;\s_n(t)=\frac{1}{\sqrt{M}}\u(\q)e^{i(\q\cdot\r_n-\om(\q)t)}
\eeq
describes the harmonic motion of the s-mode with wave-vector q.
\beq
\rho(\K,\Omega)\propto\sum_n\int dt 
e^{i\leb\K\cdot(\r_n+\s_n(t))-\Omega t\rib}\,.
\eeq
For $|\K|\sim 2\pi/\A$ and $\s_n(t)\ll \A$ we may expand
\beq
\rho(\K,\Omega)\propto\sum_n\int dt 
e^{i\leb\K\cdot(\r_n)-\Omega t\rib}\lep 1+i\K\cdot\s_n(t) +\cdots\rip
\eeq
The first term yields a finite contribution only when
\beq
\K=\k_0-\k=\G
\;\;\mbox{and}\;\;
\Omega=\omega_0-\omega=0
\eeq
which are the familiar Bragg conditions for elastic scattering.

     The second term, however, yields something new.  It only
yields a finite result when
\beq
\K\pm\q=\k_0-\k\pm\q=\G
\;\;\mbox{and}\;\;
\Omega\pm\om_s(\q)=\om_0-\om\pm\om_s(\q)=0
\eeq
When multiplied by $\hbar$ , these can be interpreted as conditions for
the conservation of (crystal) momentum and energy when the
scattering event involves the creation (destruction) of a lattice
excitation (phonon).  These processes are called Stokes and antistokes
processes, respectively, and are illustrated in Fig.~\ref{fig:stokes}.
\begin{figure}[htb]
\centerline{\includegraphics[width=0.75\textwidth,clip=true]{stokes.pdf}}
\caption[]{\em{Stokes and antistokes processes in inelastic neutron 
scattering involving the creation or absorption of a lattice phonon.}}
\label{fig:stokes}
\end{figure}

Clearly, the anti-Stokes process can only happen at finite temperatures
where real (as opposed to virtual) phonons are excited.  Thus, our classical
formalism does not correctly describe the temperature dependence of the
scattering.  Several other things are missing in the classical theory, 
including:
\begin{itemize}
\item Security in the validity of the result.
\item The effects of zero-point motion.
\item Correct temperature dependence.
\end{itemize}


\subsection{Quantum Theory of Neutron Scattering}

To address these concerns, we will do a fully quantum calculation.
Several useful references for this calculation include
\begin{itemize}
\item Ashcroft and Mermin, Appendix N, p.\ 790)
\item Callaway, p. 36--.
\item Hook and Hall (for experiment) Ch.~12 p.342-
\end{itemize}

We will imagine that the scattering shown in Fig.~\ref{fig:quantum_neutron}
occurs in a box of volume $V$.
\begin{figure}[htb]
\centerline{\includegraphics[width=0.85\textwidth,clip=true]{quantum_neutron.pdf}}
\caption[]{\em{The initial (left) and final (right) states of the
neutron and lattice during a scattering event.  The initial system state
is given by $\Psi_0=\phi_0\psi_0$, with energy $\ep_0=E_0+\hbar\om_0$
where $\om_0=k_0^2/2M$.  The final system state is given by 
$\Psi_f=\phi_f\psi_f$, with energy $\ep_f=E_f+\hbar\om_f$
where $\om_f=k_f^2/2M$. }}
\label{fig:quantum_neutron}
\end{figure}
The momentum transfer, from the neutron to the lattice is $\K=\k_0-\k_f$ and
the energy transfer which is finite for inelastic scattering is 
$\hbar\Omega=\hbar\lep\om_0-\om_f\rip$.  Again we will take the neutron-lattice 
interaction to be local
\beq
V(\r)=\sum_{r_n} V(\r-\r_n)=\frac{1}{N}\sum_{\q,n}V(\q)e^{i\q\cdot(\r-\r_n)}
=\int\frac{d^3q}{V} V_0\sum_n e^{i\q\cdot(\r-\r_n)}
\eeq
where the locality of the interaction ($V(\r-\r_n)\propto\delta(\r-\r_n)$)
indicates that $V(\q)=V(0)=V_0$.  Consistent with Aschcroft and Mermin, we 
will take
\beq
V(\r) =\frac{2\pi\hbar^2 a}{M} \frac{1}{V}\sum_{\r_n}
\int d^3qe^{i\q\cdot(\r-\r_n)}
\label{eq:ion-neutron}
\eeq
where $a$ is the scattering length, and $V_0=\frac{2\pi\hbar^2 a}{M}$ is 
chosen such that the total cross section $\sigma =4\pi a^2$.

To formulate our quantum theory, we will use Fermi's golden rule
for time dependent perturbation theory.  (This is fully
equivalent to the lowest-order Born approximation).  The
probability per unit time for a neutron to scatter from state $\k_0$ to
$\k_f$ is given by
\beqa
P&=&\frac{2\pi}{\hbar}\sum_f\delta(\ep_0-\ep_f)\left|\left<\Psi_0\left|V\right|
\Psi_f\right>\right|^2\\
&=& \frac{2\pi}{\hbar}\sum_f\delta(E_0+\hbar\om_0-E_f-\hbar\om_f)\nonumber\\
&&\left|\frac{1}{V}\int d^3r e^{i(\k-\k_0)\cdot\r}
\left<\phi_0\left|V(\r)\right|\phi_f\right>\right|^2
\eeqa
If we now substitute in the ion-neutron potential Eq.~\ref{eq:ion-neutron},
then the integral over $\r$ will yield a delta function 
$V\delta(\q+\k-\k_0)$ which allows the $\q$ integral to be evaluated
\beq
P=a^2\frac{(2\pi\hbar)^3}{(MV)^2}\sum_f\delta(E_0-E_f+\hbar\Omega)
\left|\sum_{\r_n}
\left<\phi_0\left|e^{-i\K\cdot\r_n}\right|\phi_f\right>\right|^2
\eeq

     Now, before proceeding to a calculation of the differential cross
section $\frac{d\sigma}{d\Omega d\E}$ we must be able
to convert this probability (rate) for eigenstates into a flux of neutrons
of energy $E$ and momentum $\p$.
A differential volume element of momentum space $d^3p$ contains 
$Vd^3p/(2\pi\hbar)^3$ neutron states.  While this is a natural consequence 
of the uncertainty principle, it is useful to show this in a more
quantitative sense:  Imagine a cubic volume $V = L^3$ with periodic
boundary conditions so that for any state $\Psi$ in $V$, 
\beq
\Psi(x+L,y,z)=\Psi(x,y,z)
\eeq
If we write $\Psi(\r)=\frac{1}{N}\sum_{\q}e^{i\q\cdot\r}\Psi(\q)$, 
then it must be 
that
\beq
q_x L =2\pi m\;\;\;\mbox{where $m$ is an integer}
\eeq
with similar relations for the $y$ and $z$ components.
So for each volume element of q-space $\lep\frac{2\pi}{L}\rip^3$ there is 
one such state.
In terms of states $\p=\hbar\q$, the volume of a state is
$(2\pi\hbar/L)^3$.  Thus $d^3p$
contains $Vd^3p/(2\pi\hbar)^3$ states.

     The incident neutron flux of states (velocity times density) is
\beq
\j=\frac{\hbar\k_0}{M}\left|\Psi_0\right|^2=
\frac{\hbar\k_0}{M}\left|\frac{1}{\sqrt{V}}e^{i\k_0\cdot\r}\right|^2
=\frac{\hbar\k_0}{MV}
\eeq
Then since the number of neutrons is conserved
\beq
j\frac{d\sigma}{dE d\Omega}dE d\Omega = 
\frac{\hbar k_0}{MV}\frac{d\sigma}{dE d\Omega}dE d\Omega=
PV\frac{d^3p}{(2\pi\hbar)^3} = PV\frac{p^2 dp d\Omega}{(2\pi\hbar)^3}
\eeq
And for thermal (non-relativistic) neutrons $E=p^2/2M$, so $dE=pdp/M$, and
\beq
\frac{\hbar k_0}{MV}\frac{d\sigma}{dE d\Omega}dE d\Omega=
 PV\frac{\hbar kM dE d\Omega}{(2\pi\hbar)^3}
\eeq
or
\beq
\frac{d\sigma}{dE d\Omega}=P\frac{k}{k_0}\frac{(MV)^2}{(2\pi\hbar)^3}\,.
\eeq
Substituting in the previous result for $P$
\beq
\frac{d\sigma}{dE d\Omega}=\frac{k}{k_0}\frac{(MV)^2}{(2\pi\hbar)^3}
a^2\frac{(2\pi\hbar)^3}{(MV)^2}\sum_f\delta(E_0-E_f+\hbar\Omega)
\left|\sum_{\r_n}
\left<\phi_0\left|e^{-i\K\cdot\r_n}\right|\phi_f\right>\right|^2
\eeq
or
\beq
\frac{d\sigma}{dE d\Omega}=\frac{k}{k_0}\frac{Na^2}{\hbar} S(\K,\Omega)
\eeq
where
\beq
S(\K,\Omega)=\frac{1}{N}\sum_f\delta(E_0-E_f+\hbar\Omega)
\left|\sum_{\r_n}
\left<\phi_0\left|e^{-i\K\cdot\r_n}\right|\phi_f\right>\right|^2\,.
\eeq
We may deal with the Dirac delta function by substituting
\beq
\delta(\Omega)=\int_{-\infty}^\infty \frac{dt}{2\pi} e^{i\Omega t}\,.
\eeq
so that
\beqa
S(\K,\Omega)&=&\frac{1}{N}\sum_f
\int_{-\infty}^\infty \frac{dt}{2\pi} e^{i\lep (E_0-E_f)/\hbar+\Omega\rip t}
\nonumber\\
&&\sum_{\r_n,\r_m}
\left|\left<\phi_0\left|e^{i\K\cdot\r_n}\right|\phi_f\right>\right|
\left|\left<\phi_f\left|e^{-i\K\cdot\r_m}\right|\phi_0\right>\right|\,.
\eeqa
then as 
\beq
e^{-iHt/\hbar}\left|\phi_l\right>=e^{-iE_lt/\hbar}\left|\phi_l\right>
\eeq
where $H$ is the lattice Hamiltonian, we can write this as
\beqa
S(\K,\Omega)&=&\frac{1}{N}\sum_f
\int_{-\infty}^\infty \frac{dt}{2\pi} e^{i\Omega t}
\sum_{\r_n,\r_m}
\left<\phi_0\left|e^{iHt/\hbar}e^{i\K\cdot\r_n}e^{-iHt/\hbar}\right|\phi_f
\right>\nonumber\\
&&\left<\phi_f\left|e^{-i\K\cdot\r_m}\right|\phi_0\right>\,,
\eeqa
and the argument in the first expectation value is the time-dependent 
operator $e^{i\K\cdot\r_n}$ in the Heisenberg representation
\beq
e^{i\K\cdot\r_n(t)}=e^{iHt/\hbar}e^{i\K\cdot\r_n}e^{-iHt/\hbar}\,.
\eeq
Thus,
\beqa
S(\K,\Omega)&=&\frac{1}{N}\sum_f
\int_{-\infty}^\infty \frac{dt}{2\pi} e^{i\Omega t}
\sum_{\r_n,\r_m}
\left<\phi_0\left|e^{i\K\cdot\r_n(t)}\right|\phi_f\right>
\left<\phi_f\left|e^{-i\K\cdot\r_m}\right|\phi_0\right>\nonumber\\
&=&\frac{1}{N}
\int_{-\infty}^\infty \frac{dt}{2\pi} e^{i\Omega t}
\sum_{\r_n,\r_m}
\left<\phi_0\left|e^{i\K\cdot\r_n(t)}e^{-i\K\cdot\r_m}\right|\phi_0\right>\,.
\eeqa
Now since $\r_n(t)=\r_n+\s_n(t)$ (with $\r_n$ time independent),
\beq
S(\K,\Omega)=\frac{1}{N} \sum_{n,m}
\int_{-\infty}^\infty \frac{dt}{2\pi} e^{i\lep \K\cdot(\r_n-\r_m)+\Omega t\rip}
\left<\phi_0\left|e^{i\K\cdot\s_n(t)}e^{-i\K\cdot\s_m}\right|\phi_0\right>\,.
\eeq

This formula is correct at zero temperature.  In order to describe finite 
$T$ effects (ie., anti-stokes processes involving phonon absorption) we 
must introduce a thermal average over all states
\beq
\left<\phi_0\left|A\right|\phi_0\right> \to
\left<A\right> =
\sum_le^{-\beta E_l}
\left<\phi_l\left|A\right|\phi_l\right>/\sum_le^{-\beta E_l}\,.
\eeq
With this substitution,
\beq
S(\K,\Omega)=\frac{1}{N} \sum_{n,m}
\int_{-\infty}^\infty \frac{dt}{2\pi} e^{i\lep \K\cdot(\r_n-\r_m)+\Omega t\rip}
\left< e^{i\K\cdot\s_n(t)}e^{-i\K\cdot\s_m}\right>\,.
\eeq
and $S(\K,\Omega)$ is called the dynamical structure factor.

\subsubsection{The Debye-Waller Factor}

     To simplify this relation further, recall that the
exponentiated operators within the brackets are linear functions
of the creation and annihilation operators $a^{\dag}$  and $a$.
\beq
\s_{n,\alpha}(t)=\frac{1}{\sqrt{M_\alpha N}}
\sum_{\q,s}\sqrt{\frac{\hbar}{2\om_s(\q)}} {\bf{\ep}}^s_{\alpha}(\q)
\lep a_s(\q)(t)+ a_s^{\dag}(-\q)(t)\rip e^{i\q\cdot\r_n}
\eeq
So that, in particular 
$\left<\s_{n,\alpha,i}(t)\right>= \left<s_{n,\alpha,i}(0)\right>=0$.  
Then let $A=i\K\cdot \s_{n,\alpha,i}(t)$  and $B=i\K\cdot \s_{m,\alpha,i}(0)$  
and suppose that the expectation values
of $A$ and $B$ are small.  Then
\beqa
\left< e^A e^B \right > 
&=& \left<(1+ A +\frac12 A^2+\cdots) 
(1+ B +\frac12 B^2+\cdots)\right>\nonumber\\
&\approx& \left<1+ A + B + AB +\frac12 A^2+\frac12 B^2+\cdots \right >\nonumber\\
&\approx& 1+ \frac12 \left<2AB + A^2+ B^2 \right >+\cdots\nonumber\\
&\approx& e^{\frac12 \left<2AB + A^2+ B^2 \right >}
\eeqa
This relation is in fact true to all orders, as long as $A$ and $B$
are linear functions of $a^{\dag}$  and $a$  . (c.f.\  \AM{792}, Callaway pp.\
41-48).
Thus
\beq
\left< e^{i\K\cdot\s_n(t)}e^{-i\K\cdot\s_m}\right>=
e^{-\frac12 \left< \lep\K\cdot\s_n(t)\rip^2\right>}
e^{-\frac12 \left< \lep\K\cdot\s_m\rip^2\right>}
e^{\left<  \K\cdot\s_n(t)\K\cdot\s_m\right>}\,.
\eeq
Since the Hamiltonian has no time dependence, and the lattice is
invariant under translations $r_n$
\beq
\left< e^{i\K\cdot\s_n(t)}e^{-i\K\cdot\s_m}\right>=
e^{-\left< \lep\K\cdot\s_n\rip^2\right>}
e^{\left<  \K\cdot\s_{n-m}(t)\K\cdot\s_0\right>}\,,
\eeq
where the first term is called the Debye-Waller factor $e^{-2W}$.
\beq
e^{-2W}=e^{-\left<\lep \K\cdot\s_n\rip^2\right>}\,.
\eeq
Thus letting $l=n-m$
\beq
S(\K,\Omega)= e^{-2W}\sum_{l}
\int_{-\infty}^\infty \frac{dt}{2\pi} e^{i\lep \K\cdot\r_l+\Omega t\rip}
e^{\left<  \K\cdot\s_{l}(t)\K\cdot\s_0\right>}\,.
\eeq
Here the Debye-Waller factor contains much of the crucial quantum
physics.  It is finite, even at $T=0$ due to zero-point
fluctuations, and since $\left< \K\cdot\s_n\right>^2$  will increase with 
temperature, the total strength
of the Bragg peaks will diminish with increasing $T$.  However, as
long as a crystal has long-ranged order, it will remain finite.

\subsubsection{Zero-phonon Elastic Scattering}

One may disentangle the elastic and inelastic processes by expanding the 
exponential in the equation above.
\beq
e^{\left<  \K\cdot\s_{l}(t)\K\cdot\s_0\right>}
=\sum_m\frac{1}{m!} \lep \left<  \K\cdot\s_{l}(t)\K\cdot\s_0\right>\rip^m
\eeq

If we approximate the exponential by 1, ie.\ take only the first, $m=0$ 
term, then
\beq
S_0(\K,\Omega)= e^{-2W}\sum_{l}
\int_{-\infty}^\infty \frac{dt}{2\pi} e^{i\lep \K\cdot\r_l-\Omega t\rip}
\,.
\eeq
And we recover the lowest order classical result (modified by the
Debye-Waller factor) which gives us the Bragg conditions that 
$S_0(\K,\Omega)$ is only finite when $\K=\G$ and  $\Omega=\om_0-\om_f=0$.
\beq
S_0(\K,\Omega)= e^{-2W} \delta(\Omega)N\sum_\G \delta_{\K,\G}
\,,
\eeq
\beq
\frac{d\sigma_0}{dE d\Omega}=\frac{k}{k_0}\frac{Na^2}{\hbar} 
e^{-2W} \delta(\Omega)N\sum_\G \delta_{\K,\G}
\eeq

However, now the scattering intensity is reduced by the Debye-Waller
factor $e^{-2W}$, which accounts for zero-point motion and thermal 
fluctuations.

\subsubsection{One-Phonon Inelastic Scattering}

     When $m=1$, then the scattering involves either the absorption
or creation of a phonon.  To evaluate 
\beq
S_1(\K,\Omega)= e^{-2W}\sum_{l}
\int_{-\infty}^\infty \frac{dt}{2\pi} e^{i\lep \K\cdot\r_l+\Omega t\rip}
\left<  \K\cdot\s_{l}(t)\K\cdot\s_0(0)\right>\,.
\eeq
we need
\beq
\s_{n,\alpha}(t)=\frac{1}{\sqrt{M_\alpha N}}
\sum_{\q,s}\sqrt{\frac{\hbar}{2\om_s(\q)}} {\bf{\ep}}^s_{\alpha}(\q)
\lep a_s(\q,t)+ a_s^{\dag}(-\q,t)\rip e^{i\q\cdot\r_n}
\eeq
in the Heisenberg representation, and therefore we need,
\beqa
a(\q,t) &=& e^{iHt/\hbar} a(\q) e^{-iHt/\hbar} \nonumber \\
	&=& e^{i(\om(\q) t a^{\dag}(\q) a(\q))} a(\q) 
	    e^{-i(\om(\q) t a^{\dag}(\q) a(\q)} \nonumber\\
	&=& a(\q) e^{i\lep\om(\q) t (a^{\dag}(\q) a(\q)-1)\rip}
	    e^{-i(\om(\q) t a^{\dag}(\q) a(\q)} \nonumber\\
	&=& a(\q) e^{-i\om(\q) t}
\eeqa
where we have used the fact that $(a^{\dag}a)^na=(a^{\dag}a)^{n-1}a^{\dag}aa
=(a^{\dag}a)^{n-1}(aa^{\dag}-1)a=(a^{\dag}a)^{n-1}a(a^{\dag}a-1)=
a(a^{\dag}a-1)^n$.
Similarly $a^{\dag}(\q,t)=a^{\dag}(\q) e^{i\om(\q) t}$.  Thus,
\beq
\s_{n,\alpha}(t)=\frac{1}{\sqrt{M_\alpha N}}\sum_{\q,s}
\frac{\sqrt{\hbar}e^{i\q\cdot\r_n}}{\sqrt{2\om_s(\q)}} {\bf{\ep}}^s_{\alpha}(\q)
\lep a_s(\q)e^{-i\om_s(\q)t}+ a_s^{\dag}(-\q)e^{i\om_s(\q)t}\rip 
\eeq
and
\beq
\s_{0,\alpha}(0)=\frac{1}{\sqrt{M_\alpha N}}
\sum_{\p,r}\sqrt{\frac{\hbar}{2\om_r(\p)}} {\bf{\ep}}^r_{\alpha}(\p)
\lep a_r(\p)+ a_r^{\dag}(-\p)\rip 
\eeq

Recall, we want to evaluate
\beq
S_1(\K,\Omega)= e^{-2W}\sum_{l}
\int_{-\infty}^\infty \frac{dt}{2\pi} e^{i\lep \K\cdot\r_l+\Omega t\rip}
\left<  \K\cdot\s_{l}(t)\K\cdot\s_0(0)\right>\,.
\eeq
Clearly, the only terms which survive in 
$\left<  \K\cdot\s_{l}(t)\K\cdot\s_0(0)\right>$ are those with
$r=s$ and $\p=-\q$.  Furthermore, the sum over $l$ yields a delta
function $N\sum_{\G}\delta_{\K+\q,\G}$.  Then as $\ep(\G-\k)=\ep(-\k)
=\ep^*(\k)$, and $\om(\G-\k)=\om(\k)$,
\beqa
S_1(\K,\Omega)&=& e^{-2W}\int_{-\infty}^\infty \frac{dt}{2\pi} 
e^{i \Omega t}\sum_{\q,\G,s} 
\frac{\hbar\left|\K\cdot\ep(\K)\right|^2}{2\om_s(\q)M}
 \delta_{\K+\q,\G}\\
&&\leb 
e^{-i\om_s(\q)t}\left< a_s(-\K)a_s^{\dag}(-\K)\right> +
e^{i\om_s(\q)t}\left< a_s^{\dag}(-\K)a_s(-\K)\right>
\rib\nonumber
\eeqa
The occupancy of each mode  $n(\q)$ is given by the Bose factor
\beq
\left< n(\q) \right> =\frac{1}{e^{\beta\om(\q)}-1}
\eeq
So, finally
\beqa
S_1(\K,\Omega)&=& 
e^{-2W} \sum_s \frac{\hbar}{2M\om_s(\K)} 
\left|\K\cdot{\bf{\ep}}^s(\K)\right|^2\\
&&\leb(1+n_s(\K)) \delta(-\Omega+\om_s(\K)) +
n_s(\K)\delta(\Omega+\om_s(\K))\rib
\,.\nonumber
\eeqa
For the first term, we get a contribution only when
$\Om-\om_s(\K)=\om_0-\om_f-\om_s(\K)=0$; ie., the final energy of the 
neutron is smaller than the initial energy.  The energy is lost in the 
creation of a phonon.  Note that this can happen at any temperature,
since $(1+n_s(\K))\neq 0$ at any $T$.  The second term is only finite when
$\Om+\om_s(\K)=\om_0-\om_f+\om_s(\K)=0$; ie., the final energy of the
neutron is larger than the initial energy.  The additional energy comes
from the absorption of a phonon.  Thus phonon absorption is only allowed 
at finite temperatures, and in fact, the factor  $n_s(\K)= 0$ at zero
temperature.  These terms correspond to the Stokes and anti-Stokes
processes, respectively, illustrated in Fig.~\ref{fig:stokes2}.
\begin{figure}[htb]
\centerline{\includegraphics[width=0.8\textwidth,clip=true]{stokes2.pdf}}
\caption[]{\em{Stokes and antistokes processes in inelastic neutron 
scattering involving the creation or absorption of a lattice phonon.
The antistokes process can only occur at finite-T, when $n_s(\K)\neq 0$.}}
\label{fig:stokes2}
\end{figure}

	If we were to continue our expansion of the exponential to
larger values of $m$, we would find multiple-phonon scattering processes.
However, these terms are usually of minimal contribution to the total
cross section, due to the fact that the average ionic excursion $s$
is small, and are usually neglected. 

\edo
