X-dispersionless Maxwell solver for plasma-based particle acceleration
Introduction
Plasma-based particle acceleration is a rapidly developing route towards future compact accelerators [1], [2], [3]. The reason is that plasma supports fields orders of magnitude higher than conventional accelerators [5], [6], [4]. Thus, particle acceleration can be accomplished on much shorter distances as compared with the solid-state accelerating structures. However, the plasma is a highly nonlinear medium and requires accurate and computationally efficient numerical modeling to understand and tune the acceleration process. The main workhorse for plasma simulations are Particle-in-Cell codes [7], [8], [9], [10], [11], [12] (a much longer though still incomplete list of PIC codes can be found on the web, see e.g. [13]). These provide the most appropriate description of plasma as an ensemble of particles pushed according to the relativistic equations of motion using self-consistent electromagnetic fields, which are maintained on a spatial grid [14].
From a numerical point of view, plasma-based acceleration represents a classic multi-scale problem. Here, we have the long scale of acceleration distance that can range from centimeters [15] to several meters [16], [17], and the short scale of plasma wavelength that ranges from a few micrometers to near millimeter scales. In addition, if the plasma wave is created by a laser pulse, there is additionally the laser wavelength scale in the sub-micron range. This natural scale disparity makes the simulations of plasma-based acceleration so computationally demanding.
Presently, two types of PIC codes are used to simulate the plasma-based accelerator structures: (i) universal full electromagnetic PIC codes like [7], [8], [9], [10], [11], [12], which solve the unabridged set of Maxwell equations and (ii) quasi-static PIC codes like [20], [8], [18], [19] (and many others), which analytically separate the short scale of plasma wavelength and the long propagation distance scale. The quasi-static PIC codes are proven to be both accurate and very computationally efficient when simulating beam-driven plasma-based wake field acceleration (PBWFA). Unfortunately, the quasi-static approximation for the Maxwell equations eliminates any radiation. Thus, the laser pulse driver has to be described in an envelope approximation [20]. Further, the quasi-static codes fail at simulating sharp plasma boundaries and self-trapping of particles from background plasma.
For this reason, we here consider full electromagnetic (EM) PIC codes which are usually applied for Laser Wake Field Acceleration (LWFA) in plasmas. The full EM PIC correctly describes the laser evolution even in highly nonlinear regimes. The full EM PIC codes are computationally very expensive because they do not separate the different scales.
A significant scale adjustment can be made if one makes a Lorentz transformation of the system into a reference frame moving in the direction of acceleration with a relativistic speed. This leads to the Lorentz contraction of the propagation distance with the relativistic factor , where V is the relative velocity of the reference frame. Simultaneously, the driver - and its wavelength - become longer at nearly the same factor. This so-called “Lorentz-boost” [21] evens the scale disparity and potentially gives a large computational speed up at the cost of not properly resolving backward propagating waves.
However, in “Lorentz-boosted” PIC simulations, the background plasma - both electrons and ions - is moving backward at a relativistic velocity. This moving plasma is a source of free energy that can be easily transformed into high amplitude noise fields. The major numerical mechanism for this parasitic conversion is the Cerenkov resonance [24]. The problem of most existing FDTD Maxwell solvers is that they employ the Yee lattice [25] (with a few exceptions like FBPIC [22] and INF&RNO [23]): individual components of the electromagnetic fields are located at staggered positions in space. The resulting numerical scheme includes a Courant stability restriction on the time step which leads to numerical dispersion. This results in electromagnetic waves with phase velocities below the vacuum speed of light. Thus, the relativistic particles may stay in resonance with the waves and radiate. This non-physical Cerenkov radiation plagues the Lorentz-boosted PIC simulations [26]. Moreover, even normal PIC simulations in the laboratory frame suffer from the numerical Cerenkov effect [27], [28]. Any high density bunch of relativistic particles - e.g. the accelerated witness bunch - emits Cerenkov radiation as well. This affects the bunch energy and emittance [29].
In principle, the Yee scheme can be modified - or extended - by using additional neighboring cells with the goal to tune the numerical dispersion so that the Cerenkov resonance is avoided in the zero order [30], [31]. This reduces the Cerenkov instability, but does not eliminate it. One of the reasons is that the Yee lattice itself is not Lorentz-invariant. The individual field components are located all at different positions staggered in space. In the boosted frame, the fields are Lorentz-transformed and find themselves at the wrong positions. For example, when the boosted frame moves in the X−direction, the pairs and transform one into another. Yet, they are located at different positions within the Yee lattice cell. In addition, the aliasing leads to numerical Cerenkov resonances at wavenumbers from higher Brillouin zones on the numerical grid.
The different positions of the field pairs and on the Yee lattice also cause another problem relevant to the high energy physics. When we want to simulate high current relativistic beams [32], this spatial staggering may lead to a beam numerical self-interaction. A real beam of ultra-relativistic, , particles has a small physical self-interaction due to the difference of these fields with the transverse force . Here is the unit vector in the propagation direction and is the relative difference of the particles longitudinal velocity from the speed of light c. For 50 GeV electrons with , this real difference is as small as . The transverse self-fields and of the ultra-relativistic bunch are also nearly equal with the same miniscule relative difference. However, the Yee lattice defines these fields at staggered positions in space and time. These fields must be interpolated to the same time and to individual particle positions. This interpolation leads to errors and differences between the transverse fields acting on the particle. As a consequence, the bunch self-action due to the numerical errors is many orders of magnitude larger than the real one. This results not only in the bunch numerical self-focusing/defocusing and emittance growth, but also in significant numerical bremsstrahlung and stopping - when these effects are included in the PIC code. A similar inaccuracy can occur in the interaction of laser with a co-propagating relativistic beam (see the appendix in [33]).
We conclude, the Yee lattice is not optimal for simulating high energy applications.
Section snippets
Limitations of pseudo-spectral methods
Recently, pseudo-spectral methods originally proposed by Haber et al. [34], shortly discussed in [14] and used by O. Buneman in his TRISTAN code [35] have seen a remarkable revival [10]. The seeming advantage of the spectral methods is that they are dispersionless and provide an “infinite order” of approximation, even calling the method after Haber “a pseudo-spectral analytical time-domain (PSATD) algorithm” [36].
Indeed, following Sommerfeld [37] we can write the Maxwell equations in the
The general X-dispersionless Maxwell solver
We here develop a FDTD 3D Maxwell solver that has no dispersion for plane waves propagating in vacuum in one selected direction. In plasma-based acceleration this is usually the direction of particle acceleration: the driving laser optical axis. The solver should retain its dispersionless properties not only in vacuum, but also inside dense plasmas, i.e. the optimal time step/grid step relation should not be compromised by the presence of plasma. The solver must not use spectral transformations
The three-dimensional RIP Maxwell solver in Cartesian coordinates
Let us now look at the diffraction/refraction terms. For simplicity, we use Cartesian coordinates.
We project the Yee lattice onto the plane. The grid becomes planar and has the form of Rhombi-in-Plane (RIP), as shown in Fig. 1. The pairs of transverse fields are now combined at positions according to the transport properties (23)-(24). The pair is located at the rhombi vertices . The pair is located at the rhombi vertices . The longitudinal field we
Conservation laws on the RIP grid
The RIP scheme places fields in a transverse plane as seen in Fig. 1. These field locations are perfectly suited for conservative definition of the currents, charges, field divergence and curl on the grid. The simple rule is that the trapezoidal formula must be applied in the longitudinal direction, while in the transverse direction, the usual Yee (spatial leap-frog) formula remains valid.
Dispersion and stability of the RIP scheme
We apply the plane-wave analysis to the marching equations (19), (22) and (31)-(34) with the refraction/diffraction terms (39)-(44) assuming . For simplicity, we assume uniform plasma frequency and the linear current response to the electric field . For the case of interest, , these equations become
Numerical Cerenkov instability test
As a first test, we take the numerical Cerenkov instability. We compare the standard Yee solver, the FFT-based solver (2) and the RIP solver, all implemented on the VLPL platform [8]. No artificial filtering of fields or currents is used. The initial configuration is a ellipsoidal plasma of Gaussian density profile consisting of electrons and protons moving in the X−direction with the average momentum , where α denotes the particle type (), with
Discussion
The new RIP scheme is a compact stencil FDTD Maxwell solver that removes the numerical dispersion in one selected direction. For the waves propagating in the transverse direction, it corresponds to the Yee solver. The RIP scheme is local and does not use any global spectral method. This allows for efficient parallelization via domain decomposition in all three dimensions. The computational cost of the RIP solver is comparable with that of the standard Yee solver. The RIP solver can be used for
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgements
This work has been supported in parts by BMBF project 05K16PFB and DFG PU 213/6-1 (Germany).
References (55)
- et al.
A domain decomposition method for pseudo-spectral electromagnetic simulations of plasmas
J. Comput. Phys.
(2013) - et al.
SMILEI: a collaborative, open-source, multi-purpose particle-in-cell code for plasma simulation
Comput. Phys. Commun.
(2018) - et al.
QUICKPIC: a highly efficient particle-in-cell code for modeling wakefield acceleration in plasmas
J. Comput. Phys.
(2006) - et al.
A spectral, quasi-cylindrical and dispersion-free particle-in-cell algorithm
Comput. Phys. Commun.
(2016) Numerical Cherenkov instabilities in electromagnetic particle codes
J. Comput. Phys.
(1974)- et al.
Numerical methods for instability mitigation in the modeling of laser wakefield accelerators in a Lorentz-boosted frame
J. Comput. Phys.
(2011) - et al.
Suppressing the numerical Cherenkov radiation in the Yee numerical scheme
J. Comput. Phys.
(2016) - et al.
Principles and capabilities of 3-d, E-M particle simulations
J. Comput. Phys.
(1980) - et al.
Convergence in nonlinear laser wakefield accelerators modeling in a Lorentz-boosted frame
Comput. Phys. Commun.
(2019) - et al.
Detailed analysis of the effects of stencil spatial variations with arbitrary high-order finite-difference Maxwell solver
Comput. Phys. Commun.
(2016)
Ultrahigh-order Maxwell solver with extreme scalability for electromagnetic PIC simulations of plasmas
Comput. Phys. Commun.
Comput. Phys. Commun.
Controlling the numerical Cerenkov instability in PIC simulations using a customized finite difference Maxwell solver and a local FFT based current correction
Comput. Phys. Commun.
Improved numerical Cherenkov instability suppression in the generalized PSTD PIC algorithm
Comput. Phys. Commun.
Exact charge conservation scheme for particle-in-cell simulations with an arbitrary form-factor
Comput. Phys. Commun.
A new charge conservation method in electromagnetic particle-in-cell simulations
Comput. Phys. Commun.
Rigorous charge conservation for local electromagnetic field solvers
Comput. Phys. Commun.
The virtual particle electromagnetic particle-mesh method
Comput. Phys. Commun.
Characteristics of an envelope model for laser–plasma accelerator simulation
J. Comput. Phys.
Towards an advanced linear international collider
Laser-driven plasma-wave electron accelerators
Phys. Today
Plasma accelerators
Phys. Rev. Lett.
Rev. Mod. Phys.
One-to-one direct modeling of experiments and astrophysical scenarios: pushing the envelope on kinetic plasma simulations
Plasma Phys. Control. Fusion
Particle-in-cell codes for plasma-based particle acceleration
CERN Yellow Rep.
Cited by (39)
Explicit energy-conserving modification of relativistic PIC method
2024, Journal of Computational PhysicsNumerical dispersion free in longitudinal axis for particle-in-cell simulation
2022, Journal of Computational PhysicsCitation Excerpt :To suppress numerical Cherenkov radiation, a modified Yee scheme has been proposed [18,19], but the simulated laser phase velocity and group velocity were slightly faster than the speed of light. This faster light speed can be solved using direction-splitting scheme [20,21]. Another method uses a pseudo-spectral algorithm [22–24], which solves the Maxwell equations in Fourier space contrast by using a finite-difference algorithm, and can eliminate numerical dispersion.
Undulator design for a laser-plasma-based free-electron-laser
2021, Physics ReportsCitation Excerpt :In many cases these effects can be reduced by modifications of the numerical differential operators [173,174]. In a Lorentz-boosted frame, the NCR growth in relativistic streaming plasma is very fast, and such simulations require the numerical-dispersion-free solvers, e.g. one of the pseudo-spectral analytical time-domain (PSATD) solvers [175–177], or the 1D-dispersionless rhombi-in-plane solver [178]. The presented numerical analysis gives a rather qualitative description of LPA and is not aimed to provide the benchmark parameters for the application design.
Improved modellisation of laser-particle interaction in particle-in-cell simulations
2023, Journal of Plasma PhysicsSpin-polarized <sup>3</sup>He shock waves from a solid-gas composite target at high laser intensities
2024, Plasma Physics and Controlled Fusion