Rubrics for Charge Conserving Current Mapping in Finite Element Particle in Cell Methods
Zane D. Crawford, Scott O'Connor, John Luginsland, B. Shanker
RRubrics for Charge Conserving Current Mapping in FEM-PIC
Rubrics for Charge Conserving Current Mapping in Finite ElementParticle in Cell Methods
Zane D. Crawford, a) Scott O’Connor, a) John Luginsland, and B. Shanker Department of Electrical and Computer Engineering, Michigan State University, East Lansing,MI (Dated: 29 January 2021)
Modeling of kinetic plasmas using electromagnetic particle in cell methods (EM-PIC) is a problem that is well worn, inthat methods developed have been used extensively both understanding physics and exploiting them for device design.EM-PIC tools have largely relied on finite difference methods coupled with particle representations of the distributionfunction. Refinements to ensure consistency and charge conservation have largely been an ad-hoc efforts specific tofinite difference methods. Meanwhile, solution methods for field solver have grown by leaps and bounds with significantperformance metrics compared to finite difference methods. Developing new EM-PIC computational schemes thatleverage modern field solver technology means re-examining analysis framework necessary for self-consistent EM-PIC solution. In this paper, we prescribe general rubrics for charge conservation, demonstrate how these are satisfiedin conventional finite difference PIC as well as finite element PIC, and prescribe a novel charge conserving finiteelement PIC. Our effort leverages proper mappings on to de-Rham sequences and lays a groundwork for understandingconditions that must be satisfied for consistency. Several numerical results demonstrate the applicability of these rubrics.
I. INTRODUCTION
Modeling plasmas and other charged media as sources is aneeded capability to better understand and design many real-world devices. Of the many techniques to model this phe-nomena, the particle-in-cell (PIC) method has been a popu-lar approach. A wide variety of fields depend on PIC as acomputational method, ranging from fundamental science ofspace weather, nuclear fusion, atomic processes on one hand,to high technology applications like particle accelerators, co-herent electromagnetic radiation sources, and plasma process-ing devices on the other. PIC, at a glance works as follows;charged media, represented as a collection of particles, movein a domain under the influence of fields that are due to mo-tion of these particles. This calls for self consistent solutionof Maxwell’s equations in conjunction with equation of mo-tion. Given application drivers stated earlier, there has beenextensive development of EM-PIC methods . Several ofthese have been used extensively in the design of devices .While traditional application were electrically large albeit ge-ometrically simple, there is burgeoning interest in deviceswith significantly geometric complexity but smaller footprint.As a result, there is need for EM-PIC methods with greatercapabilities .But before we delve extensively into nuances and needsof modern EM-PIC codes, it helps to understand a bit ofhistory of challenges that have been overcome in develop-ing EM-PIC methods. The principal challenge in EM-PICmodeling is self-consistency of solving Maxwell’s equationtogether with equations of motion so as to actually capturedynamics of the distribution function. Challenges lies in in-corporation of sources (and their variation over time) consis-tently into Maxwell’s equation. As a result, there has been a) Also at Department of Computational Science, Mathematics, and Engineer-ing, Michigan State University, East Lansing, MI extensive work in understanding how one may develop meth-ods that conserve charge . Early development of EM-PIC methods largely relied on finite difference time domain(FDTD-PIC), as did addressing many of the aforementionedchallenges . More recent PIC development has seen ex-tensive development mixed finite element methods (MFEM-PIC) and discontinuous-Galerkin (DG-PIC) as the EM fieldsolver .These developments are specially pertinent as EM-PICmodels that can more accurately capture geometry andphysics with greater computational efficiency would soon benecessary to modeling physics ranging from novel accelera-tors to high frequency semiconductor devices. As the com-plexity of system to be modeled increases one needs EM-PICsuites that leverage advances in both finite element and par-ticle methods. But as is to be expected, integrating the twocalls for ensuring that self-consistency is preserved. This im-plies that discrete solutions of Maxwell’s equations togetherwith Newton’s laws satisfy all of Maxwell’s equations to-gether with the equation of continuity at every time step. Thisis easier said than done. Errors can manifest themselves asspurious charges and currents. To overcome these, it is pos-sible to use procedures such as divergence cleaning to correctthe erroneous fields that arise. Approaches to correct the fieldsinclude using Lagrange multipliers , solving Poisson’s equa-tion at multiple time steps , or changing the field solver .However, there can be a significant computational cost to us-ing these methods. Of course, a more desirable alternative isto prescribe a rubric in which charge is conserved without cor-rection terms or additional equations, regardless if sources areincluded or not.Developing computational methods that are inherentlycharge conserving has been a topic of interest as long ago asthere have been EM-PIC codes . The benefits of having acharge conserving scheme are apparent; they avoid additionalwork necessary to effect divergence cleaning. Unfortunately,their development has always been closely tied to a particu-lar EM-solution technique, be in FDTD-PIC or MFEM- a r X i v : . [ phy s i c s . c o m p - ph ] J a n ubrics for Charge Conserving Current Mapping in FEM-PIC 2PIC . The latter is novel in that it develops a method thatrelies on de-Rham maps to define the proper basis sets forrepresenting the electric field and magnetic flux density (EB-MFEM-PIC). In effect, the primal quantities are defined basedon Faraday’s laws. The current mapping schemes demonstratehow they recover the results found in the early literature .As methods developed have been solver specific, significantattention has not been paid to examining the underlying math-ematical structure and why that is necessary for charge con-servation. In this work, we seek to do just that.While unusual, we start our analysis at the particle approxi-mation of the distribution function and develop a simple albeitrigorous approach to demonstrate charge continuity relation-ships between charge density and currents. The exposition inthis paper straightforward and uses the sifting property of thedelta function. This is unlike that in presented in that re-lies on an discrete approximation. As an added benefit, thisapproach lays the groundwork for developing conditions formapping particles and fields on any discrete grid, and con-sistently evolving their position and strength. In effect, wecreate a uniform framework and criterion that is necessaryto solve all of Maxwell’s equation consistently, thereby sat-isfying charge conservation as well. The rubrics prescribedare in agnostic to the method used to solve Maxwell’s equa-tions. Indeed, we demonstrate how these are applicable toFDTD and (M)FEM. We also use a these rubrics to develop anovel MFEM-PIC scheme that uses electric flux density andthe magnetic field (DH-MFEM-PIC) and validate its proper-ties including charge conservation, satisfaction of Gauss’ lawas well as comparison against quasi-analytic data.The rest of the paper is organized as follows: After definingthe EM-PIC problem in Section II and its discetization usingparticle approaches, we define the three rules in Section III inthe context of preserving the spatial and temporal connectionsbetween Gauss’ law and Ampere’s law. In Section IV, we pro-vide examples on how FDTD and two formulations of a vari-ant of FEM, the mixed finite element method (MFEM), satisfythe three conditions. One of these is novel in that it relies onthe electric flux density and the magnetic field (DH-MFEM).Next, in Section V, we provide results that show what hap-pens when the conditions are violated as well as affirm thecorrect behavior of MFEM when the conditions are met. Inaddition, several results are provided for the DH-MFEM-PICthat demonstrate necessary features for a number of examples(charge conservation, satisfaction of Gauss’ laws, correctnessagainst quasi-nalaytical data). In an Appendix, we present de-tails of the MFEM discretizations and other proofs. II. PROBLEM STATEMENT AND FORMULATION
Consider a domain Ω whose boundaries are denoted by ∂ Ω .It is assumed that domain comprises of charged species thatexist in a background medium defined by ε and µ , the per-mittivity and permeability of free space; for simplicity of theexposition, we consider only one species. It is also assumedthat there exists an electromagnetic field, both impressed andarising from motion of the charged species. Both the fields and the charged species evolve in time. The evolution of thecharged species is modeled via a modeled as a conservativephase-space distribution function (PSDF) and those of fieldsusing Maxwell’s equations. To this end, we assume that thePSDF f ( t , r , v ) must satisfy the Vlasov-Maxwell equation ∂ t f ( t , r , v ) + v · ∇ f ( t , r , v ) + qm [ E + v × B ] · ∇ v f ( t , r , v ) = . (1)The electric field E ( t , r ) and magnetic flux density B ( t , r ) sat-isfy Maxwell’s equations − ∂ t B ( t , r ) = ∇ × E ( t , r ) (2a) ∂ t D ( t , r ) = ∇ × H ( t , r ) − J ( t , r ) (2b) ∇ · B ( t , r ) = ∇ · D ( t , r ) = ρ ( t , r ) , (2d)where fields and fluxes are related through constitutive param-eters B ( t , r ) = µ H ( t , r ) (3a) D ( t , r ) = ε E ( t , r ) . (3b)The moments of the PSDF yield the the charge density ρ ( t , r ) and current density J ( t , r ) , defined as ρ ( t , r ) = q (cid:90) Ω f ( t , r , v ) d v (4a) J ( t , r ) = q (cid:90) Ω v ( t ) f ( t , r , v ) d v . (4b)It follows from the above equations that charge and currentdensities are related via the equation of continuity − ∂ t ρ ( t , r ) = ∇ · J ( t , r ) . (5)Further, assume that these equations satisfy the initial condi-tions arising from ∇ · D ( , r ) = ρ ( , r ) (6a)and ∇ · B ( , r ) = ∂ Ω may be partitioned intoseveral subdomains, Γ i , such that Γ : = ∪ i Γ i . Each Γ i has aboundary condition, either Dirichlet ( Γ D ), Neumann ( Γ N ), orimpedance ( Γ I ). The boundary conditions are defined asˆ n × E ( t , r ) = Ψ D ( t , r ) on Γ D (7a)ubrics for Charge Conserving Current Mapping in FEM-PIC 3 PushParticles MapSourcesto MeshSolve forFieldsInterpolateFields
FIG. 1: Flow of typical EM-PIC simulation.ˆ n × B ( t , r ) µ = Ψ N ( t , r ) on Γ N (7b)ˆ n × B ( t , r ) µ − Y ˆ n × ˆ n × E ( t , r ) = Ψ I ( t , r ) on Γ I (7c)where ˆ n the vector normal to the surface ∂ Ω , Y is the surface admittance defined as (cid:112) ε / µ , and Ψ D ( r , t ) , Ψ N ( r , t ) , and Ψ I ( r , t ) represent some function.These equations describe a self consistent set that mustbe solved numerically. In what follows, we describe thenumerics necessary to do so. III. SELF CONSISTENT MODELING FRAMEWORK
While we do not solve the Vlasov system directly, we fol-low the usual practice of representing the PSDF using a col-lection of basis sets and then evolving the trajectory of these inusing Newton’s law in concert with Maxwell equations. Thebasis sets are defined on a cell or simplical complex whichdiscretizes a volume Ω bounded by ∂ Ω . The mesh consistsof N g nodes, N e edges, N f faces, and N v volumes (in threedimensions).This process of modeling the charged media can be de-scribed as the EM-PIC cycle as seen Fig. 1. It begins withmapping the charged particles on the grid, forming a dis-cretized representation of ρ ( t , r ) and J ( t , r ) . These discretizedcurrents and charges act as sources in discretized Maxwellsolvers, and are used to evolve the fields for one time step. The“new” fields and flux densities are then used update the loca-tions of the charges and thereby obtain currents (and charges).The cycle then repeats. What is apparent from the the aboveare the the nuances of what ensures self consistent modelingof Maxwell and Vlasov systems discrete grids; this is the fo-cus of this Section. A. Sampling the distribution function, and equation ofcontinuity
To begin, we start with representing the PSDF using f ( t , r , v ) = N P ∑ p = S ( r − r p ( t )) δ ( v − v p ( t )) (8) where shape functions, S ( r ) , represent the spatial support ofparticles. It follows from the definitions in (4) that ρ ( t , r ) and J ( t , r ) can be represented as ρ ( t , r ) = q N p ∑ p = S ( r − r p ( t )) (9a) J ( t , r ) = q N p ∑ p = v p ( t ) S ( r − r p ( t )) . (9b)Given the discrete representation of the PSDF, it can beshown that the equation of continuity is still holds in thisdiscrete/particle setting. While this has been proven byEastwood , the approach we take is straightforward, and ex-ploits the sifting property of the delta function. Specifically,we use the definition r p ( t ) = r p ( t )+ (cid:82) t d τ v p ( τ ) and the iden-tity δ ( r − r p ( t )) = δ ( r − r p ( t )) (cid:63) s δ (cid:18) r − (cid:90) t d τ v p ( t ) (cid:19) (10)Using this identity, it can be shown that (see Appendix) ∂ t δ ( r − r p ( t )) = − ∇ · [ v p ( t ) δ ( r − r p ( t ))] (11)The equation of continuity follows via − ∂ t ρ ( t , r ) = − q N p ∑ p = ∂ t S ( r − r p ( t )) = − q N p ∑ p = S ( r ) (cid:63) ∂ t δ ( r − r p ( t ))= ∇ · (cid:34) q N p ∑ p = v p ( t ) S ( r − r p ( t )) (cid:35) = ∇ · J ( t , r ) . (12)Of course, this is not the only identity that can be derived.As shown in the Appendix, another relationship pertinent tothis paper that follows from the spectral representation of (10)is ρ ( t , r ) = ρ ( , r ) − q ∇ · N p ∑ p = (cid:90) r p ( t ) r p ( ) d ˜ r S ( r − ˜ r ( t )) (13)As shown in the Appendix, this can be transformed into a lineintegral along the path taken by each particle. Our analysis,thus far, has dwelt on the representation of the PSDF and itsmoments, viz., the charge and current. As is evident from (12)that the equation of continuity is satisfied when the PSDF isrepresented by arbitrarily shaped particles. Next, we addressmapping/measuring these charges and currents on grid (andhow they interact with Maxwell’s equations). B. Self-consistent solution of Maxwell’s Equations withactive sources
For self consistency, one needs to solve Maxwell’s equa-tions with sources defined by (9). In this pursuit, we dis-cretize Maxwell’s equations in space with a collection of cellsubrics for Charge Conserving Current Mapping in FEM-PIC 4FIG. 2: Diagram of electromagnetic quantities as differentialforms.on which we define basis functions to represent fields andsources. Additionally, a temporal discretization is chosen suchthat a consistent solution for Newton’s laws is obtained to up-date the position and trajectory of the particle. It is well knownthat de-Rham diagrams provide appropriate map betweenfields, fluxes and sources in Maxwell’s equations. In addi-tion, one needs additional framework to ensure that chargesand currents are correctly measured and integrated in time. Inwhat follows, we discuss each in turn.
1. de-Rham complex
To begin, de-Rham diagrams establish relations betweenfield, fluxes and sources. These diagrams replicate Maxwell’sequations in the discrete world; there exists extensive litera-ture on their use in developing finite elements . In a nut-shell, these permit the definition of Whitney spaces that canbe used to represent primary quantities ( φ , E , B ) or their duals( ρ , H , D ). It should be noted that the choice of fields on theprimal or dual grids are not rigid; in other words, ( ρ , H , D / J )can be defined on the primal grid and ( φ , E , B ) on the dualgrid. These relationships are shown in Figs. 2 and 3. Us-ing Whitney spaces to solve Maxwell’s equations has its an-tecedents in Bossavit’s seminal work and those of severalothers ; our presentation is along these lines. As has beenshown in these papers, these basis sets and their functionspaces are defined as follows: the 0-form defined on nodes,the 1-form defined on edges, the 2-form defined on faces, andthe 3-form defined on volumes. In this work, we choose touse the Whitney basis sets that0-form : = W ( ) ∈ H (14a)1-form : = W ( ) ( r ) ∈ H ( curl ) (14b)2-form : = W ( ) ( r ) ∈ H ( div ) (14c)3-form : = W ( ) ( r ) ∈ L (14d)The Hodge operator relates the corresponding field and fluxquantities using the correct constitutive parameters, as well asprovides a mapping between the primal and dual grids . E , J B ρ FIG. 3: Two dimensional diagram of where field and sourcequantities are defined on a simplical mesh.
2. Discrete and self consistent Gauss’ and Ampere’s Laws
In typical solutions to Maxwell equations, there is an im-plicit assumption that the solution to the two curl equationenforces Gauss’ laws together with the equation of continuity.While this is true in the continuous domain, one needs to han-dle things with care in the discrete domain. Specifically, thebasis functions used to represent or measure the fields, fluxes,currents, and charges need to be consistent with each other.Let L ∇· ◦ denote a discrete divergence operator and L ∇× ◦ denote a discrete curl operator, whose definitions are depen-dent on the basis functions. With no loss in generality, let’sdenote the basis functions, N i ( r ) and ξ i ( r ) , used to representand measure the fields and fluxes, and charges, respectively.They belong to the appropriate space of functions dependingon which quantity is being measured as determined by (14)and have local support over a simplex or cell. They also haverelationships determined by the de-Rham complex shown infigure 2. With no loss of generality, let us define the vectorsof quantities that are measured by our candidate basis func-tions. For instance, on can define d i ( t ) = (cid:104) N i ( r ) , ε E ( t , r ) (cid:105) (15a) ( L ∇× ◦ h ( t )) i = (cid:104) N i ( r ) , ∇ × µ − B ( t , r ) (cid:105) (15b) ( L ∇· ◦ d ( t )) i = (cid:104) ξ i ( r ) , ∇ · ε E ( t , r ) (cid:105) (15c) j i ( t ) = (cid:104) N i ( r ) , J ( t , r ) (cid:105) (15d)Next, let us define d ( t ) = [ d ( t ) , d ( t ) , . . . , d N γ ( t )] , h ( t ) =[ h ( t ) , h ( t ) , . . . , h N α ( t )] , and j ( t ) = [ j ( t ) , j ( t ) , . . . , j N β ( t )] ,which can be used to approximate the continuous functions,such as D ( t , r ) = ∑ N γ i = d i ( t ) N i ( r ) , where N γ is the number ofdegrees of freedom; as we will see, depending on the formula-tion it can either the number of edges or faces. This definitionis the same for N α and N β . The degree of freedom vectorsubrics for Charge Conserving Current Mapping in FEM-PIC 5may require a Hodge operator to represent the function in thecorrect space. For solutions to be consistent, the systems ofequations should be such that the divergence of Ampere’s lawyields the time derivative of Gauss’ law, which together withinitial conditions will yield Gauss’ law. While this is simplystated, there are a number of subtle nuances. In the discreteworld, it follows that the same should be valid. Specifically,the spatial basis sets and the time integrator used has to besuch that these relationships are satisfied. Note, these condi-tions must be satisfied exactly because only the curl equations,Faraday’s law and Ampere’s law, are explicitly solved. Gauss’law must be satisfied indirectly through the correct choice ofspatial basis sets and time integration.As our goal is to ensure the both Gauss’ laws and the equa-tion of continuity are satisfied, it follows that there has to be arelationship between ∇ ξ i ( r ) and N i ( r ) . We begin by consid-ering the divergence of discrete Ampere’s law; specifically, ∂ t L ∇· ◦ d ( t ) = L ∇· ◦ L ∇× ◦ h ( t ) − L ∇· ◦ j ( t ) (16)From the above, one can deduce the following: Condition 1: Consistent spatial basis sets:
The basis setsshould be such that the discrete divergence of the dis-crete rotation of the magnetic field is zero. In effect, thediscrete differential operators should maintain the sameproperties as their continuous counterpart. To serve asthe link between Ampere’s law and the continuity equa-tion, L ∇· ◦ L ∇× ◦ h ( t ) = . (17)Though the discrete differential operators can be de-fined using properties of the mesh, they must be con-sistent with the discretized Faraday’s and Gauss’ law. Condition 2: Discrete divergence:
Let ξ ( r ) denote func-tions used to measure the charge on the mesh. Then ∇ ξ i ( r ) should span the irrorational component ofboth D ( t , r ) and J ( t , r ) . In effect, span { ∇ ξ i ( r ) } = D ( t , r ) or J ( t , r ) for ∀ t . In other words, this represen-tation is complete. The need for this condition arisesfrom measuring both Gauss’ law and the equation ofcontinuity. Consider measuring Gauss’ law via ( L ∇· ◦ d ( t )) i = (cid:90) Ω d r ξ i ( r ) ρ ( r ) = (cid:90) Ω d r ξ i ( r ) ∇ · D ( t , r ) (cid:90) Ω d r ξ i ( r ) ρ ( r ) = − (cid:90) Ω d r ∇ ξ i ( r ) · D ( t , r ) + (cid:90) ∂ Ω d r ξ i ( r ) ˆ n · D ( t , r ) . (18)The preceding equation mathematically exemplifies thiscondition. An almost identical equation can be derivedfor the equation of continuity. It can be shown bound-ary terms in vanish interior to the domain and providedthe correct boundary condition is imposed on ˆ n · D ( t , r ) on ∂ Ω , Gauss’ law is implicitly satisfied together withcondition 3. Condition 3: Discrete Time Integration:
Finally, discretetime time stepping used to update the electric flux den-sity in Ampere’s equation and the path integral to used to evaluate the total current at any instant of time shouldbe consistent with each other. This will ensure that thesolution to Ampere’s law agrees with Gauss’ law. Togain more insight, we start by noting that when evolv-ing Maxwell’s equations over-time, we are in effect in-tegrating the current; however, the time integral of thecurrent is obtained via a path integral in (13). It fol-lows, that if one were to take a discrete divergence ofAmpere’s laws, Condition 1 would yield ∂ t L ∇· ◦ d ( t ) = − L ∇· ◦ j ( t ) (19)Provided Condition 2 is satisfied, it follows that choos-ing a consistent rule for integrating the charge along theparticle path and integrating Ampere’s law will ensureconservation of charge. IV. SATISFACTION OF CONDITIONS WITH VARIOUSEM-PIC METHODS
Next, we examine three different formulations of EM-PIC;FDTD-PIC, EB-MFEM-PIC, and DH-MFEM-PIC in light ofthe conditions discussed thus far. While FDTD-PIC is wellknown, EM-MFEM-PIC and DH-MFEM-PIC place dominantimportance on Faraday’s and Ampere’s laws, respectively. Asa result, the dominant variables are represented on the primalgrid, the others reside on dual grid. In what follows, we seekto show how each formulation satisfies the three conditions.
A. FDTD-PIC
FDTD-PIC is the combination of the FDTD method withPIC. It is discretized using staggered, structured grids in timeand space, as seen in Figure 4 for two dimensions. For ourdiscussion, the primal grid contains degrees of freedom forthe electric field and magnetic flux density. Likewise, the dualgrid has the degrees of freedom defined for the magnetic fieldand the electric flux density . Particles are mapped to thegrid using the 3-forms that measure the charge density. Inwhat follows, we will show that the particle weighting schemedefined by Langdon maps trivially onto the three conditionsdefined earlier.
1. Satisfaction of Condition 1
The discrete differential operators L ∇× ◦ and L ∇· ◦ on thedual grid are defined as L ∇× ◦ = (cid:40) ± n j × ˆ l i points into (out of) the face0 otherwise (20)and L ∇· ◦ = (cid:40) ± n j points out of (into) cell c i L ∇× ◦ ∈ R N f × N e , where N f is the number offaces and N e is the number of edges. The operator L ∇· ◦ ∈ R N c × N f where N c is the number of volumes.Consider an example mesh in Figure 5. Assume the edgeon the corners of the cell all point out of the page and thenormal to the edges point out of the cell. The correspondingcurl matrix definition would be [ L ∇× ◦ ] = − − −
10 0 − , (22)The divergence matrix would be [ L ∇· ◦ ] = (cid:2) (cid:3) (23)It is evident that in L ∇· ◦ L ∇× ◦ , each entry in the resultingvector is zero, satisfying (17).
2. Satisfaction of Condition 2
The basis function ξ ( r ) in (18) is a 3-form, a constant inthis formulation. Using the properties of 3-forms, (18) re-duces to (cid:90) Ω d r ξ ( r ) ∇ · X ( r ) = (cid:90) ∂ Ω d r ˆ n · X ( r ) ξ ( r ) . (24)Consider the matrix form of L div ◦ based on Figure 5 L ∇· ◦ = (cid:20) − (cid:21) . (25) FIG. 5: Discrete curl and divergence on the dual grid. Thecurl is defined using the edges ˆ l i that come out of the page.The divergence is defined using the face normals ˆ n i .A volume integral would be a summation over the cells of amesh, N v ∑ i = [ L ∇· ◦ ] i j = (cid:2) (cid:3) (26)This would effect the surface integral in (24). Therefore, thediscrete continuity equation is satisfied exactly and consis-tency condition is fully satisfied.
3. Satisfaction of Condition 3
Consider (24) in conjunction with (19). Using a leapfrogtime marching scheme, the discrete continuity equation wouldbe [ L ∇· ] d n + − d n ∆ t = − q ∆ t [ L ∇· ] (cid:90) ∂ Ω d r (cid:90) r n + r n d ˜ r ˆ n · p ( ˜ r ) , (27)where ∆ t is the time step size, p ( ˜ r ) is the path of the parti-cle. The integral computes the amount of charge that passesthrough a face. The divergence operator collects the contri-butions of each face, giving the net change of charge in eachcell, as demonstrated in figure 6. This is equivalent to the cur-rent mapping in for particles with finite shape. If the particleshape is infinitesimally small, it is equivalent to nearest-grid-point particle weighting. B. EB-MFEM-PIC Formulation
Next, we analyze EB-MFEM-PIC as defined in . This fi-nite element method can be seen as a generalization of FDTD.Like FDTD, it discretizes Maxwell’s equations using a primal-dual grid system, with the field and flux quantities defined inubrics for Charge Conserving Current Mapping in FEM-PIC 7FIG. 6: FDTD 2D particle current mapping with particles offinite size.the same manner. However, the grid itself can be unstructured,using non-cubic cells. The formulation for EB-MFEM can befound in appendix VIII C 2. We have chosen to demonstratethe validity of the arguments in Section III using brick ele-ments as it is akin to FDTD-PIC. We defer the discussion ontetrahedral elements to the next section.
1. Satisfaction of Condition 1
The discrete differential operators L ∇× ◦ and L ∇· ◦ are de-fined on the primal grid as L ∇× ◦ = (cid:40) ± n j × ˆ l i points into (out of) the face0 otherwise (28)and L ∇· ◦ = (cid:40) ± j points into (out of) node i L ∇× ◦ ∈ R N e × N f . The operator L ∇· ◦ ∈ R N n × N e where N n is the number of nodes. Consider the example meshin Figure 7. The corresponding curl matrix definition wouldbe L ∇× ◦ = (cid:2) − − (cid:3) T . (30)The divergence matrix would be L ∇· ◦ = − − − − (31)It is apparent that each entry in the resulting L ∇· ◦ L ∇× ◦ iszero, satisfying (17).
2. Satisfaction of Condition 2
The basis function for the current density spans the gradientof the basis function for the charge density. The divergence FIG. 7: FEM 2D curl and divergence mesh. The curl isdefined using the edges ˆ l i around face center ˆ n i . Thedivergence is defined with edges that point into node g i .operator L ∇· ◦ sums the edges connected to a node to satisfy(18). For perfect electrical conductor (PEC) boundary condi-tions, the current density normal to the face is zero, eliminat-ing the surface integral term. This can be seen by summing therows of (31). For non-PEC boundaries, the domain is paddedwith ghost cells where the boundary nodes are not included.The non-terminating edges then make up the surface integralterm when L ∇· ◦ is applied. Therefore, (18) is satisfied.
3. Satisfaction of Condition 3
For each cell a particle visits, the corresponding currentdensity is projected onto the edges of that cell. Using (18),the discretized continuity equation that satisfies (19) is [ L ∇· ] d n + − d n ∆ t = − q δ t [ L ∇· ] (cid:90) Ω d r (cid:90) r n + r n d ˜ r ∇ W · p ( ˜ r ) . (32)This is equivalent to particle-in-cell area weighting in FDTDand the method used in . Due to the relationship betweenthe 0-form on the primal grid and the 3-form on the dual grid,this is also equivalent to the approach defined in for finiteparticle shapes. C. DH-MFEM-PIC Formulation
The EB-MFEM-PIC formulation defined in IV B can bereformulated to such that the primal grid has the degrees offreedom associated with the magnetic field and electric fluxubrics for Charge Conserving Current Mapping in FEM-PIC 8FIG. 8: FEM 2D particle current mapping with pointparticles. The particle current is mapped to each edge thatencloses a cell along the particle path.density. The dual grid contains the degrees of freedom forthe electric field and magnetic flux density, still defined us-ing the edges and faces, respectively. The sources are repre-sented on the primal grid using 3-forms for the charge densityand 2-forms for the current density. The arguments for thisformulation satisfying the three conditions are similar to thatpresented in section IV A. Note, the principal ramification ofthis approach is that it is a natural framework for modelingsources as both the currents and electric flux density are mod-eled using the same function space, and the fields, currents andcharges lie in the primal grid. This si unlike EM-MFEM-PICformulation where one needs transformations between primaland dual spaces. The ramification of these mappings becomemore consequential in higher order schemes–a topic that willbe visited in a later paper.
1. Satisfaction of Condition 1
The discrete differential operators L ∇× ◦ and L ∇· ◦ are de-fined on the primal grid, but have the same definitions as (20)and (21). Consider the example mesh in Figure 9. Assume theedge on the corners of the cell all point out of the page and thenormal to the edges point out of the cell. The correspondingcurl matrix definition for lowest order would be [ L ∇× ◦ ] = − − − , (33)The lowest order divergence matrix would be [ L ∇· ◦ ] = (cid:2) (cid:3) (34)Analyzing L ∇· ◦ L ∇× ◦ , each entry in the resulting vector iszero, satisfying (17).
2. Satisfaction of Condition 2
For lowest order, the continuity equation reduces to (24).Consider the matrix form of L ∇· ◦ based on figure 10 L ∇· ◦ = (cid:20) − (cid:21) . (35) FIG. 9: Discrete curl mesh on the dual grid. The curl isdefined with the outward pointing edge ˆ l i .FIG. 10: FEM divergence mesh. The divergence is definedusing face normals ˆ n i that point out of cell c i .A volume integral would be a summation over the cells of amesh, N v ∑ i = [ L ∇· ◦ ] i j = (cid:2) (cid:3) (36)This would effect the surface integral in (24). Therefore, thediscrete continuity equation is satisfied exactly and consis-tency condition is fully satisfied.
3. Satisfaction of Condition 3
Again, the satisfaction of condition three relies on the cor-rect definition of current density on the faces of each cell.For consistency of the change of charge density in cells, theamount of charge that passes through the face must be mea-sured exactly. This is equivalent to what is effectively doneubrics for Charge Conserving Current Mapping in FEM-PIC 9FIG. 11: FEM particle mapping with finite shape particles.The particle is mapped to cell c and c at first and ends incell c with current mapped at on the marked faces. Theparticle shape does not affect the conservation properties ofthis method.in through different techniques. Consider the mapping ofcurrent defined by the particle in Figure 11. Here, though theparticle moves at an angle with respect to the face, the currentdensity as seen by the mesh moves in the direction of the facenormal. The solver in effect sees a charge density move fromcell c to c and the total charge density in c moves to c bytime step n + V. CONSEQUENCES OF INCONSISTENCY
In the previous sections, the three rule framework was de-fined and applied to different formulations of EM-PIC. Thereferences to previous works show how these rules encom-passes well known EM-PIC formulations. Next, we will ex-amine ramifications of violating these three rules. By exam-ining potential pitfalls and how it manifests in results, one canmore easily diagnose errors in new EM-PIC schemes.
A. Violation of Condition 1
First, consider violations of the first condition. Failing tosatisfy (17) can arise from choosing the wrong basis sets todiscretize (2). One way to do so is to choose basis sets that donot obey the de-Rahm complex. For example, if one solvesboth Ampere’s law and Faraday’s law on the same grid, stag-gered in time. Like in FDTD, we consider representing theelectric field and magnetic field; however, this time it is onthe same grid. Assume the fields are represented by 1-formsand test (2a) and (2b) with the same basis functions. The dis- f true f consistent f inconsistent TABLE I: First four rectangular cavity modes for consistentand inconsistent FEM solver in MHzcretized system is h n + = h n + − ∆ t [ L (cid:48) ∇× ◦ ] e n (37a) e n + = e n + + ∆ t ([ L (cid:48) ∇× ◦ ] h n + − j n + ) , (37b)where L (cid:48) ∇× ◦ is a modified curl operator, as the previouslydefined operator has the wrong dimensions. If this schemeis done with a structured mesh in FDTD, one would have ahigher error as the curl would be represented using forward orbackward difference rather than a central difference. Insteadof the recurrence scheme found in (51b), one would have E n + x , ( i + , j , k ) = E nx , ( i + , j , k ) + ∆ t ε ∆ y ( H n + z , ( i + , j + , k ) − H n + z , ( i + , j , k ) ) − ∆ t ε ∆ z ( H n + y , ( i + , j , k + ) − H n + y , ( i + , j , k ) ) − ∆ t ε J n + / x , ( i + , j , k ) (38)if a forward difference was used.The effects of this are more clear from an FEM perspective.Consider trying to find modes of a perfect electric conductorrectangular cavity using the generalized eigenvalue problemformed by the time harmonic form of (39) (cid:20) [ L ∗ ∇× ◦ ] − [ L ∗ ∇× ◦ ] T [ (cid:63) µ − ] (cid:21) (cid:20) be (cid:21) = Λ (cid:20) [ I ] [ (cid:63) ε ] (cid:21) (cid:20) be (cid:21) (39)where the eigenvalues in Λ = j ω . The dimensions of the cav-ity are 0 . × . × [ L ∗ ∇× ◦ ] =[ L ∇× ◦ ] and [ (cid:63) µ − ] is as defined in (39). However, in the in-consistent scheme we have defined, [ L ∗ ∇× ◦ ] i , j = (cid:90) d Ω W ( ) i · ∇ × W ( ) j (40)and [ (cid:63) µ − ] = [ (cid:63) ε ] . When the inconsistent scheme is used,spurious modes are found that do not correspond to the truemodes of the cavity as seen in table I. Using this inconsistentscheme could permit incorrect solutions into the field solver,invalidating results. B. Violation of Condition 2
Next, we consider the repercussions of violating Condition2. Consider the linear mapping scheme typically found inubrics for Charge Conserving Current Mapping in FEM-PIC 10 ∂ t ρ ∇ · J ∂ ρ + ∇ · J node 1 -0.100 -0.050 -0.150node 2 0.300 -0.850 -0.550node 3 -0.300 0.850 0.550node 4 0.100 0.050 0.150 TABLE II: Continuity equation error for linear mapping
Analytic 1 pt Quad ErrorEdge 1 1.248 1.248 -4.441e-16Edge 2 0.505 0.505 0.000Edge 3 0.248 0.248 0.000
TABLE III: Error of Analytic Integration vs Quadrature Rulefor Linear Pathearly PIC literature where the weighting function is definedas W ( r − r p ( t )) = ∏ n ∈{ x , y , z } − | r n , g − r n , p | ∆ n . (41)When the particle is mapped to the nodes as the charge den-sity, r n , g is the x , y , z − coordinate of the node. When mappedto the edges as the current density, r n , g is the x , y , z − coor-dinate of the center of the edge. It is well known that thismapping does not satisfy the continuity equation. For exam-ple, consider the particle motion in Figure 8 where the particlemoves from ( . , . ) to ( . , . ) . Table II shows the error in thecontinuity equation using this mapping. The result is an in-consistency in the solved Ampere’s law and measured Gauss’law due to spurious charge seen by the field solver as the nextmapping is done. This mapping scheme does not appropri-ately satisfy (18), as the combination of basis functions forthe current density at the node is not exactly equal to the gra-dient of the charge density basis function. Despite being in thesame direction, the basis function defined by the linear map-ping for the current density is different than the basis functionsused to measure the fields. Therefore, the de-Rham complexis not truly satisfied and the spatial link between Gauss’ lawand Ampere’s law is no longer preserved. C. Violation of Condition 3
Next, we examine how failing to correctly integrate the spa-tial integral in (19) results in a violation of Gauss’ law and thecontinuity equation due to spurious charge on the grid. Con-sider the movement of a particle on a linear path that remainsin a single cell using a 0-form basis set to measure the chargedensity. For a linear path, the path integral can be exactlycomputed either using an analytic scheme or an one pointintegration rule, as seen in Table III. This translates to ma-chine precision error in the continuity equation in Table IV.Similarly, consider the case when a particle crosses at leastone face on its path. When the path is in multiple cells, anerror occurs when the path integral is not calculated in each Q n + Q n Q n + − Q n ∆ t ∇ · J Q n + − Q n ∆ t + ∇ · J Node 1 0.284 0.459 -1.753 1.753 0.000Node 2 0.550 0.450 1.000 -1.000 -4.440e-16Node 3 0.166 0.091 0.753 -0.753 0.000
TABLE IV: Continuity Equation Error for Linear Path Q n + − Q n ∆ t + ∇ · j n + / nc Q n + − Q n + α ∆ t + Q n + α − Q n ∆ t + ∇ · j n + / cc Node 1 1.008 0.000Node 2 0.394 -4.441e-16Node 3 -0.008 2.220e-16Node 4 -1.394 -4.441e-16
TABLE V: Continuity Equation Error for Multiple Cellscell, as seen in Table V. This is an inconsistency in (19) asthough the time integral may be computed correctly, parts ofthe path integral are missing.
D. Results for Consistent DH-MFEM-PIC
Finally, we examine results when all three conditionsare satisfied. We have chosen to illustrate these us-ing DH-MFEM-PIC as results for EB-MFEM-PIC areavailable . First, we consider a particle beam in a 0.1 mPEC cylindrical cavity with a radius of 0.02 m. The particlebeam has a voltage of 7.107 kV and current of 0.25A, simu-lated by inserting ten particles per time step. The starting ve-locity of the particles entering the tube is 5 × km/s with abeam width of 0.008 m. In this test, when the three conditionsare satisfied, the beam smoothly spreads from a uniform dis-tributions. Using a DH-based EM-PIC formulation and pointparticles for particle shape, by following section IV C, a self-consistent charge conserving EM-PIC scheme was created.Point particles are used because without introducing errors,it removes the need to compute the volume of a particle pass-ing through a face, which can be extremely difficult in threedimensions, The error in the continuity equation and Gausslaw is found in Figure 13, showing that satisfying conditionstwo and three is possible with a nearest grid point algorithm inFEM. The trend of the spreading beam is compared to that ofan EB-MFEM-PIC and a FDTD based EM-PIC formulationin Figure 14.Lastly, we examine the expansion of a plasma ball. A spher-ical plasma ball is given a Gaussian distribution of ions andelectrons in a geometry with first order absorbing boundaryconditions. The experiment was simulated with the DH-FEM Analytic 1 pt Quad 1 pt error 2 pt Quad 2 pt errorEdge 1 1.253 1.250 -0.003 1.253 -6.661e-16Edge 2 0.500 0.500 0.000 0.500 2.220e-16Edge 3 0.253 0.250 0.003 0.253 0.000
TABLE VI: Error of Analytic Integration vs Quadrature Rulefor Quadratic Pathubrics for Charge Conserving Current Mapping in FEM-PIC 11FIG. 12: Particle beam comparison of DH -MFEM-PIC andFDTD-PIC. (a)(b) FIG. 13: Error in the continuity equation (13a) and Gauss’law (13b) for the particle beam case.formulation with a point particle shape. A combination of Sr + ions at 1K and electrons at 100K are placed at a densityof 5 × particles per cubic meter such that the system is ini-tially charge neutral. Further details on the experiment setupcan be found in and derivations of the analytic results in .As is evident from Figure 15, we achieve excellent agreementwith the analytic results. Again, because all the conditions aresatisfied, the continuity equation and errors in Gauss’ law arekept to machine precision. FIG. 14: Particle Beam for FDTD, EB , and DH formulationFIG. 15: Electron density of expanding adiabatic plasma ballusing DH -MFEM-PIC compared to analytic expecteddensities. VI. SUMMARY AND CONCLUSIONS
In this work, we have presented a framework to define self-consistent EM-PIC schemes. By considering the space offunctions used to represent fields, fluxes, and sources, one canensure that all the relationships between Maxwell’s equationsare preserved in a consistent manner. This has been demon-strated for FDTD and two formulations of FEM, and shownhow it generalizes and relates to well known previous work.The framework is easily related to other field solvers not men-tioned directly provided one considers the relationship of thebasis functions. Future work will apply this framework tohigher order basis functions and particle paths.
VII. ACKNOWLEDGEMENTS
This work was sponsored by the US Air Force ResearchLaboratory under contracts FA8650-19-F-1747 and FA8650-20-C-1132. This work was also supported by the Departmentof Energy Computational Science Graduate fellowship undergrant DE-FG02-97ER25308. The authors would also like tothank the HPCC Facility at Michigan State University and theDoD SMART scholarship.ubrics for Charge Conserving Current Mapping in FEM-PIC 12
VIII. APPENDIXA. Calculus on Distribution functions
The following must be understood in the sense of distri-butions. Specifically, we consider derivatives and integral ofthe delta function. These appear in (8) onwards where onerepresents the spatial variation of PSDF using S ( r − r p ( t )) = S ( r ) (cid:63) s δ ( r − r p ( t )) . Calculus on delta function is critical fordeveloping insight into implementing equations of continuityand Gauss’ law on the grid. To do so, we start with a transformto k -space to prove (11) and (13).
1. Time derivatives
Using r p ( t ) = r p + (cid:82) t d τ v p ( τ ) , the time derivative of δ ( r − r p ( t )) can be written as ∂ t δ ( r − r p ( t )) = ∂ t δ ( r − r p ) (cid:63) δ (cid:18) r − (cid:90) t d τ v p ( τ ) (cid:19) = δ ( r − r p ) (cid:63) ∂ t δ (cid:18) r − (cid:90) t d τ v p ( τ ) (cid:19) (42)Next, we expliot δ ( r − a ) = ( π ) (cid:90) ∞ − ∞ d k e j k · ( r − a ) (43)to yield ∂ t δ (cid:18) r − (cid:90) t d τ v p ( τ ) (cid:19) = ( π ) (cid:90) ∞ − ∞ d k ∂ t e j k · ( r − (cid:82) t d τ v p ( τ ) )= − ∇ · v p ( t ) δ (cid:18) r − (cid:90) t d τ v p ( τ ) (cid:19) . (44)Therefore ∂ t δ ( r − r p ( t )) = − δ ( r − r oi ) (cid:63) ∇ · v p ( t ) δ (cid:18) r − (cid:90) t d τ v p ( t ) (cid:19) = − ∇ · [ v p ( t ) δ ( r − r p ( t ))] (45) B. Time integrals
Next, starting from the time integral of (11) and using thedefinition of current density in equation (9b), (cid:90) t d τ ∇· v p ( τ ) S ( r − r p ( τ )) = S ( r ) (cid:63) s (cid:90) t d τ ∇· v p ( τ ) δ (cid:18) r − (cid:90) τ d τ (cid:48) v p ( τ (cid:48) ) (cid:19) (46) Again, using (43) and ∂ t ˜ r p ( t ) = v p ( t ) , (cid:90) t d τ ∇ · v p ( τ ) δ (cid:18) r − (cid:90) τ d τ (cid:48) v p ( τ (cid:48) ) (cid:19) = (cid:90) t d τ ( π ) ∇ · (cid:90) ∞ − ∞ d kv p ( τ ) e j k · ( r − (cid:82) τ d τ (cid:48) v p ( τ (cid:48) ) )= ∇ · ( π ) (cid:90) ∞ − ∞ d k e j k · r (cid:90) t d τ v p ( τ ) e − j k · (cid:82) τ d τ (cid:48) v p ( τ (cid:48) ) = ∇ · ( π ) (cid:90) ∞ − ∞ d k e j k · r (cid:90) r p ( t ) r p ( ) d ˜ r p e − j k · ˜ r p (47)Interchanging the order of the integrals in the above equationresults in (cid:90) t d τ ∇ · v p ( τ ) S ( r − r p ( τ )) = ∇ · (cid:90) r p ( t ) r p ( ) d ˜ r S ( r − ˜ r ) (48)This equation can be used to define the charge density in termsof the particle’s path; it can also be used to define the timeintegral of the current density. C. Maxwell Solvers
The following sections are presented purely for complete-ness and to conform with material presented earlier. Wepresent details of the three Maxwell solvers used here; FDTD,EB-MFEM, and DH-MFEM.
1. FDTD
Define the rooftop function Λ i ( r ) Λ i ( r ) = (cid:40) − | ˆ i · r | ∆ i < | ˆ i · r | ≤ ∆ i i = x , y , z . The basis functions to discretize Maxwell’sequations for FDTD are1-form : = u i Λ i + ( r − l m ) Λ i + ( r − l m ) (50a)2-form : = u i Λ i ( r − f m ) (50b)with u i = ˆ x , ˆ y , ˆ z , l m is the midpoint of each edge, and f m is themidpoint of each face. Testing Faraday’s law with the normalto cell faces and Ampere’s law is with the edge unit vectorsyields (cid:104) u i δ (( r − f m ) · ˆ i ) , ∂ t µ H ( t , r ) (cid:105) = −(cid:104) u i δ (( r − f m ) · ˆ i ) , ∇ × E ( t , r ) (cid:105) ∂ t H i ( t , f m ) = − µ − u i · ∇ × E ( t , f m )= µ − ∂ i + E i + ( t , f m ) − µ − ∂ i + E i + ( t , f m ) (51a)ubrics for Charge Conserving Current Mapping in FEM-PIC 13 (cid:104) u i δ (( r − f m ) · ˆ i ) , ∂ t ε E ( t , r ) (cid:105) = − (cid:104) u i δ (( r − l m ) · ˆ i ) , ∇ × H ( t , r ) (cid:105) − (cid:104) u i δ (( r − l m ) · ˆ i ) , J ( t , r ) (cid:105) ∂ t E i ( l m ) = − ε − u i · ∇ × H ( t , l m ) − ε − u i · J ( t , l m )= ε − ∂ i + H i + ( t , l m ) − ε − ∂ i + H i + ( t , l m ) − ε − J i ( t , l m ) . (51b)The addition on index i follows cycle { x , y , z } . This testing isequivalent to the usual presentation of the FDTD recurrencerelation. For example, consider i = x in (51a) would yield H n + x , ( i , j + , k + ) = H n − x , ( i , j + , k + ) − ∆ t µ ∆ y ( E nz , ( i , j + , k + ) − E nz , ( i , j , k + ) )+ ∆ t µ ∆ z ( E ny , ( i , j + , k + ) − E ny , ( i , j + , k ) ) . (52)Likewise, i = x in (51b) would yield E n + x , ( i + , j , k ) = E nx , ( i + , j , k ) + ∆ t ε ∆ y ( H n + z , ( i + , j + , k ) − H n + z , ( i + , j − , k ) ) − ∆ t ε ∆ z ( H n + y , ( i + , j , k + ) − H n + y , ( i + , j , k − ) ) − ∆ t ε J n + / x , ( i + , j , k ) . (53)
2. EB-MFEM Formulation
In this formulation, E ( t , r ) ≈ ∑ N e i = e i ( t ) W ( ) i ( r ) and and B ( t , r ) ≈ ∑ N f i = b i ( t ) W ( ) i ( r ) . Faraday’s law is measured withthe 2-form basis function and Ampere’s law with the 1-formbasis function to yield ∂ t b = [ L ∇× ◦ ] e (54a) ∂ t e = c [ L ∇× ◦ ] T [ (cid:63) µ − ] b − ε − j (54b)assuming the domain is freespace with the speed of light c = √ ε µ − . The Hodge operator matrix is defined as [ (cid:63) µ − ] i j = (cid:90) Ω d Ω W ( ) i ( r ) · W ( ) j ( r ) . (55)The degrees of freedom for the current are defined as j i = N p ∑ n = q n ∆ t (cid:90) r n + r n d p n · W ( ) i ( r ) (56)where p n is the particle path for the n -th particle.
3. DH-MFEM Formulation
In this formulation, H ( t , r ) ≈ ∑ N e i = h i ( t ) W ( ) i ( r ) and D ( t , r ) ≈ ∑ N f i = d i ( t ) W ( ) i ( r ) . Ampere’s law is measured with the 2-form basis function and Faraday’s law with the 1-formbasis function to yield ∂ t h = c [ L ∇× ◦ ][ (cid:63) ε − ] d (57a) ∂ t d = [ L ∇× ◦ ] h − j . (57b)The Hodge operator matrix is defined as [ (cid:63) ε − ] i j = (cid:90) Ω d Ω W ( ) i ( r ) · W ( ) j ( r ) . (58)The degrees of freedom for the current are defined as j i = N p ∑ n = q n ∆ t (cid:90) r n + r n d p n ˆ n i · J n ( r ) (59)where the particle path for the n -th particle and is only evalu-ated on the surface. J. Deca, A. Divin, G. Lapenta, B. Lembège, S. Markidis, and M. Horányi,“Electromagnetic particle-in-cell simulations of the solar wind interac-tion with lunar magnetic anomalies,” Physical review letters , 151102(2014). R. Lichters, R. E. Pfund, and J. Meyer-ter Vehn, “Lpic++. a parallel one-dimensional relativistic electromagnetic particle-in-cell code for simulat-ing laser-plasma-interaction,” Tech. Rep. (Max-Planck-Institut fuer Quan-tenoptik, 1997). S. P. Gary, S. Saito, and H. Li, “Cascade of whistler turbulence: Particle-in-cell simulations,” Geophysical research letters (2008). Y. Chen and S. E. Parker, “Electromagnetic gyrokinetic δ f particle-in-cellturbulence simulation with realistic equilibrium profiles and geometry,”Journal of Computational Physics , 839–855 (2007). P. Mardahl, H.-J. Lee, G. Penn, J. Wurtele, and N. J. Fisch, “Intense laserpulse amplification using raman backscatter in plasma channels,” PhysicsLetters A , 109–116 (2002). S. Chakrabarti, J. Martin, J. Pearson, and R. Lewis, “Developing antimat-ter containment technology: Modeling charged particle oscillations in apenning-malmberg trap,” (2003). R. A. Fonseca, L. O. Silva, F. S. Tsung, V. K. Decyk, W. Lu, C. Ren, W. B.Mori, S. Deng, S. Lee, T. Katsouleas, et al. , “Osiris: A three-dimensional,fully relativistic particle in cell code for modeling plasma based acceler-ators,” in
International Conference on Computational Science (Springer,2002) pp. 342–351. C. S. Meierbachtol, A. D. Greenwood, J. P. Verboncoeur, and B. Shanker,“Conformal electromagnetic particle in cell: A review,” IEEE Transactionson Plasma Science , 3778–3793 (2015). G. Chen, L. Chacón, and D. C. Barnes, “An energy-and charge-conserving,implicit, electrostatic particle-in-cell algorithm,” Journal of ComputationalPhysics , 7018–7036 (2011). G. Chen and L. Chacón, “An energy-and charge-conserving, nonlinearlyimplicit, electromagnetic 1d-3v vlasov–darwin particle-in-cell algorithm,”Computer Physics Communications , 2391–2402 (2014). I. V. Sokolov, “Alternating-order interpolation in a charge-conservingscheme for particle-in-cell simulations,” Computer Physics Communica-tions , 320–328 (2013). J. Villasenor and O. Buneman, “Rigorous charge conservation for lo-cal electromagnetic field solvers,” Computer Physics Communications ,306–316 (1992). A. B. Langdon, “On enforcing gauss’ law in electromagnetic particle-in-cell codes,” Computer Physics Communications , 447–450 (1992). M. Kraus, K. Kormann, P. J. Morrison, and E. Sonnendrücker, “Gempic:Geometric electromagnetic particle-in-cell methods,” Journal of PlasmaPhysics (2017). D.-Y. Na, H. Moon, Y. A. Omelchenko, and F. L. Teixeira, “Local, ex-plicit, and charge-conserving electromagnetic particle-in-cell algorithm onunstructured grids,” IEEE Transactions on Plasma Science , 1353–1362(2016). ubrics for Charge Conserving Current Mapping in FEM-PIC 14 G. Jacobs and J. S. Hesthaven, “High-order nodal discontinuous galerkinparticle-in-cell method on unstructured grids,” Journal of ComputationalPhysics , 96–121 (2006). B. Marder, “A method for incorporating gauss’ law into electromagneticpic codes,” Journal of Computational Physics , 48–55 (1987). J. P. Boris, “Relativistic plasma simulation-optimization of a hybrid code,”in
Proc. Fourth Conf. Num. Sim. Plasmas (1970) pp. 3–67. B. A. Shadwick, J. C. Bowman, and P. Morrison, “Exactly conservativeintegrators,” SIAM Journal on Applied Mathematics , 1112–1133 (1998). O. Buneman, “Dissipation of currents in ionized media,” Physical Review , 503 (1959). J. Dawson, “One-dimensional plasma model,” The Physics of Fluids ,445–459 (1962). J. W. Eastwood, “The virtual particle electromagnetic particle-meshmethod,” Computer Physics Communications , 252–266 (1991). M. C. Pinto, S. Jund, S. Salmon, and E. Sonnendrücker, “Charge-conserving fem–pic schemes on general grids,” Comptes RendusMecanique , 570–582 (2014). H. Moon, F. L. Teixeira, and Y. A. Omelchenko, “Exact charge-conservingscatter–gather algorithm for particle-in-cell simulations on unstructuredgrids: A geometric perspective,” Computer Physics Communications ,43–53 (2015). E. Tonti, “Finite formulation of the electromagnetic field,” Progress in elec-tromagnetics research , 1–44 (2001). G. A. Deschamps, “Electromagnetics and differential forms,” Proceedingsof the IEEE , 676–696 (1981). L. Demkowicz, P. Monk, L. Vardapetyan, and W. Rachowicz, “De rhamdiagram for hp finite element spaces,” Computers & Mathematics with Ap-plications , 29–38 (2000). D. Arnold, R. Falk, and R. Winther, “Finite element exterior calculus: fromhodge theory to numerical stability,” Bulletin of the American mathematicalsociety , 281–354 (2010). A. Bossavit, “Whitney forms: A class of finite elements for three-dimensional computations in electromagnetism,” IEE Proceedings A (Phys-ical Science, Measurement and Instrumentation, Management and Educa-tion, Reviews) , 493–500 (1988). R. Rieben, G. Rodrigue, and D. White, “A high order mixed vector finiteelement method for solving the time dependent maxwell equations on un-structured grids,” Journal of Computational Physics , 490–519 (2005). J.-F. Lee, R. Lee, and A. Cangellaris, “Time-domain finite-element meth-ods,” IEEE transactions on antennas and propagation , 430–442 (1997). A. Bossavit and L. Kettunen, “Yee-like schemes on staggered cellular grids:A synthesis between fit and fem approaches,” IEEE Transactions on Mag-netics , 861–867 (2000). P. Monk et al. , Finite element methods for Maxwell’s equations (OxfordUniversity Press, 2003). D. N. Arnold, R. S. Falk, and R. Winther, “Finite element exterior calculus,homological techniques, and applications,” (2006). F. L. Teixeira and W. Chew, “Lattice electromagnetic theory from a topo-logical viewpoint,” Journal of mathematical physics , 169–187 (1999). T. Tarhasaari, L. Kettunen, and A. Bossavit, “Some realizations of a dis-crete hodge operator: a reinterpretation of finite element techniques [for em field analysis],” IEEE Transactions on magnetics , 1494–1497 (1999). K. Yee, “Numerical solution of initial boundary value problems involvingmaxwell’s equations in isotropic media,” IEEE Transactions on antennasand propagation , 302–307 (1966). M. L. Stowell and D. A. White, “Discretizing transient current densi-ties in the maxwell equations,” Tech. Rep. (Lawrence Livermore NationalLab.(LLNL), Livermore, CA (United States), 2009). D.-Y. Na, F. L. Teixeira, and W. C. Chew, “Polynomial finite-size shapefunctions for electromagnetic particle-in-cell algorithms based on unstruc-tured meshes,” IEEE Journal on Multiscale and Multiphysics Computa-tional Techniques , 317–328 (2019). S. O’Connor, Z. D. Crawford, J. Verboncoeur, J. Lugisland, andB. Shanker, “A set of benchmark tests for validation of 3d particle in cellmethods,” (2021). V. Kovalev and V. Y. Bychenkov, “Analytic solutions to the vlasov equa-tions for expanding plasmas,” Physical review letters , 185004 (2003). R. Cousin, J. Larour, J. Gardelle, B. Cassany, P. Modin, P. Gouard, andP. Raymond, “Gigawatt emission from a 2.4-ghz compact magnetically in-sulated line oscillator (milo),” IEEE Transactions on Plasma Science ,1467–1475 (2007). D. N. Arnold, “Spaces of finite element differential forms,” in
Analysis andnumerics of partial differential equations (Springer, 2013) pp. 117–140. T. Z. Esirkepov, “Exact charge conservation scheme for particle-in-cell sim-ulation with an arbitrary form-factor,” Computer Physics Communications , 144–153 (2001). C.-D. Munz, P. Omnes, R. Schneider, E. Sonnendrücker, and U. Voss, “Di-vergence correction techniques for maxwell solvers based on a hyperbolicmodel,” Journal of Computational Physics , 484–511 (2000). M. C. Pinto, M. Mounier, and E. Sonnendrücker, “Handling the divergenceconstraints in maxwell and vlasov–maxwell simulations,” Applied Mathe-matics and Computation , 403–419 (2016). K. F. Warnick and P. H. Russer, “Differential forms and electromagneticfield theory,” Progress In Electromagnetics Research , 83–112 (2014). P. J. Mardahl, J. P. Verboncoeur, and C. k. Birdsall, “Progress on a 3dparticle-in-cell model of a w-band klystron,” in
IEEE Conference Record- Abstracts. 1999 IEEE International Conference on Plasma Science. 26thIEEE International Conference (Cat. No.99CH36297) (1999) pp. 161–. G. Chen, L. Chacón, L. Yin, B. J. Albright, D. J. Stark, and R. F. Bird,“A semi-implicit, energy-and charge-conserving particle-in-cell algorithmfor the relativistic vlasov-maxwell equations,” Journal of ComputationalPhysics , 109228 (2020). D. Seidel, M. Kiefer, R. Coats, T. Pointon, J. Quintenz, and W. Johnson,“The 3d, electromagnetic, particle-in-cell code, quicksilver,” Tech. Rep.(Sandia National Labs., Albuquerque, NM (USA), 1990). T. Arber, K. Bennett, C. Brady, A. Lawrence-Douglas, M. Ramsay, N. Sir-combe, P. Gillies, R. Evans, H. Schmitz, A. Bell, et al. , “Contemporaryparticle-in-cell approach to laser-plasma modelling,” Plasma Physics andControlled Fusion , 113001 (2015). J. P. Verboncoeur, “Particle simulation of plasmas: review and advances,”Plasma Physics and Controlled Fusion , A231 (2005). A. J. Christlieb, R. Krasny, J. P. Verboncoeur, J. W. Emhoff, and I. D. Boyd,“Grid-free plasma simulation techniques,” IEEE Transactions on PlasmaScience34