Advanced Solid State Physics



Magnetic effects and
Fermi surfaces


Linear response


Crystal Physics



Structural phase

Landau theory
of second order
phase transitions




Exam questions




Course notes

TUG students


Phonon contribution to thermal conductivity

Phonons transport heat from hotter regions to cooler regions. The heat flow can be calculated from the phonon dispersion relation of a crystal using the Boltzmann transport equation. The phonon dispersion relation of a crystal consists of $3N$ branches where $N$ is the number of atoms in the basis. The phonon normal modes are labeled by $(\vec{k},p)$ where $\vec{k}$ is the wave number and $p$ specifies the branch of the dispersion relation. The branches are either acoustic branches or optical branches and they each have a polarization, either transverse or longitudinal. Since phonons are bosons, the number of phonons in a $(\vec{k},p)$ mode in equilibrium is given by the Bose-Einstein distribution,

\[ \begin{equation} f_{BE} = \frac{1}{\exp\left(\frac{E-\mu}{k_B T}\right)-1}. \end{equation}\]

In thermal equilibrium, the states $\vec{k}$ and $-\vec{k}$ are equally occupied and the heat current is zero. If the system is pushed out of thermal equilibrium by applying a temperature gradient, the occupation of states $\vec{k}$ and $-\vec{k}$ can be different and the heat current density can be calculated by summing over the contributions of all phonon modes,

\[\begin{equation} \vec{j}_Q(\vec{r})=\sum_p \int D(\vec{k},p)\left(E(\vec{k},p)-\mu\right)\vec{v}_g(\vec{k},p)f(\vec{k},\vec{r}, p) d^3k. \label{eq:jQ} \end{equation}\]

Here $D(\vec{k},p)$ is the density of states, $E(\vec{k},p)$ is the energy of mode $(\vec{k},p)$, $\vec{v}_g=\frac{\nabla_{\vec{k}}E}{\hbar}$ is the group velocity of that mode and $f(\vec{k},\vec{r}, p)$ is the occupation function that tells us how many phonons occupy mode $(\vec{k},p)$ on average. Phonons are not conserved. They can be created or annihilated without a change of energy. Consequently, the chemical potential of phonons is zero, $\mu=0$ and because $\vec{j}_Q = \vec{j}_E - \mu\vec{j}_n$, the heat current is the same as the energy current,

\[\begin{equation} \vec{j}_Q(\vec{r})=\vec{j}_E(\vec{r})=\sum_p \int D(\vec{k},p)E(\vec{k},p)\vec{v}_g(\vec{k},p)f(\vec{k},\vec{r}, p) d^3k. \end{equation}\]

The allowed $\vec{k}$ states are uniformly distributed in reciprocal space so the density of states is a constant. In three dimensions the density of states is,

\[ \begin{equation} D(\vec{k},p)=\frac{1}{(2\pi)^3}. \end{equation}\]

Generally, the Boltzmann equation for the occupation function is,

\[\begin{equation} \frac{\partial f}{\partial t} = - \frac{1}{\hbar}\vec{F}_{\text{ext}}\cdot\nabla_{\vec{k}}f-\vec{v}_g\cdot\nabla_{\vec{r}}f + \frac{\partial f}{\partial t} \bigg\rvert_{collisions}. \end{equation}\]

Since we do not have to consider external forces for phonons the Boltzmann equation simplifies to

\[\begin{equation} \frac{\partial f}{\partial t} + \vec{v}_g \cdot \nabla f = \frac{\partial f}{\partial t}\bigg|_{collisions}. \label{eq:bte} \end{equation}\]

The collision term includes electron-phonon scattering, phonon-phonon scattering, phonon-defect scattering, and photon - phonon scattering. Although it is possible to use Fermi's golden rule to calculate the transition rates between the $(\vec{k},p)$ states caused by these different scattering events, a common approximation is the relaxation time approximation where the collision term is approximated as,

\[\begin{equation} \frac{\partial f}{\partial t}\bigg|_{collisions}=\frac{f_{BE} - f}{\tau}. \label{eq:rta} \end{equation}\]

Here $\tau$ is the relaxation time and describes the time it would take for the system to return to equilibrium if all external sources were removed. In this approximation the Boltzmann equation is,

\[\begin{equation} \frac{\partial f}{\partial t} + \vec{v}_g \cdot \nabla_{\vec{r}} f = \frac{f_{BE} - f}{\tau}. \label{eq:phpt} \end{equation}\]

Under steady-state conditions, $\frac{\partial f}{\partial t} = 0$, and the occupation function is,

\[\begin{equation} f(\vec{k},\vec{r},p) = f_{BE} - \tau \vec{v}_g \cdot \nabla_{\vec{r}} f. \end{equation}\]

When the system is close to equilibrium, the occupation probability $f$ will be nearly the same as the Bose-Einstein factor $f_{BE}$. To a good approximation $f$ can be replaced by $f_{BE}$ on the righthand side of the equation,

\[\begin{equation} f(\vec{k},\vec{r},p) \approx f_{BE} - \tau \vec{v}_g \cdot \nabla_{\vec{r}} f_{BE}. \label{eq:f0} \end{equation}\]

When evaluating $\nabla_{\vec{r}} f_{BE}$, we have to take into account that the temperature depends on the position $\vec{r}$. Also, we assume that the material is uniform in space, so the dispersion relation does not depend on $\vec{r}$,

\[\begin{equation} \nabla_{\vec{r}} f_{BE}(\vec{k},\vec{r}) = \frac{\partial f_{BE}}{\partial T} \nabla_{\vec{r}} T = \frac{ E \;\exp\left(\frac{E}{k_B T}\right)}{k_B T^2 \left(\exp\left(\frac{E}{k_B T}\right)-1 \right)^2 } \nabla_{\vec{r}} T = \frac{E\;\exp\left(\frac{E}{k_B T}\right)}{k_B T^2} \; f_{BE}^2 \nabla_{\vec{r}} T. \end{equation}\]

Plugging this back into equation (\ref{eq:f0}) and using $\vec{v}_g=\frac{\nabla_{\vec{k}}E}{\hbar}$ yields,

\[\begin{equation} f(\vec{k},\vec{r},p) \approx f_{BE} - \frac{\tau}{\hbar} \; \frac{E\;\exp\left(\frac{E}{k_B T}\right)}{k_B T^2} \; f^2_{BE} \nabla_{\vec{k}}E\cdot \nabla_{\vec{r}} T. \end{equation}\]

This occupation probability is substituted into (\ref{eq:jQ}) and the integration is performed over $\vec{k}$. After integrating over $\vec{k}$, the first term $f_{BE}$ yields zero since there are as many right moving states as left moving states in thermal equilibrium. The heat current is,

\[\begin{equation} \vec{j}_Q(\vec{r})= - \sum\limits_p \int \frac{\tau}{(2\pi)^3\hbar^2} \; \frac{E^2\;\exp\left(\frac{E}{k_B T}\right)}{k_B T^2} \; f^2_{BE} \nabla_{\vec{k}}E\cdot \nabla_{\vec{r}} T \nabla_{\vec{k}}E d^3k. \end{equation}\]

The heat current can be written in terms of the thermal conductivity $K$,

\[\begin{equation} \vec{j}_Q(\vec{r})= - K\nabla T, \end{equation}\]

where $K$ is a second rank tensor,

\[\begin{equation} K_{ij}= \sum\limits_p \int \frac{\tau}{(2\pi)^3\hbar^2} \; \frac{E^2\;\exp\left(\frac{E}{k_B T}\right)}{k_B T^2} \; f^2_{BE} \nabla_{\vec{k}}E\cdot \hat{e}_i \nabla_{\vec{k}}E\cdot \hat{e}_j d^3k. \end{equation}\]

For cubic crystals the thermal conductivity is a constant,

\[\begin{equation} K= \sum\limits_p \int \frac{\tau}{(2\pi)^3\hbar^2} \; \frac{E^2\;\exp\left(\frac{E}{k_B T}\right)}{k_B T^2} \; f^2_{BE} \left(\nabla_{\vec{k}}E\cdot \hat{z}\right)^2 d^3k. \end{equation}\]

Consider an acoustic phonon branch with a linear dispersion relation like in the Debye model, $E=\hbar c|\vec{k}|$ where $c$ is the speed of sound. In the Debye model, the highest phonon frequency, $\omega_D$, is related to the atomic density $n$ as $\omega_D^3 = 6\pi^2nc^3$. The Debye temperature is defined as, $\hbar \omega_D = k_BT_D$. In this model the thermal conductivity can be written as,

\[\begin{equation} K= \frac{\tau}{(2\pi)^3\hbar^2}\int\limits_0^{k_D}\int\limits_0^{\pi}\int\limits_0^{2\pi} \; \frac{\hbar^4c^4k^4\;\exp\left(\frac{\hbar ck}{k_B T}\right)}{k_B T^2\left( \exp\left(\frac{\hbar ck}{k_B T}\right)-1\right)^2} \; \sin\theta dkd\theta d\phi. \end{equation}\]

Here $k_D = \omega_D/c$ is the Debye wavenumber corresponding to the shortest phonon wavelength. The integration over the angular variables contributes a factor of $4\pi$.

\[\begin{equation} K= \frac{2\hbar^2c^4\tau}{(2\pi)^2}\int\limits_0^{k_D} \frac{k^4\;\exp\left(\frac{\hbar ck}{k_B T}\right)}{k_B T^2\left( \exp\left(\frac{\hbar ck}{k_B T}\right)-1\right)^2} \; dk. \end{equation}\]

Make a substitution $x = \frac{\hbar ck}{k_B T}$, $dx = \frac{\hbar c}{k_B T}dk$, $x_D= \frac{\hbar ck_D}{k_B T}= \frac{T_D}{T}$.

\[\begin{equation} K= \frac{2\tau}{(2\pi)^2}\frac{k_B^4T^3}{\hbar^3c}\int\limits_0^{x_D} \frac{x^4\;\exp\left(x\right)}{\left( \exp\left(x\right)-1\right)^2} \; dx. \end{equation}\]

The integral can be evaluated numerically,

$\frac{(2\pi)^2\hbar^3c}{2\tau k_B^4T_D^3}K$


The thermal conductivity due to phonons using a single
linear acoustic mode and a single relaxation time.

As $T\rightarrow 0$, $x_D \rightarrow \infty$ and the integral over $x$ can be evaluated.

$$\int\limits_0^{\infty}\frac{x^4\exp(x)}{\left(\exp(x)-1\right)^2}dx = \frac{4\pi^4}{15}$$

For low temperatures, the thermal conductivity has the form,

\[\begin{equation} K\approx \frac{2\pi^2\tau}{15}\frac{k_B^4T^3}{\hbar^3c} \qquad T < < T_D. \end{equation}\]

At high temperatures, this calculation results in a constant thermal conductivity but in usually the thermal conductivity falls like $1/T$. The problem is that a single relaxation time was assumed for all temperatures. At higher temperatures the concentration of phonons increases like $T$ and the relaxation times becomes shorter like $1/T$ due to phonon-phonon scattering. This results in a thermal conductivity of the form $K \sim 1/T$.

A very simple model for the phonon dispersion of an optical branch is that the energy is a constant. This results in a group velocity that is zero and there is no contribution to the thermal conductivity due optical phonons.

Below the thermal conductivity of some elements can be plotted below .

$K$ [W/cm K]

$T$ [K]

The thermal conductivity of silicon.[1]