Phonon BTE: Models, Methods & Applications

• Phonon BTE is a fundamental kinetic model describing non-equilibrium lattice dynamics and heat transport in crystalline solids at micro- and nanoscales.
• Modern solution strategies blend deterministic solvers, variational methods, and neural surrogates to enable efficient multiscale simulation and device optimization.
• The equation incorporates spectral, polarization, and scattering details to capture ballistic, hydrodynamic, and diffusive thermal transport regimes.

The phonon Boltzmann Transport Equation (BTE) is the fundamental kinetic equation governing non-equilibrium lattice dynamics and heat transport in crystalline solids, especially at micro- and nanoscales where standard Fourier diffusion fails. It models the space-time evolution of the phonon distribution function under the combined effects of propagation (drift), phonon-phonon and phonon-defect scattering, and boundary or source interactions. The equation incorporates full spectral, polarization, and geometric effects, and is essential for accurate analysis and design of thermal transport in semiconductors, 2D materials, alloys, and nanostructures. Modern solution strategies span direct deterministic solvers, synthetic iterative acceleration, surrogate neural models, and variational and Green’s function approaches. The BTE’s flexibility enables both forward and inverse, as well as differentiable, workflows for multiscale thermal simulation and device optimization.

1. Mathematical Formulation and Physical Interpretation

The general (time-dependent, non-gray) phonon BTE describes the evolution of the phonon occupation function $f(\mathbf{r}, \mathbf{k}, t)$:

$$\frac{\partial f(\mathbf{r}, \mathbf{k}, t)}{\partial t} + \mathbf{v}_g(\mathbf{k}) \cdot \nabla_{\mathbf{r}} f(\mathbf{r}, \mathbf{k}, t) = \left(\frac{\partial f}{\partial t}\right)_{\text{coll}} + S(\mathbf{r}, \mathbf{k}, t)$$

where $\mathbf{v}_g$ is the group velocity, and $(\partial f/\partial t)_{\text{coll}}$ is the phonon–phonon (and possibly phonon–impurity or phonon–electron) collision integral. $S$ represents external sources.

Linearization about a reference temperature $T_0$ yields the deviational BTE, both in single-mode relaxation time (RTA) and full-scattering-matrix forms:

• RTA (Single-mode):

$$\mathbf{v}_g(\mathbf{k})\cdot \nabla_{\mathbf{r}} f(\mathbf{r}, \mathbf{k}) = \frac{f^0(\mathbf{r},\mathbf{k}) - f(\mathbf{r},\mathbf{k})}{\tau(\mathbf{k})}$$

where $f^0$ is the Bose–Einstein equilibrium distribution and $\tau(\mathbf{k})$ is the mode-dependent relaxation time.

• Full matrix:

$$\mathbf{v}_n \cdot \nabla_{\mathbf{r}} g_n = \sum_{j} S_{nj}[C_j \Delta T(\mathbf{r}, t) - g_j(\mathbf{r}, t)]$$

$$\frac{\partial f(\mathbf{r}, \mathbf{k}, t)}{\partial t} + \mathbf{v}_g(\mathbf{k}) \cdot \nabla_{\mathbf{r}} f(\mathbf{r}, \mathbf{k}, t) = \left(\frac{\partial f}{\partial t}\right)_{\text{coll}} + S(\mathbf{r}, \mathbf{k}, t)$$0 is the mode-resolved deviational energy, $$\frac{\partial f(\mathbf{r}, \mathbf{k}, t)}{\partial t} + \mathbf{v}_g(\mathbf{k}) \cdot \nabla_{\mathbf{r}} f(\mathbf{r}, \mathbf{k}, t) = \left(\frac{\partial f}{\partial t}\right)_{\text{coll}} + S(\mathbf{r}, \mathbf{k}, t)$$1 is the mode heat capacity, and $$\frac{\partial f(\mathbf{r}, \mathbf{k}, t)}{\partial t} + \mathbf{v}_g(\mathbf{k}) \cdot \nabla_{\mathbf{r}} f(\mathbf{r}, \mathbf{k}, t) = \left(\frac{\partial f}{\partial t}\right)_{\text{coll}} + S(\mathbf{r}, \mathbf{k}, t)$$2 is the full linearized scattering kernel.

Physical observables are obtained by summing over all phonon modes:

• Heat flux: $$\frac{\partial f(\mathbf{r}, \mathbf{k}, t)}{\partial t} + \mathbf{v}_g(\mathbf{k}) \cdot \nabla_{\mathbf{r}} f(\mathbf{r}, \mathbf{k}, t) = \left(\frac{\partial f}{\partial t}\right)_{\text{coll}} + S(\mathbf{r}, \mathbf{k}, t)$$3
• Thermal conductivity: Derived from the linear response, with nontrivial accumulated and spectral forms in non‐gray, anisotropic, or nanostructured systems.

The BTE thus provides a direct connection from first-principles phonon properties to macroscopic thermal transport phenomena, rigorously bridging the ballistic, hydrodynamic, and diffusive regimes (Chiloyan et al., 2017, Hu et al., 2021, Chiloyan et al., 2016).

2. Collision Operators and Scattering Physics

The collision integral is the crux of phonon BTE modeling, encoding all thermalization processes:

• Full three-phonon scattering matrix $$\frac{\partial f(\mathbf{r}, \mathbf{k}, t)}{\partial t} + \mathbf{v}_g(\mathbf{k}) \cdot \nabla_{\mathbf{r}} f(\mathbf{r}, \mathbf{k}, t) = \left(\frac{\partial f}{\partial t}\right)_{\text{coll}} + S(\mathbf{r}, \mathbf{k}, t)$$4 is obtained ab initio from third-order force constants, coupling all modes via energy- and momentum-conserving (normal and Umklapp) processes. This enables hydrodynamic phenomena (second sound), spectral renormalization (relaxons), and energy conservation beyond RTA (Chiloyan et al., 2017, Romano, 2020).
• Relaxation-time approximation (RTA): Diagonalizes $$\frac{\partial f(\mathbf{r}, \mathbf{k}, t)}{\partial t} + \mathbf{v}_g(\mathbf{k}) \cdot \nabla_{\mathbf{r}} f(\mathbf{r}, \mathbf{k}, t) = \left(\frac{\partial f}{\partial t}\right)_{\text{coll}} + S(\mathbf{r}, \mathbf{k}, t)$$5 as $$\frac{\partial f(\mathbf{r}, \mathbf{k}, t)}{\partial t} + \mathbf{v}_g(\mathbf{k}) \cdot \nabla_{\mathbf{r}} f(\mathbf{r}, \mathbf{k}, t) = \left(\frac{\partial f}{\partial t}\right)_{\text{coll}} + S(\mathbf{r}, \mathbf{k}, t)$$6, which neglects momentum-conserving normal scattering and leads to inaccurate predictions in high‐Debye-temperature or 2D materials (Romano, 2020).
• Beyond Fermi’s Golden Rule (FGR): In circumstances where the linewidth from anharmonic self-energies is non-negligible compared to the mode splitting, as in strongly anharmonic or low-dimensional crystals, the standard delta-function kernel fails. Collisional broadening and generalized kinetic frameworks (derived from Kadanoff-Baym equations) yield improved, smearing-independent, and physically correct results (Lucente et al., 17 Mar 2026).
• Electron–phonon coupling: For coupled energy and charge transport, the BTE must be simultaneously solved for both systems, incorporating first-principles electron–phonon coupling. Mutual drag effects become critical for thermoelectric transport at low temperature (Protik et al., 2019).

The physical content and computational cost are set by the fidelity of the collision operator. Full-matrix treatments are mandatory for accurate second sound, hydrodynamics, and 2D/heterostructure phenomena (Chiloyan et al., 2017, Romano, 2020, Lucente et al., 17 Mar 2026).

3. Numerical Solution Methodologies

A variety of deterministic and surrogate computational strategies have been developed to solve the high-dimensional BTE:

• Discrete Ordinates and Band Discretization: The BTE is discretized over frequency/polarization (bands), angles, and spatial cells. Optimized mean-free-path (MFP) discretization schemes, such as the two-subdomain Gauss–Legendre quadrature, achieve high accuracy for non-gray BTE at minimal cost (≤10 bands for <1% error), enabling large-scale, realistic 3D simulations (Hu et al., 2021).
• Implicit Kinetic Schemes with Macroscopic Correction: Two-level solvers alternate microscopic kinetic sweeps (solving the band- and angular-resolved BTE) with macroscopic corrections (diffusion-like equations). Synthetic acceleration strategies combine a tightly coupled iteration between the BTE and a moment-based PDE (e.g., synthetic diffusion equation for temperature) for rapid convergence across all Knudsen numbers (Zhang et al., 2018, Zhang et al., 2018, Zhang et al., 2022).
• Variational and Green’s Function Approaches: For canonical geometries and experiments (e.g., transient thermal grating, TTG), variational methods yield closed-form solutions for decay rates and effective conductivities, revealing the material- and geometry-dependent suppression function for thermal transport (Chiloyan et al., 2016, Chiloyan et al., 2015). Analytical Green’s function formalisms—both in RTA and Callaway (two-τ) models—are used for multidimensional, time-dependent, and hydrodynamic regimes (Qian et al., 2024, Chiloyan et al., 2017).
• Surrogate Modeling and Differentiable Programming: Neural-augmented surrogates such as PEDS couple a low-fidelity differentiable Fourier solver with a neural network generator to interpolate between ballistic and diffusive predictions, enabling efficient inverse design over material/geometric parameter spaces (Varagnolo et al., 25 Nov 2025). Physics-informed neural networks (PINNs), including Monte Carlo sampling variants, encode the BTE and boundary conditions directly into neural losses, yielding mesh-free, parameterized solvers capable of handling multi-regime or large–$$\frac{\partial f(\mathbf{r}, \mathbf{k}, t)}{\partial t} + \mathbf{v}_g(\mathbf{k}) \cdot \nabla_{\mathbf{r}} f(\mathbf{r}, \mathbf{k}, t) = \left(\frac{\partial f}{\partial t}\right)_{\text{coll}} + S(\mathbf{r}, \mathbf{k}, t)$$7 problems (Li et al., 2021, Li et al., 2022, Lin et al., 2024).
• High-Performance Implementations: GPU-accelerated differentiable solvers (e.g., JAX-BTE) capitalize on automatic differentiation frameworks to make forward and inverse, end-to-end optimization workflows practical for complex device design and parameter fitting (Shang et al., 31 Mar 2025).

4. Boundary Conditions, Geometry, and Device Modeling

The phonon BTE is formulated with diverse boundary conditions and geometries relevant to real devices:

• Boundary conditions: Fully thermalizing (isothermal, black), specular (mirror-like), diffuse (randomizing), Robin-type (combining temperature and flux), and periodic. Accurate treatment of contact resistances and ballistic slip is essential in the ballistic and quasi-ballistic regimes (Maassen et al., 2015, Hu et al., 2021, Zhang et al., 2018).
• Nano/microstructure modeling: Real devices (FinFETs, transistors, porous media, superlattices) and experimental metrologies (TTG, TDTR) require 3D BTE solutions with full phonon dispersion, polarization, and first-principles scattering data (Shang et al., 31 Mar 2025, Hu et al., 2023, Carrete et al., 2017).
• Inverse and Data-Augmented Workflows: Differentiable programming approaches enable inverse optimization—extracting spatially resolved, geometry-dependent, or spectral transport parameters directly from measured or desired device responses (Shang et al., 31 Mar 2025, Varagnolo et al., 25 Nov 2025).

The rigorous coupling to first-principles phonon properties, along with robust band and angular discretization, is crucial for both forward prediction and device-targeted inverse design (Hu et al., 2021, Hu et al., 2023).

5. Application Regimes and Physical Phenomena

The phonon BTE accurately captures a spectrum of physical regimes and transport phenomena not accessible to macroscopic theories:

• Ballistic to Diffusive Crossover: The BTE reproduces the transition from boundary-limited ballistic transport (dominant mean free paths $$\frac{\partial f(\mathbf{r}, \mathbf{k}, t)}{\partial t} + \mathbf{v}_g(\mathbf{k}) \cdot \nabla_{\mathbf{r}} f(\mathbf{r}, \mathbf{k}, t) = \left(\frac{\partial f}{\partial t}\right)_{\text{coll}} + S(\mathbf{r}, \mathbf{k}, t)$$8) to scattering-dominated diffusion ($$\frac{\partial f(\mathbf{r}, \mathbf{k}, t)}{\partial t} + \mathbf{v}_g(\mathbf{k}) \cdot \nabla_{\mathbf{r}} f(\mathbf{r}, \mathbf{k}, t) = \left(\frac{\partial f}{\partial t}\right)_{\text{coll}} + S(\mathbf{r}, \mathbf{k}, t)$$9), including non-local temperature jumps and finite propagation speed (Maassen et al., 2015, Zhang et al., 2018).
• Hydrodynamic and Second Sound: Only full-matrix BTE or Callaway models recover hydrodynamic phenomena in phonon transport (normal processes dominate, giving rise to well-defined heat waves—second sound) (Chiloyan et al., 2017, Qian et al., 2024).
• Quasi-Equilibrium Breakdown: In strong non-equilibrium, especially after ultrafast excitation, the notion of a branch-resolved or local Bose–Einstein temperature fails (as shown by power-law decay tails and non-thermal subpopulations), demanding full kinetic resolution (Ono, 2017).
• Nanostructured and Patterned Media: Non-diffusive suppression of heat conduction due to feature-scale comparable to phonon MFPs; large-scale BTE simulation with appropriate MFP binning and realistic phonon spectra is essential for quantitative analysis (Hu et al., 2021, Hu et al., 2023).
• Coupled Electron–Phonon Transport: For comprehensive thermoelectric and low-temperature transport, the coupled solution of phonon and electron BTEs resolves mutual drag and reciprocal thermal phenomena not visible in uncoupled approaches (Protik et al., 2019).

6. Limitations, Open Challenges, and Future Directions

Despite these advances, several critical challenges remain unresolved in the phonon BTE paradigm:

• Accurate Collision Operator Construction: Standard FGR-based collisional kernels are inadequate in the presence of strong anharmonicity or for systems with flexural modes, requiring rigorous approaches derived from quantum kinetic theory and Kadanoff–Baym equations to remedy unphysical predictions and non-convergent conductivities (Lucente et al., 17 Mar 2026).
• Computational Cost: The high-dimensionality of the non-gray, full-scattering BTE renders brute-force deterministic solution impractical for “realistic” device geometries without advanced discretization, acceleration, and surrogate strategies (Zhang et al., 2018, Hu et al., 2021, Hu et al., 2023, Varagnolo et al., 25 Nov 2025).
• Transient and Strong Nonequilibrium Regimes: Most current deterministic and neural approaches are formulated for steady-state, linearized (small–$\mathbf{v}_g$0) conditions. Extending to highly transient, nonlinear, or non-equilibrium branches (e.g., after ultrafast excitation) is an active area, with emerging work in recurrent neural architectures and non-quasiparticle transport (Li et al., 2022, Ono, 2017).
• Coupling to Other Excitations: Electron–phonon and magnon–phonon coupling, as well as frequency-dependent boundary conditions, demand generalized frameworks for coupled transport and scattering (Protik et al., 2019).
• Experimental Model Extraction: Variational BTE methods are making quantitative inversion from experiment to full phonon MFP spectra more robust and systematic, but aspects of experimental validation and complex device-feature mapping require further development (Chiloyan et al., 2016, Chiloyan et al., 2015).

A plausible implication is that future progress will hinge on integrating quantum-derived collision operators, multilevel surrogate modeling, and scalable differentiable implementations, alongside deeper connection to experimental metrology and device informatics (Shang et al., 31 Mar 2025, Lucente et al., 17 Mar 2026, Varagnolo et al., 25 Nov 2025).

7. Software Frameworks and Practical Implementations

State-of-the-art open-source and proprietary solvers for the phonon BTE include:

Framework	Features	Reference
JAX-BTE	Differentiable, GPU-accelerated, forward/inverse sim.	(Shang et al., 31 Mar 2025)
GiftBTE	Deterministic, non-gray (multi-band), 3D, synthetic	(Hu et al., 2023)
almaBTE	Space-time, collision-matrix, MC & deterministic modes	(Carrete et al., 2017)
PINN/MC-PINN	Neural-BTE, parametric, mesh-free, multi-regime	(Li et al., 2021, Lin et al., 2024)
Synthetic DOM	Coupled mesoscopic/macroscopic, rapid convergence	(Zhang et al., 2018, Zhang et al., 2022)
PEDS	Differentiable surrogates, uncertainty, inverse design	(Varagnolo et al., 25 Nov 2025)

Most modern codes interface directly with first-principles phonon inputs (DFT, DFPT) and provide end-to-end data pipelines for experimental validation, device property prediction, and materials design.

---