Time Integrator Agnostic Charge Conserving Finite Element PIC
Scott O'Connor, Zane D. Crawford, O. H. Ramachandran, John Luginsland, B. Shanker
TTime Integrator Agnostic Charge Conserving Finite Element PIC
Scott O’Connor, a) Zane D. Crawford, b) O. H. Ramachandran, b) John Luginsland, and B. Shanker Department of Electrical and Computer Engineering, Michigan State University, East Lansing,MI
Developing particle-in-cell (PIC) methods using finite element basis sets, and without auxiliary divergencecleaning methods, was a long standing problem until recently. It was shown that if consistent spatial basisfunctions are used, one can indeed create a methodology that was charge conserving, albeit using a leap-frogtime stepping method. While this is a significant advance, leap frog schemes are only conditionally stableand time step sizes are closely tied to the underlying mesh. Ideally, to take full advantage of advances infinite element methods (FEMs), one needs a charge conserving PIC methodology that is agnostic to the timestepping method. This is the principal contribution of this paper. In what follows, we shall develop thismethodology, prove that both charge and Gauss’ laws are discretely satisfied at every time step, provide thenecessary details to implement this methodology for both the wave equation FEM and Maxwell Solver FEM,and finally demonstrate its efficacy on a suite of test problems. The method will be demonstrated by singleparticle evolution, non-neutral beams with space-charge, and adiabatic expansion of a neutral plasma, wherethe debye length has been resolved, and real mass ratios are used.
I. INTRODUCTION
Simulation of space charge and plasmas is critical to anumber of areas in science and engineering. These rangefrom, applications of pulsed power to particle accelera-tors to satellites and medicine . The means to do so haslargely relied on Particle-in-cell (PIC) methods. PIC hasbeen around since the 1950’s and is a popular methods ofmodeling plasma and space charge due to its simplicityand ease of use . PIC enables a self consistent solution toMaxwell’s equation and equations of motion for chargedspecies. Traditionally, PIC is based on finite differencetime domain to evolve fields . The use of regular cubicalgrids presents challenges, especially in modeling complexgeometry. Modeling curved features requires small cellsizes, and this results in a stair-stepped approximation ofthe desired geometry as well as small time steps in keep-ing with the Courant–Friedrichs–Lewy condition. Usingcut-cells has improved the geometry representation byallowing boundaries to cut across cells . Complex andfine features, as well as multi-scale objects, require theuse of a prohibitively expensive number of small cells forhigh fidelity simulations. As a result of these challenges,there has been persistent investigation into the use ofmore sophisticated field evolution techniques . A nat-ural choice is using time domain finite-element method(TDFEM) due to (a) unconditionally stable time step-ping methods, (b) ability to model complex geometries,and (c) well developed extensions to higher order (bothin representation of fields and geometry) .While TDFEM can be thought of as a panacea formodeling complex geometries, it is not so for crucialquantities that must be conserved. These include Gauss’ a) Electronic mail: [email protected]; Also at Department ofComputational Science, Mathematics, and Engineering, MichiganState University, East Lansing, MI b) Also at Department of Computational Science, Mathematics, andEngineering, Michigan State University, East Lansing, MI law and charge conservation. Indeed, developing a nu-merical scheme that implicitly conserved charge was anunsolved problem until . Prior to this development,one used divergence cleaning methods to remove spuri-ous charge accumulation . The key to realizing chargeconservation relied on (a) following the de-Rham se-quence to represent physical quantities on a mesh and(b) use explicit time stepping methods. A more recentpaper prescribes three conditions must be satisfied byself-consistent charge conserving schemes ; this asser-tion is proved and illustrated for different PIC schemes.The TDFEM-PIC method relies on Maxwell solvers, inthat one solves Maxwell’s first order equations as opposedto the wave equation, and leap-frog time stepping. Thestructure of the solver is such that one avoids a timegrowing null space corresponding to DC modes. Unfor-tunately, leap frog is only conditionally stable. As a re-sult, there is a limit on the time-step sizes that one cantake, and this closely tied to the underlying discretiza-tion. In classical TDFEM, this has been overcome us-ing Newmark-beta time stepping, which is second orderand unconditionally stable. Unfortunately, implicit timestepping poses a number of challenges to satisfaction ofconservation laws that must be satisfied and is an openproblem . This paper provides the theoretical frame-work for resolving this bottleneck.Implicit time stepping permits taking significantlylarger time steps, un-constrained by the mesh; and un-conditional stability is an added bonus. Unfortunately,as will be evident in the paper, applying these directly toTDFEM-PIC violates both Gauss’ law and the equationof continuity. In addition, in solving the field equations,one needs to evolve the locations of particles over time viaNewton’s laws. A larger time step size, implies that ad-ditional infrastructure needs to be in place to accuratelycompute all aspects of particle trajectory (including in-formation necessary to map it back on the mesh). Resolu-tion to these challenges associated implicit time steppingwith a TDFEM framework will be the main contributionof this paper. We will a r X i v : . [ phy s i c s . p l a s m - ph ] F e b
1. Develop the methods to ensure that both Gauss’law and equation of continuity is satisfied for im-plicit methods. The methods rely on insight pro-vided in Ref. .2. We will show that the proposed method is agnosticto time stepping schemes.3. We will develop methods to evolve particle param-eters (path, velocity along the path, and mappingpath to the mesh).4. Finally, we will present results validating thesemethods for both the Maxwell and wave equationTDFEM solvers.Our hope is to present the technique with sufficient lucid-ity such that they can be retrofitted with existing codes.The rest of this paper is organized as follows: In thenext Section, we present an overall rubric of implicit TD-FEM solvers (both Maxwell and wave), and why directapplication of implicit time stepping fails to conservequantities. Next, in Section III, we present details onhow these may be modified so as to conserve charge, sat-isfy Gauss’ law, and be independent of time stepping ap-proach. In addition, we present details of the methodused to evolve particle parameters. In Section IV, wepresent a number of results that validate our claims. Fi-nally, we conclude this paper in Section V outlining fu-ture directions of research. II. PRELIMINARIES
Consider a domain Ω whose boundaries are denotedby ∂ Ω. It is assumed that the domain comprise chargedspecies that exist in a background medium defined by ε and µ , the permittivity and permeability of free space,and the speed of light denoted using c = 1 / √ µ ε ; forsimplicity of the exposition, we consider only one species.It is also assumed that there exists an electromagneticfield, both impressed and arising from motion of thecharged species. Both the fields and the charged speciesevolve in time. The distribution of charge can be rep-resented by a phase space distribution function (PSDF) f ( t, r , v ) that satisfies the Vlasov equation ∂ t f ( t, r , v ) + v · ∇ f ( t, r , v )+ (1) qm [ E ( t, r ) + v × B ( t, r )] · ∇ v f ( t, r , v ) = 0 . While we do not solve this equation directly, our ap-proach is conventional in that we make a particle ap-proximation for the PSDF in (1).
A. Overview of Method
Using this PSDF, we follow the conventional definitionof the charge and current density defined as ρ ( t, r ) = q (cid:82) Ω f ( t, r , v ) d v and J ( t, r ) = q (cid:82) Ω v ( t ) f ( t, r , v ) d v as mo-ments of the PSDF. The fields, E ( t, r ) and B ( t, r ), in theVaslov equation are solutions to Maxwell’s curl equationswith the sources (charge and currents) defined earlier − ∂ B ( t, r ) ∂t = ∇ × E ( t, r ) (2a) ∂ D ( t, r ) ∂t = ∇ × H ( t, r ) − J ( t, r ) (2b)and boundary conditions. These can be either Dirichletor impedance boundary conditions on ∂ Ω D or ∂ Ω I , tobound the domain,ˆ n × E ( r , t ) = Ψ D ( r , t ) on ∂ Ω D , (3a)ˆ n × B ( r , t ) µ − Y ˆ n × ˆ n × E ( r , t ) = Ψ I ( r , t ) on ∂ Ω I . (3b)Instead of using (2), the wave equation ∇ × (cid:18) µ r ∇ × E (cid:19) + 1 c (cid:15) r ∂ E ∂t = − µ ∂ J ∂t (4)can be used instead. The magnetic field can be obtainedfrom (2a) and the impedance boundary condition is de-fined using a time derivative on (3b) and using (2a). Thefields should also satisfy Gauss’ laws ∇ · D ( t, r ) = ρ ( t, r ) (5) ∇ · B ( t, r ) = 0 (6)though they are not explicitly solved.As alluded to earlier, we use the moments of PSDFto find the fields generated and then evolve their posi-tion using Newton’s equations and Lorentz force, viz., F ( t, r ) = q ( t, r )( E ( t, r ) + v ( t, r ) × B ( t, r )), and so on, forthe duration of the simulation. Thus far, our descrip-tion has been in continuous world. To perform an actualsimulation, we would need to represent all the quantitiesinvolved in terms of functions defined on a discretizationof space and time. This is typically referred to as a parti-cle in cell (PIC) approach and is the subject of our nextdiscussion.Our starting point is the representation of both Ωand ∂ Ω in terms of a finite set of tetrahedra or a meshthat contains N s nodes, N e edges and N f faces. Onthese tetrahedra, we define basis functions that followthe de-Rham sequence, enabling us to represent fields,fluxes and sources. But before proceeding too far ahead,note that we are going to follow the usual PIC cycle; (a)map charges and currents on the mesh, (b) solve for elec-tric and magnetic fields on the mesh, (c) move particlesdue to Lorentz force and find the current due to this mo-tion, and (d) find the fields due the updated sources. Thecycle then continues.The starting point of the simulation is to define thecharge and currents due to PSDF. With no loss of gen-erality, we follow the usual procedure such that ρ ( t, r ) = q α (cid:80) N p p =1 δ ( r − r p ) and J ( t, r ) = q α (cid:80) N p p =1 v p ( t ) δ ( r − r p ).This implies that PSDF is sampled with N p shape func-tions, each being a delta function. Generalization toother shape functions is possible and is agnostic to thecrux of this paper.The electric and magnetic fields are represented us-ing Whitney basis functions . Specifically, the elec-tric fields using Whitney edge basis functions, E ( t, r ) = (cid:80) N e i =1 e i ( t ) W (1) i ( r ). The magnetic flux density is rep-resented using Whitney face basis function, B ( t, r ) = (cid:80) N f i =1 b i ( t ) W (2) i ( r ). Here, N e are the number of edgesand N f are the number of faces in the mesh. Two differ-ent approaches can be used to solve Maxwell’s equations;(a) either solve them in the coupled form or (b) solve thewave equation for the electric field and then obtain themagnetic field. To set the stage for both these solvers,we introduce the following Hodge matrix operators[ (cid:63) (cid:15) ] i,j = (cid:104) W (1) i ( r ) , ε · W (1) j ( r ) (cid:105) (7)[ (cid:63) µ − ] i,j = (cid:104) W (2) i ( r ) , µ − · W (2) j ( r ) (cid:105) , (8)the surface impedance matrix[ (cid:63) I ] i,j = (cid:104) ˆ n i × W (1) i ( r ) , µ − · ˆ n j × W (1) j ( r ) (cid:105) (9)and discrete curl operator[ ∇× ] i,j = (cid:104) ˆn i , ∇ × W (1) j ( r ) (cid:105) . (10)These matrices are used to build the semidiscreteMaxwell system (cid:20) [ I ] 00 [ (cid:63) (cid:15) ] (cid:21)(cid:124) (cid:123)(cid:122) (cid:125) ¯¯ C M (cid:20) ∂ t ¯ B∂ t ¯ E (cid:21) + (cid:20) ∇× ] c [ ∇× ] T [ (cid:63) µ − ] c [ (cid:63) I ] (cid:21)(cid:124) (cid:123)(cid:122) (cid:125) ¯¯ K M (cid:20) ¯ B ¯ E (cid:21) = (cid:20) − ¯ J(cid:15) (cid:21)(cid:124) (cid:123)(cid:122) (cid:125) ¯¯ F M (11)where the degree of freedom vectors ¯ E =[ e ( t ) , e ( t ) , . . . , e N e ( t )], ¯ B = [ b ( t ) , b ( t ) , . . . , b N f ( t )],and ¯ J = [ j ( t ) , j ( t ) , ...j N e ( t )] with j i ( t ) = (cid:104) W (1) i ( r ) , J ( t, r ) (cid:105) . For the wave equation, the sys-tem becomes[ (cid:63) (cid:15) ] (cid:124)(cid:123)(cid:122)(cid:125) ¯¯ M W ∂ t ¯ E + c [ (cid:63) I ] (cid:124)(cid:123)(cid:122)(cid:125) ¯¯ C W ∂ t ¯ E + c [ (cid:63) S ] (cid:124) (cid:123)(cid:122) (cid:125) ¯¯ K W ¯ E = − ∂ t ¯ J (12)where [ (cid:63) S ] = [ ∇× ] T [ (cid:63) µ − ][ ∇× ]. B. Unconditionally Stable Time Marching
The mixed finite element system in (2) is discretizedin time using Newmark-Beta, an unconditionally sta-ble time stepping method. This method has been ex-tensively used in for the wave equation and examined for the mixed finite element method in , allowing formuch larger time step sizes than the traditional leapfrogmethod. In this method, the fields in time are repre-sented by three temporal basis functions N n +1 − i ( t ) = (cid:88) j =0 j (cid:54) = i t − t n +1 − j t n +1 − i − t n +1 − j (13)corresponding to i ∈ [1 ,
3] and weighting function W ( t ) = t n − t ∆ t t ∈ [ t n − , t n ] t − t n ∆ t t ∈ [ t n , t n +1 ]0 otherwise . (14)This combination of basis function and weighting func-tion creates a non-disappative, unconditionally stabletime marching scheme, which can be written as recur-rence formula provided in , corresponding to parame-ters γ = 0 . β = 0 .
25. When applied to (2), thisbecomes(0 . C M + 0 . t ¯¯ K M ) ¯ X n +1 − . t ¯¯ K M ¯ X n +(0 . C M + . t ¯¯ K M ) ¯ X n − + 0 . t ¯ F n +1 M +0 . t ¯ F nM + 0 . t ¯ F n − M = 0 (15)where ¯ X i = [ ¯ B i,T , ¯ E i,T ] and ¯ F iM = [0 , − (cid:15) − ¯ J i,T ]. Like-wise, (12) becomes( ¯¯ M W + 0 . t ¯¯ C W + 0 . t ¯¯ K W ) ¯ E n +1 +( − M W − . t ¯¯ K W ) ¯ E n ( ¯¯ M W + 0 . t ¯¯ C W + . t ¯¯ K W ) ¯ E n − +0 . t ¯ F n +1 W + 0 . t ¯ F nW + 0 . t ¯ F n − W = 0 (16)However, treating the current as written in (2) will notpreserve the necessary link between Ampere’s law andGauss’ law needed to create a charge conserving scheme.This is immediately apparent after applying a discretedivergence operator to semidiscrete Ampere’s law. Afterusing the identity that [ ∇· ][ ∇× ] T = 0, this operationyields0 . ∇· ][ (cid:63) (cid:15) ] ¯ E n +1 − . ∇· ][ (cid:63) (cid:15) ] ¯ E n − = − . t [ ∇· ] ¯ J n +1 − . t [ ∇· ] ¯ J n − . t [ ∇· ] ¯ J n − . (17)When the same operator is applied to the semidiscretewave equation, one gets[ ∇· ][ (cid:63) (cid:15) ] ¯ E n +1 − ∇· ][ (cid:63) (cid:15) ] ¯ E n + [ ∇· ][ (cid:63) (cid:15) ] ¯ E n − − . t [ ∇· ] ¯ J n +1 − . t [ ∇· ] ¯ J n − . t [ ∇· ] ¯ J n − . (18)Making the substitution of ¯ ρ i = [ ∇· ][ (cid:63) (cid:15) ] ¯ E i , then it be-comes clear that neither (17) nor (18) satisfy Gauss’ lawor the continuity equation. Instead, a different treat-ment of the right hand side, the particle current density,is needed in order to create a charge conserving scheme. III. MODIFIED TDFEM-PICA. Integrator Agnostic Charge Conserving Scheme
It is apparent that, as written, time conservation failsfor both Maxwell solver and the wave equation. The rea-sons are two fold: (a) the order of time derivatives on thecurrent (on the right hand side) and those on the electricfield are off by one; (b) this requires the discrete time in-tegrator to remember initial conditions. The latter holdsthe key to solving the puzzle. Newmark time steppingschemes are, in effect, stable time integrators. The cruxof our approach is to correctly evaluate the time integralof the current. As elucidated in , the time integral ofthe current is readily obtainable, and indeed a part of thePIC scheme. Specifically, starting with the definition ofthe PDSF, − ¯ ρ n ( r ( t n )) = [ ∇· ] (cid:90) t n ¯ J n ( r ( τ )) dτ. (19)where t n = n ∆ t . As shown in this equation can berewritten as,¯ ρ n ( r ( t n )) = − [ ∇· ] (cid:90) r ( t n ) r (0) ¯ J i (˜ r ) d ˜ r . (20)Following the details in Ref. , it is immediately apparentthat for any particle p (cid:90) t n dτ v p ( τ ) δ ( r − r p ( τ )) = (cid:90) r p ( t n ) r p (0) d ˜ r δ ( r − ˜ r ) (21)Note, for each charge, its trajectory is determined bythe solution to Newton’s equations. The integrationalong a particle path can be computed to a high degree ofaccuracy. To develop a charge conserving methodology,we define ¯ G n = (cid:90) t n ¯ J n ( τ ) dτ (22)This equation is readily evaluated using (20). It followsthat instead of using ¯ J in (2) (and therefore, in (17))one can instead use ∂ t ¯ G . Discrete implementation witha Maxwell equation solver results in in the divergence ofAmpere’s law to be0 . ∇· ][ (cid:63) (cid:15) ] ¯ E n +1 − . ∇· ][ (cid:63) (cid:15) ] ¯ E n − = − . ∇· ] ¯ G n +1 + 0 . ∇· ] ¯ G n − . (23)Examining (24) term by term reveals that both sides ofthe equation are identical given that [ ∇· ] ¯ G n = − ¯ ρ n . Ina similar manner, one can use ∂ t ¯ G instead of ∂ t ¯ J in (12)to yield¯ ρ n +1 − ρ n + ¯ ρ n − = − [ ∇· ] ¯ G n +1 − ∇· ] ¯ G n − +[ ∇· ] ¯ G n − . (24) Here, we have taken the liberty of substituting, ¯ ρ n =[ ∇· ][ (cid:63) (cid:15) ] ¯ E n . At this point, we note that the proposedapproach is agnostic to the time stepping scheme (or in-tegrator) used to solve Maxwell’s equations; both theequation of continuity and Gauss’ laws are satisfied bydesign.A word of caution is in order before we proceed. While,the method developed is exact, it should be noted thatto obtain ¯ E n +1 , one needs to solve either (17) (or (18))with the appropriate substitutions for ¯ G n +1 instead of¯ J n +1 . Obviously, the solution to these sets of equationsis subject to errors that arise due to vagaries of a linearalgebraic solution (tolerances, excitation of null-spaces,etc). As a result, as will be seen in the results section, ourerrors are small but not identically zero. Next, we discussa higher order particle pusher to solve the equations ofmotion consistently. B. Particle Pusher
Using an implicit time stepping scheme has advantagesas well as challenges. The principal advantage is an un-conditionally stable time step size independent schemeas opposed to a conditionally stable scheme like leap frogwhose stability depends on the time step size. The down-side of using large time steps is that one must capturethe nuances of both the path and velocity of the parti-cle. Thus, solving the equations of motion using a Borispush , with its linear representation of the particle po-sition and velocity can introduce large errors into thesystem. Our goal is to develop a higher order scheme.As is well known, the particle positions and veloci-ties are updated by solving Newton’s equations via theLorentz force, giving us the following set of coupled firstorder ODEs for each particle, ∂ t v p ( t, r p ) = a p ( t, r p ) = q α m α ( E ( t, r p ) + v p ( t, r p ) × B ( t, r p ))(25) ∂ t r p ( t, r p ) = v p ( t, r p ) . (26)These form a pair of first order ODEs in time, and thereare a number of methods that can be applied. Our choiceis to use a higher order order Adams-Bashforth scheme.An exemplar recursion relation for v and r for a 4th orderAdam’s-Bashforth method is as follows: v n +1 p = v np + ∆ t
24 (55 a np − a n − p + 37 a n − p − a n − p )(27) r n +1 p = r np + ∆ t
24 (55 v np − v n − p + 37 v n − p − v n − p ) . (28)where ∆ t is the time step size. Given that the Newmarkscheme is second order, we choose the Adams-Bashforthscheme to be at least two orders higher so as to accommo-date a second time derivative in on ¯ G n . The path usedfor interpolating the position is a fourth order Lagrangepolynomial k = 4 + 1 is defined as, r p ( t ) = k (cid:88) j =0 r n − jp (cid:96) j ( t ) (29a) (cid:96) ( t ) = (cid:89) ≤ m ≤ km (cid:54) = j t − t n +1 − m t n +1 − j − t n +1 − m , (29b)where r p ( t ) is the position at time t and r np is the locationof particle p at the t = n ∆ t . C. Particle Path and Current Mapping
The final piece of the puzzle is mapping the path tothe underlying tesselation. In order to do so, we notethat that the integrator used to solve the equation ofmotion implicitly assumes a Lagrange polynomial inter-polant. As a result, the order of the method used mapsto order of the interpolant. This information needs to beused to find out where the particle enters and leaves thecell.Once the particle locations at each time step are knownfrom the particle push, the path through the unstruc-tured mesh needs to be found. This includes finding thelocations of where a particle enters a cell and where itleaves and is detailed in Algorithm 1. Since we are usinga higher order representation of a particle path, findingthese entry and exit points of the cell with the tetrahe-dron becomes a non-linear problem and is detail in Algo-rithm 2. Assume that we are given the normal to surfaceˆ n , and vertices of the triangle r v, , r v, , r v, . The inter-section between the trajectory r p ( t ) and the plane canbe obtained by solvingˆ n · [( r p ( t ) − r v, ) × ( r v, − r v, )] = 0 (30) Algorithm 1
Particle Path Finding Algorithm Push particle finding r p,f if r p,f is in same cell as r p,s then if All quadrature points between are in same cell then Integrate using a quadrature rule along path. Return and go onto next particle end if end if Find exit point r p,i of path in cell Integrate from r p,s to r p,i while Path not complete do Find next cell that path travels through if r f is in same cell as r p,i then if All quadrature points are in same cell then
Integrate using a quadrature rule along path.
Path is complete
Return and go onto next particle. else
Find exit point of path in cell.
Integrate from r p,is to r p,if . end if end if end while Algorithm 2
Non-Linear Bi-Section Method t s = 0, t f = 1, t h = 0 . if Any quadrature points are outside of the cell. then t f = t q end if while | t s − t f | < tol do if t h is in same cell as t s then t s = t h else t f = t h end if t h = 0 . t s + t f ) end while Note, the path r p ( t ) can be parameterized using (29).Using this parameterization, one can use a non-linear it-eration (such as Newton-Raphson) to solve (30). For con-vinience, we take a simpler approach by implementinga bi-section method that moves along the path check-ing whether candidate points are inside or outside ofthe cell. For test cases presented in this paper, thismethod converges rather robustly. Every step takesaround 47 steps to converge below a tolerance of 1 · − (0 . = 7 . · − ). Once the method converges, we thencompute the integral along each path segment in eachcell using a set of quadrature points. To illustrate this,consider Fig. 1 containing an example particle startingposition r p,s , and finishing position r p,f with the inter-section point being r p,i . The quadrature points would liealong the path between the r p,s and r p,i , then anotherset of quadrature points between r p,i and r p,f .Before we discuss results obtained using the above ap-proach, a few points are in order; to evaluate ¯ G n , (a)the integral over the path be evaluated using quadraturerules to very high precision as the order of the path isFIG. 1: Particle path for a single particles with startand location and intersection point.know; (b) when the path passes though multiple cells,the integration is broken up into pieces over each cell; (c)one can save on computational cost of by updating theintegral. IV. RESULTS
In this Section, we present a number of results demon-strating the efficacy of the proposed scheme with respectto conservation laws, as well accuracy of key steps thatare integral to the process.
A. Higher Order Particle Motion
One of the key advantages in using implicit time step-ping is the possibility of using much larger time step sizes.Unfortunately, this also implies that one needs higherorder methods to capture both the path as well as ve-locity. In this section, we demonstrate convergence ofour algorithm for particle motion using various orders ofAdams-Bashforth integrator and compare these to stan-dard non-relativistic Boris push.To do so, we set up a classic cyclotron motion testwhere a single particle was given an initial velocity ina constant magnetic field resulting in circular motion asshown in Fig. 3. The parameters are shown in TableI with a particle’s initial velocity v with a backgroundmagnetic fields B with a given mass m and charge q . Theparticle will move in a circle due to the Lorentz force asshown in Fig. 2. The relative error in both position andvelocity for various time step sizes with multiple orderof Adams-Bashforth and Boris are is shown in Fig. 3.The average error is calculated by taking the norm ofthe distance errors of each point r and dividing by thenormal of the analytic positions r a (see Ref. for details), error = || r − r a || || r a || (31)The slopes for each of the Adams-Bashforth methodsmatch its order. Boris on the other hand has a second or-der velocity update with a first order positional update.This test essentially validates out pusher as well as helps FIG. 2: Mean relative error in position forAdam-Bashforth Orders 1-5 compared with the Borispush, shown in black.FIG. 3: Mean relative error in velocity forAdam-Bashforth Orders 1-5 compared with the Borispush, shown in black.correlate error (or approximately so) in particle motionwith time step size.TABLE I: Cyclotron Motion Parameter Value B . · − ˆ z TQ − . · − Cm 9 . · − kg v · ˆ y m/s r [0 . , . , .
0] m
B. Expanding Particle Beam
Next, we consider an expanding beam test . An ex-panding particle beam is injected into a cylindrical cavitywith an initial velocity of magnitude v . As the beamtravels down the tube, the electrons repel each othercausing the beam beam to expand. This expansion ratecan be compared with other codes to validate the solu-tion. The detail of the mesh and beam parameters usedare shown in Table II.TABLE II: Expanding Particle Beam Parameters Parameter Value
Cavity Radius 20 mmCavity Length 100 mmBoundary Conditions PEC v · m/s v /c r b t ns Both the wave equation and mixed finite element tra-jectories are compared in Fig. 4 and show good agree-ment with XOOPIC (an extensively used and well val-idated quasi-2D FDTD code). We sample the electricfield half way down the tube 16 mm from the center ofthe tube. The radial field values are plotted over timeshown in Fig. (5) for simulations with different timesteps. We compare four runs with time steps of α ∆ t where α is scale factor and ∆ t = 0 .
333 ps is the largeststable step size in a leap frog time marching method forthe given mesh. Note, 2 ns corresponds to 1 transit ofthe tube. It is evident from this figure that the proposedmethod provides stable results; indeed, as is evident fromthis figure, the data at 7.5∆ t , 15∆ t and 30∆ t are almostidentical to each other, where as the one at 1485∆ t isslightly different. This points to significant gains thatcan be made with Newmark time stepping (provided themethod is charge conserving).This leads to the next argument. Shown in Fig. 6 isdata from two different methods for the same set up runusing MFEM with backward difference at ∆ t , MFEMwith Newmark at 7∆ t and the wave equation (WE) at7∆ t . As evident, all three methods conserve charge toalmost machine precision. It should be noted that bothMFEM and WE have a null space. In the case of theformer, it is fields that behave like ∇ φ ( r ), and the latter,as t ∇ φ ( r ). However, as is evident from these results, ourmapping on to these null spaces is small and behaves asexpected.To further illustrate the robustness of the method totime step sizes, in Fig. 7 we compare the satisfaction ofGauss’ law for all four time steps used in Fig. 5. Asis evident from here, charge is again conserved almostto machine precision (around 10 − for all with slight FIG. 4: Expanding particle beam macro particles in thez vs r plan. Particle locations from both mixed finiteelement methods and wave equation versions arecompared with XOOPIC beam profile.FIG. 5: Electric field values are the radial componenthalf way down the tube 16 mm from the center of thetube. Multiple simulation with different time steps areperformed.difference evolution of trajectory). C. Adiabatic Expanding Plasma
Finally, for a third validation case we simulate an adi-abatic expansion of a plasma ball with radial Gaussiandistribution in the radial direction. This case has an ana-lytic solutions and allows for good comparison and vali-dation. We change some of the parameters from the orig-inal numerical experiments such that the Debye lengthcan be fully resolved. This example is described in moredetail in . We simulate this example both MFEM andWE. For both examples we get excellent agreement inthe expansion rate with both the wave equation, Fig.FIG. 6: Discrete Gauss’s Law error per particle forNewmark-Beta mixed finite element (NM-MFEM),Newmark-Beta wave equation (WE), backwardsdifference mixed finite elements (BD-MFEM) using thecharge conservation technique provided here.FIG. 7: Discrete Gauss’s law error per particle varioustime steps using the mixed finite element methods usingNewmark time stepping.(9), and the mixed formulation, Fig. 8, when comparedwith analytic densities. V. SUMMARY
In this paper, we have presented a solution to a prob-lem that has been long-standing–charge conserving FEM-PIC methods for implicit time stepping systems withoutthe need to adopt divergence cleaning. In other words,rubrics have been developed such that conservation lawsare implicitly obeyed. Indeed, the method presented isagnostic to any time stepping scheme. We have demon-strated the efficacy of this approach for a set of testproblems, using different time step sizes and different TABLE III: Adiabatic Expanding Plasmas
Parameter Value
Mesh Radius 6mmBoundary Conditions First order ABC T ion T electron Sr + Macro-Particle Size 52012.58Min Edge Length 1.529mmMax Edge Length 6.872mm
FIG. 8: Mixed Finite Element particle beam expansiontime stepping schemes, as well as both MFEM and WEsolvers. The results reliably attest our claims. The aboveapproach opens multiple doors that will further the stateof art of FEM-PIC; these include higher order schemes inboth space and time, quasi-Helmholtz decomposition toget a better handle on null-spaces, and domain decompo-sition to effect rapid solution by parallelizing the scheme.Papers on these will be presented soon in other forums.FIG. 9: Wave Equation Adiabatic Expanding Plasma
ACKNOWLEDGMENTS
This work was supported by SMART Scholarship pro-gram. We thank the MSU Foundation for supportthrough the Strategic Partnership Grant during earlyportion of this work. This work was also supported by theDepartment of Energy Computational Science GraduateFellowship under grant DE-FG02-97ER25308. The au-thors would also like to thank the HPCC Facility, Michi-gan State University, East Lansing, MI, USA.
DATA AVAILABILITY
The data that support the findings of this study areavailable from the corresponding author upon reasonablerequest. R. Marchand, “Ptetra, a tool to simulate low orbit satellite–plasma interaction,” IEEE Transactions on Plasma Science ,217–229 (2011). R. Lemke, T. Genoni, and T. Spencer, “Three-dimensionalparticle-in-cell simulation study of a relativistic magnetron,”Physics of Plasmas , 603–613 (1999). E. Fourkal, B. Shahine, M. Ding, J. Li, T. Tajima, and C.-M.Ma, “Particle in cell simulation of laser-accelerated proton beamsfor radiation therapy,” Medical Physics , 2788–2798 (2002). C. K. Birdsall and A. B. Langdon,Plasma physics via computer simulation (CRC press, 2004). J. P. Verboncoeur, “Particle simulation of plasmas: review andadvances,” Plasma Physics and Controlled Fusion , A231(2005). C. Nieter, J. R. Cary, G. R. Werner, D. N. Smithe, and P. H.Stoltz, “Application of dey–mittra conformal boundary algo-rithm to 3d electromagnetic modeling,” Journal of Computa-tional Physics , 7902–7916 (2009). J. Squire, H. Qin, and W. M. Tang, “Geometric integrationof the vlasov-maxwell system with a variational particle-in-cellscheme,” Physics of Plasmas , 084501 (2012). P. Monk, Finite element methods for Maxwell’s equations (Ox-ford University Press, 2003). A. S. Glasser and H. Qin, “The geometric theory of chargeconservation in particle-in-cell simulations,” arXiv preprintarXiv:1910.12395 (2019). C. S. Meierbachtol, A. D. Greenwood, J. P. Verboncoeur, andB. Shanker, “Conformal electromagnetic particle in cell: A re-view,” IEEE Transactions on Plasma Science , 3778–3793(2015). J.-M. Jin, The finite element method in electromagnetics (JohnWiley & Sons, 2015). M. C. Pinto, S. Jund, S. Salmon, and E. Sonnendr¨ucker,“Charge-conserving fem–pic schemes on general grids,” ComptesRendus Mecanique , 570–582 (2014). H. Moon, F. L. Teixeira, and Y. A. Omelchenko, “Exact charge-conserving scatter–gather algorithm for particle-in-cell simula-tions on unstructured grids: A geometric perspective,” ComputerPhysics Communications , 43–53 (2015). C.-D. Munz, P. Omnes, R. Schneider, E. Sonnendr¨ucker, andU. Voss, “Divergence correction techniques for maxwell solversbased on a hyperbolic model,” Journal of Computational Physics , 484–511 (2000). Z. D. Crawford, S. O’Connor, J. Luginsland, and B. Shanker,“Rubrics for charge conserving current mapping in finite ele-ment particle in cell methods,” arXiv preprint arXiv:2101.12128(2021). G. Chen, L. Chac´on, and D. C. Barnes, “An energy-and charge-conserving, implicit, electrostatic particle-in-cell algorithm,”Journal of Computational Physics , 7018–7036 (2011). Z. Crawford, J. Li, A. Christlieb, and B. Shanker, “Uncondi-tionally stable time stepping method for mixed finite elementmaxwell solvers,” Progress In Electromagnetics Research ,17–30 (2020). O. C. Zienkiewicz, “A new look at the newmark, houbolt andother time stepping formulas. a weighted residual approach,”Earthquake Engineering & Structural Dynamics , 413–418(1977). J. P. Boris, “Relativistic plasma simulation-optimization of a hy-brid code,” in Proc. Fourth Conf. Num. Sim. Plasmas (1970) pp.3–67. S. O’Connor, Z. Crawford, J. Verboncoeur, J. Lugisland, andB. Shanker, “A set of benchmark tests for validation of 3d particlein cell methods,” arXiv preprint arXiv:2101.09299 (2021). V. Kovalev and V. Y. Bychenkov, “Analytic solutions to thevlasov equations for expanding plasmas,” Physical review letters90