A semi-implicit, energy- and charge-conserving particle-in-cell algorithm for the relativistic Vlasov-Maxwell equations
Guangye Chen, Luis Chacón, Lin Yin, Brian J. Albright, David J. Stark, Robert F. Bird
aa r X i v : . [ phy s i c s . c o m p - ph ] M a r A semi-implicit, energy- and charge-conserving particle-in-cell algorithmfor the relativistic Vlasov-Maxwell equations
G. Chen ∗ , L. Chacón, L. Yin, B. J. Albright, D. J. Stark, R. F. Bird Los Alamos National Laboratory, Los Alamos, NM 87545
Abstract
Conventional explicit electromagnetic particle-in-cell (PIC) algorithms do not conserve discrete energy ex-actly. Time-centered fully implicit PIC algorithms can conserve discrete energy exactly, but may introducelarge dispersion errors in the light-wave modes. This can lead to intolerable simulation errors where accu-rate light propagation is needed (e.g. in laser-plasma interactions). In this study, we selectively combinethe leap-frog and Crank-Nicolson methods to produce an exactly energy- and charge-conserving relativisticelectromagnetic PIC algorithm. Specifically, we employ the leap-frog method for Maxwell’s equations, andthe Crank-Nicolson method for the particle equations. The semi-implicit algorithm admits exact global en-ergy conservation, exact local charge conservation, and preserves the dispersion properties of the leap-frogmethod for the light wave. The algorithm employs a new particle pusher designed to maximize efficiencyand minimize wall-clock-time impact vs. the explicit alternative. It has been implemented in a code namediVPIC, based on the Los Alamos National Laboratory VPIC code ( 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.
Keywords: particle-in-cell, energy conservation, charge conservation, laser plasma interactions
PACS:
1. Introduction
Particle-in-cell (PIC) methods combine Lagrangian and Eulerian techniques for solving the relativisticVlasov-Maxwell equations (or a subset thereof) for kinetic plasma simulations. Specifically, the method ofcharacteristics 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 forevolving Maxwell’s equations on a mesh. The Yee scheme [1] is the most popular algorithm for integratingMaxwell’s equations, due to its simplicity and second-order accuracy. B-splines are widely used for inter-polations between the particles and the mesh. A detailed description of classical PIC methods and theiranalysis can be found in Ref. [2, 3] . ∗ Corresponding author
Email address: [email protected] (G. Chen)
Preprint submitted to Elsevier March 12, 2019 espite their numerical simplicity, PIC algorithms have been successfully applied to a large range ofplasma 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 ap-plicability of PIC for long-term multiscale simulations may be currently limited by the trustworthiness ofthe 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 ex-actly (in practice to numerical roundoff) in a discrete system (with finite timestep ∆ t and mesh size ∆ x ).Discrete charge conservation is satisfied by most PIC algorithms employed today [9, 10, 11, 12]. Discreteenergy conservation has been recently demonstrated with advances in implicit PIC algorithms, either fullyimplicit [11, 13, 14, 15, 16], or semi-implicit [17] (note that an early version of “energy-conserving” algo-rithm by Lewis [18] conserves discrete energy exactly only in the limit of ∆ t → ∼
2. The classical electromagnetic relativistic PIC method
We consider the relativistic Vlasov-Maxwell system: ∂ t f s + v · ∇ f s + q s m s ( E + v × B ) · ∇ u f s =
0, (1) ∂ B ∂ t + ∇ × E =
0, (2) ǫ ∂ E ∂ t − µ ∇ × B + j =
0. (3) ∇ · B =
0, (4) ∇ · E = ρ , (5)where f s ( x , u ) is the particle distribution function of species s in phase space, x denotes physical position, v denotes velocity, u = v γ is the proper velocity, and γ = √ + u / c is the Lorentz factor, with c the speedof light. In these equations, ǫ and µ are the vacuum permittivity and permeability respectively, E and B are the electric and magnetic field, respectively, q s and m s are the species charge and mass respectively, and j and ρ are current and charge densities, respectively, found from the distribution function as: j = ∑ s q s Z d u v f s ( r , u ) , ρ = ∑ s q s Z d u f s ( r , u ) .Equations 4 and 5 are the two involutions of Maxwell’s equations [21]. Note that conventional FDTDmethods 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 chargeconservation ∂ρ∂ t + ∇ · j = f s ( x , u ) (i.e., macro-particles): f s ( x , u ) ≈ ∑ p ∈ s w p δ ( x − x p ) δ ( u − u p ) ,where w p is the particle weight (constant in collisionless PIC simulations). The particle equations of motionread: d x p dt = v p , (7) d u p dt = q p m p ( E p + v p × B p ) . (8)3ssuming the use of a Yee mesh, the classical leap-frog scheme is written as: u n + / p − u n − / p ∆ t = q p m p E np + u n + / p + u n − / p γ np × B np ! , (9) x n + p − x np ∆ t = u n + / p γ n + / , (10) B n + / h − B n − / h ∆ t = −∇ h × E nh , (11) ǫ E n + h − E nh ∆ t = µ ∇ h × B n + / h − j n + / h , (12)where the superscript n denotes time level, the subscript p denotes a particle quantity or a field evaluated atthe particle position, ∆ t is the timestep, and the subscript h = ( i , j , k ) denotes mesh quantities and operators(using Yee finite differences). Here γ n = q + ( u n − / + q ∆ t m E n ) , and γ n + / = p + ( u n + / ) . Onepotential problem with Boris scheme is that it does not preserve the correct limit in a force-free field [26].The current density j n + / h is gathered from particles using B-spline interpolation. This scheme features aCourant-Fredrichs-Lewy (CFL) condition [27] for stability, and constrains cell sizes to be comparable to theDebye 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 particleequations. Energy conservation is achieved by advancing particles and the electric field synchronously, af-fording the scheme a semi-implicit character and demanding iteration. This iteration, however, is simple toimplement, and of rapid convergence. Simultaneously, automatic local discrete charge conservation is en-forced by the proper choice of shape interpolation functions, and a new treatment of particle cell crossingsthat does not require particles to actually stop at cell boundaries, as previously proposed in Refs. [11, 15].
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 dis-persion errors. It is well-known that the light-wave dispersion relation can be almost perfectly preservedin 1D with the Yee scheme as c ∆ t / ∆ x →
1. In multiple dimensions, the numerical dispersion varies withmodal wavelength, propagation direction, and spatial discretization [22], but it remains tolerable in practicefor many applications. In contrast, for arbitrary ∆ 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 ∆ x → π [30]. Charged particles with speed greater than the phase speed of light waves would lose muchof their energy via Cerenkov radiation [31], and may generate significant noise in the simulation. For thisreason, we keep leap-frog for the field time advance. For the analysis below, we rewrite Eq. 11 in the4ollowing equivalent way: B n + / h − B nh ∆ t /2 = −∇ h × E nh , (13) B n + h − B n + / h ∆ t /2 = −∇ h × E n + h . (14)Particles are advanced fully implicitly by the CN scheme, yielding: x n + p − x np ∆ t = v n + / p , (15) u n + p − u np ∆ t = q p m p ( E n + / p + ¯ v n + / p × B n + / p ) , (16)where we define v n + / p = u n + p + u np γ n + p + γ np , (17)¯ v n + / p = u n + p + u np s + (cid:18) u n + p + u np c (cid:19) . (18)Note that this particle pusher is similar to that in Ref. [14] in using v n + / p in the particle position update,but it employs ¯ v n + / p 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: E n + / p = ∑ h E n + / h · ¯¯ S ( x h − x p ) , (19) B n + / p = ∑ h B n + / h · ¯¯ S ( x h − x p ) . (20)Note that the field E h and B h at half timestep are obtained differently: E n + / h = E n + h + E nh B n + / h is obtained from Eq. 13. Here, x h is the grid location, x p is typically chosen to be at the center ofparticle trajectory [11], and as before the subscript h = ( i , j , k ) denotes the grid index. The specific form ofthe shape function dyad ¯¯ S and Eq. 17 are important to ensure energy and charge conservation, and will bediscussed 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 currentdensity: j n + / h = ∆ h ∑ p q p v n + / p · ¯¯ S ( x h − x p ) . (22)where ∆ 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-centeredupdate of particle positions and velocities, together with electric field, results in a coupled field-particle5ystem: x p is a function of the new-time electric field through Eqs. 15 and 16, which in turn determinesthe current that determines the field via Eq. 12. This, in turn, will require an iterative solve, which will beintroduced later (Sec. 3.6).We prove the energy conservation theorem next. Discrete energy conservation can be readily shown as follows. Multiplying Eq. 12 by E n + / h and inte-grating over space, we find: ∑ h ∆ h (cid:20) ǫ ( E n + h − E nh ) · E n + / h − ∆ t µ ( ∇ h × B n + / h ) · E n + / h + ∆ t j n + / h · E n + / h (cid:21) =
0. (23)The first term in Eq. 23 gives the change of electric energy: ǫ ∑ h ∆ h ( E n + h − E nh ) · E n + / h = ǫ ∑ h ∆ h h ( E n + h ) − ( E nh ) i ≡ W n + E − W nE ,using Eq. 21. The second term gives the change of magnetic energy: − ∆ t µ ∑ h ∆ h ( ∇ h × B n + / h ) · E n + / h = − ∆ t µ ∑ h ∆ h ( ∇ h × E n + h + E nh ) · B n + / h = µ ∑ h ∆ h h ( B n + h − B nh ) · B n + / h i = W n + B − W nB ,where we have used discrete integration by parts, Eqs. 14 and 13, and we have defined the magnetic energyas: W nB ≡ µ ∑ h ∆ h B n + / h · B n − / h , (24)using that B nh = ( B n + / h + B n − / h ) /2 (see Eq. 13). This definition is non-standard, and as we show inAppendix. 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: ∆ t ∑ h ∆ h j n + / h · E n + / h = ∆ t ∑ h E n + / h · ∑ p q p v n + / p S ( x p − x h ) = ∑ p m p v n + / p · ( u n + p − u np )= ∑ p m p c ( γ n + p − γ np ) ≡ W n + p − W np ,where we have used Eqs. 16, 17, 18, and that v n + / p · ¯ v n + / p × B n + / p =
0. (25)The energy conservation theorem sought follows: (cid:0) W E + W B + W p (cid:1)(cid:12)(cid:12) n + n = .3. Charge conservation Discrete local charge conservation requires ∂ t ρ + ∇ · j = ( i , j , k ) . Because the chargeconservation equation is linear, it is sufficient to enforce this constraint on the contributions of each particle: ( ρ p ) n + ijk − ( ρ p ) nijk ∆ t + ∇ h · j n + / p (cid:12)(cid:12)(cid:12) ijk =
0. (26)We employ first-order (trilinear) splines for the charge density: ρ np , ijk = q p S ( x i − x np ) S ( y j − y np ) S ( z k − z np ) , (27)and mixed zeroth- and first-order splines (indicated by subscripts 0 and 1 respectively) for the currentdensity: j n + / x , p , i + / , j , k = q p v n + / x S ( x i + / − x n + / p ) S n + / jk ( y p , z p ) , (28) j n + / y , p , i , j + / , k = q p v n + / y S ( y j + / − y n + / p ) S n + / ik ( z p , x p ) , (29) j n + / z , p , i , j , k + / = q p v n + / z S ( z k + / − z n + / p ) S n + / ij ( x p , y p ) , (30)where, for instance S n + / jk ( y p , z p ) ≡ " S ( y j − y n + p ) S ( z k − z n + p ) + S ( y j − y np ) S ( z k − z n + p ) + S ( y j − y n + p ) S ( z k − z np ) + S ( y j − y np ) S ( z k − z np ) ,and so on. The proof of charge conservation is as follows. We first decompose the density change into threeone-dimensional shifts (letting q p = ρ n + ijk − ρ nijk = S ( x i − x n + p ) S ( y j − y n + p ) S ( z k − z n + p ) − S ( x i − x np ) S ( y j − y np ) S ( z k − z np )= h S ( x i − x n + p ) − S ( x i − x np ) i S n + / jk ( y p , z p ) + h S ( y j − y n + p ) − S ( y j − y np ) i S n + / ik ( z p , x p ) + h S ( z k − z n + p ) − S ( z k − z np ) i S n + / ij ( x p , y p ) .Along each direction, using the particle orbit position equation, it can be shown that (see Appendix B andRef. [11]): S ( x i − x n + p ) − S ( x i − x np ) ∆ t + v n + / x S ( x i + / − x n + / p ) − S ( x i − / − x n + / p ) ∆ x =
0, (31)where x n + / p = x n + p + x np y and ˆ z directions. Equation 26 then follows. The proof requires that the particletrajectory be within a cell (without crossing a cell boundary), but is generalized in the next section. It is7orth noting that the derivation shown here results exactly the same charge and current deposition schemeas Villasenor & Buneman’s [9] (See Appendix. B). Cell crossings are generalized to ensure simultaneouscharge and energy conservation (Sec. 3.4). Moreover, the derivation can be extended to second-order shapefunctions for charge density [32].From the previous discussion, it is now clear that the shape function dyad introduced in Eqs. 19, 20 and22 must be defined as: ¯¯ S ( x ijk − x p ) = i ⊗ i S ( x i − x p ) S n + / jk ( y p , z p )+ j ⊗ j S ( y j − y p ) S n + / ik ( z p , x p )+ k ⊗ k S ( z k − z p ) S n + / ij ( x p , y p ) ,where i , j , k are Cartesian unit vectors. To ensure discrete charge and energy conservation, previous studies [11, 15] have advocated employingparticle subcycling to deal with cell crossing. In this approach, each particle substep is performed with a CNstep, which requires a nonlinear solve (e.g., by Picard iteration). While effective, this subcycling strategycan be expensive when multiple crossings occur in a single timestep (e.g., at cell corners), or when trajactoryturning 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 theoverall efficiency of the algorithm. We begin with the assumption that the trajectory is a straight line duringthe (CFL-constrained) timestep ∆ t (in this aspect, it is similar to that used in Ref. [9]). Thus, we rewrite Eq.15 and 16 as: x n + p − x np = v n + / p ∆ t , (33) u n + p − u np = N ν − ∑ ν = a ν + / p ∆ τ ν p , (34)where the superscript ν denotes a sub-step. Here a sub-step is defined as a trajectory segment within a cell,and a ν + / p = q p m p ∑ h E n + / h · ¯¯ S ( x h − x ν + / p ) is the acceleration during that segment of sub-step ν define by [ x ν p , x ν + p ] , evaluated at its center, x ν + / p = x ν + p + x ν p N ν is determined by the number of cell-crossings. If we define the trajectorysegment to be ∆ x ν p = x ν + p − x ν p = v n + / p ∆ τ ν , then the sub-timestep ∆ τ ν is given by: ∆ τ ν p = | ∆ x ν p || ∆ x p | ∆ t ,where ∆ x p = x n + p − x np . 8he definitions in Eqs. 33, 34 lead straightforwardly to discrete charge and energy conservation, asbefore [11, 15], with an orbit-averaged current density along the trajectory as:¯ j h = ∆ h ∆ t ∑ p q p N ν − ∑ ν = ∆ x ν p · ¯¯ S ( x h − x ν + / p ) . (36)It is worth pointing out the key conditions for discrete energy and charge conservation are:1. Identical shape functions are used for current deposition (Eq. 36) and electric field interpolation (Eq.19);2. The Lorentz force does no work on particles (Eq. 25).3. The shape functions are evaluated at the center of each trajectory segment (Eq. 35). As stated earlier, Maxwell’s equations feature two involutions: the solenoidal character of the magneticfield and Gauss’ law. Regarding the former, we recall that the standard Yee scheme for Maxwell’s equationsenforces the solenoidal property of magnetic field discretely, i.e., ∇ h · B =
0, (37)as long as the initial magnetic field is divergence free. This is seen by taking the discrete divergence of Eq.11, and noting that ∇ h · ∇ h × B =
0, (38)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: ǫ ∇ h · E n + h − ∇ h · E nh ∆ t + ∇ h · j n + / h = ∇ h · j n + / h = − ρ n + h − ρ nh ∆ t .Combining these two equations, we find: ∇ h · E n + h − ∇ h · E nh = ρ n + h − ρ nh ,which implies Gauss’ law at every timestep provided that it be satisfied in the beginning. The energy- and charge-conserving algorithm has been implemented in a code named iVPIC, based onthe VPIC code developed at Los Alamos National Laboratory [33]. Our iterative implementation is outlinedin Algorithm 1. Only Ampere’s equation needs iteration, as it couples the particles and the electric field9 lgorithm 1
Iterative solution of the semi-implicit conservative Vlasov-Maxwell PIC equations1. Start from state B n − / h , E nh , and n x np , v np o .2. Update magnetic field to B n + / h according to Eq. 11.3. Solve the coupled Ampere-particle system (Eqs. 12, 15, 16) for E n + h , x n + p , and v n + p by the followingiterative procedure ( k is the iteration index; lagged electric field coupling is indicated with a box).Starting with E n + = = E nh . For each iteration : x n + kp − x np = v n + / , kp ∆ t , u n + kp − u np = ∑ ν q p m p E n + kp + E np + ¯ v n + / , kp × B n + / p ∆ τ ν ,¯ j h = ∆ h ∆ t ∑ p q p v n + / , kp · ∑ ν ¯¯ S ( x h − x ν + / p ) ∆ τ ν , ǫ E n + k + h − E nh ∆ t = µ ∇ h × B n + / h − ¯ j h ,where v n + / , kp = u n + kp + u np γ n + kp + γ np , and ¯ v n + / , kp = u n + kp + u np s + (cid:18) u n + kp + u np c (cid:19) .Increment k .The number of iterations is currently specified to a fixed number.4. Update timestep counter and return to 1. 10onlinearly. However, given that we respect the CFL for stability, the coupling is not stiff numerically, and asimple Picard iteration is sufficient for fast convergence. In practice, we observe convergence to numericalround-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 tokeep old-time particle quantities available during the iteration, the particle push in our implementation iscompetitive 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 implementationand the advertised conservation properties. In particular, we consider a relativistic two-stream instability, aWeibel instability, and a stimulated Brillouin scattering (SBS) of a laser in a plasma. All results are obtainedby simulations using single (32-bit) precision IEEE floating-point rounding arithmetic, as what is originallyemployed in the VPIC code [33].
We consider a periodic relativistic system with two counter-streaming electron beams in an infinitelymassive stationary ion background. The cold beam equations can be written as ∂ n j ∂ t + ∇ · ( n j v j ) = ∂ u j ∂ t + ( v j · ∇ ) u j = qm ( E + v j × B ) ,where the n j and v j are the density and velocity of electron beam j , and u j = v j γ j . Let N j , V j , and ˜ n j , ˜ v j be the zeroth- and first-order electron density velocity for the j th beam, respectively. We assume the zeroth-order E = B =
0, and denote ˜ E the first-order electric field. If the first-order quantities vary as e i ( k · r − ω t ) , thelinearized equations in the Fourier space may be written as − ω ˜ n j + N j k · ˜ v j + k · V j ˜ n j =
0, (39) − ω ˜ u j + k · V j ˜ u j = − i qm ˜ E , (40)where we have used the electrostatic approximation, i.e., the first-order magnetic field is zero. Note thatthe first-order relation between the proper velocity and velocity is ˜u j = ˜ v j Γ j , (41)where Γ = ( − V / c ) − . The linearized Poisson equation is ǫ k · ˜E = ∑ j q j ˜ n j . (42)The dispersion relation of the system is found to be ∑ j = ω bj ( ω − k · V j ) =
1, (43)11 .000000010.000000100.000001000.000010000.000100000.00100000 0 5 10 15 20 25 W E ω p t( γ =1.0038) ivpiclinear theorynon-relativistic linear theory W E ω p t( γ =1.066) ivpiclinear theorynon-relativistic linear theory W E ω p t( γ =1.39) ivpiclinear theorynon-relativistic linear theory Figure 1: Time history of electric field energy ( W E ) of relativistic two-stream instability simulations compared with linear theories. where ω b = ω p / Γ , and ω p = p Nq / ǫ m . Equation 43 has exactly the same form of the non-relativistictwo-stream dispersion relation except for the 1/ Γ factor. For the simplest case, i.e., two beams with oppo-site velocities and equal densities ( V = − V = V , k · V = kV , N = N = N ), an analytical solution forthe maximum growth rate exists: Im [ ω ] = ω b /2, (44)with kV / ω b = √ V varying from non-relativistic to relativistic regimes. The domain size isset to be L x = π d e , where d e is electron skin depth, with N x =
32 cells, N pc =
200 particles per cell. Weemploy a time step ∆ t = ∆ x / c . The boundary conditions are periodic. Figure 1 shows the simulationresults. We see that, in the non-relativistic regime ( Γ ≃ [ ω nr ] = ω p /2. (45)As V increases, the relativistic lowering of the growth rate is more noticeable, as expected from the Γ − scaling. For Γ = ∼ We use the Weibel instability to test the long time-scale conservation properties of the algorithm. Thesimulation is performed in a periodic 1D domain of L x = d e . We employ N x =
64 cells, N pc = ∆ t = ∆ x / c . Thermal temperatures of both electrons andions are set to v th , x = c and v th , y , z = c , respectively. The cell size is about 3 Debye lengths along x at t =
0. The mass ratio is set to m i / m e = ) for both explicit and semi-implicit simulations with different number ofnonlinear iterations. It is worth noting that the original VPIC algorithm uses the method of Villasenorand Buneman [9], which has been shown to have much better energy conservation properties than thestandard momentum-conserving schemes [34]. Strictly speaking, the discrete conservation law proposed12 R e l a t i v e E ne r g y E rr o r ) steps explicit4 iterations5 iterations Figure 2: History of total energy conservation errors in long-term electron Weibel instability simulations. Energy conservation im-proves with increasing number of iterations. The insert shows the relative energy error of semi-implicit simulation with 5 iterations. in this study ensures conservation of total energy with the magnetic energy defined by W B = R d xB n + · B n − , 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 wellresolved and the CFL condition is respected. On the other hand, ˆ W B = R d xB n · B n at integer timesteps,which can be obtained from Eq. 14, is always positive and well-posed. We have found similar long-termbehavior of both definitions (i.e., conserving one measure leads to bounded errors in the other). For thesemi-implicit scheme, a few iterations have to be performed to maintain small energy errors. If a singleiteration is used, which is equivalent to a first-order forward Euler method, the corresponding energy errorwould grow even faster than the explicit one (which is second-order accurate). Energy conservation errorsdecrease rapidly with the number of iterations. Figure 2 shows the relative error of total energy defined atinteger timesteps (employing ˆ W B ) up to 10 timesteps. For this test, 5 nonlinear iterations are sufficient tomaintain conservation errors at relatively low levels ( ∼ − ) .Regarding charge conservation, Fig. 3 shows the root-mean-square of Gauss’s law evaluated on themesh with different flavors of the semi-implicit solver, all of them energy conserving, but differing in theparticle-push treatment. In particular, we consider three types of particle pusher: a single particle pushper field step, a subcycled particle push but without special treatment at cell boundaries (with two subcy-cling 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 with13 R M S ( ∇ • E - ρ ) ) steps energy-conserving, no particle sub-cyclingenergy-conserving sub-cyclingenergy-charge-conserving sub-cycling -5
0 10
Figure 3: History of errors of Gauss’s law in long-term electron Weibel instability simulations. The insert shows the root-mean-squareerror of Gauss’s equation an semi-implicit simulation using the proposed energy- and charge-conserving algorithm. the number of time steps. Particle sub-cycling improves the long-term behavior of error accumulation, butthe level of error increases very quickly at the early stage of the simulation and saturates at relatively highlevel. With the energy- and charge- conserving scheme, Gauss’s law is satisfied with high accuracy duringthe simulation, and errors remain at a very low level ( ∼ − ) . We use the iVPIC code with a plane wave source at intensity 1.25 x W/cm to simulate stimulatedRaman scattering (SRS) and stimulated Brillouin scattering (SBS) in a single laser speckle , in a under-denseplasma ( n e = n cr ) where n cr is the critical density, similar to that of Ref. [35]. The simulation is per-formed in a quasi-1D domain of size 320 × d e ), with grid points N x × N z = ×
2, and the numberof particles per cell N pc = ∆ t = ∆ t CFL , where ∆ t CFL = ∆ x / √ c in two dimensions. Thelaser is polarized along ˆ y and the simulation is in the ˆ x − ˆ z domain. We initiate the laser pulse by graduallyramping up the intensity at the x = x and periodic in z forparticles. The laser speckle is modeled as a Gaussian laser pulse polarized along the y direction. The samesetup is simulated with both VPIC and iVPIC. Figure 4 shows the backscatter reflectivity ( r = − measuredPoynting flux/Laser Poynting flux) obtained at x = R e f l e c t i v i t y time (ps) vpicivpic 0 0.1 0.2 0.3 0.4 0.5 0.6 0 5 10 15 20 25 30 35 40 45 50 R e f l e c t i v i t y time (ps) vpicivpic Figure 4: History of reflectivity in simulations with both VPIC and iVPIC (using 5 field-particle iterations). Left: instantaneousreflectivity averaged over laser period. Right: running time-averaged reflectivity (by T R T r ( t ) dt ). times (t < 7 ps) as well as the transition to SBS at late times (t > 7 ps). Although differences are observedin the timing and amplitude of the individual bursts (which are known to be sensitive to the propertiesof the noise inherent to particle simulations), excellent agreement is found in the levels of time-averagedreflectivity 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-MaxwellPIC 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 previouslyproposed fully implicit scheme [14]. It has been known that the Yee scheme conserves the discrete energyof electromagnetic fields [36, 37]. We have extended such a principle to PIC by ensuring that the energy ex-change between particles and fields is treated correctly. It requires nonlinear iteration between the electricfield update in Ampere’s equation and the particle push, but the iteration is not stiff and converges veryquickly with a simple Picard iterative scheme. To ensure simultaneous rigorous and automatic charge con-servation, we have developed a novel particle push scheme that not only deals with particle cell crossingseffectively without actually stopping the particle at cell boundaries, but also is exactly energy-conservingand directly revertible (i.e. can be solved explicitly [20]). The superior properties of the scheme have beendemonstrated 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 op-timization, each iteration in iVPIC is about 50% slower than an unvectorized explicit step, resulting in aslowdown of a factor of 6 when 4 iterations are taken. Future work will focus on the implementation ofhigher-order particle-mesh interpolations, and on the vectorization and optimization of the algorithm onmodern architectures. 15 ppendix
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 almostalways well-posed when the CFL timestep is respected, and it is a second-order accurate approximation ofthe magnetic energy defined at integer time levels. We begin by expressing B n ± around time n as: B n ± = B n ± ∆ t ∂ t B n , (A.1)to find: W B = ˆ W B − ∆ t Z d x ( ∂ t B n ) · ( ∂ t B n ) , (A.2)where we denote W B = R d xB n + · B n − , and ˆ W B = R d x ( B n ) . This proves second-order accuracy.However, a priori , the measure W B may not be positive, while ˆ 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 ω t + ikx , where ω = ω ( k ) ,the discrete Fourier transformmay be written as B ( x , t ) = N g N g − ∑ l = B ( ω l , k l ) e − i ω l t + ik l x , (A.3)where N g is the number of grid points and k l = π lL . Then, by applying the Fourier transform to W B , N g − ∑ h = B n + / ( x h ) · B n − / ( x h ) ∗ = N g N g − ∑ h = N g − ∑ l = N g − ∑ m = B ( ω l , k l ) · B ( ω k , k m ) ∗ e i ( k l − k m ) x h ℜ h e − i ( ω l t n + − ω m t n − ) i = N g N g − ∑ l = N g − ∑ m = B ( ω l , k l ) · B ( ω k , k m ) ∗ ℜ h e − i ( ω l t n + − ω m t n − ) i N g − ∑ h = e i ( k l − k m ) x h = N g N g − ∑ l = B ( ω l , k l ) · B ( ω l , k l ) ∗ cos ( ω l ∆ t ) , (A.4)where the superscript ∗ denotes complex conjugate, and ℜ [] denotes the real part of a complex number. Wehave used the orthogonal property of Fourier modes: N g − ∑ h = e i ( k l − k m ) x h = N g δ lm . (A.5)From the 1D light wave dispersion relation: (cid:20) ∆ h c ∆ t (cid:21) sin (cid:18) ω ∆ t (cid:19) = sin (cid:18) k x ∆ h (cid:19) , (A.6)and using that 2 sin ( θ /2 ) = − cos ( θ ) , we find that:cos ( ω ∆ t ) ≥ (cid:20) c ∆ t ∆ h (cid:21) cos( k x ∆ h ), (A.7)16hen c ∆ t ∆ h ≤ N g − ∑ h = B n + / ( x h ) · B n − / ( x h ) ∗ ≥ N g (cid:20) c ∆ t ∆ h (cid:21) N g − ∑ l = B ( ω l , k l ) · B ( ω l , k l ) ∗ cos( k l ∆ h ). (A.8)For the case of white noise (equal spectral content in all frequencies), it is easy to see that:1 N g N g − ∑ l = B ( ω l , k l ) · B ( ω l , k l ) ∗ cos ( k l ∆ h ) = A N g − ∑ l = cos ( k l ∆ h ) =
0, (A.9)where A is a positive constant, and provided that N g is an even number. It follows that: N g − ∑ h = B n + / ( x h ) · B n − / ( x h ) ∗ ≥
0, (A.10)i.e., the magnetic energy is positive definite if the CFL constraint ∆ tc ∆ h < ( k l ∆ h ) >
0. Therefore, as long as the field energy is not dominantlyassociated with high k (i.e. | k ∆ h | > π /2) modes, the field energy defined from B n + · B n − will bepositive.The above analysis can be straightforwardly extended to multiple dimensions. Assume a uniform grid,i.e., ∆ x = ∆ y = ∆ z = ∆ h , and ∆ tc ∆ h . √ d , where d is the number of dimensions. Using the numericaldispersion relation of the light wave in 3D, i.e., (cid:20) ∆ h c ∆ t (cid:21) sin (cid:18) ω ∆ t (cid:19) = sin (cid:18) k x ∆ h (cid:19) + sin (cid:18) k y ∆ h (cid:19) + sin (cid:18) k z ∆ h (cid:19) , (A.11)we find that: cos ( ω ∆ t ) ≥ (cid:20) c ∆ t ∆ h (cid:21) (cid:2) cos ( k x ∆ h ) + cos ( k y ∆ h ) + cos ( k z ∆ h ) (cid:3) . (A.12)Clearly, the field energy given by Eq. 24 is well-posed if it is well-posed in each direction. 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 thesame 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 ofthe two schemes.Assume that the particle trajectory is from ( x np , y np , z np ) to ( x n + p , y n + p , z n + p ) , which lies within a singlecell. We employed B-splines which are written as (letting ∆ x = ∆ y = ∆ z =
1, and x i = y j = z k =
0, withoutloss of generality): S ( x i + / − x p ) =
1, (B.1) S ( x i − x p ) = − x p , (B.2) S ( x i + − x p ) = x p , (B.3)17nd 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 + / ),we write the equation for the node ( i ) and ( i + ) respectively aslhs = ( − x n + p ) − ( − x np ) ∆ t + v n + / x − =
0, (B.4)lhs = x n + p − x np ∆ t + v n + / x − =
0, (B.5)where we have used Eq. 15. A formal derivation can be found in Ref. [11].Adopting similar notations as in Ref. [9], we denote ∆ x p = x n + p − x np , ∆ y p = y n + p − y np , ∆ z p = z n + p − z np ,and ¯ ξ = ( x n + p + x np ) /2, ¯ η = ( y n + p + y np ) /2, ¯ ζ = ( z n + p + z np ) /2, Eqs. 27-30 for the node ( i + j + k + ) are written as ρ n + p , i + j + k + = q p ( ¯ ξ + ∆ x p )( ¯ η + ∆ y p )( ¯ ζ + ∆ z p ) , (B.6) ρ np , i + j + k + = q p ( ¯ ξ − ∆ x p )( ¯ η − ∆ y p )( ¯ ζ − ∆ z p ) , (B.7) j n + / x , p , i + / , j + k + = q p ( ∆ x p ∆ t ¯ η ¯ ζ + ∆ y p ∆ z p ) , (B.8) j n + / y , p , i + j + / , k + = q p ( ∆ y p ∆ t ¯ ζ ¯ ξ + ∆ z p ∆ x p ) , (B.9) j n + / z , p , i + j + k + / = q p ( ∆ z p ∆ t ¯ ξ ¯ η + ∆ x p ∆ y p ) , (B.10) j n + / x , p , i + / , j + k + =
0, (B.11) j n + / y , p , i + j + / , k + =
0, (B.12) j n + / z , p , i + j + k + / =
0, (B.13)and their substitutions into Eq. 26 yield ∆ x p ¯ η ¯ ζ + ∆ y p ¯ ζ ¯ ξ + ∆ z p ¯ ξ ¯ η = ( ¯ ξ + ∆ x p )( ¯ η + ∆ y p )( ¯ ζ + ∆ z p ) − ( ¯ ξ − ∆ x p )( ¯ η − ∆ y p )( ¯ ζ − ∆ z p ) , (B.14)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’sscheme, i.e., we split the trajectory into segments, each of which lies in one cell, and the above treatmentapplies. 18 eferences [1] K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations inisotropic media,” IEEE Transactions on antennas and propagation , vol. 14, no. 3, pp. 302–307, 1966.[2] C. K. Birdsall and A. B. Langdon,
Plasma physics via computer simulation . CRC Press, 2004.[3] R. W. Hockney and J. W. Eastwood,
Computer simulation using particles . CRC Press, 1988.[4] J. Dawson, “One-dimensional plasma model,”
The Physics of Fluids , vol. 5, no. 4, pp. 445–459, 1962.[5] A. Hasegawa and H. Okuda, “One-dimensional plasma model in the presence of a magnetic field,”
The Physics of Fluids , vol. 11, no. 9, pp. 1995–2003, 1968.[6] W. Lee, “Gyrokinetic approach in particle simulation,”
The Physics of Fluids , vol. 26, no. 2, pp. 556–562,1983.[7] C. Huang, V. K. Decyk, C. Ren, M. Zhou, W. Lu, W. B. Mori, J. H. Cooley, T. M. Antonsen Jr, andT. Katsouleas, “Quickpic: A highly efficient particle-in-cell code for modeling wakefield accelerationin plasmas,”
Journal of Computational Physics , vol. 217, no. 2, pp. 658–679, 2006.[8] A. B. Langdon and B. F. Lasinski, “Electromagnetic and relativistic plasma simulation models,”
Meth-ods in Computational Physics , vol. 16, pp. 327–366, 1976.[9] J. Villasenor and O. Buneman, “Rigorous charge conservation for local electromagnetic field solvers,”
Computer Physics Communications , vol. 69, no. 2-3, pp. 306–316, 1992.[10] T. Z. Esirkepov, “Exact charge conservation scheme for particle-in-cell simulation with an arbitraryform-factor,”
Computer Physics Communications , vol. 135, no. 2, pp. 144–153, 2001.[11] G. Chen, L. Chacón, and D. C. Barnes, “An energy-and charge-conserving, implicit, electrostaticparticle-in-cell algorithm,”
Journal of Computational Physics , vol. 230, no. 18, pp. 7018–7036, 2011.[12] M. C. Pinto, S. Jund, S. Salmon, and E. Sonnendrücker, “Charge-conserving FEM–PIC schemes ongeneral grids,”
Comptes Rendus Mecanique , vol. 342, no. 10-11, pp. 570–582, 2014.[13] S. Markidis and G. Lapenta, “The energy conserving particle-in-cell method,”
Journal of ComputationalPhysics , vol. 230, no. 18, pp. 7037–7052, 2011.[14] G. Lapenta and S. Markidis, “Particle acceleration and energy conservation in particle in cell simula-tions,”
Physics of Plasmas , vol. 18, no. 7, p. 072101, 2011.[15] G. Chen and L. Chacón, “A multi-dimensional, energy-and charge-conserving, nonlinearly implicit,electromagnetic Vlasov–Darwin particle-in-cell algorithm,”
Computer Physics Communications , vol. 197,pp. 73–87, 2015. 1916] L. Chacón and G. Chen, “A curvilinear, fully implicit, conservative electromagnetic PIC algorithm inmultiple dimensions,”
Journal of Computational Physics , vol. 316, pp. 578–597, 2016.[17] G. Lapenta, “Exactly energy conserving semi-implicit particle in cell formulation,”
Journal of Computa-tional Physics , vol. 334, pp. 349–366, 2017.[18] H. R. Lewis, “Energy-conserving numerical approximations for Vlasov plasmas,”
Journal of Computa-tional Physics , vol. 6, no. 1, pp. 136–141, 1970.[19] Y. Chen and G. Toth, “Gauss’s law satisfying energy-conserving semi-implicit particle-in-cell method,” arXiv preprint arXiv:1808.05745 , 2018.[20] A. V. Higuera and J. R. Cary, “Structure-preserving second-order integration of relativistic chargedparticle trajectories in electromagnetic fields,”
Physics of Plasmas , vol. 24, no. 5, p. 052104, 2017.[21] W. M. Seiler,
Involution: The formal theory of differential equations and its applications in computer algebra ,vol. 24. Springer Science & Business Media, 2009.[22] A. Taflove and S. C. Hagness,
Computational electrodynamics: the finite-difference time-domain method .Artech house, 2005.[23] B. Marder, “A method for incorporating Gauss’ law into electromagnetic pic codes,”
Journal of Compu-tational Physics , vol. 68, no. 1, pp. 48–55, 1987.[24] A. B. Langdon, “On enforcing Gauss’ law in electromagnetic particle-in-cell codes,”
Computer PhysicsCommunications , vol. 70, no. 3, pp. 447–450, 1992.[25] C.-D. Munz, P. Omnes, R. Schneider, E. Sonnendrücker, and U. Voss, “Divergence correction tech-niques for Maxwell solvers based on a hyperbolic model,”
Journal of Computational Physics , vol. 161,no. 2, pp. 484–511, 2000.[26] J.-L. Vay, “Simulation of beams or plasmas crossing at relativistic velocity,”
Physics of Plasmas , vol. 15,no. 5, p. 056701, 2008.[27] R. Courant, K. Friedrichs, and H. Lewy, “On the partial difference equations of mathematical physics,”
IBM journal of Research and Development , vol. 11, no. 2, pp. 215–234, 1967.[28] C. Sun and C. Trueman, “Unconditionally stable crank-nicolson scheme for solving two-dimensionalMaxwell’s equations,”
Electronics Letters , vol. 39, no. 7, pp. 595–597, 2003.[29] G. Sun and C. Trueman, “Unconditionally-stable FDTD method based on Crank-Nicolson scheme forsolving three-dimensional Maxwell equations,”
Electronics Letters , vol. 40, no. 10, pp. 589–590, 2004.2030] G. Chen and L. Chacón, “An energy-and charge-conserving, nonlinearly implicit, electromagnetic1d-3v vlasov–darwin particle-in-cell algorithm,”
Computer Physics Communications , vol. 185, no. 10,pp. 2391–2402, 2014.[31] J. D. Jackson,
Classical electrodynamics . John Wiley & Sons, 2012.[32] A. Stanier, L. Chacón, and G. Chen, “A fully implicit, conservative, non-linear, electromagnetic hybridparticle-ion/fluid-electron algorithm,”
Journal of Computational Physics , vol. 376, pp. 597–616, 2019.[33] K. J. Bowers, B. Albright, L. Yin, B. Bergen, and T. Kwan, “Ultrahigh performance three-dimensionalelectromagnetic relativistic kinetic plasma simulation,”
Physics of Plasmas , vol. 15, no. 5, p. 055703,2008.[34] A. Pukhov, “Three-dimensional electromagnetic relativistic particle-in-cell code VLPL (Virtual LaserPlasma Lab),”
Journal of Plasma Physics , vol. 61, no. 3, pp. 425–433, 1999.[35] B. Albright, L. Yin, K. Bowers, and B. Bergen, “Multi-dimensional dynamics of stimulated brillouinscattering in a laser speckle: Ion acoustic wave bowing, breakup, and laser-seeded two-ion-wave de-cay,”
Physics of Plasmas , vol. 23, no. 3, p. 032703, 2016.[36] R. Schuhmann and T. Weiland, “Conservation of discrete energy and related laws in the finite integra-tion technique,”
Progress In Electromagnetics Research , vol. 32, pp. 301–316, 2001.[37] F. Edelvik, R. Schuhmann, and T. Weiland, “A general stability analysis of FIT/FDTD applied to lossydielectrics and lumped elements,”