4.1 The finite element electromagnetic Poisson structure

To start, let us consider electromagnetic fields on the continuous manifold $\varOmega$, using the temporal gauge wherein the electric potential vanishes, ${\phi =0}$. The configuration space for such fields is the set ${Q=\{ {\boldsymbol {A}}~|~ {\boldsymbol {A}}\in \varLambda ^{1}(\varOmega )\}}$ of possible vector potentials, defined as differential ${\text {1-forms}}$ over $\varOmega$. To find a variable conjugate to $ {\boldsymbol {A}}( {\boldsymbol {x}})$, we compute the following variational derivative of the electromagnetic Lagrangian expressed in Gaussian units,

(4.1a,b)\begin{equation} L_\text{EM}=\frac{1}{8{\rm \pi}}\int\mathrm{d}{{\boldsymbol{x}}}\left(\left| {-\frac{1}{c}{\dot{\boldsymbol{A}}}} \right|^{2}- \left| {\mathrm{d}{\boldsymbol{A}}} \right|^{2}\right),\quad {\boldsymbol{Y}}=\frac{\delta L_\text{EM}}{\delta{\dot{\boldsymbol{A}}}}=\frac{{\dot{\boldsymbol{A}}}}{4{\rm \pi} c^{2}}. \end{equation}

Clearly, ${ {\boldsymbol {Y}}\in \varLambda ^{1}(\varOmega )}$ is also a ${\text {1-form}}$ over $\varOmega$ corresponding to negative the electric field, ${ {\boldsymbol {Y}}=- {\boldsymbol {E}}/4{\rm \pi} c}$. As in (3.6), ${\left | {\alpha } \right |^{2}=(\alpha,\alpha )_p}$ in (4.1a) denotes the standard inner product on $\mathbb {R}^{3}$ for ${p=1,2}$.

The full phase space is then given by the cotangent bundle ${T^{*}Q=\{ {\boldsymbol {A}}, {\boldsymbol {Y}}\}}$ whose canonical symplectic structure is defined by the Poisson bracket (Marsden & Weinstein Reference Marsden and Weinstein1982)

(4.2)\begin{equation} \{F,G\}=\int\left(\frac{\delta F}{\delta{\boldsymbol{A}}}\boldsymbol{\cdot}\frac{\delta G}{\delta{\boldsymbol{Y}}}-\frac{\delta G}{\delta{\boldsymbol{A}}}\boldsymbol{\cdot}\frac{\delta F}{\delta{\boldsymbol{Y}}}\right)\mathrm{d}{\boldsymbol{x}}. \end{equation}

Here, ${F[ {\boldsymbol {A}}, {\boldsymbol {Y}}]}$ and ${G[ {\boldsymbol {A}}, {\boldsymbol {Y}}]}$ are arbitrary functionals on ${T^{*}Q}$.

We now map this geometric description of fields on $\varOmega$ to its triangulation $\mathcal {T}_h$. On $\mathcal {T}_h$, the fields $ {\boldsymbol {A}}$ and $ {\boldsymbol {Y}}$ can be defined by their expansion in the corresponding basis for finite element ${\text {1-forms}}$

(4.3)\begin{equation} \left.\begin{gathered} {\boldsymbol{A}}(t,{\boldsymbol{x}})={\boldsymbol{a}}(t)\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{1}({\boldsymbol{x}})\\ {\boldsymbol{Y}}(t,{\boldsymbol{x}})={\boldsymbol{y}}(t)\boldsymbol{\cdot}\mathbb{M}_1^{{-}1}\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{1}({\boldsymbol{x}}). \end{gathered}\right\} \end{equation}

Here, ${ { \boldsymbol {\varLambda }}^{1}}$ is an ${N_1\times 1}$ vector of basis elements and ${ {\boldsymbol {a}}, {\boldsymbol {y}}\in \mathbb {R}^{N_1}}$ denote coefficients, as in (3.2). $ {\boldsymbol {a}}$ and $ {\boldsymbol {y}}$ are identified as dynamical variables by explicitly notating their time dependence. The inverse factor of the 1-form mass matrix $\mathbb {M}_1$ in (4.3) follows from computing the conjugate momentum of $ {\boldsymbol {a}}$ in the following discretization of $L_\textrm {EM}$:

(4.4)\begin{align} L_\text{EM}& =\frac{1}{8{\rm \pi}}\int_{\left| {\mathcal{T}_h} \right|}\mathrm{d}{{\boldsymbol{x}}} \left(\left| {-\frac{1}{c}\dot{{\boldsymbol{a}}}\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{1}} \right|^{2}-\left| {{\boldsymbol{a}} \boldsymbol{\cdot}\mathrm{d}{ \boldsymbol{\varLambda}}^{1}} \right|^{2}\right)\nonumber\\ & =\frac{1}{8{\rm \pi}}\int_{\left| {\mathcal{T}_h} \right|}\mathrm{d}{{\boldsymbol{x}}} \left(\frac{1}{c^{2}}\dot{{\boldsymbol{a}}}\boldsymbol{\cdot}\mathbb{M}_1\boldsymbol{\cdot}\dot{{\boldsymbol{a}}}- {\boldsymbol{a}}\boldsymbol{\cdot}\mathbb{C}^{T}\mathbb{M}_2\mathbb{C}\boldsymbol{\cdot}{\boldsymbol{a}}\right),\quad {\boldsymbol{y}}=\frac{\delta L_\text{EM}}{\delta\dot{{\boldsymbol{a}}}}=\frac{\mathbb{M}_1\dot{{\boldsymbol{a}}}}{4{\rm \pi} c^{2}}, \end{align}

where we have substituted ${\mathrm {d} { \boldsymbol {\varLambda }}^{1}=\mathbb {C}^{T} { \boldsymbol {\varLambda }}^{2}}$ from table 1 and used mass matrices as defined in (3.5). Comparing (4.1a,b) and (4.4) yields the desired expansion for $ {\boldsymbol {Y}}$ in (4.3). (This inverse factor of $\mathbb {M}_1$ will further ensure that the Poisson bracket of (4.9) is canonical.)

To discretize the Poisson bracket of (4.2), we must first express ${\delta F/\delta {\boldsymbol {A}}}$ in terms of ${\partial F/\partial {\boldsymbol {a}}}$ for a discrete variation ${\delta {\boldsymbol {A}}=\delta {\boldsymbol {a}}\boldsymbol {\cdot } { \boldsymbol {\varLambda }}^{1}}$. Since variational derivatives are valued dually to their variations, following Kraus et al. (Reference Kraus, Kormann, Morrison and Sonnendrücker2017), it is appropriate to require that

(4.5)\begin{equation} {\left\langle\frac{\delta F}{\delta{\boldsymbol{A}}},\delta{\boldsymbol{a}}\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{1}\right\rangle}_{L^{2}\varLambda^{1}}={\left\langle\frac{\partial F}{\partial{\boldsymbol{a}}},\delta{\boldsymbol{a}}\right\rangle}_{\mathbb{R}^{N_1}}. \end{equation}

Here, ${\langle \cdot,\cdot \rangle _{L^{2}\varLambda ^{1}}=\int _{\left | {\mathcal {T}_h} \right |}\mathrm {d} {\boldsymbol {x}}(\cdot,\cdot )_{p=1}}$ denotes the ${L^{2}\varLambda ^{1}}$ inner product, and ${\langle \cdot,\cdot \rangle _{\mathbb {R}^{N_1}}}$ the standard inner product on $\mathbb {R}^{N_1}$. In particular, setting $\delta a_j=\delta _{ij}$ for some fixed, arbitrary ${i\in [1,N_1]}$ and using ${(\cdot,\cdot )_p}$ as defined in (3.6), (4.5) yields

(4.6)\begin{equation} \int_{\left| {\mathcal{T}_h} \right|}\mathrm{d}{\boldsymbol{x}}\left(\frac{\delta F}{\delta{\boldsymbol{A}}},\varLambda^{1}_i({\boldsymbol{x}})\right)_{p=1}=\frac{\partial F}{\partial a_i}. \end{equation}

To solve for ${\delta F/\delta {\boldsymbol {A}}}$, we expand ${\delta F/\delta {\boldsymbol {A}}= {\boldsymbol {f}}\boldsymbol {\cdot } { \boldsymbol {\varLambda }}^{1}}$ in the ${ { \boldsymbol {\varLambda }}^{1}}$ basis for some ${ {\boldsymbol {f}}\in \mathbb {R}^{N_1}}$, such that evaluating (4.6) implies ${ {\boldsymbol {f}}=\mathbb {M}_1^{-1}\boldsymbol{\cdot}\partial F/\partial {\boldsymbol {a}}}$. Therefore,

(4.7)\begin{equation} \frac{\delta F}{\delta{\boldsymbol{A}}}=\frac{\partial F}{\partial{\boldsymbol{a}}}\boldsymbol{\cdot}\mathbb{M}_1^{{-}1}\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{1}. \end{equation}

The discrete variation ${\delta {\boldsymbol {Y}}=\delta {\boldsymbol {y}}\boldsymbol {\cdot }\mathbb {M}_1^{-1}\boldsymbol {\cdot } { \boldsymbol {\varLambda }}^{1}}$ establishes a similar result for $ {\boldsymbol {Y}}$, namely

(4.8)\begin{equation} \frac{\delta F}{\delta{\boldsymbol{Y}}}=\frac{\partial F}{\partial{\boldsymbol{y}}}\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{1}. \end{equation}

Thus, to derive the discrete Poisson bracket, we substitute (4.7)–(4.8) into (4.2) and integrate over $\left | {\mathcal {T}_h} \right |$ to find

(4.9)\begin{equation} \{F,G\}=\frac{\partial {F}}{\partial {{\boldsymbol{a}}}}\boldsymbol{\cdot}\frac{\partial {G}}{\partial {{\boldsymbol{y}}}}-\frac{\partial {G}}{\partial {{\boldsymbol{a}}}}\boldsymbol{\cdot}\frac{\partial {F}}{\partial {{\boldsymbol{y}}}}. \end{equation}

True to its continuous counterpart in (4.2), the discrete Poisson bracket of (4.9) is in canonical symplectic form (i.e. Darboux coordinates).

4.2 The finite element Vlasov–Maxwell system

The full Poisson structure of the discrete Vlasov–Maxwell system readily follows from the foregoing analysis of its electromagnetic subsystem. To describe a system of $L$ particles, we let ${ {\boldsymbol {X}}_\ell, {\boldsymbol {P}}_\ell \in \mathbb {R}^{3}}$ ${\ell \in [1,L]}$ denote the position and momentum of the $\ell \textrm {th}$ particle and let ${m_\ell }$ and ${q_\ell }$ denote its mass and charge. Particles may then be characterized by a Klimontovich distribution, ${f( {\boldsymbol {x}}, {\boldsymbol {p}})=\sum _\ell \delta ( {\boldsymbol {x}}- {\boldsymbol {X}}_\ell )\delta ( {\boldsymbol {p}}- {\boldsymbol {P}}_\ell )}$. Particle phase space is defined as usual with a canonical bracket on position ${ {\boldsymbol {X}}_\ell }$ and momentum ${ {\boldsymbol {P}}_\ell }$.

Combining (4.9) with the Poisson bracket for these $L$ particles, therefore, the discrete Vlasov–Maxwell Poisson structure is given by

(4.10)\begin{equation} \{F,G\}=\frac{\partial {F}}{\partial {{\boldsymbol{a}}}}\boldsymbol{\cdot}\frac{\partial {G}}{\partial {{\boldsymbol{y}}}}-\frac{\partial {G}}{\partial {{\boldsymbol{a}}}}\boldsymbol{\cdot}\frac{\partial {F}}{\partial {{\boldsymbol{y}}}}+\sum_{\ell=1}^{L}\left(\frac{\partial F}{\partial{\boldsymbol{X}}_\ell}\boldsymbol{\cdot}\frac{\partial G}{\partial{\boldsymbol{P}}_\ell}-\frac{\partial G}{\partial{\boldsymbol{X}}_\ell}\boldsymbol{\cdot}\frac{\partial F}{\partial{\boldsymbol{P}}_\ell}\right). \end{equation}

Here, $F$ and $G$ are arbitrary functionals on the discrete Poisson manifold ${(M_d,\{\cdot,\cdot \})}$, where each point ${m_d\in M_d}$ is defined by the data

(4.11)\begin{align} m_d& =({\boldsymbol{a}},{\boldsymbol{y}},{\boldsymbol{X}}_1,\dots,{\boldsymbol{X}}_L,{\boldsymbol{P}}_1,\dots,{\boldsymbol{P}}_L)\nonumber\\ & \in\mathbb{R}^{N_1}\times\mathbb{R}^{N_1}\times\mathbb{R}^{3L}\times\mathbb{R}^{3L}. \end{align}

We now define a Hamiltonian ${H_\textrm {VM}=H_\textrm {VM}[ {\boldsymbol {a}}, {\boldsymbol {y}}, {\boldsymbol {X}}_i, {\boldsymbol {P}}_i]}$ on $M_d$, given in Gaussian units by

(4.12)\begin{equation} \left.\begin{gathered} H_\text{VM}=H_\text{EM}+H_\text{Kinetic}\\ \text{where}\quad H_\text{EM}=\frac{1}{8{\rm \pi}}\int_{\left| {\mathcal{T}_h} \right|} \mathrm{d}{\boldsymbol{x}}\left(\left| {-4{\rm \pi} c{\boldsymbol{Y}}} \right|^{2}+\left| {\mathrm{d}{\boldsymbol{A}}} \right|^{2}\right)\\ \qquad\qquad\qquad\qquad\quad=\frac{1}{8{\rm \pi}}\left((4{\rm \pi} c)^{2}{\boldsymbol{y}}\boldsymbol{\cdot}\mathbb{M}_1^{{-}1}\boldsymbol{\cdot} {\boldsymbol{y}}+{\boldsymbol{a}}\boldsymbol{\cdot}\mathbb{C}^{T}\mathbb{M}_2\mathbb{C}\boldsymbol{\cdot}{\boldsymbol{a}}\right)\\ H_\text{Kinetic}=\sum_{\ell=1}^{L}\frac{1}{2m_\ell}\left|{\boldsymbol{P}}_\ell-\frac{q_\ell}{c}{\boldsymbol{A}}({\boldsymbol{X}}_\ell)\right|^{2} \end{gathered}\right\}, \end{equation}

where ${ {\boldsymbol {A}}( {\boldsymbol {X}}_\ell )= {\boldsymbol {a}}\boldsymbol {\cdot } { \boldsymbol {\varLambda }}^{1}( {\boldsymbol {X}}_\ell )}$. In $H_\textrm {EM}$, we have substituted (3.5) and (4.3). The difference in $H_\textrm {Kinetic}$ is taken componentwise, i.e. ${P_{\ell \mu }-(q_\ell /c) {\boldsymbol {a}}\boldsymbol {\cdot } { \boldsymbol {\varLambda }}^{1}( {\boldsymbol {X}}_\ell )_\mu }$ for ${\mu \in \{1,2,3\}}$. Hereafter, we denote components of $ {\boldsymbol {P}}_\ell$ by $P_{\ell \mu }$ and those of $ {\boldsymbol {X}}_\ell$ by $X_\ell ^{\mu }$.

4.3 Gauge structure

We now examine the gauge structure of this discrete Vlasov–Maxwell system and derive its corresponding momentum map. We first note that, because ${\mathbb {C}\mathbb {G}=0}$, $H_\textrm {VM}$ of (4.12) is invariant under any gauge transformation ${\varPhi _{\exp ( {\boldsymbol {s}})}:M_d\rightarrow M_d}$ of the form

(4.13)\begin{equation} \varPhi_{\exp({\boldsymbol{s}})}\left(\begin{matrix} {\boldsymbol{a}}\\ {\boldsymbol{y}}\\ {\boldsymbol{X}}_\ell\\ {\boldsymbol{P}}_\ell \end{matrix}\right)=\left(\begin{matrix} {\boldsymbol{a}}+\mathbb{G}{\boldsymbol{s}}\\ {\boldsymbol{y}}\\ {\boldsymbol{X}}_\ell\\ {\boldsymbol{P}}_\ell+\frac{q_\ell}{c}\mathbb{G}{\boldsymbol{s}}\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{1}({\boldsymbol{X}}_\ell)\end{matrix}\right) \end{equation}

$\forall ~{\ell \in [1,L]}$, ${ {\boldsymbol {s}}\in \mathbb {R}^{N_0}}$. Since ${\varPhi _{\exp ( {\boldsymbol {s}})}\circ \varPhi _{\exp ( {\boldsymbol {t}})}=\varPhi _{\exp ( {\boldsymbol {s}}+ {\boldsymbol {t}})}}$, such transformations form an abelian group. Evidently, they are also generated by vector fields ${X_{ {\boldsymbol {s}}}\in \mathfrak {X}(M_d)}$ of the form

(4.14)\begin{equation} X_{{\boldsymbol{s}}}=\left.\frac{\mathrm{d}}{\mathrm{d}\epsilon}\right|_{\epsilon=0}\varPhi_{\exp(\epsilon{\boldsymbol{s}})}= \mathbb{G}_{jk}s_k\left[\partial_{a_j}+\sum_{\ell=1}^{L}\frac{q_\ell}{c}\varLambda^{1}_j({\boldsymbol{X}}_\ell)_\mu\partial_{P_{\ell\mu}}\right], \end{equation}

expressed using Einstein summation convention.

We may check whether $\varPhi$ of (4.13) is a canonical transformation, that is, whether it preserves the Poisson bracket of (4.10)

(4.15)\begin{equation} \{F,G\}\circ\varPhi_{\exp({\boldsymbol{s}})}=\{F\circ\varPhi_{\exp({\boldsymbol{s}})},G\circ\varPhi_{\exp({\boldsymbol{s}})}\} \end{equation}

$\forall \ {\boldsymbol {s}}$. It suffices to check this condition infinitesimally, i.e. whether

(4.16)\begin{equation} X_{{\boldsymbol{s}}}(\{F,G\})=\{X_{{\boldsymbol{s}}}(F),G\}+\{F,X_{{\boldsymbol{s}}}(G)\}. \end{equation}

After cancelling terms by the equality of mixed partials, verifying (4.16) reduces to checking that

(4.17)\begin{align} 0& =\sum_{\ell=1}^{L}\frac{q_\ell}{c}\mathbb{G}_{jk}s_k\left[\partial_{X_\ell^{\nu}}\varLambda^{1}_j({\boldsymbol{X}}_\ell)_\mu-\partial_{X_\ell^{\mu}}\varLambda^{1}_j({\boldsymbol{X}}_\ell)_\nu\right] (\partial_{P_{\ell\mu}}F)(\partial_{P_{\ell\nu}}G)\nonumber\\ & =\sum_{\ell=1}^{L}\frac{q_\ell}{c}\left[\mathbb{G}{\boldsymbol{s}}\boldsymbol{\cdot}\mathrm{d}{ \boldsymbol{\varLambda}}^{1}({\boldsymbol{X}}_\ell)_{\nu\mu}\right] (\partial_{P_{\ell\mu}}F)(\partial_{P_{\ell\nu}}G). \end{align}

Each term in this sum is seen to vanish, however, since

(4.18)\begin{equation} \mathbb{G}{\boldsymbol{s}}\boldsymbol{\cdot}\mathrm{d}{ \boldsymbol{\varLambda}}^{1}={\boldsymbol{s}}\boldsymbol{\cdot}\mathbb{G}^{T}\mathrm{d} { \boldsymbol{\varLambda}}^{1}={\boldsymbol{s}}\boldsymbol{\cdot}\mathrm{d}(\mathbb{G}^{T}{ \boldsymbol{\varLambda}}^{1})={\boldsymbol{s}}\boldsymbol{\cdot}\mathrm{d}\mathrm{d}{ \boldsymbol{\varLambda}}^{0}=0. \end{equation}

Thus, $X_{ {\boldsymbol {s}}}$ generates a canonical group action, as desired.

We now find the momentum map $\mu$ of this canonical group action by solving for its generating functions. In particular, given any ${ {\boldsymbol {s}}\in \mathbb {R}^{N_0}}$, we seek a generating function $\mu _{ {\boldsymbol {s}}}$ such that ${X_{ {\boldsymbol {s}}}=\{\cdot,\mu _{ {\boldsymbol {s}}}\}}$ as derivations, as in (2.2). By comparing (4.10) and (4.14), we see that ${X_{ {\boldsymbol {s}}}=\{\cdot,\mu _{ {\boldsymbol {s}}}\}}$ holds if and only if

(4.19)\begin{equation} \left.\begin{gathered} \frac{\partial\mu_{{\boldsymbol{s}}}}{\partial{\boldsymbol{a}}}=0\\ \frac{\partial\mu_{{\boldsymbol{s}}}}{\partial{\boldsymbol{y}}}=\mathbb{G}{\boldsymbol{s}}\\ \frac{\partial\mu_{{\boldsymbol{s}}}}{\partial{\boldsymbol{X}}_\ell}={-}\frac{q_\ell}{c}\mathbb{G}{\boldsymbol{s}}\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{1}({\boldsymbol{X}}_\ell)\\ \frac{\partial\mu_{{\boldsymbol{s}}}}{\partial{\boldsymbol{P}}_\ell}=0\quad \forall\ {\ell\in[1,L]}. \end{gathered}\right\} \end{equation}

Since ${\mathbb {G}^{T} { \boldsymbol {\varLambda }}^{1}=\mathrm {d} { \boldsymbol {\varLambda }}^{0}}$ is an exact form, this linear system of partial differential equations is readily solved by

(4.20)\begin{equation} \mu_{{\boldsymbol{s}}}=\mathbb{G}{\boldsymbol{s}}\boldsymbol{\cdot}{\boldsymbol{y}}-\sum_{\ell=1}^{L}\frac{q_\ell}{c}{\boldsymbol{s}}\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{0}({\boldsymbol{X}}_\ell). \end{equation}

The momentum map $\mu$ characterizing all generating functions $\{\mu _{ {\boldsymbol {s}}}\}$ is defined by requiring that ${\mu \boldsymbol {\cdot } {\boldsymbol {s}}=\mu _{ {\boldsymbol {s}}}}\ \forall \ {\boldsymbol {s}}$. Therefore, $\mu$ is given by

(4.21)\begin{equation} \mu=\mathbb{G}^{T}{\boldsymbol{y}}-\sum_{\ell=1}^{L}\frac{q_\ell}{c}{ \boldsymbol{\varLambda}}^{0}({\boldsymbol{X}}_\ell). \end{equation}

Setting ${\mu =0}$, (4.21) is a discrete form of Gauss’ law,

(4.22)\begin{equation} 0=(\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{E}}-4{\rm \pi}\rho)/4{\rm \pi} c, \end{equation}

where ${ {\boldsymbol {E}}=-4{\rm \pi} c {\boldsymbol {Y}}}$ and ${\rho =\sum _\ell q_\ell { \boldsymbol {\varLambda }}^{0}( {\boldsymbol {X}}_\ell )}$.

Any non-zero value of $\mu$ indicates that the divergence of the electric field is not entirely accounted for by the dynamical particles labelled ${\ell \in [1,L]}$. Since ${\dot {\mu }=0}$, such a non-zero $\mu$ acts as a fixed, external, non-dynamical background charge that persists throughout a simulation, and in a manner that remains entirely structure preserving. In particular, we may regard $\mu$ as representing an external charge density,

(4.23)\begin{equation} \mu\sim\rho_\text{ext}/c. \end{equation}

As we shall demonstrate, (4.23) can be useful in establishing precise initial conditions in a PIC simulation.

5.1 Continuous-time equations of motion

Let us now derive equations of motion via the Hamiltonian of (4.12) and the Poisson bracket of (4.10). We find

(5.1)\begin{equation} \left.\begin{gathered} \dot{X}_\ell^{\mu}=\{X_\ell^{\mu},H_\text{VM}\}=\frac{1}{m_\ell}\left(P_{\ell\mu}-\frac{q_\ell}{c}{\boldsymbol{A}}({\boldsymbol{X}}_\ell)_\mu\right)\\ \dot{P}_{\ell\mu}=\{P_{\ell\mu},H_\text{VM}\}=\frac{q_\ell}{c}\dot{X}_\ell^{\nu}\partial_{X_\ell^{\mu}}{\boldsymbol{A}}({\boldsymbol{X}}_\ell)_\nu\\ \dot{{\boldsymbol{a}}}=\{{\boldsymbol{a}},H_\text{VM}\}=4{\rm \pi} c^{2}\mathbb{M}_1^{{-}1}\boldsymbol{\cdot}{\boldsymbol{y}}\\ \dot{{\boldsymbol{y}}}=\{{\boldsymbol{y}},H_\text{VM}\}={-}\frac{1}{4{\rm \pi}}\mathbb{C}^{T}\mathbb{M}_2\mathbb{C} \boldsymbol{\cdot}{\boldsymbol{a}}+\sum_{\ell=1}^{L}\frac{q_\ell}{c}\dot{X}_\ell^{\mu}{ \boldsymbol{\varLambda}}^{1}({\boldsymbol{X}}_\ell)_\mu \end{gathered}\right\}, \end{equation}

where ${ {\boldsymbol {A}}( {\boldsymbol {X}}_\ell )_\mu = {\boldsymbol {a}}\boldsymbol {\cdot } { \boldsymbol {\varLambda }}^{1}( {\boldsymbol {X}}_\ell )_\mu }$ using the component notation of (3.3).

We remark that these equations of motion allow us to reexpress the conservation of the momentum map, ${\dot {\mu }=0}$, as a discrete, local charge conservation law in conservative form. In particular, we may take the time derivative $\dot {\mu }$ of (4.21) and substitute $\dot { {\boldsymbol {y}}}$ of (5.1), noting that ${\mathbb {C}\mathbb {G}=0}$ to find

(5.2)\begin{equation} (\mathbb{G}^{T}{\boldsymbol{j}}-\dot{\rho})/c=0, \end{equation}

where ${\rho =\sum _\ell q_\ell { \boldsymbol {\varLambda }}^{0}( {\boldsymbol {X}}_\ell )}$ and ${ {\boldsymbol {j}}=\sum _\ell q_\ell \dot {X}^{\mu }_\ell { \boldsymbol {\varLambda }}^{1}( {\boldsymbol {X}}_\ell )_\mu }$. As in (4.21), we note a correspondence between the discrete and continuous operators ${\mathbb {G}^{T}\sim (-\boldsymbol {\nabla }\boldsymbol {\cdot })}$. Thus, (5.2) is a discrete equivalent of the charge conservation law ${\partial _t{\rho }+\boldsymbol {\nabla }\boldsymbol {\cdot } {\boldsymbol {j}}=0}$, as seen in Kraus et al. (Reference Kraus, Kormann, Morrison and Sonnendrücker2017).

Equation (5.1) is sometimes referred to as a semi-discrete system, since it describes a discretely defined system evolving in continuous time. We now proceed to a fully discrete, algorithmic system by defining a GCSM.

5.2 Discrete-time equations of motion via a gauge-compatible splitting

Using the Vlasov–Maxwell splitting discovered in He et al. (Reference He, Qin, Sun, Xiao, Zhang and Liu2015) and adapted to canonical coordinates in Glasser & Qin (Reference Glasser and Qin2020), we split $H_\textrm {VM}$ of (4.12) into five sub-Hamiltonians, as follows:

(5.3)\begin{equation} \left.\begin{gathered} H_\text{VM}=H_{{\boldsymbol{A}}}+H_{{\boldsymbol{Y}}}+H^{x}_\text{Kinetic}+H^{y}_\text{Kinetic}+H^{z}_\text{Kinetic}\\ \text{where}\quad H_{{\boldsymbol{A}}}=\frac{1}{8{\rm \pi}}{\boldsymbol{a}}\boldsymbol{\cdot}\mathbb{C}^{T}\mathbb{M}_2\mathbb{C}\boldsymbol{\cdot}{\boldsymbol{a}}\\ H_{{\boldsymbol{Y}}}=\frac{1}{8{\rm \pi}}(4{\rm \pi} c)^{2}{\boldsymbol{y}}\boldsymbol{\cdot}\mathbb{M}_1^{{-}1}\boldsymbol{\cdot}{\boldsymbol{y}}\\ H^{\alpha}_\text{Kinetic}=\sum_{\ell=1}^{L}\frac{1}{2m_\ell} \left(P_{\ell\alpha}-\frac{q_\ell}{c}{\boldsymbol{A}}({\boldsymbol{X}}_\ell)_\alpha\right)^{2} \end{gathered}\right\} \end{equation}

$\forall \ {\alpha \in \{x,y,z\}}$. We immediately observe that each sub-Hamiltonian remains invariant under the gauge transformation $\varPhi _{\exp ( {\boldsymbol {s}})}$ defined in (4.13).

Let us now examine the equations of motion of each subsystem, omitting equations for the subsystems’ static degrees of freedom

(5.4)\begin{equation} \left.\begin{gathered} H_{{\boldsymbol{A}}}:\ \dot{{\boldsymbol{y}}}={-}\frac{1}{4{\rm \pi}}\mathbb{C}^{T}\mathbb{M}_2\mathbb{C}\boldsymbol{\cdot}{\boldsymbol{a}}\\ H_{{\boldsymbol{Y}}}:\ \dot{{\boldsymbol{a}}}=4{\rm \pi} c^{2}\mathbb{M}_1^{{-}1}\boldsymbol{\cdot}{\boldsymbol{y}}\\ H_{\text{Kinetic}}^{\alpha}:\ \begin{cases} \dot{X}_\ell^{\alpha}=\dfrac{1}{m_\ell}\left(P_{\ell\alpha}-\dfrac{q_\ell}{c}{\boldsymbol{A}}({\boldsymbol{X}}_\ell)_\alpha\right)\\ \dot{P}_{\ell\mu}=\dfrac{q_\ell}{c}\dot{X}_\ell^{\alpha}\partial_{X_\ell^{\mu}}{\boldsymbol{A}}({\boldsymbol{X}}_\ell)_\alpha\\ \dot{{\boldsymbol{y}}}=\sum\limits_{\ell=1}^{L}\dfrac{q_\ell}{c}\dot{X}_\ell^{\alpha}{ \boldsymbol{\varLambda}}^{1}({\boldsymbol{X}}_\ell)_\alpha. \end{cases} \end{gathered}\!\!\!\!\!\!\!\right\}\end{equation}

To clarify the notation in (5.4), we emphasize that, in the $H_{\textrm {Kinetic}}^{\alpha }$ subsystem, ${\dot {X}_\ell ^{\mu }=0}$ for ${\mu \neq \alpha }$. (Here, $\alpha$ is regarded as fixed while $\mu$ ranges over all $\{{x,y,z}\}$ indices.) Thus, the equations of motion for ${X}_\ell ^{\mu \neq \alpha }$ are omitted above.

Furthermore, it follows from a simple calculation that ${\ddot {X}_\ell ^{\alpha }=0}$ in $H_{\textrm {Kinetic}}^{\alpha }$ so that ${\dot {X}_\ell ^{\alpha }}$ is constant during the evolution of each subsystem. As a result, all subsystems above are exactly integrable. Equation (5.4) therefore defines a GCSM.

More concretely, an evolution over the timestep $[t,t+\varDelta ]$ in each subsystem is fully specified by

(5.5)\begin{equation} \left.\begin{gathered} H_{{\boldsymbol{A}}}:\ {\boldsymbol{y}}(t+\varDelta)={\boldsymbol{y}}(t)-\frac{\varDelta}{4{\rm \pi}}\mathbb{C}^{T}\mathbb{M}_2\mathbb{C}\boldsymbol{\cdot}{\boldsymbol{a}}(t)\\ H_{{\boldsymbol{Y}}}:\ {\boldsymbol{a}}(t+\varDelta)={\boldsymbol{a}}(t)+\varDelta 4{\rm \pi} c^{2}\mathbb{M}_1^{{-}1}\boldsymbol{\cdot}{\boldsymbol{y}}(t)\\ H_{\text{Kinetic}}^{\alpha}:\ \begin{cases} X_\ell^{\alpha}(t+\delta)=X_\ell^{\alpha}(t)+\dfrac{\delta}{m_\ell} \left(P_{\ell\alpha}(t)-\dfrac{q_\ell}{c}{\boldsymbol{a}}(t)\cdot{ \boldsymbol{\varLambda}}^{1}({\boldsymbol{X}}_\ell(t))_\alpha\right)\\ P_{\ell\mu}(t+\varDelta)=P_{\ell\mu}(t)+\dfrac{q_\ell}{c}\dot{X}_\ell^{\alpha}(t){\boldsymbol{a}}(t) \boldsymbol{\cdot}\displaystyle\mathop{\int}\limits_t^{t+\varDelta}\mathrm{d} t'\partial_{X_\ell^{\mu}(t')} { \boldsymbol{\varLambda}}^{1}({\boldsymbol{X}}_\ell(t'))_\alpha\\ {\boldsymbol{y}}(t+\varDelta)={\boldsymbol{y}}(t)+\sum\limits_{\ell=1}^{L}\dfrac{q_\ell}{c} \dot{X}_\ell^{\alpha}(t)\displaystyle\mathop{\int}\limits_t^{t+\varDelta} \mathrm{d} t'{ \boldsymbol{\varLambda}}^{1}({\boldsymbol{X}}_\ell(t'))_\alpha. \end{cases} \end{gathered}\!\!\!\!\!\!\!\right\} \end{equation}

In (5.5), $t$ is a fixed initial time and ${\delta \in [0,\varDelta ]}$ parametrizes the particle trajectory ${ {\boldsymbol {X}}_\ell (t)\rightarrow {\boldsymbol {X}}_\ell (t+\varDelta )}$ during one timestep of the ${H_{\textrm {Kinetic}}^{\alpha }}$ subsystem, which forms a straight line segment in the $\hat {\alpha }$ direction. Since $ { \boldsymbol {\varLambda }}^{1}$ is comprised of piecewise polynomial differential $\textrm {1-forms}$, $ { \boldsymbol {\varLambda }}^{1}$ and its derivatives are integrable in closed form along the straight path ${ {\boldsymbol {X}}_\ell (t+\delta )}$. Thus, (5.5) defines a symplectic algorithm – specifically a GCSM – that can be computed exactly.

6.1 Landau damping

Using Whitney form finite elements, a 650-cell domain $\mathcal {T}_h$ with periodic boundaries is constructed. Each cell is assigned width ${w_x=2.4\times 10^{-2}\ \textrm {cm}}$, and ${26\times 10^{6}}$ electrons are simulated (40 000 per cell, when unperturbed). With electron temperature at ${T_e=5}$ keV, the set-up has Debye length ${\lambda _D=1.0\textrm { cm}}$ and plasma frequency ${\omega _p=3.0\times 10^{9}\ \textrm {rad}\ \textrm {s}^{-1}}$, roughly mirroring physical parameters of Xiao et al. (Reference Xiao, Qin, Liu, He, Zhang and Sun2015).

We assume the electric field to have initial perturbation at ${t=0}$

(6.1)\begin{equation} {\boldsymbol{Y}}={-}\frac{{\boldsymbol{E}}}{4{\rm \pi} c}={-}\frac{E_0\cos(kx)}{4{\rm \pi} c}\hat{{\boldsymbol{x}}}, \end{equation}

where ${E_0=1.2\ \textrm {statV}\ \textrm {cm}^{-1}}$ and ${k\lambda _D=0.8}$. To construct this perturbation, we first project the continuous 1-form field $ {\boldsymbol {Y}}_\textrm {cont}$ onto its Whitney form approximation, i.e. ${ {\boldsymbol {Y}}_\textrm {cont}\xrightarrow {{\rm \pi} _h} {\boldsymbol {Y}}= {\boldsymbol {y}}\boldsymbol {\cdot }\mathbb {M}_1^{-1}\boldsymbol {\cdot } { \boldsymbol {\varLambda }}^{1}}$. In particular, we solve for $ {\boldsymbol {y}}$ in (4.3) such that the integrals of ${ {\boldsymbol {Y}}_\textrm {cont}}$ and $ {\boldsymbol {Y}}$ agree on edges of the discretized domain. This procedure yields a sinusoidal $\mu _\textrm {field}=\mathbb {G}^{T} {\boldsymbol {y}}$ as depicted in figure 2.

Figure 2. The terms of (4.21) are plotted over the simulation domain at time ${t=0}$, characterizing initial conditions by the momentum map ${\mu =\mu _\textrm {field}+\mu _\textrm {particle}}$.

From the particle side, electron momenta are initialized in phase space by randomly selecting their velocities $\dot { {\boldsymbol {X}}}_\ell$ from a Maxwellian distribution of temperature ${T_e=5}$ keV. Taking ${ {\boldsymbol {A}}=0}$ at ${t=0}$, initial canonical momenta are therefore given by ${ {\boldsymbol {P}}_\ell =m_\ell \dot { {\boldsymbol {X}}} _\ell }$. Electron positions are then initialized to be consistent with Gauss’ law. More precisely, we demand that particle positions, in combination with the electric field perturbation, ensure the constancy of the total momentum map $\mu$. Following (4.23), a non-zero constant momentum map can be understood as a fixed, homogeneous, ion background charge, ${\mu \sim \rho _\textrm {ext}/c}$. Indeed, in a Landau damping simulation that takes only electrons to be dynamical, such a constant background charge constitutes the desired Gauss’ law remainder.

We therefore optimize electron positions $\{ {\boldsymbol {X}}_\ell \}$ to ensure ${\mu _\textrm {particle}=\sum _\ell ({\left | {e} \right |}/{c}) { \boldsymbol {\varLambda }}^{0}( {\boldsymbol {X}}_\ell )}$ closely satisfies

(6.2)\begin{equation} \mu=\mu_\text{field}+\mu_\text{particle}\approx N_\text{ppc}\left| {e} \right|/c, \end{equation}

where ${\mu _\textrm {field}=\mathbb {G}^{T} {\boldsymbol {y}}}$ and ${N_\textrm {ppc}=40\,000}$ denotes the number of particles per cell. In this way, the positive background charge, characterized by $\mu$, is homogeneous across the simulation domain and enforces quasineutrality with the dynamical electrons.

This optimization of ${\mu _\textrm {particle}}$ is carried out in two stages. First, electrons are randomly selected via rejection sampling (von Neumann Reference von Neumann1951) from the distribution

(6.3)\begin{equation} n_e(x)=n_0\left[1+\frac{kE_0\sin(kx)}{4{\rm \pi}\left| {e} \right|n_0}\right], \end{equation}

which satisfies ${\boldsymbol {\nabla }\boldsymbol {\cdot } {\boldsymbol {E}}=4{\rm \pi} \left | {e} \right |(n_0-n_e)}$. The electron positions ${\{ {\boldsymbol {X}}_\ell \}}$ are then further optimized using Nesterov accelerated gradient descent (Nesterov Reference Nesterov1983) to minimize the objective function

(6.4)\begin{equation} \operatorname{arg\,min}\limits_{\{{\boldsymbol{X}}_\ell\}}\left|\mu-\frac{N_\text{ppc}\left| {e} \right|}{c}\right|^{2}. \end{equation}

The resulting momentum map, plotted in red in figure 2, thus defines a background charge that is homogeneous to a high degree of accuracy.

Note that an alternative initialization, undertaken for example by Kraus et al. (Reference Kraus, Kormann, Morrison and Sonnendrücker2017), reverses the above procedure by randomly generating electron positions first and then solving for $ {\boldsymbol {y}}$ to satisfy Gauss’ law. This alternative is more straightforward computationally than the procedure described above, but it may afford less precision in the specification of the initial electric field perturbation, whenever such precision is desirable.

With initialized fields and electrons, the simulation is then evolved using a first-order Lie–Trotter splitting (Trotter Reference Trotter1959) derived from (5.5), in particular,

(6.5)\begin{equation} {\boldsymbol{u}}(t+\varDelta)=\exp(\Delta H_\text{Kinetic}^{x})\exp(\Delta H_{{\boldsymbol{Y}}})\exp(\Delta H_{{\boldsymbol{A}}}){\boldsymbol{u}}(t), \end{equation}

where ${ {\boldsymbol {u}}(t)}$ denotes the simulation state at time $t$ – i.e. ${ {\boldsymbol {u}}=m_d\in M_d}$, a point in phase space as defined in (4.11).

In figure 3, we first examine the evolution of the momentum map throughout the simulation domain. We note that, while $\mu _\textrm {particle}$ (multicolour) exhibits an oscillation and decay consistent with Landau damping, $\mu$ (grey) remains constant over time, consistent with the conservation of the momentum map.

Figure 3. With their initial conditions as depicted in figure 2, the total momentum map $\mu$ (grey) is compared with $\mu _\textrm {particle}$ (multicolour) as the two functions evolve over time. Whereas $\mu _\textrm {particle}$ exhibits a decaying sinusoid consistent with Landau damping, $\mu$ remains constant to machine precision. The momentum map $\mu$ constitutes a physical representation of the fixed positive background charge implicit in the simulation.

To compare this simulation with theory, the evolution of the (normalized) electric field is plotted in figure 4 (which may be compared with figure 2 of Xiao et al. Reference Xiao, Qin, Liu, He, Zhang and Sun2015). The results agree with a theoretical expectation of (i) electrostatic Langmuir wave oscillation at a frequency ${\omega _p=3.0\times 10^{9}\ \textrm {rad}\ \textrm {s}^{-1}}$, and (ii) Landau damping at a decay rate ${\omega _i={\omega _p}/{\kappa ^{3}}\sqrt {({{\rm \pi} }/{8})}\exp (-{(1+3\kappa ^{2})}/{2\kappa ^{2}})=3.9\times 10^{8} \textrm {rad}\ \textrm {s}^{-1}}$, where ${\kappa =k\lambda _D}$. Furthermore, as is characteristic of a symplectic algorithm, the error in the total energy, measured by $H_\textrm {VM}$ of (4.12), is well bounded throughout the simulation. This error is plotted in figure 5.

Figure 4. The evolution of an electrostatic wave over time is simulated with a first-order Lie–Trotter splitting (Trotter Reference Trotter1959) of (5.5). The blue time series denotes the (normalized) log modulus of the electric field ${ {\boldsymbol {E}}=-4{\rm \pi} c {\boldsymbol {Y}}}$, where ${\left | { {\boldsymbol {E}}} \right |}$ is computed over the simulation domain by the ${L^{2}\varLambda ^{1}}$ norm. The theoretical Landau damping rate of the wave in a Maxwellian plasma is depicted as a red line, decaying at a rate of ${\omega _i={\omega _p}/{\kappa ^{3}}\sqrt {({{\rm \pi} }/{8})}\exp (-{(1+3\kappa ^{2})}/{2\kappa ^{2}})}$ for ${\kappa =k\lambda _D}$.

Figure 5. The log error in the total energy of a first-order Lie–Trotter splitting (Trotter Reference Trotter1959) Landau damping simulation.

Having demonstrated the canonical finite element formalism's efficacy in simulating an electrostatic problem, we next consider the electromagnetic simulation of the Weibel instability. Before that, however, we introduce the ‘canonical 1X2P phase space’ in which we will conduct our simulation.

6.2 Canonical 1X2P phase space

Let us consider the restriction of phase space to one spatial dimension and two dimensions of canonical momentum. Despite its computational efficiency, the 1X2P setting is capable of simulating a number of non-trivial electromagnetic problems. Our development here of 1X2P phase space essentially adapts for canonical coordinates the techniques of Crouseilles et al. (Reference Crouseilles, Einkemmer and Faou2015) and Kraus et al. (Reference Kraus, Kormann, Morrison and Sonnendrücker2017).

To characterize 1X2P, we must reflect the lack of $y$ dependence in the finite element expansions of our fields. In this setting, we therefore define $ {\boldsymbol {A}}$ and $ {\boldsymbol {Y}}$ by analogy with (4.3) as

(6.6)\begin{equation} \left.\begin{gathered} {\boldsymbol{A}}={\boldsymbol{a}}_x\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{1}(x)+{\boldsymbol{a}}_y\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{0}(x)\\ {\boldsymbol{Y}}={\boldsymbol{y}}_x\boldsymbol{\cdot}\mathbb{M}_1^{{-}1}\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{1}(x)+{\boldsymbol{y}}_y\boldsymbol{\cdot} \mathbb{M}_0^{{-}1}\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{0}(x). \end{gathered}\right\} \end{equation}

Here, ${ {\boldsymbol {a}}_x}$ denotes a vector of coefficients that pair only with the finite element 1-form basis ${ { \boldsymbol {\varLambda }}^{1}(x)}$ – a basis in 1X2P which has only $x$-components and whose spatial dependence is only one-dimensional. We further note that $ {\boldsymbol {a}}_y$ is defined as the coefficients of a 0-form. Since there is no way to associate 1-forms with $\hat {y}$-directed edges in the 1X2P setting, 0-forms are the most natural representation of the $y$ components of $ {\boldsymbol {A}}$. This becomes especially clear when examining finite elements restricted to have spatial $x$-dependence. The expansion of $ {\boldsymbol {Y}}$ in (6.6) follows similarly, with appropriate mass matrices computed by integrals over the one-dimensional domain.

Turning now to the characterization of particle phase space, we retain two components of the canonical momentum, $P_{\ell x}$ and $P_{\ell y}$. Only one component of spatial dependence is retained, which we shall denote $X_\ell$.

The appropriate Hamiltonian is then computed by restricting (4.12) as follows:

(6.7)\begin{align} H_\text{1X2P}& =\frac{1}{8{\rm \pi}}\left((4{\rm \pi} c)^{2}\left[{\boldsymbol{y}}_x\boldsymbol{\cdot}\mathbb{M}_1^{{-}1}\boldsymbol{\cdot}{\boldsymbol{y}}_x+{\boldsymbol{y}}_y \boldsymbol{\cdot}\mathbb{M}_0^{{-}1}\boldsymbol{\cdot}{\boldsymbol{y}}_y\right]+{\boldsymbol{a}}_y\boldsymbol{\cdot}\mathbb{G}^{T}\mathbb{M}_1\mathbb{G}\boldsymbol{\cdot}{\boldsymbol{a}}_y\right) \nonumber\\ & \quad +\sum_{\ell=1}^{L}\frac{1}{2m_\ell}\left(\left|P_{\ell x}-\frac{q_\ell}{c}{\boldsymbol{a}}_x\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{1}(X_\ell)\right|^{2}+\left|P_{\ell y}-\frac{q_\ell}{c}{\boldsymbol{a}}_y\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{0}(X_\ell)\right|^{2}\right). \end{align}

We note that the term ${ {\boldsymbol {a}}_x\boldsymbol {\cdot }\mathbb {C}^{T}\mathbb {M}_2\mathbb {C}\boldsymbol {\cdot } {\boldsymbol {a}}_x}$ is absent; while it would ordinarily arise from the magnetic energy $|\mathrm {d} {\boldsymbol {A}}|^{2}$, $\mathrm {d}$ annihilates 1-forms in the 1X2P setting. Indeed, the magnetic field in the 1X2P formalism is given simply by ${ {\boldsymbol {B}}=\mathrm {d} {\boldsymbol {A}}= {\boldsymbol {a}}_y\boldsymbol {\cdot }\mathrm {d} { \boldsymbol {\varLambda }}^{0}=\mathbb {G} {\boldsymbol {a}}_y\boldsymbol {\cdot } { \boldsymbol {\varLambda }}^{1}}$.

The Poisson bracket of the restricted phase space follows by simply omitting the inapplicable terms of (4.10). In particular

(6.8)\begin{equation} \{F,G\}_\text{1X2P}=\left(\frac{\partial {F}}{\partial {{\boldsymbol{a}}_x}}\boldsymbol{\cdot}\frac{\partial {G}}{\partial {{\boldsymbol{y}}_x}}+\frac{\partial {F}}{\partial {{\boldsymbol{a}}_y}} \boldsymbol{\cdot}\frac{\partial {G}}{\partial {{\boldsymbol{y}}_y}}+\sum_{\ell=1}^{L}\frac{\partial F}{\partial X_\ell}\frac{\partial G}{\partial P_{\ell x}}\right)-(F\leftrightarrow G). \end{equation}

We observe that, despite its simplicity, the 1X2P phase space retains the gauge structure of the original problem. In particular we note its gauge symmetry

(6.9)\begin{equation} \varPhi_{\exp({\boldsymbol{s}})}\left(\begin{matrix} {\boldsymbol{a}}_x\\ P_{\ell x} \end{matrix}\right)=\left(\begin{matrix} {\boldsymbol{a}}_x+\mathbb{G}{\boldsymbol{s}}_x\\ P_{\ell x}+\dfrac{q_\ell}{c}\mathbb{G}{\boldsymbol{s}}_x\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{1}(X_\ell)_x\end{matrix}\right), \end{equation}

and the resulting momentum map

(6.10)\begin{equation} \mu_\text{1X2P}=\mathbb{G}^{T}{\boldsymbol{y}}_x-\sum_{\ell=1}^{L}\frac{q_\ell}{c}{ \boldsymbol{\varLambda}}^{0}(X_\ell). \end{equation}

Finally, we derive the equations of motion in the 1X2P setting, which are straightforward reexpressions of their full phase space counterparts in (5.1)

(6.11)\begin{equation} \left.\begin{gathered} \dot{X}_\ell=\{X_\ell,H_\text{1X2P}\}=\frac{1}{m_\ell}\Big(P_{\ell x}-\frac{q_\ell}{c}{\boldsymbol{a}}_x\cdot{ \boldsymbol{\varLambda}}^{1}(X_\ell)\Big)\\ \dot{P}_{\ell x}=\{P_{\ell x},H_\text{1X2P}\}=\frac{q_\ell}{c}\left[\dot{X}_\ell{\boldsymbol{a}}_x\cdot\left(\partial_{X_\ell} { \boldsymbol{\varLambda}}^{1}(X_\ell)\right)+\dot{Y}_\ell{\boldsymbol{a}}_y\boldsymbol{\cdot} \left(\partial_{X_\ell}{ \boldsymbol{\varLambda}}^{0}(X_\ell)\right)\right] \\ \dot{{\boldsymbol{a}}}_x=\{{\boldsymbol{a}}_x,H_\text{1X2P}\}=4{\rm \pi} c^{2}\mathbb{M}_1^{{-}1}\boldsymbol{\cdot}{\boldsymbol{y}}_x\\ \dot{{\boldsymbol{a}}}_y=\{{\boldsymbol{a}}_y,H_\text{1X2P}\}=4{\rm \pi} c^{2}\mathbb{M}_0^{{-}1}\boldsymbol{\cdot}{\boldsymbol{y}}_y\\ \dot{{\boldsymbol{y}}}_x=\{{\boldsymbol{y}}_x,H_\text{1X2P}\}=\sum_{\ell=1}^{L}\frac{q_\ell}{c}\dot{X}_\ell{ \boldsymbol{\varLambda}}^{1}(X_\ell)\\ \dot{{\boldsymbol{y}}}_y=\{{\boldsymbol{y}}_y,H_\text{1X2P}\}=\sum_{\ell=1}^{L}\frac{q_\ell}{c} \dot{Y}_\ell{ \boldsymbol{\varLambda}}^{0}(X_\ell)-\frac{1}{4{\rm \pi}}\mathbb{G}^{T}\mathbb{M}_1\mathbb{G}\boldsymbol{\cdot}{\boldsymbol{a}}_y. \end{gathered}\right\} \end{equation}

Here, with an abuse of notation that ignores the non-existence of $Y_\ell$, we define

(6.12)\begin{equation} \dot{Y}_\ell=\frac{1}{m_\ell}\left(P_{\ell y}-{\boldsymbol{a}}_y\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{0}(X_\ell)\right). \end{equation}

We observe that $P_{\ell y}$ is conserved in (6.11), which can be viewed as a consequence of the independence from $Y_\ell$ of the 1X2P formulation.

We further note that a gauge-compatible splitting of $H_\textrm {1X2P}$ is readily defined by the following four sub-Hamiltonians, in agreement with (5.3)

(6.13)\begin{equation} \left.\begin{gathered} H_\text{1X2P}=H_{{\boldsymbol{A}}}+H_{{\boldsymbol{Y}}}+H^{x}_\text{Kinetic}+H^{y}_\text{Kinetic}\\ \text{where}\quad H_{{\boldsymbol{A}}}=\frac{1}{8{\rm \pi}}{\boldsymbol{a}}_y\boldsymbol{\cdot}\mathbb{G}^{T}\mathbb{M}_1\mathbb{G}\boldsymbol{\cdot}{\boldsymbol{a}}_y\\ H_{{\boldsymbol{Y}}}=\frac{1}{8{\rm \pi}}\left((4{\rm \pi} c)^{2}\left[{\boldsymbol{y}}_x\boldsymbol{\cdot}\mathbb{M}_1^{{-}1} \boldsymbol{\cdot}{\boldsymbol{y}}_x+{\boldsymbol{y}}_y\boldsymbol{\cdot}\mathbb{M}_0^{{-}1}\boldsymbol{\cdot}{\boldsymbol{y}}_y\right]\right)\\ H^{x}_\text{Kinetic}=\sum_{\ell=1}^{L}\frac{1}{2m_\ell} \left|P_{\ell x}-\frac{q_\ell}{c}{\boldsymbol{a}}_x\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{1}(X_\ell)\right|^{2}\\ H^{y}_\text{Kinetic}=\sum_{\ell=1}^{L}\frac{1}{2m_\ell} \left|P_{\ell y}-\frac{q_\ell}{c}{\boldsymbol{a}}_y\boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{0}(X_\ell)\right|^{2}. \end{gathered}\right\} \end{equation}

Using the Poisson bracket ${\{\cdot,\cdot \}_\textrm {1X2P}}$ of (6.8), we thereby derive the following subsystem equations of motion:

(6.14)\begin{align} \begin{array}{cc} H_{{\boldsymbol{A}}}:\begin{cases} \dot{{\boldsymbol{y}}}_y={-}\dfrac{1}{4{\rm \pi}}\mathbb{G}^{T}\mathbb{M}_1\mathbb{G}\boldsymbol{\cdot}{\boldsymbol{a}}_y \end{cases} & H_{{\boldsymbol{Y}}}: \begin{cases} \dot{{\boldsymbol{a}}}_x=4{\rm \pi} c^{2}\mathbb{M}_1^{{-}1}\boldsymbol{\cdot}{\boldsymbol{y}}_x\\ \dot{{\boldsymbol{a}}}_y=4{\rm \pi} c^{2}\mathbb{M}_0^{{-}1}\boldsymbol{\cdot}{\boldsymbol{y}}_y \end{cases}\\ H_{\text{Kinetic}}^{x}:\begin{cases} \dot{X}_\ell=\dfrac{1}{m_\ell}\left(P_{\ell x}-\dfrac{q_\ell}{c}{\boldsymbol{a}}_x \boldsymbol{\cdot}{ \boldsymbol{\varLambda}}^{1}({\boldsymbol{X}}_\ell)\right)\\ \dot{P}_{\ell x}=\dfrac{q_\ell}{c}\dot{X}_\ell{\boldsymbol{a}}_x\boldsymbol{\cdot} \left(\partial_{X_\ell}{ \boldsymbol{\varLambda}}^{1}(X_\ell)\right)\\ \dot{{\boldsymbol{y}}}_x=\mathop{\sum}\limits_{\ell=1}^{L}\dfrac{q_\ell}{c}\dot{X}_\ell{ \boldsymbol{\varLambda}}^{1}(X_\ell) \end{cases} & H_{\text{Kinetic}}^{y}: \begin{cases} \dot{P}_{\ell x}=\dfrac{q_\ell}{c}\dot{Y}_\ell{\boldsymbol{a}}_y\boldsymbol{\cdot} \left(\partial_{X_\ell}{ \boldsymbol{\varLambda}}^{0}(X_\ell)\right)\\ \dot{{\boldsymbol{y}}}_y=\sum\limits_{\ell=1}^{L}\dfrac{q_\ell}{c}\dot{Y}_\ell{ \boldsymbol{\varLambda}}^{0}(X_\ell). \end{cases} \end{array} \end{align}

Since these equations of motion are – as (5.5) – exactly solvable, and since the sub-Hamiltonians of (6.13) preserve the gauge symmetry of (6.9), it is clear that (6.14) defines a GCSM that exactly preserves ${\mu _\textrm {1X2P}}$.

6.3 Weibel instability

We now apply the 1X2P formulation defined above to the simulation of the Weibel instability. We initialize the simulation in close agreement with the parametrization of the Weibel instability simulations of Crouseilles et al. (Reference Crouseilles, Einkemmer and Faou2015) and Kraus et al. (Reference Kraus, Kormann, Morrison and Sonnendrücker2017). Specifically, in a periodic domain of 64 cells we simulate ${N_e=100\,032}$ electrons. Continuing to work in Gaussian units, we take perturbation wavenumber ${k=1.25({\omega _p}/{c})}$ and simulation domain length $L=2{\rm \pi} /k$. We seed a magnetic perturbation by defining the vector potential to be

(6.15)\begin{equation} {\boldsymbol{A}}=\frac{B_0}{k}\sin(kx)\hat{y}, \end{equation}

where ${B_0=-5.7\ \textrm {mG}}$. We sample initial electron velocities from the anisotropic distribution

(6.16)\begin{equation} f_e(x,v_x,v_y)=\frac{1}{2{\rm \pi}\sigma_x\sigma_y}\exp\left(-\left[\frac{v_x^{2}}{2\sigma_x^{2}}+\frac{v_y^{2}}{2\sigma_y^{2}}\right]\right), \end{equation}

where ${\sigma _x=0.02c/\sqrt {2}}$ and ${\sigma _y=\sqrt {12}\sigma _x}$. To initialize the initial electron distribution to be independent of $x$ (i.e. spatially uniform) as above, we evenly space electrons throughout the simulation domain at separation of $L/N_e$. $ {\boldsymbol {Y}}$ is correspondingly initialized to zero. $ {\boldsymbol {A}}$ and $ {\boldsymbol {Y}}$ are again modelled using Whitney (0- and 1-) forms.

Finally, simulating the Weibel instability with a Lie–Trotter splitting of stepsize $1/40\omega _p$ yields the evolution plotted in figure 6. The magnetic field growth rate is in strong agreement with the analytic prediction (Weibel Reference Weibel1959),

(6.17)\begin{equation} \gamma=\frac{\omega_pk\sigma_y}{\sqrt{\omega_p^{2}+c^{2}k^{2}}}. \end{equation}

Lastly, the total energy is also well bounded over the lifetime of the Weibel instability simulation, as depicted in figure 7.

Figure 6. The growth rate of the magnetic field energy closely approximates the analytic model of the Weibel instability. As the magnetic field ramps up, the anisotropy of the electron velocity is reduced. This plot may be compared with figure 1 of Kraus et al. (Reference Kraus, Kormann, Morrison and Sonnendrücker2017).

Figure 7. The log error in the total energy of a first-order Lie–Trotter splitting (Trotter Reference Trotter1959) Weibel instability simulation (Weibel Reference Weibel1959).