Modeling of Relativistic Plasmas with a Conservative Discontinuous Galerkin Method

Abstract

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.

Figures (5)

• Figure 1: Electric-field screening by continuously injected electron--positron plasma: comparison between Gkeyll and TRISTAN-MP v2 simulations. (a) Spectrum of electric-field fluctuations at late times, highlighting an excess of short-wavelength oscillations in the TRISTAN-MP v2 simulations (dashed lines) in comparison to Gkeyll runs (solid lines). (b) Time evolution of the electric-field energy density; the inset highlights the early-time behavior. At late times, TRISTAN-MP v2 (dashed lines) exhibits growth of the electric field, whereas Gkeyll (solid lines) shows physically consistent decay. (c) Electron and positron distribution functions, showing peaks at large $u$ that indicate efficient acceleration at early times in the unshielded field, and a broad peak near $u \sim 0$ corresponding to plasma injected at late times after the field becomes marginally screened.
• Figure 2: 2x3v Gkeyll simulation of relativistic magnetic reconnection in pair plasma with magnetization parameter $\sigma = 1$. The top row shows the electric-current distribution, while the bottom rows display particle distribution functions sampled from different regions of the current layer. A reference power-law spectrum, $dN/d\gamma \propto \gamma^{-p}$ with $p = 4$, is overplotted to facilitate comparison with previous PIC studies. A Maxwell–Jüttner energy spectrum, constructed using the average plasma-frame density, bulk velocity, and temperature, is also shown to highlight the efficient production of high-energy particles. The Gkeyll energy spectrum is plotted using absolute values to suppress small gaps arising from negative values of the distribution function caused by discretization errors.
• Figure 3: Comparison of growth rates of the two-stream (a) and filamentation (b) instabilities between simulations and warm relativistic linear theory as a function of the unstable-mode wavenumber. In both cases, the fluid-frame moments of the initial distribution are $n_0 = 1.0$ (the sum of the densities of both beams) and $T_0/m_0c^2 = 0.04$. The initial beam velocities are $v_b/c = 0.99$ for the two-stream case and $v_b/c = 0.9$ for the filamentation case. Simulations are performed using both uniform and non-uniform grids in velocity space.
• Figure 4: Evolution of the two-stream instability from initial conditions (a) into the initial nonlinear saturation (b) and long-time evolution (c) in phase space. The distribution function is normalized to its peak value at each time for clarity of visualization. Panel (d) shows the evolution of the energy in the electrostatic, $E_x$, component. We observe an initially slower-growing box-scale mode, followed by the emergence and dominance of the fastest-growing mode in the system. The relative change in total energy (e) demonstrates that energy conservation is governed by the time-step size, exhibiting convergence consistent with $(\Delta t)^3$, as expected for the third-order strong-stability-preserving Runge–Kutta method used for temporal discretization. Panel (f) shows the integrated $L^2$ energy, which decreases monotonically across resolutions, indicating $L^2$ stability.
• Figure 5: Same as Fig. \ref{['fig:TS_evolution']}, but for the filamentation instability in the full three-dimensional phase space (panel (d) shows the evolution of the energy in the $B_z$ component).