An \mathcal{O}(N) Maxwell solver with improved numerical dispersion properties

TL;DR

This paper introduces an (N) Maxwell solver based on finite element methods that enhances numerical dispersion properties in relativistic PIC simulations, suitable for parallel computing and reducing instabilities.

Contribution

It presents a novel finite element Maxwell solver with (N) complexity that improves dispersion accuracy and suppresses numerical instabilities in high-speed plasma simulations.

Findings

Reduces Numerical Cherenkov instability growth rate

Achieves (N) computational complexity

Enhances dispersion properties in relativistic PIC simulations

Abstract

A Maxwell solver derived from finite element method with \mathcal{O}(N) computing cost is developed to improve the numerical dispersion properties in relativistic particle-in-cell (PIC) simulations. The correction of the dispersion relation of the electromagnetic wave is achieved using the neighboring cells via an iteration scheme without decomposing into Fourier modes. The local nature of the communication is ideally suited to massively parallel computer architectures. This Maxwell solver constrains the Numerical Cherenkov instability (NCI) for the ultra-relativistic drifting pair plasma in x direction to large wave vectors for two dimensional grid. The growth rate of NCI is suppressed by using the low pass filtering.

Tables

Table 1. Table 1: Parameters for test problem of a drifting pair plasma. The definition of electron skin depth is d e = c ω p = c / 4 ​ π ​ q 2 ​ n e γ ​ m e subscript 𝑑 𝑒 𝑐 subscript 𝜔 𝑝 𝑐 4 𝜋 superscript 𝑞 2 subscript 𝑛 𝑒 𝛾 subscript 𝑚 𝑒 d_{e}=\frac{c}{\omega_{p}}=c/\sqrt{\frac{4\pi q^{2}n_{e}}{\gamma m_{e}}} .

domain size	$L_x=16d_e$, $L_y=8d_e$
boundary condition	periodic, both in $x$ and $y$
number of cells	$N_x=256$, $N_y=128$
pseudo-particles per cell	$N_{\mathrm{PPC}}=512$ (256 for each species)
drift Lorentz factor	$\gamma =1000$
temperature	$k_B{T}_e=k_B{T}_i=0.01m_e{c}^2$
time step	$\Delta t/\Delta t_{\mathrm{CFL}}=0.1,0.2,\dots 1.0$
particle shape	1th, 2nd, 4th order B-spline

Equations

$$\frac{1}{c}{D}_t{B}_x$$

$$\frac{1}{c}{D}_t{B}_y$$

$$\frac{1}{c}{D}_t{B}_z$$

$$\frac{1}{c}{D}_t({\alpha }_{yz}{E}_x)$$

$$\frac{1}{c}{D}_t({\alpha }_{zx}{E}_y)$$

$$\frac{1}{c}{D}_t({\alpha }_{xy}{E}_z)$$

$$D_x({\alpha }_{yz}{E}_x)+D_y({\alpha }_{zx}{E}_y)+D_z({\alpha }_{xy}{E}_z)$$

$$D_x{B}_x+D_y{B}_y+D_z{B}_z$$

$$D_x{f}_{k,l,m,n}$$

$${\alpha }_x{f}_{k,l,m,n}$$

$$E_x^{k+1/2,l,m,n},E_y^{k,l+1/2,m,n},E_z^{k,l,m+1/2,n},$$

$$B_x^{k,l+1/2,m+1/2,n+1/2},B_y^{k+1/2,l,m+1/2,n+1/2},$$

$$B_z^{k+1/2,l+1/2,m,n+1/2},$$

$$j_x^{k+1/2,l,m,n},j_y^{k,l+1/2,m,n},j_z^{k,l,m+1/2,n}$$

$${\rho }^{k,l,m,n+1/2}$$

$${\alpha }_x$$

$$D_x$$

$$D_t$$

$$[\overrightarrow{k}]\times \overrightarrow{E}$$

$$[\overrightarrow{k}]\times ({\mathcal{A}}_B\overrightarrow{B}$$

$$[k{]}_i$$

$$[\omega ]$$

$${\mathcal{A}}_B$$

$${\mathcal{A}}_E$$

$$[\omega {]}^2{\mathcal{A}}_E\overrightarrow{E}+c^2[\overrightarrow{k}]\times {\mathcal{A}}_B([\overrightarrow{k}]\times \overrightarrow{E})=-4\pi i[\omega ]\overrightarrow{ȷ}$$

$${ϵ}_{xx}^{(\mathrm{vac})}$$

$${ϵ}_{yy}^{(\mathrm{vac})}$$

$${ϵ}_{zz}^{(\mathrm{vac})}$$

$${ϵ}_{xy}^{(\mathrm{vac})}$$

$${ϵ}_{yz}^{(\mathrm{vac})}$$

$${ϵ}_{zx}^{(\mathrm{vac})}$$

$$\det {\overleftrightarrow{ϵ}}^{(\mathrm{vac})}$$

$$\displaystyle\qquad-\alpha_{x}^{-1}[k]_{x}^{2}+\alpha_{y}^{-1}[k]_{y}^{2}+\alpha_{z}^{-1}[k]_{z}^{2}\big{)}^{2}$$

$${\alpha }_t^{-1}\frac{[\omega {]}^2}{c^2}={\alpha }_x^{-1}[k{]}_x^2+{\alpha }_y^{-1}[k{]}_y^2+{\alpha }_z^{-1}[k{]}_z^2$$

$$G(\theta )=\frac{{\sin }^2\theta }{1-\frac{2}{3}{\sin }^2\theta }$$

$$\frac{G(\frac{\omega \Delta t}{2})}{c^2\Delta t^2}=\frac{G(\frac{k_x\Delta x}{2})}{\Delta x^2}+\frac{G(\frac{k_y\Delta y}{2})}{\Delta y^2}+\frac{G(\frac{k_z\Delta z}{2})}{\Delta z^2}$$

$$c\Delta t_{3D}\le c\Delta t_{3D,\mathrm{CFL}}=(1/\Delta x^2+1/\Delta y^2+1/\Delta z^2{)}^{-1/2}$$

$$\omega /k$$

$$\displaystyle\phantom{=c\bigg{[}1}+\mathcal{O}(k^{4})\bigg{]}$$

$$\omega /k$$

$$\displaystyle\phantom{=c\bigg{[}1}+\mathcal{O}(k^{4})\bigg{]}$$

$$\frac{G(\frac{\omega \Delta t}{2})}{c^2\Delta t^2}=\frac{G(\frac{k_x\Delta x}{2})}{\Delta x^2}+\frac{G(\frac{k_y\Delta y}{2})}{\Delta y^2}$$

Full text

1 Introduction

Particle-in-cell (PIC) method[2] is widely used for the simulations of plasma dynamics ranging from Laser Plasma Accelerators (LPAs) to collisionless astrophysical problems. The traditional Yee scheme[15] introduces a resonance between the beam and electromagnetic (EM) wave, which results in numerical Cherenkov instability (NCI) [5]. The resonance is a numerical artifact, which is due to two fundamental properties in the numerical schemes: (1) In the Maxwell solver, the dispersion relation of the vacuum electromagnetic waves deviates from the physical one, i.e. $\omega=ck$ (we call $\omega=ck$ luminous, $\omega>ck$ superluminous, $\omega<ck$ subluminous). (2) The dispersion of the drifting plasma has contribution of aliases[9]. For example, in Yee solver the dispersion is superluminous along any axis but luminous along diagonals and the main beam will intersect at some direction between these two directions for low $k$. This numerical artifact has significant impact on the quality and physical interpretation of relativistic PIC simulations. The numerical instabilities near $k=0$, once generated, are difficult to be differentiated from the physics modes, which are also near $k=0$. A numerical scheme without numerical resonance near $k=0$ or suppressing the NCI growth rate resolves the numerical artifact. Several approaches have been proposed to improve the numerical dispersion properties[14, 3, 7, 8, 6, 13, 16, 10], including the optimal choice of time step and spectral solvers.

One approach to minimize the NCI growth rate is achieved by optimal choice of time step, which was found for both energy conserving and momentum conserving scheme by using the analytical formula for the asymptotic growth rate[14]. The asymptotic NCI growth rate is useful for optimizing the numerical dispersion properties of different interpolation schemes and Maxwell solvers.

Another approach to improve the numerical dispersion is to use the spatial Fourier transform based methods, such as Pseudo-Spectral Time Domain (PSTD) algorithms[7, 8, 6]. Those methods can constrain numerical resonance between EM mode and beams to large wave vectors, i.e. no resonance near $k=0$. The numerical instabilities can then be suppressed by applying low-pass filtering[13]. The non-local nature of communication cost and the $\mathcal{O}(N\log N)$ computational cost of the spatial Fourier transform algorithms make those methods challenging to scale on massively parallel computers. In the generalized PSTD algorithm[6], the components of the wavevectors can be replaced by the Fourier transforms of finite difference approximations to spatial derivatives on a grid, which reduces the communication cost to local and computational cost to $\mathcal{O}(N)$. The mixed FD-FFT solver use 1D FFT only[16] or high order FDTD in one direction only[10]. However, the mixed solver nature in different directions requires correction due to the loss of charge conservation.

In this paper, we aim to develop an alternative numerical scheme on multi-dimensional grid with superluminous numerical dispersion in a large range of $k$ including $k=0$, and with $\mathcal{O}(N)$ computing cost and only local communication. The discretized set of Maxwell equations without lumping, i.e. keeping the averaging operators, from Eastwood 1991[4] is found to satisfy all our requirements and has been implemented in this work using the EPOCH code[1]. If the charge conserving current deposition is used, the set of equations also conserve the Gauss’s law for electric field, there is no need for divergence cleaning of electric field or correction of the current. The divergence-free nature of the magnetic field is also preserved. In principle, the numerical scheme we develop here can be implemented in any Yee-grid based PIC code and generalized to complex geometries. Our numerical scheme modifies Ampere’s equations, instead of modifying the Faraday equation in the finite-difference time-domain (FDTD) method[3].

This paper is organized as follows. The algorithm for solving Maxwell equations and the dispersion properties of the numerical electromagnetic waves in vacuum are described in Sec 2. The NCI growth rate for drifting pair plasma is discussed in Sec 3, following the analytical technique in Xu 2013[14]. The numerical experiments are carried out to confirm the two dimensional EM dispersion relation and the NCI growth rate in Sec 4.

2 Maxwell Solver

In Eastwood [4], the coupled relativistic Vlasov-Maxwell set of equations are derived using finite elements in both space and time. Their finite element derivations are suitable for complex geometries. In this paper, we solve the Vlasov-Maxwell equation on uniform rectangular grid by using Faraday and Ampere’s equations from Eastwood, without modifying other modules from the usual finite difference PIC, i.e., force interpolation, current deposition and particle pusher. The discretized Faraday equations from Eastwood 1991 can be generalized to 3D

[TABLE]

and the Ampere’s equations are

[TABLE]

and the Gauss’s law for electric field and magnetic field

[TABLE]

where difference operators $D_{x}$, $D_{y}$, $D_{z}$, $D_{t}$ and averaging operators $\alpha_{x}$, $\alpha_{y}$, $\alpha_{z}$, $\alpha_{t}$ respectively act on the spatial indices $k$, $l$, $m$ and time step index $n$:

[TABLE]

and similarly for $D_{y}$ and $\alpha_{y}$ on index $l$, and for $D_{z}$ and $\alpha_{z}$ on index $m$, and for $D_{t}$ and $\alpha_{t}$ on the index $n$. The product of two averaging operators $\alpha_{y}$ and $\alpha_{z}$ is abbreviated as $\alpha_{yz}$, similarly for averaging operators in other directions. All the indices can be integers or half integers, depending on the staggered directions of the grid and time step. As derived from finite element method[4], the averaging stencil is $(\frac{1}{6},\frac{2}{3},\frac{1}{6})$ for spatial grid in all three spatial dimensions and also for time step. The variables for electric field $\overrightarrow{E}$, magnetic field $\overrightarrow{B}$ and current density $\overrightarrow{\jmath}$ are defined either at grid point or has a half-grid offset, on integer or half-integer time step, as same as in the Yee scheme. The whole set of discrete variables is

[TABLE]

where $k$, $l$, $m$, $n$ are integers. By lumping the averaging operators, $\alpha_{x},\alpha_{y},\alpha_{z},\alpha_{t}\to 1$ in Eqs(1) to Eq(6), we recover the discretized equations for Yee solver. In the numerical scheme we use in this work, we keep those averaging operators instead of lumping them in order to achieve the desired numerical dispersion properties.

By performing Fourier transform on Eqs(1) to (6), assuming that all the variables have $e^{i(k_{x}x+k_{y}y+k_{z}z-\omega t)}$ form, the difference and averaging operators can be written in terms of frequency $\omega$ or wave vector $\overrightarrow{k}$

[TABLE]

and similarly for $\alpha_{t},\alpha_{y},\alpha_{z},D_{y},D_{z}$, where $\Delta x$ is the grid size in $x$ direction, $\Delta t$ is the time step. The Eqs(1) to (6) in Fourier space can be written as

[TABLE]

where

[TABLE]

[TABLE]

By multiplying $[\omega]$ on Eq(15) and using Eq(14) to eliminate $\overrightarrow{B}$, we get

[TABLE]

We define the vacuum part of dielectric tensor by $\overleftrightarrow{\epsilon}^{(\mathrm{vac})}\cdot\overrightarrow{E}=[\omega]^{2}\mathscr{A}_{E}\overrightarrow{E}+c^{2}[\overrightarrow{k}]\times\mathscr{A}_{B}([\overrightarrow{k}]\times\overrightarrow{E})$, then the matrix elements of $\overleftrightarrow{\epsilon}^{(\mathrm{vac})}$ can be written down

[TABLE]

By lumping the averaging operators, i.e. letting $\mathscr{A}_{E}\to\mathbb{1}$ and $\mathscr{A}_{B}\to\mathbb{1}$ in Eq(20), the vacuum dielectric tensor for Yee solver is recovered.

It is derived in Appendix A that Eq(7) is conserved automatically if it is fulfilled at the initial condition and a charge conserving current deposition scheme is used. And Eq(8) is conserved if it is fulfilled at the initial condition.

2.1 Dispersion relation for EM waves in vacuum

In vacuum we have zero current $\overrightarrow{j}=0$, the dispersion relation of EM waves can be obtained by calculating the determinant of the dielectric tensor $\overleftrightarrow{\epsilon}^{(\mathrm{vac})}$

[TABLE]

The $\omega=0$ mode represents the background field while the other root of $\omega$ represents the EM mode. Thus we obtain the vacuum EM wave dispersion relation

[TABLE]

We define the function $G(\theta)$ to be

[TABLE]

Then Eq(28) becomes

[TABLE]

For numerical stability, the frequency $\omega$ must be real which, by Eq(30) implies the same time step constraint as Yee solver

[TABLE]

The detail derivation of Eq(31) is in Appendix B. For $k\to 0$, we expand the phase speed $\omega/k$ to second order in $k$

[TABLE]

From Eq(32) we derive in Appendix B that the EM mode for sufficiently small $k$ always has phase speed larger than or equal to $c$, i.e. the speed of light, so the resonance between EM mode and unaliased plasma beam with any physical speed $v<c$ cannot be located near $k=0$.

For a two-dimensional grid, we can write down the dispersion relation of the EM mode, the CFL limit of time step and the expansion of phase speed in the similar way

[TABLE]

[TABLE]

[TABLE]

For the the $\Delta x=\Delta y$ case, we plot the phase velocity in $k$ space for different values of time step in Fig. 1.

In the 2D or 3D case for $\Delta t=\Delta t_{\mathrm{CFL}}$, the phase speed is constant and equals to $c$ along the diagonal of the first Brillouin zone, i.e. $|k_{x}|\Delta x=|k_{y}|\Delta y=|k_{z}|\Delta z$ in 3D case and $|k_{x}|\Delta x=|k_{y}|\Delta y$ in 2D case. For $\Delta t\neq\Delta t_{\mathrm{CFL}}$ or off-diagonal direction, the second order term in Eq(32) and (33) is always larger than zero, thus $\omega/k$ is larger than $c$ for sufficiently small $k$.

2.2 Numerical iteration scheme

2.2.1 Equations for three dimensional case

Updating $\overrightarrow{B}$ is fully explicit by using Eqs(1) to (3). However, Eqs(4) to (6) contain the time averaging operator on the L.H.S., which require the field at a previous step and a future step. This requirement can be eliminated by substituting Eq(1) to (3) into the averaging operators. For example, the time averaging of $B_{z}$ can be written as

[TABLE]

where we used Eq(3) in the derivation. Similarly we can write down the time averaging of $B_{y}$ as

[TABLE]

Substituting Eq(36) and (37) into the time averaging terms in Eq(4), we obtain

[TABLE]

where $n$ is the index for $n$th time step. If we define the auxiliary variable $E_{x}^{\prime}$ such that

[TABLE]

then we can simplify Eq(38) to be

[TABLE]

The equations for the auxiliary variables are similar in $y$ and $z$ directions

[TABLE]

And the change of electric field vector from $n$th step to $(n+1)$th step $\Delta\overrightarrow{E}=(\Delta E_{x},\Delta E_{y},\Delta E_{z})$ and the change of the auxiliary vector $\Delta\overrightarrow{E}^{\prime}=(\Delta E_{x}^{\prime},\Delta E_{y}^{\prime},\Delta E_{z}^{\prime})$ satisfy the linear transform

[TABLE]

The matrix elements of $A$ can then be written down

[TABLE]

[TABLE]

where we define the dimensionless operator $\bar{D}_{x}$ to be the value difference between the right cell and left cell in $x$ direction

[TABLE]

and similarly for $y$ and $z$ directions. $\bar{D}$ with multiple indexes represents the abbreviation of multiplication, e.g. $\bar{D}_{xx}=\bar{D}_{x}\bar{D}_{x}$ and $\bar{D}_{xxyy}=\bar{D}_{x}\bar{D}_{x}\bar{D}_{y}\bar{D}_{y}$.

To update the electric field, we first compute the change of auxiliary field explicitly by using Eqs(40) to (42). Then the change of the electric field can be obtained by solving Eq(43) and the electric field is updated simply by the increment $\overrightarrow{E}^{n+1}=\overrightarrow{E}^{n}+\Delta\overrightarrow{E}$. To solve for $\Delta\overrightarrow{E}$ by Eq(43), we propose to use the iteration method listed in Alg 1. We prove in Appendix C that the infinity norm of matrix $\overleftrightarrow{A}$ is always less unity in 1D and 2D case. For 3D case, the infinity norm of $A$ is less than unity if $\Delta t\leq\frac{\min(\Delta x,\Delta y,\Delta z)}{\sqrt{2}c}$. In the 3D and $\Delta x=\Delta y=\Delta z$ case the infinity norm of matrix $A$ is always less than $1$ because $\Delta t\leq\frac{\Delta x}{\sqrt{3}c}$. Thus the iteration in Alg 1 always converges to the exact solution. Further more the infinity norm or eigenvalues of matrix $\overleftrightarrow{A}$ does not depend on the number of cells on the mesh, which implies that the speed of convergence is not dependent on the size of the problem.

For several reasons, we propose to use a predetermined number for the number of iterations instead of determine $|\Delta\overrightarrow{E}_{(m+1)}-\Delta\overrightarrow{E}_{(m)}|<\varepsilon$ during the iterations, (1) calculating the norm $|\Delta\overrightarrow{E}_{(m+1)}-\Delta\overrightarrow{E}_{(m)}|$ and determine $|\Delta\overrightarrow{E}_{(m+1)}-\Delta\overrightarrow{E}_{(m)}|<\varepsilon$ requires a global reduction and broadcasting among processors (2) keeping fixed number of iterations can keep the consistency of the error of dispersion relation throughout all time steps (3) fixed number of iterations is sufficient for reducing the error, e.g. Fig. 2(a) shows that for $\Delta x=\Delta y$ and $\Delta t=\Delta t_{\mathrm{CFL}}$ case, 20 iterations are sufficient to keep the low error for numerical dispersion relation, especially for the modes near $k=0$. The precondition for the matrix $\overleftrightarrow{A}$ or relaxation method can be used to accelerate the convergence and will be subjects of future reports.

The computational cost of the solver is $\mathcal{O}(N)$ where $N$ is the number of cells. For each iteration, the new quantity $\Delta\overrightarrow{E}_{(m+1)}$ (subscript in parentheses for iteration step, not time step or spatial grid) only depends on the old quantity $\Delta\overrightarrow{E}_{(m)}$ at the nearest neighbor and next nearest neighbor cells, so that in the domain decomposed mesh, the communication is local, i.e. communication only needed to get the quantities at the adjacent cells. Our Maxwell solver needs slightly more memory than Yee solver due to the auxiliary field variables. Comparing to the global FFT-based solvers, where the computing cost is $\mathcal{O}(N\log N)$ and the requirement for communication is non-local, our algorithm is less computational expensive and more scalable. In Fig. 2(b), we show the weak-scaling performance of our solver by testing without particles, for 100 time steps. Each processor resolves a $256\times 256$ grid, regardless of what number of processors $N_{\mathrm{proc}}$ is used. The runtime only increases from 14.32 seconds for $N_{\mathrm{proc}}=64$ to 15.75 seconds for $N_{\mathrm{proc}}=14641$, which implies that the slowdown of the speed is only $9\%$ and the scaling of the solver is close to $\mathcal{O}(N)$.

2.2.2 Equations for two dimensional case

The Maxwell equations on two dimensional grid in $x-y$ plane can be written down by setting $D_{z}\to 0$ and $\alpha_{z}\to 1$ in Eqs(1) to (3) and (40) to (42). We obtain

[TABLE]

The matrix elements of $\overleftrightarrow{A}$ on two dimensional grid in $x-y$ plane can be written down by setting $\bar{D}_{z}\to 0$ in Eqs(44) to (49). We obtain

[TABLE]

Using Alg 1, $E_{x}$ and $E_{y}$ need to be solved simultaneously because $A_{xy}\neq 0$, but $E_{z}$ can be solved independent of $E_{x}$ and $E_{y}$.

2.2.3 Equations for one dimensional case

The Maxwell equations on one dimensional grid in $x$ direction can be written down by setting $D_{y}\to 0$ and $\alpha_{y}\to 1$ in Eqs(51) to (56). We obtain

[TABLE]

The matrix elements of $\overleftrightarrow{A}$ on one dimensional grid in $x$ direction can be written down by setting $\bar{D}_{y}\to 0$ in Eqs(57) to (61). We obtain

[TABLE]

Using Alg 1, $E_{y}$ and $E_{z}$ can be solved independently and there is no need to solve for $E_{x}$ because $E_{x}=E_{x}^{\prime}$.

3 NCI growth rate

For a cold plasma drifting with ultra-relativistic velocity $v_{x}\to c$, we derive in Appendix E the NCI growth rate following Xu 2013[14]. In two dimensional grid in $x-y$ plane, the asymptotic expression for NCI growth rate is

[TABLE]

where $\omega_{p}=\sqrt{\frac{4\pi q^{2}n_{e}}{\gamma m_{e}}}$ is the relativistic electron plasma frequency. The notation $S$ with subscript is for the interpolation functions for the corresponding component of $\overrightarrow{E}$, $\overrightarrow{B}$ and $\overrightarrow{J}$. The averaging operators in Fourier space is given in Eq(11) and $[k]_{i}$ in Eq(16). The $\xi_{0}$, $\xi_{1}$, $\xi_{2}$, $\xi_{3}$ are defined in Eqs(132) to (135). The growth rate is only nonzero for the condition that $(\omega^{\prime},k_{x}^{\prime},k_{y}^{\prime})$ sits near the EM modes and beam modes

[TABLE]

For our momentum conserving scheme we have

[TABLE]

where

[TABLE]

Thus

[TABLE]

4 Numerical verifications

We implement our Maxwell solver in the PIC code EPOCH 2D[1], where we use the default Boris particle push and the options for particle shapes in the code. We use $m=20$ iterations for solving Eq(43).

4.1 2D dispersion relation for electromagnetic wave

To initialize the simulation, we only put a point source in one cell of the two dimensional domain, so that the spatial Fourier transform of field in the initial frame is nonzero for all wavevectors. In the simulations, we use $\Delta x=\Delta y$, and number of cells $N_{x}\times N_{y}=128\times 128$. The time step is $\Delta t=0.1\Delta t_{\mathrm{CFL}},0.5\Delta t_{\mathrm{CFL}},0.7\Delta t_{\mathrm{CFL}},1.0\Delta t_{\mathrm{CFL}}$. And the simulation is run for $t_{\mathrm{max}}=256\frac{\Delta x}{c}$. No particles are loaded in this test. In the post-processing, the discrete Fourier transform is performed for electric field in both space and time to test the numerical dispersion relation. To make the data periodic in time before the discrete Fourier transform, we apply the Hann window function

[TABLE]

The numerical value of the frequency for each $(k_{x},k_{y}$) is calculated to be the weighted value of $\omega$ in the discrete Fourier transform, i.e,

[TABLE]

where $\tilde{E}_{z}(\omega,k_{x,}k_{y})$ is the discrete Fourier transform of $E_{z}(t,x,y)$, $\omega$ is the discretized frequency in the discrete Fourier transform, and $(k_{x},k_{y})$ is the discretized wavevector. Since $E_{z}(t,x,y)$ is always real-valued, the transform in $t\to\omega$ satisfies $\tilde{E}_{z}(\omega)=\tilde{E}_{z}^{*}(-\omega)$, thus it is sufficient to calculate the $\omega\geq 0$ components. The results for phase velocities in the numerical tests are plotted in Fig. 1 (d) to (f) against the phase velocities calculated from Eq(34). The results of numerical experiments are consistent with Eq(34) for $\Delta t=0.5\Delta t_{\mathrm{CFL}},0.7\Delta t_{\mathrm{CFL}},1.0\Delta t_{\mathrm{CFL}}$, which suggests that the number of iterations is sufficient. The EM wave at high $\overrightarrow{k}$ becomes more superluminal, i.e. larger $v_{ph}$, for smaller time steps.

4.2 Location of NCI in $\overrightarrow{k}$ space and growth

rate

We carry out numerical tests to verify the location of NCI in $\overrightarrow{k}$ space and compare the growth rate in the test runs with the analytical expression in Eq(77). In our test runs, a pair plasma is initialize with a drifting speed $\beta$ in $x$ direction. And the Lorentz factor is $\gamma=\frac{1}{\sqrt{1-\beta^{2}}}$. The numerical parameters of the test runs are listed in Table 1. We use the method in [17] to load the particles with relativistic shifted-Maxwellian distribution, which takes care of the spatial part of the transformation $d^{3}x\to d^{3}x^{\prime}$, and the acceptance efficiencies are 50% for generic cases and 100% for symmetric distributions. We use momentum conserving scheme for field interpolation and current deposition in our simulations.

The locations of NCI resonance modes in $\overrightarrow{k}$ space are computed by solving Eq(71) and (72) for the main beam $\nu_{x}=0$ and first order aliasing beam $\nu_{x}=\pm 1$, as shown in Fig. 3 (a), where $\nu_{x}$ is the index for aliasing beam as in Eq(121). The smallest value of $|k_{x}|$ in the NCI modes is plotted in Fig. 3(e), showing that if $\Delta t<\frac{\Delta t_{\mathrm{CFL}}}{\sqrt{2}}$ only $\nu_{x}=\pm 1$ resonance exists, if $\frac{\Delta t_{\mathrm{CFL}}}{\sqrt{2}}<\Delta t<\Delta t_{\mathrm{CFL}}$ both $\nu_{x}=0$ and $\nu_{x}=\pm 1$ resonance exist. However, the simulation shows that the dominant NCI mode is from the $\nu_{x}=\pm 1$ resonance which indicates that an effective method to suppress such resonance is desired. The location of NCI modes is bounded by $\frac{2.06}{\Delta x}<|k_{x}|<\frac{\pi}{\Delta x}$ for all the valid time step $\Delta t\leq\Delta t_{\mathrm{CFL}}$. The NCI modes in the numerical experiments are visualized in Fig. 3 (b)-(d) by plotting the spatial Fourier transform of the out-of-plane magnetic field. The locations of NCI in $\overrightarrow{k}$ space are consistent between the theory and the numerical experiments as shown in Fig. 3 (f)-(h). The growth rate in Eq(77) is essentially zero on $k_{x}$ axis, which is also seen in Fig. 3 (b)-(d). There is no higher order aliasing $|\nu_{x}|\geq 2$, as solved by using Eq(71) and (72).

In the test runs for drifting pair plasma, the fastest growing NCI mode will dominate in the growth of the total field energy, thus the growth rate for the amplitude of the fastest growing NCI mode is half of the growth rate for the total field energy for WKB approximation. The growth rate of the fastest growing NCI mode can be calculated by maximizing Eq(77) over $k_{x}$ or $k_{y}$. We calculate the growth rate of the fastest growing NCI mode, i.e. the most unstable NCI mode, as shown in the curves in Fig. 3(i). In the same figure, we also show the growth rate calculated from the numerical experiments by calculating half of the growth rate of total field energy. Using high order B-spline particle shape significantly reduces the NCI growth rate, although more computational cost is required. The trend of the NCI growth rate are consistent between the analytical expression and the numerical results.

For the $\Delta t=0.5\Delta t_{\mathrm{CFL}}$ case, we use the low pass filtering with our solver and make the comparison with Yee solver, as shown in Fig. 4. The growth rate of NCI is significantly suppressed as shown in Fig. 4(a) and (c).

5 Conclusions and discussions

We demonstrate that the Maxwell solver from Eastwood[4] has several good numerical properties that allow for the efficient mitigation of NCI. The dispersion relation is always superluminal in a large range of $k$ including $k=0$, thus the NCI resonance are separated from the physical modes in $k$ space. There is only high $k$ NCI mode, i.e. $\frac{2.06}{\Delta x}<|k_{x}|<\frac{\pi}{\Delta x}$ for $\nu_{x}=0$ and $\nu_{x}=\pm 1$ NCI. Those high $k$ modes can be suppressed by applying a low-pass filter on the current during the simulation. An alternative approach to deal with the high $k$ mode is to post-process the fields and currents with a low-pass filter after the simulation. Although our algorithm requires solving linear equations using an iteration method, it is local, of $\mathcal{O}(N)$ computational cost and scalable for multi-dimensional simulation applications, which has advantage over the methods based on spatial Fourier transform. This solver also possess conservation properties such that no correction scheme is needed. The examples in this paper are using 2D grid, but our implementation and discussions on the numerical properties can be generalized to 3D grid. The averaging stencil we used in this work is $(\frac{1}{6},\frac{2}{3},\frac{1}{6})$. Another stencil can be potentially found to further improve the numerical properties. The rate of convergence may also be improved by using more optimized iteration scheme. These further improvements and generalizations will be subjects of future reports.

The Maxwell solver we developed was derived from finite element method in Eastwood 1991[4]. The particle pusher we use is the Boris push, which conserves the phase space volume for particles. If we derive the particle pusher using the finite element method in Eastwood 1991, then the particle pusher should be implicit. By the combination of that implicit particle pusher and the Maxwell solver we used in this work, the numerical solution of all the particle and field variables should strictly satisfy the minimal action principle and has better conservation properties than traditional PIC methods. As we learned from Eastwood 1991, the numerical conservation laws can be well defined and without error terms in the finite element method. With finite element methods, we can potentially construct numerical schemes which is structure preserving for both particles and fields in the relativistic case, as an alternative approach to the method derived from Hamiltonian for the non-relativistic case[12]. Further more, the finite element PIC is suitable for complex geometries.

6 Acknowledgements

Research presented in this paper was supported by the Center for Space and Earth Science(CSES) program and Laboratory Directed Research and Development(LDRD) program of Los Alamos National Laboratory(LANL). The simulations were performed with LANL Institutional Computing which is supported by the U.S. Department of Energy National Nuclear Security Administration under Contract No. 89233218CNA000001, and with the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation(NSF) grant number ACI-1548562.

Appendix A Conservation of the divergence of the fields

Starting from the L.H.S of Eq(8) and apply $D_{t}$ on it, and using Eqs(1) to (3), we obtain

[TABLE]

Thus Eq(8) is conserved during the evolution as long as it is conserved initially.

The time evolution of the difference between the L.H.S of Eq(7) can be calculated by using Eqs(4) to (6) and the continuity equation for charge conserving deposition scheme

[TABLE]

We obtain

[TABLE]

where the last step requires the use of a charge-conserving current deposition scheme. Thus Eq(7) is conserved during the evolution as long as it is conserved initially and the charge conserving deposition scheme is used. If one does not use a charge-conserving current deposition scheme, then the correction for the the Gauss’s law for electric field is needed, otherwise the finite grid instability can grow due to the loss of charge conservation[11].

Appendix B Derivation for the limit of time step and the superluminal property

of dispersion relation

For real value of $\omega$, the function $G(\frac{\omega\Delta t}{2})$ satisfies $0\leq G(\frac{\omega\Delta t}{2})\leq 3.$ From Eq(30) we need for real $\omega$ that

[TABLE]

So

[TABLE]

And $\Delta t$ needs to be less or equal to the minimum value on the R.H.S, which is reached when $G(\frac{k_{y}\Delta x}{2})=G(\frac{k_{y}\Delta x}{2})=G(\frac{k_{y}\Delta x}{2})=3$, thus we get Eq(31)

[TABLE]

The 2D case can be derived in the same way.

From (32) we have

[TABLE]

Using Eq(31) we have

[TABLE]

Using Cauchy inequality $(a_{1}^{2}+a_{2}^{2}+a_{3}^{2})(b_{1}^{2}+b_{2}^{2}+b_{3}^{2})\geq(a_{1}b_{1}+a_{2}b_{2}+a_{3}b_{3})^{2}$ and letting $a_{1}=k_{x}^{2}\Delta x$, $a_{2}=k_{y}^{2}\Delta y$, $a_{3}=k_{z}^{2}\Delta z$, $b_{1}=\frac{1}{\Delta x}$, $b_{2}=\frac{1}{\Delta y}$, $b_{3}=\frac{1}{\Delta z}$, then the R.H.S of Eq(87) is greater than or equal to zero, thus we have $\omega/k\geq c$. For the L.H.S. of Eq(87) to be equal to zero, we need $c\Delta t=\frac{1}{(1/\Delta x^{2}+1/\Delta y^{2}+1/\Delta z^{2})^{1/2}}$ and $\frac{a_{1}}{b_{1}}=\frac{a_{2}}{b_{2}}=\frac{a_{3}}{b_{3}}$, i.e. $|k_{x}|\Delta x|=|k_{y}|\Delta y=|k_{z}|\Delta z$. The 2D case can be derived in the same way.

Appendix C Derivation of the inequality for infinite norm of matrix $\overleftrightarrow{A}$

The infinite norm of matrix $\overleftrightarrow{A}$ is simply the maximum absolute row sum of the matrix. The absolute $x$ row sum is $\sum_{j=x,y,z}|A_{xj}|=|A_{xx}|+|A_{xy}|+|A_{xz}|$. The operators $\bar{D}_{x},\bar{D}_{y}$, etc are diagonalized in Fourier space, thus it is convenient to analyze the properties of matrix $\overleftrightarrow{A}$ in Fourier space. In Fourier space $\bar{D}_{x}=2i\sin(k_{x}\Delta x/2)$ satisfies $|\bar{D}_{x}|\leq 2$. Similarly, $|\bar{D}_{y}|\leq 2$ and $|\bar{D}_{z}|\leq 2$. And $\bar{D}_{xx}=-4\sin^{2}(k_{x}\Delta x/2)$ satisfy $-4\leq\bar{D}_{xx}\leq 0$. Similarly, $-4\leq\bar{D}_{yy}\leq 0$ and $-4\leq\bar{D}_{zz}\leq 0$. Thus

[TABLE]

We use Eq(31), which makes $1-\frac{c^{2}\Delta t^{2}}{\Delta y^{2}}>0$, $1-\frac{c^{2}\Delta t^{2}}{\Delta z^{2}}>0$, $1-\frac{c^{2}\Delta t^{2}}{\Delta y^{2}}-\frac{c^{2}\Delta t^{2}}{\Delta z^{2}}>0$ and $\frac{1-\frac{c^{2}\Delta t^{2}}{\Delta y^{2}}}{1-\frac{c^{2}\Delta t^{2}}{\Delta y^{2}}-\frac{c^{2}\Delta t^{2}}{\Delta z^{2}}}>1.$ The maximum of $A_{xx}$ is 0 when $\bar{D}_{yy}=\bar{D}_{zz}=0$, thus $|A_{xx}|=-A_{xx}$. Using \bigg{|}\bar{D}_{zx}\bigg{|}\leq 4, \bigg{|}\bar{D}_{xy}\bigg{|}\leq 4 we have

[TABLE]

Thus for the absolute $x$ row sum

[TABLE]

The R.H.S. of Eq(91) is monotonic in $\bar{D}_{yy}$ and $\bar{D}_{zz}$ as derived in Appendix C.1, and reaches maximum when $\bar{D}_{yy}=\bar{D}_{zz}=-2$ , thus

[TABLE]

When $\Delta t\leq\frac{\Delta x}{\sqrt{2}c}$ we have $\sum_{j=x,y,z}|A_{xj}|<1$. Similarly, we have $\sum_{j=x,y,z}|A_{yj}|<1$ if $\Delta t\leq\frac{\Delta y}{\sqrt{2}c}$ and we have $\sum_{j=x,y,z}|A_{zj}|<1$ if $\Delta t\leq\frac{\Delta z}{\sqrt{2}c}$. Thus the infinite norm of matrix $A$ is less than $1$ if $\Delta t\leq\frac{\min(\Delta x,\Delta y,\Delta z)}{\sqrt{2}c}$. In the $\Delta x=\Delta y=\Delta z$ case the infinite norm of matrix $A$ is always less than $1$ because $\Delta t\leq\frac{\Delta x}{\sqrt{3}c}$.

For 2D case, from Eqs(57) to (61) we have

[TABLE]

Similarly we can prove $\sum_{j=x,y}|A_{yj}|<1$. And

[TABLE]

For 1D case we have

[TABLE]

To calculate the speed of convergence of the iteration, we need to compute the eigenvalues of matrix $\overleftrightarrow{A}$. We derive in Appendix D how the error of the dispersion relation depends on the number of iterations. For symmetric matrix $\overleftrightarrow{A}$, the absolute value of eigenvalue $|\lambda|$ is always smaller or equal to the infinity norm, thus $|\lambda|<1$ and $(\mathbb{1}-\overleftrightarrow{A})$ is invertible if the infinity norm is less than 1.

C.1 Derivation of the inequality for the infinite norm of $\overleftrightarrow{A}$

in 3D

The R.H.S of Eq(91) can be written as

[TABLE]

where the constants $a_{y1}$, $a_{z1}$, $a_{xy}$, $a_{y2}$, $a_{z2}$ are

[TABLE]

and the variables $d_{y}=-\bar{D}_{yy}/6$ and $d_{z}=-\bar{D}_{zz}/6$ satisfy $0\leq d_{y}\leq\frac{2}{3}$ and $0\leq d_{z}\leq\frac{2}{3}$. By take the derivative of $f$ w.r.t $d_{y}$ and $d_{z}$ we can check the monotonicity, e.g.

[TABLE]

Since \big{(}1-\frac{c^{2}\Delta t^{2}}{\Delta y^{2}}-\frac{c^{2}\Delta t^{2}}{\Delta z^{2}}\big{)}>0, $\frac{\partial f}{\partial d_{y}}$ reaches minimum at $d_{z}=\frac{2}{3}$

[TABLE]

Thus $f$ increases as $d_{y}$ increases. Similarly $f$ increases as $d_{z}$ increases.

Appendix D Speed of convergence and error term for iteration

We constrain the derivation to the case where the infinity norm of matrix $\overleftrightarrow{A}$ is less than 1, which is always true in 1D or 2D. Using finite number of iterations for Alg 1, we have the error term, which is the difference between the increment electric field we get after $m$ iterations and after infinite number of iterations

[TABLE]

where we used in the last step

[TABLE]

which is true for the exact solution $\Delta\overrightarrow{E}_{(\infty)}$ in Eq(43). Thus by mathematical induction from Eq(100) we have

[TABLE]

where we used $\Delta\overrightarrow{E}_{(0)}=\Delta\overrightarrow{E}^{\prime}$ in the last step. Multiply Eq(102) on the left by $(\mathbb{1}-\overleftrightarrow{A})$ and use Eq(101), we obtain

[TABLE]

Thus

[TABLE]

Because the eigenvalue of matrix $\overleftrightarrow{A}$ has $|\lambda|<1$, Eq(104) always converge to $(\mathbb{1}-\overleftrightarrow{A})^{-1}\cdot\Delta\overrightarrow{E}^{\prime}$ with $m\to\infty$. Since Eq(104) holds for each time step, we have

[TABLE]

From the Fourier transform of Eqs(40) to (42), we obtain

[TABLE]

where

[TABLE]

is the spatial averaging operator and $\omega_{(m)}$ is the frequency of electromagnetic wave with $m$ iterations. In Eqs(1) to 3, we use $\overrightarrow{E}_{(m)}$ as an approximation for $\overrightarrow{E}_{(\infty)}$ to update $\overrightarrow{B}$, thus

[TABLE]

From Eqs(105) and (106) and (108) we can eliminate $\overrightarrow{B}$ and $\overrightarrow{E}^{\prime}$ and get the equation for $\overrightarrow{E}_{(m)}$

[TABLE]

For $m\to\infty$

[TABLE]

If $\overrightarrow{E}_{(m)}$ is a eigenvector of $\overleftrightarrow{A}$ with eigenvalue $\lambda$, then from Eq(105) we know that both $\overrightarrow{E}_{(\infty)}$ and $\overrightarrow{E}_{(0)}$ are eigenvectors of $\overleftrightarrow{A}$ with eigenvalue $\lambda$. By comparing Eq(109) and Eq(110) we obtain

[TABLE]

where $\omega_{(\infty)}$ is the frequency of electromagnetic wave with infinite number of iterations, which is given by Eq(28). We can rewrite Eq(111) as

[TABLE]

With sufficient number of iterations, $\omega_{(m)}$ is close to $\omega_{(\infty)}$, thus

[TABLE]

For each set of $(k_{x},k_{y},k_{z})$, $\overleftrightarrow{A}$ has three eigenvalues $\lambda_{1},\lambda_{2},\lambda_{3}$. And $\lambda_{1},\lambda_{2},\lambda_{3}$ are functions of $(k_{x},k_{y},k_{z})$. Let $\lambda_{3}$ be the eigenvalue with the largest absolute value among $\lambda_{1},\lambda_{2},\lambda_{3}$, i.e. $|\lambda_{1}|\leq|\lambda_{2}|\leq|\lambda_{3}|$, then with sufficient number of iterations the dominant error term for dispersion relation is

[TABLE]

Appendix E Asymptotic NCI growth rate

Following the derivation in Xu 2013[14], we derive the asymptotic expression for NCI growth rate in 2D. We use the relativistic electron plasma frequency $\omega_{p}=\sqrt{\frac{4\pi q^{2}n_{e}}{\gamma m_{e}}}$ instead of the expression $\omega_{p}=\sqrt{\frac{4\pi q^{2}n_{e}}{m_{e}}}$. The definition of $\omega_{p}$ with $\gamma$ has physical meaning of the oscillation frequency in cold drifting plasma. Eq(10) in Xu 2013 and Eq(20) in this paper are the same for the R.H.S, but on the L.H.S of Eq(20) in this paper, the vacuum dielectric tensor are different from Xu 2013 and is written down in Eqs(21) to (26). In this paper we have $[\overrightarrow{k}]_{E}=[\overrightarrow{k}]_{B}=[\overrightarrow{k}]$ written down in Eq(16). The dielectric tensor for the current part of drifting plasma is same as that in Xu 2013. We ignore high order terms in $1/\gamma$, in $v_{0}\to c$ limit

[TABLE]

where

[TABLE]

are the aliased frequency and wavevector, the interpolation functions have dummy variables with aliasing wavevector. We expand $\omega^{\prime}$ around the beam resonance $\omega^{\prime}=k_{x}^{\prime}$ and write $\omega^{\prime}=k_{x}^{\prime}+\delta\omega^{\prime}$, where $\delta\omega^{\prime}$ is a small term. Denoting $Q_{ij}$ as the numerator $\epsilon_{ij}^{(J)}$ we can write $\epsilon_{ij}=\epsilon_{ij}^{(\mathrm{vac})}+\epsilon_{ij}^{(J)}$ as

[TABLE]

Using $\det(\epsilon)=0$ and only keeping $x$ and $y$ terms, we can obtain

[TABLE]

where

[TABLE]

Now we use the condition that $(\omega^{\prime},k_{x}^{\prime})$ sits near the EM modes

[TABLE]

and expand the finite difference operator $[\omega]=\xi_{0}+\delta\omega^{\prime}\xi_{1}$ , and $\alpha_{t}=\xi_{2}+\delta\omega^{\prime}\xi_{3}$ where

[TABLE]

[TABLE]

and then we further expand $[\omega]$ to first order in $A_{1}$, since this term is sensitive near the EM mode, while only keeping the zeroth order of $[\omega]$ in $B_{1}$ and $C_{1}$

[TABLE]

We then obtain a cubic equation, we can drop the small $B_{1}$ term

[TABLE]

So the growth rate

[TABLE]