Phonon Boltzmann Transport Equation

Updated 28 December 2025

Phonon Boltzmann Transport Equation is a kinetic framework describing heat flow through lattice vibrations in solids, encompassing ballistic, diffusive, and hydrodynamic regimes.
Its formulation leverages linearization and first-principles inputs to resolve mode-dependent phenomena such as thermal boundary resistance, second sound, and nonlocal effects.
Advanced numerical techniques—including discrete ordinates, variational methods, and machine learning surrogates—enhance simulation accuracy across multiscale, nonlocal heat transport conditions.

The phonon Boltzmann Transport Equation (BTE) is a fundamental kinetic framework for modeling heat transport by lattice vibrations in solid-state systems. The BTE provides a first-principles description of phonon evolution in position, momentum, energy, and polarization spaces, bridging ballistic, hydrodynamic, and diffusive regimes. It enables prediction and analysis of various phenomena such as size effects, thermal boundary resistance, second sound, and thermal conductivity in nanostructured and bulk materials. The equation is indispensable for both theoretical studies and computational simulations involving multiscale and nonlocal heat flow.

1. Mathematical Formulation and Structure

The general form of the linearized frequency- and direction-resolved phonon BTE under the relaxation-time approximation (RTA) is

$\frac{\partial f(\mathbf{r},\omega,\mathbf{s},t)}{\partial t} + \mathbf{v}_g(\omega) \cdot \nabla_{\mathbf{r}} f = -\frac{f - f_0(T(\mathbf{r},t))}{\tau(\omega)},$

where $f$ is the non-equilibrium phonon distribution, $f_0$ is the Bose-Einstein equilibrium at temperature $T$ , $\mathbf{v}_g(\omega)$ is group velocity, and $\tau(\omega)$ the relaxation time. After linearization and in energy form, key variables are the energy density $e$ and the heat flux $\mathbf{q}$ via

$e(\mathbf{r},\omega,\mathbf{s},t) = \hbar \omega D(\omega) [f - f_0(T_{ref})],\quad \mathbf{q}(\mathbf{r}) = \int \mathbf{v}_g e\,d\omega\,d\Omega,$

with density of states $D(\omega)$ .

Spectral and mode-resolved versions of the BTE treat each branch and polarization independently, allowing $p$ and $\lambda$ indices for polarization/wavevector. Full first-principles implementations sample the Brillouin zone and include ab initio collision operators (Jain et al., 29 May 2025, 1804.01729). Non-gray approaches utilize frequency- and polarization-dependent $v$ , $\tau$ , and $C$ (Zhang et al., 2018, Hu et al., 2021, Hu et al., 2023).

2. Physical Regimes and Limiting Behaviors

2.1 Ballistic versus Diffusive Transport

In the ballistic limit ( $L\ll\Lambda$ , $\tau\to\infty$ ), the collision term vanishes, and phonons propagate at their group velocity, leading to wave-like nonequilibrium evolution (Maassen et al., 2015). The diffusive limit ( $L\gg\Lambda$ , frequent collisions) recovers Fourier’s law with conductivity $\kappa=C_V v \Lambda/2$ in the gray model, and more generally

$\kappa = \frac{1}{V} \sum_{\lambda} \hbar \omega_\lambda v_\lambda^2 \tau_\lambda \frac{\partial f^0_\lambda}{\partial T}$

for mode-resolved implementations (Jain et al., 29 May 2025, 1804.01729).

2.2 Hydrodynamic and Collective Effects

When momentum-conserving Normal processes dominate over Umklapp/resistive processes, phonon flow becomes hydrodynamic. The appropriate model is the Callaway dual-relaxation formulation: $\frac{\partial f}{\partial t} + \mathbf{v}_g \cdot \nabla f = -\frac{f-f^N_{eq}}{\tau_N} - \frac{f-f^R_{eq}}{\tau_R}$ yielding Guyer-Krumhansl equations in the hydrodynamic regime, supporting second sound and Poiseuille flow (Guo et al., 2020, Qian et al., 26 Oct 2024).

2.3 Transition Criteria

Transport regime is dictated by the Knudsen number ( $\mathrm{Kn}=\Lambda/L$ ), with boundary/surface scattering, interface transmission (Kapitza resistance), and spatial scale $L$ tuning the crossover (Péraud et al., 2015, Ordonez-Miranda et al., 2015, Chiloyan et al., 2017, Romano, 2020).

3. Numerical Solution Methodologies

3.1 Discrete Ordinates, FVM, and Synthetic Iteration

Deterministic solvers (e.g., DOM + FVM) discretize frequency, direction, space (and optionally time), solving for $e(\omega, p, s, r, t)$ at each point (Hu et al., 2023). Synthetic iterative schemes couple microscopic BTE solvers with macroscopic diffusion steps, achieving rapid convergence in near-diffusive systems while retaining ballistic accuracy (Zhang et al., 2018, Zhang et al., 2022).

Band discretization can be optimized using mean free path domain splitting and Gauss-Legendre quadrature, reducing band number to $N<10$ for 1% accuracy (Hu et al., 2021).

3.2 Variational and Green's Function Methods

Variational techniques exploit trial (often Fourier-like) solutions to minimize PDE residuals, yielding closed-form expressions for multidimensional decay rates and conductivity including suppression functions capturing nonlocality and material dependence (Chiloyan et al., 2015, Chiloyan et al., 2016). Full Green’s function methodologies invert the BTE with the full scattering matrix, enabling modeling of ultrafast experiments and nonlocal phenomena; these are essential for hydrodynamic/collective effects (Chiloyan et al., 2017, Qian et al., 26 Oct 2024, Romano, 2020).

3.3 Machine Learning and Surrogates

Physics-informed neural networks (PINNs and MC-PINNs) encode the BTE and boundary conditions as neural-network loss functions in high-dimensional variable spaces (space, frequency, direction, time, polarization), bypassing grid-based discretization and directly learning device-parametric solutions (Li et al., 2021, Li et al., 2022, Lin et al., 20 Aug 2024). Hybrid approaches combine fast, interpretable low-fidelity PDE solvers with deep surrogates and active learning, reducing the required simulation data for fast conductivity prediction and design (Varagnolo et al., 25 Nov 2025).

4. Boundary Conditions and Interface Modeling

Boundary conditions in BTE solvers range from thermalizing (isothermal), adiabatic (specular or diffuse reflection), periodic, and partial transmission/reflection (for interfaces with thermal boundary resistance) (Péraud et al., 2015, Hu et al., 2023, Gamba et al., 2020). Thermal interface resistance is modeled via frequency-dependent reflection/transmission coefficients, which can be inferred via PDE-constrained inverse problems using adjoint methods and stochastic optimization (Gamba et al., 2020). Asymptotic analysis yields matched jump-type boundary conditions and universal boundary-layer corrections quantified in terms of material and interface properties (Péraud et al., 2015).

5. First-Principles and Ab Initio Implementation

Contemporary computational schemes input phonon mode frequencies, group velocities, and collision matrices from density functional theory (DFT) or direct supercell calculations, treating both three-phonon and four-phonon scattering, multi-channel transport (Wigner/BTE), and temperature-dependent interatomic force constants (Jain et al., 29 May 2025, 1804.01729). Iterative linear algebra and variational minimization schemes (e.g., ShengBTE, almaBTE, fourPhonon) are deployed for the full linearized BTE (Jain et al., 29 May 2025).

6. Experimental Interpretation and Inverse Approaches

The spectral BTE enables quantitative analysis of transient thermal grating (TTG), time-domain thermoreflectance (TDTR), and other pump-probe thermal experiments. Variational and inversion methods relate measured effective thermal conductivity and decay rates to underlying mean free path distributions, enabling extraction of the bulk and size-dependent phonon thermal conductivity accumulation function (Chiloyan et al., 2015, Chiloyan et al., 2016).

Inverse problems also arise in reconstructing interface transmission coefficients from surface temperature data, utilizing adjoint-based Fréchet derivatives and stochastic gradient descent with rigorous regularity guarantees (Gamba et al., 2020).

7. Practical Applications and Toolchains

Deterministic BTE solvers such as GiftBTE provide scalable, parameter-free simulation capability for nanostructure thermal conductivity, laser heating, and transistor self-heating. Synthetic iterative schemes and advanced spatial/angle discretization achieve computational efficiency and accuracy suitable for practical device analysis (Hu et al., 2023). PINN and MC-PINN frameworks are increasingly adopted for device-level design where high dimensionality and parametric flexibility are required (Li et al., 2021, Lin et al., 20 Aug 2024).

Physics-enhanced deep surrogates incorporate physical inductive biases with data-driven corrections and active learning, achieving fast and robust conductivity prediction for design optimization across ballistic and diffusive regimes (Varagnolo et al., 25 Nov 2025).

The phonon Boltzmann Transport Equation is the definitive mesoscale framework for the theory and simulation of heat flow in materials, enabling physically rigorous prediction from first principles. Its modern solution methodologies—both deterministic and learning-based—are at the forefront of material design, experiment interpretation, and computational prediction for micro- and nanoscale thermal engineering.