\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 5: Thermal Properties of Crystal Lattices}
\author{Debye}
\input{../defs}
\def\baselinestretch{1.4}

\bdo
\maketitle
%\large
\tableofcontents
\pagebreak
%\Large

In the previous chapter, we have shown that the motion of a harmonic crystal 
can be described by a set of decoupled harmonic oscillators.
\beq
H=\frac12\sum_{\k,s} \left|\P_s(\k)\right|^2
+ \om_s^2(\k) \left|\Q_s(\k)\right|^2
\eeq
At a given temperature T, the occupancy of a given mode is
\beq
\left<n_s(\k)\right> =\frac{1}{e^{\beta \om_s(\k)}-1}
\eeq
In this chapter, we will apply this information to calculate the
thermodynamic properties of the ionic lattice, in addition to addressing
questions regarding its long-range order in the presence of lattice
vibrations (i.e. do phonons destroy the order).  In order to evaluate
the different formulas for these quantities, we will first discuss two
matters of formal convenience.

\section{Formalism}

To evaluate some of these properties we can use the virial
theorem and, integrals over the density of states.  

\subsection{The Virial Theorem}

Consider the Hamiltonian for a quantum system $H(\x,\p)$, where $\x$ and
$\p$ are the canonically conjugate variables.  Then the expectation value 
of any function of these canonically conjugate variables $f(\x,\p)$ in a 
stationary state (eigenstate) is constant in time.  Consider
\beqa
\frac{d}{dt}\left<\x\cdot\p\right>&=&
\frac{i}{\hbar}\left<\leb H,\x\cdot\p\rib\right>=
\frac{i}{\hbar}\left< H \x\cdot\p - \x\cdot\p H \right>\nonumber\\
&=&
\frac{iE}{\hbar}\left< \x\cdot\p - \x\cdot\p \right>=0
\eeqa
where $E$ is the eigenenergy of the stationary state.  Let
\beq
H=\frac{p^2}{2m} + V(\x)\,,
\eeq
then
\beqa
0&=& \left<\leb\frac{p^2}{2m} + V(\x),\x\cdot\p\rib\right>\nonumber\\
 &=& \left<
\leb\frac{p^2}{2m},\x\cdot\p\rib +
\leb V(\x),\x\cdot\p\rib
\right>\nonumber\\
 &=& \left<
\frac{1}{2m}\leb p^2,\x\rib\cdot\p +
\x\cdot\leb V(\x),\p\rib
\right>
\eeqa
Then as $\leb p^2,\x\rib = \p\leb p,\x\rib + \leb p,\x\rib\p = -2i\hbar\p$
and $\p=-i\hbar\grad_x$,
\beq 0=
\left< 
\frac{-i\hbar}{m} p^2 +
i\hbar\x\cdot \grad_xV(\x)
\right>
\eeq
or, the Virial theorem:
\beq
2\left<T\right>=\left<\x\cdot \grad_xV(\x)\right>
\eeq
Now, lets apply this to a harmonic oscillator where $V=\frac12 m\om^2x^2$
and $T=p^2/2m$, $<H>=<T>+<V>=\hbar\om(n+\frac12)$ and .  We get
\beq
\left<T\right>=\left<V(\x)\right>=\frac12 \hbar\om(n+\frac12)
\eeq

Lets now apply this to find the RMS excursion of a lattice site in
an elemental lattice $r=1$
\beqa
\left<s^2\right> 
&=& \frac{1}{N} \sum_{n,i}\left<s^2_{n,i}\right>\nonumber\\
&=& \frac{1}{N} \left<\sum_{n,i}\frac{1}{NM} 
\sum_{\q,s,\k,r} Q_s(\q) \ep^s_i(\q) e^{i\q\cdot\r_n}
Q_r(\k) \ep^r_i(\k) e^{i\k\cdot\r_n} \right>\nonumber\\
&=& \frac{1}{N} \left<\sum_{n,i}\frac{1}{NM} 
\sum_{\q,s,\k,r} Q_s(\q) \ep^s_i(\q) e^{i\q\cdot\r_n}
Q_r(-\k) \ep^r_i(-\k) e^{-i\k\cdot\r_n} \right>\nonumber\\
&=& \frac{1}{N} \left<\sum_{n,i}\frac{1}{NM} 
\sum_{\q,s,\k,r} Q_s(\q) \ep^s_i(\q) e^{i\q\cdot\r_n}
Q_r^*(\k) \ep^{*r}_i(\k) e^{-i\k\cdot\r_n} \right>\nonumber\\
&=& \frac{1}{NM}\sum_{\q,s}\left< \left|Q_s(\q)\right|^2\right>\nonumber\\
&=& \frac{1}{NM}\sum_{\q,s}\frac{2}{\om_s^2(\q)}
\left<\frac12 \om_s^2(\q) \left|Q_s(\q)\right|^2\right>\nonumber\\
&=& \frac{1}{NM}\sum_{\q,s}\frac{2}{\om_s^2(\q)}
\frac12 \hbar \om_s(\q)\lep n_s(\q)+\frac12\rip\nonumber\\
\left<s^2\right> &=& \frac{1}{NM}\sum_{\q,s}\frac{\hbar}{\om_s(\q)}
\lep n_s(\q)+\frac12\rip
\eeqa

This integral, which must be finite in order for the system to
have long-range order, is still difficult to perform.  However,
the integral may be written as a function of $\om_s(\q)$ only.
\beq
\left<s^2\right> = \frac{1}{NM}\sum_{\q,s}\frac{\hbar}{\om_s(\q)}
\lep \frac{1}{e^{\beta\om_s(\q)}-1}+\frac12\rip
\eeq
It would be convenient therefore, to introduce a density of phonon states
\beq
Z(\om) =\frac{1}{N}\sum_{\q,s}\delta\lep\om-\om_s(\q)\rip
\eeq
so that
\beq
\left<s^2\right> = \frac{\hbar}{M}\int d\om Z(\om)
\frac{1}{\om} \lep n(\om) +\frac12\rip
\eeq

\subsection{The Phonon Density of States}

     In addition to the calculation of $\left<s^2\right>$, the density of 
states $Z(\om)$ is also useful in the calculation of $E=<H>$, the partition
function, and the related thermodynamic properties

\begin{wrapfigure}{l}{0.5\textwidth}
\includegraphics[width=0.49\textwidth,clip=true]{firstbz.pdf}
\protect\caption{\emph{First Brillouin zone of the square lattice}\label{fig:firstbz}}
\end{wrapfigure}
In order to calculate
\beq
Z(\om) =\frac{1}{N}\sum_{\q,s}\delta\lep\om-\om_s(\q)\rip
\eeq
we must first better define the sum over $\q$.  As we discussed last
chapter for a 3-d cubic system, we will assume that we have a periodic
finite lattice of $N$ basis points, or $N^{1/3}$ in each of the principle
lattice directions $\aone$, $\atwo$, and $\athree$.  Then, we want the 
Fourier representation to respect the periodic boundary conditions (pbc), so
\beq
e^{i\q\cdot\lep\r+N^{1/3}(\aone+\atwo+\athree)\rip}=e^{i\q\cdot\r}
\eeq
This means that $q_i=2\pi m/N^{1/3}$, where $m$ is an integer, and $i$
indicates one of the principle lattice directions.  In addition, we only 
want unique values of $q_i$, so we will choose those within the first 
Brillouin zone, so that
\beq
\G\cdot\q \leq \frac12 G^2
\eeq
The size of this region is the same as that of a unit cell of the
reciprocal lattice $\gone\cdot\lep\gtwo\times\gthree\rip$.  Since 
there are $N$ states in this region, the density of $\q$ states is
\beq
\frac{N}{\gone\cdot\lep\gtwo\times\gthree\rip} =
\frac{NV_c}{(2\pi)^3}=\frac{V}{(2\pi)^3}
\eeq
where $V_c$ is the volume of a Bravais lattice cell 
($\aone\cdot\atwo\times\athree$), and $V$ is the lattice volume.

\begin{wrapfigure}{r}{0.5\textwidth}
\includegraphics[width=0.49\textwidth,clip=true]{dosgrid.pdf}
\protect\caption{\emph{States in q-space.  $S_\om$ is the surface of constant 
$\om=\om_s(\q)$, so that $d^3q=dS_wdq_\perp=
\frac{dS_\om d\om}{\grad_q\om_s(\q)}$.}\label{fig:dosgrid}}
\end{wrapfigure}
Clearly as $N\to\infty$ the density increases until a continuum of states
is formed (all that we need here is that the spacing between
q-states be much smaller than any physically relevant value of
$\q$). The number of states in a frequency interval $d\om$ is then
given by the volume of q-space between the surfaces defined by
$\om=\om_s(\q)$ and $\om=\om_s(\q)+d\om$ multiplied by $V/(2\pi)^3$
\beqa
Z(\om)d\om &=& \frac{V}{(2\pi)^3}\int_\om^{\om+d\om} d^3q\\
&=& \frac{V}{(2\pi)^3} d\om \sum_s \int d^3q \delta\lep\om-\om_s(\q)\rip\nonumber
\eeqa
As shown in Fig.~\ref{fig:dosgrid}, $d\om=\grad_q\om_s(\q)dq_\perp$, and
\beq
d^3q=dS_wdq_\perp=\frac{dS_\om d\om}{\grad_q\om_s(\q)}\,,
\eeq
where $S_\om$ is the surface in q-space of constant $\om=\om_s(\q)$.
Then
\beq
Z(\om)d\om = \frac{V}{(2\pi)^3} d\om \sum_s \int_{\om=\om_s(\q)} \frac{dS_\om}{\grad_q\om_s(\q)}
\eeq
Thus the density of states is high in regions where the
dispersion is flat so that $\grad_q\om_s(\q)$ is small.
\begin{figure}
\begin{center}
\includegraphics[width=0.85\textwidth,clip=true]{1d.pdf}
\protect\caption{\emph{Linear harmonic chain.  The dispersion $\om(q)$ is linear
near $q=0$, and flat near $q=\pm \pi/a$ .  Thus, the density of states (DOS) is 
flat near $\om=0$ corresponding to the acoustic mode (for which
$\grad_q\om(q)=$constant), and divergent near $\om=\om_0$ corresponding to the 
peak of the dispersion (where $\grad_q\om(q)=0$) }\label{fig:1d}
}
\end{center}
\end{figure}

As an example, consider the 1-d Harmonic chain shown in Fig.~\ref{fig:1d}.
The phonon dispersion must be symmetric around $q=0$ and $q=pi/a$ due to the 
point group symmetries.  In addition, there must be an acoustic mode as
shown in Chapter 4.  As a result, the dispersion is linear near $q=0$ and 
has a peak at the one boundary.  Since the peak is flat there is a corresponding 
peak in the phonon DOS.  In fact, any point within the Brillouin zone for which
$\grad_q\om(q)=0$   (cusp, maxima, minima) will yield an integrable singularity 
in the DOS.

\section{Models of Lattice Dispersion}
\subsection{The Debye Model}
     For most thermodynamic properties, we are interested in the
modes $\hbar\om(\q)\sim k_B T$ which are low frequency modes in general.  
\begin{figure}
\begin{center}
\includegraphics[width=0.85\textwidth,clip=true]{dispersion.pdf}
\protect\caption{\emph{Dispersion for the diatomic linear chain.  In the Debye model,
we replace the acoustic mode by a purely linear mode and ignore any optical modes.}\label{fig:debye}
}
\end{center}
\end{figure}
From a very general set of (symmetry) constraints we have
argued that all interacting lattices in which the total energy is
invariant to an overall arbitrary rigid shift in the location of the
lattice must have at least one acoustic
mode, where for small $\om_s(\q)=c_s |\q|$.  Thus, for the thermodynamic
properties of the lattice, we care predominantly about the 
limit $\om(\q)\to 0$.  This physics is rather accurately described by the 
Debye model.


In the Debye model, we will assume that all modes are acoustic (elastic), so
that $\om_s(\q)=c_s |\q|$ for all $s$ and $\q$, then $\grad_q\om_s(\q)=c_s$
for all $s$ and $\q$, and
\beqa
Z(\om)
&=&\frac{V}{(2\pi)^3}\sum_s \int \frac{dS_\om}{\grad_q\om_s(\q)}\nonumber\\
&=&\frac{V}{(2\pi)^3}\sum_s\int \frac{dS_\om}{c_s}
\eeqa
The surface integral may be evaluated, and yields a constant 
$\int dS_\om= S_s$ for each branch.  Typically $c_s$ is different for 
different modes.  However, we will assume that the system is isotropic, so
$c_s=c$. If the dispersion is isotopic, 
then the surface of constant $\om_s(\q)$ is just a sphere, so the surface 
integral is trivial
\beq
S_s(\om=\om(\q))=\left\{\begin{array}{c} 2\;\;\mbox{for $d=1$}\\
2\pi q= 2\pi\om/c\;\;\mbox{for $d=2$}\\
4\pi q^2=4\pi \om^2/c^2 \;\;\mbox{for $d=3$}
\end{array}\right.
\eeq
then since the number of modes $=d$
\beq
Z(\om)=\frac{V}{(2\pi)^3}\left\{\begin{array}{c} 2/c\;\;\mbox{for $d=1$}\\
2\pi q= 4\pi\om/c^2\;\;\mbox{for $d=2$}\\
4\pi q^2=12\pi \om^2/c^3 \;\;\mbox{for $d=3$}
\end{array}\right\}\;\;\;0<\om<\om_D
\label{eq:debyedos}
\eeq
Note that since the total number of states is finite, we have
introduced a cutoff $\om_D$ on the frequency.

\subsection{The Einstein Model}

     ``Real'' two dimensional systems, i.e., a monolayer of gas (He) deposited 
on an atomically perfect surface (Vycor), may be better described by an 
Einstein model where each atom oscillates with a frequency $\om_0$ and does 
not interact with its neighbors.  The model is dispersionless $\om(\q)=\om_0$, 
and the DOS for this system is a delta function $Z(\om)=c\delta(\om-\om_0)$.
Note that it does not have an acoustic mode; however, this is not in 
violation of the discussion in the last chapter.  Why?
\begin{figure}
\begin{center}
\includegraphics[width=0.75\textwidth,clip=true]{einstein.pdf}
\protect\caption{\emph{Helium on a Vycor surface.  Each He is attracted
weakly to the surface by a van der Waals attraction and sits in a local
minimum of the surface lattice potential.}\label{fig:einstein}
}
\end{center}
\end{figure}

\section{Thermodynamics of Crystal Lattices}

We are now in a situation to calculate many of the thermodynamic
properties of crystal lattices.  However before addressing such
questions as the lattice energy free energy and specific heat we
should see if our model has long-range order... ie., is it
consistent with our initial assumptions.

\subsection{Mermin-Wagner Theory, Long-Range Order}

For simplicity, we will work on an elemental lattice model.
We may define long-ranged order (LRO) as a finite value of
\beqa
\langle s^2 \rangle &=&\frac{\hbar}{MN}\sum_{\q,s}\frac{1}{\om_s(\q)}
\lep \frac{1}{e^{\beta\om_s(\q)}-1} + \frac12\rip\nonumber\\
&=&
\frac{\hbar}{2MN} \sum_{\q,s}\frac{{\rm{sinh}}\lep\beta\om_s(\q)/2\rip}
{\om_s(\q){\rm{cosh}}\lep\beta\om_s(\q)/2\rip}
\eeqa
Since we expect all lattices to melt for some high temperature, we are
interested only in the $T\to 0$ limit.  Given the factor of 
$\frac{1}{\om_s(\q)}$ in the summand, we are most interested in acoustic 
modes since they are the ones which will cause a divergence.
\beq
\lim_{\beta\to\infty}\langle s^2 \rangle =
\frac{\hbar}{2MN} \lim_{\beta\to\infty} 
\sum_{\q,s}\frac{1}{\om_s(\q)}\lep \frac{1}{\beta\om_s(\q)} + \frac12\rip
\eeq
The low frequency modes are most important so a Debye
model may be used
\beq
\lim_{\beta\to\infty}\langle s^2 \rangle \approx
\frac{\hbar}{2MN}
\int_0^{\om_D} d\om Z(\om) 
\frac{1}{\om}\lep \frac{1}{\beta\om} + \frac12\rip\,,
\eeq
where $Z(\om)$ is the same as was defined above in Eq.~\ref{eq:debyedos}.
\beq
Z(\om)=\frac{V}{(2\pi)^3}\left\{\begin{array}{c} 2/c\;\;\mbox{for $d=1$}\\
2\pi q= 4\pi\om/c^2\;\;\mbox{for $d=2$}\\
4\pi q^2=12\pi \om^2/c^3 \;\;\mbox{for $d=3$}
\end{array}\right\}\;\;\;0<\om<\om_D
\eeq

Thus in 3-d the DOS always cancels the $1/\om$ singularity but in two
dimensions the singularity is only cancelled when $T=0$ ($\beta=\infty$), 
and in one dimension $\langle s^2 \rangle=\infty$  for all $T$.
\begin{table}[htb]
\begin{center}
\begin{tabular}{|l|l|l|l|}\hline
$\langle s^2 \rangle$	& $d=1$		& $d=2$		& $d=3$ \\ \hline
$T=0$	&$\infty$	& finite	& finite \\
$T\neq 0$&$\infty$	& $\infty$	& finite \\\hline
\end{tabular}
\end{center}
\caption[]{\em{$\langle s^2 \rangle$ for lattices of different dimension,
assuming the presence of an acoustic mode.}}
\label{lattice-MW}
\end{table}
This is a specific case of the Mermin-Wagner Theorem.  We should emphasize 
that the result $\langle s^2 \rangle=\infty$  does not mean that our
theory has failed.  The harmonic approximation requires that the
near-neighbor strains must be small, not the displacements.

Physically, it is easy to understand why one-dimensional systems do
not have long range order, since as you go along the chain, the displacements
of the atoms can accumulate to produce a very large rms displacement.
In higher dimensional systems, the displacements in any direction are 
constrained by the neighbors in orthogonal directions.  
\begin{figure}
\begin{center}
\includegraphics[width=0.75\textwidth,clip=true]{1ddis.pdf}
\protect\caption{\emph{Random fluctuaions of atoms in a 1-d lattice may accumulate
to produce a very large displacements of the atoms.}\label{fig:1ddis}
}
\end{center}
\end{figure}
``Real'' two dimensional systems, i.e., a monolayer of gas deposited on an atomically 
perfect surface, do have long-range order even at finite temperatures due to 
the surface potential (corrugation of the surface).  These may be better 
described by an Einstein model where each atom oscillates with a frequency 
$\om_0$ and does not interact with its neighbors.  The DOS for this system 
is a delta function as described above.  For such a DOS, 
$\langle s^2 \rangle$ is always finite.  You will explore this physics, in 
much more detail, in your homework.


\subsection{Thermodynamics}

     We will assume that our system is in equilibrium with a heat
bath at temperature $T$.  This system is described by the canonical
ensemble, and may be justified by dividing an infinite system into
a finite number of smaller subsystems.  Each subsystem is expected to
interact weakly with the remaining system which also acts as the subsystems
heat bath.  The probability that any state in the subsystem is occupied is
given by
\beq
P\lep\left\{n_s(\k)\right\}\rip \propto 
e^{-\beta E\lep\left\{n_s(\k)\right\}\rip}
\eeq
Thus the partition function is given by
\beqa
{\cal{Z}} &=& \sum_{\left\{n_s(\k)\right\} }
e^{-\beta E\lep\left\{n_s(\k)\right\}\rip}\nonumber\\
&=& \sum_{\left\{n_s(\k)\right\}}
e^{-\beta \sum_{\k,s}\hbar\om_s(\k)\lep n_s(\k) + \frac12\rip}\nonumber\\
&=& \prod_{s,\k} {\cal{Z}}_s(\k)
\eeqa
where ${\cal{Z}}_s(\k)$ is the partition function for the mode
$s,\k$; i.e.\ the modes are independent and decouple.
\beqa
{\cal{Z}}_s(\k)
&=&\sum_n e^{-\beta \hbar\om_s(\k)\lep n_s(\k) + \frac12\rip }\nonumber\\
&=& e^{-\beta \hbar\om_s(\k)/2 }
\sum_n e^{-\beta \hbar\om_s(\k)\lep n_s(\k)\rip }\nonumber\\
&=& \frac{ e^{-\beta \hbar\om_s(\k)/2 }}{1-e^{-\beta \hbar\om_s(\k) }}\nonumber\\
&=& \frac{1}{2\sinh\lep\beta\om_s(\k)/2\rip}
\eeqa

The free energy is given by
\beq
{\cal{F}} = -k_BT\ln\lep{\cal{Z}}\rip=
k_BT\sum_{\k,s}\ln\lep2\sinh\lep\beta\hbar\om_s(\k)/2\rip\rip
\eeq
Since $d\CE=Td\CS-PdV$ and $d\CF=Td\CS -PdV -Td\CS-\CS dT$, the entropy is
\beq
{\cal{S}}= -\lep\pde{{\cal{F}}}{T}\rip_V\,,
\eeq
and system energy is then given by
\beq
{\cal{E}}={\cal{F}}+T{\cal{S}} ={\cal{F}}-T\lep\pde{{\cal{F}}}{T}\rip_V
\eeq
where constant volume $V$ is guaranteed by the harmonic approximation
(since $<s>=0$).
\beq
{\cal{E}}=\sum_{\k,s}\frac12 \hbar\om_s(\k)
{\rm{coth}}\lep\beta\om_s(\k)/2\rip
\eeq
The specific heat is then given by
\beq
C=\lep\der{{\cal{E}}}{T}\rip_V=
k_B\sum_{\k,s} \lep\beta\hbar\om_s(\k)\rip^2
{\rm{csch}}^2\lep\beta\om_s(\k)/2\rip
\eeq
where ${\rm{csch}}\lep x\rip=1/{\rm{sinh}}\lep x\rip$

Consider the specific heat of our 3-dimensional Debye model.
\beqa
C&=&k_B\int_0^{\om_D} d\om {\cal{Z}}(\om) \lep\beta\om/2\rip^2
{\rm{csch}}^2\lep\beta\om/2\rip\nonumber\\
&=& k_B\int_0^{\om_D} d\om \lep\frac{12V\pi\om^2}{(2\pi c)^3}\rip \lep\beta\om/2\rip^2
{\rm{csch}}^2\lep\beta\om/2\rip
\eeqa
Where the Debye frequency $\om_D$ is determined by the requirement that
\beq
3rN=\int_0^{\om_D} d\om {\cal{Z}}(\om)= 
\int_0^{\om_D} d\om \lep\frac{12V\pi\om^2}{(2\pi c)^3}\rip\,,
\eeq
or $V/(2\pi c)^3=3rN/(4\om_D^3)$.

Clearly the integral for C is a mess, except in the high and low
T limits.  At high temperatures $\beta\hbar\om_D/2\ll 1$,
\beq
C\approx k_B\int_0^{\om_D} d\om {\cal{Z}}(\om) = 3Nrk_B
\eeq
This is the well known classical result (equipartition theorem) which 
attributes $(1/2) k_B$ of the specific heat to each quadratic degree of 
freedom.  Here for each element of the basis we have 6 quadratic degrees 
of freedom (three translational, and three momenta).

At low temperatures, $\beta\hbar\om/2\gg 1$,
\beq
{\rm{csch}}^2\lep\beta\om/2\rip \approx 2e^{-\beta\om/2}
\eeq
Thus, at low $T$, only the low frequency modes contribute, so
the upper bound of integration may be extended to $\infty$
\beq
C\approx 12\pi k_B \frac{3rN}{4\om_D^2} \int_0^\infty d\om \om^2 
\lep\frac{\beta\hbar\om}{2}\rip^2 2e^{-\beta\om_s(\k)/2}\,.
\eeq
If we make the change of variables $x=\beta\hbar\om/2$, we get
\beqa
C&\approx& \frac{9k_BrN\pi}{2}\lep\frac1{\om_D\beta\hbar}\rip^3
\int_0^\infty dx x^4 e^{-x}\\
 &\approx & \frac{9k_BrN\pi}{2}\lep\frac1{\om_D\beta\hbar}\rip^3 24
\eeqa
Then, if we identify the Debye temperature $\theta_D=\hbar\om_D/k_B$,
we get
\beq
C\approx 96 \pi rNk_B\lep\frac{T}{\theta_D}\rip^3
\eeq
{\em{$C\propto T^3$ at low temperature is the characteristic 
signature of low-energy phonon excitations}}.

\subsection{Thermal Expansion, the Gruneisen Parameter}

Consider a cubic system of linear dimension $L$.  If unconstrained, we expect
that the volume of this system will change with temperature (generally expand
with increasing $T$, but not always. cf.\ ice or Si).  We define the
coefficient of free expansion ($P=0$) as
\beq
\alpha_L=\frac{1}{L}\der{L}{T} \;\;\;\mbox{or}\;\;\;
\alpha_V=3\alpha_L=\frac{1}{V}\der{V}{T}\,.
\eeq
Of course, this measurement only makes sense in equilibrium.
\beq
P=-\lep\der{{\cal{F}}}{V}\rip_T=0
\eeq

	As mentioned earlier, since $<s>=0$ in the harmonic approximation, 
a harmonic crystal does not expand when heated.  Of course, real crystals
do, so that lack of thermal expansion of a harmonic crystal can be considered
a limitation of the harmonic theory.  To address this limitation, we
can make a quasiharmonic approximation.  Consider a more general potential
between the ions, of the form
\beq
V(x)=bx+cx^3+\frac12 m\om^2x^2
\eeq
and let's see if any of these terms will produce a temperature dependent
displacement.  The last term is the usual harmonic term, which we have
already shown does not produce a T-dependent $<x>$.  Also the first term 
does not have the desired effect!  It corresponds to a 
temperature-independent shift in the oscillator, as can be seen by completing
the square
\beq
\frac12 m\om^2x^2 + bx = \frac12\lep x + \frac{b}{m\om^2}\rip^2
-\frac{b^2}{2m\om^2}\;.
\eeq
Then $<x>=-\frac{b}{m\om^2}$, independent of the temperature; that is,
assuming that $b$ is temperature-independent.  What we need is a temperature
dependent coefficient $b$!
\begin{figure}
\begin{center}
\includegraphics[width=0.75\textwidth,clip=true]{cubic.pdf}
\protect\caption{\emph{Plot of the potential $V(x)=\frac12 mx^2+cx^3$ when $m=\omega=1$ 
and $c=0.0,\,0.1$.  The average position of a particle $<x>$ in the anharmonic 
potential, $c=0.1$, will shift to the left as the energy (temperature) is 
increased; whereas, that in the harmonic potential, $c=0$, is fixed $<x>=0$. }\label{fig:cubic}
}
\end{center}
\end{figure}

     The cubic term has the desired effect.  As can be seen in Fig~\ref{fig:cubic},
as the average energy (temperature) of a particle trapped in a cubic
potential increases, the mean position of the particle shifts.
However, it also destroys the solubility of the model.  To get around
this, approximate the cubic term with a mean-field decomposition.
\beq
cx^3\approx c\eta x\left<x^2\right>+c(1-\eta)x^2\left<x\right>
\eeq
and treat these two terms separately (the new parameter $\eta$ is to be 
determined self-consistently, usually by minimizing the free energy
with respect to $\eta$).  The first term yields the needed temperature 
dependent shift of $<x>$
\beq
\frac12 m\om^2x^2 + c\eta x\left<x^2\right>=
\frac12 m\om^2\lep x + \frac{c\eta<x^2>}{m\om^2}\rip^2 -
\frac{c^2<x^2>^2}{2m\om^2}
\eeq
so that $<x>=-\frac{c\eta<x^2>}{m\om^2}$.  Clearly the renormalization of 
the equilibrium position of the harmonic oscillator will be temperature 
dependent.  While the second term, $(1-\eta)x^2<x>$, yields a shift in 
the frequency $\om\to\om'$
\beq
\om'=\om\lep 1+\frac{2(1-\eta)<x>}{m\om^2}\rip^\frac12
\eeq
which is a function of the equilibrium position.
Thus a mean-field description of the cubic term is consistent with the 
observed physics.

     In what follows, we will approximate the effect of the anharmonic cubic 
term as a shift in the equilibrium position of the lattice (and hence the 
lattice potential) and a change of $\om$ to $\om'$; however, we imagine that 
the energy levels remain of the form
\beq
E_n=\hbar\om'(<x>)(n+\frac12),
\eeq
and that $<x>$ varies with temperature,
consistent with the mean-field approximation just described. 

     To proceed, imagine the cube of cubic system to be made up
of oscillators which are independent.  Since the final result
can be formulated as a sum over these independent modes, consider
only one.  In equilibrium, where $P=-\lep\der{ {\cal{F}} }{V}\rip_T=0$,
the free energy of one of the modes is
\beq
{\cal{F}}=\Phi + \frac12\hbar\om + k_BT\ln\lep 1-e^{-\beta\hbar\om}\rip
\eeq
and (following the notation of Ibach and L\"uth), let the lattice potential
\beq
\Phi=\Phi_0+\frac12f(a-a_0)^2 + \cdots
\eeq
where $f$ is the spring constant.  Then
\beq
0=P=\lep\der{ {\cal{F}} }{a}\rip_T=f(a-a_0)+ 
\frac{1}{\om}\pde{\om}{a}\lep \frac12 \hbar\om-\frac{\hbar\om}{1-e^{-\beta\hbar\om}}\rip,
\eeq
If we identify the last term in parenthesis as $\ep(\om,T)$, and
solve for $a$, then
\beq
a=a_0 -\frac{1}{\om f}\ep(\om,T)\pde{\om}{a}
\eeq

Since we now know $a(T)$ for a single mode, we may calculate the linear
expansion coefficient for this mode
\beq
\alpha_L=\frac{1}{a_0}\der{a}{T} = -\frac{1}{a^2_0 f}\pde{\ln w}{\ln a}
\pde{\ep(\om,T)}{T}
\eeq
To generalize this to a solid let $\alpha_L\to\alpha_V$ (as discussed above)
and $a_0^2f\to V^2\der{P}{V}=V\kappa$ ($\kappa$ is the bulk modulus) and 
sum over all modes the modes $\k,s$
\beq
\alpha_V=\frac{1}{V}\der{V}{T}=
\frac{1}{\kappa V}\sum_{\k,s} -\pde{\ln\om_s(\k)}{\ln V}
\pde{\ep\lep \om_s(\k),T\rip}{T}\,.
\eeq
Clearly (due to the factor of $\pde{\ep}{T}$), $\alpha_V$ will have 
a behavior similar to that of the specific heat ($\alpha_V\sim T^3$ for
low $T$, and $\alpha_V=$constant for high $T$). In addition, for many 
lattices, the Gruneisen number
\beq
\gamma=\pde{\ln\om_s(\k)}{\ln V}
\eeq
shows a weak dependence upon $s,\k$, and may be replaced by its
average, called the Gruneisen parameter
\beq
\left<\gamma\right>=\left<\pde{\ln\om_s(\k)}{\ln V}\right>\,,
\eeq
typically on the order of two.

Before proceeding to the next section, I would like to reexamine
the cubic term in a crystal where
\beqa
\sum_ls_l^3 &=& \frac{1}{(2MN)^{3/2}}
\sum_{l,\k,\q,\p} \frac{\hbar^{3/2}}{\sqrt{\om(\q)\om(\k)\om(\p)}}
e^{i(\p+\q+\k-\G)\cdot\r_l}  \\
&&\lep a(\k)+a^{\dag}(-\k) \rip
\lep a(\p)+a^{\dag}(-\p) \rip
\lep a(\q)+a^{\dag}(-\q) \rip\,.\nonumber
\eeqa
The sum over $l$ yields a delta function $\delta_{\p+\q+\k,\G}$ (ie., 
crystal momentum conservation).  
\begin{figure}
\begin{center}
\includegraphics[width=0.75\textwidth,clip=true]{threephon.pdf}
\protect\caption{\emph{Three-phonon processes resulting from cubic terms in the
inter-ion potential. Six other three-phonon processes are possible. }\label{fig:threephon}
}
\end{center}
\end{figure}
Physically, these processes correspond 
to phonon decay in which a phonon can decompose into two others.  As we 
shall see, such anharmonic processes are crucial to the calculation of the 
thermal conductivity, $\kappa$, of crystals.

\subsection{Thermal Conductivity}

     Metals predominately carry heat with free electrons, and are considered 
to be good conductors.  Insulators, which lack free electrons, predominantly 
carry heat with lattice vibrations -- phonons.  Nevertheless, some very hard 
insulating crystals have very high thermal conductivities - diamond C
\begin{table}[htb]
\begin{center}
\begin{tabular}{|l|l|l|}\hline
material/T	& 273.2K	& 298.2K 	\\ \hline
C		& 26.2		& 23.2		\\
Cu		& 4.03		& 4.01		\\ \hline
\end{tabular}
\end{center}
\caption[]{\em{The thermal conductivities of copper and diamond (CRC) (
in $\mu$Ohm-cm).}}
\label{thermcon}
\end{table}
which is often highly temperature dependent.  However, most
insulators are not good thermal conductors.  The square of the thermal conductivity also
figures into the figure of merit for thermoelectrics\cite{thermoelectrics}.  This subsection
will be devoted to understanding what makes stiff crystals like diamond
such good conductors of heat.

     The thermal conductivity $\kappa$ is measured by setting up a small 
steady thermal gradient across the material, then
\beq
Q=-\kappa\grad T
\eeq
where $Q$ is the thermal current density; i.e., the energy times the
density times the velocity.  If the thermal current is in the x-direction,
then
\beq
Q_x=\frac{1}{V}\sum_{\q,s} \hbar \om_s(\q)\langle n_s(\q)\rangle v_{sx}(\q)
\eeq
where the group velocity is given by $v_{sx}(\q)=\pde{\om_{s}(\q)}{q_x}$.
Since we assume  $\grad T$ is small, we will only look at the linear
response of the system where $\langle n_s(\q)\rangle$ deviates little from 
its equilibrium value $\langle n_s(\q)\rangle^0$.  Furthermore since  
$\om_s(\q)=\om_s(-\q)$, 
\beq
v_{sx}(-\q)=\pde{\om_{s}(-\q)}{-q_x}= -v_{sx}(\q)
\eeq
Thus as $\langle n_s(-\q)\rangle^0=\langle n_s(\q)\rangle^0$
\beq
Q_x^0=
\frac{1}{V}\sum_{\q,s} \hbar \om_s(\q)\langle n_s(\q)\rangle^0 v_{sx}(\q)=0
\eeq
since the sum is over all $\q$ in the B.Z. Thus if we expand
\beq
\langle n_s(\q)\rangle=\langle n_s(\q)\rangle^0+\langle n_s(\q)\rangle^1+\cdots
\eeq
we get
\beq
Q_x\approx
\frac{1}{V}\sum_{\q,s} \hbar \om_s(\q)\langle n_s(\q)\rangle^1 v_{sx}(\q)
\eeq
since we presumably already know $\om_s(\k)$, the calculation of $Q$ and 
hence $\kappa$ reduces to the evaluation of the linear change in $<n>$.

Within a region, $<n>$ can change in two ways.  Either phonons can diffuse 
into the region, or they can decay through an anharmonic (cubic) term into 
other modes.
\begin{figure}
\begin{center}
\includegraphics[width=0.75\textwidth,clip=true]{phonregion.pdf}
\protect\caption{\emph{Change of phonon density within a trapazoidal region.
$\langle n_s(\q)\rangle$ can change either by phonon decay or by
phonon diffusion into and out of the region. }\label{fig:phonregion}
}
\end{center}
\end{figure}
so
\beq
\der{\langle n\rangle}{t} =
\left.\pde{\langle n\rangle}{t}\right|_{\rm{diffusion}} +
\left.\pde{\langle n\rangle}{t}\right|_{\rm{decay}}
\eeq
However $\der{\langle n\rangle}{t}=0$  since we are in a steady state.  
The decay process is usually described by a relaxation time  $\tau$
(or a mean-free path $l=v\tau$ )
\beq
\left.\pde{\langle n\rangle}{t}\right|_{\rm{decay}}
=-\frac{<n>-<n>^0}{\tau}=\approx -\frac{<n>^1}{\tau}
\eeq
The diffusion part of $\der{\langle n\rangle}{t}$ is addressed pictorially
in Fig.~\ref{fig:phondiff}.
\begin{figure}
\begin{center}
\includegraphics[width=0.75\textwidth,clip=true]{phondiff.pdf}
\protect\caption{\emph{Phonon diffusion.  In time $\Delta t$, all the phonons in the
left, source region, will travel into the region of interest on the right,
while those on the right region will all travel out in time $\Delta t$.
Thus, $\Delta n/\Delta t= (n_{left}-n_{right})/\Delta t$.}\label{fig:phondiff}
}
\end{center}
\end{figure}
Formally,
\beqa
\left.\pde{\langle n\rangle}{t}\right|_{\rm{diffusion}}
&\approx &
\frac{\langle n(x-v_x\Delta t)\rangle-\langle n(x)\rangle}{\Delta t}
\approx
-\pde{\langle n(x)\rangle}{x} v_x\\
&\approx & -v_x \pde{\langle n(x)\rangle}{T}\pde{T}{x}
\approx  -v_x \pde{\lep\langle n(x)\rangle^0+ \langle n(x)\rangle^1\rip}{T}\pde{T}{x}
\nonumber
\eeqa
Keeping only the lowest order term,
\beq
\left.\pde{\langle n\rangle}{t}\right|_{\rm{diffusion}}
\approx -v_x \pde{\langle n(x)\rangle^0}{T}\pde{T}{x}
\eeq
Then as $\der{\langle n\rangle}{t}=0$,
\beq
\left.\pde{\langle n\rangle}{t}\right|_{\rm{diffusion}}
=
-\left.\pde{\langle n\rangle}{t}\right|_{\rm{decay}}
\eeq
or
\beq
\langle n\rangle^1 = - v_x\tau_s(\q) \pde{\langle n(x)\rangle^0}{T}\pde{T}{x}\,.
\eeq
Thus
\beq
Q_x\approx
-\frac{1}{V}\sum_{\q,s} \hbar \om_s(\q) v_{sx}^2(\q)
\tau_s(\q) \pde{\langle n(x)\rangle^0}{T}\pde{T}{x}\,,
\eeq 
and since $Q=-\kappa\grad T$ 
\beq
\kappa \approx
\frac{1}{V}\sum_{\q,s} \hbar \om_s(\q) v_{sx}^2(\q)
\tau \pde{\langle n(x)\rangle^0}{T}\,.
\eeq 

From this relationship we can learn several things.  First since
$\kappa\sim v_{sx}^2(\q)$, phonons near the zone boundary or optical modes 
with small $\v_s(\q)=\grad_q\om_s(\q)$ contribute little to the thermal conductivity.  Also, stiff materials, with very fast speed of the acoustic
modes $v_{sx}(\q)\approx c$  will have a large $\kappa$. Second, since 
$\kappa\sim\tau_s(\q)$, and $l_s(\q)=v_{sx}(\q)\tau_s(\q)$,
$\kappa$ will be small for materials with short 
mean-free paths.  The mean-free path is effected by defects, anharmonic 
Umklapp processes, etc.  We will explore this effect, especially its 
temperature dependence, in more detail.

     At low $T$, only low-energy, acoustic, modes can be excited (those with 
$\hbar\om_s(\q)\sim k_B T$).  These modes have
\beq
v_s(q)=c_s
\eeq
In addition, since the momentum of these modes $q\ll G$, we only have to
worry about anharmonic processes which do not involve a reciprocal lattice 
vector $\G$ in lattice momentum conservation.  Consider one of the three-phonon
anharmonic processes of phonon decay shown in Fig.~\ref{fig:threephon} (with
$\G=0$).  For these processes $Q\sim\hbar\om c$ so the thermal 
current is not disturbed by anharmonic processes.
Thus the anharmonic terms at low T do not affect the mean-free path, so 
the thermal resistivity (the inverse of the conductivity) is dominated 
by scattering from impurities in the bulk and surface imperfections
at low temperatures.

     At high $T$  momentum conservation in an anharmonic process
may involve a reciprocal lattice vector $\G$ if the $q_1$ of an excited
mode is large enough and there exists a sufficiently small $\G$ so that
$q_1>G/2$ (c.f.\ Fig.~\ref{fig:umklapp}).  
\begin{figure}
\begin{center}
\includegraphics[width=0.75\textwidth,clip=true]{umklapp.pdf}
\protect\caption{\emph{Umklapp processes involve a reciprocal lattice vector
$\G$ in lattice momentum conservation.  They are possible whenever
$q_1 > G/2$, for some $\G$, and involve a virtual reversal of the
momentum and heat carried by the phonons (far right).}\label{fig:umklapp}
}
\end{center}
\end{figure}
This is called an Umklapp process, and it involves a very large
change in the heat current (almost a reversal).  Thus the mean-free path
$l$ and $\kappa$ are very much smaller for high temperatures where $\q_1$
can be larger than half the smallest $\G$.

     So what about diamond?  It is very hard and very stiff, so the sound
velocities $c_s$ are
large, and so thermally excited modes for which $k_B T\sim\hbar\om\sim\hbar cq$
involve small $\q_1$  for which Umklapp processes are irrelevant. 
Second $\kappa\sim c^2$ which is large.  Thus $\kappa$ for diamond is huge!


\begin{thebibliography}{00}

%% \bibitem must have the following form:
%%   \bibitem{key}...
%%

\bibitem{thermoelectrics} https://en.wikipedia.org/wiki/Thermoelectric{\textunderscore}materials

\end{thebibliography}

\edo
