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.

Plasma accelerators provide high accelerating gradients and promise very compact next-generation particle accelerators as drivers for brilliant light sources1,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 emittance9 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 relation9,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,13pseudo-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.

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 z is given by a multiplication with ikz, 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 Δ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 scheme22,23 that is intrinsically stable when modeling relativistically drifting plasmas with uniform velocity.

The PSATD field propagation equations in Cartesian coordinates are

Ẽn+1=CẼn+iSk̂×B̃nSckJ̃n+1/2+ik̂k[(SckΔt1)ρ̃n+1+(CSckΔt)ρ̃n],
(1a)
B̃n+1=CB̃niSk̂×Ẽn+i1Cckk̂×J̃n+1/2,
(1b)

where k is the wave vector of length kkx2+ky2+kz2 with k̂=k/k, F̃ 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Δt) and S=sin(kcΔ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 form21 and performing a first-order Taylor expansion in Δ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 Fr, Fθ, and Fz of a field F are connected to their spectral components F̃+,m,F̃,m, and F̃z by the combination of a Hankel transform (along 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) by23 

F̃=ikS̃(F̃+,mF̃,mF̃z,m)=(kS̃m/2kS̃m/2ikzS̃m),
(2a)
F̃=ik×Ṽ(F̃+,mF̃,mF̃z,m)=(kzṼ+,mikṼz,m/2kzṼ,mikṼz,m/2ikṼ+,m+ikṼ,m,),
(2b)

where S̃ is any scalar field and Ṽ 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 2n of a function F at a discrete point p on the spatial grid is25 

(DzF)p=j=1nαn,jFp+jFpj2jΔzαn,j=(1)j+12(n!)2(nj)!(n+j)!.
(3)

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

(DzF)p=kzi[kz]F̃kzeikzpΔz.
(4)

Here, [kz] is the modified wavenumber given as

[kz]=j=1nαn,jsin(kzjΔz)jΔz,
(5)

with Δz being the grid spacing and F̃ 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 kz, kx, and ky are replaced by [kz],[kx],and[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.

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 kz, kx, and ky by [kz],[kx],and[ky], including in the expression of kkx2+ky2+kz2, which in turn impacts the expressions of the coefficients 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]̂S 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 kz 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 Δr=Δz=1µm. It is then propagated for one timestep Δt=Δ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) 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 α* derived from the modified PSATD operations given by

α*=|F1(C)2|+|F1(Sk)2|+|F1(S[kz])2|,
(6)

where F1 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=0.5µm1. 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 ΔNz cells. In order to show how the stencil extent Δ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.

FIG. 1.

Expansion of an electromagnetic unit pulse after one PIC cycle with modified kz (top). The field pushed with the unmodified PSATD method (green line) shows non-negligible field contributions along the entire simulation grid, while the fields modeled with the order 32 PSATD (black line) decrease by 16 orders of magnitude over ΔNz cells and reach the machine precision level. The number of cells over which the fields decrease (bottom) depends on the stencil order. The dashed lines show the stencil coefficients derived from the modified PSATD equations.

FIG. 1.

Expansion of an electromagnetic unit pulse after one PIC cycle with modified kz (top). The field pushed with the unmodified PSATD method (green line) shows non-negligible field contributions along the entire simulation grid, while the fields modeled with the order 32 PSATD (black line) decrease by 16 orders of magnitude over ΔNz cells and reach the machine precision level. The number of cells over which the fields decrease (bottom) depends on the stencil order. The dashed lines show the stencil coefficients derived from the modified PSATD equations.

Close modal

In a domain decomposition scheme, the fields would only need to be exchanged between neighboring processes in guard regions with the size of Δ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 Δ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 coordinates21 and a local FFT based current correction on staggered grids15 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.

In general, the phase velocity vϕ of an electromagnetic mode with wave vector k and angular frequency ω is vϕ2=ω2/k2. The allowed electromagnetic modes in a vacuum are described by the dispersion relation

ω2=c2k2.
(7)

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

ω2=c2(kr2+[kz]2).
(8)

Here, a spurious numerical dispersion is introduced, where the phase velocity in a vacuum decreases for modes with increasing kz. 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

vϕ,z=vbeam,
(9)

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

vbeam=ckr2+[kz]2kz.
(10)

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 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, the phase velocity vϕ is always greater than vbeam (i.e., no dashed line), and no NCR is present.

FIG. 2.

Dispersion relation for the arbitrary order PSATD solver (solid lines) on the kz-Axis (kr = 0). Finite stencils introduce a spurious dispersion relation, which leads to a reduced phase velocity vϕ,z=ω/kz of modes with large kz and allows unphysical intersections (box) with the mode corresponding to the velocity of a relativistic particle (dashed line), here γbeam=3. Particle beams can then travel with the same velocity as electromagnetic modes, which leads to coherent emission of spurious radiation, referred to as NCR.

FIG. 2.

Dispersion relation for the arbitrary order PSATD solver (solid lines) on the kz-Axis (kr = 0). Finite stencils introduce a spurious dispersion relation, which leads to a reduced phase velocity vϕ,z=ω/kz of modes with large kz and allows unphysical intersections (box) with the mode corresponding to the velocity of a relativistic particle (dashed line), here γbeam=3. Particle beams can then travel with the same velocity as electromagnetic modes, which leads to coherent emission of spurious radiation, referred to as NCR.

Close modal
FIG. 3.

Spatial Fourier transforms of the accelerating field for simulations with different PSATD stencil orders. The simulations with the finite order show resonances around the modes with a phase velocity matching the velocity of the electron beam (dashed lines) given by Eq. (10). The magnitude of NCR decreases for increasing stencil order. The NCR-free infinite order PSATD simulation is shown for reference.

FIG. 3.

Spatial Fourier transforms of the accelerating field for simulations with different PSATD stencil orders. The simulations with the finite order show resonances around the modes with a phase velocity matching the velocity of the electron beam (dashed lines) given by Eq. (10). The magnitude of NCR decreases for increasing stencil order. The NCR-free infinite order PSATD simulation is shown for reference.

Close modal

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.

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 setups6,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×1017cm3 is simulated. A non-linear wakefield is driven by a 1 GeV Gaussian beam with 180 pC charge, 8.4µm rms length, and 5µm rms width, compare Fig. 4(a). The 25 MeV Gaussian witness bunch of 50 pC charge, 2µm rms length, and 1µm rms width is positioned in the accelerating region at the back of the bubble.

FIG. 4.

(a) Simulation setup featuring a driver and witness beam after propagating 35mm. Colors correspond to the accelerating field, the charge density is given in gray (arb. units), and ζ=zct are the co-moving coordinates with ζ = 0 at the center of the driver beam. The short PSATD4 stencil leads to distortions of the longitudinal field around the witness beam (b) that are not present for an infinite order stencil (c).

FIG. 4.

(a) Simulation setup featuring a driver and witness beam after propagating 35mm. Colors correspond to the accelerating field, the charge density is given in gray (arb. units), and ζ=zct are the co-moving coordinates with ζ = 0 at the center of the driver beam. The short PSATD4 stencil leads to distortions of the longitudinal field around the witness beam (b) that are not present for an infinite order stencil (c).

Close modal

The simulation box has a resolution of 5 cells/µm longitudinally and 1.5 cells/µm transversally, with 16 particles per cell distributed as 2×2×4 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 Δt=Δ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 Ez 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 Ez. 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.9 

FIG. 5.

Longitudinal phase space of the witness beam after propagation through 35mm of plasma. ξ is the internal bunch coordinate with ξ = 0 at the center of the witness beam. For decreasing stencil order, it is increasingly spoiled by a high frequency oscillation. This causes a growth of projected and slice energy spread (slice length of 0.5μm here).

FIG. 5.

Longitudinal phase space of the witness beam after propagation through 35mm of plasma. ξ is the internal bunch coordinate with ξ = 0 at the center of the witness beam. For decreasing stencil order, it is increasingly spoiled by a high frequency oscillation. This causes a growth of projected and slice energy spread (slice length of 0.5μm here).

Close modal

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 Ez in combination with beam loading, so that an s-shape like structure is imprinted on the phase space also with PSATD. 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 PSATD16 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 Warp27 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 (and arbitrary order) solver, which, if used for this problem, would also produce the correct physical results as the FBPIC PSATD algorithm.

FIG. 6.

Top: Growth of the witness beam's slice energy spread during the propagation through 35mm of plasma. The slice energy spread is averaged over all the slices within ±3 standard deviations around the bunch center. The solid and dashed lines show data obtained with the PSATD and Yee scheme, respectively. The PSATD16 and PSATD lines almost overlay each other. Bottom: Errors in emittance (blue), slice energy spread (black), and longitudinal field (red crosses) after 35mm. For a stencil order increasing from 4 to 256, the error introduced by NCR in the electric field decreases by 8 orders of magnitude. Consequently, the error of the final emittance and slice energy spread quickly approaches zero.

FIG. 6.

Top: Growth of the witness beam's slice energy spread during the propagation through 35mm of plasma. The slice energy spread is averaged over all the slices within ±3 standard deviations around the bunch center. The solid and dashed lines show data obtained with the PSATD and Yee scheme, respectively. The PSATD16 and PSATD lines almost overlay each other. Bottom: Errors in emittance (blue), slice energy spread (black), and longitudinal field (red crosses) after 35mm. For a stencil order increasing from 4 to 256, the error introduced by NCR in the electric field decreases by 8 orders of magnitude. Consequently, the error of the final emittance and slice energy spread quickly approaches zero.

Close modal

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

Errσγ=1NsliceiNslice(σγ,f,iσγ,,iσγ,,i)2,
(11)
Errϵ=ϵfϵϵ,
(12)
ErrEz=1NzNrNz,Nr(FFT2(Ez)fFFT2(Ez)FFT2(Ez))2,
(13)

respectively. The subscripts f and indicate the quantities from finite and infinite order simulations. For the calculation of Errσγ,Nslice slices of length 0.5 µm within ±3σz of the witness beam are considered. For ErrEz, the sum is calculated over all the cells of the spatial Fourier transform of Ez (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. The beam's energy spread is even more sensitive to potential NCR than the beam emittance. Errϵ and Errσγ 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σγ<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).

FIG. 7.

Starting from the parameters used in Figs. 3–6, either the charge, the propagation length, or the longitudinal grid resolution is scanned, and the required stencil order to fulfill Errσγ<0.05 is shown. Top: Required stencil order depending on the propagation length for a beam charge of 50 pC and a resolution of 5 cells/μm. Bottom: Required stencil order depending on the beam charge (b) with a resolution of 5 cells/μm and depending on the longitudinal grid resolution (c) with a beam charge of 50 pC, both evaluated after 20 mm of propagation in plasma.

FIG. 7.

Starting from the parameters used in Figs. 3–6, either the charge, the propagation length, or the longitudinal grid resolution is scanned, and the required stencil order to fulfill Errσγ<0.05 is shown. Top: Required stencil order depending on the propagation length for a beam charge of 50 pC and a resolution of 5 cells/μm. Bottom: Required stencil order depending on the beam charge (b) with a resolution of 5 cells/μm and depending on the longitudinal grid resolution (c) with a beam charge of 50 pC, both evaluated after 20 mm of propagation in plasma.

Close modal

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×1019cm3.

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.

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×1019 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.

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.

1.
A. R.
Maier
,
A.
Meseck
,
S.
Reiche
,
C. B.
Schroeder
,
T.
Seggebrock
, and
F.
Grüner
,
Phys. Rev. X
2
,
031019
(
2012
).
2.
Z.
Huang
,
Y.
Ding
, and
C. B.
Schroeder
,
Phys. Rev. Lett.
109
,
204801
(
2012
).
3.
C. B.
Schroeder
,
E.
Esarey
,
C. G. R.
Geddes
,
C.
Benedetti
, and
W. P.
Leemans
,
Phys. Rev. Spec. Top. Accel. Beams
13
,
101301
(
2010
).
4.
W. P.
Leemans
,
B.
Nagler
,
A. J.
Gonsalves
,
C.
Toth
,
K.
Nakamura
,
C. G. R.
Geddes
,
E.
Esarey
,
C. B.
Schroeder
, and
S. M.
Hooker
,
Nat. Phys.
2
,
696
(
2006
).
5.
W. P.
Leemans
,
A. J.
Gonsalves
,
H.-S.
Mao
,
K.
Nakamura
,
C.
Benedetti
,
C. B.
Schroeder
,
C.
Tóth
,
J.
Daniels
,
D. E.
Mittelberger
,
S. S.
Bulanov
,
J.-L.
Vay
,
C. G. R.
Geddes
, and
E.
Esarey
,
Phys. Rev. Lett.
113
,
245002
(
2014
).
6.
M.
Litos
,
E.
Adli
,
W.
An
,
C. I.
Clarke
,
C. E.
Clayton
,
S.
Corde
,
J. P.
Delahaye
,
R. J.
England
,
A. S.
Fisher
,
J.
Frederico
,
S.
Gessner
,
S. Z.
Green
,
M. J.
Hogan
,
C.
Joshi
,
W.
Lu
,
K. A.
Marsh
,
W. B.
Mori
,
P.
Muggli
,
N.
Vafaei-Najafabadi
,
D.
Walz
,
G.
White
,
Z.
Wu
,
V.
Yakimenko
, and
G.
Yocky
,
Nature
515
,
92
(
2014
).
7.
X.
Wang
,
R.
Zgadzaj
,
N.
Fazel
,
Z.
Li
,
S. A.
Yi
,
X.
Zhang
,
W.
Henderson
,
Y. Y.
Chang
,
R.
Korzekwa
,
H. E.
Tsai
,
C. H.
Pai
,
H.
Quevedo
,
G.
Dyer
,
E.
Gaul
,
M.
Martinez
,
A. C.
Bernstein
,
T.
Borger
,
M.
Spinks
,
M.
Donovan
,
V.
Khudik
,
G.
Shvets
,
T.
Ditmire
, and
M. C.
Downer
, “
Quasi-monoenergetic laser-plasma acceleration of electrons to 2 GeV
,”
Nat. Commun.
4
(published online 2013).
8.
B. B.
Godfrey
,
J. Comput. Phys.
15
,
504
(
1974
).
9.
R.
Lehe
,
A.
Lifschitz
,
C.
Thaury
,
V.
Malka
, and
X.
Davoine
,
Phys. Rev. Spec. Top. Accel. Beams
16
,
021301
(
2013
).
10.
B. M.
Cowan
,
D. L.
Bruhwiler
,
J. R.
Cary
,
E.
Cormier-Michel
, and
C. G. R.
Geddes
,
Phys. Rev. Spec. Top. Accel. Beams
16
,
041303
(
2013
).
11.
A. D.
Greenwood
,
K. L.
Cartwright
,
J. W.
Luginsland
, and
E. A.
Baca
,
J. Comput. Phys.
201
,
665
(
2004
).
12.
J.-L.
Vay
,
C.
Geddes
,
E.
Cormier-Michel
, and
D.
Grote
,
J. Comput. Phys.
230
,
5908
(
2011
).
13.
J.-L.
Vay
and
A.
Arefiev
,
AIP Conf. Proc.
1777
,
030002
(
2016
).
14.
H.
Vincenti
and
J.-L.
Vay
,
Comput. Phys. Commun.
200
,
147
(
2016
).
15.
F.
Li
,
P.
Yu
,
X.
Xu
,
F.
Fiuza
,
V. K.
Decyk
,
T.
Dalichaouch
,
A.
Davidson
,
A.
Tableman
,
W.
An
,
F. S.
Tsung
,
R. A.
Fonseca
,
W.
Lu
, and
W. B.
Mori
,
Comput. Phys. Commun.
214
,
6
17
(
2017
).
16.
R.
Lehe
,
M.
Kirchen
,
I. A.
Andriyash
,
B. B.
Godfrey
, and
J.-L.
Vay
,
Comput. Phys. Commun.
203
,
66
(
2016
).
17.
R.
Lehe
,
M.
Kirchen
,
S.
Jalas
, and
dornmai
, fbpic/fbpic: 0.2.0 [Data set], Zenodo (
2017
).
18.
I.
Haber
,
R.
Lee
,
H.
Klein
, and
J.
Boris
, in
Proceedings of Sixth Conference on Numerical Simulation Plasmas, Berkeley, CA
(
1973
), pp.
46
48
.
19.
R.
Lehe
,
C.
Thaury
,
E.
Guillaume
,
A.
Lifschitz
, and
V.
Malka
,
Phys. Rev. Spec. Top. Accel. Beams
17
,
121301
(
2014
).
21.
J.-L.
Vay
,
I.
Haber
, and
B. B.
Godfrey
,
J. Comput. Phys.
243
,
260
(
2013
).
22.
M.
Kirchen
,
R.
Lehe
,
B. B.
Godfrey
,
I.
Dornmair
,
S.
Jalas
,
K.
Peters
,
J.-L.
Vay
, and
A. R.
Maier
,
Phys. Plasmas
23
,
100704
(
2016
).
23.
R.
Lehe
,
M.
Kirchen
,
B. B.
Godfrey
,
A. R.
Maier
, and
J.-L.
Vay
,
Phys. Rev. E
94
,
053305
(
2016
).
24.
A.
Lifschitz
,
X.
Davoine
,
E.
Lefebvre
,
J.
Faure
,
C.
Rechatin
, and
V.
Malka
,
J. Comput. Phys.
228
,
1803
(
2009
).
26.
J.
Grebenyuk
,
A. M.
de la Ossa
,
T.
Mehrling
, and
J.
Osterhoff
, in
Proceedings of the first European Advanced Accelerator Concepts Workshop
(
2013
).
J.
Grebenyuk
,
A. M.
de la Ossa
,
T.
Mehrling
, and
J.
Osterhoff
, [
Nucl. Instrum. Methods Phys. Res. Sec. A
740
,
246
(
2014
)].
27.
J.-L.
Vay
,
D. P.
Grote
,
R. H.
Cohen
, and
A.
Friedman
,
Comput. Sci. Discovery
5
,
014019
(
2012
).
28.
S.
Jalas
,
I.
Dornmair
,
R.
Lehe
,
H.
Vincenti
,
J.-L.
Vay
,
M.
Kirchen
, and
A. R.
Maier
, Accurate modeling of plasma acceleration with arbitrary order pseudo-spectral particle-in-cell methods [Data set], Zenodo (
2017
).