A second-order structure- and positivity-preserving convex limiting method for the Vlasov equations

1 Introduction

• High-dimensional Vlasov advection makes accurate, efficient computation difficult, while simultaneously preserving mass, energy, momentum, field-divergence constraints, positivity, and asymptotic properties remains challenging.The difficulty is especially acute beyond classical 1D1V and 1D2V phase spaces.Source: §1 · p. 1
• To the authors’ knowledge, no continuous finite element approximation of the Vlasov system had rigorously established positivity preservation.The paper identifies compact representation and memory efficiency as attractive properties of continuous finite elements for kinetic plasma simulations.Source: §1 · p. 2
• High-order positivity-preserving finite element schemes require a robust first-order method, a high-order discretization, and conservative limiting, but the first-order component had been a key obstacle.Graph-based viscosity formulations supplied a continuous finite element analogue of Lax–Friedrichs; residual-based viscosity and convex limiting [24] are among approaches used for high-order schemes.Source: §1 · p. 2
• Residual-based viscosity is challenging to extend to high dimensions because projection adds computational overhead and backward differentiation can require storing multiple solution states.Replacing the time derivative with a spatial flux approximation may still leave an expensive projection problem.Source: §1 · pp. 2–3
• Artificial viscosity, mass lumping, and convex limiting can violate discrete continuity and thereby compromise discrete Gauss’s law; the paper proposes weak-form divergence cleaning to enforce that constraint.Divergence cleaning removes spurious divergence errors from computed fields.Source: §1 · p. 3
• The proposed method combines low-order maximum-principle-enforcing viscosity with high-order residual-based viscosity using convex limiting, and its tensor-product residual construction avoids storing solutions at multiple time levels.The residual formulation is designed for tensor-product finite element spaces, where memory efficiency matters for high-dimensional phase-space simulations.Source: §1 · p. 3
• The work demonstrates a provably positivity- and structure-preserving continuous finite element framework for Vlasov-type problems, realized with piecewise multilinear elements and second-order phase-space accuracy.The authors present this as a computationally efficient alternative to DG for high-dimensional settings; the second-order restriction is not intrinsic to the limiting strategy.Source: §1 · p. 3

2 Preliminaries

• The preliminaries introduce the Vlasov–Maxwell equations and a tensor-product finite element method for approximating the Vlasov equation.Source: §2 · p. 3

2.1 Governing equations

• The governing system couples the Vlasov equation to Maxwell’s equations; the Vlasov–Poisson simplification disregards the magnetic field and obtains the electric field from Poisson’s equation.The text defines the species distribution function and charge density, and assumes unit values for the speed of light, permittivity, and permeability for simplicity.Source: §2.1 · p. 4

2.2 Finite element approximations

• The phase-space finite element space is built as a tensor product of continuous piecewise bilinear spaces over configuration space and velocity space.The phase-space mesh is likewise formed from the configuration-space and velocity-space meshes.Source: §2.2 · p. 4

2.3 Tensor-product finite element method

• The tensor-product Galerkin formulation seeks a time-dependent finite element approximation and yields a system for its nodal degrees of freedom.The basis functions on phase space are products of configuration-space and velocity-space basis functions.Source: §2.2 · pp. 4–5, §2.3 · p. 5

3.1 Low-order scheme

• With the CFL condition and 0 < cfl ≤ 1, the low-order scheme satisfies a discrete maximum principle and preserves nonnegativity.The step-size restriction is given by CFL condition (11).Source: §3.1 · pp. 6–7

3.2 High-order scheme

• Standard residual projection is prohibitively expensive for tensor-product discretizations; the proposed full-Vlasov residual aims to retain robustness while reducing artificial diffusion.The earlier lower-dimensional residual approach may accumulate diffusion in both physical and velocity dimensions.Source: §3.2 · p. 8

3.2.1 A mass-lumped residual approximation

• The mass-lumped residual approximation replaces the discrete time derivative with a spatial approximation and uses spatial operators, making it suitable for high-dimensional tensor-product discretizations.The governing equation gives ∂tf = −β · ∇f; the auxiliary space uses the same polynomial space on a coarser mesh.Source: §3.2.1 · p. 8

3.2.2 Residual-viscosity method

• The residual approximation feeds the high-order residual-viscosity construction and its finite-element scheme, with the time step set by CFL condition (11).The section introduces the viscosity matrix before giving the high-order finite-element scheme.Source: §3.2.2 · pp. 8–9

3.3 Convex limiting

• The convex limiter conserves total mass; with relaxed extrema, positivity holds for 0 ≤ r_h ≤ 1, although the strict local maximum principle is no longer implied.The unrelaxed limiter preserves the discrete maximum principle.Source: §3.3 · pp. 10–11

4.1 Semi-discretization of the Maxwell’s equations

• The semidiscrete Maxwell fields use finite element spaces drawn from the de Rham complex, with E_h and B_h assigned to its corresponding electromagnetic spaces.The formulation seeks E_h and B_h for a given current density J_h.Source: §4.1 · p. 12
• The magnetic divergence constraint is preserved by the semidiscrete scheme if it holds at t = 0.The proof uses the discrete de Rham structure, for which the divergence of the curl term vanishes.Source: §4.1 · p. 12

4.2 A divergence cleaning technique

• Artificial viscosity and nonlinear limiting can disrupt the continuity equation, so the earlier correction strategy is nontrivial to extend to this method.The paper instead applies divergence cleaning after each time step.Source: §4.2 · p. 13
• After divergence cleaning, the electric field satisfies Gauss’s law in the weak form.The cleaning procedure is applied after each time step.Source: §4.2 · p. 13

4.3 Time discretization and summary of the algorithms

• SSP-RK3 is written as convex combinations of Forward Euler steps and inherits their positivity property under the same CFL condition.The same time integrator couples the Vlasov and Maxwell equations; the Maxwell update occurs at each time level.Source: §4.3 · pp. 13–14

5 Numerical experiments

• The experiments cover Vlasov–Poisson and Vlasov–Maxwell benchmarks in 1D1V, 1D2V, and 2D2V settings, generally using a CFL number of 0.3.The streaming Weibel test uses a smaller CFL number of 0.2.Source: §5 · p. 15, §5.2 · pp. 18–19

5.1.1 Landau damping

• Residual-viscosity errors are slightly larger on coarser residual spaces, but convergence rates remain optimal in all tested cases.With the convex limiter and relaxation, L1- and L2-norm errors remain comparable to the unlimited results and converge at second order.Source: §5.1.1 · pp. 15–16, Table 1 · p. 16, Table 2 · p. 17
• The convex limiter confines strong Landau-damping extrema and achieves positivity, whereas residual viscosity alone does not guarantee positivity.The Galerkin minimum reaches approximately −0.14; the limited distribution is similar to the unlimited residual-viscosity solution but spans a smaller range.Source: §5.1.1 · p. 16, Figure 1 · p. 17

5.1.2 Two-stream instability

• For k = 0.5 and α = 0.01, the two-stream distribution strictly preserves positivity at t = 10, 15, and 20, with a diagonal-cut lower bound slightly above zero.Across different mesh resolutions, the overall distribution remains qualitatively similar. Using coarser spaces V2h, V4h, and V8h produces residuals of larger magnitude, corresponding to increased numerical diffusion.Source: §5.1.2 · p. 18, Figure 3 · p. 18

5.2 1D2V Vlasov–Maxwell equations

• Divergence cleaning significantly reduces Gauss’s-law constraint violations, while the electromagnetic-field growth rate agrees with 0.03 reported in [28].A small Gauss-error drift persists because finite-element error in the Poisson solve accumulates over time.Source: §5.2 · p. 19, Figure 6 · p. 20, Figure 5 · p. 20

5.3 2D2V diocotron instability

• In the 2D2V diocotron simulation, charge density remains strictly positive and resembles the result reported in [33, Fig. 9].The run uses 10000 P1 elements in the physical domain and 202 Q1 elements in velocity space.Source: §5.3 · p. 21, Figure 7 · p. 21

• The continuous finite element scheme combines a low-order discrete-maximum-principle scheme with a high-order scheme through convex limiting.Source: §6 · p. 21
• The scheme preserves positivity of the distribution function and charge density, and achieves second-order accuracy for smooth solutions with first-order polynomial approximations.Source: §6 · p. 21
• The method has been applied to four-dimensional problems, but is computationally expensive for higher polynomial degrees and full six-dimensional simulations.Source: §6 · p. 21