Stable discrete representation of relativistically drifting plasmas

Representing the electrodynamics of relativistically drifting particle ensembles in discrete, co-propagating Galilean coordinates enables the derivation of a Particle-in-Cell algorithm that is intrinsically free of the Numerical Cherenkov Instability, for plasmas flowing at a uniform velocity. Application of the method is shown by modeling plasma accelerators in a Lorentz-transformed optimal frame of reference.

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 (PIC) methods are commonly used to self-consistently solve the electromagnetic interaction of particle ensembles [1][2][3][4]. 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 [5] in the optimal frame of reference [6] or astrophysical plasma interactions [7]. In those cases, the applicability of the to-date electromagnetic PIC algorithms is fundamentally limited by the Numerical Cherenkov Instability (NCI) [8][9][10][11], 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. Maxwell's and the 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 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 Cherenkovlike radiation to be emitted. Although many algorithms, such as pseudo-spectral solvers [12], 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][14][15] or damping [16] 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][18][19], by scaling the deposited currents with a frequency-dependent factor [20,21], or by artificially modifying the physical electromagnetic dispersion relation [22][23][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 arXiv:1608.00215v1 [physics.plasm-ph] 31 Jul 2016 and the continuity equation becomes (∂ t − v gal · ∇ ) ρ + ∇ · j = 0. Here, ∇ denotes the spatial derivative with respect to the Galilean coordinates x . For v gal = 0, these equations reduce to their well-known original form.
Using the Pseudo-Spectral Analytical Time Domain (PSATD) framework [12], 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, hence constant over one time step in the Galilean coordinates x . The resulting Galilean-PSATD equations for the advance of the spectral field components,Ê andB, from time step n∆t to (n + 1)∆t are then given by (see [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 v p in a co-propagating set of coordinates x with v gal = v p . 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 v p .
FIG. 1. Schematic drawing illustrating the Galilean concept. Without applying a Galilean coordinate transformation to the Particle-In-Cell equations (Standard ), a plasma flowing with velocity vp in z (represented by a single particle) would propagate a distance vp∆t with respect to the numerical grid (represented by a single cell) during one time step ∆t. However, in a Galilean transformed discrete space x with vgal = (0, 0, v gal = vp) the plasma particles remain static with respect to the discrete grid nodes, which themselves propagate a distance z + v gal ∆t in the original coordinate system x.
The algorithm is implemented in the Warp code [26], for Cartesian coordinates, as well as in the recently developed quasi-cylindrical [27] code Fbpic [28]. In [25] 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 [5] 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 [6]. 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 counterpropagate with a relativistic velocity −βc. Thereby, a speed-up by orders of magnitude can be achieved that scales as ∝ γ 2 boost , with the maximum speed-up typically limited to ≈ 2γ 2 wake , 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 laser pulse with wavelength λ = 800 nm, peak normalized vector potential a 0 = 1.5, pulse length cτ = 8 µm and waist w 0 = 30 µm that propagates through a matched plasma guiding channel with an on-axis electron density n e = 1 · 10 18 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 40 cells per µm 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.
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 counterpropagates 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. The lower half shows the same simulation, mirrored along the x = 0 axis, but conducted with the standard PSATD solver. Here, a fast growing, virulent NCI can be observed. 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 Lorentzboosted 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 [29]. In order to validate the accuracy of our new method, results from the stable Lorentz-boosted frame simulation are compared to 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 Lorentzboosted 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 backtransformation 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 0 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, caus-ing 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 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 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.
Future research will cover the applicability of the Galilean scheme to other solvers, the parallelization based on domain decomposition [30] with arbitrary-order spectral methods [4,24,31] and the potential generalization to support arbitrary relativistic plasma flows. For example, the new method could directly be extended to model collisionless astrophysical shocks [7] 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. Work at LBNL was funded by the Director, Office of Science, Office of High Energy Physics, U.S. Dept. of Energy under Contract No. DE-AC02-05CH11231, including from the Laboratory Directed Research and Development (LDRD) funding from Berkeley Lab.