Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

PHY452H1S Basic Statistical Mechanics. Problem Set 7: BEC and phonons

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


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

Question: Bose-Einstein condensation (BEC) in one and two dimensions

Obtain the density of states N(\epsilon) in one and two dimensions for a particle with an energy-momentum relation

\begin{aligned}E_\mathbf{k} = \frac{\hbar^2 \mathbf{k}^2}{2 m}.\end{aligned} \hspace{\stretch{1}}(1.1)

Using this, show that for particles whose number is conserved the BEC transition temperature vanishes in these cases – so we can always pick a chemical potential \mu < 0 which preserves a constant density at any temperature.


We’d like to evaluate

\begin{aligned}N_d(\epsilon) \equiv\sum_\mathbf{k}\delta(\epsilon - \epsilon_\mathbf{k})\approx\frac{L^d}{(2 \pi)^d} \int d^d \mathbf{k} \delta\left( \epsilon - \frac{\hbar^2 k^2}{2 m} \right),\end{aligned} \hspace{\stretch{1}}(1.2)

We’ll use

\begin{aligned}\delta(g(x)) = \sum_{x_0} \frac{\delta(x - x_0)}{\left\lvert {g'(x_0)} \right\rvert},\end{aligned} \hspace{\stretch{1}}(1.3)

where the roots of g(x) are x_0. With

\begin{aligned}g(k) = \epsilon - \frac{\hbar^2 k^2}{2 m},\end{aligned} \hspace{\stretch{1}}(1.4)

the roots k^{*} of g(k) = 0 are

\begin{aligned}k^{*} = \pm \sqrt{\frac{2 m \epsilon }{\hbar^2}}.\end{aligned} \hspace{\stretch{1}}(1.5)

The derivative of g(k) evaluated at these roots are

\begin{aligned}g'(k^{*}) &= -\frac{\hbar^2 k^{*}}{m} \\ &= \mp \frac{\hbar^2}{m}\frac{\sqrt{2 m \epsilon}}{ \hbar } \\ &= \mp \frac{\hbar \sqrt{2 m \epsilon} }{m}.\end{aligned} \hspace{\stretch{1}}(1.6)

In 2D, we can evaluate over a shell in k space

\begin{aligned}N_2(\epsilon) &= \frac{A}{(2 \pi)^2} \int_0^\infty 2 \pi k dk\left( \delta \left( k - k^{*}  \right) + \delta \left( k + k^{*}  \right)  \right)\frac{m}{\hbar \sqrt{2 m \epsilon} } \\ &= \frac{A}{2 \pi} \not{{k^{*}}}\frac{m}{\hbar^2 \not{{k^{*}}} }\end{aligned} \hspace{\stretch{1}}(1.7)


\begin{aligned}\boxed{N_2(\epsilon) = \frac{2 \pi A m}{h^2}.}\end{aligned} \hspace{\stretch{1}}(1.8)

In 1D we have

\begin{aligned}N_1(\epsilon) &= \frac{L}{2 \pi} \int_{-\infty}^\infty dk\left( \delta \left( k - k^{*}  \right) + \delta \left( k + k^{*}  \right)  \right)\frac{m}{\hbar \sqrt{2 m \epsilon} } \\ &= \frac{2 L}{2 \pi} \frac{m}{\hbar \sqrt{2 m \epsilon} }.\end{aligned} \hspace{\stretch{1}}(1.9)

Observe that this time for 1D, unlike in 2D when we used a radial shell in k space, we have contributions from both the delta function roots. Our end result is

\begin{aligned}\boxed{N_1(\epsilon) =\frac{2 L}{h} \sqrt{\frac{m}{2 \epsilon}}.}\end{aligned} \hspace{\stretch{1}}(1.10)

To consider the question of the BEC temperature, we’ll need to calculate the density. For the 2D case we have

\begin{aligned}\rho = \frac{N}{A} &= \frac{1}{A} A \int \frac{d^2 \mathbf{k}}{(2 \pi)^2} f(e_\mathbf{k}) \\ &= \frac{1}{A} \frac{2 \pi A m}{h^2}\int_0^\infty d\epsilon \frac{1}{{ z^{-1} e^{\beta \epsilon} -1 }} \\ &= \frac{2 \pi m}{h^2 \beta}\int_0^\infty dx \frac{1}{{ z^{-1} e^{x} -1 }} \\ &= -\frac{2 \pi m k_{\mathrm{B}} T}{h^2} \ln (1 - z) \\ &= -\frac{1}{{\lambda^2}} \ln (1 - z).\end{aligned} \hspace{\stretch{1}}(1.11)

Recall for the 3D case that we had an upper bound as z \rightarrow 1. We don’t have that for this 2D density, so for any value of k_{\mathrm{B}} T > 0, a corresponding value of z can be found. That is

\begin{aligned}z &= 1 - e^{-\rho \lambda^2} \\ &= 1 - e^{-\rho h^4/(2 \pi m k_{\mathrm{B}} T)^2}.\end{aligned} \hspace{\stretch{1}}(1.1.12)

For the 1D case we have

\begin{aligned}\rho &= \frac{N}{L} \\ &= \frac{1}{L} L \int \frac{dk}{2 \pi} f(e_\mathbf{k}) \\ &= \frac{1}{L} \frac{2 L}{h} \sqrt{\frac{m}{2}}\int_0^\infty d\epsilon \frac{1}{{\sqrt{\epsilon}}}\frac{1}{{ z^{-1} e^{\beta \epsilon} -1 }} \\ &= \frac{1}{{h}} \sqrt{\frac{2 m}{\beta}} \int_0^\infty \frac{x^{1/2 - 1}}{z^{-1} e^x - 1} \\ &= \frac{1}{{h}} \sqrt{\frac{2 m}{\beta}} \Gamma(1/2) f^-_{1/2}(z),\end{aligned} \hspace{\stretch{1}}(1.1.12)


\begin{aligned}\rho= \frac{1}{{\lambda}} f^-_{1/2}(z).\end{aligned} \hspace{\stretch{1}}(1.1.12)

See fig. 1.1 for plots of f^-_\nu(z) for \nu \in \{1/2, 1, 3/2\}, the respective results for the 1D, 2D and 3D densities respectively.

Fig 1.1: Density integrals for 1D, 2D and 3D cases

We’ve found that f^-_{1/2}(z) is also unbounded as z \rightarrow 1, so while we cannot invert this easily as in the 2D case, we can at least say that there will be some z for any value of k_{\mathrm{B}} T > 0 that allows the density (and thus the number of particles) to remain fixed.

Question: Estimating the BEC transition temperature

Find data for the atomic mass of liquid {}^4 He and its density at ambient atmospheric pressure and hence estimate its BEC temperature assuming interactions are unimportant (even though this assumption is a very bad one!).

For dilute atomic gases of the sort used in Professor
Thywissen’s lab
, one typically has a cloud of 10^6 atoms confined to an approximate cubic region with linear dimension 1 \mu\,m. Find the density – it is pretty low, so interactions can be assumed to be extremely weak. Assuming these are {}^{87} Rb atoms, estimate the BEC transition temperature.


With an atomic weight of 4.0026, the mass in grams for one atom of Helium is

\begin{aligned}4.0026 \,\text{amu} \times \frac{\text{g}}{6.022 \times 10^{23} \text{amu}} &= 6.64 \times 10^{-24} \text{g} \\ &= 6.64 \times 10^{-27} \text{kg}.\end{aligned} \hspace{\stretch{1}}(1.15)

With the density of liquid He-4, at 5.2K (boiling point): 125 grams per liter, the number density is

\begin{aligned}\rho &= \frac{\text{mass}}{\text{volume}} \times \frac{1}{{\text{mass of one He atom}}} \\ &= \frac{125 \text{g}}{10^{-3} m^3} \times \frac{1}{{6.64 \times 10^{-24} g}} \\ &= \frac{125 \text{g}}{10^{-3} m^3} \times \frac{1}{{6.64 \times 10^{-24} g}} \\ &= 1.88 \times 10^{28} m^{-3}\end{aligned} \hspace{\stretch{1}}(1.16)

In class the T_{\mathrm{BEC}} was found to be

\begin{aligned}T_{\mathrm{BEC}} &= \frac{1}{k_{\mathrm{B}}} \left( \frac{\rho}{\zeta(3/2)}  \right)^{2/3} \frac{ 2 \pi \hbar^2}{M} \\ &= \frac{1}{{1.3806488 \times 10^{-23} m^2 kg/s^2/K}} \left( \frac{\rho}{ 2.61238 }  \right)^{2/3} \frac{ 2 \pi (1.05457173 \times 10^{-34} m^2 kg / s)^2}{M} \\ &= 2.66824 \times 10^{-45} \frac{\rho^{2/3}}{M} K.\end{aligned} \hspace{\stretch{1}}(1.17)

So for liquid helium we have

\begin{aligned}T_{\mathrm{BEC}} &= 2.66824 \times 10^{-45} \left( 1.88 \times 10^{28}  \right)^{2/3} \frac{1}{{ 6.64 \times 10^{-27} }} K \\ &= 2.84 K.\end{aligned} \hspace{\stretch{1}}(1.18)

The number density for the gas in Thywissen’s lab is

\begin{aligned}\rho &= \frac{10^6}{(10^{-6} \text{m})^3} \\ &= 10^{24} m^{-3}.\end{aligned} \hspace{\stretch{1}}(1.1.19)

The mass of an atom of {}^{87} Rb is

\begin{aligned}86.90 \,\text{amu} \times \frac{10^{-3} \text{kg}}{6.022 \times 10^{23} \text{amu}} = 1.443 \times 10^{-25} \text{kg},\end{aligned} \hspace{\stretch{1}}(1.1.19)

which gives us

\begin{aligned}T_{\mathrm{BEC}} &= 2.66824 \times 10^{-45} \left( 10^{24}  \right)^{2/3} \frac{1}{{ 1.443 \times 10^{-25} }} K \\ &= 1.85 \times 10^{-4} K.\end{aligned} \hspace{\stretch{1}}(1.1.19)

Question: Phonons in two dimensions

Consider phonons (quanta of lattice vibrations) which obey a dispersion relation

\begin{aligned}E_\mathbf{k} = \hbar v \left\lvert {\mathbf{k}} \right\rvert\end{aligned} \hspace{\stretch{1}}(1.1.22)

for small momenta \left\lvert {\mathbf{k}} \right\rvert, where v is the speed of sound. Assuming a two-dimensional crystal, phonons only propagate along the plane containing the atoms. Find the specific heat of this crystal due to phonons at low temperature. Recall that phonons are not conserved, so there is no chemical potential associated with maintaining a fixed phonon density.

The energy density of the system is

\begin{aligned}\frac{E}{V} &= \int \frac{d^2 \mathbf{k}}{(2 \pi)^2} \frac{\epsilon}{ e^{\beta \epsilon} - 1 } \\ &= \int d\epsilon \frac{N(\epsilon)}{V} \frac{\epsilon}{ e^{\beta \epsilon} - 1 }.\end{aligned} \hspace{\stretch{1}}(1.23)

For the density of states we have

\begin{aligned}\frac{N(\epsilon) }{V} &= \int \frac{d^2 \mathbf{k}}{(2 \pi)^2} \delta( \epsilon - \epsilon_\mathbf{k} ) \\ &= \frac{1}{{(2 \pi)^2}} 2 \pi \int_0^\infty k dk \delta( \epsilon - \hbar v k ) \\ &= \frac{1}{{2 \pi}} \int_0^\infty k dk \delta \left( k - \frac{\epsilon}{\hbar v}  \right) \frac{1}{{\hbar v}} \\ &= \frac{1}{{2 \pi}} \frac{\epsilon}{(\hbar v)^2}.\end{aligned} \hspace{\stretch{1}}(1.24)

Plugging back into the energy density we have

\begin{aligned}\frac{E}{V} &= \frac{2 \pi}{(\hbar v)^2}\int_0^\infty d\epsilon \frac{\epsilon^2}{ e^{\beta \epsilon} - 1 } \\ &= \frac{\pi \left( k_{\mathrm{B}} T \right)^3 }{(\hbar v)^2}\zeta(3),\end{aligned} \hspace{\stretch{1}}(1.25)

where \zeta(3) \approx 2.40411. Taking derivatives we have

\begin{aligned}\boxed{C_V = \frac{dE}{dT} = V\frac{3 \pi k_{\mathrm{B}}^3 T^2 }{(\hbar v)^2}\zeta(3).}\end{aligned} \hspace{\stretch{1}}(1.1.26)


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: