Scalable spectral solver in Galilean coordinates for eliminating the numerical Cherenkov instability in particle-in-cell simulations of streaming plasmas

Discretizing Maxwell’s equations in Galilean (comoving) coordinates allows the derivation of a pseudospectral solver that eliminates the numerical Cherenkov instability for electromagnetic particle-in-cell simulations of relativistic plasmas ﬂowing at a uniform velocity. Here we generalize this solver by incorporating spatial derivatives of arbitrary order, thereby enabling efﬁcient parallelization by domain decomposition. This allows scaling of the algorithm to many distributed compute units. We derive the numerical dispersion relation of the algorithm and present a comprehensive theoretical stability analysis. The method is applied to simulations of plasma acceleration in a Lorentz-boosted frame of reference.


I. INTRODUCTION
The occurrence of the numerical Cherenkov instability (NCI) [1][2][3][4] can severely limit the applicability of the particlein-cell (PIC) method [5,6] for simulations of relativistic beams or plasmas. The instability arises from the unphysical coupling of relativistic particles to numerically distorted electromagnetic modes and aliases thereof.
To zeroth order, numerical errors in the field solver can artificially reduce the phase velocity of electromagnetic waves, which causes relativistic particles to emit numerical Cherenkov radiation (NCR). NCR is most commonly observed in finite-difference time-domain (FDTD) solvers [7], which introduce spurious dispersion when approximating the derivatives in Maxwell's equations by taking finite differences of the fields in space and time. In order to mitigate NCR, specialized FDTD solvers with improved numerical dispersion have been developed [8][9][10][11].
A more general solution is given by the class of pseudospectral time-domain (PSTD) solvers [12,13], which represent the spatial derivatives in the frequency domain. In the special case of the pseudospectral analytical time-domain (PSATD) solver [14,15], Maxwell's equations are integrated analytically in time under the assumption of constant currents. The PSATD solver correctly models the dispersion of electromagnetic waves and thereby mitigates the problem of NCR in simulations of relativistic beams [16,17]. It is, however, not * manuel.kirchen@desy.de generally free of the NCI and remains unstable for simulations of relativistic plasmas [18].
In recent work [19,20], we extended the PSATD solver to eliminate the NCI for simulations of streaming plasmas. As shown in Fig. 1, the fundamental idea is to apply a simple Galilean coordinate transformation to the PIC equations, such that a moving plasma appears to be at rest with respect to the static numerical grid. Following discretization and analytical integration under the assumption of henceforth comoving currents, the Galilean-PSATD solver becomes intrinsically free of the NCI for a relativistic plasma flowing at a uniform velocity. Unlike previous strategies to suppress the NCI [18,[21][22][23][24][25][26], the Galilean scheme does not require artificial corrections or numerical modifications of the underlying solver. Furthermore, it is applicable in both Cartesian and cylindrical geometry. Potential applications of the method include simulations of plasma accelerators in an optimal, Lorentz-boosted frame of reference [27][28][29][30] or relativistic collisionless shocks [31][32][33] in astrophysics.
However, in its original formulation [19,20], the solver is not scalable to many distributed compute units. In PIC codes, scalability is typically achieved by spatially decomposing the simulation box into smaller domains that are computed in parallel. When using FDTD solvers, continuity of the physics can be satisfied by exchanging particles and fields only at the local domain boundaries. In contrast, the global field transformations in pseudospectral solvers require, in principle, costly communication across the entire numerical grid.
However, it has been shown that domain decomposition is still possible with the PSATD solver, by using local transforms and large boundary layers [15]. In addition, by incorporating a finite-difference operator into the spectral field solver equations [34], an arbitrary-order PSATD solver can be derived. At the cost of reducing the accuracy of the spatial derivatives to that of an FDTD solver, the evolution of the fields on the grid becomes more local. Although this makes a simulation again prone to spurious dispersion and NCR, the order can be tuned to balance accuracy and locality enabling scalability to many compute units [35], while correctly modeling the physics [17,36].
In this article, we show that this parallelization strategy can also be applied to the Galilean-PSATD solver, while retaining its unique stability properties. For a beam or plasma comoving with the Galilean coordinates, the arbitrary-order Galilean-PSATD solver remains free of NCR and eliminates the NCI. We thereby remove the limitation of the original solver to efficiently scale to many compute units and lay the foundation for intrinsically stable and massively parallel PIC simulations of relativistically streaming plasmas.
In Sec. II, we derive the arbitrary-order Galilean-PSATD solver, verify its spatial locality, and discuss the vacuum dispersion relation. A theoretical stability analysis is presented in Sec. III, proving the solver's insusceptibility to the NCI for a relativistically streaming plasma, which is verified by two-dimensional (2D) Cartesian PIC simulations with WARPX [37]. Finally, in Sec. IV we apply the parallelized solver to multi-GPU Lorentz-boosted frame simulations of laserplasma acceleration using the quasi-3D PIC code FBPIC [16].

A. Derivation of the discretized Maxwell equations
As explained in previous work [19,20], in the Galilean where ∇ denotes a spatial derivative with respect to the Galilean coordinates x . In the finite-order PSATD scheme, the spatial derivatives ∇ are replaced by discretized derivativeŝ ∇ on a grid, which are accurate to order n order and correspond, in Fourier space, to a multiplication by a modified wave vector [k] [17,34]. As explained in detail in Sec. II B, we thereby introduce spatial locality to the evolution of the fields on the grid, which is the prerequisite for efficient parallelization by domain decomposition. Maxwell's equations are then transformed into Fourier space and integrated analytically over one time step. More explicitly, the representation in Fourier space of Maxwell's equations with discretized derivatives is where the Fourier transform is defined byF (k, t ) = Just like in the case of the infinite-order Galilean scheme [19,20], these equations can be integrated analytically under the assumption that the currentĴ is constant over one time step, and that the fieldsÊ,B,Ĵ , andρ satisfy the conservation equations: (3) The result of this analytical integration iŝ where we use the shorthand notationsÊ n ≡Ê (k, n t ),B n ≡ B(k, n t ), and In addition, the discretized continuity equation reads Because the current J d from a direct deposition scheme (from the particles to the grid) does not automatically satisfy the continuity equation, this current can be corrected in Fourier space:Ĵ Note that the only difference to the infinite-order Galilean-PSATD solver [19,20] is that the wave vector k has been replaced by the modified wave vector [k] of order n order , including in the definitions of the coefficients θ , ν, χ 1 , χ 2 , and χ 3 . In fact, for n order → ∞, [k] → k and Eqs. (2a) to (7b) then equal to the original formulation of the solver as in Ref. [20]. In addition, in the limit v gal = 0, the above expressions reduce to the standard PSATD equations [14,15]. On the other hand, for v gal = v plasma , the solver becomes intrinsically free of the NCI for a relativistically drifting plasma.
For the infinite-order case, a qualitative explanation for the absence of the instability was found to be threefold [19,20]: First, as the solver correctly models the electromagnetic dispersion relation, the NCI is mitigated to zeroth order (no NCR). Second, the immobility of the plasma with respect to the numerical grid suppresses the coupling of spatial aliases of the fields to the particles. Third, the analytical integration naturally takes into account the comoving currents in the direction of the plasma motion and thereby eliminates remaining sources of the instability. As to why those criteria-and in particular the first one-hold true even in the finite-order case will be explained in detail throughout this article.
Apart from its resilience against the NCI, the arbitraryorder Galilean-PSATD solver possesses several other advantageous properties. Due to the analytical integration in time, the solver is not bound to a CFL condition-in principle allowing for arbitrarily long time steps. The fields are centered in space and-except for the currents J-also in time. This avoids interpolation errors common to FDTD solvers or the simpler PSTD solver, where the E and B fields are staggered in time. Another advantage of the Galilean scheme is its independence of geometry. Analogous equations in cylindrical geometry with azimuthal mode decomposition can be readily derived from the equations listed above [20], although in this case finite-order derivatives are possible in the axial direction only. The solver has been implemented in the PIC frameworks WARP [38] and WARPX [37] in 2D and 3D Cartesian coordinates and in the quasi-3D code FBPIC [16].

B. Locality of the Galilean-PSATD solver with a finite-order stencil
To achieve locality of the field solver, we introduced the discrete spatial derivatives∇ in Eqs. (1a) and (1b), representing finite differences of arbitrary order n order . More explicitly, for a scheme with centered fields, the discrete derivative D of a function F at point p and along an axis z on the spatial grid can be written as [39] Because Eq. (8a) only takes into account the 2n ≡ n order adjacent grid points for calculating the spatial derivative, the accuracy of this operation is inversely proportional to its locality. The equivalent spectral representation of this operator can be defined as [17,34] withF the spectral transform of F and where we introduced the modified wave number where [k z ] → k z for n order → ∞. From Eq. (9) we see that the derivative of finite order n order in real space can be reproduced in spectral space by replacing the wave vector k with its modified counterpart [k]. With such a modification the simpler PSTD solver, which is only secondorder accurate in time, becomes mathematically equivalent to an arbitrary-order FDTD solver [34]. As a result, the evolution of the fields on the grid becomes strictly local. In case of the (Galilean-)PSATD solver, however, the effect of the discrete derivatives on the locality is more complex, because the modified wave vector also appears in the additional coefficients C and S and-if |v gal | = 0-also in θ , ν, χ 1 , χ 2 , and χ 3 [Eqs. (5a) to (5d)]. These additional coefficients stem from the analytical integration in time, which is equivalent to taking infinitesimal steps with the PSTD solver. In fact, it was shown that the PSTD solver can be derived from the PSATD solver by performing a first-order Taylor expansion in t [15]. As a result, the effective imprint of the (Galilean-)PSATD solver is wider than that of a PSTD (or FDTD) solver of the same order.
As shown in previous work [17], we can study the locality by applying an inverse Fourier transform to the set of coefficients for advancing the fields. By that we get an equivalent to the stencil of an FDTD solver, that is the weightings of the finite differences on the spatial grid. In cylindrical coordinates and where we use finite-order derivatives in z only, the stencil coefficients (11) can be numerically calculated for arbitrary values of the transverse wave number k ⊥ , t, and v gal = v gal u z from a sum of inverse Fourier transforms F −1 along z.
Alternatively, the stencil can be retrieved directly from a simulation. We use FBPIC to initialize a δ-peak in the E and B fields and propagate it for a single time step t = z/c. The resulting amplitude and reach of the fields across the grid is a direct measure for the stencil function. (Note that the propagation of this δ-peak should not be confused with that of a real physical signal, since it does not satisfy ∇ · E = 0 and ∇ · B = 0 at t = 0. Instead, this is a simple empirical way to obtain the stencil from a simulation.) Figure 2 compares the stencils from the numerical evaluation [Eq. (11)] and the simulation. As expected, the infiniteorder stencil affects the fields globally and therefore extents over the entire grid. Contrarily, a stencil of n order = 32 quickly falls off to machine precision at a distance N z z from the point of origin. While symmetric for v gal = 0, the stencil leans toward the direction of the grid velocity for v gal = c. Hence, the stencil becomes slightly wider but the solver remains local for |v gal | > 0. In practice, when spatially decomposing the simulation box into domains that are computed in parallel, each subgrid is extended by the stencil extent N z . These additional cells serve as a boundary layer to exchange the fields between neighboring domains at each time step, thereby ensuring continuity of the physics across the global grid. Using significantly fewer boundary cells would cause spurious distortions in the electromagnetic fields. Efficient scaling is only possible if N z is small compared to the total grid size, thus n order should be low, which decreases the accuracy. While highorder stencils are needed for some applications [36], it has been shown before [17], and we show again in Sec. IV, that relatively low orders are already sufficient to correctly model other applications such as plasma acceleration.
It should be noted that the charge conserving correction of the currents J, Eq. (7a), remains nonlocal even for a finite order. However, and in stark contrast to the E and B fields, truncating its effective stencil function at N z typically has no negative effect as it only causes very small variations in the currents. By preference, it is also possible to use a charge conserving current deposition that preserves locality. Although initially derived for 3D Cartesian coordinates only [15], this deposition scheme has recently been extended to cylindrical coordinates and will be published elsewhere. For the results in this article, we restrict ourselves to using the simple current correction from Eq. (7a).

C. Asymmetry of the numerical vacuum dispersion relation in Galilean coordinates
In the previous section we showed that the finite-order Galilean-PSATD solver is more local than its infinite-order counterpart, but its accuracy is lower (by construction) which causes the vacuum dispersion relation to be altered.
Considering Eqs. (4a) and (4b) in vacuum (i.e.,Ĵ = 0 and ρ = 0), with electromagnetic modes of the form e ik·x−iωt = e ik·x −i(ω−k·v gal )t leads, after some algebra, to the numerical dispersion relation in vacuum describing the propagation (in both directions) of electromagnetic waves on the discretized grid. Figure 3 shows the one-dimensional dispersion relation for different values of v gal in a single direction z (i.e., along k z ). As expected, in the trivial case with v gal = 0, and for analytic dispersion (n order = ∞), it is a linear relation with slope c for both forward-and backward-propagating waves. As the solver becomes more local (n order = 8, 32, 128), highfrequency waves slow down and travel at a phase velocity v = ω/k that is smaller than the vacuum speed of light.
In simulations, these slow electromagnetic modes can couple resonantly to relativistic particles, leading to the excitation of unphysical radiation known as numerical Cherenkov radiation (NCR). Frequencies supported by relativistic particles traveling at a velocity v 0 along z can be described by an artificial beam mode k z v 0 . NCR occurs at intersections of this beam mode with the distorted electromagnetic mode ω [Eq. (12)], i.e., if the main NCI resonance condition is fulfilled. Naturally, electromagnetic waves in vacuum would obey ω = ck > k z v 0 and thus never cross this beam mode. For finite stencils, however, the unphysical resonance condition is fulfilled whenever c[k] = k z v 0 . However, once the comoving velocity v gal is increased, the numerical dispersion relation of the Galilean-PSATD solver becomes asymmetric. Increasing v gal improves the dispersion of forward-propagating waves (i.e., those having ω and k z of the same sign) but worsens the dispersion of backward-propagating waves. For v gal = v 0 , there will be no intersection between the electromagnetic mode and the beam mode. Inserting the vacuum dispersion relation, Eq. (12), into the resonance condition, Eq. (13), gives c[k] − [k z ]v 0 = 0, a condition which is never met. Moreover, at the limit v gal = c, waves propagating in the direction of the comoving frame exhibit perfect physical dispersion regardless of the solver's order.
This unique property of the numerical dispersion relation is a direct consequence of the Galilean transformation and a fundamental prerequisite for the stability of the solver. For v gal = v 0 , the asymmetry of the dispersion relation effectively suppresses any NCR in the direction of the beam or plasma and-in combination with the other stability propertieseliminates the NCI even in the finite-order case, as shown in Sec. III. The dispersion of backward-propagating waves, on the other hand, worsens compared to the reference case (v gal = 0), which has to be taken into account in practice. For example, it could be necessary to use a higher order n order of the solver to mitigate problems from NCR in simulations that include a backward-propagating relativistic beam. Nevertheless, we show in Sec. IV that relatively low orders are still sufficient to not degrade such a counterpropagating beam in Lorentz-boosted frame simulations of plasma acceleration.

A. Numerical dispersion relation of a cold relativistic plasma
In Sec. II, we defined the finite-order Galilean-PSATD solver and showed that the locality of the field evolution makes it suitable for parallelization. In addition, based on the asymmetrical features of the vacuum dispersion relation, we derived a first, intuitive understanding for why the Galilean scheme retains its stability properties under such a finite-order approximation of the fields. To get the whole picture, however, we can derive a theoretical dispersion relation that includes the particle feedback of a moving plasma on the electromagnetic fields and that considers all of the approximations of the numerical scheme [2][3][4]18].
We therefore extend the derivation of the dispersion relation that we originally presented in Ref. [20] to include the modified wave vector [k] and summarize the main differences in Appendix. The resulting expressions, Eq. (14) with Eqs. (15a) to (16b) represent the numerical dispersion relation of the arbitrary-order Galilean-PSATD solver, relating ω and k, in the presence of a neutral, uniform plasma flowing at v 0 = v 0 u z (with Lorentz factor γ 0 ) on a periodic 2D Cartesian grid, Note that we only consider the case where the Galilean transformation is along the direction of the moving plasma, i.e., v gal = v gal u z , and that we use the following shorthand notations:  In fact, and as we have discussed in detail in Ref. [20], this rather complicated equation takes into account all numerical properties of the underlying algorithm, such as the finite temporal and spatial resolution, the discretization on the grid, and the particularities of the PIC cycle (e.g., the particle field gathering, charge and current deposition, current correction, etc.). In the limit of infinite resolution and for any value of v gal , Eq. (14) supports two physical modes, the relativistic electromagnetic mode ω 2 = c 2 k 2 + ω 2 p /γ 0 and the relativistic plasma mode ω = k z v 0 ± ω p γ −3/2 0 [20]. For finite resolutions, on the other hand, it recovers all of the unphysical couplings between distorted numerical modes that can lead to the numerical Cherenkov instability.
Especially the coefficients ξ are of interest, as they represent the numerical plasma response, taking into account the particle shape factorŜ (k), the spatial field smoothinĝ T (k), and a sum over all spatial aliases K m (K m = m x 2π x u x + m z 2π z u z and k m = k + K m ). The latter arise from the mismatch of sampling the continuous particle distribution on a discrete grid and play a key role in the occurrence of the NCI. It should again be noted that for infinite order [k] → k and Eq. (14) then equals to the original dispersion relation as derived in Ref. [20]. Furthermore, setting v gal = 0 recovers the dispersion relation of the standard PSATD solver that would be equivalent to the one formerly derived in Ref. [18].

B. Evaluation of the numerical Cherenkov instability growth rate
In the following section, we begin with an analytical evaluation of the dispersion relation, which closely follows [18] and reveals some of the sources of NCI growth. This analysis, however, is limited and Eq. (14) has to be solved numerically to obtain the NCI growth rate across all frequencies. Consequently, the results of this numerical evaluation are compared to actual 2D PIC simulations with WARPX. In the following, we consider three variants of the solver and use the same test case as in Refs. [3,18,20] with the exact parameters listed in Table I. For the sake of simplicity, we restrict ourselves to using a finite stencil in z only, thus [k] = [k z ] 2 + k 2 x . Nevertheless, the main findings, such as the elimination of the NCI, are equally valid when using finite-order stencils along multiple axes.
To begin with the analytical evaluation, it is useful to discard the plasma related terms of the dispersion relation and thereby recover the vacuum dispersion relation, i.e., the very first term of Eq. (14), Although the vacuum dispersion relation was already given by Eq. (12) and discussed in Sec. II, it now becomes evident that the above equation supports an infinite number of temporal aliases as seen by the particles, with n the temporal alias number. Even for the standard PSATD solver with infinite-order accuracy (PSATD ∞ ), this causes a distortion of the electromagnetic mode at large k if the time step t is too large. This is because particles sample the fields at discrete times and at maximum see a frequency of ±π/ t. If the actual frequency of the fields supported by the grid exceeds this value, then particles see a temporal alias of the electromagnetic mode. This can be seen in Fig. 4(a), where the electromagnetic mode (solid lines) of the PSATD ∞ solver is shown for a fixed transverse wave number k x = (π/2) x −1 . In our example, c t = 1.2 z > π/k max , where k max = π ( z −2 + x −2 ) 1/2 is the maximum wave number of the grid. Hence, particles see a temporal alias of the fields (gray lines) and ω decreases once the frequencies supported by the grid, c[k] exceed the maximum frequency supported by the time step, π c/ (1.2 z). Conversely, in the finite-order case (PSATD 8 ), the main electromagnetic mode is always distorted and eventually starts decreasing at high k z . These distorted electromagnetic modes can couple to spurious plasma or beam modes that emerge from the resonance condition in the plasma coefficients ξ , i.e., whenever the sine-term in the denominator of these expressions becomes zero. Here m z is the spatial alias number in the direction of the moving plasma, and n is the index of the temporal alias. Not surprisingly, for m z = 0 and n = 0 we get the NCR resonance condition that was already mentioned in Sec. II C, Eq. (13). The m z = 0 beam mode and the first two spatial aliases m z = [−1, +1] are shown in Fig. 4. Intersections of these beam modes-or their temporal aliases-with the electromagnetic mode can support resonant growth of the NCI.
As expected, the m z , n = (0, 0) resonance only occurs for finite orders (PSATD 8 ), while the standard PSATD solver (PSATD ∞ ) solely suffers from aliased resonances. By setting v gal = v 0 in Eq. (20) all resonances with spatial aliases m z = 0 vanish [ Fig. 4(c)], which is a direct consequence of the numerical grid following the moving plasma in the Galilean frame. In addition, the asymmetrical dispersion relation prevents the m z = 0 beam mode to intersect with the main electromagnetic mode in case of the finite-order Galilean-PSATD solver (Galilean-PSATD 8 ).
Despite that, Eq. (20) suggests that an NCI resonance with m z = 0 and a temporal alias n = 0 of the fields could still occur for any order of the Galilean-PSATD solver. For example, the m z , n = (0, 1) intersection in Fig. 4(a) does not FIG. 4. Distorted electromagnetic modes (black solid lines) and temporal aliases (gray solid lines) for forward-and backward-(±ω t) propagating waves and spurious beam modes (dashed lines) for different variants of the solver (see Table I vanish for v gal = v 0 . However, as we will learn from the numerical analysis in the following, the ν-dependent terms in the coefficients χ 5 and χ 5 seem to cancel this resonance for v gal = v 0 and also prevent the occurrence of an otherwise typical nonresonant growth of the NCI for m z = 0. These additional terms are a direct consequence of the analytical integration when deriving the PSATD equations and come from the fact that the currents are assumed to be constant over one time step in the comoving frame, which is equivalent to assume the currents to be comoving with the plasma in the original frame. Again, those final observations have been made by numerically solving Eq. (14), because it is not easily possible to extract more information from the analytical expressions. Figure 5 shows the theoretical NCI growth rates Im(ω) for the three solver cases and compares them to WARPX simulations. The possible curves of resonant growth are shown for the different beam modes and are labeled according to their spatial and temporal alias number. The k x locations of intersecting modes can be derived as a function of k z by inserting the vacuum dispersion relation Eq. (12) into the resonance condition Eq. (20) and solve it for k x [18]: The standard PSATD ∞ case suffers from nonresonant growth of the NCI over a wide range of frequencies and in addition shows resonant growth along the m z , n = (0, 1) and m z , n = (−1, 0) resonance lines. A very similar pattern of growth can be seen for the finite-order PSATD 8 case, although here three lines of resonant growth, m z , n = (0, 0), (0,1), and (−1, −1) are visible, with the m z , n = (0, 0) resonance being the dominant one. In contrast, numerically solving the dispersion relation for the Galilean-PSATD 8 case predicts absolutely no growth across all frequencies, proving the absence of the NCI for a uniformly streaming plasma. The estimated growth rates from the PIC simulations validate the above results and show excellent quantitative agreement for all three cases. Note that the visible residual background can be attributed to numerical noise in the simulation.
Although not shown here, another useful conclusion can be drawn by slightly modifying the dispersion relation before numerically solving it. Explicitly replacing [k z ] by k z in ω and ν removes the asymmetric properties of the dispersion relation and thereby reintroduces the m z , n = (0, 0) resonance for the otherwise stable case with v gal = v 0 . Indeed, when numerically solving this modified dispersion relation, NCI growth appears only along the m z , n = (0, 0) resonance line. This again highlights the self-consistency of the Galilean approach. Despite being a fundamental ingredient to cancel the resonance, the asymmetry of the vacuum dispersion relation was not artificially introduced to the scheme but rather it is a natural consequence of deriving the PIC equations in the comoving coordinates.
In summary, the stability analysis unfolded the different sources of NCI growth and we demonstrated how the Galilean coordinate transformation introduces unique properties that lead to their compensation. Consequently, we proved the elimination of the numerical Cherenkov instability for a finite, localized stencil in the case of a cold relativistic plasma moving uniformly in one direction. In the next section, the practical usefulness of the solver is validated by applying it to simulations of laser-plasma acceleration using the Lorentzboosted frame technique.

IV. APPLICATION TO PLASMA ACCELERATION
In a plasma wakefield accelerator [40], an intense driver, either a short laser or particle beam, excites a trailing plasmadensity wave, sustaining electric fields that are orders of magnitude higher than in radiofrequency-based accelerators. Particle-in-cell simulations of plasma accelerators are costly, because spatially resolving the small scale objects, such as the laser wavelength, limits the propagation distance within one  14) and simulated growth rates are estimated from the difference of the Fourier transformed fields between ω p,r t 90 and 120. Mode intersections that support resonant growth of the NCI are calculated from Eq. (21) and are overlaid as solid lines and labeled according to their spatial (m z ) and temporal (n) alias number. time step. Hence, many iterations are required to follow the interaction through the longer plasma.
A much more efficient ratio of scales can be achieved by transforming the entire simulation from the laboratory frame to a relativistic frame of reference moving in the direction of the driver [27], i.e., where γ = (1 − β 2 ) −1/2 and β = v boost /c is the speed of the relativistic frame in the z direction. In this Lorentz-boosted frame, the length scales of copropagating relativistic objects, such as the laser pulse or the plasma wave, are elongated by a factor γ (1 + β ). In contrast, objects that have been initially at rest, i.e., the plasma, are contracted by a factor γ and counterpropagate at a relativistic velocity v plasma = −βc in the opposite direction. As a consequence, the total simulation length decreases and the propagation distance per time step increases, reducing the overall number of iterations by (1 + β ) 2 γ 2 [28].
Although the Lorentz-boosted frame technique enables orders of magnitude speedups, its applicability is limited because the counterstreaming plasma eventually gives rise to an NCI. However, as we have originally shown in Ref. [19], we can use the infinite-order Galilean-PSATD solver to perform the simulation in a coordinate system that is comoving with the relativistic plasma, hence v gal = −βc, and thereby eliminate the NCI. In this section, we extend this result to the finiteorder Galilean-PSATD solver. This extension significantly increases the practical relevance of the Galilean scheme as it allows to perform large-scale, parallelized Lorentz-boosted frame simulations of plasma accelerators while still eliminating the NCI.
In the following, we show FBPIC simulations of a nonlinear plasma wakefield in the blow-out regime that is used to double the energy of a 1-GeV electron beam with 100 pC of charge over a distance of 30 mm. The plasma wave is resonantly driven by a laser pulse that is kept in focus by a matched plasma guiding channel. The setup has been simulated in the laboratory frame and in a Lorentz-boosted frame with γ boost = 10. Simulation parameters are summarized in Table II. Except for the infinite-order simulation, the computation is parallelized to eight graphics processing units by spatial domain decomposition in the axial direction.   Table  I). (a) Galilean-PSATD 8 ; (b) PSATD ∞ ; (c) PSATD 8 . Simulation parameters are listed in Table II. A destructive NCI develops unless setting v gal = v plasma . Figure 6 compares the stability of the Lorentz-boosted frame simulation for the three different variants of the solver in accordance with Sec. III. At the time shown, the laser pulse, propagating to the right, has already reached the end of the counterstreaming plasma. The trailing plasma wave and the accelerated electron beam are visible as modulations of the charge density ρ . While the onset of an NCI is visible in case of the standard PSATD ∞ solver, a strong NCI has already developed for the finite-order PSATD 8 solver. The distorted dispersion relation leads to a higher growth rate of the instability in this case. In contrast, the Galilean-PSATD 8 solver remains completely stable for the entire simulation distance, while the only difference here is the use of v gal = v plasma . Although the excitation of a strong nonlinear plasma wave causes deviations from the initially uniform plasma flow, no signs of an NCI become apparent. The high-density electron bunch that moves opposite to the Galilean frame is not degraded by NCR or other detrimental numerical effects. Despite the presence of a distorted vacuum dispersion relation in this direction, the accuracy of the solver for the relatively low order n order = 8 is already sufficient to mitigate NCR [17] in this use case.
The results of this finite-order Lorentz-boosted frame simulation are back-transformed to the laboratory frame and benchmarked against a high accuracy (n order = 128) reference simulation with γ boost = 1, i.e., a laboratory-frame simulation. Figure 7 shows a direct comparison of the longitudinal wakefield E z and the charge density ρ toward the end of the simulation after 25 mm of propagation. The guided laser expels almost all electrons from the axis and forms a trailing plasma cavity sustaining tens of GeV/m accelerating gradients. Within the blow-out, the high charge electron bunch drives its own wakefield and absorbs energy from the wake. The beam loading flattens the accelerating fields along the bunch and thereby minimizes the accumulation of linearly correlated energy spread. The on-axis accelerating field agrees well between both simulations with only small deviations near the region of highest density at the back of the first plasma oscillation where the electron trajectories collapse.
A comparison of the evolution of characteristic laser and electron beam properties is shown in Fig. 8. The laser is initially focused 250 μm behind the entrance of the plasma channel into a 1-mm-long sinusoidal density upramp. As it reaches relativistic intensities the nonlinear interaction with the plasma causes the pulse to self-modulate. As a consequence, the pulse length τ decreases and the transverse waist w oscillates along the propagation distance. The electron bunch is initially focused at the end of the upramp and situated at the back of the first plasma oscillation, 32.5 μm behind the centroid of the laser. While slowly dephasing in the direction of the driver, it gets accelerated to E kin 2.1 GeV and accumulates a relative energy spread of σ E 1.4%. A mismatch between the transverse beam phase space and the strong focusing forces leads to betatron decoherence causing a slight growth of the emittance that eventually saturates to ε n 1.3 mm mrad. Both the nonlinear laser plasma interaction and the electron beam dynamics are correctly reproduced in the Lorentz-boosted frame. At the end of the propagation distance, the shown quantities deviate by less than a percent between both simulations. The agreement is confirmed by a direct comparison of the longitudinal and transverse phase spaces at z prop = 30 mm.
In summary, the results prove the applicability of the finite-order Galilean-PSATD solver for typical simulations of plasma acceleration in the Lorentz-boosted frame. For the example case considered here, the required number of PIC iterations is reduced by two orders of magnitude, while the use of a local stencil of order n order = 8 allows for efficient parallelization by domain decomposition.

V. CONCLUSION AND DISCUSSION
In this article, we extended the formulation of the Galilean-PSATD solver [19,20] by incorporating spatial derivatives of finite order for the electromagnetic fields. This adds locality to the field solver [17,34] and allows to scale the algorithm to many distributed compute units via domain decomposition [15,35]. While the modification leads to errors in the vacuum dispersion, the solver retains its unique stability properties that lead to the elimination of the numerical Cherenkov instability for a relativistically streaming plasma. The natural discretization in Galilean coordinates introduces an asymmetry to the numerically distorted dispersion relation which cancels the main NCI resonance for a beam or plasma propagating at the comoving grid velocity.
The benefits of this asymmetry come at the cost of reducing the accuracy in the direction opposite to the Galilean frame. However, the tunability of the solver always allows to balance accuracy and locality to preserve the correct physics [17,36]. The proposed solver leads to no increase in computational cost compared to the standard PSATD solver, apart from a minor increase in communication overhead when parallelized. Variations of the finite-order Galilean-PSATD solver will be subject of future studies. For example, adding spatial staggering of the electromagnetic fields could further improve the scalability.
The solver was applied to practically relevant simulations of plasma acceleration in an optimal Lorentz-boosted frame of reference [27]. Although the excitation of a plasma wave leads to deviations from the uniform streaming velocity, we demonstrated that the solver remains stable and accurate. Furthermore, even when using low orders, a relativistic highdensity beam that counterpropagated with respect to the Galilean frame did not show any signs of an NCI and was not degraded by NCR.
Nevertheless, large deviations from the bulk plasma velocity or simulations of counterstreaming plasmas, such as studied in astrophysics, can again lead to an instability. Likewise, a dense relativistic beam that counterpropagates with respect to the Galilean frame can still degrade from NCR when using very low orders and grid resolutions. In these situations, however, the simultaneous combination of multiple comoving frames in more than one direction is a possible solution.
Apart from suppressing NCR, simulations including relativistic beams could further benefit from a Galilean frame. Recent work [41] showed that even for an ideal numerical dispersion and in the absence of Cherenkov like fields, simulations of relativistic beams can suffer from erroneous space charge fields. However, by modeling such systems in a frame comoving with the beam, it should be possible to recover the correct physical self-fields with the arbitrary-order Galilean-PSATD solver.
In this section, we derive the theoretical dispersion relation for a relativistic flowing plasma, with the finite-order Galilean scheme. The derivation closely follows that of Appendix B in Ref. [20], which derived the theoretical dispersion relation for the infinite-order case. Thus, in this section, we only emphasize the differences due to the finite-order scheme.
As a reminder, Appendix B in Ref. [20] considered the discretized Vlasov-Maxwell system for a uniform, flowing plasma, with small perturbations (in the fields E and B, and in the distribution function f ). The perturbations were then decomposed into Fourier eigenmodes, and a dispersion relation connecting ω and k for each eigenmode was obtained. Importantly, all numerical effects (e.g., particle shape, discretization on a grid, etc.) were taken into account in the derivation.