Modeling of Relativistic Plasmas with a Conservative Discontinuous Galerkin Method

We present a new method for solving the relativistic Vlasov–Maxwell system of equations, applicable to a wide range of extreme high-energy-density astrophysical and laboratory environments. The method directly discretizes the kinetic equation on a high-dimensional phase-space grid using a discontinuous Galerkin finite element approach, yielding a high-order, conservative numerical scheme that is free from the Poisson noise inherent to traditional Monte-Carlo methods. A novel and flexible velocity-space mapping technique enables the efficient treatment of the wide range of energy scales characteristic of relativistic plasmas, including QED pair-production discharges, instabilities in strongly magnetized plasmas surrounding neutron stars, and relativistic magnetic reconnection. Our noise-free approach is capable of providing unique insight into plasma dynamics, enabling detailed analysis of electromagnetic emission and fine-scale phase-space structure.

💡 Research Summary

The authors present a high‑order, conservative discontinuous Galerkin (DG) scheme for directly solving the relativistic Vlasov–Maxwell system on a phase‑space grid. By introducing a flexible four‑velocity mapping u(η₁,η₂,η₃) together with its Jacobian J_u, the method can stretch the velocity grid to cover many orders of magnitude in particle energy while preserving the structure of the differential operators. The kinetic equation is multiplied by test functions and integrated by parts, yielding a weak form that contains both volume and surface integrals. Choosing the test function to be the kinetic energy (m|v|²/2 in the non‑relativistic limit or mc²γ in the relativistic case) guarantees that the surface contributions telescope and cancel when summed over all cells, thereby ensuring exact discrete energy conservation in the semi‑discrete (time‑continuous) limit.

A continuous sub‑space Wₚ⁰,h is defined so that the discrete Lorentz factor γ_h and the relativistic current J_h can be represented with the same basis functions. The current is expressed through the transformed gradient operator ∇_U, which naturally couples to Maxwell’s equations and yields the correct electromagnetic work term ∫E·J dx. The authors prove phase‑space incompressibility (∇·α=0) and L² stability for the DG discretization, and they demonstrate numerically that total energy drift remains below 10⁻⁸ over long integrations.

Implementation is carried out within the open‑source Gkeyll framework. Each phase‑space cell uses polynomial bases of at least quadratic order in velocity space, enabling the representation of non‑polynomial quantities such as γ. Time integration employs strong‑stability‑preserving Runge‑Kutta schemes, and the mapping parameters (u_min, u_max, N_u) can be tuned to resolve ultra‑relativistic particles (γ up to 10⁶) without excessive grid refinement. The code exploits GPU acceleration and local communication patterns to keep the computational cost manageable despite the high dimensionality (up to 3‑D configuration × 4‑D velocity).

Two benchmark problems illustrate the method’s advantages. First, an electric‑field screening test mimicking pair‑production discharges in pulsar polar caps shows that the DG solver reproduces the expected decay of electric‑field energy and yields smooth electric‑field spectra, whereas a comparable PIC simulation (TRISTAN‑MP) suffers from excess high‑frequency noise and spurious field growth. The particle distribution functions from the DG run display clear acceleration peaks and a later low‑velocity plasma component, matching theoretical expectations. Second, relativistic magnetic reconnection in a strongly magnetized pair plasma is simulated; the DG results capture fine structures in the particle phase space and produce clean electromagnetic emission spectra, again free of Monte‑Carlo noise.

Although the grid‑based approach requires more memory and operations than PIC, the authors argue that for problems where noise contaminates diagnostics—such as precise radiation spectra, wave‑particle interactions, or the growth of micro‑instabilities—the DG method offers a superior cost‑to‑accuracy ratio, especially on modern multi‑core and GPU architectures.

In conclusion, the paper delivers a robust, energy‑conserving DG framework for relativistic kinetic plasmas, validates it against physically relevant scenarios, and makes the implementation publicly available. Future directions include adaptive velocity‑space refinement, incorporation of collisional and radiation‑reaction operators, and scaling to full 3‑D astrophysical simulations.