X-dispersionless Maxwell solver for plasma-based particle acceleration

https://doi.org/10.1016/j.jcp.2020.109622Get rights and content

Highlights

  • The new Maxwell solver completely removes the numerical dispersion along one selected direction.

  • It accurately simulates the electromagnetic pulse down to its phase

  • the new solver removes greatly reduces the numerical Cerenkov instability completely.

  • The new solver is local and thus is infinitely scalable via simple domain decomposition.

Abstract

A semi-implicit finite difference time domain (FDTD) numerical Maxwell solver is developed for full electromagnetic Particle-in-Cell (PIC) codes for the simulations of plasma-based acceleration. The solver projects the volumetric Yee lattice into planes transverse to a selected axis (the particle acceleration direction). The scheme - by design - removes the numerical dispersion of electromagnetic waves running parallel the selected axis. The fields locations in the transverse plane are selected so that the scheme is Lorentz-invariant for relativistic transformations along the selected axis. The solver results in “Galilean shift” of transverse fields by exactly one cell per time step. This eases greatly the problem of numerical Cerenkov instability (NCI). The fields positions build rhombi in plane (RIP) patterns. The RIP scheme uses a compact local stencil that makes it perfectly suitable for massively parallel processing via domain decomposition along all three dimensions. No global/local spectral methods are involved.

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 γ=1/1V2/c2, 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 Ey,Bz and Ez,By 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 Ey,Bz and Ez,By 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, γ1, particles has a small physical self-interaction due to the difference of these fields with the transverse force q(E+β||e||×B). Here e|| is the unit vector in the propagation direction and 1β||=1v||/c1/2γ2 is the relative difference of the particles longitudinal velocity v|| from the speed of light c. For 50 GeV electrons with γ105, this real difference is as small as 1β||51011. The transverse self-fields E and B 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 (Y,Z) 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 Ey,Bz is located at the rhombi vertices (i,j+1/2,k). The pair Ez,By is located at the rhombi vertices (i,j,k+1/2). The longitudinal field Ex 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 F=F˜exp(iωt+ikr). For simplicity, we assume uniform plasma frequency ωp2=4πne2/γ and the linear current response to the electric field 2cΔsinωτ2J˜=iq2nE˜. For the case of interest, cτ=hx=Δ, these equations become2ΔsinωΔ2ccoskxΔ2E˜y=2ΔsinkxΔ2cosωΔ2cB˜z+2hzsinkzhz2coskxΔ2B˜x+ωp2Δ2csinωΔ2cE˜y2ΔsinωΔ2ccoskxΔ2E˜z=2ΔsinkxΔ2cosωΔ2cB˜y2hysinkyhy2

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 n=n0exp(r2/σ2) consisting of electrons and protons moving in the X−direction with the average momentum <p0>/mαc=(px0,0,0), where α denotes the particle type (α=e,p), with mp/me=1846

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)

  • H. Vincenti et al.

    Ultrahigh-order Maxwell solver with extreme scalability for electromagnetic PIC simulations of plasmas

    Comput. Phys. Commun.

    (July 2018)
  • B.B. Godfrey et al.

    Comput. Phys. Commun.

    (2015)
  • Fei Li et al.

    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.

    (2017)
  • Brendan B. Godfrey et al.

    Improved numerical Cherenkov instability suppression in the generalized PSTD PIC algorithm

    Comput. Phys. Commun.

    (2015)
  • T.Z. Esirkepov

    Exact charge conservation scheme for particle-in-cell simulations with an arbitrary form-factor

    Comput. Phys. Commun.

    (2001)
  • T.Y. Umeda et al.

    A new charge conservation method in electromagnetic particle-in-cell simulations

    Comput. Phys. Commun.

    (2003)
  • John Villasenor et al.

    Rigorous charge conservation for local electromagnetic field solvers

    Comput. Phys. Commun.

    (1992)
  • James W. Eastwood

    The virtual particle electromagnetic particle-mesh method

    Comput. Phys. Commun.

    (1991)
  • Benjamin M. Cowan et al.

    Characteristics of an envelope model for laser–plasma accelerator simulation

    J. Comput. Phys.

    (2011)
  • Towards an advanced linear international collider

  • W. Leemans et al.

    Laser-driven plasma-wave electron accelerators

    Phys. Today

    (2009)
  • C. Joshi et al.

    Plasma accelerators

  • T. Tajima et al.

    Phys. Rev. Lett.

    (1979)
  • E. Esarey et al.

    Rev. Mod. Phys.

    (2009)
  • R.A. Fonseca et al.

    One-to-one direct modeling of experiments and astrophysical scenarios: pushing the envelope on kinetic plasma simulations

    Plasma Phys. Control. Fusion

    (2008)
  • A. Pukhov

    Particle-in-cell codes for plasma-based particle acceleration

    CERN Yellow Rep.

    (2016)
  • Cited by (39)

    • Numerical dispersion free in longitudinal axis for particle-in-cell simulation

      2022, Journal of Computational Physics
      Citation 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 Reports
      Citation 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.

    View all citing articles on Scopus
    View full text