High-Order Phase Space Discretizations

• High-Order Phase Space Discretizations are numerical schemes that approximate evolution equations with arbitrarily high accuracy while preserving geometric structures and invariants such as energy and symplectic form.
• They incorporate methods like discrete gradient techniques, high-order discontinuous Galerkin schemes, and specialized preconditioners to tackle complex problems in Hamiltonian systems, plasma physics, and quantum dynamics.
• These approaches enable robust simulations in curved and high-dimensional geometries, ensuring long-term stability, bound preservation, and efficient error control.

High-order phase space discretizations refer to numerical schemes that approximate evolution equations in phase space (usually including both position and momentum/velocity—or other conjugate variables) to arbitrary high order in time and/or space, with a particular emphasis on geometric structure, conservation properties, and accuracy. In complex dynamical and physical systems—ranging from Hamiltonian mechanics and plasma physics to transport in curved geometries and quantum phase-space dynamics—the development of robust, high-order, and physically faithful phase space discretizations is a foundational challenge. The literature provides a diverse toolkit, from discrete gradient methods and spectral elements to sophisticated finite volume, discontinuous Galerkin, and structure-preserving integrators tailored for specific applications.

1. Geometric and Invariant-Preserving Discretizations

A central principle in phase space discretization is the preservation of geometric structure and invariants (such as energy, symplectic form, or measure). Discrete gradient methods are a prominent approach for Hamiltonian systems, ensuring the exact conservation of the first integral.

For one-dimensional Hamiltonian systems

$\dot{x} = p, \qquad \dot{p} = -V'(x)$

with Hamiltonian $E = \frac{1}{2}p^2 + V(x)$, the classic discrete gradient method (modified midpoint rule) advances $(x_n, p_n)$ via: $$\frac{x_{n+1} - x_n}{h} = \frac{1}{2}(p_{n+1} + p_n),\qquad \frac{p_{n+1} - p_n}{h} = -\frac{V(x_{n+1}) - V(x_n)}{x_{n+1} - x_n}$$ These integrators exactly preserve the energy but are only second order. The "GR–N" method improves accuracy by replacing the fixed time step $h$ with a series expansion: $\delta = h + \sum_{k=3}^N a_k(x, p) h^k$ where coefficients $a_k(x, p)$ are chosen so that the numerical scheme matches the Taylor expansion of the exact flow to order $N$. For instance, $a_3 = \frac{1}{12} V''(x)$ and $a_4 = \frac{1}{24} p V'''(x)$. This produces integrators of arbitrarily high order that still enjoy exact conservation of the Hamiltonian and demonstrate exceptional long-term stability and phase space fidelity, even for large timesteps (Cieśliński et al., 2010).

2. High-Order Discretizations in Complex Geometries

High-order phase space discretizations are crucial for solving kinetic, transport, and geometric PDEs in non-Cartesian or curved domains:

• Curvilinear Coordinates and Conservation: In high-energy astrophysics and radiation transport, phase space advection must respect nontrivial geometry (e.g., spherical or cylindrical coordinates, or even general relativity's curved spacetimes). High-order discontinuous Galerkin (DG) methods, combined with upwind numerical fluxes and strong-stability-preserving Runge-Kutta (SSP-RK) time integration, deliver both high accuracy and strict bound-preservation ($0 \leq f \leq 1$ for distribution functions). Preservation of the divergence-free property of the phase-space velocity is essential for maintaining upper bounds, and the divergence form is handled to machine precision with rigorous proofs underpinning the bound-preserving techniques. Local limiters (e.g., Zhang–Shu scaling) further guarantee that the numerical solution remains within physical bounds, and carefully chosen CFL conditions ensure numerical stability (Endeve et al., 2014).
• Mapped Multiblock Grids for Plasma Edge Simulation: For gyrokinetic simulations in tokamak geometry, high-order finite-volume discretizations on mapped multiblock (MMB) grids permit accurate phase space discretization even across the complex separatrix/X-point topology. High-order product rules and exact treatment of divergence-free velocities (rewritten via skew-symmetric tensors and recast in surface/edge-integral form) enable robust and efficient simulation across heterogeneous geometries, maintaining conservation properties and minimizing numerical dissipation (Dorr et al., 2017).

3. Temporal and Space-Time Accuracy

Discrete Gradient Methods of Arbitrarily High Order

Recent developments generalize discrete gradient schemes to arbitrary order through order theory based on formal B-series and P-series expansions. The Average Vector Field (AVF) discrete gradient,

$$H_{\text{AVF}}(x, y) = \int_0^1 \nabla H((1-\xi)x + \xi y) d\xi$$

is especially significant due to its symmetry and the vanishing of skew-symmetric error terms. Order conditions are systematized through matching series coefficients for every rooted tree up to order $p$, yielding explicit construction of third-, fourth-, or even higher-order integrators (bypassing the need for higher-order tensors). Such methods are shown to outperform classical Runge–Kutta or Gauss–Legendre schemes in long-time invariant preservation for mechanical and learned (e.g., neural network) Hamiltonian systems (Eidnes, 2020).

Space-Time Discretizations with Error Control

Combining high-order conforming Galerkin discretizations in space with high-order (possibly discontinuous) Galerkin time discretizations allows for robust and adaptive phase space schemes, with rigorously guaranteed a posteriori error estimators independent of both spatial and temporal polynomial degree. Flux reconstructions and Radau-type reconstructions enable direct error control in unified norms (e.g., $L^2(H^1) \cap H^1(H^{-1})$), and the approach naturally accommodates arbitrary mesh refinement or coarsening between time steps—crucial for adaptive high-dimensional phase space problems (Ern et al., 2016).

High-Order Semi-Lagrangian and Exponential Integrators

High-order semi-Lagrangian methods and exponential integrators (e.g., DIRK–CF scheme) provide an efficient way to resolve convection-dominated dynamics (notably in incompressible Navier–Stokes) without severe CFL restrictions. These schemes leverage spectral element spatial discretization and explicit characteristic tracing to attain high accuracy even in the presence of index-2 DAEs (after projection to divergence-free subspaces), outperforming standard projection or IMEX methods in both accuracy and computational efficiency for high-order phase space flows (Celledoni et al., 2012).

4. Structure-Preserving and Volume-Preserving Schemes

Numerical methods that preserve volume (incompressibility of phase space) are vital for simulating Hamiltonian and kinetic systems over long times:

• Volume-Preserving Splitting for Relativistic Particles: For the Lorentz force in time-dependent electromagnetic fields, splitting the extended phase space equations into explicitly solvable volume-preserving subflows (translations, electric kicks, and magnetic rotations) and utilizing processing techniques enables explicit, high-order, volume-preserving integrators. These methods are time-symmetric and conserve invariants such as energy and angular momentum, with careful stability analysis revealing significantly enlarged permissible timestep domains for high-order processed schemes (He et al., 2016).
• Measure-Preserving Schemes for Time-Reparametrized Hamiltonian and Nonholonomic Systems: For systems that are Hamiltonian up to time rescaling (e.g., Hamiltonizable nonholonomic systems), discretizations constructed via backward error analysis, and high-order expansions of the modified/altered Hamiltonian, yield formal invariance of the smooth phase-space measure. This guarantees long-time qualitative fidelity, even if classical symplecticity is absent (García-Naranjo et al., 2020).

5. Bound- and Maximum Principle–Preserving Design

In many phase space applications, ensuring that the numerical solution respects physical bounds (e.g., positivity, upper/lower limits) is essential. This is achieved via:

• Bound-Preserving Limiters: High-order DG schemes and WENO finite volume methods are combined with local (subcell) or global monolithic convex (GMC) flux limiters, ensuring updated values remain within admissible bounds for explicit or implicit (DIRK) high-order time schemes. The two-tiered limiting strategy—first applying a high-order method, then postprocessing to enforce maximum/minimum bounds—ensures no step size restriction in fully implicit settings and guarantees optimal accuracy is preserved in smooth regions (Endeve et al., 2014, Kuzmin et al., 2020, Luna et al., 2021).
• Positivity with Optimization Techniques: For degenerate parabolic equations and gradient flows, pairing high-order LDG discretization with time-implicit DIRK schemes and enforcing positivity/mass conservation via a KKT (Karush–Kuhn–Tucker) limiter (solved as a variational inequality with Lagrange multipliers) ensures not only positivity but also unconditional entropy dissipation, matching the qualitative long-time behavior of the continuous system (Yan et al., 2023).

6. Efficient Solver Strategies for High-Order Phase Space Discretizations

High-order phase space approximations sharply increase computational demands due to large numbers of degrees of freedom and nontrivial matrix structure. Recent work advances scalable solvers and preconditioners:

• Monolithic Multigrid with $ph$-Coarsening: For high-order stationary Stokes discretizations, multigrid frameworks that coarsen both the mesh ($h$) and approximation order ($p$) demonstrate substantial reductions in setup and solve times, especially in the high-order regime. Taylor–Hood elements benefit particularly, while block-factorization techniques may be advantageous for mixed systems (e.g., Scott–Vogelius elements) (Voronin et al., 9 Jul 2024).
• Auxiliary and Fictitious Space Preconditioners in H(div): Interior penalty DG methods in H(div) (used for exactly divergence-free discretizations) are efficiently preconditioned using subspace correction (with vertex patches), fictitious space, and auxiliary space strategies, resulting in robust solvers with condition numbers weakly dependent on mesh parameter $h$, polynomial degree $p$, and penalty parameter $\eta$, as substantiated on both Cartesian and unstructured meshes (Pazner, 21 Nov 2024).
• Space-Time Multigrid for Tensor-Product Discretizations: In fully coupled space-time finite element frameworks (treating time as another geometrical dimension), space-time multigrid with tensor-product decompositions and cell-wise additive Schwarz smoothers achieve scalable, high-throughput solutions for both parabolic and hyperbolic PDEs, supporting high polynomial degree and demonstrating empirical throughput exceeding a billion DoFs per second on extreme-scale HPC problems (Margenberg et al., 8 Aug 2024).

7. Applications and Forward Directions

High-order phase space discretizations are integral not only to classical mechanics and plasma physics but extend to:

• Quantum Mechanics and Quasiprobability Representations: Discrete phase space methods constructed via the discrete Heisenberg–Weyl group and their connection to continuous $SU(N)$ quantization kernels provide a unified language for quantum-classical correspondence and simulation of strongly interacting quantum systems. Discrete positive-P representations help mitigate diverging trajectories in stochastic simulations (Žunkovič, 2015).
• Wave Transport and Direction-Preserving Discretizations: In acoustics and electromagnetic wave transport, Petrov–Galerkin boundary integral discretizations using directional subsets of globally defined direction sets preserve the correct propagation of phase-space densities across composite domains. The ability to interpolate between discrete nodal directions (using piecewise-constant or piecewise-linear test functions in the momentum coordinate) enables high-order accuracy in complex engineering environments (Chappell et al., 2023).

Summary Table: Pillars of High-Order Phase Space Discretization

Pillar	Methodologies/Papers	Key Benefits
Invariant Preservation	(Cieśliński et al., 2010, Eidnes, 2020)	Exact energy, symplectic or measure preservation
High-Order/Curved Geometry	(Endeve et al., 2014, Dorr et al., 2017)	Accuracy in non-Cartesian, multi-block geometries
Structure-Preserving Volume/Measure	(He et al., 2016, García-Naranjo et al., 2020)	Long-run stability, physical fidelity
Bound/Maximum Principle Preservation	(Kuzmin et al., 2020, Luna et al., 2021)	Positivity, bounds, maximum principle in phase space
Efficient High-Order Solvers	(Voronin et al., 9 Jul 2024, Pazner, 21 Nov 2024)	Scalability for large DoF, high-$p$
Quantum/Transport Applications	(Žunkovič, 2015, Chappell et al., 2023)	Mixed discrete/continuous phase space, directionality

High-order phase space discretizations are thus characterized by their ability to achieve arbitrarily high accuracy while preserving essential geometric, algebraic, or physical structures—including energy, invariants, bounds, and conservation laws—within arbitrarily complex spatial and temporal domains. Recent advances enable application to stiff, high-dimensional, or geometrically intricate problems, with emerging solver and preconditioner frameworks ensuring computational feasibility at scale.