A semi-implicit, energy- and charge-conserving particle-in-cell algorithm for the relativistic Vlasov-Maxwell equations

TL;DR

This paper introduces a semi-implicit particle-in-cell algorithm that conserves energy and charge exactly while maintaining accurate light-wave dispersion, improving simulations of relativistic plasma phenomena.

Contribution

It combines leap-frog and Crank-Nicolson methods to create an energy- and charge-conserving PIC algorithm with accurate light-wave dispersion properties.

Findings

Exact energy and charge conservation achieved.

Preserves light-wave dispersion properties.

Validated with relativistic plasma instability simulations.

Abstract

Conventional explicit electromagnetic particle-in-cell (PIC) algorithms do not conserve discrete energy exactly. Time-centered fully implicit PIC algorithms can conserve discrete energy exactly, but may introduce large dispersion errors in the light-wave modes. This can lead to intolerable simulation errors where accurate light propagation is needed (e.g. in laser-plasma interactions). In this study, we selectively combine the leap-frog and Crank-Nicolson methods to produce an exactly energy- and charge-conserving relativistic electromagnetic PIC algorithm. Specifically, we employ the leap-frog method for Maxwell's equations, and the Crank-Nicolson method for the particle equations. The semi-implicit algorithm admits exact global energy conservation, exact local charge conservation, and preserves the dispersion properties of the leap-frog method for the light wave. The algorithm employs a new particle pusher designed to maximize efficiency and minimize wall-clock-time impact vs. the explicit alternative. It has been implemented in a code named iVPIC, based on the Los Alamos National Laboratory VPIC code (\url{https://github.com/losalamos/vpic}). We present numerical results that demonstrate the properties of the scheme with sample test problems: relativistic two-stream instability, Weibel instability, and laser-plasma instabilities.

Equations

$${\partial }_t{f}_s+\mathbf{v}\cdot \nabla f_s+\frac{q_s}{m_s}(\mathbf{E}+\mathbf{v}\times \mathbf{B})\cdot {\nabla }_{\mathbf{u}}{f}_s$$

$$\frac{\partial \mathbf{B}}{\partial t}+\nabla \times \mathbf{E}$$

$${ϵ}_0\frac{\partial \mathbf{E}}{\partial t}-\frac{1}{{\mu }_0}\nabla \times \mathbf{B}+\mathbf{j}$$

$$\nabla \cdot \mathbf{B}$$

$$\nabla \cdot \mathbf{E}$$

$$\mathbf{j}$$

$$\rho$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \mathbf{j}=0$$

$$f_s(\mathbf{x},\mathbf{u})\approx \sum _{p\in s}{w}_p\delta (\mathbf{x}-{\mathbf{x}}_p)\delta (\mathbf{u}-{\mathbf{u}}_p),$$

$$\frac{d{\mathbf{x}}_p}{dt}$$

$$\frac{d{\mathbf{u}}_p}{dt}$$

$$\frac{{\mathbf{u}}_p^{n+\text{nicefrac}12}-{\mathbf{u}}_p^{n-\text{nicefrac}12}}{\Delta t}$$

$$\frac{{\mathbf{x}}_p^{n+1}-{\mathbf{x}}_p^n}{\Delta t}$$

$$\frac{{\mathbf{B}}_h^{n+\text{nicefrac}12}-{\mathbf{B}}_h^{n-\text{nicefrac}12}}{\Delta t}$$

$${ϵ}_0\frac{{\mathbf{E}}_h^{n+1}-{\mathbf{E}}_h^n}{\Delta t}$$

$$\frac{{\mathbf{B}}_h^{n+\text{nicefrac}12}-{\mathbf{B}}_h^n}{\Delta t/2}$$

$$\frac{{\mathbf{B}}_h^{n+1}-{\mathbf{B}}_h^{n+\text{nicefrac}12}}{\Delta t/2}$$

$$\frac{{\mathbf{x}}_p^{n+1}-{\mathbf{x}}_p^n}{\Delta t}$$

$$\frac{{\mathbf{u}}_p^{n+1}-{\mathbf{u}}_p^n}{\Delta t}$$

$${\mathbf{v}}_p^{n+\text{nicefrac}12}$$

$${\bar{\mathbf{v}}}_p^{n+\text{nicefrac}12}$$

$${\mathbf{E}}_p^{n+\text{nicefrac}12}$$

$${\mathbf{B}}_p^{n+\text{nicefrac}12}$$

$${\mathbf{E}}_h^{n+\text{nicefrac}12}=\frac{{\mathbf{E}}_h^{n+1}+{\mathbf{E}}_h^n}{2}$$

$${\mathbf{j}}_h^{n+\text{nicefrac}12}=\frac{1}{{\mathbf{\Delta }}_h}\sum _p{q}_p{\mathbf{v}}_p^{n+\text{nicefrac}12}\cdot \overline{\stackrel{ˉ}{\mathbf{S}}}({\mathbf{x}}_h-{\mathbf{x}}_p).$$

$$\sum _h{\mathbf{\Delta }}_h\left[{ϵ}_0({\mathbf{E}}_h^{n+1}-{\mathbf{E}}_h^n)\cdot {\mathbf{E}}_h^{n+\text{nicefrac}12}-\frac{\Delta t}{{\mu }_0}({\nabla }_h\times {\mathbf{B}}_h^{n+\text{nicefrac}12})\cdot {\mathbf{E}}_h^{n+\text{nicefrac}12}+\Delta t{\mathbf{j}}_h^{n+\text{nicefrac}12}\cdot {\mathbf{E}}_h^{n+\text{nicefrac}12}\right]=0.$$

$$\frac{{ϵ}_0}{2}\sum _h{\mathbf{\Delta }}_h({\mathbf{E}}_h^{n+1}-{\mathbf{E}}_h^n)\cdot {\mathbf{E}}_h^{n+\text{nicefrac}12}=\frac{{ϵ}_0}{2}\sum _h{\mathbf{\Delta }}_h\left[({\mathbf{E}}_h^{n+1}{)}^2-({\mathbf{E}}_h^n{)}^2\right]\equiv W_E^{n+1}-W_E^n,$$

$$-\frac{\Delta t}{{\mu }_0}\sum _h{\mathbf{\Delta }}_h({\nabla }_h\times {\mathbf{B}}_h^{n+\text{nicefrac}12})\cdot {\mathbf{E}}_h^{n+\text{nicefrac}12}$$

$$W_B^n\equiv \frac{1}{2{\mu }_0}\sum _h{\mathbf{\Delta }}_h{\mathbf{B}}_h^{n+\text{nicefrac}12}\cdot {\mathbf{B}}_h^{n-\text{nicefrac}12},$$

$$\Delta t\sum _h{\mathbf{\Delta }}_h{\mathbf{j}}_h^{n+\text{nicefrac}12}\cdot {\mathbf{E}}_h^{n+\text{nicefrac}12}$$

$${\mathbf{v}}_p^{n+\text{nicefrac}12}\cdot {\bar{\mathbf{v}}}_p^{n+\text{nicefrac}12}\times {\mathbf{B}}_p^{n+\text{nicefrac}12}=0.$$

$${\left(W_E+W_B+W_p\right)|}_n^{n+1}=0.$$

$$\frac{({\rho }_p{)}_{ijk}^{n+1}-({\rho }_p{)}_{ijk}^n}{\Delta t}+{{\nabla }_h\cdot {\mathbf{j}}_p^{n+\text{nicefrac}12}|}_{ijk}=0.$$

$${\rho }_{p,ijk}^n=q_p{S}_1(x_i-x_p^n)S_1(y_j-y_p^n)S_1(z_k-z_p^n),$$

$$j_{x,p,i+\text{nicefrac}12,j,k}^{n+\text{nicefrac}12}$$

$$j_{y,p,i,j+\text{nicefrac}12,k}^{n+\text{nicefrac}12}$$

$$j_{z,p,i,j,k+\text{nicefrac}12}^{n+\text{nicefrac}12}$$

$${\mathbb{S}}_{11,jk}^{n+\text{nicefrac}12}(y_p,z_p)\equiv \frac{1}{3}$$

Full text

1 Introduction

Particle-in-cell (PIC) methods combine Lagrangian and Eulerian techniques for solving the relativistic Vlasov-Maxwell equations (or a subset thereof) for kinetic plasma simulations. Specifically, the method of characteristics is employed to solve Vlasov equation, evolving macro-particles according to the Newton (or Lorentz) equations of motion. Finite-difference time-domain (FDTD) methods are commonly used for evolving Maxwell’s equations on a mesh. The Yee scheme [1] is the most popular algorithm for integrating Maxwell’s equations, due to its simplicity and second-order accuracy. B-splines are widely used for interpolations between the particles and the mesh. A detailed description of classical PIC methods and their analysis can be found in Ref. [2, 3] .

Despite their numerical simplicity, PIC algorithms have been successfully applied to a large range of plasma descriptions, including electrostatic [4], low-frequency electromagnetic (Darwin [5]), gyrokinetic [6], quasi-static [7], and the most general electromagnetic and relativistic systems [8]. However, the applicability of PIC for long-term multiscale simulations may be currently limited by the trustworthiness of the simulation, as measured by accumulated errors of conservation properties such as charge, momentum, and energy. Hereafter, we will adopt the notion of “discrete conservation” as conservation satisfied exactly (in practice to numerical roundoff) in a discrete system (with finite timestep $\Delta t$ and mesh size $\Delta x$). Discrete charge conservation is satisfied by most PIC algorithms employed today [9, 10, 11, 12]. Discrete energy conservation has been recently demonstrated with advances in implicit PIC algorithms, either fully implicit [11, 13, 14, 15, 16], or semi-implicit [17] (note that an early version of “energy-conserving” algorithm by Lewis [18] conserves discrete energy exactly only in the limit of $\Delta t\rightarrow 0$). Discrete momentum conservation has only been realized in the simplest electrostatic models integrated by explicit schemes [2, Ch.8.6].

In this study, we seek to develop an energy- and charge-conserving scheme for the relativistic Vlasov-Maxwell system. Compared to the fully implicit approach proposed by Lapenta and Markidis [5, 14], which is based on the Crank-Nicolson (CN) scheme, our approach is less dispersive for light wave modes and consequently less subject to numerical Cerenkov (or Cherenkov) radiation issues. For the CN scheme solving Maxwell’s equations, the numerical Cerenkov radiation can be so problematic that, in practice, some numerical dissipation is needed at the cost of discrete energy conservation [13]. It is also worth noting that the methods developed by Lapenta and his group do not conserve charge (see Ref. [19] for an improvement on this issue). Energy- and charge-conserving PIC algorithms based on fully implicit approaches have been developed, but so far only for electrostatic [11] and Vlasov-Darwin systems (e.g. [15]).

Based on a semi-implicit approach, the proposed energy- and charge-conserving PIC algorithm has relatively low light-wave dispersion (w.r.t. the CN scheme) for electromagnetic relativistic plasma simulations. The algorithm combines the conventional leap-frog scheme for the field equations with a CN update of the particle equations of motion. The CN scheme combines the position update of Ref. [14] and velocity update of Ref. [20]. Energy conservation is achieved by implicitly coupling the electric field advance and particle push, while charge conservation rests on a novel treatment of particle cell crossings that, unlike previous studies [11, 15], demands no particle stopping at cell boundaries. The leap-frog field update retains a Courant-Friedrics-Lewy (CFL) time-step stability constraint, which is necessary for stability. However, respecting the CFL constraint keeps the numerical light-wave dispersion low, and ensures relatively fast convergence ($\sim 5$ Picard iterations) of the electric-field/particle-push coupling. This approach is mostly motivated by simulations of laser-plasma interactions, where light wave and relativistic effects dominate. However, the algorithm can be used for non-relativistic simulations as well.

In what follows, Sec. 2 briefly reviews the classical electromagnetic relativistic PIC algorithm, and introduces notation. Section 3 introduces the new algorithmic elements that enable the exact discrete energy and charge conservation. Section 4 presents some numerical experiments demonstrating the correctness and long-term conservation properties of the scheme. We close with some concluding remarks in Sec. 5.

2 The classical electromagnetic relativistic PIC method

We consider the relativistic Vlasov-Maxwell system:

[TABLE]

where $f_{s}(\mathbf{x},\mathbf{u})$ is the particle distribution function of species $s$ in phase space, $\mathbf{x}$ denotes physical position, $\mathbf{v}$ denotes velocity, $\mathbf{u}=\mathbf{v}\gamma$ is the proper velocity, and $\gamma=\sqrt{1+u^{2}/c^{2}}$ is the Lorentz factor, with $c$ the speed of light. In these equations, $\epsilon_{0}$ and $\mu_{0}$ are the vacuum permittivity and permeability respectively, $\mathbf{E}$ and $\mathbf{B}$ are the electric and magnetic field, respectively, $q_{s}$ and $m_{s}$ are the species charge and mass respectively, and $\mathbf{j}$ and $\rho$ are current and charge densities, respectively, found from the distribution function as:

[TABLE]

Equations 4 and 5 are the two involutions of Maxwell’s equations [21]. Note that conventional FDTD methods typically advance Eqs. 2 and 3 only [22]. Equation 4 is automatically satisfied by Yee’s scheme, and Eq. 5 may be enforced by specially designed divergence-cleaning methods [23, 24, 25] when charge conservation

[TABLE]

is not automatically satisfied.

Equation 1 is discretized by samples of $f_{s}(\mathbf{x},\mathbf{u})$ (i.e., macro-particles):

[TABLE]

where $w_{p}$ is the particle weight (constant in collisionless PIC simulations). The particle equations of motion read:

[TABLE]

Assuming the use of a Yee mesh, the classical leap-frog scheme is written as:

[TABLE]

where the superscript $n$ denotes time level, the subscript $p$ denotes a particle quantity or a field evaluated at the particle position, $\Delta t$ is the timestep, and the subscript $h=(i,j,k)$ denotes mesh quantities and operators (using Yee finite differences). Here $\gamma^{n}=\sqrt{1+(\mathbf{u}^{n-\nicefrac{{1}}{{2}}}+\frac{q\Delta t}{2m}\mathbf{E}^{n})^{2}}$, and $\gamma^{n+\nicefrac{{1}}{{2}}}=\sqrt{1+(\mathbf{u}^{n+\nicefrac{{1}}{{2}}})^{2}}$. One potential problem with Boris scheme is that it does not preserve the correct limit in a force-free field [26]. The current density $\mathbf{j}_{h}^{n+\nicefrac{{1}}{{2}}}$ is gathered from particles using B-spline interpolation. This scheme features a Courant-Fredrichs-Lewy (CFL) condition [27] for stability, and constrains cell sizes to be comparable to the Debye length to suppress finite-grid instabilities [2].

3 The energy- and charge-conserving electromagnetic relativistic PIC

method

We derive a discrete global energy and local charge conserving scheme for the relativistic Vlasov-Maxwell system by combining the leap-frog time advance for Maxwell equations and CN for the particle equations. Energy conservation is achieved by advancing particles and the electric field synchronously, affording the scheme a semi-implicit character and demanding iteration. This iteration, however, is simple to implement, and of rapid convergence. Simultaneously, automatic local discrete charge conservation is enforced by the proper choice of shape interpolation functions, and a new treatment of particle cell crossings that does not require particles to actually stop at cell boundaries, as previously proposed in Refs. [11, 15].

3.1 Conservative semi-implicit Vlasov-Maxwell PIC algorithm

The proposed algorithm employs the leap-frogged Maxwell update in Eqs. 11, 12, coupled with a time-centered (CN) particle push (given below). Leap-frog is selected to ensure relatively low numerical dispersion errors. It is well-known that the light-wave dispersion relation can be almost perfectly preserved in 1D with the Yee scheme as $c\Delta t/\Delta x\rightarrow 1$. In multiple dimensions, the numerical dispersion varies with modal wavelength, propagation direction, and spatial discretization [22], but it remains tolerable in practice for many applications. In contrast, for arbitrary $\Delta t$, the CN scheme applied to Maxwell equations [28, 29] distorts the light wave dispersion relation, and slows down the light wave phase speed significantly as $k\Delta x\rightarrow\pi$ [30]. Charged particles with speed greater than the phase speed of light waves would lose much of their energy via Cerenkov radiation [31], and may generate significant noise in the simulation. For this reason, we keep leap-frog for the field time advance. For the analysis below, we rewrite Eq. 11 in the following equivalent way:

[TABLE]

Particles are advanced fully implicitly by the CN scheme, yielding:

[TABLE]

where we define

[TABLE]

Note that this particle pusher is similar to that in Ref. [14] in using $\mathbf{v}_{p}^{n+\nicefrac{{1}}{{2}}}$ in the particle position update, but it employs $\bar{\mathbf{v}}_{p}^{n+\nicefrac{{1}}{{2}}}$ for velocity the update, as in Ref. [20]. This push can be made energy- and charge-conserving, as shown in the following sections. The field interpolations are given by:

[TABLE]

Note that the field $\mathbf{E}_{h}$ and $\mathbf{B}_{h}$ at half timestep are obtained differently:

[TABLE]

and $\mathbf{B}_{h}^{n+\nicefrac{{1}}{{2}}}$ is obtained from Eq. 13. Here, $\mathbf{x}_{h}$ is the grid location, $\mathbf{x}_{p}$ is typically chosen to be at the center of particle trajectory [11], and as before the subscript $h=(i,j,k)$ denotes the grid index. The specific form of the shape function dyad $\bar{\bar{\mathbf{S}}}$ and Eq. 17 are important to ensure energy and charge conservation, and will be discussed in Secs. 3.2 and 3.3. We employ a direct inversion of Eq. 16, as described in Ref. [20].

The conservative Vlasov-Maxwell PIC algorithm is closed with the following definition for the current density:

[TABLE]

where $\bm{\Delta}_{h}$ is the cell volume. Note that, we have used identical shape functions for the electric field (Eq. 19) and the current density (Eq. 22) to ensure exact energy conservation [11, 15]. Also note that time-centered update of particle positions and velocities, together with electric field, results in a coupled field-particle system: $\mathbf{x}_{p}$ is a function of the new-time electric field through Eqs. 15 and 16, which in turn determines the current that determines the field via Eq. 12. This, in turn, will require an iterative solve, which will be introduced later (Sec. 3.6).

We prove the energy conservation theorem next.

3.2 Energy conservation

Discrete energy conservation can be readily shown as follows. Multiplying Eq. 12 by $\mathbf{E}_{h}^{n+\nicefrac{{1}}{{2}}}$ and integrating over space, we find:

[TABLE]

The first term in Eq. 23 gives the change of electric energy:

[TABLE]

using Eq. 21. The second term gives the change of magnetic energy:

[TABLE]

where we have used discrete integration by parts, Eqs. 14 and 13, and we have defined the magnetic energy as:

[TABLE]

using that $\mathbf{B}_{h}^{n}=(\mathbf{B}_{h}^{n+\nicefrac{{1}}{{2}}}+\mathbf{B}_{h}^{n-\nicefrac{{1}}{{2}}})/2$ (see Eq. 13). This definition is non-standard, and as we show in Appendix. A, it is almost always well posed as long as the CFL condition is respected.

The last term in Eq. 23 equals the change in kinetic energy:

[TABLE]

where we have used Eqs. 16, 17, 18, and that

[TABLE]

The energy conservation theorem sought follows:

[TABLE]

3.3 Charge conservation

Discrete local charge conservation requires $\partial_{t}\rho+\nabla\cdot\mathbf{j}=0$ at every grid point $(i,j,k)$. Because the charge conservation equation is linear, it is sufficient to enforce this constraint on the contributions of each particle:

[TABLE]

We employ first-order (trilinear) splines for the charge density:

[TABLE]

and mixed zeroth- and first-order splines (indicated by subscripts 0 and 1 respectively) for the current density:

[TABLE]

where, for instance

[TABLE]

and so on. The proof of charge conservation is as follows. We first decompose the density change into three one-dimensional shifts (letting $q_{p}=1$):

[TABLE]

Along each direction, using the particle orbit position equation, it can be shown that (see Appendix B and Ref. [11]):

[TABLE]

where

[TABLE]

and similarly for the$\hat{y}$ and $\hat{z}$ directions. Equation 26 then follows. The proof requires that the particle trajectory be within a cell (without crossing a cell boundary), but is generalized in the next section. It is worth noting that the derivation shown here results exactly the same charge and current deposition scheme as Villasenor & Buneman’s [9] (See Appendix. B). Cell crossings are generalized to ensure simultaneous charge and energy conservation (Sec. 3.4). Moreover, the derivation can be extended to second-order shape functions for charge density [32].

From the previous discussion, it is now clear that the shape function dyad introduced in Eqs. 19, 20 and 22 must be defined as:

[TABLE]

where $\mathbf{i}$, $\mathbf{j}$, $\mathbf{k}$ are Cartesian unit vectors.

3.4 Particle cell crossing

To ensure discrete charge and energy conservation, previous studies [11, 15] have advocated employing particle subcycling to deal with cell crossing. In this approach, each particle substep is performed with a CN step, which requires a nonlinear solve (e.g., by Picard iteration). While effective, this subcycling strategy can be expensive when multiple crossings occur in a single timestep (e.g., at cell corners), or when trajactory turning points occur near cell faces and the nonlinear iteration converges slowly.

Here, we introduce a simpler approach to avoid subcycling at cell crossings, which is critical for the overall efficiency of the algorithm. We begin with the assumption that the trajectory is a straight line during the (CFL-constrained) timestep $\Delta t$ (in this aspect, it is similar to that used in Ref. [9]). Thus, we rewrite Eq. 15 and 16 as:

[TABLE]

where the superscript $\nu$ denotes a sub-step. Here a sub-step is defined as a trajectory segment within a cell, and $\mathbf{a}_{p}^{\nu+\nicefrac{{1}}{{2}}}=\frac{q_{p}}{m_{p}}\sum_{h}\mathbf{E}_{h}^{n+\nicefrac{{1}}{{2}}}\cdot\bar{\bar{\mathbf{S}}}(\mathbf{x}_{h}-\mathbf{x}_{p}^{\nu+\nicefrac{{1}}{{2}}})$ is the acceleration during that segment of sub-step $\nu$ define by $[\mathbf{x}_{p}^{\nu},\mathbf{x}_{p}^{\nu+1}]$, evaluated at its center,

[TABLE]

The number of sub-steps $N_{\nu}$ is determined by the number of cell-crossings. If we define the trajectory segment to be $\Delta\mathbf{x}_{p}^{\nu}=\mathbf{x}_{p}^{\nu+1}-\mathbf{x}_{p}^{\nu}=\mathbf{v}_{p}^{n+\nicefrac{{1}}{{2}}}\Delta\tau^{\nu}$, then the sub-timestep $\Delta\tau^{\nu}$ is given by:

[TABLE]

where $\Delta\mathbf{x}_{p}=\mathbf{x}_{p}^{n+1}-\mathbf{x}_{p}^{n}$.

The definitions in Eqs. 33, 34 lead straightforwardly to discrete charge and energy conservation, as before [11, 15], with an orbit-averaged current density along the trajectory as:

[TABLE]

It is worth pointing out the key conditions for discrete energy and charge conservation are:

Identical shape functions are used for current deposition (Eq. 36) and electric field interpolation (Eq. 19); 2. 2.

The Lorentz force does no work on particles (Eq. 25). 3. 3.

The shape functions are evaluated at the center of each trajectory segment (Eq. 35).

3.5 Preservation of the Maxwell involutions

As stated earlier, Maxwell’s equations feature two involutions: the solenoidal character of the magnetic field and Gauss’ law. Regarding the former, we recall that the standard Yee scheme for Maxwell’s equations enforces the solenoidal property of magnetic field discretely, i.e.,

[TABLE]

as long as the initial magnetic field is divergence free. This is seen by taking the discrete divergence of Eq. 11, and noting that

[TABLE]

discretely by the Yee scheme. In practice, the solenoidal property is satisfied to numerical round-off levels.

Gauss’ law is enforced in our scheme by the exact charge conservation property of the particle update, as follows. Taking the divergence of Ampere’s law (Eq. 12) and utilizing Eq. 38, there results:

[TABLE]

The conservation of charge independently enforces:

[TABLE]

Combining these two equations, we find:

[TABLE]

which implies Gauss’ law at every timestep provided that it be satisfied in the beginning.

3.6 Iterative algorithm

The energy- and charge-conserving algorithm has been implemented in a code named iVPIC, based on the VPIC code developed at Los Alamos National Laboratory [33]. Our iterative implementation is outlined in Algorithm 1. Only Ampere’s equation needs iteration, as it couples the particles and the electric field nonlinearly. However, given that we respect the CFL for stability, the coupling is not stiff numerically, and a simple Picard iteration is sufficient for fast convergence. In practice, we observe convergence to numerical round-off (in single precision) in around 5 iterations. Note that particles are pushed once per iteration, at a cost comparable to an explicit particle push. Therefore, except for the additional storage needed to keep old-time particle quantities available during the iteration, the particle push in our implementation is competitive with a state-of-the-art explicit push.

4 Numerical results

We exercise iVPIC on some simple test problems that demonstrate the correctness of the implementation and the advertised conservation properties. In particular, we consider a relativistic two-stream instability, a Weibel instability, and a stimulated Brillouin scattering (SBS) of a laser in a plasma. All results are obtained by simulations using single (32-bit) precision IEEE floating-point rounding arithmetic, as what is originally employed in the VPIC code [33].

4.1 Two-stream instability

We consider a periodic relativistic system with two counter-streaming electron beams in an infinitely massive stationary ion background. The cold beam equations can be written as

[TABLE]

where the $n_{j}$ and $\mathbf{v}_{j}$ are the density and velocity of electron beam $j$, and $\mathbf{u}_{j}=\mathbf{v}_{j}\gamma_{j}$. Let $N_{j}$, $\mathbf{V}_{j}$, and $\tilde{n}_{j}$, $\tilde{\mathbf{v}}_{j}$ be the zeroth- and first-order electron density velocity for the $j^{th}$ beam, respectively. We assume the zeroth-order $\mathbf{E}=\mathbf{B}=0$, and denote $\tilde{\mathbf{E}}$ the first-order electric field. If the first-order quantities vary as $e^{i(\mathbf{k}\cdot\mathbf{r}-\omega t)}$, the linearized equations in the Fourier space may be written as

[TABLE]

where we have used the electrostatic approximation, i.e., the first-order magnetic field is zero. Note that the first-order relation between the proper velocity and velocity is

[TABLE]

where $\Gamma=(1-V^{2}/c^{2})^{-1/2}$. The linearized Poisson equation is

[TABLE]

The dispersion relation of the system is found to be

[TABLE]

where $\omega_{b}=\omega_{p}/\Gamma^{3}$, and $\omega_{p}=\sqrt{Nq/\epsilon_{0}m}$. Equation 43 has exactly the same form of the non-relativistic two-stream dispersion relation except for the $1/\Gamma^{3}$ factor. For the simplest case, i.e., two beams with opposite velocities and equal densities ($V_{1}=-V_{2}=V_{0}$, $\mathbf{k}\cdot\mathbf{V}=kV_{0}$, $N_{1}=N_{2}=N_{0}$), an analytical solution for the maximum growth rate exists:

[TABLE]

with $kV_{0}/\omega_{b}=\sqrt{3}/2$.

We set up simulations with $V_{0}$ varying from non-relativistic to relativistic regimes. The domain size is set to be $L_{x}=2\pi d_{e}$, where $d_{e}$ is electron skin depth, with $N_{x}=32$ cells, $N_{pc}=200$ particles per cell. We employ a time step $\Delta t=0.99\Delta x/c$. The boundary conditions are periodic. Figure 1 shows the simulation results. We see that, in the non-relativistic regime ($\Gamma\simeq 1$), the simulation linear growth rate matches well with both Eq. 44 and the non-relativistic one, given by:

[TABLE]

As $V_{0}$ increases, the relativistic lowering of the growth rate is more noticeable, as expected from the $\Gamma^{-3}$ scaling. For $\Gamma=1.39$, the relativistic growth rate decreases by a factor of $\sim 0.37$. In both weakly and strongly relativistic simulations, the change of the growth rate is very well captured.

4.2 Weibel instability

We use the Weibel instability to test the long time-scale conservation properties of the algorithm. The simulation is performed in a periodic 1D domain of $L_{x}=10d_{e}$ . We employ $N_{x}=64$ cells, $N_{pc}=200$ number of particles per cell, and a time step $\Delta t=0.99\Delta x/c$. Thermal temperatures of both electrons and ions are set to $v_{th,x}=0.1c$ and $v_{th,y,z}=0.3c$, respectively. The cell size is about 3 Debye lengths along $x$ at $t=0$. The mass ratio is set to $m_{i}/m_{e}=1836$. The boundary conditions are periodic.

We consider energy conservation first. Fig. 2 shows the total energy error history as a function of the number of timesteps (up to $10^{7}$) for both explicit and semi-implicit simulations with different number of nonlinear iterations. It is worth noting that the original VPIC algorithm uses the method of Villasenor and Buneman [9], which has been shown to have much better energy conservation properties than the standard momentum-conserving schemes [34]. Strictly speaking, the discrete conservation law proposed in this study ensures conservation of total energy with the magnetic energy defined by $W_{B}=\frac{1}{2}\int d\mathbf{x}\mathbf{B}^{n+1/2}\cdot\mathbf{B}^{n-1/2}$, which is not always well-posed (see Appendix A for detailed discussion about $W_{B}$). However, $W_{B}$ is almost always well-posed except in pathological cases, as long as the light wave mode of interest are well resolved and the CFL condition is respected. On the other hand, $\hat{W}_{B}=\frac{1}{2}\int d\mathbf{x}\mathbf{B}^{n}\cdot\mathbf{B}^{n}$ at integer timesteps, which can be obtained from Eq. 14, is always positive and well-posed. We have found similar long-term behavior of both definitions (i.e., conserving one measure leads to bounded errors in the other). For the semi-implicit scheme, a few iterations have to be performed to maintain small energy errors. If a single iteration is used, which is equivalent to a first-order forward Euler method, the corresponding energy error would grow even faster than the explicit one (which is second-order accurate). Energy conservation errors decrease rapidly with the number of iterations. Figure 2 shows the relative error of total energy defined at integer timesteps (employing $\hat{W}_{B}$) up to $10^{7}$ timesteps. For this test, 5 nonlinear iterations are sufficient to maintain conservation errors at relatively low levels ($\sim 10^{-4})$.

Regarding charge conservation, Fig. 3 shows the root-mean-square of Gauss’s law evaluated on the mesh with different flavors of the semi-implicit solver, all of them energy conserving, but differing in the particle-push treatment. In particular, we consider three types of particle pusher: a single particle push per field step, a subcycled particle push but without special treatment at cell boundaries (with two subcycling steps), and the proposed single-step charge-conserving particle push. As discussed earlier (Sec. 3.5), Gauss’s law should be always satisfied everywhere if it is satisfied at $t=0$. Clearly, for the non-charge-conserving single-step particle-push case, Gauss’s law is violated and the errors accumulate secularly with the number of time steps. Particle sub-cycling improves the long-term behavior of error accumulation, but the level of error increases very quickly at the early stage of the simulation and saturates at relatively high level. With the energy- and charge- conserving scheme, Gauss’s law is satisfied with high accuracy during the simulation, and errors remain at a very low level ($\thicksim 10^{-5}$) .

4.3 Laser-plasma instabilities

We use the iVPIC code with a plane wave source at intensity $1.25x10^{15}\>\mathrm{W}/\mathrm{cm}^{2}$to simulate stimulated Raman scattering (SRS) and stimulated Brillouin scattering (SBS) in a single laser speckle , in a under-dense plasma ($n_{e}=0.05n_{cr}$) where $n_{cr}$ is the critical density, similar to that of Ref. [35]. The simulation is performed in a quasi-1D domain of size $320\times 0.11$ ($d_{e}^{2}$), with grid points $N_{x}\times N_{z}=7776\times 2$, and the number of particles per cell $N_{pc}=512$. timestep $\Delta t=0.98\,\Delta t_{CFL}$, where $\Delta t_{CFL}=\Delta x/\sqrt{2}c$ in two dimensions. The laser is polarized along $\hat{y}$ and the simulation is in the $\hat{x}-\hat{z}$ domain. We initiate the laser pulse by gradually ramping up the intensity at the $x=0$ boundary, and the laser propagates into the plasma. We employ absorbing boundary conditions for fields, and absorbing boundary conditions in $x$ and periodic in $z$ for particles. The laser speckle is modeled as a Gaussian laser pulse polarized along the $y$ direction. The same setup is simulated with both VPIC and iVPIC. Figure 4 shows the backscatter reflectivity ($r=1-$measured Poynting flux/Laser Poynting flux) obtained at $x=0$ from the quasi-1D simulation as a function of time. Except for some differences in the fluctuations that are sensitive to the thermal noise and random number seeds used, the instantaneous and running time-averaged the reflectivities agree well between VPIC and iVPIC. The two simulations reproduce both the bursty nature of SRS in the electron trapping regime at early times (t < 7 ps) as well as the transition to SBS at late times (t > 7 ps). Although differences are observed in the timing and amplitude of the individual bursts (which are known to be sensitive to the properties of the noise inherent to particle simulations), excellent agreement is found in the levels of time-averaged reflectivity between VPIC and iVPIC as well as the timing of the transition from SRS to SBS.

5 Discussion and summary

We have developed the first low-dispersion, energy- and charge-conserving relativistic Vlasov-Maxwell PIC algorithm. It is semi-implicit, leap-frogging electromagnetic fields and time-centering the particle push. The light wave dispersion errors are the same as the standard explicit PIC [2], and lower than the previously proposed fully implicit scheme [14]. It has been known that the Yee scheme conserves the discrete energy of electromagnetic fields [36, 37]. We have extended such a principle to PIC by ensuring that the energy exchange between particles and fields is treated correctly. It requires nonlinear iteration between the electric field update in Ampere’s equation and the particle push, but the iteration is not stiff and converges very quickly with a simple Picard iterative scheme. To ensure simultaneous rigorous and automatic charge conservation, we have developed a novel particle push scheme that not only deals with particle cell crossings effectively without actually stopping the particle at cell boundaries, but also is exactly energy-conserving and directly revertible (i.e. can be solved explicitly [20]). The superior properties of the scheme have been demonstrated on three numerical examples of varying complexity and dimensionality, from electrostatic (two-stream) to fully electromagnetic (Weibel, SBS). At this point, without vectorization and careful optimization, each iteration in iVPIC is about 50% slower than an unvectorized explicit step, resulting in a slowdown of a factor of 6 when 4 iterations are taken. Future work will focus on the implementation of higher-order particle-mesh interpolations, and on the vectorization and optimization of the algorithm on modern architectures.

Appendix

Appendix A Well-posedness of magnetic energy definition

We address the question of whether the magnetic energy definition in Eq. 24 is well-posed or not, i.e., whether the energy is positive and provide a physically meaningful measure. We show that it is almost always well-posed when the CFL timestep is respected, and it is a second-order accurate approximation of the magnetic energy defined at integer time levels. We begin by expressing $\mathbf{B}^{n\pm 1/2}$ around time $n$ as:

[TABLE]

to find:

[TABLE]

where we denote $W_{B}=\frac{1}{2}\int d\mathbf{x}\mathbf{B}^{n+1/2}\cdot\mathbf{B}^{n-1/2}$, and $\hat{W}_{B}=\frac{1}{2}\int d\mathbf{x}(B^{n})^{2}$. This proves second-order accuracy. However, a priori, the measure $W_{B}$ may not be positive, while $\hat{W}_{B}$ is. For simplicity, we first show that $W_{B}$ is almost always well-posed in a 1D analysis.

Assuming that the magnetic field varies as $e^{-i\omega t+ikx}$, where $\omega=\omega(k)$,the discrete Fourier transform may be written as

[TABLE]

where $N_{g}$ is the number of grid points and $k_{l}=\frac{2\pi l}{L}$. Then, by applying the Fourier transform to $W_{B}$,

[TABLE]

where the superscript $*$ denotes complex conjugate, and $\Re[]$ denotes the real part of a complex number. We have used the orthogonal property of Fourier modes:

[TABLE]

From the 1D light wave dispersion relation:

[TABLE]

and using that $2\sin^{2}(\theta/2)=1-\cos(\theta)$, we find that:

[TABLE]

when $\frac{c\Delta t}{\Delta_{h}}\leq 1$ . Therefore:

[TABLE]

For the case of white noise (equal spectral content in all frequencies), it is easy to see that:

[TABLE]

where $A$ is a positive constant, and provided that $N_{g}$ is an even number. It follows that:

[TABLE]

i.e., the magnetic energy is positive definite if the CFL constraint $\frac{\Delta tc}{\Delta_{h}}<1$ is respected. In general, simulations will have varying spectral content with frequency, but low-frequency, well-resolved modes will carry the most information, and therefore the white noise analysis is conservative. For a mode with more than four grid points per wavelength, $\cos(k_{l}\Delta_{h})>0$. Therefore, as long as the field energy is not dominantly associated with high $k$ (i.e. $\left|k\Delta_{h}\right|>\pi/2$) modes, the field energy defined from $\mathbf{B}^{n+1/2}\cdot\mathbf{B}^{n-1/2}$ will be positive.

The above analysis can be straightforwardly extended to multiple dimensions. Assume a uniform grid, i.e., $\Delta x=\Delta y=\Delta z=\Delta_{h}$, and $\frac{\Delta tc}{\Delta_{h}}\lesssim\frac{1}{\sqrt{d}}$, where $d$ is the number of dimensions. Using the numerical dispersion relation of the light wave in 3D, i.e.,

[TABLE]

we find that:

[TABLE]

Clearly, the field energy given by Eq. 24 is well-posed if it is well-posed in each direction.

Appendix B Equivalence of Villasenor and Buneman’s charge conservation scheme

and the proposed one

The proposed current and charge density deposition scheme and cell-crossing scheme are exactly the same as Villasenor and Buneman’s [9]. The derivation presented in Sec. 3.3 appears to be new, however, and can be easily extended to higher-order splines [32]. Here we explicitly demonstrate the equivalence of the two schemes.

Assume that the particle trajectory is from $(x_{p}^{n},y_{p}^{n},z_{p}^{n})$ to $(x_{p}^{n+1},y_{p}^{n+1},z_{p}^{n+1})$, which lies within a single cell. We employed B-splines which are written as (letting $\Delta x=\Delta y=\Delta z=1$, and $x_{i}=y_{j}=z_{k}=0$, without loss of generality):

[TABLE]

and similarly in $y$ and $z$ directions.

We first show that Eq. 31 is trivally satisfied. Assuming that the particle is in a cell centered at ($i+\nicefrac{{1}}{{2}}$), we write the equation for the node ($i$) and $(i+1)$ respectively as

[TABLE]

where we have used Eq. 15. A formal derivation can be found in Ref. [11].

Adopting similar notations as in Ref. [9], we denote $\Delta x_{p}=x_{p}^{n+1}-x_{p}^{n}$, $\Delta y_{p}=y_{p}^{n+1}-y_{p}^{n}$, $\Delta z_{p}=z_{p}^{n+1}-z_{p}^{n}$, and $\bar{\xi}=(x_{p}^{n+1}+x_{p}^{n})/2$, $\bar{\eta}=(y_{p}^{n+1}+y_{p}^{n})/2$, $\bar{\zeta}=(z_{p}^{n+1}+z_{p}^{n})/2$, Eqs. 27-30 for the node $(i+1,j+1,k+1)$ are written as

[TABLE]

and their substitutions into Eq. 26 yield

[TABLE]

which is exactly the same as Eq. 38 of Ref. [9] except for the subscript $p$ adopted in this study.

For trajectories across cell boundaries, the treatment is exactly the same as Villasenor and Buneman’s scheme, i.e., we split the trajectory into segments, each of which lies in one cell, and the above treatment applies.