Theory: Boltzmann Transport Equation

In the phonon picture, thermal transport can be modeled in the limit of small anharmonicity with perturbation theory. Given a small temperature gradient \(\nabla T\), the phonon heat flux is given by the phonon energy and phonon group velocity \(\mathbf{v}_{\lambda}\),

\[ \mathbf{J} = \frac{1}{V}\sum_{\lambda} \tilde{n}_{\lambda} \hbar \omega_{\lambda} \mathbf{v}_{\lambda} \]

where the phonon energy is calculated by \(\tilde{n}_{\lambda} \hbar \omega_{\lambda}\), with non-equilibrium phonon occupation number \(\tilde{n}_{\lambda}\). This non-equilibrium occupation number can be described as a linear response to the temperature gradient,

\[ \tilde{n}_{\lambda} \approx n_{\lambda} - \mathbf{v}_{\lambda}\tau_{\lambda}\frac{dn_{\lambda}}{dT} \nabla T \]

From Fourier's law, \(\mathbf{J} = \mathbf{\kappa} \nabla T\) and the phonon heat capacity \(c_{\lambda} = \hbar \omega_{\lambda}\frac{dn_{\lambda}}{dT}\), we can write the phonon heat conductivity as,

\[ \mathbf{\kappa}_{\alpha \beta} = \frac{1}{V}\sum_{\lambda} c_{\lambda}\mathbf{v}_{\lambda}^{\alpha}\mathbf{v}_{\lambda}^{\beta}\tau_{\lambda} \]

where \(\alpha\) and \(\beta\) are the Cartesian coordinates, and \(\tau_{\lambda}\) is the phonon lifetime. In Boltzmann transport theory, the phonon lifetime is estimated from the phonon-phonon scattering event, which can be calculated from perturbation theory.

In the harmonic approximation, phonon modes are described by eigenvalues \(\omega_s^2(\mathbf{q})\) and eigenvectors \(\boldsymbol{e}_s(\mathbf{q})\) of the dynamical matrix \(D(\mathbf{q})\). Due to the orthogonality of eigenvectors, the vibrational modes do not couple with each other. To study the thermal transport, one has to consider phonon-phonon interactions, which can be described using perturbation theory by treating third or higher-order force constants as the perturbation to phonons,

\[\begin{split} E \left(\{\mathbf{R}^0 + \Delta \mathbf{R}\}\right) \approx E\left(\{\mathbf{R}^0\}\right) & + \cancel{ \sum\limits_{I} \left.\frac{\partial E}{\partial \mathbf{R}_I}\right\vert_{\mathbf{R}^0} \Delta\mathbf{R}_{I} } \qquad (1)\\ & + \frac{1}{2} \sum\limits_{I,J} \left.\frac{\partial^2 E}{\partial \mathbf{R}_I\partial \mathbf{R}_J}\right\vert_{\mathbf{R}^0} \Delta\mathbf{R}_{I}\Delta\mathbf{R}_{J} \\ & + \frac{1}{6} \sum\limits_{I,J,K} \left.\frac{\partial^3 E}{\partial \mathbf{R}_I\partial \mathbf{R}_J\partial \mathbf{R}_K}\right\vert_{\mathbf{R}^0} \Delta\mathbf{R}_{I}\Delta\mathbf{R}_{J}\Delta\mathbf{R}_{K} \\ & + \mathcal{O}(\Delta\mathbf{R}^4) \end{split}\]

Here we go to the third-order derivative tensor of the potential energy surface, or 3rd-order force constants \(\mathbf{\Phi}_{IJK} = \frac{\partial^3 E}{\partial \mathbf{R}_I \partial\mathbf{R}_J \partial\mathbf{R}_K}\). In order to calculate the 3rd-order force constant, one can again use the finite displacements method,

\[\begin{split} \mathbf{\Phi}_{IJK} = & \left.\frac{\partial^3 E}{\partial \mathbf{R}_I\partial \mathbf{R}_J\partial \mathbf{R}_K}\right\vert_{\mathbf{R}^0} \\ = & - \left.\frac{\partial^2 }{\partial \mathbf{R}_I\partial \mathbf{R}_J} \mathbf{F}_K\right\vert_{\mathbf{R}^0} \\ \approx & - \frac{ \mathbf{F}_K(\mathbf{R}_I^0 + \varepsilon_I \mathbf{d}_I , \mathbf{R}_J^0 + \varepsilon_J \mathbf{d}_J)}{\varepsilon_I \varepsilon_J} \quad ,\qquad (2) \end{split}\]

which uses the forces \(\mathbf{F}_K\) of atom \(K\) at the displacements of atoms \(\mathbf{R}_I\) and \(\mathbf{R}_J\). Due to the displacements of the first atoms, the crystal symmetry is reduced to the point symmetry, and the number of the displaced supercells will increase. However, the third-order force constants also decay fast for non-ionic crystals, so a cutoff ion pair distance can be applied to reduce the computational cost.

With the 3rd-order force constants we can calculate the phonon lifetimes. The phonon lifetime \(\tau_\lambda\) of mode \(\lambda\) due to the phonon-phonon scattering is related to the imaginary part of the phonon self energy (\(\Sigma = \Delta + i\Gamma\)), by \(1/\tau_{\lambda} = 2\Gamma_{\lambda}\).

The phonon self energy is calculated by,

\[\begin{split} \Gamma_{\lambda} = & \frac{\hbar \pi}{16} \sum_{\lambda' \lambda''} |\mathbf{\Phi}_{\lambda \lambda' \lambda''} |^2 [(n_{\lambda'} + n_{\lambda''} + 1)\delta(\omega_{\lambda} - \omega_{\lambda'} - \omega_{\lambda''}) \\ + & 2(n_{\lambda'} - n_{\lambda''})\delta(\omega_{\lambda} - \omega_{\lambda'} + \omega_{\lambda''})] \end{split}\]

where \(n_{\lambda}\) is the phonon occupation number at equilibrium. The sum is over momentum conservation with \(\mathbf{q} + \mathbf{q}' + \mathbf{q}'' = \mathbf{G}\), and the delta function ensures the energy conservation. The three-phonon matrix is given by,

\[ \mathbf{\Phi_{\lambda \lambda' \lambda''}} = \sum_{IJK} \frac{ \mathbf{e}^{I}_{\lambda} \mathbf{e}^{J}_{\lambda'} \mathbf{e}^{K}_{\lambda''}}{\sqrt{M_I M_J M_K} \sqrt{\omega_\lambda \omega_{\lambda'} \omega_{\lambda''}} } \mathbf{\Phi}_{IJK}e^{(i\mathbf{q}\cdot\mathbf{R}_I + \mathbf{q}'\cdot\mathbf{R}_J + i\mathbf{q}''\cdot\mathbf{R}_K)} \]

where \(M_I\) is the mass of atom \(I\), \(\mathbf{e}_{\lambda}^{I}\) is the eigenvector of phonon mode \(\lambda\) and atom \(I\).