Numerical instability due to relativistic plasma drift in EM-PIC simulations
The numerical instability observed in the Electromagnetic-Particle-in-cell (EM-PIC) simulations with a plasma drifting with relativistic velocities is studied using both theory and computer simulations. We derive the numerical dispersion relation for a cold plasma drifting with a relativistic velocity and find an instability attributed to the coupling between the beam modes of the drifting plasma and the electromagnetic modes in the system. The characteristic pattern of the instability in Fourier space for various simulation setups and Maxwell Equation solvers are explored by solving the corresponding numerical dispersion relations. Furthermore, based upon these characteristic patterns we derive an asymptotic expression for the instability growth rate. The asymptotic expression greatly speeds up the calculation of instability growth rate and makes the parameter scan for minimal growth rate feasible even for full three dimensions. The results are compared against simulation results and good agreement is found. These results can be used as a guide to develop possible approaches to mitigate the instability. We examine the use of a spectral solver and show that such a solver when combined with a low pass filter with a cutoff value of $|\vec{k}|$ essentially eliminates the instability while not modifying modes of physical interest. The use of spectral solver also provides minimal errors to electromagnetic modes in the lowest Brillouin zones.
💡 Research Summary
This paper addresses a long‑standing problem in electromagnetic particle‑in‑cell (EM‑PIC) simulations: the emergence of a numerical instability when a plasma drifts at relativistic velocities. The authors begin by deriving a general numerical dispersion relation for a cold plasma moving with a drift velocity β = v/c. By discretizing both Maxwell’s equations and the particle push, they obtain a relation that couples the beam mode of the drifting plasma (ω ≈ k·v) with the discrete electromagnetic (EM) modes of the grid (ω² = c²k² + ω_p²). The coupling produces complex frequency solutions, i.e., a growth rate γ > 0, whenever the beam and EM branches intersect in k‑space.
The analysis shows that the instability is most pronounced near the edges of the first Brillouin zone, and its strength depends sensitively on the grid spacing Δx, the time step Δt, and the Courant number S = cΔt/Δx. Larger S values shift the EM dispersion upward, increasing the overlap with the beam branch and thereby amplifying γ. The authors systematically explore 1‑D, 2‑D, and 3‑D configurations, confirming that the same mechanism operates in all dimensions.
To make the evaluation of γ practical for large‑scale simulations, they derive an asymptotic expression:
γ ≈ (π/2)·(ω_p/γ_b)·(Δt/Δx)·|sin(kΔx/2)|·|1 − S·cos(kΔx/2)|,
where γ_b is the Lorentz factor of the drifting plasma. This formula captures the dependence on wave‑number k and grid parameters and can be evaluated orders of magnitude faster than solving the full numerical dispersion relation. Using this expression, the authors perform exhaustive three‑dimensional parameter scans, identifying optimal Δt/Δx ratios and Courant numbers that minimize the growth rate while preserving physical fidelity.
The paper then compares several Maxwell solvers. Traditional Yee‑FDTD and higher‑order staggered schemes exhibit the predicted instability across a broad range of k. In contrast, a pseudo‑spectral (FFT‑based) solver reproduces the continuous EM dispersion with negligible phase error, dramatically reducing the overlap with the beam mode. Nevertheless, even the spectral solver can admit high‑k noise. To eliminate this residual source, the authors apply a low‑pass filter that removes all Fourier components with |k| > k_c, where k_c is chosen just below the first Brillouin‑zone boundary. Simulations demonstrate that this combination (spectral solver + low‑pass filter) suppresses the instability to machine‑precision levels without affecting physically relevant low‑k modes such as laser‑plasma interactions or wakefield generation.
Extensive validation is performed using both 2‑D and 3‑D relativistic drift scenarios. Measured growth rates match the asymptotic formula and the full numerical dispersion solution within a few percent. Moreover, the filtered spectral solver yields EM wave phase and group velocities that agree with analytical values to better than 10⁻⁴, confirming that the mitigation strategy does not compromise accuracy.
Finally, the authors provide practical guidelines: (1) compute the expected growth rate using the asymptotic expression as soon as drift parameters are known; (2) if γ exceeds a tolerable threshold, reduce the Courant number or adjust Δt/Δx; (3) whenever possible, adopt a pseudo‑spectral Maxwell solver; (4) apply a low‑pass filter with cutoff just below the Brillouin‑zone edge to eliminate remaining high‑k coupling. Following these steps enables stable, accurate EM‑PIC simulations of relativistically drifting plasmas, opening the door to reliable modeling of high‑energy density physics, laser‑wakefield acceleration, and astrophysical jet dynamics.