2_debye_model.md

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.

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]

fit = curve_fit(c_debye, T, c, 500)
T_D = fit[0][0]

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_debye(temp, T_D), label=f'Debye model, $T_D={T_D:.5}K$')
ax.plot(temp, c_einstein(temp, T_E), label=f'Einstein model, $T_E={T_E:.5}K$')
ax.set_ylim(bottom=0, top=3)
ax.set_xlim(0, 215)
ax.set_xlabel('$T(K)$')
ax.set_ylabel(r'$C/k_B$')
ax.legend(loc='lower right');

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:

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