Particle in Cell (PIC) simulations are a widely used tool for the investigation of both laser- and beam-driven plasma acceleration. It is a known issue that the beam quality can be artificially degraded by numerical Cherenkov radiation (NCR) resulting primarily from an incorrectly modeled dispersion relation. Pseudo-spectral solvers featuring infinite order stencils can strongly reduce NCR—or even suppress it—and are therefore well suited to correctly model the beam properties. For efficient parallelization of the PIC algorithm, however, localized solvers are inevitable. Arbitrary order pseudo-spectral methods provide this needed locality. Yet, these methods can again be prone to NCR. Here, we show that acceptably low solver orders are sufficient to correctly model the physics of interest, while allowing for parallel computation by domain decomposition.

## INTRODUCTION

Plasma accelerators provide high accelerating gradients and promise very compact next-generation particle accelerators as drivers for brilliant light sources^{1,2} and high energy physics.^{3} Using laser or beam drivers, the acceleration of electron beams by multiple GeV over few-cm distances has been demonstrated in experiments.^{4–7} The improvement of beam quality, especially in terms of energy spread, remains crucial for applications. Due to the non-linear nature of the beam-plasma interaction, Particle-in-Cell (PIC) simulations are a vital tool to model the beam injection and transport in the plasma. PIC algorithms self-consistently solve the interaction of spatially discretized electromagnetic fields with charged particles defined in a continuous phase space. However, these codes can suffer from numerical Cherenkov radiation (NCR)^{8} that is caused primarily by incorrect modeling of the electromagnetic dispersion relation. Mitigation of this effect is crucial as it gives rise to an unphysical degradation of beam quality, which leads to wrong predictions when using those results as input for studies on applications such as compact Free-Electron Lasers. For example, NCR has been found to artificially decrease the beam quality in terms of transverse beam emittance^{9} for finite-difference time domain (FDTD) solvers, such as the commonly used Yee scheme. Efforts to circumvent this issue include the modification of the FDTD dispersion relation^{9,10} and digital filtering.^{11,12} In contrast to FDTD solvers, pseudo-spectral solvers, which solve Maxwell's equations in spectral space, strongly reduce NCR—and sometimes even suppress it.

In practice, PIC simulations can easily demand thousands of hours of computation time. Therefore, implementations that are efficient and scalable for parallel production are unavoidable. The parallelization of PIC codes is typically done by decomposing the simulation volume into domains, which are evaluated by individual processes. For FDTD algorithms, efficient scaling to several thousand parallel processes is possible. Pseudo-spectral solvers, however, act globally on the grid due to the needed Fourier transform, which limits efficient parallelization. An approach to overcome this limitation are *arbitrary order* pseudo-spectral solvers,^{13–15} which only locally affect the fields when solving Maxwell's equations. On the downside, this—again—introduces spurious numerical dispersion and can make a simulation prone to NCR.

In the following, we investigate numerical Cherenkov radiation and its implications on the physics of a simulation for arbitrary order pseudo-spectral solvers using the spectral, quasi-cylindrical PIC algorithm implemented in the code FBPIC.^{16,17}

This paper is structured as follows. We will first summarize the concept of pseudo-spectral solvers, with an emphasis on the arbitrary order,^{13} *pseudo-spectral analytic time domain* (PSATD),^{18} followed by a discussion of the main characteristics of NCR. In the last part, we numerically analyze the mitigation of NCR in the arbitrary order PSATD scheme. As a figure of merit for the physical accuracy of the solver, we use the effect of NCR on the beam emittance and the slice energy spread. With parameter scans, we show the dependence of NCR on physical and numerical parameters of the simulation.

## FINITE ORDER PSEUDO-SPECTRAL SOLVERS

FDTD solvers are a widely used method in PIC codes. These algorithms use finite difference operators to approximate derivatives in time and space when numerically solving Maxwell's equations, which introduces a spurious numerical dispersion. This approximation can be done with arbitrary precision by increasing the used stencil order, i.e., taking contributions from more surrounding grid nodes into account. Yet, this has the drawback of also increasing the computational demand of the calculation. Further, common FDTD algorithms represent the electromagnetic fields on a staggered grid, where the electric and magnetic fields are not defined at the same points in time and space. This is known to cause errors when modeling the direct interaction between an electron beam and a co-propagating laser.^{19}

The use of finite differences in space can be avoided by transforming Maxwell's equations into spectral space. In that case, the differential operator $\u2202z$ is given by a multiplication with *ik _{z}*, which can be numerically evaluated without approximations. The resulting solver class is referred to as the

*pseudo-spectral time domain*(PSTD),

^{20}and it describes the limit of an FDTD solver with an infinite stencil order.

^{14}Yet, this solver still applies finite difference approximations in time, and thereby, like FDTD solvers, its timestep is limited by a Courant condition. It is subject to spurious numerical dispersion, but unlike FDTD schemes, the erroneous dispersion caused by the PSTD scheme results in modes traveling faster than the speed of light.

^{21}

Under the assumption of constant currents ** J** over one timestep $\Delta t$, this constraint can be overcome by integrating Maxwell's equations analytically in time. With a precise representation of derivatives in time and space, the resulting

*pseudo-spectral analytic time domain*(PSATD) scheme is not subject to a Courant condition and is free of spurious numerical dispersion. Furthermore, it is currently the only solver that can be combined with the Galilean scheme

^{22,23}that is intrinsically stable when modeling relativistically drifting plasmas with uniform velocity.

The PSATD field propagation equations in Cartesian coordinates are

where ** k** is the wave vector of length $k\u2261kx2+ky2+kz2$ with $k\u0302=k/k$, $F\u0303$ is the Fourier transform of a field

*F*, with

*E*and

*B*being the electric and magnetic field, and

*ρ*and

*J*the charge and current density. Further, $C=cos(kc\Delta t)$ and $S=sin(kc\Delta t$). With pseudo-spectral methods, the fields are naturally centered in space and time, eliminating a source of interpolation errors due to the usual staggering of the electromagnetic fields in the popular Yee FDTD solver. The equations of the PSTD scheme can be retrieved from those of the PSATD scheme, by writing the PSATD equations in an equivalent staggered form

^{21}and performing a first-order Taylor expansion in $\Delta t$.

In the following, we will use the PSATD scheme in a quasi-cylindrical geometry as derived in Ref. 16. In this geometry, the fields are decomposed into azimuthal modes *m*, which can then be represented on separate 2D grids.^{24} In practice, as the contributions from most of the modes become zero for near-symmetrical systems, only a few of these have to be taken into account. Consequently, this geometry allows the modeling of near-symmetrical problems at a drastically reduced computational effort compared to full 3D simulations. In the quasi-cylindrical PSATD scheme, the real-space components *F _{r}*, $F\theta $, and

*F*of a field

_{z}**are connected to their spectral components $F\u0303+,m,\u2009F\u0303\u2212,m$, and $F\u0303z$ by the combination of a Hankel transform (along**

*F**r*) and a Fourier transform (along

*z*). The equations of the quasi-cylindrical PSATD can be derived from the Cartesian PSATD scheme by replacing the differential operators in Eqs. (1a) and (1b) by

^{23}

where $S\u0303$ is any scalar field and $V\u0303$ any vector field. This representation of the PSATD scheme preserves the advantages of the Cartesian PSATD scheme, i.e., an ideal dispersion relation, a centered representation of the fields and no Courant condition for the timestep and combines them with the reduced computational demand associated with a quasi-cylindrical geometry.

The propagation of the fields in spectral space in Eqs. (1a) and (1b) corresponds to a global operation in spatial space, where information is transferred over the entire simulation grid. However, for domain decomposition, a locality of the solver operations is necessary. In order to achieve an arbitrary order local stencil, the solver equations are modified.

The formulations of the modified pseudo-spectral scheme can be derived from the expression of a finite difference with the arbitrary order. For a centered scheme, the discrete derivative with stencil order $2\u2009n$ of a function *F* at a discrete point *p* on the spatial grid is^{25}

To get a spectral method with a finite stencil derivative, consider the spectral representation of this operator

Here, $[kz]$ is the modified wavenumber given as

with $\Delta z$ being the grid spacing and $F\u0303$ being the spectral form of *F*. This means that the equivalent to a finite stencil in real space can be achieved by using a modified PSTD algorithm, whereby *k _{z}*,

*k*, and

_{x}*k*are replaced by $[kz],\u2009[kx],\u2009and\u2009[ky]$ in the discretized Maxwell equations in spectral space. Because this modified PSTD algorithm is mathematically equivalent to an FDTD scheme (with high-order spatial derivatives), it has the same locality.

_{y}Similar to the PSTD algorithm, in order to obtain a PSATD scheme with improved locality, we follow identical steps. More precisely, the discretized Maxwell equations Eqs. (1a) and (1b) are modified by replacing *k _{z}*,

*k*, and

_{x}*k*by $[kz],\u2009[kx],\u2009and\u2009[ky]$, including in the expression of $k\u2261kx2+ky2+kz2$, which in turn impacts the expressions of the coefficients

_{y}*C*and

*S*. Consequently, unlike the modified PSTD algorithm, the modified PSATD scheme is not strictly local. In order to evaluate its locality, we get the corresponding stencil coefficients on the spatial grid, by applying an inverse Fourier transform to the operators used to advance the fields in Eqs. (1a) and (1b), i.e., $[k]\u0302S$ and

*C*.

Next, we will use this scheme in quasi-cylindrical geometry with the PIC framework FBPIC to verify the locality of the solver. In FBPIC, only *k _{z}* is modified because the Fourier transform is only used along the

*z*-axis. Therefore, the corresponding solver stencil remains infinite in the radial direction, and the locality is only gained in the longitudinal direction. To examine this, a

*δ*-peak in E and B is initialized on the axis in the center of the simulation box with a grid resolution of $\Delta r=\Delta z=1\u2009\xb5m$. It is then propagated for one timestep $\Delta t=\Delta z/c$. Fig. 1 (top) shows the expansion of the pulse for a stencil order of 32 compared to the standard PSATD. When using the unmodified PSATD ($PSATD\u221e$) scheme, the simulation shows field contributions along the entire axis, whereas the fields modeled with the finite order PSATD quickly decrease in amplitude and reach the machine precision level. The expansion of the unit pulse strictly follows the stencil coefficients $\alpha *$ derived from the modified PSATD operations given by

where $F\u22121$ denotes an inverse Fourier transform along the *z*-axis. Here, these coefficients are calculated numerically, and instead of an analytical Fourier transform, an FFT was used. The coefficients are evaluated for $k\u22a5=0.5\u2009\xb5m\u22121$. The stencil coefficients as well as the unit pulse reach machine precision related noise after decreasing by 16 orders of magnitude over a range of $\Delta Nz$ cells. In order to show how the stencil extent $\Delta Nz$ scales with the arbitrary order PSATD, the same simulation is done with varying stencil orders. Fig. 1 (bottom) shows that the stencil extent behaves non-linearly to the applied stencil order. While for high orders, the extent is similar to the stencil order, the reach of the solver tends to be bigger than the order in low order cases.

In a domain decomposition scheme, the fields would only need to be exchanged between neighboring processes in guard regions with the size of $\Delta Nz$ cells as this is the maximum distance information can travel within one timestep. These guard cells hold copies of the fields in the neighboring domains. At each timestep, the fields are evaluated independently in each domain, and the information is only exchanged in these cells. When significantly less than $\Delta Nz$ guard cells are used, truncation of the stencil can lead to spurious fields.^{14}

Note that for parallelization also, a localized current correction or a charge conserving current deposition is needed. For that, an exact current deposition in k-space for Cartesian coordinates^{21} and a local FFT based current correction on staggered grids^{15} have been shown.

To summarize, in high order FDTD schemes, the computational demand scales with the stencil order as an extension of the stencil to more neighboring cells calls for more individual numerical operations. In contrast, the stencil in spatial space for the arbitrary order PSATD is merely determined by the modification of the *k* values in spectral space. Apart from this, the scheme is identical to the infinite order PSATD. Therefore, the computational demand of the solver is independent of the stencil order and consequently also of the precision of the field solver. However, the improved locality of these schemes comes along with approximations on the integration of the electromagnetic fields, which also result in spurious numerical dispersion. The implications of this will be discussed in the following.

## NUMERICAL CHERENKOV RADIATION

In general, the phase velocity $v\varphi $ of an electromagnetic mode with wave vector ** k** and angular frequency

*ω*is $v\varphi 2=\omega 2/k2$. The allowed electromagnetic modes in a vacuum are described by the dispersion relation

When modifying ** k** in the longitudinal direction, the dispersion relation is

Here, a spurious numerical dispersion is introduced, where the phase velocity in a vacuum decreases for modes with increasing *k _{z}*. Fig. 2 shows the distorted dispersion relation caused by the arbitrary order PSATD solver compared to the accurate case (Eq. (7)) given by the unmodified, infinite order PSATD method. Due to the distortion of the dispersion relation, relativistic particles can travel with the same velocity as some electromagnetic modes. This can lead to coherent emission of spurious radiation referred to as numerical Cherenkov radiation. Consequently, the affected modes obey

where $vbeam$ is the velocity of the charged particles in the simulation. Using Eq. (8), this can be rewritten as

This condition describes the modes, where the spurious electromagnetic dispersion relation intersects with the beam mode corresponding to the velocity of a charged particle. This effect is illustrated in Fig. 3. The dashed lines indicate the modes fulfilling Eq. (10) and thus correspond to the intersection points of the dispersion relation and the beam mode as depicted in Fig. 2. These intersection points denote the modes where NCR is excited. They are overlaid with the spatial Fourier transform of the longitudinal electric field from FBPIC simulations, with parameters as discussed in the section Accurate Modeling of Plasma Accelerators in Finite Order PSATD. The panels of Fig. 3 show the results from simulations done with different stencil orders. They all show excited modes that can be attributed to the physics in the simulation of a beam driven plasma accelerator, e.g., the plasma wave or betatron radiation. However, compared to the unmodified $PSATD\u221e$ scheme simulation, additional excited modes are visible around the modes prone to NCR. The excitations are especially strong for low solver orders. They decrease for higher stencil orders as the resonant modes of Eq. (10) move away from the modes that can be attributed to the physics of the simulation. Eventually, for the $PSATD\u221e$, the phase velocity $v\varphi $ is always greater than $vbeam$ (i.e., no dashed line), and no NCR is present.

The analysis of the Fourier transformed electric field consequently can be used as an indication for the presence and the magnitude of NCR. Next, we will show the influence of these unphysical modes on the properties of an electron beam.

## ACCURATE MODELING OF PLASMA ACCELERATORS IN THE FINITE ORDER PSATD

To show that already low stencil orders are sufficient to mitigate errors due to NCR in the context of plasma-based acceleration, PIC simulations with the spectral, quasi-cylindrical code FBPIC for beam-driven plasma accelerator parameters are presented. Typical setups^{6,26} feature high bunch charges and consequently high currents, as they aim for applications in future brilliant light sources and high energy physics.

The propagation of cold, i.e., zero emittance, monoenergetic driver and witness electron beams through a plasma of density $1\xd71017\u2009cm\u22123$ is simulated. A non-linear wakefield is driven by a 1 GeV Gaussian beam with 180 pC charge, $8.4\u2009\xb5m$ rms length, and $5\u2009\xb5m$ rms width, compare Fig. 4(a). The 25 MeV Gaussian witness bunch of 50 pC charge, $2\u2009\xb5m$ rms length, and $1\u2009\xb5m$ rms width is positioned in the accelerating region at the back of the bubble.

The simulation box has a resolution of 5 cells/$\xb5m$ longitudinally and 1.5 cells/$\xb5m$ transversally, with 16 particles per cell distributed as $2\xd72\xd74$ in longitudinal, radial, and azimuthal directions, respectively. A total of 900 × 300 cells in two azimuthal modes is used. While the pseudo-spectral solver of FBPIC is technically not limited by a CFL condition, a timestep of $\Delta t=\Delta z/c$ was used. For the same physical parameters, the stencil order of the arbitrary order PSATD is varied. The NCR-free simulation with an infinite order stencil is used as a reference in the following.

Both Figs. 4 and 5 show simulation results after a propagation through 35 mm of plasma. In Figs. 4(b) and 4(c), the accelerating field *E _{z}* and the charge density

*ρ*in the vicinity of the witness bunch can be seen. In the case of a low stencil order, strong distortions are visible in

*E*. The distortions are the spatial space equivalent to resonances visible for the $PSATD4$ in Fig. 3. These spurious fields affect the witness beam shape, which corresponds to an increase of beam emittance, as it has been reported for 3D FDTD codes.

_{z}^{9}

We also find that the longitudinal phase space is affected, as shown in Fig. 5, which compares the phase spaces modeled with several stencil orders. The projected energy spread increases during propagation through the plasma due to the slope of *E _{z}* in combination with beam loading, so that an s-shape like structure is imprinted on the phase space also with $PSATD\u221e$. However, for decreasing stencil order, NCR causes an increasingly violent high frequency oscillation on the longitudinal phase space, leading to an unphysical growth of both slice and projected energy spread. It can be observed that the perturbation is especially pronounced behind the beam center, where the current density is the highest, while the head remains mostly unaffected.

Fig. 6 (top) shows the evolution of the slice energy spread during the propagation through the plasma for different solver orders. For finite order stencils, the energy spread grows non-linearly with *z*. While $PSATD4$ and $PSATD8$ still deviate significantly from the Cherenkov free infinite order case, the impact of NCR on the beam in the PSATD_{16} simulation is barely visible even after a long propagation distance. The dashed lines correspond to simulations done with the FDTD Yee scheme in 3D Cartesian and quasi-cylindrical geometry. These are performed with the PIC code Warp^{27} using the same parameters and resolution as FBPIC. With the Yee scheme, the growth rate of spurious slice energy spread is even greater than for the low order PSATD case. Note that the Warp framework also features a $PSATD\u221e$ (and arbitrary order) solver, which, if used for this problem, would also produce the correct physical results as the FBPIC $PSATD\u221e$ algorithm.

Fig. 6 (bottom) shows the error of the slice energy spread, the transverse beam emittance, and the longitudinal electric field after a propagation length of 35 mm depending on the stencil order. Here, the errors are defined as

respectively. The subscripts *f* and $\u221e$ indicate the quantities from finite and infinite order simulations. For the calculation of $Err\sigma \gamma ,\u2009Nslice$ slices of length 0.5 µm within $\xb13\sigma z$ of the witness beam are considered. For $ErrEz$, the sum is calculated over all the cells of the spatial Fourier transform of *E _{z}* (compare Fig. 3). This quantity can be regarded as a measure for the strength of the spurious fields caused by NCR.

As can be expected from Fig. 3, Fig. 6 (bottom) confirms that the unphysical contributions of NCR to the field, $ErrEz$, depend strongly on the stencil order and decrease by 8 orders of magnitude as the stencil order increases from 4 to 256. With the strength of NCR decreasing, the beam quality in terms of energy spread and emittance also converges toward the correct values given by the reference case with a PSATD$\u2009\u221e$. The beam's energy spread is even more sensitive to potential NCR than the beam emittance. $Err\u03f5$ and $Err\sigma \gamma $ are below 0.05 for the stencil orders of 12 and 32, respectively. In the following, this will also be used as a criterion for sufficient accuracy.

The stencil order necessary to meet this criterion depends on the physical and numerical parameters of the simulation. Fig. 7 is obtained from scans over the beam charge, the propagation length as well as the longitudinal grid resolution and shows the required stencil order to achieve $Err\sigma \gamma <0.05$. Fig. 7 (top) shows the required stencil order depending on the propagation length in plasma for the same parameters used in Figs. 3–6. For a reduced propagation length, accurate modeling is possible with lower stencil orders. This is a result of the difference in the growth rate of the spurious slice energy spread seen in Fig. 6 (top).

NCR also scales with the bunch charge density and the longitudinal grid resolution. The connection of the required stencil order with these parameters is shown in Fig. 7 (bottom). As the computational demand increases quadratically with the grid resolution, this scan and for consistency also the charge scan are evaluated after a propagation length of 20 mm. Clearly, high charges call for increased stencil orders (see Fig. 7(b)). However, note that in the present case, due to the small beam dimensions, 50 pC already corresponds to a high peak current of 3 kA and a peak electron density of $1\xd71019\u2009cm\u22123$.

The required stencil order can be reduced by increasing the longitudinal grid resolution (see Fig. 7(c)). However, for the propagation length and charge considered here, it is necessary to apply a stencil order of at least 8 even when using a resolution of 20 cells/*μ*m.

## CONCLUSION

In conclusion, NCR stemming from erroneous modeling of the dispersion relation leads to beam phase space degradation, which manifests itself in a spurious growth of slice energy spread and emittance. This effect is present in both the FDTD Yee scheme and the finite order PSATD solvers. Errors from NCR rapidly decrease for higher stencil orders, eventually showing the same results as the NCR-free infinite order PSATD solver.

In the case discussed here, which is inspired by typical beam driven plasma acceleration parameters, a stencil order of 32 effectively suppresses spurious beam quality degradation. We have chosen a long propagation distance of 35 mm and a high witness beam charge of 50 pC corresponding to an electron density of $1\xd71019$ cm^{−3}. This is a conservative example in the sense that our results are applicable to many setups in the field of plasma acceleration that use more moderate parameters. This is confirmed by parameter scans that show that the constraints on the required stencil order are relaxed for a lower beam charge and a shorter propagation length.

In practice, this means that with around 60 guard cells between neighboring domains, parallelization is possible while retaining the physics of the considered problem.

Analogous to the frequently used Yee scheme, growth of NCR in the finite order stencil PSATD scheme can also be reduced by increasing the resolution of the simulation grid. However, this is an inefficient solution especially in schemes where the timestep is limited by a CFL-condition and hence directly linked to the grid resolution. Yet, this means that in simulations with a high spatial resolution, the effects of NCR are less severe. This, for example, would be the case for laser plasma acceleration, where a high resolution is required in any case to resolve the laser wavelength. Therefore, in these cases, even lower order stencils than suggested here are sufficient to model the physics without artifacts from NCR.

The arbitrary order PSATD scheme preserves the benefits of pseudo-spectral solvers, e.g., the analytic integration in time or a centered definition of the electric and magnetic fields. It further allows the independent adjustment of the precision of the electromagnetic dispersion relation and the resolution of the simulation. This way the needed spatial resolution is governed by the physical problem and not by the mitigation of NCR.

Input scripts to reproduce the presented FBPIC and Warp simulations can be found in Ref. 28.

## ACKNOWLEDGMENTS

We gratefully acknowledge the computing time provided on the supercomputer JURECA under Project No. HHH20. 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 the Berkeley Lab.