Can a Grid Beat the Noise in Extreme-Plasma Simulations?

Inside a pulsar's magnetosphere, a rotation-induced electric field pulls electrons from the neutron star surface and accelerates them to ultra-relativistic energies. Those electrons emit high-energy curvature photons that convert into electron–positron pairs in the strong magnetic field. The resulting pair plasma loads and screens the accelerating electric field, while continuous pair injection can excite coherent plasma waves—some electromagnetic, a possible route to pulsar radio emission. Simulating this cleanly is hard because the traditional particle-in-cell, or PIC, method uses a finite number of macroparticles, and finite sampling brings Poisson noise, the statistical jitter that can contaminate the weak electromagnetic signals researchers want to study.

A new arXiv preprint by James Juno and colleagues presents a grid-based solver for the relativistic Vlasov–Maxwell equations—the equations for a collisionless plasma of charged particles moving near light speed. The authors report that their method is conservative, free of the Poisson noise inherent to Monte Carlo/PIC methods, and flexible enough for the wide energy ranges of neutron-star magnetospheres and high-energy-density laboratory plasmas.

Poisson noise in pair-discharge simulations

Traditional PIC evolves macroparticles and samples them to estimate charge and current. It has underpinned work on extreme astrophysical and laser-plasma systems, but a finite sample always carries Poisson noise. Adding more particles reduces that noise only slowly, so weak signals can be buried. In pulsar polar caps, pair production loads the plasma and screens the parallel electric field; continuous pair injection excites plasma waves, some of which can become coherent radio emission. As pair density rises, PIC noise can contaminate the electromagnetic emission signal.

A mapped phase-space grid instead of macroparticles

The new solver directly evolves the particle distribution function on a grid in phase space—the combined space of position and velocity. The scheme is a discontinuous Galerkin finite element method: it represents the distribution as piecewise polynomials and is designed to conserve energy. A key innovation is a mapped velocity grid. Relativistic plasmas span a huge range of energies, up to Lorentz factors of thousands in the pair-discharge test; a uniform velocity grid would need enormous resolution to cover that range. The mapping stretches and compresses velocity coordinates so the grid can cover a wide dynamic range efficiently. Because there are no random macroparticles, the method avoids Poisson noise. That lets it resolve fine phase-space structure and weak electromagnetic fluctuations that PIC simulations can lose.

Finding 1: Screening a pulsar-like electric field

The authors test the solver on a benchmark inspired by pulsar polar caps: an electron–positron plasma in a strong electric field with fresh pairs injected continuously. They compare Gkeyll, their open-source simulation framework, with the established PIC code TRISTAN-MP v2. In the PIC runs, the late-time electric-field spectrum has an excess of short-wavelength oscillations, and the electric-field energy grows at late times. In Gkeyll, those small-scale modes are absent, and the electric field decays in a physically consistent way. The grid-based run also produces a smooth particle distribution over four orders of magnitude in four-velocity. Its cost is comparable to a PIC run with a fiducial 1,000 particles per cell. The practical payoff, the authors write, is that parameter scans can avoid limiting the pair multiplicity to keep PIC noise under control.

Finding 2: Reconnection's nonthermal spectra, resolved locally

The second test is relativistic magnetic reconnection, in which magnetic field lines rearrange into a lower-energy state and release magnetic energy to the plasma. At moderate to high magnetization—when magnetic energy density is comparable to or greater than the plasma rest-mass energy density—reconnection produces substantial nonthermal particle acceleration. The Gkeyll run uses two spatial dimensions and three velocity dimensions (2x3v), a magnetization parameter σ=1, and a small out-of-plane guide field of 0.01 B0. It recovers the expected power-law spectrum dN/dγ ∝ γ^-4: the number of particles per Lorentz factor falls as the fourth power of γ. Crucially, it recovers that spectrum in local spatial subvolumes rather than by integrating over the whole simulation box to accumulate statistics. That enables spatially resolved phase-space diagnostics, which PIC studies often cannot afford because they typically rely on global volume integration.

What the conservative DG solver does not yet do

This is a numerical methods advance, not a new telescope observation. Its benchmarks are simulations compared with linear theory and with another code. The authors are explicit about limitations. The solver does not enforce exact charge conservation or strict positivity of the distribution function, though they report negligible violations in these tests that do not affect energy conservation or stability. The pair-discharge simulations also neglect photon propagation and pair conversion, effects that future versions could include. The current velocity grid assumes independent mappings in each velocity coordinate. Future work may use spherical coordinates in velocity space, which could naturally resolve total energy, or coordinates aligned with the local magnetic field. The authors also see a path toward a general-relativistic Vlasov solver for compact objects. The paper is a preprint and has not yet been through peer review.

From pulsar screening to magnetar radio emission

The new grid-based approach is complementary to particle simulations: rather than adding more and more particles, it removes sampling noise. If it holds up in broader use, it could help researchers trace how pulsars and magnetars generate radio emission, how reconnection turns magnetic energy into energetic particles, and how collisionless relativistic plasmas behave where the usual slow, warm-plasma approximations no longer apply. The authors note that the solver is already being used in work on a potential radio emission mechanism from quiescent magnetars and that longer studies of pair-production discharges are underway.