Variational Framework for Structure-Preserving Electromagnetic Particle-in-Cell Methods

In this article we apply a discrete action principle for the Vlasov–Maxwell equations in a structure-preserving particle-field discretization framework. In this framework the finite-dimensional electromagnetic potentials and fields are represented in a discrete de Rham sequence involving general finite element spaces, and the particle-field coupling is represented by a set of projection operators that commute with the differential operators. With a minimal number of assumptions which allow for a variety of finite elements and shape functions for the particles, we show that the resulting variational scheme has a general discrete Poisson structure and thus leads to a semi-discrete Hamiltonian system. By introducing discrete interior products we derive a second type of space discretization which is momentum preserving, based on the same finite elements and shape functions. We illustrate our method by applying it to spline finite elements, and to a new spectral discretization where the particle-field coupling relies on discrete Fourier transforms.


Introduction
Since the early days of Particle-in-Cell (PIC) schemes, plasma physicists have devised variational algorithms based on least action principles to preserve key invariants such as the total energy and Gauss's laws [25,26,14].In parallel, a Hamiltonian structure of the Vlasov-Maxwell equations has been proposed, that involves a non-canonical Poisson bracket [30,34,28].Although the first methods were developed for finite difference field solvers, many improvements have been made and in the last decade several schemes have been proposed that rely on the de Rham structure of the Maxwell equations [6,20] to guarantee an exact preservation of proper discrete Gauss laws for general Finite Element PIC methods on general meshes [12], later extended to variational PIC schemes in e.g.[33] and [15,32], where it was shown that variational spectral methods also preserve the total momentum of the plasma.
Following these ideas a Geometric Electromagnetic PIC (GEMPIC) method based on spline finite elements has been proposed in [23], that possess a Hamiltonian structure relying on a discrete Poisson bracket.Coupled with Hamiltonian splitting methods [18,13,19], this approach leads to fully discrete schemes that preserve a modified energy, discrete Gauss laws, and the Poisson structure of the semi-discrete problem, including its associated Casimir invariants [23].
In this article, we extend these constructions to a flexible and general setting that allows for arbitrary structure-preserving discretizations of the electromagnetic fields and a variety of particle-field coupling operators which in particular includes almost arbitrary smoothing shape functions.By applying a discrete action principle we rigorously derive a variational system of discrete Vlasov-Maxwell equations, and we show that it has a non-canonical Poisson structure.This approach allows for instance to derive numerical Maxwell solvers with a strong Ampère and Gauss equation, and also extends the strong Faraday solver of [23] to more general particlefield coupling schemes.Another direct application is the design of variational spectral particle methods, where the Maxwell equations are solved in discrete Fourier spaces.
The outline of the paper is as follows.In Section 2 we first present the commuting de Rham complex that serves as the basis of our discrete derivation.This setting is now common in the structure-preserving (mimetic) discretization of Maxwell equations, and has been thoroughly studied in the Finite Element Exterior Calculus (FEEC) literature.For the Vlasov-Maxwell equations it describes how the particle-field coupling operators are connected with the differential operators involved in the discrete Maxwell equations.Then, we derive a variational particle discretization of the Vlasov-Maxwell system in a strong Ampère formulation from a discrete action principle, and analyze its main conservation properties together with its discrete Poisson structure.In Section 3, we present a variant of our method that preserves exactly the Gauss laws and the total momentum.In Section 4 a matrix form of the equations is carefully detailed, which also allows to derive a matrix form of the discrete Poisson bracket.In Section 5, we show how our analysis extends to a more general setting, and easily applies to the case of strong Faraday solvers.A detailed application to the case of structure-preserving Spline and Fourier discretizations is then presented in Section 6, with particle-field coupling operators based on geometric degrees of freedom which amount to discrete Fourier transforms in the spectral case.In Section 7, we conclude with preliminary numerical experiments that validate our approach and we compare the results obtained by various space discretizations that fit into our general framework, including different Maxwell solvers and different orders of particle smoothing.

Variational particle-field discretization 2.1 Maxwell equations and particle trajectories
A kinetic description of the dynamics of a plasma in an electromagnetic field (E, B) models the particles of species s by a distribution function f s in phase-space that evolves according to the Vlasov equation where m s and q s denote the mass and charge of the particle species s.The self-consistent fields evolve according to Maxwell's equations which are coupled to the Vlasov equation through the charge and current densities, We refer to e.g.[17,4] for a detailed presentation of these equations.As the Vlasov equation is a conservative transport equation, the distribution function f s is constant over time along the characteristic trajectories for that species, which are solution to the characteristic ODEs d dt In particle methods the distribution function is often represented by a collection of N macroparticles with phase-space positions (X p , V p )(t) and weights w p , of the form where S is a shape function that can either by the Dirac δ distribution or some smoothing kernel, depending on the particular configuration of the particle method.Starting from a collection of initial positions (X 0 p , V 0 p ), p = 1, . . ., N , the weights are initialized so as to provide a good approximation to the initial density f 0 s , and the particle positions are evolved according to some discrete characteristic equation, in order to approximate the trajectories (4).For the solution of Maxwell's equations, a grid-based solver is commonly used.

Structure of Maxwell's equations and finite element exterior calculus
As has been evidenced by several key contributions in the last decades [5,6,20], the Maxwell equations (2) possess a geometric structure where a central role is played by de Rham sequence In order to derive structure-preserving schemes we will follow the framework of Finite Element Exterior Calculus (FEEC) developed in e.g.[29,20,1,2,8,11].A central feature of these approaches is to involve a discretization that preserves the sequence (6) at the discrete level, and that admits a sequence of projection operators Π 0 , . . ., Π 3 mapping infinite-dimensional function spaces into discrete ones: In our framework, it is these operators Π ℓ , together with some shape (smoothing) functions S, that will encode the coupling mechanism between the particles and the discrete fields.Here the top row contains the infinite-dimensional domain spaces V ℓ of the operators Π ℓ , which are in general proper subsets of the natural Hilbert spaces involved in the sequence (6), and the bottom row consists of general discrete spaces such as finite-element or spectral spaces, see e.g.Section 6.
A key ingredient in our variational derivation will be that the operators Π ℓ make the diagram commuting.In practice many choices can be made for these operators and the associated finiteelement spaces where the fields are discretized.Each choice will result in a different coupling mechanism between the particles and the fields, but all of them will lead to Hamiltonian systems, provided the following property holds.Assumption 1.The operators Π ℓ : V ℓ → V ℓ h are such that: • the diagram (7) commutes, i.e., we have • the domain spaces V ℓ are translation invariant function (or distribution) spaces, in the Since the commuting projection operators will be applied to particle shape functions, we also need to specify when these shapes are admissible.
Definition 1 (admissible shape functions).A shape function S is said to be admissible for a given sequence of operators Π ℓ if it belongs to the domain spaces V 0 and V 3 of Π 0 and Π 3 , and if for any e ∈ Ê 3 , eS belongs to the domains V 1 and V 2 of Π 1 and Π 2 .
Remark 1.In practice, the translation invariance assumption corresponds to defining the projection operators on domain spaces V ℓ characterized by some homogeneous regularity over Ê 3 , which simplifies the notion of an admissible shape function S. In special cases where one works with localized or heterogeneous domain spaces, some additional care may need to be taken to guarantee that the projection operators can be applied on the shape functions.

Discretizing the Ampère or Faraday equations in strong form
In the article [23] the discretization ansatz was to consider fields in the spaces with φh and Ãh denoting discrete representations of the scalar and vector potentials, and this has led to an approximation of Ampère's and Faraday's laws in weak and strong form, respectively.Although the analysis presented here readily applies to the ansatz (11), it also covers the dual choice which leads to a new discrete model involving a strong Ampère law and a weak Faraday law.Throughout this article we will thus focus on this new ansatz (12), and describe in Section 5 how our results apply to the 'strong Faraday' ansatz (11).
In both cases, the discrete equations in weak form will involve the discrete adjoints to the strong differential operators, i.e., , and ψ h ∈ V 0 h .These discrete operators may be seen as the discrete Riesz representants of the differential operators in distribution's sense.

Discrete Action principle
We now derive a general geometric electromagnetic particle method where, following the ansatz (12), the Ampère equation is discretized in a strong sense.Here the coupling mechanism is essentially encoded in the abstract operators Π ℓ that are only assumed to satisfy the commuting diagram properties, see Assumption 1, and in the shape function S that must be admissible in the sense of Definition 1.
To do so we follow a discrete variational principle in the spirit of [30,33,22,15], based on Low's Lagrangian functional for the Vlasov-Maxwell equations [27], Here the curves X = X(t; x 0 , v 0 ), X ′ = X ′ (t; x 0 , v 0 ), V = V (t; x 0 , v 0 ) depend on time and on the initial conditions, and we recall that in a variational derivation they represent independent variables of the functional, in particular the prime symbol does not stand for a derivative.We also note that a different set of characteristics is associated to each particle species, which has been left implicit here for notational simplicity.Formally, the Vlasov-Maxwell equations can be derived as the Euler-Lagrange equations associated with this Lagrangian, as shown in [27].Here we will carefully apply this principle at the discrete level, starting from the discrete Lagrangian functional This Lagrangian is a function of discrete variables,

are arbitrary collections of trajectories and
h are arbitrary finite element potential fields.In (15) the dependence on t is implicit, and again we recall that the prime symbol does not mean a derivative, as all these functions are independent in the variational derivation.Finally the coupling potentials are defined as where S Xp (x) = S(x − X p ) denotes the shape function centered on a particle.We note that L h is formally derived from the continuous functional ( 14) by (i) replacing the initial density f s by its Dirac approximation in (5), i.e., f δ s,N (t 0 , x 0 , v 0 ) = N p=1 w p δ(x 0 − X 0 p )δ(v 0 − V 0 p ), (ii) using trajectories satisfying (X, V )(t; X 0 p , V 0 p ) = (X p , V p )(t), (iii) potential fields in the discrete (finite element) spaces, (iv) weak discrete differentials (13) instead of the exact ones, and finally (v) the coupling fields (16) defined with admissible shape (smoothing) functions in (5).In the case of several species, each density f s is approximated by a different set of discrete particles, so that we actually have s = s(p) in (15).For this reason it will be convenient to denote in the sequel particle masses and charges by m p := w p m s(p) and q p := w p q s(p) , for p = 1, . . ., N.
The discrete Action functional is then defined as and following a discrete action principle we look for generalized trajectories that form an extremum of S h .We already point out that the resulting equations will only involve the fields hence they will be gauge-independent.Formally, extremality conditions for S h are associated to the Euler-Lagrange equations of the discrete Lagrangian functional (15).Thus we look for X N , V N , φ h , and A h such that the following relations hold for all t ∈ [0, T ], with functional Gateaux derivatives evaluated at (X N , For the variations with respect to V N , we compute for an arbitrary VN = ( Vp ) p=1,...,N , so that (20) gives Using the coupling potentials (16), we compute for the variations with respect to for an arbitrary XN = ( Xp ) p=1,...,N .We then write Equation ( 21) for a variation of a single particle 1 ≤ p ≤ N along the unit basis vector e α ∈ Ê 3 for some dimension 1 ≤ α ≤ 3. Thus we take Xp ′ = δ p ′ ,p e α , which gives where we have used the commuting diagram property ( 8)-( 10) of the operators Π ℓ , and the definition (19) of the fields in the last equality.Using the linearity of the projection operator we rewrite the magnetic rotation term as with a coupling magnetic field defined at the particle position as Defining similarly the coupling electric field by we arrive at a velocity equation of the form Turning to the variations with respect to A h , using again (16) we compute so that Equation (23) gives The latter can be rewritten only in terms of the fields (19) and the particle current defined as and ( 17), as where we have used again the definition of the weak operators (13).Since both −∂ t E h + curl B h and Π 2 J S N belong to V 2 h , and (31) holds for all Ǎh ∈ V 2 h , it leads to an Ampère equation in strong form, In turn, a weak Faraday equation involving the discrete curl (13), follows from the definition of the fields (19): Indeed, for all Bh ∈ V 1 h we have which amounts to (33), by using the fact that Ω grad w φ h • curl Bh = Ω φ h div curl Bh = 0.For the variations with respect to φ h we use once more (16) and compute for an arbitrary φh ∈ V 3 h , so that (22) gives Using the field E h defined in (19) and noting that (34) must hold for all φh ∈ V 3 h , we arrive at a Gauss law in strong form, see again (5), (17).Finally a discrete magnetic Gauss law, this time in weak form, follows again from the definition (19) of B h = curl w A h , writing that

The variational equations
Gathering the findings of the variational derivation just detailed, we obtain a system of semidiscrete equations where the fields h are governed by the discrete Ampère and Faraday equations with a weak curl w : (13), and particles follow the trajectory equations with coupling fields defined by ( 27)-( 28), namely where (e 1 , e 2 , e 3 ) is an orthonormal basis of Ê 3 .These evolution equations are completed with two discrete Gauss laws, div with ρ S N = N p=1 q p S Xp and the weak divergence operator div w : V 1 h → V 0 h defined by (13).We note that here the first Gauss law has been derived from the variational principle (considering variations in the electric potential), whereas the second one follows from the definition of the magnetic field.

Derivation of a discrete Hamiltonian and an associated Poisson bracket
In this section we describe how the above variational equations can be associated with a discrete Poisson bracket.Following Hamilton's method [18, Sec.VI.1.2],we observe that our discrete Lagrangian has two nonzero conjugate momenta given by ( 25) and ( 30), which we may identify with their Riesz representant in the proper spaces.Assuming that the discrete solution satisfies the variational equations ( 37)-(40), we have Using the form of the coupling potential ( 16) and the variational Gauss law (40) we have so that the resulting Hamiltonian can be reformulated as a function of the fields (19), namely By construction this Hamiltonian is preserved by any solution satisfying the Euler-Lagrange equations ( 20)- (23).Following [3, Sec.40-A], a discrete Poisson bracket {F h , G h } can then be associated to the evolution equations ( 37)-( 39), such that holds for an arbitrary functional F h of the discrete solution.To identify this bracket we may simply consider linear functionals of the form defined by and G h = H h .Since the Poisson bracket should be a bilinear antisymmetric expression of the derivatives of its respective functionals, which read (upon identification with their proper discrete Riesz representant) and observing by linearity of F that (42) just amounts to the evolution equations ( 37)-(39) written in weak forms, with XN , VN , Ěh , Bh as test fields, we verify that (42) holds with the following discrete bracket where we remind that the coupling magnetic field B S (X p ) is defined in (39) and involves the projection operator Π 1 .We observe that this field plays the role of a parameter of the bracket, as do the shape functions centered on the particle positions, S Xp .A different role is played by the electric coupling terms, which enter the bracket through the product of V -E derivatives.
Below we will verify that this bracket is a (non-canonical) Poisson bracket in the sense of [18, Def.VII.2.4], in particular it satisfies the Jacobi identity.We note that other brackets involving different coupling fields B S would still be antisymmetric, and hence also energy-preserving.As the different projection operators are connected by the commuting diagram properties which have been used in several steps of the least action principle derivation, such brackets would probably not be variational, but they could maybe still satisfy the Jacobi identity.

Semi-discrete conservation properties of the variational system
One major property of the above derivation is that the resulting semi-discrete system has a Poisson structure, under the very general assumption that the diagram ( 7) is commuting.This result, whose proof will be given in Section 4.3, implies in particular that the evolution equations ( 37)-( 39) preserve all the functionals F h such that {F h , H h } = 0, which includes the Hamiltonian itself, F h = H h , but also all the Casimirs of the bracket (43) which are the functionals C h such that {C h , G h } = 0 for all G h , and new Casimirs may be derived using the Jacobi identity, see e.g.[18].An important example is provided by the functionals associated to an arbitrary φh ∈ V 3 h .The fact that they are Casimirs will be verified just below, and it implies that the discrete Gauss law div E h = Π 3 ρ S N is preserved by our equations.Theorem 1 will be most conveniently proven on a matrix form of the equations, which we will describe in Section 4. However a few basic conservation properties can be proven with a direct argument.Theorem 2. Under the conditions of Theorem 1, the evolution equations (37)-(39) preserve the discrete Hamiltonian (41) as well as the variational Gauss laws (40).
Remark 2 (weak Gauss law).Similarly as for the GEMPIC method [23], the magnetic Gauss law plays the role of a pseudo-Casimir, in the sense that its conservation is actually needed to establish that the evolution system has a discrete Hamiltonian structure.With a strong-Ampère ansatz (12), we observe that this divergence-free constraint is only preserved in a weak sense, see (36).Although this may seem very weak, we will see below that it is the natural discrete invariant that provides a Poisson structure for the resulting Hamiltonian system.
Proof.The preservation of the magnetic Gauss law readily follows from the weak Faraday equation in (37), indeed we have for all ϕ h ∈ V 0 h , using again the definition (13) of the weak curl operator.Turning to the electric Gauss law, we use d dt X p (t) = V p to compute for an arbitrary smooth function ψ which shows that the continuity equation always holds in distribution's sense, independently of the discrete particle trajectories.Taking next the divergence of the discrete Ampère equation in (37), the commuting diagram property (10) (which holds thanks to the admissibility of S) allows us to write where the last equality follows from (45) and from the time-invariance of the operator Π 3 .Integrating over time this shows that the electric Gauss law is indeed preserved.Another argument consists of verifying that any functional of the form (44) is indeed a Casimir.To do so we compute that the (Riesz representants of the) functional derivatives of C h read As for the derivatives δC h δVp and δC h δB h , they vanish.For the discrete bracket (43) we thus find Here the first term vanishes for arbitrary vectors δG h δVp ∈ Ê 3 , by using the commuting diagram property and the definition of the weak gradient operator.As for the second term, a discrete integration by parts yields Ω grad w φh N is an invariant of the evolution system.Finally to verify the energy conservation, we may simply observe that the bracket (43) is antisymmetric, so that F h = H h is an obvious invariant of (42).A more pedestrian argument is to first compute using (37) where we have used the adjoint definition of curl w , and then, using the trajectory equations (38)-(39), d dt which shows that the discrete energy (41) is indeed constant over time.

Generic Gauss and momentum preserving schemes
Similarly as for the method in [23], the semi-discrete scheme derived above is in general not momentum-preserving.However it is possible to describe a general variant that preserves both the Gauss laws and a discrete momentum.This modified scheme comes at the price of losing the discrete Hamiltonian (Poisson) structure and the conservation of energy, but it may be preferred for problems where momentum preservation is critical.

Particle-field coupling with discrete interior products
Our momentum-preserving schemes rely on discrete interior products of the form which involve the continuous interior products ı ℓ eα : V ℓ+1 → V ℓ associated with a canonical unit vector e α , α ∈ 1, 3 , namely and where A h,α is a linear approximation operator, such that the operators I ℓ eα map every discrete space to its predecessor in the sequence, as stated in (46).
As a key property, denoting by d 0 = grad, d 1 = curl and d 2 = div, we require that the associated discrete Lie derivatives, defined as Specifically, the momentum preserving properties will rely on the following relations and

Gauss and momentum-preserving schemes
Using the discrete interior products described above, we obtain the following result.
Theorem 3. The scheme obtained by coupling the discrete Maxwell equations (37)-(38) with the modified particle equations with coupling fields defined as preserves the discrete Gauss laws (40), as well as the discrete momentum Remark 3. Given the form (46)-(47) of I 1 eα and the linearity of A h,α , we have which makes clear how (52) approximates the exact momentum along e α .Similarly, we have which shows that the discrete magnetic force involved in (50) is indeed an approximation of the "natural" term V p × B h (X p ). However it is not possible in general to write R S (B h , X p , V p ) as a product of the form V p × B S (X p ) for some field B S , because the approximation operators A h,α involved in the trajectory equation depend a priori on the component α of the latter.
Proof.We first observe that the arguments used in the proof of Theorem 2 for the conservation of the discrete Gauss laws did not rely on the particle trajectory equation, hence they are still valid for the modified scheme.Turning to the discrete momentum, we compute using (50) Using next (37) we write where we have used the definition of the weak curl operator in the second equality, the relations (48)-(49) in the third one and the preservation of the discrete (weak and strong) Gauss laws in the last one.

Interior products based on directional averaging on tensor-product spaces
In this section we show that a simple construction based on directional averaging allows to design momentum-preserving schemes when the compatible sequence involves tensor-product spaces of the form where the univariate spaces Í α h , Î α h , form an exact sequence along each dimension α ∈ 1, 3 , Lemma 1. Assume that the univariate sequences (56) are exact, with spaces Í α h invariant over translations of ±h α , 1 ≤ α ≤ 3. Then the discrete interior products I ℓ eα = A h,α ı ℓ eα defined by composing the exact interior products (47) with the directional averaging operator, and similarly for α = 2, 3, map V ℓ+1 h to V ℓ h .Furthermore, they satisfy the relations (48)-(49).
Proof.Let us show that 47) and the tensor-product structure (54)-(56).The exact sequence property (56) then allows us to write Λ 1,α kα = ∂ α Γ 0,α kα for some Γ 0,α kα ∈ Í α h , which yields which belongs to V 0 h , according to (54) and the discrete translation invariance.The argument for the other spaces is similar.Turning to (48)-(49) we next observe that the directional averaging operators are of the form A h,α G = µ α * G with a symmetric measure µ α (−x) = µ α (x).Thus, for all α, β, and any function G.This allows to write a proof that is formally the same as for the continuous interior product (47).Thus, using that (curl which proves (48).The relation (49) follows by a similar argument.
4 The semi-discrete Hamiltonian system as a system of ordinary differential equations In this section, we express the variational particle method (37)-(38) as a system of ordinary differential equations.This will allow us to introduce some useful notation for our general framework, and to verify the Hamiltonian structure of the semi-discrete system.

Commuting diagrams with degrees of freedom
One practical approach to build commuting projection operators is to introduce one additional layer in the diagram (7), consisting of coefficient spaces C ℓ = Ê N ℓ corresponding to the choice of specific bases for the finite-dimensional spaces V ℓ h with dimension N ℓ .This approach is somehow parallel to the geometric construction of [24] where commuting de Rham complexes are described for differential forms.As we consider here a a finite element setting, we will follow similar principles but our construction does not involve differential forms.
In this diagram the main novel ingredient is the degrees of freedom σ ℓ = (σ ℓ i ) 1≤i≤N ℓ , which must be unisolvent for the finite-dimensional spaces V ℓ h in the usual sense that they must be one-to-one when restricted to these spaces.The spaces V ℓ then denote the domains of these degrees of freedom, and as above we consider a conforming discretization in the sense that V ℓ h ⊂ V ℓ .The other discrete entities can then be determined from the degrees of freedom.
• The "interpolation" operators I ℓ are characterized by the right-inverse property σ ℓ I ℓ g = g for all g ∈ C ℓ .In particular, the basis functions Λ ℓ i ∈ V ℓ h defined by the usual duality relations correspond to Λ ℓ i = I ℓ e ℓ i where e ℓ i = (δ i,j ) 1≤j≤N ℓ is a canonical basis vector of C ℓ .It is sometimes convenient to stack the basis functions into colum vectors Λ ℓ = (Λ ℓ i ) 1≤i≤N ℓ , and to use a matrix notation for stacked functionals evaluated on vectors of functions.With this convention, the duality relation (59) reads • The matrices D ℓ ∈ Ê N ℓ+1 ×N ℓ correspond to the differential operators d 0 = grad, d 1 = curl and d 2 = div in the respective bases, namely • The projection operators are defined as and they are characterized by the relations indeed we have This setting proves particularly useful in practice, as it allows to restate the commuting diagram properties (7) as a linear relation between degrees of freedom.
Proof.The proof is a matter of elementary computations.For instance, (65) yields where we have used twice the characterization (63).

The semi-discrete Hamiltonian system in matrix form
The introduction of a third layer in the commuting diagram offers the possibility to rewrite the semi-discrete scheme (37)-(38) as a system of ordinary differential equations in matrix form.
To do so we collect all the dynamic variables in a global vector where the (column) block-vectors ) N collect all the particle positions and velocities as in Section 2.4, while the vectors the coefficients of the electric and magnetic fields in their respective bases.Using these degrees of freedom, we observe that the coupling fields (39) read where M ℓ is the standard finite-element mass matrix in the corresponding basis of V ℓ h , ℓ = 1, 2, The value of the coupling fields at the particle positions may then be expressed as block-vectors, where S ℓ (X) ∈ (Ê 3 ) N ×N ℓ denotes the matrix with (3 × 1) blocks We finally let r(b) = (e α × e β ) • b 1≤α,β≤3 ∈ Ê 3×3 be the rotation matrix and we denote by R(b(X)) ∈ (Ê 3×3 ) N ×N the block-diagonal rotation matrix with blocks Then the particle trajectory equations (38)-(39), can be written in the block-matrix form where W q m = diag( qp mp : p ∈ 1, N ) is the diagonal weighting matrix carrying the particles charge to mass ratios, and where we have denoted the block-diagonal rotation matrix associated with the coupling magnetic field.Observe that its diagonal blocks read R 1 (X, B) p,p = σ 1 (e α × e β S Xp ) ⊤ M 1 B 1≤α,β≤3 .
Turning to the field equations (37), we see that the strong Ampère equation can be expressed directly on the degrees of freedom σ 2 .From the characterization of the projection operator (63) we have , S 2 (X) the matrix defined in (68) and W q the diagonal weighting matrix carrying the particles charges.Finally the weak Faraday equation is tested against the basis functions Λ 1 i .By definition of the weak curl operator (13) this yields with M 1 and M 2 the mass matrices recalled in (66).Finally, rewriting the discrete Hamiltonian H h (X N , V N , E h , B h ) as a function of the array variables with W m the diagonal weighting matrix carrying the particle masses, see (41), we obtain for the corresponding derivatives which allows us to rewrite the equations (37)-(39) as a system of ODEs with a structure matrix given by In particular, System (77) may be rewritten in the form of a Poisson system dU dt = {U, H}(U) with a discrete bracket defined as {F, G}(U) Note that this is just the matrix form of the discrete bracket {F h , G h } given in (43), where F h and G h are the same functionals as F and G but seen as functions of the finite element fields E h , B h .Compared with the Poisson matrix found in [23, Eq. (4.29)], we observe that the main difference lies in the fact that the particle-field coupling blocks now involve the degrees of freedom of the smoothed particles through the matrices S 2 and S 1 which involve the generic commuting diagram operators Π 2 and Π 1 , see ( 68) and (73).In particular, the similarity of both matrices allows us to easily verify the Poisson structure of the semi-discrete system (37)-(39).

Proof of Theorem 1
Since we have rewritten our equations in a matrix form, it suffices to show that J is a Poisson matrix in the sense of [18, Def.VII.2.4], i.e., that it is skew-symmetric and it satisfies the matrix Jacobi identity.This will show that (80) is a (non-canonical) Poisson bracket and that (79), namely (77), is a Poisson system.
Using that weighting matrices like W q m are diagonal, and that R 1 (X, B) is skew-symmetric, we easily verify that J = −J T .To verify the matrix Jacobi identity, we then observe that J has the same form as the one involved in the original GEMPIC scheme, see [23,Eq. (4.29)], with CM −1 1 and Λ 1 (X)M −1 1 replaced by (M 1 ) −1 C ⊤ and S 2 (X), respectively (the mass and curl matrices being defined for different spaces, due to the different ansatz in the fields).We also note that R 1 plays the role of the magnetic rotation matrix B in [23], with smoothed coupling terms as already observed.In particular, we may follow the same reasoning to verify that it satisfies the Jacobi identity, which amounts to verifying that the analog of Eqs.(4.34) and (4.38) hold in our case.Using the block-diagonal matrix R 1 (X, B) defined by (73), and taking (p, α), (p, β), and (p, γ) as multi-indices corresponding to b, c, and d, Equation (4.34) becomes (for q p = 0) Using the expression R 1 (X, B) (p,α),(p,β) = σ 1 (e α × e β S Xp ) ⊤ M 1 B seen above, this amounts to By antisymmetry, we see that the function in parentheses vanishes if two of the components coincide, so that we may assume w.l.o.g. that (α, β, γ) = (1, 2, 3).Then this function is just ∇S Xp and the above equation amounts to where we have used the commuting diagram property and the admissibility of the shape function S. The desired equality (81) then follows from the discrete magnetic Gauss law, see (40).
The second equality to verify is the analog of Equation (4.37) from [23], which reads here (given the above matrix correspondence and correcting a typo on the sign of the right-hand side) By antisymmetry of R 1 , we see that both sides vanish for α = β, so let us assume w.l.o.g. that (α, β) = (1, 2).Then R 1 (X, B) (p,α),(p,β) = σ 1 (e 3 S Xp ) ⊤ M 1 B and by differentiating these entries and those of the matrix S 2 , see (68), the equality becomes In vector terms this writes σ 2 curl(e 3 S Xp ) = Cσ 1 (e 3 S Xp ), which directly follows from the commuting diagram property as seen in Lemma 2. Thus (82) holds, which shows that J satisfies the Jacobi identity and is indeed a Poisson matrix.

Propagation in time
Based on its Poisson structure, geometric time propagation schemes can be derived for our variational system in the same way in [23].More precisely, a variational integrator can be derived from a Hamiltonian splitting, that yields a scheme that is explicit in time.We refer to [23,Sec. 5.1] where the resulting equations are detailed for the weak Ampère case and delta shape functions.This kind of splitting has originally been proposed for the Vlasov-Maxwell system in [19,35] as a Hamiltonian splitting and later been constructed from a fully discrete action principle in [36].On the other hand, energy-conserving time propagators can be derived by an antisymmetric splitting of the Poisson matrix combined with a suitable discrete-gradient time propagation of the substeps as explained in [21].In our numerical experiments, we consider the energy-and Gauss-conserving discrete-gradient method from [21], which demonstrates the best the conservation properties of the phase-space discretization, and the Hamiltonian splitting from [23] due to its simplicity.

Generalization and application to the strong Faraday model
Before turning to the description of particular discretizations of Maxwell's equations, it may be useful to pause for a moment and make some comments on the above findings.In our variational derivation we have explicitly required that (7) was a commuting diagram, and by doing so we have made two implicitly assumptions: first, we have considered that the discrete sequence involved strong differential operators, which corresponds to a conforming discretization.Second, we have referred to the operators Π ℓ as projection operators.Although these are standard properties to assume, they played no particular role in our analysis, be it in the variational derivation of Section 2, or in the proof of its Hamiltonian structure.In particular, our results directly apply to a more general setting of the form where the discrete differential operators grad, curl, div no longer need to coincide with the exact ones (in particular, the discrete spaces V ℓ h need not be conforming in H 1 , H(curl) and H(div)), and the Π ℓ no longer need to be projection operators.In this generalized setting the only assumptions are that: (i) the solid diagram in (83) commutes, (ii) the lower discrete differential are adjoint to the upper ones in the sense of ( 13), namely Ω φh div * Fh = − Ω ( grad φh ) • Fh must hold for all φh ∈ V 0 h and Fh ∈ V 1 h , and so on.
Our variational derivation then applies verbatim, starting from the discrete Lagrangian h and coupling potentials defined as in (16).The resulting variational equations, analog to (37)-(40), read with coupling fields defined similarly as in (39).Our analysis then shows that these general equations preserve both the corresponding discrete Gauss laws and the Hamiltonian, and that they have a discrete Poisson structure.This allows to extend our results to a wider range of discrete settings, including the structure-preserving DG-type Conga discretizations developed in [10,11] where both E and B are represented in broken finite element spaces.Our results also apply to the discrete ansatz (11) corresponding to a strong Faraday equation.For this case we may consider a conforming (strong) discretization of the form (7), and set so that the ansatz (11) takes a form similar to the one (12) considered above, namely A commuting diagram (83) involving the spaces (85) can then be obtained as follows: define the commuting (upper) discrete differential operators as the weak operators (13), i.e.
and for the projection operators Π ℓ simply take the L 2 projections on the discrete spaces, The commutation property is indeed easily verified: For the grad operator, using the embedding div : h and the characterization of L 2 projections, we can write , which shows that Π 1 grad = grad Π 0 holds on V 0 (which may be taken here as H 1 (Ω)).The same argument also applies for the operators curl and div.With this construction one recovers the Hamiltonian particle method of [23], with general shape functions.The discrete Poisson matrix thus takes the same form, with particle-field coupling terms encoded in block matrices Λ ℓ (X) which extend the corresponding matrices in [23] to the case of a general shape function S.
6 Application to tensor-product spline and Fourier field solvers In this section, we apply the above method to the case of tensor-product finite element spaces defined on cartesian domains.Following the interpolation / histopolation approach of [16,24], we review a general method for designing commuting diagrams, which is based on geometric degrees of freedom that can then be associated to finite element spaces of various types.In this article, we detail two applications, one using splines and another one using truncated Fourier spaces.

Geometric degrees of freedom with commuting properties
Let us equip the cartesian domain Ω = [0, L] 3 with a tensor-product grid using M α nodes along each dimension α, On this mesh, we consider evaluation functionals defined on the various geometric elements: • point evaluations on the nodes • edge integrals along some dimension 1 ≤ α ≤ 3, • face integrals normal to some dimension 1 ≤ α ≤ 3, • and cell integrals where we have denoted by [a, b] the convex hull of a ∪ b.A set of "geometric" degrees of freedom can then be derived from these local functionals: If these degrees of freedom are associated to spaces and for which they are unisolvent, then they define a unique set of dual basis functions Λℓ i according to (59), which may also be called "geometric": for the space V 0 h for example these basis functions correspond to the interpolatory basis associated with the nodes x m , for the space V 3 h they correspond to histopolation basis functions, and for the intermediate spaces they involve a combination of both.A key property of this construction is the following.Lemma 3. The degrees of freedom defined by (93) are well-defined on the domains per , where we have denoted by L 1 per the space of L-periodic and locally L 1 functions, and by anisotropic Sobolev spaces of W s,1 type.Moreover if the σℓ are unisolvent on the spaces V ℓ h , then the resulting projection operators Πℓ characterized by the relations (63), namely satisfy the commuting diagram property Proof.The fact that these degrees of freedom are well-defined on the above domains follows from standard Sobolev inequalities, see e.g.[7,Rem. 13].The commuting diagram properties are then easy to verify by applying the Stokes formula and Lemma 2. For the gradient for instance, we consider some ϕ ∈ V 0 and compute According to Lemma 2, this specifies the gradient matrix D0 ∈ Ê N 1 ×N 0 such that σ1 (grad ϕ) = D0 σ0 (ϕ) and also implies grad Π0 ϕ = Π1 grad ϕ.The same argument works for the other operators.
In the construction above, we see that the commuting properties rely only on the geometric nature of the degrees of freedom, and not on the tensor-product structure of the grid.However, this tensor-product structure allows us to specify the form of the differential matrices.Setting ϕ = Λ0 k in the proof of Lemma 3, we find indeed the following representation of D0 where I Mα is the identity matrix of size M α × M α , d α is a univariate differential matrix and the Kronecker matrix product is defined as (c ⊗ b ⊗ a) m,n = c m 3 ,n 3 b m 2 ,n 2 a m 1 ,n 1 .In the same way, we find where O M denotes the zero square matrix of size In practice, the basis functions Λℓ i defined by the geometric degrees of freedom according to (59) may not be the most convenient to use, either because they have no simple expression, or because some other basis Λ ℓ i has better locality properties, or leads to simpler discrete Maxwell equations.One then needs to determine the coefficients of the geometric projections in this new practical basis, which amounts to finding degrees of freedom σ ℓ i that are dual to the practical basis functions and lead to the same projection operator Π ℓ = Πℓ as the geometric ones.Using the stacked vector notation introduced in Section 4.1 for the geometric basis Λℓ and the practical basis Λ ℓ , these new degrees of freedom σ ℓ are characterized by the relations Introducing the matrix which gives a practical formula for computing the coefficients of the geometric projections in the practical basis.Accordingly, the differential matrices in this new basis read Note that K 0 is a Vandermonde matrix when Λ 0 is a monomial basis.For this reason the matrices K ℓ are sometimes referred to as a generalized Vandermonde matrices.

Compatible finite elements based on B-splines
Compatible finite elements based on splines on a Cartesian grid have been studied by Buffa, Sangalli, Vázquez and co-authors, see e.g.[8,9], and in [23] they have been used to implement the strong Faraday GEMPIC formulation.Here we describe how spline spaces can be used in conjunction with the geometric degrees of freedom described in Section 6.1.
For simplicity, we consider periodic boundaries and regular knot sequences with M α knots per dimension.Denoting by N p α,k the univariate B-spline of degree p along x α , associated with the knots (kh α , . . ., (k + p + 1)h α ) where h α = L Mα , see e.g.[31], the first space in the sequence consists of tensor-product splines of multi-variate degree (p 1 , p 2 , p 3 ), namely and the full sequence reads The fact that this is indeed a sequence follows from the well-known relation Introducing for convenience the scaled B-splines along x α , yields a particularly simple formula for the derivative operator in the corresponding basis.In particular, it makes it convenient to equip the vector-valued spaces V 1 h , V 2 h with the basis functions and the last, scalar-valued space V 3 h with In practice, B-splines are appealing because of their minimal support property, however they are not dual to the geometric degrees of freedom defined in (93) so that new degrees of freedom must be computed as described at the end of Section 6.1.For the nodal degrees of freedom, the change of basis matrix reads and a common choice of interpolation nodes x m consists of Greville points, which coincide with the knot sequence for regular splines of odd degrees, and with their midpoints for even degrees.More generally, we observe that K 0 is invertible as long as the degrees of freedom σ 0 are unisolvent, which holds iff the grid satisfies the spline interpolation condition, see e.g.[31,Th. 4.61].Using the tensor-product structure and the locality of the B-splines, we see that K 0 is the Kronecker product of three banded matrices, which are also circulant for regular Greville points.Moreover, as B-splines satisfy by construction see ( 99), (100), we have hence the matrix block K 1,α = K 1 (α,m),(α,k) m,k coincides with K 0 , the other blocks of K 1 being clearly zero.Similarly we find that K 3 and K 2,α also coincide with K 0 , so that (with obvious notation) From relation (99) we also see that the one-dimensional derivative matrices -and hence, every D ℓ -are the same as for the geometric basis.As for the three-dimensional mass matrices, they are the Kronecker product of the one-dimensional mass matrices which are a circulant matrices with 2p α + 1 non-zero entries per row in each dimension.Finally we note that the discrete interior products (46) based on the directional averaging operator (57) may be evaluated using the relation (102), writing e.g.

Compatible finite elements based on Fourier spaces
With periodic boundary conditions, another option is to consider a sequence of compatible finite elements made of discrete Fourier spaces.Such spectral elements are very common in particle solvers, with particle-field interaction usually based on discrete Fourier transforms and FFT algorithms.Here we describe a coupling based on the geometric degrees of freedom described in Section 6.1.To match the dimensions of the grid, we consider spaces with M α = 2K α + 1 modes per dimension, of the form , where we have denoted −K, K = 3 α=1 −K α , K α , and These discrete spaces clearly form a de Rham sequence, as the derivative of a Fourier mode is the same mode up to a complex scaling factor.One interesting feature of the canonical modal basis is that it leads to diagonal Maxwell equations.Indeed the differential matrices D ℓ have the same simple block and Kroneckerproduct structure as (96)-( 98), here with diagonal one-dimensional derivative matrices and the mass matrices are all diagonal due to the orthogonality of the basis functions, with M ℓ = L 3 I N ℓ for the chosen normalization.However, as the modal basis is not dual to the geometric degrees of freedom from Section 6.1, we need to determine the proper change of basis formulas in order to apply the geometric interpolation-histopolation projections Πℓ , as we did for the B-splines in the previous section.To do so, it is convenient to consider regular interpolation nodes x m = (m 1 h 1 , m 2 h 2 , m 3 h 3 ), with h α = L Mα .The nodal change of basis matrix reads then which is a standard DFT matrix as well as its inverse, where we remind that M = M 1 M 2 M 3 , see (94).The interpolation operator in the modal basis then takes the well-known form
For the other projections in the sequence we proceed similarly as in Section 6.2, noting that mαhα In particular, writing T α := diag T α k,k = T α kα : k ∈ −K, K we find K 1,α = K 0 T α for the matrix block K 1,α = K 1 (α,m),(α,k) m,k , and similarly K 2,α = K 0 T α−1 T α+1 and K 3 = K 0 T 1 T 2 T 3 .The expression of the different degrees of freedom in the modal Fourier basis reads then where we note that all the T α matrices are clearly diagonal and invertible.To apply the discrete interior products (46) based on directional averaging (57), we finally need to evaluate for all α ∈ 1, 3 and k ∈ −K, K .

Numerical illustration in reduced phase space
In this section, we will show some numerical results obtained with the proposed schemes in a reduced phase space.All results are obtained with an implementation of the strong Ampère scheme within the SeLaLib library.We study the variational semi-discretization as derived in Section 2-which is energy conserving-as well as the momentum-preserving semi-discretization as derived in Section 3.For the basis of the finite element field solver, both splines and Fourier modes are considered.The shape function is chosen to be a B-spline of varying degree.
As for the time discretization, we compare a Hamiltonian splitting scheme for both space discretization methods, see Section 4.4.Only when considering the conservation properties we also provide results for the variational scheme with an energy-conserving discrete gradient time discretization.We use a time step of ∆t = 0.05, the linear solvers use a tolerance of 10 −15 and the nonlinear iterations in the discrete gradient method have a tolerance of 10 −12 .

Physical model
For the numerical study we consider a reduced phase space with one periodic spatial and one or two velocity dimensions, namely Moreover, we simulate an electron distribution in a neutralizing ion background, which differs from the multi-species Vlasov-Maxwell system in that the average current is substracted from the total one in order for the model to be momentum preserving.In particular, the reduced Maxwell system then reads In some cases this model will be further reduced to 1d1v phase space by skipping v 2 , E 2 and B 3 , so that the equation for E 1 above remains as the only field equation.As a first test case, we consider the Weibel instability in 1d2v phase-space as studied in [23] with an initial value of As a second test case, we consider the two-stream instability in 1d1v phase-space with initial value with parameters ǫ = 0.001 and k = 0.2.The initial field E 1 is again determined from Gauss' law.For this test case, the reference solution is also produced with a Fourier solver and a piecewise affine spline as shape function, but the grid resolution is reduced to 31 cells (and 15 modes) while the particle number is increased to 5 • 10 6 .Note that this test case requires a lot more particles to produce qualitative results compared to the Weibel test case.In Sections 7.2 to 7.4 below we study the influence of different numerical parameters using the relevant energy curves for these two test cases, namely the magnetic and electric energy, plotted in Figures 1 and 2 respectively.In Section 7.5 we finally compare the long-time conservation properties of the schemes, looking at different error curves shown in Figure 3.

Influence of the shape function
We first study the influence of the shape function.Here, we expect two counteracting effects: On the one hand, a higher degree of the shape function yields smoother data for the field solver which can yield better results.On the other hand, higher order smoothing kernels smear out the influence of particles which yields a damping.This latter effect is clearly seen in the simulations with M = 7 cells (grid points) of Figures 1a and 2a.For this coarse resolution, low order splines give rather good results whereas higher order shapes lead to a visible damping in the instability growth rate for both test cases.Increasing the number of cells to M = 15 while keeping the number of particles constant as in Figures 1b and 2b, we observe both effects: In this case, the degree one spline yields too noisy data for the field solver, while a degree of e.g.seven yields too high damping and an intermediate degree of four yields rather accurate results.Our results also show that when increasing also the number of particles, the choice of the shape function is of lesser importance (cf.Figures 1d and 2d).

Influence of the space semi-discretization
In Figures 1c-1d and 2c-2d we next compare the variational scheme presented in Section 2 with the momentum-preserving variant from Section 3.Here we use the spectral finite element solver and a Hamiltonian splitting time discretization.With this configuration, the momentumpreserving scheme yields clearly worse results for the coarse resolution runs (in Figures 1c and  2c), as the instability growth rate is damped similarly as with higher order shape functions.With increased resolution (i.e., using twice as many cells and four times as many particles for both test cases), we find that both schemes yield rather good results for various orders of the shape function (in Figures 1d and 2d).Finally, we see in Figure 1f that the long-time accuracy of the variational semi-discretization can be significantly better than that of the momentumpreserving one: here the Weibel instability is run with a small number of particles and we find a qualitatively wrong behavior for the momentum-preserving scheme using a piecewise affine shape function, where other schemes perform correctly.In Figure 2f a similar comparison is done with the two-stream instability, using a higher particle resolution as required for this test case to produce qualitatively correct results.The long-time behavior is then found to be qualitatively good for the different schemes and shapes.

Influence of the finite element solver
In Figures 1e and 2e we then compare the different field solvers, namely the spectral solver and finite element solvers based on splines of degree one to three.Using a piecewise affine spline for the shape function and low resolution runs we find that the accuracy of the low order fem solver is of bad quality and it improves for higher orders and for the spectral solver.This observation holds for the two test cases.

Conservation properties
We now compare the conservation properties of the various methods.For this we consider long times simulations with both the variational and the momentum-preserving discretizations.For the time stepping, we consider in both cases a Hamiltonian splitting as before but we also provide the solution with an energy-conserving discrete gradient propagator for the variational scheme to show that the semi-discretization is indeed energy-conserving.Figures 3a and 3b show the relative error in energy conservation for the various runs.We can see that the energy is conserved up to the tolerance of the linear solvers for the variational scheme with an energyconserving discrete gradient time propagator.If we use the Hamiltonian splitting instead, there is an energy error but its behavior is oscillatory and decreases with decreased time step.This is the typical behavior for such Poisson integrators.Finally, we see that the energy error is larger for the momentum-preserving scheme, in particular for the low order shape function with a low particle resolution.For the variational scheme, on the other hand, the energy error does not depend on the shape function.
Figure 3c and 3d show the error in momentum for the various methods.We can see that the momentum-preserving scheme indeed preserves momentum up to machine precision.On the other hand, for the variational scheme the error in momentum increases as soon as the nonlinear phase of the simulations starts and later flattens out at a certain level.As this error level seems to be rather independent of the propagator, and is smaller for higher order shape functions, we conjecture that it is dominated by the error in the spatial semi-discretization.
Finally the error in Gauss' law as a function of time is shown in Figures 3e and 3f for the two test cases, respectively.We can see that all scheme preserve Gauss' law to machine precision.

Theorem 1 .
If the operators Π ℓ satisfy Assumption 1 and if the shape function S is admissible in the sense of Definition 1, then the discrete bracket (43) is a (non-canonical) Poisson bracket and the semi-discrete equations (37)-(39) are a Poisson system in the sense of [18, Def.VII.2.4].