    from matplotlib import pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

import numpy as np
from scipy.optimize import curve_fit
from scipy.integrate import quad

from common import draw_classic_axes, configure_plotting

configure_plotting()

Solutions for lecture 2 exercises¶

Warm-up exercises¶

For low $T$, $T_D/T \rightarrow \infty$. The heat capacity is then given by:

\[ C \overset{\mathrm{low \: T}}{\approx} 9Nk_{\mathrm{B}}\left(\frac{T}{T_{D}}\right)^3\int_0^{\infty}\frac{x^4{\mathrm{e}}^x}{({\mathrm{e}}^x-1)^2}{\mathrm{d}}x. \]
See the plot below (shown for $T_{D,1} < T_{D,2}$).
The polarization describes the direction of motion of the atoms in the wave with respect to the direction in which the wave travels. In 3D, there are three possible polarizations.
The integral can be expressed as $$ \int \mathrm{d}k_x \, \mathrm{d}k_y = \int_0^{2\pi} \mathrm{d}\phi \int_0^\infty k \, \mathrm{d}k. $$
The Debye frequency $\omega_D$ is the frequency of the vibrational mode with the highest eigenfrequency. It has a corresponding Debye temperature $T_D = \hbar\omega_D/k_B$, which sets the temperature scale above which essentially all vibrational modes in the system are thermally populated.
The wavelength is of the order of the interatomic spacing:

\[ \lambda = (\frac{4}{3}\pi)^{1/3} a. \]

fig, ax = plt.subplots()
T = np.linspace(0.1, 3)
T_D = [1,2]
ax.plot(T, (T/T_D[0])**3, 'b-', label = r'$T_{D,1}$')
ax.plot(T, (T/T_D[1])**3, 'r-', label = r'$T_{D,2}$')
ax.set_ylim([0,3])
ax.set_xlim([0,3])
ax.set_xlabel('$T$')
ax.set_xticks([0])
ax.set_xticklabels(['$0$'])
ax.set_ylabel('$C$')
ax.set_yticks([0])
ax.set_yticklabels(['$0$'])
ax.legend();

svg

Exercise 1: Deriving the density of states for the linear dispersion relation of the Debye model¶

$\omega = v_s|\mathbf{k}|$
The distance between nearest-neighbour points in $\mathbf{k}$-space is $2\pi/L$. The density of $\mathbf{k}$-points in 1, 2, and 3 dimensions is $L/(2\pi)$, $L^2/(2\pi)^2$, and $L^3/(2\pi)^3$ respectively.
Express the number of states between frequencies $0<\omega<\omega_0$ as an integral over k-space. Do so for 1D, 2D and 3D. Do not forget the possible polarizations. We assume that in $d$ dimensions there are $d_p$ polarizations.

\[\begin{align*} N_\text{states, 1D} & = \frac{L}{2\pi} \int_0^{k_0} 2 \, \mathrm{d}k \\ N_\text{states, 2D} & = 2 \left(\frac{L}{2\pi}\right)^2 \int_0^{k_0} 2\pi k \, \mathrm{d}k \\ N_\text{states, 3D} & = 3 \left(\frac{L}{2\pi}\right)^3 \int_0^{k_0} 4\pi k^2 \, \mathrm{d}k \end{align*}\]
We use $k = |\mathbf{k}| = \omega/v_s$ and $\mathrm{d}k = \mathrm{d}\omega/v_s$ to get

\[\begin{align*} N_\text{states, 1D} & = \frac{L}{2\pi} \int_0^{\omega_0} 2 \frac{1}{v_s} \, \mathrm{d}\omega := \int_0^{\omega_0} g_{1D}(\omega) \, \mathrm{d}\omega \\ N_\text{states, 2D} & = 2 \left(\frac{L}{2\pi}\right)^2 \int_0^{\omega_0} 2\pi \frac{\omega}{v_s^2} \, \mathrm{d}\omega := \int_0^{\omega_0} g_{2D}(\omega) \, \mathrm{d}\omega \\ N_\text{states, 3D} & = 3 \left(\frac{L}{2\pi}\right)^3 \int_0^{\omega_0} 4\pi \frac{\omega^2}{v_s^3} \, \mathrm{d}\omega := \int_0^{\omega_0} g_{3D}(\omega) \, \mathrm{d}\omega \end{align*}\]

The integral boundaries set the frequency region in which you calculate the density of states.
The density of states is the number of states per unit frequency. It therefore has units of inverse frequency.

Exercise 2: Debye model in 2D¶

The energy stored in the vibrational modes of a two-dimensional Debye solid is:

\[\begin{align*} E & = 2 \left(\frac{L}{2\pi}\right)^2 \int \left(n_B(\omega(\mathbf{k}))+\frac{1}{2}\right)\hbar\omega(\mathbf{k}) \, \mathrm{d}^2\mathbf{k} \\ & = \frac{L^2}{\pi v_s^2\hbar^2\beta^3}\int_{0}^{\beta\hbar\omega_D}\frac{x^2}{e^{x} - 1}\,\mathrm{d}x + E_0 \end{align*}\]
The high-$T$ limit implies $\beta \rightarrow 0$. Therefore, $n_B \approx k_B T/\hbar\omega$, and the integral becomes particularly illuminating:

\[ E = \int_{0}^{\omega_D} \hbar\omega n_B(\omega) g(\omega) d\omega \approx \int_{0}^{\omega_D} k_B T g(\omega) d\omega = N_\text{modes} k_B T \]

where we neglected the zero-point energy. In 2D, we have $N_\text{modes} = 2N_\text{atoms}$, so we recover the 2D Dulong-Petit law $C_v = \mathrm{d}E/\mathrm{d}T = 2 k_B$ per atom.
In the low temperature limit, the high-energy modes are not excited so we can safely let the upper boundary of the integral go to infinity. For convenience, we write $g(\omega) = \alpha \omega$, with $\alpha = \frac{L^2}{\pi v_s^2}$. We get

\[ E = \int_0^{\omega_D}\frac{\hbar\omega g(\omega)}{e^{\hbar\omega/k_B T}- 1}\,\mathrm{d}\omega \approx \alpha\frac{k_B^3 T^3}{\hbar^2}\int_0^\infty\frac{x^2}{e^x-1}\,\mathrm{d}x \]

From this we find $C_v = \mathrm{d}E/\mathrm{d}T = K T^2$, with

\[ K = 3\alpha\frac{k_B^3}{\hbar^2}\int_0^\infty\frac{x^2}{e^x-1}\,\mathrm{d}x \]

Exercise 3: Longitudinal and transverse vibrations with different sound velocities¶

The key idea is that the total energy in the individual harmonic oscillators (the vibrational modes) is the sum of the energies in the individual oscillators. There is, however, a small ambiguity in how to impose the Debye cutoff once the two branches have different sound velocities. The exercise does not specify which convention to use, so both of the following are acceptable.

In both cases,

\[\begin{align*} g_\parallel(\omega) & = \frac{L^3}{2\pi^2} \frac{\omega^2}{v_\parallel^3}, \\ g_\perp(\omega) & = 2 \frac{L^3}{2\pi^2} \frac{\omega^2}{v_\perp^3}, \end{align*}\]

where $g_\perp$ already includes the two transverse polarizations.

The simplest interpretation is to keep a single Debye frequency for the combined density of states, exactly as in the lecture:

\[ E = \int_0^{\omega_D}\frac{\hbar \omega \left[g_\parallel(\omega) + g_\perp(\omega)\right]}{e^{\beta\hbar\omega} - 1}d\omega + E_Z. \]

Here $\omega_D$ is chosen such that $\int_0^{\omega_D}\left[g_\parallel(\omega) + g_\perp(\omega)\right]d\omega = 3N$. This gives

\[ E = \frac{L^3 k_B^4 T^4}{2\pi^2\hbar^3}\left(\frac{2}{v_\perp^3} + \frac{1}{v_\parallel^3}\right)\int_{0}^{\beta\hbar\omega_D}\frac{x^3}{e^{x} - 1}dx. \]

A somewhat more physical interpretation is to use a common cutoff in $\mathbf{k}$-space, since the normal modes are really counted by allowed $\mathbf{k}$-values. Then the different branches have different Debye frequencies,

\[ \omega_{D,\parallel} = v_\parallel k_D, \qquad \omega_{D,\perp} = v_\perp k_D, \]

with $k_D$ chosen so that the total number of modes is $3N$. The energy is then written as

\[ E = \int_0^{\omega_{D,\parallel}}\frac{\hbar \omega g_\parallel(\omega)}{e^{\beta\hbar\omega} - 1}d\omega + \int_0^{\omega_{D,\perp}}\frac{\hbar \omega g_\perp(\omega)}{e^{\beta\hbar\omega} - 1}d\omega + E_Z, \]

or equivalently

\[ E = \frac{L^3 k_B^4 T^4}{2\pi^2\hbar^3}\left[ \frac{1}{v_\parallel^3}\int_{0}^{\beta\hbar\omega_{D,\parallel}}\frac{x^3}{e^{x} - 1}dx + \frac{2}{v_\perp^3}\int_{0}^{\beta\hbar\omega_{D,\perp}}\frac{x^3}{e^{x} - 1}dx \right]. \]

This second choice is somewhat more physical, but for the limiting behaviors asked in the exercise both conventions lead to the same final answers.
In the high-$T$ limit, we have $\beta \rightarrow 0$. Therefore, $n_B \approx k_B T/\hbar\omega$, and in either convention the integral for $E$ becomes:

\[ E \approx k_B T \, N_\text{modes} + E_Z \]

and we are left with the Dulong-Petit law $C = 3N_\text{atoms} k_B$.
In the low temperature limit, we can let the upper integral boundary go to infinity as in exercise 1. This is valid for either a single $\omega_D$ or separate $\omega_{D,\parallel}$ and $\omega_{D,\perp}$, provided the temperature is well below the relevant Debye scale(s). We then get

\[ C \approx \frac{2\pi^2 k_B^4 L^3}{15\hbar^3}\left(\frac{2}{v_\perp^3} + \frac{1}{v_\parallel^3}\right)T^3 \]

where we used $\int_{0}^{\infty}\frac{x^3}{e^{x} - 1}dx = \frac{\pi^4}{15}$. So the ambiguity in the cutoff prescription only affects the intermediate-temperature interpolation, not the high-$T$ and low-$T$ limits requested here.

Exercise 4: Anisotropic sound velocities¶

In this case, the velocity depends on the direction. Note however that, in contrast with the previous exercise, the polarization does not affect the dispersion of the waves. We get

\[\begin{align*} E =& 3\left(\frac{L}{2\pi}\right)^3 \int \frac{\hbar\omega(\mathbf{k})}{e^{\beta\hbar\omega(\mathbf{k})}-1} \, \mathrm{d}k_x \, \mathrm{d}k_y \, \mathrm{d}k_z + E_Z \\ =& 3\left(\frac{L}{2\pi}\right)^3\frac{1}{v_x v_y v_z} \int \frac{\hbar\kappa}{e^{\beta\hbar\kappa} - 1} \, \mathrm{d}\kappa_x \, \mathrm{d}\kappa_y \, \mathrm{d}\kappa_z + E_Z, \end{align*}\]

where we made the substitutions $\kappa_x = k_x v_x,\kappa_y = k_y v_y, \kappa_z = k_z v_z$ so that $\omega = \kappa$. Going to spherical coordinates:

\[\begin{align*} E &= \frac{3 L^3}{2\pi^2}\frac{1}{v_x v_y v_z}\int_{0}^{\kappa_D} \frac{\hbar\kappa^3}{e^{\beta\hbar\kappa} - 1} \, \mathrm{d}\kappa + E_Z \end{align*}\]

To calculate the specific heat $C = \frac{dE}{dT}$, let's differentiate this expression directly:

\[\begin{align*} C &= \frac{\mathrm{d}E}{\mathrm{d}T} = \frac{3 L^3}{2\pi^2}\frac{1}{v_x v_y v_z}\int_{0}^{\kappa_D} \frac{\mathrm{d}}{\mathrm{d}T}\left(\frac{\hbar\kappa^3}{e^{\beta\hbar\kappa} - 1}\right) \, \mathrm{d}\kappa \\ &= \frac{3 L^3}{2\pi^2}\frac{1}{v_x v_y v_z}\int_{0}^{\kappa_D} \hbar\kappa^3 \left(\frac{\partial}{\partial\beta}\frac{1}{e^{\beta\hbar\kappa} - 1}\right) \frac{\mathrm{d}\beta}{\mathrm{d}T} \, \mathrm{d}\kappa \\ &= \frac{3 L^3}{2\pi^2}\frac{1}{v_x v_y v_z}\int_{0}^{\kappa_D} \hbar\kappa^3 \left(-\frac{\hbar\kappa e^{\beta\hbar\kappa}}{(e^{\beta\hbar\kappa}-1)^2}\right) \left(-\frac{\beta}{T}\right) \, \mathrm{d}\kappa \\ &= \frac{3 L^3}{2\pi^2}\frac{1}{v_x v_y v_z} \cdot \frac{\hbar^2\beta}{T}\int_{0}^{\kappa_D} \frac{\kappa^4 e^{\beta\hbar\kappa}}{(e^{\beta\hbar\kappa} - 1)^2} \, \mathrm{d}\kappa. \end{align*}\]

Making the substitution $x = \beta\hbar\kappa$ (so $d\kappa = dx/(\beta\hbar)$) yields

\[\begin{align*} C &= \frac{3 L^3}{2\pi^2}\frac{1}{v_x v_y v_z} \cdot \frac{\hbar^2\beta}{T} \cdot \frac{1}{(\beta\hbar)^5} \int_{0}^{\beta\hbar\kappa_D} \frac{x^4 e^x}{(e^x - 1)^2} \, \mathrm{d}x \\ &= \frac{3 L^3}{2\pi^2}\frac{1}{v_x v_y v_z} \cdot \frac{1}{T\beta^4\hbar^3} \int_{0}^{\beta\hbar\kappa_D} \frac{x^4 e^x}{(e^x - 1)^2} \, \mathrm{d}x. \end{align*}\]

Using $\beta = 1/(k_B T)$ so that $1/(T\beta^4) = k_B^4 T^3$, we obtain the final result

\[\begin{align*} C &= \frac{3 L^3 k_B^4 T^3}{2\pi^2\hbar^3}\frac{1}{v_x v_y v_z} \int_{0}^{\beta\hbar\kappa_D} \frac{x^4 e^x}{(e^x - 1)^2} \, \mathrm{d}x. \end{align*}\]

This is the same form as the isotropic result but with the replacement $1/v^3 \to 1/(v_x v_y v_z)$.