Low-Rank IMEX Schemes Overview

• Low-rank IMEX schemes are time integration methods that combine dynamical low-rank approximation with implicit–explicit discretizations to efficiently resolve stiff kinetic problems.
• They reduce high-dimensional computational costs via a separable low-rank ansatz and maintain stability through QR orthonormalizations and projector-splitting strategies.
• These methods exhibit asymptotic-preserving properties, ensuring accurate limiting behavior in stiff regimes while preserving key physical characteristics.

Low-rank IMEX (implicit–explicit) schemes are time integration methods designed for kinetic equations where the dominant dynamics exhibit tensor structure and stiffness arises, e.g., from collision terms. In the context of the Vlasov–Poisson–Fokker–Planck (VPFP) system, such schemes combine dynamical low-rank (DLR) approximation—where the high-dimensional solution $f(x,v,t)$ is approximated in a matrix product form—with IMEX integrators capable of stably and efficiently resolving stiff source terms. These methods exhibit efficient computational scaling, facilitate the preservation of physical properties, and—when appropriately designed—possess asymptotic-preserving (AP) features that guarantee accurate limiting behavior as stiffness parameters vanish (Zhang et al., 17 Jan 2026).

1. Dynamical Low-Rank Ansatz and Problem Structure

Central to the approach is the low-rank ansatz for the solution:

$$f(x,v,t) \approx \sum_{i,j=1}^r X_i(x,t)\,S_{ij}(t)\,V_j(v,t),$$

where $\{X_i\} \subset L^2(\Omega_x)$ and $\{V_j\} \subset L^2(\mathbb{R}^{d_v})$ are orthonormal in $x$ and $v$, respectively, and $S(t)\in\mathbb{R}^{r\times r}$ is a coefficient matrix (see eq (2.6), (Zhang et al., 17 Jan 2026)). The orthonormality conditions $\langle X_i,X_{i'}\rangle_x = \delta_{ii'}$ and $\langle V_j,V_{j'}\rangle_v = \delta_{jj'}$ guarantee an efficient and stable manifold representation. This parametrization reduces the degrees of freedom from $N_x N_v$ to $r(N_x + N_v + r)$, for modest $r$.

2. IMEX Time Discretizations: First- and Second-Order Low-Rank Schemes

The time-stepping is organized by applying a hybrid IMEX treatment to the low-rank manifold. For reference, the underlying full-tensor IMEX step for the discretized equation is:

$$(f^{n+1}-f^n)/\Delta t = \mathcal{A}_h(f^n) + (1/\epsilon) \mathcal{L}_h[E^{n+1}](f^{n+1}),$$

where $\mathcal{A}_h$ implements advection and $\mathcal{L}_h$ is the Fokker–Planck collision operator.

First-Order Low-Rank IMEX (Algorithm 1, eqs (4.7)–(4.9)):

• K-step: Backward Euler (implicit in $\mathcal{L}_h$, explicit in $\mathcal{A}_h$), updating coefficients along $V_j$, uses a density-prediction update.
• QR orthonormalization: Ensures the updated factors form an orthonormal basis.
• S-step: Forward Euler for the coefficient matrix $S$, fully explicit.
• L-step: Backward Euler (implicit in $\mathcal{L}_h$, explicit in $\mathcal{A}_h$), analogous to the K-step but updates along $X_i$.
• QR orthonormalization: On completion, reconstructs the orthonormal basis for velocity factors.

The scheme applies a Lie–Trotter K–S–L splitting—projecting onto the tangent space of the fixed-rank manifold at each step (see Section 2.2, eqs (2.8)–(2.10)).

Second-Order Low-Rank IMEX (Algorithm 2):

The second-order scheme employs a Strang-type projector splitting:

$$K(\tfrac{\Delta t}{2}) \rightarrow S(\tfrac{\Delta t}{2}) \rightarrow L(\Delta t) \rightarrow S(\tfrac{\Delta t}{2}) \rightarrow K(\tfrac{\Delta t}{2}),$$

where in each substep a two-stage IMEX Runge–Kutta update analogous to Filbet–Jin (2010) is executed (see eqs (4.10)–(4.11), Section 4.2).

3. Separable Discretization of the Fokker–Planck Operator

Spatial discretization of the Fokker–Planck operator is structured to enable efficient low-rank projections. In 1D1V, the operator is rendered via a three-point stencil (eq (3.4)), admitting the form:

$$[\mathcal{L}_h(f)]_{p,q} = (1/\Delta v^2)\big[ \alpha_{q+1/2}\beta_p f_{p,q+1} - \left(\frac{1}{\alpha_{q+1/2}\beta_p} + \alpha_{q-1/2}\beta_p \right) f_{p,q} + \frac{1}{\alpha_{q-1/2}\beta_p} f_{p,q-1} \big],$$

with $\alpha_{q+1/2} = \exp( \Delta v (v_q+v_{q+1})/4 )$ and $\beta_p = \exp( -\Delta v E_p/2 )$ (eq (3.7)). This allows $\mathcal{L}_h$ to be factorized as $$\beta\otimes T^{(\alpha)} + \beta^{-1}\otimes T^{(1/\alpha)}$$, i.e., as sums of tensor products. Consequently, the contractions required for the K, S, L updates incur computational cost $O(r^2N_v)$ or $O(r^2N_x)$, rather than $O(N_xN_v)$, amplifying computational efficiency while maintaining accuracy (Section 3.1–3.2).

4. Projector-Splitting Integrators and Manifold Geometry

The integration within the rank-$r$ manifold is governed by the Dirac–Frenkel variational principle, which projects time derivatives onto the tangent space:

$$P = (I \otimes P_V) - (P_X \otimes P_V) + (P_X \otimes I),$$

producing three evolution subflows (the aforementioned K, S, L steps). Each subflow solves a projected evolution equation, and successive QR factorizations enforce orthonormality and rank maintenance in the basis sets (Section 2.2). Both first- and second-order low-rank IMEX schemes utilize this projector-splitting structure, ensuring compatibility with the underlying geometric constraints of the DLR manifold.

5. Asymptotic-Preserving Analysis

A central feature of the first-order low-rank IMEX scheme is its asymptotic-preserving (AP) property in a small field-fluctuation regime. The proof relies on a micro–macro decomposition $f = \rho M + g$ with $\langle g, 1 \rangle_v^h = 0$ (eq (5.1)), coercivity of the discrete Fokker–Planck operator ($$-\langle\mathcal{L}_h(f),g/M\rangle_{xv}^h \geq \gamma_h\|g/M\|_{xv}^{h,2}$$, Proposition 5.1), and estimates on the projector fluctuation ($$\|(I-P_X)\mathcal{L}_h(f)\|_{xv}^h \leq \kappa_h\|f\|_{xv}^h$$, Proposition 5.2). Under the small field-fluctuation assumption $\kappa_h M_{\text{max}} < \gamma_h$ (Assumption 5.1), the main theorem establishes:

$$\|\mathcal{L}_h^{n+1}(f^{n+1})\|_{xv}^h \leq \Theta( \epsilon\|G^{n+1}\| + \kappa_h M_{\text{max}} \sqrt{L_v} \|\rho^{n+1}\|_x^h )$$

and

$$\|f^{n+1} - \rho^{n+1} M^{n+1}\|_{xv}^h = O(\epsilon) + O(\kappa_h),$$

demonstrating that the error remains controlled as $\epsilon\to 0$ and for modest $\kappa_h$ (Theorem 5.3, Section 5).

6. Algorithmic, Implementation, and Practical Considerations

Hybrid IMEX integration is achieved by applying backward Euler (L-stable) treatment for stiff collision terms in K and L steps and forward Euler in the S-step to mitigate instability, consistent with the approaches in [54], [55] of (Zhang et al., 17 Jan 2026). Each K-step solves $N_x$ decoupled $r\times r$ linear systems. Each L-step requires inversion of a block-tridiagonal system of dimension $rN_v$ with $r\times r$ blocks. QR factorizations after each substep maintain orthogonality at $O(r^2N_x)$ or $O(r^2N_v)$ computational cost. Mass conservation is only approximate due to rank truncation; a “density-prediction” update addresses this. The numerical rank remains fixed in presented schemes; adaptive strategies (e.g., SVD-truncated splitting) are not implemented but are compatible with this structure (cf. [62], [63]).

7. Empirical Performance and Regime Robustness

Numerical experiments confirm the schemes' stability, robustness, and AP property across kinetic, fluid, and mixed regimes, achieving high fidelity at modest ranks. The approach balances accuracy with significant computational savings due to the separable spatial–velocity discretization and low-rank structure, as demonstrated in Section 6 of (Zhang et al., 17 Jan 2026).

---

Key equations, algorithmic steps, and proofs referenced above can be found in (Zhang et al., 17 Jan 2026).