\documentclass[11pt,fleqn]{article}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% template that does not use Revtex4
%%% but allows special fonts
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%% Please use this template.
%%% Edit it using e.g. Notepad
%%% Ignore the header (do not change it)
%%% Process the file in the Latex site 
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Page setup
\topmargin -1.5cm      
\oddsidemargin -0.04cm   
\evensidemargin -0.04cm  
\textwidth 16.59cm
\textheight 24cm 
\setlength{\parindent}{0cm} 
\setlength{\parskip}{0cm} 


% Fonts
\usepackage{latexsym}
\usepackage{amsmath} 
\usepackage{amssymb} 
\usepackage{bm}
\usepackage[pdftex]{graphicx}

% Math symbols I
\newcommand{\sinc}{\mbox{sinc}}
\newcommand{\const}{\mbox{const}}
\newcommand{\trc}{\mbox{trace}}
\newcommand{\intt}{\int\!\!\!\!\int }
\newcommand{\ointt}{\int\!\!\!\!\int\!\!\!\!\!\circ\ }
\newcommand{\ar}{\mathsf r}
\newcommand{\im}{\mbox{Im}}
\newcommand{\re}{\mbox{Re}}

% Math symbols II
\newcommand{\eexp}{\mbox{e}^}
\newcommand{\bra}{\left\langle}
\newcommand{\ket}{\right\rangle}

% Mass symbol
\newcommand{\mass}{\mathsf{m}} 
\newcommand{\Mass}{\mathsf{M}} 

% More math commands
\newcommand{\tbox}[1]{\mbox{\tiny #1}}
\newcommand{\bmsf}[1]{\bm{\mathsf{#1}}} 
\newcommand{\amatrix}[1]{\begin{matrix} #1 \end{matrix}} 
\newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}}

% Other commands
\newcommand{\hide}[1]{}
\newcommand{\drawline}{\begin{picture}(500,1)\line(1,0){500}\end{picture}}
\newcommand{\bitem}{$\bullet$ \ \ \ }
\newcommand{\Cn}[1]{\begin{center} #1 \end{center}}
\newcommand{\mpg}[2][1.0\hsize]{\begin{minipage}[b]{#1}{#2}\end{minipage}}
\newcommand{\Dn}{\vspace*{3mm}}

% Figures
\newcommand{\putgraph}[2][0.30\hsize]{\includegraphics[width=#1]{#2}}

% heading
\newcommand{\exnumber}[1]{\newcommand{\exnum}{#1}}
\newcommand{\heading}[1]{\begin{center} {\Large {\bf Ex\exnum:} #1} \end{center}}
\newcommand{\auname}[1]{\begin{center} {\bf Submitted by:} #1 \end{center}}


\begin{document} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\exnumber{2215}
\heading{Heat capacity of solids} 
\auname{Asaf Szulc}


{\bf The problem:}

\Dn

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Consider a piece of solid whose low laying 
excitations are bosonic modes that have 
spectral density $g(\omega)=C\omega^{\alpha-1}$
up to a cutoff frequency ${\omega_c}$, 
as in the well-known Debye model (items 1-5). 
Similar description applies for magnetic materials (item 6).
%
%
In items 7-8 assume that the solid is a "glass", 
whose low laying excitations are like two level entities 
that have a spectral density $g(\omega)$. 


\begin{itemize}

\item[(1)] 
Write a general expression for the energy $E(T)$ of the system.
This expression may involve a numerical prefactor that is defined 
by an $\alpha$ dependent definite integral.   

\item[(2)] 
Write a general expression for the heat capacity $C(T)$.

\item[(3)] 
Write a general expression for the variance $\text{Var}(r)$ 
of an atom that reside inside the solid.

\item[(4)]
Determine what are $\alpha$ and $C$ and $\omega_c$ 
for a piece of solid that consists of $N$ atoms that occupy a 
volume $L^d$ in $d=1,2,3$ dimensions, 
assuming a dispersion relation $\omega=c|k|$, as for "phonons".  

\item[(5)]
Write explicitly what are $C(T)$ and $\text{Var}(r)$ 
for ${d=1,2,3}$. Be careful with the evaluation of $\text{Var}(r)$.  
In all cases consider both low temperatures ($T\ll \omega_c$), 
and high temperatures ($T\gg \omega_c$). 

\item[(6)]
Point out what would be $\alpha$ if the low laying excitations 
had a dispersion relation $\omega=a|k|^2$ as for "magnons".   

\item[(7)]
What is the heat capacity of a "glass" whose two level 
entities have excitation energies $\omega = \Delta$,
where $\Delta$ has a uniform distribution with density $C$.    
 
\item[(8)]
What is the heat capacity of a "glass" whose two level 
entities have excitation energies $\omega = \omega_c \exp(-\Delta)$, 
where the barrier $\Delta>0$ has a uniform distribution with density $D$.   

\end{itemize}

\Dn\Dn


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf The solution:}

\Dn

(1) 
Taking the bosonic occupation function with $\mu=0$ and the spectral density $g(\omega)=C\omega^{\alpha-1}$, then integrating up to a cutoff frequency ${\omega_c}$ yields:
\begin{equation}
	E = \int_0^\infty{g(\epsilon)\epsilon d\epsilon f(\epsilon)} = 
	\int_0^{\omega_c}{g(\omega)\omega d\omega \frac{1}{e^{\omega/T}-1}} = 
	C \int_0^{\omega_c} {\frac{\omega^{\alpha} d\omega}{e^{\omega/T}-1}}
\end{equation}

%%%%%%%%%%%%
\Dn

(2) 
The heat capacity is the derivative of the energy with respect to the temperature:
\begin{equation}
	\label{heat_capacity}
	C(T) = \frac{dE}{dT} = C \int_0^{\omega_c} {\frac{e^{\omega/T}}{(e^{\omega/T}-1)^2} \left(\frac{\omega}{T}\right)^2 \omega^{\alpha -1} d\omega} \equiv
	CT^{\alpha} F\left(\frac{\omega_c}{T}\right)
\end{equation}

%%%%%%%%%%%%
\Dn

(3)
The Hamiltonian of our problem is:
\begin{equation}
	{\cal H} = \sum_{i=1}^N {\frac{p_i^2}{2m}} +\sum_{\lbrace ij\rbrace (nn)}K_{ij}(x_i-x_j)^2 \quad \longmapsto \quad
	{\cal H} = \sum_k {\frac{P_k^2}{2m}} +\sum_k {\frac{1}{2}m\omega_k^2Q_k^2}
\end{equation}
where we were lucky enough that the basis in which the Hamiltonian is diagonalized is the Discrete Fourier Transform of $x_i$, $p_i$ given by:
\begin{equation}
	Q_k = \frac{1}{\sqrt{N}} \sum_l {e^{ikl}x_l}\ ; \quad\quad 
	P_k = \frac{1}{\sqrt{N}} \sum_l {e^{-ikl}p_l}
\end{equation}
The Hamiltonian was written for a $1D$ case, however, the generalization to other dimensions is trivial.
Now, as we expressed the position as a sum of independent random variables, we can write:
\begin{equation}
	Var(r) = \sum_k Var(Q_k) = \sum_k {\frac{3\omega f(\omega)}{m\omega^2}} = 
	\frac{3C}{m} \int_{\frac{2\pi c}{L}}^{\omega_c} {\frac{\omega^{\alpha-2}d\omega}{e^{\omega/T}-1}}
\end{equation}
where we implemented the equipartition theorem for low temperatures ($\omega f(\omega)/2$ for each degree of freedom), then we used the given density of states for the transformation to an integral and defined minimum and maximum cutoff frequencies.
%%%%%%%%%%%%
\Dn

(4) 
Using the connection, $E = \omega = c|k|$, it is convenient to calculate the density of states in k-space.
States with equal energy lie on a d-dimensional spherical shell in k-space, divided by the volume of
a unit cell gives the density of states:
\begin{equation}
	\label{density}
	g(k)dk = \left(\frac{d}{2^d}\right) \frac{\Omega_d k^{d-1}dk}{\left(\frac{\pi}{L}\right)^d} 
	\quad \Longrightarrow \quad
	g(\omega)d\omega = \frac{L^{d}d}{(2\pi c)^{d}} \Omega_{d} \omega^{d -1}d\omega
\end{equation}
where $\Omega_{d} = \{2, 2\pi, 4\pi\}$ for $d = \{1,2,3\}$ and the pre-factor is due to the consideration of 'd' different modes of the lattice (one for each degree of freedom) and only positive values of 'k'.\\

And now to business, first we will plug eqn.(\ref{density}) in the number of states equation. For a d-dimensional crystal with $N$ atoms there are $Nd$ normal modes so we can write:
\begin{equation}
	\label{states}
	Nd = \int_0^{\omega_c} {g(\omega)d\omega} = 
	\frac{\Omega_{d}d}{(2\pi c)^{d}} L^{d} \int_0^{\omega_c} {\omega^{d -1}d\omega}
\end{equation}
second, in the equation for the heat capacity to get:
\begin{equation}
	\label{heat}
	C(T) = \frac{d}{dT}\left(\int_0^{\omega_c}{g(\omega)\omega d\omega \frac{1}{e^{\omega/T}-1}} \right) = 
	\frac{\Omega_{d}d}{(2\pi c)^{d}} L^{d} \int_0^{\omega_c} {\frac{e^{\omega/T}}{(e^{\omega/T}-1)^2} \left(\frac{\omega}{T}\right)^2 \omega^{d -1} d\omega}
\end{equation}
Solving eqn.(\ref{states}) and comparing eqn.(\ref{heat}) to eqn.(\ref{heat_capacity}) leads to the final result:
\begin{equation}
	\omega_c = \frac{2\pi c}{L}\left(\frac{Nd}{\Omega_{d}}\right)^{1/d} \ ; \quad\quad
	C = \frac{\Omega_{d}d}{(2\pi c)^{d}} L^{d} = \frac{Nd^2}{\omega_c^d} \ ; \quad\quad 
	\alpha = d
\end{equation}


%%%%%%%%%%%%
\Dn

(5)
In this section, throughout our entire calculation for the low temperature limit, $T \ll \omega_c$, we take $x_c =
\omega_c/T \rightarrow \infty$ and integrating $F(x_c)$ by parts. For the high temperature limit, $T \gg \omega_c$, we take $x \rightarrow 0$ $\Rightarrow$ $e^x \approx 1 + x$ and we keep only the lowest order in $x$. Thus we get:\\

%%% ---------- 1D ---------- %%%
For $1D$ ($\alpha=1$, $\Omega_{\alpha}=2$):\\
\begin{equation}
	C(T) = \frac{NT}{\omega_c} F(x_c) = \left\{ 
	\begin{array}{l l}
	\frac{\pi}{3}\frac{LT}{c} & \quad T \ll \omega_c\\
	N & \quad T \gg \omega_c
	\end{array} \right.
\end{equation}

\begin{equation}
	Var(r) = \frac{3L}{m\pi c} \int_{\frac{2\pi c}{L}}^{\omega_c} {\frac{d\omega}{\omega(e^{\omega/T}-1)}} = \left\{ 
	\begin{array}{l l}
	0 & \quad T \ll \omega_c\\
	\frac{3}{2\pi^2}\left(\frac{T}{mc^2}\right)L^2 & \quad T \gg \omega_c
	\end{array} \right.
\end{equation}

%%% ---------- 2D ---------- %%%
For $2D$ ($\alpha=2$, $\Omega_{\alpha}=2\pi$):\\
\begin{equation}
	C(T) = \frac{4N}{x_c^2} F(x_c) = \left\{ 
	\begin{array}{l l}
	7.212\cdot\frac{L^2T^2}{\pi c^2} & \quad T \ll \omega_c\\
	2N & \quad T \gg \omega_c
	\end{array} \right.
\end{equation}

\begin{equation}
	Var(r) = \frac{3L^2}{m\pi c^2} \int_{\frac{2\pi c}{L}}^{\omega_c} {\frac{d\omega}{(e^{\omega/T}-1)}} = \left\{ 
	\begin{array}{l l}
	0 & \quad T \ll \omega_c\\
	\frac{3ln(N)}{2\pi}\left(\frac{T}{mc^2}\right)L^2 & \quad T \gg \omega_c
	\end{array} \right.
\end{equation}
\\
%%% ---------- 3D ---------- %%%
For $3D$ ($\alpha=3$, $\Omega_{\alpha}=4\pi$):\\
\begin{equation}
	C(T) = \frac{9N}{x_c^3} F(x_c) = \left\{ 
	\begin{array}{l l}
	\frac{2\pi^2L^3}{5c^3}T^3 & \quad T \ll \omega_c \quad\quad (Debye\ T^3\ law)\\
	3N & \quad T \gg \omega_c \quad\quad (Dulong\ Petit\ law)
	\end{array} \right.
\end{equation}

\begin{equation}
	Var(r) = \frac{9L^3}{2m\pi^2 c^3} \int_{\frac{2\pi c}{L}}^{\omega_c} {\frac{\omega d\omega}{(e^{\omega/T}-1)}} = \left\{ 
	\begin{array}{l l}
	0 & \quad T \ll \omega_c\\
	\frac{9}{\pi}\left(\frac{3N}{4\pi}\right)^{\frac{1}{3}} \left(\frac{T}{mc^2}\right)L^2  & \quad T \gg \omega_c
	\end{array} \right.
\end{equation}

%%%%%%%%%%%%
\Dn

(6)
The calculation of the new density of states is similar to the previous problem. Therefore, we continue by plugging the new dispersion relation $\omega = a|k|^2$ into eqn.(\ref{density}) and get:
\begin{equation}
	g(\omega)d\omega = \Omega_{d} L^{d}\frac{md}{(2\pi)^{d}} \left(\frac{\omega}{a}\right)^{d/2 - 1} \frac{d\omega}{2ma} 
	\quad \Longrightarrow \quad 
	g(\omega) = C\omega^{d/2 - 1} \equiv C\omega^{\alpha-1}
\end{equation}
hence $\alpha = d/2$.

%%%%%%%%%%%%
\Dn

(7)
The partition function for a two level system is:
\begin{equation}
	Z(\beta) = 1+e^{-\beta\omega}
\end{equation}
leading to the same calculation but with "+" instead of "-" sign. Plugging in $g(\Delta) = C$ to get:
\begin{equation}
	C(T) = \int_0^{\infty} {\frac{g(\Delta)e^{\Delta /T}}{(e^{\Delta /T} +1)^2} \left(\frac{\Delta}{T}\right)^2 d\Delta} = 
	-CT \int_0^{\infty} {\frac{x^2e^{x}}{(e^{x} +1)^2} dx} = 
	\frac{\pi^2}{6}CT
\end{equation} 
where we defined $x = -\Delta/T$ and integrated by parts.
%%%%%%%%%%%%
\Dn

(8)
Here we use the given connection, $\omega = \omega_c e^{-\Delta}$ and then introduce another change of variables,
$x = -\omega_c e^{-\Delta}/T$ , to get:
\begin{equation}
	C(T) = \int_0^{\infty} {\left(\frac{\omega_c e^{-\Delta}}{T}\right)^2\frac{g(\Delta)e^{-\frac{\omega_c}{T}e^{-\Delta}}}{(e^{-\frac{\omega_c}{T}e^{-\Delta}}+1)^2}d\Delta} = 
	 \int_0^{-x_c} {\frac{Dxe^xdx}{(e^x+1)^2}} = \left\{
	\begin{array}{l l}
	ln(2)D & \quad T \ll \omega_c\\
	\frac{\omega_c^2}{8T^2}D  & \quad T \gg \omega_c
	\end{array} \right.
\end{equation} 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document} 
