On Geometric Fourier Particle In Cell Methods
Martin Campos Pinto, Jakob Ameres, Katharina Kormann, Eric Sonnendrücker
OOn Geometric Fourier Particle In Cell Methods
Martin Campos Pinto , Jakob Ameres , Katharina Kormann , and EricSonnendr¨ucker Max-Planck-Institut f¨ur Plasmaphysik, Garching, Germany Technische Universit¨at M¨unchen, Zentrum Mathematik, Garching, Germany
February 4, 2021
Abstract
In this article we describe a unifying framework for variational electromagnetic particleschemes of spectral type, and we propose a novel spectral Particle-In-Cell (PIC) schemethat preserves a discrete Hamiltonian structure. Our work is based on a new abstractvariational derivation of particle schemes which builds on a de Rham complex where Low’sLagrangian is discretized using a particle approximation of the distribution function. In thisframework, which extends the recent Finite Element based Geometric Electromagnetic PIC(GEMPIC) method to a wide variety of field solvers, the discretization of the electromagneticpotentials and fields is represented by a de Rham sequence of compatible spaces, and theparticle-field coupling procedure is described by approximation operators that commutewith the differential operators involved in the sequence. In particular, for spectral Maxwellsolvers the choice of truncated L projections using continuous Fourier transform coefficientsfor the commuting approximation operators yields the gridless Particle-in-Fourier method,whereas spectral Particle-in-Cell methods are obtained by using discrete Fourier transformcoefficients computed from a grid. By introducing a new sequence of spectral pseudo-differential approximation operators, we then obtain a novel variational spectral PIC methodwith discrete Hamiltonian structure that we call Fourier-GEMPIC. Fully discrete schemesare then derived using a Hamiltonian splitting procedure, leading to explicit time steps thatpreserve the Gauss laws and the discrete Poisson bracket associated with the Hamiltonianstructure. These explicit steps are found to share many similarities with a standard spectralPIC method that appears as a Gauss and momentum-preserving variant of the variationalmethod. As arbitrary filters are allowed in our framework, we also discuss aliasing errorsand study a natural back-filtering procedure to mitigate the damping caused by anti-aliasingsmoothing particle shapes. Representing the electromagnetic fields in truncated Fourier spaces has been a standard practicein plasma simulation, from the early particle schemes [26, 21, 3] to more recent parallel high-performance codes [37, 11, 16, 30].In these spectral particle methods, two main approaches exist. The simplest consists ofcoupling the particles to the fields via continuous Fourier transforms, which leads to a gridlessmethod sometimes called
Particle-in-Fourier (PIF) [30, 1, 29]. This method, which may includesmoothing techniques through low pass filters or smoothing shape functions, is naturally chargeand energy preserving. More importantly, it can be derived from a variational principle [13,33, 38], allowing for numerical schemes with very good stability on long time ranges. Using1 a r X i v : . [ phy s i c s . c o m p - ph ] F e b he translational invariance of the discrete Fourier spaces, it has also been shown to preservemomentum.Another approach that is often preferred for simulations where many modes are needed,such as turbulence in tokamak plasmas, is to use a grid for the coupling, and discrete Fouriertransforms (DFT). This leads to spectral Particle-in-Cell (PIC) methods [24, 25, 11, 16], whichwe may call Fourier-PIC here. Smooth (continuous) particle shapes are then necessary for thecoupling to be well-defined, which in Fourier space corresponds to a low-pass filtering. Usinga DFT grid allows to localize the coupling, but it also causes dispersion and aliasing errors[24, 25, 3], which in turn have been shown to lead to numerical instabilities of various types,including grid heating [4], electrostatic finite grid instabilities [22] and numerical Cherenkovinstabilities [15, 40, 16]. In most cases these methods are momentum preserving and in somecases they are also charge preserving, see e.g. [16]. However, it seems that until very recently aproper derivation of variational spectral PIC methods was still missing.This has been done in part in a new article [9] where a flexible yet rigorous method isproposed to design variational, gauge-free electromagnetic particle schemes with a Hamiltonianstructure relying on a non-canonical discrete Poisson bracket. Formulated in a general frame-work with a minimal set of properties, this approach essentially extends the recent GeometricElectromagnetic PIC (GEMPIC) scheme [23], based on spline finite elements and point particles,to a wide range of field solvers and particle coupling techniques. In particular, one applicationdescribed in [9] consists of a new Hamiltonian spectral PIC method, where the particle-fieldcoupling is done through a DFT grid. However, as this method relies on particle-field operatorsdefined through geometric degrees of freedom, the resulting deposition algorithms involve vol-ume integrals of the shaped particles as they travel through space, which somehow deviates fromstandard PIC scheme that are based on pointwise evaluations of the particle shape functions.In this article we thus propose a novel variational spectral PIC method called Fourier-Gempic, that has a discrete Hamiltonian structure and relies on particle-field coupling tech-niques very similar to those encountered in standard PIC schemes. Our method is obtainedby applying the abstract derivation of [9] to truncated Fourier spaces, combined with novelparticle-field coupling operators of pseudo-differential type. By observing that our variationalderivation combined with continuous ( L -orthogonal) projection operators leads to the grid-less PIF method, we show that this framework actually unifies the formulation of variationalspectral particle methods, where the standard PIC scheme appears as a momentum-conservingvariant, albeit non-variational.By applying a Hamiltonian splitting technique in the spirit of [10, 19] we are able to proposefully discrete schemes for the three different methods, with explicit time steps that provide anadditional insight into their differences and similarities. For the Hamiltonian GEMPIF andFourier-GEMPIC methods, these fully discrete schemes preserve the total energy within thetime splitting error as guaranteed by backward error analysis [17]. For the PIF and standardPIC methods, they preserve the total momentum. All of them preserve the Gauss laws tomachine accuracy.In addition to a Hamiltonian structure which guarantees very good stability properties onlong time ranges, the Fourier-Gempic method has the ability to include various shape functionsand filter coefficients in Fourier space. In particular, we show that a natural back-filteringmechanism can be associated to the usual low-pass filtering effect of high order spline shapes inorder to strongly reduce aliasing errors inherent in the DFT, without damping relevant modesin the computational range.The outline is as follows. In Section 2 we recall the two main coupling approaches for Fourier-2article methods, namely Particle-in-Fourier (PIF) and Fourier-PIC, we discuss aliasing errorsand consider a simple back-filtering technique to mitigate the smoothing effect of anti-aliasingsplines. In Section 3 we then present the general form of a variational spectral scheme asderived in [9], as well as that of the momentum-conserving variant, and apply it to two par-ticular sets of particle-field approximation operators: using L projections which correspondto continuous Fourier coefficients, this leads to the gridless PIF method, which coincides withits momentum-conserving variant. Using a novel class of pseudo-differential DFT operators,we obtain a new Fourier-PIC method, whose Hamiltonian structure is guaranteed by the com-muting de Rham diagram property of the new DFT approximation operators and the analysisfrom [9]. In Section 4 we then describe fully discrete schemes of arbitrary orders, obtained byapplying a Hamiltonian splitting procedure. Both in the gridless and PIC cases, we provideexplicit formulas for the different steps of the discrete schemes. The standard Fourier-PICcoupling is then found to coincide with the momentum-conserving variant of the HamiltonianFourier-GEMPIC method, while the latter differs in the computation of the pushing fields. InSection 5 we assess the basic numerical properties of the different methods. Their accuracy iscompared on standard test cases, as well as their long-time stability and conservation proper-ties. It is shown that all Hamiltonian methods are very stable in energy and momentum, andthat back-filtered methods are very accurate for approximating the fundamental growth anddamping rates, even for low-resolution runs. Electromagnetic particle models formally consist of Maxwell’s equations for the fields, ∂ t E ( x , t ) − curl B ( x , t ) = − J N ( x , t ) (1) ∂ t B ( x , t ) + curl E ( x , t ) = 0 (2)div E ( x , t ) = ρ N ( x , t ) (3)div B ( x , t ) = 0 , (4)coupled with discrete particles with positions X p and velocities V p , p = 1 , . . . , N , subject toLorentz force trajectoriesdd t X p ( t ) = V p ( t ) , dd t V p ( t ) = q p m p (cid:0) E ( X p ( t ) , t ) + V p ( t ) × B ( X p ( t ) , t ) (cid:1) . (5)Here q p and m p represent the charge and mass of the p -th numerical particle, and the currentand charge densities, which are the sources in Maxwell’s equations, are obtained by summingeach particle contribution, ρ N ( t, x ) = (cid:88) p =1 ··· N q p δ ( x − X p ( t )) and J N ( t, x ) = (cid:88) p =1 ··· N q p V p ( t ) δ ( x − X p ( t )) . (6)In the limit N → ∞ of infinitely many particles, this model approximates the kinetic Vlasovequation [36] where the plasma is represented by a continuous phase-space density function f ( t, x , v ), and the choice of Dirac densities in (6) corresponds to a pointwise evaluation of thecontinuous charge and current densities ρ ( t, x ) = (cid:82) qf ( t, x , v ) d v , J ( t, x ) = (cid:82) q v f ( t, x , v ) d v .In spectral solvers, electromagnetic fields are represented by truncated Fourier expansions3f the form E K ( t, x ) = (cid:88) k ∈ (cid:74) − K,K (cid:75) E k ( t )e π k · x L B K ( t, x ) = (cid:88) k ∈ (cid:74) − K,K (cid:75) B k ( t )e π k · x L for x ∈ [0 , L ] , (7)here with 2 K + 1 modes per dimension, and the sources need to be properly represented in thesame truncated Fourier spaces. The discrete Maxwell’s equations then take the form − ∂ t E k + (cid:16) π k L (cid:17) × B k = J S k ,∂ t B k + (cid:16) π k L (cid:17) × E k = 0 for k ∈ (cid:74) − K, K (cid:75) (8)and the trajectory equations for the particles read dd t X p = V p , dd t V p = q p m p ( E S ( X p ) + V p × B S ( X p )) for p ∈ (cid:74) , N (cid:75) . (9)In (8)–(9) we have denoted the coupling terms by • J S k the Fourier coefficients of the current density seen by the discrete field, • E S and B S , the electromagnetic field seen by the particles.These coupling terms involve a shape function S which is primarily used to define an auxiliaryparticle current density J SN ( x ) := J N ∗ S ( x ) = N (cid:88) p =1 q p V p S ( x − X p ) (10)from which the Fourier coefficients J S k are derived. This shape S may either be the Dirac massor some smooth approximation of it such as a B-spline, in which case J SN corresponds to aconvolution smoothing of J N = J δN .How the coupling terms are precisely defined will then characterize the numerical methodat this semi-discrete level. Spectral particle methods can essentially be divided in two classes:gridless Particle-in-Fourier methods where the particles are directly coupled to the fields, and
Particle-in-Cell methods that use an intermediate grid and a Discrete Fourier Transform (DFT)to localize the coupling steps. In the remainder of this section we recall the main features ofthese methods and discuss aliasing errors and anti-aliasing techniques.
The simplest option consists of a gridless coupling as in e.g. [26, 37, 11, 13], which leads to amethod sometimes called
Particle-in-Fourier [30, 1, 29]. Here, the coupling current terms J S k are simply obtained as the Fourier coefficients of J SN . Letting F k ( G ) := (cid:16) L (cid:17) (cid:90) [0 ,L ] G ( x )e − π k · x L d x (11)4enote the k -th (continuous) Fourier coefficient of an arbitrary function G , this gives J S k := F k ( J SN ) = (cid:16) L (cid:17) (cid:88) p =1 ··· N q p V p (cid:90) [0 ,L ] S ( x − X p )e − π k · x L d x . (12)The coupling fields are then defined by continuous convolution products E S ( x ) = (cid:90) [0 ,L ] E K ( ˜ x ) S ( ˜ x − x ) d ˜ x , B S ( x ) = (cid:90) [0 ,L ] B K ( ˜ x ) S ( ˜ x − x ) d ˜ x (13)evaluated at the particle positions x = X p . It will sometimes be convenient to rewrite thisgridless coupling in terms of the Fourier coefficients of the function S X p ( x ) = S ( x − X p ), as J S k = (cid:88) p =1 ··· N q p V p F k ( S X p ) (14)and E S ( X p ) = L (cid:88) k ∈ (cid:74) − K,K (cid:75) E k F k ( S X p ) , B S ( X p ) = L (cid:88) k ∈ (cid:74) − K,K (cid:75) B k F k ( S X p ) , (15)where F denotes the complex conjugate of F . We note that here a Dirac shape can be used, S = δ , since Fourier coefficients are well defined for Dirac distributions [14], F k ( δ X p ) = (cid:16) L (cid:17) e − π k · X pL and this is indeed a standard option in gridless methods [30, 1]. Here we keep the possibility ofusing arbitrary shapes, for the sake of generality and clarity of exposition.In [13, 33] it was shown that this semi-discrete system can be derived from a discretevariational principle, and that it preserves exactly the charge (namely the Gauss laws), as wellas the total energy and momentum of the system. However, the coupling is global: everyparticle contributes directly to every Fourier mode, and vice versa. For problems involving alarge number of Fourier modes, this leads to a computational complexity of O ( N K ) which isprohibitive for simulations using a large number N of particles. A standard approach [24, 25, 21, 3] consists of using an intermediate grid with M points perdimension, M ≥ K + 1, and discrete Fourier transforms (DFT), which is more efficient forsimulations where a large number of Fourier modes are needed. This approach is sometimesreferred to as spectral or pseudo-spectral PIC [11, 16]. Denoting by F M, k ( G ) := (cid:16) M (cid:17) (cid:88) m ∈ (cid:74) ,M (cid:75) G ( m h )e − π k · m M with h = LM , (16)the discrete
Fourier coefficients associated with this grid, the current source is then defined as J S k := F M, k ( J SN ) (17)with a smoothed current density given again by J SN = J N ∗ S , see (10). In practice this amountsto first depositing this current on the grid as in a standard Particle-in-Cell method, J pic m := J SN ( m h ) = (cid:88) p =1 ··· N q p V p S ( m h − X p ) for m ∈ (cid:74) , M (cid:75) J S k = (cid:16) M (cid:17) (cid:88) m ∈ (cid:74) ,M (cid:75) J pic m e − π k · m M for k ∈ (cid:74) − K, K (cid:75) . The pushing fields are then defined by a discrete convolution, E S ( x ) := h (cid:88) m ∈ (cid:74) ,M (cid:75) E K ( m h ) S ( m h − x ) , B S ( x ) := h (cid:88) m ∈ (cid:74) ,M (cid:75) B K ( m h ) S ( m h − x ) (18)also evaluated at the particle positions x = X p . In practice the steps are similar, in a transposedorder: the field values (7) are first computed on the grid, which corresponds to an inverse DFT E pic m := E K ( m h ) = (cid:88) k ∈ (cid:74) − K,K (cid:75) E k e π k · m M B pic m := B K ( m h ) = (cid:88) k ∈ (cid:74) − K,K (cid:75) B k e π k · m M for m ∈ (cid:74) , M (cid:75) , and they are gathered on the particles with the shape function S , as in a standard PIC method E S ( X p ) = h (cid:88) m ∈ (cid:74) ,M (cid:75) E pic m S ( m h − X p ) B S ( X p ) = h (cid:88) m ∈ (cid:74) ,M (cid:75) B pic m S ( m h − X p ) . We note that these coupling terms can be rewritten in a form similar to (14)–(15), now withthe discrete Fourier coefficients of the function S X p ( x ) = S ( x − X p ). Indeed we have J S k = (cid:88) p =1 ··· N q p V p F M, k ( S X p ) (19)and E S ( X p ) = L (cid:88) k ∈ (cid:74) − K,K (cid:75) E k F M, k ( S X p ) , B S ( X p ) = L (cid:88) k ∈ (cid:74) − K,K (cid:75) B k F M, k ( S X p ) . (20)For these terms to be well-defined, we see that S must now be at least continuous. A commonchoice is to take (the periodic extension of) a tensor-product B-spline of degree κ ≥ S ( x ) := (cid:88) r ∈ Z S hκ ( x + r L ) where S hκ ( x ) := (cid:16) h (cid:17) (cid:89) α ∈ (cid:74) , (cid:75) ˆ S κ (cid:16) x α h (cid:17) , (21)with cardinal univariate B-splines defined on the reference grid as ˆ S ( x ) := [ − , ] ( x ) andˆ S κ ( x ) := ˆ S ∗ ˆ S κ − ( x ) = (cid:90) − ˆ S κ − ( x − y ) d y for κ ≥ . As these shape functions have localized supports, supp( S ) = (cid:2) − h (cid:0) κ +12 (cid:1) , h (cid:0) κ +12 (cid:1)(cid:3) , the particlesonly interact with the ( κ + 1) neighbouring grid nodes, which makes the deposition/gathering6teps local. Moreover if M is a power of two, then the DFTs can be efficiently performed withan FFT algorithm, leading to a computational complexity of O ( N ( κ + 1) ) + O ( M log( M ))that is more affordable in simulations involving a large number of particles and Fourier modes.In [29] it is observed that combining DFT couplings with the filtering method of [34] providesan efficient approximation of the gridless method, corresponding to a nonequispaced fast Fouriertransform [32]. This technique will be revisited in Section 2.4 as a natural back-filtering method.More generally, we note that Fourier filtering is commonly used in modern PIC codes in orderto reduce the statistical noise inherent to particle approximations, see e.g. [28, 18].Fourier-PIC coupling often leads to momentum-preserving schemes. In some cases they havebeen shown to be also charge-preserving, see e.g. [16] where the DFT current deposition is seenas a spectral adaptation of the classical Esirkepov method [12]. However, it does not preservethe energy and in general it cannot be derived from a variational principle. As is well known (see e.g. [3, Sec. 8-7]), smooth particle shapes have a low-pass filtering effect inFourier space. This is most easily seen in the gridless case, where the coupling fields (12)–(13)defined by continuous convolution products satisfy J S k = σ k F k ( J N ) for k ∈ (cid:74) − K, K (cid:75) (22)and (cid:40) F k ( E S ) = σ k E k F k ( B S ) = σ k B k for k ∈ Z (23)where we have denoted σ k := L F k ( S ) and set E k := B k := 0 for | k | ∞ > K . Notice that σ k = σ k for symmetric shapes. As smoother functions are associated with faster decreasingspectra, we can clearly see the filtering effect of smooth shapes. Specifically, with the Diracshape S = δ we have σ k = 1 for all k , hence no filtering. With a B-spline (21) of degree κ ∈ N and scale h = L/M , we have σ k = L F k ( S ) = (cid:89) α =1 (cid:18) sinc (cid:16) πk α µ (2 K + 1) (cid:17)(cid:19) κ +1 with µ := M K + 1 ≥ . (24)Here sinc( θ ) := θ sin θ and µ is the oversampling parameter . Since | πk α µ (2 K +1) | < π µ ≤ π for k ∈ (cid:74) − K, K (cid:75) , this makes explicit how modes | k | ∞ ≈ K are damped with “smoother” splines,namely higher degrees and coarser grids.In the Fourier-PIC case a similar filtering effect can be observed, but an additional phe-nomenon enters into play. Indeed the use of a grid leads to the well-known aliasing effect [5], an M -periodization of the discrete Fourier coefficients through the superposition of high frequencymodes, i.e. F M, k ( G ) = (cid:88) (cid:96) ∈ Z F k + (cid:96) M ( G ) for any G ∈ C . Applying this equality to J SN and using again that F k ( J SN ) = σ k F k ( J N ), we find that in theFourier-PIC case the modes of the coupling current (17) read J S k = (cid:88) (cid:96) ∈ Z σ k + (cid:96) M F k + (cid:96) M ( J N ) for k ∈ (cid:74) − K, K (cid:75) . (25)7ere the aliases are the modes corresponding to (cid:96) (cid:54) = 0, all outside the main range (cid:74) − K, K (cid:75) ,since M ≥ K + 1. For the pushing fields (18) the discrete convolution leads to a dual aliasingphenomenon, of the form F k ( E S ) = σ k (cid:88) (cid:96) ∈ Z E k + (cid:96) M F k ( B S ) = σ k (cid:88) (cid:96) ∈ Z B k + (cid:96) M for k ∈ Z . (26)Again, the aliases consist of the (cid:96) (cid:54) = 0 terms, which now correspond to the main modes of E K and B K contributing to higher frequencies of the coupling field. Indeed, using that E k = 0 for k / ∈ (cid:74) − K, K (cid:75) we may rewrite (26) as F k (cid:48) ( E S ) = (cid:40) σ k + (cid:96) M E k for k (cid:48) = k + (cid:96) M ∈ (cid:74) − K, K (cid:75) + M Z k (cid:48) / ∈ (cid:74) − K, K (cid:75) + M Z (27)and similarly for B S .The repercussions of aliasing in numerical simulations have been studied since the earlydays of computational plasma modelling, either through linearized dispersion analysis or fullynonlinear studies [24, 31, 15, 25]. By introducing spurious modes which can then be coupledin the nonlinear models, aliasing has often been recognized as the source of many issues in thesimulations, including grid heating [4], finite grid instabilities [22] and numerical Cherenkovinstabilities [40, 16]. It is clear from (25)–(27) that smooth shapes with a low-pass filtering effect, such as B-splines,may be used for anti-aliasing purposes. However, by filtering also some frequencies withinthe computational range (cid:74) − K, K (cid:75) , they can lead to the overdamping of relevant modes, inparticular for low-resolution discretizations.In order to mitigate the aliasing errors and thus reduce the associated instabilities, a success-ful approach has consisted in associating the anti-aliasing properties of smooth spline shapeswith additional ad-hoc filters. In [16] for instance, the authors show that many instabilitiescan be strongly reduced by using filters determined so as to reduce specific growing modes inthe dispersion relations. And in [29], the particular DFT coupling that is proposed to reducealiasing is based on the nonequispaced fast Fourier transforms (NFFT) [34, 32] which preciselyinvolves filter coefficients that match the low-pass filter effect of the smoothing splines.Here we propose to interpret these filtering techniques as an effective back-filtering method.Indeed, it is easily seen from (25) that a simple solution for the overdamping issue consists ofdividing each deposited current mode with the corresponding shape filter coefficient, leading toa new current defined as J S, bf k := 1 σ k F M, k ( J SN ) . (28)For the pushing fields the idea is the same but we see from (26) that the back-filtering needsto be applied on the original field, rather than on the coupling terms. This leads to setting E S, bf ( x ) := h (cid:88) m ∈ (cid:74) ,M (cid:75) E bf K ( m h ) S ( m h − x ) , B S, bf ( x ) := h (cid:88) m ∈ (cid:74) ,M (cid:75) B bf K ( m h ) S ( m h − x ) (29)8ith back-filtered fields defined as E bf K ( x ) := (cid:88) k ∈ (cid:74) − K,K (cid:75) (cid:16) σ k (cid:17) E k e π k · x L B bf K ( x ) := (cid:88) k ∈ (cid:74) − K,K (cid:75) (cid:16) σ k (cid:17) B k e π k · x L . (30)With this coupling, formulas (25)–(26) become J S, bf k = F k ( J N ) + (cid:88) (cid:96) (cid:54) =0 σ k + (cid:96) M σ k F k + (cid:96) M ( J N ) for k ∈ (cid:74) − K, K (cid:75) (31)and F k ( E S, bf ) = E k + (cid:88) (cid:96) (cid:54) =0 (cid:16) σ k σ k + (cid:96) M (cid:17) E k + (cid:96) M F k ( B S, bf ) = B k + (cid:88) (cid:96) (cid:54) =0 (cid:16) σ k σ k + (cid:96) M (cid:17) B k + (cid:96) M for k ∈ Z (32)where we have separated each contribution into its main mode ( (cid:96) = 0) and the filtered aliases.And again, using explicitly that E k = 0 for k / ∈ (cid:74) − K, K (cid:75) we can rewrite (32) as F k (cid:48) ( E S, bf ) = (cid:16) σ k + (cid:96) M σ k (cid:17) E k if k (cid:48) = k + (cid:96) M ∈ (cid:74) − K, K (cid:75) + M Z k (cid:48) / ∈ (cid:74) − K, K (cid:75) + M Z (33)and similarly for B S, bf . For B-splines we can see from (24) that σ k is far from 0 in the range k ∈ (cid:74) − K, K (cid:75) . Thus, back-filtering allows to reduce the amplitude of all the aliased modesin a similar proportion as with a standard filtering, but without damping any mode in thecomputational range. In this section we follow the variational structure-preserving discretization framework of [9] andapply it to discrete Fourier spaces. This essentially allows us to extend the Finite Elementspline GEMPIC method from [23] to spectral Maxwell solvers.
We remind that a central feature of this framework is to preserve the de Rham sequence of R , H ( R ) ∇−−→ H (curl; R ) curl −−−−−→ H (div; R ) div −−−−→ L ( R ) (34)at the discrete level, and to admit a sequence of projection operators Π , . . . , Π mappinginfinite-dimensional function spaces into the discrete ones, such that the following diagramcommutes: V V V V V K V K V K V K Π ∇ Π ∇ Π Π curlcurl divdiv (35)9e point out that such commuting de Rham diagrams are a key tool in Finite Element ExteriorCalculus (FEEC), see e.g. [6, 20, 2, 7, 8]. In our framework, it is these operators Π (cid:96) , togetherwith the shape functions S , that will encode the coupling mechanism between the particles andthe discrete fields. The bottom row thus consists of truncated Fourier spaces V K := V K := Span (cid:0) { Λ k : k ∈ (cid:74) − K, K (cid:75) } (cid:1) V K := V K := Span (cid:0) { Λ k e α : k ∈ (cid:74) − K, K (cid:75) , α ∈ (cid:74) , (cid:75) } (cid:1) where e α is the unit basis vector in the α dimension, andΛ k ( x ) := e π k · x L , x ∈ [0 , L ] is the Fourier basis function of index k ∈ Z . Using the fact that the discrete spaces V K and V K coincide in this spectral setting, we find that the weak discrete differential operators associatedwith the sequence (35) coincide with the strong (continuous) ones. As a consequence, thevariational method derived in [9] takes the following form: the field equations read (cid:40) − ∂ t E K + curl B K = Π J SN ∂ t B K + curl E K = 0 with Π J SN = (cid:88) p =1 ··· N q p Π ( V p S X p ) (36)with B K , E K in V K = V K , and the particle trajectories read dd t X p = V p dd t V p = q p m p (cid:0) E S ( X p ) + V p × B S ( X p ) (cid:1) with E Sα ( X p ) = (cid:90) E α ( x )(Π α S X p )( x ) d x B Sα ( X p ) = (cid:90) B α ( x )(Π α S X p )( x ) d x (37)for p ∈ (cid:74) , N (cid:75) and α ∈ (cid:74) , (cid:75) , with E α and B α the α -components of E K and B K .As for the operators Π (cid:96) , several choices can then be made that lead to a commuting diagramand each choice results in a different coupling mechanism between the particles and the fields.In [9] a set of projections was presented that is based on interpolation and hispolation, which inpractice amounts to performing discrete Fourier transforms (DFT) on a grid with M = 2 K + 1nodes but also involves surface and volume integrals of the particle shapes. In this article wewill study two types of approximation operators: L projection operators corresponding tocontinuous Fourier transforms, as recalled in Section 3.4, and new pseudo-differential operatorsbased on discrete Fourier transforms, that we present in Section 3.5. These choices will thenrespectively lead to two variational schemes with Hamiltonian structure, namely the gridlessParticle-in-Fourier method, and a new spectral Fourier-GEMPIC method. Remark 1.
Although we often follow the common usage and refer to Π (cid:96) as commuting projec-tion operators, we emphasize that we do not require them to be actual projections in the sensethat one would have Π (cid:96) = I on V (cid:96)K . Indeed this property is not needed for the Hamiltonianstructure of the resulting schemes [9], and by relaxing it we can directly extend our analysis tocoupling methods that involve filtering or back-filtering mechanisms. In [9], a variant of the abstract system (36)–(37) was also proposed, that is a priori not Hamilto-nian but preserves both the Gauss laws and a discrete momentum. This modified system involves10he same operators from the general commuting diagram (35), and some interior products cou-pled with a dimension-dependent approximation operator A h,α . Specifically, the momentum-preserving variant consists of the same field solver (36) as above, and of a modified particlepusher where the Lorentz term F S ( X p , V p ) = E S ( X p ) + V p × B S ( X p ) from (37) is replaced by˜ F S ( X p , V p ) := ˜ E S ( X p ) + ˜ R S ( X p , V p )with components given by ˜ E Sα ( X p ) = (cid:90) E α ( x )( A h,α Π S X p )( x ) d x ˜ R Sα ( X p , V p ) = (cid:90) B K ( x ) · (cid:0) A h,α ( e α × Π ( V p S X p ) (cid:1) ( x ) d x . (38)In the general construction of [9] the operator A h,α was defined as a directional averaging ona grid with M points, h = L/M , which in Fourier space amounts to a diagonal filtering ofthe form ( A h,α G ) k = sinc (cid:0) πk α M (cid:1) G k . This allows for the method to be well posed in generalpolynomial or spline finite element settings. With spectral solvers however, this averaging is notneeded and one may simply take the identity operator, ( A h,α G ) k = G k . We may then rewritethese pushing fields in terms of Fourier coefficients, for a clearer comparison. The fields in theHamiltonian pusher (37) take the form E Sα ( X p ) = L (cid:88) k E α, k (Π α S X p ) k R Sα ( X p , V p ) = (cid:88) ν = ± νV p,α + ν B α − ν ( X p ) = L (cid:88) k (cid:88) ν = ± νB α − ν, k V p,α + ν (Π α − ν S X p ) k , (39)whereas the momentum-preserving ones read ˜ E Sα ( X p ) = L (cid:88) k E α, k (Π S X p ) k ˜ R Sα ( X p , V p ) = L (cid:88) k (cid:88) ν = ± νB α − ν, k V p,α + ν (Π α + ν S X p ) k . (40)In particular, we observe that the two schemes differ in the particle-field coupling operatorsinvolved in the pushing fields. After performing a convenient time discretization, we will see inSection 4.6 that this latter formulation will result in fully discrete PIC schemes with standardparticle-field coupling terms, as described in Section 2.2. In [9] we have shown that the semi-discrete equations (36)–(37) have a discrete Hamiltonianstructure, which in particular implies that they preserve the total energy and any discreteCasimir functional. In the particular setting of truncated Fourier spaces, a few basic conserva-tion properties can be proven with a direct argument.
Theorem 1.
The variational spectral particle method (36) – (37) preserves the strong Gauss laws (cid:40) div E K = Π ρ SN div B K = 0 (41)11 s well as the discrete energy H ( t ) = 12 N (cid:88) p =1 m p | V p ( t ) | + 12 (cid:90) [0 ,L ] (cid:16) | E K ( t, x ) | + | B K ( t, x ) | (cid:17) d x . (42) If in addition the projection operators satisfy Π α = Π α = Π for ≤ α ≤ , (43) then it also preserves the discrete momentum P ( t ) = N (cid:88) p =1 m p V p ( t ) + (cid:90) [0 ,L ] E K ( t, x ) × B K ( t, x ) d x . (44) Proof.
The preservation of the magnetic Gauss law readily follows from the strong Faradayequation in (36). As for the electric Gauss law, we compute for an arbitrary smooth function ϕ dd t (cid:90) ρ SN ( t, x ) ϕ ( x ) d x = N (cid:88) p =1 q p (cid:90) S ( ˜ x ) V p · ∇ ϕ ( ˜ x + X p ) d ˜ x = (cid:90) J SN ( t, x ) · ∇ ϕ ( x ) d x which shows that the continuity equation ∂ t ρ SN + div J SN = 0 (45)always holds in distribution’s sense. Taking next the divergence of the discrete Amp`ere equationin (36) the commuting diagram property (35) allows us to write ∂ t div E K = − div Π J SN = − Π div J SN = ∂ t Π ρ SN where the last equality follows from (45) and from the time-invariance of the operator Π .Integrating over time this shows that the electric Gauss law is indeed preserved. To show theenergy conservation, we next compute using (36)dd t (cid:16) (cid:90) | E K | + | B K | (cid:17) = (cid:90) E K · ( ∇ × B K − Π J SN ) − B K · ∇ × E K = − (cid:90) E K · Π J SN and, using (37),dd t (cid:16) N (cid:88) p =1 m p | V p | (cid:17) = N (cid:88) p =1 q p (cid:90) Π ( V p S X p ) · E K = (cid:90) E K · Π J SN , so that H is indeed constant over time. We finally turn to the momentum conservation andassume that (43) holds. Then the coupling fields take the form G S ( X p ) = (cid:90) G K ( x )(Π S X p )( x ) d x with G = E or B , and the deposited current reads Π J SN = (cid:80) p q p Π ( V p S X p ) = (cid:80) p q p V p Π ( S X p ) . Using (37) wecomputedd t N (cid:88) p =1 m p V p = N (cid:88) p =1 q p (cid:90) (cid:0) E K Π S X p + V p × B K Π S X p (cid:1) = (cid:90) E K Π ρ SN + (cid:90) (Π J SN ) × B K , (cid:82) Ω G × (curl G ) = (cid:82) Ω ∇ ( G ) − ( G ·∇ ) G = (cid:82) Ω (div G ) G valid for an arbitrary function G , allows us to computedd t (cid:90) Ω E K × B K = (cid:90) Ω (curl B K − Π J SN ) × B K − (cid:90) Ω E K × (curl E K )= − (cid:90) Ω (div B K ) B K − (cid:90) Ω Π J SN × B K − (cid:90) Ω (div E K ) E K = − (cid:90) Ω (Π J SN ) × B K − (cid:90) Ω (Π ρ SK ) E K where in the last equality we have used the Gauss laws (41). The result follows. (cid:3) Because truncated Fourier spaces have the particular property that they are stable under spacedifferentiation, L projection operators can be used for the commuting diagram. This choiceessentially corresponds to a gridless Particle-in-Fourier coupling described above, and in thisarticle we will refer to the resulting method as a Geometric Electromagnetic Particle-in-Fourier (GEMPIF) method, to emphasize its natural expression in the general GEMPIC framework. Inorder to account for general filtering and back-filtering mechanisms, we consider an arbitrarycollection of Hermitian coefficients γ k = γ − k ∈ C , and setΠ := Π α := Π α := Π for α ∈ (cid:74) , (cid:75) , (46)where Π is the operator that maps a function to its γ -filtered Fourier series of rank K ,Π G := (cid:88) k ∈ (cid:74) − K,K (cid:75) γ k F k ( G )Λ k , (47)with F k the continuous Fourier coefficient operator, see (11). We observe that with unit filters γ k = 1 the operators in (46) coincide with the L projection P K on V = V , characterized by (cid:90) [0 ,L ] P K G ( x )Λ k ( x ) d x = (cid:90) [0 ,L ] G ( x )Λ k ( x ) d x for k ∈ (cid:74) − K, K (cid:75) . For general filter coefficients we have the following result.
Lemma 1.
The filtered L projection operators Π (cid:96) defined by (46) – (47) satisfy the commutingdiagram property (35) in a distribution sense, with periodic distributions as domain spaces: V = V = D (cid:48) per and V = V = ( D (cid:48) per ) . Proof.
We may consider that γ k = 1, as the general case follows easily. The L projections overtruncated Fourier spaces classically extend to periodic distributions [14], writing e.g. F k (cid:0) ( ∂ α ) a δ X p (cid:1) = ( − a L (cid:90) [0 ,L ] δ X p ( ∂ α ) a Λ k ( x ) d x = (cid:16) L (cid:17) (cid:16) πk α L (cid:17) a e − π k · X pL for α ∈ (cid:74) , (cid:75) and a ∈ N . The commuting diagram property is then derived from the symmetryof the discrete de Rham sequence. For instance, using that div V = div V ⊂ V = V allowsto write for all F ∈ V (cid:90) ( ∇ Π G ) · F = − (cid:90) (Π G ) div F = − (cid:90) G div F = (cid:90) ( ∇ G ) · F = (cid:90) (Π ∇ G ) · F G ∈ D (cid:48) per . Since ∇ V ⊂ V this shows that ∇ Π G = Π ∇ G . Theother relations curl Π G = Π curl G and div Π G = Π div G are proven in the same way. (cid:3) In our framework, spectral PIC methods are obtained with commuting projection operators thatinvolve discrete Fourier transforms on a finite grid. In [9] we have described a set of projectionswhich rely on geometric degrees of freedom, namely nodal interpolations for V and edge, face,and volume “histopolations” for V , V , and V , respectively. As a result the current depositioninvolves face integrals which need to be integrated over the particle trajectories, which somehowdeviates from standard PIC methods where common deposition procedures are based on pointevaluations of the particle shapes. For this reason we consider here an alternate constructionbased on a new sequence of projection operators, obtained by combining DFT coefficients on agrid with M ≥ K + 1 nodes as in Section 2.2 with standard derivatives and anti-derivativesin Fourier variables. As we will see, these new projections will lead to deposition methods thatonly involve pointwise evaluations of the particle shapes. Following the terminology of [23],we will refer to the resulting methods as Fourier-Geometric Electromagnetic Particle-in-Cell (Fourier-GEMPIC) methods. Beginning with the space V , we letΠ G = (cid:88) k ∈ (cid:74) − K,K (cid:75) (Π G ) k Λ k (48)be defined on C the space of continuous, L -periodic functions, by its coefficients(Π G ) k := γ k ˜ F M, k ( G ) := γ k ˜ F M,k ⊗ ˜ F M,k ⊗ ˜ F M,k ( G ) . (49)Here, the values γ k = γ − k ∈ C are again Hermitian filters, and the univariate operators ˜ F M,k α are defined as˜ F M,k α ( G ) := (cid:40) L (cid:82) L G ( x α ) d x α if k α = 0 F M,k α ( G ) = M (cid:80) Mm α =1 G ( m α h )e − πkαmαM if k α (cid:54) = 0 . (50)We note that for γ = 1 the operator Π is a conservative discrete Fourier transform, indeed (cid:90) [0 ,L ] Π G ( x ) d x = L (Π G ) = L ˜ F M, ( G ) = (cid:90) [0 ,L ] G ( x ) d x . For the vector-valued V we then setΠ G = (cid:88) α =1 (cid:88) k ∈ (cid:74) − K,K (cid:75) (Π α G α ) k Λ k e α with (Π α G α ) k := γ k ( ˆ D k ,α ) − ˜ F M, k ( ˜ D k ,α G α ) , with pseudo-differential operators defined asˆ D k ,α c k := (cid:40) c k if k α = 0 πk α L c k if k α (cid:54) = 0 and ˜ D k ,α G := (cid:40) G if k α = 0 ∂ α G if k α (cid:54) = 0 (51)for any c k ∈ C and any function G in the anisotropic regularity space C , × C , × C , ,where we have denoted C ,α := { G ∈ C : ∂ α G ∈ C } . (52)14imilarly for the vector-valued space V we define (using a circular convention α ≡ α + 3 forthe dimension indices) Π G = (cid:88) α =1 (cid:88) k ∈ (cid:74) − K,K (cid:75) (Π α G α ) k Λ k e α with (Π α G α ) k := γ k ( ˆ D k ,α +1 ˆ D k ,α +2 ) − ˜ F M, k ( ˜ D k ,α +1 ˜ D k ,α +2 G α )for any function G in the anisotropic regularity space C , , × C , , × C , , , where C ,α,β := { G ∈ C : ∂ α ∂ β G ∈ C } . (53)Finally for the scalar-valued space V we define Π G = (cid:88) k ∈ (cid:74) − K,K (cid:75) (Π G ) k Λ k with (Π G ) k := γ k ( ˆ D k , ˆ D k , ˆ D k , ) − ˜ F M, k ( ˜ D k , ˜ D k , ˜ D k , G ) , for all G in the anisotropic regularity space C , , , := { G ∈ C : ∂ ∂ ∂ G ∈ C } . (54) Lemma 2.
The above pseudo-differential operators Π (cid:96) satisfy the commuting diagram property (35) with domains defined as V = C , , , , V = C , , × C , , × C , , , V = C , × C , × C , , V = C . Moreover if γ k = 1 , they are projection operators on their respective Fourier spaces.Proof. The projection property can be checked by direct computation. To verify the commutingdiagram property we consider again the case γ k = 1, as the general case follows easily. We beginwith the last relation and observe that˜ F M, k ( ∂ α G α ) = πk α L (Π G ) k ,α holds for every dimension α ∈ (cid:74) , (cid:75) and mode k ∈ Z : this follows from the definition of Π if k α (cid:54) = 0, and from that of the conservative Fourier transform ˜ F M, k if k α = 0. We then have, foran arbitrary k ,(Π div G ) k = (cid:88) α (Π ( ∂ α G α )) k = (cid:88) α ˜ F M, k ( ∂ α G α ) = (cid:88) α πk α L (Π G ) k ,α = (div Π G ) k which shows that Π div = div Π holds on V . Similarly, for all k and α , we have(Π curl G ) k ,α = (cid:0) Π α ( ∂ α +1 G α +2 − ∂ α +2 G α +1 ) (cid:1) k = ( ˆ D k ,α ) − ˜ F M, k (cid:0) ˜ D k ,α ( ∂ α +1 G α +2 − ∂ α +2 G α +1 ) (cid:1) = πk α +1 L (Π G ) k ,α +2 − πk α +2 L (Π G ) k ,α +1 = (curl Π G ) k ,α so that Π curl = curl Π holds on V . Finally we compute, again for all k , α ,(Π ∇ G ) k ,α = (Π α ∂ α G ) k = ( ˆ D k ,α +1 ˆ D k ,α +2 ) − ˜ F M, k ( ˜ D k ,α +1 ˜ D k ,α +2 ∂ α G )= πk α L (Π G ) k = ( ∇ Π G ) k ,α which shows that Π ∇ = ∇ Π holds on V and ends the proof. (cid:3) .6 GEMPIF and Fourier-GEMPIC methods Now that we have defined two sets of projection operators Π (cid:96) with commuting diagram prop-erties, we may specify the spectral variational scheme (36)–(37). As announced above, theGeometric Particle-in-Fourier (GEMPIF) method corresponds to the case where the operatorsΠ (cid:96) are defined as the L projections of Section 3.4. Then any periodic distribution S ∈ D (cid:48) per isadmissible and the coupling terms take the form J S k = γ k (cid:88) p =1 ··· N q p V p F k ( S X p ) = γ k F k ( S ) (cid:88) p =1 ··· N q p V p e − π k · X pL (55)and E Sα ( X p ) = L (cid:88) k ∈ (cid:74) − K,K (cid:75) γ k E α, k F k ( S X p ) = L (cid:88) k ∈ (cid:74) − K,K (cid:75) γ k F k ( S ) E α, k e π k · X pL B Sα ( X p ) = L (cid:88) k ∈ (cid:74) − K,K (cid:75) γ k B α, k F k ( S X p ) = L (cid:88) k ∈ (cid:74) − K,K (cid:75) γ k F k ( S ) B α, k e π k · X pL . (56)Up to the arbitrary filter coefficients γ k , this corresponds to the gridless coupling of Section 2.1.With the pseudo-differential operators defined in Section 3.5, the coupling terms involvemodified Fourier coefficients such as (Π α S X p ) k = ( ˆ D k ,α ) − ˜ F M, k ( ˜ D k ,α S X p ), see (51), and similarcoefficients for Π α S X p , involving derivatives along dimensions α + 1 and α −
1. The resultingcoupling terms read then J Sα, k = γ k (cid:16) πk α L (cid:17) − (cid:16) M (cid:17) (cid:88) p, m q p v p,α ∂ α S ( m h − X p )e − π k · m M (57)and E Sα ( X p ) = γ k h (cid:88) k , m E α, k (cid:16) − πk α L (cid:17) − ∂ α S ( m h − X p )e π k · m M B Sα ( X p ) = γ k h (cid:88) k , m B α, k (cid:16) − πk α +1 L (cid:17) − (cid:16) − πk α − L (cid:17) − ∂ α +1 ∂ α − S ( m h − X p )e π k · m M (58)with the modifications specified in (49)–(51) when k α = 0 or k α ± = 0. We observe that in orderto be admissible, shape functions need to be in the mixed C space (54). For tensor-productB-splines this corresponds to using at least quadratic splines in each dimension.Although these terms involve a DFT grid, they differ from the standard spectral PIC cou-pling terms recalled in Section 2.2, which are not associated with a variational principle. Forthe purpose of comparison we rewrite the latter in the case of a general γ -filtering, J Sα, k = γ k (cid:16) M (cid:17) (cid:88) p, m q p v p,α S ( m h − X p )e − π k · m M (59)and E Sα ( X p ) = γ k h (cid:88) k , m E α, k S ( m h − X p )e π k · m M B Sα ( X p ) = γ k h (cid:88) k , m B α, k S ( m h − X p )e π k · m M . (60)Below we will see how these differences translate in a fully discrete setting.16 Fully discrete schemes based on Hamiltonian time splitting
In this section we specify some fully discrete schemes for the variational semi-discrete systemsdescribed above. As these schemes are derived from a Hamiltonian splitting procedure, it willbe convenient to first rewrite the general equations (36)–(37) in a matrix form.
Similarly as in [9, 23], we gather the degrees of freedom of the discrete solution into multi-indexarrays. Particle unknowns will be written as X = ( X p,α ) ( p,α ) ∈ (cid:74) ,N (cid:75) × (cid:74) , (cid:75) , V = ( V p,α ) ( p,α ) ∈ (cid:74) ,N (cid:75) × (cid:74) , (cid:75) ∈ R N , (with an implicit time-dependance) and electro-magnetic field coefficients will be denoted by E = ( E α, k ,σ ) ( α, k ,σ ) ∈ (cid:74) , (cid:75) × (cid:74) − K,K (cid:75) × (cid:74) , (cid:75) , B = ( B α, k ,σ ) ( α, k ,σ ) ∈ (cid:74) , (cid:75) × (cid:74) − K,K (cid:75) × (cid:74) , (cid:75) ∈ R K +1) where each component-wise Fourier coefficient is decomposed into its real and imaginary parts, G α, k = G α, k , + i G α, k , with G = E or B. For notational reasons it is convenient to see these arrays as column vectors (with an arbi-trary ordering of the multi-indices), and to gather them into a global array of time-dependentunknowns, U = XVEB . (61)We may then rewrite the equations of the geometric Fourier-particle method (36)–(37) in termsof these coefficients. Using the fact that the vector-valued operators Π (cid:96) are defined component-wise, in the sense that their component along any dimension α ∈ (cid:74) , (cid:75) reads(Π (cid:96) ( G )) α = Π (cid:96)α ( G α ) (62)for some scalar operators Π (cid:96)α , we may rewrite the particle trajectories as dd t X p,α = V p,α dd t V p,α = L q p m p (cid:88) k (cid:16) E α, k (Π α S X p ) k ± V p,α ± B α ∓ , k (Π α ∓ S X p ) k (cid:17) = L q p m p (cid:88) k ,σ (cid:16) E α, k ,σ (Π α S X p ) k ,σ ± V p,α ± B α ∓ , k ,σ (Π α ∓ S X p ) k ,σ (cid:17) (63)for all p ∈ (cid:74) , N (cid:75) and α ∈ (cid:74) , (cid:75) . Here we have used the fact that the right-hand side is real,and we recall the circular convention ( α ≡ α + 3) for the dimension indices. In terms of theabove arrays, this gives (cid:40) dd t X = V dd t V = W qm (cid:0) S ( X ) M E + R ( X , B ) V (cid:1) where W qm = diag( q p m p : p ∈ (cid:74) , N (cid:75) ) is the diagonal weighting matrix carrying the particlescharge to mass ratios, S ( X ) is the Π Fourier-particle coupling matrix, S ( X ) ( p,β ) , ( α, k ,σ ) = δ α,β (Π α S X p ) k ,σ , (cid:96) = L I N (cid:96) is the (diagonal) finite element “mass” matrix for the Fourier basis in V (cid:96)K , and R ( X , B ) is the skew-symmetric matrix corresponding to magnetic rotation, R ( X , B ) ( p,α ) , ( p (cid:48) ,β ) = δ p,p (cid:48) (cid:88) k ,σ L (cid:0) δ α +1 ,β B α − , k ,σ (Π α − S X p ) k ,σ − δ α − ,β B α +1 , k ,σ (Π α +1 S X p ) k ,σ (cid:1) = δ p,p (cid:48) (cid:88) ν = ± δ α + ν,β (cid:90) νB α − ν ( x )(Π α − ν S X p )( x ) d x . (64)As for the Maxwell solver, we first rewrite (36) in terms of (complex) Fourier coefficients, namely dd t E α, k = (cid:88) ν = ± ν πL k α + ν B α − ν, k − (cid:88) p q p v p,α Π α ( S X p ) k dd t B α, k = − (cid:88) ν = ± ν πL k α + ν E α − ν, k (65)for k ∈ (cid:74) − K, K (cid:75) and α ∈ (cid:74) , (cid:75) . In terms of the above arrays, this rewrites as (cid:40) dd t E = C B − S ( X ) T W q V dd t B = − C E (66)with C the symmetric matrix of the curl operator, which writes here C ( α, k ,σ ) , ( β, (cid:96) ,τ ) = δ k , (cid:96) (1 − δ σ,τ )( − σ +1 πL (cid:0) δ α − ,β k α +1 − δ α +1 ,β k α − (cid:1) . Writing next the discrete Hamiltonian (42) as a function of the array variable (61), H ( U ) = V (cid:62) W m V + E (cid:62) M E + B (cid:62) M B (67)we find for the corresponding derivatives ∇ U H ( U ) = ∇ X H∇ V H∇ E H∇ B H ( U ) = W m V M E M B and this allows us to rewrite the abstract spectral particle method (36)–(37) in the form of anon-canonical Hamiltonian system dd t U = J ( U ) ∇ U H ( U ) (68)with J ( U ) = J ( X , B ) = W m − W m W qm R ( X , B ) W m W qm S ( X ) 00 − S ( X ) T W qm C ( M ) − − ( M ) − C . (69)In [9] we have shown that J is a Poisson matrix in the sense of [17, Def. VII.2.4], i.e., it isskew-symmetric and it satisfies the matrix Jacobi identity. In particular, System (68) may berewritten in the usual form dd t U = { U , H} with a discrete Poisson bracket given by {F , G} = ( ∇ U F ) T J ∇ U G . .2 Hamiltonian splitting time discretization Following [17, 10], we now apply a splitting procedure to the semi-discrete Hamiltonian sys-tem (68). This will provide us with a series of Hamiltonian structure-preserving schemes ofvarious orders in time, for the variational Fourier-particle equations (36)–(37).Similarly as in [19, 23], we split the kinetic part along the three different dimensions, butkeep together the electric and magnetic parts as we are solving the Maxwell equations in Fourierspaces. This leads us to the following Hamiltonian splitting, H = (cid:88) α =1 H v α + H EB (70)with H v α ( U ) := N (cid:88) p =1 m p | V p,α | and H EB ( U ) := 12 (cid:90) Ω (cid:16) | E K ( x ) | + | B K ( x ) | (cid:17) d x where we remind that U carries the time dependent coefficients of the full solution, see (61).This splitting has two key properties. It leads to split steps that can all be solved exactly, andit preserve the fact that J ( X , B ) is a Poisson matrix, see Theorem 2. As a consequence, weknow that any combination of the split steps will provide a Hamiltonian time scheme whichpreserves the Casimir invariants and the total energy up to some constant time discretizationerror, see e.g. [17]. Specifically, we decompose System (68) into the subsystems, dd t U ( t ) = J ( U ) ∇ U H v α ( U ) for α ∈ (cid:74) , (cid:75) (71)and dd t U ( t ) = J ( U ) ∇ U H EB ( U ) . (72)Denoting by ϕ τ,v α and ϕ τ,EB the corresponding solution flow maps, we can use standard com-position methods to obtain time integrators of various orders, as described e.g. in [23]: eithera first order Lie-Trotter schemeΦ ∆ t,L = ϕ ∆ t,EB ◦ ϕ ∆ t,v ◦ ϕ ∆ t,v ◦ ϕ ∆ t,v or a second-order Strang scheme Φ ∆ t,S = ϕ ∆ t/ ,L ◦ ϕ ∗ ∆ t/ ,L (73)where ϕ ∗ τ,L := ϕ − − τ,L denotes the adjoint flow. Since each split flow is the exact solution of anautonomous system they are all self-adjoint, i.e. symmetric, which yieldsΦ ∆ t,S = ϕ ∆ t/ ,EB ◦ ϕ ∆ t/ ,v ◦ ϕ ∆ t/ ,v ◦ ϕ ∆ t,v ◦ ϕ ∆ t/ ,v ◦ ϕ ∆ t/ ,v ◦ ϕ ∆ t/ ,EB . Similarly we can use a fourth-order Suzuki-Yoshida scheme ϕ ∆ t,S = ϕ γ ∆ t,S ◦ ϕ γ ∆ t,S ◦ ϕ γ ∆ t,S (74)with γ = 1 / (2 − / ) and γ = − / / (2 − / ), or higher-order composition methods, seee.g. [27, 17].In the sections below we specify the resulting equations for each split step, and we providean explicit solution for the GEMPIF and Fourier-GEMPIC methods described in Section 3.6.19 .3 Discrete Hamiltonian subsystems Before giving the solutions we detail the subsystems (71) and (72).
Kinetic H v α subsystems. For each dimension α ∈ (cid:74) , (cid:75) , System (71) reads dd t X = V [ α ]dd t V = W qm R ( X , B ) V [ α ]dd t E = − S ( X ) T W q V [ α ]dd t B = 0 i.e., dd t X p = V [ α ] p dd t V p = q p m p V [ α ] p × B S ( X p ) dd t E K = − Π ( J S, [ α ] N ) dd t B K = 0(for all p ), where we have denoted V [ α ] p = e α V p,α and similarly for V [ α ] , J S, [ α ] N . Expressedin scalar coefficients and using the definition (37) of B S , this gives (again with the circularconvention α ≡ α + 3 on dimension indices) dd t X p,α = V p,α dd t X p,α + ν = 0 dd t V p,α = 0 dd t V p,α + ν = − ν q p m p L (cid:88) k V p,α B α − ν, k (Π α − ν S X p ) k dd t E α, k = − (cid:88) p q p V p,α (Π α S X p ) k dd t E α + ν, k = 0 dd t B k = 0 (75)for all p ∈ (cid:74) , N (cid:75) , k ∈ (cid:74) − K, K (cid:75) and ν = ± Electromagnetic H EB subsystems. For the electromagnetic part, System (72) reads dd t X = 0 dd t V = W qm S ( X ) M E dd t E = C B dd t B = − C E i.e., dd t X p = 0 dd t V p = q p m p E S ( X p ) dd t E K = curl B K dd t B K = − curl E K for all p ∈ (cid:74) , N (cid:75) . Using scalar coefficients and the definition (37) of E S , this reads dd t X p = 0 dd t V p,α = q p m p L (cid:88) k E α, k (Π α S X p ) k dd t E k = π k L × B k dd t B k = − π k L × E k (76)for all α ∈ (cid:74) , (cid:75) , p ∈ (cid:74) , N (cid:75) and k ∈ (cid:74) − K, K (cid:75) . This splitting enjoys the following properties.20 heorem 2. All the split steps above preserve the discrete Gauss laws (41) , namely div E K = Π ρ N = N (cid:88) p =1 q p Π S X p and div B K = 0 (77) and the fact that J ( X , B ) is a Poisson bracket. In particular, any combination of the individualflows ϕ τ,v α and ϕ τ,EB preserves the discrete Hamiltonian structure. If in addition the operators Π (cid:96) satisfy (43) , then the discrete momentum (44) is also preserved exactly.Proof. By applying the computations from the proof of Theorem 1 to any of the split steps, oneverifies that the discrete Gauss laws are preserved, as well as the discrete momentum in thecase where (43) holds. In particular the Gauss law div B K = 0 is preserved, so that Theorem 1from [9] applies and this shows that J defines a discrete Poisson bracket at each split step. Thestandard theory of Hamiltonian splitting schemes then applies, see [17]. (cid:3) In the GEMPIF method the operators Π (cid:96) are defined as L projections with general filtercoefficients γ k , see Section 3.4. In particular we have(Π α S X p ) k = (Π α S X p ) k = (cid:16) L (cid:17) γ k (cid:90) Ω S ( x − X p )e − π k · x L d x = γ k F k ( S )e − π k · X pL (78)for any direction α ∈ (cid:74) , (cid:75) . This allows us to give explicit solutions for each split subsystem. Directional kinetic steps.Lemma 3.
In the GEMPIF method, the exact solution ϕ τ,v α : U → U ( τ ) of the kinetic splitstep (75) in a direction α ∈ (cid:74) , (cid:75) is given by the explicit expressions ϕ τ,v α : X p,α ( τ ) = X p,α + τ V p,α X p,α + ν ( τ ) = X p,α + ν V p,α ( τ ) = V p,α V p,α + ν ( τ ) = V p,α + ν − ντ q p m p L (cid:88) k ∈ (cid:74) − K,K (cid:75) B α − ν, k (cid:16) γ k F k ( S ) ˆ V p,α, k ( τ ) (cid:17) E α, k ( τ ) = E α, k − τ γ k F k ( S ) (cid:88) p =1 ··· N q p ˆ V p,α, k ( τ ) E α + ν, k ( τ ) = E α + ν, k B K ( τ ) = B K . for all p ∈ (cid:74) , N (cid:75) , k ∈ (cid:74) − K, K (cid:75) and ν = ± , with ˆ V p,α, k ( τ ) := V p,α e − π k · X pL if k α = 0 (cid:16) − πk α τL (cid:17) − (cid:16) e − πkατVp,αL − (cid:17) e − π k · X pL if k α (cid:54) = 0 . (79)21 roof. In (75), the velocity equations readdd t V p,α + ν = − ν q p m p L (cid:88) k B α − ν, k (cid:16) γ k F k ( S ) V p,α e − π k · X pL (cid:17) and for the electric field we havedd t E k ,α = − γ k F k ( S ) (cid:88) p q p V p,α e − π k · X pL . Since the only time-varying term in the right hand sides is dd t X p ( t ) = V [ α ] p = V p,α e α , weintegrate (cid:90) τ V p,α e − π k · X p ( t ) L d t = τ ˆ V p,α, k ( τ )with the expression in (79), which proves the lemma. (cid:3) Electromagnetic step.
In (76) the source-free Maxwell equations have an explicit solution(see e.g. [35]), which allows to solve also for the particles. The resulting flow takes the form ϕ τ,EB : X p ( τ ) = X p V p,α ( τ ) = V p,α + q p m p L (cid:88) k ∈ (cid:74) − K,K (cid:75) (cid:16) (cid:90) τ E α, k ( t ) d t (cid:17) (Π α S X p ) k E k ( τ ) = E k + ˆ k × (cid:16) (1 − c ( τ ))(ˆ k × E k ) + i s ( τ ) B k (cid:17) B k ( τ ) = B k + ˆ k × (cid:16) (1 − c ( τ ))(ˆ k × B k ) − i s ( τ ) E k (cid:17) (80)for all p ∈ (cid:74) , N (cid:75) , α ∈ (cid:74) , (cid:75) and k ∈ (cid:74) − K, K (cid:75) . Here the projected shaped particle is given by(78), and the integrated electric field reads (cid:90) τ E k ( t ) d t = τ E k + ˆ k × (cid:16)(cid:0) τ − L π | k | s ( τ ) (cid:1) (ˆ k × E k ) + i L π | k | (1 − c ( τ )) B k (cid:17) (81)where we have set c ( τ ) = cos (cid:0) π | k | τL (cid:1) , s ( τ ) = sin (cid:0) π | k | τL (cid:1) with | k | = ( k · k ) , (82)and ˆ k := k / | k | if k (cid:54) = 0, otherwise ˆ k := 0. We now consider the Fourier-GEMPIC method defined by the pseudo-differential DFT operatorsfrom Section 3.5.
Directional kinetic steps.
We have the following result.22 emma 4.
For the operators Π (cid:96) defined as in Section 3.5, the exact solution ϕ τ,v α : U → U ( τ ) to the kinetic split step (75) in a direction α ∈ (cid:74) , (cid:75) is given by the explicit expressions ϕ τ,v α : X p,α ( τ ) = X p,α + τ V p,α X p,α + ν ( τ ) = X p,α + ν V p,α ( τ ) = V p,α V p,α + ν ( τ ) = V p,α + ν − ντ q p m p L (cid:88) k ∈ (cid:74) − K,K (cid:75) B α − ν, k ˆ V ,νp,α, k ( τ ) E α, k ( τ ) = E α, k − τ (cid:88) p =1 ··· N q p ˆ V p,α, k ( τ ) E α + ν, k ( τ ) = E α + ν, k B k ( τ ) = B k (83) for all p ∈ (cid:74) , N (cid:75) , k ∈ (cid:74) − K, K (cid:75) and ν = ± , with ˆ V ,νp,α, k ( τ ) := V p,α (Π α + ν S X p ) k if k α = 0 − (cid:16) πk α L (cid:17) − τ (cid:2) (Π α + ν S X p ) k (cid:3) τ if k α (cid:54) = 0 (84) and ˆ V p,α, k ( τ ) := V p,α (Π S X p ) k if k α = 0 − (cid:16) πk α L (cid:17) − τ (cid:2) (Π S X p ) k (cid:3) τ if k α (cid:54) = 0 . (85) Proof.
We remind the velocity equation from (75),dd t V p,α + ν = − ν q p m p L (cid:88) k B α − ν, k V p,α (Π α − ν S X p ) k and the electric one, dd t E α, k = − (cid:88) p q p V p,α (Π α S X p ) k . Here by definition of the pseudo-differential projection operators, we have(Π α − ν S X p ) k = γ k ˆ D − k ,α ˆ D − k ,α + ν ˜ F M, k ( ˜ D k ,α ˜ D k ,α + ν S X p ) (86)and (Π α S X p ) k = γ k ˆ D − k ,α ˜ F M, k ( ˜ D k ,α S X p ) . (87)We begin with the term involving Π as the computations are simpler. We first observe that if k α = 0, then V p,α (Π α S X p ) k = γ k V p,α ˜ F M, k ( S X p ) = V p,α (Π S X p ) k e α .Next if k α (cid:54) = 0, we compute V p,α (Π α S X p ) k = γ k (cid:16) πk α L (cid:17) − V p,α ˜ F M, k ( ∂ α S X p )= γ k (cid:16) πk α L (cid:17) − V p,α M (cid:88) m ∂ α S ( m h − X p ( t ))e − π k · m M = − γ k (cid:16) πk α L (cid:17) − M (cid:88) m dd t (cid:16) S ( m h − X p ( t )) (cid:17) e − π k · m M = − (cid:16) πk α L (cid:17) − dd t (Π S X p ) k . It follows that we can integrate exactly (cid:82) τ V p,α (Π α S X p ) k d t = τ ˆ V p,α, k ( τ ) with the expressionfrom (85), which provides the update of E α, k . Turning to the terms involving Π we observethat if k α = 0, then V p,α (Π α − ν S X p ) k = γ k V p,α ˆ D − k ,α + ν ˜ F M, k ( ˜ D k ,α + ν S X p ) = V p,α (Π α + ν S X p ) k and again this quantity is constant, for the same reason as above. Now if k α (cid:54) = 0, then using dd t X p ( t ) = V [ α ] p = V p,α e α we compute V p,α (Π α − ν S X p ) k = γ k (cid:16) πk α L (cid:17) − V p,α ˆ D − k ,α + ν ˜ F M, k ( ∂ α ˜ D k ,α + ν S X p )= γ k (cid:16) πk α L (cid:17) − ˆ D − k ,α + ν V p,α M (cid:88) m ∂ α ˜ D k ,α + ν S ( m h − X p ( t ))e − π k · m M = − γ k (cid:16) πk α L (cid:17) − ˆ D − k ,α + ν M (cid:88) m dd t (cid:16) ˜ D k ,α + ν S ( m h − X p ( t )) (cid:17) e − π k · m M = − (cid:16) πk α L (cid:17) − dd t (Π α + ν S X p ) k . It follows that (cid:82) τ V p,α (Π α − ν S X p ) k d t = τ ˆ V ,νp,α, k ( τ ) , now with the expression from (84). Thisgives the update for V p,α + ν . (cid:3) Electromagnetic step.
The exact solution flow ϕ τ,EB is again given by (80)–(82), up toreplacing the Π projected term in (80) by its specific expression, see (87). By applying the same splitting method as above to the momentum and Gauss preserving schemepresented in Section 3.2, we obtain a standard Fourier-PIC coupling as described in Section 2.2,as follows. For the kinetic subsystem (75) along α , the modification from (39) to (40) appliesto α ← α + ν and amounts to replacing the operator Π α − ν by Π α . For the electromagnetic24ubsystem (76) it just amounts to replacing Π α by Π . The modified equations read then dd t X p,α = V p,α dd t X p,α + ν = 0 dd t V p,α = 0 dd t V p,α + ν = − ν q p m p L (cid:88) k V p,α B α − ν, k (Π α S X p ) k dd t E α, k = − (cid:88) p q p V p,α (Π α S X p ) k dd t E α + ν, k = 0 dd t B k = 0 and dd t X p = 0 dd t V p = q p m p L (cid:88) k E k (Π S X p ) k dd t E k = π k L × B k dd t B k = − π k L × E k , and explicit solutions are obtained by computing as in the proof of Lemma 4. For the kineticsubsystem along α , the flow ϕ τ,v α is given by ϕ τ,v α : X p,α ( τ ) = X p,α + τ V p,α X p,α + ν ( τ ) = X p,α + ν V p,α ( τ ) = V p,α V p,α + ν ( τ ) = V p,α + ν − ντ q p m p L (cid:88) k ∈ (cid:74) − K,K (cid:75) B α − ν, k ˆ V p,α, k ( τ ) E α, k ( τ ) = E α, k − τ (cid:88) p q p ˆ V p,α, k ( τ ) E α + ν, k ( τ ) = E α + ν, k B K ( τ ) = B K where we remind that ˆ V p,α, k ( τ ) is defined in (85). For the electromagnetic step the flow ϕ τ,EB corresponds to (80)–(82) with a modified velocity equation, namely V p,α ( τ ) = V p,α + q p m p L (cid:88) k ∈ (cid:74) − K,K (cid:75) (cid:16) (cid:90) τ E α, k ( t ) d t (cid:17) (Π S X p ) k . (88)In particular, we see that the fully discrete steps all involve the operator Π , see (49), whichcorresponds to the standard DFT coupling described in Section 2.2, up to the use of exactintegrals (50) along zero-modes. This scheme does not have a Hamiltonian structure, but itsatisfies some important conservation properties. Lemma 5.
The Fourier-PIC scheme described above preserves the discrete Gauss laws (41) and the discrete momentum (44) .Proof.
The proof is a matter of elementary computations, similar to the ones involved in theproofs of Theorem 1. One key observation is that the above modifications in the velocity equa-tions have no effect on the Gauss laws being preserved, and they precisely lead to an exactmomentum preservation. (cid:3) .7 Summary of the proposed methods In the above sections we have derived geometric Fourier-particle methods following the DiscreteAction Principle formalized in [9], where the coupling between the fields and the particles isrepresented by abstract projection operators on the truncated Fourier spaces that satisfy acommuting diagram property. This results in semi-discrete schemes that preserve the Gausslaws and have a discrete Hamiltonian structure relying on a non-canonical Poisson bracket.When the coupling operators are defined as L projections, the coupling is gridless andessentially relies on continuous Fourier coefficients. The resulting method preserves the totalmomentum in addition to the charge and energy, and corresponds to the Particle-In-Fourier(PIF) approach. Here it is called GEMPIF to emphasize its geometric nature.When the coupling operators are defined as pseudo-differential DFT projections, the cou-pling involves a grid and essentially relies on discrete Fourier coefficients. The resulting methodcan be seen as a variant of standard spectral PIC methods and is called Fourier-GEMPIC.Fully discrete schemes of various orders in time have then been constructed by a Hamiltoniansplitting procedure. These schemes are Poisson maps, in particular they preserve the discretePoisson bracket and a modified energy that approximates the exact one to the time order ofthe splitting.Finally we have observed that Fourier-PIC schemes with standard coupling terms can be ob-tained as a variant of the above Fourier-GEMPIC method. These schemes are not Hamiltonian,but they preserve exactly the Gauss laws and the total momentum.Since our approach handles general shape functions S and arbitrary filter coefficients γ k ,we may use high order shapes with anti-aliasing properties in the Fourier-GEMPIC methods,and the ad-hoc back-filtering procedure (28)–(30) to avoid the damping of relevant modes inthe computational range. In this section we present some numerical results obtained with several of the methods describedabove. Specifically, we will compare • the GEMPIF scheme described in Section 4.4 with point shape functions ( S = δ ), • the Fourier-GEMPIC scheme described in Section 4.5, involving a DFT grid with M pointsand B-spline shapes of degree κ . To assess the benefits of back-filtering, we will use twoversions of this scheme: a plain smoothed Fourier SF-GEMPIC method corresponding to γ k = 1, and a back-filtered Fourier BFF-GEMPIC method that involves filter coefficients γ k := 1 σ k = (cid:89) α =1 (cid:18) sinc (cid:16) πk α µ (2 K + 1) (cid:17)(cid:19) − ( κ +1) with µ := M K + 1 ≥ µ is the oversampling parameter, see (24). • the Fourier-PIC scheme described in Section 4.6, which also uses a DFT grid with M points and B-spline shapes of degree κ . Similarly as above, this scheme will be used intwo versions, namely a basic smoothed Fourier SF-PIC method corresponding to γ k = 1,and a back-filtered Fourier BFF-PIC method involving the coefficients (89).26 .1 Periodic plasma test-cases
For our numerical experiments we consider a periodic one-species model in 1D2V and 1D1Vsimilarly as in [9], with zero-mean current to preserve the total momentum, i.e. ∂ t E = − J + 1 L (cid:90) L J , ∂ t E + ∂ x B = − J + 1 L (cid:90) L J , ∂ t B + ∂ x E = 0 , and standard plasma test-cases: classical Landau damping test-cases corresponding to f ( x, v ) = (1 + (cid:15) cos( k x )) 1 √ π exp (cid:16) − v (cid:17) with k = 0 . (cid:15) = 0 . .
01 for the strong and weak damping, a standard
Weibelinstability [39, 23] where f ( x, v , v ) = 12 πσ σ exp (cid:16) − (cid:16) v σ + v σ (cid:17)(cid:17) and where the initial magnetic field is B ( t = 0 , x ) = β cos( k x ), with σ = 0 . / √ σ = √ σ , β = 10 − and k = 1 .
25, a two-stream instability f ( x, v ) = (1 + (cid:15) cos( k x )) 12 √ π (cid:16) exp (cid:0) − ( v + 2 . (cid:1) − exp (cid:0) − ( v − . (cid:1)(cid:17) with (cid:15) = 5 · − and k = 0 .
2, and a bump-on-tail instability corresponding to f ( x, v ) = (1 + (cid:15) cos( k x )) 1 √ π (cid:16) δ exp (cid:0) − v (cid:1) + 2(1 − δ ) exp (cid:0) − v − . (cid:1)(cid:17) with δ = 0 . (cid:15) = 5 · − and k = 0 .
3. For these test-cases the domain size is L = π k , exceptfor the last one where we take L = π k . The initial B field is zero unless stated otherwise,and the initial E field is computed from the particles by solving the periodic Poisson equation.In practice we run all these cases using a 1D2V implementation of the above methods. Whenthe initial distribution is only 1D1V we sample the particles following a Maxwellian in v withthermal velocity σ = 1. In order to compare the qualitative response of the different schemes, we first plot the relevantenergy curves for the above test-cases using rather coarse numerical parameters such as K , M and κ , but a number of particles N high enough to match the reference damping or growthrates.In Figure 1 and 2 we plot the energy curves for the above test cases, using K = 3 Fouriermodes, splines of degree κ = 3 and DFT grids with M = 8 (left) or M = 16 points (right). Forthe strong (top) and weak (bottom) Landau damping test-cases in Figure 1 we use N = 5 · and N = 10 particles, respectively. For the Weibel (top), two-stream (center) and bump-on-tail instabilities in Figure 2 we use N = 10 , N = 5 · and N = 10 particles, respectively.The gray curves show reference runs, obtained using the GEMPIF scheme with 5 · particlesand K = 15 Fourier modes. In all these runs, the time scheme is given by a fourth-order(Suzuki-Yoshida) composition method (74) with a time step of ∆ t = 0 . µ = 16 / ≈ κ = 5 for some of the above test cases. For the two filteredmethods, this has no visible effect: they are still on top of the gridless method (which involvesno splines). For the two unfiltered methods this actually degrades the results, as the energycurves are more distant from the reference one than with κ = 3.28 a) strong damping with an M = 8 grid (b) strong damping with an M = 16 grid(c) weak damping with an M = 8 grid (d) weak damping with an M = 16 grid Figure 1: Landau damping test-cases: E field energy curves, using K = 3 Fourier modes, N = 5 · particles for the strong damping (top) and N = 10 for the weak damping (bottom).Gridded (GEMPIC and PIC) methods use B-spline shapes of degree κ = 3 and DFT grids with M = 8 (left) or M = 16 points (right). 29 a) Weibel instability with an M = 8 grid (b) Weibel instability with an M = 16 grid(c) two stream instability with an M = 8 grid (d) two stream instability with an M = 16 grid(e) bump-on-tail instability with an M = 8 grid (f) bump-on-tail instability with an M = 16 grid Figure 2: Energy behavior for instability test-cases: B field energy for the Weibel instability(top) using N = 10 particles, and E field energy for the two-stream (center) and bump-on-tailinstabilities (bottom), using N = 5 · and N = 10 particles, respectively. All runs use K = 3Fourier modes, and gridded (GEMPIC and PIC) methods use B-spline shapes of degree κ = 3and DFT grids with M = 8 (left) or M = 16 points (right).30 a) strong damping with an M = 8 grid (b) strong damping with an M = 16 grid(c) weak damping with an M = 8 grid (d) weak damping with an M = 16 grid(e) bump-on-tail instability with an M = 8 grid (f) bump-on-tail instability with an M = 16 grid Figure 3: Landau damping and bump-on-tail instability test-cases, using the same parametersas in Figures 1 and 2 and B-spline shapes of degree κ = 5.31 .3 Long-time conservation properties In order to assess the long-time conservation properties of the different methods, we show inFigure 4 several errors for the Weibel instability test-case on a time range ten times longer thanabove. The plotted errors are in energy conservation,Err n H = |H n − H |H with H n = 12 N (cid:88) p =1 m p | V np | + 12 (cid:90) [0 ,L ] (cid:16) | E nK ( x ) | + | B nK ( x ) | (cid:17) d x (90)momentum conservation,Err n P = (cid:107)P n − P (cid:107) (cid:107)P (cid:107) with P n = N (cid:88) p =1 m p V np + (cid:90) [0 ,L ] E nK ( x ) × B nK ( x ) d x (91)and Gauss law Err nG = (cid:90) [0 ,L ] (cid:12)(cid:12)(cid:12)(cid:16) div E nK − N (cid:88) p =1 q p Π S X np (cid:17) ( x ) (cid:12)(cid:12)(cid:12) d x. (92)For each method we use the same numerical parameters as in Figure 2, and low-resolution runswith spline shapes of degree κ = 3 and κ = 7.For the energy conservation we observe overall a good behavior for all the methods. Herethe main observation is that a clear separation is visible between the three geometric methodsand the two non-geometric ones. While the former show a very good stability of the energyconservation over these long-time ranges, the latter sometimes display sensible growth in theenergy error, characteristic of a less stable behavior. Another interesting difference is that whenthe DFT grid is refined (right plots), the energy errors made by the non-geometric methods getreduced, while those of the geometric methods hardly change. This behavior is consistent withthe backward error analysis which states that Hamiltonian splitting time integrators of order q preserve exactly a modified energy which approximates the exact one with the same order.As for the momentum conservation, the first observation is that the curves confirm theexact conservation of the GEMPIF method and the two non-geometric schemes (SF-PIC andBFF-PIC). For the Fourier-GEMPIC methods the conservation is only approximate, howeverwe observe a very good stability over these long-time simulations. We also observe that theaccuracy to which the momentum is preserved is improved by refining grid or increasing thespline order.Finally the Gauss error curves confirm the exact charge conserving property of all methods.32 a) Energy errors with shapes of degree κ = 3 (b) Energy errors with shapes of degree κ = 7(c) Momentum errors with shapes of degree κ = 3 (d) Momentum errors with shapes of degree κ = 7(e) Gauss law errors with shapes of degree κ = 3 (f) Gauss law errors with shapes of degree κ = 7 Figure 4: Long-time conservation properties. Energy conservation errors (top), momentumconservation errors (middle) and Gauss law errors (bottom) are shown for Weibel instabilityruns using N = 10 particles and K = 3 Fourier modes. Gridded (Fourier-Gempic and Fourier-PIC) methods use a DFT grid with M = 16 points corresponding to an oversampling parameterof µ = M K +1 ≈ . κ = 3 (left) or κ = 7 (right). Similar curveshave been observed for higher resolution runs with analogous oversampling parameter, such as K = 13 and M = 64. Compared to the runs in Fig. 2, the simulation range is ten times longer.33 .4 Convergence studies We finally study how the conservation of energy and momentum is improved by refining thenumerical parameters.In Figure 5 and 6 we plot the time-averaged energy errors N t (cid:80) N t − n =0 Err n H , for the fivemethods (using the same color key as above) as a function of the time step, for various DFTgrids and spline degrees. Results shown in Figure 5 are obtained with the second order Strangscheme (73), while those in Figure 6 correspond to the fourth-order Suzuki-Yoshida scheme(74). For an easier comparison, both figures use the same scale.For the non-geometric methods we observe an improvement in the energy errors when thetime step decreases, but the convergence is limited by the resolution of the grid and the degreeof the splines. This artifact is not present with the geometric methods, where the convergenceof the energy errors holds almost independently of the grid resolution and spline degree. Thisconfirms the behavior already observed in the long-time runs, see Figure 4. It is also in strongagreement with the backward error analysis, which predicts a convergence of the energy errorsof the same order as the time scheme. Here the runs correspond to the Weibel instability with amoderate time range ( T = 500), but results obtained with other test cases have showed similarbehavior.In Figure 7 we plot the time-averaged momentum errors N t (cid:80) N t − n =0 Err n P , for the five methods(again with the same color key), as a function of the ratio 2 K/M ≈ /µ , for various splinedegrees. For the GEMPIF and Fourier-PIC methods the error is close to machine accuracy,since the methods are exactly momentum conserving. For the Fourier-GEMPIC methods weobserve that in every case they converge to 0, with a rate close to O ( K/M ) κ +1 . Here the test-case is again the Weibel instability with K = 3 modes as in Section 5.2, but the same convergencebehavior was observed with other test-cases and a higher number of Fourier modes. This work has been carried out within the framework of the EUROfusion Consortium and hasreceived funding from the Euratom research and training programme 2014-2018 and 2019-2020under grant agreement No 633053. The views and opinions expressed herein do not necessarilyreflect those of the European Commission. 34 a) M = 8 and κ = 3 (b) M = 16 and κ = 3(c) M = 8 and κ = 5 (d) M = 16 and κ = 5 Figure 5: Energy error convergence: time-averaged energy conservation errors as a function ofthe time step ∆ t , for a Weibel instability test-case with K = 3 Fourier modes, N = 5 · particles and various grids and spline degrees, as indicated. Results reported here use a second-order (Strang) time scheme. 35 a) M = 8 and κ = 3 (b) M = 16 and κ = 3(c) M = 8 and κ = 5 (d) M = 16 and κ = 5 Figure 6: Energy error convergence curves, using the same parameters as in Figure 5, and afourth-order (Suzuki-Yoshida) time scheme. 36 a) Spline degree κ = 3 (b) Spline degree κ = 5(c) Spline degree κ = 9 (d) Spline degree κ = 11 Figure 7: Momentum error convergence: time-averaged momentum conservation errors as afunction of the ratio 2
K/M , using the Weibel instability test-case with K = 3 Fourier modes, N = 5 · particles and spline degrees as indicated. Results reported here use a fourth-order(Suzuki-Yoshida) time scheme. 37 eferences [1] Jakob Ameres. Stochastic and Spectral Particle Methods for Plasma Physics . Dissertation,Technische Universit¨at M¨unchen, M¨unchen, 2018.[2] Douglas N. Arnold, Richard S. Falk, and Ragnar Winther. Finite Element Exterior Cal-culus: from Hodge theory to numerical stability.
Bulletin of the American MathematicalSociety , 47:281–354, 2010. doi: 10.1090/S0273-0979-10-01278-4.[3] C K Birdsall and A.B. Langdon.
Plasma physics via computer simulation . Adam Hilger,IOP Publishing, 1991.[4] Charles K. Birdsall and Neil Maron. Plasma self-heating and saturation due to numericalinstabilities.
J. Comput. Physics , 36(1):1–19, June 1980.[5] R B Blackman and J W Tukey.
The Measurement of Power Spectra . Dover, 1958.[6] Alain Bossavit.
Computational electromagnetism: variational formulations, complemen-tarity, edge elements . Academic Press, 1998.[7] Annalisa Buffa, Giancarlo Sangalli, and Rafael V´azquez. Isogeometric analysis in elec-tromagnetics: B-splines approximation.
Computer Methods in Applied Mechanics andEngineering , 199(17):1143–1152, 2010. doi: 10.1016/j.cma.2009.12.002.[8] Martin Campos Pinto and Eric Sonnendr¨ucker. Compatible Maxwell solvers with particlesI: conforming and non-conforming 2d schemes with a strong Ampere law.
The SMAIjournal of computational mathematics , 3:53–89, 2017. doi: 10.5802/smai-jcm.20. URL smai-jcm.centre-mersenne.org/item/SMAI-JCM_2017__3__53_0/ .[9] Martin Campos Pinto, Katharina Kormann, and Eric Sonnendr¨ucker. Variational frame-work for structure-preserving electromagnetic Particle-In-Cell methods. arXiv preprintarXiv:2101.09247 , 2021.[10] Nicolas Crouseilles, Lukas Einkemmer, and Erwan Faou. Hamiltonian splitting for theVlasov–Maxwell equations.
Journal of Computational Physics , 283:224–240, 2015. doi:10.1016/j.jcp.2014.11.029.[11] Viktor K Decyk. Description of Spectral Particle-in-Cell Codes from the UPIC Framework.
Presentation at ISSS-10 , 2011. URL https://picksc.idre.ucla.edu/wp-content/uploads/2015/05/UPICModels.pdf .[12] T Zh Esirkepov. Exact charge conservation scheme for particle-in-cell simulation with anarbitrary form-factor.
Computer Physics Communications , 135(2):144–153, 2001.[13] Evstati G. Evstatiev and Bradley A. Shadwick. Variational formulation of particle algo-rithms for kinetic plasma simulations.
Journal of Computational Physics , 245:376–398,2013. doi: 10.1016/j.jcp.2013.03.006.[14] Claude Gasquet and Patrick Witomski.
Fourier Analysis and Applications , volume 30 of
Texts in Applied Mathematics . Springer New York, New York, NY, 1999.[15] Brendan B Godfrey. Numerical Cherenkov instabilities in electromagnetic particle codes.
Journal of Computational Physics , 15(4):504–521, August 1974.3816] Brendan B Godfrey, Jean-Luc Vay, and Irving Haber. Numerical stability analysis of thepseudo-spectral analytical time-domain PIC algorithm.
Journal of Computational Physics ,258:689–704, February 2014.[17] Ernst Hairer, Christian Lubich, and Gerhard Wanner.
Geometric Numerical Integration .Springer, 2006.[18] R Hatzky, R Kleiber, A K¨onies, A Mishchenko, M Borchardt, A Bottino, and E. Son-nendr¨ucker. Reduction of the statistical error in electromagnetic gyrokinetic particle-in-cell simulations.
Journal of Plasma Physics , 85(1):905850112, 2019. doi: 10.1017/s0022377819000096.[19] Yang He, Hong Qin, Yajuan Sun, Jianyuan Xiao, Ruili Zhang, and Jian Liu. Hamiltonianintegration methods for Vlasov–Maxwell equations.
Physics of Plasmas , 22:124503, 2015.doi: 10.1063/1.4938034.[20] Ralf Hiptmair. Finite elements in computational electromagnetism.
Acta Numerica , 11:237–339, 2002.[21] R.W. Hockney and J.W. Eastwood.
Computer simulation using particles . Taylor & Francis,Inc, Bristol, PA, USA, 1988.[22] C-K Huang, Yong Zeng, Ying Wang, Michael D Meyers, Sunghwan Yi, and Brian J Al-bright. Finite grid instability and spectral fidelity of the electrostatic particle-in-cell algo-rithm.
Computer Physics Communications , 207:123–135, 2016.[23] Michael Kraus, Katharina Kormann, Philip J. Morrison, and Eric Sonnendr¨ucker. GEM-PIC: Geometric electromagnetic particle-in-cell methods.
Journal of Plasma Physics , 83(4), 2017.[24] A Bruce Langdon. Effects of the spatial grid in simulation plasmas.
Journal of Computa-tional Physics , 6(2):247–267, October 1970.[25] A Bruce Langdon. Kinetic theory for fluctuations and noise in computer simulation ofplasma.
Physics of Fluids , 22(1):163–10, 1979.[26] A Bruce Langdon and Charles K. Birdsall. Theory of Plasma Simulation Using Finite-SizeParticles.
Physics of Fluids , 13(8):2115–2122, August 1970.[27] Robert I. McLachlan and G. Reinout W. Quispel. Splitting methods.
Acta Numerica , 11:341–434, 2002. doi: 10.1017/S0962492902000053.[28] B F McMillan, S Jolliet, A Bottino, P Angelino, T M Tran, and L Villard. Rapid Fourierspace solution of linear partial integro-differential equations in toroidal magnetic confine-ment geometries.
Computer Physics Communications , 181(4):715 – 719, 04 2010. doi:10.1016/j.cpc.2009.12.001.[29] Matthew S Mitchell, Matthew T Miecnikowski, Gregory Beylkin, and Scott E Parker.Efficient Fourier Basis Particle Simulation.
Journal of Computational Physics , pages 837–847, August 2019. 3930] N Ohana, A Jocksch, E Lanti, TM Tran, S Brunner, C Gheller, F Hariri, and L Villard.Towards the optimization of a gyrokinetic Particle-In-Cell (PIC) code on large-scale hybridarchitectures. In
Journal of Physics: Conference Series , volume 775, page 012010. IOPPublishing, 2016.[31] Okuda, Hideo. Nonphysical noises and instabilities in plasma simulation due to a spatialgrid.
J. Comput. Physics , 10(3):475–486, December 1972.[32] Gerlind Plonka, Daniel Potts, Gabriele Steidl, and Manfred Tasche. Fast Fourier Trans-forms for Nonequispaced Data. In
Numerical Fourier Analysis , pages 377–419. Birkh¨auser,Cham, Cham, 2018.[33] Bradley A. Shadwick, Alexander B. Stamm, and Evstati G. Evstatiev. Variational formu-lation of macro-particle plasma simulation algorithms.
Physics of Plasmas , 21(5):055708,2014. doi: 10.1063/1.4874338.[34] Gabriele Steidl. A note on fast Fourier transforms for nonequispaced grids.
Advances inComputational Mathematics , 9(3-4):337–352, 1998.[35] Jean-Luc Vay, Irving Haber, and Brendan B Godfrey. A domain decomposition method forpseudo-spectral electromagnetic simulations of plasmas.
Journal of Computational Physics ,243:260–268, June 2013.[36] Jr H D Victory and Edward J Allen. The Convergence Theory of Particle-In-Cell Methodsfor Multidimensional Vlasov–Poisson Systems.
SIAM Journal on Numerical Analysis , 28(5):1207 – 1241, 1991. doi: 10.1137/0728065. URL http://epubs.siam.org/doi/abs/10.1137/0728065 .[37] G Vlad, S Briguglio, G Fogaccia, and B Di Martino. Gridless finite-size-particle plasmasimulation.
Computer physics communications , 134(1):58–77, 2001.[38] Stephen D Webb. A spectral canonical electrostatic algorithm.
Plasma Physics and Con-trolled Fusion , 58(3):034007, 2016.[39] Erich S. Weibel. Spontaneously growing transverse waves in a plasma due to an anisotropicvelocity distribution.
Physical Review Letters , 2:83–84, 1959. doi: 10.1103/PhysRevLett.2.83.[40] Xinlu Xu, Peicheng Yu, Samual F Martins, Frank S Tsung, Viktor K Decyk, Jorge Vieira,Ricardo A Fonseca, Wei Lu, Luis O Silva, and Warren B Mori. Numerical instability dueto relativistic plasma drift in EM-PIC simulations.