Stat_Mech.tex

\documentclass{article}
\usepackage{amsthm}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{wrapfig}
\usepackage{tabu}
\usepackage{hyperref}
\graphicspath {{Images/}}

\makeatletter
\def\lecture{\@ifnextchar[{\@lectureWith}{\@lectureWithout}}
\def\@lectureWith[#1]{\medbreak\refstepcounter{section}%
  \renewcommand{\leftmark}{Lecture \thesection}
  \noindent{\addcontentsline{toc}{section}{Lecture \thesection: #1\@addpunct{.}}%
  \sectionfont Lecture \thesection. #1\@addpunct{.}}\medbreak}
\def\@lectureWithout{\medbreak\refstepcounter{section}%
  \renewcommand{\leftmark}{Lecture \thesection}
  \noindent{\addcontentsline{toc}{section}{Lecture \thesection.}%
  \sectionfont Lecture \thesection.}\medbreak}
\makeatother
\newcommand{\pardif}[2]{\frac{\delta#1}{\delta#2}}
\newcommand{\secpardif}[2]{\frac{\delta^{2}#1}{\delta#2^{2}}}
\newcommand{\thermdif}[4]{\frac{\delta#1}{\delta#2}\vert_{#3,#4}}
\newcommand{\secthermdif}[4]{\frac{\delta^{2}#1}{\delta#2^{2}}\vert_{#3,#4}}
\newcommand{\dif}[2]{\frac{d#1}{d#2}}
\newcommand{\bltz}{k_{B}}
\newcommand{\sumser}[2]{\sum\limits_{#1}^{#2}}
\newcommand{\dtk}{\frac{d^{3}k}{2\pi^{3}}}	
\newcommand{\ebk}{e^{-\beta\epsilon_{k}}}
\newcommand{\qmom}{\frac{\hbar^{2}}{2m}}
\newcommand{\brad}{R_{\nu}(T)}
\newcommand{\braddiff}{R_{\nu}(T)d\nu}
\newcommand{\mbb}[1]{\mathbb{#1}}
\newcommand{\aras}{\hat{a}^{\alpha^{\dagger}}(\mbb{k})}
\newcommand{\alow}{\hat{a}^{\alpha}(\mbb{k})}

\title{Lectures on Statistical Mechanics and Thermodynamics}
\author{Paolo Bedaque}
\begin{document}
\maketitle

\part{Foundations of Statistical Mechanics}

%\lecture
The object of Statistical Mechanics is to study systems with a large number of degrees of freedom by focusing solely on the macroscopic properties of the system.  

\begin{figure}
	\includegraphics[width=0.5\linewidth]{Images/mic_to_mac.png}
	\caption{Multiple microscopic states make up a single macroscopic state of a system.}
\end{figure}

Why, you may ask, are we only focusing on these macroscopic properties?  Surely, given enough information we could model any system by following its individual members!  This is not necessarily the case because the sheer volume of information necessary to do such a calculation is in fact an impediment to performing such a calculation.  Let's illustrate this using a simple example of a system consisting of 20 particles in the ground state.  Suppose, through divine intervention or otherwise, we have obtained the complete wave equation for each particle in the system.  Before we begin to calculate any properties of this system we must first store this information.  

\begin{figure}
\centering
	\begin{minipage}{0,49\textwidth}
		\frame{\includegraphics[width=0.5\linewidth]{Images/fig_1_1.png}}
		\caption{Wavefunction of our system in dimension $x_{1}$}
	\end{minipage}
%	\begin{minipage}{0.49\textwidth}
%		\frame{\includegraphics[\width=\linewidth]]{imagefile}}
%		\caption{test}
%	\end{minipage}
\label{fig:information_graph}
\end{figure}

As shown in figure 2, every particle wavefunction will be reasonably approximated by 10 points in each dimension.  Taking the 3-dimensional case, we find that the amount of memory we would need would be $(10)^{3N}*8\simeq10^{61}$ bytes.  We include the factor of 8 since each real number will be represented by 8 bytes.  Assuming we store all of this data on terabyte drives that have a storage capacity of $\frac{10^{12} bytes}{0.4 kgs}$, we find that the total mass of all of the drives required to store the information describing this 20 particle system to be on the order of $10^{49}$ kilograms, which is much much more than the mass of the Milky Way galaxy!  Clearly, it is impossible to store this much information, much less manipulate it.  Therefore, we employ statistics to deal with the system at a macroscopic level. 

\section{Hamiltonian Dynamics}

Consider a hypothetical system consisting of some number of particles which can be described by their individual positions, $q_{i}$ and their momentums, $p_{i}$.  We can indicate the paths these particles travel on a {\bf{Phase Space}} diagram.  

%Insert phase space figure

This diagram displays the evolution of the systems constituent particles as their positions and momentums change.  The evolution of this system is governed by the Hamilton Equations:  $\dif{q_{1}}{t} = \pardif{\mathcal{H}}{p_{i}}$ and $\dif{p_{i}}{t}=-\pardif{\mathcal{H}}{q_{i}}$  The evolution of $p$ and $q$ are clearly deterministic, therefore we can express these values as time-dependant functions.  $q_{i}(t), p_{i}(t) \rightarrow q_{i}(t=0)=q_{i}^{0}$ and $p_{i}(t=0)=p_{i}^0$

There are a couple of tricks we can employ to idealize the measurement process and simplify manipulations of these functions.

\section{Idealization of the Measurement Process}

Let's start by defining some value $\overline{f}$, which is the result of measuring function $f(p_{i},q_{i})$.  

$$\overline{f}=\lim_{T\to\infty}\frac{1}{T}\int_{\delta}^{\delta+T}dtf(q(t),p(t))=\widetilde{f}(q_{i}^0,p_{i}^0)$$

$T=$ duration of the measurement, $\delta=$ time the measurement starts, $\widetilde{f}=$ some function of the initial conditions.

This expression is an idealization of an arbitrary measurement.  In reality, the time T is much larger than the microscopic scales in which the functions of $q(t)$ and $p(t)$ vary.  For example, suppose the function $f$ is the pressure on the wall of a gas containing vessel.  The quantity $f(q(t),p(t))$ flucuates like the diagram below.

%Insert pressure fluctuation diagram

However, except in specialized cases, the measuring apparatus is not accurate enough the capture every tiny fluctuation of the force on the wall of the vessel.  Instead, it measures an average of $f(q(t),p(t))$ over a finite period of time:

%Insert modified pressure fluctuation diagram

Taking the $\lim_{T\to\infty}$ is a mathematical idealization of this averaging out of function fluctuations.  Returning to our analysis of the hypothetical measurement $\overline{f}$, we see that the function $\widetilde{f}(q_{i}^0,p_{i}^0)$ is not really a function of the initial condition.  Rather, it is a function of the \textit{trajectory} of the system evolution.  This may not be immediately apparent, but can be illustrated by choosing another initial point on the trajectory of our function $\widetilde{f}=(q_{i}^0,p_{i}^0)$.  As $T\to\infty$, we see that our new initial points $q_{i}^{0'}$ and $p_{i}^{0'}$ yield the same value of $\widetilde{f}$ as $q_{i}^0$ and $p_{i}^0$.  i.e., $\widetilde{f}(q_{i}^0,p_{i}^0)=\widetilde{f}(q_{i}^{0'},p_{i}^{0'})$  This is because the difference between the two trajectories is neglible compared to the infinite extended trajectories as $T\to\infty$.

%Insert trajectory phase diagram

Now, it could be possible that two different trajectories give different values for $\overline{f}$.

%Insert IC Trajectory phase diagram

However, it turns out that, for systems we are interested in, that it does not happen.  These systems have the property that any individual trajectory will pass through almost \textbf{every} point in phase space.  This is called the \textbf{Ergodic Property}.  

A couple of important things to keep in mind.  The assumption that the trajectory passes through every point in the phase space is contingent on the the assumption that these points have the same value of energy (and any other conserved quantities) as the initial point.  After all, these values are conserved by the Hamiltonian Flow.  The phrase "almost every point" also has a specific meaning, in that it is "almost every point" in the measure theory sense.

It is really difficult to prove that a realistic system has the ergodic property.  Only a few systems have actually been proven rigorously to be ergodic (for example, the Sinai Stadium).  Moving forward, we will operate on the assumption that the systems we are interested are, in fact, ergodic.    

To expand on the concept of ergodicity, it is very easy to find systems that are not ergodic.  Any conserved value restricts the trajectories to lie on a submanifold of the phase space.  If the only conserved quantities are the ones resulting from the standard symmetries (i.e. energy, momeentum, etc.) we can question the validity of ergodicity within the restricted subspace with constant values of these quantities.  There are systems with so many other conserved quantities that the subspace is one dimensional (integrable systems), thereby invalidating the ergodicity of the system.  The Integrable/Ergodic classification is also non-binary, there are intermediate types of systems that exist between the two extremes, such as mixing systems.

\section{Introduction to Thermodynamics}

At its core level, many thermodynamic problems can be phrased as a comparison between two states of a system.  The pre-evolutionary equilibrium state and the post-evolutionary equilibrium state.

%Insert Thermoproblem breakdown diagram.

Since every microstate of the system is equally probably, the final macroscopic state will be the one which corresponds to the largest number of microstates:

$$\Gamma(E_{1}',V_{1}',N_{1}',E_{2}',V_{2}',N_{2}')=\Gamma_{1}(E_{1}',V_{1}',N_{1}')\Gamma_{2}(E_{2}',V_{2}',N_{2}') $$

Where $\Gamma(E_{1}',V_{1}',N_{1}',E_{2}',V_{2}',N_{2}')$ is the number of states of the composite system, $\Gamma_{1}(E_{1}',V_{1}',N_{1}')$ is the number of microstates of subsystem 1, and $\Gamma_{2}(E_{2}',V_{2}',N_{2}')$ is the number of microstates of subsystem 2.

Let's look at the case where the wall between the two subsystems allows the exchange of energy but not of volume or particle number.  In this case we find:

\begin{multline}
E_{1}+E_{2}=E_{1}'+E_{2}' \\
V_{1}=V_{1}', V_{2}=V_{2}' \\
N_{1}=N_{1}', N_{2}=N_{2}' \\
\Gamma(E_{1}',E_{2}')=\Gamma_{1}(E_{1}')\Gamma_{2}(E_{2}')
\end{multline}

This implies that

$$\Gamma_{Max, E_{total} fixed}\Leftrightarrow\frac{1}{\Gamma_{1}}\dif{\Gamma_{1}(E_{1}')}{E_{1}'}-\frac{1}{\Gamma_{2}(E_{2}')}\dif{\Gamma_{2}(E_{2}')}{E_{2}'}=0$$  

So, knowing the total number of initial microstates of the two subsystems allows us to predict the final state of the overall system.  This is summed up in a very clean and useful fashion in the definition of \textbf{entropy}:

$$S=\bltz\ln(\Gamma)$$
   
Notice that entropy is extensive as long as the subsystems are large enough that the boundary effects are negligible.

\textit{How can Entropy grow?}

Within the constraints of the problem, there is an apparent paradox in staying that the entropy is maximized.  That's because, using Hamiltonian equations, we can show S is a constant.

Define $\rho(q,p,t=0)$ to be an initial distribution of identical systems (this is called an \textbf{ensemble}.)  As the Hamiltonian Flow evolves in time every element of the ensemble is carried with it and the distribution evolves to $\rho(q,p,t)=\rho(q(t),p(t))$.

%Insert flow evolution diagram

Because the number of systems in the ensemble is fixed $\rho$ obeys the continuity equation:

$$\dif{\rho}{t}+\nabla*(\rho\nabla)=0$$
or
\begin{equation}
\begin{split}
0=\pardif{\rho}{t}+\sum_{i=1}^{3N}[\pardif{(\rho\dot{q_{i}})}{q_{i}}+\pardif{\rho\dot{p_{i}}}{p_{i}}] \\
	=\pardif{\rho}{t}+
\end{split}	 
\end{equation}
%Finish Equation

This result is known as the \textbf{Liouville Theorem}.  It implies that the density of elements on the ensemble does not change as it evolves along the Hamiltonian Flow (of course, it does change at any fixed position in phase space).  As the number of systems in the ensemble is fixed, the Liouville Theorem implies that the volume of the phase space occupied by the ensemble is fixed.  

This means that the entropy, being the log of the volume, cannot change either.  So how can the entropy grow if the microscopic equations of motion show that it does not?  A resolution to this paradox is to notice that, for systems in which statistical mechanics holds, the Hamiltonian Flow takes nice looking, civilized ensembles and turns them into convoluted shapes:

%Insert convoluted flow diagram

The average of smooth observables (i.e. the value of $f(q,p)$) does not distinguish between the average performed with the true $\rho(t)$ or the "coarse grained" $\widetilde{\rho}(t)$, as shown below.

%Insert convoluted flow diagram 2

$$\overline{f}=\int_{}^{}dqdpf(q,p)\rho(q,p,t)\simeq\int_{}^{}dqdpf(q,p)\widetilde{\rho}(q,p,t)$$

For all practical purposes, the entropy grows so long as we only look at observables that are very smooth.  Macroscopic observables don't care about the precise poition of the particles and tend to be smooth in phase space.  Systems whose Hamiltonian Flow have this property are called "\textbf{Mixing}."  It is easy to show that ergodic systems are mixing, but keep in mind that mixing systems are not necessarily ergodic.    

\part{Thermodynamics}

From the general property of entropy we can derive a massive amount of information about a system.  In fact, if the entropy as a function of extensive variables is known, everything else about the thermodynamic behavior of the system can be found!  The definition of entropy is reiterated below:

$$S=\bltz\ln(\Gamma)$$

We include the Boltzmann Constant, $\bltz$, for historical reasons.  However, if we judiciously choose our units (by measuring temperature in units of energy, as we should) we can set $\bltz=1$.  A good thing to remember in this regard is that $10,000K\simeq1eV$.  

Since entropy is so important, let's review a few properties satisfied by the entropy:

\begin{itemize}
	\item S as a function of the \underline{extensive} variables contains all the information about a system that you need.
	\item S is itself extensive: $S(\lambda*E,\lambda*V,\lambda*N)=\lambda*S(E,V,N)$
	\item $\pardif{S}{E}\geq0$  We will eventually see that this implies $T\geq0$
	\item $S(E,V,N)$ is a concave function.  We will eventually see that this implies stability: \\
	$S(E+\Delta*E,V+\Delta*V,N+\Delta*N)+S(E-\Delta*E,V-\Delta*V,N-\Delta*N)\leq2S(E,V,N)$\\ 
		%insert plot of S showing concavity
	\item $S=0$ at $\frac{1}{\pardif{S}{E}}=0 (T=0)$  This property depends on how we normalize the number of microstates, $\Gamma$.  Quantum mechanics determines the right normalization that implies $S(T=0)-0$
	\item $S$ is maximized in equilibrium within the constraints set by the physical setup.  This is a direct consequence of the fact that all microstates are equally probable.
\end{itemize}

Inverting the entropy function $S(E,V,N)$ allows us to obtain a function for the energy of the system in terms of entropy: $E(S,V,N)$.  If $S,V,N$ are changed slowly (i.e. \textbf{"adiabatically"}) so that the system stays in equilibrium at all intermediate states we find:

$$dE=\thermdif{E}{S}{V}{N}dE+\thermdif{S}{V}{E}{N}dV+\thermdif{S}{N}{E}{V}dN$$

where

$$\thermdif{S}{E}{V}{N}=\frac{1}{T}, \thermdif{S}{V}{E}{N}=\frac{P}{T}, \thermdif{S}{N}{E}{V}=\frac{\mu}{T}$$

To derive the first relation, we used $\thermdif{E}{S}{V}{N}=\frac{1}{den}$

%Finish equation of state section

It's important to note that S is defined for equilibrium macroscopic states only.  So dS or dE refer to changes slow enough such that the system continuously remains in equilibrium during its transition from its initial state to its final state.  i.e., at any given time during the transition the system is in an equilibrium state.  For instance, if the volume is changed suddenly and then the system is left to relax to equilibrium the change in entropy is not given by $\frac{p}{T}dV$ (it will, in fact, be larger). 

From these relations we can derive the Equations of State:

\begin{tabu} to 0.8\textwidth { | X[l] | X[r] | }
	\hline
	$T=\thermdif{E}{S}{V}{N}=T(S,V,N)$ & $\frac{1}{T}=\thermdif{S}{E}{V}{N}=\frac{1}{T(E,V,N)}$ \\
	\hline
	$-P=\thermdif{E}{V}{S}{N}=P(S,V,N)$ & $-\frac{P}{T}=\thermdif{S}{V}{E}{N}=-\frac{P(E,V,N)}{T(E,V,N)}$ \\
	\hline
	$\mu=\thermdif{E}{N}{S}{V}=\mu(S,V,N)$ & $\frac{\mu}{T}=\thermdif{S}{N}{E}{V}=\frac{\mu(E,V,N)}{T(E,V,N)}$ \\
	\hline
\end{tabu}

These equations of state are particularly important because they are frequently accessible to experiments.  Notice that the same quantities as functions of other variables are \underline{not} equations of state and may not contain the same information. 

As mentioned earlier, all thermodynamical problems can be reduced to a transition between two states of a system:

%Insert diagram

Depending on the properties of the wall between the two subsystems, we find that:

\begin{itemize}
	\item Wall allows Energy exchange $\rightarrow T_{1}' = T_{2}'$
	\item Wall allows Volume exchange $\rightarrow V_{1}' = V_{2}'$
	\item Wall allows Particle exchange $\rightarrow \mu_{1}' = \mu_{2}'$
\end{itemize}

\section{Non-Relativistic, Non-Interacting, Monoatomic Gas}

Let us compute the entropy, S, for one example:

$$\mathcal{H}(\mathbb{r}_{i},\mathbb{p}_{i})=\sumser{i=1}{N}\frac{\mathbb{p}_{i}^{2}}{2m}, (0<\mathbb{r}_{i}<L)$$

\begin{equation}
\begin{split}
\Gamma(E,V,N)=\int\prod_{i=1}^{N}d\mathbb{r}_{i}\mathbb{p}_{i}\delta(E-sumser{i=1}{N}\frac{\mathbb{p}_{i}^{2}}{2m}) \\ 
	=[\int_{0}^{L}dx]^{3N}\int_{-\infty}^{\infty}dp_{ix}dp_{iy}dp_{iz}....dp_{nx}dp_{ny}dp_{nz}*2m\delta(2mE-\sumser{i=1}{3N}p_{i}^2) \\
	=V^{N}2mS^{3N}(R=\sqrt{2mE}) \\
	=2mV^{N}\frac{2\pi^{\frac{3N}{2}}}{\Gamma(\frac{3N}{2})}(2mE)^{\frac{3N-1}{2}} 
\end{split}
S=\bltz\ln{\Gamma}=\bltz[N\ln{V}+\frac{3N-1}{2}\ln{E}-\ln{\Gamma(\frac{3N}{2})}+\frac{3N}{2}\ln{\sum}+Constant]\overrightarrow{N\rightarrow\infty}\bltz[N\ln{V}+\frac{3N}{2}\ln{E}-\frac{3N}{2}\ln{N}+\frac{3}{2}N\ln{\prod}]=N\bltz\ln{\frac{VE^{\frac{3}{2}}}{N^{\frac{3}{2}}}}
\end{equation}
The above expression violates extensivity (2) and the Neust Postulate (5).  The violation of the Neust Postualte is not surprising: the behavior of a gas at T=0 is dictated by quantum mechanics, not classical mechanics (the atoms wavelength is much larger than the separation between them).  What is surprising is that quantum mechanics connects the calculation above even at large T.  What is missing from the classical calculation above is that the gas is composed of \underline{identical} particles (bosons or fermions, it doesn't matter at high T).  So the correct way of counting states is:

$$S=\bltz\ln{\Gamma}=\bltz\int\prod_{i=1}^{N}\frac{dq_{i}dp_{i}}{N!h^N}\delta(E-\sumser{i=1}{N}\frac{\mathbb{p}_{i}^2}{2m})$$

The denominator factor of $N!$ avoids double counting of microstates with the same particle position and the factor of $h^{N}$ becomes important when considering the classical limit of quantum gases.

As the number of particles in the system becomes very large we find that the expression becomes

$$S=\bltz*N\ln{\frac{VE^{\frac{3}{2}}}{N^{\frac{5}{2}}}}$$

From this entropy we can find the equations of state:

\begin{tabu} to 0.8\textwidth { | X[l] | X[r] | }
	\hline
	$\frac{1}{T}=\thermdif{S}{E}{V}{N}=\frac{3}{2}\frac{\bltz*N}{E}$ & $E=\frac{3}{2}\bltz*N*T$ \\
	\hline
	$-\frac{P}{T}=\thermdif{S}{V}{E}{N}=-\frac{\bltz*N}{V}$ & $PV=\bltz*NT$ \\
	\hline
	$\frac{\mu}{T}=\thermdif{S}{N}{E}{V}=\bltz\ln{\frac{VE^{\frac{3}{2}}}{N^{\frac{5}{2}}}}-\frac{5}{2}\bltz$ & $ $ \\
	\hline
\end{tabu}

Note that we found the Boyle-Mariott law without even discussing atoms colliding against walls or other kinetic theory considerations.  There's a lot behind the assumption of ergodicity.  

This leads us to....

\section{Euler Relation}

$$S(\lambda*E,\lambda*V,\lambda*N)=\lambda*S(E,V,N)$$

Taking the derivative of this expression in terms of $\lambda$ yields

$$E()\pardif{S}{E}=\frac{1}{T})+V(\pardif{S}{V}=-\frac{P}{T})+N(\pardif{S}{N}=\frac{\mu}{T})=S$$

Which can be re-arranged to yield the following equations:

$$E=TS-PV+\mu*N$$
$$S=\frac{E}{T}-\frac{PV}{T}+\frac{\mu*N}{T}$$

\subsection{Recovering the Entropy or Energy from an Equation of State}

If three equations of state are availible:

$$T(S,V,N)=\thermdif{E}{S}{V}{N}$$
$$-P(S,V,N)=\thermdif{E}{V}{S}{N}$$
$$\mu(S,V,N)=\thermdif{E}{N}{S}{V}$$

$$\rightarrow E(S,V,N)=T(S,V,N)=T(S,V,N)S-P(S,V,N)V+\mu(S,V,N)N$$

If only two equations of state are availible: $E(S,V,N)=NE(\frac{S}{N},\frac{V}{N},1)$

The equality above is due to extensivity.

$$T=\thermdif{E}{S}{V}{N}=T(S,V,N)=T(\frac{S}{N},\frac{V}{N},1)$$

$$-P=\thermdif{E}{V}{S}{N}=P(S,V,N)=P(\frac{S}{N},\frac{V}{N},1)$$

Note that there is no N!.  T and P are intensive.

$$\downarrow$$

$$E(\frac{S}{V},\frac{V}{N},1)=\int_{0}^{\frac{S}{N}}d(\frac{S}{N})T(\frac{S}{N},\frac{V}{N},1)+E(0,\frac{V}{N},1)=\int_{0}^{\frac{S}{N}}d(\frac{S}{N})T(\frac{S}{N},\frac{V}{N},1)-\int_{0}^{\frac{V}{N}}d(\frac{V}{N})P(0,\frac{V}{N},1)+E(0,0,1)$$

In the above expressions, $T(\frac{S}{N},\frac{V}{N},1)=\thermdif{E}{S}{V}{N}$ and $E(0,0,1)=E(\frac{S}{N},\frac{V}{N},1)$.  The latter equality can be found up to one constant independant of S,V, or N.

Now, if only one equation of state is availible, you cannot recover the entropy or energy.  There is simply not enough information in a single equation of state to extract these system properties.  

\subsection{EXAMPLE: Non-relativistic, Monoatomic, Non-interacting Gas}

$$S=\bltz*N\ln{\frac{VE^{\frac{3}{2}}}{N^{\frac{5}{2}}}}, E=(\frac{N^{\frac{5}{2}}}{V}e^{\frac{S}{\bltz*N}})^\frac{2}{3}=\frac{N^{\frac{5}{2}}}{V^{\frac{2}{3}}}*e^{\frac{2S}{3\bltz*N}}$$

$$T=\thermdif{E}{S}{V}{N}=\frac{2}{3\bltz*N}\frac{N^\frac{5}{3}}{V^\frac{2}{3}}e^\frac{2S}{3\bltz*N}$$

$$-P=\thermdif{E}{V}{S}{N}=-\frac{2}{3}\frac{N^\frac{5}{3}}{V^\frac{5}{3}}e^\frac{2S}{3\bltz*N}$$

$$E(\frac{S}{N},\frac{V}{N},1)=\int_{0}^{\frac{S}{N}}d(\frac{S}{N})\frac{2}{3\bltz}(\frac{N}{V})^{\frac{2}{3}}e^{\frac{2S}{3\bltz*N}}-\int_{0}^{\frac{V}{N}}d(\frac{V}{N})\frac{2}{3}(\frac{N}{V})^{\frac{5}{3}}+E(0,0,1)$$

$$E(\frac{S}{N},\frac{V}{N},1)=\frac{2}{3\bltz}(\frac{N}{V})^{\frac{2}{3}}\frac{3\bltz*N}{2}[e^{\frac{2S}{3\bltz*N}}-1]-\frac{2}{3}(\frac{-3}{2})(\frac{N}{V})^\frac{2}{3}$$

$$E(\frac{S}{N},\frac{V}{N},1)=\frac{N^{\frac{5}{3}}}{V^{\frac{1}{3}}}e^{\frac{2S}{3\bltz*N}}$$

\section{Minimum Energy Principle}

Before formally stating the minimum energy principle, let's begin with an analogy.  We know that a circle is the figure with the largest area for a given perimeter.  Of course, it is also the the figure with the smallest perimeter for a given area.  Even if a given circle was built by someone with a fixed length string trying to enclose the largest area it will still be true that the final circle minimizes the perimeter for the given area.  

The same thing is true for the maximum entropy principle.  After some constraints are changed the final state will be the one with the largest entropy for a given \textbf{total} energy (the constraint that the total energy is unchanged is always there).  But the final equilibrium state can also be characterized as the state with the smallest energy for a fixed entropy:

%Insert entropy surface diagrams

Convince yourself that S as a function of total E and subsystem energy $E_{1}^{1}$ has the general shape depicted above.  The minimum energy principle follows from the picture above.  

Now let's see how the minimum energy principle works in our archetypal example   :

%Insert thermo diagram

$$E=E_{1}(S_{1}^{1},V_{1},N_{1})+E_{2}(S_{2}^{'},V_{2},N_{2})$$

Where $S_{2}^{'}=S-S_{1}^{'}$.

Therefore, we find that the minimum energy, E, for fixed entropy, S, is:

$$\pardif{E}{S_{1}^{'}}=0 \rightarrow \pardif{E_{1}}{S_{1}^{'}}-\pardif{E_{2}}{S_{2}^{'}}=0 \rightarrow T_{1}=T_{2}$$

Which confirms what we expected: In its final state our example subsystem pair reaches an equilibrium temperature.

\section{Efficiency of Engines}

In our thermodynamic terminology, an 'engine' is a machine that extracts work from the the energy difference between two systems in contact.  To clarify, the engine takes energy from the 'hot' source, uses some of that energy to perform work, and ejects the remaining energy into the 'cold' sink.  In many cases, the 'cold' sink is the atmosphere: one example where this is true is a car engine, which throws off heat as it performs the work of moving the car.  

In designing our thermodynamic engine, we would like it to perform this energy transfer over and over again in a closed cycle.  I.e. we want the system to end up in the same macroscopic state as it started.  For our purposes right now we will assume that all changes in the system caused by the energy transfers occur slowly (adiabatically) so, that at every point in the cycle the machine is in equilibrium and its macroscopic state is specified by E,V,N or S,V,N.  It turns out, that for thermodynamic calculations it is better to use the variables T,S,N instead.  In our engine we want to keep N fixed by having impermeable walls separating the machine from the rest of the Universe.  This machine cycle can be represented by a closed loop in the TxS space:

%Insert engine cycle diagram

Let us start by considering one specific cycle, the \textbf{Carnot Cycle}:

%Insert Carnot Cycle/Piston diagrams

%Insert Carnot Cycle explaination diagram

For a fixed N: $dE=TdS-PdV$ so that 

$$\Delta Q_{1}=\int_{S_{1}}^{S_{2}}TdS=T_{h}(S_{2}-S_{1})$$
$$\Delta Q_{2}=\int_{S_{2}}^{S_{1}}TdS=-T_{c}(S_{2}-S_{1})$$

%Insert Carnot explaination diagram

Since after a full cycle the energy of the system goes back to what it was in the beginning we find:

$$\Delta Q_{1}+\Delta Q_{2}=\Delta W=$$ the area inside the cycles in the TxS plane.

In the expression above, $\Delta Q_{1}+\Delta Q_{2}$ is the net heat absorbed by the system and $\Delta W$ is the work done \textbf{by} the system.  

%Insert Carnot cycle deltaW diagram

For practical reasons, we would like to get as much work as possible from a given amount of heat absorbed.  The efficiency of the engine is defined as"

$$\epsilon=\frac{\Delta W}{\Delta Q_{1}}=\frac{(T_{h}-T_{c})(S_{2}-S_{1})}{T_{h}(S_{2}-S_{1})}=1-\frac{T_{c}}{T_{h}}<1$$

Now, one very important thing to note about this expression is that the total efficiency is \textit{less than} not \textit{less than or equal to} one because some of the energy is wasted and dispersed into the cold sink.  This begs the question: is there a way that we could change the cycle to obtain a larger efficiency? 

%Insert diagrams 

Going back to our TxS plane, we would want to maximize the area inside of the cycle ($\Delta W$) and minimize the area below the top curve ($\Delta Q_{1}$).  It's clear, for instance, that the cycle on the left can have its efficiency improved by changing it as shown below.  We end up then with the Carnot cycle, which is the best that can be done.  Thus, we have the universal result:

$$\epsilon\leq1-\frac{T_{c}}{T_{h}}$$  

Notice that we only used the fact that an ideal gas was used inside the engine to describe the different phases of the cycle in words.  The main result, $E\leq1-\frac{T_{c}}{T_{h}}$ does not depend on that and is valid universally.  

Finally, let us revisit the assumption of adiabaticity, namely, that all phases of the cycle were infinitely slow.  In each phase of the cycle the machine foes from the one initial state A to a final state B exchanging a certain amount of heat $\Delta Q$ (with a source at temperature T) and performing some work $\Delta W$ (positive or negative).  I'll assume the work is done on a reversible work source, namely, a system which can exchange work but not heat (a paper weight on top of a piston is an example).  Similarly, heat is exchanged with a reversible heat source.

%Insert heat exchange diagram

For given initial and final states A and B (and consequently, given changes in energy $E_{B}-E_{A}$ and entropy) the maximum amount of work will be delivered if $\Delta Q$ is a minimum as $E_{B}-E_{A}=-(\Delta W+\Delta Q)$.  The process with the smallest $\Delta Q$ will correspond to the smallest possible increase of the entropy of the universe (system+heat source, the work source has no entropy).  How much is the minimum amount of heat to be transferred to the heat source?  As the total entropy of the universe must increase (otherwise the system would stay at A) we have 

$$S_{B}-S_{A}+T_{r}\Delta Q\geq0$$

Where $T_{r}$ is the temperature of the heat source.  That means that the minimum amount of heat "wasted" on the heat source happens when the total entropy of the universe is unchanged as is given by 

$$\Delta Q=\frac{S_{A}-S_{B}}{T_{r}}$$

If more heat than $\frac{S_{A}-S_{B}}{T_{r}}$ is delivered to the heat source the universe's entropy will increase and less work will be delivered to the work source.  A similar analysis applies when the system is absorbing heat or work.  The main point is that \textbf{reversible processes}, those where the entropy of the universe does not change, lead to the maximum efficiency.  Thus, our bound on the efficiency of an engine obtained assuming adiabaticity, is not violated by non-adiabatic processes.  Notice that the adiabaticity condition means that the engine works really slowly.  Practical Machines, which need to do their job fast, are not adiabatic and much less efficient than the theoretical limit would expect (roughly a factor of 3 less efficient).

\section{Thermodynamic Stability}

Previously we determined the equilibrium condition between two subsystems: 

$$T_{1}=T_{2}$$ for a diathermal wall

$$P_{1}=P_{2}$$ for a moveable wall

$$\mu_{1}=\mu_{2}$$ for a porous wall

using the maximum entropy principle.  But we only imposed the $S=S_{1}+S_{2}$ was a critical point: it could have been a maximum, a minimum, or a saddle point.  We'll now see what the consequences are of demanding that the final equilibrium state is indeed a maximum.  

Let's consider first the case of two subsystems separated by a diathermal wall.  The condition that $S=S_{1}+S+{2}$ is a maximum implies:

%insert final state of equilibrium diagram

$$\dif{S}{E_{1}^{'}}=0\rightarrow\pardif{S_{1}}{E_{1}}=\pardif{S_{2}}{E_{2}}$$
$$\frac{d^{2}S}{dE_{1}^{2}}<0\rightarrow\frac{\delta^{2}S_{1}}{\delta E_{1}^{2}}+\frac{\delta^{2}S_{2}}{\delta E_{2}^{2}}<0$$

The quantity $\frac{\delta^{2}S}{\delta E^{2}}$ has a physical interpretation:

$$\secthermdif{S}{E}{V}{N}=\thermdif{\frac{1}{T}}{E}{V}{N}=-\frac{1}{T^{2}}\thermdif{T}{E}{V}{N}$$

$$=-\frac{1}{T^{2}}\frac{1}{\thermdif{E}{T}{V}{N}}=-\frac{1}{T^{2}N}\frac{1}{c_{v}}$$

Where $c_{v}$ is the system's specific heat at a constant volume.  From the above we see that the stability condition is then:

$$-\frac{1}{N_{1}T^{2}}\frac{1}{c_{v}^{1}}-\frac{1}{N_{2}T^{2}}\frac{1}{c_{v}^{2}}<0$$
or
$$\frac{1}{N_{1}c_{v}^{1}}+\frac{1}{N_{2}c_{v}^{2}}>0$$

This is true because $T_{1}=T_{2}$ in equilibrium.  By taking one subsystem to be much larger we conclude that each one of $c_{v}^{1}$ and $c_{v}^{2}$ should be positive.  This result is physically obvious.  Suppose $c_{v}^{1}<0$.  In this case, the system would get \textit{colder} by absorbing heat.  This would, of course, cause the system to absorb even more heat in a runaway cycle.  Clearly, this is an unstable situation.  

Next, let's consider the case where both E and V are exchanged between the subsystems.  We have:

$$0=\pardif{S}{E_{1}^{1}}=\pardif{S_{1}}{E_{1}}-\pardif{S_{2}}{E_{2}}$$
$$0=\pardif{S}{V_{1}^{1}}=\pardif{S_{1}}{V_{1}}-\pardif{S_{2}}{V_{2}}$$

$$0>d^{2}S=(dE_{1}^{1}dV_{1}^{1}) \left[ \begin{smallmatrix} \secpardif{S_{1}}{(E_{1}^{1})}&\frac{\delta^{2}S_{1}}{\delta E_{1}^{1}\delta V_{1}^{1}}\\ \frac{\delta^{2}S_{1}}{\delta E_{1}^{1}\delta V_{1}^{1}}&\secpardif{S_{1}}{(V_{1}^{1})} \end{smallmatrix} \right] 
\begin{bmatrix}
	dE_{1}^{1} \\
	dV_{1}^{1} 
\end{bmatrix}
+1\leftrightarrow2$$

We need that last term, 2, to be negative for any $dE_{1}^{1},dV_{1}^{1}$.  For the second variation $d^{2}S$ to be negative for any $E_{1}^{1},V_{1}^{1}$ it's necessary that the eigenvalues of 

$$\left[ \begin{smallmatrix} \secpardif{S_{1}}{(E_{1}^{2})}&\frac{\delta^{2}S_{1}}{\delta E_{1}^{1}\delta V_{1}^{1}}\\ \frac{\delta^{2}S_{1}}{\delta E_{1}^{1}\delta V_{1}^{1}}&\secpardif{S_{1}}{(V_{1}^{1})} \end{smallmatrix} \right]$$

to be negative.  So let's work out a general condition for this to be true.

$$\left| \begin{smallmatrix} a-\lambda&b\\ b&c-\lambda \end{smallmatrix} \right|=0$$

where

$$\lambda^{2}-(a+c)\lambda+ac-b^{2}=0$$
$$\lambda=\frac{a+c}{2}\pm\sqrt{\frac{a+c}{2}^{2}-(ac-b^{2})}<0$$
$$a+c<0 \& ac-b^{2}>0$$
$$a<0 \& c-\frac{b^{2}}{a}<0$$

In this case, $a=\secthermdif{S}{E}{V}{N}=\thermdif{\frac{1}{T}}{E}{V}{N}=-\frac{1}{T^{2}}\thermdif{T}{E}{V}{N}=-\frac{1}{T^{2}N}\frac{1}{c_{v}}<0$.  Note that the last part of the expression is the same condition as before, namely that $c_{v}>0$.

$$0>c-\frac{b^{2}}{a}=\secthermdif{S}{V}{E}{N}-(\frac{\delta^{2}S}{\delta E\delta V})^{2}\frac{1}{\secthermdif{S}{E}{V}{N}}$$

%Insert Diagram

Although this last expression is rather messy, we can rewrite in a more clear way with a transparent physical meaning.

$$\thermdif{\frac{P}{T}}{V}{T}{ }=\nabla(\frac{P}{T})*\mathbb{n}$$
$$\thermdif{\frac{P}{T}}{V}{T}{ }=\thermdif{\frac{P}{T}}{V}{E}{ }-\thermdif{T}{V}{E}{ }\frac{1}{\thermdif{T}{E}{V}{ }}\thermdif{\frac{P}{T}}{E}{V}{ }$$
$$\thermdif{\frac{P}{T}}{V}{T}{ }=\secpardif{S}{V}+T^{2}\frac{\delta^{2}S}{\delta E\delta V}\vert_{E}\frac{1}{-T^{2}\secthermdif{S}{E}{V}{ }}\frac{\delta^{2}S}{\delta V\delta E}$$  

So we find $\thermdif{\frac{P}{T}}{V}{T}{ }=\frac{1}{T}\thermdif{P}{V}{T}{ }<0$.  The $\frac{1}{T}\thermdif{P}{V}{T}{ }<0$ term can be rewritten to yield:

$$-\frac{1}{VT}[-\frac{1}{V}\thermdif{V}{P}{T}{ }]^{-1}=-\frac{1}{VT}\frac{1}{R_{T}}$$
	
We have defined a new quantity above, $R_{T}$.  The quantity is the \textbf{Isothermal Compressibility} of the gas.  The stability condition we defined earlier implies that the isothermal compressiblity is positive, however, this result is not immediately obvious.  A negative compressibility would mean that the system expands as the pressure increases.  This behaviour would increase the pressure of the other system leading to a runaway expansion of our observed system.

\part{Canonical Ensemble}

%Insert Canonical ensemble resevoir diagram

For many situations, it is much easier to treat a system as if it is in contact with an infinite heat reservoir.  This important property of this reservoir is that it is unaffected by any energy it absorbs due to its size.  In order to derive an ensemble to describe this system, let's start with the microcanonical ensemble:

$$\rho(q,p)~\Gamma(E-\mathcal{H}(q,p))$$

Where $q,p$ are the coordinates of the subsystem only, $\Gamma$ is the function which represents the number of microstates of the reservoir with energy $E-\mathcal{H}(q,p)$.  The terms within the gamma function represent the total energy of the reservoir (E) and the subsystem energy ($\mathcal{H}$).  Continuing our derivation:

$$\rho(q,p)~e^{\frac{1}{\bltz}*S_{R}(E-\mathcal{H}(q,p))}$$
Where $S_{R}$ is the reservoir's entropy.
$$\rho(q,p)\cong e^{\frac{1}{\bltz}S_{R}(E)-\frac{1}{\bltz}\mathcal{H}(q,p)\pardif{S_{R}}{E}}$$
In the expression above, we see that $\frac{1}{\bltz}S_{R}(E)$ is a constant independent of the coordinants $(q,p)$, and the $\pardif{S_{R}}{E}\rightarrow\frac{1}{T_{R}}=\frac{1}{T}$ term shows the equilibrium equivalence of the reservoir temperature ($T_{R}$) and the subsystem temperature ($T$).  Thus:

$$\rho(q,p)~e^{-\frac{\mathcal{H}(q,p)}{\bltz T}}$$

Therefore, we see that $\rho(q,p)$ is independent of the reservoir's properties.  This expression illustrates what is known as the \textbf{Canonical Ensemble}.

$$\rho(q,p)=\frac{e^{-\beta\mathcal{H}(q,p)}}{Z}$$
$$Z=\int dqdpe^{-\beta\mathcal{H}(q,p)}$$
$$\beta=\frac{1}{\bltz T}$$

\section{Equivalence with the Microcanonical Ensemble}

To see the equivalence of these two ensembles, let's start with the partition function:

$$Z=\int dqdpe^{-\beta\mathcal{H}(q,p)}=\int dE\Gamma(E)e^{-\beta E}$$

Where $dE\Gamma(E)$ is the number of microstates with energy between $E$ and $E+dE$.  Therefore

$$Z=\int dE e^{\frac{S(E)}{\bltz}-\beta E}$$

As we see in the above expression, $S(E)$ is the microcanonical ensemble and the entire exponent of $e$ is an extensive quantity which is $~N$ \textit{very} concentrated around the maximum of the integrand.

%Insert integral diagram

$$Z\cong\Delta e^{-\beta(E-TS(E))}\vert_{E=E*}$$

The energy constraint insists that the energy minimizes the value $E-TS(E)$.  This value has a particular significance and therefore has a specific name...

\section{Helmholtz Free Energy and the Partition Function}

As stated before, the quantity $E-TS(E,V,N)\vert_{E=E*}$ has a specific name: \textbf{Helmholtz Free Energy}, and is denoted by the function $F(T,V,N)$.  Before we further discuss the importance of this quantity, let's first explore the significance of the normalization factor Z in the Canonical Ensemble.  Z also has a specific name, the \textbf{Partition Function}.  It is extremely useful in its own right:

$$\pardif{\ln{Z}}{\beta}=\frac{1}{T}\pardif{Z}{\beta}=\frac{1}{Z}\int dqdp(-\mathcal{H})e^{-\beta\mathcal{H}(q,p)}=<-E>$$

Where $<-E>$ is the average of the energy ion the canonical ensemble.

$$\secpardif{\ln{Z}}{\beta}=\frac{\delta}{\delta\beta}\frac{1}{Z}\pardif{Z}{\beta}=-\frac{1}{Z^{2}}(\pardif{Z}{\beta})^{2}+\frac{1}{Z}\secpardif{Z}{\beta}$$

$$\secpardif{\ln{Z}}{\beta}=-<E>^{2}+\frac{1}{Z}\int dqdp\mathcal{H}^{2}(q,p)e^{-\beta\mathcal{H}(q,p)}$$

$$\secpardif{\ln{Z}}{\beta}=-<E>^{2}+<E^{2}>$$

$$\secpardif{\ln{Z}}{\beta}=<(E-<E>)^{2}>$$

However, keep in mind that 

$$\frac{\sqrt{<(E-<E>)^{2}}}{<E>}~\frac{\sqrt{N}}{N}~\frac{1}{\sqrt{N}}\rightarrow_{N>>1}0$$

Where $\sqrt{<(E-<E>)^{2}}$ is the fluctuation in energy and the $<E>$ is the average energy.

This shows that in the thermodynamic limit ($N\rightarrow\infty$) the energy does not fluctuate and only the microstates with an energy equal to the average value contributes significantly when computing averages using the canonical ensemble.  In other words, the canonical ensemble is equivalent to the microcanonical ensemble.

%Insert E/\rho diagram 

There are two specific situations where the canonical ensemble is used:

\begin{itemize}
	\item The system is coupled to a heat reservoir
	\item The system is isolated but it is easier to compute averages in the canonical ensemble.  This ends up being a very common situation
\end{itemize}

Previously we have argued that the averages of observables computed in the canonical ensemble are the same (in the Thermodynamic Limit) as the microcanonical ensemble.  Now, we would like to demonstrate that the information contained in $S(E,V,N)$ (or, equivalently, E(S,V,N)) is also contained in the Helmholtz Free Energy, $F(T,V,N)$(notice the arguments T,V,N; \textbf{NOT} S,V,N!).  To do this, we'll have to utilize a mathematical algorithm known as the Legendre Transform.

\subsection{Legendre Transforms}

Suppose $f(x)$ is not known but $f(p)$, where $p=\dif{f}{x}$, \textit{is} known.  Can we recover the function $f(x)$ from $f(p)$?  \textbf{NO!}  This is clearly illustrated by the figure below:

%Insert L-T graph proof

As an analytical example, let's take a look at a couple of functions and see if we can perform this recovery:

$$f(x)=x^{2}, p=\dif{f}{x}=2x \rightarrow f(p)=(\frac{p}{2})^{2}$$
$$f(x)=(x-a)^{2}, p=\dif{f}{x}-2(x-a) \rightarrow f(p)=(\frac{p}{2})^{2}$$

From two different starting functions we extracted two different values of p, however, $f(p)$ in both cases is the same.  Clearly we cannot recover a unique initial function from the information contained solely in $f(p)$.  

Now, for single variable functions this constant shift is not a big deal.  But for functions of several variables we end up seeing shifts and alterations that are functions of other variables, which greatly complicates the situation.  

Even though we cannot recover the initial function from the $f(p)$ values shown above, is it possible to construct some $f(p)$ where it \textit{is} possible to reconstruct the initial function?  It turns out the answer is yes!  A function of $p=\dif{f}{x}$ containing all of the information needed to reconstruct $f(x)$ is the intercept of the tangent with the vertical axis:

%Insert L-T Intercept diagrams

$\rightarrow p=\frac{f-\tilde{f}}{x}\rightarrow\tilde{f}(p)=px-f(x)\vert_{x such that \frac{f}{x}=p }$

Returning to thermodynamics, we see from its definition that $-F(T,V,N)$ is the Legendre Transform of $E(S,V,N)$ in relation to the variable $S$.  Indeed, the Legendre Transform of $E$ is: 

$$E(S,V,N)\rightarrow\tilde{E}(T,V,N)=TS-E(S,V,N)\vert_{S=S* such that \pardif{E}{S}=T}$$
$$T=\pardif{E}{S}$$

From the definition of the Helmholtz Free Energy:

$$F(T,V,N)=E-TS\vert_{E=E*, such that 0=\frac{\delta}{\delta E}(E-TS)=1-T\pardif{S}{E}\vert_{E=E*}}$$

But the value $S*$ for which $\thermdif{E}{S}{S*}{ }=T$ corresponds to $E^{*}=E(S^{*},V,N)$, so

$$F(T,V,N)=-\tilde{E}(E,V,N)$$

Notice that:

$$\thermdif{F}{T}{V}{N}=\frac{\delta}{\delta t}[(E(S,V,N)-TS)\vert_{S such that \thermdif{E}{S}{S*}{ }=T}]\vert_{V,N}$$
$$\thermdif{F}{T}{V}{N}=\thermdif{E}{S}{V}{N}\thermdif{S*}{T}{V}{N}-S^{*}-T\thermdif{S^{*}}{T}{V}{N}$$

(Remember that $\thermdif{E}{S}{V}{N}=T$!)

$$\thermdif{F}{T}{V}{N}=-S^{*}$$

For expedience, from now on we will be dropping the * superscript.

$$\thermdif{F}{V}{T}{N}=\frac{\delta}{\delta V}[(E(S,V,N)-TS)\vert_{S such that \thermdif{E}{S}{S^{*}}{ }}]\vert_{T,N}$$
$$=\thermdif{E}{V}{S}{N}+\thermdif{E}{S}{V}{N}\pardif{S^{*}}{V}-T\pardif{S^{*}}{V}=\thermdif{E}{V}{S}{N}=-P$$
$$\thermdif{F}{N}{T}{V}=\thermdif{E}{N}{S}{V}+\thermdif{E}{S}{V}{N}\pardif{S^{*}}{N}-T\pardif{S^{*}}{N}=\thermdif{E}{N}{V}{N}=\mu$$

Thus, $dF=-SdT-PdV+\mu dN$  Also, as a consequence of Euler's relations, $E=TS-PV+\mu N\rightarrow F=-PV+\mu N$

\subsubsection{EXAMPLE: Classical, Ideal, Non-Relativistic, Monoatomic Gas}

$$\mathcal{H}(\mathbb{r}_{i},\mathbb{p}_{i})=\sumser{i=1}{N}\frac{\mathbb{p}_{i}^{2}}{2m}$$
$$Z=\int\prod_{i=1}^{N}\frac{d^{3}\mathbb{r}_{i}d^{3}\mathbb{p}_{i}}{N!}e^{-\beta\sum_{i=1}^{N}\frac{\mathbb{p}_{i}^{2}}{2m}}$$
$$=[\int d^{3}\mathbb{r}\int d^{3}\mathbb{p}e^{-\beta\frac{\mathbb{p}^{2}}{2m}}]^{3}\frac{1}{N!}$$

(Remember that $\int d^{3}\mathbb{r}$ is simply $V$)
 
$$=V^{N}\left[\int_{-\infty}^{\infty}dpe^{-\frac{p^{2}}{2m\bltz T}}\right]^{3N}\frac{1}{N!}$$

$$=V^{N}(\frac{2m\bltz T}{\pi})^{\frac{3N}{2}}\frac{1}{N!}$$

Therefore, in the $N>>1$ limit,

$$F\cong-\bltz T\ln{Z}=-N\bltz T[\ln(V(\frac{2m\bltz T}{\pi})^{\frac{3}{2}})-\ln{N}]$$

$$F\cong-N\bltz T\ln{\frac{V}{N}(\frac{2m\bltz T}{\pi})^{\frac{3}{2}}}$$

This shows that $F$ is extensive, which it should be.  Now, we take the Legendre Transform of the Helmholtz Free Energy to find $E(S,V,N)$ from $F(T,V,N)$:     

$$\pardif{F}{T}=-S=-N\bltz\ln{\frac{V}{N}(\frac{2m\bltz T}{\pi})^{\frac{3}{2}}}\frac{-3N\bltz}{2}$$
$$\rightarrow\bltz T=\frac{\pi}{2m}\frac{N^{\frac{2}{3}}}{V^{\frac{2}{3}}}e^{\frac{2}{3}\frac{S}{N\bltz}-1}$$
$$E=F+TS=-N\bltz T\ln{\frac{V}{N}(\frac{2m\bltz T}{\pi})^{\frac{3}{2}}}+\frac{\pi}{2m}(\frac{N}{V})^{\frac{2}{3}}e^{\frac{2S}{3N\bltz}-1}S$$
$$E=\frac{3}{2}\frac{\pi}{2m}\frac{N^{\frac{5}{3}}}{V^{\frac{2}{3}}}e^{\frac{2S}{3N\bltz}-1}$$

Which is the same value that we calculated before using the microcanonical ensemble.

\subsubsection{EXAMPLE: Classical model of a Paramagnet}

%Insert magnet model diagram

$$\mathcal{H}(\theta_{i},L_{i})=\sumser{i=1}{N}\frac{L_{i}^{2}}{2I}-H\mu\cos{\theta_{i}}$$

$$Z=\int\prod_{i=1}^{N}d\theta_{i}dL_{i}e^{-\beta\sumser{i=1}{N}[\frac{L_{i}^{2}}{2I}-H\mu\cos{\theta_{i}}]}$$
$$=[\int_{0}^{2\pi}d\theta e^{\beta H\mu\cos{\theta}}]^{N}[\int_{-\infty}^{\infty}dLe^{-\frac{\beta L^{2}}{2I}}]^{N}$$

Here we see that $[\int_{0}^{2\pi}d\theta e^{\beta H\mu\cos{\theta}}]^{N}$ is simply the $I_{0}$ Bessel Function, specifically $I_{0}(\beta H\mu)$.  The second term, $[\int_{-\infty}^{\infty}dLe^{-\frac{\beta L^{2}}{2I}}]^{N}$, is approximately the value of $\sqrt{T}$.  Therefore,

$$Z~[I_{0}(\beta H\mu)\sqrt{T}]^{N}$$

Thus we can easily express the Helmholtz Free Energy of the system as 

$$F=-\bltz T\ln{Z}=-\bltz T\ln{\sqrt{T}I_{0}(\beta H\mu)}$$

Now, let's derive the value of an observable.  For example, the average magnetic moment per site in our system: $<M>$

$$\frac{1}{N}<M>=<\frac{1}{N}\sumser{i=1}{N}\mu\cos{\theta_{i}}>=\frac{1}{Z}\int d\theta_{i}dL_{i}\frac{\mu}{N}\sum\cos{\theta_{i}}e^{-\beta\mathcal{H}}$$
$$=\frac{1}{\beta N}\frac{\delta}{\delta H}\ln{Z}=-\pardif{F}{H}$$
$$=\frac{I_{1}(\frac{H\mu}{\bltz T})}{I_{0}(\frac{H\mu}{\bltz T})}$$

The above expression yields two different values based on the limit of $H\mu$.  

$$H\mu<<\bltz T\rightarrow \frac{H\mu T}{2}$$
$$H\mu>>\bltz T\rightarrow sign(H)$$

%Insert <M> graph diagram

\subsubsection{EXAMPLE: Classical Harmonic Oscillators and the Equipartition Theorem}

Systems of linked harmonic oscillators are extremely important in modelling many physical phenomena.  Let's examine such a system using the equipartition theorem:

In this case, the Hamiltonian that describes the system is: $\mathcal{q_{i},p_{i}}=\sumser{i=1}{N}[\frac{p_{i}^{2}}{2m}+\frac{m\omega^{2}}{2}q_{i}^{2}]$

Therefore, we can use this to write the system's partition function:

$$Z=\int\prod_{i=1}^{N}dq_{i}dp_{i}e^{-\beta\sumser{i=1}{N}[\frac{p_{i}^{2}}{2m}+\frac{m\omega^{2}}{2}q_{i}^{2}]}$$
$$=[\int_{-\infty}^{\infty}dqe^{-\frac{\beta m\omega^{2}}{2}q^{2}}]^{N}[\int_{-\infty}^{\infty}dpe^{-\frac{\beta p^{2}}{2m}}]^{N}$$

The first term is approximately $\frac{1}{\sqrt{m\omega^{2}\bltz T}}$ and the second term is approximately $\sqrt{2m\bltz T}$.

Many of these variables cancel or reduce, leaving us with

$$Z~(\frac{\bltz T}{\omega})^{N}=(\beta\omega)^{-N}$$

Solving for the average energy yields:

$$<E>=-\frac{1}{Z}\pardif{Z}{\beta}=-\pardif{\ln{Z}}{\beta}=\frac{N}{\beta}=N\bltz T$$

This result is important, as it shows in our classical model that the average energy of any one harmonic oscillator at temperature T is $\bltz T$, entirely independent of $\omega$.  Oscillators with a large $\omega$ (stiff spring constant) will oscillate less, the oscillators with a small $\omega$ (loose spring constant) will oscillate more, but all will have the same average energy:

%Insert oscillator energy graphs.

If we had worked this example for free particles in one dimension instead of oscillators we would have found that $<E>=\frac{1}{2}\bltz T$ (see appendix).  The fact that many physical systems can be written in terms of coupled harmonic oscillators (normal modes) implies the \textbf{Equipartition Theorem}.  This theorem is extremely powerful and is very useful in solving classical systems, so it is stated separately below:

\textbf{Equipartition Theorem}: Each quadratic term in the Hamiltonian contributes $\frac{\bltz T}{2}$ to the system's average energy ($<E>$).

\subsubsection{EXAMPLE: Classical, Non-relativistic, Diatomic Gas}

Assuming the vibrational mode can be treated as classical, we have:

$$\mathcal{H}=\sumser{i=1}{N}\frac{\mathbb{p_{i}^{2}}}{2m}+\frac{p_{i}^{2}}{2m}+\frac{m\omega^{2}}{2}x_{i}^{2}$$

Disecting this hamiltonian, we see that the first term is a momentum vector.  This represents the center of mass momentum, so it contributes 3 quatratic terms per molecule.  The second and third terms together represent the molecule's vibrational energy, where $x_{i}$ is the distance between the atoms in the molecule.  The vibrational energy contributes 2 quadractic terms per molecule.  Since we know this information, it is now trivial to construct an expression for the average energy using the Equipartition Theorem:

$$\frac{1}{N}<E>=(3+2)\frac{1}{2}\bltz T=\frac{5}{2}\bltz T$$

\section{Quantum Statistical Mechanics}

So far we have been assuming that the systems we have been working with behaved Classically at the microscopic level.  Now let's rework our previous formalism to describe systems which exist within the realm of Quantum Mechanics.  There are two ways to accomplish this:  We can trust analogies and guess the form of most expressions or we can generalize the concept of the ensemble for quantum systems.  The second option is much more powerful, however, the first is quicker and will be enough for the majority of the systems we will be dealing with.

\begin{tabu} to 0.8\textwidth { | X[r] | X[c] | X[c] | }
	\hline
	 & Classical & Quantum \\
	\hline
	Observables: & Function in Phase Space $\rightarrow f(q,p)$ & Operator in Hilbert Space $\rightarrow\hat{f}$ \\
	\hline
	Averages in the Canonical Ensemble: & $\frac{1}{Z}\int dqdpe^{-\beta H}f(q,p), Z=\int dqdpe^{-\beta H}$ & $\frac{1}{Z}Tr[e^{-\beta\hat{H}}\hat{f}], Z=Tr[e^{-\beta\hat{H}}]$\\
	\hline
	Partition Function: & $Z=\int dE\Gamma(E)e^{-\beta E}$ & $Z=\sum_{n}<n|e^{-\beta\hat{H}}|n>=\sum_{n}e^{-\beta E_{n}}\cong\int dE\Gamma(E)e^{-\beta E}$\\
	\hline
\end{tabu}  

A couple of notes regarding the quantum partition function: $n$ is are the eigenstates of $\hat{H}$ and $E_{n}$ are the eigenvalues.  $dE\Gamma(E)$ is the number of quantum states between E and E+dE.  The integral over these states approximates the partition function only for very dense spectra.  

\subsection{The Density Matrix}

Although this concept won't actually be used in this course, it is a basic thing that any good physicist worth their salt should know.  Suppose we know the quantum state of a system at a given time, for example

$$|\Psi>$$

We can compute the expectation value of an observable (represented by operator $\hat{A}$ in this case) by

$$<A>=<\Psi|\hat{A}|\Psi>$$

Where $<A>$ is the average over results obtained by many identically prepared systems.  An interesting observation of this statement is that we can get different results from identical measurements on identically prepared systems.  This uncertainty is due to the nature of Quantum Mechanics.  In particular, if a system is prepared in the state $|\Psi>=\alpha|\Psi_{1}>+\beta|\Psi_{2}>$, $<A>$ is given by:

$$<A>=<\Psi|\hat{A}|\Psi>=|\alpha|^{2}<\Psi_{1}|\hat{A}|\Psi_{1}>+|\beta|^{2}<\Psi_{2}|\hat{A}|\Psi_{2}>+\alpha^{*}\beta<\Psi_{1}|\hat{A}|\Psi_{2}>+\alpha\beta^{*}<\Psi_{2}|\hat{A}|\Psi_{1}>$$

$\alpha^{*}\beta<\Psi_{1}|\hat{A}|\Psi_{2}>+\alpha\beta^{*}<\Psi_{2}|\hat{A}|\Psi_{1}>$ are the interference terms, so they may have either sign.  Now, suppose that eh state $|\Psi>$ is not known for sure, not because of any uncertainty imposed by quantum mechanics but just because of practical considerations.  Let us state that there is a probability $p_{1}$ that the system is in state $|\Psi_{1}>$ and probability $p_{2}$ that it is in state $|\Psi_{2}>$.  The actual wavefunction is either $|\Psi_{1}>$ or $|\Psi_{2}>$.  Therefore, we expect that the expectation value of observable $A$ is:

$$<A>=p_{1}<\Psi_{1}|\hat{A}|\Psi_{1}>+p_{2}<\Psi_{2}|\hat{A}|\Psi_{2}>=Tr[(p_{1}|\Psi_{1}><\Psi_{1}|+p_{2}|\Psi_{2}><\Psi_{2}|)\hat{A}]=Tr[\hat{\rho}\hat{A}]$$

The operator $\hat{\rho}$ is called the \textbf{Density Matrix} and it is the quantum analogue of the phase space distribution $\rho(q_{i},p_{i})$.  Of course, this definition generalizes to the case where the classical uncertainty extends to several states:

$$\hat{\rho}=\sum_{n}p_{n}|n><n|$$

From this definition a few properties directly follow:

\begin{itemize}
	\item $Tr\hat{\rho}$ (since $\sum_{n}p_{n}=1$)
	\item $\hat{\rho}^{\dagger}=\hat{\rho}$ ($\rho$ is hermitian, since the analogue of $\rho(q_{i},p_{i})$ must be real)
\end{itemize}

If $\hat{\rho}=|\Psi><\Psi|$, that is, if we know the quantum state of a system with 100\% confidence, we say that the system is in a \textit{pure} state.  Otherwise we say the system is in a \textit{mixed} state.  For a pure state $\hat{\rho}^{2}=|\Psi><\Psi|\Psi><\Psi|=|\Psi><\Psi|=\hat{\rho}$, so $\hat{\rho}^{2}=\hat{\rho} \leftrightarrow \hat{\rho}$ is a pure state.

Notice the difference between $|\Psi>=\frac{1}{\sqrt{2}}|\Psi_{1}>+\frac{1}{\sqrt{2}}|\Psi_{2}>$ and $\hat{\rho}=\frac{1}{2}|\Psi_{1}><\Psi_{1}|+\frac{1}{2}|\Psi_{2}><\Psi_{2}|$.  Both could be described as a 50-50 combination of $|\Psi_{1}>$ and $|\Psi_{2}>$ but they represent completely different physical situations.  

\subsection{EXAMPLE: Quantum Harmonic Oscillator}

We begin, as always, with the Hamiltonian.

$$H=\frac{\hat{p}^{2}}{2m}+\frac{m\omega^{2}}{2}\hat{x}^{2}\rightarrow E_{n}=(n+\frac{1}{2})\hbar\omega,|n>, n=0,1,...$$

$$Z=\sum_{n_{1}=0}^{\infty}\sum_{n_{2}=0}^{\infty}\sum_{n_{3}=0}^{\infty}....\sum_{n_{N}=0}^{\infty}e^{\beta\hbar\omega[(n_{1}+\frac{1}{2})+(n_{2}+\frac{1}{2})+....+(n_{N}+\frac{1}{2})]}$$
$$Z=[\sum_{n=0}^{\infty}e^{-\beta\hbar\omega(n+\frac{1}{2})}]^{N}$$
$$Z=[e^{-\frac{\beta\hbar\omega}{2}}\frac{1}{1-e^{-\beta\hbar\omega}}]^{N}$$
$$Z=[\frac{1}{e^{\frac{\beta\hbar\omega}{2}}-e^{-\frac{\beta\hbar\omega}{2}}}]^{N}$$
$$Z=[\frac{1}{2\sinh{\frac{\beta\hbar\omega}{2}}}]^{N}$$

$$\frac{1}{N}<E>=-\frac{1}{ZN}\pardif{Z}{\beta}=-\frac{1}{N}\pardif{\ln{Z}}{\beta}=\frac{1}{N}\pardif{N\ln{2\sinh{\frac{\beta\hbar\omega}{2}}}}{\beta}$$
$$=\frac{\cosh{\frac{\beta\hbar\omega}{2}}}{\sinh{\frac{\beta\hbar\omega}{2}}}=\frac{\hbar\omega}{2}\coth{\frac{\beta\hbar\omega}{2}}$$

Similar to the classical situation, this expression again yields two different solutions based on the limit of $\bltz T$.  

$$\bltz T>>\hbar\omega \rightarrow \frac{\hbar\omega}{2}+\bltz T$$
$$\bltz T<<\hbar\omega \rightarrow \frac{\hbar\omega}{2}+\hbar\omega e^{-\frac{\hbar\omega}{\bltz T}}$$

Note that the first result is the classical result $\bltz T$ modified by the addition of zero-point energy ($\frac{\hbar\omega}{2}$).  The second result includes the same zero-point energy, but also includes "Boltzmann Supressed" thermal energy.  This arises because in the $\bltz T<<\hbar\omega$ limit the high frequency modes are "frozen".  This runs contrary to the classical equipartition theory.

\subsection{EXAMPLE: Specific Heat of Molecular Gases}

Molecules have vibrational and rotational degrees of freedom in addition to the standard translational ones.  For instance, for $H_{2}$:

$$\mathcal{H}=\sumser{i=1}{N}\frac{\mathbb{P}_{i}^{2}}{2M}+\frac{\mathbb{P}_{i}^{2}}{2m_{red}}+m_{red}\frac{\omega^{2}}{2}x_{i}^{2}+\frac{\mathbb{L}_{i}^{2}}{2I}$$

Breaking this down by terms, we have: $\frac{\mathbb{P}_{i}^{2}}{2M}$, which is the center of mass kinetic energy ($M=2m$); $\frac{\mathbb{P}_{i}^{2}}{2m_{red}}$, which is the vibrational kinetic energy ($\frac{1}{m_{red}}=\frac{1}{m}+\frac{1}{m}$); $m_{red}\frac{\omega^{2}}{2}x_{i}^{2}$, which is the vibrational potential energy ($x_{i}$ is the separation between atoms); and $\frac{\mathbb{L}_{i}^{2}}{2I}$ which is the rotational energy.

Classically, we express the partition function thusly:

$Z=\int\prod_{i=1}^{N}d^{3}\mathbb{R}_{i}d^{3}\mathbb{P}_{i}dp_{i}dx_{i}d^{2}\mathbb{L}_{i}e^{-\beta\mathcal{H}}$

Note that there is no rotation around the molecule's axis.  This is because such rotation would result in a state that is indistinguishable from the initial state and therefore have no effect on our system.

If all of these degrees of freedom are classical,

%Where are the rest of the notes?

\part{Quantum Gases}

Previously, using the grand canonical ensemble we found the grand canonical potential of a quantum ideal gas of identical particles:

$$\Omega(T,V,\mu)=\pm g\bltz TV\int\frac{d^{3}k}{(2\pi)^{3}}\ln(1\mp e^{-\beta(\frac{\hbar^{2}k^{2}}{2m}-\mu)})$$

In this expression, we nicely sum up the main components of the canonical potential for both bosonic and fermionic gases.  To clarify, bosons behave according to the upper operators and fermions utilize the lower operators (i.e. for bosons there is a plus sign outside the log and a minus sign inside, and the converse is true for fermions).  We also rely on $g$, which is the degeneracy of single particle (orbital) states, $2s+1$ for particles of spin s.  Within the log, the value $\frac{\hbar^{2}k^{2}}{2m}$ is the energy $E_{k}$.  This is the energy of one particle with momentum of $\hbar k$.  

Now, $N=-\pardif{\Omega}{\mu}$ since $-\pardif{\Omega}{\mu}=\frac{\delta}{\delta \mu}\frac{1}{\beta}\ln{Tr[e^{-\beta(\hat{H}-\mu\hat{N})}]}$

$$=\frac{1}{\beta}\frac{Tr[\beta e^{-\beta(\hat{H}-\mu\hat{N})}\hat{N}]}{Tr[e^{-\beta(\hat{H}-\mu\hat{N})}]}=Tr[\hat{\rho}\hat{N}]=<N>$$
 
Using this information we can find the particle density:

$$\frac{N}{V}=-\pardif{\Omega}{\mu}=\mp\bltz Tg\int \frac{d^{3}k}{(2\pi^{3})}\frac{\mp\beta e^{-\beta(E_{k}-\mu)}}{1\mp e^{-\beta(E_{k}-\mu)}}$$

$$=g\int\frac{d^{3}k}{(2\pi^{3})}\frac{1}{e^{\beta(E_{k}-\mu)}\mp1}$$

Similarly,

$$S=-\pardif{\Omega}{T}=\mp\bltz gV\int\dtk\ln{1\mp\ebk}+\bltz gV\int\dtk\frac{\frac{1}{\bltz T^{2}}\ebk\epsilon_{k}}{1\mp\ebk}$$

Where we have defined a new variable $\epsilon_{k}=E_{k}-\mu$.  Continuing:

$$=\mp\bltz gV\int\dtk[\ln{1\mp\ebk}\mp n(k)\ln{e^{\beta\epsilon_{k}}}]$$

The terms we are integrating over are equivalent to $\frac{1}{1\pm n(k)}$ and $\frac{1\pm n(k)}{n(k)}$.  Therefore,

$$=\mp gV\int\dtk[-\ln{1\pm n(k)}\mp n(k)\ln{1\pm n(k)}\pm n(k)\ln{n(k)}]$$
$$=\bltz gV\int\dtk[(n(k)\pm 1)\ln{1\pm n(k)}-n(k)\ln{n(k)}]$$
$$=\pm gV\int\dtk[(1\pm n(k))\ln{1\pm n(k)}\mp n(k)\ln{n(k)}]$$

We know $\Omega=E-TS-\mu N$, and, by the Euler relation, $E=TS-PV+\mu N$, so $\Omega=-PV$.  That means that 

$$E=\Omega+TS+\mu N=\pm V\int\dtk[\ln{1\mp\ebk}+(1+n(k))\ln{1\pm n(k)}\mp n(k)\ln{n(k)}\pm\frac{\mu}{\bltz T}n(k)]$$
$$=\pm\bltz TV\int\dtk[\ln{\frac{1}{1\pm n(k)}}+(1\pm n(k))\ln(1\pm n(k))\mp n(k)\ln{n(k)}]$$

Notice that the expression $(1\pm n(k))\ln(1\pm n(k))\mp n(k)\ln{n(k)}=\pm\ln{\frac{1\pm n(k)}{n(k)}}\pm\mu\beta n(k)$.

Therefore,

$$\pm n(k)\ln{e^{\beta\epsilon_{k}}\pm1\mp1}=\pm n(k)\beta\epsilon_{k}$$
$$=\pm\bltz TV\int\dtk\pm n(k)\beta(\epsilon_{k}+\beta\mu)$$
$$=V\int\dtk n(k)E_{k}$$

As expected, the energy density is given by the sum over all momenta of the energy of each particle with that momentum multiplied by how many particle have that value of momentum.  

\section{Classical Limit}

Let's consider the regime where $e^{\beta\mu}<<1$ $(-\mu<<\bltz)$.  This implies

$$n(k)=\frac{1}{e^{\beta(E_{k}-\mu)}\mp 1}\cong e^{\beta\mu}e^{-\beta E_{k}}<<1$$

Which means our system has a low occupation number.  One thing to note is that $n(k)$ is simply $e^{\beta\mu}$ multiplied by the Maxwell-Boltzmann Distribution.  Continuing our analysis, let's calculate the expected particle density:

$$\frac{N}{V}=g\int\dtk n(k)\cong gVe^{\beta\mu}\int_{0}^{\infty}dk\frac{4\pi k^{2}}{8\pi^{3}}e^{-\frac{\beta\hbar^{2}k^{2}}{2m}}$$

Where $\int_{0}^{\infty}dk\frac{4\pi k^{2}}{8\pi^{3}}e^{-\frac{\beta\hbar^{2}k^{2}}{2m}}=\frac{1}{8\pi^{3}(\frac{2\pi m\bltz T}{\hbar^{2}})^{\frac{3}{2}}}$.  To simplify further, we recall the the Thermal De Broglie Wavelength is $\lambda_{T}=\sqrt{\frac{\hbar^{2}}{2\pi m\bltz T}}$, meaning we can write our particle density expression in the following form:
	
$$\frac{N}{V}=ge^{\beta\mu}\frac{1}{\lambda_{T}^{3}}$$

The Helmholtz Free Energy is given by:

$$F=E-TS=\Omega+\mu N\cong\pm\bltz T gV\int\dtk\ln{1\mp\ebk}+\mu N$$

Let's stop here for a moment and look at the above expression.  We know that $\ln{1\mp\ebk}\cong\mp e^{\beta\mu}e^{-\beta E_{k}}$ in the limit $-\mu>>T$  This means that in this limit we can approximate the first term, $\pm\bltz T gV\int\dtk\ln{1\mp\ebk}$, as $TN$.  In our classical case this is much less than the value of the second term $\mu N$.  

Continuing our analysis, we find:

$$F\cong N\ln{\frac{(2\pi)^{3}\lambda_{T}^{3}\mu}{gV}}\frac{1}{\beta}=N\mu$$

We obtained this value by solving for N in our previous derivation of the particle density.  Therefore, we find:

$$F\cong-\bltz TN[\ln{\frac{VT^{\frac{3}{2}}}{N}+Const.}]$$

Which is the Classical Result for Helmholtz Free Energy.  Note that in the Classical calculation the $\ln{\frac{1}{N}}$ factor comes from the $\frac{1}{N!}$ "fudge factor" we used previously.

In summary, we have derived that the conditions for classicality can be rephrased as:

$$e^{\beta\mu}<<1\leftrightarrow\frac{N}{V}\lambda_{T}^{3}<<1\leftrightarrow(\frac{\lambda_{T}}{d})^{3}<<1$$

Where $d$ is the interparticle distance.  All of these things are true for \textit{Dilute Gases}, so we come to the conclusion that \textbf{Dilute Gases are Classical.}

The expressions for S,F,... are similar for bosons and fermions.  The only difference is just a few signs but those signs can completely change the physics so it behooves us to look at these two cases separately.

\section{Ideal Fermi Gas}

For Fermi Gases, the occupation number is defined as follows:

$$n(k)=\frac{1}{e^{\beta(E_{k}-\mu)}+1}$$

%Insert Fermi Momentum Diagram

The value of $k$ where $n$ goes from $\cong 1$ to $\cong 0$:

$$E_{k}-\mu=0\rightarrow k=\frac{\sqrt{2m\mu}}{\hbar}\equiv k_{f}$$

$k_{f}$ is a special term which we will call the \textbf{Fermi Momentum}.

Examining this term, let's find the range of $k$ which $n$ goes from $\cong1$ to $\cong0$:

$$E_{k}-\mu\cong\bltz T\rightarrow\frac{\hbar^{2}}{2m}(k^{2}-k_{f}^{2})\cong\bltz T\rightarrow k-k_{f}\cong\frac{m}{\hbar k_{f}}\frac{\bltz T}{\hbar}$$

We are able to make this leap because $(k^{2}-k_{f}^{2})=(k+k_{f})(k-k_{f})\cong2k_{f}(k-k_{f})$ and $k\cong k_{f}$.  The term $\frac{\hbar k_{f}}{m}$ is another special quantity known as the \textbf{Fermi Velocity}.

Defined in words, the Fermi Momentum and the Fermi Velocity are the momentum and group velocity of a fermion at the Fermi surface.  

We know that the dilute ($e^{\beta\mu}<<1$) regime is classical.  So now we'll concentrate on the opposite regime, i.e. a degenerate Fermi gas.

\subsubsection{The Fermi Surface and the Sommerfeld Expansion}

%Insert Fermi Sphere diagram

In order to extract the physics from the formulae above we will need integrals of the form

$$\int_{0}^{\infty}dkk^{2\alpha}n(k)$$

in the low temperature (high $\beta$) limit.  The low T expansion is a little tricky and is known as the \textbf{Sommerfeld Expansion}.  

$$I(\alpha)=\int_{0}^{\infty}dkk^{2\alpha}n(k)=\int_{0}^{\infty}dE\frac{m}{\hbar^{2}}k^{2\alpha-1}n(E)$$

Noting that here $k^{2\alpha-1}=(\frac{2mE}{\hbar^{2}})^{2\alpha-1}$ and $E=\frac{\hbar^{2}k^{2}}{2m}$, $dE=\frac{\hbar^{2}k}{2m}dk$.

Therefore,

$$I(\alpha)=(\frac{2m}{\hbar^{2}})^{\alpha-\frac{1}{2}}\frac{m}{\hbar^{2}}\int_{0}^{\infty}dEE^{\alpha-\frac{1}{2}}n(E)$$

Where $\frac{1}{e^{\beta(E-\mu)+1}}$.

The \underline{difference} between $n(E)$ and its zero temperature limit is concentrated around the Fermi Surface.  The following steps make use of this observation.

%Insert T=0 graph

$$\int_{0}^{\infty}dEE^{\alpha-\frac{1}{2}}n(E)=n(E)\frac{E^{\alpha+\frac{1}{2}}}{\alpha+\frac{1}{2}}\vert_{0}^{\infty}-\int_{0}^{\infty}dE\frac{E^{\alpha+\frac{1}{2}}}{\alpha+\frac{1}{2}}\frac{dn(E)}{dE}$$

$$=\beta\int_{0}^{\infty}\frac{dE}{\alpha+\frac{1}{2}}[\mu^{\alpha+\frac{1}{2}}+\frac{(E-\mu)^{\alpha+\frac{1}{2}}}{\beta}\mu^{\alpha-\frac{1}{2}}+\frac{(E-\mu)^{2}\mu^{\alpha-\frac{3}{2}}}{2\beta^{2}}(\alpha+\frac{1}{2})(\alpha-\frac{1}{2})+...]*\frac{e^{\beta(E-\mu)}}{(e^{\beta(E-\mu)^{2}}+1)^{2}}$$

$$E^{\alpha+\frac{1}{2}}=\mu^{\alpha+\frac{1}{2}}+(E-\mu)(\alpha+\frac{1}{2})\mu^{\alpha-\frac{1}{2}}+\frac{1}{2}(E-\mu)^{2}(\alpha+\frac{1}{2})(\alpha-\frac{1}{2})\mu^{\alpha-\frac{3}{2}}+...$$

$$=\frac{1}{\alpha+\frac{1}{2}}\int_{-\beta\mu}^{\infty}dx[\mu^{\alpha+\frac{1}{2}}+x\mu^{\alpha-\frac{1}{2}}(\alpha+\frac{1}{2})+x^{2}\frac{\mu^{\alpha-\frac{3}{2}}}{2\beta^{2}}(\alpha+\frac{1}{2})(\alpha-\frac{1}{2})+...]\frac{e^{x}}{(e^{x}+1)^{2}}$$

Where $\beta(E-\mu)=x$ and $\beta dE=dx$.

Now, the integrand is exponentially small between $x=-\beta\mu$ and $x=-\infty$.  So we can state

$$\cong\frac{1}{\alpha+\frac{1}{2}}\int_{-\infty}^{\infty}dx[\mu^{\alpha+\frac{1}{2}}+x(\alpha+\frac{1}{2})\mu^{\alpha-\frac{1}{2}}+x^{2}\frac{\mu^{\alpha-\frac{3}{2}}}{2\beta^{2}}(\alpha+\frac{1}{2})(\alpha-\frac{1}{2})+....]*\frac{e^{x}}{(e^{x}+1)^{2}}
$$

But $\int_{-\infty}^{\infty}dx\frac{x^{n}e^{x}}{(e^{x}+1)^{2}}=0$ for odd n, $\int_{-\infty}^{\infty}dx\frac{e^{x}}{(e^{x}+1)^{2}}=1$, $\int_{-\infty}^{\infty}dx\frac{x^{2}e^{x}}{(e^{x}+1)^{2}}=\frac{\pi^{2}}{3}$, ...

so

$$\int_{0}^{\infty}dEE^{\alpha-\frac{1}{2}}n(E)\cong\frac{1}{\alpha+\frac{1}{2}}[\mu^{\alpha+\frac{1}{2}}+\frac{\pi^{2}}{6}(\alpha^{2}-\frac{1}{4})\bltz^{2}T^{2}+\theta((\frac{\bltz T}{\mu})^{4})]$$

%Insert diagram

Using the Sommerfeld Expansion we find:

$$N=gV\int\dtk n(k)=\frac{gV}{2\pi^{2}}\int_{0}^{\infty}dkn(k)-\frac{gV}{2\pi^{2}}\frac{\sqrt{2m}}{\hbar}\frac{m}{\hbar^{2}}\int_{0}^{\infty}dEE^{\frac{1}{2}}n(E)$$

$\alpha=1$ in the Sommerfeld formula.  The integral $ \int_{0}^{\infty}dEE^{\frac{1}{2}}n(E)$ is equivalent to $\frac{1}{\frac{3}{2}[\mu^{\frac{3}{2}}+\frac{\pi^{2}}{6}\frac{3}{4}\frac{(\bltz T)^{2}}{n^{\frac{1}{2}}}]}$

$$\cong gV\frac{\sqrt{2}}{3\pi^{2}}(\frac{m\mu}{\hbar^{2}})^{\frac{3}{2}}[1+\pi^{2}(\frac{\bltz T}{\mu})^{2}+...]$$

Where $\frac{\sqrt{2}}{3\pi^{2}}=\frac{k_{f}^{3}}{6\pi^{2}}$.

$$E=gV\int\dtk E_{k}n(k)=gV\frac{\hbar^{2}}{2\pi^{2}}\frac{1}{2m}\int_{0}^{\infty}dkk^{4}n(k)$$
$$=gV\frac{1}{2\pi^{2}}\qmom(\frac{2m}{\hbar^{2}})^{\frac{3}{2}}\frac{m}{\hbar^{2}}\int_{0}^{\infty}dEE^{\frac{3}{2}}n(E)$$

$\alpha=2$ in the Sommerfeld formula.  $\int_{0}^{\infty}dEE^{\frac{3}{2}}n(E)$ is equivalent to $\frac{1}{\frac{5}{2}}\mu^{\frac{5}{2}}[1+\frac{\pi^{2}}{6}(\frac{\bltz T}{\mu})^{2}\frac{15}{4}+...]$.

$$\cong gV\qmom(\frac{m\mu}{\hbar^{2}})^{\frac{5}{2}}\frac{2^{\frac{3}{2}}}{5\pi^{2}}[1+5\pi^{2}(\frac{\bltz T}{\mu})^{2}+...]$$

In the above expression the '1' term in the summation is the Fermi Sphere and the subsequent terms are the particle hole excitations.  We can invert the $\alpha=1$ expression to find $\mu$ as a function of $N$ up to higher orders in $\frac{\bltz T}{\mu}$:

$$(\frac{m\mu}{\hbar^{2}})^{\frac{3}{2}}=\frac{N}{gV}\frac{3\pi^{2}}{\sqrt{2}}\frac{1}{1+\frac{\pi^{2}}{8}(\frac{\bltz T}{\mu})^{2}}$$

$$\mu^{\frac{3}{2}}\cong\frac{\hbar^{3}}{m^{\frac{3}{2}}}(\frac{N}{gV}\frac{3\pi^{2}}{\sqrt{2}})^{\frac{2}{3}}[1-\frac{\pi^{2}}{8}(\frac{\bltz T}{\mu})^{2}]^{\frac{2}{3}}$$

$\frac{\hbar^{3}}{m^{\frac{3}{2}}}(\frac{N}{gV}\frac{3\pi^{2}}{\sqrt{2}})^{\frac{2}{3}}$ is equivalent to $[\qmom(\frac{N}{V}\frac{6\pi^{2}}{g}^{\frac{2}{3}})]^{\frac{3}{2}}$.  Within the brackets, 1 is the leading term and $\frac{\pi^{2}}{8}(\frac{\bltz T}{\mu})^{2}]^{\frac{2}{3}}$ is the $\theta(\frac{\bltz T}{\mu})$ correction term.

Therefore,

$$\mu^{\frac{3}{2}}\cong\epsilon_{f}^{\frac{3}{2}}[1-\frac{\pi^{2}}{12}(\frac{\bltz}{\epsilon_{f}})^{2}]$$

Within the brackets, we see the difference between $\epsilon_{f}$ and $\mu$ is higher order in $\frac{\bltz T}{\epsilon_{f}}$.

$$E\cong\qmom gV(\frac{m\mu}{\hbar^{2}})^{\frac{5}{2}}\frac{2^{\frac{3}{2}}}{5\pi^{2}}[1+\frac{5\pi^{2}}{8}(\frac{\bltz T}{\mu})^{2}]$$

$$\cong\frac{3}{5}N\epsilon_{f}[1+\frac{5}{12}\pi^{2}(\frac{T}{T_{f}})^{2}+...]$$

Again, within the brackets we have the '1' term representing the Fermi Sphere, the subsequent term representing particle-hole excitations which depend on $T_{f}$: the Fermi Temperature.

The Specific Heat at Constant Volume is:

$$c_{v}=\frac{1}{N}\thermdif{E}{T}{V}{N}=\frac{3}{5}\epsilon_{f}\frac{5\pi^{2}}{12}\frac{2T}{T_{f}^{2}}=\frac{\pi^{2}}{2}\frac{\bltz T}{T_{f}}$$

One thing that is important to note is that the specific heat is linear in T, as opposed to the classical $c_{v}~\frac{3N}{2}\bltz$.

\section{Ideal Bose Gas}

For an ideal Bose gas, we have the following occupation number expression:

$$n(k)=\frac{1}{e^{\beta(E_{k}-\mu)}-1}$$

%insert graph

This expression diverges at $E_{k}=\mu$ if $\mu>0$.  That's because there is no energy cost in adding bosons to the system.  Therefore, the chemical potential must be negative.

This means that 

$$N=gV\int\dtk\frac{1}{\frac{1}{z}e^{\beta E_{k}}-1}=gV\frac{1}{2\pi^{2}}\frac{m}{\beta\hbar^{2}}\int_{0}^{\infty}dx\sqrt{\frac{2m}{\beta\hbar^{2}}}\frac{\sqrt{x}}{\frac{1}{z}e^{x}-1}$$

Where $\frac{\beta\hbar^{2}k^{2}}{2m}=x$ and $\frac{\beta\hbar^{2}k}{m}dk=dx$.  The variable $z$ represents a special property of the system known as the \textbf{fugacity}.  In this case we know that $0<z=e^{\beta\mu}<1$.

$$N=gV\frac{1}{4\pi^{2}}(\frac{2m\bltz T}{\hbar^{2}})^{\frac{3}{2}}\frac{\sqrt{\pi}}{2}g_{\frac{3}{2}}(z)$$

with $$g_{\alpha}(z)=\frac{1}{\Gamma(\alpha)}\int_{0}^{\infty}dx\frac{x^{\alpha-1}}{\frac{1}{z}e^{x}-1}$$

$$=\frac{1}{\Gamma(\alpha)}\int_{0}^{\infty}dxx^{\alpha-1}\sumser{n=1}{\infty}e^{-nx}z^{n}$$

$$=\sumser{n=1}{\infty}\frac{z^{n}}{n^{\alpha}}$$

so $g_{\frac{3}{2}}(1)=\sumser{n=1}{\infty}\frac{1}{n^{\frac{3}{2}}}=\zeta(\frac{3}{2})\cong2.61$  As a reminder, $\zeta$ is the Riemann Zeta Function.  

%Insert graph

$$N=gV(\frac{m\bltz T}{2\pi\hbar^{2}})^{\frac{3}{2}}g_{\frac{3}{2}(z)}$$

We can see that $(\frac{m\bltz T}{2\pi\hbar^{2}})^{\frac{3}{2}}=\frac{1}{\lambda_{T}^{3}}$  It seems that if $\frac{N}{V}$and T are chosen so that $\frac{N}{gV}\lambda_{T}^{3}>\zeta(\frac{3}{2})$.  However, there is no such value of z that would satisfy this equation.  The mistake we made was in approximating the sum over single particle states by an integral.  This approximation is valid if and only if n(k) is smooth.  But the occupation number is singular at $\mathbb{k}=0$: a macroscopic fraction of the particles occupy the $\mathbb{k}=0$ state and explains the difference between N and $\frac{gV}{\lambda_{T}^{3}}g_{\frac{3}{2}}(1).$  This is obvious at $T=0$.  All particles want (and can) to fall to the $\mathbb{k}=0$ state:
	
$$N=gV\int\dtk\frac{1}{\frac{1}{z}e^{\beta E_{k}}-1}+N_{0}$$

Where $N_{0}$ is the number of particles in the $\mathbb{k}=0$ mode. 

$$=\frac{gV}{\lambda_{T}^{3}}g_{\frac{3}{2}}(z)+N_{0}$$

As long as $\frac{N}{gV}<\frac{\zeta(\frac{3}{2})}{\lambda_{t}^{3}}$, $N_{0}$ vanishes.  If $\frac{N}{gV}>\frac{\zeta(\frac{3}{2})}{\lambda_{T}^{3}}$, $N_{0}$ is given by $N_{0}=N-\frac{gV}{\lambda_{T}^{3}}\zeta(\frac{3}{2})$.  The transition happens at 

$$\frac{1}{\lambda_{T}^{3}}=(\frac{m\bltz T_{c}}{\hbar^{2}})^{\frac{3}{2}}=\frac{N}{gV}\frac{1}{\zeta(\frac{3}{2})}$$

or

$$\bltz T_{c}=\frac{\hbar^{2}}{m}n^{\frac{2}{3}}\frac{1}{\zeta(\frac{3}{2})^{\frac{2}{3}}}$$
	
Where $T_{c}$ is the \textbf{Bose Condensation Temperature} and $n^{\frac{2}{3}}$ is the density.  Above $T_{c}$, $N_{0}=0$ and $N=\frac{gV}{\lambda_{T}^{3}}g_{\frac{3}{2}}(z)$ is valid.  Notice that Bose Condensation occurs when the thermal wavelength is larger than the typical interparticle distance $d\cong(\frac{N}{gV})^{\frac{1}{3}}$:

$$\lambda_{T}\leq d$$

Implies \textbf{Bose-Einstein Condensation}.  When the particle wavelength begin to overlap, patterns of interference cause the individual particle waves to sync up.  When this occurs, some very strange properties arise.  Thermodynamic quantities have different expressions above and below $T_{c}$:

$$\frac{N_{0}}{N}=1-\frac{gV}{N}\frac{\zeta(\frac{3}{2})}{\lambda_{T}^{3}}, 0=1-\frac{gV}{N}\frac{\zeta(\frac{3}{2})}{\lambda_{T}^{3}}$$
$$\downarrow$$
$$1-(\frac{T}{T_{c}})^{\frac{3}{2}}, T<T_{c}$$
$$0, T>T_{c}$$

$$E=gV\int\dtk\frac{E_{k}}{\frac{1}{z}e^{\beta E_{k}}-1}=gV\frac{1}{2\pi^{2}}\frac{\hbar^{2}}{2m}\int_{0}^{\infty}dk\frac{k^{4}}{\frac{1}{z}e^{\beta\qmom k^{2}}-1}$$

$$\frac{\beta\hbar^{2}k^{2}}{2m}=x$$
$$\frac{\beta\hbar^{2}}{m}kdk=dx$$

$$=gV\frac{1}{2\pi^{2}}\qmom\frac{m}{\beta\hbar^{2}}(\frac{2m}{\hbar^{2}})^{\frac{3}{2}}\int_{0}^{\infty}dx\frac{x^{\frac{3}{2}}}{\frac{1}{z}e^{x}-1}$$
$$=\frac{gV}{\beta}(\frac{m\bltz T}{2\pi\hbar^{2}})^{\frac{3}{2}}\frac{3}{2}\frac{4}{3\sqrt{\pi}}\int_{0}^{\infty}dx\frac{x^{\frac{3}{2}}}{\frac{1}{z}e^{x}-1}$$

Note that $\frac{4}{3\sqrt{\pi}}\int_{0}^{\infty}dx\frac{x^{\frac{3}{2}}}{\frac{1}{z}e^{x}-1}$ is $g_{\frac{5}{2}}(z)$.

$$=\frac{3}{2}\frac{gV\bltz T}{\lambda_{T}^{3}}g_{\frac{5}{2}}(z), T>T_{c}$$
$$=\frac{3}{2}\frac{gV\bltz T}{\lambda_{T}^{3}}g_{\frac{5}{2}}(1), T<T_{c}$$

The specific heat at low temperatures is

$$c_{v}=\frac{1}{N}\thermdif{E}{T}{V}{N}=\frac{3}{2}\frac{gV}{N}g_{\frac{5}{2}}(1)(\frac{m}{2\pi\bltz})^{\frac{3}{2}}(\bltz T)^{\frac{3}{2}}\frac{5}{2}$$
$$=\frac{15}{4}\frac{gV}{N}\frac{g_{\frac{5}{2}}(1)}{\lambda_{T}^{3}}$$
$$~T^{\frac{3}{2}}$$

As expected, since the number of particles outside of the condensate is $~T^{\frac{3}{2}}$; the condensate does not carry energy.

An experimental realization of a Bose-Einstein Condensation requires low T, weak interactions, and low mass (in order to make $\lambda_{T}$ larger).  Most substances freeze solid well before they reach $T_{c}$ due to their properties.  However, Helium (4He, to be precise, 3He is a fermion) as a noble gas has small interactions and a small mass and has been forced into a condensate state in experimental settings.  Some of the properties derived above are \textit{qualitatively} correct for liquid 4He but the interactions are still significant so the ideal gas approximation is not quantitatively correct.  4He around $T\cong T_{c}$ is a liquid, its atoms are almost touching each other and the interactions are frequent.  Still, Bose-Einstein condensation does occur in 4He and leads to a bunch of very unusual properties (zero viscosity, extreme heat conduction, etc.) called \textbf{superfluidity}.  To properly study superfluidity we need better tools ("second quantization").  Dilute cold gases of alkali atoms (Na, K,...) were also shown to exhibit Bose-Einstein condensation at extreme low temperatures achievable in atomic traps.  These gases are very close to being ideal.

\section{Phonons}

Although it may seem as though phonons are unrelated to bosonic gases, their behavior is actually very similar.  Consider a lattice of ions (one dimensional for simplicity) connected to each nearest neighbor by a spring:

%insert diagram

Where $x_{n}(t)$ is the \textit{deviation} from the spring equilibrium position of the $n^{th}$ ion.  The classical hamiltonian for such a system would be:

$$\mathcal{H}=\sumser{n=1}{N}[\frac{P_{n}^{2}}{2m}+\frac{m\omega^{2}}{2}(x_{n}-x_{n-1})^{2}]$$

The difference between $X_{n}$ and $x_{n-1}$ is the stretching/compression of the spring.  Notice how the hamiltonian is simply the hamiltonian for a system of coupled harmonic oscillators.  

Again, for simplicity, let's assume $x_{0}=x_{n}$ (periodic boundry conditions).  Change the variables to express the normal modes:

$$x_{n}(t)=\frac{1}{N}\sum_{k}e^{ikan}Q_{k}(t)$$
$$p_{n}(t)=\frac{1}{N}\sum_{k}e^{ikan}P_{k}(t)$$

With $a$ being the lattice spacing and k in the range (1st Brillouin Zone):

%Insert diagram

Let's rewrite $\mathcal{H}$ in terms of the new coordinates $Q_{n},P_{n}$:

$$\sumser{n=1}{N}P_{n}^{2}=\frac{1}{N^{2}}\sumser{n=1}{N}\sum_{k}\sum_{k'}e^{i(k+k')an}P_{k}P_{k'}$$

Notice that the expression $ \sumser{n=1}{N}e^{i(k+k')an}$ is simply the delta function multiplied by N, i.e. $N\delta_{k_{1}-k'}$.

$$\sumser{n=1}{N}=\frac{1}{N}\sum_{k}P_{k}P_{-k}=\frac{1}{N}\sum_{k}P_{k}P_{k}^{*}$$

$$\sumser{n=1}{N}(x_{n}-x_{n-1})^{2}=\sumser{n=1}{N}(x_{n}^{2}+x_{n-1}^{2}-x_{n}x_{n-1}-x_{n+1}x_{n})$$

$$=\sumser{n=1}{N}\frac{1}{N^{2}}\sum_{k_{1}k'}Q_{k}Q_{k'}(2e^{ian(k+k')}-e^{iakn+iak'(n-1)}-e^{iak(n+1)+iak'n})$$

$$\sum_{k}Q_{k}Q_{-k}(2-2\cos{ka})=\sum_{k}Q_{k}Q_{k}^{*}4\sin^{2}(\frac{ka}{2})$$

Therefore we can write the hamiltonian as

$$\mathcal{H}=\frac{1}{N}\sum_{k}[\frac{P_{k}^{*}P_{k}}{2m}+\frac{m\omega^{2}4\sin^{2}(\frac{ka}{2})}{2}Q_{k}^{*}Q_{k}]=\frac{1}{N}\sum_{k}[\frac{P_{k}^{*}P_{k}}{2m}+\omega_{k}^{2}Q_{k}^{*}Q_{k}]$$

As we see above, in our system we have 2N \textit{uncoupled} oscillators (N from real part, and N from imaginary) but only N are independent because $Q_{k}^{*}=Q_{-k}$

Classically, each oscillator contributes with $\int dP_{k}dQ_{k}e^{-\beta(\frac{P_{k}^{2}}{2m}+\frac{m}{2}\omega_{k}^{2}Q_{k}^{2})}$ to the partition function.  Quantum mechanically we have:

$$Z=Tr[e^{-\beta\hat{H}}]=\sum_{n_{1},n_{2},n_{3},...}<n_{1}n_{2}...\vert e^{-\beta\hat{H}}\vert n_{1}n_{2}...>$$

As before, Z is the canonical partition function, and $n_{1}...$ is the excitation level of each oscillator.  Continuing:

$$=\sumser{n_{1}=0}{\infty}e^{-\beta\hbar\omega_{k_{1}}(n_{1}+\frac{1}{2})}\sumser{n_{2}=0}{\infty}e^{-\beta\hbar\omega_{k_{2}}(n_{2}+\frac{1}{2})}....$$

$$=\frac{e^{-\beta\hbar\frac{\omega_{k1}}{2}}}{1-e^{-\beta\hbar\omega_{k1}}}\frac{e^{-\beta\hbar\frac{\omega_{k2}}{2}}}{1-e^{-\beta\hbar\omega_{k2}}}....$$

$$F=-\bltz T\ln{z}=-\bltz T\sum_{k}\ln{\frac{e^{-\beta\hbar\frac{\omega_{k}}{2}}}{1-e^{-\beta\hbar\omega_{k}}}}$$

$$=\sum_{\mathbb{k}}\frac{\omega_{k}}{2}+\bltz T\sum_{\mathbb{k}}\ln{1-e^{-\beta\hbar\omega_{k}}}$$

Where the first term is the T-independant Zero-point energy and the second term is a gas of bosonic "quasi-particles" with dispersion relation $\omega_{\mathbb{k}}=2\omega\vert\sin{\frac{ka}{2}}\vert$ and $\mu=0$.

%insert diagram

The state |0000...|...00> has all the properties we expect of a single boson state with a dispersion relation which, at small momenta, is $\omega_{k}=Vk$.  These particles can be created or destroyed (that is, the harmonic oscillator state can be raised or lowered) by heating/cooling/shaking the system.  Thus, their number is not conserved and is not an independent variable but is determined by T,V, etc...  Consequently, there is no $\mu$ cooresponding to the number of phonons.  Had we included higher powers of $X_{n}$ in the hamiltonian the 1,2,3,... phonon states (states like |11000...> or |20000...>) would no longer be eigenstates of $\hat{H}$ so |t=0>=|110000...> could evolve in time to, say, |t>0>=|100100...> or |1111000...>.  In other words, two phonons could scatter from each other and change their momenta or even create other phonons:

%insert scattering diagram

Differently from "regular" particles, there is a maximum momentum $k=\frac{\pi}{a}$ for phonons (actually, we don't know if that is true; how can we be sure there is no maximum momenta for an electron?).  


In three dimensions the dispersion relation is complicated and dependant on the crystal structure.  It's, in general, not rotationally symmetric and, since atoms can move in three directions there are three modes of oscillation for every momentum $\mathbb{k}$ one longitudinal and two transverse:

%Insert mode diagrams

All three modes obey $\omega=ck$ (c is the speed of sound, not light) for small $\mathbb{k}$ but the transverse and longitudinal phonons but with different values of c for the two kinds of modes.

Let us now compute the specific heat of a solid due to lattice vibrations(that is, phonons) by making \textbf{The Debye Approximation}:

\begin{itemize}
	\item $\omega=ck$ for all 3 modes with the same c
	\item Brillouin Zone (range of allowed $\mathbb{k}$, in the cubic lattice case $-\frac{\pi}{a}<k_{x,y,z}<\frac{\pi}{a}$) is approximated by a sphere with radius $k=k_{b}$ given by: $3N=V\int^{k_0}\dtk$ where we have 3N atom degrees of freedom and $V\int^{k_D}\dtk$ fourier modes.
\end{itemize} 

Note that the number of degrees of freedom is the same whether we use $\mathbb{r}_{n}$ or $\mathbb{Q}_{\mathbb{k}}$ to describe them.  Therefore, we have:

$$c_{v}=\frac{1}{N}\thermdif{E}{T}{V}{N}=\frac{1}{N}V\frac{\delta}{\delta T}\int^{\bltz}\dtk3\frac{\bar{\omega_{k}}}{e^{\beta\hbar\omega_{k}}-1}$$

$$=\frac{V}{N}3\frac{1}{\bltz T^{2}}\int^{k_{0}}\dtk\frac{(\hbar\omega_{k})^{2}}{(e^{\beta\hbar\omega_{k}}-1)^{2}}e^{\beta\hbar\omega_{k}}$$

$$=3\frac{V}{N}\frac{1}{2\pi^{2}}\frac{(\hbar c)^{2}}{\bltz T^{2}}(\frac{\bltz T}{\hbar c})^{5}\int_{0}^{\beta\hbar ck_{D}}dx\frac{x^{4}}{(e^{x}-1)^{2}}e^{x}$$

Where $x=\frac{\hbar c}{\bltz T}k$.  

For $\bltz T << \hbar ck_{D}$:

$$c_{v}=3\frac{V}{N}\bltz\frac{(\bltz T)^{2}}{(\hbar c)^{3}}\frac{1}{2\pi^{2}}\int_{0}^{\infty}dx\frac{x^{4}e^{x}}{(e^{x}-1)^{2}}=3\frac{V}{N}\bltz\frac{(\bltz T)^{2}}{(\hbar c)^{3}}\frac{1}{2\pi^{2}}\frac{4\pi^{4}}{15}$$

$$c_{v}=\frac{V}{N}\frac{2\pi^{2}}{5}(\frac{\bltz T}{\hbar c})^{3}$$

The $~T^{3}$ low temperature behavior is generic for phonons.  Note:  In lattices with two atoms per cell there are vibration modes where the distance between these atoms oscillates.  Their spectrum is gapped (that is, $\omega_{k}\rightarrow\omega_{0}\neq0$ as $k\rightarrow0$) and, since the atom displacement generates an electric dipole, they couple strongly to electromagnetic waves.  They are called \textit{optical phonons}, as opposed to the $\omega\rightarrow ck$ modes that are called \textit{acoustic phonons}.

\section{Photons and Blackbody Radiation}

Material bodies emit electromagnetic radiation when heated.  In general, the amount and frequencies of the radiation depend on the material properties, such as its size, shape, makeup, etc.  A lot of them, however, are close to one idealized limit: \textbf{The Blackbody}.  A blackbody is an object that absorbs (instead of reflecting) all radiation falling onto it.  All blackbodies, regardless of composition, radiate the same way for any given temperature.  There's a simple argument showing that:

%insert blackbody diagram

This shows that the energy emitted by any blackbody of temperature T (per time, area, etc) is the same, regardless of its size, shape, or composition!  It's a universal function $R_{\nu}(T)$:

$$R_{\nu}(T)d\nu$$

Which describes the power emitted per area between frequencies $\nu$ and $\nu+d\nu$.  A practical way of building a blackbody is to build a box with a hole.  Almost all radiation falling on the hole bounces back inside the internal walls and is eventually absorbed.  Thus, the hold is an excellent approximation to a blackbody.  Theoretically, this is also the best way of analysing blackbody radiation.  Since $\brad$ is universal we might as well consider the simplest black body we can imagine:

%insert simple blackbody diagram.

The electromagnetic waves inside the box are in equilibrium with the walls the same way that the molecules inside a box reach thermal equilibrium.  So, let us apply the laws of statistical mechanics to the waves inside the box.  We need to recast Maxwell's Equations in a hamiltonian formalism to apply what we know of statistical mechanics.  For that we use the variables $\mathbb{A}(\mathbb{r},t)$ instead of $\mathbb{E}(\mathbb{r},t)$ and $\mbb{B}(\mbb{r},t)$.  We first notice that $\mbb{A}$'s differing by a \textit{Gauge Transformation} ($\mbb{A}\rightarrow\mbb{A}+\nabla X$) which corresponds to the same $\mbb{E},\mbb{B}$ (the physical observables).  We don't want to double count, that is, include two $\mbb{A}$'s in the partition function corresponding to the same physics, so we impose the following condition:

$$\nabla\mbb{A}=0$$
$$\theta=0$$
$$\mbb{E}=-\pardif{\mbb{A}}{t}-\nabla\theta=-\pardif{\mbb{A}}{t}$$
$$\mbb{B}=\nabla X\mbb{A}$$

Where $\theta$ is the electric potential and $\mbb{A}$ is the magnetic potential.  The lagrangian for the electromagnetic field (without charges or constants) is:

$$\mathcal{L}=\frac{1}{8\pi}\int d^{3}\mbb{r}[(\frac{\delta\mbb{A}}{c\delta t})^{2}-(\nabla X\mbb{A})^{2}]$$

Extremeizing this lagrangian leads to $\frac{1}{c^{2}}\secpardif{\mbb{A}}{t}-\nabla^{2}\mbb{A}=0$, which is equivalent to Maxwell's equations.  Notice that there are an infinite number of possible values of $\mbb{r}$ and, consequently, an infinite number of variables $\mbb{A}(\mbb{r})$ (compare that with $q_{i}$, i=1,...,n) that we had before.  From $\mathcal{L}$ we can find the hamiltonian:

$$\mbb{\Pi}(\mbb{r})=\pardif{\mathcal{H}}{\mbb{A}(\mbb{r})}=\frac{1}{4\pi c^{2}}\pardif{\mbb{A}}{t}$$

$$\mathcal{H}=\int d^{3}\mbb{r}[\mbb{\Pi}(\mbb{r},t)\pardif{\mbb{A}(\mbb{r},t)}{t}-\mathcal{L}]$$
$$\mathcal{H}=\frac{1}{8\pi}\int d^{3}\mbb{r}[(4\pi^{2})^{2}\mbb{\Pi}^{2}(\mbb{r},t)+(\nabla X\mbb{A})^{2}]$$

Note that $(4\pi^{2})^{2}\mbb{\Pi}^{2}(\mbb{r},t)=\mbb{E}^{2}(\mbb{r},t)$ and $(\nabla X\mbb{A})^{2}=\mbb{B}^{2}(\mbb{r},t)$

Compare the equation above to $\mathcal{H}=\sumser{i=1}{N}p_{i}\dot{q}_{i}-\mathcal{L}$.  The hamiltonian is quadratic in the momenta $\mbb{\Pi}$ and $\mbb{A}$ so we should be able to find the normal coordinates which transform it into decoupled harmonic oscillators.  We change variables to $\mbb{Q}_{k}$, $\mbb{P}_{k}$ defined by:

$$\mbb{A}(\mbb{r},t)=\sqrt{\frac{4\pi c}{V}}\sum_{\mbb{k}}\mbb{Q}(\mbb{k},t)e^{i\mbb{k}\mbb{r}}$$
$$\mbb{\Pi}(\mbb{r},t)=\frac{1}{\sqrt{4\pi c^{2}V}}\sum_{\mbb{k}}\mbb{P}(\mbb{k},t)e^{i\mbb{k}\mbb{r}}$$
$$V=L^{3}$$
$$\mbb{k}=\frac{2\pi}{L}(n_{x},n_{y},n_{z})$$
$$n_{x},n_{y},n_{z}=0,\pm1,\pm2,...$$
$$\mbb{Q}^{*}(\mbb{k})=\mbb{Q}(-\mbb{k})$$
$$\mbb{P}^{*}(\mbb{k})=\mbb{P}(-\mbb{k})$$

Since $\mbb{A}$ and $\mbb{\Pi}$ are real.  $\nabla\mbb{A}=0\rightarrow\mbb{Q}(\mbb{k})\mbb{k}=0$, $\mbb{P}(\mbb{k})\mbb{k}=0$

$$\mathcal{H}=\frac{1}{2}\sum_{\mbb{k}}[\mbb{P}(\mbb{k},t)\mbb{P}(-\mbb{k},t)+c^{2}k^{2}\mbb{Q}(\mbb{k})\mbb{Q}(-\mbb{k},t)]$$
$$c^{2}k^{2}=\omega_{k}^{2}$$

This whole expression is the sum of decoupled harmonic oscillators, very much like the phonon hamiltonian.  We can now start doing statistical mechanics!  Classically, the energy density is given by the equipartition theorem:

$$E=\sum_{\mbb{k}}\frac{1}{2}\bltz T*2*2\cong\infty$$

The two '2' terms represent the kinetic+potential energy contributions and the two polarizations per mode, respectively.  The expression is approximately infinity since there are an infinite number of allowed values of $\mbb{k}$

The divergence arises because, contrary to the phonon case, there is no cutoff at high $\mbb{k}$.  Still, the contribution from each frequency range is finite:

$$\frac{E}{V}=\frac{2\bltz T}{2\pi^{2}}\int_{0}^{\infty}dk k^{2}=2\bltz T\int\dtk=\int_{0}^{\infty}d\nu\frac{\bltz T}{\pi^{2}}\frac{(2\pi)^{3}}{c^{3}}\nu^{2}$$

Which is, of course, the Rayleigh-Jeans formula.  The failure of this formula at large frequencies is what precipitated the advent of quantum theory.  

Now that we have the classical case out of the way, we will now work out the quantum theory.  We start by assuming that the coordinates $\mbb{A}(\mbb{r})$ and momenta $\mbb{\Pi}(\mbb{r})$ become operators satisfying the usual commutation rules.  That implies that $\mbb{Q}$ and $\mbb{P}$ become operators satisfying:

$$[\hat{Q}_{i}(\mbb{k}),\hat{P}_{j}(\mbb{q}^{'})]=i\hbar\delta_{ij}\delta_{\mbb{k1-k'}}$$

All the other commutators including Q's and P's vanish.  We now define:

$$\hat{\alpha}^{\alpha}(\mbb{k})=[\sqrt{\frac{\omega_{k}}{2\hbar}}\hat{\mbb{Q}}({\mbb{k}})-\frac{i}{\sqrt{2\hbar\omega_{k}}}\hat{\mbb{P}}(\mbb{k})]\mbb{\epsilon}_{\mbb{k}}^{\alpha}$$
$$\hat{\alpha}^{\dagger\alpha}(\mbb{k})=[\frac{\omega_{k}}{2\hbar}\hat{\mbb{Q}}({\mbb{-k}})-\frac{i}{\sqrt{2\hbar\omega_{k}}}\hat{\mbb{P}}(\mbb{-k})]\mbb{\epsilon}_{\mbb{k}}^{\alpha}$$

then,

$$[\hat{a}^{\alpha}(\mbb{k}),\hat{a}^{\beta\dagger}(\mbb{k^{'}})]=-i\sqrt{\frac{\omega_{k}}{2\hbar}}\frac{1}{\sqrt{2\hbar\omega_{k'}}}[\hat{Q}_{i}(\mbb{k}),\hat{P}_{j}(\mbb{k'})]\epsilon^{\alpha}_{i}(\mbb{k})\epsilon_{j}^{\beta}(\mbb{k'})-i\sqrt{\frac{\omega_{k'}}{2\hbar}}\frac{1}{\sqrt{2\hbar\omega_{k}}}[\hat{P}_{j}(\mbb{k}),\hat{Q}_{i}(\mbb{k'})]\epsilon_{j}^{\alpha}\epsilon^{\beta}_{i}(\mbb{k'})=\delta_{\alpha\beta}\delta_{ij}$$

$$[\alow,\hat{a}^{\beta}(\mbb{k'})]=[\aras,\hat{a}^{\beta^{\dagger}}(\mbb{k'})]=0$$

Thus, $\aras$ is the raising operator of the harmonic oscillator made with momentum $\mbb{k}$ and polarization $\alpha$.  In other words, $\aras$ creates a photon with momentum $\mbb{k}$ and polarization $\alpha$.  Similarly, $\alow$ destroys a photon.  The hamiltonian can be written in terms of creation/destruction operators:

$$\hat{H}=\frac{1}{2}\sum_{\alpha=1,2}V\int\dtk[-\frac{\hbar\omega_{k}}{2}(\alow-\hat{a}^{\alpha^{\dagger}}(-\mbb{k}))(\hat{a}^{\alpha}(-\mbb{k})-\aras)+\omega_{k}^{2}\frac{\hbar}{2\omega_{k}}(\alow+\hat{a}^{\alpha^{\dagger}}(-\mbb{k}))(\hat{a}^{\alpha}(-\mbb{k})+\aras)]$$

$$=\frac{V}{2}\sum_{\alpha=1,2}\int\dtk[\frac{\hbar\omega_{k}}{2}2\alow\aras+\frac{\hbar\omega_{k}}{2}2\aras\alow]$$
$$=V\sum_{\alpha=1,2}\int\dtk[\aras\alow+\frac{1}{2}]\hbar\omega_{k}$$

We have one oscillator per value of $\mbb{k}$ and $\alpha$.  The calculation of the partition function now follows the previous example of the phonons:

$$Z=Tr[e^{-\beta\hat{H}}]=\Pi_{\mbb{k},\alpha}\sumser{n_{k}=0}{\infty}e^{-\beta\hbar\omega(n_{k}+\frac{1}{2})}$$
$$=\Pi_{\mbb{k},\alpha}\frac{e^{-\beta\hbar\frac{\omega_{k}}{2}}}{1-e^{-\beta\hbar\omega_{k}}}=\Pi_{\mbb{k},\alpha}\frac{2}{\sinh{\beta\hbar\frac{\omega_{k}}{2}}}$$

And consequently,

$$F=-\bltz T\ln{Z}=\bltz T\sum_{\mbb{k},\alpha}\ln{\frac{1-e^{-\beta\hbar\omega_{k}}}{e^{-\beta\hbar\frac{\omega_{k}}{2}}}}$$

$$=\bltz T\sum_{\mbb{k},\alpha}[\ln{1-e^{-\beta\hbar\omega_{k}}+\bltz T\beta\hbar\omega_{k}}]$$

The term $\bltz T\beta\hbar\omega_{k}$ is T-independent, zero-point energy.  This time this term is infinite because there is no cutoff in $\mbb{K}$.  It is a contribution to the cosmological constant.

$$\cong2\bltz TV\int\dtk\ln{1-e^{-\beta\hbar\omega_{k}}}+const$$

Here we have 2 polarizations multiplied by the free energy of an ideal, massless gas of bosons plus a constant.  Now we can find the energy density:

$$\frac{E}{V}=2\int\dtk\frac{\hbar ck}{e^{\beta\hbar ck}-1}=\frac{2}{(\beta\hbar c)^{3}}\frac{1}{2\pi^{2}}\frac{1}{\beta}\int_{0}^{\infty}dx\frac{x^{3}}{e^{x}-1}$$
$$x=\beta\hbar ck$$
$$\int_{0}^{\infty}dx\frac{x^{3}}{e^{x}-1}=\frac{\pi^{4}}{15}$$

$$\frac{E}{V}\frac{\pi^{2}(\bltz T)^{4}}{15\hbar^{3}c^{3}}$$

The above result is commonly known as the \textbf{Stefan-Boltzman Law}.

Also,

$$\frac{E}{V}=2\frac{1}{2\pi^{2}}\int_{0}^{\infty}dkk^{2}\frac{\hbar ck}{e^{\beta\hbar ck}-1}$$
$$=\frac{\hbar c}{\hbar^{2}}(\frac{2\pi}{c})^{4}\int_{0}^{\infty}d\nu\frac{\nu^{3}}{e^{\beta\hbar\nu}-1}$$
$$=\int_{0}^{\infty}d\nu\frac{16\pi^{2}\hbar}{c^{3}}\frac{\nu^{2}}{e^{\beta\hbar\nu}-1}$$

$$\frac{16\pi^{2}\hbar}{c^{3}}\frac{\nu^{2}}{e^{\beta\hbar\nu}-1}=u_{T}(\nu)$$

$u_{T}(\nu)$ is Planck's Law, which defines the energy per volume in the frequency range ($\nu,\nu+d\nu$).

\part{Non-Ideal Systems}

In this part we will explore approximation methods that allow us to approximately solve various non-ideal systems.

\section{Cluster Expansion}

\begin{itemize}
	\item Method to approximately compute the thermodynamic properties of gases
	\item Results in an expansion in powers of the density.  Thus, it works well at low density.  It's a well controlled approximation in the sense that it becomes better and better as the density becomes smaller.
	\item It's an example of a diagrammatic expansion.  Many features generalize to many other diagrammatic expansions.
\end{itemize}

The goal is to compute 

$$Z=\int\prod_{i=1}^{N}\frac{d^{3}\mbb{r}_{i}d^{3}\mbb{p}_{i}}{N!h^{3N}}e^{-\beta[\sumser{i=1}{N}\qmom+\sum_{i<j}V(\mbb{r}_{i}-\mbb{r}_{j})]}$$

We include the N! denominator to account for indistinguishability and the $h^{3N}$ because it makes Z dimensionless.  Both of these modifications are justified by the quantum calculation.

 



\end{document}