Exponential Runge-Kutta methods for stiff kinetic equations

1 Introduction

• Explicit exponential Runge–Kutta integrators are proposed for homogeneous Boltzmann and related kinetic equations.The approach is based on classical time splitting and decomposition of the collision gain into equilibrium and non-equilibrium parts.Source: §1 · pp. 1–2
• The methods work uniformly over a wide range of Knudsen numbers and avoid nonlinear systems even in stiff regimes.The paper derives sufficient conditions for unconditional asymptotic stability and asymptotic preservation; methods up to order 4 are exact for BGK-type equations.Source: §1 · p. 2
• For homogeneous equations, the Maxwellian truncation criterion of is equivalent to the penalization method of; stiff-regime splitting accuracy is outside this paper’s scope.The authors identify the optimal truncation criterion as an open problem.Source: §1 · pp. 1–2
• The approach extends to large stiff ODE systems with dissipative relaxation operators and conserved quantities, beyond the Boltzmann prototype.Examples include discrete kinetic theory, shallow water equations, granular gases, and traffic flows.Source: §1 · p. 2
• The decomposition separates a stiff linear relaxation term from a nonlinear part, allowing exponential methods to solve the linear component exactly and treat the nonlinear component explicitly.The linearization is on the asymptotically large time scale, rather than the short time scale used by standard linearization techniques.Source: §1 · p. 2

2 The Boltzmann equation

• The Boltzmann collision operator preserves local collision invariants and satisfies an entropy inequality, while its equilibria are local Maxwellians.Integrating against the invariants yields non-closed conservation laws for the macroscopic moments.Source: §2 · p. 3
• As the mean free path tends to zero, the distribution approaches a local Maxwellian and the macroscopic variables satisfy the compressible Euler system.For small nonzero Knudsen numbers, the Chapman–Enskog expansion gives the compressible Navier–Stokes equations as a second-order approximation in ε.Source: §2 · p. 4

2.1 Operator splitting and stability definitions

• The analysis uses time splitting between free transport and collisions; Lie–Trotter splitting is first order, while Strang splitting is second order.Splitting methods above order two cannot avoid negative time steps in both transport and collision, which may compromise positivity and Monte Carlo use.Source: §2.1 · pp. 4–5
• The paper restricts its analysis to operator-splitting schemes; extension to non-splitting schemes is under study.Accuracy properties of splitting methods in stiff regimes are explicitly outside the manuscript’s scope.Source: §2.1 · pp. 4–5
• A time discretization is asymptotic preserving when, as ε→0, it becomes a consistent discretization of the reduced system independently of initial data and step size.This definition does not require preservation of time accuracy in the stiff limit; that stronger property is called asymptotic accuracy.Source: §2.1 · p. 4
• For the split scheme, an AP relaxation step projects onto the local Maxwellian, after which the transport step gives a kinetic approximation to Euler.Thus AP of the relaxation step suffices for AP of the full splitting, with analogous results for higher-order splitting methods.Source: §2.1 · p. 5
• Unconditional asymptotic stability requires numerical solutions from equal-moment initial data to converge together as the step count grows, independently of step size.The section also defines unconditional contractivity and entropy monotonicity; higher-order unconditionally stable Runge–Kutta schemes do not satisfy the entropy inequality without a step restriction.Source: §2.1 · pp. 5–6

3 Exponential methods

• Exponential methods avoid costly fully implicit collision solves by separating a stiff dissipative linear relaxation term from the stiff nonlinear deviation term and treating the linear part exactly.The decomposition is useful when initial data are far from local equilibrium, where the problem is stiff as a whole.Source: §3 · pp. 6–7

3.1 Explicit exponential Runge-Kutta schemes

• Explicit exponential Runge–Kutta schemes apply an underlying explicit Runge–Kutta method to the integrating-factor equation; integrating-factor and exponential-time-differencing variants are described.For integrating-factor schemes, the accuracy order matches that of the underlying method.Source: §3.1 · pp. 7–8

3.2 Time Relaxed methods

• Time-relaxed schemes arise by truncating a Taylor expansion of the transformed equation and inherit asymptotic preservation, nonnegativity, contractivity, and entropic stability under the stated density and convex-weight properties.They rely on bilinearity of P(f, f), unlike the general exponential Runge–Kutta family.Source: §3.2 · pp. 9–10

3.3 Main properties

• Under coefficient conditions, general exponential Runge–Kutta schemes are unconditionally contractive and stable in the d2 metric, and asymptotic preserving; for Maxwellian molecules, the contractivity argument uses the Boltzmann equation’s nonexpansion in that metric.For integrating-factor schemes, a nonnegative-coefficient underlying Runge–Kutta method of matching stage count and order gives unconditional contractivity and asymptotic stability; the stated coefficient condition also ensures asymptotic preservation.Source: §3.3 · pp. 10–12
• Nonnegative exponential coefficients and weights yield unconditional positivity and entropy behavior; for integrating-factor methods, sufficient underlying Runge–Kutta conditions also give convexity, which supports Monte Carlo evaluation.The underlying-method conditions are restrictive; the cited examples include second-order midpoint and Runge methods and third-order Heun, but not classical fourth-order Runge–Kutta.Source: §3.3 · pp. 13–14

3.4 Generalizations and implementation

• The parameter μ used to estimate the collision-operator spectrum can be chosen at each time step to improve the estimate, with choices discussed for variable hard spheres.Source: §3.4 · p. 14
• Positivity-enforcing μ choices can overestimate the collision-operator spectrum, especially in Monte Carlo simulations; more accurate μa or μs choices do not ensure P(f, f) ≥ 0.Their stability may require a time-step restriction.Source: §3.4 · p. 15
• As suggested in, μ can estimate the spectral radius of Q linearized around the Maxwellian M, but this estimate is easily computed only in deterministic settings.Source: §3.4 · p. 15
• For first-order IF methods, μ ≥ μp − ε/∆t ensures unconditional positivity, while contractivity and asymptotic stability hold under a weaker restriction.The passage illustrates the condition for ν = 1; it omits the weaker restriction’s expression.Source: §3.4 · pp. 15–16
• Higher-order IF schemes are asymptotic preserving provided their underlying method satisfies Proposition 4.Source: §3.4 · p. 16

4 Numerical results

• The tests apply first-, second-, and third-order integration-factor schemes to homogeneous fourth-moment evolution and an inhomogeneous shock problem.A large particle count and statistical averaging are used to focus on time-discretization error; all tested methods satisfy Proposition 5 and can use Monte Carlo collision-integral evaluation.Source: §4 · p. 16

4.1 Homogeneous relaxation

• The fourth-moment L2 error shows the expected convergence rate for all schemes, while hard spheres have larger error than Maxwell molecules at matched μΔt/ε.For hard spheres, μ is an upper bound for the collision cross section, and the test does not reduce possible overestimation of μ.Source: Figure 1 · p. 17, §4.1 · p. 18

4.2 A shock wave computation

• At the smallest tested Knudsen number, both schemes have larger errors but capture departure from equilibrium, with the second-order scheme reproducing the profile with sufficient accuracy.The Knudsen number is one tenth of the time step in this case.Source: §4.2 · p. 18
• As the Knudsen number decreases, the first-order scheme develops discrepancies from the reference sooner, while the second-order scheme remains in good agreement at the halved value.Both schemes agree well with the reference at ε = 10^-3.Source: §4.2 · p. 18

• The schemes generalize earlier methods and share the fundamental property of asymptotic preservation.Source: §5 · p. 18
• The schemes were developed for dilute-gas Boltzmann equations, but extend straightforwardly to collisional kinetic equations with a smooth equilibrium state.Source: §5 · p. 18
• Using the ES-BGK equilibrium function instead of the standard Maxwellian offers a more accurate penalization of collision dynamics.Source: §5 · p. 18
• Extending the schemes to the non-splitting case is left for future research.Source: §5 · pp. 18, 20