Relativistic extension of a charge-conservative finite element solver for time-dependent Maxwell-Vlasov equations

Particle-in-cell (PIC) algorithms^3–7 have been a very successful tool in many scientific and engineering applications such as electron accelerators,^8–10 laser-plasma interactions,^6,11–14 astrophysics,^15,16 vacuum electronic devices,^17–19 and semiconductor devices.^20–23 In many cases, the particles of interest are often in the relativistic regime and the relevant physical phenomena need to be described by taking into account fully relativistic effects. For example, runaway electrons can cause plasma discharge disruptions in fusion devices. The interaction of runaway electron beams with plasma turbulence requires full electromagnetic treatment in relativistic regimes. Relativistic plasma waves propagating in energetic electron-positron (pair) plasmas have also been of interest in astrophysics (pulsar atmospheres). Relativistic PIC algorithms can be found in a variety of Refs. 9, 11–14, 20, and 24–27.

Historically, PIC simulations have been carried out using finite-difference field solvers based on regular (structured) grids, such as the finite-difference time-domain (FDTD) algorithm. However, regular structured grids are less suited for representing arbitrary geometries with slanted or curved boundaries, leading to a “staircased” representation of such geometries. To enable modeling of more complex geometries, locally conformal finite-difference schemes are sometimes employed.^28,29 Although some geometric flexibility is gained by this strategy, fundamental challenges remain when employing locally conformal schemes such as the lack of stable and systematic mesh refinement strategies. In addition, new challenges are introduced, for instance, the lack of energy conservation and slowing down of simulations due to a more stringent stability condition on the time step.

It is widely recognized that arbitrary geometries are best represented by irregular (unstructured) meshes. The finite element method (FEM) is naturally suited for such meshes, but for many years, optimal integration of FEM solvers in PIC algorithms was inhibited by violation of the charge continuity equation in discretizations based on irregular meshes. The lack of charge conservation leads to incorrect modeling of the underlying physics and biased electron trajectories. There have been a number of available strategies to enforce charge-conservation on irregular meshes such as projection methods or hyperbolic cleaning techniques; however, these approaches are ad-hoc and not entirely satisfactory as they require a time-consuming global Poisson solver at each time step or introduce a correction potential that may subtly alter the physics.

In recent years, a first-principles solution to this problem has been achieved by a number of works which share in common a representation of dynamical variables on the mesh (electromagnetic fields, currents, and charges) by Whitney forms^1,30–33 and a consistent transfer of information between the particle positions/velocities and the mesh variables (during the scatter/gather steps of the PIC algorithm) by means of Whitney forms as well. In the exterior calculus framework,^1,31,34,35 Whitney forms can be viewed as a consistent discrete representation of the differential forms of various degrees representing the dynamical variables on irregular meshes. This has paved a way for the consistent integration of PIC algorithms with full-wave finite element time-domain (FETD) field solvers² on such meshes.

In this paper, a charge conserving FETD PIC algorithm previously developed for time-dependent Maxwell-Vlasov equations on irregular meshes^1,2 is extended to the relativistic regime. In particular, we integrate Boris,^36,57 Vay,³⁷ and Higuera–Cary³⁸ relativistic pushers in the conservative PIC-FETD algorithm for solving time-dependent Maxwell-Vlasov equations and provide a brief comparison among them. Several examples such as particle cyclotron motion, harmonic particle oscillation in the Lorentz-boosted frame, and relativistic Bernstein modes in magnetized charge-neutral (pair) plasmas are presented for validation. We adopt MKS units throughout this work.

We consider Maxwell's equations

coupled to the Vlasov equation governing the phase space distribution function $f(r,v,t)$ of collisionless particles

where q and m₀ are the particle charge and rest mass, respectively, $γ−2=1−v2/c2$⁠, with c being the speed of light, and we have assumed a single species for simplicity.

The Maxwell-Vlasov system can be efficiently solved numerically using particle-in-cell (PIC) algorithms whereby $f(r,v,t)$ is represented by a collection of superparticles (i.e., a coarse-graining of the phase space distribution). As is well known, time-dependent PIC algorithms consist of four procedures at each time step: the field update incorporating Maxwell's equations, the gather step transferring information from the mesh to the particle positions, the particle push incorporating the Lorentz force law and Newton's equation of motion, and the scatter step transferring particle dynamics information to the mesh. As Maxwell's equations are already relativistic, extension of the FETD PIC algorithm to the relativistic regime should focus on the particle push. Because of this, we will focus below primarily on this aspect of the algorithm. For completeness, we also describe the other three steps as they integrate into the full algorithm, albeit more briefly. Further details about the field solver, scatter, and gather steps can be found in Refs. 1 and 2.

On an irregular mesh, the electric field intensity E(r,t) and the magnetic flux density B(r,t) can be expanded by using vector proxies of Whitney forms as^39–43 

where N_n, N_e, and N_f are the total number of free nodes, edges, and faces in the mesh and $Wi1(r)$ and $Wi2(r)$ are the vector proxies of Whitney 1- and 2-forms^39,41 associated (1:1) with edges and faces of the mesh, respectively.⁴⁴ The electric current density and charge density on the mesh can be likewise expressed as

where $Wi0$ is a Whitney 0-form associated with the grid nodes. In what follows, the dynamical degrees of freedom $ei(t), bi(t), ii(t)$⁠, and $qi(t)$ are grouped into column vectors denoted as e, b, i, and q. By applying the generalized Stokes' theorem and the Galerkin method to obtain a spatial discretization of Maxwell's equations and by performing a leap-frog discretization in time, the field update algorithm for Eqs. (1) and (2) is obtained as^1,2,43

with the discrete version of the divergence constraints $∇·B=0$ and $∇·ϵ0E=ρ$ being automatically² fulfilled for all time steps, i.e.,

where the superscript n denotes the time-step index, and C and S are the incidence matrices representing the discrete exterior derivative operator or, equivalently, the discrete curl and discrete divergence operators^34,35,45,46 distilled from their metric structure. Due to their metric-free character of the incidence matrices, all their elements are integers in the set ${−1,0,1}$⁠. The matrices C and S refer to the primal mesh, while $C̃$ and $S̃$ refer to the dual mesh.^34,35,45,46 It can be shown^45,46 that $C̃=CT$ and $S̃=ST$⁠. In addition, $[⋆μ−1]$ and $[⋆ϵ]$ in Eq. (9) are diagonally dominant symmetric positive definite matrices representing the discrete Hodge star operator,^35,41,43 which encodes the metric aspects of the mesh. It should be noted that the inverse matrix $[⋆ϵ]−1$ present in Eq. (9) is not computed directly because this would lead to a dense system update. To avoid the need for a linear solve at every time step, a sparse approximate inverse (SPAI) is precomputed in a parallel fashion and with tunable accuracy. The accuracy is controlled by a parameter associated with the sparsity level of the approximate inverse, which yields exponential convergence to the exact inverse (in the sense of the Frobenius norm).^2,43 Details of the parallel SPAI implementation can be found in Ref. 43.

In the gather step, the field values are interpolated using the corresponding Whitney forms in Eqs. (4) and (5) computed at particle positions, i.e.,

where $rp$ is the position of the p-th particle.

In the particle update step, the particle mass is modified to account for relativistic effects such that

where $up=γpvp, vp$ is the velocity of the p-th particle, and γ_p is its relativistic factor defined as $γp−2=1−|vp|2/c2$⁠. Using the central-differences to approximate the time derivatives, Eqs. (14) and (15) are discretized as

where $v¯p$ is the mean particle velocity between the $n±12$ time steps, which can also be approximated as $u¯p/γ¯p$ with $u¯p=γ¯pv¯p$⁠.

In the non-relativistic case, $γp→1, v¯p$ can be chosen based on the midpoint rule, viz., $v¯pn=vpn=(vpn+12+vpn−12)/2$⁠, to obtain updated phase coordinates explicitly. In this case, the (non-relativistic) Boris algorithm is typically used not only due to its computationally efficient velocity update obtained by separating irrotational (electric) and rotational (magnetic) forces but also because of its long-term numerical stability. The latter property essentially means that, in spite of not being symplectic, the non-relativistic Boris algorithm preserves phase-space volume such that it provides energy conservation bounded within a finite interval. Note that every symplectic integrator guarantees phase-space volume-preservation but not vice-versa. In contrast, in the relativistic regime, $v¯p$ should be carefully determined to accurately model the kinetics of high-energy particles. Next, we examine in detail three different relativistic pushers proposed by Boris, Vay, and Higuera-Cary.

The main tenet of the relativistic Boris pusher is basically similar to the non-relativistic-Boris-pusher, viz., separation of irrotational and rotational forces.⁷ Importantly, it averages $v¯p$ as

where $v(u)=u/1+|u|2/c2, ϵpn=αEpn$⁠, and $α=qΔt/2m0$⁠. The particle velocity update in the relativistic Boris pusher follows the procedure below⁷ 

where $uB−,uB′$⁠, and $uB+$ are the auxiliary vectors and the subscript B refers to the Boris algorithm. In addition, $tB=βpn/γ¯p,B, sB=2tB/(1+|tB|2)$⁠, and $βpn=αBpn$⁠. The factor $γ¯p,B$ is computed as

and to obtain $Bpn$⁠, we set

Note that the relativistic Boris pusher has two variants: with and without correction. The relativistic Boris pusher without correction uses $tB$ as defined above. The relativistic Boris pusher with correction uses $tB=(βp/|βp|)tan (|βpn|/γ¯p,B)$ instead.

The separation of two different forces can be easily observed by substituting Eqs. (19) and (22) into Eq. (17), which results in

In Eq. (25), the effect of $Epn$ is completely removed, so that only magnetic rotation is left.

The relativistic Boris pusher preserves volumes in the phase-space because the determinant of Jacobian for time-update map $ψBn:(rpn,upn−12)→(rpn+1,upn+12)$ equals one.^36,38 To verify that, we express determinant of the Jacobian for $ψBn$ as

where from (16) we have that $∂rpn+1/∂rpn=I¯¯+∂upn+12/∂xpn$⁠, with $I¯¯$ being the 3 × 3 identity matrix. If we assume electromagnetic fields to be uniform along the particle trajectory during one time step, $∂upn+12/∂xpn=0$⁠. In addition, it is clear that $∂xpn+1/∂xpn=I¯¯$⁠. Substituting $upn+12$ from (17) into (16) and taking a derivative w.r.t. $∂upn−12$ of the resulting equation, we obtain $∂xpn+1/∂upn−12=Δt(∂upn+12/∂upn−12)$⁠. As a result, the determinant in (26) simply takes the form $|∂upn+12/∂upn−12|$⁠. This derivative can be computed by splitting one time-update into two half-time-updates and evaluating each serially as follows:

But since³⁸ 

it follows that $|∂upn+12/∂upn−12|=1$⁠, and therefore, the relativistic-Boris-pusher is phase-space volume-preserving. This means that energy conservation is attained for long time simulations except for residual errors stemming from numerical artifacts such as finite-precision roundoff. However, the main disadvantage of the relativistic-Boris-pusher is that it cannot accurately capture the particle acceleration by electric forces. This is because magnetic rotation fails to consider the varying relativistic factor due to electric field effects during the magnetic rotation.

The Vay pusher corrects trajectories of relativistic particles experiencing electric fields by averaging $v¯p$ as

The particle velocity update follows the procedure below³⁷ 

where $tV=βpn/γpn+12, sV=2tV/(1+|tV|2), uV*=uV′·βpn/c$⁠, and $γV′=1+|uV′|′2/c2$⁠. The Vay pusher correctly models energetic particle motion under electric forces with Lorentz (relativistic) invariance. In other words, the relation between the particle's trajectory observed in (relativistic) moving and laboratory frames satisfies the Lorentz transformation.

The Higuera–Cary pusher provides an accurate treatment of electric forces by approximating the average velocity as³⁸ 

The particle velocity update is basically similar to Boris algorithm except for the relativistic factor as

Therefore, the Higuera-Cary pusher is also phase-space volume-preserving.

In this step, node-based charge density values and edge-based current density values are determined on the mesh from the positions and velocities of the particles. In order to achieve exact charge conservation on the mesh, the charge density is obtained by evaluating node-based Whitney 0-forms at each particle position and the current density is obtained by evaluating (integrating) edge-based Whitney 1-forms along each particle trajectory.^1,2 This correspondence is consistent with the discrete version of the continuity equation and applies to both non-relativistic and relativistic regimes equally.

The charge Q of the p-th particle is distributed (scattered) to nearby nodes using the Whitney 0-form such that

where the subscript i indicates the nodal index and $Wi0$ is the Whitney 0-form, equivalent to barycentric coordinates or point $rp$ with respect to node i, denoted as $λi(rp)$⁠. Likewise, the current along an edge ij (indexed here by its two end nodes i and j) produced by the movement of the p-th particle with charge Q from $rp,s$ to $rp,f$ during $Δt$ is determined by evaluating the line integral of the Whitney 1-form $Wij1(rp)$ (associated with the edge ij) along a path from $rp,s$ to $rp,f$⁠, i.e.,

where we have used the shorthand notation $λis=λi(rp,s),λif=λi(rp,f)$⁠, and similarly for λ_j. It should be stressed that Eq. (38) is an exact expression and no numerical quadratures are necessary to evaluate the trajectory integral.

In order to verify charge conservation, we examine the discrete continuity equation

which is an equality relating node-based quantities and it should remain valid on every node of the mesh. By considering an arbitrary node labeled as node i, the first term in Eq. (39) can be expressed as

On the other hand, assuming for simplicity that two edges are connected to node i, labeled as ik and ij (the generalization for more edges is straightforward), the second term in Eq. (39) can be expanded as

where we have used Eq. (38) and $λi(rp)+λj(rp)+λk(rp)=1$⁠, a basic property of the barycentric coordinates on a triangle with nodes i, j, and k. Since the sum of Eqs. (40) and (41) is equal to zero, the continuity equation is satisfied. Gauss' law can be also be verified in a geometrical fashion that is illustrated in Ref. 1.

To validate our relativistic PIC algorithm, we first examine a synchrocyclotron example. As electrons usually have velocities near to the speed of light in this case, progressively more energy needs to be delivered to accelerate them due to relativistic effects. The relativistic mass increase results in a lower orbital (cyclotron) frequency. Therefore, the driving RF electric field should have variable frequencies matching this relativistic cyclotron frequency. Figure 1(b) illustrates a computational mesh used for the simulation where the centripetal force from external magnets is present in the red region leading to a circular motion of electrons and a longitudinal RF electric field is present in the vertical blue strip leading to a periodic electron acceleration.

Figure 2(a) shows electron cyclotron motion in the non-relativistic regime, where the relativistic factor is assumed to be one. An electron is injected at (x,y) = (0.52, 0.5) m with an initial velocity of $|v0|=v0=1×107$ m/s. The static magnetic force is determined to be $Bz=m0v0/qr=2.84281×10−3$ T for an initial orbital radius r of 0.02 m and the RF electric force is set to be E_x = 2 × 10⁵ V/m. The thickness of the vertical strip in which the longitudinal electric field is present is 0.02 m. As can be seen, the spacing of two adjacent orbits becomes successively smaller due to the increasing velocity. Figure 2(b) shows the trajectory of the electron with the same initial conditions where we use the PIC algorithm with the relativistic Boris pusher with correction discussed in Sec. II C 1. In this case, the frequency of the RF electric force is set to be constant (79.6 MHz), which results in an unsynchronized phase between particle velocity and electric force and a mixed trajectory. In Fig. 2(c), the frequency of the RF electric field is matched to the synchrocyclotron frequency given by $f=qB/2πγm0$⁠. This frequency is shown as a function of the number of time steps in Fig. 3. Therefore, the in-phase acceleration is maintained at all times and a circular trajectory is observed at higher energies. Note that the total distance over which an electron moves in this case is shorter than that for the non-relativistic case because of its smaller speed caused by the relativistic mass. Figure 4 shows electron velocity magnitudes in the three cases. It is clearly observed that the deceleration occurs near 4000 time steps for the second case. Also, the third case shows slightly smaller magnitudes than the first one due to the relativistic mass.

Tables I–III provide a verification of Gauss' law. The amount of charge on arbitrarily selected mesh nodes is recorded at different time steps. As the rightmost columns in these tables show, the normalized residual due to the discrete version of Gauss' law is near the double precision floor (<10⁻¹⁵).

In order to compare how accurately the three different kinds of relativistic particle pushers capture relativistic E × B drift motions, we consider a harmonic oscillatory motion of a positron in the Lorentz-boosted frame with γ_f = 2 such as in Ref. 37. Initial parameters of the harmonic motion are transformed via the Lorentz transformation into the moving frame along $ŷ$⁠, and PIC simulations are performed in the moving frame. At the end of the simulation, we re-transform the phase coordinates from the moving frame back to the laboratory frame by using the inverse Lorentz transformation. We compare the resultant trajectories obtained with three different particle pushers and analytical predictions in Fig. 5. As we discussed in Sec. III C, the relativistic Boris pusher (without or with correction) cannot correctly capture relativistic E × B motion; on the other hand, results obtained with the Vay pusher and Higuera-Cary pusher accurately match the analytical prediction.

The non-relativistic electron Bernstein mode^47,48 is a purely electrostatic plasma wave propagating normal to a stationary magnetic field. It has been mainly explored in magnetic plasma confinement fusion as a promising alternative to conventional electron cyclotron electromagnetic waves such as the ordinary (O) or extraordinary (X) modes which have frequency cutoffs associated with plasma density.⁴⁹ The Bernstein mode is free from the density cut-off, and as a result, it is able to reach the core of over-dense plasmas in Tokamak devices and heat the plasma electrons effectively. It is well known that conventional non-relativistic Bernstein waves are present at harmonics of electron cyclotron resonances.⁴⁸ Figure 6 shows the dispersion relation in terms of the normalized frequency $ω̂=ω/ωc$ and the normalized transverse wavenumber $k̂⊥=kxc/ωc$ of non-relativistic electron Bernstein and X modes (see the colormap) and compares PIC simulation results with analytical predictions (dashed red line). In this case, $ωc=9.0×1011 rad/s,B0=5.13ẑ T,ωp=8.7×1011, n0=2.4×1020 m−3$⁠, and the initial isotropic speed distribution (equilibrium state) obeys a Maxwellian distribution with v_th = 0.07c. The simulation parameters are chosen similar to Ref. 50.

Recently, analytical works have been done to characterize the behaviors of Bernstein modes in relativistic electron-positron pair-plasmas.^51–54 There are several features in this case distinguishing the classical wave from the relativistic one: (1) The classical Maxwellian distribution (equilibrium state) is modified to the Maxwell-Boltzmann-Jüttner distribution (relativistic Maxwellian), (2) the mobility of positively charged particles is identical to that of negatively charged particles, and (3) the conventional dispersion relations are significantly transformed to undamped or damped closed curve shapes. In this example, we use our FETD PIC algorithm to perform simulations of Bernstein modes propagating in relativistic magnetized pair-plasmas and compare the results with analytical predictions.

To derive analytical dispersion relations of magnetized plasma waves,^48,51–54 we obtain a complex permittivity tensor, $ϵ¯¯$⁠, associated with plasma currents. First of all, we consider a small perturbation imposed on equilibrium magnetized pair-plasmas with parameters of

where f_s is a distribution function represented in the phase space for the species s, E is the electric field intensity, B is the magnetic flux density, and subscriptions of 0 and 1 denote the equilibrium and perturbed quantities, respectively. Note that the perturbed electromagnetic fields are proportional to $ei(k·r−ωt)$⁠. Equilibrium relativistic electron-positron pair-plasmas are typically described by Maxwell-Boltzmann-Jüttner (relativistic Maxwellian) distribution which is given by

where $η=m0ckBT$⁠, k_B is the Boltzmann constant, T is the kinetic temperature, $K2(·)$ is the modified Bessel function of the second kind, and $γ=(1+p2c2)−1/2$⁠. The evolution of the distribution function is governed by the Vlasov equation. Its first-order approximation takes the form

Substituting Eq. (45) into Eq. (46) and integrating Eq. (46) over time, solutions for the perturbed distribution function, f_1,s, can be obtained. Then, plasma currents are calculated based on f_1,s as

where $σ¯¯$ is a conductivity tensor from which we obtain the complex permittivity associated with plasma currents as

We are interested in longitudinal electrostatic plasma waves propagating in the x-direction. Thus, $(ω,kx)$ curves yielding zeros of ϵ_xx form the dispersion relations for Bernstein modes in a magnetized relativistic pair-plasma. The expression for ϵ_xx can be written as⁵¹ 

where $F2 3(12,1;32,1−a,1+a;−b2)$ is a hypergeometric function defined as

$ω̂p=ωp/ωc, p̂=p/(m0c) β̂=k̂⊥p̂$⁠, and $Jν(·)$ denotes the Bessel function of the first kind for ν. One can find details in Ref. 51 on how to numerically compute the integral in Eq. (49), which exhibits singularities at harmonics of the (rest) cyclotron frequencies.

We consider the case of $ω̂p=3$ and η = 20. Other parameters are specified as follows: $B0=5ẑ T,ωc=8.7941×1011 rad/s,ωp=2.6382×1012 rad/s,ne=2.1870×1021 m−3$⁠, and electron or positron density, $n0=2.1870×1021$⁠. The Debye length, λ_D, is equal to 2.55 × 10⁻⁵m, and the characteristic (relativistic) gyroradius, r_g, becomes 7.92 × 10⁻⁵m. An irregular mesh with triangular elements of size l_x × l_y is constructed with l_y = r_g and l_x = 1000 × r_g, and the average mesh element size is comparable to λ_D. The number of nodes, edges, and faces in the mesh is 7252, 18 123, and 10 872, respectively. The left and right boundaries of the mesh are terminated by perfectly matched layers (PMLs)^55,56 to mimic open boundaries. Periodic boundary conditions (PBCs) are applied at the top and bottom boundaries. In order to obtain the dispersion relation in (k_x, ω) for Bernstein waves propagating along x, we spatially sample the electric field along the x direction at each time-step and perform a Fourier transform on the resulting dataset in (x, t). The PML treatment of the left and right boundaries not only reduces unwanted reflections but also avoids aliasing effects in the dispersion relation caused by PBC with low sampling rates and the presence of image sources. Whenever particles meet a PBC boundary wall, they are removed and reassigned the corresponding relative positions on the other PBC boundary wall and the same momentum. The total number of superparticles for e and p species representing 1.7222 × 10¹⁰ electrons and positrons, respectively, is set to $Nsp,e+Nsp,p=4×105$⁠. Note that $Nsp,e$ and $Nsp,p$ are identical, and our simulation is initialized so that the average number of superparticles per grid element is around 40. This number is chosen as an attempt to provide a good trade-off between simulation speed and the rate of plasma self-heating. Initially, superparticles are uniformly distributed on the mesh in a pairwise fashion, i.e., each species-e superparticle is collocated with another species-p superparticle. This arrangement produces zero initial electric field. Based on the initial Maxwell-Boltzmann-Jüttner distribution [Fig. 7(a)], superparticles for species e are then launched in random 2D directions. The corresponding p superparticles are simultaneously launched with the same speed in the opposite direction [see Fig. 7(b)]. Figure 7(a) compares the Maxwell-Boltzmann-Jüttner distribution and classical Maxwellian distribution. It can be seen that the classical Maxwellian velocity distribution starts to gradually deviate from the relativistic one above the most probable velocity, which is about v/c = 0.21. We simulate up to 100 000 of time-steps and employ Δ t = 0.01 ps, which is equivalent to a Courant factor of 0.2. Then, we perform space and time Fourier transforms of sampled data to obtain the dispersion relations.

Figure 8 illustrates the dispersion relations for relativistic Bernstein waves with η = 1/20, compared to analytical predictions. It is observed in Fig. 8 that there are solutions of curved shape between every two neighboring harmonics of the (rest) cyclotron frequency. Our PIC simulations capture this feature quite well, which distinguishes the relativistic Bernstein wave from the classical one as predicted by theory. It is also interesting to note that every upper curve of each solution is considerably weaker since, as pointed out in Ref. 53, the damping coefficient for the upper curve is larger than that for the lower curve. In addition, there are two kinds of stationary modes at around $ω̂=0.9$ and $ω̂=4.1$⁠. As shown in Ref. 51, these stationary modes are nearly dependent on ω_p and the gap and location frequency between two stationary modes becomes smaller as η decreases (ultra-relativistic). It should be noted that even though PIC simulations are somewhat noisy, some physical damping is certainly present (especially in warm plasmas and for high $k̂x$⁠). Introducing a spectral filtering scheme to reduce the noise may lead to heating and artifacts. It is highly desirable to compare numerical results with analytical predictions here since although there are many PIC codes available, there is little understanding of how well PIC codes describe classical plasma effects, like Bernstein modes. It is expected that the most efficient way to reduce the numerical noise is to use high-order particle gather/scatter procedures. We leave these issues to a future work.

In order to check charge conservation, Fig. 9 illustrates normalized residuals for the discrete continuity equation (DCE) [Fig. 9(a)] and the discrete Gauss law (DGL) [Fig. 9(b)] across all mesh nodes at three different time-steps, n = 10 000, 20 000, and 30 000. The normalized residuals for DCE and DGL at the k-th node, NR_DCE and NR_DGL, are defined as

for $k=1,2,…,Nn$⁠. It is observed in Figs. 9(a) and 9(b) that the normalized residuals are near the double precision floor, which again indicates that no spurious charges are deposited at the nodes.

A finite element time-domain particle-in-cell algorithm is presented for the relativistic Maxwell-Vlasov equations. The key feature of the algorithm consists in combining a new charge-conserving scatter/gather scheme for irregular meshes with relativistic particle pushers for efficient plasma simulations. Three different relativistic particle pushers are compared and several numerical examples are provided for illustrative purposes.

This work was supported by NSF under Grant No. ECCS-1305838, by OSC under Grant Nos. PAS-0110 and PAS-0061, and by the OSU Presidential Fellowship program.