• 336,988

# Posts Tagged ‘average energy’

## 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.

## Relativisitic Fermi gas

Posted by peeterjoot on April 20, 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.)]

## Question: Relativisitic Fermi gas ([1], pr 9.3)

Consider a relativisitic gas of $N$ particles of spin $1/2$ obeying Fermi statistics, enclosed in volume $V$, at absolute zero. The energy-momentum relation is

\begin{aligned}\epsilon = \sqrt{(p c)^2 + \epsilon_0^2},\end{aligned} \hspace{\stretch{1}}(1.1)

where $\epsilon_0 = m c^2$, and $m$ is the rest mass.

Find the Fermi energy at density $n$.

With the pressure $P$ defined as the average force per unit area exerted on a perfectly-reflecting wall of the container.

Set up expressions for this in the form of an integral.

Define the internal energy $U$ as the average $\epsilon - \epsilon_0$.
Set up expressions for this in the form of an integral.

Show that $P V = 2 U/3$ at low densities, and $P V = U/3$ at high densities. State the criteria for low and high densities.

There may exist a gas of neutrinos (and/or antineutrinos) in the cosmos. (Neutrinos are massless Fermions of spin $1/2$.) Calculate the Fermi energy (in eV) of such a gas, assuming a density of one particle per $\text{cm}^3$.

Attempt exact evaluation of the various integrals.

We’ve found [3] that the density of states associated with a 3D relativisitic system is

\begin{aligned}\mathcal{D}(\epsilon) = \frac{4 \pi V}{(c h)^3} \epsilon \sqrt{\epsilon^2 -\epsilon_0^2},\end{aligned} \hspace{\stretch{1}}(1.0.2)

For a given density $n$, we can find the Fermi energy in the same way as we did for the non-relativisitic energies, with the exception that we have to integrate from a lowest energy of $\epsilon_0$ instead of $0$ (the energy at $\mathbf{p} = 0$). That is

\begin{aligned}n &= \frac{N}{V} \\ &= \left( 2 \frac{1}{{2}} + 1 \right)\frac{4 \pi}{(c h)^3} \int_{\epsilon_0}^{\epsilon_{\mathrm{F}}}d\epsilon \epsilon \sqrt{ \epsilon^2 -\epsilon_0^2} \\ &= \frac{8 \pi}{(c h)^3}\frac{1}{{3}} {\left.{{\left( x^2 - \epsilon_0^2 \right)^{3/2}}}\right\vert}_{{\epsilon_0}}^{{\epsilon_{\mathrm{F}}}} \\ &= \frac{8 \pi}{3 (c h)^3}\left( \epsilon_{\mathrm{F}}^2 - \epsilon_0^2 \right)^{3/2}.\end{aligned} \hspace{\stretch{1}}(1.0.2)

Solving for $\epsilon_{\mathrm{F}}/\epsilon_0$ we have

\begin{aligned}\frac{\epsilon_{\mathrm{F}}}{\epsilon_0} =\sqrt{\left( \frac{3 (c h)^3 n}{8 \pi \epsilon_0^3} \right)^{2/3}+ 1}.\end{aligned} \hspace{\stretch{1}}(1.0.2)

We’ll see the constant factor above a number of times below and designate it

\begin{aligned}n_0 = \frac{8 \pi}{3} \left( \frac{\epsilon_0}{c h} \right)^3,\end{aligned} \hspace{\stretch{1}}(1.0.2)

so that the Fermi energy is

\begin{aligned}\frac{\epsilon_{\mathrm{F}}}{\epsilon_0} =\sqrt{\left( \frac{n}{n_0} \right)^{2/3}+ 1}.\end{aligned} \hspace{\stretch{1}}(1.0.2)

For the pressure calculation, let’s suppose that we have a configuration with a plane in the $x,y$ orientation as in fig. 1.1.

Fig 1.1: Pressure against x,y oriented plane

It’s argued in [4] section 6.4 that the pressure for such a configuration is

\begin{aligned}P = n \int p_z u_z f(\mathbf{u}) d^3 \mathbf{u},\end{aligned} \hspace{\stretch{1}}(1.7)

where $n$ is the number density and $f(\mathbf{u})$ is a normalized distribution function for the velocities. The velocity and momentum components are related by the Hamiltonian equations. From the Hamiltonian eq. 1.1 we find \footnote{ Observe that by squaring and summing one can show that this is equivalent to the standard relativisitic momentum $p_x = \frac{m v_x}{\sqrt{ 1 - \mathbf{u}^2/c^2}}$.} (for the x-component which is representative)

\begin{aligned}u_x \\ &= \frac{\partial {\epsilon}}{\partial {p_x}} \\ &= \frac{\partial {}}{\partial {p_x}}\sqrt{(p c)^2 +\epsilon_0^2} \\ &= \frac{ p_x c^2 }{\sqrt{(p c)^2 +\epsilon_0^2}}.\end{aligned} \hspace{\stretch{1}}(1.8)

For $\alpha \in \{1, 2, 3\}$ we can summarize these velocity-momentum relationships as

\begin{aligned}\frac{u_\alpha}{c} = \frac{ c p_\alpha }{ \epsilon }.\end{aligned} \hspace{\stretch{1}}(1.9)

Should we attempt to calculate the pressure with this parameterization of the velocity space we end up with convergence problems, and can’t express the results in terms of $f^+_\nu(z)$. Let’s try instead with a distribution over momentum space

\begin{aligned}P=n \int \frac{(c p_z)^2}{\epsilon} f(c \mathbf{p}) d^3 (c \mathbf{p}).\end{aligned} \hspace{\stretch{1}}(1.10)

Here the momenta have been scaled to have units of energy since we want to express this integral in terms of energy in the end. Our normalized distribution function is

\begin{aligned}f(c \mathbf{p})\propto \frac{\frac{1}{{ z^{-1} e^{\beta \epsilon} + 1 }}}{\int \frac{1}{{ z^{-1} e^{\beta \epsilon} + 1 }} d^3 (c \mathbf{p})},\end{aligned} \hspace{\stretch{1}}(1.11)

but before evaluating anything, we first want to change our integration variable from momentum to energy. In spherical coordinates our volume element takes the form

\begin{aligned}d^3 (c \mathbf{p}) &= 2 \pi (c p)^2 d (c p) \sin\theta d\theta \\ &= 2 \pi (c p)^2 \frac{d (c p)}{d \epsilon} d \epsilon \sin\theta d\theta.\end{aligned} \hspace{\stretch{1}}(1.12)

Implicit derivatives of

\begin{aligned}c^2 p^2 = \epsilon^2 - \epsilon_0^2,\end{aligned} \hspace{\stretch{1}}(1.13)

gives us

\begin{aligned}\frac{d (c p)}{d\epsilon}= \frac{\epsilon}{c p}=\frac{\epsilon}{\sqrt{\epsilon^2 -\epsilon_0^2}}.\end{aligned} \hspace{\stretch{1}}(1.0.14)

Our momentum volume element becomes

\begin{aligned}d^3 (c \mathbf{p}) \\ &= 2 \pi (c p)^2 \frac{\epsilon}{\sqrt{\epsilon^2 - \epsilon_0^2 }}d \epsilon \sin\theta d\theta \\ &= 2 \pi \left( \epsilon^2 - \epsilon_0^2 \right)\frac{\epsilon}{\sqrt{\epsilon^2 - \epsilon_0^2 }}d \epsilon \sin\theta d\theta \\ &= 2 \pi \epsilon \sqrt{ \epsilon^2 - \epsilon_0^2} d \epsilon \sin\theta d\theta.\end{aligned} \hspace{\stretch{1}}(1.0.14)

For our distribution function, we can now write

\begin{aligned}f(c \mathbf{p}) d^3 (c \mathbf{p})= C \frac{\epsilon \sqrt{ \epsilon^2 - \epsilon_0^2} d \epsilon }{ z^{-1} e^{\beta \epsilon} + 1 }\frac{ 2 \pi \sin\theta d\theta }{ 4 \pi \epsilon_0^3 },\end{aligned} \hspace{\stretch{1}}(1.0.14)

where $C$ is determined by the requirement $\int f(c \mathbf{p}) d^3 (c \mathbf{p}) = 1$

\begin{aligned}C^{-1} = \int_{0}^\infty \frac{(y + 1)\sqrt{ (y + 1)^2 - 1} dy }{ z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1 }.\end{aligned} \hspace{\stretch{1}}(1.0.14)

The z component of our momentum can be written in spherical coordinates as

\begin{aligned}(c p_z)^2= (c p)^2 \cos^2\theta= \left( \epsilon^2 - \epsilon_0^2 \right)\cos^2\theta,\end{aligned} \hspace{\stretch{1}}(1.0.18)

Noting that

\begin{aligned}\int_0^\pi \cos^2\theta \sin\theta d\theta =-\int_0^\pi \cos^2\theta d(\cos\theta)= \frac{2}{3},\end{aligned} \hspace{\stretch{1}}(1.0.19)

all the bits come together as

\begin{aligned}P &= \frac{C n}{3 \epsilon_0^3 } \int_{\epsilon_0}^\infty\left( \epsilon^2 - \epsilon_0^2 \right)^{3/2} \frac{1}{{ z^{-1} e^{\beta \epsilon} + 1 }} d \epsilon \\ &= \frac{n \epsilon_0}{3} \int_{1}^\infty\left( x^2 - 1 \right)^{3/2} \frac{1}{{ z^{-1} e^{\beta \epsilon_0 x} + 1 }} dx.\end{aligned} \hspace{\stretch{1}}(1.0.19)

Letting $y = x - 1$, this is

\begin{aligned}P= \frac{C n \epsilon_0}{3} \int_{0}^\infty \frac{ \left( (y + 1)^2 - 1 \right)^{3/2} } { z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1 } dy.\end{aligned} \hspace{\stretch{1}}(1.0.19)

We could conceivable expand the numerators of each of these integrals in power series, which could then be evaluated as a sum of $f^+_\nu(z e^{-\beta \epsilon_0})$ terms.

Note that above the Fermi energy $n$ also has an integral representation

\begin{aligned}n &= \left(2\left( \frac{1}{{2}} \right) + 1\right)\int_{\epsilon_0}^\infty d\epsilon \mathcal{D}(\epsilon) \frac{1}{{ z^{-1} e^{\beta \epsilon} + 1}} \\ &= \frac{8 \pi}{(c h)^3} \int_{\epsilon_0}^\infty d\epsilon\frac{\epsilon \sqrt{\epsilon^2 - \epsilon_0^2} }{ z^{-1} e^{\beta \epsilon} + 1} \\ &= \frac{8 \pi \epsilon_0^3}{(c h)^3} \int_{0}^\infty dy\frac{(y + 1)\sqrt{(y + 1)^2 - 1} }{ z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1},\end{aligned} \hspace{\stretch{1}}(1.0.19)

or

\begin{aligned}\boxed{n = \frac{3 n_0}{C}.}\end{aligned} \hspace{\stretch{1}}(1.0.23)

Observe that we can use this result to remove the dependence of pressure on this constant $C$

\begin{aligned}\boxed{\frac{P}{n_0 \epsilon_0}= \int_{0}^\infty dy \frac{ \left( (y + 1)^2 - 1 \right)^{3/2} } { z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1 }.}\end{aligned} \hspace{\stretch{1}}(1.0.24)

Now for the average energy difference from the rest energy $\epsilon_0$

\begin{aligned}U &= \left\langle{{\epsilon - \epsilon_0}}\right\rangle \\ &= \int_{\epsilon_0}^\infty d\epsilon \mathcal{D}(\epsilon) f(\epsilon) (\epsilon - \epsilon_0) \\ &= \frac{8 \pi V}{(c h)^3}\int_{\epsilon_0}^\infty d\epsilon \frac{ \epsilon(\epsilon - \epsilon_0) \sqrt{ \epsilon^2 - \epsilon_0 } }{ z^{-1} e^{\beta \epsilon} + 1} \\ &= \frac{8 \pi V \epsilon_0^4}{(c h)^3}\int_{0}^\infty dy\frac{ y ( y - 1 ) \sqrt{ (y + 1)^2 - 1 }}{ z^{-1} e^{\beta \epsilon} + 1}.\end{aligned} \hspace{\stretch{1}}(1.0.24)

So the average energy density difference from the rest energy, relative to the rest energy, is

\begin{aligned}\boxed{\frac{\left\langle{{\epsilon - \epsilon_0}}\right\rangle}{V \epsilon_0} =3 n_0 \int_{0}^\infty dy \frac { y (y + 1)\sqrt{(y + 1)^2 - 1} } { z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1 }.}\end{aligned} \hspace{\stretch{1}}(1.0.26)

From eq. 1.0.24 and eq. 1.0.26 we have

\begin{aligned}\frac{1}{{n_0}} &= 3 \frac{V \epsilon_0} {\left\langle{{\epsilon - \epsilon_0}}\right\rangle} \int_{0}^\infty \frac { y (y + 1)\sqrt{(y + 1)^2 - 1} dy } { z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1 } \\ &= \frac{\epsilon_0}{P} \int_{0}^\infty \frac{ \left( (y + 1)^2 - 1 \right)^{3/2} } { z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1 } dy,\end{aligned} \hspace{\stretch{1}}(1.0.26)

or

\begin{aligned}P V =\frac{U}{3}\frac{ \int_{0}^\infty \frac{ \left( (y + 1)^2 - 1 \right)^{3/2} } { z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1 } dy}{ \int_{0}^\infty \frac { y (y + 1)\sqrt{(y + 1)^2 - 1} dy } { z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1 }}.\end{aligned} \hspace{\stretch{1}}(1.0.26)

This ratio of integrals is supposed to resolve to 1 and 2 in the low and high density limits. To consider this let’s perform one final non-dimensionalization, writing

\begin{aligned}\begin{aligned} \\ x &= \beta \epsilon_0 y \\ \theta &= \frac{1}{{\beta \epsilon_0}} = \frac{k_{\mathrm{B}} T}{\epsilon_0} \\ \bar{\mu} &= \mu - \epsilon_0 \\ \bar{z} &= e^{\beta \bar{\mu}}.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.29)

The density, pressure, and energy take the form

\begin{aligned}\frac{n}{n_0}= 3 \theta\int_{0}^\infty dx\frac{(\theta x + 1)\sqrt{(\theta x + 1)^2 - 1} }{ \bar{z}^{-1} e^{x} + 1}\end{aligned} \hspace{\stretch{1}}(1.0.30a)

\begin{aligned}\frac{P}{n_0 \epsilon_0}= \theta \int_{0}^\infty dx \frac{ \left( (\theta x + 1)^2 - 1 \right)^{3/2} } { \bar{z}^{-1} e^{x} + 1 }\end{aligned} \hspace{\stretch{1}}(1.0.30b)

\begin{aligned}\frac{\left\langle{{\epsilon - \epsilon_0}}\right\rangle}{V \epsilon_0 n_0} =3 \theta^2 \int_{0}^\infty dx \frac { x (\theta x + 1)\sqrt{(\theta x + 1)^2 - 1} } { \bar{z}^{-1} e^{x} + 1 }.\end{aligned} \hspace{\stretch{1}}(1.0.30c)

We can rewrite the square roots in the number density and energy density expressions by expanding out the completion of the square

\begin{aligned}(1 + \theta x) \sqrt{ (1 + \theta x)^2 - 1}=(1 + \theta x) \sqrt{ 1 + \theta x + 1 }\sqrt{ 1 + \theta x - 1 }= \sqrt{2 \theta} x^{1/2} (1 + \theta x) \sqrt{ 1 + \frac{\theta x}{2}},\end{aligned} \hspace{\stretch{1}}(1.0.30c)

Expanding the distribution about $\bar{z} e^{-x} = 0$, we have

\begin{aligned}\frac{1}{ \bar{z}^{-1} e^{x} + 1}=\frac{\bar{z} e^{-x}}{ 1 + \bar{z} e^{-x}}=z e^{-x} \sum_{s = 0}^\infty (-1)^s \left( \bar{z} e^{-x} \right)^s,\end{aligned} \hspace{\stretch{1}}(1.0.32)

allowing us to write, in the low density limit with respect to $\bar{z}$

\begin{aligned}\frac{n}{n_0}= 3 \sqrt{2}\theta^{3/2} \sum_{s=0}^\infty(-1)^s\bar{z}^{s + 1}\int_{0}^\infty dx x^{1/2}(1 + \theta x) \sqrt{ 1 + \frac{\theta x}{2}} e^{-x(1 + s)} \end{aligned} \hspace{\stretch{1}}(1.0.33a)

\begin{aligned}\frac{P}{n_0 \epsilon_0}= \theta\sum_{s=0}^\infty(-1)^s\bar{z}^{s + 1} \int_{0}^\infty dx\left( (\theta x + 1)^2 - 1 \right)^{3/2} e^{-x(1 + s)} \end{aligned} \hspace{\stretch{1}}(1.0.33b)

\begin{aligned}\frac{\left\langle{{\epsilon - \epsilon_0}}\right\rangle}{V \epsilon_0 n_0} =3 \sqrt{2} \theta^{5/2} \sum_{s=0}^\infty(-1)^s\bar{z}^{s + 1} \int_{0}^\infty dx x^{3/2} (1 + \theta x) \sqrt{ 1 + \frac{\theta x}{2}} e^{-x(1 + s)} .\end{aligned} \hspace{\stretch{1}}(1.0.33c)

Low density result

An exact integration of the various integrals above is possible in terms of special functions. However, that attempt (included below) introduced an erroneous extra factor of $\theta$. Given that this end result was obtained by tossing all but the lowest order terms in $\theta$ and $\bar{z}$, let’s try that right from the get go.

For the pressure we have an integrand containing a factor

\begin{aligned}\left( (\theta x + 1)^2 -1 \right)^{3/2}&= \left( \theta x + 1 - 1 \right)^{3/2}\left( \theta x + 1 + 1 \right)^{3/2} \\ &= \theta^{3/2} x^{3/2} 2^{3/2} \left( 1 + \frac{\theta x}{2} \right)^{3/2} \\ &= 2 \sqrt{2} \theta^{3/2} x^{3/2} \left( 1 + \frac{\theta x}{2} \right)^{3/2}\approx2 \sqrt{2} \theta^{3/2} x^{3/2} \end{aligned} \hspace{\stretch{1}}(1.0.33c)

Our pressure, to lowest order in $\theta$ and $\bar{z}$ is then

\begin{aligned}\frac{P}{\epsilon_0 n_0} = 2 \sqrt{2} \theta^{5/2} \bar{z} \int_0^\infty x^{3/2} e^{-x} dx= 2 \sqrt{2} \theta^{5/2} \bar{z} \Gamma(5/2).\end{aligned} \hspace{\stretch{1}}(1.0.33c)

Our energy density to lowest order in $\theta$ and $\bar{z}$ from eq. 1.0.33c is

\begin{aligned}\frac{U}{V \epsilon_0 n_0} &= 3 \sqrt{2} \theta^{5/2} \bar{z} \int_{0}^\infty dx x^{3/2} e^{-x} \\ &= 3 \sqrt{2} \theta^{5/2} \bar{z} \Gamma(5/2).\end{aligned} \hspace{\stretch{1}}(1.0.33c)

Comparing these, we have

\begin{aligned}\frac{1}{{\epsilon_0 n_0\sqrt{2} \theta^{5/2} \bar{z} \Gamma(5/2)}} &= 3 \frac{V}{U} \\ &= \frac{2}{P},\end{aligned} \hspace{\stretch{1}}(1.0.37)

or in this low density limit

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

High density limit

For the high density limit write $\bar{z} = e^y$, so that the distribution takes the form

\begin{aligned}f(\bar{z}) &= \frac{1}{ \bar{z}^{-1} e^{x} + 1} \\ &= \frac{1}{ e^{x - y} + 1}.\end{aligned} \hspace{\stretch{1}}(1.0.39)

This can be approximated by a step function, so that

\begin{aligned}\frac{P}{n_0 \epsilon_0}\approx \int_{0}^y \theta dx\left( (\theta x + 1)^2 - 1 \right)^{3/2} \end{aligned} \hspace{\stretch{1}}(1.0.40a)

\begin{aligned}\frac{U}{V \epsilon_0 n_0} \approx3 \int_{0}^\infty \theta dx \theta x (\theta x + 1)\sqrt{(\theta x + 1)^2 - 1} \end{aligned} \hspace{\stretch{1}}(1.0.40b)

With a change of variables $u = \theta x + 1$, we have

\begin{aligned}\begin{aligned}\frac{P}{n_0 \epsilon_0} &\approx \int_{1}^{\theta y + 1x} du\left( u^2 - 1 \right)^{3/2} \\ &=\frac{1}{8} \left((2 \theta y (\theta y+2)-3) \sqrt{\theta y (\theta y+2)} (\theta y+1)+3 \ln \left(\theta y+\sqrt{\theta y (\theta y+2)}+1\right)\right) \\ &\approx\frac{1}{4} \left( \theta \ln \bar{z} \right)^4\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.41a)

\begin{aligned}\begin{aligned}\frac{U}{V \epsilon_0 n_0} &\approx3 \int_{1}^{\theta y + 1x} (u^2 - u)\sqrt{u^2 - 1} \\ &=\frac{3}{24} \left(\sqrt{\theta y (\theta y+2)} (\theta y (2 \theta y (3 \theta y+5)-1)+3)-3 \left(\ln \left(\theta y+\sqrt{\theta y (\theta y+2)}+1\right)\right)\right) \\ &\approx\frac{3}{4} \left( \theta \ln \bar{z} \right)^4\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.41b)

Comparing both we have

\begin{aligned}\frac{4}{\epsilon_0 n_0 \left( \theta \ln \bar{z} \right) } = \frac{1}{{P}} = \frac{3 V}{U},\end{aligned} \hspace{\stretch{1}}(1.0.42)

or

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

\begin{aligned}{\left.{{\epsilon_{\mathrm{F}}}}\right\vert}_{{n = 1/(0.01)^3}} = 6.12402 \times 10^{-35} \text{J} \times 6.24150934 \times 10^{18} \frac{\text{eV}}{\text{J}} = 3.82231 \times 10^{-16} \text{eV}\end{aligned} \hspace{\stretch{1}}(1.0.43)

Wow. That’s pretty low!

Pressure integral

Of these the pressure integral is yields directly to Mathematica

\begin{aligned}\begin{aligned} \int_{0}^\infty & dx\left( (\theta x + 1)^2 - 1 \right)^{3/2} e^{-x(1 + s)} \\ &=\frac{3 \theta e^{(s+1)/\theta}}{(s + 1)^2} K_2\left( \frac{s+1}{\theta } \right) \\ &=\frac{3 \sqrt{\frac{\pi }{2}} \theta ^{3/2}}{(s+1)^{5/2}}+\frac{45 \sqrt{\frac{\pi }{2}} \theta ^{5/2}}{8 (s+1)^{7/2}}+\frac{315 \sqrt{\frac{\pi }{2}} \theta ^{7/2}}{128 (s+1)^{9/2}}-\frac{945 \sqrt{\frac{\pi }{2}} \theta ^{9/2}}{1024 (s+1)^{11/2}}+\frac{31185 \sqrt{\frac{\pi }{2}} \theta ^{11/2}}{32768 (s+1)^{13/2}} + \cdots\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.45)

where $K_2(z)$ is a modified Bessel function [5] of the second kind as plotted in fig. 1.2.

Fig 1.2: Modified Bessel function of the second kind

Plugging this into the series for the pressure, we have

\begin{aligned}\frac{P}{n_0 \epsilon_0}= 3 \left( \frac{k_{\mathrm{B}} T}{\epsilon_0} \right)^2\sum_{s=0}^\infty(-1)^s\frac{\left( \bar{z} e^{\epsilon_0/k_{\mathrm{B}} T} \right)^{s + 1}}{(s + 1)^2}K_2\left( (s+1) \epsilon_0/k_{\mathrm{B}} T \right).\end{aligned} \hspace{\stretch{1}}(1.0.46)

Plotting the summands $3 (-1)^s \frac{\theta^2}{(s + 1)^2} \left( \bar{z} e^{ 1/\theta} \right)^{s + 1} K_2\left((s+1)/\theta\right)$ for $\bar{z} = 1$ in fig. 1.4 shows that this mix of exponential Bessel and quadratic terms decreases with $s$.

Plotting this sum in fig. 1.3 numerically to 10 terms, shows that we have a function that appears roughly polynomial in $\bar{z}$ and $\theta$.

Fig 1.3: Pressure to ten terms in z and theta

Fig 1.4: Pressure summands

For small $\bar{z}$ it can be seen graphically that there is very little contribution from anything but the $s = 0$ term of this sum. An expansion in series for a few terms in $\bar{z}$ and $\theta$ gives us

\begin{aligned}\begin{aligned}\frac{P}{\epsilon_0 n_0}&=\sqrt{\pi} \theta^{5/2} \left(\frac{3 \bar{z}}{\sqrt{2}}-\frac{3 \bar{z}^2}{8}+\frac{\bar{z}^3}{3 \sqrt{6}}-\frac{3 \bar{z}^4}{32 \sqrt{2}}+\frac{3 \bar{z}^5}{25 \sqrt{10}}\right) \\ &+\sqrt{\pi} \theta^{7/2} \left(\frac{45 \bar{z}}{8 \sqrt{2}}\right) -\frac{45 \bar{z}^2}{128}+\frac{5 \bar{z}^3}{24 \sqrt{6}}-\frac{45 \bar{z}^4}{1024 \sqrt{2}}+\frac{9 \bar{z}^5}{200 \sqrt{10}}\\ &+\sqrt{\pi} \theta^{9/2} \left(\frac{315 \bar{z}}{128 \sqrt{2}}-\frac{315 \bar{z}^2}{4096}+\frac{35 \bar{z}^3}{1152 \sqrt{6}}-\frac{315 \bar{z}^4}{65536 \sqrt{2}}+\frac{63 \bar{z}^5}{16000 \sqrt{10}}\right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.47)

This allows a $k_{\mathrm{B}} T \ll m c^2$ and $\bar{z} \ll 1$ approximation of the pressure

\begin{aligned}\frac{P}{\epsilon_0 n_0} = \frac{3}{2} \sqrt{2 \pi} \bar{z} \theta^{5/2}.\end{aligned} \hspace{\stretch{1}}(1.0.48)

Number density integral

For the number density, it appears that we can evaluate the integral using integration from parts applied to eq. 1.0.30.30

\begin{aligned}\frac{n}{n_0}= \theta\int_{0}^\infty dx\frac{3 (\theta x + 1)\sqrt{(\theta x + 1)^2 - 1} }{ \bar{z}^{-1} e^{x} + 1}=\theta\int_{0}^\infty dx\left( \frac{d}{dx} \left( (\theta x + 1)^2 - 1 \right) ^{3/2} \right)\frac{1}{ \bar{z}^{-1} e^{x} + 1}={\left.{{\theta\left( (\theta x + 1)^2 - 1 \right)^{3/2}\frac{1}{ \bar{z}^{-1} e^{x} + 1}}}\right\vert}_{{0}}^{{\infty}}-\theta\int_{0}^\infty dx\left( (\theta x + 1)^2 - 1 \right)^{3/2}\frac{ -\bar{z}^{-1} e^{x} }{ \left( \bar{z}^{-1} e^{x} + 1 \right)^2}=\theta\int_{0}^\infty dx\left( (\theta x + 1)^2 - 1 \right)^{3/2}\frac{ \bar{z} e^{-x} }{ \left( 1 + \bar{z} e^{-x} \right)^2}.\end{aligned} \hspace{\stretch{1}}(1.0.48)

Expanding in series, gives us

\begin{aligned}\frac{n}{n_0}=\theta\sum_{s = 0}^\infty\binom{-2}{s}\bar{z}^{s + 1} \int_{0}^\infty dx\left( (\theta x + 1)^2 - 1 \right)^{3/2} e^{-x(s + 1)}=3 \theta^2\sum_{s = 0}^\infty\binom{-2}{s}\frac{\left( \bar{z} e^{1/\theta} \right)^{s + 1}}{(s + 1)^2}K_2\left( \frac{s+1}{\theta } \right).\end{aligned} \hspace{\stretch{1}}(1.0.48)

Here the binomial coefficient has the meaning given in the definitions of \statmechchapcite{nonIntegralBinomialSeries}, where for negative integral values of $b$ we have

\begin{aligned}\binom{b}{s}\equiv(-1)^s \frac{-b}{-b + s} \binom{-b+s}{-b}.\end{aligned} \hspace{\stretch{1}}(1.0.51)

Expanding in series to a couple of orders in $\theta$ and $\bar{z}$ we have

\begin{aligned}\frac{n}{n_0} = \frac{\sqrt{2 \pi}}{36} \theta^{1/2} \left(\left(2 \sqrt{3} \bar{z} - 9/\sqrt{2} \right) \bar{z} +18 \right) \bar{z}+\frac{5 \sqrt{ 2 \pi}}{576} \theta^{3/2} \left(\left(4 \sqrt{3} \bar{z} - 27/\sqrt{2}\right) \bar{z} +108 \right) \bar{z}+ \cdots\end{aligned} \hspace{\stretch{1}}(1.0.52)

To first order in $\theta$ and $\bar{z}$ this is

\begin{aligned}\frac{n}{n_0} = \frac{1}{{2}} \sqrt{ 2 \pi } \bar{z} \theta^{1/2},\end{aligned} \hspace{\stretch{1}}(1.0.53)

which allows a relation to pressure

\begin{aligned}P V = 3 N (k_{\mathrm{B}} T)^2 /\epsilon_0.\end{aligned} \hspace{\stretch{1}}(1.0.54)

It’s kind of odd seeming that this is quadratic in temperature. Is there an error?

Energy integral

Starting from eq. 1.0.30c and integrating by parts we have

\begin{aligned}\frac{\left\langle{{\epsilon - \epsilon_0}}\right\rangle}{V \epsilon_0 n_0} &= 3 \theta^2 \int_{0}^\infty dx \frac { x (\theta x + 1)\sqrt{(\theta x + 1)^2 - 1} } { \bar{z}^{-1} e^{x} + 1 } \\ &= -\theta^2 \int_{0}^\infty dx\left( (\theta x + 1)^2 - 1 \right)^{3/2}\frac{d}{dx} \left( \frac{x} { \bar{z}^{-1} e^{x} + 1 } \right) \\ &= -\theta^2 \int_{0}^\infty dx\left( (\theta x + 1)^2 - 1 \right)^{3/2}\left( \frac{1} { \bar{z}^{-1} e^{x} + 1 } - \frac{x \bar{z}^{-1} e^{x} } { \left( \bar{z}^{-1} e^{x} + 1 \right)^2 } \right) \\ &= \theta^2 \int_{0}^\infty dx \left( (\theta x + 1)^2 - 1 \right)^{3/2} \frac{ (x - 1)\bar{z}^{-1} e^{x} - 1} { \left( \bar{z}^{-1} e^{x} + 1 \right)^2 } \\ &= \theta^2 \int_{0}^\infty dx \left( (\theta x + 1)^2 - 1 \right)^{3/2} \frac{ (x - 1)\bar{z} e^{-x} - \bar{z}^2 e^{-2 x}} { \left( 1 + \bar{z} e^{-x} \right)^2 } \\ &= \theta^2\sum_{s=0}^\infty \binom{-2}{s} \int_{0}^\infty dx \left( (\theta x + 1)^2 - 1 \right)^{3/2} \left( (x - 1)\bar{z} e^{-x} - \bar{z}^2 e^{-2 x} \right) (\bar{z} e^{-x})^s \\ &= \theta^2\sum_{s=0}^\infty \binom{-2}{s} \bar{z}^{s + 1} \int_{0}^\infty dx \left( (\theta x + 1)^2 - 1 \right)^{3/2} \left( (x - 1) e^{-x(s + 1)} - \bar{z} e^{-x(s + 2)} \right).\end{aligned} \hspace{\stretch{1}}(1.0.54)

The integral with the factor of $x$ doesn’t have a nice closed form as before (if you consider the $K_2$ a nice closed form), but instead evaluates to a confluent hypergeometric function [6]. That integral is

\begin{aligned}\int_0^{\infty } x \left((\theta x+1)^2-1\right)^{3/2} e^{-x (1+s)} dx = \frac{15 \sqrt{\pi } \theta^3 U\left(-\frac{3}{2},-4,\frac{2 (s+1)}{\theta }\right)}{8 (s+1)^5},\end{aligned} \hspace{\stretch{1}}(1.0.54)

and looks like fig. 1.5. Series expansion shows that this hypergeometricU function has a $\theta^{3/2}$ singularity at the origin

Fig 1.5: Plot of HypergeometricU, and with theta^5 scaling

\begin{aligned}U\left(-\frac{3}{2},-4,\frac{2 (s+1)}{\theta }\right)=\frac{2 \sqrt{2} \sqrt{s+1} s+2 \sqrt{2} \sqrt{s+1}}{\theta^{3/2}}+\frac{21 \sqrt{s+1}}{2 \sqrt{2} \sqrt{\theta }}+ \cdots\end{aligned} \hspace{\stretch{1}}(1.57)

so our multiplication by $\theta^5$ brings us to zero as seen in the plot. Evaluating the complete integral yields the unholy mess

\begin{aligned}\frac{\left\langle{{\epsilon - \epsilon_0}}\right\rangle}{V \epsilon_0 n_0} &= \sum_{s=0}^\infty \theta^2 (-1)^s (s+1) \bar{z}^{s+1} \Bigl( \\ &\frac{105 \sqrt{\pi } \theta^3 U\left(-\frac{1}{2},-4,\frac{2 (s+1)}{\theta }\right)}{16 (s+1)^5} \\ &- \frac{3 \sqrt{\pi } \theta^2 U\left(-\frac{1}{2},-2,\frac{2 (s+1)}{\theta }\right)}{2 (s+1)^3} \\ &- \frac{3 \sqrt{\pi } \theta^2 \bar{z} U\left(-\frac{1}{2},-2,\frac{2 (s+2)}{\theta }\right)}{2 (s+2)^3} \\ &+\frac{(\theta -2) (-3 \theta +2 s+2) e^{\frac{s+1}{\theta }} K_2\left(\frac{s+1}{\theta }\right)}{\theta (s+1)^2} \\ &-\frac{2 (\theta -2) e^{\frac{s+1}{\theta }} K_1\left(\frac{s+1}{\theta }\right)}{\theta (s+1)} \\ &+\frac{\bar{z} (-3 \theta +2 s+4) e^{\frac{s+2}{\theta }} K_2\left(\frac{s+2}{\theta }\right)}{(s+2)^2} \\ &-\frac{2 \bar{z} e^{\frac{s+2}{\theta }} K_1\left(\frac{s+2}{\theta }\right)}{s+2} \Bigr),\end{aligned} \hspace{\stretch{1}}(1.58)

to first order in $\bar{z}$ and $\theta$ this is

\begin{aligned}\frac{\left\langle{{\epsilon - \epsilon_0}}\right\rangle}{V \epsilon_0 n_0} =\frac{9}{4} \sqrt{2 \pi} \bar{z} \theta^{7/2}.\end{aligned} \hspace{\stretch{1}}(1.59)

Comparing pressure and energy we have for low densities (where $\bar{z} \approx 0$)

\begin{aligned}\frac{1}{{\epsilon_0 n_0 \sqrt{2 \pi} \bar{z} \theta^{5/2}}} = \frac{3}{2} \frac{1}{{P}} = \frac{9}{4} \theta \frac{V}{U},\end{aligned} \hspace{\stretch{1}}(1.0.60)

or

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

It appears that I’ve picked up an extra factor of $\theta$ somewhere, but at least I’ve got the $2/3$ low density expression. Given that I’ve Taylor expanded everything anyways around $\bar{z}$ and $\theta$ this could likely have been done right from the get go, instead of dragging along the messy geometric integrals. Reworking this part of this problem like that was done above.

# References

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

[2] Peeter Joot. Basic statistical mechanics., chapter {Non integral binomial coefficient}. \natexlab{a}. URL http://sites.google.com/site/peeterjoot2/math2013/phy452.pdf.

[3] Peeter Joot. Basic statistical mechanics., chapter {Relativisitic density of states}. \natexlab{b}. URL http://sites.google.com/site/peeterjoot2/math2013/phy452.pdf.

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

[5] Wolfram. BesselK, \natexlab{a}. URL http://reference.wolfram.com/mathematica/ref/BesselK.html. [Online; accessed 11-April-2013].

[6] Wolfram. HyperGeometricU, \natexlab{b}. URL http://reference.wolfram.com/mathematica/ref/HypergeometricU.html. [Online; accessed 17-April-2013].

## 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

This is an ungraded set of answers to the problems posed.

## 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.

## Pathria chapter 4 diatomic molecule problem

Posted by peeterjoot on March 18, 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.)]

## Question: Diatomic molecule ([1] pr 4.7)

Consider a classical system of non-interacting, diatomic molecules enclosed in a box of volume $V$ at temperature $T$. The Hamiltonian of a single molecule is given by

\begin{aligned}H(\mathbf{r}_1, \mathbf{r}_2, \mathbf{p}_1, \mathbf{p}_2) = \frac{1}{{2m}} \left( \mathbf{p}_1^2 + \mathbf{p}_2^2 \right)+\frac{1}{{2}} K \left\lvert {\mathbf{r}_1 - \mathbf{r}_2} \right\rvert^2.\end{aligned} \hspace{\stretch{1}}(1.0.1)

Study the thermodynamics of this system, including the dependence of the quantity $\left\langle{{r_{12}^2}}\right\rangle$ on $T$.

Partition function
First consider the partition function for a single diatomic pair

\begin{aligned}Z_1 &= \frac{1}{{h^6}} \int d^6 \mathbf{p} d^6 \mathbf{r} e^{-\beta \frac{ \mathbf{p}_1^2 + \mathbf{p}_2^2 }{2m}} e^{-\beta K\frac{ \left\lvert {\mathbf{r}_1 - \mathbf{r}_2} \right\rvert^2 }{2}} \\ &= \frac{1}{{h^6}} \left( \frac{2 \pi m}{\beta} \right)^{6/2}\int d^3 \mathbf{r}_1 d^3 \mathbf{r}_2 e^{-\beta K\frac{ \left\lvert {\mathbf{r}_1 - \mathbf{r}_2} \right\rvert^2 }{2}}\end{aligned} \hspace{\stretch{1}}(1.0.2)

Now we can make a change of variables to simplify the exponential. Let’s write

\begin{aligned}\mathbf{u} = \mathbf{r}_1 - \mathbf{r}_2\end{aligned} \hspace{\stretch{1}}(1.0.3a)

\begin{aligned}\mathbf{v} = \mathbf{r}_2,\end{aligned} \hspace{\stretch{1}}(1.0.3b)

or

\begin{aligned}\mathbf{r}_2 = \mathbf{v}\end{aligned} \hspace{\stretch{1}}(1.0.4a)

\begin{aligned}\mathbf{r}_1=\mathbf{u} + \mathbf{v}.\end{aligned} \hspace{\stretch{1}}(1.0.4b)

Our volume element is

\begin{aligned}d^3 \mathbf{r}_1 d^3 \mathbf{r}_2 = d^3 \mathbf{u} d^3 \mathbf{v} \frac{\partial(\mathbf{r}_1, \mathbf{r}_2)}{\partial(\mathbf{u}, \mathbf{v})}.\end{aligned} \hspace{\stretch{1}}(1.0.5)

It wasn’t obvious to me that this change of variables preserves the volume element, but a quick Jacobian calculation shows this to be the case

\begin{aligned}\frac{\partial(\mathbf{r}_1, \mathbf{r}_2)}{\partial(\mathbf{u}, \mathbf{v})} &= \begin{vmatrix}\partial r_{11}/\partial u_1 & \partial r_{11}/\partial u_2 &\partial r_{11}/\partial u_3 &\partial r_{11}/\partial v_1 &\partial r_{11}/\partial v_2 &\partial r_{11}/\partial v_3 \\ \partial r_{12}/\partial u_1 & \partial r_{12}/\partial u_2 &\partial r_{12}/\partial u_3 &\partial r_{12}/\partial v_1 &\partial r_{12}/\partial v_2 &\partial r_{12}/\partial v_3 \\ \partial r_{13}/\partial u_1 & \partial r_{13}/\partial u_2 &\partial r_{13}/\partial u_3 &\partial r_{13}/\partial v_1 &\partial r_{13}/\partial v_2 &\partial r_{13}/\partial v_3 \\ \partial r_{21}/\partial u_1 & \partial r_{21}/\partial u_2 &\partial r_{21}/\partial u_3 &\partial r_{21}/\partial v_1 &\partial r_{21}/\partial v_2 &\partial r_{21}/\partial v_3 \\ \partial r_{22}/\partial u_1 & \partial r_{22}/\partial u_2 &\partial r_{22}/\partial u_3 &\partial r_{22}/\partial v_1 &\partial r_{22}/\partial v_2 &\partial r_{22}/\partial v_3 \\ \partial r_{23}/\partial u_1 & \partial r_{23}/\partial u_2 &\partial r_{23}/\partial u_3 &\partial r_{23}/\partial v_1 &\partial r_{23}/\partial v_2 &\partial r_{23}/\partial v_3 \end{vmatrix} \\ &= \begin{vmatrix}\partial r_{11}/\partial u_1 & \partial r_{11}/\partial u_2 &\partial r_{11}/\partial u_3 &\partial r_{11}/\partial v_1 &\partial r_{11}/\partial v_2 &\partial r_{11}/\partial v_3 \\ \partial r_{12}/\partial u_1 & \partial r_{12}/\partial u_2 &\partial r_{12}/\partial u_3 &\partial r_{12}/\partial v_1 &\partial r_{12}/\partial v_2 &\partial r_{12}/\partial v_3 \\ \partial r_{13}/\partial u_1 & \partial r_{13}/\partial u_2 &\partial r_{13}/\partial u_3 &\partial r_{13}/\partial v_1 &\partial r_{13}/\partial v_2 &\partial r_{13}/\partial v_3 \\ 0 & 0 & 0 &\partial r_{21}/\partial v_1 &\partial r_{21}/\partial v_2 &\partial r_{21}/\partial v_3 \\ 0 & 0 & 0 &\partial r_{22}/\partial v_1 &\partial r_{22}/\partial v_2 &\partial r_{22}/\partial v_3 \\ 0 & 0 & 0 &\partial r_{23}/\partial v_1 &\partial r_{23}/\partial v_2 &\partial r_{23}/\partial v_3 \end{vmatrix} \\ &= 1.\end{aligned} \hspace{\stretch{1}}(1.0.6)

Our remaining integral can now be evaluated

\begin{aligned}\int d^3 \mathbf{r}_1 d^3 \mathbf{r}_2 e^{-\beta K\frac{ \left\lvert {\mathbf{r}_1 - \mathbf{r}_2} \right\rvert^2 }{2}} &= \int d^3 \mathbf{u} d^3 \mathbf{v} e^{-\beta K \left\lvert {\mathbf{u}} \right\rvert^2 /2 } \\ &= V \int d^3 \mathbf{u} e^{-\beta K \left\lvert {\mathbf{u}} \right\rvert^2 /2 } \\ &= V \int d^3 \mathbf{u} e^{-\beta K \left\lvert {\mathbf{u}} \right\rvert^2 /2 } \\ &= V \left( \frac{ 2 \pi }{ K \beta } \right)^{3/2}.\end{aligned} \hspace{\stretch{1}}(1.0.7)

Our partition function is now completely evaluated

\begin{aligned}Z_1 = V\frac{1}{{h^6}} \left( \frac{2 \pi m}{\beta} \right)^{3}\left( \frac{ 2 \pi }{ K \beta } \right)^{3/2}.\end{aligned} \hspace{\stretch{1}}(1.0.8)

As a function of $V$ and $T$ as in the text, we write

\begin{aligned}Z_1 = V f(T)\end{aligned} \hspace{\stretch{1}}(1.0.9a)

\begin{aligned}f(T) = \left( \frac{m }{h^2 } \sqrt{\frac{(2\pi)^3}{K}} \right)^3\left( k_{\mathrm{B}} T \right)^{9/2}.\end{aligned} \hspace{\stretch{1}}(1.0.9b)

Gibbs sum

Our Gibbs sum, summing over the number of molecules (not atoms), is

\begin{aligned}Z_{\mathrm{G}} &= \sum_{N_r = 0}^\infty \frac{z^{N_r}}{N_r!} Z_1^{N_r} \\ &= e^{ z V f(T) },\end{aligned} \hspace{\stretch{1}}(1.0.10)

or

\begin{aligned}q &= \ln Z_{\mathrm{G}} \\ &= z V f(T) \\ &= P V \beta.\end{aligned} \hspace{\stretch{1}}(1.0.11)

The fact that we can sum this as an exponential series so nicely looks like it’s one of the main advantages to this grand partition function (Gibbs sum). We can avoid any of the large $N!$ approximations that we have to use when the number of particles is explicitly fixed.

Pressure

The pressure follows

\begin{aligned}P &= z f(T) k_{\mathrm{B}} T \\ &= e^{\mu/k_{\mathrm{B}} T}\left( \frac{m }{h^2 } \sqrt{\frac{(2\pi)^3}{K}} \right)^3\left( k_{\mathrm{B}} T \right)^{11/2}.\end{aligned} \hspace{\stretch{1}}(1.0.12)

Average energy

\begin{aligned}\left\langle{{H}}\right\rangle &= -\frac{\partial {q}}{\partial {\beta}} \\ &= - z V \frac{9}{2} \frac{f(T)}{T} \frac{\partial {T}}{\partial {\beta}} \\ &= z V \frac{9}{2} \frac{f(T)}{T^3} \frac{1}{{k_{\mathrm{B}}}},\end{aligned} \hspace{\stretch{1}}(1.0.13)

or

\begin{aligned}\left\langle{{H}}\right\rangle = e^{\mu/k_{\mathrm{B}} T} V \frac{9}{2} k_{\mathrm{B}}^2 \left( \frac{m }{h^2 } \sqrt{\frac{(2\pi)^3}{K}} \right)^3\left( k_{\mathrm{B}} T \right)^{3/2}.\end{aligned} \hspace{\stretch{1}}(1.0.14)

Average occupancy

\begin{aligned}\left\langle{{N}}\right\rangle &= z \frac{\partial {}}{\partial {z}} \ln Z_{\mathrm{G}} \\ &= z \frac{\partial {}}{\partial {z}} \left( z V f(T) \right) \\ &= z V f(T)\end{aligned} \hspace{\stretch{1}}(1.0.15)

but this is just $q$, or

\begin{aligned}\left\langle{{N}}\right\rangle &= e^{\mu/k_{\mathrm{B}} T} V\left( \frac{m }{h^2 } \sqrt{\frac{(2\pi)^3}{K}} \right)^3\left( k_{\mathrm{B}} T \right)^{9/2}.\end{aligned} \hspace{\stretch{1}}(1.0.16)

Free energy

\begin{aligned}F &= - k_{\mathrm{B}} T \ln \frac{ Z_{\mathrm{G}} }{z^N} \\ &= - k_{\mathrm{B}} T \left( q - N \ln z \right) \\ &= N k_{\mathrm{B}} T \beta \mu - k_{\mathrm{B}} T q \\ &= z V f(T) \mu - k_{\mathrm{B}} T z V f(T) \\ &= z V f(T) \left( \mu - k_{\mathrm{B}} T \right)\end{aligned} \hspace{\stretch{1}}(1.0.17)

\begin{aligned}F = e^{\mu/k_{\mathrm{B}} T} V \left( \mu - k_{\mathrm{B}} T \right)\left( \frac{m }{h^2 } \sqrt{\frac{(2\pi)^3}{K}} \right)^3\left( k_{\mathrm{B}} T \right)^{9/2}.\end{aligned} \hspace{\stretch{1}}(1.0.18)

Entropy

\begin{aligned}S &= \frac{U - F}{T} \\ &= \frac{V}{T} e^{\mu/k_{\mathrm{B}} T} \left( k_{\mathrm{B}} T \right)^{3/2}\left( \frac{m }{h^2 } \sqrt{\frac{(2\pi)^3}{K}} \right)^3\left( \frac{9}{2} k_{\mathrm{B}}^2 - \left( \mu - k_{\mathrm{B}} T \right) \left( k_{\mathrm{B}} T \right)^3 \right).\end{aligned} \hspace{\stretch{1}}(1.0.19)

Expectation of atomic separation

The momentum portions of the average will just cancel out, leaving just

\begin{aligned}\left\langle{r_{12}^2}\right\rangle &= \frac{\int d^3 \mathbf{r}_1 d^3 \mathbf{r}_2 \left( \mathbf{r}_1 - \mathbf{r}_2 \right)^2 e^{-\beta K \left( \mathbf{r}_1 - \mathbf{r}_2 \right)^2 /2 }}{\int d^3 \mathbf{r}_1 d^3 \mathbf{r}_2 e^{-\beta K \left( \mathbf{r}_1 - \mathbf{r}_2 \right)^2 /2 }} \\ &= \frac{ \int d^3 \mathbf{u} \mathbf{u}^2 e^{-\beta K \mathbf{u}^2 /2 }}{\int d^3 \mathbf{u} e^{-\beta K \mathbf{u}^2 /2 }} \\ &= \frac{\int da db dc \left( a^2 + b^2 + c^2 \right) e^{-\beta K \left( a^2 + b^2 + c^2 \right) /2}}{\int e^{-\beta K \left( a^2 + b^2 + c^2 \right)/2}} \\ &= 3 \frac{\int da a^2 e^{-\beta K a^2/2}\int db dc e^{-\beta K \left( b^2 + c^2 \right) /2}}{\int e^{-\beta K \left( a^2 + b^2 + c^2 \right)/2 }} \\ &= 3 \frac{\int da a^2 e^{-\beta K a^2/2}}{\int e^{-\beta K a^2/2}}\end{aligned} \hspace{\stretch{1}}(1.0.20)

Expanding the numerator by parts we have

\begin{aligned}\int da a^2 e^{-\beta K a^2/2} \\ &= \int a d\frac{ e^{-\beta K a^2/2}}{- 2 \beta K/2} \\ &= \frac{1}{\beta K}\int e^{-\beta K a^2/2}.\end{aligned} \hspace{\stretch{1}}(1.0.21)

This gives us

\begin{aligned}\boxed{\left\langle r_{12}^2 \right\rangle = \frac{3}{\beta K} = \frac{3 k_{\mathrm{B}} T}{K}.}\end{aligned} \hspace{\stretch{1}}(1.0.22)

# References

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

## probability forms of entropy

Posted by peeterjoot on March 16, 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.)]

## Question: Entropy as probability

[1] points out that entropy can be written as

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

where

\begin{aligned}P_i = \frac{e^{-\beta E_i}}{Z}\end{aligned} \hspace{\stretch{1}}(1.0.2a)

\begin{aligned}Z = \sum_i e^{-\beta E_i}.\end{aligned} \hspace{\stretch{1}}(1.0.2b)

Show that this follows from the free energy $F = U - T S = -k_{\mathrm{B}} \ln Z$.

In terms of the free and average energies, we have

\begin{aligned}\frac{S}{k_{\mathrm{B}}} &= \frac{U - F}{k_{\mathrm{B}} T} \\ &= \beta \left( -\frac{\partial {\ln Z}}{\partial {\beta}} \right) - \beta \left( -k_{\mathrm{B}} T \ln Z \right) \\ &= \frac{\sum_i \beta E_i e^{-\beta E_i}}{Z} +\ln Z \\ &= -\sum_i P_i \ln e^{-\beta E_i} + \sum_i P_i \ln Z \\ &= -\sum_i P_i \ln \frac{e^{-\beta E_i}}{Z} P_i \\ &= -\sum_i P_i \ln P_i.\end{aligned} \hspace{\stretch{1}}(1.0.3)

## Question: Entropy in terms of grand partition probabilites ( [2] pr 4.1)

Generalize \cref{pr:entropyProbabilityForm:1} to the grand canonical scheme, where we have

\begin{aligned}P_{r, s} = \frac{e^{-\alpha N_r - \beta E_s}}{Z_{\mathrm{G}}}\end{aligned} \hspace{\stretch{1}}(1.0.4a)

\begin{aligned}Z_{\mathrm{G}} = \sum_{r,s} e^{-\alpha N_r - \beta E_s}\end{aligned} \hspace{\stretch{1}}(1.0.4b)

\begin{aligned}z = e^{-\alpha} = e^{\mu \beta}\end{aligned} \hspace{\stretch{1}}(1.0.4c)

\begin{aligned}q = \ln Z_{\mathrm{G}},\end{aligned} \hspace{\stretch{1}}(1.0.4d)

and show

\begin{aligned}S = - k_{\mathrm{B}} \sum_{r,s} P_{r,s} \ln P_{r,s}.\end{aligned} \hspace{\stretch{1}}(1.0.5)

With

\begin{aligned}\beta P V = q,\end{aligned} \hspace{\stretch{1}}(1.0.6)

the free energy takes the form

\begin{aligned}F = N \mu - P V = N \mu - q/\beta,\end{aligned} \hspace{\stretch{1}}(1.0.7)

so that the entropy (scaled by $k_{\mathrm{B}}$) leads us to the desired result

\begin{aligned}\frac{S}{k_{\mathrm{B}}} &= \beta U - N \mu \beta + q/(\beta k_{\mathrm{B}} T) \\ &= -\beta \frac{\partial {q}}{\partial {\beta}} - z \mu \beta \frac{\partial {q}}{\partial {z}} + q \\ &= \frac{1}{{Z_{\mathrm{G}}}}\sum_{r, s}\left( -\beta (-E_s) - \mu \beta N_r \right) e^{-\alpha N_r - \beta E_s}+ \ln Z_{\mathrm{G}} \\ &= \sum_{r, s} \ln e^{ \alpha N_r + \beta E_s } P_{r,s} + \left( \sum_{r, s} P_{r, s} \right)\ln Z_{\mathrm{G}} \\ &= -\sum_{r, s} \ln \frac{e^{ -\alpha N_r - \beta E_s }}{Z_{\mathrm{G}}} P_{r,s} \\ &= -\sum_{r, s} P_{r, s} \ln P_{r, s}\end{aligned} \hspace{\stretch{1}}(1.0.8)

# References

[1] E.A. Jackson. Equilibrium statistical mechanics. Dover Pubns, 2000.

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

## Quantum anharmonic oscillator

Posted by peeterjoot on March 13, 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.)]

## Question: Quantum anharmonic oscillator ([1] pr 3.30)

The energy levels of a quantum-mechanical, one-dimensional, anharmonic oscillator may be approximated as

\begin{aligned}\epsilon_n = \left( n + \frac{1}{{2}} \right) \hbar \omega - x \left( n + \frac{1}{{2}} \right)^2 \hbar \omega\qquad n = 0, 1, 2, \cdots\end{aligned} \hspace{\stretch{1}}(1.0.1)

The parameter $x$, usually $\ll 1$, represents the degree of anharmonicity. Show that, to the first order in $x$ and the fourth order in $u \equiv \hbar \omega/k_{\mathrm{B}} T$, the specific heat of a system of $N$ such oscillators is given by

\begin{aligned}C = N k_{\mathrm{B}} \left( \left( 1 - \frac{1}{{12}} u^2 + \frac{1}{{240}} u^4 \right) + 4 x \left( \frac{1}{{u}} + \frac{1}{{80}} u^3 \right) \right)\end{aligned} \hspace{\stretch{1}}(1.0.2)

We can expand the partition function in a first order Taylor series about $x = 0$, then evaluate the sums

\begin{aligned}Z_1 &= \sum_{n = 0}^\infty \exp\left( -\beta \left( n + \frac{1}{{2}} \right) \hbar \omega + \beta x \left( n + \frac{1}{{2}} \right)^2 \hbar \omega \right) \\ &= \sum_{n = 0}^\infty \exp\left( - \left( n + \frac{1}{{2}} \right) u + x \left( n + \frac{1}{{2}} \right)^2 u \right)\approx\sum_{n = 0}^\infty e^{ - \left( n + \frac{1}{{2}} \right) u }\left( 1 + x u \left( n + \frac{1}{{2}} \right)^2 \right).\end{aligned} \hspace{\stretch{1}}(1.0.3)

The quadratic sum can be evaluated indirectly as it can be expressed as a derivative

\begin{aligned}Z_1 &= \left( 1 + x u \frac{d^2}{du^2} \right)\sum_{n = 0}^\infty e^{ - \left( n + \frac{1}{{2}} \right) u } \\ &= \left( 1 + x u \frac{d^2}{du^2} \right)e^{-u/2}\sum_{n = 0}^\infty e^{ - n u } \\ &= \left( 1 + x u \frac{d^2}{du^2} \right)e^{-u/2}\frac{ 1 }{1 - e^{-u} } \\ &= \left( 1 + x u \frac{d^2}{du^2} \right)\frac{ 1 }{e^{u/2} - e^{-u/2} } \\ &= \left( 1 + x u \frac{d^2}{du^2} \right)\frac{ 1 }{2 \sinh(u/2)}\end{aligned} \hspace{\stretch{1}}(1.0.4)

Finally, evaluation of the derivatives gives us

\begin{aligned}Z_1=\frac{ 1 }{\sinh(u/2)}\left( 1 + x u \frac{2 \coth^2(u/2) - 1}{8} \right).\end{aligned} \hspace{\stretch{1}}(1.0.5)

Now we’d like to compute the specific heat in terms of derivatives of $u$. First, for the average energy

\begin{aligned}\left\langle{{H}}\right\rangle &= -N \frac{\partial {}}{\partial {\beta}} \ln Z_1 \\ &= -N \hbar \omega \frac{\partial {}}{\partial {u}} \ln Z_1.\end{aligned} \hspace{\stretch{1}}(1.0.6)

The specific heat is

\begin{aligned}C_{\mathrm{V}} &= \frac{\partial {\left\langle{{H}}\right\rangle}}{\partial {T}} \\ &= \frac{\partial {u}}{\partial {T}} \frac{\partial {\left\langle{{H}}\right\rangle}}{\partial {u}} \\ &= \frac{\hbar \omega}{k_{\mathrm{B}}} \frac{\partial {(1/T)}}{\partial {T}} \frac{\partial {\left\langle{{H}}\right\rangle}}{\partial {u}} \\ &= -\frac{\hbar \omega}{k_{\mathrm{B}} T^2} \frac{\partial {\left\langle{{H}}\right\rangle}}{\partial {u}} \\ &= -\frac{k_{\mathrm{B}} u^2}{\hbar \omega} \frac{\partial {\left\langle{{H}}\right\rangle}}{\partial {u}},\end{aligned} \hspace{\stretch{1}}(1.0.7)

or

\begin{aligned}C_{\mathrm{V}} = N k_{\mathrm{B}} u^2 \frac{\partial^2 {{}}}{\partial {{u}}^2} \ln Z_1.\end{aligned} \hspace{\stretch{1}}(1.0.8)

Actually computing that is messy algebra (See \nbref{pathria_3_30.nb}), and the result isn’t particularily interesting looking. The plot fig. 1.1 is interesting though and shows negative heat capacities near zero and a funny little jog near $C_{\mathrm{V}} = 0$.

Fig 1.1: Quantum anharmonic heat capacity

Also confirmed in the Mathematica notebook is equation eq. 1.0.2, which follows by first doing a first order series expansion in $x$, then a subsequent series expansion in $u$.

# References

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

## Midterm II reflection, take II, with approximate anharmonic oscillator solution

Posted by peeterjoot on March 11, 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.)]

## Question: Perturbation of classical harmonic oscillator (2013 midterm II p2)

Consider a single particle perturbation of a classical simple harmonic oscillator Hamiltonian

\begin{aligned}H = \frac{1}{{2}} m \omega^2 \left( {x^2 + y^2} \right) + \frac{1}{{2 m}} \left( {p_x^2 + p_y^2} \right) + a x^4 + by^6\end{aligned} \hspace{\stretch{1}}(1.0.12)

Calculate the canonical partition function, mean energy and specific heat of this system.

This problem can be attempted in two ways, the first of which was how I did it on the midterm, differentiating under the integral sign, leaving the integrals in exact form, but not evaluated explicitly in any way.

Alternately, by Taylor expanding around $c = 0$ and $d = 0$ with those as the variables in the Taylor expansion (as now done in the Pathria 3.29 problem), we can form a solution in short order. Given my low midterm mark, it seems very likely that this was what was expected.

Performing a two variable Taylor expansion of $Z$, about $(c, d) = (0, 0)$ we have

\begin{aligned}Z \approx\frac{2 \pi m}{\beta}\int dx dye^{- \beta m \omega^2 x^2/2}e^{- \beta m \omega^2 y^2/2}\left( 1 - \beta a x^4 - \beta b y^6 \right)=\frac{2 \pi m}{\beta}\frac{ 2 \pi}{\beta m \omega^2}\left( 1 - \beta a \frac{3!!}{(\beta m \omega^2)^2} - \beta b \frac{5!!}{(\beta m \omega^2)^3} \right),\end{aligned} \hspace{\stretch{1}}(1.0.22)

or

\begin{aligned}\boxed{Z \approx\frac{(2 \pi/\omega)^2}{\beta^2}\left( 1 - \frac{3 a }{\beta (m \omega^2)^2} - \frac{15 b }{\beta^2 (m \omega^2)^3} \right).}\end{aligned} \hspace{\stretch{1}}(1.0.23)

Now we can calculate the average energy

\begin{aligned}\left\langle{{H}}\right\rangle = - \frac{\partial {}}{\partial {\beta}}\ln Z= - \frac{\partial {}}{\partial {\beta}}\left( -2 \ln \beta + \ln \left( 1 - \frac{3 a }{\beta (m \omega^2)^2} - \frac{15 b }{\beta^2 (m \omega^2)^3} \right) \right)=\frac{2 \beta}-\frac{ \frac{3 a }{\beta^2 (m \omega^2)^2}+ \frac{30 b }{\beta^3 (m \omega^2)^3}}{ 1 - \frac{3 a }{\beta (m \omega^2)^2} - \frac{15 b }{\beta^2 (m \omega^2)^3}}.\end{aligned} \hspace{\stretch{1}}(1.0.24)

Dropping the $c$, $d$ terms of the denominator above, we have

\begin{aligned}\boxed{\left\langle{{H}}\right\rangle=\frac{2 \beta}- \frac{3 a }{\beta^2 (m \omega^2)^2}- \frac{30 b }{\beta^3 (m \omega^2)^3}.}\end{aligned} \hspace{\stretch{1}}(1.0.25)

The heat capacity follows immediately

\begin{aligned}\boxed{C_{\mathrm{V}} = \frac{1}{{k_{\mathrm{B}}}} \frac{\partial {\left\langle{{H}}\right\rangle}}{\partial {T}}= 2 - \frac{6 a k_{\mathrm{B}} T}{(m \omega^2)^2} - \frac{90 k_{\mathrm{B}}^2 T^2 b }{(m \omega^2)^3}.}\end{aligned} \hspace{\stretch{1}}(1.0.26)

## Heat capacity of perturbed harmonic oscillator

Posted by peeterjoot on March 11, 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.)]

This problem was suggested as prep for the second midterm, but I spent too much time on my problem set. That’s pretty unfortunate since this showed exactly the approach that was expected for the second midterm problem. Not hard, just not obvious in the heat of the moment how to do that Taylor expansion.

## Question: Anharmonic oscillator ([1] pr 3.29)

The potential energy of a one-dimensional, anharmonic oscillator may be written as

\begin{aligned}V(q) = c q^2 - g q^3 - f q^4,\end{aligned} \hspace{\stretch{1}}(1.1)

where $c$, $g$, and $f$ are positive constant; quite generally, $g$ and $f$ may be assumed to be very small in value.

Show that the leading contribution of anharmonic terms to the heat cpacity of the oscillator, assumed classical, is given by

\begin{aligned}\frac{3}{2} k_{\mathrm{B}}^2 \left( {\frac{f}{c^2} + \frac{5}{4} \frac{g^2}{c^3} } \right) T,\end{aligned} \hspace{\stretch{1}}(1.2)

To the same order, show that the mean value of the position coordinate $q$ is given by

\begin{aligned}\frac{3}{4} \frac{g k_{\mathrm{B}} T}{c^2}.\end{aligned} \hspace{\stretch{1}}(1.3)

Our partition function is

\begin{aligned}Z = \int dp dq e^{-\beta p^2/2m} e^{ -\beta \left( { c q^2 - g q^3 - f q^4} \right) }= \sqrt{\frac{2 \pi m}{\beta} }\int dq e^{ -\beta \left( { c q^2 - g q^3 - f q^4} \right) }\end{aligned} \hspace{\stretch{1}}(1.4)

How to expand this wasn’t immediately clear to me (as it wasn’t on the midterm either). We can’t Taylor expand in $q$, because there’s no single position $q$ that is of interest to expand around (we are integrating over all $q$). What we can do though is Taylor expand about the values $f$ and $g$, which are assumed to be small. Here’s the two variable Taylor expansion of this perturbed harmonic oscillator exponential. With

\begin{aligned}A(f, g) = e^{ -\beta \left( { c q^2 - g q^3 - f q^4} \right) }\end{aligned} \hspace{\stretch{1}}(1.5)

The expansion to second order is

\begin{aligned}A(f, g) &= A(0, 0) + f {\left. \frac{\partial A}{\partial f} \right\vert}_{f = 0} + g {\left. \frac{\partial A}{\partial g} \right\vert}_{g = 0} + \frac{1}{2} f^2 {\left. \frac{\partial^2 A}{\partial f^2} \right\vert}_{f = 0} + \frac{1}{2} g^2 {\left. \frac{\partial^2 A}{\partial g^2} \right\vert}_{g = 0} + f g {\left. \frac{\partial^2 A}{\partial g \partial f} \right\vert}_{f, g = 0} + \cdots \\ &= e^{ -\beta c q^2 } \left( 1 + g \beta q^3 + f \beta q^4+ \frac{1}{2} g^2 \left( \beta q^3 \right)^2 + \frac{1}{2} f^2 \left( \beta q^4 \right)^2 + f g \left( \beta q^3 \right) \left( \beta q^4 \right) + \cdots \right) \\ &= e^{ -\beta c q^2 } \left( 1 + g \beta q^3 + f \beta q^4+ \frac{1}{2} g^2 \beta^2 q^6+ f g \beta^2 q^7+ \frac{1}{2} f^2 \beta^2 q^8 + \cdots \right) \end{aligned} \hspace{\stretch{1}}(1.6)

This can now be integrated by parts, where any odd powers are killed. For even powers we have

\begin{aligned}\int q^{2 N} e^{-a q^2} dq &= \int q^{2 N - 1} d \frac{ e^{-a q^2} }{-2 a} \\ &= \frac{2 N - 1}{2a} \int q^{2 (N - 1)} e^{-a q^2} dq \\ &= \frac{(2 N - 1)!!}{(2a)^N} \sqrt{\frac{\pi}{a}}.\end{aligned} \hspace{\stretch{1}}(1.7)

This gives us

\begin{aligned}Z &= \sqrt{ \frac{\pi}{\beta c} }\sqrt{\frac{2 \pi m}{\beta} }\left( {1 + f \beta \frac{3!!}{(2 \beta c)^2}+ \frac{1}{{2}} g^2 \beta^2 \frac{5!!}{(2 \beta c)^3}+ \frac{1}{{2}} f^2 \beta^2 \frac{7!!}{(2 \beta c)^4}+ \cdots} \right) \\ &= \frac{\pi}{\beta}\sqrt{ \frac{2 m}{c} }\left( {1 + \frac{3 f}{4 c^2 \beta}+ \frac{15 g^2}{16 c^3 \beta}+ \frac{105 f^2}{32 c^4 \beta}+ \cdots} \right).\end{aligned} \hspace{\stretch{1}}(1.8)

Retaining only the first two terms of the expansion, we have

\begin{aligned}\boxed{Z \approx\frac{\pi}{\beta}\sqrt{ \frac{2 m}{c} }\left( {1 + \frac{3 f}{4 c^2 \beta}+ \frac{15 g^2}{16 c^3 \beta}} \right).}\end{aligned} \hspace{\stretch{1}}(1.9)

Our average energy, in this approximation

\begin{aligned}\left\langle{{H}}\right\rangle &= -\frac{\partial {}}{\partial {\beta}} \ln Z \\ &\approx -\frac{\partial {}}{\partial {\beta}}\left( {- \ln \beta + \ln \left( 1 + \frac{3 f}{4 c^2 \beta}+ \frac{15 g^2}{16 c^3 \beta} \right) } \right) \\ &= \frac{1}{{\beta}} + \frac{1}{{\beta^2}} \frac{\frac{3 f}{4 c^2} + \frac{15 g^2}{16 c^3} }{1 + \frac{1}{{\beta}}\left( {\frac{3 f}{4 c^2} + \frac{15 g^2}{16 c^3} } \right)} \\ &= k_{\mathrm{B}} T + k_{\mathrm{B}}^2 T^2\left( { \frac{3 f}{4 c^2} + \frac{15 g^2}{16 c^3} } \right)\left( {1 - k_{\mathrm{B}} T \left( \frac{3 f}{4 c^2} + \frac{15 g^2}{16 c^3} \right) + \cdots} \right).\end{aligned} \hspace{\stretch{1}}(1.10)

So to first order in $T$ our specific heat is

\begin{aligned}C_{\mathrm{V}} = \frac{1}{{k_{\mathrm{B}}}} \frac{\partial {\left\langle{{H}}\right\rangle}}{\partial {T}}\approx1+ 2 k_{\mathrm{B}} T \left( {\frac{3 f}{4 c^2} + \frac{15 g^2}{16 c^3} } \right)\end{aligned} \hspace{\stretch{1}}(1.11)

or

\begin{aligned}\boxed{C_{\mathrm{V}} = 1+ k_{\mathrm{B}} T \left( {\frac{3 f}{2 c^2} + \frac{15 g^2}{8 c^3} } \right)+ \cdots}\end{aligned} \hspace{\stretch{1}}(1.12)

For the coordinate

\begin{aligned}\left\langle{{q}}\right\rangle &= \sqrt{\frac{2 \pi m}{\beta} }\frac{1}{{Z}}\int q e^{ -\beta \left( { c q^2 - g q^3 - f q^4} \right) } dq \\ &= \sqrt{\frac{2 \pi m}{\beta} }\frac{ \int q e^{ -\beta c q^2 } \left( { 1 + g \beta q^3 + f \beta q^4 } \right) dq}{ \frac{\pi}{\beta} \sqrt{ \frac{2 m}{c} } \left( { 1 + \frac{3 f}{4 c^2 \beta} + \frac{15 g^2}{16 c^3 \beta} } \right)} \\ &\approx \sqrt{\frac{\not{{2 \pi m}}}{\not{{\beta}}} }g \not{{\beta}} \frac{3!!}{(2 \beta c)^2} \sqrt{\frac{\not{{\pi}}}{\not{{\beta}} \not{{c}}}}\frac{1}{{ \frac{\not{{\pi}}}{\beta} \sqrt{ \frac{\not{{2 m}}}{\not{{c}}} }}},\end{aligned} \hspace{\stretch{1}}(1.13)

or

\begin{aligned}\boxed{\left\langle{{q}}\right\rangle \approx\frac{3 g k_{\mathrm{B}} T}{4 c^2}.}\end{aligned} \hspace{\stretch{1}}(1.14)

Compare this to the expectation of the coordinate for an unperturbed harmonic oscillator

\begin{aligned}\left\langle{{q}}\right\rangle = \frac{\int q e^{ -\beta c q^2}}{\int e^{-\beta c q^2}} = 0.\end{aligned} \hspace{\stretch{1}}(1.15)

We now have a temperature dependence to the expecation of the coordinate that we didn’t have for the harmonic oscillator.

# References

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

## Midterm II reflection

Posted by peeterjoot on March 10, 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.)]

Here’s some reflection about this Thursday’s midterm, redoing the problems without the mad scramble. I don’t think my results are too different from what I did in the midterm, even doing them casually now, but I’ll have to see after grading if these solutions are good.

## Question: Magnetic field spin level splitting (2013 midterm II p1)

A particle with spin $S$ has $2 S + 1$ states $-S, -S + 1, \cdots S-1, S$. When exposed to a magnetic field, state splitting results in energy $E_m = \hbar m B$. Calculate the partition function, and use this to find the temperature specific magnetization. A “sum the geometric series” hint was given.

Our partition function is

\begin{aligned}Z &= \sum_{m = -S}^S e^{-\hbar \beta m B} \\ &= e^{-\hbar \beta S B}\sum_{m = -S}^S e^{-\hbar \beta (m + S) B} \\ &= e^{\hbar \beta S B}\sum_{n = 0}^{2 S} e^{-\hbar \beta n B}.\end{aligned} \hspace{\stretch{1}}(1.0.1)

Writing

\begin{aligned}a = e^{-\hbar \beta B},\end{aligned} \hspace{\stretch{1}}(1.0.2)

that is

\begin{aligned}Z &= a^{-S}\sum_{n = 0}^{2 S} a^n \\ &= a^{-S} \frac{ a^{2 S + 1} - 1 }{a - 1} \\ &= \frac{ a^{S + 1} - a^{-S} }{a - 1} \\ &= \frac{ a^{S + 1/2} - a^{-S - 1/2} }{a^{1/2} - a^{-1/2}}.\end{aligned} \hspace{\stretch{1}}(1.0.3)

Substitution of $a$ gives us

\begin{aligned}\boxed{Z = \frac{ \sinh( \hbar \beta B (S + 1/2) ) }{ \sinh( \hbar \beta B /2 ) }.}\end{aligned} \hspace{\stretch{1}}(1.0.4)

To calculate the magnetization $M$, I used

\begin{aligned}M = -\left\langle{{H}}\right\rangle/B.\end{aligned} \hspace{\stretch{1}}(1.0.5)

As [1] defines magnetization for a spin system. It was pointed out to me after the test that magnetization was defined differently in class as

\begin{aligned}\mu = \frac{\partial {B}}{\partial {F}}.\end{aligned} \hspace{\stretch{1}}(1.0.6)

These are, up to a sign, identical, at least in this case, since we have $\beta$ and $B$ travelling together in the partition function. In terms of the average energy

\begin{aligned}M &= -\frac{\left\langle{{H}}\right\rangle}{B} \\ &= \frac{1}{{B}} \frac{\partial {}}{\partial {\beta}} \ln Z(\beta B) \\ &= \frac{1}{{Z B}} \frac{\partial {}}{\partial {\beta}}Z(\beta B) \\ &= \frac{1}{{Z}} \frac{\partial {}}{\partial {(\beta B)}} Z(\beta B)\end{aligned} \hspace{\stretch{1}}(1.0.7)

Compare this to the in-class definition of magnetization

\begin{aligned}\mu &= \frac{\partial {F}}{\partial {B}} \\ &= \frac{\partial {}}{\partial {B}} \left( - k_{\mathrm{B}} T \ln Z(\beta B) \right) \\ &= -\frac{\partial {}}{\partial {B}} \frac{\ln Z (\beta B)}{\beta} \\ &= -\frac{1}{{\beta Z}} \frac{\partial {}}{\partial {B}} Z(\beta B) \\ &= -\frac{1}{{Z}} \frac{\partial {}}{\partial {(\beta B)}} Z(\beta B).\end{aligned} \hspace{\stretch{1}}(1.0.8)

For this derivative we have

\begin{aligned}\frac{\partial {}}{\partial {(\beta B)}} \ln Z \\ &= \frac{\partial {}}{\partial {(\beta B)}} \ln \frac{ \sinh( \hbar \beta B (S + 1/2) ) }{ \sinh( \hbar \beta B /2 ) } \\ &= \frac{\partial {}}{\partial {(\beta B)}} \left( \ln \sinh( \hbar \beta B (S + 1/2) ) - \ln \sinh( \hbar \beta B /2 ) \right) \\ &= \frac{\hbar }{2}\left( (2 S + 1) \coth( \hbar \beta B (S + 1/2) ) - \coth( \hbar \beta B /2 ) \right).\end{aligned} \hspace{\stretch{1}}(1.0.9)

This gives us

\begin{aligned}\mu &= -\frac{1}{{Z}} \frac{\hbar }{2}\left( (2 S + 1) \coth( \hbar \beta B (S + 1/2) ) - \coth( \hbar \beta B /2 ) \right) \\ &= -\frac{ \sinh( \hbar \beta B /2 ) }{ \sinh( \hbar \beta B (S + 1/2) ) }\frac{\hbar }{2}\left( (2 S + 1) \coth( \hbar \beta B (S + 1/2) ) - \coth( \hbar \beta B /2 ) \right)\end{aligned} \hspace{\stretch{1}}(1.0.10)

After some simplification (done offline in \nbref{midtermTwoQ1FinalSimplificationMu.nb}) we get

\begin{aligned}\boxed{\mu = \hbar \frac{(s+1) \sinh(\hbar \beta B s)-s \sinh(\hbar \beta B (s+1))}{\cosh(\hbar \beta B(2 s+1)) - 1}.}\end{aligned} \hspace{\stretch{1}}(1.0.11)

I got something like this on the midterm, but recall doing it somehow much differently.

## Question: Pertubation of classical harmonic oscillator (2013 midterm II p2)

Consider a single particle perturbation of a classical simple harmonic oscillator Hamiltonian

\begin{aligned}H = \frac{1}{{2}} m \omega^2 \left( x^2 + y^2 \right) + \frac{1}{{2 m}} \left( p_x^2 + p_y^2 \right) + a x^4 + b y^6\end{aligned} \hspace{\stretch{1}}(1.0.12)

Calculate the canonical partition function, mean energy and specific heat of this system.

There were some instructions about the form to put the integrals in.

The canonical partition function is

\begin{aligned}Z &= \int dx dy dp_x dp_y e^{-\beta H} \\ &= \int dx e^{-\beta \left( \frac{1}{{2}} m \omega^2 x^2 + a x^4 \right)}\int dy e^{-\beta \left( \frac{1}{{2}} m \omega^2 y^2 + b y^6 \right)}\int dp_x dp_y e^{-\beta p_x^2/2 m} e^{-\beta p_y^2/2 m}.\end{aligned} \hspace{\stretch{1}}(1.0.13)

With

\begin{aligned}u = \sqrt{\frac{\beta}{2m}} p_x\end{aligned} \hspace{\stretch{1}}(1.0.14a)

\begin{aligned}v = \sqrt{\frac{\beta}{2m}} p_y,\end{aligned} \hspace{\stretch{1}}(1.0.14b)

the momentum integrals are

\begin{aligned}\int dp_x dp_y e^{-\beta p_x^2/2 m} e^{-\beta p_y^2/2 m} \\ &= \frac{2m}{\beta}\int du du e^{- u^2 - v^2} \\ &= \frac{m}{\beta}2 \pi\int 2 r dr e^{- r^2} \\ &= \frac{2 \pi m}{\beta}.\end{aligned} \hspace{\stretch{1}}(1.0.15)

Writing

\begin{aligned}f(x) = \frac{1}{{2}} m \omega^2 x^2 + a x^4\end{aligned} \hspace{\stretch{1}}(1.0.16a)

\begin{aligned}g(x) = \frac{1}{{2}} m \omega^2 y^2 + b y^4,\end{aligned} \hspace{\stretch{1}}(1.0.16b)

we have

\begin{aligned}\boxed{Z = \frac{2 \pi m}{\beta}\int dx e^{- \beta f(x)}\int dy e^{- \beta g(y)}.}\end{aligned} \hspace{\stretch{1}}(1.0.17)

The mean energy is

\begin{aligned}\left\langle{{H}}\right\rangle &= \frac{\int H e^{-\beta H}}{\int e^{-\beta H}} \\ &= -\frac{\partial {}}{\partial {\beta}} \ln \int e^{-\beta H} \\ &= \frac{\partial {}}{\partial {\beta}} \left( \ln \beta -\ln \int dx e^{- \beta f(x)} -\ln \int dy e^{- \beta g(y)} \right) \\ &= \frac{1}{{\beta}} + \frac{\int dx f(x) e^{- \beta f(x)}}{\int dx e^{- \beta f(x)}}+ \frac{\int dy g(y) e^{- \beta g(y)}}{\int dy e^{- \beta g(y)}}.\end{aligned} \hspace{\stretch{1}}(1.0.18)

The specific heat follows by differentiating once more

\begin{aligned}C_{\mathrm{V}} &= \frac{\partial {\left\langle{{H}}\right\rangle}}{\partial {T}} \\ &= \frac{\partial {\beta}}{\partial {T}}\frac{\partial {\left\langle{{H}}\right\rangle}}{\partial {\beta}} \\ &= -\frac{1}{{k_{\mathrm{B}} T^2}}\frac{\partial {\left\langle{{H}}\right\rangle}}{\partial {\beta}} \\ &= -k_{\mathrm{B}} \beta^2\frac{\partial {\left\langle{{H}}\right\rangle}}{\partial {\beta}} \\ &= - k_{\mathrm{B}} \beta^2\left( -\frac{1}{{\beta^2}} + \frac{\partial {}}{\partial {\beta}} \left( \frac{ \int dx f(x) e^{- \beta f(x)} } { \int dx e^{- \beta f(x)} } + \frac{ \int dy g(y) e^{- \beta g(y)} } { \int dy e^{- \beta g(y)} } \right) \right).\end{aligned} \hspace{\stretch{1}}(1.0.19)

Differentiating the integral terms we have, for example,

\begin{aligned}\frac{\partial {}}{\partial {\beta}} \frac{\int dx f(x) e^{- \beta f(x)}}{\int dx e^{- \beta f(x)}}=-\frac{\int dx f^2(x) e^{- \beta f(x)}}{\int dx e^{- \beta f(x)}}-\left( \frac{ \int dx f(x) e^{- \beta f(x)} } { \int dx e^{- \beta f(x)} } \right)^2,\end{aligned} \hspace{\stretch{1}}(1.0.20)

so that the specific heat is

\begin{aligned}\boxed{C_{\mathrm{V}} =k_{\mathrm{B}} \left(1 + \frac{\int dx f^2(x) e^{- \beta f(x)}}{\int dx e^{- \beta f(x)}}+\left( \frac{ \int dx f(x) e^{- \beta f(x)} } { \int dx e^{- \beta f(x)} } \right)^2+ \frac{\int dy g^2(y) e^{- \beta g(y)}}{\int dy e^{- \beta g(y)}}+\left( \frac{ \int dy g(y) e^{- \beta g(y)} } { \int dy e^{- \beta g(y)} } \right)^2\right).}\end{aligned} \hspace{\stretch{1}}(1.0.21)

That’s as far as I took this problem. There was a discussion after the midterm with Eric about Taylor expansion of these integrals. That’s not something that I tried.

# References

[1] C. Kittel and H. Kroemer. Thermal physics. WH Freeman, 1980.