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

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

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

2. See plot below (shown for $T_{D,1} < T_{D,2}$)

3. The polarization describes the direction of the motion of the atoms in the wave with respect to the direction in which the wave travels. In 3D, there are only 3 different polarizations possible.

4. The integral can be expressed as $$\int dk_x dk_y = \int k dk d\phi$$

5. The Debye frequency $\omega_D$ is the frequency of the vibrational mode with the highest eigenfrequency. It has corresponding Debye temperature $T_D = \hbar\omega_D/k_B$, which is the temperature above which all the vibrational modes in the system become excited

6. 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();

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

1. $\omega = v_s|\mathbf{k}|$
2. 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.

3. 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} & = 1_p \frac{L}{2\pi} \int_0^{k_0} 2 dk \\N_\text{states, 2D} & = 2_p \left(\frac{L}{2\pi}\right)^2 \int_0^{k_0} 2\pi k dk \\N_\text{states, 3D} & = 3_p \left(\frac{L}{2\pi}\right)^3 \int_0^{k_0} 4\pi k^2 dk\end{align*}$$

4. We use $k=\mathbf{k}| = \omega/v_s$ and $dk = d\omega/v_s$ to get

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

   The integral boundaries set the frequency region in which you calculate the density of states.

5. The density of states is the number of states per unit frequency. It has units of 1 over frequency

Exercise 2: Debye model in 2D

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

   $$\begin{align*}E & = 2_p \frac{L^2}{4\pi^2}\int_0^{\infty}(n_B(\omega(\mathbf{k}))+\frac{1}{2})\hbar\omega(\mathbf{k}) d\mathbf{k} \\& = \frac{L^2}{\pi v^2\hbar^2\beta^3}\int_{0}^{\beta\hbar\omega_D}\frac{x^2}{e^{x} - 1}dx + E_0\end{align*}$$

2. 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} = 2_p N_\text{atoms}$, so that we recover the 2D law of Dulong–Petit $C_v = dE/dT = 2 k_B$ per atom.

3. 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}d\omega \approx \alpha\frac{k^3 T^3}{\hbar^2}\int_0^\infty\frac{x^2}{e^x-1}dx$$

   From which we find $C_v = dE/dT = K T^2$, with

   $$K = 3\alpha\frac{k^3}{\hbar^2}\int_0^\infty\frac{x^2}{e^x-1}dx$$

Exercise 3: Longitudinal and transverse vibrations with different sound velocities

1. 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: $E = \int_0^{\omega_D}\frac{\hbar \omega g(\omega)}{e^{\beta\hbar\omega} - 1}d\omega + E_Z$, where $g(\omega) = g_\parallel(\omega) + g_\perp(\omega)$. Using

   $$\begin{align*}g_\parallel & = 1_p \frac{L^3}{2\pi^2} \frac{\omega^2}{v_\parallel^3}, \\g_\perp & = 2_p \frac{L^3}{2\pi^2} \frac{\omega^2}{v_\perp^3},\end{align*}$$

   we get

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

2. As in exercise 1, in the high-T limit, we have $\beta \rightarrow 0. Therefore, $n_B \approx k_B T/\hbar\omega$ and the integral for $E$ becomes:

   $$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 +E_Z = N_\text{modes} k_B T + E_Z$$

   and we are left with the Dulong-Petit law $C = 3N_\text{atoms} k_B$.

3. In the low temperature limit, we can let the upper integral boundary go to infinity as in exercise 1. This yields

   $$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}$.

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_p\left(\frac{L}{2\pi}\right)^3 \int \frac{\hbar\omega(\mathbf{k})}{e^{\beta\hbar\omega(\mathbf{k})}-1} dk_x dk_y dk_z + E_Z \\=& 3_p\left(\frac{L}{2\pi}\right)^3\frac{1}{v_x v_y v_z} \int \frac{\hbar\kappa}{e^{\beta\hbar\kappa} - 1}d\kappa_x d\kappa_y 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_p 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} 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{dE}{dT} = \frac{3_p L^3}{2\pi^2}\frac{1}{v_x v_y v_z}\int_{0}^{\kappa_D} \frac{d}{dT}\left(\frac{\hbar\kappa^3}{e^{\beta\hbar\kappa} - 1}\right) d\kappa \\&= \frac{3_p L^3}{2\pi^2}\frac{1}{v_x v_y v_z}\int_{0}^{\kappa_D} \frac{\hbar\kappa^3 \cdot \hbar\kappa \cdot e^{\beta\hbar\kappa}}{(e^{\beta\hbar\kappa} - 1)^2} \cdot \frac{d\beta}{dT} d\kappa\end{align*}$$

Since $\beta = \frac{1}{k_B T}$, we have $\frac{d\beta}{dT} = -\frac{1}{k_B T^2} = -\frac{\beta}{T}$. Substituting this:

$$\begin{align*}C &= - \frac{3_p L^3}{2\pi^2}\frac{1}{v_x v_y v_z}\int_{0}^{\kappa_D} \frac{\hbar^2\kappa^4 e^{\beta\hbar\kappa}}{(e^{\beta\hbar\kappa} - 1)^2} \cdot \frac{\beta}{T} d\kappa \\&= - \frac{3_p 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} d\kappa\end{align*}$$

Now, we can make the substitution $x = \beta\hbar\kappa$, which gives $d\kappa = \frac{dx}{\beta\hbar}$ and changes our limits to $\int_{0}^{\beta\hbar\kappa_D}$:

$$\begin{align*}C &= - \frac{3_p L^3}{2\pi^2}\frac{1}{v_x v_y v_z} \cdot \frac{\hbar^2\beta}{T} \cdot \int_{0}^{\beta\hbar\kappa_D} \left(\frac{x}{\beta\hbar}\right)^4 \frac{e^x}{(e^x - 1)^2} \cdot \frac{dx}{\beta\hbar} \\&= - \frac{3_p 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} dx \\&= - \frac{3_p 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} dx\end{align*}$$

Since $\frac{1}{T\beta^4} = \frac{k_B^4 T^3}{1} = k_B^4 T^3$, we get:

$$\begin{align*}C &= -3_p \frac{3_p L^3}{2\pi^2}\frac{1}{v_x v_y v_z} \cdot \frac{k_B^4 T^3}{\hbar^3} \int_{0}^{\beta\hbar\kappa_D} \frac{x^4 e^x}{(e^x - 1)^2} dx \\&= 3_p \frac{3_p 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} dx\end{align*}$$

This is similar to the result with isotropic linear dispersion, with the difference being the factor $1/(v_x v_y v_z)$ instead of $1/v^3$.