Semi-implicit Lax-Wendroff kinetic scheme for electron-phonon coupling

A semi-implicit Lax-Wendroff scheme is developed for electron-phonon coupling process in metals based on the two-temperature kinetic equations. The core of this method is to integrate the evolution information of physical equations into the numerical modeling process, which leads to that the time step or cell size is not limited by the relaxation time and mean free path. Specifically, the finite difference method is used to solve the kinetic model again when reconstructing the interfacial distribution function, through which the particle migration, scattering and electron-phonon coupling processes are coupled together within a single time step. Numerical tests demonstrate that this method could efficiently capture electron-phonon coupling or heat conduction processes from the ballistic to diffusive regimes. It provides a new tool for describing electron-phonon coupling or thermal management in microelectronic devices.

💡 Research Summary

This paper presents a novel numerical method for simulating electron‑phonon coupling in metals, based on the two‑temperature Boltzmann transport equations (BTE). Traditional two‑temperature models (TTM) and many existing kinetic solvers (e.g., lattice Boltzmann, explicit discrete ordinate, Monte‑Carlo) either assume diffusive transport or suffer from severe time‑step restrictions tied to the carrier relaxation time τ and mean free path λ. To overcome these limitations, the authors develop a semi‑implicit Lax‑Wendroff kinetic scheme that integrates the evolution of the physical equations directly into the discretization process.

The governing equations for electrons and phonons are written as ∂f_e/∂t + v_e·∇_x f_e = (f_eq_e – f_e)/τ_e – G/(4π)(T_e – T_p)
∂f_p/∂t + v_p·∇_x f_p = (f_eq_p – f_p)/τ_p + G/(4π)(T_e – T_p) ,
where f denotes the distribution function, v the group velocity, τ the relaxation time, C the specific heat, and G the electron‑phonon coupling constant. Macroscopic quantities (energy density, heat flux, temperature) are obtained by angular moments of f.

The core of the new scheme is a two‑stage time integration. First, the distribution at the half‑time step (tⁿ⁺¹/₂) is computed at cell interfaces using a combination of forward‑Euler for the source terms and backward‑Euler for the advection terms, together with a second‑order spatial interpolation for the gradients. This yields a semi‑implicit update that is second‑order accurate in both space and time. The half‑step values are then used to update the macroscopic fields at the interfaces, from which the equilibrium distribution functions are recomputed. Finally, the full‑step values (tⁿ⁺¹) at cell centers are obtained by applying the same Lax‑Wendroff discretization, now using the freshly updated interface values. Because the relaxation terms are treated implicitly, the allowable time step is no longer limited by τ, and the spatial grid can be coarser than λ without loss of accuracy.

The authors validate the method with quasi‑1D heat‑conduction simulations in gold (Au) films of varying thickness (10 nm to 10 µm). Physical parameters are λ_e = 33 nm, λ_p = 1.5 nm, τ_e = 0.0243 ps, τ_p = 0.679 ps. Isothermal boundaries are imposed at the two film surfaces. Results show a clear transition: for 10 nm films the phonon temperature profile is highly nonlinear, indicating ballistic transport; for 10 µm films the profile becomes nearly linear, consistent with Fourier diffusion. The semi‑implicit Lax‑Wendroff results match those from TTM, explicit discrete ordinate method (DOM), and discrete unified gas‑kinetic scheme (DUGKS) across the whole Knudsen‑number range, while allowing much larger Δt and Δx. Convergence is achieved when the relative temperature change falls below 10⁻⁸.

Key advantages of the proposed scheme are:

1. Relaxation‑time‑independent time step – the method remains stable for CFL numbers up to ~0.4 regardless of τ, enabling efficient simulation of ultra‑fast processes.
2. Unified treatment of migration, scattering, and coupling – all three mechanisms are incorporated within a single time step, preserving the physical coupling between electrons and phonons.
3. Second‑order accuracy – both spatial and temporal discretizations achieve O(Δx², Δt²) error, which is confirmed by the smooth transition between ballistic and diffusive regimes.
4. Computational efficiency – the additional cost of solving the interface distribution is modest; overall CPU time is reduced compared with explicit DOM or Monte‑Carlo for the same accuracy.

The paper also discusses limitations. The current implementation assumes a single relaxation time and a single group velocity for each carrier type, which neglects the multi‑band nature of electrons and the spectrum of phonon modes in real metals. Extending the method to three‑dimensional, anisotropic geometries will require more sophisticated angular quadrature and interface reconstruction. Moreover, the coupling constant G is treated as a material constant; incorporating its temperature or carrier‑density dependence would improve realism.

In conclusion, the semi‑implicit Lax‑Wendroff kinetic scheme offers a powerful, scalable tool for multiscale heat‑transfer analysis in micro‑ and nano‑electronic devices. By removing the restrictive time‑step constraints of traditional kinetic solvers while retaining full kinetic fidelity, it bridges the gap between ballistic and diffusive heat transport. Future work integrating multi‑band electron models, full phonon dispersion, and temperature‑dependent coupling will further enhance its applicability to next‑generation thermal management challenges.