X-Dispersionless Maxwell solver for plasma-based particle acceleration

TL;DR

This paper introduces a semi-implicit Maxwell solver for plasma-based particle acceleration that eliminates numerical dispersion and Cerenkov instability, enabling efficient, Lorentz-invariant simulations suitable for parallel computing.

Contribution

A novel semi-implicit Maxwell solver projecting onto transverse planes to remove dispersion and instability in plasma acceleration simulations.

Findings

Eliminates numerical dispersion of electromagnetic waves.

Removes numerical Cerenkov instability (NCI).

Enables efficient parallel computation with local stencil.

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) that makes the scheme quasi-one-dimensional. The scheme 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 removes 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.

Equations

$$\frac{\partial \hat{\mathbf{F}}}{\partial t}$$

$${\hat{\mathbf{F}}}^{n+1}$$

$${\widehat{F_{||}}}^{n+1}$$

$${\hat{\mathbf{F}}}_{\perp }^{n+1}$$

$$\frac{1}{c}\frac{\partial E_x}{\partial t}$$

$$\frac{1}{c}\frac{\partial E_y}{\partial t}$$

$$\frac{1}{c}\frac{\partial E_z}{\partial t}$$

$$\frac{1}{c}\frac{\partial B_x}{\partial t}$$

$$\frac{1}{c}\frac{\partial B_y}{\partial t}$$

$$\frac{1}{c}\frac{\partial B_z}{\partial t}$$

$${\Gamma }_x$$

$${\Gamma }_y$$

$${\Gamma }_z$$

$${\Phi }_x$$

$${\Phi }_y$$

$${\Phi }_z$$

$$\frac{E_{y(i+1)}^{n+1}+E_{y(i)}^{n+1}-E_{y(i+1)}^n-E_{y(i)}^n}{2c\tau }$$

$$\frac{E_{z(i+1)}^{n+1}+E_{z(i)}^{n+1}-E_{z(i+1)}^n-E_{z(i)}^n}{2c\tau }$$

$$\frac{E_{x(i)}^{n+1}-E_{x(i)}^n}{c\tau }$$

$$\frac{B_{y(i)}^{n+1}+B_{y(i+1)}^{n+1}-B_{y(i)}^n-B_{y(i+1)}^n}{2c\tau }$$

$$\frac{B_{z(i)}^{n+1}+B_{z(i+1)}^{n+1}-B_{z(i)}^n-B_{z(i+1)}^n}{2c\tau }$$

$$\frac{B_{x(i)}^{n+1}-B_{x(i)}^n}{c\tau }$$

$$T_{y(i)}^{+(n+1)}=E_{y(i)}^{n+1}+B_{z(i)}^{n+1}$$

$$T_{z(i)}^{+(n+1)}=E_{z(i)}^{n+1}+B_{y(i)}^{n+1}$$

$$T_{y(i)}^{+(n+1)}$$

$$T_{z(i)}^{+(n+1)}$$

$$T_{y(i)}^{-(n+1)}=E_{y(i)}^{n+1}-B_{z(i)}^{n+1}$$

$$T_{z(i)}^{-(n+1)}=E_{z(i)}^{n+1}-B_{y(i)}^{n+1}$$

$$T_{y(i)}^{-(n+1)}$$

$$T_{z(i)}^{-(n+1)}$$

$$E_{y(i)}^{n+1}$$

$$E_{z(i)}^{n+1}$$

$$B_{y(i)}^{n+1}$$

$$B_{z(i)}^{n+1}$$

$$E_{y(i)}^{(n+1)}$$

$$E_{z(i)}^{(n+1)}$$

$$B_{y(i)}^{(n+1)}$$

$$B_{z(i)}^{(n+1)}$$

$${\Gamma }_{y(i+1/2,j+1/2,k)}^{n+1/2}$$

$${\Gamma }_{z(i+1/2,j,k+1/2)}^{n+1/2}$$

$${\Gamma }_{x(i,j,k)}^{n+1/2}$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Plasma-based particle acceleration is a rapidly developing route towards future compact accelerators [1, 2, 3, 4]. The reason is that plasma supports fields orders of magnitude higher than conventional accelerators [5, 6]. Thus, particle acceleration can be accomplished on much shorter distances as compared with the solid-state accelerating strucutures. 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 longerthough 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 aproximation 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 $\gamma=1/\sqrt{1-V^{2}/c^{2}}$, 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 $E_{y},B_{z}$ and $E_{z},B_{y}$ 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 $E_{y},B_{z}$ and $E_{z},B_{y}$ 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, $\gamma\gg 1$, particles has a small physical self-interaction due to the difference of these fields with the transverce force $q\left(\mathbf{E_{\perp}+\beta_{||}}\mathbf{e_{||}}\times\mathbf{B}_{\perp}\right)$. Here $\mathbf{e_{||}}$ is the unit vector in the propagation direction and $1-\beta_{||}=1-v_{||}/c\approx 1/2\gamma^{2}$ is the relative difference of the particles longitudinal velocity $v_{||}$ from the speed of light $c$. For $50\,$GeV electrons with $\gamma\approx 10^{5}$, this real difference is as small as $1-\beta_{||}\approx 5\cdot 10^{-11}$. The transverse self-fields ${\bf E}_{\perp}$ and ${\bf B}_{\perp}$ 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 transvers 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.

2 Limitations of pseudo-spectral methods

Recently, pseudo-spectral methods origninally 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 Fourier space as

[TABLE]

where $\hat{\mathbf{J}}=FFT[\mathbf{J}]$ is the Fourier image of the real current while $\hat{\mathbf{F}}=FFT[\mathbf{F}]$ is the Fourier image of the complex electromagnetic field $\mathbf{F}=\mathbf{E}+i\mathbf{B}$. It is straightforward to show that the numerical scheme advancing the fields from the time step $n$ to $n+1$ in the form

[TABLE]

is dispersionless in vacuum and provides second order approximation for the plasma currents. Here, $C_{\mathbf{k}}=\cos(ck\tau)$, $\tilde{C}_{\mathbf{k}}=\cos(ck\tau/2)$, $\mathbf{S}_{\mathbf{k}}=\mathbf{e_{k}}\sin(ck\tau)$, $\tilde{\mathbf{S}}_{\mathbf{k}}=\mathbf{e_{k}}\sin(ck\tau/2)$, $\tau$ is the time step, $k=|\mathbf{k}|$, $\mathbf{e_{k}}=\mathbf{k}/k$.

The FFT-based solvers are intrinsically global. This means, they need information about fields in the full simulation domain to update the local field at a particular point in space. This contradicts the causality principle of the special relativity: only fields within $c\tau$ distance from the space point may cause the local fields to change. The propagator (2) explicitly separates the fields and currents into the propagating transverse fields and non-propagating longitudinal fields.

Indeed, the longitudinal part of the field is ${\bf F}_{||}^{n}=\left(\mathbf{e_{k}}\cdot\hat{\mathbf{F}}^{n}\right)\mathbf{e_{k}}$, and the transverse part is ${\bf F}_{\perp}^{n}={\bf F}^{n}-\left(\mathbf{e_{k}}\cdot\hat{\mathbf{F}}^{n}\right)\mathbf{e_{k}}$, so that ${\bf F}={\bf F_{||}}+{\bf F_{\perp}}$. The same is valid for the current ${\bf J}={\bf J_{||}}+{\bf J_{\perp}}$. The Eq. (2) projected onto the vector ${\bf e_{k}}$ is

[TABLE]

The transverse field components are updated according to

[TABLE]

We see that the transverse fields (4) propagate with the speed of light. The longitudinal component (3) does not propagate anywhere. Taking divergence of (3) we arrive at the Poisson equation.

Computationally, one can use spectral algorithms that are ”local” to one simulation sub-domain [38, 39]. In this case, an ultra-high order finite differences scheme can be designed. The resulting convolution is then effciently computed within one sub-domain with the help of a spectral transformation. The resulting schemes show excellent parallel scalability [39]

Although the spectral solver (2) removes the numerical dispersion, it does not remove aliasing errors and the numerical Cerenkov instability in pseudo-spectral codes persists, even when at a lower rate [40, 41]. In an effort to remove the Cerenkov instability in pseudo-spectral codes, filtering currents of the most unstable modes is often applied [42, 43]. This artificial filtering, however, may lead to additional unphysical effects in the pseudo-spectral simulations.

Another approach is elimination of numerical Cherenkov instability in flowing-plasma particle-in-cell simulations by using Galilean coordinates [44]. This approach removes relative motion between the numerical grid and the streaming plasma: the grid cells flow together with the plasma. Unfortunately, this trick works only in one direction and does not help removing numerical Cerenkov emission of the high current bunch being accelerated in the opposite direction.

We conclude that pseudo-spectral methods are far from ideal candidates for PIC simulations and that a better FDTD method is required. In this work, a new FDTD solver is presented that does not employ spectral transformations and yet has the unique property of having no numerical dispersion along one selected spatial axis. The positions of the transverse field pairs $\left(E_{y},B_{z}\right)$ $\left(E_{z},B_{y}\right)$ are colocated in the RIP scheme. This is Lorentz-invariant and greatly improves the accuracy in calculating the transverse force acting on a relativistic particle moving along the $X-$axis.

3 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 and should have a compact local stencil. This is the pre-requisite for efficient parallelization via domain decomposition. In short, we develop an efficient Maxwell solver for full three-dimensional problems where one axis is distinguished from the two others (e.g. the laser- or beam-propagation axis).

We select the $X-$direction for dispersionless propagation. For electromagnetic waves propagating in $X$, we have the Maxwell equations

[TABLE]

[TABLE]

Here, the vector

[TABLE]

combines the vacuum diffraction $\mathbf{E}$ and the medium response (currents) $\mathbf{J}$, while

[TABLE]

is the vacuum diffraction operator for $\mathbf{B}$.

We use a semi-implicit trapezoidal (sometimes called “implicit midpoint”) scheme for the discretization of the transverse fields on a 3D grid. We write here explicitly the $i-$index along the $X-$axis only as the scheme can be easily generalized for arbitrary transverse geometries (e.g. Cartesian, or cylindrical, etc.):

[TABLE]

[TABLE]

Here, $\tau$ is the time step and $h_{x}$ is the spatial grid step in the $X-$direction.

These equations (17)-(22) build a system of coupled linear equations relating the updated fields at the time step $n+1$ with already known fields at the time steps $n$ and $n+1/2$. Although this implicit system of linear equations can generally be solved using a fast matrix inversion method (the system has a sparse matrix), we will be interested in the special case $c\tau=h_{x}=\Delta$. In this particular case, the inversion is straightforward.

First, we add Eqs. (17)+(21) and (18)+(20). to obtain transport components

[TABLE]

or simply

[TABLE]

Then, we substract the same Eqs. (17)-(21) and (18)-(20) to obtain

[TABLE]

or

[TABLE]

These are the marching equations. The transport components $T_{y,z}^{+/-}$ must be shifted one cell in the corresponding direction and the diffraction/refraction terms be correctly added.

This marching has a form of “Galiliean field shift” exactly by single cell per time step. Thus, instead of shifting the grid following the relativistic plasma [44], the RIP solver shifts the transverse fields so that the relativistic particle sees the same fields when it enters the new cell. In one-dimensional geometry, the new algorithm defaults to the well known advective algorithm introduced by Birdsall and Langdon [14]. In [24], and as reported also in [14], on the stability of various electromagnetic PIC schemes, it is stated that ”the improved stability associated with the advective differencing schemes is due not so much to the dispersionless vacuum transport of the fields, per se, as to the less conventional methods of determining the mesh current usually employed with advective differencing”.

For the fields, we get

[TABLE]

[TABLE]

or simply

[TABLE]

[TABLE]

4 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 (3)-(3). The pair $\mathbf{E_{y}},\mathbf{B_{z}}$ is located at the rhombi vertices $(i,j+1/2,k)$. The pair $\mathbf{E_{z}},\mathbf{B_{y}}$ is located at the rhombi vertices $(i,j,k+1/2)$. The longitudinal field $\mathbf{E_{x}}$ we place at point $(i,j,k)$ which is the center of the full integer rhombus. The longitudinal field $\mathbf{B_{x}}$ we place at the center of the half integer rhombus $(i,j+1/2,k+1/2)$. The grid is shown in Fig.1.

Then, the diffraction/refraction terms at the half time step will be:

[TABLE]

[TABLE]

Similar formulas are obtained for the fields at the half-time steps, where all fields and currents are shifted by a half time step.

The use of the transport vectors $\mathbf{T}_{\perp}$ makes the boundary conditions in the $x-$direction trivial. One sets the inbound T vectors equal to the incident laser pulse and outbound T vectors to zero at the boundaries. This procedure absorbs waves normally incident on the boundaries exactly.

It seems that we have to maintain two sets of fields for each time step: fields at the full step and at the half step. The particles however, can be pushed just once per time step. For the particle push we use the symplectic semi-implicit mid-point scheme of Higuera and Hary [45] (pushers of Boris [46] and Vay [47] produce hardly discernible results) at the full time step:

[TABLE]

where $\mathbf{p}=(\mathbf{p}_{\alpha}^{n+1/2}+\mathbf{p}_{\alpha}^{n-1/2})/2$ and $\gamma=\sqrt{1+p^{2}/m^{2}c^{2}}$.

These momenta are used to generate currents $\mathbf{j}^{n+1/2}$ at the half time steps. Currents at the full time step required to push the half-time step fields can be obtained by simple averaging on the grid

[TABLE]

To ensure the Lorentz invariance and charge conservation of the scheme, the current components are defined within the cell at the same positions as the corresponding $\mathbf{E}-$field components.

5 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.

5.1 Generalized rigorous charge conservation

The numerical continuity equation on the RIP grid has the form

[TABLE]

The charge conservation on the grid can be enforced in various ways. One can correct currents and solve an elliptic problem after each time step [14]. One can atomatically calculate currents inside the cell co that the charge is locally conserved [48]. One can move particles in a zig-zag along the axises [49]. The charge-conserving closure is not unique and an infinite number of other schemes can be easily generated like charge splitting curvy trajectory particle motion, etc. All these schemes generate different curl currents on the grid and thus have different noise properties.

However, the only true 2-nd order accurate current closure is the rigorous charge conservation method introduced originally by Villacenor and Buneman [50]. This scheme assumes the straight particle trajectory during the time step. All other methods fail to do so. The Esirkepov scheme [48] coincides with the Villacenor and Buneman scheme *identically * as long as the particle stays inside one cell during the time step. It gives different results, however, as soon as the particle crosses boundaries.

We use a generalized rigorous charge conservation (GRCC) method based on [50], compare also [51]. It is not limited to the Cartesian geometry and is valid for any particle shape. Let us suppose, we have selected a form-factor $\mathbf{w}$ for the current deposition by the numerical macro-particles. The macro-particle $\alpha$ will then induce an instantaneous current $\mathbf{j(t)}$ on the grid with components:

[TABLE]

where $\mathbf{r}^{\alpha}(t)$ is the instantaneous particle position inside the cell. Depending on the form-factor, the particle may induce instantaneous currents at many grid cells. We write expressions for one cell only, as the others are analogous.

The particle starts its motion at the time step $t^{n}$ at the position $\mathbf{r}^{\alpha}(t^{n})$ and finishes at the time step $t^{n+1}=t^{n}+\Delta/c$ at the position $\mathbf{r}^{\alpha}(t^{n+1})$. The straight particle trajectory is parameterized as $\mathbf{r}^{\alpha}(t)=\mathbf{r}^{\alpha}(t^{n})+(t-t^{n})/(t^{n+1}-t^{n})\mathbf{\delta r}^{\alpha}$, where $\mathbf{\delta r}^{\alpha}=\mathbf{r}^{\alpha}(t^{n+1})-\mathbf{r}^{\alpha}(t^{n})$. The current induced by the particle on the grid during the time step is then

[TABLE]

When the macroparticle crosses the cell boundaries, the integral (49) must be split along the straight particle trajectory in parts, where the particle center belongs to one particular cell.

For most popular particle shapes (box, triangle, quadratic, spline, etc.), the integration in (49) is done analytically (using any symbolic integration software) and the function $\mathbf{W^{\alpha}}(\mathbf{r}^{\alpha}(t^{n}),\mathbf{\delta r}^{\alpha})$ is easily coded. Villacenor and Buneman did it explicitly for the case of Cartesian grid and a square box particle shape.

This described GRCC algorithm preserves the discretized Maxwell-Gauss equations automatically by the RIP scheme.

5.2 Divergence conservation

In the same manner we define the curl and divergence of the fields. For example, $\nabla\cdot\mathbf{B}$ is defined as

[TABLE]

at the middle cell position $\left(i+1/2,j+1/2,k+1/2\right)$. We average fields along the $X-$axis according to the trapezoidal rule while the usual leap-frog Yee rule is applied along the transverse coordinates. It is straightforward to check that the RIP scheme preserves $\nabla\cdot\mathbf{B}$ defined in this way.

The fields at the half-time steps are required to calculate the diffraction terms only. Without diffraction, the need to maintain the additional set of fields at half-time steps vanishes and the RIP scheme becomes identical to the standard 1D PIC scheme [14], which is the workhorse of 1D plasma simulations due to its excellent stability and accuracy.

6 Dispersion and stability of the RIP scheme

We apply the plane-wave analysis to the marching equations (19), (22) and (3)-(3) with the refraction/diffraction terms (4)-(4) assuming $\mathbf{F}=\mathbf{\tilde{F}}\exp\left(-i\omega t+i\mathbf{k}\mathbf{r}\right)$. For simplicity, we assume uniform plasma frequency $\omega_{p}^{2}=4\pi ne^{2}/\gamma$ and the linear current response to the electric field $\frac{2c}{\Delta}\sin\frac{\omega\tau}{2}\tilde{\mathbf{J}}=iq^{2}n\tilde{\mathbf{E}}$. For the case of interest, $c\tau=h_{x}=\Delta$, these equations become

[TABLE]

[TABLE]

The dispersion relation in vacuum ($\omega_{p}=0)$ is rather simple:

[TABLE]

The stability condition in vacuum is

[TABLE]

In the presence of plasmas, it is modified to

[TABLE]

The RIP scheme combines dispersionless properties of the standard 1D solver along the $X-$axis with the Yee dispersion for waves running in the transverse direction. Indeed, setting $k_{y}=k_{z}=0$ in the dispersion relation (57), we immediately obtain $\omega=ck_{x}$ and the phase velocity

[TABLE]

for plane waves propagating in the $X-$direction.

Conversely, setting $k_{x}=0$, we obtain the usual 2D Yee dispersion relation for waves propagating in the transverse direction

[TABLE]

with all its known advantages and drawbacks.

Mention that the NDF scheme introduced in [52] has a different stability condition:

[TABLE]

so that one can set $c\tau=h_{x}$ only in vacuum. The presence of plasma, $\omega_{p}>0$, one has to choose $c\tau<h_{x}$ and the dispersionless properties of the NDF scheme are compromised.

7 Numerical tests of the RIP Maxwel solver

7.1 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=n_{0}\exp\left(-r^{2}/\sigma^{2}\right)$ consisting of electrons and protons moving in the $X-$direction with the average momentum $<\mathbf{p_{0}}>/m_{\alpha}c=\left(p_{x0},0,0\right)$, where $\alpha$ denotes the particle type ($\alpha=e,p$), with $m_{p}/m_{e}=1846$. To seed the instability, the electrons have a small initial temperature $<\left(\mathbf{p}_{0}-<\mathbf{p}_{0}>\right)^{2}>=\sigma_{p}^{2}$. In relativistically normalized units, the simulation parameters are: the peak plasma density is $n_{0}=1$ with the corresponding non-relativistic plasma frequency $\omega_{p}=\sqrt{4\pi n_{0}e^{2}/m_{e}}$. The initial particle momenta $p_{x0}=-10$ and $\sigma_{p}=10^{-4}$. The grid steps were $h_{y}=h_{z}=1.88\,c/\omega_{p}$ and $h_{x}=c\tau=\Delta=0.63\,c/\omega_{p}$. As a diagnostics for the comparison, we selected the growth of the maximum local field intensity $I=\mathbf{E^{2}}+\mathbf{B}^{2}$ on the grid. The results are shown in Fig.2

We see that the fluctuating fields in simulations using the Yee scheme grow to the non-linear saturation within a few plasma oscillations. This is because the Yee solver is exposed to the first order Cerenkov resonance. The FFT-based solver is dispersionless and avoids the first order Cerenkov resonance. Still, the second order aliasing of the spectral FFT solver leads to the numerical Cerenkov instabilty, though at a lower growth rate as compared with the standard Yee solver.

In contrast, the RIP solver is free from NCI. The noise in the RIP scheme remains many orders of magnitude lower over a long simulation time of $t=100\cdot 2\pi/\omega_{p}$. The very slow growth of the noise fields here has nothing to do with the Cerenkov resonance, but is the unavoidable “numerical heating” always present in PIC codes.

Finally, we do another simulation with a stationary plasma, $<\mathbf{p_{0}}>/m_{\alpha}c=\left(0,0,0\right)$, while keeping all other parameters the same. We observe here no numerical heating at all. Intensity of fluctuating noise fields remains constant over many hundreds of plasma periods here. The higher absolute level of the noise for the streaming plasma is the natural consequence of the larger initial noise current source in this case. Fig.2 demonstrates clearly that the RIP scheme is much less subject to Cerenkov instability for plasmas drifting along the selected axis.

We stress here that the FFT-based method (2) is identical to the RIP solver for waves running in the $X-$direction. Yet, we observe a quite different behaviour with respect to the numerical Cerenkov instability. The reason for this difference has to be studied further.

We mention here that the numerical Cerenkov instability of uniformly streaming plasma can be alleviated by using a co-moving grid as proposed by Lehe et al. [44]. The method exploits a Galilean transformation to a grid in which the background plasma does not stream through the cell boundaries. Yet, the dense bunch of accelerated particles moves in the opposite direction at twice the light speed relative to this grid and is fully exposed to the Cerenkov resonance.

7.2 Laser-driven plasma bubble

As the second numerical test, we select laser-plasma particle acceleration in the bubble regime [53]. A circularly polarized laser pulse with initial vector potential $\mathbf{A}=\Re\left[a(\xi,\mathbf{r_{\perp}})(\mathbf{e_{y}}+i\mathbf{e_{z}})\exp(ik\xi)\right]$ is used. Here, $\xi=x-ct$ and the envelope shape has been selected as a ellipsoidal Gaussian $a(\xi,\mathbf{r_{\perp}})=a_{0}\exp(-\xi^{2}/\sigma_{||}^{2}-r_{\perp}^{2}/\sigma_{\perp}^{2})$ with the amplitude $a_{0}=5$ , length $\sigma_{||}=5\lambda$ and radius $\sigma_{\perp}=5\lambda$, where the laser wavelength $\lambda=2\pi/k$. The plasma consisting of electrons and protons has an initial density $n=0.01\,n_{c}$, where $n_{c}=m_{e}\omega^{2}/4\pi e^{2}$ is the critical density. At the plasma boundary, the density increases lineraly from $n=0$ to $n=0.01\,n_{c}$ over a length $L=38\,\lambda$. The simulation results after an acceleration distance of $L_{a}=300\,\lambda$ are shown in Fig.3. The simulation box has the size $40\lambda\times 40\lambda\times 40\lambda$. The grid steps are $h_{x}=0.05\lambda$, $h_{y}=0.25\lambda$, $h_{z}=0.25\lambda$ and the time step is $\tau=0.045\lambda/c$ in the Yee simulation, and $\tau=h_{x}/c$ in the RIP simulation.

We see that the trapped electron bunch of the bubble has a fine longitudinal structure in the Yee simulation. At the same time, the bubble accelerating field $E_{x}$ is rippled with the short-wavelength radiation emitted by the relativistic electrons due to the numerical Cerenkov resonance. This numerical emission is clearly seen in Fig. 3(c) as the bow-like short wavelength radiation emanating from the dense electron bunch. The RIP simulation shows a rather smooth electron bunch and no signatures of numerical Cerenkov emission. The $E_{y}-$field of the relativistic electron bunch has a clean quasi-static form: it is not bow-shaped, but perpendicular to the bunch. Further, a small additional numerical dephasing can be observed at the leading edge of the bubble.

To check the RIP scheme convergence, we did an additional simulation with rough resolution. We doubled the longitudinal grid step and the time step to $h_{x}=c\tau=0.1\lambda$, so that we have only 10 cells per laser wavelength. The results are shown in the last row in Fig. 3. One observes little difference from the higher resolution simulation, shown in the middle row in Fig. 3. Compare the phase of the laser pulse seen in the electron density perturbations in frames (f) and (i). The RIP simulation even at this rough resolution accurately describes the laser phase. The Yee solver gives a completely wrong laser phase as seen in the frame (b). At the resolution of just 10 cells per wavelength (or even smaller), the RIP solver is a good alternative to full electromagnetic PIC codes that employ the envelope approximation [23, 54] when the laser pulse is short.

7.3 Ultra-low transverse emittance beam propagation in vacuum

Finally, we check how well the RIP scheme preserves the beam emittance. The future colliders and XFEL light sources must have beams with ultra-low transverse emittance. Emittances in the range of a few nanometers, or even picometers have to be realized. There are several approaches, how such beams can be generated using the conventional accelerators. Plasma-based acceleration also might reach such ultra-low beam emittances, using, e.g. the Trojan-horse injection [55]. Thus, very accurate simulation methods are required, where the emittance is preserved at picometer levels.

Lehe et al [29] have shown that the Yee Maxwell solver has problems with emittance conservation and suggested an ”improved” one that shows a better emittance preservation. There are several sources of emittance growth in PIC simulations: numerical heating, numerical Cerenkov instability and the wrong field interpolation to the particle position due the staggered mesh. Lehe et al. modified the Yee solver to remove the zero-order Cerenkov resonance. This improved the emittance preservation [29].

In this sub-section, we simulate an electron bunch propagation in vacuum. The electron bunch has initial energy of 10 GeV and normalized transverse emittance of $\varepsilon_{y}=\varepsilon_{z}=1~{}$nm. The bunch has a Gaussian shape with $\sigma_{x}=\sigma_{y}=\sigma_{z}=1.41~{}\mu$m. We simulate two cases, where the bunch carries either $0.6~{}$nC or $6.0~{}$nC charge. The simulations are done using the Yee or RIP solver for total propagation distance of $1~{}$mm. The grid steps were $h_{x}=0.1~{}\mu$m, $h_{y}=h_{z}=0.15~{}\mu$m. The time step for the Yee solver was $c\tau=0.7h_{x}$.

The emittance evolution is shown in Fig. 4. The Yee solver shows a very fast inital jump of the bunch emittance due to the incorrect field interpolation to the particle positions on the staggered Yee mesh. The jump is higher for the higher bunch charge. Later, the bunch becomes unstable due to the numerical Cerenkov instability and the emittance grows steadily. After the full propagation distance of $1~{}$mm, the final numerical emittance grew in the Yee simulations to $8.5~{}\mu$m for the high current bunch of $6~{}$nC and to $0.8~{}\mu$m for the low current bunch of $0.6~{}$nC.

The RIP solver shows an excellent conservation of the emittance. For the high-current bunch of $6~{}$nC, the transverse emittance grew only by $10~{}$pm, from $1~{}nm$ to $1.01~{}$nm after the full $1~{}mm$ propagation distance. This miniscule emittance growth is physical: it is due to the Coulomb explosion of the high current bunch.

The bunch dynamics is shown in Fig. 5. It shows the normalized bunch density $n_{e}/n_{c}$ and the tranverse force acting on the electrons $F_{y}=e(E_{y}-B_{z})/mc\omega$. Here, $\omega=2\pi c/\lambda$, $\lambda=1~{}\mu$m and $n_{c}=m\omega^{2}/4\pi e^{2}$. The frames (a) and (b) are taken for the RIP solver after the full $1~{}$mm propagation distance, while the frames (c) and (d) are taken for the Yee solver after $0.1~{}$mm propagation. We take here the high current case of $6~{}$nC charge bunch.

The Yee solver fails with the transverse force by many orders of magnitude. We also see the development of NCI in the transverse force causing self-modulation of the bunch tail. Apparently, the Yee solver is not the best choice when one wants to simulate low emittance bunches.

The RIP solver accurately reproduces the transverse force due to the self-interaction down to the machine precision. The normalized transverse force is at the level of $F_{y}=3.7\cdot 10^{-7}$ in this case. Thus, the RIP Maxwell solver is perfectly suited to simulate bunches with sub-nanometer emittances.

8 Discussion

The new RIP scheme is a compact stencil FDTD Maxwell solver that removes the numerical dipersion 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 costs of the RIP solver is comparable with that of the standard Yee solver. The RIP solver can be used for simulations of quasi-1D physics problems like laser wake field acceleration. This RIP marching algorithm has a form of “Galiliean field shift” exactly by single cell per time step. Thus, instead of shifting the grid following the relativistic plasma [44], the RIP solver shifts the transverse fields so that the relativistic particle sees the same fields when it enters the new cell. Apparently, this procedure greatly reduces the numerical Cerenkov instability.

Acknowledgements

This work has been supported in parts by BMBF and DFG (Germany).