 ## 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$ $T/T_D$

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