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
Approximate integrals with a small parameter. To check yourself, answer how \(\int_0^1 \sqrt{\exp(-\alpha x)} dx\) depends on \(\alpha\) when \(\alpha \to \infty\).

Learning goals

After this lecture you will be able to:

Describe the concept of reciprocal space and allowed wave vectors
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 measured heat capacity 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();

svg

We see that 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 scale with temperature at low \(T\)?

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

The Debye model¶

The key simplification of the Einstein model is to consider the atoms as independent quantum harmonic oscillators. Instead of considering each atom as an independent harmonic oscillator, Peter Debye considered the sound waves in a material - the collective motion of atoms - as independent harmonic oscillators.

Sound waves¶

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

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

where \(\mathbf{\delta 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).

The sound waves are characterized by their frequency, wavevector, and polarization. The frequency of the sound wave is related to its wavevector \(\mathbf{k}\) through the dispersion relation

\[ \omega = v_s|\mathbf{k}|, \]

where \(v_s\) is the sound velocity of a material.

As discussed in the previous lecture, the quantum mechanical excitations of a harmonic oscillator are called phonons, and the expected number of phonons in the oscillator at temperature \(T\) is given by the Bose-Einstein distribution \(n_B(\beta \hbar \omega(\mathbf{k}))\). Instead of having \(3N\) oscillators with the same frequency \(\omega_0\) as in the Einstein model, we now have \(3N\) oscillators (the vibrational modes) with frequencies depending on \(\textbf{k}\) through the dispersion relation \(\omega(\mathbf{k}) = v_s|\mathbf{k}|\). Apart from this crucial difference with the Einstein model, the calculation of the expectation value of the energy stored in the vibrational modes proceeds in the same way: The expected value of the total energy stored in the oscillators (which, from now on, we will simply denote as the total energy \(E\)) is given by the sum of the energy stored in the individual oscillators. These oscillators are characterized by their wavevector \(\mathbf{k}\):

\[\begin{align*} E &= 3 \sum_\mathbf{k} \left(\frac{1}{2}\hbar\omega(\mathbf{k})+\hbar \omega(\mathbf{k}) n_{B}(\beta \hbar \omega(\mathbf{k}))\right)\\ &= 3 \sum_\mathbf{k} \left(\frac{1}{2}\hbar\omega(\mathbf{k})+\frac{\hbar\omega(\mathbf{k})}{ e^{\hbar\omega(\mathbf{k})/{k_BT}}-1}\right). \end{align*}\]

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

Where does the factor 3 come 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:

The heat capacity \(C\) is a macroscopic property: it should not depend on the material's shape and only be proportional to its volume. Therefore, we consider a material with a simple shape to make the calculation of \(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 \(\mathbf{\delta r}\) is periodic inside the material. Let us consider a translation by \(L\) in the \(x\)-direction

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

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

\[ \mathbf{\delta r}_0 e^{i(\mathbf{k} \cdot \mathbf{r}+k_xL-\omega t)} = \mathbf{\delta r}_0 e^{i(\mathbf{k} \cdot \mathbf{r}-\omega t)}, \]

or after canceling common prefactors

\[ e^{i k_x L} = e^{i 0} = 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)\text{, where }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 \({\textbf k}\)-value per volume \(\left(\frac{2\pi}{L}\right)^3\) in \(k\)-space.

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

\[ \sum_\mathbf{k} \approx \frac{L^3}{(2\pi)^3}\int \textrm{d} \textbf{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 to rewrite the total energy as an integral:

\[\begin{align*} E &= 3 \sum_\mathbf{k} \left(\frac{1}{2}\hbar\omega(\mathbf{k})+\frac{\hbar\omega(\mathbf{k})}{ e^{\hbar\omega(\mathbf{k})/k_B T}-1}\right) \\ &\approx 3 \left(\frac{L}{2\pi}\right)^3 \int \textrm{d} \textbf{k} \left(\frac{1}{2}\hbar\omega(\mathbf{k})+\frac{\hbar\omega(\mathbf{k})}{e^{\hbar\omega(\mathbf{k})/k_B T}-1}\right). \end{align*}\]

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

Density of states¶

Because the dispersion is \(\omega(\mathbf{k}) = v_s |\mathbf{k}|\), the integrand of the total energy depends only on the length \(|\mathbf{k}|\) of the wavevector. Therefore, it is convenient to rewrite the integral in spherical coordinates.

\[ \int \textrm{d} \textbf{k} \to \int\limits_0^{2\pi}\textrm{d}\varphi\int\limits_0^{\pi} \sin(\theta) \: \textrm{d}\theta \int\limits_0^\infty k^2 \textrm{d}k = 4 \pi \int \limits_{0}^{\infty} k^2 \textrm{d} k. \]

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

\[\begin{align*} E &= 3 \frac{L^3}{(2\pi)^3}\int_0^\infty 4 \pi k^2 \left(\frac{1}{2}\hbar\omega(k)+\frac{\hbar\omega(k)}{ {e}^{\hbar\omega(k)/{k_BT}}-1}\right) \textrm{d}k \\ &= 3 \frac{L^3}{(2\pi)^3}\frac{4 \pi}{v_s^3}\int_0^\infty \omega^2 \left(\frac{1}{2}\hbar\omega+\frac{\hbar\omega}{ e^{\hbar\omega/{k_BT}}-1}\right) \textrm{d}\omega. \end{align*}\]

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 \(\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 + \textrm{d} \omega\). In other words, the number of modes in this interval equals \(g(\omega) \textrm{d} \omega\).

With this definition, the integral becomes

\[ E = \int\limits_0^{\infty}\textrm{d} \omega \left(\frac{1}{2}\hbar\omega+\frac{\hbar\omega}{ e^{\hbar\omega/{k_BT}}-1}\right)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 \(\textbf{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_sk\).
\(v_s^{-3}\) is from the linear dispersion relation \(\omega = v_sk\).

So in our case, due to the spherical symmetry, \(g(\omega)\textrm{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 \(\hbar \omega\) do not depend on temperature, we split up the integral of the total energy to temperature-dependent and temperature-independent parts:

\[ E = E_{\textrm Z} + \frac{3L^3}{2\pi^2 v_s^3}\int\limits_0^\infty \left(\frac{\hbar\omega}{ e^{\hbar\omega/k_BT}-1}\right)\omega^2 {\textrm{d}}\omega. \]

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

The integral depends on the temperature through the \(e^{\hbar\omega/k_BT}\) term. In order evaluate the integral, we substitute \(x\equiv\frac{\hbar\omega}{k_BT}\) and remove the temperature dependence of the integrand:

\[\begin{align*} E &= E_Z + \frac{3L^3}{2\pi^2 v_s^3}\frac{\left(k_BT\right)^4}{\hbar^3}\int\limits_0^\infty\frac{x^3}{ {\textrm e}^x-1} {\textrm{d}}x \\ &= E_Z + \frac{3L^3}{2\pi^2 v_s^3}\frac{\left(k_BT\right)^4}{\hbar^3}\frac{\pi^4}{15}, \end{align*}\]

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{ {\textrm{d}}E}{ {\textrm{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!

At temperature \(T\), only phonon modes with an energy below the thermal energy \(E_{\textrm{T}} = k_B T\) become thermally excited. These modes have \(\hbar \omega(\mathbf{k}) \lesssim k_B T\).
By substituting the dispersion relation into the above inequality, we conclude that these modes have wave vectors \(|\mathbf{k}| \lesssim k_B T /\hbar v_s\). Therefore the number of excited modes is proportional to the volume of a sphere \(V_{\textbf{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{align*} N_\textrm{modes} &= V_{\textbf{k}} \left(\frac{L}{2\pi}\right)^3\\ &\sim \left( |\mathbf{k}| L \right)^3\\ &\sim (k_B T L/\hbar v_s)^3 \end{align*}\]

where we have substituted \(|\mathbf{k}| \simeq k_B T /\hbar v_s\) and left out all numerical factors.
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{align*} C &= N_{\textrm{modes}} k_B \\ &\propto k_B (k_B T L/\hbar v_s)^3. \end{align*}\]

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 \(3Nk_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_{\textrm{modes}} = \int_0^{\infty} \textrm{d} \omega g(\omega) \to \infty \]

As a result, we now incorrectly predict that the heat capacity also goes to infinity \(C \propto N_{\textrm{modes}} k_B \to \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) = \left\{ \begin{array}{ll} \frac{3L^3\omega^2}{2\pi^2v_s^3} & \omega \leq \omega_D \\ 0 & \omega>\omega_D \end{array} \right. \]

Let us now compute \(\omega_D\). We know that for a 3D system with \(N\) atoms has to have exactly \(3N\) phonon modes.

\[\begin{align*} 3N &= \int_0^{\omega_D} g(\omega)d\omega = \frac{3L^3}{2 \pi^2 v_s^3} \int_0^{\omega_D}\omega^2 {\textrm{d}}\omega\\ &= \frac{L^3\omega_D^3}{2\pi^2v_s^3}, \end{align*}\]

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×L×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{\hbar \omega}{e^{\hbar\omega/k_BT}-1} \textrm{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=9Nk_B\left(\frac{T}{T_{D}}\right)^3\int_0^{T_{D}/T}\frac{x^4{\textrm e}^x}{({\textrm e}^x-1)^2}{\textrm{d}}x, \]

where we made use of the substitution \(x\equiv\frac{\hbar\omega}{k_BT}\) and defined the Debye temperature \(T_{D}\equiv\frac{\hbar\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_BT \gg \hbar\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');

svg

Debye model clearly wins!

Conclusions¶

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}|\).
The phonon modes have a constant density of \((L/2\pi)^3\) in the reciprocal / \(k\)-space.
The total energy and heat capacity are obtained by integrating the contribution of the individual modes over \(k\)-space.
The density of states \(g(\omega)\) is the number of states per frequency. With a dispersion relation \(ω = v_s|\mathbf{k}|\), \(g(\omega)\) is proportional to \(\omega^2\) for a 3D bosonic system.
At low temperatures the phonon heat capacity is proportional to \(T^3\).
Phonon modes only exist up until the Debye frequency \(\omega_D\), after which there are no modes in the system.

Exercises¶

Exercises with an asterisk (*) are considered to be at the essential/basic level

Warm-up exercises*¶

Express the heat capacity of the 3D Debye model in the low-\(T\) limit in terms of \(T_D\).
Make a sketch of the heat capacity in the low-\(T\) limit for two different Debye temperatures.
Discuss the concept of polarization and what directions it can have for a wave in three dimensions.
Express the two-dimensional integral over k-space \(\int\mathrm{d}k_x\mathrm{d}k_y\) in terms of polar coordinates.
The Einstein model has a material-dependent frequency \(\omega_0 = k_\mathrm{B} T_E/\hbar\) 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? What does this parameter represent?
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*: Deriving the density of states for the linear dispersion relation of the Debye model¶

In this lecture, we found that the linear dispersion considered in the Debye model yields a density of states \(g(\omega)\propto \omega^2\) (in three dimensions). In this exercise, we will practice this important derivation again and extend it to 1D and 2D.

Write down the dispersion relation of the vibrational modes in the Debye model.
Assume periodic boundary conditions. What is the distance between nearest-neighbour points in \(\mathbf{k}\)-space? What is the density of \(\mathbf{k}\)-points in 1, 2, and 3 dimensions?
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.
Transform these integrals into integrals over frequency for 1D, 2D and 3D using the dispersion relation. Indicate the integral boundaries. Extract the density of states. Are the integral boundaries important for the result?
Discuss the concept of density of states.

Exercise 2*: Debye model in 2D¶

Here, we analyze the vibrational energy and heat capacity of a two-dimensional Debye solid as a function of temperature.

Formulate an integral expression for the energy stored in the vibrational modes of a two-dimensional Debye solid as a function of \(T\).
Calculate the heat capacity in the high \(T\) limit.
At low \(T\), show that \(C_V=KT^{n}\). Find \(n\). Express \(K\) as an indefinite integral (similarly to what was done during the lecture)³.

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

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

In the lecture, we derived the low-temperature Debye heat capacity assuming that all the vibrational modes have the same sound velocity \(v_s\). In reality, 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 |\mathbf{k}|\), and two transverse modes with \(\omega = v_\bot |\mathbf{k}|\)

Write down the total energy of the thermally excited phonons in this material (hint: use the same reasoning as in the Lithium exercise).
Verify that at high \(T\) you reproduce the Dulong-Petit law.
Compute the behavior of the heat capacity at low \(T\).

Exercise 4: Anisotropic sound velocities¶

(adapted from ex 2.6b of "The Oxford Solid State Basics" by S.Simon) We consider a 3D solid in which the sound velocity is anisotropic. I.e., the sound velocity depends on the direction in which the wave travels: (\(v_x \neq v_y \neq v_z\)). As a result, the dispersion is \(\omega = \sqrt{v_x^2 k_x^2 + v_y^2 k_y^2 + v_z^2 k_z^2}\).

Express the thermal energy stored in the vibrational modes of the system as an integral over k-space.
Make a variable substitution to restore the spherical symmetry of the integrand. Express the integral in spherical coordinates.
Formulate an integral expression for the heat capacity. Make a variable substitution to extract the temperature of the heat capacity in the low-\(T\) limit (as done in the lecture notes). Is the scaling with temperature different from the case of isotropic sound velocities?

Data is taken from C. Kittel, Solid State Physics, 2ed Wiley (1956). ↩
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. ↩
This integral evaluates to the famous Riemann zeta function (See Chapter 2.3 of the book for more details). ↩
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). ↩