Neural network learns Boltzmann collision dynamics with decomposed relaxation structure

Solving the Boltzmann equation exactly is one of the most computationally punishing tasks in computational physics. The equation describes a gas through the evolution of its molecular velocity distribution, and its collision operator couples every molecular velocity to every other through binary interactions. When a gas departs from local equilibrium, as in rarefied flows or regions with steep gradients, this operator determines how the distribution relaxes as it is transported through physical space. Directly evaluating the collision integral over a discretized velocity space requires multiple integrations over relative velocity and scattering direction, repeated at every physical grid point and every time step, which is why reduced collision models have long been essential to time-dependent kinetic simulations.

A team of researchers has now introduced a relaxation-decomposed network, or RDN, that learns the local collision response directly from distribution data while keeping the recognizable structure of a classical relaxation model. The work, published in Results in Engineering, replaces the expensive evaluation of the Boltzmann collision integral with a neural approximation that can be embedded in a kinetic solver. The approach follows the Bhatnagar–Gross–Krook (BGK) tradition, which replaces the full collision operator with relaxation toward a local Maxwellian at a single scalar relaxation time. That compact form makes collision evaluation inexpensive, but a single relaxation rate offers limited control over how different molecular velocities respond, and extensions of relaxation models must still be assessed through quantities that resolve nonequilibrium transport rather than density, velocity, and temperature alone.

The new architecture splits the collision response into two separately learned branches. A residual distribution branch takes the deviation between the current distribution and the local Maxwellian as input and generates five state-dependent coefficients that build an exponential correction of the Maxwellian. A second residual branch takes the same deviation, together with a prescribed reference relaxation time, and outputs a positive relaxation time at each discrete velocity node. The two branches therefore modify the distribution difference and its rate of relaxation through separate learned maps. Because the exponential function is strictly positive and the exponential linear unit used in the relaxation branch keeps the correction above a critical threshold, positivity of both the parameterized distribution and the relaxation times is enforced directly by the model structure rather than by the size of the training error.

The mathematical analysis distinguishes what each branch contributes near equilibrium. Expanding the collision operator around a collision fixed point shows that the linearized response is determined jointly by the Jacobian of the corrected distribution map and the nodewise relaxation frequencies at the fixed point. At first order, the relaxation-time branch acts only through those fixed-point frequencies, because its input-dependent variation multiplies a distribution difference that vanishes at equilibrium. At second order, however, the curvature of the distribution map and the input-dependent modulation of the relaxation frequencies contribute separately to the nonlinear response. The authors note that exact equilibrium consistency, in which the learned collision response vanishes when the input is the discrete Maxwellian, holds only if the distribution network satisfies a zero-input condition that a biased network need not meet, and the training objective contains no separate equilibrium-residual term to guarantee it.

Training uses reference collision responses generated by a deterministic kinetic solver. For a one-dimensional smooth wave on a periodic domain, full Boltzmann trajectories and their local collision responses were computed with the DGFS solver on a physical grid of 1280 points and a molecular-velocity domain represented by a uniform tensor-product grid of 25 cubed nodes. Twenty training trajectories were sampled with Knudsen numbers drawn logarithmically between 0.01 and 1, with density and temperature baselines between 0.8 and 1.2 and perturbation amplitudes between 0.1 and 0.7, yielding 10,000 training samples. Each sample contains the discrete distribution, the BGK reference relaxation time, and the reference collision response at the same time and physical location. The loss combines the squared error of the collision response with soft penalties on mass and momentum conservation residuals, while energy is supervised only through the full collision target.

The wave test allows a strict examination of the learned collision term. Across Knudsen numbers of 0.01, 0.1, and 1.0, RDN reproduces the locations of the maxima and minima of the reference density, velocity, and temperature profiles and preserves their phase relationships, including a streamwise velocity that develops from an initially quiescent state through a nonuniform pressure field. Deviations concentrate around the extrema where spatial gradients are largest. The learned relaxation time varies with both physical position and molecular velocity, producing alternating bands of faster and slower relaxation in velocity space, precisely the degree of freedom absent from the scalar BGK relaxation time. The mass, momentum, and energy residuals of the learned response remain below the BGK residuals over most of the domain and are comparable to or smaller than those of the Shakhov model in the displayed state, although the energy residual stays finite because no energy-invariant projection is applied.

Ablation experiments confirm that both branches are needed. When only the corrected-distribution branch is active, the relative collision-response error is approximately 0.13; when only the velocity-dependent relaxation-time branch is active, it is approximately 0.15. Combining the two branches reduces this error to approximately 6.7 times ten to the minus three. The one-step distribution error follows the same ordering, falling from roughly 2.4 and 2.9 times ten to the minus three for the single-branch variants to roughly 1.3 times ten to the minus four for the full model. The two branches therefore do not provide interchangeable corrections: the distribution branch changes the state toward which the local distribution relaxes, while the relaxation-time branch controls the velocity-dependent rate of that relaxation.

The benchmarks grow progressively harder. In a Sod shock tube, RDN reproduces the principal plateaus, extrema, and propagation directions of the reference solution at three reported times, but the discriminating information lies in the narrow transition layers. Heat flux and normal-stress differences, which directly measure nonequilibrium transport through the shock layer, expose discrepancies barely visible in the primary fields, with differences concentrated in peak heights and local widths. In a three-dimensional lid-driven cavity at a Knudsen number of 1, the relative L2 error against DSMC reference data is 0.45 percent for density, 6.98 percent for the shear stress component Pxy, and 17.33 percent for Pyz, a secondary stress component associated with weaker three-dimensional shear. The correlations remain high, at 0.9971 for density and Pxy and 0.9856 for Pyz, indicating that the model captures the shape and location of the stress structures more accurately than their amplitudes at the most demanding points.

The final benchmark extends the assessment to a time-dependent argon–krypton mixture in oscillatory Couette flow, with wall velocities of 50 and 500 meters per second and species-pair Knudsen numbers near or above one. RDN follows the species-resolved wall-parallel velocity profiles for both gases at four phases of the forcing cycle, reproducing profile curvature, interior extrema, and the phase-dependent reversal of the flow, with the remaining deviations concentrated near the driven wall where local gradients are strongest. Because a separate network is trained for each problem, the experiments assess the learned model under the stated conditions rather than transfer between flow configurations. Across all four problems, the authors conclude that the approach supports relaxation-decomposed learning as a structured approximation of collision dynamics within the tested configurations, while identifying exact conservation enforcement, accuracy in strong-gradient moments, and three-dimensional stress prediction as the most direct targets for further development.

Subject of Research: A relaxation-decomposed neural network that learns nonequilibrium Boltzmann collision responses for kinetic gas-flow simulations

Article Title: Learning nonequilibrium Boltzmann collision dynamics with a relaxation-decomposed network

Article References: Feng, M., Shan, M., Kuai, L., Yang, Y., Yin, C., & Han, Q. (2026). Learning nonequilibrium Boltzmann collision dynamics with a relaxation-decomposed network. Results in Engineering, 32, Article 113740. https://doi.org/10.1016/j.rineng.2026.113740

Image Credits: AI Generated

DOI: 10.1016/j.rineng.2026.113740

Keywords: Boltzmann equation, collision operator, neural network, BGK model, rarefied gas dynamics, nonequilibrium moments, Sod shock tube, lid-driven cavity, Couette flow, machine learning, kinetic theory, conservation laws