A new computational solver substantially reduced the work needed to reach a converged solution in modeled heat-transport problems for suspended monolayer graphene. Across 56 benchmark points where both methods converged, the solver cut the iteration count by a median factor of 4.26 and wall-clock time by 2.46. The largest recorded gains were 55.5 times fewer iterations and 37.3 times less wall time.
The study concerns the stationary, linearized, non-gray Callaway phonon Boltzmann transport equation, a detailed model in which different heat-carrying vibration modes can have different relaxation times. The proposed general synthetic iterative scheme, or GSIS, was compared with conventional source iteration, known as CIS, using the same numerical setting.
The problem grows harder at larger scales
The central numerical problem is slow convergence at long wavelengths, meaning errors that vary gradually across the modeled material. Fourier analysis found that the CIS contraction factor, the fraction of an error left after each iteration, approaches 1 as the characteristic length increases. The analysis also found that reducing that remaining error by a fixed factor requires quadratically more iterations.
GSIS behaved differently in the tested graphene calculation. Its spectral radius, another measure of how much error survives an iteration, reached about 0.93 near a characteristic length of 0.08 micrometres, then approached 0.409 at large lengths. The corresponding large-length calculation required about 26 iterations.
An extended sweep of the iso4 configuration tested lengths through 1 millimetre. GSIS continued to converge across that range with iteration counts in the thousands. CIS reached its preset limit of 2 × 10^5 iterations at about 113 micrometres without meeting the stopping criterion, while GSIS converged in approximately 2.4 × 10^3 iterations. The resulting cap-to-GSIS ratio exceeded 80, although that figure is a lower bound because the CIS run had not converged.
The speed came without a visible change in the answer
Faster convergence would matter only if the accelerated method retained the computed result. At all 56 jointly converged finite-domain points, CIS and GSIS overlapped at the plotting scale. Their relative effective-conductivity difference had a median of 3.38 × 10−6 and a maximum of 4.39 × 10−5. The study attributes these small discrepancies to stopping at a finite tolerance rather than to a difference at the exact fixed point.
A separate analytic bulk calculation reproduced the cited reference conductivity curve for the Callaway model on the spectral grid used in the study. That comparison supports the internal consistency of the implemented material model and its numerical quadrature, but it was not independent experimental validation.
A two-level numerical design
GSIS combines a kinetic sweep with a separate solve for synthetic macroscopic equations. Its construction uses exact energy and quasi-momentum conservation, first-order Chapman–Enskog constitutive relations, and high-order terms taken from the kinetic solution. The discretization couples a nodal discontinuous Galerkin method, which represents the solution with local polynomials, to a hybridizable discontinuous Galerkin solve.
The numerical contraction study used suspended monolayer graphene and retained its LA, TA and ZA acoustic phonon branches. The main benchmark calculations used eight Gauss–Legendre frequency nodes per branch, polynomial degree 3, a 20 by 20 spatial resolution and 800 triangular elements. CIS and GSIS ran on identical meshes with the same stopping tolerance, 10−7, and the same iteration cap, 2 × 10^5.
The biggest gains appeared in difficult cases
The largest iteration reduction in the benchmark suite occurred in Case B-1 at a characteristic length of 100 micrometres, where the iteration-count speed-up was 55.5 times and the wall-time speed-up was 35.5 times. The largest wall-time gain occurred in Case B-2 at 10 micrometres and 600 kelvin, reaching 37.3 times.
The implementation was also tested on two unstructured configurations, including a cavity with an internal curved boundary and diffuse non-thermalizing boundaries. The two methods produced the same resolved temperature distribution on the plotted scale in those tests.
Evidence remains limited to the tested model
The reported timing advantage is not a universal hardware figure. Wall-clock ratios depend on the implementation, hardware, thread allocation and the extra hybridizable discontinuous Galerkin solve, so iteration counts are the more portable comparison. The long-range sweep also does not show that capped CIS runs would never converge; it shows only that they failed to meet the prescribed tolerance within the stated cap.
The asymptotic conclusions depend on a nonsingular fixed-point increment problem and a stated spatial-discretization requirement. Within those assumptions, the analysis reported compatibility between the synthetic and kinetic fields and recovery of a Guyer–Krumhansl-like hydrodynamic equation and Fourier’s heat-conduction law in the corresponding limits.
The document is an arXiv version 1 preprint dated 28 August 2026. The supplied material reports no funding source.
Paper data and sources
Original title: A Synthetic Iterative Scheme for Non-Gray Phonon Boltzmann Transport Equation with Dual Relaxation Times
Authors: Dingtao Shen, Jia Liu, Wei Su
Journal/Repository: arXiv
Status: Preprint, not yet peer-reviewed
First online: 2026-08-28
DOI: Not available
Original paper · Full text