• 358,603

# Posts Tagged ‘Fermion’

## A final pre-exam update of my notes compilation for ‘PHY452H1S Basic Statistical Mechanics’, Taught by Prof. Arun Paramekanti

Posted by peeterjoot on April 22, 2013

Here’s my third update of my notes compilation for this course, including all of the following:

April 21, 2013 Fermi function expansion for thermodynamic quantities

April 20, 2013 Relativistic Fermi Gas

April 10, 2013 Non integral binomial coefficient

April 10, 2013 energy distribution around mean energy

April 09, 2013 Velocity volume element to momentum volume element

April 04, 2013 Phonon modes

April 03, 2013 BEC and phonons

April 03, 2013 Max entropy, fugacity, and Fermi gas

April 02, 2013 Bosons

April 02, 2013 Relativisitic density of states

March 28, 2013 Bosons

plus everything detailed in the description of my previous update and before.

## Fermi-Dirac function expansion for thermodynamic quantities

Posted by peeterjoot on April 21, 2013

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

In section 8.1 of [1] are some Fermi-Dirac \index{Fermi-Dirac function} expansions for $P$, $N$, and $U$. Let’s work through these in detail.

Our starting point is the relations

\begin{aligned}P V \beta = \ln Z_{\mathrm{G}} = \sum \ln \left( 1 + z e^{-\beta \epsilon} \right)\end{aligned} \hspace{\stretch{1}}(1.0.1a)

\begin{aligned}N = \sum \frac{1}{{ z^{-1} e^{\beta \epsilon} + 1 }}.\end{aligned} \hspace{\stretch{1}}(1.0.1b)

Recap. Density of states

We’ll employ the 3D non-relativisitic density of states

\begin{aligned}\mathcal{D}(\epsilon) &= \sum_\mathbf{k} \delta(\epsilon - \epsilon_\mathbf{k}) \\ &\sim V \int \frac{d^3 \mathbf{k}}{(2 \pi)^3}\delta(\epsilon - \epsilon_\mathbf{k}) \\ &= \frac{4 \pi V}{(2 \pi)^3}\int dk k^2 \delta\left( \epsilon - \frac{\hbar^2 k^2}{2 m} \right) \\ &= \frac{4 \pi V}{(2 \pi)^3}\int dk k^2 \frac{ \delta\left( k - \sqrt{2 m \epsilon}/\hbar \right)}{ \frac{\hbar^2}{m} \frac{\sqrt{2 m \epsilon}}{\hbar}} \\ &= \frac{2 V}{(2 \pi)^2 }\frac{m}{\hbar^2}\sqrt{\frac{2 m \epsilon}{\hbar^2}},\end{aligned} \hspace{\stretch{1}}(1.0.1b)

or

\begin{aligned}\boxed{\mathcal{D}(\epsilon)=\frac{V}{(2 \pi)^2 }\left( \frac{2 m}{\hbar^2} \right)^{3/2}\epsilon^{1/2}.}\end{aligned} \hspace{\stretch{1}}(1.0.1b)

Density

Now let’s make our integral approximation of the sum for $N$. That is

\begin{aligned}N &= g \int d\epsilon \mathcal{D}(\epsilon) \frac{1}{{ z^{-1} e^{\beta \epsilon} + 1 }} \\ &= g \frac{V}{(2 \pi)^2 }\left( \frac{2 m}{\hbar^2} \right)^{3/2}\int_0^\infty d\epsilon \frac{\epsilon^{1/2}}{ z^{-1} e^{\beta \epsilon} + 1 } \\ &= g \frac{V}{(2 \pi)^2 \beta^{3/2}}\left( \frac{2 m}{\hbar^2} \right)^{3/2}\int_0^\infty du \frac{u^{1/2}}{ z^{-1} e^{u} + 1 } \\ &= g \frac{V}{(2 \pi)^2 \beta^{3/2}}\left( \frac{2 m}{\hbar^2} \right)^{3/2}\Gamma(3/2) f_{3/2}(z) \\ &= g \frac{V}{(2 \pi)^2 \beta^{3/2}}\frac{\left( 2 m k_{\mathrm{B}} T \right)^{3/2}}{\hbar^3}\frac{1}{{2}} \sqrt{\pi}f_{3/2}(z)\\ &= g V \not{{2}} \pi\frac{\left( 2 m k_{\mathrm{B}} T \right)^{3/2}}{h^3}\frac{1}{{\not{{2}}}} \sqrt{\pi}f_{3/2}(z),\end{aligned} \hspace{\stretch{1}}(1.0.1b)

or

\begin{aligned}\frac{N}{V} = g \frac{\left( 2 \pi m k_{\mathrm{B}} T \right)^{3/2}}{h^3}f_{3/2}(z).\end{aligned} \hspace{\stretch{1}}(1.0.5)

With

\begin{aligned}\lambda = \frac{h}{\sqrt{ 2 \pi m k_{\mathrm{B}} T }},\end{aligned} \hspace{\stretch{1}}(1.0.6)

this gives us the desired density result from the text

\begin{aligned}\boxed{\frac{N}{V}=\frac{g}{\lambda^3} f_{3/2}(z).}\end{aligned} \hspace{\stretch{1}}(1.0.7)

Pressure

For the pressure, we can do the same, but have to integrate by parts

\begin{aligned}P V \beta &= g \sum \ln \left( 1 + z e^{-\beta \epsilon} \right) \\ &\sim g \frac{V}{(2 \pi)^2 }\left( \frac{2 m}{\hbar^2} \right)^{3/2}\int_0^\infty d\epsilon \epsilon^{1/2} \ln \left( 1 + z e^{-\beta \epsilon} \right) \\ &= - g \frac{V}{(2 \pi)^2 }\left( \frac{2 m}{\hbar^2} \right)^{3/2}\int_0^\infty d\epsilon \frac{2}{3} \epsilon^{3/2} \frac{-\beta z e^{-\beta \epsilon} }{ 1 + z e^{-\beta \epsilon} } \\ &= g\frac{V}{(2 \pi)^2 }\left( \frac{2 m}{\hbar^2} \right)^{3/2}\frac{2}{3} \frac{1}{{\beta^{3/2}}}\int_0^\infty dx\frac{x^{3/2}}{z^{-1} e^{x} + 1 } \\ &= g\frac{2}{3} 2 \pi V\frac{\left( 2 m k_{\mathrm{B}} T \right)^{3/2}}{h^3 }\Gamma(5/2)f_{5/2}(z) \\ &= g\frac{2}{3} 2 \pi V\frac{\left( 2 m k_{\mathrm{B}} T \right)^{3/2}}{h^3 }\frac{3}{2} \frac{1}{2} \sqrt{\pi}f_{5/2}(z) \\ &= g V\frac{\left( 2 \pi m k_{\mathrm{B}} T \right)^{3/2}}{h^3 }f_{5/2}(z),\end{aligned} \hspace{\stretch{1}}(1.0.7)

or

\begin{aligned}\boxed{P \beta = \frac{g}{\lambda^3} f_{5/2}(z).}\end{aligned} \hspace{\stretch{1}}(1.0.9)

Energy

The average energy is the last thermodynamic quantity to come very easily. We have

\begin{aligned}U &= - \frac{\partial {}}{\partial {\beta}} \ln Z_{\mathrm{G}} \\ &= - \frac{\partial {T}}{\partial {\beta}} \frac{\partial {}}{\partial {T}} \ln Z_{\mathrm{G}} \\ &= - \frac{\partial {(1/k_{\mathrm{B}} T)}}{\partial {\beta}} \frac{\partial {}}{\partial {T}} P V \beta \\ &= \frac{1}{{k_{\mathrm{B}} \beta^2}}\frac{\partial {}}{\partial {T}} \frac{g V}{\lambda^3} f_{5/2}(z) \\ &= g V k_{\mathrm{B}} T^2f_{5/2}(z)\frac{\partial {}}{\partial {T}} \frac{\left( 2 \pi m k_{\mathrm{B}} T \right)^{3/2}}{h^3} \\ &= \frac{3}{2} \frac{g V k_{\mathrm{B}} T}{\lambda^3}f_{5/2}(z).\end{aligned} \hspace{\stretch{1}}(1.0.9)

From eq. 1.0.7, we have

\begin{aligned}\frac{g V}{\lambda^3} = \frac{N}{f_{3/2}(z) },\end{aligned} \hspace{\stretch{1}}(1.0.11)

so the energy takes the form

\begin{aligned}\boxed{U = \frac{3}{2} N k_{\mathrm{B}} T \frac{f_{5/2}(z)}{f_{3/2}(z) }.}\end{aligned} \hspace{\stretch{1}}(1.0.11)

We can compare this to the ratio of pressure to density

\begin{aligned}\frac{P \beta}{n} = \frac{f_{5/2}(z)}{f_{3/2}(z) },\end{aligned} \hspace{\stretch{1}}(1.0.11)

to find

\begin{aligned}U= \frac{3}{2} N k_{\mathrm{B}} T \frac{P V \beta}{N}= \frac{3}{2} P V,\end{aligned} \hspace{\stretch{1}}(1.0.11)

or

\begin{aligned}\boxed{P V = \frac{2}{3} U.}\end{aligned} \hspace{\stretch{1}}(1.0.11)

# References

[1] RK Pathria. Statistical mechanics. Butterworth Heinemann, Oxford, UK, 1996.

## PHY452H1S Basic Statistical Mechanics. Problem Set 6: Max entropy, fugacity, and Fermi gas

Posted by peeterjoot on April 3, 2013

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

# Disclaimer

## Question: Maximum entropy principle

Consider the “Gibbs entropy”

\begin{aligned}S = - k_{\mathrm{B}} \sum_i p_i \ln p_i\end{aligned} \hspace{\stretch{1}}(1.1)

where $p_i$ is the equilibrium probability of occurrence of a microstate $i$ in the ensemble.

For a microcanonical ensemble with $\Omega$ configurations (each having the same energy), assigning an equal probability $p_i= 1/\Omega$ to each microstate leads to $S = k_{\mathrm{B}} \ln \Omega$. Show that this result follows from maximizing the Gibbs entropy with respect to the parameters $p_i$ subject to the constraint of

\begin{aligned}\sum_i p_i = 1\end{aligned} \hspace{\stretch{1}}(1.2)

(for $p_i$ to be meaningful as probabilities). In order to do the minimization with this constraint, use the method of Lagrange multipliers – first, do an unconstrained minimization of the function

\begin{aligned}S - \alpha \sum_i p_i,\end{aligned} \hspace{\stretch{1}}(1.3)

then fix $\alpha$ by demanding that the constraint be satisfied.

For a canonical ensemble (no constraint on total energy, but all microstates having the same number of particles $N$), maximize the Gibbs entropy with respect to the parameters $p_i$ subject to the constraint of

\begin{aligned}\sum_i p_i = 1,\end{aligned} \hspace{\stretch{1}}(1.4)

(for $p_i$ to be meaningful as probabilities) and with a given fixed average energy

\begin{aligned}\left\langle{{E}}\right\rangle = \sum_i E_i p_i,\end{aligned} \hspace{\stretch{1}}(1.5)

where $E_i$ is the energy of microstate $i$. Use the method of Lagrange multipliers, doing an unconstrained minimization of the function

\begin{aligned}S - \alpha \sum_i p_i - \beta \sum_i E_i p_i,\end{aligned} \hspace{\stretch{1}}(1.6)

then fix $\alpha, \beta$ by demanding that the constraint be satisfied. What is the resulting $p_i$?

For a grand canonical ensemble (no constraint on total energy, or the number of particles), maximize the Gibbs entropy with respect to the parameters $p_i$ subject to the constraint of

\begin{aligned}\sum_i p_i = 1,\end{aligned} \hspace{\stretch{1}}(1.7)

(for $p_i$ to be meaningful as probabilities) and with a given fixed average energy

\begin{aligned}\left\langle{{E}}\right\rangle = \sum_i E_i p_i,\end{aligned} \hspace{\stretch{1}}(1.8)

and a given fixed average particle number

\begin{aligned}\left\langle{{N}}\right\rangle = \sum_i N_i p_i.\end{aligned} \hspace{\stretch{1}}(1.9)

Here $E_i, N_i$ represent the energy and number of particles in microstate $i$. Use the method of Lagrange multipliers, doing an unconstrained minimization of the function

\begin{aligned}S - \alpha \sum_i p_i - \beta \sum_i E_i p_i - \gamma \sum_i N_i p_i,\end{aligned} \hspace{\stretch{1}}(1.10)

then fix $\alpha, \beta, \gamma$ by demanding that the constrains be satisfied. What is the resulting $p_i$?

Writing

\begin{aligned}f = S - \alpha \sum_{j = 1}^\Omega p_j,= -\sum_{j = 1}^\Omega p_j \left( k_{\mathrm{B}} \ln p_j + \alpha \right),\end{aligned} \hspace{\stretch{1}}(1.11)

our unconstrained minimization requires

\begin{aligned}0 = \frac{\partial {f}}{\partial {p_i}}= -\left( k_{\mathrm{B}} \left( \ln p_i + 1 \right) + \alpha \right).\end{aligned} \hspace{\stretch{1}}(1.12)

Solving for $p_i$ we have

\begin{aligned}p_i = e^{-\alpha/k_{\mathrm{B}} - 1}.\end{aligned} \hspace{\stretch{1}}(1.13)

The probabilities for each state are constant. To fix that constant we employ our constraint

\begin{aligned}1 = \sum_{j = 1}^\Omega p_j= \sum_{j = 1}^\Omega e^{-\alpha/k_{\mathrm{B}} - 1}= \Omega e^{-\alpha/k_{\mathrm{B}} - 1},\end{aligned} \hspace{\stretch{1}}(1.14)

or

\begin{aligned}\alpha/k_{\mathrm{B}} + 1 = \ln \Omega.\end{aligned} \hspace{\stretch{1}}(1.15)

Inserting eq. 1.15 fixes the probability, giving us the first of the expected results

\begin{aligned}\boxed{p_i = e^{-\ln \Omega} = \frac{1}{{\Omega}}.}\end{aligned} \hspace{\stretch{1}}(1.0.16)

Using this we our Gibbs entropy can be summed easily

\begin{aligned}S &= -k_{\mathrm{B}} \sum_{j = 1}^\Omega p_j \ln p_j \\ &= -k_{\mathrm{B}} \sum_{j = 1}^\Omega \frac{1}{{\Omega}} \ln \frac{1}{{\Omega}} \\ &= -k_{\mathrm{B}} \frac{\Omega}{\Omega} \left( -\ln \Omega \right),\end{aligned} \hspace{\stretch{1}}(1.0.16)

or

\begin{aligned}\boxed{S = k_{\mathrm{B}} \ln \Omega.}\end{aligned} \hspace{\stretch{1}}(1.0.16)

For the “action” like quantity that we want to minimize, let’s write

\begin{aligned}f = S - \alpha \sum_j p_j - \beta \sum_j E_j p_j,\end{aligned} \hspace{\stretch{1}}(1.0.16)

for which we seek $\alpha$, $\beta$ such that

\begin{aligned}0 &= \frac{\partial {f}}{\partial {p_i}} \\ &= -\frac{\partial {}}{\partial {p_i}}\sum_j p_j\left( k_{\mathrm{B}} \ln p_j + \alpha + \beta E_j \right) \\ &= -k_{\mathrm{B}} (\ln p_i + 1) - \alpha - \beta E_i,\end{aligned} \hspace{\stretch{1}}(1.0.16)

or

\begin{aligned}p_i = \exp\left( - \left( \alpha - \beta E_i \right) /k_{\mathrm{B}} - 1 \right).\end{aligned} \hspace{\stretch{1}}(1.0.16)

Our probability constraint is

\begin{aligned}1 &= \sum_j \exp\left( - \left( \alpha - \beta E_j \right) /k_{\mathrm{B}} - 1 \right) \\ &= \exp\left( - \alpha/k_{\mathrm{B}} - 1 \right)\sum_j \exp\left( - \beta E_j/k_{\mathrm{B}} \right),\end{aligned} \hspace{\stretch{1}}(1.0.16)

or

\begin{aligned}\exp\left( \alpha/k_{\mathrm{B}} + 1 \right)=\sum_j \exp\left( - \beta E_j/k_{\mathrm{B}} \right).\end{aligned} \hspace{\stretch{1}}(1.0.16)

Taking logs we have

\begin{aligned}\alpha/k_{\mathrm{B}} + 1 = \ln \sum_j \exp\left( - \beta E_j/k_{\mathrm{B}} \right).\end{aligned} \hspace{\stretch{1}}(1.0.16)

We could continue to solve for $\alpha$ explicitly but don’t care any more than this. Plugging back into the probability eq. 1.0.16 obtained from the unconstrained minimization we have

\begin{aligned}p_i = \exp\left( -\ln \sum_j \exp\left( - \beta E_j/k_{\mathrm{B}} \right) \right)\exp\left( - \beta E_i/k_{\mathrm{B}} \right),\end{aligned} \hspace{\stretch{1}}(1.0.16)

or

\begin{aligned}\boxed{p_i = \frac{ \exp\left( - \beta E_i/k_{\mathrm{B}} \right)}{ \sum_j \exp\left( - \beta E_j/k_{\mathrm{B}} \right)}.}\end{aligned} \hspace{\stretch{1}}(1.0.16)

To determine $\beta$ we must look implicitly to the energy constraint, which is

\begin{aligned}\left\langle{{E}}\right\rangle &= \sum_i E_i p_i \\ &= \sum_iE_i\left( \frac{ \exp\left( - \beta E_i/k_{\mathrm{B}} \right) } { \sum_j \exp\left( - \beta E_j/k_{\mathrm{B}} \right) } \right),\end{aligned} \hspace{\stretch{1}}(1.0.16)

or

\begin{aligned}\boxed{\left\langle{{E}}\right\rangle = \frac{ \sum_i E_i \exp\left( - \beta E_i/k_{\mathrm{B}} \right)}{ \sum_j \exp\left( - \beta E_j/k_{\mathrm{B}} \right)}.}\end{aligned} \hspace{\stretch{1}}(1.0.16)

The constraint $\beta$ ($=1/T$) is given implicitly by this energy constraint.

Again write

\begin{aligned}f = S - \alpha \sum_j p_j - \beta \sum_j E_j p_j - \gamma \sum_j N_j p_j.\end{aligned} \hspace{\stretch{1}}(1.0.16)

The unconstrained minimization requires

\begin{aligned}0 &= \frac{\partial {f}}{\partial {p_i}} \\ &= -\frac{\partial {}}{\partial {p_i}}\left( k_{\mathrm{B}} (\ln p_i + 1) + \alpha + \beta E_i + \gamma N_i \right),\end{aligned} \hspace{\stretch{1}}(1.0.16)

or

\begin{aligned}p_i = \exp\left( -\alpha/k_{\mathrm{B}} - 1 \right) \exp\left( -(\beta E_i + \gamma N_i)/k_{\mathrm{B}} \right).\end{aligned} \hspace{\stretch{1}}(1.0.16)

The unit probability constraint requires

\begin{aligned}1 &= \sum_j p_j \\ &= \exp\left( -\alpha/k_{\mathrm{B}} - 1 \right) \sum_j\exp\left( -(\beta E_j + \gamma N_j)/k_{\mathrm{B}} \right),\end{aligned} \hspace{\stretch{1}}(1.0.16)

or

\begin{aligned}\exp\left( -\alpha/k_{\mathrm{B}} - 1 \right) =\frac{1}{{\sum_j\exp\left( -(\beta E_j + \gamma N_j)/k_{\mathrm{B}} \right)}}.\end{aligned} \hspace{\stretch{1}}(1.0.16)

Our probability is then

\begin{aligned}\boxed{p_i = \frac{\exp\left( -(\beta E_i + \gamma N_i)/k_{\mathrm{B}} \right)}{\sum_j\exp\left( -(\beta E_j + \gamma N_j)/k_{\mathrm{B}} \right)}.}\end{aligned} \hspace{\stretch{1}}(1.0.16)

The average energy $\left\langle{{E}}\right\rangle = \sum_j p_j E_j$ and average number of particles $\left\langle{{N}}\right\rangle = \sum_j p_j N_j$ are given by

\begin{aligned}\left\langle{{E}}\right\rangle = \frac{E_i\exp\left( -(\beta E_i + \gamma N_i)/k_{\mathrm{B}} \right)}{\sum_j\exp\left( -(\beta E_j + \gamma N_j)/k_{\mathrm{B}} \right)}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

\begin{aligned}\left\langle{{N}}\right\rangle = \frac{N_i\exp\left( -(\beta E_i + \gamma N_i)/k_{\mathrm{B}} \right)}{\sum_j\exp\left( -(\beta E_j + \gamma N_j)/k_{\mathrm{B}} \right)}.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

The values $\beta$ and $\gamma$ are fixed implicitly by requiring simultaneous solutions of these equations.

## Question: Fugacity expansion ([3] Pathria, Appendix D, E)

The theory of the ideal Fermi or Bose gases often involves integrals of the form

\begin{aligned}f_\nu^\pm(z) = \frac{1}{{\Gamma(\nu)}} \int_0^\infty dx \frac{x^{\nu - 1}}{z^{-1} e^x \pm 1}\end{aligned} \hspace{\stretch{1}}(1.36)

where

\begin{aligned}\Gamma(\nu) = \int_0^\infty dy y^{\nu-1} e^{-y}\end{aligned} \hspace{\stretch{1}}(1.37)

denotes the gamma function.

Obtain the behavior of $f_\nu^\pm(z)$ for $z \rightarrow 0$ keeping the two leading terms in the expansion.

For Fermions, obtain the behavior of $f_\nu^\pm(z)$ for $z \rightarrow \infty$ again keeping the two leading terms.

For Bosons, we must have $z \le 1$ (why?), obtain the leading term of $f_\nu^-(z)$ for $z \rightarrow 1$.

For $z \rightarrow 0$ we can rewrite the integrand in a form that allows for series expansion

\begin{aligned}\frac{x^{\nu - 1}}{z^{-1} e^x \pm 1} &= \frac{z e^{-x} x^{\nu - 1}}{1 \pm z e^{-x}} \\ &= z e^{-x} x^{\nu - 1}\left( 1 \mp z e^{-x} + (z e^{-x})^2 \mp (z e^{-x})^3 + \cdots \right)\end{aligned} \hspace{\stretch{1}}(1.38)

For the $k$th power of $z e^{-x}$ in this series our integral is

\begin{aligned}\int_0^\infty dx z e^{-x} x^{\nu - 1} (z e^{-x})^k &= z^{k+1}\int_0^\infty dx x^{\nu - 1} e^{-(k + 1) x} \\ &= \frac{z^{k+1}}{(k+1)^\nu}\int_0^\infty du u^{\nu - 1} e^{- u} \\ &= \frac{z^{k+1}}{(k+1)^\nu} \Gamma(\nu)\end{aligned} \hspace{\stretch{1}}(1.39)

Putting everything back together we have for small $z$

\begin{aligned}\boxed{f_\nu^\pm(z) =z\mp\frac{z^{2}}{2^\nu}+\frac{z^{3}}{3^\nu}\mp\frac{z^{4}}{4^\nu}+ \cdots}\end{aligned} \hspace{\stretch{1}}(1.40)

We’ll expand $\Gamma(\nu) f_\nu^+(e^y)$ about $z = e^y$, writing

\begin{aligned}F_\nu(e^y) &= \Gamma(\nu) f_\nu^+(e^y) \\ &= \int_0^\infty dx \frac{x^{\nu - 1}}{e^{x - y} + 1} \\ &= \int_0^y dx \frac{x^{\nu - 1}}{e^{x - y} + 1}+\int_y^\infty dx \frac{x^{\nu - 1}}{e^{x - y} + 1}.\end{aligned} \hspace{\stretch{1}}(1.41)

The integral has been split into two since the behavior of the exponential in the denominator is quite different in the $x y$ ranges. Observe that in the first integral we have

\begin{aligned}\frac{1}{{2}} \le \frac{1}{e^{x - y} + 1} \le \frac{1}{{1 + e^{-y}}}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

Since this term is of order 1, let’s consider the difference of this from $1$, writing

\begin{aligned}\frac{1}{e^{x - y} + 1} = 1 + u,\end{aligned} \hspace{\stretch{1}}(1.0.42)

or

\begin{aligned}u = \frac{1}{e^{x - y} + 1} - 1 &= \frac{1 -(e^{x - y} + 1)}{e^{x - y} + 1} \\ &= \frac{-e^{x - y} }{e^{x - y} + 1} \\ &= -\frac{1}{1 + e^{y - x}}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

This gives us

\begin{aligned}F_\nu(e^y) &= \int_0^y dx x^{\nu - 1} \left( 1 - \frac{ 1 } { 1 + e^{y - x} } \right)+\int_y^\infty dx \frac{x^{\nu - 1}}{e^{x - y} + 1} \\ &= \frac{y^\nu}{\nu}-\int_0^y dx \frac{ x^{\nu - 1} } { 1 + e^{y - x} }+\int_y^\infty dx \frac{x^{\nu - 1}}{e^{x - y} + 1}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

Now let’s make a change of variables $a = y - x$ in the first integral and $b = x - y$ in the second. This gives

\begin{aligned}F_\nu(e^y) = \frac{y^\nu}{\nu}-\int_0^\infty da \frac{ (y - a)^{\nu - 1} } { 1 + e^{a} }+\int_0^\infty db \frac{(y + b)^{\nu - 1}}{e^{b} + 1}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

As $a$ gets large in the first integral the integrand is approximately $e^{-a} (y-a)^{\nu -1}$. The exponential dominates this integrand. Since we are considering large $y$, we can approximate the upper bound of the integral by extending it to $\infty$. Also expanding in series we have

\begin{aligned}F_\nu(e^y) &\approx \frac{y^\nu}{\nu}+\int_0^\infty da \frac{ (y + a)^{\nu - 1} -(y - a)^{\nu - 1} } { 1 + e^{a} } \\ &= \frac{y^\nu}{\nu}+\int_0^\infty da \frac{1}{{e^a + 1}}\left( \left( \frac{1}{{0!}} y^{\nu-1} a^0 + \frac{1}{{1!}} (\nu-1) y^{\nu-2} a^1 + \frac{1}{{2!}} (\nu-1) (\nu-2) y^{\nu-3} a^2 + \cdots \right) - \left( \frac{1}{{0!}} y^{\nu-1} (-a)^0 + \frac{1}{{1!}} (\nu-1) y^{\nu-2} (-a)^1 + \frac{1}{{2!}} (\nu-1) (\nu-2) y^{\nu-3} (-a)^2 + \cdots \right) \right) \\ &= \frac{y^\nu}{\nu}+ 2\int_0^\infty da \frac{1}{{e^a + 1}} \left( \frac{1}{{1!}} (\nu-1) y^{\nu-2} a^1 + \frac{1}{{3!}} (\nu-1) (\nu-2) (\nu - 3)y^{\nu-4} a^3 + \cdots \right) \\ &= \frac{y^\nu}{\nu}+ 2\sum_{j = 1, 3, 5, \cdots} \frac{y^{\nu - 1 - j}}{j!} \left( \prod_{k = 1}^j (\nu-k) \right)\int_0^\infty da \frac{a^j}{e^a + 1} \\ &= \frac{y^\nu}{\nu}+ 2\sum_{j = 1, 3, 5, \cdots} \frac{y^{\nu - 1 - j}}{j!} \frac{ \Gamma(\nu) } {\Gamma(\nu - j)}\int_0^\infty da \frac{a^j}{e^a + 1}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

For the remaining integral, Mathematica gives

\begin{aligned}\int_0^\infty da \frac{a^j}{e^a + 1}=\left( 1 - 2^{-j} \right) j! \zeta (j+1),\end{aligned} \hspace{\stretch{1}}(1.0.42)

where for $s > 1$

\begin{aligned}\zeta(s) = \sum_{k=1}^{\infty} k^{-s}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

This gives

\begin{aligned}F_\nu(e^y) \approx \frac{y^\nu}{\nu}+ 2\sum_{j = 1, 3, 5, \cdots} y^{\nu - 1 - j}\frac{ \Gamma(\nu) } {\Gamma(\nu - j)}\left( 1 - 2^{-j} \right) \zeta(j+1),\end{aligned} \hspace{\stretch{1}}(1.0.42)

or

\begin{aligned}f_\nu^+(e^y) &\approx y^\nu\left( \frac{1}{\nu \Gamma(\nu)} + 2 \sum_{j = 1, 3, 5, \cdots} \frac{ 1 } {\Gamma(\nu - j) y^{j + 1} } \left( 1 - 2^{-j} \right) \zeta(j+1) \right) \\ &= \frac{y^\nu}{\Gamma(\nu + 1)}\left( 1 + 2 \sum_{j = 1, 3, 5, \cdots} \frac{ \Gamma(\nu + 1) } {\Gamma(\nu - j) } \left( 1 - 2^{-j} \right) \frac{\zeta(j+1)}{ y^{j + 1} } \right),\end{aligned} \hspace{\stretch{1}}(1.0.42)

or

\begin{aligned}\boxed{f_\nu^+(e^y) \approx \frac{y^\nu}{\Gamma(\nu + 1)}\left( 1 + 2 \nu \sum_{j = 1, 3, 5, \cdots} (\nu-1) \cdots(\nu - j) \left( 1 - 2^{-j} \right) \frac{\zeta(j+1)}{ y^{j + 1} } \right).}\end{aligned} \hspace{\stretch{1}}(1.0.42)

Evaluating the numerical portions explicitly, with

\begin{aligned}c(j) = 2 \left(1-2^{-j}\right) \zeta (j+1),\end{aligned} \hspace{\stretch{1}}(1.0.42)

\begin{aligned}\begin{aligned}c(1) &= \frac{\pi^2}{6} \\ c(3) &= \frac{7 \pi^4}{360} \\ c(5) &= \frac{31 \pi^6}{15120} \\ c(7) &= \frac{127 \pi^8}{604800},\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.42)

so to two terms ($j = 1, 3$), we have

\begin{aligned}\boxed{f_\nu^+(e^y) \approx \frac{y^\nu}{\Gamma(\nu + 1)}\left( 1 + \nu(\nu-1) \frac{\pi^2}{6 y^{2}} + \nu(\nu-1)(\nu-2)(\nu -3) \frac{7 \pi^4}{360 y^4} \right).}\end{aligned} \hspace{\stretch{1}}(1.0.42)

In order for the Boson occupation numbers to be non-singular we require $\mu$ less than all $\epsilon$. If that lowest energy level is set to zero, this is equivalent to $z < 1$. Given this restriction, a $z = e^{-\alpha}$ substitution is convenient for investigation of the $z \rightarrow 1$ case. Following the text, we'll write

\begin{aligned}G_\nu(e^{-\alpha})=\Gamma(\nu)f_\nu^-(e^{-\alpha}) = \int_0^\infty dx \frac{x^{\nu - 1}}{e^{x + \alpha} - 1}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

For $\nu = 1$, this is integrable

\begin{aligned}\frac{d}{dx} \ln\left( 1 - e^{-x - \alpha} \right) &= \frac{e^{-x - \alpha}}{ 1 - e^{-x - \alpha} } \\ &= \frac{1}{{ e^{x + \alpha} - 1}},\end{aligned} \hspace{\stretch{1}}(1.0.42)

so that

\begin{aligned}G_1(e^{-\alpha}) &= \int_0^\infty dx \frac{1}{e^{x + \alpha} - 1} \\ &= {\ln \left( 1 - e^{-x - \alpha} \right)}_{0}^{\infty} \\ &= \ln 1 - \ln \left( 1 - e^{- \alpha} \right) \\ &= -\ln \left( 1 - e^{- \alpha} \right).\end{aligned} \hspace{\stretch{1}}(1.0.42)

Taylor expanding $1 - e^{-\alpha}$ we have

\begin{aligned}1 - e^{-\alpha} = 1 - \left( 1 - \alpha + \alpha^2/2 - \cdots \right).\end{aligned} \hspace{\stretch{1}}(1.0.42)

Noting that $\Gamma(1) = 1$, we have for the limit

\begin{aligned}\lim_{\alpha \rightarrow 0} G_1(e^{-\alpha}) \rightarrow - \ln \alpha,\end{aligned} \hspace{\stretch{1}}(1.0.42)

or

\begin{aligned}\lim_{z\rightarrow 1} f_\nu^-(z)= -\ln (-\ln z).\end{aligned} \hspace{\stretch{1}}(1.0.42)

For values of $\nu \ne 1$, the denominator is

\begin{aligned}e^{\alpha + x} - 1 = (\alpha + x) + (\alpha + x)^2/2 + \cdots\end{aligned} \hspace{\stretch{1}}(1.0.42)

To first order this gives us

\begin{aligned}f_\nu^-( e^{-\alpha} ) \approx \frac{1}{{\Gamma(\nu)}} \int_0^\infty dx \frac{1}{x + \alpha}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

Of this integral Mathematica says it can be evaluated for $0 < \nu < 1$, and has the value

\begin{aligned}\frac{1}{{\Gamma(\nu)}} \int_0^\infty dx \frac{1}{x + \alpha}=\frac{\pi}{\sin(\pi\nu)} \frac{1}{\alpha^{1 - \nu} \Gamma (\nu )}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

From [1] 6.1.17 we find

\begin{aligned}\Gamma(z) \Gamma(1-z) = \frac{\pi}{\sin(\pi z)},\end{aligned} \hspace{\stretch{1}}(1.0.42)

with which we can write

\begin{aligned}\boxed{f_\nu^-( e^{-\alpha} ) \approx \frac{ \Gamma(1 - \nu)}{ \alpha^{1 - \nu} }.}\end{aligned} \hspace{\stretch{1}}(1.0.42)

## Question: Nuclear matter ([2], prob 9.2)

Consider a heavy nucleus of mass number $A$. i.e., having $A$ total nucleons including neutrons and protons. Assume that the number of neutrons and protons is equal, and recall that each of them has spin-$1/2$ (so possessing two spin states). Treating these nucleons as a free ideal Fermi gas of uniform density contained in a radius $R = r_0 A^{1/3}$, where $r_0 = 1.4 \times 10^{-13} \text{cm}$, calculate the Fermi energy and the average energy per nucleon in MeV.

Our nucleon particle density is

\begin{aligned}\rho &= \frac{N}{V} \\ &= \frac{A}{\frac{4 \pi}{3} R^3} \\ &= \frac{3 A}{4 \pi r_0^3 A} \\ &= \frac{3}{4 \pi r_0^3} \\ &= \frac{3}{4 \pi (1.4 \times 10^{-13} \text{cm})^3} \\ &= 8.7 \times 10^{37} (\text{cm})^{-3} \\ &= 8.7 \times 10^{43} (\text{m})^{-3}\end{aligned} \hspace{\stretch{1}}(1.67)

With $m$ for the mass of either the proton or the neutron, and $\rho_m = \rho_p = \rho/2$, the Fermi energy for these particles is

\begin{aligned}\epsilon_{\mathrm{F}} = \frac{\hbar^2}{2m} \left( \frac{6 \pi (\rho/2)}{2 S + 1} \right)^{2/3},\end{aligned} \hspace{\stretch{1}}(1.68)

With $S = 1/2$, and $2 S + 1 = 2(1/2) + 1 = 2$ for either the proton or the neutron, this is

\begin{aligned}\epsilon_{\mathrm{F}} = \frac{\hbar^2}{2 m} \left( \frac{3 \pi^2 \rho}{2} \right)^{2/3}.\end{aligned} \hspace{\stretch{1}}(1.69)

\begin{aligned}\begin{aligned}\hbar &= 1.05 \times 10^{-34} \,\text{m^2 kg s^{-1}} \\ m &= 1.67 \times 10^{-27} \,\text{kg}\end{aligned}.\end{aligned} \hspace{\stretch{1}}(1.70)

This gives us

\begin{aligned}\epsilon_{\mathrm{F}} &= \frac{(1.05 \times 10^{-34})^2}{2 \times 1.67 \times 10^{-27}} \left( \frac{3 \pi^2 }{2} \frac{8.7 \times 10^{43} }{2} \right)^{2/3}\text{m}^4 \frac{\text{kg}^2}{s^2} \frac{1}{{\text{kg}}} \frac{1}{{\text{m}^2}} \\ &= 3.9 \times 10^{-12} \,\text{J} \times \left( 6.241509 \times 10^{12} \frac{\text{MeV}}{J} \right) \approx 24 \text{MeV}\end{aligned} \hspace{\stretch{1}}(1.71)

In lecture 16

we found that the total average energy for a Fermion gas of $N$ particles was

\begin{aligned}E = \frac{3}{5} N \epsilon_{\mathrm{F}},\end{aligned} \hspace{\stretch{1}}(1.72)

so the average energy per nucleon is approximately

\begin{aligned}\frac{3}{5} \epsilon_{\mathrm{F}} \approx 15 \,\text{MeV}.\end{aligned} \hspace{\stretch{1}}(1.73)

## Question: Neutron star ([2], prob 9.5)

Model a neutron star as an ideal Fermi gas of neutrons at $T = 0$ moving in the gravitational field of a heavy point mass $M$ at the center. Show that the pressure $P$ obeys the equation

\begin{aligned}\frac{dP}{dr} = - \gamma M \frac{\rho(r)}{r^2},\end{aligned} \hspace{\stretch{1}}(1.74)

where $\gamma$ is the gravitational constant, $r$ is the distance from the center, and $\rho(r)$ is the density which only depends on distance from the center.

In the grand canonical scheme the pressure for a Fermion system is given by

\begin{aligned}\beta P V &= \ln Z_{\mathrm{G}} \\ &= \ln \prod_\epsilon \sum_{n = 0}^1 \left( z e^{-\beta \epsilon} \right)^n \\ &= \sum_\epsilon \ln \left( 1 + z e^{-\beta \epsilon} \right).\end{aligned} \hspace{\stretch{1}}(1.75)

The kinetic energy of the particle is adjusted by the gravitational potential

\begin{aligned}\epsilon &= \epsilon_\mathbf{k}- \frac{\gamma m M}{r} \\ &= \frac{\hbar^2 \mathbf{k}^2}{2m}- \frac{\gamma m M}{r}.\end{aligned} \hspace{\stretch{1}}(1.76)

Differentiating eq. 1.75 with respect to the radius, we have

\begin{aligned}\beta V \frac{\partial {P}}{\partial {r}} &= -\beta \frac{\partial {\epsilon}}{\partial {r}}\sum_\epsilon \frac{z e^{-\beta \epsilon}}{ 1 + z e^{-\beta \epsilon} } \\ &= -\beta\left( \frac{\gamma m M}{r^2} \right)\sum_\epsilon \frac{1}{ z^{-1} e^{\beta \epsilon} + 1} \\ &= -\beta\left( \frac{\gamma m M}{r^2} \right)\left\langle{{N}}\right\rangle.\end{aligned} \hspace{\stretch{1}}(1.77)

Noting that $\left\langle{{N}}\right\rangle m/V$ is the average density of the particles, presumed radial, we have

\begin{aligned}\boxed{\frac{\partial P}{\partial r} &= -\frac{\gamma M}{r^2} \frac{m \left\langle N \right\rangle}{V} \\ &= -\frac{\gamma M}{r^2} \rho(r).}\end{aligned} \hspace{\stretch{1}}(1.0.78)

# References

[1] M. Abramowitz and I.A. Stegun. \emph{Handbook of mathematical functions with formulas, graphs, and mathematical tables}, volume 55. Dover publications, 1964.

[2] Kerson Huang. Introduction to statistical physics. CRC Press, 2001.

[3] RK Pathria. Statistical mechanics. Butterworth Heinemann, Oxford, UK, 1996.

## An updated compilation of notes, for ‘PHY452H1S Basic Statistical Mechanics’, Taught by Prof. Arun Paramekanti

Posted by peeterjoot on March 27, 2013

Here’s my second update of my notes compilation for this course, including all of the following:

March 27, 2013 Fermi gas

March 26, 2013 Fermi gas thermodynamics

March 26, 2013 Fermi gas thermodynamics

March 23, 2013 Relativisitic generalization of statistical mechanics

March 21, 2013 Kittel Zipper problem

March 18, 2013 Pathria chapter 4 diatomic molecule problem

March 17, 2013 Gibbs sum for a two level system

March 16, 2013 open system variance of N

March 16, 2013 probability forms of entropy

March 14, 2013 Grand Canonical/Fermion-Bosons

March 13, 2013 Quantum anharmonic oscillator

March 12, 2013 Grand canonical ensemble

March 11, 2013 Heat capacity of perturbed harmonic oscillator

March 10, 2013 Langevin small approximation

March 10, 2013 Addition of two one half spins

March 10, 2013 Midterm II reflection

March 07, 2013 Thermodynamic identities

March 06, 2013 Temperature

March 05, 2013 Interacting spin

plus everything detailed in the description of my first update and before.

## PHY452H1S Basic Statistical Mechanics. Lecture 16: Fermi gas. Taught by Prof. Arun Paramekanti

Posted by peeterjoot on March 27, 2013

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

# Disclaimer

Peeter’s lecture notes from class. May not be entirely coherent.

# Fermi gas

Review

Continuing a discussion of [1] section 8.1 content.

We found

\begin{aligned}n_{\mathbf{k}} = \frac{1}{{e^{\beta(\epsilon_k - \mu)} + 1}}\end{aligned} \hspace{\stretch{1}}(1.2.1)

With no spin

\begin{aligned}\int n_\mathbf{k} \times \frac{d^3 k}{(2\pi)^3} = \rho\end{aligned} \hspace{\stretch{1}}(1.2.2)

Fig 1.1: Occupancy at low temperature limit

Fig 1.2: Volume integral over momentum up to Fermi energy limit

\begin{aligned}\epsilon_{\mathrm{F}} = \frac{\hbar^2 k_{\mathrm{F}}^2}{2m}\end{aligned} \hspace{\stretch{1}}(1.2.3)

gives

\begin{aligned}k_{\mathrm{F}} = (6 \pi^2 \rho)^{1/3}\end{aligned} \hspace{\stretch{1}}(1.2.4)

\begin{aligned}\sum_\mathbf{k} n_\mathbf{k} = N\end{aligned} \hspace{\stretch{1}}(1.2.5)

\begin{aligned}\mathbf{k} = \frac{2\pi}{L}(n_x, n_y, n_z)\end{aligned} \hspace{\stretch{1}}(1.2.6)

This is for periodic boundary conditions \footnote{I filled in details in the last lecture using a particle in a box, whereas this periodic condition was intended. We see that both achieve the same result}, where

\begin{aligned}\Psi(x + L) = \Psi(x)\end{aligned} \hspace{\stretch{1}}(1.2.7)

Moving on

\begin{aligned}\sum_{k_x} n(\mathbf{k}) = \sum_{p_x} \Delta p_x n(\mathbf{k})\end{aligned} \hspace{\stretch{1}}(1.2.8)

with

\begin{aligned}\Delta k_x = \frac{2 \pi}{L} \Delta p_x\end{aligned} \hspace{\stretch{1}}(1.2.9)

this gives

\begin{aligned}\sum_{k_x} n(\mathbf{k}) = \sum_{n_x} \frac{L}{2\pi} \Delta k_x \rightarrow \frac{L}{2\pi} \int d k_x\end{aligned} \hspace{\stretch{1}}(1.2.10)

Over all dimensions

\begin{aligned}\sum_{\mathbf{k}} n_\mathbf{k} = \left( \frac{L}{2\pi} \right)^3 \left( \int d^3 \mathbf{k} \right)n(\mathbf{k})=N\end{aligned} \hspace{\stretch{1}}(1.2.11)

so that

\begin{aligned}\rho = \int \frac{d^3 \mathbf{k}}{(2 \pi)^3}\end{aligned} \hspace{\stretch{1}}(1.2.12)

Again

\begin{aligned}k_{\mathrm{F}} = (6 \pi^2 \rho)^{1/3}\end{aligned} \hspace{\stretch{1}}(1.2.13)

## Example: Spin considerations

{example:basicStatMechLecture16:1}{

\begin{aligned}\sum_{\mathbf{k}, m_s} = N\end{aligned} \hspace{\stretch{1}}(1.2.14)

\begin{aligned}\sum_{\mathbf{k}, m_s} \frac{1}{{e^{\beta(\epsilon_k - \mu)} + 1}} = (2 S + 1)\left( \int \frac{d^3 \mathbf{k}}{(2 \pi)^3} n(\mathbf{k}) \right)L^3\end{aligned} \hspace{\stretch{1}}(1.2.15)

This gives us

\begin{aligned}k_{\mathrm{F}} = \left( \frac{ 6 \pi^2 \rho }{2 S + 1} \right)^{1/3}\end{aligned} \hspace{\stretch{1}}(1.2.16)

and again

\begin{aligned}\epsilon_{\mathrm{F}} = \frac{\hbar^2 k_{\mathrm{F}}^2}{2m}\end{aligned} \hspace{\stretch{1}}(1.2.17)

}

High Temperatures

Now we want to look at the at higher temperature range, where the occupancy may look like fig. 1.3

Fig 1.3: Occupancy at higher temperatures

\begin{aligned}\mu(T = 0) = \epsilon_{\mathrm{F}}\end{aligned} \hspace{\stretch{1}}(1.2.18)

\begin{aligned}\mu(T \rightarrow \infty) \rightarrow - \infty\end{aligned} \hspace{\stretch{1}}(1.2.19)

so that for large $T$ we have

\begin{aligned}\frac{1}{{e^{\beta(\epsilon_k - \mu)} + 1}} \rightarrow e^{-\beta(\epsilon_k - \mu)}\end{aligned} \hspace{\stretch{1}}(1.2.20)

\begin{aligned}\rho &= \int \frac{d^3 \mathbf{k}}{(2 \pi)^3} e^{\beta \mu} e^{-\beta \epsilon_k} \\ &= e^{\beta \mu} \int \frac{d^3 \mathbf{k}}{(2 \pi)^3} e^{-\beta \epsilon_k} \\ &= e^{\beta \mu} \int dk \frac{4 \pi k^2}{(2 \pi)^3} e^{-\beta \hbar^2 k^2/2m}.\end{aligned} \hspace{\stretch{1}}(1.2.21)

Mathematica (or integration by parts) tells us that

\begin{aligned}\frac{1}{{(2 \pi)^3}} \int 4 \pi^2 k^2 dk e^{-a k^2} = \frac{1}{{(4 \pi a )^{3/2}}},\end{aligned} \hspace{\stretch{1}}(1.2.22)

so we have

\begin{aligned}\rho &= e^{\beta \mu} \left( \frac{2m}{ 4 \pi \beta \hbar^2} \right)^{3/2} \\ &= e^{\beta \mu} \left( \frac{2 m k_{\mathrm{B}} T 4 \pi^2 }{ 4 \pi h^2} \right)^{3/2} \\ &= e^{\beta \mu} \left( \frac{2 m k_{\mathrm{B}} T \pi }{ h^2} \right)^{3/2}\end{aligned} \hspace{\stretch{1}}(1.2.23)

Introducing $\lambda$ for the thermal de Broglie wavelength, $\lambda^3 \sim T^{-3/2}$

\begin{aligned}\lambda \equiv \frac{h}{\sqrt{2 \pi m k_{\mathrm{B}} T}},\end{aligned} \hspace{\stretch{1}}(1.2.24)

we have

\begin{aligned}\rho = e^{\beta \mu} \frac{1}{{\lambda^3}}.\end{aligned} \hspace{\stretch{1}}(1.2.25)

Does it make any sense to have density as a function of temperature? An inappropriately extended to low temperatures plot of the density is found in fig. 1.4 for a few arbitrarily chosen numerical values of the chemical potential $\mu$, where we see that it drops to zero with temperature. I suppose that makes sense if we are not holding volume constant.

Fig 1.4: Density as a function of temperature

We can write

\begin{aligned}\boxed{e^{\beta \mu} = \left( \rho \lambda^3 \right)}\end{aligned} \hspace{\stretch{1}}(1.2.26)

\begin{aligned}\frac{\mu}{k_{\mathrm{B}} T} = \ln \left( \rho \lambda^3 \right)\sim -\frac{3}{2} \ln T\end{aligned} \hspace{\stretch{1}}(1.2.27)

or (taking $\rho$ (and/or volume?) as a constant) we have for large temperatures

\begin{aligned}\mu \propto -T \ln T\end{aligned} \hspace{\stretch{1}}(1.2.28)

The chemical potential is plotted in fig. 1.5, whereas this $- k_{\mathrm{B}} T \ln k_{\mathrm{B}} T$ function is plotted in fig. 1.6. The contributions to $\mu$ from the $k_{\mathrm{B}} T \ln (\rho h^3 (2 \pi m)^{-3/2})$ term are dropped for the high temperature approximation.

Fig 1.5: Chemical potential over degenerate to classical range

Fig 1.6: High temp approximation of chemical potential, extended back to T = 0

Pressure

\begin{aligned}P = - \frac{\partial {E}}{\partial {V}}\end{aligned} \hspace{\stretch{1}}(1.2.29)

For a classical ideal gas as in fig. 1.7 we have

Fig 1.7: Ideal gas pressure vs volume

\begin{aligned}P = \rho k_{\mathrm{B}} T\end{aligned} \hspace{\stretch{1}}(1.2.30)

For a Fermi gas at $T = 0$ we have

\begin{aligned}E &= \sum_\mathbf{k} \epsilon_k n_k \\ &= \sum_\mathbf{k} \epsilon_k \Theta(\mu_0 - \epsilon_k) \\ &= \frac{V}{(2\pi)^3} \int_{\epsilon_k < \mu_0} \frac{\hbar^2 \mathbf{k}^2}{2 m} d^3 \mathbf{k} \\ &= \frac{V}{(2\pi)^3} \int_0^{k_{\mathrm{F}}} \frac{\hbar^2 \mathbf{k}^2}{2 m} d^3 \mathbf{k} \\ &= \frac{V}{(2\pi)^3} \frac{\hbar^2}{2 m} \int_0^{k_{\mathrm{F}}} k^2 4 \pi k^2 d k\propto k_{\mathrm{F}}^5\end{aligned} \hspace{\stretch{1}}(1.2.31)

Specifically,

\begin{aligned}E(T = 0) = V \times \frac{3}{5} \underbrace{\epsilon_{\mathrm{F}}}_{\sim k_{\mathrm{F}}^2}\underbrace{\rho}_{\sim k_{\mathrm{F}}^3}\end{aligned} \hspace{\stretch{1}}(1.2.32)

or

\begin{aligned}\frac{E}{N} = \frac{3}{5} \epsilon_{\mathrm{F}}\end{aligned} \hspace{\stretch{1}}(1.2.33)

\begin{aligned}E = \frac{3}{5} N \frac{\hbar^2}{2 m} \left( 6 \pi^2 \frac{N}{V} \right)^{2/3} = a V^{-2/3},\end{aligned} \hspace{\stretch{1}}(1.2.34)

so that

\begin{aligned}\frac{\partial {E}}{\partial {V}} = -\frac{2}{3} a V^{-5/3}.\end{aligned} \hspace{\stretch{1}}(1.2.35)

\begin{aligned}P &= -\frac{\partial {E}}{\partial {V}} \\ &= \frac{2}{3} \left( a V^{-2/3} \right)V^{-1} \\ &= \frac{2}{3} \frac{E}{V} \\ &= \frac{2}{3} \left( \frac{3}{5} \epsilon_{\mathrm{F}} \rho \right) \\ &= \frac{2}{5} \epsilon_{\mathrm{F}} \rho.\end{aligned} \hspace{\stretch{1}}(1.2.36)

We see that the pressure ends up deviating from the classical result at low temperatures, as sketched in fig. 1.8. This low temperature limit for the pressure $2 \epsilon_{\mathrm{F}} \rho/5$ is called the degeneracy pressure.

Fig 1.8: Fermi degeneracy pressure

# References

[1] RK Pathria. Statistical mechanics. Butterworth Heinemann, Oxford, UK, 1996.

## PHY452H1S Basic Statistical Mechanics. Lecture 15: Grand Canonical/Fermion-Bosons. Taught by Prof. Arun Paramekanti

Posted by peeterjoot on March 14, 2013

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

# Disclaimer

Peeter’s lecture notes from class. May not be entirely coherent.

# Grand Canonical/Fermion-Bosons

Was mentioned that three dimensions confines us to looking at either Fermions or Bosons, and that two dimensions is a rich subject (interchange of two particles isn’t the same as one particle cycling around the other ending up in the same place — how is that different than a particle cycling around another in a two dimensional space?)

Definitions

1. Fermion. Antisymmetric under exchange. $n_k = 0, 1$
2. Boson. Symmetric under exchange. $n_k = 0, 1, 2, \cdots$

In either case our energies are

\begin{aligned}\epsilon_k = \frac{\hbar^2 k^2}{2m},\end{aligned} \hspace{\stretch{1}}(1.2.1)

For Fermions we’ll have occupation filling of the form fig. 1.1, where there can be only one particle at any given site (an energy level for that value of momentum). For Bosonic systems as in fig. 1.2, we don’t have a restriction of only one particle for each state, and can have any given number of particles for each value of momentum.

Fig 1.1: Fermionic energy level filling for free particle in a box

Fig 1.2: Bosonic free particle in a box energy level filling

Our Hamiltonian is

\begin{aligned}H = \sum_k \hat{n}_k \epsilon_k,\end{aligned} \hspace{\stretch{1}}(1.2.2)

where we have a number operator

\begin{aligned}N = \sum \hat{n}_k,\end{aligned} \hspace{\stretch{1}}(1.2.3)

such that

\begin{aligned}\left[{N},{H}\right] = 0.\end{aligned} \hspace{\stretch{1}}(1.2.4)

\begin{aligned}Z_{\mathrm{G}} = \sum_{N=0}^\infty e^{\beta \mu N}\sum_{n_k, \sum n_k = N} e^{-\beta \sum_k n_k \epsilon_k}.\end{aligned} \hspace{\stretch{1}}(1.2.5)

While the second sum is constrained, because we are summing over all $n_k$, this is essentially an unconstrained sum, so we can write

\begin{aligned}Z_{\mathrm{G}} &= \sum_{n_k}e^{\beta \mu \sum_k n_k}e^{-\beta \sum_k n_k \epsilon_k} \\ &= \sum_{n_k} \left( \prod_k e^{-\beta(\epsilon_k - \mu) n_k} \right) \\ &= \prod_{n} \left( \sum_{n_k} e^{-\beta(\epsilon_k - \mu) n_k} \right).\end{aligned} \hspace{\stretch{1}}(1.2.6)

Fermions

\begin{aligned}\sum_{n_k = 0}^1 e^{-\beta(\epsilon_k - \mu) n_k} = 1 + e^{-\beta(\epsilon_k - \mu)}\end{aligned} \hspace{\stretch{1}}(1.2.7)

Bosons

\begin{aligned}\sum_{n_k = 0}^\infty e^{-\beta(\epsilon_k - \mu) n_k} = \frac{1}{{1 - e^{-\beta(\epsilon_k - \mu)}}}\end{aligned} \hspace{\stretch{1}}(1.2.8)

($\epsilon_k - \mu \ge 0$).

Our grand partition functions are then

\begin{aligned}Z_{\mathrm{G}}^f = \prod_k \left( 1 + e^{-\beta(\epsilon_k - \mu)} \right)\end{aligned} \hspace{\stretch{1}}(1.0.9a)

\begin{aligned}Z_{\mathrm{G}}^b = \prod_k \frac{1}{{ 1 - e^{-\beta(\epsilon_k - \mu)} }}\end{aligned} \hspace{\stretch{1}}(1.0.9b)

We can use these to compute the average number of particles

\begin{aligned}\left\langle{{n_k^f}}\right\rangle = \frac{1 \times 0 + e^{-\beta(\epsilon_k - \mu)} \times 1}{ 1 + e^{-\beta(\epsilon_k - \mu)} }=\frac{1}{{ 1 + e^{-\beta(\epsilon_k - \mu)} }}\end{aligned} \hspace{\stretch{1}}(1.0.10)

\begin{aligned}\left\langle{{n_k^b}}\right\rangle = \frac{1 \times 0 + e^{-\beta(\epsilon_k - \mu)} \times 1+e^{-2 \beta(\epsilon_k - \mu)} \times 2 + \cdots}{ 1+e^{-\beta(\epsilon_k - \mu)} +e^{-2 \beta(\epsilon_k - \mu)} }\end{aligned} \hspace{\stretch{1}}(1.0.11)

This chemical potential over temperature exponential

\begin{aligned}e^{\beta \mu} \equiv z,\end{aligned} \hspace{\stretch{1}}(1.0.12)

is called the fugacity. The denominator has the form

\begin{aligned}D = 1 + z e^{-\beta \epsilon_k}+ z^2 e^{-2 \beta \epsilon_k},\end{aligned} \hspace{\stretch{1}}(1.0.13)

so we see that

\begin{aligned}z \frac{\partial {D}}{\partial {z}} = z e^{-\beta \epsilon_k}+ 2 z^2 e^{-2 \beta \epsilon_k}+ 3 z^3 e^{-3 \beta \epsilon_k}+ \cdots\end{aligned} \hspace{\stretch{1}}(1.0.14)

Thus the numerator is

\begin{aligned}N = z \frac{\partial {D}}{\partial {z}},\end{aligned} \hspace{\stretch{1}}(1.0.15)

and

\begin{aligned}\left\langle{{n_k^b}}\right\rangle &= \frac{z \frac{\partial {D_k}}{\partial {z}} }{D_k} \\ &= z \frac{\partial {}}{\partial {z}} \ln D_k \\ &= \cdots \\ &= \frac{1}{{ e^{\beta(\epsilon_k - \mu)} - 1}}\end{aligned} \hspace{\stretch{1}}(1.0.16)

What is the density $\rho$?

For Fermions

\begin{aligned}\rho = \frac{N}{V} =\frac{1}{{V}} \sum_{\mathbf{k}}\frac{1}{{ e^{\beta(\epsilon_\mathbf{k} - \mu)} + 1}}\end{aligned} \hspace{\stretch{1}}(1.0.17)

Using a “particle in a box” quantization where $k_\alpha = 2 \pi m_\alpha/L$, in a $d$-dimensional space, we can approximate this as

\begin{aligned}\boxed{\rho = \int \frac{d^d k}{(2 \pi)^d}\frac{1}{{ e^{\beta(\epsilon_k - \mu)} - 1}}.}\end{aligned} \hspace{\stretch{1}}(1.0.18)

This integral is actually difficult to evaluate. For $T \rightarrow 0$ ($\beta \rightarrow \infty$, where

\begin{aligned}n_k = \Theta(\mu - \epsilon_k).\end{aligned} \hspace{\stretch{1}}(1.0.19)

This is illustrated in, where we also show the smearing that occurs as temperature increases fig. 1.3.

Fig 1.3: Occupation numbers for different energies

With

\begin{aligned}E_{\mathrm{F}} = \mu(T = 0),\end{aligned} \hspace{\stretch{1}}(1.0.20)

we want to ask what is the radius of the ball for which

\begin{aligned}\epsilon_k = E_{\mathrm{F}}\end{aligned} \hspace{\stretch{1}}(1.0.21)

or

\begin{aligned}E_{\mathrm{F}} = \frac{\hbar^2 k_{\mathrm{F}}^2}{2m},\end{aligned} \hspace{\stretch{1}}(1.0.22)

so that

\begin{aligned}k_{\mathrm{F}} = \sqrt{\frac{2 m E_{\mathrm{F}}}{\hbar^2}},\end{aligned} \hspace{\stretch{1}}(1.0.23)

so that our density where $\epsilon_k = \mu$ is

\begin{aligned}\rho &= \int_{k \le k_{\mathrm{F}}} \frac{d^3 k}{(2 \pi)^3} \times 1 \\ &= \frac{1}{{(2\pi)^3}} 4 \pi \int^{k_{\mathrm{F}}} k^2 dk \\ &= \frac{4 \pi}{3} k_{\mathrm{F}}^3 \frac{1}{{(2 \pi)^3}},\end{aligned} \hspace{\stretch{1}}(1.0.24)

so that

\begin{aligned}k_{\mathrm{F}} = (6 \pi^2 \rho)^{1/3},\end{aligned} \hspace{\stretch{1}}(1.0.25)

Our chemical potential at zero temperature is then

\begin{aligned}\mu(T = 0) = \frac{\hbar^2}{2m} (6 \pi^2 \rho)^{2/3}.\end{aligned} \hspace{\stretch{1}}(1.0.26)

\begin{aligned}\rho^{-1/3} = \mbox{interparticle spacing}.\end{aligned} \hspace{\stretch{1}}(1.0.27)

We can convince ourself that the chemical potential must have the form fig. 1.4.

Fig 1.4: Large negative chemical potential at high temperatures

Given large negative chemical potential at high temperatures our number distribution will have the form

\begin{aligned}\left\langle{{n_k}}\right\rangle = e^{-\beta (\epsilon_k - \mu)} \propto e^{-\beta \epsilon_k}\end{aligned} \hspace{\stretch{1}}(1.0.28)

We see that the classical Boltzmann distribution is recovered for high temperatures.

We can also calculate the chemical potential at high temperatures. We’ll find that this has the form

\begin{aligned}e^{\beta \mu} = \frac{4}{3} \rho \lambda_T^3,\end{aligned} \hspace{\stretch{1}}(1.0.29)

where this quantity $\lambda_T$ is called the Thermal de Broglie wavelength.

\begin{aligned}\lambda_T = \sqrt{\frac{ 2 \pi \hbar^2}{m k_{\mathrm{B}} T}}.\end{aligned} \hspace{\stretch{1}}(1.0.30)