```python
from matplotlib import pyplot
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()
```
_(based on chapter 2.2 of the book)_
!!! summary "Learning goals"
After this lecture you will be able to:
- Describe the concept of reciprocal space and allowed momenta
- Write down the total energy of phonons given the temperature and the dispersion relation
- Estimate heat capacity due to phonons in high temperature and low temperature regimes of the Debye model
### Debye model
The Einstein model explained the experimental data quite well, but still slightly underestimated the observed values of $C$. Apparently the "each atom is an oscillator"-idea is too simplistic.
Peter Debye (1884 – 1966) suggested to instead consider _normal modes_: sound waves that propagate through a solid with velocity $v=\omega/k$, where $k=2\pi/\lambda$ is the _wave number_. Instead of just multiplying $\bar{\varepsilon}$ by $3N$, we now use:
$$E=\int\limits_0^\infty\left(\frac{1}{2}\hbar\omega+\frac{\hbar\omega}{ {\rm e}^{\hbar\omega/{k_{\rm B}T}}-1}\right)g(\omega){\rm d}\omega$$
$g(\omega)$ is the _density of states_: the number of normal modes found at each position along the $\omega$-axis. How do we calculate $g(\omega)$?
#### Reciprocal space, periodic boundary conditions
Each normal mode can be described by a _wave vector_ ${\bf k}$. A wave vector represents a point in _reciprocal space_ or _k-space_. We can find $g(\omega)$ by counting the number of normal modes in k-space and then converting those to $\omega$.
How far apart are the normal modes in k-space? This is determined by the boundaries of the solid. One way to treat the boundaries is by using _fixed boundary conditions_ (like a guitar string), resulting in modes $k=\pi/L$, $2\pi/L$, $3\pi/L$ etc., where $L$ is the length of the solid.
In this course, however, we will exclusively use _periodic boundary conditions_, where one edge of the solid should connect seamlessly to the opposite edge. This results in:
$$k=\ ..., \frac{-4\pi}{L}, \frac{-2\pi}{L}, 0, \frac{2\pi}{L}, \frac{4\pi}{L}, ...$$
The three-dimensional wave vector ${\bf k}$ can be any threefold permutation of these values. All possible values of ${\bf k}$ then form a grid in k-space:
There is one allowed ${\bf k}$ per $\left(\frac{2\pi}{L}\right)^3$. The number of ${\bf k}$-values inside a sphere with radius $k$:
$$N=\frac{\frac{4}{3}\pi k^3}{\left(\frac{2\pi}{L}\right)^3}=\frac{Vk^3}{6\pi^2}$$
This means for the density of states $g(k)$ as a function of $k$:
$$g(k)=\frac{ {\rm d}N}{ {\rm d}k}=\frac{Vk^2}{2\pi^2}$$
The density of states $g(\omega)$ in frequency space then becomes:
$$g(\omega)=g(k)\frac{ {\rm d}k}{ {\rm d}\omega}=\frac{Vk^2}{2\pi^2}\frac{ {\rm d}k}{ {\rm d}\omega}$$
In general, ${\rm d}k/{\rm d}\omega$ can be difficult to calculate; we will see more of this later. But going back to the Debye model for now, where we assume simple sound waves with $v=\omega/k$, this reduces to $g(\omega)=V\omega^2/2\pi^2v^3$. The total energy then becomes:
$$ E=E_{\rm Z}+\frac{3V}{2\pi^2 v_{\rm s}^3}\int\limits_0^\infty\left(\frac{\hbar\omega}{ {\rm e}^{\hbar\omega/k_{\rm B}T}-1}\right)\omega^2{\rm d}\omega$$
Here, the factor 3 comes from the fact that every wave has three polarizations (two transversal, one longitudinal). The term $E_{\rm Z}$ goes to infinity through integration. This is no problem, as it doesn't count towards the heat capacity.
Substitute $x\equiv\frac{\hbar\omega}{k_{\rm B}T}$:
$$\Rightarrow E=E_{\rm Z}+\frac{3V}{2\pi^2 v_{\rm s}^3}\frac{\left(k_{\rm B}T\right)^4}{\hbar^3}\int\limits_0^\infty\frac{x^3}{ {\rm e}^x-1}{\rm d}x$$
The integral on the right is a constant, $\left(\frac{\pi^4}{15}\right)$ $\Rightarrow$ $C=\frac{ {\rm d}E}{ {\rm d}T}\propto T^3$.
#### Debye's interpolation for medium T
The above approximation works very well at low temperature. But at high temperature, $C$ should of course settle at $3k_{\rm B}$ (the Dulong-Petit value). The reason why the model breaks down, is that it assumes that there is an infinite number of harmonic oscillators up to infinite frequency.
Debye proposed an approximation: all phonons are acoustic (i.e. constant sound velocity) until a certain cut-off frequency, beyond which there are no phonons.
$$
g(\omega) = \left\{
\begin{array}{ll}
\frac{3V\omega^2}{2\pi^2v_{\rm s}^3} & \omega<\omega_{\rm D} \\
0 & \omega>\omega_{\rm D}
\end{array}
\right.
$$
What determines the _Debye frequency_ $\omega_{\rm D}$?
$$
\int_0^{\omega_{\rm D}}g(\omega){\rm d}\omega=\frac{V\omega_{\rm D}^3}{2\pi^2v_{\rm s}^3}=3N,
$$
where $N$ is the number of atom, so $3N$ the number of degrees of freedom.
Substitute in $E$, differentiate to $T$:
$$
\Rightarrow C=9Nk_{\rm B}\left(\frac{T}{\Theta_{\rm D}}\right)^3\int_0^{\Theta_{\rm D}/T}\frac{x^4{\rm e}^x}{({\rm e}^x-1)^2}{\rm d}x,
$$
where $x=\frac{\hbar\omega}{k_{\rm B}T}$ and $\Theta_{\rm D}\equiv\frac{\hbar\omega_{\rm D}}{k_{\rm B}}$, the _Debye temperature_.
```python
pyplot.rcParams['axes.titlepad'] = 20
T = np.array([1.35,2.,3.,4.,5.,6.,7.,8.,10.,12.,14.,16.,20.,28.56,36.16,47.09,55.88,65.19,74.56,83.91,103.14,124.2,144.38,166.78,190.17,205.3])
c = np.array([0.,0.,0.,0.,0.,0.,0.0719648,0.1075288,0.2100368,0.364008,0.573208,0.866088,1.648496,4.242576,7.07096,10.8784,13.47248,15.60632,17.27992,18.6188,20.33424,21.63128,22.46808,23.05384,23.47224,23.68144])
c *= 3/24.945 #24.954 is 3Nk_B
def c_einstein(T, T_E):
x = T_E / T
return 3 * x**2 * np.exp(x) / (np.exp(x) - 1)**2
def integrand(y):
return y**4 * np.exp(y) / (np.exp(y) - 1)**2
@np.vectorize
def c_debye(T, T_D):
x = T / T_D
return 9 * x**3 * quad(integrand, 0, 1/x)[0]
temp = np.linspace(1, 215, 100)
fit = curve_fit(c_einstein, T, c, 500)
T_E = fit[0][0]
#delta_T_E = np.sqrt(fit[1][0, 0])
#print(f"T_E = {T_E:.5} ± {delta_T_E:.3} K")
fit = curve_fit(c_debye, T, c, 500)
T_D = fit[0][0]
#delta_T_D = np.sqrt(fit[1][0, 0])
#print(f"T_D = {T_D:.5} ± {delta_T_D:.3} K")
fig, ax = pyplot.subplots()
ax.scatter(T, c)
ax.set_title('Heat capacity of silver compared to the Debye and Einstein models')
ax.plot(temp, c_einstein(temp, T_E), label='Einstein model')
ax.plot(temp, c_debye(temp, T_D), label='Debye model')
ax.set_ylim(bottom=0, top=3.4)
ax.set_xlim(0, 215)
ax.set_xlabel('$T(K)$')
ax.set_ylabel(r'$C/k_B$')
ax.set_xticks([T_E, T_D])
ax.set_xticklabels(['$T_E$','$T_D$'])
ax.set_yticks([3])
ax.set_yticklabels(['$3$'])
ax.legend(loc='upper left')
pyplot.hlines([3], 0, 215, linestyles='dashed')
pyplot.vlines([T_E,T_D], 0, 3.4, linestyles='dashed')
#draw_classic_axes(ax, xlabeloffset=0.3)
```
## Exercises
### Exercise 1: Debye model: concepts
1. Describe the concept of k-space. What momenta are allowed in a 2D system with dimensions $L\times L$?
2. The probability to find an atom of a 1D solid that originally had a position $x$ at a displacement $\delta x$ is shown on this plot:
```python
def psi_squared(delta_x, x):
return delta_x**2 * np.exp(-delta_x**2) * np.sin(4*np.pi*x)**2
x = np.linspace(0, 1, 200)
delta_x = np.linspace(-2, 2, 200)
pyplot.imshow(psi_squared(delta_x.reshape((-1, 1)), x.reshape((1, -1))), cmap='gist_heat_r', extent=(0, 3, -1, 1))
pyplot.ylabel(r'$\delta x$')
pyplot.xlabel(r'$x$')
pyplot.xticks((0, 3), ('$0$', '$L$'))
pyplot.yticks((), ())
cbar = pyplot.colorbar()
cbar.set_ticks(())
cbar.set_label(r'$|\psi^2|$')
```
Describe how many phonons in which $k$-state this solid has.
Explain your answer.
??? hint
There are $n=2$ phonons in the state with $k=4\pi/L$ and $n=2$ phonons in a state with $k=-4\pi/L$.
3. Explain the concept of density of states.
4. Calculate the density of states $g(\omega)$ for a 3D, 2D and 1D systems with linear dispersion $\omega=vk$.
### Exercise 2: Debye model in 2D
1. State the assumptions of the Debye theory.
2. Determine the energy of a two-dimensional solid as a function of $T$ using the Debye approximation (the integral can't be solved analytically).
3. Calculate the heat capacity in the limit of high $T$ (hint: it goes to a constant).
4. At low $T$, show that $C_V=KT^{n}$. Find $n$. Express $K$ as a definite integral.
### Exercise 3: Different phonon modes
During the lecture we derived the low-temperature heat capacity assuming that all the phonons have the same sound velocity $v$.
In reality the longitudinal and transverse modes have different sound velocities (see [Wikipedia](https://en.wikipedia.org/wiki/Sound#Longitudinal_and_transverse_waves) for an illustration of different sound wave types).
Assume that there are two types of excitations:
* One longitudinal mode with $\omega = v_\parallel |k|$
* Two transverse modes with $\omega = v_\bot |k|$
1. Write down the total energy of phonons in this material *(hint: use the same reasoning as in the [Lithium exercise](1_einstein_model.md#exercise-4-total-heat-capacity-of-a-diatomic-material))*.
2. Verify that at high $T$ you reproduce the Dulong-Petit law.
3. Compute the behavior of heat capacity at low $T$.
### Exarcise 4: Anisotropic sound velocities
Suppose now that the velocity is anisotropic ($v_x \neq v_y \neq v_z$) and $\omega = \sqrt{v_x^2 k_x^2 + v_y^2 k_y^2 + v_z^2 k_z^2}$.
How does this change the Debye result for the heat capacity?
??? hint
Write down the total energy as an integral over $k$, then change the integration variables so that the spherical symmetry of the integrand is restored.