SHARP: A Spatially Higher-order, Relativistic Particle-in-Cell Code

TL;DR

SHARP introduces a higher-order relativistic PIC algorithm that significantly reduces numerical heating and improves long-term plasma simulation accuracy by conserving charge and momentum exactly.

Contribution

The paper presents a novel fifth-order accurate relativistic PIC algorithm that enhances simulation fidelity and energy conservation compared to traditional second-order schemes.

Findings

Accurately reproduces kinetic features of linear plasma regimes.

Reduces numerical heating and energy non-conservation.

Demonstrates the importance of combined particle number increase and cell size reduction for convergence.

Abstract

Numerical heating in particle-in-cell (PIC) codes currently precludes the accurate simulation of cold, relativistic plasma over long periods, severely limiting their applications in astrophysical environments. We present a spatially higher-order accurate relativistic PIC algorithm in one spatial dimension, which conserves charge and momentum exactly. We utilize the smoothness implied by the usage of higher-order interpolation functions to achieve a spatially higher-order accurate algorithm (up to fifth order). We validate our algorithm against several test problems -- thermal stability of stationary plasma, stability of linear plasma waves, and two-stream instability in the relativistic and non-relativistic regimes. Comparing our simulations to exact solutions of the dispersion relations, we demonstrate that SHARP can quantitatively reproduce important kinetic features of the linear regime. Our simulations have a superior ability to control energy non-conservation and avoid numerical heating in comparison to common second-order schemes. We provide a natural definition for convergence of a general PIC algorithm: the complement of physical modes captured by the simulation, i.e., those that lie above the Poisson noise, must grow commensurately with the resolution. This implies that it is necessary to simultaneously increase the number of particles per cell and decrease the cell size. We demonstrate that traditional ways for testing for convergence fail, leading to plateauing of the energy error. This new PIC code enables us to faithfully study the long-term evolution of plasma problems that require absolute control of the energy and momentum conservation.

Tables

Table 1. Table 1: Relativistic Two-stream Simulation Parameters.

$u_b/c$	$v_b/c$	$N_c$ aaNumber of cells.	$N_{\mathrm{pc}}$ bbNumber of computation particles per cell.	${\Gamma }_m/{\omega }_p$ ccTheoretical predictions for the maximum growth rate.	$k_{m}c/{\omega }_p$ ddFastest growing wave mode predicted theoretically.	${\hat{\omega }}_m^1$ and ${\hat{\omega }}_m^2$ ee$\hat{\omega }=\omega /{\omega }_p$: other solutions given by Equation (57) at the fastest growing wave mode $k_m$.	${\hat{\omega }}_m^3$ ee$\hat{\omega }=\omega /{\omega }_p$: other solutions given by Equation (57) at the fastest growing wave mode $k_m$.
1	0.707	2441	491.6	0.2102	0.5147	$\pm 0.8142$	$-0.2102i$
2	0.894	6138	488.76	0.1057	0.2047	$\pm 0.4095$	$-0.1057i$
3	0.948	10948	493.24	0.0629	0.1148	$\pm 0.2435$	$-0.0629i$
4	0.970	16668	497.96	0.0422	0.0754	$\pm 0.1636$	$-0.0422i$

Equations

$${\partial }_t{f}_s(x,u,t)+\frac{u}{\gamma }{\partial }_x{f}_s(x,u,t)+\frac{q_{s}E(x,t)}{m_s}{\partial }_u{f}_s(x,u,t)=0,$$

$${\partial }_{x}E(x,t)-\frac{\rho (x,t)}{{ϵ}_0}=0\text{and}{\partial }_{t}E(x,t)+\frac{j(x,t)}{{ϵ}_0}=0,$$

$$\rho (x,t)$$

$$j(x,t)$$

$$f_s^K(x,u)=\sum _{i_s}^{N_s^p}\delta (x-x_{i_s})\delta (u-u_{i_s}),$$

$$f_s(x,u)=w\sum _{i_s}^{N_s}S(x-x_{i_s})\delta (u-u_{i_s})\approx \int d{x}^{\prime }{f}_s^K(x^{\prime },u)S(x^{\prime },x),$$

$$\int S(x,x_{i_s})dx=1.$$

$$\frac{k_B{T}_s}{M_s{c}^2}=\frac{\int du(\gamma -1)f_s}{\int du{f}_s}=\frac{\int du(\gamma -1)f_s^K}{\int du{f}_s^K}=\frac{k_B{T}_s^p}{m_s{c}^2}.$$

$$\frac{d{x}_{i_s}}{dt}$$

$$E_{i_s}$$

$$\rho (x,t)$$

$$j(x,t)$$

$$S(x,x_i)\to S^1\left(\frac{|x-x_i|}{\Delta x}\right)=\frac{1}{\Delta x}\{\begin{array}{l}1,\text{If }|x-x_i|<\Delta x/2,\\ \\ 0,\text{otherwise. }\end{array}$$

$$E_{k+1}^n-E_k^n$$

$${\rho }_{k+1/2}^n$$

$$W^m\left(\frac{x_k-x_{i_s}^n}{\Delta x}\right)$$

$$E_k^n=E_0^n+\frac{\Delta x}{{ϵ}_0}\sum _{j=0}^{k-1}{\rho }_{j+1/2}^n.$$

$${\partial }_t{E}_{\mathrm{tot}}$$

$$E_{\mathrm{tot}}^n$$

$$\sum _{k=0}^{N_c-1}k\left(E_{k+1}^n-E_k^n\right)=\sum _{k=1}^{N_c}(k-1)E_k^n-\sum _{k=1}^{N_c-1}k{E}_k^n=N_c{E}_{N_c}^n-E_{\mathrm{tot}}^n.$$

$$E_0^n=\frac{E_{\mathrm{tot}}^n}{N_c}+\frac{\Delta x}{{ϵ}_0{N}_c}\sum _{k=0}^{N_c-1}k{\rho }_{k+1/2}^n.$$

$$E_{i_s}^n$$

$$I_k$$

$$E_{k+\frac{1}{2}}$$

$$E_{i_s}^n$$

$$u_{i_s}^{n+1/2}$$

$$x_{i_s}^{n+1}$$

$$\begin{array}{r}{t}_0=\sqrt{m_0{ϵ}_0/(q_0^2{n}_0)},{\mathbb{E}}_0=\sqrt{n_0{m}_0{c}^2/{ϵ}_0},\\ {\rho }_0=q_0{n}_0,j_0={\rho }_{0}c,x_0=c{t}_0.\end{array}$$

$${\omega }_0^2\equiv t_0^{-2}=\frac{q_0^2{n}_0}{m_0{ϵ}_0}={\omega }_p^2=\sum _s\frac{Q_s^2{n}_s}{M_s{ϵ}_0}\to n_0=\sum _s\frac{{\bar{Q}}_s^2{n}_s}{{\bar{M}}_s},$$

$$n_0\Delta x=\sum _s\frac{{\bar{Q}}_s^2{n}_s\Delta x}{{\bar{M}}_s}=\frac{1}{N_c}\sum _s\frac{{\bar{Q}}_s^2{N}_s}{{\bar{M}}_s}.$$

Full text

1 Introduction

The PIC method is a very powerful numerical tool to study the evolution of plasmas, it is used to model plasmas ranging from laboratory experiments to astrophysical environments. First proposed in one spatial dimension (1D) by Buneman (1959) and Dawson (1962), the general idea of this algorithm is straightforward: it follows the trajectory of particles with $N$-body methods, while solving Maxwell’s equation on a Eulerian grid. The communication between grid points and particles is achieved through interpolation. The general loop (described in Figure 1) consists of first interpolating particle positions and velocities to a spatial grid to solve for the resulting charge and current density. Maxwell’s equations are then solved on the grid with these source terms to find the self-consistent electromagnetic (electrostatic in 1D) fields. Fields are then interpolated back to the particle positions to calculate the Lorentz force, and hence, acceleration, to forward evolve the particles in time using a so-called pusher. This reduces the number of computational operations from $\sim O(N^{2})$ (such as in the case of $N$-body methods) to $\sim O(N)$, where $N$ is the number of particles in the simulation. This also results in eliminating all wave modes in the electromagnetic fields on scales smaller than the cell size on the grid upon which they are computed.

A major test of the accuracy and fidelity of different PIC schemes is their ability to preserve conserved quantities such as energy, momentum, and charge. Often, this is required to accurately study subdominant, relativistic populations that typically arise in astrophysical contexts. Examples include nonthermal particle populations accelerated at shocks and reconnection events (Spitkovsky, 2008; Lyutikov et al., 2016), propagation of cosmic rays (Riquelme & Spitkovsky, 2009), interaction of accretion disks and coronae (Miller & Stone, 1999), and TeV blazar driven beam instabilities (Broderick et al., 2012). In the latter, this problem is especially severe, with the beams being both numerically and energetically subdominant while being highly relativistic; even a small degree of heating in the background can impact or overwhelm the evolution of beam plasma instabilities.

Direct interpolation of particle data to construct charge and current densities on the grid, in general, leads to a violation of charge conservation. However, modern algorithms typically use charge conserving methods to perform such step while maintaining the charge conservation (e.g., Eastwood, 1991; Villasenor & Buneman, 1992; Esirkepov, 2001; Umeda et al., 2003). Energy and momentum conservation on the other hand appear to be mutually exclusive (see, e.g., Brackbill, 2016). Due to their importance, several schemes have been developed to ameliorate their non-conservation using different underlying methodologies.

Recently introduced implicit methods (Chen et al., 2011; Lapenta & Markidis, 2011; Markidis & Lapenta, 2011) can in theory preserve total energy exactly (though in practice may not), while violating momentum conservation. For computationally simpler explicit, momentum conserving schemes, energy conservation is improved by filtering the deposited grid moments (charge and current densities), as is done in the TRISTAN-MP code (Buneman et al., 1993; Spitkovsky, 2005). However, filtering, when used with a momentum conserving scheme, leads to a violation of momentum conservation and non-vanishing self-forces (see Appendix A.3). In lieu of filtering, energy conservation is also improved by using higher-order interpolation functions, which is the approach that we will adopt below for both forward- (from particles to fields) and back-interpolation (from fields to particles) steps. This has the added advantage that momentum conservation is maintained.

In this paper, we describe an implementation of the PIC algorithm in 1D that uses high-order spline functions (up to the fifth order) for the forward- and back-interpolation steps of the algorithm. We couple this with an exact Poisson solver and a second-order symplectic integrator, i.e., leap frog, to produce a second-order accurate code called SHARP-1D. SHARP-1D displays superior energy-conservation properties while conserving the momentum exactly. The smoothness coming from the usage of high-order interpolation functions is utilized to construct an up to fifth-order spatially accurate algorithm.111Due to the usage of lower order interpolation functions, current algorithms perform the back-interpolation step with second-order spatial accuracy (e.g., Lapenta & Markidis, 2011; Haugbølle & Frederiksen, 2013; Brackbill, 2016; Lapenta, 2016). Despite the high spatial order accuracy, our code remains second-order accurate because it is limited by the accuracy of the particle pusher. This will be addressed in future work.

The paper is organized as follows. In Section 2, we describe the basic equations of the PIC method, our choice of discretization of the equations, and discuss sources for the numerical error. After discussing the order of accuracy of the solution, we describe our choice of code units and the implementation of SHARP-1D. In Section 3, we discuss the conservation properties of different PIC algorithms. In Section 4, we demonstrate the different capabilities of SHARP-1D by validating it against several test problems: thermal stability of plasma, plasma oscillation frequency and linear-Landau damping of standing plasma waves, and two-stream instabilities in both relativistic and non-relativistic regimes. We compare several of the results of SHARP-1D to the results of TRISTAN-MP simulations in Section 5. Finally, we study the convergence properties for our algorithm in Section 6 and discuss the performance of SHARP-1D in Section 7. We conclude in Section 8.

2 The PIC method

The evolution of the particles that comprise a plasma is described by the Boltzmann and Maxwell’s equations. In the absence of collisions, the particles are described by the Vlasov equation, which in one spatial dimension is

[TABLE]

where $s$ denotes a particle species, characterized by its charge $q_{s}$ and mass $m_{s}$, $u=\gamma v$ is the spatial component of the four-velocity, $\gamma=\sqrt{1+(u/c)^{2}}$ is the Lorentz factor, $E(x,t)$ is the electric field, and $f_{s}(x,u,t)$ is the phase-space distribution functions of particles of species $s$.

In one dimension, Maxwell’s equations imply that the magnetic field is constant and along the direction of the particle motion. Thus, it impacts the evolution of neither the particle nor the electric field, and we take it to be zero henceforth. Therefore, Maxwell’s equations reduce to

[TABLE]

where $\rho$ and $j$ are the charge and current densities, respectively. 222Here, we use SI units, and $\epsilon_{0}$ is the permittivity of the vacuum. To convert to CGS units, $\epsilon_{0}$ needs to be replaced by $1/4\pi$. Ultimately, we convert to a system of numerical units, obviating the distinction between these (Section 2.8). This set of equations is closed by the following equations for $\rho$ and $j$:

[TABLE]

2.1 Smoothing Phase-space Distribution Function

The distribution function for point-like particles is given by the Klimontovich distribution function:

[TABLE]

where $N_{s}^{p}$ is the number of physical (point-like) particles and $\delta(x)$ is the Dirac delta function.

Direct usage of this singular distribution function is impractical for two reasons. First, the number of particles to be simulated is too large to be tractable, and thus, second, this distribution function would result in an overwhelming shot noise due to the finite number of particles used in practice (Lipatov, 2002). Both can be mitigated in simulations by using a smoothed approximation for the distribution function. Thus, we approximate the distribution function by

[TABLE]

where $N_{s}$ is the number of macro-particles (defined below), and $w=N_{s}^{p}/N_{s}$ is the number of physical particles that a macro-particle represents.

The macro-particles (also called computational particles333Both macro-particles and particles will be interchangeably used to mean the same thing: particles used in simulations. When we refer to physical particles, it will be explicitly specified.) have a “shape” $S(x,x_{i_{s}})$, i.e., a smoothed localized charge distribution that has $w$ physical particles centered at $x_{i_{s}}$, with

[TABLE]

The charge and the mass for these computational particles are $Q_{s}=wq_{s}$ and $M_{s}=wm_{s}$, respectively.

The macro-particles have the same plasma frequency $\omega_{P}$ as the physical particles they represent in simulations: they have same charge-to-mass ratio and hence $\omega_{p}^{2}=q_{s}^{2}N^{p}_{s}/(Vm_{s}\epsilon_{0})=Q_{s}^{2}N_{s}/(VM_{s}\epsilon_{0})$, where $V$ is the volume. Both macro-particles and physical particles also have the same normalized temperature $\theta_{s}=k_{B}T_{s}/M_{s}c^{2}=k_{B}T_{s}^{p}/m_{s}c^{2}$, where $T_{s}$ and $T_{s}^{p}$ are the effective temperatures of the macro-particles and physical particles, respectively. This arises from the assumption that the macro-particles are monolithic, and thus have the same velocity distribution as the underlying physical particles, and thus

[TABLE]

Inserting the phase-space distribution function in Equation (5) into the Vlasov–Maxwell system (Equations (1) and (2)), we obtain from the first two moments the following equations of motion for the macro-particle of species $s$:

[TABLE]

where $E(x)$ is the solution of Maxwell’s equations (2) with moments given by

[TABLE]

where $v_{i_{s}}=u_{i_{s}}/\gamma_{i_{s}}$. Below we explain how implicit discretization of such a system of equations is achieved in our code.

2.2 Spatial grid

For a system of macro-particles in a periodic box (line) of length $L$, we divide our domain into $N_{c}$ cells each of size $\Delta x=L/N_{c}$. Assuming $k\in\{0,1,...,N_{c}-1\}$, we define the $k$th cell ($c_{k}$) centered at $x_{k+1/2}\equiv x_{k}+\Delta x/2$ as $x\in[x_{k},x_{k+1})$, where $x_{k}\equiv\Delta x\hskip 2.84544ptk$.

We adopt spline functions extending over a number of grid cells as the shape function of the macro-particles. Therefore, the distribution of physical particles inside these macro-particles is symmetric around their center and extends over a number of computational cells depending on the order $m$ of the spline functions used. For instance, $m=1$ is a top-hat distribution (shape) given by

[TABLE]

We also define, $t^{n}=n\times\Delta t$, $E^{n}\equiv E(t^{n})$, i.e., superscript $n$ denotes the $n$th time step for the particular quantity. We choose the time step $\Delta t$ such that $c\Delta t\leq\Delta x$, to obey the Courant–Friedrichs–Lewy (CFL) stability condition in 1D (Courant et al., 1967).

2.3 Charge and current deposition

To obtain a discrete set of equations that governs the evolution of such macro-particles, we begin by integrating the first equation in (2) over the $k$th cell $c_{k}$:

[TABLE]

Here, $\rho^{n}_{k+1/2}$ is the average charge density inside $c_{k}$ at $t=t^{n}$:

[TABLE]

where

[TABLE]

defines the weighed contribution of a macro-particle at $x_{i_{s}}$ to the average charge density of the $k$th cell.

The explicit forms for both shape $S^{m}$ and weight $W^{m}$ functions that we use in our code are given in Appendix B. In Figure 2, we plot different weight functions.

2.4 Solving Maxwell’s equations

Equation (13) gives the change in the electric field at cell edges exactly. However, a complete solution also requires the boundary condition, $E_{k=0}$. Therefore, we rewrite Equation (13) as

[TABLE]

We see that $E^{n}_{k}$ inherits the error of $E^{n}_{0}$. To find $E^{n}_{0}$, we first find the sum of cell-edges field $E^{n}_{\rm tot}\equiv\sum_{k}E^{n}_{k}$. The second equation in (2) can be re-written (for $m>0$) as

[TABLE]

The first equality uses $\sum_{k}\Delta xS^{m}(x_{k},x_{i})=\sum_{k}W^{m-1}(x_{k},x_{i})=1$, which is a property of spline functions implemented in our code (see Table B). Note that $j_{\rm tot}$ can be calculated exactly at each time-step using the macro-particles’ velocities. If the plasma macro-particles have a total current, then $\partial_{t}E\neq 0$, i.e., setting $\partial_{t}E$ to zero will artificially add a constant electric field on the grid or equivalently a counter current.

Equation (17) is approximated to second, third, and fourth orders of accuracy, respectively, as

[TABLE]

Proceeding in a similar way to generate higher-order accurate, asymmetric estimates for $E^{n}_{\text{tot}}$ results in numerically unstable approximations. Thus, we stop at fourth-order accurate method given in Equation (20).

To find $E_{0}$ from $E_{\rm tot}$, we multiply Equation (13) by the index $k$, and then sum over all cells. The left-hand side is given by

[TABLE]

Then, using the periodicity of $E_{k}$, i.e., $E_{N_{c}}=E_{0}$, we can write

[TABLE]

Thus, for a given set of particle data $\{x^{n}_{i_{s}},v^{n-1/2}_{i_{s}}\}$, we are able to find the electric field at the edges of the cells, $E^{n}_{k}$. Importantly, the only source of error is the error introduced in finding $E^{n}_{\rm tot}$.

The error, for fixed $\Delta t/\Delta x$ (motivated by the CFL condition), is of order $O(\Delta x^{4})$ if we use Equation (18) to update $E_{\rm tot}^{n}$. The error drops to $O(\Delta x^{5})$ or $O(\Delta x^{6})$ with Equations (19) or (20), respectively. It is also important here to note that an $O(\Delta x^{k})$ error in the force would introduce an $O(\Delta x^{k+1})$ error in the updated data of the macro-particles.

2.5 Back-interpolation: Force on Macro-particles

Having determined the electric field using Equations (13), (14), and (18) or (19) or (20) to different orders of accuracy from the macro-particle positions and velocities, we now calculate the force on the individual macro-particles. The effective electric field that acts on the macro-particle (Equation (9)) can be determined from the electric field at each cell face, $E_{k}^{n}$ as follows

[TABLE]

In Appendix A.2, we show that the approximation in Equation (23) leads to exact momentum conservation and vanishing self-forces. Another possible approximation of $E^{n}_{i_{s}}$, that has the same order of accuracy is given by

[TABLE]

The approximation in Equation (24) is used, for instance, in Haugbølle & Frederiksen (2013). It generally leads to a violation of momentum conservation and unphysical self-forces as shown in Appendix A.1. The order of accuracy for the back-interpolation step is typically second order (e.g., Lapenta & Markidis, 2011; Haugbølle & Frederiksen, 2013; Brackbill, 2016; Lapenta, 2016) because of the use of lower order interpolation functions.

Using higher-order interpolation functions implies a smoother representation of the phase-space distribution function. As a consequence, we can assume a smoother representation of the electric fields444An order $m$ interpolation function means that $f$ is $m+1$ times spatially differentiable.. Therefore, since we implemented up to fifth-order spline interpolation function, this allows constructing up to fifth-order accurate back-interpolation. To derive a higher-order accurate method, we approximate the integration $I_{k}=\int_{c_{k}}E^{n}(x)S(x,x^{n}_{i_{s}})dx$ using Simpson’s rule:

[TABLE]

where we define $S^{m}_{k,i_{s}}\equiv S^{m}(x_{k},x^{n}_{i_{s}})$.

To utilize Equation (25), we need to approximate $E^{n}_{k+\frac{1}{2}}$ in terms of cell-edge fields $E^{n}_{k}$ as follows

[TABLE]

Using Equation (25) and since $\Delta xS^{m}\in[0,1]$ (see Table B), $O(\Delta x^{4})$ error in $E_{k+\frac{1}{2}}$ implies $O(\Delta x^{4})$ error order in $E_{i_{s}}$ and $O(\Delta x^{6})$ error in $E_{k+\frac{1}{2}}$ implies $O(\Delta x^{5})$ error order in $E_{i_{s}}$ (because of the error order in Equation (25)). Therefore, the approximate electric field $E^{n}_{i_{s}}=\sum_{k}I_{k}$ on a macro-particle at $x_{i_{s}}$, using the periodicity of $E_{k}$, can be expressed as

[TABLE]

To summarize, here we showed how we can find the forces on the individual macro-particles for a given discretized field on a grid. The force error can be of order $O(\Delta x^{2})$, $O(\Delta x^{4})$, or $O(\Delta x^{5})$ if we use Equation (23), (29), or (29), respectively. Given the numerical error in finding $E_{k}^{n}$, so far we have shown how, for a given set of particle data $\{x^{n}_{i_{s}},v^{n-1/2}_{i_{s}}\}$, we can find the forces on such particles $F^{n}_{i_{s}}$. The error can be of order $O(\Delta x^{2})$, $O(\Delta x^{4})$, or $O(\Delta x^{5})$, by employing consecutively higher-order equations in finding $E^{n}_{0}$.

2.6 Pusher: particle update

To push the individual particles, we use a leapfrog scheme to discretize the equations of motions for the particles,

[TABLE]

The code assumes that the initial positions of the particles are provided at $t=0$, but the initial momenta are given at $t=-\Delta t/2$. It also assumes that at $t=-\Delta t$ the sum of electric field at cell edges was zero, i.e., $E_{\rm tot}^{n=-1}=0$.

2.7 Error sources in our algorithm

Here, we summarize the different sources of error in the algorithm presented above.

The use of discretized equations to compute $E_{\rm tot}$, yielding up to $O(\Delta t^{5})$ accurate schemes (Equations (18)-(20)). 2. 2.

The interpolation of the field from the grid to the macro-particle to calculate the force $F_{i_{s}}$. This can be done at $O(\Delta x)$ in such a fashion that total momentum is conserved exactly, or at higher accuracy, $O(\Delta x^{3})$ and $O(\Delta x^{4})$ at the cost of (slightly) violating momentum conservation. 3. 3.

The updating of the particle positions, which is currently performed at $O(\Delta t^{3})$.

For fixed $\Delta t/\Delta x$, as implied by the CFL condition, the above imply that our code is fundamentally second-order accurate, limited by the particle pusher. That is, despite improving energy and momentum conservation, the order of the interpolation function does not set the order of accuracy of the overall scheme. Implementing a higher-order symplectic integrator within our scheme would improve the convergence order; we leave this point for future work. Nevertheless, as we will show in Sections 4-6, the improved spatial order produces substantial practical enhancements in the code performance.

2.8 Normalized equations

Using the fiducial units ($c,n_{0},q_{0},m_{0}$), we define the following scales

[TABLE]

We then define our dimensionless variables as $\bar{t}=t/t_{0}$, $dt=\Delta t/t_{0}$, $\bar{x}=x/x_{0}$, $h=\Delta x/x_{0}$, $\bar{u}=u/c$, $\bar{Q}_{s}=Q_{s}/q_{0}$, $\bar{M}_{s}=M_{s}/m_{0}$, $\bar{\rho}=\rho/\rho_{0}$, $\bar{j}=j/j_{0}$ and $\bar{E}=E/\mathbb{E}_{0}$. We also identify $t_{0}$ as the time scale of the plasma frequency of the entire plasma that includes contributions from all species. This defines $n_{0}$ as follows

[TABLE]

where $n_{s}$ is the number density of the macro-particles of species $s$ and

[TABLE]

In the case where all species have the same mass and charge, $n_{0}\Delta x=(\sum_{s}N_{s})/N_{c}=N_{\rm t}/N_{c}$. In such a case, we chose our fiducial units so that $\bar{Q}_{s}^{2}=1$, $\bar{M}_{s}=1$, therefore, $n_{0}=\sum_{s}n_{s}$.

In terms of the dimensionless variables defined above, the equations that are solved by the code can be written as follows

[TABLE]

The general loop in the code is then

[TABLE]

A schematic representation of this loop is shown in Figure 1.

2.9 Implementation

SHARP-1D is implemented in C++, and is massively parallelized using MPI. The parallelization is done by distributing macro-particles on different processors, while reserving the first processor (rank = 0) to compute the electric field on the grid and to manage outputs.

3 Conserved quantities in PIC

As we mentioned in the introduction, a major test of the accuracy and fidelity of PIC schemes is their ability to preserve conserved quantities such as energy, momentum, and charge. Below we discuss such conservation laws in PIC schemes and how well they are respected when different methodologies are used. We begin by showing that our algorithm is charge conserving. Then, we compare the energy and momentum conservation properties when implicit and explicit techniques are employed555Typically discretization in PIC is done using either explicit or implicit schemes. In explicit schemes (such as the algorithm presented here), particle data are used first to calculate the fields on the grid and then particle advancement in time is carried out using these fields. On the other hand, when implicit discretization is employed, the equations for fields on the grid and evolution equations of particles have to be solved simultaneously in order to evolve forward in time. . We also show that the usage of higher-order interpolation leads to a decrease in the aliasing, which improves energy conservation while maintaining exact momentum conservation.

3.1 Charge conservation

In traditional implementations of PIC, the direct interpolation from particles to grid points, in order to calculate the grid charge and current densities, leads to violation of the continuity equation on the grid. The reason for that is calculating the current density requires the knowledge of both particles positions and velocities at the half-time step and approximating the particles positions at half-time step leads to an error in the calculated current density when particles cross cell-boundaries.

Recently, several methods were proposed where the calculation of the current density on the grid from the particles is done such that the continuity equation is satisfied on the grid at all times (Eastwood, 1991; Villasenor & Buneman, 1992; Esirkepov, 2001; Umeda et al., 2003). In the presented algorithm, we locally obey the discretized continuity equation at all times (i.e., we use a charge conserving scheme). The discretized current density666In fact, our presented algorithm does not require to calculate the current density for solving Poisson’s equation. coincides with the current density proposed in Esirkepov (2001).

To see this, we integrate the continuity equation $\partial_{t}\rho(x,t)+\partial_{x}J_{x}(x,t)=0$ over a cell of size $\Delta x$,

[TABLE]

Therefore, we can write

[TABLE]

Using the second equation in (2), the current density at cell edges can be expressed as follows

[TABLE]

Then using Equations (42) and (13), we can write

[TABLE]

Therefore, our scheme obeys exactly the second-order accurate continuity equation.

3.2 Energy and momentum conservation

In general, for PIC schemes, energy and momentum conservation appear to be mutually exclusive (see, e.g., Brackbill, 2016). Momentum non-conservation typically comes from non-vanishing self-forces and errors in the interaction forces – an example of such a case is shown in Appendix A.1. These non-physical forces can produce macroscopic non-physical instabilities (Langdon, 1973). Energy non-conservation can also produce dramatic changes in the evolution of the plasmas. Since such an error has a secular (unbounded) growth, energy non-conservation imposes a serious limitation on the ability to study the nonlinear phenomena which occurs in long time scales (compared to $\omega_{p}^{-1}$).

For instance, Lapenta & Markidis (2011) demonstrated that in the two-stream instability, the errors in the energy are disproportionately distributed to the fast particles. In particular, they demonstrated that even though the per-particle violation in energy conservation is small, the disproportionate distribution of energy non-conservation leads to errors in the distribution of particles, especially at the high energy end. This is important for particle acceleration in relativistic situations (Lapenta & Markidis, 2011). Also, results that probe the long-term behavior of particle distribution functions starting from linearly unstable conditions such as tenuous beam instabilities or particle acceleration are subject to these issues (e.g. Sironi & Giannios, 2014; Ardaneh et al., 2015; Park et al., 2015).

Traditional explicit algorithms lead to numerical increase in the total energy (numerical heating) while conserving the total momentum exactly (Birdsall & Langdon, 1991; Hockney & Eastwood, 1988). For explicit schemes, an energy conserving algorithm was developed in Lewis (1970). However, the total energy is conserved only in the limit of $\Delta t\rightarrow 0$. In practice, there will be numerical heating of the plasma because of the finite timestep.

On the other hand, traditional implicit algorithms tend to decrease the total energy numerically (numerical cooling), while violating the total momentum conservation (Brackbill & J. Forslund, 1985). Recently, implicit algorithms that, in principle, conserve the total energy exactly, while still violating momentum conservation, were introduced: for non-relativistic/classical plasmas by Markidis & Lapenta (2011) and Chen et al. (2011) and relativistic plasmas by Lapenta & Markidis (2011). However, in practice, these algorithms use the Jacobian-Free Newton Krylov (JFNK) method to solve the full implicit system, introducing an error that depends on the accuracy of the Newton or Picard iteration. This leads to violation in energy conservation that can be controlled by increasing the accuracy of such methods.

One major source of energy non-conservation is the coupling between the wave modes resolved by the grid and their aliases. Aliases are wave modes that differ by an integer number of $2\pi/\Delta x$. The reason for such coupling is that continuous particle data (which support wave modes that include aliases of wave modes resolved by the grid) are used in the construction of phase-space moments at a discrete set control of points on the physical grid. Therefore, for momentum conserving schemes, we improve the energy conservation by decreasing the coupling of the wave modes resolved by the grid with their aliases.

One way to decrease the effect of aliasing is done by filtering the deposited grid moments (charge and current densities), this is used, for instance, in TRISTAN (Buneman et al., 1993) and its parallel version TRISTAN-MP (Spitkovsky, 2005). Filtering, however, when used with a momentum conserving scheme results in violating the momentum conservation and non-vanishing self-forces (an example for such case is presented in Appendix A.3).

Alternatively, energy conservation is naturally improved when using higher-order interpolation functions, as employed here. In Fourier space, using higher-order interpolation function is qualitatively equivalent to low-pass filtering. However, for instance, filtering when used with $W^{0}$ produces larger energy errors in comparison to the errors produced when first-order interpolation, $W^{1}$, is used (cf. Section 8.7 of Birdsall & Langdon 1991). That is, higher-order shape functions considerably decrease the impact of aliasing, resulting in an improved energy conservation for the same reasons as filtering, while simultaneously maintaining higher accuracy evolution of the underlying system. This includes, potentially, exact momentum conservation. As shown explicitly in Appendix C, the discretization of the plasma into macro-particles with shape functions of order $m$ produces a spectral smearing with a width that scales as $(k\Delta x)^{-m}$, producing a corresponding an exponential decrease in aliasing for the resolved modes with shape-function order.

4 Validation

We begin the assessment of the numerical algorithm presented in Section 2 with the physical validation of simulation results against known results. This represents the first in two critical numerical tests, the second being convergence, which is treated in Section 6.

Specifically, we present comparisons with the analytical or semi-analytical results on the following test problems:

thermal stability of a uniform plasma (Section 4.1), 2. 2.

standing plasma waves – Plasma oscillations and linear Landau damping (Section 4.2), and 3. 3.

two-stream instabilities – non-relativistic and relativistic (Section 4.3).

Because validation essentially consists of quantitative comparisons with known results, we are currently limited to phenomena in the linear regime. This appears to be a wide-spread difficulty within plasma simulations. Nevertheless, weak validation in the nonlinear regime can be found when we discuss code comparisons in Section 5. Here, we assume that all plasma species have the same mass, and, up to a sign, the same charge. Thus, we take as our fiducial units, $q^{2}_{0}=Q^{2}_{s}$ and $m_{0}=M_{s}$, which implies that $n_{0}=\sum_{s}n_{s}$.

4.1 Thermal stability of plasma

In the absence of microscopic radiative processes (e.g., Compton scattering, bremsstrahlung), a uniform, thermal plasma should not evolve. In practice, even a uniform thermal plasma will numerically heat. This arises mainly as a consequence of aliasing as explained in Section 3.2.

There are two temperature scales that are often relevant for the numerical heating of plasma simulated with variants of the PIC algorithm. The first is the numerical Debye temperature,

[TABLE]

which is the temperature at which the Debye length is equal to square of the cell size in code units. Note that this is a purely numerical quantity that defines those temperatures below which the discretization of Maxwell’s equations no longer resolves the Debye length. Typically, for $\theta<\theta_{D}$ second-order accurate codes will exhibit virulent numerical heating until $\theta\approx\theta_{D}$ (see, e.g., Birdsall & Maron (1980)). Thus, $\theta_{D}$ presents a key numerical limitation on the classes of plasmas that have been simulated to date. However, as shown in Section 3.2, implementing higher-order spatial interpolation decreases the effect of aliasing and hence considerably decrease such heating as shown in Figure 3.

The second is the Poisson temperature, set by the Poisson fluctuations in the reconstruction of the particle distribution. This is set by equating the average potential and kinetic energies of randomly distributed particles (Appendix D.3), though in the cold and hot limits (i.e., non-relativistic and relativistic velocity dispersion limits, respectively) this reduces to

[TABLE]

Temperatures below $\theta_{P}$ are not well defined numerically. Note that the ordering of $\theta_{P}$ and $\theta_{D}$ is not fixed, though we will consider cases when $\theta_{P}<\theta_{D}$ exclusively.

Here, we study the heating due to different approximations in PIC algorithms. To this end, we start all of our simulations with temperatures higher than $\theta_{P}$ and study the evolution of the plasma temperatures. We perform a series of simulations for two populations of negatively and positively charged macro-particles in a periodic box of length $\bar{L}=5$ and cell size of $h=0.1$. The total number of macro-particles $N_{\rm t}=2\times 10^{5}$. Therefore, $\theta_{P}\approx 1.04\times 10^{-5}$ and $\theta_{D}=10^{-2}$.

Each simulation is started at a different initial temperature: $\theta=10^{-1},10^{-2},10^{-3},10^{-4}$. In all simulations presented here, we use second-order accurate back-interpolation, i.e., Equation (23).

Since for all simulations here, $\theta\ll 1$, we start with an initial distribution function that is given by

[TABLE]

For the various initial temperatures, we show in Figure 3 the evolution of the plasma temperatures (top panel) and that of the fractional energy error of the plasma (bottom panel). We show results when first-order interpolation $W^{1}$ (left), third-order $W^{3}$ (middle), and fifth-order $W^{5}$ (right) are used in deposition and back-interpolation steps.

As expected, for the first-order interpolation (the scheme that is most commonly employed in existing PIC codes), uncontrolled heating is observed for all $\theta\leq\theta_{D}$. This heating subsides when $\theta\approx 3\theta_{D}$, requiring between one and two cells per Debye length. Hence, first order PIC algorithms face severe computational requirements to resolve cold plasmas.

However, using higher-order spatial interpolation significantly reduces the temperature at which uncontrolled numerical heating occurs. By fifth order, temperatures four orders of magnitude smaller than $\theta_{D}$, and only an order of magnitude larger than $\theta_{P}$, can be resolved for millions of plasma timescales. That is, high-order spatial interpolation extends the range of temperatures and timescales that can be simulated.

The marked improvement of SHARP-1D is a direct result of the corresponding improvement in energy conservation. The bottom panels of Figure 3 show the evolution of the growth in the energy of each simulation; in all cases, the unphysical heating can be fully attributed to the failure to conserve energy. However, the fractional energy non-conservation is improved by nearly three orders of magnitude for each simulated temperature as the spatial interpolation order is increased from $W^{1}$ to $W^{5}$. The net result of higher spatial order is, therefore, the ability to run simulations orders of magnitude longer with orders of magnitude lower resolutions.

4.2 Stability of standing linear plasma waves

We now turn to the stability and evolution of plasmas with linear perturbations, specifically, standing waves. Key validation tests are the reproduction of oscillation frequencies, dispersion relations, and linear Landau damping rates. We begin with a discussion of the anticipated values followed by quantitative comparisons of standing wave evolution.

4.2.1 Linear Dispersion Relations and Growth Rates

As shown in Appendix E, the linear dispersion relation for thermal plasmas are conveniently expressed in terms of a handful of dimensionless quantities:

[TABLE]

[TABLE]

[TABLE]

where $v_{\rm th}=\theta^{1/2}c$ is the thermal velocity dispersion and $k_{D}$ is the wavenumber associated with the Debye length. The linear dispersion relation for a non-relativistic thermal population of uniformly distributed electrons with a fixed positively charged background, i.e., infinitely massive ions, is given by (see Appendix E for more details)

[TABLE]

where Erfi is the complex error function defined as $\text{Erfi}(v_{p})=-i\text{ Erf}(iv_{p})$. For a given $\hat{k}$, the roots $\hat{\omega}_{j}=\omega_{j}/\omega_{p}$ can be found by solving Equation (50) numerically; the real and imaginary values for $\omega_{j}$ yield the oscillating frequency and the growing/damping rate of the mode with wavenumber $k$, respectively, at a given thermal velocity $v_{\rm th}$.

Approximate expressions for the roots of Equation (50) are often obtained for the limit $\hat{k}\ll 1$. The most common of these, and the standard expression found in most textbooks (e.g., Boyd & Sanderson, 2003, and referred to here as “Standard”) is

[TABLE]

where $\omega_{r}$ and $\omega_{i}$ are the real and imaginary components of $\omega$, respectively. A more accurate expression, derived by a higher-order approximation to Equation (50) is given by McKinstrie et al. (1999, referred to as “Extended”).

[TABLE]

Both of these approximates are shown in comparison to the full numerical solution (“Numerical”) in Figure 4. It is immediately evident that the regime of applicability of both the Standard and Extended approximations is limited to $\hat{k}<0.25$, with the numerical consequence that neither are quantitatively accurate in the rapid damping regime, where numerical validation experiments can most easily be performed.

For convenience, we provide below a numerical fitting formula based on the formulation of McKinstrie et al. (1999) to the full numerical solution for $\hat{k}\in[0,0.6]$ that is good to 4% throughout and better than 0.5% above $\hat{k}=0.3$

[TABLE]

In the following subsections, we report a series of simulations to test the code in the regimes $\omega_{i}\sim 0$ (undamped modes) and $\omega_{i}\neq 0$ (damped modes). We use a fixed neutralizing background and negatively charged plasma macro-particles whose initial distribution function is given by

[TABLE]

In all simulations, we use $\theta=10^{-3}$, $\alpha=10^{-2}$, and fifth-order spatial interpolation, i.e., $W^{5}$. The damping rates and oscillation frequencies in different simulations of this section are also shown in Figure 4.

{deluxetable*}

ccccc ccc ccc cc

\tablewidth18.6cm \tablecaptionStanding plasma waves simulation parameters and results.

\tablehead $\hat{k}$ & $N_{\rm t}$ $\bar{\lambda}$ \tablenotemarka $\bar{\lambda}/h$ \tablenotemarkb $\bar{L}/\bar{\lambda}$ $\omega_{r}/\omega_{p}$ \tablenotemarkc
$\omega^{sim}_{r}/\omega_{p}$ \tablenotemarkd
$R_{r}$ \tablenotemarke
$\omega_{i}/\omega_{p}$ \tablenotemarkc
$\omega^{sim}_{i}/\omega_{p}$ \tablenotemarkd
$R_{i}$\tablenotemarke
$\omega_{p}\hskip 0.56917ptt_{\rm min}$ \tablenotemarkf $\omega_{p}\hskip 0.56917ptt_{\rm max}$ \tablenotemarkf

\startdata $0.01$ $10^{9}$ $19.869$ $1761$ $2$ $1.00015$

$0.99925\pm 3\times 10^{-7}$ $0.9991$    –       –    – [math] $10^{3}$

$0.35$ $5\times 10^{8}$ $0.568$ $67$ $12$ $1.22095$ $1.21989\pm 2\times 10^{-5}$ $0.9991$ $-0.03432$ $-0.0344\pm 4\times 10^{-5}$ $1.0036$ $1.7$ $30.84$

$0.40$ $5\times 10^{8}$ $0.497$ $62$ $15$ $1.28506$ $1.28145\pm 4\times 10^{-5}$ $0.9972$ $-0.06613$ $-0.0687\pm 8\times 10^{-5}$ 1.0389$$1.7 $16.80$

$0.45$ $5\times 10^{8}$ $0.442$ $68$ $19$ $1.35025$ $1.34617\pm 1\times 10^{-5}$ $0.9969$ $-0.10629$ $-0.1066\pm 2\times 10^{-5}$ 1.00326$$1.6 $16.65$

$0.50$ $5\times 10^{8}$ $0.397$ $68$ $19$ $1.41566$ $1.40786\pm 2\times 10^{-5}$

$0.9945$ $-0.15336$ $-0.1546\pm 6\times 10^{-5}$ 1.0081$$1.6 $8.77$ \enddata\tablenotetextaWavelength of the initially excited mode in code units. \tablenotetextbNumber of cells used to resolve the initially excited wavelength. \tablenotetextcNumerical solution of Equation (50). \tablenotetextdOscillation frequencies and damping rates found by fitting the changes in simulations. The error in the fitting parameters corresponds to 99% confidence level.

\tablenotetexte$R_{r}=\omega^{\text{sim}}_{r}/\omega_{r}$ and $R_{i}=\omega^{\text{sim}}_{i}/\omega_{i}$, where $\omega^{\text{sim}}_{r}$ and $\omega^{\text{sim}}_{i}$ are the oscillation frequency and damping rate obtained from fitting simulation results, respectively. An example for such a fit is shown in Figure 5 for the case of $\hat{k}=0.01$ and in Figure 7 for the case of $\hat{k}=0.45$ . \tablenotetextfTo fit our simulation results, we specify a time range [$t_{\rm min},t_{\rm max}$] in the simulation over which we carry out the fitting.

4.2.2 Plasma oscillations

At small $\hat{k}$, the linear Landau damping rate is vanishingly small, i.e., $\omega_{i}\approx 0$. As a result, a linear perturbation should oscillate providing quantitative tests in the form of the oscillation frequency and evolution of the mode amplitude. We initialize the simulation with an excited mode such that $\hat{k}=0.01$. The theoretical predictions (numerical solution of Equation (50)) are $\omega_{r}/\omega_{p}=1.00015$ and $\omega_{i}=-4.7\times 10^{-2167}$. The rest of our simulation parameters are given in Table 4.

Figure 5 shows the fitted values for the electric field of the initially excited oscillation modes. The oscillation frequency of the initially excited mode is found to be within 0.09% of the theoretically predicted oscillation frequency (this is measured by fitting the oscillation of such a mode in the simulation over about 159 oscillation periods). The oscillation frequency and amplitude are still in excellent agreement with the theoretical prediction until the end of the simulation time $t=10^{3}\omega_{p}^{-1}.$

In Figure 6, we show the evolution of the averaged energy in the excited modes (over five plasma periods). We find that most of the power stays in the excited mode for the whole simulation period. Coupled with the degree of energy conservation during the simulation, this implies that no more that 0.1% of the initial energy in the mode leaks into other degrees of freedom of the plasma (e.g., other modes or heating of the plasma).

4.2.3 Linear Landau damping rates

For larger $\hat{k}$, corresponding to comparatively larger wavenumbers, the damping rates become large. By $\hat{k}=0.35$ the wave should damp by one $e$-fold in about six wave oscillation periods. Again, this provides a number of quantitative tests of SHARP-1D: oscillation frequencies and damping rates. Thus, here we report on simulations at large $\hat{k}$; the values of parameters in these simulations and their results are summarized in Table 4. The damping rates and oscillation frequencies in different simulations are also shown in Figure 4.

In all simulations, the oscillation frequencies and damping rates are found to be within 0.5% and 0.8% of the theoretical predictions of linear theory, respectively. The two columns $R_{r}$ and $R_{i}$ of Table 4 report the ratio between the oscillation frequency and damping rate in simulations to their theoretically predicted values, obtained by numerically solving Equation (50). The damping rate and the oscillation frequency of the simulation are obtained by fitting the evolution of the Fourier component of the electric field that corresponds to the initially excited wave mode. In Figure 7, we show an example (for $\hat{k}=0.45$ simulation) of the fitting carried out to find oscillation frequencies and damping rates of different simulations.

A sudden drop in the mode energy at the first period can be observed in Figure 7. This is a characteristic feature of linear Landau simulations, i.e., when damping rate is comparable to the plasma frequency. Such a feature is also present when other simulation methods are used to simulate the evolution of such modes in a thermal plasma (e.g., Besse & Sonnendrücker, 2003; Rossmanith & Seal, 2011).

4.3 Two-stream Instability

We now consider the quantitative accuracy with which SHARP-1D can reproduce a dynamical instability – in 1D the primary example is the two-stream instability. This provides an opportunity to also assess the relativistic performance of the code, through the simulation of relativistic beams. Thus, we will consider two limiting regimes: non-relativistic ($v_{b}\ll c$) and relativistic ($u_{b}\gg c$), where $v_{b}$ is the speed of the streams and $u_{b}=v_{b}/\sqrt{1-(v_{b}/c)^{2}}$ is the specific momentum.

As with linear Landau damping, we will begin with a general discussion of the anticipated instability properties and then move onto quantitative comparisons.

4.3.1 Instability Growth Rates

In the non-relativistic regime, i.e., the streams are moving with non-relativistic speeds and have non-relativistic temperatures, the linear dispersion relation for two counter streams in such limit is given by (see Appendix E for more details)

[TABLE]

where $z_{b}=\bar{v}_{b}/\sqrt{\theta}$, $\bar{v}_{b}=v_{b}/c$, $\hat{k}=kc\sqrt{\theta}/\omega_{p}$ and $v_{p}=\omega/kc\sqrt{\theta}$. In general, this must be numerically solved to obtain the mode frequencies. The solution generally consist of two oscillatory modes, a growing mode, and several damping modes. We present some of these solutions in Table 10. In Figure 8, we show the numerical solutions for the maximum growth rates (and the mode growing with such a rate) as a function of the stream speed and its temperature.

In the relativistic limit, the beam velocity distribution exhibits a narrow peak very close to the speed of light $c$, and is thus well described by the cold-plasma limit, i.e., $\theta=0$. Within this limit, the linear dispersion relation for two relativistic counter streaming $e^{+}$-$e^{-}$ populations (with speed $v_{b}$) is given by

[TABLE]

In this case, it is possible to obtain analytic solutions:

[TABLE]

where $\hat{k}_{b}=kv_{b}/\omega_{p}$. When $\gamma_{b}^{3}\hat{k}_{b}^{2}<1$, these again correspond to two oscillating modes, a growing mode and a damping mode. The positive imaginary root is maximized at the wavenumber

[TABLE]

at which the growth rate is

[TABLE]

where $\Gamma_{m}=\Im(\omega_{m})$ is the growth rate of that wavenumber.

The anticipated growth rates for both the relativistic and non-relativistic regimes are shown in Figure 9. At low beam velocities the finite temperatures of the beams suppress the growth rates appreciably relative to the cold-plasma limit, highlighting the importance of numerically solving the dispersion relation in the non-relativistic regime. In both cases (non-relativistic and relativistic streams) investigated below, we use a fifth-order interpolation function $W^{5}$.

4.3.2 Non-relativistic Two-stream Simulations

The non-relativistic two-stream simulations are initialized with an initial co-moving temperature of the streams at $\theta=10^{-4}$. We use $\Delta x\sim 0.001c/\omega_{p}$ and $L/\lambda_{m}=20$, where $\lambda_{m}=2\pi/k_{m}$ is the fastest growing wavelength. We start with a uniform distribution of macro-particles. Therefore, the initial distribution function is given by

[TABLE]

Other simulation parameters and the theoretical prediction are given in Table 10.

Generally, these simulations exhibit excellent quantitative agreement with the results of the linear theory. This is clearly evident in Figure 10, which shows the evolution of the amplitude of the most rapidly growing mode (based on linear theory) in the non-relativistic simulations as a function of linear growth times. These should be compared to the solid black line, which shows the expected exponential growth, i.e., $e^{\Gamma_{m}t}$; the correspondence lasts over 4-8 $e$-folding times, i.e., 2-3.5 orders of magnitude, ending when the instability saturates non-linearly.

{deluxetable*}

ccccc cccc

\tablewidth18.6cm \tabletypesize \tablecaption Non-relativistic two-stream simulation parameters.

\tablehead $v_{b}/c$ $N_{c}$ \tablenotemarka $N_{\rm pc}$ \tablenotemarkb $\Gamma_{m}/\omega_{p}$ \tablenotemarkc $k_{m}c/\omega_{p}$ \tablenotemarkd $\hat{\omega}_{m}^{1}$ and $\hat{\omega}_{m}^{2}$ \tablenotemarke
$\hat{\omega}_{m}^{3}$ \tablenotemarke
$\hat{\omega}_{m}^{4}$ \tablenotemarke
$\hat{\omega}_{m}^{5}$ \tablenotemarke

\startdata 0.02 5046 900. 0.168555 24.9 $\pm 1.40842-0.0161625i$ $-0.485447i$ $-0.915823i$ $-1.32416i$

0.025 5324 900. 0.243771 23.6 $\pm 1.47658-0.0123756i$ $-0.342545i$ $-0.673817i$ $-0.987125i$

0.03 6013 899.85 0.284945 20.9 $\pm 1.48041-0.00479742i$ $-0.249662i$ $-0.48542i$ $-0.715823i$

0.035 6904 900. 0.307734 18.2 $\pm 1.46288-0.00103408i$ $-0.213058i$ $-0.356584i$ $-0.524922i$

0.04 7903 900. 0.321003 15.9 $\pm 1.44312-0.00011898i$ $-0.247524i$ $-0.391793i$ $-0.527878i$

0.05 9974 899.91 0.334645 12.6 $\pm 1.41697-3.246\times 10^{-7}i$ $-0.340762i$ $-0.388706i$ $-0.494053i$

0.06 12083 900. 0.34111 10.4 $\pm 1.40108-1.294\times 10^{-10}i$ $-0.341415i$ $-0.467936i$ $-0.496439i$

0.08 16110 900. 0.346906 7.8 $\pm 1.39159-3.816\times 10^{-19}i$ $-0.346906i$ $-0.396172i$ $-0.524619i$

0.09 18212 900. 0.348379 6.9 $\pm 1.38628-1.148\times 10^{-24}i$ $-0.348379i$ $-0.548489i$ $-0.556374i$ \enddata\tablenotetextaNumber of cells. \tablenotetextbNumber of computation particles per cell. \tablenotetextcTheoretical predictions for the maximum growth rate. \tablenotetextdFastest growing wave mode predicted theoretically. \tablenotetexte$\hat{\omega}=\omega/\omega_{p}$: other solutions of Equation (55) at the fastest growing wave mode $k_{m}$.

The initial oscillations correspond to other solutions of the dispersion relation, i.e., they are fully described by the linear analysis of the two-stream instability: in addition to the ultimately dominant exponentially growing mode, the linear dispersion relation admits other damping and oscillatory modes.

The inset in Figure 10 shows a reconstruction of the $v_{b}/c=0.02$ simulation, which shows prominent oscillations at the beginning, in terms of the full complement of plasma modes shown in Table 10. Thus, the simulation quantitatively reproduces all of the anticipated linear features. A similar exercise is possible with the remaining non-relativistic simulations as well.

4.3.3 Relativistic Two-stream Simulations

We also performed a series of simulations for streams moving with relativistic speeds $u_{b}=\gamma_{b}v_{b}$. In all simulations, the initial co-moving temperature of streams is $\theta=3\times 10^{-3}$, $\Delta x\sim 0.05c/\omega_{p}$ and $L/\lambda_{m}=10$, where $\lambda_{m}=2\pi/k_{m}$ is the fastest growing wavelength. We start with uniformly distributed macro-particles, i.e., the instability here also grows from numerical noise. The initial distribution function is given by a combination of Maxwell–Jüttner distributions:

[TABLE]

where $K_{1}$ is the Bessel function of the first kind, and $\theta$ is the temperature in the co-moving frame of each beam. All other simulations parameters, along with theoretical predictions, are shown in Table 1.

Figure 11 shows again an excellent quantitative agreement between the growth rates of the fastest growing Fourier component of the electric field in different simulations and the theoretical predictions (solid-black curve). As in the non-relativistic case, the simulation exhibits exponential growth with the anticipated growth rate over three to four $e$-folding timescales. As for the non-relativistic case, the initial oscillations can be identified with the oscillatory, non-growing modes (described by and below Equation (57)). A fit for $u_{b}/c=4.0$ that includes the oscillatory components is shown as an inset inside Figure 11.

5 Comparison with TRISTAN-MP

TRISTAN-MP is a publicly available PIC code to study plasma physics relevant for astrophysical problems (e.g., Spitkovsky, 2008; Riquelme & Spitkovsky, 2009; Philippov et al., 2015). Here, we compare some of the results we obtain from SHARP-1D with those obtained using TRISTAN-MP. In all cases, the same initial data are used. Generally, we find a substantial improvement in the ability to conserve energy and avoid numerical heating in SHARP simulations, and a good agreement in short timescale phenomena for which energy non-conservation is not substantial.

In our test problems, the performance of SHARP compares very favorably to TRISTAN-MP, typically running roughly an order of magnitude faster. We caution, however, that this may not be an entirely fair comparison since TRISTAN-MP, as a 3D code, may not be optimized for 1D problems.

5.1 Numerical Particulars of TRISTAN-MP

In addition to the initial conditions, TRISTAN-MP has a number of specific numerical parameters that impact its performance. While we are unable to perform an exhaustive analysis of each, we did explore the result of varying a handful of these. TRISTAN-MP is a 3D and 2D PIC code. Here, we use the 2D version with one or two cells in one of the spatial dimensions to run it in an effective 1D setup, which enables a fair comparison to SHARP-1D.

5.1.1 Filtering

TRISTAN-MP provides the ability to low-pass filter the deposited grid moments, e.g., current densities on the grid, damping high-frequency noise prior to using them to solve Maxwell’s equations. This reduces the coupling between the wave modes resolved on the grid with their aliases leading to improvements in the momentum and energy conservation of the algorithm. In TRISTAN-MP, filtering is accomplished with a three-cell stencil that generates a weighted average between the charge current density in a given cell with its neighbors. This operation may be repeated as many times as desired, smoothing the moments on progressively larger scales.

It is not a priori clear how many passes of the three-cell filtering operation are optimal in a given problem and we experimented with a number of different choices for the comparison problems presented here. We find that after a small number of filtering passes, typically three to five, the qualitative improvement is only moderate for $\theta\lesssim\theta_{D}$ and negligible for $\theta>\theta_{D}$.

This modest improvement comes with the additional computational cost, set by the addition of a substantial number of transverse grid cells required by many filtering passes. Thus, when comparing the numerical heating in SHARP-1D (when $W^{5}$ is used), we use four filtering passes with a 2D simulation box for TRISTAN-MP that is only two cells wide in the $x$-direction. All other comparisons employ only three filtering passes or less with a 2D simulation box for TRISTAN-MP that is only one cell wide in the $x$-direction.

5.1.2 Electromagnetic Mode Speeds

To suppress the numerical Cerenkov instability in TRISTAN-MP it is possible to independently set the ratio of the propagation speed of transverse electromagnetic modes to the speed of light. This is implemented explicitly via an additional numerical coefficient in Maxwell’s equations. Typically, this is set near to unity, e.g., 1.025. However, in the 1D electrostatic case, the numerical Cerenkov instability does not exist, and we have verified that this factor does not qualitatively change any of the results from the TRISTAN-MP simulations.

5.2 Thermal stability and energy conservation

In Figure 12 we compare the evolution of the temperature and energy error in a pair of simulations described in Section 4.1. These are chosen such that in one case the Debye length is resolved by the grid cell (red curves) and when it is not resolved (blue curves). In both cases, TRISTAN-MP (dashed lines) exhibits a significantly larger violation of energy conservation, differing only in the timescale over which this occurs. When the Debye length is resolved, the numerical heating occurs more slowly, becoming untenable only after $3.5\times 10^{6}$ plasma timescales. On the other hand, when the Debye length is not resolved, the numerical heating dominates the initial thermal energy almost instantly. In both cases by $4\times 10^{6}\omega_{p}^{-1}$, both simulations have generated similar relative degrees of numerical heating, i.e., the ratio of the energy errors to the original thermal energy of the plasma.

In comparison, SHARP-1D (solid lines) reduces the numerical heating rate drastically. When the Debye length is resolved (high temperature) the factional errors are fixed near $10^{-5}$ throughout the simulation. Lower temperature plasmas exhibit similar absolute heating rates, and therefore the relative heating for cold plasmas appears larger. However, even when the Debye length is unresolved by an order of magnitude, the plasma continues to be well modeled. The origin of the improvement in the numerical heating is the improvement in the order of interpolation.

5.3 Stability of standing linear plasma waves

Here, we compare the evolution of a standing plasma wave, where the linear Landau damping can be ignored, i.e., the ability of both TRISTAN-MP and SHARP-1D to maintain a small amplitude oscillating wave mode.

The simulation setup is similar to that in Section 4.2: a fixed uniform background of ions with thermal electrons in a box with size $L=40\hskip 1.42271ptc/\omega_{p}$. The electrons are initially uniformally distributed. The initially excited mode is added through a velocity perturbation: we first initialize electron velocities using Equation (46) with $\theta=10^{-3}$, then add a position dependent velocity perturbation to individual particles’ velocity by adding $\beta\cos(2\pi x/\lambda)$ to their velocities, where the initially excited wavelength $\lambda=20\hskip 1.42271ptc/\omega_{p}$, with $\beta=0.01\lambda/2\pi$, i.e., after about $0.25$ of a plasma period this will introduce a density perturbation with an amplitude of 0.01.777Unlike the simulations in Section 4.2, here we initialize a perturbation in the velocity, which is more easily done with the native initialization routines of TRISTAN-MP. We have verified that there are no significant differences when the mode is initialized as a density perturbation, and hence do so when we later compare a Landau damped wave. The cell size is $\Delta x=0.01\hskip 1.42271ptc/\omega_{p}$ and since $kv_{\rm th}/\omega_{p}=0.009934$, linear Landau damping can be ignored. In all simulations described here, we fix the number of electrons per cell to $N_{pc}=1250$. We also note that in all simulations the Debye length is well resolved, i.e., $\lambda_{D}=3.162\hskip 1.42271pt\Delta x$.

In Figure 13 we show the long-term evolution of the energy in this isolated wave mode. In both the SHARP-1D (solid curves) and TRISTAN-MP (dashed curves) simulations, the square of the amplitude ($|\tilde{E}_{k}|^{2}$) exhibits very small, long-timescale oscillations. In the SHARP-1D simulations, these are confined to within 0.8% of the initial value over the entire simulation. In contrast, the TRISTAN-MP simulations also exhibit a secular growth in the mode amplitude, leading to an approximate energy increase of 10% by $8\times 10^{4}\omega^{-1}_{p}$.

This behavior is independent of the interpolation order of SHARP-1D or number of smoothing filters employed in TRISTAN-MP. Even when employing $W^{1}$ with SHARP-1D, the mode amplitude continues to execute only small oscillations about the fixed value, accurately reproducing the expectation for the linear evolution of the mode. Because the Debye length is well resolved in this case, filtering improves the energy conservation only slightly in TRISTAN-MP, again, making little difference to the mode evolution.

The origin of the unphysical growth in the mode in the TRISTAN-MP simulations is unclear. The heating of the background is insufficient to appreciably Landau damp or excite the mode. We have run additional simulations with SHARP-1D employing a momentum non-conserving scheme in SHARP-1D with $W^{1}$, i.e., using Equation (24) (which is the back-interpolation scheme used in TRISTAN-MP), finding similar results. That is, excluding both the differences in order and the back-interpolation of the fields as the source of the secular growth in the mode energy. The only untested distinction remaining between SHARP-1D and TRISTAN-MP is the way in which the electric fields are updated.

Next, we compare the evolution of a shorter wavelength mode ($kv_{\rm th}/\omega_{p}=0.45$) that has a high linear Landau damping rate. The simulation’s setup is exactly the same as in Section 4.2.3, but with a lower number of macro-particles ($N_{\rm t}=5\times 10^{7}$). The time scale for damping is much smaller than the time needed for the energy non-conservation to affect the evolution of such modes. Thus, the result from both TRISTAN-MP (dashed blue curve) and SHARP-1D (solid red line) match exactly as shown in Figure 14. As noted before, the initial drop in the wave amplitude is also present when other simulation methods are used (e.g., Besse & Sonnendrücker, 2003; Rossmanith & Seal, 2011).

5.4 Two-stream Instability

We now compare the performance on dynamical instabilities, i.e., the two-stream instability as described in Section 4.3, for non-relativistic and relativistic streams. Here, again the timescale, on which such instability grows, is much shorter than the time scale needed for non-energy conservation to affect the evolution. Therefore, we find the same same linear evolution in both codes.

Figure 15 shows a comparison for two of the non-relativistic stream velocities reported in Section 4.3.2. As before, time is measured in $e$-foldings of the most unstable mode. In Figure 16 comparisons for two relativistic stream simulations reported in Section 4.3.3 are shown. Again the same linear evolution of the instability is found in both codes.

6 Convergence

At the end of Section 2.2, we discussed the accuracy and different errors introduced in our numerical scheme. We saw that the dominant error is of order $O(h^{3})$, arising from the order of particle pusher. Here, we assess the convergence of SHARP-1D and demonstrate that the numerical error decreases as expected. In particular, we develop a general criterion for convergence studies of PIC simulations by requiring that the ratio of the energy in the shortest wave mode to the energy in the Poisson noise of simulation to be, at least, fixed.

We then present a test case where the error in the total energy of plasmas is used as the measure of error: such a test shows that the definition for convergence motivated above leads to a decrease in the error at the expected rate and typical methods to test for convergence fail: a slower decrease in the error is observed as resolutions increase leading to a plateau in the error, where increasing the resolutions no longer lead to a decrease in the error.

6.1 “Resolution” in PIC Algorithms

In general, for a PIC-type algorithm, three notions of resolution are relevant for simulations.

Spatial resolution of the grid, i.e., $h$. This also determines the temporal resolution. 2. 2.

Momentum resolution, set by the number of particles used to construct the charge and current density at each cell, i.e., the number of particles per cell, $N_{\rm pc}$. 3. 3.

Spectral resolution, set by the size of the “spectral-cell,” which, for each spatial-dimension, is given by of $2\pi/L$, where $L$ is the box-size.

The third is rarely discussed in PIC simulations and arises when a physical phenomena (in linear or nonlinear regimes) has a narrow spectral support. In such cases, higher resolution simulations will require increasing the three types of resolutions simultaneously. Here we will focus on the first two, leaving a complete discussion of the third for future work.

Generally, it is necessary to increase all relevant resolutions simultaneously to study algorithmic convergence. As we will see below, this requires increasing the number of particles per cell, rather than fixing it as it is typically done. A similar requirement was found for smooth-particle hydrodynamics (SPH) simulations, where the convergence also requires increasing the number of fluid-particles within the smoothing volume of each particle to study convergence (Zhu et al., 2015).

6.2 Definition of Convergence – Equivalent Simulation

The notions of both spatial resolution (i.e., $h$) and momentum resolution (i.e., $N_{\rm pc}$) place different constraints on the range of underlying wave modes that can effectively be simulated. Thus, some care must be taken to ensure that as these resolutions are increased simultaneously, the underlying wave complement of the physical system resolves ever smaller scales.

The discrete nature of the macro-particles places a floor on the amplitude of a mode that can be effectively resolved.888This is what sets the thermal floor delineated by $\theta_{P}$ in Equation (45). The average potential energy of the particle distribution, or Poisson noise, is

[TABLE]

where recall that $f_{m}$ is a coefficient that depends on the spatial order of the algorithm. Any mode with an energy less than $\mathscr{E}_{\rm noise}^{m}$ is effectively unresolvable.

For comparison, we compute the energy in a single plasma mode that can be revolved on the grid999The wavelength $\lambda$ is resolved by the grid, if $\lambda/L$ is an integer.. For a single mode in plasma with charge density given by

[TABLE]

the electric field is given by

[TABLE]

where $A$ is the amplitude of the initial perturbation for a mode with wavelength $\lambda$. Hence, the total electric field energy is given by

[TABLE]

Averaging over a full-period and assuming that $\left\langle E^{2}_{0}\right\rangle=0$, the averaged potential energy in the wave mode is then

[TABLE]

Whether or not a mode can be resolved is then determined by the ratio of $\left\langle\mathscr{E}\right\rangle$ to $\mathscr{E}_{\rm noise}^{m}$,

[TABLE]

When $r>1$, the mode is resolved on the grid, while when $r<1$ it is dominated by the Poisson noise in the simulation and rapidly randomized. Because $r\propto\lambda^{2}$ this statement is also a function of wavelength, with the smallest wavelength modes being the most marginal. That is, the Poisson noise limit, $\mathscr{E}_{\rm noise}^{m}$, sets a minimum mode wavelength, $\lambda_{\rm min}$ the simulation can resolve, independent of the spatial resolution of the grid. Improving the fidelity of the simulation requires, therefore, concurrent increases in spatial resolution (e.g., $h$), momentum resolution (e.g., $N_{\rm pc}$), and spectral resolution (e.g., $L$).

Explicitly, requiring that modes on the smallest spatial scales are resolved, i.e., $\lambda_{\rm min}\propto h$, then translates into the requirement that $hN_{pc}/L$ is, at least, fixed.101010Such requirement (a fixed ratio of the energy in the shortest wavelength mode to the Poisson energy) implies that the ratio of the energy in a specific mode to the energy in the Poisson noise of simulation increases as $\eta^{2}$ with improving resolution. When spectral resolution is not important, this requires that if the spatial resolution increases by a factor $\eta$ then $N_{\rm pc}$ must grow by a similar factor, i.e.,

[TABLE]

This is unsurprising – convergence requires simultaneous and equal increases in the spatial and momentum resolutions. It does mean, however, that convergence studies are numerically demanding, as they scale as $\eta^{3}$, even in 1D. This is, of course, exactly the factor one would find in a purely Eulerian scheme for solving the 1D Boltzmann equation, which is similar to 2D hydrodynamics. It does make clear, however, that the inherent randomness of the particle description does nothing to improve the convergence characteristics.

6.3 SHARP-1D Convergence – An Explicit Example

We now provide an explicit example of convergence testing, as described in the previous section, using SHARP-1D. To illustrate both the convergence of SHARP-1D under this definition, and equally importantly, the lack of convergence under separate definitions often employed, we do this for an extreme range of $\eta$, extending over two orders of magnitude.

We begin with a fiducial simulation, which defines $\eta=1$. This is comprised of a population of electrons, with total number of macro-particles of $N_{p}=8950$, and a fixed neutralizing background. We use a box with normalized length $\bar{L}=39.96175$ and the initial normalized temperature for electrons $\theta_{ini}=10^{-3}$. We start with a single excited mode with amplitude $A=10^{-2}$ and wavelength of $\bar{\lambda}=\bar{L}/2$. Therefore, $\hat{k}=2\pi\sqrt{\theta}/\bar{\lambda}=0.0099441$, i.e., the linear perturbation should oscillate without damping during the entire simulation time, $T=100\omega_{p}^{-1}$.

For our fiducial simulation, we set the cell size, $\Delta x$, such that $1/h=c/(\Delta x\omega_{p})=8.958$, i.e., $N_{c}=358$ and $N_{pc}=N_{p}/N_{c}=25$. For all simulations in this section, we use fifth-order interpolation ($W^{5}$). Note that in this simulation the box is sufficiently large to spectrally resolve all relevant features of the dispersion relation, and thus we do not consider it further here.

The accuracy measure we employ is the normalized maximum error over the duration of the simulation:

[TABLE]

where $\Delta\mathscr{E}$ is the energy change in the total energy, $\mathscr{E}_{\rm th}$ is the initial thermal energy, i.e., excluding rest mass energy. Using various definitions of the energy error, i.e., average error, result in qualitatively identical results. Note that this is not the only accuracy measure we might use; others include the amplitude or phase of the wave, or the ability to reproduce other known solutions. It does have the property that it is fundamentally well understood (energy is conserved), not explicitly conserved by the code (like momentum), and universally defined.

Figure 17 shows the impact of increasing the various relevant notions of resolution independently and together. Increasing either the spatial or momentum resolution independently leads to a plateau in $\epsilon_{m}$. The location and magnitude of this plateau depends upon the value of the non-converging resolution, i.e., when converging in spatial resolution, $N_{pc}$, or when converging in momentum resolution, $h$. This is qualitatively distinct from the case when both resolutions are increased simultaneously, for which no plateau is evident over two orders of magnitude.

Quantitatively, using our definition of convergence, SHARP-1D converges as $\eta^{-2.8}\propto h^{2.8}$. This is very similar to the anticipated $h^{3}$, with the implication that the algorithm performance is well understood. Because this is ultimately set by the currently second-order symplectic integration employed in the particle pusher, implementing higher-order spatial interpolation will not improve this convergence rate. However, the value of the higher-order spatial interpolation appears in the magnitude of the energy error (i.e., its coefficient). This is clearly evident in Figure 18, which shows the error in the $\eta=1$ simulation for different interpolation orders. Nevertheless, achieving the full benefits of the higher-order spatial implementations will require implementing an appropriate particle pusher.

7 Performance of SHARP-1D

To quantify the increase in the computational cost of using higher-order interpolation functions, we ran a simulation using 8950 macro-particles on five processors for all implemented orders. In Figure 19, we show the relative increase in computational time for both deposition and back-interpolation steps after running each simulation up to $t~{}\omega_{p}=$ $10^{3}$ (red), $10^{5}$ (green), and $10^{6}$ (black). The computational cost per update of using $W^{5}$ is $2.28$ times larger than $W^{1}$. We have verified that this is independent of the number of macro-particles and the number of processors.

The advantages of using higher-order interpolation are problem dependent. However, if we use the error in the total energy as a measure of accuracy, we can attempt to quantify the difference by computing the relative computational cost of simulations with different interpolation orders holding the level of accuracy fixed. Figure 18 shows that for a simulation that runs until time $t=10^{6}~{}\omega_{p}^{-1}$, the energy error, $\epsilon_{m}$, is smaller by a factor of $10^{3}$ when $W^{5}$ is used instead of $W^{1}$.

To achieve a similar accuracy using $W^{1}$, i.e., decreasing $\epsilon_{m}$ by a factor of $10^{3}$, $\eta$ needs to be increased by $10^{3/2.8}\sim 11.8$ (where we employed the scaling in the right-hand panel of Figure 17). Consequently, both $N_{\rm pc}$ and $N_{\rm c}$ each have to increase by $11.8$, increasing the number of steps by the same factor. Therefore, using $W^{1}$, the computational cost increases by a factor of $11.8^{3}\sim 1640$, $\sim 730$ times that required by $W^{5}$. That is, to achieve the same level of accuracy, a simulation that uses $W^{5}$ is about $730$ faster than a simulation that uses $W^{1}$ with improved resolutions.111111 For reference, on Intel Xeon 3.47 GHZ CPUs, the computational time to evolve 8950 macro-particles on five processors for 22396417 steps is 1732 seconds per processor when using $W^{1}$.

SHARP-1D exhibits a near linear strong scaling, i.e., for fixed problem size, with the number of processors employed, $N_{\rm pr}$. By varying $N_{\rm pr}$ between 20 and 300, we find

[TABLE]

8 Conclusion

In this paper, we present a self-consistent discretization for the governing equations of plasma made of macro-particles in 1D (i.e., the Vlasov–Poisson equations) implemented in the SHARP-1D code. It employs a self-consistent force on such macro-particles that is accurate up to fifth order and provides an essential step toward higher-order accurate PIC schemes. The over all accuracy of the algorithm is, however, limited by the accuracy of the particle pusher which is still a second order accurate symplectic method (leap-frog).

SHARP-1D conserves momentum exactly, and despite its second-order accuracy, when higher-order interpolation functions are used, better energy conservation and lower numerical heating is evident. SHARP-1D simulations of a thermal plasma, whose Debye length is 10 times smaller than the cell size and which employ spatial interpolation accurate to fifth order only have an energy error, which is better than 1% of the initial thermal energy. Moreover, it shows a negligible numerical heating over a very long time (up to millions of inverse plasma frequencies, see Figure 3).

We present a validation of SHARP-1D against some test problems: the thermal stability of plasmas, the stability of linear plasma modes, and the two-stream instability in the relativistic and non-relativistic regimes.

To perform such validation tests, we determine the correct modes of thermal plasmas (oscillation frequencies and damping and growing rates) by solving the corresponding linear dispersion relations numerically. This is done for thermal plasmas that are both stationary and counter streaming. For convenience, we provide a fit to the oscillation frequencies and the damping rates in the linear regime of thermal plasmas up to $\hat{k}=kv_{\rm th}/\omega_{p}=0.6$.

In all test problems, SHARP-1D demonstrates the ability to reproduce kinetic effects of the linear regime both qualitatively and quantitatively. This includes reproducing the correct oscillation frequencies and damping rates for different modes of the thermal plasma, and also reproducing all oscillating, growing, and damping modes in counter streaming plasmas in the relativistic and non-relativistic regimes.

Results from SHARP-1D in both relativistic and non-relativistic regimes are contrasted with results from TRISTAN-MP. A substantial improvement in the ability to conserve energy and control numerical heating is shown when SHARP-1D is used.

Importantly, the improved performance due to higher spatial order does not come at the cost of increased execution time; to achieve the same level of accuracy, we have shown that, for SHARP-1D, a simulation with $W^{5}$ is almost three orders of magnitude faster than a simulation with $W^{1}$ and improved resolutions.

Finally, we develop a general criterion for convergence studies of PIC simulations by requiring that the ratio of the energy in the shortest wave mode to the energy in the Poisson noise of the simulation to be at lease fixed.

An example study for such a convergence test is presented, where the decrease in the energy error for plasma, as different relevant resolutions are increased, is used as a measure for convergence. Both the number of particles per cell and the spatial resolution of the grid are crucial resolution elements: increasing only one relevant resolution results not only in slower decrease of the error, but also in a plateau where the error does not decrease any longer as such resolution increases. Faster decrease of the error without any plateau is achieved when all relevant resolutions are increased simultaneously for such plasma.

This new PIC code provides a new avenue that enables the faithful study of the long-term evolution of plasma problems (in one dimension) that require absolute control of the energy and momentum conservation. Those include, e.g., the oblique instability driven by the highly anisotropic TeV pair beams that emerge from TeV gamma-rays that propagate from blazars to us or interactions of relativistic plasma components with a non-relativistic background plasma over long time scales (Broderick et al., 2012).

Acknowledgments

We would like to thank Anatoly Spitkovsky for providing access to TRISTAN-MP, substantial guidance in its use, and a number of constructive comments on the manuscript. M.S. and A.E.B. receive financial support from the Perimeter Institute for Theoretical Physics and the Natural Sciences and Engineering Research Council of Canada through a Discovery Grant. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Research and Innovation. P.C. gratefully acknowledges support from the NASA ATP program through NASA grant NNX13AH43G, and the NSF through grant AST-1255469. C.P. gratefully acknowledges support by the European Research Council through ERC-CoG grant CRAGSMAN-646955 and by the Klaus Tschira Foundation. E.P. gratefully acknowledges support by the Kavli Foundation. Support for A.L. was provided by an Alfred P. Sloan Research Fellowship, NASA ATP Grant NNX14AH35G, and NSF Collaborative Research Grant #1411920 and CAREER grant #1455342.

Appendix A Momentum conservation

For the purpose of comparison, we start by calculating the correct interaction force for point particle in 1D, by using

[TABLE]

Therefore, the correct interaction term in 1D for point-particles ($m=0$) (in code units) is given by

[TABLE]

On the other hand, the force on a macro-particle, with charge $q_{\alpha}$ and centered at $x_{\alpha}\epsilon[0,L]$ on a periodic box, is given by $\bar{F}_{\alpha}\equiv F_{\alpha}/(q_{0}\mathbb{E}_{0})=\bar{Q}_{\alpha}\bar{E}_{\alpha}$, where

[TABLE]

By defining

[TABLE]

The exact equations for the fields on the grid edges are given by Equations (35,36), these can be written as

[TABLE]

The solution of Equation (75) can be expressed as

[TABLE]

where $A_{jk}$ is anti-symmetric matrix given by

[TABLE]

Two possible, second-order accurate, approximations for Equation (73) can be written as follows

[TABLE]

A.1 Non-momentum Conserving Second-order scheme: approximation (78)

If we use (76) and net-charge neutrality, i.e., $\sum_{k}\bar{\rho}_{k}=\sum_{k}\bar{\rho}_{k+\frac{1}{2}}=0$, then the total force on all macro-particles is given by

[TABLE]

Therefore, the net-force on the system does not depend on the choice of $E_{0}$ (because of charge neutrality), and since using higher interpolation functions makes the variation in the interpolated density smoother, it decreases both terms in Equation (80), i.e., using higher-order interpolation improves the momentum conservation.

A.1.1 Origin of Momentum Non-conservation: Self-forces and wrong interaction forces

To see the origin of momentum non-conservation, we examine the interpolated force in the case in which there are only two macro-particles. Using $\bar{Q}_{1}+\bar{Q}_{2}=0$, $\sum_{k}W^{m}_{k,i_{s}}=1$, and

[TABLE]

the force on macro-particle at $x_{1}$ is given by

[TABLE]

Therefore, the non-vanishing self-force, $\bar{F}_{\text{self}}$, and the numerical interaction force, $\bar{F}_{\text{int}}$, are given by

[TABLE]

Therefor, the self-force ($\sim h/N_{pc}=L/N_{p}$) here vanishes only in the limit of infinite number of macro-particles $N_{p}\rightarrow\infty$.

A.2 Momentum Conserving Second-order Scheme: approximation (79)

If we use (76) and net-charge neutrality, i.e., $\sum_{k}\bar{\rho}_{k}=\sum_{k}\bar{\rho}_{k+1/2}=0$, and

[TABLE]

The net-force is, then, given by

[TABLE]

Therefore, the net-force on the system, again, does not depend on the choice of $E_{0}$ (because of charge neutrality), and it is always exactly zero, therefore, the momentum is exactly conserved. Hence we call this a momentum conserving scheme.

A.2.1 Vanishing of the Self-forces exactly

To see how self-forces vanish, we examine the interpolated force in the case in which there are only two macro-particles. Using $\bar{Q}_{1}+\bar{Q}_{2}=0$, $\sum_{k}W^{m}_{k,i_{s}}=1$ and

[TABLE]

the force on macro-particle at $x_{1}$ is given by

[TABLE]

Therefore, the vanishing self-force, $\bar{F}^{\text{self}}$, and the numerical interaction force, $\bar{F}^{\text{int}}$, are given by

[TABLE]

We can see from comparing Equations (72) and (86) that the numerically calculated interaction force is shorter in range, and that the usage of higher-order interpolation (larger macro-particles) results in smoothing the interaction forces.

A.3 Filtering for Momentum Conserving Scheme: Non-vanishing Self-forces

Here, we show that filtering of the grid charge densities, e.g. 1D equivalence of filtering implemented in TRISTAN-MP, lead to a non-vanishing self-forces and a violation of momentum conservation.

After one-filter, the filtered charge densities $\bar{\rho}^{\rm f}_{k}$ are given by

[TABLE]

Therefore, for a momentum conserving scheme, e.g., Equation (79), if we used (83) and replaced $\rho_{k+1/2}$ by $\rho_{k+1/2}^{\rm f}$, the net force is given by

[TABLE]

Which shows that filtering the grid deposited moment (grid charge density) leads to violation of momentum conservation in otherwise a momentum conserving scheme.

To see the origin of such violation, we look, as we did before, at the interpolated force in case of having only two macro-particles, using $\bar{Q}_{1}+\bar{Q}_{2}=0$, $\sum_{k}W^{m}_{k,i_{s}}=1$ and

[TABLE]

the force on macro-particle at $x_{1}$ is given by

[TABLE]

Therefore, filtering leads to wrong interaction forces and a non-vanishing self-force given by

[TABLE]

Appendix B Shape and weight functions explicit form

If we define $y=(x-x_{i})/\Delta x=(\bar{x}-\bar{x}_{i})/h$, the shape functions, $S^{m}(y)$, assumed for the macro-particles and the corresponding weight functions, $W^{m}(y)$, used in the interpolation steps in the code are given in Table B. {deluxetable*}p0.5cmp8cmp8cm \tabletypesize \tablecolumns3 \tablewidth0pt \tablecaptionShape and weight functions implemented in the code.

\tablehead $m$ $S^{m}(y)$ $W^{m}(y)$ \startdata

0

$\delta(y)$

$\begin{cases}1&0\leq|y|\leq\dfrac{1}{2}\\ 0&\text{Otherwise}\end{cases}$

1

$\dfrac{1}{\Delta x}\begin{cases}1&0\leq|y|\leq\dfrac{1}{2}\\ \\ 0&\text{Otherwise}\end{cases}$

$\begin{cases}1-|y|&0\leq|y|<1\\ \\ 0&\text{Otherwise}\end{cases}$

2

$\dfrac{1}{\Delta x}\begin{cases}1-|y|&0\leq|y|<1\\ \\ 0&\text{Otherwise}\end{cases}$

$\begin{cases}\dfrac{3}{4}-y^{2}&0<|y|<\dfrac{1}{2}\\ \\ \dfrac{1}{8}(3-2|y|)^{2}&\dfrac{1}{2}\leq|y|<\dfrac{3}{2}\\ \\ 0&\text{Otherwise}\end{cases}$

3

$\dfrac{1}{\Delta x}\begin{cases}\dfrac{3}{4}-y^{2}&0<|y|<\dfrac{1}{2}\\ \\ \dfrac{1}{8}(3-2|y|)^{2}&\dfrac{1}{2}\leq|y|<\dfrac{3}{2}\\ \\ 0&\text{Otherwise}\end{cases}$

$\begin{cases}\dfrac{2}{3}-y^{2}+|y|^{3}/2&0<|y|<1\\ \\ \dfrac{1}{6}(2-|y|)^{3}&1\leq|y|<2\\ \\ 0&\text{Otherwise}\end{cases}$

4

$\dfrac{1}{\Delta x}\begin{cases}\dfrac{2}{3}-y^{2}+|y|^{3}/2&0<|y|<1\\ \\ \dfrac{1}{6}(2-|y|)^{3}&1\leq|y|<2\\ \\ 0&\text{Otherwise}\end{cases}$

$\begin{cases}\dfrac{115}{192}-\dfrac{5y^{2}}{8}+\dfrac{y^{4}}{4}&0<|y|<\dfrac{1}{2}\\ \\ \dfrac{1}{96}\left[55+20|y|-120y^{2}+80|y|^{3}-16y^{4}\right]&\dfrac{1}{2}\leq|y|<\dfrac{3}{2}\\ \\ \dfrac{1}{384}(5-2|y|)^{4}&\dfrac{3}{2}\leq|y|<\dfrac{5}{2}\\ \\ 0&\text{Otherwise}\end{cases}$

5

$\displaystyle\dfrac{1}{\Delta x}\begin{cases}\dfrac{115}{192}-\dfrac{5y^{2}}{8}+\dfrac{y^{4}}{4}&0<|y|<\dfrac{1}{2}\\ \\ \dfrac{1}{96}\left[55+20|y|-120y^{2}+80|y|^{3}-16y^{4}\right]&\dfrac{1}{2}\leq|y|<\dfrac{3}{2}\\ \\ \dfrac{1}{384}(5-2|y|)^{4}&\dfrac{3}{2}\leq|y|<\dfrac{5}{2}\\ \\ 0&\text{Otherwise}\end{cases}$

$\displaystyle\begin{cases}\dfrac{11}{20}-\dfrac{y^{2}}{2}+\dfrac{y^{4}}{4}-\dfrac{|y|^{5}}{12}&0\leq|y|\leq 1\\ \\ \dfrac{17}{40}+\dfrac{5|y|}{8}-\dfrac{7y^{2}}{4}+\dfrac{5|y|^{3}}{4}-\dfrac{3y^{4}}{8}+\dfrac{|y|^{5}}{24}&1<|y|<2\\ \\ \dfrac{1}{120}(3-|y|)^{5}&2\leq|y|<3\\ \\ 0&\text{Otherwise}\end{cases}$ \enddata

Appendix C Aliasing

The Fourier components of grid quantities, in 1D, $g_{k}$ are such that $g_{k}=g_{k+pk_{g}}$, where $p$ is some integer and $k_{g}=2\pi/\Delta x$ is the wave mode associated with the cell size $\Delta x$ on that physical grid. Therefore, for a continuous particles number density $n(x)$, the Fourier component of the grid charge density is given by (Birdsall & Langdon, 1991)

[TABLE]

where $\tilde{S}(k)$ is Fourier transform of our interpolation function and $\tilde{n}(k)$ is Fourier transform of $n(x)$. Therefore, all aliases of $k$ (wave modes that differ from $k$ by integer number of $k_{g}$) contribute when grid quantities are calculated. Clearly, this will feedback on the particle quantities, when the grid quantities are used to calculate the force on the particles to evolve them. The strength of the coupling between aliases (the source of this error) depends on how fast $\tilde{S}^{m}(k)$ falls off for large $k$, as can be seen in Equation (90).

The Fourier transform of our interpolation functions (spline functions of order $m$, see Table B) is given by

[TABLE]

Therefore, Using higher-order interpolation functions (larger $m$) in our code leads to a decrease in the strength of the couplings between grid wave modes and their aliases), which results in improvements in energy conservation as seen in Section 4.1.

Appendix D Poisson noise

Here we calculate the noise when a finite number of computational particles are used to represent a uniform distribution function. We calculate the total energy density due to such noise in Appendix D.1 and then find the power spectrum for such noise in Appendix D.2. In Appendix D.3, we calculate the temperature, $\theta_{p}$, set by the energy in such noise.

D.1 Average Potential Energy From Uniformly Distributed Macro-particles

Using Equation (10) and the first equation in (2), we can write the electric field associated with plasma particles on a periodic box of length $L$, i.e., $x_{i_{s}}\in[0,L)$ as follows

[TABLE]

Where, $q^{m}(x,x_{i_{s}})=\int_{0}^{x}dx^{{}^{\prime}}S^{m}(x^{{}^{\prime}},x_{i_{s}})$. The periodicity of the box implies that the plasma is neutral. Therefore,

[TABLE]

The spatial averaging of $q^{m}$, for uniformly distributed macro-particles, is

[TABLE]

For such macro-particles the average of the electric field is zero:

[TABLE]

However, due to the finite number of macro-particles, the average potential energy is non-zero, to calculate such energy we need to calculate

[TABLE]

For the shape functions implemented in SHARP-1D (their explicit forms are given in Appendix B), the integral in (LABEL:ex2) is given by

[TABLE]

Therefore,

[TABLE]

The average electrostatic potential energy due to the finite number of macro-particles is, then, given by

[TABLE]

Here, $N_{c}=L/\Delta x$ is the number of macro-cells. If we assume that all plasma species have the same mass, and absolute value of charge, we then make the choice of our fiducial units as, $q^{2}_{0}=Q^{2}_{s}$ and $m_{0}=M_{s}$ (that implies $n_{0}=\sum_{s}n_{s}$). Therefore,

[TABLE]

Equation (100) shows that using higher-order shape functions decreases the noise coming from the fact that we are using a finite number of macro-particles. The decrease that we gain in the potential energy noise is $f_{m}\bar{L}^{2}/2N_{c}$. For a given box size, this improvement is lowered, if we increase the number of cells $N_{c}$ because it means a decreases in cell size, which means also a decrease in the size of the macro-particles. On the other hand, if we increase the number of cells while keeping the cell size fixed, i.e., by increasing the box-size $\bar{L}$, that improvement due to using higher-order interpolation functions increases.

D.2 Spectrum of the Poison Noise

To find the spectrum of such Poisson noise, we average the Fourier transform of the grid electric fields. Using Equation (92), the Fourier components of the electric field are given by

[TABLE]

By defining

[TABLE]

and using the fact that $E(x)$ is a real valued function, the absolute value for such Fourier components are given by

[TABLE]

Averaging over a periodic box and assuming the macro-particles are uniformly distributed, we can write

[TABLE]

By using Equations (102, 94)

[TABLE]

Therefore, the averaged magnitude for the Fourier components can be written as

[TABLE]

If all plasma species have the same mass, and absolute value of charges, such average can be written in code units as

[TABLE]

D.3 Heating Due to Noise

If the energy due to Poisson noise is converted to heat that puts a floor in the temperature PIC scheme can simulate. Here we estimate such a temperature floor $\theta_{p}$. The energy due to Poisson noise is calculate in Appendix D.1. It is given by

[TABLE]

where $f_{m}$ is defined in (97). This noise is due to the finite number of macro-particles used in the simulations. If this energy is converted to thermal energy, it would lead to heating of the plasma up to a temperature $\theta_{p}$. If the plasmas are at thermal equilibrium, the momentum distribution is given by Maxwell–Jüttner distribution, and hence the temperature of plasmas is related to the kinetic energy, $\mathscr{\bar{K}}$, as follows

[TABLE]

Here, $\theta_{s}=k_{B}T_{s}/m_{0}c^{2}$ is the normalized temperature of species $s$ with $N_{s}$ of macro-particles, $N$ is the total number of macro-particles from all species and $K_{0}$, $K_{1}$ are the Bessel functions of zeroth and first kind respectively. Therefore,

[TABLE]

Hence, if a plasma of macro-particles starts with a temperatures below $\theta_{p}$ the Poisson noise will non-physically heat such plasma.

Appendix E Dispersion relation For Non-relativistic warm plasma

For non-relativistic ($\gamma^{3}\approx 1$) warm plasma, i.e., $0<\theta\ll 1$, we can write

[TABLE]

If we assume no net current in the plasmas, i.e., the momentum distribution of all species is such that $\sum_{s}Q_{s}\int vf^{s}_{0}(v)dv=0$, then the linear dispersion relation of uniformly distributed plasma is given by

[TABLE]

where $v_{p}=\hat{\omega}/\hat{k}$, $\hat{\omega}=\omega/\omega_{p}$, $\hat{k}=kc\sqrt{\theta}/\omega_{p}$, i.e., $v_{p}=\omega/kc\sqrt{\theta}$, $\omega^{2}_{p}=\sum_{s}\omega^{2}_{ps}$, $\omega^{2}_{ps}=Q_{s}^{2}n_{s}/\epsilon_{0}M_{s}$, and $\bar{v}=v/c$.

If we assume that $\Im(\hat{\omega})>0$, i.e., $\Im(\omega)>0$, then

[TABLE]

where $z\equiv\bar{v}/\sqrt{\theta}$. Extending the definition of $\chi_{s}(v_{p})$ to the entire complex plane can be done as follows (Brambilla, 1998)

[TABLE]

where $\hat{\omega}_{s}\equiv\omega_{ps}/\omega_{p}$. Therefore, for $\Im(v_{p})\neq 0$, we can then write

[TABLE]

where Erfi is the complex error function, which is defined as $\text{Erfi}(v_{p})=-i\text{ Erf}(iv_{p})$.

E.1 Standing Linear Plasma Waves

In the case of thermal electrons with fixed neutralizing background ($\hat{\omega}_{s}=1$), the dispersion relation is then given by

[TABLE]

E.2 Two-stream Instability

In the case of two population of thermal electrons (both have the same number density), propagating in two opposite directions with speed $v_{b}$, with fixed neutralizing background, therefore $\hat{\omega}^{2}_{s}=1/2$, and the linear dispersion relation is then given by ($z_{b}\equiv v_{b}/c\sqrt{\theta}$)

[TABLE]