Numerical instability due to relativistic plasma drift in EM-PIC simulations

We mainly follow the notation in Ref. [4] to derive the numerical dispersion relation for a cold plasma drifting with relativistic speeds. We note that the multi-dimensional analysis in [4] solves for $\phi$ and $\vec{A}$ and is not valid for relativistic drifts, and the 1D analysis in [3] predicts growth in 1D. On the other hand, our analysis includes relativistic mass effects, is valid in multi-dimensions, and it predicts instability only in multi-dimensions (in agreement with our simulations). Since most EM-PIC codes now in use solve for the electric field $\vec{E}$ and magnetic field $\vec{B}$ directly (with finite difference or spectral solvers), we derive a numerical dispersion relation directly using these two quantities. Gaussian units will be used; in addition, particle mass and velocity will be normalized to electron mass and the speed of light.

For a multi-dimensional simulation setup in Cartesian coordinates, the EM field interpolated on a particle can be expressed as

where $m$ is the time index and $\vec{n}$ is the grid index; $\overleftrightarrow{S}$ is the interpolation dyadic used to obtain the appropriate field at $\vec{x}$ and $t=m\Delta t$; $\vec{E}_{m,\vec{n}}$ and $\vec{B}_{m,\vec{n}}$ stands for the electromagnetic forces at time grid index $m$ and space grid index $\vec{n}$. For momentum conserving field interpolation $\overleftrightarrow{S_{E}}$ and $\overleftrightarrow{S_{B}}$ are equal and are scalar functions times the unit dyadic while for energy conserving field interpolation $\overleftrightarrow{S_{E}}$ and $\overleftrightarrow{S_{B}}$ are not equal in each direction (the interpolation dyadics for $\vec{E}$, $\vec{B}$, and $\vec{j}$ are given in A). The momentum change of the particle is related to the change in the distribution function of the plasma by the linearized Vlasov equation

where $\vec{p}$ is the particle momentum, and $\gamma$ is the particle Lorentz factor. After Fourier transforming, the Vlasov Equation becomes

Note that $\vec{E}$ and $\vec{B}$ are defined at the discrete grid position and discrete time step, so its Fourier transform in $(w,k)$ is periodic, i.e.,

where

Note that when the EM field are staggered (such as on a Yee mesh), there is an addition $(-1)^{\sum_{i}\nu_{i}}$ term in $\vec{E}(\omega^{\prime},\vec{k}^{\prime})$ and $\vec{B}(\omega^{\prime},\vec{k}^{\prime})$, where $\hat{i}$ is summed over the directions for which the EM fields are half-grid offset [5]; we absorb these additional terms into $\overleftrightarrow{S_{E}}$ and $\overleftrightarrow{S_{B}}$ to keep Eq. (2) correct (see A).

Replacing $(w,\vec{k})$ with $(w^{\prime},\vec{k}^{\prime})$ in Eq. (1), and using Eq. (2.1), we obtain

The current density $\vec{j}$ due to the movement of the particles can be expressed as

where $\overleftrightarrow{S_{j}}(\vec{x}^{\prime}-\vec{x})$ is the dyadic for the current deposit. After Fourier transforming we obtain

We can now proceed in the normal way to obtain a dispersion relation. We start from Faraday’s and Ampere’s Law,

which upon Fourier transforming gives,

$[k]_{E}$ and $[k]_{B}$ are the finite difference operator of the corresponding Maxwell solver schemes for $\vec{E}$ and $\vec{B}$ fields. We follow the notation in Ref. [4], and use $[\cdot]$ exclusively to indicate the finite difference operator. Applying $[\vec{k}]_{B}\times$ to both sides of Eq. (6), we end up with the coupled wave equation for $\vec{E}$ and $\vec{j}$,

Using Eq. (7) and (5), we could obtain

and if we normalize the distribution function such that $f_{0}=n_{0}f^{n}_{0}$, use the definition of plasma frequency

and use the expression for the distribution function in Eq. (4), we finally obtain

which is a generalized dispersion relation for a plasma of finite size particles drifting on a grid. We note that the use of additional smoothers and filters can be incorporated into the dispersion relation by adding additional $S_{SM}(\vec{k}^{\prime})$ terms outside the summation over Brillouin zones (essentially it multiplies the $\omega^{2}_{p}$ term).

We next examine the dispersion relation in the limit of a cold plasma including the possibility that the drift is near the speed of light.

The 1D, 2D, and 3D dispersion relations show that a drifting plasma leads to beam modes in the dispersion relation that are associated with longitudinal oscillations. A beam mode roughly satisfies the dispersion relation

In addition, a drifting plasma also supports transverse EM waves that are described by the dispersion relation

which for the high gammas are the dispersion relation on the grid for transverse EM modes in vacuum. In 3D, as the three components of the electric field are coupled, we expect to find instabilities near the intersections of the EM modes and beam modes in ($k_{1}$, $k_{2}$, and $k_{3}$) space. However, in 2D $E_{3}$ is de-coupled from $E_{1}$ and $E_{2}$; thus the instability (coupling) can only occur in $k_{1}$ and $k_{2}$ space.

In 1D, all the three components of the electric field are de-coupled so instability can only occur if either the longitudinal beam or transverse EM mode are numerically unstable themselves. We next concentrate on the 2D case since it is easier than the 3D but it still has the possibility of numerically unstable modes.

As we show next, we observe in the simulations and in the solutions to our dispersion relation that in fact the unstable wave numbers and frequencies lie at the intersection of $\omega^{\prime}-k^{\prime}_{1}v_{0}=0$ (which we refer to as a vacuum beam mode) and the vacuum dispersion relation for EM waves. We will use this observation to derive asymptotic growth rates for the instability in section 5.

Without loss of generality, we use the results in section 2 to study the numerical instability induced by the relativistic plasma drift in a 2D system. According to the dispersion relation in 2D, we expect to observe instability in $E_{2}$ (and $B_{3}$). By calculating the maximum imaginary part of $\omega$ for real values of ($k_{1},k_{2})$ for Eq. (20), we can obtain the characteristic pattern of the instability in Fourier space, as well as the growth rate of the instability. We can also plot the real part of $\omega$ for $k_{1},k_{2}$. These results can be used later to compare with the simulation results.

The dispersion relation is general and can be used to examine different choices in Maxwell Equation solvers, in differences between energy and momentum conserving field interpolation, in differences between charge conserving and direct current deposition schemes, and the use of smoothing and low pass filters. In this paper, we are emphasizing that our dispersion relation agrees well with the simulation results for cases studied, that we can predict the region of unstable modes by plotting where the beam and EM modes intersect in $\vec{k}$ and $\omega$ space; that we can obtain an asymptotic expression for the growth in 3D which agrees well with the simulation for various finite difference solvers (including the values of $\Delta t$ tht minimize the growth rate); and the advantages of using a spectral solver from the point of view of eliminating the instability and not attempting to carry out a comprehensive survey of all available choices listed above.

We illustrate the instability using a 2D case with the standard Yee solver [15]. We choose the grid parameters and time step that satisfies the Courant Condition [2] to eliminate known numerical instability from the EM modes. We use the parameters in Tab. 1, and substitute the finite difference operator for the Yee solver (see A) into the 2D dispersion relation. We assume linear (area) interpolation, momentum conserving field interpolation, and a charge conserving current deposition (see Appendix A).

After obtaining all the roots $(\omega,k_{1},k_{2})$, we plot the dependence of the growth rate in the $(\omega_{r},k_{1})$ space [figure 1 (c)], as well as in the $(k_{1},k_{2})$ space [figure 1 (d)]. It is evident that all the instabilities are near the main or aliased beam modes. Since the terms with $|\mu|\leq 1$, and $|\nu_{1}|\leq 1$ are the most important, we neglect higher order terms when solving Eq. (20). Higher order $\mu$ and $\vec{\nu}$ terms can be included in the summation if needed. These additional terms lead to additional unstable modes in $(k_{1},k_{2})$ space with lower growth rates as well as to small modifications to the growth rate and location of the original modes. A plot in $(k_{1},k_{2})$ space with more terms included are presented in figure 1 (b), in which we use the asymptotic expression Eq. (5.1) for the growth rate (see section 5 for more details).

While the results in figure 1 (c) and (d) are numerically calculated from Eq. (20), the location of the unstable modes can also be conveniently predicted by plotting the intersection of the EM modes Eq. (22) and beam modes Eq. (21) in $(k_{1},k_{2},\omega_{r})$ space. This is shown in a 3D plot [figure 1 (a)]. By examining the unstable pattern in $(k_{1},k_{2})$ space we see that the central part of pattern comes from the intersections of the EM modes and main beam mode ($\mu=0$ and $\nu=0$), while the part at the four corners can be identified from the intersections of the EM modes and first order spatial aliasing beam modes ($\mu=0$ and $\nu_{1}=\pm 1$). As we argue in section 5, a key to mitigating the instability is to manipulate the instability pattern through a careful choice of the Maxwell Equation solver. Making a plot in $(k_{1},k_{2},\omega_{r})$ of the intersection of the EM and beam modes for various solvers becomes a useful method for examining where the unstable modes reside without having to solve the full dispersion relation.

To compare with the results in section 3.1, we conducted simulation studies in the 2D system using the EM-PIC code OSIRIS [13]. In these simulations, a neutral plasma with both the ion and electrons drifting in $x_{1}$ at the same relativistic velocity of $\gamma=50.0$ is initialized throughout the entire simulation box. Periodic boundary conditions for fields and particles are used. Other crucial parameters for the simulation setup are identical to the theoretical study in section 3.1.

As is shown in figure 2 (a), the total EM energy starts to grow violently as the plasma drifts relativistically. The exponential growth indicates that a numerical instability is occurring. In addition, the EM field energy in $E_{2}$ and $B_{3}$ and that in $E_{3}$ and $B_{2}$ are shown separately. As predicted by the 2D dispersion relation the $E_{3}$ and $B_{2}$ modes are stable and do not grow. The pattern of $E_{2}$ at $t=100~\omega^{-1}_{p}$ is plotted in figure 1 (e) and (f), and good agreement for the location and relative amplitude of the unstable modes is obtained when compared against the theoretical prediction [figure 1 (c) and (d)].

The EM energy grows with a lower rate after $t=110~\omega^{-1}_{p}$ [figure 2 (a)]. The plasma density in this regime is highly modulated by the EM fields. The first order perturbation in plasma electron density [figure 2 (b)] shows a similar pattern as for $E_{2}$ [figure 2 (d)] , which confirms they are coupling in the system. Note that no exponential energy growth can be seen in the $E_{3}$ field [figure 2 (c)]

From the simulation we find that for later times after the instability has evolved into a nonlinear state, the same pattern in $(k_{1},k_{2})$ space as that of the linear regime still exists. This indicates that the instability will remain near the intersections of the EM modes and beam modes and that both the linear and nonlinear growth can be mitigated through eliminating or controlling the intersections.

We also carried out a numerical investigation of the 1D dispersion relation Eq. (2.2.2), and (18) using the same simulation parameters as in Tab. 1 (with the 1D Courant condition). This confirms that there is no numerical instability under these parameters which is expected since $E_{1}$ is de-coupled from $E_{2}$ and $E_{3}$ in Eq. (2.2.2) and each mode is itself stable.

We have done the same simulation studies on the use of other finite difference solvers beside the Yee solver, including the Karkkainen [8, 16] and 4th-order solver [17]. Good agreement between theory and simulations was also found [18]. Some results are shown in figure 5 and discussed in section 5.2.

In section 3, we obtained the instability pattern and growth rate by numerically solving the dispersion relation equation, which is feasible on a modern laptop computer in 2D scenario, but much more difficult in 3D. As mentioned above, we observed the highest growth rate at the intersections of the beam modes and EM modes. Taking advantage of this observation, we here derive an asymptotic expression for the solutions near the beam mode $\omega^{\prime}=k^{\prime}_{1}v_{0}$. These expressions will not only speed up the instability pattern analysis in 3D, but also provide more insights into the dependence of instability pattern and growth rate to the grid sizes and time step used in simulation.

We denote the terms which are summing over $\mu$ and $\vec{\nu}$ in $\epsilon_{ij}$ as $Q_{ij}$, i.e.

We expand $\omega^{\prime}$ around the beam mode $\omega^{\prime}=k^{\prime}_{1}v_{0}$, and write $\omega^{\prime}=k^{\prime}_{1}v_{0}+\delta\omega^{\prime}$, where $\delta\omega^{\prime}$ is a small term. In addition, we will use the relativistic limit $v_{0}\rightarrow 1$. We will also truncate the summation over $\nu_{2}$ and $\nu_{3}$, keeping only the $\nu_{2}=\nu_{3}=0$ terms. Using $\det(\overleftrightarrow{\epsilon})=0$, and dropping terms of higher order of $(\omega^{2}/\gamma)^{2}$, $(\omega^{2}/\gamma)^{3}$, $\ldots$, we can obtain

where

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