• 360,180

# Posts Tagged ‘phonon’

Posted by peeterjoot on January 19, 2014

## Question: Quadratic Debye phonons (2013 midterm pr B2)

Assume a quadratic dispersion relation for the longitudinal and transverse modes

\begin{aligned}\omega = \left\{\begin{array}{l}b_{\mathrm{L}} q^2 \\ b_{\mathrm{T}} q^2\end{array}\right..\end{aligned} \hspace{\stretch{1}}(1.2)

### Part a

Find the density of states.

### Part b

Find the Debye frequency.

### Part c

In terms of $k_{\mathrm{B}} \Theta = \hbar \omega_{\mathrm{D}}$, and

\begin{aligned}\mathcal{I} = \int_0^\infty \frac{y^{5/2} e^{y} dy}{\left( { e^y - 1} \right)^2 },\end{aligned} \hspace{\stretch{1}}(1.2)

find the specific heat for $k_{\mathrm{B}} T \ll \hbar \omega_{\mathrm{D}}$.

### Part d

Find the specific heat for $k_{\mathrm{B}} T \gg \hbar \omega_{\mathrm{D}}$.

### Part a

Working straight from the definition

\begin{aligned}Z(\omega) &= \frac{V}{(2 \pi)^3 } \sum_{L, T} \int \frac{df_\omega}{ \left\lvert { \boldsymbol{\nabla}_\mathbf{q} \omega } \right\rvert } \\ &= \frac{V}{(2 \pi)^3 } \left( { {\left.{{\frac{4 \pi q^2}{2 b_{\mathrm{L}} q} }}\right\vert}_{{\mathrm{L}}} + {\left.{{\frac{2 \times 4 \pi q^2}{2 b_{\mathrm{T}} q} }}\right\vert}_{{\mathrm{T}}} } \right) \\ &= \frac{V}{4 \pi^2 } \left( { \frac{q_{\mathrm{L}}}{b_{\mathrm{L}}} + \frac{2 q_{\mathrm{T}}}{b_{\mathrm{T}}} } \right).\end{aligned} \hspace{\stretch{1}}(1.3)

With $q_{\mathrm{L}} = \sqrt{\omega/b_{\mathrm{L}}}$ and $q_{\mathrm{T}} = \sqrt{\omega/b_{\mathrm{T}}}$, this is

\begin{aligned}Z(\omega) = \frac{V}{4 \pi^2 } \left( { \frac{1}{b_{\mathrm{L}}^{3/2}} + \frac{2}{b_{\mathrm{T}}^{3/2}} } \right)\sqrt{\omega}\end{aligned} \hspace{\stretch{1}}(1.4)

### Part b

The Debye frequency was given implicitly by

\begin{aligned}\int_0^{\omega_{\mathrm{D}}} Z(\omega) d\omega = 3 r N,\end{aligned} \hspace{\stretch{1}}(1.5)

which gives

\begin{aligned}3 r N=\frac{2}{3} \frac{V}{4 \pi^2 } \left( { \frac{1}{b_{\mathrm{L}}^{3/2}} + \frac{2}{b_{\mathrm{T}}^{3/2}} } \right)\omega_{\mathrm{D}}^{3/2}=\frac{V}{6 \pi^2 } \left( { \frac{1}{b_{\mathrm{L}}^{3/2}} + \frac{2}{b_{\mathrm{T}}^{3/2}} } \right)\omega_{\mathrm{D}}^{3/2}\end{aligned} \hspace{\stretch{1}}(1.6)

### Part c

Assuming a Bose distribution and ignoring the zero point energy, which has no temperature dependence, the specific heat, the temperature derivative of the energy density, is

\begin{aligned}C_{\mathrm{V}} &= \frac{d}{d T} \frac{1}{{V}} \int Z(\omega) \frac{\hbar \omega}{ e^{\hbar \omega/ k_{\mathrm{B}} T } - 1} d\omega \\ &= \frac{1}{{V}} \frac{d}{d T} \int Z(\omega) \frac{\hbar \omega}{ \hbar \omega/ k_{\mathrm{B}} T + \frac{1}{{2}}( \hbar \omega/k_{\mathrm{B}} T)^2 + \cdots } d\omega \\ &\approx \frac{1}{{V}} \frac{d}{d T} \int Z(\omega) k_{\mathrm{B}} T d\omega \\ &= \frac{1}{{V}} k_{\mathrm{B}} 3 r N.\end{aligned} \hspace{\stretch{1}}(1.7)

### Part d

First note that the density of states can be written

\begin{aligned}Z(\omega) = \frac{9 r N}{ 2 \omega_{\mathrm{D}}^{3/2} } \omega^{1/2},\end{aligned} \hspace{\stretch{1}}(1.8)

for a specific heat of

\begin{aligned}C_{\mathrm{V}} &= \frac{d}{d T} \frac{1}{{V}} \int_0^\infty \frac{9 r N}{ 2 \omega_{\mathrm{D}}^{3/2} } \omega^{1/2} \frac{\hbar \omega}{ e^{\hbar \omega/ k_{\mathrm{B}} T } - 1} d\omega \\ &= \frac{9 r N}{ 2 V \omega_{\mathrm{D}}^{3/2} } \int_0^\infty d\omega \omega^{1/2} \frac{d}{d T} \frac{\hbar \omega}{ e^{\hbar \omega/ k_{\mathrm{B}} T } - 1} \\ &= \frac{9 r N}{ 2 V \omega_{\mathrm{D}}^{3/2} } \int_0^\infty d\omega \omega^{1/2} \frac{-\hbar \omega}{ \left( {e^{\hbar \omega/ k_{\mathrm{B}} T } - 1} \right)^2 } e^{\hbar \omega/k_{\mathrm{B}} T} \hbar \omega/k_{\mathrm{B}} \left( {-\frac{1}{{T^2}}} \right) \\ &= \frac{9 r N k_{\mathrm{B}} }{ 2 V \omega_{\mathrm{D}}^{3/2} } \left( { \frac{ k_{\mathrm{B}} T}{\hbar} } \right)^{3/2}\int_0^\infty d \frac{\hbar \omega}{k_{\mathrm{B}} T} \left( {\frac{\hbar \omega}{k_{\mathrm{B}} T}} \right)^{1/2} \frac{1}{ \left( {e^{\hbar \omega/ k_{\mathrm{B}} T } - 1} \right)^2 } e^{\hbar \omega/k_{\mathrm{B}} T} \left( { \frac{\hbar \omega}{k_{\mathrm{B}} T} } \right)^2 \\ &= \frac{9 r N k_{\mathrm{B}} }{ 2 V \omega_{\mathrm{D}}^{3/2} } \left( { \frac{ k_{\mathrm{B}} T}{\hbar} } \right)^{3/2}\int_0^\infty dy \frac{y^{5/2} e^y }{ \left( {e^y - 1} \right)^2 } \\& = \frac{9 r N k_{\mathrm{B}} }{ 2 V } \left( { \frac{ T}{\Theta} } \right)^{3/2} \mathcal{I}.\end{aligned} \hspace{\stretch{1}}(1.9)

## One atom basis phonons in 2D

Posted by peeterjoot on January 19, 2014

Let’s tackle a problem like the 2D problem of the final exam, but first more generally. Instead of a square lattice consider the lattice with the geometry illustrated in fig. 1.1.

Fig 1.1: Oblique one atom basis

Here, $\mathbf{a}$ and $\mathbf{b}$ are the vector differences between the equilibrium positions separating the masses along the $K_1$ and $K_2$ interaction directions respectively. The equilibrium spacing for the cross coupling harmonic forces are

\begin{aligned}\begin{aligned}\mathbf{r} &= (\mathbf{b} + \mathbf{a})/2 \\ \mathbf{s} &= (\mathbf{b} - \mathbf{a})/2.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.1)

Based on previous calculations, we can write the equations of motion by inspection

\begin{aligned}\begin{aligned}m \dot{d}{\mathbf{u}}_\mathbf{n} = &-K_1 \text{Proj}_{\hat{\mathbf{a}}} \sum_\pm \left( { \mathbf{u}_\mathbf{n} - \mathbf{u}_{\mathbf{n} \pm(1, 0)}} \right)^2 \\ &-K_2 \text{Proj}_{\hat{\mathbf{b}}} \sum_\pm \left( { \mathbf{u}_\mathbf{n} - \mathbf{u}_{\mathbf{n} \pm(0, 1)}} \right)^2 \\ &-K_3 \text{Proj}_{\hat{\mathbf{r}}} \sum_\pm \left( { \mathbf{u}_\mathbf{n} - \mathbf{u}_{\mathbf{n} \pm(1, 1)}} \right)^2 \\ &-K_4 \text{Proj}_{\hat{\mathbf{s}}} \sum_\pm \left( { \mathbf{u}_\mathbf{n} - \mathbf{u}_{\mathbf{n} \pm(1, -1)}} \right)^2.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.2)

Inserting the trial solution

\begin{aligned}\mathbf{u}_\mathbf{n} = \frac{1}{{\sqrt{m}}} \boldsymbol{\epsilon}(\mathbf{q}) e^{i( \mathbf{r}_\mathbf{n} \cdot \mathbf{q} - \omega t) },\end{aligned} \hspace{\stretch{1}}(1.3)

and using the matrix form for the projection operators, we have

\begin{aligned}\begin{aligned}\omega^2 \boldsymbol{\epsilon} &=\frac{K_1}{m} \hat{\mathbf{a}} \hat{\mathbf{a}}^\text{T} \boldsymbol{\epsilon}\sum_\pm\left( { 1 - e^{\pm i \mathbf{a} \cdot \mathbf{q}} } \right) \\ & +\frac{K_2}{m} \hat{\mathbf{b}} \hat{\mathbf{b}}^\text{T} \boldsymbol{\epsilon}\sum_\pm\left( { 1 - e^{\pm i \mathbf{b} \cdot \mathbf{q}} } \right) \\ & +\frac{K_3}{m} \hat{\mathbf{b}} \hat{\mathbf{b}}^\text{T} \boldsymbol{\epsilon}\sum_\pm\left( { 1 - e^{\pm i (\mathbf{b} + \mathbf{a}) \cdot \mathbf{q}} } \right) \\ & +\frac{K_3}{m} \hat{\mathbf{b}} \hat{\mathbf{b}}^\text{T} \boldsymbol{\epsilon}\sum_\pm\left( { 1 - e^{\pm i (\mathbf{b} - \mathbf{a}) \cdot \mathbf{q}} } \right) \\ &=\frac{4 K_1}{m} \hat{\mathbf{a}} \hat{\mathbf{a}}^\text{T} \boldsymbol{\epsilon} \sin^2\left( { \mathbf{a} \cdot \mathbf{q}/2 } \right)+\frac{4 K_2}{m} \hat{\mathbf{b}} \hat{\mathbf{b}}^\text{T} \boldsymbol{\epsilon} \sin^2\left( { \mathbf{b} \cdot \mathbf{q}/2 } \right) \\ &+\frac{4 K_3}{m} \hat{\mathbf{r}} \hat{\mathbf{r}}^\text{T} \boldsymbol{\epsilon} \sin^2\left( { (\mathbf{b} + \mathbf{a}) \cdot \mathbf{q}/2 } \right)+\frac{4 K_4}{m} \hat{\mathbf{s}} \hat{\mathbf{s}}^\text{T} \boldsymbol{\epsilon} \sin^2\left( { (\mathbf{b} - \mathbf{a}) \cdot \mathbf{q}/2 } \right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.4)

This fully specifies our eigenvalue problem. Writing

\begin{aligned}\begin{aligned}S_1 &= \sin^2\left( { \mathbf{a} \cdot \mathbf{q}/2 } \right) \\ S_2 &= \sin^2\left( { \mathbf{b} \cdot \mathbf{q}/2 } \right) \\ S_3 &= \sin^2\left( { (\mathbf{b} + \mathbf{a}) \cdot \mathbf{q}/2 } \right) \\ S_4 &= \sin^2\left( { (\mathbf{b} - \mathbf{a}) \cdot \mathbf{q}/2 } \right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.5.5)

\begin{aligned}\boxed{A = \frac{4}{m}\left( { K_1 S_1 \hat{\mathbf{a}} \hat{\mathbf{a}}^\text{T} + K_2 S_2 \hat{\mathbf{b}} \hat{\mathbf{b}}^\text{T} + K_3 S_3 \hat{\mathbf{r}} \hat{\mathbf{r}}^\text{T} + K_4 S_4 \hat{\mathbf{s}} \hat{\mathbf{s}}^\text{T}} \right),}\end{aligned} \hspace{\stretch{1}}(1.0.5.5)

we wish to solve

\begin{aligned}A \boldsymbol{\epsilon} = \omega^2 \boldsymbol{\epsilon} = \lambda \boldsymbol{\epsilon}.\end{aligned} \hspace{\stretch{1}}(1.0.6)

Neglecting the specifics of the matrix at hand, consider a generic two by two matrix

\begin{aligned}A = \begin{bmatrix}a & b \\ c & d\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.6)

for which the characteristic equation is

\begin{aligned}0 &= \begin{vmatrix}\lambda - a & - b \\ -c & \lambda -d \end{vmatrix} \\ &= (\lambda - a)(\lambda - d) - b c \\ &= \lambda^2 - (a + d) \lambda + a d - b c \\ &= \lambda^2 - (Tr A) \lambda + \left\lvert {A} \right\rvert \\ &= \left( {\lambda - \frac{Tr A}{2}} \right)^2- \left( {\frac{Tr A}{2}} \right)^2 + \left\lvert {A} \right\rvert.\end{aligned} \hspace{\stretch{1}}(1.0.6)

So our angular frequencies are given by

\begin{aligned}\omega^2 = \frac{1}{{2}} \left( { Tr A \pm \sqrt{ \left(Tr A\right)^2 - 4 \left\lvert {A} \right\rvert }} \right).\end{aligned} \hspace{\stretch{1}}(1.0.6)

The square root can be simplified slightly

\begin{aligned}\left( {Tr A} \right)^2 - 4 \left\lvert {A} \right\rvert \\ &= (a + d)^2 -4 (a d - b c) \\ &= a^2 + d^2 + 2 a d - 4 a d + 4 b c \\ &= (a - d)^2 + 4 b c,\end{aligned} \hspace{\stretch{1}}(1.0.6)

so that, finally, the dispersion relation is

\begin{aligned}\boxed{\omega^2 = \frac{1}{{2}} \left( { d + a \pm \sqrt{ (d - a)^2 + 4 b c } } \right),}\end{aligned} \hspace{\stretch{1}}(1.0.6)

Our eigenvectors will be given by

\begin{aligned}0 = (\lambda - a) \boldsymbol{\epsilon}_1 - b\boldsymbol{\epsilon}_2,\end{aligned} \hspace{\stretch{1}}(1.0.6)

or

\begin{aligned}\boldsymbol{\epsilon}_1 \propto \frac{b}{\lambda - a}\boldsymbol{\epsilon}_2.\end{aligned} \hspace{\stretch{1}}(1.0.6)

So, our eigenvectors, the vectoral components of our atomic displacements, are

\begin{aligned}\boldsymbol{\epsilon} \propto\begin{bmatrix}b \\ \omega^2 - a\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.6)

or

\begin{aligned}\boxed{\boldsymbol{\epsilon} \propto\begin{bmatrix}2 b \\ d - a \pm \sqrt{ (d - a)^2 + 4 b c }\end{bmatrix}.}\end{aligned} \hspace{\stretch{1}}(1.0.6)

Square lattice

There is not too much to gain by expanding out the projection operators explicitly in general. However, let’s do this for the specific case of a square lattice (as on the exam problem). In that case, our projection operators are

\begin{aligned}\hat{\mathbf{a}} \hat{\mathbf{a}}^\text{T} = \begin{bmatrix}1 \\ 0\end{bmatrix}\begin{bmatrix}1 & 0\end{bmatrix}=\begin{bmatrix}1 & 0 \\ 0 & 0\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.16a)

\begin{aligned}\hat{\mathbf{b}} \hat{\mathbf{b}}^\text{T} = \begin{bmatrix}0\\ 1 \end{bmatrix}\begin{bmatrix}0 &1 \end{bmatrix}=\begin{bmatrix}0 & 0 \\ 0 & 1\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.16b)

\begin{aligned}\hat{\mathbf{r}} \hat{\mathbf{r}}^\text{T} = \frac{1}{{2}}\begin{bmatrix}1 \\ 1 \end{bmatrix}\begin{bmatrix}1 &1 \end{bmatrix}=\frac{1}{{2}}\begin{bmatrix}1 & 1 \\ 1 & 1\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.16c)

\begin{aligned}\hat{\mathbf{s}} \hat{\mathbf{s}}^\text{T} = \frac{1}{{2}}\begin{bmatrix}-1 \\ 1 \end{bmatrix}\begin{bmatrix}-1 &1 \end{bmatrix}=\frac{1}{{2}}\begin{bmatrix}1 & -1 \\ -1 & 1\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.16d)

\begin{aligned}\begin{aligned}S_1 &= \sin^2\left( { \mathbf{a} \cdot \mathbf{q} } \right) \\ S_2 &= \sin^2\left( { \mathbf{b} \cdot \mathbf{q} } \right) \\ S_3 &= \sin^2\left( { (\mathbf{b} + \mathbf{a}) \cdot \mathbf{q} } \right) \\ S_4 &= \sin^2\left( { (\mathbf{b} - \mathbf{a}) \cdot \mathbf{q} } \right),\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.16d)

Our matrix is

\begin{aligned}A = \frac{2}{m}\begin{bmatrix}2 K_1 S_1 + K_3 S_3 + K_4 S_4 & K_3 S_3 - K_4 S_4 \\ K_3 S_3 - K_4 S_4 & 2 K_2 S_2 + K_3 S_3 + K_4 S_4\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.16d)

where, specifically, the squared sines for this geometry are

\begin{aligned}S_1 = \sin^2 \left( { \mathbf{a} \cdot \mathbf{q}/2 } \right) = \sin^2 \left( { a q_x/2} \right)\end{aligned} \hspace{\stretch{1}}(1.0.19a)

\begin{aligned}S_2 = \sin^2 \left( { \mathbf{b} \cdot \mathbf{q}/2 } \right) = \sin^2 \left( { a q_y/2} \right)\end{aligned} \hspace{\stretch{1}}(1.0.19b)

\begin{aligned}S_3 = \sin^2 \left( { (\mathbf{b} + \mathbf{a}) \cdot \mathbf{q}/2 } \right) = \sin^2 \left( { a (q_x + q_y)/2} \right)\end{aligned} \hspace{\stretch{1}}(1.0.19c)

\begin{aligned}S_4 = \sin^2 \left( { (\mathbf{b} - \mathbf{a}) \cdot \mathbf{q}/2 } \right) = \sin^2 \left( { a (q_y - q_x)/2} \right).\end{aligned} \hspace{\stretch{1}}(1.0.19d)

Using eq. 1.0.6, the dispersion relation and eigenvectors are

\begin{aligned}\omega^2 = \frac{2}{m} \left( { \sum_i K_i S_i \pm \sqrt{ (K_2 S_2 - K_1 S_1)^2 + (K_3 S_3 - K_4 S_4)^2 } } \right)\end{aligned} \hspace{\stretch{1}}(1.0.20.20)

\begin{aligned}\boldsymbol{\epsilon} \propto\begin{bmatrix}K_3 S_3 - K_4 S_4 \\ K_2 S_2 - K_1 S_1 \pm \sqrt{ (K_2 S_2 - K_1 S_1)^2 + (K_3 S_3 - K_4 S_4)^2 } \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.20.20)

This calculation is confirmed in oneAtomBasisPhononSquareLatticeEigensystem.nb. Mathematica calculates an alternate form (equivalent to using a zero dot product for the second row), of

\begin{aligned}\boldsymbol{\epsilon} \propto\begin{bmatrix}K_1 S_1 - K_2 S_2 \pm \sqrt{ (K_2 S_2 - K_1 S_1)^2 + (K_3 S_3 - K_4 S_4)^2 } \\ K_3 S_3 - K_4 S_4 \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.20.20)

Either way, we see that $K_3 S_3 - K_4 S_4 = 0$ leads to only horizontal or vertical motion.

With the exam criteria

In the specific case that we had on the exam where $K_1 = K_2$ and $K_3 = K_4$, these are

\begin{aligned}\omega^2 = \frac{2}{m} \left( { K_1 (S_1 + S_2) + K_3(S_3 + S_4) \pm \sqrt{ K_1^2 (S_2 - S_1)^2 + K_3^2 (S_3 - S_4)^2 } } \right)\end{aligned} \hspace{\stretch{1}}(1.0.22.22)

\begin{aligned}\boldsymbol{\epsilon} \propto\begin{bmatrix}K_3 \left( { S_3 - S_4 } \right) \\ K_1 \left( { (S_1 - S_2) \pm \sqrt{ (S_2 - S_1)^2 + \left( \frac{K_3}{K_1} \right)^2 (S_3 - S_4)^2 } } \right)\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.22.22)

For horizontal and vertical motion we need $S_3 = S_4$, or for a $2 \pi \times \text{integer}$ difference in the absolute values of the sine arguments

\begin{aligned}\pm ( a (q_x + q_y) /2 ) = a (q_y - q_y) /2 + 2 \pi n.\end{aligned} \hspace{\stretch{1}}(1.0.22.22)

That is, one of

\begin{aligned}\begin{aligned}q_x &= \frac{2 \pi}{a} n \\ q_y &= \frac{2 \pi}{a} n\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.22.22)

In the first BZ, that is one of $q_x = 0$ or $q_y = 0$.

System in rotated coordinates

On the exam, where we were asked to solve for motion along the cross directions explicitly, there was a strong hint to consider a rotated (by $\pi/4$) coordinate system.

The rotated the lattice basis vectors are $\mathbf{a} = a \mathbf{e}_1, \mathbf{b} = a \mathbf{e}_2$, and the projection matrices. Writing $\hat{\mathbf{r}} = \mathbf{f}_1$ and $\hat{\mathbf{s}} = \mathbf{f}_2$, where $\mathbf{f}_1 = (\mathbf{e}_1 + \mathbf{e}_2)/\sqrt{2}, \mathbf{f}_2 = (\mathbf{e}_2 - \mathbf{e}_1)/\sqrt{2}$, or $\mathbf{e}_1 = (\mathbf{f}_1 - \mathbf{f}_2)/\sqrt{2}, \mathbf{e}_2 = (\mathbf{f}_1 + \mathbf{f}_2)/\sqrt{2}$. In the $\{\mathbf{f}_1, \mathbf{f}_2\}$ basis the projection matrices are

\begin{aligned}\hat{\mathbf{a}} \hat{\mathbf{a}}^\text{T} = \frac{1}{{2}}\begin{bmatrix}1 \\ -1\end{bmatrix}\begin{bmatrix}1 & -1\end{bmatrix}= \frac{1}{{2}} \begin{bmatrix}1 & -1 \\ -1 & 1\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.25a)

\begin{aligned}\hat{\mathbf{b}} \hat{\mathbf{b}}^\text{T} = \frac{1}{{2}}\begin{bmatrix}1 \\ 1\end{bmatrix}\begin{bmatrix}1 & 1\end{bmatrix}= \frac{1}{{2}} \begin{bmatrix}1 & 1 \\ 1 & 1\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.25b)

\begin{aligned}\hat{\mathbf{r}} \hat{\mathbf{r}}^\text{T} = \begin{bmatrix}1 & 0 \\ 0 & 0\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.25c)

\begin{aligned}\hat{\mathbf{s}} \hat{\mathbf{s}}^\text{T} = \begin{bmatrix}0 & 0 \\ 0 & 1\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.25d)

The dot products that show up in the squared sines are

\begin{aligned}\mathbf{a} \cdot \mathbf{q}=a \frac{1}{{\sqrt{2}}} (\mathbf{f}_1 - \mathbf{f}_2) \cdot (\mathbf{f}_1 k_u + \mathbf{f}_2 k_v)=\frac{a}{\sqrt{2}} (k_u - k_v)\end{aligned} \hspace{\stretch{1}}(1.0.26a)

\begin{aligned}\mathbf{b} \cdot \mathbf{q}=a \frac{1}{{\sqrt{2}}} (\mathbf{f}_1 + \mathbf{f}_2) \cdot (\mathbf{f}_1 k_u + \mathbf{f}_2 k_v)=\frac{a}{\sqrt{2}} (k_u + k_v)\end{aligned} \hspace{\stretch{1}}(1.0.26b)

\begin{aligned}(\mathbf{a} + \mathbf{b}) \cdot \mathbf{q} = \sqrt{2} a k_u \end{aligned} \hspace{\stretch{1}}(1.0.26c)

\begin{aligned}(\mathbf{b} - \mathbf{a}) \cdot \mathbf{q} = \sqrt{2} a k_v \end{aligned} \hspace{\stretch{1}}(1.0.26d)

So that in this basis

\begin{aligned}\begin{aligned}S_1 &= \sin^2 \left( { \frac{a}{\sqrt{2}} (k_u - k_v) } \right) \\ S_2 &= \sin^2 \left( { \frac{a}{\sqrt{2}} (k_u + k_v) } \right) \\ S_3 &= \sin^2 \left( { \sqrt{2} a k_u } \right) \\ S_4 &= \sin^2 \left( { \sqrt{2} a k_v } \right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.26d)

With the rotated projection operators eq. 1.0.5.5 takes the form

\begin{aligned}A = \frac{2}{m}\begin{bmatrix}K_1 S_1 + K_2 S_2 + 2 K_3 S_3 & K_2 S_2 - K_1 S_1 \\ K_2 S_2 - K_1 S_1 & K_1 S_1 + K_2 S_2 + 2 K_4 S_4\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.26d)

This clearly differs from eq. 1.0.16d, and results in a different expression for the eigenvectors, but the same as eq. 1.0.20.20 for the angular frequencies.

\begin{aligned}\boldsymbol{\epsilon} \propto\begin{bmatrix}K_2 S_2 - K_1 S_1 \\ K_4 S_4 - K_3 S_3 \mp \sqrt{ (K_2 S_2 - K_1 S_1)^2 + (K_3 S_3 - K_4 S_4)^2 }\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.26d)

or, equivalently

\begin{aligned}\boldsymbol{\epsilon} \propto\begin{bmatrix}K_4 S_4 - K_3 S_3 \mp \sqrt{ (K_2 S_2 - K_1 S_1)^2 + (K_3 S_3 - K_4 S_4)^2 } \\ K_1 S_1 - K_2 S_2 \\ \end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.26d)

For the $K_1 = K_2$ and $K_3 = K_4$ case of the exam, this is

\begin{aligned}\boldsymbol{\epsilon} \propto\begin{bmatrix}K_1 (S_2 - S_1 ) \\ K_3 \left( { S_4 - S_3 \mp \sqrt{ \left( \frac{K_1}{K_3} \right)^2 (S_2 - S_1)^2 + (S_3 - S_4)^2 } } \right)\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.26d)

Similar to the horizontal coordinate system, we see that we have motion along the diagonals when

\begin{aligned}\pm \frac{a}{\sqrt{2}} (k_u - k_v) = \frac{a}{\sqrt{2}} (k_u + k_v) + 2 \pi n,\end{aligned} \hspace{\stretch{1}}(1.0.26d)

or one of

\begin{aligned}\begin{aligned}k_u &= \sqrt{2} \frac{\pi}{a} n \\ k_v &= \sqrt{2} \frac{\pi}{a} n\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.26d)

Stability?

The exam asked why the cross coupling is required for stability. Clearly we have more complex interaction. The constant $\omega$ surfaces will also be more complex. However, I still don’t have a good intuition what exactly was sought after for that part of the question.

Numerical computations

A Manipulate allowing for choice of the spring constants and lattice orientation, as shown in fig. 1.2, is available in phy487/oneAtomBasisPhonon.nb. This interface also provides a numerical calculation of the distribution relation as shown in fig. 1.3, and provides an animation of the normal modes for any given selection of $\mathbf{q}$ and $\omega(\mathbf{q})$ (not shown).

Fig 1.2: 2D Single atom basis Manipulate interface

Fig 1.3: Sample distribution relation for 2D single atom basis.

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

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

# Disclaimer

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

or

\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)

or

\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)

## PHY452H1S Basic Statistical Mechanics. Lecture 20: Bosons. Taught by Prof. Arun Paramekanti

Posted by peeterjoot on April 2, 2013

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

# Disclaimer

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

# Bosons

In order to maintain a conservation of particles in a Bose condensate as we decrease temperature, we are forced to change the chemical potential to compensate. This is illustrated in fig. 1.1.

Fig 1.1: Chemical potential in Bose condensation region

Bose condensatation occurs for $T < T_{\mathrm{BEC}}$. At this point our number density becomes (except at $\mathbf{k} = 0$)

\begin{aligned}n(\mathbf{k}) = \frac{1}{{e^{\beta \epsilon_\mathbf{k}} - 1}}.\end{aligned} \hspace{\stretch{1}}(1.2.1)

Except for $\mathbf{k} = 0$, $n(\mathbf{k})$ is well defined, and not described by this distribution. We are forced to say that

\begin{aligned}N = N_0 + \sum_{\mathbf{k} \ne 0} n(\mathbf{k}) = N_0 + V\int \frac{d^3 \mathbf{k}}{(2 \pi)^3} \frac{1}{{ e^{\beta \epsilon_\mathbf{k}} - 1 }}.\end{aligned} \hspace{\stretch{1}}(1.2.1)

Introducing the density of states, our density is

\begin{aligned}\rho = \rho_0 + \int_0^\infty d\epsilon \frac{N(\epsilon)}{e^{\beta \epsilon} - 1 },\end{aligned} \hspace{\stretch{1}}(1.2.3)

where

\begin{aligned}N(\epsilon) = \frac{1}{{4 \pi^2}} \left( \frac{2m}{\hbar} \right)^{3/2} \epsilon^{1/2}.\end{aligned} \hspace{\stretch{1}}(1.2.4)

We worked out last time that

\begin{aligned}\rho = \rho_0 + \rho \left( \frac{T}{T_{\mathrm{BEC}}} \right)^{3/2},\end{aligned} \hspace{\stretch{1}}(1.2.4)

or

\begin{aligned}\rho_0 = \rho \left( 1 - \left( \frac{T}{T_{\mathrm{BEC}}} \right) ^{3/2} \right).\end{aligned} \hspace{\stretch{1}}(1.2.6)

This is plotted in fig. 1.2.

Fig 1.2: Density variation with temperature for Bosons

\begin{aligned}\rho_0 = \frac{N_{\mathbf{k} = 0}}{V}.\end{aligned} \hspace{\stretch{1}}(1.7)

For $T \ge T_{\mathrm{BEC}}$, we have $\rho_0 = 0$. This condensation temperature is

\begin{aligned}T_{\mathrm{BEC}} \propto \rho^{2/3}.\end{aligned} \hspace{\stretch{1}}(1.8)

This is plotted in fig. 1.3.

Fig 1.3: Temperature vs pressure demarkation by T_BEC curve

There is a line for each density that marks the boundary temperature for which we have or do not have this condensation phenomina where $\mathbf{k} = 0$ states start filling up.

Specific heat: $T < T_{\mathrm{BEC}}$

\begin{aligned}\frac{E}{V} &= \int \frac{d^3 \mathbf{k}}{(2 \pi)^3} \frac{1}{{ e^{\beta \hbar^2 k^2/2m} - 1}}\frac{\hbar^2 k^2}{2m} \\ &= \int_0^\infty d\epsilon N(\epsilon) \frac{1}{{ e^{\beta \epsilon} - 1 }} \epsilon \\ &\propto \int_0^\infty d\epsilon \frac{\epsilon^{3/2}}{ e^{\beta \epsilon} - 1 } \\ &\propto \left( k_{\mathrm{B}} T \right)^{5/2},\end{aligned} \hspace{\stretch{1}}(1.9)

so that

\begin{aligned}\frac{C}{V} \propto \left( k_{\mathrm{B}} T \right)^{3/2}.\end{aligned} \hspace{\stretch{1}}(1.10)

Compare this to the classical and Fermionic specific heat as plotted in fig. 1.4.

Fig 1.4: Specific heat for Bosons, Fermions, and classical ideal gases

One can measure the specific heat in this Bose condensation phenomina for materials such as Helium-4 (spin 0). However, it turns out that Helium-4 is actually quite far from an ideal Bose gas.

Photon gas

A system that is much closer to an ideal Bose gas is that of a gas of photons. To a large extent, photons do not interact with each other. This allows us to calculate black body phenomina and the low temperature (cosmic) background radiation in the universe.

An important distinction between a photon sea and some of these other systems is that the photon number is actually not fixed.

Photon numbers are not “conserved”.

If a photon interacts with an atom, it can impart energy and disappear. An excited atom can emit a photon and change its energy level. In a thermodynamic system we can generally expect that introducing heat will generate more photons, whereas a cold sink will tend to generate fewer photons.

We have a few special details that distinguish photons that we’ll have to consider.

1. spin 1.
2. massless, moving at the speed of light.
3. have two polarization states.

Because we do not have a constraint on the number of particles, we essentially have no chemical potential, even in the grand canonical scheme.

Writing

\begin{aligned}\lambda = \left\{\begin{array}{l l}+1 & \quad \mbox{Right circular polarization} \\ -1 & \quad \mbox{Left circular polarization}\end{array}\right.\end{aligned} \hspace{\stretch{1}}(1.11)

Our number density, since we have no chemical potential, is of the form

\begin{aligned}n_{\mathbf{k}, \lambda}= \frac{1}{{e^{\beta \epsilon_{\mathbf{k}, \lambda}} - 1 }},\end{aligned} \hspace{\stretch{1}}(1.12)

Observe that the average number of photons in this system is temperature dependent. Because this chemical potential is not there, it can be quite easy to work out a number of the thermodynamic results.

Photon average energy density

We’ll now calculate the average energy density of the photons. The energy of a single photon is

\begin{aligned}\epsilon_{\mathbf{k}, \lambda} = \hbar c k = \hbar \omega,\end{aligned} \hspace{\stretch{1}}(1.2.13)

so that the average energy density is

\begin{aligned}\frac{E}{V} &= \sum_{\mathbf{k}, \lambda} \frac{1}{{ e^{ \beta \epsilon_\mathbf{k}} - 1}} \epsilon_\mathbf{k}\rightarrow\underbrace{2}_{\text{number of polarizations}}\int \frac{d^3 \mathbf{k}}{(2 \pi)^3}\frac{ \hbar c k}{ e^{ \beta \epsilon_\mathbf{k}} - 1} \\ &= 2 \int_0^\infty d\epsilon \underbrace{\frac{1}{{(2 \pi)^3}} 4 \pi \frac{\epsilon^2}{(\hbar c)^3} }_{\text{Photon density of states}}\frac{\epsilon}{e^{\beta \epsilon} - 1} \\ &= \frac{1}{{\pi^2}} \frac{1}{{ (\hbar c)^3 }} \int_0^\infty d\epsilon \frac{\epsilon^3}{e^{\beta \epsilon} - 1}\end{aligned} \hspace{\stretch{1}}(1.2.13)

Mathematica tells us that this integral is

\begin{aligned}\int_0^\infty d\epsilon \frac{\epsilon^3}{e^{\beta \epsilon} - 1} =\frac{\pi ^4}{15 \beta ^4},\end{aligned} \hspace{\stretch{1}}(1.2.13)

for an end result of

\begin{aligned}\frac{E}{V} =\frac{\pi^2}{15} \frac{1}{{(\hbar c)^3}} \left( k_{\mathrm{B}} T \right)^4.\end{aligned} \hspace{\stretch{1}}(1.2.13)

Phonons and other systems

There is a very similar phenomina in matter. We can discuss lattice vibrations in a solid. These are called phonon modes, and will have the same distribution function where the only difference is that the speed of light is replaced by the speed of the sound wave in the solid. Once we understand the photon system, we are able to look at other Bose distributions such as these phonon systems. We’ll touch on this very briefly next time.