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.
In a diverse array of plasma environments, from the magnetospheres of pulsars (Philippov and Kramer 2022) to runaway electrons in fusion reactors (Breizman et al. 2019) to ultra-intense laser-plasma experiments (Gonoskov et al. 2022), we observe plasmas which are not only hot enough and diffuse enough to be collisionless, i.e. described by kinetic theory, but also relativistic. These plasmas are some combination of relativistically hot, with temperatures exceeding the rest mass energy of the particular particles, T s > m s c 2 , and moving at relativistic velocities, v ∼ c where v is either a bulk velocity of the plasma or a particular particle velocity. Modeling these plasmas thus requires solving the Vlasov, or Boltzmann, equation for the evolution of the particle distribution function with the necessary modifications to include relativity,
where u is the spatial component of the four-velocity, u = γv, γ = 1 + |u| 2 /c 2 is the particle Lorentz boost factor, q s and m s are the particle charge and rest mass respectively, E and B are the electric and magnetic fields respectively, C[f ] includes all discrete particle effects such as collisions and radiation reaction, and S is an explicit source term due to, e.g., pair production. Eq. ( 1), for any number of plasma species of interest, is then coupled to Maxwell’s equations via (four-) velocity moments of the distribution function to obtain the charge and current density used to evolve the electric and magnetic fields and close the system of equations.
Traditionally, the numerical solution of Eq. ( 1) is done with Monte Carlo methods due to the high dimensionality of the Vlasov equation, reducing the six dimensional solve to evolving “macroparticles” in (up to) three spatial dimensions, and sampling these particles to construct the distribution function and deposit charges and currents on a grid to couple to a solver for Maxwell’s equations-typically referred to as the particle-in-cell, or PIC method. This approach has proved historically fruitful, revealing new insights into a diverse array of plasma phenomena, and there now exist a number of production PIC codes for modeling extreme astrophysical and laser-plasma systems (Fonseca et al. 2002;Spitkovsky 2005;Fonseca et al. 2008;Bowers et al. 2009;Cerutti et al. 2013;Germaschewski et al. 2016;Parfrey et al. 2019;Hakobyan et al. 2024;Hakobyan et al. 2025;Galishnikova et al. 2025). However, Monte Carlo methods inevitably suffer from Poisson noise due to the statistics from using a finite number of particles. Unfortunately, this noise decreases only very slowly, ∝ N ppc , where N ppc is the number of particles per grid cell (Birdsall and Langdon 1990), and thus in some applications where we are interested in precise phase space dynamics and detailed diagnostics of electromagnetic emission, the standard PIC technique can be inadequate for understanding the plasma’s evolution.
The alternative, directly discretizing Eq. ( 1) on a phase space grid and numerically integrating the high dimensional partial differential equation for the physical system of interest, would eliminate sources of error from Monte Carlo methods such as noise, but at increased computational cost. This increased computational cost has proved a significant deterrent within the field of plasma physics, especially for the relativistic variant of the kinetic equation. As computational capabilities have advanced with new hardware and the development of novel algorithms, a number of production codes for both the continuum discretization of the multi-species Vlasov-Maxwell system equations (Vencels et al. 2016;Juno et al. 2018;Roytershteyn and Delzanno 2018;Hakim and Juno 2020;Allmann-Rahn et al. 2022;Grünwald et al. 2025) and hybrid-kinetic reductions relevant when the electron kinetics can be ignored (Valentini et al. 2007;Kempf et al. 2013;Palmroth et al. 2018) have been developed, but only for low energy density space, astrophysical, and laboratory plasmas. While extensions of numerical methods for kinetic equations such as the Lattice Boltzmann method have been applied to relativistic neutral Boltzmann solvers (Ambruş et al. 2022) and reduced geometry relativistic Vlasov-Maxwell solvers for targeted applications (Wettervik et al. 2017;Ye and Chen 2025), extending the bevy of algorithmic developments to relativistic plasmas more generally has proved nontrivial for two principal reasons.
Firstly, significant breakthroughs in robustness were obtained with conservative grid-based methods for the Vlasov equation, but these breakthroughs relied on particular choices of basis functions in finite element (Cheng et al. 2014;Juno et al. 2018;Kormann et al. 2025) or spectral (Delzanno 2015;Koshkarov et al. 2021;Pagliantini et al. 2023) methods that allowed one to prove conservation relations such as conservation of energy were implicitly contained in the spatial discretization of the Vlasov equation, and there exists no obvious basis choice when we add
This content is AI-processed based on open access ArXiv data.