In microscale and nanoscale heat conduction, size effects mainly involve two key length scales: phonon wavelength and mean free path (MFP).
The First-Principles Phonon Boltzmann Transport Equation (PBTE) method is essentially a combination of the following three techniques:
First-principles (ab initio) methods refer to techniques that solve the electronic structure and atomic interactions of materials within the framework of quantum mechanics, using minimal empirical parameters. They are primarily based on numerical solutions to the Schrödinger equation, yielding properties such as ground-state electron density, band structure, total energy, and other physical quantities. In solid-state physics and materials science, first-principles methods are typically based on Density Functional Theory (DFT). As the full many-body Schrödinger equation is difficult to solve directly, DFT provides an efficient way to handle many-body problems via the electron density and is widely applied to metals, semiconductors, insulators, molecules, and surfaces.
The core concepts of DFT originate from the Hohenberg-Kohn theorem and the Kohn-Sham equation:
Where:
In practice, one must choose an exchange-correlation functional (e.g., LDA, GGA, or hybrid functionals) and adopt a method to treat real and reciprocal space (e.g., pseudopotentials, PAW, APW+lo, etc.).
A crucial step in thermal conductivity calculations is determining the interatomic force constants (IFCs) of the material. These constants are used to construct the dynamical matrix and calculate lattice dynamical properties such as phonon dispersion and phonon lifetimes. IFCs can be obtained from first-principles calculations using two main methods:
In the DFPT framework, small periodic perturbations are applied to the periodic structure, and the linearized Kohn-Sham equations are solved to obtain the first-, second-, or third-order responses of the potential energy with respect to atomic displacements.
In this method, small displacements (e.g., 0.01 Å) are applied to atoms in a supercell, and the resulting changes in energy or forces are used to compute second-order and higher-order IFCs:
Both methods rely on DFT-level force or energy calculations. Since phonon properties are sensitive to calculation accuracy, careful choices of exchange-correlation functional, plane-wave cutoff energy, and k-point mesh are necessary to ensure accurate IFCs.
In harmonic lattice dynamics, the second-order IFCs are used to obtain the phonon dispersion relation $\omega_\lambda(\mathbf{q})$. Once dispersion is known, specific heat of each phonon mode $\lambda$ can be computed. The phonon group velocity is defined as the gradient of the frequency with respect to wave vector: $\mathbf{v}\lambda = \nabla\mathbf{q} \omega_\lambda$ The relaxation time $\tau_\lambda$ is obtained through anharmonic lattice dynamics, which involves both second- and higher-order IFCs.
For a periodic crystal with small displacements around equilibrium, the total potential energy $U$ can be expanded as a Taylor series:
\[U = U_0 + \frac{1}{2!} \sum_{ij}\sum_{\alpha\beta} \Phi_{ij}^{\alpha\beta} u_i^\alpha u_j^\beta + \frac{1}{3!} \sum_{ijk}\sum_{\alpha\beta\gamma} \Psi_{ijk}^{\alpha\beta\gamma} u_i^\alpha u_j^\beta u_k^\gamma + \mathcal{O}(u^4)\]Where:
Due to force equilibrium, first-order terms are absent. Neglecting higher-order terms leads to the harmonic approximation.
If atom $i$ belongs to the $b$-th atom in the $l$-th unit cell and $j$ is in the $l’$-th cell, the equation of motion becomes:
\[m_b\frac{d^2 u_{lb}^\alpha (t)}{d t^2} = -\sum_{l'b',\beta} \Phi_{lb,l'b'}^{\alpha\beta} u_{l'b'}^{\beta} (t)\]Assuming plane wave solutions:
\[u_{lb}^\alpha (t) = \frac{1}{\sqrt{m_b}}\Lambda_\lambda e_{b,\lambda}^\alpha e^{i(\mathbf q\cdot\mathbf R_l - \omega_\lambda t)}\]Substituting gives the eigenvalue equation:
\[\omega_\lambda^2 \mathbf e _{b,\lambda} = \mathbf D_{bb'}^{\alpha\beta}(\textbf q) \textbf e_{b,\lambda}\]Where the dynamical matrix is:
\[\mathbf D_{bb'}^{\alpha\beta}(\textbf q) =\frac{1}{\sqrt{m_bm_{b'}}}\sum_{l'}\Phi_{0b,l'b'}^{\alpha\beta} e^{i\mathbf q\cdot(\mathbf R_{l'} - \mathbf R_0)}\]Solving gives the phonon dispersion $\omega_\lambda(\mathbf{q})$ and corresponding eigenvectors ${e_{b,\lambda}^\alpha}$.
In real materials, phonons have finite lifetimes due to various scattering mechanisms. Within the BTE framework, these are expressed as scattering rates $\Gamma_\lambda$ or equivalently relaxation times $\tau_\lambda$.
In nonmetals, the dominant scattering is phonon-phonon interactions, i.e., anharmonic lattice effects. Quantum mechanically, the crystal Hamiltonian is divided into harmonic and anharmonic parts. Third-order scattering (three-phonon processes) includes combination and splitting processes. Energy and momentum conservation must be satisfied (including Umklapp processes):
\[\Gamma_{\lambda\lambda'\lambda''}^{\pm} = \frac{\hbar\pi}{4}\begin{Bmatrix}n_{\lambda'}-n_{\lambda''} \\ n_{\lambda'} + n_{\lambda''}+1 \end{Bmatrix} \frac{\delta(\omega_\lambda\pm\omega_{\lambda'} - \omega_{\lambda''})}{\omega_{\lambda}\omega_{\lambda'}\omega_{\lambda''}}|V_{\lambda\lambda'\lambda''}^{\pm}|^2\Delta_{\mathbf{G}, \,\mathbf{q}\pm \mathbf{q}' - \mathbf{q}''}\]Caused by mass or bonding perturbations from impurities (e.g., isotopes, defects):
\[\Gamma_{\lambda\lambda'} = \frac{\pi \omega_\lambda^2}{2}\delta(\omega_\lambda - \omega_{\lambda'})\sum_bg(b) |\mathbf e_{b,\lambda}^*\cdot \mathbf e_{b,\lambda'}|^2\]Where $g(b) = \sum_s f_s(b)(1 - \frac{m_{b,s}}{\overline{m}_b})^2$ and $\overline{m}b = \sum_s f_s(b) m{b,s}$.
According to Fourier’s law, the thermal conductivity $\kappa$ characterizes the ability of a material to conduct heat:
\[\mathbf{J} = -\kappa \nabla T\]Where $\mathbf{J}$ is the heat flux vector and $\nabla T$ is the temperature gradient. In anisotropic materials, $\kappa$ is a tensor.
To predict $\kappa$ from phonon properties, one must use the BTE to describe the non-equilibrium phonon distribution under a small temperature gradient. Assuming steady-state and small gradients, the linearized BTE is:
\[-\mathbf{v}_\lambda \nabla T \frac{\partial n_\lambda^0}{\partial T} = \frac{n_\lambda'}{\tau_\lambda}\]Where:
The Bose–Einstein distribution:
\[n_\lambda^0 = \frac{1}{\exp\left(\frac{\hbar \omega_\lambda}{k_B T}\right) - 1}\]The heat flux contributed by phonons:
\[\mathbf{J} = \frac{1}{V} \sum_\lambda \hbar\omega_\lambda \mathbf{v}_\lambda n_\lambda\]Where $V = N_0 \cdot \Omega$ is the total volume ($\Omega$ = unit cell volume, $N_0$ = number of $\mathbf{q}$-points). Comparing with Fourier’s law gives the thermal conductivity tensor:
\[\kappa^{\alpha\beta} = \frac{1}{V}\sum_\lambda \hbar\omega_\lambda \frac{\partial n_\lambda^0}{\partial T} v_\lambda^\alpha v_\lambda^\beta \tau_\lambda = \sum_\lambda c_\lambda v_\lambda^\alpha v_\lambda^\beta \tau_\lambda\]Where $c_\lambda = \frac{\hbar\omega_\lambda}{V}\frac{\partial n_\lambda^0}{\partial T}$ is the mode heat capacity. For isotropic systems:
\[\kappa = \frac{1}{3V}\sum_\lambda \hbar\omega_\lambda \frac{\partial n_\lambda^0}{\partial T} |\mathbf{v}_\lambda|^2 \tau_\lambda\]The SMRTA assumes that all phonon modes except mode $\lambda$ remain in equilibrium:
\[\begin{cases} n_\lambda = n_\lambda^0 + n_\lambda' \\ n_{\lambda'} = n_{\lambda'}^0 \\ n_{\lambda''} = n_{\lambda''}^0 \end{cases}\]The SMRTA relaxation time becomes:
\[\frac{1}{\tau_\lambda^0} = \sum_{\lambda'\lambda''}^+ \Gamma_{\lambda\lambda'\lambda''}^+ + \sum_{\lambda'\lambda''}^- \frac{1}{2} \Gamma_{\lambda\lambda'\lambda''}^- + \sum_{\lambda'}\Gamma_{\lambda\lambda'}\]Where superscript 0 indicates this is the zero-order approximation.
To overcome SMRTA limitations, the Full Iterative Method solves the BTE self-consistently:
\[\begin{cases} n_\lambda = n_\lambda^0 + n_\lambda' \\ n_{\lambda'} = n_{\lambda'}^0 + n_{\lambda'}' \\ n_{\lambda''} = n_{\lambda''}^0 + n_{\lambda''}' \end{cases}\]Initial Guess: Use $\tau_\lambda^0$ from SMRTA.
Self-Consistent Iteration: Form a coupled system of equations for all $n_\lambda’$ and iterate until convergence.
The relaxation time becomes:
\[\tau_\lambda = \tau_\lambda^0 + \tau_\lambda^0 \Delta_\lambda\] \[\Delta_\lambda = \sum_{\lambda'\lambda''}^+ \Gamma_{\lambda\lambda'\lambda''}^+ (\xi_{\lambda\lambda''}\tau_{\lambda''} - \xi_{\lambda\lambda'}\tau_{\lambda'}) \\ + \sum_{\lambda'\lambda''}^- \frac{1}{2} \Gamma_{\lambda\lambda'\lambda''}^-(\xi_{\lambda\lambda''}\tau_{\lambda''} + \xi_{\lambda\lambda'}\tau_{\lambda'}) \\+ \sum_{\lambda'}\Gamma_{\lambda\lambda'}\xi_{\lambda\lambda'}\tau_{\lambda'}\]This method includes more detailed interactions and better captures phonon redistribution caused by Normal processes, especially in high thermal conductivity or low-temperature systems. However, it is also more computationally intensive.
The general workflow for calculating phonon thermal conductivity includes:
IFCs are the basis for lattice dynamics and thermal transport. They can be obtained via:
Two main extraction methods:
IFCs theoretically exist between all atomic pairs. In practice, a cutoff radius is applied. The impact of different cutoffs should be tested.
Also, due to numerical noise and truncation, computed IFCs may violate translational invariance or crystal symmetry. These must be corrected to ensure accurate thermal conductivity.
With second- and third-order IFCs, combine anharmonic lattice dynamics and BTE using:
This method requires no fitting parameters, only the initial atomic structure. It is highly predictive and widely validated against experiments.
Despite its strength, this method is still sensitive to:
Nevertheless, PBTE is considered one of the most reliable methods for predicting lattice thermal conductivity. It not only yields total κ, but also mode-resolved contributions, enabling studies of interfacial thermal conductance and nanoscale transport.
Several open-source packages support this workflow:
These integrate well with major first-principles codes like:
They enable an end-to-end automated workflow from structure optimization to IFC extraction and thermal conductivity prediction.