A Conservative Discontinuous Galerkin Algorithm for Particle Kinetics on Smooth Manifolds

This paper presents a novel, conservative discontinuous Galerkin algorithm for simulating particle kinetics on smooth manifolds that utilizes Hamiltonian formulations to exactly conserve density and energy, incorporates a BGK collision operator with an iterative scheme for preserving collisional invariants, and demonstrates its efficacy through various test cases including rotating geometries and shock problems.

The Gist

Imagine you are trying to simulate how a swarm of tiny, invisible bees moves around inside a complex, curved room. Maybe the room is shaped like a perfect sphere, or maybe it's a wobbly, saddle-shaped surface. In the real world, these bees (particles) don't just fly in straight lines; they follow the curves of the room, and sometimes they bump into each other.

This paper presents a new, highly accurate computer program designed to track these bees without making mistakes or adding "fake" noise to the simulation. Here is how the authors did it, explained in everyday terms:

1. The Map and the Compass (Hamiltonian Systems)

To tell the bees where to go, the authors use a special kind of map called a Hamiltonian. Think of this as a master rulebook that tells every bee exactly how to move based on the shape of the room.

• The "Canonical" Rulebook: The authors found a special way to write these rules (using "canonical coordinates") that makes the math incredibly clean and efficient. It's like having a compass that always points true north, no matter how twisty the path gets. This method ensures that the total number of bees and their total energy never magically appear or disappear during the simulation.
• The "Non-Canonical" Rulebook: Sometimes, the "perfect" compass is hard to use because the room is too weirdly shaped. The authors also created a backup set of rules (non-canonical) that is a bit messier but works better for specific shapes, like a polar map where distances get squished near the center.

2. The Digital Tiles (Discontinuous Galerkin)

Instead of trying to draw the whole room as one giant, smooth picture, the authors chop the room into millions of tiny, separate tiles.

• Imagine a mosaic. Each tile has its own little drawing of how the bees are moving inside it.
• The magic of their method is that they can talk to the neighbors on the edges of these tiles to make sure the bees flow smoothly from one tile to the next.
• Why this is cool: Because they use these tiles, they can use very high-resolution math (like a super-high-definition camera) without needing a supercomputer the size of a city. It's efficient and precise.

3. The "Bump" and the "Bounce" (Collisions)

In the real world, bees bump into each other. The authors added a special "bump" mechanic to their simulation.

• The BGK Operator: This is a simplified way to model collisions. Imagine that if the bees get too chaotic, this mechanic gently nudges them back toward a calm, organized state (like a teacher calming down a noisy classroom).
• The Safety Net: They built a special "iterative" loop (a check-and-fix cycle) into the code. After every bump, the computer checks: "Did we accidentally lose a bee? Did we create extra energy?" If the answer is yes, the loop fixes it immediately. This ensures the simulation stays physically honest.

4. Spinning Rooms (Rotation)

The authors also tested what happens if the room itself is spinning, like a merry-go-round.

• They showed that by tweaking the "rulebook" (the Hamiltonian) just a little bit, they could account for the spinning. This is crucial for simulating things like gas swirling around a spinning black hole or a neutron star.
• They proved that even with the spinning, their method still conserves energy and particle counts perfectly.

5. The Tests (Did it work?)

To prove their new program works, they ran three famous "stress tests":

• The Sod Shock: They created a scenario where a wall of gas suddenly breaks, creating a shockwave. They showed that their computer simulation matched the exact mathematical answer perfectly, even when the gas was bumping into itself a lot (fluid limit) or not at all (collisionless limit).
• The Kelvin-Helmholtz Instability: They simulated two streams of gas sliding past each other on a sphere and a saddle shape. This usually creates beautiful, swirling "cat's eye" patterns. Their simulation captured these swirls with incredible detail, showing exactly how the gas behaves without the "static" or "graininess" that plagues other methods.
• The Rotating Sphere: They tracked a single "blob" of gas moving on a spinning sphere. The blob followed the exact path predicted by physics, including the weird curves caused by the spin (Coriolis force).

The Bottom Line

The authors have built a new, robust tool for simulating how particles move on curved surfaces.

• It's conservative: It never loses or gains energy or particles by mistake.
• It's quiet: Unlike other methods that are "noisy" (like static on a radio), this one gives a clean, clear picture of the physics.
• It's flexible: It works on flat floors, curved spheres, and spinning worlds.

The paper concludes by saying this tool is a stepping stone. While they tested it on non-relativistic (non-light-speed) scenarios, the same mathematical foundation can eventually be used to simulate the extreme gravity around black holes and neutron stars, helping us understand the universe's most violent environments.

Technical Summary

Technical Summary: A Conservative Discontinuous Galerkin Algorithm for Particle Kinetics on Smooth Manifolds

Problem Statement
The study of kinetic theory in strong gravitational fields, such as those surrounding black holes or neutron stars, presents a fundamental challenge. While General Relativistic Magnetohydrodynamics (GRMHD) is the standard for modeling such systems, it relies on fluid closures that often fail to capture non-thermal effects or violate causality in relativistic regimes. Conversely, General Relativistic Particle-in-Cell (GRPIC) methods, while flexible, suffer from intrinsic statistical noise that can artificially thermalize distributions and obscure delicate phase-space structures. Directly discretizing the Vlasov equation in 6D phase space offers a noise-free alternative, but existing continuum methods are often computationally expensive and complex to implement on curved manifolds. This paper addresses the need for a high-order, conservative, and noise-free numerical algorithm capable of simulating particle kinetics on smooth manifolds, serving as a stepping stone toward general relativistic kinetic theory.

Methodology
The authors develop a Discontinuous Galerkin (DG) scheme for the kinetic equation on $d$-dimensional Riemannian manifolds. The core of the methodology involves:

1. Hamiltonian Formulation: The motion of free-streaming (collisionless) particles is formulated using Hamiltonian mechanics. The authors present two distinct formulations:

   • Canonical Formulation: Uses covariant momentum components ($p_i$) as coordinates. The Poisson bracket is the standard canonical form, and the Hamiltonian is $H = \frac{1}{2}g^{ij}p_i p_j$. This formulation yields geodesic motion without explicitly requiring Christoffel symbols in the update equations.
   • Non-Canonical Formulation: Uses contravariant momentum components ($p^i$). This leads to a complex, non-canonical Poisson bracket where the metric and its derivatives appear explicitly. This approach is introduced to handle coordinate singularities (e.g., near $r=0$ in polar coordinates) where canonical momentum resolution becomes inefficient.

2. Discontinuous Galerkin Discretization:

   • The phase space is discretized into cells using piecewise polynomial spaces (Serendipity basis).
   • The scheme enforces the continuity of the discretely represented Hamiltonian and Poisson tensor components across cell interfaces.
   • A Lax flux is used for the surface terms, allowing for either central fluxes (for $L^2$ conservation) or upwind fluxes (for monotonic $L^2$ decay).
   • Time advancement is performed using a Strong Stability Preserving (SSP) third-order Runge-Kutta scheme.

3. Collisional Physics:

   • A Bhatnagar-Gross-Krook (BGK) collision operator is employed to model relaxation toward local thermodynamic equilibrium.
   • To avoid the severe time-step restrictions of explicit collisional terms (especially in the fluid limit where collision frequency $\nu \to \infty$), the authors use a splitting strategy. The advection term is advanced explicitly, while the collision term is treated implicitly using a first-order Euler step.
   • Iterative Correction: A critical component is an iterative Picard algorithm that corrects the discretely projected equilibrium distribution ($f_M$) to ensure its moments (density, momentum, energy) exactly match the moments of the distribution function. This ensures strict conservation of collisional invariants.

4. Rotating Manifolds: The framework incorporates rotation (e.g., a uniformly rotating sphere) by modifying the Hamiltonian to include a Coriolis-like term ($-\Omega p_\phi$) while maintaining the canonical structure.

Key Contributions

• Conservative DG Scheme on Manifolds: The paper constructs a DG scheme that exactly conserves total particle number and total energy for time-independent Hamiltonians.
• Stability Proofs: The authors prove that for the canonical formulation, the scheme conserves the $L^2$ norm of the distribution function when using central fluxes and monotonically decays it when using upwind fluxes. This relies on the proof that the discrete phase space is incompressible in the canonical limit.
• Asymptotic Preserving Property: The scheme is shown to be asymptotic preserving; in the high-collisionality limit, the kinetic solver naturally recovers the Euler equations (fluid limit) and satisfies Rankine-Hugoniot conditions without ad-hoc fluid closures.
• Non-Canonical Extension: The paper derives a non-canonical bracket using normalized momentum coordinates to mitigate geometric dependencies in momentum space resolution, offering a trade-off between computational efficiency and provable $L^2$ stability.

Results and Benchmarks
The algorithm was implemented in the Gkeyll code and tested on several problems:

• Sod Shock Tests: In both 1D flat space and 2D annular disk geometries, the kinetic solver reproduced the exact Riemann solution in the fluid limit ($\nu = 15,000$) and correctly modeled free-streaming in the collisionless limit. The scheme demonstrated machine-precision conservation of total particles and energy, and exact conservation of angular momentum in axisymmetric cases.
• Kelvin-Helmholtz Instability (KHI): Simulations on the surface of a sphere and a hyperbolic surface successfully reproduced the KHI vortex formation in the fluid limit. The continuum approach allowed for the extraction of the deviation from equilibrium ($f - f_M$) with high accuracy ($O(10^{-5})$), enabling the calculation of transport coefficients.
• Rotating Sphere: A test case involving a Gaussian bump on a uniformly rotating sphere demonstrated agreement between the collisionless kinetic solver and single-particle geodesic trajectories, correctly capturing Coriolis and centrifugal effects.
• Convergence: The scheme showed stable convergence without spurious oscillations around shocks, attributed to the $L^2$ stability of the canonical bracket scheme.

Significance and Future Outlook
The paper establishes a robust, first-principles numerical framework for kinetic theory on curved manifolds. Its significance lies in providing a noise-free, conservative alternative to PIC methods that can resolve detailed phase-space structures, which is essential for understanding transport and instabilities in complex geometries.

The authors explicitly state that this work is a precursor to developing tools for General Relativity (GR). While the current implementation handles 2D manifolds embedded in 3D Euclidean space, the structural similarity between the canonical DG scheme on manifolds and the GR Hamiltonian suggests that the method can be extended to 4D spacetime. The authors note that a full GR implementation will require further mathematical development, particularly regarding the use of tetrad frames to handle collisions and intrinsic curvature, which will be presented in future work. The current approach provides a "natural extension" for modeling systems ranging from fusion plasmas to astrophysical environments like accretion disks and pulsar magnetospheres.