Structure and structure-preserving algorithms for plasma physics

Hamiltonian and action principle (HAP) formulations of plasma physics are reviewed for the purpose of explaining structure preserving numerical algorithms. Geometric structures associated with and emergent from HAP formulations are discussed. These include conservative integration, which exactly conserves invariants, symplectic integration, which exactly preserves the Hamiltonian geometric structure, and other Hamiltonian integration techniques. Basic ideas of variational integration and Poisson integration, which can preserve noncanonical Hamiltonian structure, are discussed. Metriplectic integration, which preserves the structure of conservative systems with both Hamiltonian and dissipative parts, is proposed. Two kinds of simulated annealing, a relaxation technique for obtaining equilibrium states, are reviewed: one that uses metriplectic dynamics, which maximizes an entropy at fixed energy, the other that uses double bracket dynamics, which preserves Casimir invariants. Throughout, applications to plasma systems are emphasized. The paper culminates with a discussion of GEMPIC [Kraus et al., arXiv:1609.03053v1 [math.NA] (2016)], a particle in cell code that incorporates Hamiltonian and geometrical structure preserving properties.


I. INTRODUCTION
The purpose of this article is to give an overview of computational algorithms designed to preserve structure of dynamical systems, and to review some recent progress in the development of structure preserving algorithms for plasma simulation.
The first question, of course, is what is meant by structure. Clearly all computational scientists attempt to preserve the solution of their systems, the ultimate structure, which would mean this article would be about the proud history of many decades of computational plasma physics research that attempts to do just that. Thus it is necessary to rein in the purview, and this is done by concentrating on properties associated with or that emerge from the Hamiltonian and action principle (HAP) structure of plasma models, structure that ultimately originates from the fundamental electromagnetic interaction. Of interest are algorithms of relevance to plasma science that preserve various HAP properties that occur in finite systems, ordinary differential equations (ODEs), e.g., those that describe particle orbits or magnetic field lines, and in infinite systems of partial differential equations (PDEs), plasma field theories, e.g., various plasma fluid models and kinetic theories. Clearly ever since the advent of computation, researchers have attempted to design methods of structure preservation (see e.g. Refs. 1 and 2 for early examples), and there has been continual progress over the years (see e.g. Refs. 3 and 4 for more modern surveys). The perspective taken here is to describe in broad brush strokes the types of structure that can be preserved, illustrated with selected examples; a) Electronic Mail: morrison@physics.utexas.edu; Homepage: https://web2.ph.utexas.edu/˜morrison/ however, global value judgements and comparisons of the efficacy of various techniques will be avoided. Many references are cited, but constructing a complete bibliography would be a daunting if not impossible task, even with the given limitation of scope. Rather, a representative although somewhat biased selection of literature is given, but one that is hoped will provide a gateway into the field. The article also has the goal of highlighting the community of plasma researchers working in this area, as exemplified by workshops organized by Hong Qin, Eric Sonnendrücker, and myself: • Geometric Algorithms and Methods for Plasma Physics Workshop (GAMPP) May 13-15, 2014 Hefei, China (Hong Qin) • GAMPP II, September 12-16, 2016 in Garching, Germany (PJM, Hong Qin, and Eric Sonnendrücker) • Mini-Conference at this meeting: New Developments in Algorithms and Verification of Gyrokinetic Simulations (Amitava Bhattacharjee and Eric Sonnendrücker) It is possible there will be a GAMPP III in 2018.
The HAP structure of classical plasma physics originates from the relativistic action principle for a collection of charged particles interacting self-consistently with the electromagnetic fields they generate. This system has the dynamical variables (q i (t), φ(x, t), A(x, t)), where q i is the position of the ith particle, φ is the electrostatic potential, and A the vector potential. The Lorentz co-variant action is given by where (1) represents the relativistic kinetic energy, (2) provides the coupling between the particle dynamics and the fields, and (3) is the pure field contribution. Here s is a sum over species, N is the number of particles of each species, and species indices on N , q i , e, and m, are suppressed to avoid clutter. In (3) E and B are shorthand notation for the usual expressions in terms of the scalar and vector potentials Faraday's law and ∇ · B = 0 follow automatically from these expressions. Extremization of (1)-(3) yields Maxwell's equations coupled to particle dynamics. We will do this extremization by the equivalent procedure of functional differentiation, which we carefully define here because it will be used extensively later on. For the variation of S with respect to φ we consider the functional derivative δS/δφ, which is defined as follows: where δφ is an arbitrary function of its arguments, subject to admissible boundary conditions that assure the vanishing of terms obtained by integrations by parts. In (5), after δφ is isolated one can extract δS/δφ(x, t) = 0 by removing the integrals and δφ. The argument (x, t) of δS/δφ(x, t) is displayed so one knows which integrals to remove. (See Ref. 5 for an extended discussion of functional differentiation.) For the action of (1)-(3) the functional derivative δS/δφ defined by (5) yields Poisson's equation as expected. Similarly, the Ampère-Maxwell law follows from from δS/δA(x, t), where the current density is given by Lastly, the variation on particles is given by δS/δq i (t), and this produces the relativistic version of Newton's second law. This functional derivative involves integrating by parts in time and performing the integration over d 3 x so as to evaluate the Lorentz force on the particle positions. Thus the coupling term provides both the sources in Maxwell's equations and the force for the particle dynamics. The lion's share of plasma physics is embodied in the action principle of (1)-(3), so from one lofty point of view the discipline is complete. However, this naive point of view misses the beautiful emergent collective phenomena of plasma science that has required considerable effort to unravel. The unraveling has entailed various approximations, reductions resulting in various fluid and kinetic models that elucidate various plasma phenomena. Generally speaking, the reduced models that originate from the action S contain both dissipative and nondissipative processes, and the action principle origin can be obscured. Nevertheless, HAP properties bubble to the surface as part of essentially all plasma models, ranging from magnetic field line behavior, to MHD and more general magnetofluid dynamics, to various kinetic theories, and even the BBGKY hierarchy itself. We refer the reader to Refs. 5-13, which are reviews or contain a significant review component of HAP structure.
The paper is organized as follows. Section II discusses conservative integration, schemes that are designed to exactly conserve constants of motion. Next, Sec. III treats symplectic integrators that preserve canonical Hamiltonian structure, which is followed by Sec. IV where various other forms of canonical and noncanonical Hamiltonian integration are described, including variational integration and Poisson integration. Section V describes metriplectic integration, a proposal to build structure preserving integrators for systems that are conservative with both Hamiltonian and dissipative parts. Simulated annealing, a relaxation method that preserves structure while obtaining equilibrium states, is considered in Sec. VI. In Sec. VII GEMPIC algorithms for the Maxwell-Vlasov (MV) system are described. These are particle in cell (PIC) codes that preserve geometric structure, including the noncanonical Hamiltonian structure of the MV system, its Casimir invariants, and the geometry of Maxwell's equations. Here many of the features of previous sections are incorporated: it is shown how to obtain a semidiscrete Hamiltonian reduction of the MV system, yielding a finite-dimensional Hamiltonian system with noncanonical Poisson, which is then integrated with a Poisson integrator. Finally, the paper concludes with Sec. VIII. result in non-HAP reductions. Examples of this kind of reduction include truncation of the BBGKY hierarchy leading, to the Landau-like collision operators, and the derivation of quasilinear theory. For these models the Hamiltonian structure is lost, but conservation of energy and possibly other invariants, may be maintained. Sometimes one is aware that energy survives an approximation, by direct calculations after the fact or by construction, or that an unknown Hamiltonian structure actually exists. For some such systems, the long time states may be governed to a large degree by the preservation of a few invariants, whether or not the system is Hamiltonian.
The failure of an integrator to preserve an invariant is demonstrated by the phenomenon of energy drift. To understand this consider a system with a conserved energy or part of an energy given by a kinetic energy expression of the form where m is a particle mass and v its velocity. Let v calculated be the value of the velocity produced by some algorithm and δv the error, both after a single time step. Thus v calculated = v exact + δv and after a single time step Over many time steps one may expect the cross term v exact · δv to average to zero, while the term |δv| 2 will systematically give rise to a monotonically changing E. Indeed, for some standard integrators, energy may tend to increase dramatically. Over the years many researchers have used various techniques for enforcing the conservation of invariants. If one begins with a PDE then a common mode of development is to first obtain a semidiscrete approximation, an energy conserving set of differential equations, and then a next step would be to construct an integrator that enforces energy conservation. There are many examples of this is procedure; we mention spectral reduction 14 that was designed to give accurate spectral and cascades for turbulence described by the 2D Euler equation. For this system a semidiscrete reduction that exactly conserves energy and enstrophy was obtained, as desired for accurate spectra, and the resulting set of ordinary differential equations was then solved by a conservative integra-

A. Conservation Structure -Sets of Invariants
A set of M first order ordinary differential equations (ODEs), written aṡ where dot denotes total time differentiation and V is some vector field, possesses a conserved energy function where in the first equality of (10), which follows from the chain rule, repeated summation over the index a is assumed. For such systems withz being the initial value of the dependent variable z = (z 1 , z 2 , . . . , z M ). Similarly, any additional invariant I must satisfy V i ∂I/∂z i = 0. For a partial (integro) differential equation (PDE) we write where χ(µ, t) is a multicomponent field variable dependent on a 'spatial' variable µ of any dimension and time. For such PDEs, subscripts will denote partial differentiation. Systems of the form of (12) possess an invariant energy functional E[χ] iḟ where δE/δχ(µ) denotes the functional derivative as described in Sec. I, here defined in terms of the pairing where the first variation δE[χ; δχ], a linear operation on δχ, is the Fréchet derivative. The ' · ' of (13) is a shorthand for summation over the components of the field χ should it have more than one. The goal of a CI is to have the following exactly satisfied: withχ being the initial value of χ. Thus, the goal of a conservative integrator is to construct a numerical algorithm to maintain exactly (11) for an ODE or (15) for a PDE. Or, if there are other known invariants, to also preserve these as desired.

B. Conservative Integrators
Although various means have been used to preserve invariants, a systematic procedure amenable to various types of integrators was developed and applied in Refs. 14, 25, and 26. The procedure of these references can be used to adapt common integrators, such as Euler's method, predictor-corrector, etc., so as to exactly conserve an invariant or a family of invariants.
The procedure of Ref. 25 is quite simple. For a set of ODEs, first one modifies the system by adding a 'correction' vector field V C to (9), i.e., this system becomeṡ Next, it is required that and this must be at a rate fast enough so as to preserve the order of the desired integrator. Finally, one chooses V C so has to enforce exactly at each time step. When (18) is only a single constraint for a system of equations, there is freedom on how to enforce this; if a family of invariants exists, then one has more equations like (18), but there is still latitude for enforcing the family.
Several physical examples are solved with this type of CI in Refs. 14, 25, and 26; in particular, Euler's method is applied to Euler's rigid body equations, as a simple proof of principle example. Another example, one that illustrates well the role invariants can play in answering a question about the long time behavior of a system is the famous Kepler problem. Although, this problem is integrable, with an easily evaluated solution, it serves as a illustrative example. The computational question of interest is to obtain the path traced out by the orbit over a long period of time, a question that is critically dependent on invariance.
As usual, the Kepler problem can be reduced to motion in a plane, governed by the following three ODEs: where (r, θ) are polar coordinates, the gravitational potential φ = −K/r with constant K, is the magnitude of the angular momentum, and m is the reduced mass. This system conserves the energy plus an extra invariant, the Runge-Lenz, The Runge-Lenz serves to lock together the radial and angular frequencies for all values of the energy -as a consequence all bounded orbits (E < 0) are ellipses independent of chosen initial conditions. With the initial conditionsr,θ = 0, andv r = 0, A only has the xcomponent Figure 1 depicts calculated and known orbits for second order predictor corrector (PC) with and without the CI constraints (see Ref. 25 for details). In panel (a) plain PC is used and we see the calculated orbit deviate significantly from the known exact elliptical path. In panel (b) energy is conserved, but not the RL vector resulting in a wobbling precession. Panel (c) shows the case where both energy and the RL vector are conserved, giving rise to near perfect alignment of the calculated and know path. In panel (d) we have run the integrator with only energy conservation for a long time depicting the expected annular region sampled because of precession. Similarly Runge's RK4 gives poor performance in conserving Runge's vector. Thus, if one is interested in the planet path, straight integrations may give poor performance, while a CI can be a valuable tool.
Often physical systems are conservative, yet possess vector fields that have Hamiltonian and non-Hamiltonian components. For example, this is the case for the Vlasov equation with the Landau-Fokker-Planck collision operator. A formalism that describes the nature of this ubiquitous occurrence, called metriplectic dynamics in Ref. 27, will be discussed in Sec. V where metriplectic integrators (MIs) will be proposed. Before discussing MIs, a combination of a CI and a HI, various forms of HIs will be discussed in Secs. III and IV.

III. SYMPLECTIC INTEGRATION (SI)
Whereas the geometric content of CI is fairly minimal, i.e., the preservation of a few constraints restricting the dynamics to submanifolds, SI is both global and local in nature and amounts to preserving the entire symplectic geometry of phase space, a consequence of the Hamiltonian form. Symplectic geometry allows one to measure sizes of particular two-dimensional and higher even dimensional objects, akin to the measurement of lengths and angles of Riemannian geometry.
The idea behind SI is old: a written account occurred in 1956 in a report of de Vogelaere 1 and was widely known and used in the accelerator physics community. In fact, prior to this paper, symplectic algorithms were used unbeknownst to their originators, viz. the method used by Verlet and Störmer and the leapfrog method. A flurry of papers on symplectic integrators appeared in the 1980s by various authors, notably Ref. 28. A general and interesting historical account from an accelerator physics point of view is given in Ref. 29. In the plasma physics community SI in the infinitedimensional context has been employed, where it has been used to integrate semidiscrete Hamiltonian equations for describing the two-stream instability 30,31 and for understanding facets of the nonlinear plasma evolution in the single and multiwave models 32-34 of the bump-on-tail instability.

A. Canonical Hamiltonian Structure
The concept of SI is a consequence of Hamilton's equations, a set of 2N first order differential equations determined entirely by the Hamiltonian function H(q, p) depending on the canonical coordinates q and momenta p. Sometimes it is convenient to let z = (q, p) and rewrite (24) in tensorial forṁ where the canonical Poisson matrix and bracket are given, respectively, by with 0 N being an N ×N matrix of zeros, I N being the N × N identity matrix, and f and g being arbitrary functions of z.
An important idea of Hamiltonian mechanics is to effect special coordinate transformations to obtain solutions or approximate solutions. Upon representing a coordinate change as z = z(z), 2N functions of 2N variables, with the inversez =z(z) Eqs. (25) becomė where we see that H transforms as a scalar and the Poisson matrix J as a contravariant 2-tensor, If the second equality of (27) is satisfied, then the transformation is called a canonical transformation (CT) (or symplectomorphism in mathematics) and the equations in the new coordinates have clearly identified new coordinates and momentaz = (q,p) and, consequently, the form of Hamilton's equations is preserved: A main theorem of Hamiltonian dynamics is that the time advance map is itself a CT; i.e., if z = z(z, t) is the solution of (25) at time t with initial conditionz, then for any fixed t the coordinate change z ↔z is a CT.
To understand the importance of this main theorem, consider three nearby trajectories, z r (t), z r (t)+δz(t), and z r (t) + δz(t), as depicted in Fig. 2 and let us investigate the following area-like quantity: where A simple calculation for our three Hamiltonian trajectories reveals Thus, δ 2 A, whatever it means, is preserved by the Hamiltonian dynamics. To understand its meaning, consider the case where N = 1. Examination of Fig. 3 reveals that δ 2 A measures the area of a parallelogram, i.e., |δ 2 A| = |δpδq − δqδp| = |δz × δz|. For N > 1 this quantity is a sum over such areas, δ 2 A = i δp i δq i − δq i δp i , and is called the first Poincaré invariant. In more modern notation we would call this the symplectic two-form and write it as ω c = i dp i ∧ dq i . In tensorial notation the canonical symplectic 2-form = (canonical Poisson matrix) −1 , i.e., ω c αγ J γβ c = δ β α . There are consequences of (28): upon summing (integrating) over the infinitesimal areas δ 2 A one sees that extended suitable two-dimensional ribbons in phase space are preserved under Hamiltonian dynamics. Since an initial area in any q − p subplane is preserved, one can form a 2N -dimensional box, i.e., a volume that will be preserved. This is of course the famous Liouville theorem of Hamiltonian mechanics, which is a consequence of the preservation of δ 2 A but not vice versa. In addition to regions of dimension 2 and 2N , Hamiltonian dynamics preserves everything in between, i.e., all the Poincare invariants which are surfaces of dimension 2, 4, 6, ... 2N.
Another family of invariants, intimately related to δ 2 A, are the loop integrals where γ is any closed curve in phase space. By a generalization of Stokes theorem, these loop integrals and the areas δ 2 A are related. Thus, Hamiltonian dynamics has a phase space circulation theorem, similar to that of ordinary fluids.
It should be borne in mind that if Hamiltonian dynamics did not have the properties described above, the nature of phase space would be quite different. For example, surfaces of section of magnetic field lines would look very different, and lack the similarity we see. Also, non-SI of Hamiltonian systems can give rise to spurious damping or growth, a problem of concern when long time integration is desired.
For infinite-dimensional canonical Hamiltonian systems with conjugate fields χ = (ψ, π), (24) becomes where H is now a Hamiltonian functional like the kinetic energy in a fluid d 3 xρ |v| 2 /2, partial derivatives are replaced by functional derivatives, and the Poisson bracket operating on two functionals F and G is where the fields χ depend on µ as well as time.

B. Symplectic Integrators
Various approaches to symplectic integration have been followed. Consider a simple early example, starting from (9) with a possibly time dependent vector field V . Suppose one could find a time dependent coordinate change that takes this system intȱ z =V ≡ 0 : (32) i.e., in the new coordinates the solution is constant in time,z(z 0 , t) =z 0 . Upon denoting such a timedependent coordinate change by z = Φ(z, t), the initial conditions of the problem in the two system of coordinates are related by z 0 = Φ(z 0 , t 0 ), with t 0 being the initial time. If z(z 0 , t 0 , t) were the solution to the original problem (9), then the solutions to (9) and (32) would in general be related by Φ z(z 0 , t 0 , t), t = z(z 0 , t 0 , t). But since the solution to (32) is constant in time, we have and the coordinate change is seen to solve the problem. All of the above is moot unless one has a means of obtaining Φ, but this can be attempted for Hamiltonian systems by making use of mixed variable generating functions, which not only provides an avenue for approximating Φ but does so in terms of a canonical transformation. For clarity consider the case of a single degree of freedom and seek a canonical transformation z = (q, p) ↔z = (q,p). This can be done by introducing a mixed variable generating function of type three, F 3 (q, p, t), and giving the transformation as follows: Because the transformation has explicit time dependence, there is a new Hamiltonian given bȳ and if this new Hamiltonian can be made constant or zero, then a Hamiltonian version of (33) applies. Consider the simple case where H = p 2 /2 + V (q, t), with a time-dependent potential V , and choose F 3 = −q p − t(p 2 /2+V (q, 0)). By (35) this gives rise to the new Hamil-tonianH where V q = ∂V /∂q. Equation (36) gives upon expansion for small t,H (37) is O(t) the same is true for Hamilton's differential equations, with the solutions Thus, the map is symplectic and accurate to first order, which means that if we take t to be the step size ∆t, then the Taylor series in ∆t of the approximate solution matches the actual to first order. While the above example is rather simple, it serves to reveal a path for obtaining higher order methods, such as the 3rd order method described in Ref. 28 and the higher order methods of Refs. 35-40. One feature of transformations of given by (34) is that these relations for general Hamiltonians are implicit and approximate solution of these formulas by e.g. Newton's method may lead to a violation of the symplectic property. A similar difficulty arises with Lie transform approaches (see e.g. Ref. 41), where the canonical transformation is generated by a series as follows: Truncation of this series again in general results in a transformation that is not exactly canonical. A way around this problem is to make an approximation using Cremona maps (see e.g. Ref. 42). Using the notation {·, g} = X g for a Hamiltonian vector field one incorporates the Campbell-Baker-Hausdorff theorem to write where the H i are special polynomials for which the infinite series truncates, thus guaranteeing exact symplecticity.
To summarize, every time step of a symplectic integrator is a CT and, consequently, the geometric structure described in Sec. III A is exactly preserved. However, it is not the case that energy is exacly conserved as it is for CI. Instead energy is 'shadowed', which means there is a nearby Hamiltonian that is exactly conserved. In this way one can bound the deviations of the energy, in accordance with the order of the scheme, and avoid energy drift (see e.g. Ref. 4). A criticism of symplectic integration is that making the time step adaptable can be difficult. One approach that overcomes this difficulty and allows for adaptive time steps is that of Ref. 43.
For Hamiltonian field theories, PDEs that have canonical Poisson brackets of the form of (31), semidiscrete Hamiltonian reductions are straightforward. One can simply expand the fields (ψ, π) in terms of any complete set of basis functions and then truncate to obtain a finitedimensional canonical Poisson bracket for ODEs, as is typical with spectral methods. In this way it is easy to show that the truncated system has finite-dimensional Hamiltonian form. (See Refs. 44 and 45 for an early comparison of canonical vs. noncanonical Hamiltonian reductions.) Once the semidiscrete Hamiltonian reduction is at hand, one can turn to SI as a possibility for the time advancement of the ODEs.

IV. HAMILTONIAN INTEGRATION (HI)
Although SI as described in Sec. III is the most common type of integration that preserves Hamiltonian form, other methods and approaches exist. In this section we discuss some of these.

A. Variational Integration (VI)
The main idea behind VI is to make approximations in an action functional, giving rise to time step equations as Euler-Lagrange equations. In this way one can obtain symplectic integrators that may also preserve other invariants such as momenta.

Variational Structure -Canonical
Consider the canonical case where we start from Hamilton's variational principle of mechanics. Recall, Lagrange's equations arise from extremization of the following functional over paths: where L is the Lagrangian. Thus, δS/δq = 0 yields With the Lagrangian L = m|q| 2 /2−φ(q, t), Lagrange's equations are equivalent to Newton's second law.
For convex Lagrangians, the Legendre transform makes the identification of Lagrange's N second order equations with Hamiltonian's 2N first order equations. Recall, one defines the canonical momentum by p i = ∂L/∂q i and if then by the inverse function theorem one can solve for all the momentaq i as a function of the p i and construct Then with the Hamilton of (44), Hamilton's equations of (24) are equivalent to Lagrange's equations (42). When W = 0 one cannot solve for all the momentȧ q i as a function of the p i . For this more complicated situation one must employ Dirac constraint theory (see e.g. Refs. 46 and 47), which we will briefly consider in Sec. IV A 4.

Variational Integrators -Canonical
For finite-dimensional systems, the strategy for obtaining a variational integrator is to discretize the time integration of the action and then vary. As a simple example, assume t 1 − t 0 is small and use the trapezoidal rule to approximate the action as where a linear trajectory approximation has been used. Next, we add up to make a discrete action as follows: Finally, upon variation, which now is merely differentiation, one obtains a symplectic map for each time step. There has been much work on VI with an entry point into the literature being Ref. 48. Historically, the idea to discretize actions is old and has been used for various purposes in Hamiltonian dynamical systems theory, both theoretically and computationally. Theory work proposing VI appeared in the 1980s (e.g. Ref. 49), but actual working codes for the more difficult problem of finding periodic orbits were given in Refs. 50-52, where use as an ODE integrator was also proposed.
Variational integrators for PDEs proceed by discretizing in the spacial variable as well as time. An example of this for MHD was given in Ref. 53 where discrete exterior calculus was used to discretize the label of Newcomb's MHD Lagrangian. 54

Variational Structure -Noncanonical
Noncanonical variational structure occurs when the Legendre transform fails. Consider the action functional of the form of (41) with the Lagrangian L =q a A a (q) − φ(q, t) , a = 1, 2, . . . , M which would produce the canonical momenta p a = ∂L/∂q a = A a . For this case one cannot solve for any of theq a in terms of the p a because the Lagrangian is totally degenerate, i.e., by (43) W ≡ 0 with rank zero. However, for this totally degenerate case the Euler-Lagrange equations are with and one has a 'half-sized' Hamiltonian system in noncanonical variables with φ as the Hamiltonian. This example is of interest in plasma physics because guiding center orbits are governed by Littlejohn's Lagrangian 55 (see also Ref. 56) that has the following form: where A(x, t) is the vector potential.

Variational Integrators -Noncanonical
One can use Dirac's constraint procedure 57 to obtain a canonical Hamiltonian system (see Ref. 58 where this is done for (49)) but for the purpose of VI one can directly discretize the action with (49)  The particular case of the Maxwell-Vlasov system will be studied in Sec. VII.

Noncanonical Hamiltonian Structure
The noncanonical generalization of the Hamiltonian form of (25) is given bẏ where the following noncanonical Poisson bracket replaces (26): For a bracket of the form (51) to be a good Poisson bracket it must have the following properties for all functions f, g, h: • antisymmetry: • Jacobi identity: where cyc means cyclic permutation over f gh in the first expression and over abc in the second.
This noncanonical generalization of the Hamiltonian mechanics is reasonable because of an old theorem due to Darboux, which states that if det J = 0 then there exists a coordinate change that (at least locally) brings J into the canonical form J c of (26). Recalling the J transform as a rank 2 contravariant tensor, this canonizing transformationz ↔ z would satisfy However, the more interesting case is the one studied by Sophus Lie where det J = 0. This case is degenerate and gives rise to Casimir invariants (Lie's distinguished functions), which are constants of motion for any possible Hamiltonian that satisfy Because of the degeneracy, there is no coordinate transformation to canonical form; however, a theorem known to Lie (see e.g. Refs. 69 and 70) which we call the Lie-Darboux theorem states there is a transformation to the following degenerate canonical form: Instigated in a major way by the noncanonical Poisson brackets for plasma models, manifolds with the addition of degenerate Poisson bracket structure, known as Poisson manifolds, have now been widely studied (see e.g. Ref. 71). The local structure of a Poisson manifold is depicted in Fig. 5, where it is seen that the phase space is foliation by the level sets of the Casimir invariants. For an M -dimensional system there exist M −2N Casimir invariants, and an orbit that initially lies on such a surface defined by the level sets of the initial Casimir invariants remains there. These surfaces, called symplectic leaves, have dimension 2N and the phase space is generically foliated by them.
Lie-Poisson brackets are a special form of noncanonical Poisson bracket that typically appear in matter models where F and G are arbitrary functionals, the Poisson matrix becomes a Poisson operator J , and the fields χ = (χ 1 , . . . , χ N ) are not divided into canonical pairs (ψ, π) as in (31). In general the Poisson operator may depend on the fields χ(µ, t), but most importantly J must still satisfy infinite-dimensional versions of (52) and (53), the latter of which may be challenging to show (see Ref. 6). A example of an infinite-dimensional noncanonical Lie-Poisson bracket is that for the Vlasov-Poisson system, 6,45,72 where f (x, v, t) is the distribution function and the Poisson operator for this case is given by Here the value of the electron charge and mass have been set to unity and consequently the Hamiltonian of the system, the total energy, is given by where|∇φ| 2 is a shorthand for the quadratic Green's function expression in terms of f . Using δH/δf = v 2 /2 − φ(f ; x) = E, the particle energy, the Vlasov-Poisson system is expressed in Hamiltonian form as In Sec. VII we will see that the Poisson bracket of (58) is a part of the more complicated full Vlasov-Maxwell Hamiltonian structure.

Poisson Integrators
Symplectic integrators preserve, step-by-step, the canonical Poisson matrix J c (cf. (27)). Thus, it would seem natural to consider transformationsz ↔ z that preserve the form of a given noncanonical Poisson matrix, i.e. transformations for which Such transformations are called Poisson maps. If an integrator is constructed to maintain (60) step-by-step, then up to coordinate changes there is a preserved canonical J c and one has all the structure of SI maintained. To see this, suppose (60)  For noncanonical Hamiltonian field theories, PDEs that have canonical Poisson brackets of the form of (57) with J depending on χ, semidiscrete Hamiltonian reductions are not straightforward. In fact, prior to the publication of Ref. 68 considerable effort was spent on trying to find a finite-dimensional Hamiltonian projection of MHD by expansion in a complete set of basis functions followed by various truncations, but it was learned that such a procedure results in the failure of the Jacobi identity. This is more easily seen for the simpler case of (58), which was tried in Ref. 44 by expansion in a Fourier basis -there it was noted for this projection that the process of truncation destroys the Jacobi identity. However, it was also noted there that expansion in terms of particle degrees of freedom results in a good Poisson bracket. This procedure will be used in Sec. VII where we will discuss GEMPIC, which is a kind of PI in the PDE context for the full Maxwell-Vlasov system.

V. METRIPLECTIC INTEGRATION (MI)
As noted in Sec. II reductions of the general HAP formalism of (1)-(3) and its associated Hamiltonian description may give rise to dissipation, which generally results in vector fields that have Hamiltonian and non-Hamiltonian or dissipative components. Generally speaking, dissipation takes one of two forms in continuum matter models that describe media, such as fluids and plasmas. The first form is exemplified by the Navier-Stokes equation where viscosity removes energy from the system, while the second form is exemplified by the fluid theory that conserves energy but allows for viscous heating, thermal diffusion, and entropy production by using a general equation of heat transfer. 77 Transport equations with collision operators, such as the Boltzmann, Landau-Fokker-Planck, or gyrokinetic 78 collision operators, are of the second form because they conserve mass, momentum, and energy, but produce entropy. Metriplectic dynamics 27,79 is general formalism that embodies systems like these that in a real dynamical sense embody both the first and second laws of thermodynamics, i.e., have energy conservation and entropy production.  85). There has been much subsequent work, notably by Grmela and collaborators, 86 and myself and others. [87][88][89][90][91] A metriplectic vector field has two parts, a part that is Hamiltonian and a part that is a degenerate gradient or metric flow. Consider first a metric flow on a finitedimensional phase space manifold, which in coordinates has the forṁ z a = g ab ∂S ∂z b = (z a , S) , a = 1, 2, . . . , M where S is an 'entropy' function and the symmetric bracket (f, g) = (g, f ), defined on arbitrary phase space functions f and g, is defined by Metriplectic dynamics requires two things of the matrix g: (i) that it be positive semi-definite so that the dynamics satisfies dS dt = (S, S) ≥ 0 , which can be used to build-in asymptotic stability, i.e. an 'H-theorem', and (ii) degeneracy so as to conserve an energy H that will act as a Hamiltonian, Paired with the gradient flow is a noncanonical Hamiltonian flow of the form of (50) with Poisson bracket (51), and the two together define a metriplectic vector field generated by some function F as follows: The metric and Hamiltonian components are then mated by requiring F = H +S with the entropy S selected from the set of Casimir invariants of the noncanonical Poisson bracket. Given this structure we have the following: • a 1st Law: Thus metriplectic dynamics is a dynamical paradigm that embodies both the first and second laws of thermodynamics. A finite-dimensional example of a metriplectic system based on the Hamiltonian structure of the free rigid body was given in Ref. 27. For this example, the energy of the rigid body is conserved, but the magnitude of the angular momentum, which acts as the entropy, monotonically changes until the system relaxes to one of its stable rotations about a principal axis.
A general infinite-dimensional symmetric bracket has the form where for metriplectic dynamics G ij is chosen to guarantee positive semidefiniteness (F, F ) ≥ 0, the symmetry (F, G) = (G, F ), and to build-in the degeneracy condition (H, F ) = 0 for all functionals F . A most important class of metriplectic systems are transport equations of the form where the first two terms on the right-hand side of (65), the Hamiltonian Vlasov terms for electrons, can be generated by the Poisson bracket of (58), while the metriplectic formalism is completed by using the following to generate the collision operator, ∂f /∂t| c : 27 .
In (66) z = (x, v), i, j = 1, 2, 3, and it remains to tailor this expression to fit a specific collision operator. Because Casimir invariants are candidate entropies, and the specific choice of the entropy density function s will determine the state to which the system relaxes, i.e., the state of thermal equilibrium. If the tensor T ij is given by with the functions M and s related by the following entropy compatibility condition: then it only remains to determine the tensor w ij (z, z ). Any choice that satisfies the symmetry conditions will assure conservation of mass, momentum, and energy. Finally, according to the metriplectic prescription, the system will relax to the state determined by extremization of F = H + C, i.e., where recall E = δH/δf . If the solution to (68) is stable, which is assured by the Kruskal-Oberman 92 monotonicity stability condition, then at least on a formal level dS/dt ≥ 0 will cause S to increase until the monotonic state f = (ds/df ) −1 (−E) is achieved. Thus, this collision operator can be tailored so that the system relaxes to any distribution function that is monotonic in the energy. Proving the relaxation property is detailed, but it can be shown by paralleling the arguments of Ref. 93. The Landau collision operator is obtained by choosing the kernel In this way, one can construct a collision operator that relaxes to any equilibrium state monotonic in the energy.

B. Metriplectic Integrators
MI should at once have a PI component for the Hamiltonian part of its vector field, preserving the symplectic nature of phase space and the Casimir foliation, while it should have a CI component for the dissipative part that exactly conserves invariants. Various ideas come to mind on how to achieve this, e.g. by splitting the time step. At present we know of no published metriplectic integrators, finite or infinite. However, very recently progress has been made and examples of MI are on the horizon. 95,96 VI. SIMULATED ANNEALING (SA) By simulated annealing we mean a numerical relaxation procedure where structure is used to obtain physical equilibrium states by constructing usually nonphysical dynamics that possess certain geometrical constraints. In this section we discuss two kinds of SA: one based on the metriplectic dynamics discussed in Sec. V and the other called double bracket dynamics (introduced in Ref. 97), which is also based on Hamiltonian structure.

A. Symmetric Bracket Structure and SA
Both metriplectic and double bracket SA use a symmetric bracket to achieve relaxation. Such brackets can be constructed in a variety of ways, giving rise to a bracket of the form of (64) or its finite-dimensional analog. Since such brackets are interesting independent of their SA application, we explore their construction in some detail.
One way of building in degeneracy is with brackets of the following form: where K ij is positive semidefinite and { , } is any Poisson bracket. This form, given in Ref. 98, has a very general geometric significance. 88 To see this suppose P is any phase space manifold that has both Hamiltonian and Riemannian structure. Because P has Riemannian structure, given any two vector fields X 1,2 ∈ X(P ) the following is defined: If the two vector fields are Hamiltonian, in terms of either a degenerate or a nondegenerate Poisson matrix, say X F , X G , then we natually have a symmetric bracket given by which is an abstract way of writing (69). Thus, natural symmetric brackets exist for Kähler manifolds, which naturally have metriplectic or double bracket flows. If P is a manifold with any degenerate Poisson bracket, then the associated Casimir invariants C must satisfy (F, C) ≡ 0 for all F . One possibility is the case where { , } is a Lie-Poisson bracket. Another possibility is to build a Dirac bracket with desired degeneracies, as was done in Ref. 98 for SA. Given any bracket, canonical or not, one can build in Casimir invariants by Dirac's construction. This is easily done in finite dimensions for the case of building in invariants, say D i . For this case one calculates the matrix Ω ij = {D i , D j }, which must be invertible, whence the Dirac bracket is given by where it is easily seen that  (69). Similarly, one could use antisymmetric brackets with more slots. A triple bracket will be used for metriplictic SA in Sec. VI B 1.
Metriplictic SA proceeds just as in Sec. V, by employing a metriplectic dynamical system to numerically effect the variational principle δF = δ(H + S) = 0, and finds the state consistent with the initial energy or other invariants. On the other hand, double bracket SA extremizes H subject to constancy of all Casimir invariants. We will give example of these in Secs. VI B 1 and VI B 2, respectively.
B. Examples of Simulated Annealing

Metriplectic SA
As a numerical tool, metriplectic SA is in its infancy. The work described in this section is joint work with Glenn Flierl 104,105 and is preliminary in nature. Additional preliminary yet promising results have been obtained with and by collaborators at IPP in Garching 106 where this approach is being explored for obtaining MHD equilibria. It is important to clarify that metriplectic SA differs from the MI of Sec. V B in purpose: metriplectic SA is a tool for finding equilibrium states via, in general, unphysical dynamics, while MI is an integrator that preserves the MI geometric structure.
The example at hand concerns quasigeostrophic flow with topography, as governed by the following PDE for the potential vorticity ζ(x, y, t): where in the present section [f, g] = f x g y − f y g x and the stream function ψ satisfies ζ = ∇ 2 ψ + T , with T (x) modeling bottom topography. (See e.g. Ref. 107.) Observe, this system also has an interpretation for guiding center plasmas and is closely related to the Hasegawa-Mima equation for plasma drift waves. 108 Our goal is to find an asymptotic time-independent state consistent with conservation of energy. We seek this state by a procedure that extremizes the enstrophy, in particular, we will construct a metriplectic system that minimizes (Note, with a flip of a sign it can be made to maximize Z, which is not important since a physical process is not being modeled.) We construct our symmetric bracket by means of the following antisymmetric triple bracket: 102,103 where A, B, C are arbitrary functionals and we use the notation A ζ := δA/δζ. That (74) is completely antisymmetric is easily shown for the periodic boundary conditions we use by integration by parts. From (74) Now we choose H = − d 2 x ψζ/2, the physical quasigeostrophic energy, and suppose giving rise to which coincidentally is an equation of the form of that in Ref. 109, which was argued to be of physical origin.
As an example we consider a ridge described by T = e −x 2 /2 and integrate forward using a pseudospectral code until the system approaches a relaxed state. Figure 6 shows an initial dumbbell shaped concentration of vorticity being sheared out along the ridge as the artificial time (determined for convenience by the parameter γ) progresses.

Double Bracket SA
The original double bracket formulation of Ref. 97 was generalized in two ways in Ref. 98: by introducing a smoothing kernel and by imposing constraints using Dirac brackets. This allowed a much wider class of states to be found. In subsequent work this formalism has been applied to reduced MHD. 106,[110][111][112] Here, in Fig. 8, we display a result from Ref. 98 for Euler's equation for ideal vorticity dynamics in two dimensions. In this example we seek a rotating state for the scalar vorticity that has twofold symmetry. Double bracket SA is used to find this as an equilibrium in a rotating frame of reference. We refer the reader to Ref. 98 for details and many examples.   (c) t = 120.  126 GEMPIC is a functioning Maxwell-Vlasov PIC code that has been tested and is available in the SelabLib library. 127 Many additional references to the relevant literature can be found in Ref. 124.
GEMPIC is a culmination that uses many of the ideas about structure presented in previous sections of this paper. A semidiscrete reduction of the full Maxwell-Vlasov system is obtained, either from its variational principle or by reduction of its noncanonical Poisson bracket, resulting in a finite-dimensional noncanonical Hamiltonian the-ory. The finite-dimensional theory exactly conserves invariants, including finite-dimensional versions of the energy and Casimir invariants, and a symplectic method, specifically a PI splitting method, is devised for its numerical solution. Thus GEMPIC brings together many ideas presented here, but it also uses ideas for the solution of Maxwell's equation based on its geometric structure.
In Sec. VII A we review the HAP structure of the Maxwell-Vlasov system, while in Sec. VII B we describe its semidiscretization and resultant Hamiltonian structure. Obtaining this structure is complicated, so only a sketch is given, with a goal of giving the flavor of the calculations of Ref. 124. The section ends with Sec.VII C, a discussion of the splitting method.

A. Maxwell-Vlasov HAP Structure
An action principle for the Maxwell-Vlasov system was proposed by Low, 128 which can be viewed as a continuum or smoothed out version of the nonrelativistic limit of the discrete particle action principle of (1)-(3). Low's action treats the particles by means of a Lagrangian or material continuum variable labeled by its initial condition, which replaces the discrete particle index i. For simplicity we suppress entirely the sum over species with its presence being implicit. Additional material on action principles for the Maxwell-Vlasov system can be found in Refs. 7 and 129.
Suppose at t = 0 a 'particle' is located at every point z := (x,v) of phase space. Thus, in the action of (1)-(3) we replace the discrete label by this continuum label, q i → q(z, t), and we replace the sums as follows: where dz = d 3x d 3v and f 0 can be viewed as a probability or number density attached to each point of phase space. This procedure results in the following action: where now all variables are fields, the particle phase space field q(z, t), and the electromagnetic fields φ(x, t) and A(x, t).
The functional derivative δS[q, φ, A]/δq(z, t) produces the partial differential equation for Lagrangian particle orbits. This equation can be shown to be equivalent to the Vlasov equation if we define the usual Eulerian expression for the distribution function as f (x, v, t) := f 0 (z), wherez is the initial condition of the particle located at z at time t, i.e.,z is determined by inverting z = (x, v) = (q(z, t),q(z, t)). This inversion of the map z ↔ z is possible by uniqueness and, moreover, it has unit Jacobian.
As with the discrete particle action of (1)-(3), Faraday's law and ∇·B = 0 follow because of the introduction of the potentials, φ and A with the relations (4), and the remaining two Maxwell equations follow from δS/δφ(x, t) and δS/δA(x, t), but now with the sources Thus with the interpretation above, Low's action has the Maxwell-Vlasov system as its Euler-Lagrange equations, which for completeness we record below, Note that (77) reverts to (1)-(3) if we suppose that initially particles are located at N isolated points, i.e. upon using q i (t) := q(z i , t), and setting all the 'weights' w i = 1. The counterpart of the Low action is the purely Eulerian gauge-free noncanonical Hamiltonian theory conceived of in Ref. 72 that uses the natural field variables (f, E, B). This theory uses the conserved Maxwell-Vlasov-energy as its Hamiltonian, together with the following noncanonical Poisson bracket where [f, g] = ∇f ·∂ v g −∇g ·∂ v f . With the noncanonical bracket of (86) and the Hamiltonian (85), one obtains the Vlasov-Maxwell system of (80), (81), and (82) Thus the Jacobi identity is satisfied for arbitrary functionals F, G, H defined on divergence-free magnetic fields. The constraints of (83), satisfy {F, C} = 0 pointwise for all F . Thus, C E is a Casimir invariant, while C B could be called a semi-Casimir, because it is intertwined with the satisfaction of the Jacobi identity per (88). The relativistic Vlasov-Maxwell theory follows upon replacing m|v| 2 /2 in the first term of the Hamiltonian (85) with −mc 2 1 − |v| 2 /c 2 , and the theory can be written in manifestly covariant form, 132,133 but this will not be pursued further here.

B. Discrete HAP Maxwell-Vlasov
Let us now consider semidiscrete reductions of the Maxwell-Vlasov system, systems of ODEs that are designed to approximate the dynamics. Many approaches for this have been pursued, but that most natural to a physicist is to start from the action principle, because this provides a direct avenue for maintaining the exact conservation of energy and other invariants in the reduced theory. Many such approximations have been made in plasma physics, e.g. that of Refs. 129, 134, and 135, where the latter two were specifically designed for obtaining computable models, the single-and multi-wave models of plasma physics. 136 However, specifically for the purpose of finding a semidiscrete reduction for solving the Maxwell-Vlasov system, the early work of Ralph Lewis on variational PIC stands out. [137][138][139] His approach to PIC did not take hold, most likely because NGP (nearestgrid-point) was used for the force which proved to be prohibitively noisy for computers of the day.
The idea of Lewis was to insert expansions for the fields and particles in terms of sums over basis functions into (77), truncate the sums, integrate overz, and then vary as in Hamilton's principle of mechanics to obtain sets of ODEs that exactly conserve invariants. The ODEs obtained by this procedure describe both particle and field degrees of freedom with their coupling provided automatically by the coupling terms of (77). This variational approach of Lewis was revived in more recent works, 113,114,117,118 where the use of better basis functions and/or the introduction of shape functions was seen to dramatically reduce noise and other sources of error. For example, upon inserting Fourier series for the potentials, into the field term of (77) and representing the particles by a shape function S as follows: where d 3 x S = 1, yields an action of the form Then, extremization gives Lagrange's equations in terms of the Lagrangian L in the usual manner. By construction, these equations are conservative and have Hamiltonian structure. GEMPIC differs from the above procedure by using advances in the integration of Maxwell's equations that exploit its intrinsic geometry, viz., that E is a 1-form and B is a 2-form. Because of its desirable properties, the early finite difference discretization of Yee 140 was generalized to a complete discrete differential calculus in later work, 141 and these ideas can be incorporated into a Maxwell-Vlasov PIC code. However, GEMPIC adapts Finite Element Exterior Calculus (FEEC), a finite element or spline approach that is described in detail in Refs. 142 and 143. The upshot of this approach is that one obtains discrete versions of the grad, div, and curl operations, in terms of matrices, with properties like curl grad =0 maintained. An important aspect of this approach is that there are now many finite element spaces, splines etc. that achieve this. In fact, the Fourier bases of (90) is a special case of the general framework.
To see how this progresses, consider the expansions where Λ 1,2 i are the 1-form and 2-form bases with the desired properties. Given these properties, Eqs. (4) take the discrete matrix forms where e and b have as components the e i and b i , ϕ and a are the amplitudes for expansions analogous to (93) for the potentials φ and A, and G and C are the discrete gradient and curl operators. The discrete form of Maxwell's equations become the following ODEs in matrix form: where D is the discrete divergence operator, j and are suitable n-tuples representing the charge and current densities (3-forms and 2-forms) defined in terms of the particle degrees of freedom, means transpose, and the matrices M 1,2 are defined by the pairing of the bases elements Λ 1,2 which in general are not orthogonal. This discretization of Maxwell's equations is coupled with a particle discretization like (91) to obtain a finitedimensional system. That this resulting system possesses noncanonical Hamiltonian form, follows upon projecting the Poisson bracket of (86). This is done by using functional chain rule expressions such as those given in Refs. 44 and 45, e.g., with v i =q being the canonical momentum conjugate to q i . To see the details of this calculation we refer the reader to Ref. 124 where (X, V) denote the 6N particle degrees of freedom and z in index form, z a , has a = 1, 2, . . . , M with M including the particle and e and b degrees of freedom.
with I 3×3 denoting the 3 × 3 identity matrix and M p and M q being diagonal matrices with elements m i w i and e i w i , respectively. The matrices B(X, b) and Λ 1 (X), depend on components of z and are given in terms of Λ 1,2,3 i with Λ 3 i being the 3-form basis elements (see Ref. 124 for details). Although complicated, in the end J is just a matrix of functions of the dynamical variables, which was shown directly in Refs. 122 and 124 to satisfy the Jacobi identity, and thus with a Hamiltonian we have a finite, yet very large, noncanonical Hamiltonian system of the form of (50).

C. Hamiltonian Splitting
The discretized form of the Hamiltonian of (85) has three parts,Ĥ =Ĥ p +Ĥ E +Ĥ B , witĥ This decomposition is the basis of the splitting method developed in Refs. 119, 144, and 145. The essential idea is that each of the subsystemṡ can be solved exactly and each is a Poisson map. Thus, their composition is a Poisson map, and we can construct a time step that exactly conserves the Poisson bracket. This is an example of the PI discussed in Sec. IV B. The integrator is constructed as a Lie series as in (40). For example a first order integrator is given by z(∆t) = e ∆t X E e ∆t X B e ∆t Xp 1 e ∆t Xp 2 z(0), where the Xs are Hamiltonian vector fields, e.g., For different arrangements of exponentials, one can obtain a second order integrator -we direct the reader to the references for details.

VIII. CONCLUSIONS
In this paper HAP formulations of plasma dynamical systems have been reviewed, and methods for integrating differential equations that preserve such structure have been surveyed. In Sec. II CI schemes that exactly conserve constants of motion were discussed. Then in Sec. III SI, which preserves canonical Hamiltonian structure, viz. the areas of (28), and the associated Poincaré invariants and loop integrals of (29), was described. This was followed by two other kinds of HI, the VI of Sec. IV A, based on approximating variational principles, and the PI of Sec. IV B for Hamiltonian systems with the noncanonical Poisson bracket formulation that is typical of plasma models. Next MI was proposed in Sec. V, an integration scheme that would preserve the structure of systems that are conservative with both Hamiltonian and dissipative parts. Simulated annealing, the relaxation method for obtaining equilibrium states was described in Sec. VI. Examples of metriplectic SA and double bracket SA were given. Finally, in Sec. VII the GEMPIC algorithm for the Maxwell-Vlasov (MV) system are described. The GEM-PIC code is a fitting finale, for it is a culmination of the work of many researchers that exactly preserves geometric structure described in previous sections of this paper.
Given the many time-honored computational methods, it is important and obvious to note that structure preserving methods like those described here may be unnecessary. However, for some problems it may be easy and inexpensive to try one of the methods described here, and depending on the problem being addressed it might prove to be superior, but no structure preserving integrator is going to be a panacea for all computational ills. For some problems such techniques are useful, and may even be essential to obtain an accurate solution.
There are many avenues for future work, both of a general nature and specific to plasma physics. Since suggesting the idea of MI, several researchers have already taken up the challenge and have made progress. It appears in the PDE context that MI will be useful for handling collision operators. Application to drift kinetic or gyrokinetic theories is underway 146 , but because of the different variational and Hamiltonian structure of these theories (see e.g. Refs. 147 and 148), modifications are necessary. Lastly we mention that Hamiltonian reductions like that of GEMPIC and others 149 open the field of discrete gyrokinetics. Conventional gyrokinetics reduces an infinite-dimensional theory, the MV system, to another reduced infinite-dimensional theory. Given a discretization like GEMPIC one could proceed in two novel ways: i) start from the finite-dimensional Hamiltonian theory of the discretization and then use the considerable lore of more than a century of rigorous Hamiltonian perturbation methods to arrive at a reduced faithful system. This is facilitated by the adaptation of canonical perturbations methods to noncanonical 150,151 or ii) use multiscale methods (e.g. Ref. 152) on the large but finite GEMPIC-like system, thereby circumventing the usual gyrokinetic expansions by incorporating modern day numerical tools. Neither of these approaches are trivial, and their efficacy remains to be determined.