GEMPIC: Geometric ElectroMagnetic Particle-In-Cell Methods

1 Introduction

• The Vlasov–Maxwell system couples kinetic particle dynamics with Maxwell’s equations and possesses variational, Hamiltonian, and Casimir structures that discretizations should preserve.These structures imply conserved quantities associated with symmetries and Poisson-bracket degeneracy.Source: §1 · p. 3
• Maxwell discretizations must preserve the deRham-complex compatibility conditions, div B = 0, and Gauss’ law, div E = ρ, to obtain stable methods.The relevant vector-calculus identities are curl grad = 0 and div curl = 0.Source: §1 · p. 3
• Earlier geometric PIC approaches preserved selected energy, charge, or gauge properties, but existing bracket discretizations did not establish all requirements for a Hamiltonian semi-discrete system.In particular, prior work either lacked the Jacobi-identity proof or omitted Casimir conservation.Source: §1 · p. 4
• The authors unify preceding ideas in a flexible FEEC framework for structure-preserving Vlasov–Maxwell discretization.The framework uses compatible Finite Element spaces rather than requiring one particular basis.Source: §1 · pp. 4–5
• The semi-discrete noncanonical Poisson structure preserves anti-symmetry, the Jacobi identity, and Casimir invariants, thereby remaining Hamiltonian.These guarantees hold for any compatible Finite Element basis, including the spline spaces discussed later.Source: §1 · p. 5
• Poisson time integrators based on Hamiltonian splitting and higher-order compositions extend these conservation properties to the fully discrete scheme.The splitting approach provides explicit time stepping and supports higher-order methods.Source: §1 · p. 5
• Preserving Hamiltonian structure targets the full hierarchy of Poincaré integral invariants, symmetries, local conservation laws, and Casimir invariants, not merely global energy or momentum.Global energy conservation alone may leave large finite-dimensional systems with physically incorrect solutions.Source: §1 · p. 5
• The paper proceeds from bracket theory and FEEC foundations to discretization, splitting, long-time analysis, and numerical tests of nonlinear Landau damping and the Weibel instability.The discrete analysis includes Jacobi identity, Casimir invariants, and discrete Gauss’ law, while experiments use one-dimensional-two-velocity splines.Source: §1 · p. 5

2 The Vlasov–Maxwell System

• The Vlasov–Maxwell system couples species distribution functions to electric and magnetic fields through charge and current sources, subject to initial and boundary conditions.Numerical approximations must maintain compatibility between discrete sources and divergence and curl operators to preserve Maxwell constraints.Source: §2 · p. 6

2.1 Noncanonical Hamiltonian Structure

• Its noncanonical Poisson bracket is antisymmetric, satisfies Leibniz’ rule and the Jacobi identity, and generates the Vlasov and Maxwell equations from the particle-field energy Hamiltonian.Preserving this structure discretely requires compatible finite-element spaces, provided naturally by Finite Element Exterior Calculus, to maintain continuity, Gauss’ law, and div B = 0.Source: §2.1 · pp. 7–8

2.2 Invariants

• The system conserves Casimirs associated with phase-space volume, Gauss’ law, and magnetic divergence, while momentum maps are conserved when they generate Hamiltonian-preserving symmetries.The magnetic-divergence functional is a pseudo-Casimir whose bracket properties require div B = 0; momentum conservation additionally depends on satisfying Poisson’s equation.Source: §2.2 · pp. 8–9

3 Finite Element Exterior Calculus

• FEEC systematically analyzes stability and convergence of mixed finite-element discretizations using geometrical and topological structure, with Maxwell’s equations as a central example.Source: §3 · p. 10

3.1 Maxwell’s Equations

• Interpreting Maxwell fields as differential forms determines a natural finite-element discretization: primal equations are treated strongly, while Ampère’s equation and Gauss’ law use weak dual formulations.E and H are 1-forms, D and B are 2-forms, while charge and current are represented in the finite-element setting as 0- and 1-forms.Source: §3.1 · p. 10

3.2 Finite Element Spaces of Differential Forms

• Compatible spaces V0, V1, V2, V3 reproduce the differential-form complex and commuting projections, producing stable discretizations across choices including Nédélec, Raviart–Thomas, mimetic, and spline elements.Spline finite elements are selected for their smoothness properties, which can reduce noise in PIC codes.Source: §3.2 · pp. 11–12

3.3 Finite Element Discretization of Maxwell’s Equations

• The compatible discretization represents curl, gradient, and divergence through discrete matrices and mass matrices, while satisfying div B = 0 strongly, a property important for the Poisson bracket’s Jacobi identity.Source: §3.3 · pp. 12, 13, 15

3.4 Example: B-Spline Finite Elements

• Tensor-product B-spline bases realize the differential-form spaces, naturally represent gradient, curl, and divergence actions, and yield a semi-discrete magnetic field satisfying div Bh = 0.The construction assumes periodic boundary conditions when rearranging the spline sums.Source: §3.4 · pp. 16–18

4 Discretization of the Hamiltonian Structure

• The discretization replaces the continuum distribution with Np particle characteristics and restricts the electromagnetic fields to finite-element subspaces V1 and V2, producing variables X, V, e, and b.The FEEC spaces assign φh and ρh to V0, Ah, Eh, and Jh to V1, and Bh to V2.Source: §4 · p. 18

4.1 Discretization of the Functional Field Derivatives

• Functionals of the finite-element fields become functions of their coefficient vectors, with functional derivatives represented through the finite-element mass matrices and dual bases.The same procedure is applied to the remaining field-derivative terms before discretizing the distribution function.Source: §4.1 · pp. 18, 20

4.2 Discretization of the Functional Particle Derivatives

• Assuming particles with masses, signed charges, weights, positions, and velocities, functionals of the distribution become functions F(X,V), yielding the discrete kinetic bracket.The kinetic contribution is obtained by replacing the species distribution with its particle approximation and inserting the resulting particle derivatives.Source: §4.2 · pp. 20–21

4.3 Discrete Poisson Bracket

• The discrete bracket is expressed in matrix form using the curl matrix, particle evaluation matrices, and diagonal mass and charge matrices for the coupled particle–field dynamics.Its Poisson matrix depends on particle positions and magnetic degrees of freedom and is explicitly anti-symmetric.Source: §4.3 · pp. 21–22

4.3 Discrete Poisson Bracket / 4.4 Jacobi Identity

• The semi-discrete Poisson bracket is anti-symmetric and satisfies the Jacobi identity when the discrete magnetic field is divergence-free and the finite-element spaces form a de Rham complex.The magnetic field must be closed but need not be exact, while the curl-basis compatibility condition completes the Jacobi verification.Source: §4.3 · p. 22, §4.4 · pp. 24–25

4.5 Discrete Hamiltonian and Equations of Motion

• The Hamiltonian is discretized in the particle and finite-element variables, and Hamilton’s equations evolve particle positions, velocities, electric-field coefficients, and magnetic-field coefficients.The mass and charge matrices provide compact notation for the Poisson structure and need not be explicitly constructed or inverted.Source: §4.5 · pp. 25–26

4.6 Discrete Gauss’ Law

• The discrete equations conserve Gauss’ law and include a discrete continuity equation, although conservation of the relation in the fully discrete system depends on the time-stepping scheme.The derivation uses the complex identity CG = 0 and defines discrete charge and current densities from the particles.Source: §4.6 · p. 26

4.7 Discrete Casimir Invariants

• The semi-discrete Casimirs conserve the discrete magnetic divergence and Gauss’ law, and remain valid under the weaker complex property rather than requiring an exact sequence.An exact sequence is assumed to identify the invariants, but the final verification requires only Im D⊤ ⊆ Ker C⊤ and grad Λ0(X) = Λ1(X)G.Source: §4.7 · p. 27

5 Hamiltonian Splitting

• The discrete equations split into particle, electric-field, and magnetic-field subsystems, and composing their exact Poisson maps yields Poisson structure-preserving integrators.Each subsystem flow is a Poisson map, and compositions of Poisson maps remain Poisson.Source: §5 · p. 28

5.1 Solution of the Subsystems

• Exact subsystem integration uses electric- and magnetic-field flows plus further particle sub-splitting when general magnetic fields prevent direct integration.Only one particle-position component changes per substep, making the required straight-line trajectory integrals tractable; costs are comparable to charge-conserving methods.Source: §5.1 · pp. 29–30

5.2 Splitting Methods

• Lie–Trotter provides first-order integration, symmetric Strang provides second-order integration, and tuned or composed schemes extend the method to smaller-error and fourth-order integrators.The free parameter α can reduce the error constant, with α = 0.1932 giving particularly small error.Source: §5.2 · p. 31

5.3 Backward Error Analysis

• The Lie–Trotter integrator preserves the modified energy H̃ = H + ∆t H̃1 to O(∆t^2), whereas the original energy H is preserved only to O(∆t).The modified Hamiltonian is derived through Baker–Campbell–Hausdorff expansions of the subsystem Hamiltonian vector fields.Source: §5.3 · pp. 32–33

6 Example: Vlasov–Maxwell in 1D2V

• Under the 1D2V assumptions, the fields and distribution depend on x, v1, and v2, with div B = 0 manifest.The reduction sets x = (x, 0, 0), v = (v1, v2, 0), E = (E1, E2, 0), and B = (0, 0, B3).Source: §6 · p. 33

6.1 Noncanonical Hamiltonian Structure

• The reduced noncanonical bracket and Hamiltonian generate the reduced Vlasov–Maxwell equations, preserving the Hamiltonian formulation.The reduced bracket with the reduced Hamiltonian generates equations (6.4), (6.5), and (6.6).Source: §6.1 · p. 34

6.2 Reduced Jacobi Identity

• The reduced bracket satisfies the Jacobi identity because reduced functionals are closed under the full bracket.Closure means the bracket of two reduced functionals is again reduced; reductions omitting magnetic-field dependence generally lack closure and Jacobi identity.Source: §6.2 · pp. 34–35

6.3 Discrete deRham Complex in one Dimension

• The one-dimensional discrete de Rham complex assigns E1 and B3 to one-form spaces and E2 to a zero-form space.The field solver represents E_h with degree-p splines and D_h, B_h with degree-(p−1) splines.Source: §6.3 · p. 35

6.4 Discrete Poisson Bracket

• The reduced system has a discrete Poisson bracket that also satisfies the Jacobi identity.The identity follows either by direct proof, analogous to the full bracket, or by the reduced-functional closure argument.Source: §6.4 · pp. 35–36

6.5 Discrete Hamiltonian and Equations of Motion

• The discrete Hamiltonian produces equations of motion in direct correspondence with the reduced continuous equations.The reduced degrees of freedom are u = (X, V1, V2, d, e, b).Source: §6.5 · p. 36

6.6 Hamiltonian Splitting

• Hamiltonian splitting solves the discrete dynamics through separate subflows for HD + HE, HB, Hp1, and Hp2.The HD + HE subflow updates particle momenta using the discrete electric-field couplings, while the other subflows are solved separately.Source: §6.6 · pp. 36–37

7 Numerical Experiments

• The experiments compare Hamiltonian splitting with Boris–Yee for 1d2v electron dynamics, using degree-three 0-form splines and Sobol particle loading with antithetic sampling.The comparisons target conservation of total energy and Casimirs across multiple test cases.Source: §7 · p. 37

7.1 Weibel Instability

• 0.02784 growth rate is verified for the Weibel instability, while Hamiltonian splitting preserves Gauss’ law each step and higher orders reduce energy error.The simulation used 100,000 particles, 32 grid points, degree-3 and degree-2 splines, and ∆t = 0.05.Source: §7.1 · p. 38

7.2 Streaming Weibel Instability

• The streaming Weibel test verifies the 0.03 growth rate and again shows exact Gauss’ law conservation with Hamiltonian splitting plus similar energy-error behavior across schemes.The simulation used 20,000,000 particles, 128 grid points, degree-3 and degree-2 splines, and ∆t = 0.01; earlier E1 growth is attributed to particle noise.Source: §7.2 · p. 40

7.3 Strong Landau Damping

• For strong Landau damping, fitted damping and growth rates are in good agreement with other codes, while fourth-order methods show excellent energy conservation over time.Source: §7.3 · p. 42

7.4 Backward Error Analysis

• Backward error analysis confirms the theoretical convergence rates for total and corrected energy errors across time steps ∆t = 0.01, 0.02, 0.05.The Lie–Trotter Hamiltonian error is order ∆t, while the first-order corrected modified Hamiltonian error is order (∆t)2.Source: §7.4 · pp. 42–43

7.5 Momentum Conservation

• 1.6×10−16 maximum deviation demonstrates that Strang splitting conserves the tested momentum relation to very high accuracy, although momentum is not conserved exactly.The reported simulation used 20,000,000 particles, 128 grid points, and ∆t = 0.01; momentum error remains small during the linear phase.Source: §7.5 · pp. 43, 47

• GEMPIC preserves the Poisson structure, Casimir invariants, Gauss’ law, and div B = 0 through geometric spatial discretization and Hamiltonian splitting.The method is a rare genuine Poisson integrator; energy is not exact, but backward error analysis shows its error is independent of degrees of freedom, particle count, and time steps.Source: §8 · p. 47
• Numerical experiments in 1d2v with spline electromagnetic-field bases verified the method’s favourable conservation and stability properties.These conservation laws arise from preserving the equations’ underlying geometric structure rather than being manually imposed.Source: §8 · p. 47
• The implementation has cost comparable to charge-conserving methods, uses mixed-degree field bases and exact line integrals, and extends conceptually to one, two, and three dimensions.It is somewhat more expensive than non-charge-conserving methods, which can suffer spurious instabilities; a three-dimensional implementation and method comparison were deferred to later work.Source: §8 · p. 47
• Future work includes deRham-compatible finite element spaces, gyrokinetic and relativistic Vlasov–Maxwell systems, alternative energy-preserving time stepping, and collisions.Non-Hamiltonian effects such as collisions may be represented microscopically by stochastic Hamiltonian processes or macroscopically by metriplectic brackets.Source: §8 · pp. 47–48

Discretization scheme

• The alternative Boris–Yee scheme is formulated with conforming Finite Elements and staggered variables evaluated at integer and half-integer time levels.The staggering uses X, d, and e at t_n+1/2, while V and b are evaluated at t_n.Source: A The Boris–Yee scheme · p. 48

Time stepping

• The Vlasov–Maxwell system is propagated from the previous staggered state through a sequence of field and particle updates over each time step.The listed updates include computing b^n and then d^n+1/2 according to the scheme’s equations.Source: A The Boris–Yee scheme · pp. 48–49

Initialization

• Initialization samples X^0 and V^0, sets the initial fields, solves Poisson’s equation for d^0, and computes the first half-step variables using a half time step.The first half-step substitutes b^0 and V^0 for values at time Δt/4.Source: A The Boris–Yee scheme · p. 49
• Although the startup half-step introduces an error of order (Δt)^2, occurring only once, the overall Boris–Yee scheme remains second-order accurate.The stated accuracy concerns the complete scheme after the one-time initialization error.Source: A The Boris–Yee scheme · p. 49