Debye model

    from matplotlib import pyplot
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()
  

(based on chapter 2.2 of the book)

Expected prior knowledge

Before the start of this lecture, you should be able to:

• Recall how atoms are modeled in the Einstein model
• Derive the heat capacity of a solid within the Einstein model
• Describe how the frequency of a sound wave depends on the wavenumber $k=2\pi /\lambda$, with $\lambda$ the wavelength.
• Express a volume integral in spherical coordinates

Learning goals

After this lecture you will be able to:

• Describe the concept of reciprocal space and allowed momenta
• Describe the concept of a dispersion relation
• Derive the total number and energy of phonons in an object given the temperature and dispersion relation
• Estimate the heat capacity due to phonons in the high- and low-temperature regimes of the Debye model

Deficiency of the Einstein model

In the previous lecture, we observed that the Einstein model explained the heat capacity of solids quite well. However, we can see that something goes wrong if we compare the heat capacity predicted by the Einstein model to the that of silver¹:

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

temp = np.linspace(1, 215, 100)

fit = curve_fit(c_einstein, T, c, 500)
T_E = fit[0][0]

fig, ax = pyplot.subplots()
ax.scatter(T, c, label='Silver')
ax.plot(temp, c_einstein(temp, T_E), label=f'Einstein model, $T_E={T_E:.5}K$')
ax.text(T_E+3, 1.5, r'$T=T_E$', ha='left', color='r');
ax.plot([T_E, T_E], [0, 3], 'r--')
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.set_title(r'Heat capacity of silver and a fit of the Einstein model')
ax.legend();

The low-temperature heat capacity of silver is underestimated by the Einstein model. This underestimation is not obvious at first, but as we will see, this subtle difference is due to a profound physical phenomenon.

How does $C$ predicted by the Einstein model behave at low $T$?

When $T\to 0$, $T_E/T\to \mathrm{\infty }$. Therefore neglecting $1$ in the denominator we get $C\propto {\left(\frac{T_E}{T}\right)}^2{e}^{-T_E/T}$, and the heat capacity should be exponentially small.

Debye model

The key simplification of the Einstein model is to consider the atoms as independent quantum harmonic oscillators. Instead of independent oscillators, Peter Debye considered the collective motion of atoms as sound waves.

Sound waves

A sound wave is a collective motion of atoms through a solid. The displacement $\delta \mathbf{r}$ of an atom at position $\mathbf{r}$ and time $t$ is discribed by

$$\delta \mathbf{r}={\delta \mathbf{r}}_0{e}^{i(\mathbf{k}\cdot \mathbf{r}-\omega t)},$$

where ${\delta \mathbf{r}}_0$ is the amplitude of the wave and $\mathbf{k}=(k_x,k_y,k_z)$ the wave vector. The wavelength $\lambda$ is related to the wavevector $\mathbf{k}$ though $\lambda =2\pi /|\mathbf{k}|$.

The wave depends on time only through the factor $e^{-i\omega t}$. Therefore these waves are normal modes: oscillations of a system in which all parts of the system oscillate with the same frequency and fixed phase relation.

In addition to direction of the wave $k$, each sound wave has another degree of freedom: the direction in which the atoms themselves move or the wave polarization. Per wavevector $\mathbf{k}$ there are three modes in a 3D solid: two tranverse (perpendicular to $\mathbf{k}$) and one longitudinal mode (parallel to $\mathbf{k}$).

The space containing all possible values of $\mathbf{k}$ is called the $k$-space (also named the reciprocal space).

Despite having a position-dependent phase, each sound waves is an independent harmonic oscillator. The quantum mechanical excitations of this harmonic oscillator motion are called phonons—the particles of sound. Phonons are bosons and therefore their statistics is described by the Bose-Einstein distribution $n_B(\hslash \omega (\mathbf{k}))$. Debye used the description of phonons to model the heat capacity of solids.

The frequency of these phonons depends on its wavevector $\mathbf{k}$ through the dispersion relation

$$\omega =v_s|\mathbf{k}|,$$ where $v_s$ is the sound velocity of a material.

To summarize, instead of having $3N$ oscillators with the same frequency ${\omega }_0$, we now have $3N$ possible phonon modes with frequencies depending on $\mathbf{\text{k}}$ through the dispersion relation $\omega (\mathbf{k})=v_s|\mathbf{k}|$. The expected value of the total energy (which we, for simplicity, from now on will denote as the total energy) is given by the sum over the energy of all possible phonon modes characterized by a wavevector $\mathbf{k}$:

$$\begin{array}{rl}⟨E⟩& =3\sum _{\mathbf{k}}(\frac{1}{2}\hslash \omega (\mathbf{k})+\hslash \omega (\mathbf{k})n_{BE}(\hslash \omega (\mathbf{k})))\\ & =3\sum _{\mathbf{k}}(\frac{1}{2}\hslash \omega (\mathbf{k})+\frac{\hslash \omega (\mathbf{k})}{e^{\hslash \omega (\mathbf{k})/k_{B}T}-1}).\end{array}$$

Here we used that the expected occupation number is $n_B(\hslash \omega (\mathbf{k}))$.

Where does the factor 3 originate from?

The factor 3 originates from the three possible normal mode polarizations (directions in which the atoms can move) per wavevector $\mathbf{k}$.

We still have several open questions:

• Normal modes depend on the material's shape. What impact does this have on the heat capacity?
• Which $\mathbf{k}$ are possible and which are not?
• If all $\mathbf{k}$ are possible, shouldn't $E$ be infinite?

Periodic boundary conditions

We can answer all the above questions by realizing the following:

$C$ is a macroscopic property: it should not depend on the material's shape and should only be proportional to its volume. Therefore, we consider a material with a simple shape to make the calculation for $C$ easier.

The easiest option people have invented so far is a box of size $V=L^3$ with periodic boundary conditions².

Periodic boundary conditions require that the atomic displacement $\delta \mathbf{r}$ is periodic inside the material. Let us consider a translation by $L$ in the $x$-direction

$$\delta \mathbf{r}(\mathbf{r}+L\hat{\mathbf{x}})=\delta \mathbf{r}(\mathbf{r}).$$

Substituting the plane wave definition into the periodicity requirement, we get

$${\delta \mathbf{r}}_0{e}^{i(\mathbf{k}\cdot \mathbf{r}+k_{x}L-\omega t)}={\delta \mathbf{r}}_0{e}^{i(\mathbf{k}\cdot \mathbf{r}-\omega t)},$$

or after canceling common prefactors

$$e^{i{k}_{x}L}=e^{i0}=1.$$

We see that in order to satisfy the periodic boundary conditions, $k_x=n_x\frac{2\pi }{L}$, for $n_x\in \mathbb{Z}$. The same condition holds for the $y$- and $z$-direction. Hence the allowed values for $\mathbf{k}$ that satisfy the periodic boundary conditions are given by

$$\mathbf{k}=\frac{2\pi }{L}(n_x,n_y,n_z),\{n_x,n_y,n_z\}\in \mathbb{Z}.$$

We observe something very important: the periodic boundary conditions discretize $k$-space. The allowed values of $\mathbf{k}$ form a regular grid in $k$-space.

We observe that in 3D, there is exactly one allowed $\mathbf{\text{k}}$-value per volume ${\left(\frac{2\pi }{L}\right)}^3$ in $k$-space.

When we consider larger and larger box sizes, $L\to \mathrm{\infty }$, the volume per allowed mode becomes smaller and smaller. This implies that for a large enough $L$, we can approximate the sum over $\mathbf{\text{k}}$ as an integral:

$$\sum _{\mathbf{k}}\approx \frac{L^3}{(2\pi {)}^3}\int \text{d}\mathbf{\text{k}}$$

Integral over $k$-space

This conversion from a sum over a discrete grid of $k$-space states to a volume integral is one of the extremely commonly used ideas in solid state physics: it provides us a way to count all the possible waves.

We can use this approximation and rewrite the total energy as an integral:

$$\begin{array}{rlr}⟨E⟩& =3\sum _{\mathbf{k}}(\frac{1}{2}\hslash \omega (\mathbf{k})+\frac{\hslash \omega (\mathbf{k})}{e^{\hslash \omega (\mathbf{k})/k_{B}T}-1})& \approx 3\frac{L^3}{(2\pi {)}^3}\int \text{d}\mathbf{\text{k}}(\frac{1}{2}\hslash \omega (\mathbf{k})+\frac{\hslash \omega (\mathbf{k})}{e^{\hslash \omega (\mathbf{k})/k_{B}T}-1})\end{array}$$.

Where $\omega (\mathbf{k})=v_s|\mathbf{k}|$ is the dispersion relation, and the integral goes over the 3-dimensional $k$-space.

Density of states

The integrand of the total energy depends only on $|\mathbf{k}|$. Because of this symmetry, the integrand is convenient to rewrite in spherical coordinates.

$$\int \text{d}\mathbf{\text{k}}\to \underset{0}{\overset{2\pi }{\int }}\text{d}\phi \underset{0}{\overset{\pi }{\int }}\sin (\theta )\text{d}\theta \underset{0}{\overset{\mathrm{\infty }}{\int }}{k}^2\text{d}k=4\pi \underset{0}{\overset{\mathrm{\infty }}{\int }}{k}^2\text{d}k.$$

Performing the change of variables, we obtain the expression for the total energy in spherical coordinates is

$$\begin{array}{rl}⟨E⟩& =3\frac{L^3}{(2\pi {)}^3}{\int }_0^{\mathrm{\infty }}4\pi k^2(\frac{1}{2}\hslash \omega (k)+\frac{\hslash \omega (k)}{e^{\hslash \omega (k)/k_{B}T}-1})\text{d}k\\ & =3\frac{L^3}{(2\pi {)}^3}\frac{4\pi }{v_s^3}{\int }_0^{\mathrm{\infty }}{\omega }^2(\frac{1}{2}\hslash \omega +\frac{\hslash \omega }{e^{\hslash \omega /k_{B}T}-1})\text{d}\omega .\end{array}$$

We utilized the dispersion relation $\omega (k)=v_s|k|$ and omitted the absolute value of $k$ due to the integral over $k$ only running from 0 to $\mathrm{\infty }$ after conversion to spherical coordinates. The integral above can be split up into two factors. The factor inside the brackets describes the average energy of a phonon mode with frequency $\omega$. The other factor is the density of states $g(\omega )$.

Definition: density of states

The density of states $g(\omega )$ counts the total number of available normal modes inside a frequency interval $\omega +\text{d}\omega$. In other words, the number of modes in this interval equals $g(\omega )\text{d}\omega$.

With this definition, the integral becomes

$$E=\underset{0}{\overset{\mathrm{\infty }}{\int }}\text{d}\omega (\frac{1}{2}\hslash \omega +\frac{\hslash \omega }{e^{\hslash \omega /k_{B}T}-1})g(\omega ),$$

with

$$g(\omega )=3{\left(\frac{L}{2\pi }\right)}^3\frac{4\pi {\omega }^2}{v_s^3}.$$

We interpret the integral above as follows: we multiply the number of modes $g(\omega )$ by the average energy of a single mode at a given frequency $\omega$ and integrate over all frequencies.

Let us separate $g(\omega )$ into a product of individual factors:

• $3$ comes from the number of possible polarizations in 3D (two transversal, one longitudinal).
• $(\frac{L}{2\pi }{)}^3$ is the density of $\mathbf{\text{k}}$ points in $k$-space.
• $4\pi$ is the area of a unit sphere.
• ${\omega }^2$ is due to the area of a sphere in $k$-space being proportional to its squared radius $k^2$ and by having a linear dispersion relation $\omega =v_{s}k$.
• $v_s^{-3}$ is from the linear dispersion relation $\omega =v_{s}k$.

So in our case, due to the spherical symmetry, $g(\omega )\text{d}\omega$ can be obtained by calculating the density of states of a volume element $dV=4\pi k^{2}dk$ in $k$-space and substituting the dispersion relation $\omega (k)$.

Low temperatures

Because $g(\omega )$ and $\hslash \omega$ do not depend on temperature, we split up the integral of the total energy to temperature-dependent and temperature-independent parts:

$$⟨E⟩=E_{\text{Z}}+\frac{3L^3}{2{\pi }^2{v}_s^3}\underset{0}{\overset{\mathrm{\infty }}{\int }}\left(\frac{\hslash \omega }{e^{\hslash \omega /k_{B}T}-1}\right){\omega }^2\text{d}\omega .$$

The term $E_{\text{Z}}$ is the temperature-independent zero-point energy of the phonon modes. Despite $E_{\text{Z}}$ diverging towards infinity, does not contribute to $C$.

The integral depends on the temperature through the $e^{\hslash \omega /k_{B}T}$ term. In order evaluate the integral, we substitute $x\equiv \frac{\hslash \omega }{k_{B}T}$ and remove the temperature dependence of the integrand:

$$\begin{array}{rl}⟨E⟩& =E_Z+\frac{3L^3}{2{\pi }^2{v}_s^3}\frac{{\left(k_{B}T\right)}^4}{{\hslash }^3}\underset{0}{\overset{\mathrm{\infty }}{\int }}\frac{x^3}{{\text{e}}^x-1}\text{d}x\\ & =E_Z+\frac{3L^3}{2{\pi }^2{v}_s^3}\frac{{\left(k_{B}T\right)}^4}{{\hslash }^3}\frac{{\pi }^4}{15},\end{array}$$

where we used the fact that the integral is equal to $\frac{{\pi }^4}{15}$⁴. As we can see, the energy scales as $T^4$. Therefore we can conclude that

$$C=\frac{\text{d}E}{\text{d}T}\propto T^3.$$

We recover the empirical $T^3$ dependence of $C$ at low temperatures!

Can we understand this without going through the integration? Turns out we can!

1. At temperature $T$, only phonon modes with an energy below the thermal energy $E_{\text{T}}=k_{B}T$ become thermally excited. These modes have $\hslash \omega (\mathbf{k})\lesssim k_{B}T$.

2. By substituting the dispersion relation into the above inequality, we conclude that these modes have wave vectors $|\mathbf{k}|\lesssim k_{B}T/\hslash v_s$. Therefore the number of excited modes is proportional to the volume of a sphere $V_{\mathbf{\text{k}}}=\frac{4\pi }{3}|\mathbf{k}{|}^3$, multiplied by the density of modes in $k$-space, ${\left(\frac{L}{2\pi }\right)}^3$. Thus the total number of excited modes is

   $$\begin{array}{rl}{N}_{\text{modes}}& =V_{\mathbf{\text{k}}}{\left(\frac{L}{2\pi }\right)}^3\\ & \sim {\left(|\mathbf{k}|L\right)}^3\\ & \sim (k_{B}TL/\hslash v_s{)}^3\end{array}$$

   where we have substituted $|\mathbf{k}|\simeq k_{B}T/\hslash v_s$ and left out all numerical factors.

3. When thermally excited, the motion of these modes resembles that of classical harmonic oscillators. Therefore, each mode contributes $k_B$ to the heat capacity (Equipartition theorem). As a result, the heat capacity is

   $$\begin{array}{rl}C& =N_{\text{modes}}{k}_B\\ & \propto k_B(k_{B}TL/\hslash v_s{)}^3.\end{array}$$

Debye's interpolation for medium $T$

We observed that the above approximation yields a correct scaling of the heat capacity at low temperatures. We also know that for a 3D material with $N$ atoms, $C$ should converge to $3N{k}_B$ at high temperatures (the law of Dulong–Petit). However at high temperatures, phonon modes with all values of $\omega$ become thermally excited, and the number of these modes tends towards infinity:

$$N_{\text{modes}}={\int }_0^{\mathrm{\infty }}\text{d}\omega g(\omega )\to \mathrm{\infty }$$

As a result, we now incorrectly predict that the heat capacity also goes to infinity $C\propto N_{\text{modes}}{k}_B\to \mathrm{\infty }$. Hence, the model breaks down for high temperatures.

To fix this problem Debye realised that there should be as many phonon modes in the system as there are degrees of freedom. In a 3D material with $N$ atoms, there are a total of $3N$ normal modes, and not more.

In view of this fact, Debye proposed a fix to the problem: assume that there is a maximal frequency ${\omega }_D$ (Debye frequency), beyond which there are no phonons. We have no good justification for this assumption yet, but it is reasonable because the atoms certainly cannot move with infinite frequency.

$$g(\omega )=\{\begin{array}{ll}\frac{3L^3{\omega }^2}{2{\pi }^2{v}_s^3}& \omega \le {\omega }_D\\ 0& \omega >{\omega }_D\end{array}$$

Let us now compute ${\omega }_D$. We know that for a 3D system with $N$ atoms has to have exactly $3N$ phonon modes.

$$\begin{array}{rl}3N& ={\int }_0^{{\omega }_D}g(\omega )d\omega =\frac{3L^3}{2{\pi }^2{v}_s^3}{\int }_0^{{\omega }_D}{\omega }^2\text{d}\omega \\ & =\frac{L^3{\omega }_D^3}{2{\pi }^2{v}_s^3},\end{array}$$

which gives us

$${\omega }_D=v_s(6{\pi }^{2}N{)}^{1/3}/L.$$

Both $N$ and $L$ are arbitrary, however we are considering an $L\times L\times L$ box with $N$ atoms, so $L/N^{1/3}$ is the distance between neighboring atoms, and therefore ${\omega }_D$ does not depend on the box size.

Using the corrected expression for the total energy that includes the high frequency cut-off, the total energy without the zero-point motion part is

$$⟨E⟩={\int }_0^{{\omega }_D}g(\omega )\frac{\hslash \omega }{e^{\hslash \omega /k_{B}T}-1}\text{d}\omega .$$

We now substitute the previously calculated density of states $g(\omega )$ and the Debye frequency ${\omega }_D$ into the equation of the total energy. Then the heat capacity yields

$$C=9N{k}_B{\left(\frac{T}{T_D}\right)}^3{\int }_0^{T_D/T}\frac{x^4{\text{e}}^x}{({\text{e}}^x-1{)}^2}\text{d}x,$$

where we made use of the substitution $x\equiv \frac{\hslash \omega }{k_{B}T}$ and defined the Debye temperature $T_D\equiv \frac{\hslash {\omega }_D}{k_B}$. Similar to all the isolated harmonic oscillators in the Einstein model becoming thermally excited when $T\gtrsim T_E$, when $T\gtrsim T_D$, all the phonon modes in a Debye solid become thermally excited. The number of phonons in each mode will keep on increasing with $T$ as described by the Bose-Einstein distribution, scaling linearly with $T$ when $k_{B}T\gg \hslash \omega$.

Below is once again the plot of the measured heat capacity of silver fitted by the Einstein model and the Debye model.

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]

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 fitted by 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');

Debye model clearly wins!

Conclusions

1. The Debye model assumes that atoms in materials move in a collective fashion, described by quantized normal modes with a dispersion relation $\omega =v_s|\mathbf{k}|$.
2. The phonon modes have a constant density of $(L/2\pi {)}^3$ in the reciprocal / $k$-space.
3. The total energy and heat capacity are obtained by integrating the contribution of the individual modes over $k$-space.
4. The density of states $g(\omega )$ is the number of states per frequency. With a dispersion relation $\omega =v_s|\mathbf{k}|$, $g(\omega )$ is proportional to ${\omega }^2$ for a 3D bosonic system.
5. At low temperatures the phonon heat capacity is proportional to $T^3$.
6. Phonon modes only exist up until the Debye frequency ${\omega }_D$, after which there are no modes in the system.

Exercises

Quick warm-up exercises

1. Express the heat capacity in the low-$T$ limit in terms of $T_D$.
2. Make a sketch of the heat capacity in the low-$T$ limit for two different Debye temperatures.
3. Why are there only 3 polarizations when there are 6 degrees of freedom in three-dimensions for an oscillator?
4. Express the two-dimensional integral $\int \mathrm{d}{k}_x\mathrm{d}{k}_y$ in terms of polar coordinates. You can assume rotational symmetry.
5. The Einstein model has a material-dependent frequency ${\omega }_0=k_{\mathrm{B}}{T}_E/\hslash$ of the quantum harmonic oscillators as a free fitting parameter. What is the material-dependent parameter that plays a similar role in the Debye model?
6. Derive an expression for the shortest possible wavelength in the Debye model it in terms of the interatomic distance $a$. Hint: assume that the number of atoms is given by $N=V/a^3$. Discuss if the answer is reasonable.

Exercise 1: Debye model: concepts

Consider the probability to find an atom of a 1D solid that originally had a position $x$ at a displacement $\delta x$ shown below:

def psi_squared(delta_x, x):
    factor = np.sin(4*np.pi*x)**2 + .001
    return delta_x**2 * np.exp(-delta_x**2 / factor) / factor

x = np.linspace(0, 1, 200)
delta_x = np.linspace(-2, 2, 200)

# Now to plotting
pyplot.figure()
ax = pyplot.gca()
im = ax.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((), ())
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
cbar = pyplot.colorbar(im, cax=cax)
cbar.set_ticks(())
cbar.set_label(r'$|\psi|^2$')

1. Describe which $k$-states contain a phonon. Explain your answer.

   Hint

   There are two $k$-states which contain a phonon.

2. Describe the concept of $k$-space. What momenta are allowed in a 2D system with dimensions $L\times L$?

3. Explain the concept of density of states.
4. Calculate the phonon density of states $g(\omega )$ of a 3D, 2D and 1D solid with linear dispersion $\omega =v_s|\mathbf{k}|$.

Exercise 2: Debye model in 2D

1. State the assumptions of the Debye model.
2. Determine the energy of a two-dimensional solid as a function of $T$ using the Debye approximation. You do not have to solve the integral.
3. Calculate the heat capacity in the high $T$ limit.
4. At low $T$, show that $C_V=K{T}^n$. Find $n$. Express $K$ as an indefinite integral (similarly to what done during the lecture)³.

Exercise 3: Different phonon modes

(adapted from ex 2.6a of "The Oxford Solid State Basics" by S.Simon)

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_{\mathrm{\perp }}|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$.

Exercise 4: Anisotropic sound velocities

(adapted from ex 2.6b of "The Oxford Solid State Basics" by S.Simon)

Suppose now that the velocity is anisotropic ($v_x\ne v_y\ne 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.

---

1. Data is taken from C. Kittel, Solid State Physics, 2ed Wiley (1956). ↩

2. An alternative way to treat the boundaries is by using fixed boundary conditions (like a guitar string), resulting in standing waves with $k=\pi /L$, $2\pi /L$, $3\pi /L$, …. This gives the same answer, but it usually more ugly, and takes more work. ↩

3. This integral evaluates to the famous Riemann zeta function (See Chapter 2.3 of the book for more details). ↩

4. The integrand can be solved by reducing it to the Riemann zeta functions and then solving the remaining new integral (see page 12 of the book). ↩