Superluminal phase velocity approach for suppression of Numerical Cherenkov Instability in Maxwell solver

A modification of a standard finite-difference time-domain (FDTD) Maxwell solver is proposed based on the inclusion of the terms of the higher order in the finite differences approximation. By selecting the coefficients of this hybrid scheme the numeric dispersion of the electromagnetic waves can be fine-tuned. It is shown that a particular choice of the coefficients leads to the phase velocity of all the waves being slightly larger than speed of light. Results of the 3D Particles-In-Cells (PIC) simulations demonstrate that the proposed scheme of the Maxwell solver is efficient to mitigate the Numerical Cherenkov Instability (NCI) when simulating ultrarelativistic charged beams.


Introduction
The upcoming FACET-II facility [1] is expected to deliver extremely bright electron and positron beams. This will bring the ability to conduct unprecedented experiments on beam interaction with solids, gamma generation, probing strong field quantum electrodynamics (QED) and more. To guide these future experiments a rigorous search for the optimal configurations should be carried out. While analytical examination of the mentioned processes can give a hint what these configurations should look like, the numeric simulations are expected to provide the quantitative estimations. The state-of-the-art method for coupled modeling of plasmas and electromagnetic fields is the Particles-In-Cells (PIC) method. Although there is a general problem of the scheme usually used to solve the Maxwell's equations -FDTD. The problem is that the dispersion of the waves which exist in this scheme differs from the vacuum dispersion ω = kc. In particular, there exist waves which phase velocity is smaller than the speed of light. Thus the ultrarelativistic charged particles can satisfy the Cherenkov synchronism condition and resonantly excite such waves in vacuum -the effect called the Numerical Cherenkov Instability (NCI). Several ways of mitigating this effect were proposed, such as using high-order approximation of the spatial derivatives [2,3,4], applying filters [5], or resorting to the spectral approach and solving the Maxwell's equation in the Fourier domain [6], which all have their own advantages and use cases. We propose another scheme of mitigating the NCI which can be easily implemented in the existing PIC codes and which is mostly suitable for the simulation of the ultrarelativistic beams interacting with plasmas. The key feature of our scheme is that in contrary to the common FDTD scheme, all the waves propagate with the superluminal velocity.

Hybrid scheme development
We propose the following alteration of the usual FDTD scheme are the finite-difference operator and averaging operator correspondingly. Here α is one of x, y, z and j α are the nodes of the grid.
Seeking for the solution of the Eqs.
(1) -(6) in the form E, B ∝ e −iωt+ikr we get the following dispersion relation where A α = 1 ∆α sin kα∆α 2 (a 0α + a 1α cos k α ∆α). Analysis of this relation shows that increase of the coefficients a 1α leads to the increase of the magnitude of the phase velocity of the waves propagating along the coordinate axes v ph,α = ω/k α for all values of the k α ∆ α in the fundamental Brillouin zone |k α ∆α ≤ π|. Hence we can choose the coefficients a 1α in such a way that the phase velocity of the waves is larger than the speed of light inside and is equal to the speed of light on the boundaries of the fundamental Brillouin zone. The coefficients that satisfy this condition are found as follows In that case the maximum value of the A α,max = sin (π∆t/2∆α) /∆t and thus the stability condition of proposed scheme is the following sin 2 π 2 ∆t ∆x + sin 2 π 2 ∆t ∆y + sin 2 π 2 ∆t ∆z < 1.
We compared the dispersion relations of proposed hybrid scheme with the typical schemes used in PIC codes: FDTD and NDF [7]. The results of the comparison are shown in the Fig. 1. As described above in the proposed scheme the phase velocity of all waves in the fundamental Brillouin zone is greater than the speed of light, in contrast to the FDTD scheme, where all the waves propagate slower than the speed of light, and the NDFX scheme, where only the waves propagating along one coordinate axis (in this case, the x-axis), have a phase velocity equal to the speed of light. The difference between the group velocity and the speed of light, which is present in all schemes, is of the same order of magnitude for all three schemes.

PIC tests
We also compared the results of the identical numerical simulations of the propagation of an ultrarelativistic electron beam in vacuum using different schemes. The comparison results are shown in Fig. 2. The NCI grows most rapidly in the FDTD scheme. In the NDFX scheme, despite the exact vacuum dispersion of the waves propagating along the axis of the beam propagation, instability is also present. It occurs firstly due to excitation of the waves propagating at a small angle to the axis which have phase velocity smaller than the speed of light and secondly due to the intersection of the dispersion of the beam and electromagnetic waves outside the fundamental Brillouin zone (the so-called frequency aliasing [8]). Our scheme gets rid of the first factor completely as all the waves are superluminal. As for the aliasing in general the growth rate of the NCI decreases with increase of the order of the Brillouin zone, in which the beam and the wave dispersions intersect. As seen from the results of the PIC simulations the impact of the aliasing is not present in our scheme for long enough time (> 0.12 ps). The simulations were performed using the QUILL code [9]. The simulation the cell size is 0.003µm × 0.01µm × 0.01µm. The time step ∆t was set to 0.6∆x/c for the proposed scheme and the FDTD scheme and it was set to ∆x/c for the NDFX scheme.

Disccussion and conclusion
As the proposed scheme only modifies the stencil of the FDTD scheme, it can be easily implemented in the existing PIC code, compared to, for example, the recently proposed RIP scheme [10] which requires the code to be restructured quite significantly. The deviation of the dispersion relation from the vacuum one in the proposed scheme compared to the FDTD or the NDF schemes is of the same order of magnitude and is also mostly significant in the high frequency area, so if the fields of interest are resolved well on the grid that deviation can be negligibly small. Although if one would be interested in the long-term evolution of the laser pulses the total error can be significant. Thus we argue the scheme is most suitable for simulating dense charged beams interaction with stationary targets, laser pulses or other beams expected to be experimentally realized at FACET-II facility in future, as well as for simulating laser-plasma interaction where the bunches of relativistic charged particles are produced.  Figure 2. Results of the 3D PIC simulations of the ultrarelativistic electron beam propagation in the vacuum. Top row: distribution of the y component of the electric field; bottom row: distribution of the electron density at the initial moment (a) and after 120 fs calculated with the proposed scheme (b), the NDFX scheme (c) and FDTD scheme (d).
In conclusion, we developed an easy to implement scheme for numeric solution of the Maxwell's equations. A particular choice of the parameters of the scheme allows to alter the dispersion relation of the electromagnetic waves in such a way that their phase velocity exceed the speed of light in the fundamental Brillouin zone. Thus the NCI is strongly mitigated in that scheme which is verified by the results of the 3D-PIC simulations.