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

Describing complex physics beyond analytical theories requires numerical modeling of the underlying equations in discrete space. In plasma physics, astrophysics, or accelerator physics, Particle-In-Cell1–4 (PIC) methods are commonly used to self-consistently solve the electromagnetic interaction of particle ensembles. A PIC algorithm iteratively solves Maxwell's equations on a discrete grid with particles following the equations of motion in a continuous space.

Some of the physical systems accessible with the PIC method feature plasmas drifting at relativistic velocities, for example, when modeling plasma-based particle accelerators5 in the optimal frame of reference6 or astrophysical plasma interactions.7 In those cases, the applicability of the to-date electromagnetic PIC algorithms is fundamentally limited by the numerical Cherenkov instability8–11 (NCI), which either falsifies the numerical results or causes virulent growth of unphysical waves.

Here, we present a novel discrete formulation of the fundamental kinetic equations of plasmas, i.e., the Maxwell's and Newton-Lorentz equations, that represents the physics in a moving Galilean frame of reference and thereby is intrinsically free of the NCI for plasmas drifting at uniform relativistic velocities.

The NCI originates from the coupling of distorted electromagnetic modes with spurious particle modes. Distortions of the electromagnetic field modes are caused by numerical inaccuracies of the discretized field-solving algorithm. Spurious spatial and temporal aliases of the physical particle modes result from the numerical mismatch of sampling the continuously distributed particle quantities to the discrete field grid. To the first order, for example, numerical Cherenkov radiation (NCR) can occur, if the dispersion relation of the electromagnetic waves is numerically distorted. In this case, particles moving at relativistic velocities vp ≈ c couple resonantly to electromagnetic waves of high frequency, which propagate at a spurious phase velocity vΦ < vp < c, causing Cherenkov-like radiation to be emitted. Although many algorithms, such as pseudo-spectral solvers,12 do not suffer from NCR, higher order NCI effects severely limit the stable modeling of relativistic plasmas.

So far, no electromagnetic, fully explicit PIC algorithm is intrinsically free of NCI, even for the simple case of a plasma drifting at a uniform relativistic velocity. Previously developed suppression strategies can limit the NCI growth rate, thereby retaining the physical meaningfulness of a simulation. For example, wide-band smoothing13–15 or damping16 of the currents or electromagnetic fields can hinder the development of the instability. Coupling of unphysical modes can be mitigated by slightly changing the ratio of the electric and magnetic fields as seen by the particles,17–19 by scaling the deposited currents with a frequency-dependent factor,20,21 or by artificially modifying the physical electromagnetic dispersion relation.22–24 Yet, all of these techniques rely on numerical methods that potentially alter the physics and could affect the results obtained with the algorithm.

In contrast, the method presented in this paper inherently eliminates the NCI for a relativistically drifting plasma, as opposed to suppressing its growth by the measures described above. From a heuristic point of view, the main difference between modeling a plasma at rest, showing no NCI, and a relativistically drifting plasma is that the particles move with respect to the static numerical grid. Thus, intuitively, by mathematically representing the underlying discrete equations such that this discrepancy in relative movement is eliminated, the NCI should be suppressed.

This is achieved by applying a Galilean coordinate transformation of the form

x=xvgalt
(1)

to the frame of reference in which a plasma is moving at a relativistic velocity. Consequently, the equations of motion and Maxwell's equations transform to

dxdt=pγmvgal,dpdt=q(E+pγm×B),(tvgal·)B=×E,1c2(tvgal·)E=×Bμ0j,

and the continuity equation becomes (tvgal·)ρ+·j=0. Here, denotes the spatial derivative with respect to the Galilean coordinates x. For vgal=0, these equations reduce to their well-known original form.

Using the Pseudo-Spectral Analytical Time Domain12 (PSATD) framework, the last two equations are transformed to Fourier space and can then be integrated analytically in time. As the quantities are only known at discrete times in a PIC algorithm, the time evolution of ρ and j needs to be explicitly taken into account during integration. Typically, the currents are assumed to be constant over one time step Δt in the original coordinates x. A key difference of our new scheme is that we assume the currents to be co-moving with respect to the original coordinates x, and hence constant over one time step in the Galilean coordinates x. The resulting Galilean-PSATD equations for the advance of the spectral field components, Ê and B̂, from time step nΔt to (n+1)Δt are then given by (see Ref. 25 for a derivation)

B̂n+1=θ2CB̂nθ2Sckik×Ên+θχ1ϵ0c2k2ik×Ĵn+1/2,Ên+1=θ2CÊn+θ2Skcik×B̂n+iνθχ1θ2Sϵ0ckĴn+1/21ϵ0k2(χ2ρ̂n+1θ2χ3ρ̂n)ik,C=cos(ckΔt),S=sin(ckΔt),k=|k|,ν=k·vgalck,θ=eik·vgalΔt/2,θ*=eik·vgalΔt/2,χ1=11ν2(θ*Cθ+iνθS),χ2=χ1θ(1C)θ*θ,χ3=χ1θ*(1C)θ*θ,

where k is the wavevector. The currents Ĵ at time (n+1/2)Δt and the charge density ρ̂ at time nΔt and (n+1)Δt are generated by the particles and deposited to the grid nodes before being transformed to Fourier space. Subsequently, the updated fields are transformed back to real space and interpolated to the particles, which are then advanced in time using the Galilean transformed equations of motion.

This algorithm allows to model a plasma moving at vp in a co-propagating set of coordinates x with vgal=vp. As shown in Fig. 1, the flowing plasma particles now remain static with respect to the numerical grid. Because of this and the co-moving current assumption, the NCI is completely eliminated for particles streaming at the velocity vp.

FIG. 1.

Schematic drawing illustrating the Galilean concept. Without applying a Galilean coordinate transformation to the Particle-In-Cell equations (Standard), a plasma flowing with velocity vp in z (represented by a single particle) would propagate a distance vpΔt with respect to the numerical grid (represented by a single cell) during one time step Δt. However, in a Galilean transformed discrete space x with vgal=(0,0,vgal=vp), the plasma particles remain static with respect to the discrete grid nodes, which themselves propagate a distance z + vgalΔt in the original coordinate system x.

FIG. 1.

Schematic drawing illustrating the Galilean concept. Without applying a Galilean coordinate transformation to the Particle-In-Cell equations (Standard), a plasma flowing with velocity vp in z (represented by a single particle) would propagate a distance vpΔt with respect to the numerical grid (represented by a single cell) during one time step Δt. However, in a Galilean transformed discrete space x with vgal=(0,0,vgal=vp), the plasma particles remain static with respect to the discrete grid nodes, which themselves propagate a distance z + vgalΔt in the original coordinate system x.

Close modal

The algorithm is implemented in the WARP26 code, for Cartesian coordinates, as well as in the recently developed quasi-cylindrical27 code FBPIC.28 In the study of Lehe et al.,25 we also present an analytical derivation of the dispersion relation and conduct a detailed empirical and theoretical stability analysis for uniform relativistically flowing plasmas. Here, we restrict ourselves to presenting the general concept and the practical demonstration of the stability and accuracy of our new method with a direct application. In the following, Lorentz-boosted frame simulations of plasma acceleration with FBPIC are presented.

Plasma-based accelerators5 can sustain high field gradients, allowing for the acceleration of charged particles within distances shorter by orders of magnitude compared to conventional accelerators. In a plasma-wakefield accelerator, an intense driver beam (a high intensity laser pulse or particle bunch) propagates through an underdense plasma and induces a charge separation on the sub-mm scale. This leads to the excitation of a trailing density wave carrying large electric fields, suitable for the acceleration of electron bunches to high energies.

The natural frame of reference for PIC simulations of plasma accelerators is the laboratory frame. In this frame of reference, the physical objects of small scale, i.e., the laser or particle beam, propagate at relativistic velocities in a single direction while interacting with a large scale object that is static, i.e., the plasma. A Lorentz transformation in the propagation direction of the driver beam then relaxes the requirements on the spatial resolution while contracting the required simulation distance. In this Lorentz-boosted frame,6 the co-propagating quantities, e.g., the laser or the plasma wavelength, are elongated by γ(1 + β), whereas the previously static lengths, such as the plasma, are contracted by γ and counter-propagate with a relativistic velocity –βc. Thereby, a speed-up by orders of magnitude can be achieved that scales as γboost2, with the maximum speed-up typically limited to 2γwake2, i.e., the phase velocity of the wake, in the case of laser-plasma acceleration.

In the following, we show simulations of a non-linear plasma wave driven by a Gaussian laser pulse with wavelength λ = 800 nm, peak normalized vector potential a0 = 1.5, pulse length  = 8 µm, and waist w0 = 30 µm ( and w0 correspond to 2σ of the longitudinal and transverse laser intensity distribution, respectively) that propagates through a matched plasma guiding channel with an on-axis electron density ne=1·1018cm3. In the generated wakefields, a 1 pC electron bunch of size σz = 1 µm, σr = 2 µm, and normalized emittance ϵn = 0.5 mm mrad, located at the back of the first wave bucket, is accelerated from 100 MeV to 687 MeV within a propagation distance of zprop ≈ 14.3 mm. The resolution of the simulation is 32 cells per λ in the longitudinal and 2 cells per µm in the transverse direction. Third order particle shapes are used with 24 particles per cell. The time step is set to Δt = Δz/c. Please note that the simulation box and the physical quantities in the propagation direction are equally elongated in the Lorentz-boosted frame. Therefore, the number of cells in the box, as well as the longitudinal resolution with respect to the laser wavelength λ, is the same in the lab and the Lorentz-boosted frame simulation.

As described above, the occurrence of the NCI, caused by the counter-streaming relativistic plasma, can hinder the application of the Lorentz-boosted frame method for simulations of plasma-wakefield accelerators. With our new method, however, such a simulation is modeled in a Galilean transformed coordinate system that counter-propagates to the Lorentz-boosted frame with the velocity vgal = –βc in the direction of the boosted plasma. With respect to the numerical grid, the background plasma is thus static, whereas the elongated quantities, such as the laser pulse and the electron bunch, propagate with a velocity increased by the same amount with respect to the grid.

Fig. 2 shows the charge density obtained in a Lorentz boosted frame (γboost = 13) with boosted longitudinal coordinate zboost = γ(zvt). The upper-half of the plot shows the results of a simulation with the Galilean-PSATD solver, whereas the lower-half shows the corresponding results of a standard PSATD simulation. Here, a fast growing NCI can be observed. In contrast, the same simulation remains completely stable when modeled in the Galilean transformed discrete space. We emphasize that all numerical parameters are the same in these simulations, except for the difference in using vgal = –βc instead of vgal = 0 for the Galilean-PSATD equations. Thus, the absence of the instability results solely from the Galilean transformation of the underlying discrete equations. Even though the electron bunch and the grid move in opposite directions, we do not observe any NCI around the bunch. This can be explained by the fact that the electron bunch has a density that is much lower than the plasma, as it is elongated in the Lorentz-boosted frame. Moreover, due to its non-zero charge, it is probably much less affected by higher order numerical Cherenkov effects. Likewise, in a laboratory frame simulation, a relativistic electron bunch does typically not lead to an instability, as long as NCR is suppressed.29 

FIG. 2.

Charge density ρ obtained from a Lorentz-boosted frame simulation (γboost = 13) of a non-linear laser-plasma wave. At the time step shown, a part of the laser pulse, propagating to the right, has already left the plasma, which is flowing to the left. The upper half corresponds to a Galilean-PSATD simulation with vgal = –βc, showing no instability. The lower half shows the same simulation, mirrored along the x = 0 axis, but conducted with the standard PSATD solver. Here, a fast growing, virulent NCI can be observed.

FIG. 2.

Charge density ρ obtained from a Lorentz-boosted frame simulation (γboost = 13) of a non-linear laser-plasma wave. At the time step shown, a part of the laser pulse, propagating to the right, has already left the plasma, which is flowing to the left. The upper half corresponds to a Galilean-PSATD simulation with vgal = –βc, showing no instability. The lower half shows the same simulation, mirrored along the x = 0 axis, but conducted with the standard PSATD solver. Here, a fast growing, virulent NCI can be observed.

Close modal

In order to validate the accuracy of our new method, results from the stable Lorentz-boosted frame simulation are compared with a laboratory frame simulation. Fig. 3 shows the electric fields at the end of the acceleration distance. The upper half of the plots shows the results of the reference simulation (γboost = 1), whereas the lower half shows the corresponding back-transformed results of the Lorentz-boosted frame simulation (γboost = 13). Both the longitudinal fields Ez and the transverse fields Ey show no differences. Note that the results in the Lorentz-boosted frame are obtained within only a few thousand time steps, whereas the lab frame simulation takes more than half a million time steps to complete. We achieve a speed-up of ≈287 (≈92% of the optimal speed-up) with FBPIC, where the only overhead is the on-the-fly back-transformation of data to the laboratory frame.

FIG. 3.

Comparison of the accelerating fields (Ez fields and on-axis lineout) and the focusing fields (Ey fields and off-axis, y = 11.25 µm, lineout). The upper half of the plots shows the results of a laboratory frame simulation (solid line) with γboost = 1. The lower half shows the back-transformed results of a Lorentz-boosted frame simulation (dashed line) with γboost = 13, mirrored along the x, y = 0 axis.

FIG. 3.

Comparison of the accelerating fields (Ez fields and on-axis lineout) and the focusing fields (Ey fields and off-axis, y = 11.25 µm, lineout). The upper half of the plots shows the results of a laboratory frame simulation (solid line) with γboost = 1. The lower half shows the back-transformed results of a Lorentz-boosted frame simulation (dashed line) with γboost = 13, mirrored along the x, y = 0 axis.

Close modal

Furthermore, we compare characteristic bunch and laser parameters to demonstrate that the physics is preserved in the Lorentz-boosted frame. Fig. 4 shows the evolution of the laser waist w0 and the pulse duration τ, as well as the kinetic energy Ekin, the rms energy spread σE, and the normalized emittance ϵn of the accelerated electron bunch. During the propagation through the plasma guiding channel, the laser pulse self-focuses transversely and the pulse duration shortens due to the relativistic interaction with the plasma. The electron bunch is initially situated at the minimum of the accelerating field and slips towards the laser pulse during the propagation. It is accelerated to 687 MeV, while accumulating an rms energy spread of ≈11.5%, due to the slope of the accelerating field. As the bunch enters the plasma, strong transverse fields act on it abruptly, causing transverse oscillations of the bunch size and growth of the emittance ϵn to around 1.4 mm mrad. In direct comparison with the laboratory simulation, all the quantities shown differ only on the sub-percent level at the end of the propagation distance, which resembles a remarkable precision.

FIG. 4.

Comparison of the laser and electron bunch evolution between the laboratory frame (solid line) and the Lorentz-boosted frame (dashed line) simulation. The upper plot shows the pulse duration τ (blue) and the laser waist w0 (red), and the lower plot shows the kinetic Energy Ekin (red), the rms energy spread σE (gray area), and the normalized emittance ϵn (blue) over the complete acceleration distance of zprop ≈ 14.3 mm.

FIG. 4.

Comparison of the laser and electron bunch evolution between the laboratory frame (solid line) and the Lorentz-boosted frame (dashed line) simulation. The upper plot shows the pulse duration τ (blue) and the laser waist w0 (red), and the lower plot shows the kinetic Energy Ekin (red), the rms energy spread σE (gray area), and the normalized emittance ϵn (blue) over the complete acceleration distance of zprop ≈ 14.3 mm.

Close modal

In conclusion, we have proposed a novel discrete formulation of the fundamental kinetic equations of plasmas in Galilean transformed coordinates. To the best of our knowledge, we thereby derived the first electromagnetic, fully explicit and self-consistent PIC representation that is intrinsically free of the NCI for plasmas flowing at a uniform velocity. Our concept is not reliant on otherwise inevitable numerical corrections and, unlike most of the previous NCI suppression strategies, it is independent of the specific geometry. This allows us to combine the accuracy and efficiency of a spectral, quasi-3D PIC algorithm with the superior stability properties of the presented method. Applying the Galilean scheme to simulations of plasma accelerators in the Lorentz-boosted frame yields excellent agreement, while achieving a close-to-optimal speed-up of more than two orders of magnitude in practice.

The Galilean-PSATD solver satisfies a combination of specific requirements that cause the complete absence of NCI for a plasma drifting at uniform relativistic velocity, which can be summarized qualitatively as follows: It is free of spurious numerical dispersion; thus, first order NCR effects are suppressed. Furthermore, the analytic integration in the time domain in combination with the co-moving current assumption explicitly retains stability properties that are otherwise lost when directly discretizing the Maxwell equations. Finally, the suppressed movement of the particles with respect to the numerical grid eliminates aliasing effects due to sampling the particle fields to the grid and vice versa. The Galilean scheme is therefore not directly applicable to other solvers, such as standard Finite-Difference Time Domain (FDTD) solvers. It is unclear whether it would be possible to derive an implementation using other solvers (e.g., FDTD) that would retain the advantages obtained when integrating with the Galilean-PSATD solver. This requires further studies.

Nevertheless, it is possible to apply arbitrary-order spectral methods4,24,30 to the Galilean-PSATD solver, which enables the efficient parallelization based on spatial domain decomposition31 employing domain-local Fourier transforms of the fields. These topics will be covered by future research, as well as the potential generalization to support arbitrary relativistic plasma flows. For example, the new method could directly be extended to model collisionless astrophysical shocks7 involving two plasmas, by employing separate numerical grids for each plasma using different Galilean transformed coordinates. Taking advantage of the superposition principle, only the electromagnetic fields would be shared between those individual grids.

We gratefully acknowledge the computing time provided on the supercomputer JURECA under project HHH20 and on the PHYSnet cluster of the University of Hamburg. For verification of the method in Cartesian geometry with Warp, this research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. 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 Berkeley Lab.

1.
O.
Buneman
,
C. W.
Barnes
,
J. C.
Green
, and
D. E.
Nielsen
,
J. Comput. Phys.
38
,
1
(
1980
).
2.
J. M.
Dawson
,
Rev. Mod. Phys.
55
,
403
(
1983
).
3.
R. W.
Hockney
and
J. W.
Eastwood
,
Computer Simulation Using Particles
(
CRC Press
,
1988
).
4.
C. K.
Birdsall
and
A. B.
Langdon
,
Plasma Physics via Computer Simulation
(
CRC Press
,
2004
).
5.
E.
Esarey
,
C. B.
Schroeder
, and
W. P.
Leemans
,
Rev. Mod. Phys.
81
,
1229
(
2009
).
7.
L.
Sironi
,
A.
Spitkovsky
, and
J.
Arons
,
Astrophys. J.
771
,
54
(
2013
).
8.
B. B.
Godfrey
,
J. Comput. Phys.
15
,
504
(
1974
).
9.
B. B.
Godfrey
,
J. Comput. Phys.
19
,
58
(
1975
).
10.
B. B.
Godfrey
and
J.-L.
Vay
,
J. Comput. Phys.
248
,
33
(
2013
).
11.
X.
Xu
,
P.
Yu
,
S. F.
Martins
,
F. S.
Tsung
,
V. K.
Decyk
,
J.
Vieira
,
R. A.
Fonseca
,
W.
Lu
,
L. O.
Silva
, and
W. B.
Mori
,
Comput. Phys. Commun.
184
,
2503
(
2013
).
12.
I.
Haber
,
R.
Lee
,
H.
Klein
, and
J.
Boris
,
Proc. Sixth Conf. on Num. Sim. Plasmas, Berkeley, CA
(
1973
), pp.
46
48
.
13.
J.-L.
Vay
,
C. G. R.
Geddes
,
C.
Benedetti
,
D. L.
Bruhwiler
,
E.
Cormier-Michel
,
B. M.
Cowan
,
J. R.
Cary
, and
D. P.
Grote
,
AIP Conf. Proc.
1299
,
244
249
(
2010
).
14.
J. L.
Vay
,
C. G. R.
Geddes
,
E.
Cormier-Michel
, and
D. P.
Grote
,
J. Comput. Phys.
230
,
5908
(
2011
).
15.
J.-L.
Vay
,
C. G. R.
Geddes
,
E.
Cormier-Michel
, and
D. P.
Grote
,
Phys. Plasmas (1994-present)
18
,
030701
(
2011
).
16.
S. F.
Martins
,
R. A.
Fonseca
,
L. O.
Silva
,
W.
Lu
, and
W. B.
Mori
,
Comput. Phys. Commun.
181
,
869
(
2010
).
17.
B. B.
Godfrey
and
J.-L.
Vay
,
J. Comput. Phys.
267
,
1
(
2014
).
18.
B. B.
Godfrey
, e-print arXiv:1408.1146 [physics].
19.
B. B.
Godfrey
and
J.-L.
Vay
,
Comput. Phys. Commun.
196
,
221
(
2015
).
20.
B. B.
Godfrey
,
J.-L.
Vay
, and
I.
Haber
,
J. Comput. Phys.
258
,
689
(
2014
).
21.
B. B.
Godfrey
,
J. L.
Vay
, and
I.
Haber
,
IEEE Trans. Plasma Sci.
42
,
1339
(
2014
).
22.
P.
Yu
,
X.
Xu
,
V. K.
Decyk
,
F.
Fiuza
,
J.
Vieira
,
F. S.
Tsung
,
R. A.
Fonseca
,
W.
Lu
,
L. O.
Silva
, and
W. B.
Mori
,
Comput. Phys. Commun.
192
,
32
(
2015
).
23.
P.
Yu
,
X.
Xu
,
A.
Tableman
,
V. K.
Decyk
,
F. S.
Tsung
,
F.
Fiuza
,
A.
Davidson
,
J.
Vieira
,
R. A.
Fonseca
,
W.
Lu
,
L. O.
Silva
, and
W. B.
Mori
,
Comput. Phys. Commun.
197
,
144
(
2015
).
24.
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
, e-print arXiv:1605.01496 [physics].
25.
R.
Lehe
,
M.
Kirchen
,
B. B.
Godfrey
,
A. R.
Maier
, and
J.-L.
Vay
, e-print arXiv:1608.00227 [physics].
26.
J.-L.
Vay
,
D. P.
Grote
,
R. H.
Cohen
, and
A.
Friedman
,
Comput. Sci. Discovery
5
,
014019
(
2012
).
27.
A. F.
Lifschitz
,
X.
Davoine
,
E.
Lefebvre
,
J.
Faure
,
C.
Rechatin
, and
V.
Malka
,
J. Comput. Phys.
228
,
1803
(
2009
).
28.
R.
Lehe
,
M.
Kirchen
,
I. A.
Andriyash
,
B. B.
Godfrey
, and
J.-L.
Vay
,
Comput. Phys. Commun.
203
,
66
(
2016
).
29.
R.
Lehe
,
A.
Lifschitz
,
C.
Thaury
,
V.
Malka
, and
X.
Davoine
,
Phys. Rev. Spec. Top. – Accel. Beams
16
,
021301
(
2013
).
30.
H.
Vincenti
and
J. L.
Vay
,
Comput. Phys. Commun.
200
,
147
(
2016
).
31.
J.-L.
Vay
,
I.
Haber
, and
B. B.
Godfrey
,
J. Comput. Phys.
243
,
260
(
2013
).