Stable discrete representation of relativistically drifting plasmas

Describing complex physics beyond analytical theories requires numerical modeling of the underlying equations in discrete space. In plasma physics, astrophysics, or accelerator physics, Particle-In-Cell^1–4 (PIC) methods are commonly used to self-consistently solve the electromagnetic interaction of particle ensembles. A PIC algorithm iteratively solves Maxwell's equations on a discrete grid with particles following the equations of motion in a continuous space.

Some of the physical systems accessible with the PIC method feature plasmas drifting at relativistic velocities, for example, when modeling plasma-based particle accelerators⁵ in the optimal frame of reference⁶ or astrophysical plasma interactions.⁷ In those cases, the applicability of the to-date electromagnetic PIC algorithms is fundamentally limited by the numerical Cherenkov instability^8–11 (NCI), which either falsifies the numerical results or causes virulent growth of unphysical waves.

Here, we present a novel discrete formulation of the fundamental kinetic equations of plasmas, i.e., the Maxwell's and Newton-Lorentz equations, that represents the physics in a moving Galilean frame of reference and thereby is intrinsically free of the NCI for plasmas drifting at uniform relativistic velocities.

The NCI originates from the coupling of distorted electromagnetic modes with spurious particle modes. Distortions of the electromagnetic field modes are caused by numerical inaccuracies of the discretized field-solving algorithm. Spurious spatial and temporal aliases of the physical particle modes result from the numerical mismatch of sampling the continuously distributed particle quantities to the discrete field grid. To the first order, for example, numerical Cherenkov radiation (NCR) can occur, if the dispersion relation of the electromagnetic waves is numerically distorted. In this case, particles moving at relativistic velocities v_p ≈ c couple resonantly to electromagnetic waves of high frequency, which propagate at a spurious phase velocity v_Φ < v_p < c, causing Cherenkov-like radiation to be emitted. Although many algorithms, such as pseudo-spectral solvers,¹² do not suffer from NCR, higher order NCI effects severely limit the stable modeling of relativistic plasmas.

So far, no electromagnetic, fully explicit PIC algorithm is intrinsically free of NCI, even for the simple case of a plasma drifting at a uniform relativistic velocity. Previously developed suppression strategies can limit the NCI growth rate, thereby retaining the physical meaningfulness of a simulation. For example, wide-band smoothing^13–15 or damping¹⁶ of the currents or electromagnetic fields can hinder the development of the instability. Coupling of unphysical modes can be mitigated by slightly changing the ratio of the electric and magnetic fields as seen by the particles,^17–19 by scaling the deposited currents with a frequency-dependent factor,^20,21 or by artificially modifying the physical electromagnetic dispersion relation.^22–24 Yet, all of these techniques rely on numerical methods that potentially alter the physics and could affect the results obtained with the algorithm.

In contrast, the method presented in this paper inherently eliminates the NCI for a relativistically drifting plasma, as opposed to suppressing its growth by the measures described above. From a heuristic point of view, the main difference between modeling a plasma at rest, showing no NCI, and a relativistically drifting plasma is that the particles move with respect to the static numerical grid. Thus, intuitively, by mathematically representing the underlying discrete equations such that this discrepancy in relative movement is eliminated, the NCI should be suppressed.

This is achieved by applying a Galilean coordinate transformation of the form

to the frame of reference in which a plasma is moving at a relativistic velocity. Consequently, the equations of motion and Maxwell's equations transform to

and the continuity equation becomes $(∂t−vgal·∇′)ρ+∇′·j=0$⁠. Here, $∇′$ denotes the spatial derivative with respect to the Galilean coordinates $x′$⁠. For $vgal=0$⁠, these equations reduce to their well-known original form.

Using the Pseudo-Spectral Analytical Time Domain¹² (PSATD) framework, the last two equations are transformed to Fourier space and can then be integrated analytically in time. As the quantities are only known at discrete times in a PIC algorithm, the time evolution of ρ and j needs to be explicitly taken into account during integration. Typically, the currents are assumed to be constant over one time step Δt in the original coordinates x. A key difference of our new scheme is that we assume the currents to be co-moving with respect to the original coordinates x, and hence constant over one time step in the Galilean coordinates $x′$⁠. The resulting Galilean-PSATD equations for the advance of the spectral field components, $Ê$ and $B̂$⁠, from time step nΔt to $(n+1)Δt$ are then given by (see Ref. 25 for a derivation)

where k is the wavevector. The currents $Ĵ$ at time $(n+1/2)Δt$ and the charge density $ρ̂$ at time nΔt and $(n+1)Δt$ are generated by the particles and deposited to the grid nodes before being transformed to Fourier space. Subsequently, the updated fields are transformed back to real space and interpolated to the particles, which are then advanced in time using the Galilean transformed equations of motion.

This algorithm allows to model a plasma moving at $vp$ in a co-propagating set of coordinates $x′$ with $vgal=vp$⁠. As shown in Fig. 1, the flowing plasma particles now remain static with respect to the numerical grid. Because of this and the co-moving current assumption, the NCI is completely eliminated for particles streaming at the velocity $vp$⁠.

The algorithm is implemented in the WARP²⁶ code, for Cartesian coordinates, as well as in the recently developed quasi-cylindrical²⁷ code FBPIC.²⁸ In the study of Lehe et al.,²⁵ we also present an analytical derivation of the dispersion relation and conduct a detailed empirical and theoretical stability analysis for uniform relativistically flowing plasmas. Here, we restrict ourselves to presenting the general concept and the practical demonstration of the stability and accuracy of our new method with a direct application. In the following, Lorentz-boosted frame simulations of plasma acceleration with FBPIC are presented.

Plasma-based accelerators⁵ can sustain high field gradients, allowing for the acceleration of charged particles within distances shorter by orders of magnitude compared to conventional accelerators. In a plasma-wakefield accelerator, an intense driver beam (a high intensity laser pulse or particle bunch) propagates through an underdense plasma and induces a charge separation on the sub-mm scale. This leads to the excitation of a trailing density wave carrying large electric fields, suitable for the acceleration of electron bunches to high energies.

The natural frame of reference for PIC simulations of plasma accelerators is the laboratory frame. In this frame of reference, the physical objects of small scale, i.e., the laser or particle beam, propagate at relativistic velocities in a single direction while interacting with a large scale object that is static, i.e., the plasma. A Lorentz transformation in the propagation direction of the driver beam then relaxes the requirements on the spatial resolution while contracting the required simulation distance. In this Lorentz-boosted frame,⁶ the co-propagating quantities, e.g., the laser or the plasma wavelength, are elongated by γ(1 + β), whereas the previously static lengths, such as the plasma, are contracted by γ and counter-propagate with a relativistic velocity –βc. Thereby, a speed-up by orders of magnitude can be achieved that scales as $∝γboost2$⁠, with the maximum speed-up typically limited to $≈2γwake2$⁠, i.e., the phase velocity of the wake, in the case of laser-plasma acceleration.

In the following, we show simulations of a non-linear plasma wave driven by a Gaussian laser pulse with wavelength λ = 800 nm, peak normalized vector potential a₀ = 1.5, pulse length cτ = 8 µm, and waist w₀ = 30 µm (cτ and w₀ correspond to 2σ of the longitudinal and transverse laser intensity distribution, respectively) that propagates through a matched plasma guiding channel with an on-axis electron density $ne=1·1018 cm−3$⁠. In the generated wakefields, a 1 pC electron bunch of size σ_z = 1 µm, σ_r = 2 µm, and normalized emittance ϵ_n = 0.5 mm mrad, located at the back of the first wave bucket, is accelerated from 100 MeV to 687 MeV within a propagation distance of z_prop ≈ 14.3 mm. The resolution of the simulation is 32 cells per λ in the longitudinal and 2 cells per µm in the transverse direction. Third order particle shapes are used with 24 particles per cell. The time step is set to Δt = Δz/c. Please note that the simulation box and the physical quantities in the propagation direction are equally elongated in the Lorentz-boosted frame. Therefore, the number of cells in the box, as well as the longitudinal resolution with respect to the laser wavelength λ, is the same in the lab and the Lorentz-boosted frame simulation.

As described above, the occurrence of the NCI, caused by the counter-streaming relativistic plasma, can hinder the application of the Lorentz-boosted frame method for simulations of plasma-wakefield accelerators. With our new method, however, such a simulation is modeled in a Galilean transformed coordinate system that counter-propagates to the Lorentz-boosted frame with the velocity v_gal = –βc in the direction of the boosted plasma. With respect to the numerical grid, the background plasma is thus static, whereas the elongated quantities, such as the laser pulse and the electron bunch, propagate with a velocity increased by the same amount with respect to the grid.

Fig. 2 shows the charge density obtained in a Lorentz boosted frame (γ_boost = 13) with boosted longitudinal coordinate z_boost = γ(z – vt). The upper-half of the plot shows the results of a simulation with the Galilean-PSATD solver, whereas the lower-half shows the corresponding results of a standard PSATD simulation. Here, a fast growing NCI can be observed. In contrast, the same simulation remains completely stable when modeled in the Galilean transformed discrete space. We emphasize that all numerical parameters are the same in these simulations, except for the difference in using v_gal = –βc instead of v_gal = 0 for the Galilean-PSATD equations. Thus, the absence of the instability results solely from the Galilean transformation of the underlying discrete equations. Even though the electron bunch and the grid move in opposite directions, we do not observe any NCI around the bunch. This can be explained by the fact that the electron bunch has a density that is much lower than the plasma, as it is elongated in the Lorentz-boosted frame. Moreover, due to its non-zero charge, it is probably much less affected by higher order numerical Cherenkov effects. Likewise, in a laboratory frame simulation, a relativistic electron bunch does typically not lead to an instability, as long as NCR is suppressed.²⁹ 

In order to validate the accuracy of our new method, results from the stable Lorentz-boosted frame simulation are compared with a laboratory frame simulation. Fig. 3 shows the electric fields at the end of the acceleration distance. The upper half of the plots shows the results of the reference simulation (γ_boost = 1), whereas the lower half shows the corresponding back-transformed results of the Lorentz-boosted frame simulation (γ_boost = 13). Both the longitudinal fields E_z and the transverse fields E_y show no differences. Note that the results in the Lorentz-boosted frame are obtained within only a few thousand time steps, whereas the lab frame simulation takes more than half a million time steps to complete. We achieve a speed-up of ≈287 (≈92% of the optimal speed-up) with FBPIC, where the only overhead is the on-the-fly back-transformation of data to the laboratory frame.

Furthermore, we compare characteristic bunch and laser parameters to demonstrate that the physics is preserved in the Lorentz-boosted frame. Fig. 4 shows the evolution of the laser waist w₀ and the pulse duration τ, as well as the kinetic energy E_kin, the rms energy spread σ_E, and the normalized emittance ϵ_n of the accelerated electron bunch. During the propagation through the plasma guiding channel, the laser pulse self-focuses transversely and the pulse duration shortens due to the relativistic interaction with the plasma. The electron bunch is initially situated at the minimum of the accelerating field and slips towards the laser pulse during the propagation. It is accelerated to 687 MeV, while accumulating an rms energy spread of ≈11.5%, due to the slope of the accelerating field. As the bunch enters the plasma, strong transverse fields act on it abruptly, causing transverse oscillations of the bunch size and growth of the emittance ϵ_n to around 1.4 mm mrad. In direct comparison with the laboratory simulation, all the quantities shown differ only on the sub-percent level at the end of the propagation distance, which resembles a remarkable precision.

In conclusion, we have proposed a novel discrete formulation of the fundamental kinetic equations of plasmas in Galilean transformed coordinates. To the best of our knowledge, we thereby derived the first electromagnetic, fully explicit and self-consistent PIC representation that is intrinsically free of the NCI for plasmas flowing at a uniform velocity. Our concept is not reliant on otherwise inevitable numerical corrections and, unlike most of the previous NCI suppression strategies, it is independent of the specific geometry. This allows us to combine the accuracy and efficiency of a spectral, quasi-3D PIC algorithm with the superior stability properties of the presented method. Applying the Galilean scheme to simulations of plasma accelerators in the Lorentz-boosted frame yields excellent agreement, while achieving a close-to-optimal speed-up of more than two orders of magnitude in practice.

The Galilean-PSATD solver satisfies a combination of specific requirements that cause the complete absence of NCI for a plasma drifting at uniform relativistic velocity, which can be summarized qualitatively as follows: It is free of spurious numerical dispersion; thus, first order NCR effects are suppressed. Furthermore, the analytic integration in the time domain in combination with the co-moving current assumption explicitly retains stability properties that are otherwise lost when directly discretizing the Maxwell equations. Finally, the suppressed movement of the particles with respect to the numerical grid eliminates aliasing effects due to sampling the particle fields to the grid and vice versa. The Galilean scheme is therefore not directly applicable to other solvers, such as standard Finite-Difference Time Domain (FDTD) solvers. It is unclear whether it would be possible to derive an implementation using other solvers (e.g., FDTD) that would retain the advantages obtained when integrating with the Galilean-PSATD solver. This requires further studies.

Nevertheless, it is possible to apply arbitrary-order spectral methods^4,24,30 to the Galilean-PSATD solver, which enables the efficient parallelization based on spatial domain decomposition³¹ employing domain-local Fourier transforms of the fields. These topics will be covered by future research, as well as the potential generalization to support arbitrary relativistic plasma flows. For example, the new method could directly be extended to model collisionless astrophysical shocks⁷ involving two plasmas, by employing separate numerical grids for each plasma using different Galilean transformed coordinates. Taking advantage of the superposition principle, only the electromagnetic fields would be shared between those individual grids.

We gratefully acknowledge the computing time provided on the supercomputer JURECA under project HHH20 and on the PHYSnet cluster of the University of Hamburg. For verification of the method in Cartesian geometry with Warp, this research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Work at LBNL was funded by the Director, Office of Science, Office of High Energy Physics, U.S. Department of Energy under Contract No. DE-AC02-05CH11231, including the Laboratory Directed Research and Development (LDRD) funding from Berkeley Lab.