3.1 The continuous space–time Klimontovich–Maxwell model

We now use the preceding N2T procedure to systematically derive a charge conservation law for the continuous space–time Klimontovich–Maxwell system. This system specializes a Vlasov–Maxwell system to the following distribution function defined by $N$ point particles

(3.1)$$\begin{eqnarray}f(t,\boldsymbol{x},\boldsymbol{v})=\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D6FF}^{(3)}(\boldsymbol{x}-\boldsymbol{X}_{j}(t))\unicode[STIX]{x1D6FF}^{(3)}(\boldsymbol{v}-\dot{\boldsymbol{X}}_{j}(t)).\end{eqnarray}$$

The Klimontovich–Maxwell system is accordingly described by the following action:

(3.2)$$\begin{eqnarray}\displaystyle S=\int \text{d}^{4}x~{\mathcal{L}}[\unicode[STIX]{x1D719},\boldsymbol{A},\boldsymbol{X}_{i}] & = & \displaystyle \int \text{d}^{4}x\left[\frac{1}{2}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}+\unicode[STIX]{x2202}_{t}\boldsymbol{A})^{2}-\frac{1}{2}(\unicode[STIX]{x1D735}\times \boldsymbol{A})^{2}\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D6FF}_{j}\boldsymbol{\cdot }\left(\frac{1}{2}m_{j}\dot{\boldsymbol{X}}_{j}^{2}-q_{j}\unicode[STIX]{x1D719}+q_{j}\boldsymbol{A}\boldsymbol{\cdot }\dot{\boldsymbol{X}}_{j}\right)\right].\end{eqnarray}$$

Here, $\boldsymbol{A}=\boldsymbol{A}(t,\boldsymbol{x})$ is the vector potential, $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}(t,\boldsymbol{x})$ is the electric potential, $\boldsymbol{X}_{i}=\boldsymbol{X}_{i}(t)$ are particle positions and particle mass and charge are denoted by $m_{i}$ and $q_{i}$, respectively. We have also used the following shorthand for the delta function:

(3.3)$$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{j}:=\unicode[STIX]{x1D6FF}^{(3)}(\boldsymbol{x}-\boldsymbol{X}_{j}(t)).\end{eqnarray}$$

We apply Euler operators to derive the Euler–Lagrange equations of each field

(3.4)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D60C}_{\unicode[STIX]{x1D719}}({\mathcal{L}})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{E}-\unicode[STIX]{x1D70C},\\ \displaystyle \unicode[STIX]{x1D60C}_{\boldsymbol{A}}({\mathcal{L}})=\unicode[STIX]{x2202}_{t}\boldsymbol{E}-\unicode[STIX]{x1D735}\times \boldsymbol{B}+\boldsymbol{J},\\ \displaystyle \unicode[STIX]{x1D60C}_{\boldsymbol{X}_{i}}({\mathcal{L}})=\unicode[STIX]{x1D6FF}_{i}\boldsymbol{\cdot }[-m_{i}\ddot{\boldsymbol{X}}_{i}+q_{i}(\boldsymbol{E}+\dot{\boldsymbol{X}}_{i}\times \boldsymbol{B})],\end{array}\right\}\end{eqnarray}$$

where we have used the distributional derivative

(3.5)$$\begin{eqnarray}\int f(\unicode[STIX]{x1D702})\unicode[STIX]{x1D6FF}^{\prime }(\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D702}=-\int f^{\prime }(\unicode[STIX]{x1D702})\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D702},\end{eqnarray}$$

with $\unicode[STIX]{x1D702}\in \{t,\boldsymbol{x}\}$, and where

(3.6)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \boldsymbol{E}(t,\boldsymbol{x}):=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}(t,\boldsymbol{x})-\unicode[STIX]{x2202}_{t}\boldsymbol{A}(t,\boldsymbol{x}),\\ \boldsymbol{B}(t,\boldsymbol{x}):=\unicode[STIX]{x1D735}\times \boldsymbol{A}(t,\boldsymbol{x}),\\ \displaystyle \unicode[STIX]{x1D70C}(t,\boldsymbol{x}):=\mathop{\sum }_{j=1}^{N}q_{j}\unicode[STIX]{x1D6FF}_{j},\\ \displaystyle \boldsymbol{J}(t,\boldsymbol{x}):=\mathop{\sum }_{j=1}^{N}q_{j}\dot{\boldsymbol{X}}_{j}(t)\unicode[STIX]{x1D6FF}_{j}.\end{array}\right\}\end{eqnarray}$$

As noted in (2.2), an Euler operator $\unicode[STIX]{x1D60C}_{u}$ for an arbitrary field $u$ is essentially defined to accommodate integration by parts, such as that in (3.5). In particular, total derivatives in ${\mathcal{L}}$ – e.g. $(f\unicode[STIX]{x1D6FF})^{\prime }$ – that contribute to boundary terms of the action integral $S=\int {\mathcal{L}}\,\text{d}^{4}x$ – e.g. $f\unicode[STIX]{x1D6FF}|_{-\infty }^{\infty }$ – lie in the kernel of $\unicode[STIX]{x1D60C}_{u}$, such that $\unicode[STIX]{x1D60C}_{u}({\mathcal{L}}+\text{Div}\,\unicode[STIX]{x1D6FE})=\unicode[STIX]{x1D60C}_{u}({\mathcal{L}})$. Indeed, the operator relation $\unicode[STIX]{x1D60C}_{u}\circ \text{Div}=0$ always holds (see the ‘variational complex’ of Olver (Reference Olver1986)), where $\text{Div}$ denotes a divergence.

We now note that the action of (3.2) is invariant under the following gauge transformation:

(3.7)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D719}\rightarrow \unicode[STIX]{x1D719}^{\prime }=\unicode[STIX]{x1D719}+\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D706},\\ \boldsymbol{A}\rightarrow \boldsymbol{A}^{\prime }=\boldsymbol{A}-\unicode[STIX]{x1D735}\unicode[STIX]{x1D706}.\end{array}\right\}\end{eqnarray}$$

In particular, the electromagnetic terms of the Lagrangian are invariant, while the coupled particle terms pick up a divergence – namely, ${\mathcal{L}}\rightarrow {\mathcal{L}}+\unicode[STIX]{x2202}_{\unicode[STIX]{x1D707}}\unicode[STIX]{x1D6FE}^{\unicode[STIX]{x1D707}}$, where $\unicode[STIX]{x1D6FE}^{\unicode[STIX]{x1D707}}=-\sum _{j}q_{j}\unicode[STIX]{x1D6FF}_{j}\unicode[STIX]{x1D706}\cdot (1,\dot{\boldsymbol{X}}_{j})$ – that vanishes on the boundary of $S$. The vector field corresponding to this transformation – equivalent to (2.9) – is given by

(3.8)$$\begin{eqnarray}\boldsymbol{v}_{\unicode[STIX]{x1D706}}=\mathop{\sum }_{\unicode[STIX]{x1D6FC}}Q^{\unicode[STIX]{x1D6FC}}[\unicode[STIX]{x1D706}]\unicode[STIX]{x2202}_{u^{\unicode[STIX]{x1D6FC}}}=(\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D706})\unicode[STIX]{x2202}_{\unicode[STIX]{x1D719}}-(\unicode[STIX]{x1D735}\unicode[STIX]{x1D706})\boldsymbol{\cdot }\unicode[STIX]{x2202}_{\boldsymbol{A}}\end{eqnarray}$$

for an arbitrary smooth function $\unicode[STIX]{x1D706}(x)$.

Finally, given EOM in (3.4) and the characteristics of our gauge symmetry in (3.8), we may derive the differential identity of N2T using the construction of (2.6)

(3.9)$$\begin{eqnarray}\displaystyle \mathop{\sum }_{\unicode[STIX]{x1D6FC}}Q^{\unicode[STIX]{x1D6FC}}[\unicode[STIX]{x1D706}]\unicode[STIX]{x1D60C}_{u^{\unicode[STIX]{x1D6FC}}}({\mathcal{L}}) & = & \displaystyle Q^{\unicode[STIX]{x1D719}}[\unicode[STIX]{x1D706}]\boldsymbol{\cdot }\unicode[STIX]{x1D60C}_{\unicode[STIX]{x1D719}}({\mathcal{L}})+Q^{\boldsymbol{A}}[\unicode[STIX]{x1D706}]\boldsymbol{\cdot }\unicode[STIX]{x1D60C}_{\boldsymbol{A}}({\mathcal{L}})\nonumber\\ \displaystyle & = & \displaystyle (\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D706})[\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{E}-\unicode[STIX]{x1D70C}]-\unicode[STIX]{x1D735}\unicode[STIX]{x1D706}\boldsymbol{\cdot }[\unicode[STIX]{x2202}_{t}\boldsymbol{E}-\unicode[STIX]{x1D735}\times \boldsymbol{B}+\boldsymbol{J}]\end{eqnarray}$$

such that

(3.10)$$\begin{eqnarray}\displaystyle 0 & = & \displaystyle \unicode[STIX]{x1D60C}_{\unicode[STIX]{x1D706}}\left[\mathop{\sum }_{\unicode[STIX]{x1D6FC}}Q^{\unicode[STIX]{x1D6FC}}[\unicode[STIX]{x1D706}]\unicode[STIX]{x1D60C}_{u^{\unicode[STIX]{x1D6FC}}}({\mathcal{L}})\right]\nonumber\\ \displaystyle & = & \displaystyle -\unicode[STIX]{x2202}_{t}[\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{E}-\unicode[STIX]{x1D70C}]+\unicode[STIX]{x1D735}\boldsymbol{\cdot }[\unicode[STIX]{x2202}_{t}\boldsymbol{E}-\unicode[STIX]{x1D735}\times \boldsymbol{B}+\boldsymbol{J}]\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D70C}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{J}.\end{eqnarray}$$

In the final line, we have noted the equality of mixed partials and the vanishing divergence of the curl.

The N2T differential identity arising from the Klimontovich–Maxwell Lagrangian’s local gauge symmetry evidently discovers the charge conservation law itself. By construction, this conservation law must hold off-shell and identically; in particular, equation (3.10) does not require the equations of motion in order to hold true. It is a trivial conservation law – also referred to as a ‘strong’ or ‘improper’ conservation law (Brading & Brown Reference Brading and Brown2000) – an often-overlooked fact that is immediately verified upon examining the definitions of $\unicode[STIX]{x1D70C}$ and $\boldsymbol{J}$ in (3.6).

3.2 The geometric PIC method of Squire et al. (Reference Squire, Qin and Tang2012)

We now derive an analogous charge conservation law for the discrete, gauge-symmetric Vlasov–Maxwell PIC method defined by Squire et al. (Reference Squire, Qin and Tang2012). In this PIC scheme, space–time is discretized by a $d$-dimensional spatial simplicial complex (comprised of triangles in two dimensions or tetrahedra in three dimensions) whose structure is held constant throughout a uniformly discretized time. The time dimension may be envisaged as forming temporal edges that extend orthogonally from the spatial simplices, as in a triangular prism. We denote this ($d$+1)-dimensional prismal complex $P_{C}$. We use DEC (Desbrun et al. Reference Desbrun, Hirani, Leok and Marsden2005) to define fields on $P_{C}$ that are single valued on its $k$-cells (or their circumcentric duals) for $0\leqslant k\leqslant d+1$. In the present paper, we shall assume a spatial dimensionality $d=3$, such that $P_{C}$ is four-dimensional, with three-dimensional spatial tetrahedra comprising each time slice.

We first review some elements of DEC formalism that are necessary in the present study. It will be useful to distinguish the spatial edges from the temporal edges of $P_{C}$, so we denote a vertex of $P_{C}$ by , where $i$ is the spatial index of the vertex and $n$ is its temporal index. A discrete $0$-form $\unicode[STIX]{x1D6FC}$ is then defined by its values at each vertex, and a discrete $1$-form $\unicode[STIX]{x1D6FD}$ by its values on each edge

(3.11)

Here, we have expressed the discrete forms (equivalently, cochains) $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ in terms of their cochain bases, where $\unicode[STIX]{x1D6E5}_{n}^{i}$ is an element of the 0-cochain basis that maps  to 1 and all other vertices to 0; $\unicode[STIX]{x1D6E5}_{n}^{ij}$ is similarly an element of the 1-cochain basis that maps the oriented edge $[ij]$ to 1 and all others to 0. (A temporal edge is understood to be oriented in the positive time direction, and its cochain is denoted $\unicode[STIX]{x1D6E5}_{n-1/2}^{i}$.) Discrete $k$-forms of higher degree may be constructed with cochain bases in essentially the same way. The formalism of cochain bases will prove especially useful when we derive EOM for the dynamical fields on $P_{C}$ – that is, when we define a DEC Euler operator.

Let us denote the set of all vertices in $P_{C}$ by $\{v\}$, the set of spatial and temporal edges by $\{e\}=\{e_{s}\}\sqcup \{e_{t}\}$, and the set of spatial and ‘spatio-temporal’ faces by $\{\,f\}=\{\,f_{s}\}\sqcup \{\,f_{t}\}$. The DEC exterior derivative $\text{d}$, satisfying $\text{d}^{2}=0$, may be defined (Elcott & Schröder Reference Elcott and Schröder2005; Desbrun, Kanso & Tong Reference Desbrun, Kanso and Tong2006) by a matrix multiplication in the cochain basis. For a 1-form $\unicode[STIX]{x1D6FD}$, we see this as follows:

(3.12)$$\begin{eqnarray}\text{d}\unicode[STIX]{x1D6FD}=\text{d}(\unicode[STIX]{x1D6FD}_{e}\unicode[STIX]{x1D6E5}^{e})=\unicode[STIX]{x1D6FD}_{e}\,\text{d}\unicode[STIX]{x1D6E5}^{e}=\unicode[STIX]{x1D6FD}_{e}W_{f}^{e}\unicode[STIX]{x1D6E5}^{f},\end{eqnarray}$$

where the matrix entry $W_{f}^{e}$ stores the weight – $\{\pm 1,0\}$ – of the 1-cochain $\unicode[STIX]{x1D6E5}^{e}$ in the 2-cochain $\unicode[STIX]{x1D6E5}^{f}$. (We recall that the boundary operator on chains – $\unicode[STIX]{x2202}$ – is similarly determined by $W_{e}^{f}=(W_{f}^{e})^{T}$.) We adopt the Einstein summation convention in (3.12) and hereafter for prismal complex indices: $\{v\}$, $\{e\}$, and $\{\,f\}$.

For example, the electromagnetic gauge field $A$ – a discrete 1-form defined on all edges of $P_{C}$ – neatly splits in into an electric potential $\unicode[STIX]{x1D719}_{n-1/2}^{i}:=-A_{n-1/2}^{i}$ and a vector potential $\boldsymbol{A}_{n}^{ij}:=A_{n}^{ij}$, as follows:

(3.13)$$\begin{eqnarray}\displaystyle A & = & \displaystyle -\mathop{\sum }_{i,n}\unicode[STIX]{x1D719}_{n-1/2}^{i}\unicode[STIX]{x1D6E5}_{n-1/2}^{i}+\mathop{\sum }_{[ij],n}\boldsymbol{A}_{n}^{ij}\unicode[STIX]{x1D6E5}_{n}^{ij}\nonumber\\ \displaystyle & = & \displaystyle -\unicode[STIX]{x1D719}_{e_{t}}\unicode[STIX]{x1D6E5}^{e_{t}}+\boldsymbol{A}_{e_{s}}\unicode[STIX]{x1D6E5}^{e_{s}}.\end{eqnarray}$$

Using (3.12), we may correspondingly express $\text{d}A$ as

(3.14)$$\begin{eqnarray}\displaystyle \text{d}A & = & \displaystyle -\unicode[STIX]{x1D719}_{e_{t}}\,\text{d}\unicode[STIX]{x1D6E5}^{e_{t}}+\boldsymbol{A}_{e_{s}}\,\text{d}\unicode[STIX]{x1D6E5}^{e_{s}},\nonumber\\ \displaystyle & = & \displaystyle (-\unicode[STIX]{x1D719}_{e_{t}}W_{f_{t}}^{e_{t}}+\boldsymbol{A}_{e_{s}}W_{f_{t}}^{e_{s}})\unicode[STIX]{x1D6E5}^{f_{t}}+\boldsymbol{A}_{e_{s}}W_{f_{s}}^{e_{s}}\unicode[STIX]{x1D6E5}^{f_{s}}\nonumber\\ \displaystyle & = & \displaystyle E\wedge \text{d}t+B\nonumber\\ \displaystyle & = & \displaystyle F,\end{eqnarray}$$

where we have made use of the 1-form $\text{d}t:=\sum _{e_{t}}\unicode[STIX]{x1D6E5}^{e_{t}}$ and wedge product to implicitly define the spatial 1- and 2-forms $E$ and $B$, respectively, and the Faraday 2-form $F$ (Stern et al. Reference Stern, Tong, Desbrun, Marsden, Chang, Holm, Patrick and Ratiu2015). (As a note of caution, we emphasize that the preceding vector potential $\boldsymbol{A}$ is a single number on each spatial edge, and its bold notation is only suggestive. On the other hand, the Whitney interpolant of $\boldsymbol{A}$ will coordinatize $\mathbb{R}^{3}$ and thereby extend the single valued $\boldsymbol{A}$ from the spatial edges of $P_{C}$ to a 3-component vector field, as we shall see.)

The map from primal $k$-forms to dual $(4-k)$-forms on $P_{C}$ is effected via the metric-dependent Hodge star operator, $\star$. The Hodge star is defined (Abraham, Marsden & Ratiu Reference Abraham, Marsden and Ratiu1988; Desbrun et al. Reference Desbrun, Hirani, Leok and Marsden2005) such that the symmetric, metric-induced inner product $(\cdot ,\cdot )$ on $k$-forms satisfies

(3.15)$$\begin{eqnarray}(\unicode[STIX]{x1D714},\unicode[STIX]{x1D708})\unicode[STIX]{x1D707}=\unicode[STIX]{x1D714}\wedge \star \unicode[STIX]{x1D708}\end{eqnarray}$$

for primal $k$-forms $\unicode[STIX]{x1D714}$ and $\unicode[STIX]{x1D708}$ and volume top form $\unicode[STIX]{x1D707}$. For a $k$-form $\unicode[STIX]{x1D714}$ on an $n$-dimensional cell complex ($n=4$ on $P_{C}$), it can therefore be shown that

(3.16)$$\begin{eqnarray}\star (\star \unicode[STIX]{x1D714})=(-1)^{k(n-k)+\text{Ind}(g)}\unicode[STIX]{x1D714}.\end{eqnarray}$$

Here, $\text{Ind}(g)=\#\{\text{eig}[g]<0\}$ represents the index of the metric $g$. In the present context, we adopt a $(-+++)$ convention for our Lorentzian metric $g=\unicode[STIX]{x1D702}$, such that $\text{Ind}(g)=1$ on $P_{C}$. The dual $(4-k)$-form $\star \unicode[STIX]{x1D6FC}$ is thus defined on a dual chain $(\star \unicode[STIX]{x1D70E})^{4-k}$ as follows:

(3.17)$$\begin{eqnarray}\langle \star \unicode[STIX]{x1D6FC},\star \unicode[STIX]{x1D70E}\rangle =\unicode[STIX]{x1D716}(\unicode[STIX]{x1D70E})\frac{|\star \,\unicode[STIX]{x1D70E}|}{|\unicode[STIX]{x1D70E}|}\langle \unicode[STIX]{x1D6FC},\unicode[STIX]{x1D70E}\rangle ,\end{eqnarray}$$

where

(3.18)$$\begin{eqnarray}\unicode[STIX]{x1D716}(\unicode[STIX]{x1D70E})=\left\{\begin{array}{@{}ll@{}}+1\quad & \text{if }\unicode[STIX]{x1D70E}\text{ is entirely spacelike}\\ -1\quad & \text{otherwise}.\end{array}\right.\end{eqnarray}$$

Here, $|\unicode[STIX]{x1D70E}^{k}|$ denotes the $k$-volume of the $k$-dimensional $\unicode[STIX]{x1D70E}^{k}$ (where $|\unicode[STIX]{x1D70E}^{0}|=1$ for a single vertex), and $\langle \cdot ,\cdot \rangle$ denotes a $k$-cochain evaluated on a $k$-chain in (3.17).

Integration by parts on $P_{C}$ – which is necessary for the derivation of EOM, as in (2.2) – may be facilitated via the codifferential operator, $\unicode[STIX]{x1D6FF}$. Up to boundary contributions, $\unicode[STIX]{x1D6FF}$ is a formal adjoint to $\text{d}$, that is

(3.19)$$\begin{eqnarray}(\text{d}\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})\unicode[STIX]{x1D707}=(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FD})\unicode[STIX]{x1D707}+\text{d}(\unicode[STIX]{x1D6FC}\wedge \star \unicode[STIX]{x1D6FD}).\end{eqnarray}$$

When acting on a $k$-form defined on an $n$-dimensional complex, $\unicode[STIX]{x1D6FF}$ is given by (Abraham et al. Reference Abraham, Marsden and Ratiu1988; Desbrun et al. Reference Desbrun, Hirani, Leok and Marsden2005)

(3.20)$$\begin{eqnarray}\unicode[STIX]{x1D6FF}=(-1)^{n(k-1)+1+\text{Ind}(g)}\star \,\text{d}\star .\end{eqnarray}$$

We observe that, whereas $\text{d}$ maps a $k$-form to a $(k+1)$-form, $\unicode[STIX]{x1D6FF}$ reduces its degree to a $(k-1)$ form.

Having briefly reviewed relevant elements of DEC, we may now restate the discrete action of Squire et al. (Reference Squire, Qin and Tang2012), defined on $P_{C}$

(3.21)$$\begin{eqnarray}\displaystyle S & = & \displaystyle \mathop{\sum }_{V_{\unicode[STIX]{x1D70E}^{2}}}-\frac{1}{2}\,\text{d}A\wedge \star \,\text{d}A+\mathop{\sum }_{p,n}\left\{\frac{h}{2}m_{p}\left|\frac{\boldsymbol{X}_{n+1/2}^{p}-\boldsymbol{X}_{n-1/2}^{p}}{h}\right|^{2}-q_{p}\mathop{\sum }_{i}\unicode[STIX]{x1D719}_{n-1/2}^{i}\unicode[STIX]{x1D711}^{i}(\boldsymbol{X}_{n-1/2}^{p})\right.\nonumber\\ \displaystyle & & \displaystyle +\left.q_{p}\left(\frac{\boldsymbol{X}_{n+1/2}^{p}-\boldsymbol{X}_{n-1/2}^{p}}{h}\right)\boldsymbol{\cdot }\mathop{\sum }_{[ij]}\boldsymbol{A}_{n}^{ij}\int _{t_{n-1/2}}^{t_{n+1/2}}\,\text{d}t\,\unicode[STIX]{x1D711}^{ij}(\boldsymbol{X}^{p}(t))\right\}.\end{eqnarray}$$

In (3.21), we have denoted a sum over support volumes $V_{\unicode[STIX]{x1D70E}^{2}}$ for the primal-dual 4-form $\text{d}A\wedge \star \text{d}A$, with $A$ defined as in (3.13). $V_{\unicode[STIX]{x1D70E}^{2}}$ represents the convex hull of the 2-chain $\unicode[STIX]{x1D70E}^{2}$ and its dual $\star \unicode[STIX]{x1D70E}^{2}$ on which $\langle \text{d}A,\unicode[STIX]{x1D70E}^{2}\rangle$ and $\langle \star \text{d}A,\star \unicode[STIX]{x1D70E}^{2}\rangle$ are respectively defined. The symbol $h$ denotes the time step, $n$ the time index and $p$ the particle index. $\boldsymbol{X}^{p}(t)$ is defined as the constant velocity path between the particle’s staggered-time positions $\boldsymbol{X}_{n-1/2}^{p}$ and $\boldsymbol{X}_{n+1/2}^{p}$. In particular, particle paths are chosen to have straight line trajectories between the staggered times $t\in [(n-1/2)h,(n+1/2)h]$, $\forall n\in \mathbb{Z}$.

The Whitney 0-form $\unicode[STIX]{x1D711}^{i}(x)$ and 1-form $\unicode[STIX]{x1D711}^{ij}(x)$ interpolate $\unicode[STIX]{x1D719}$ and $\boldsymbol{A}$ to an arbitrary point $x\in P_{C}$ (Bossavit Reference Bossavit1988). In effect, $\unicode[STIX]{x1D711}^{i}$ and $\unicode[STIX]{x1D711}^{ij}$ complete the spatial components of the cochain bases $\unicode[STIX]{x1D6E5}_{n}^{i}$ and $\unicode[STIX]{x1D6E5}_{n}^{ij}$ adopted in (3.11) by extending DEC forms to the convex hull of $P_{C}$. In the continuous space–time of the Klimontovich–Maxwell system, the everywhere-defined gauge fields $\unicode[STIX]{x1D719}(t,\boldsymbol{x})$ and $\boldsymbol{A}(t,\boldsymbol{x})$ were ‘attached’ to point particles by the delta function, $\unicode[STIX]{x1D6FF}^{(3)}(\boldsymbol{x}-\boldsymbol{X}_{j}(t))$. In the prismal complex $P_{C}$, Whitney forms play this role by interpolating the gauge fields to the locations of point particles. Likewise, while we continue to avoid ascribing any geometric notion to point particles themselves, we see that Whitney forms on $P_{C}$ attach geometry to the charge densities and currents of the particles, as did the delta function in (3.6). For example, the spatial dot product in (3.21) composes $\boldsymbol{X}_{n+1/2}^{p}$ with the Whitney-interpolated 3-component vector field $\boldsymbol{A}_{n}^{ij}\unicode[STIX]{x1D711}^{ij}(\boldsymbol{X}^{p}(t))$, (where $\boldsymbol{A}_{n}^{ij}$ represents a single number and $\unicode[STIX]{x1D711}^{ij}$ a 3-component vector).

We now follow the continuous space–time N2T procedure of (3.4)–(3.10) by examining the equations of motion and gauge symmetry of $S$ in (3.21). As we have already seen, the differential structure of space–time has been replaced in the discrete setting by the prismal complex $P_{C}$, its DEC operators and Whitney forms. To derive the Euler–Lagrange equations of $S$, therefore, we must define an Euler operator for fields defined on this discrete structure. By analogy with (2.3), such an operator must implement a variational derivative on the space of fields defined on $P_{C}$.

Consider, for example, a $k$-form $\unicode[STIX]{x1D6FC}$ defined by its expansion in $k$-cochain basis elements: $\unicode[STIX]{x1D6FC}=\sum _{\unicode[STIX]{x1D70E}}\unicode[STIX]{x1D6FC}_{\unicode[STIX]{x1D70E}}\unicode[STIX]{x1D6E5}^{\unicode[STIX]{x1D70E}}$, where $\unicode[STIX]{x1D70E}$ ranges over all $k$-cells on $P_{C}$. Since each component $\unicode[STIX]{x1D6FC}_{\unicode[STIX]{x1D70E}}$ of $\unicode[STIX]{x1D6FC}$ can be varied independently, the variational derivative of $\unicode[STIX]{x1D6FC}_{\unicode[STIX]{x1D70E}}$ takes the simple form of a partial derivative. We therefore define the Euler operator $\unicode[STIX]{x1D60C}_{\unicode[STIX]{x1D6FC}_{\unicode[STIX]{x1D70E}}}$ on the action $S$ as follows:

(3.22)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D60C}_{\unicode[STIX]{x1D6FC}_{\unicode[STIX]{x1D70E}}}(S) & := & \displaystyle \frac{\unicode[STIX]{x2202}S}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}_{\unicode[STIX]{x1D70E}}}\nonumber\\ \displaystyle & = & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}_{\unicode[STIX]{x1D70E}}}\sum L.\end{eqnarray}$$

(We note that in a discrete setting, variational derivatives are made to act on the entire action $S$, rather than on the Lagrangian, because discrete Lagrangians are necessarily non-local.) As usual we will assume that all fields and their variations have compact support, such that any divergence term in $L$ – which contributes to $S=\sum L$ only at the boundary – vanishes under $\unicode[STIX]{x1D60C}_{\unicode[STIX]{x1D6FC}_{\unicode[STIX]{x1D70E}}}$. Equation (3.22) is the DEC counterpart to the continuous Euler operator defined in (2.3), and is now applied to derive our EOM.

To calculate $\unicode[STIX]{x1D60C}_{\unicode[STIX]{x1D719}_{e_{t}}}(S)$ and $\unicode[STIX]{x1D60C}_{\boldsymbol{A}_{e_{s}}}(S)$, we first re-express $\text{d}A\wedge \star \text{d}A$ in $S$ using (3.15)–(3.20) and the invariance of $S$ under the addition of a divergence ($L\rightarrow L+\text{d}\unicode[STIX]{x1D6FE}$)

(3.23)

where ${\approx}$ indicates equality up to an (ignorable) divergence. Then, using (3.13) and noting the symmetry of the intermediate expression $(\text{d}A,\text{d}A)\unicode[STIX]{x1D707}$ above, we apply the Euler operator of (3.22) to derive the EOM for $A$ as follows:

(3.24a)$$\begin{eqnarray}\displaystyle & \displaystyle 0=\unicode[STIX]{x1D60C}_{\unicode[STIX]{x1D719}_{e_{t}}}(S)=\unicode[STIX]{x1D6E5}^{e_{t}}\wedge \text{d}\star \,\text{d}A-\unicode[STIX]{x1D70C}^{e_{t}}, & \displaystyle\end{eqnarray}$$

(3.24b)$$\begin{eqnarray}\displaystyle & \displaystyle 0=\unicode[STIX]{x1D60C}_{\boldsymbol{A}_{e_{s}}}(S)=-\unicode[STIX]{x1D6E5}^{e_{s}}\wedge \text{d}\star \,\text{d}A+J^{e_{s}}, & \displaystyle\end{eqnarray}$$

(3.25a)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}^{e_{t}}:=\mathop{\sum }_{p}q_{p}\unicode[STIX]{x1D711}^{i}(\boldsymbol{X}_{t(e_{t})}^{p}), & \displaystyle\end{eqnarray}$$

(3.25b)$$\begin{eqnarray}\displaystyle & \displaystyle J^{e_{s}}:=\mathop{\sum }_{p}q_{p}\left(\frac{\boldsymbol{X}_{t_{f}(e_{s})}^{p}-\boldsymbol{X}_{t_{i}(e_{s})}^{p}}{h}\right)\boldsymbol{\cdot }\int _{t_{i}(e_{s})}^{t_{f}(e_{s})}\text{d}t\,\unicode[STIX]{x1D711}^{e_{s}}(\boldsymbol{X}^{p}(t)). & \displaystyle\end{eqnarray}$$

$i$

$e_{t}$

$\boldsymbol{X}_{t(e_{t})}^{p}$

$p$

$e_{t}$

$\boldsymbol{X}_{t_{i}(e_{s})}^{p}$

$\boldsymbol{X}_{t_{f}(e_{s})}^{p}$

$t=n$

$e_{s}$

It is worth pausing to interpret these EOM. We first observe that the primal-dual wedge product in (3.24a) is only non-vanishing on the spatial $(\star \unicode[STIX]{x1D6E5}^{e_{t}})$ component of $\text{d}\star \,\text{d}A$. This follows from the definition of the primal-dual wedge product, which is only non-zero on the convex hulls of a cell and its dual: $CH(\unicode[STIX]{x1D70E},\star \unicode[STIX]{x1D70E})$. Reading off from (3.14), therefore, equation (3.24a) becomes

(3.26)$$\begin{eqnarray}\unicode[STIX]{x1D70C}^{e_{t}}=\unicode[STIX]{x1D6E5}^{e_{t}}\wedge \text{d}\star (E\wedge \text{d}t)=\text{d}D\wedge \unicode[STIX]{x1D6E5}^{e_{t}},\end{eqnarray}$$

Gauss’s law for the electric displacement dual 2-form $D$, as expected. An analogous interpretation of (3.24b) yields a discrete Ampère–Maxwell law. We have omitted the $\unicode[STIX]{x1D60C}_{\boldsymbol{X}^{p}}(S)$ EOM for particle trajectories, as they will not be necessary for the derivation of charge conservation via N2T – just as they were unnecessary in (3.9)–(3.10). These implicit time-step particle EOM are derived in Squire et al. (Reference Squire, Qin and Tang2012).

Having derived our field equations of motion, we must now examine the gauge symmetry of the action $S$ in (3.21). In particular, $S$ is invariant under the local gauge transformation $A\rightarrow A-\text{d}f$, defined by

(3.27)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D719}_{n+1/2}^{i}\rightarrow \unicode[STIX]{x1D719}_{n+1/2}^{i}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}_{n+1/2}^{i}=\unicode[STIX]{x1D719}_{n+1/2}^{i}+(f_{n+1}^{i}-f_{n}^{i}),\\ \boldsymbol{A}_{n}^{ij}\rightarrow \boldsymbol{A}_{n}^{ij}+\unicode[STIX]{x1D6FF}\boldsymbol{A}_{n}^{ij}=\boldsymbol{A}_{n}^{ij}-(f_{n}^{j}-f_{n}^{i}),\end{array}\right\}\end{eqnarray}$$

where $f=f_{v}\unicode[STIX]{x1D6E5}^{v}$ is an arbitrary primal 0-form on $P_{C}$.

After all, the electromagnetic term of $S$ – $\text{d}A\wedge \star \text{d}A$ – is clearly invariant under $A\rightarrow A-\text{d}f$, since $\text{d}^{2}=0$. Furthermore, as noted in Squire et al. (Reference Squire, Qin and Tang2012), the gauge invariance of the particle terms of $S$ follows from a defining property of Whitney interpolation: $\text{d}_{\text{c}}((\unicode[STIX]{x1D6FC})_{\text{interp}})=(\text{d}_{\text{d}}\unicode[STIX]{x1D6FC})_{\text{interp}}$, where $\text{d}_{\text{c}}$ and $\text{d}_{\text{d}}$ denote continuous and discrete exterior derivatives, respectively and $(\cdot )_{\text{interp}}$ denotes Whitney interpolation. Equation (3.27) therefore transforms $L\rightarrow L+\text{d}_{c}\unicode[STIX]{x1D6FE}$, adding a divergence term analogous to the transformation of (3.7) for the continuous space–time Vlasov–Maxwell system.

Following (3.8), the gauge transformation of (3.27) is seen to be generated by a vector field

(3.28)$$\begin{eqnarray}\boldsymbol{v}_{f}=\mathop{\sum }_{\unicode[STIX]{x1D6FC}}Q^{\unicode[STIX]{x1D6FC}}[f]\unicode[STIX]{x2202}_{u^{\unicode[STIX]{x1D6FC}}}=\mathop{\sum }_{e_{t}}(\text{d}_{e_{t}}\,f)\unicode[STIX]{x2202}_{\unicode[STIX]{x1D719}_{e_{t}}}-\mathop{\sum }_{e_{s}}(\text{d}_{e_{s}}\,f)\unicode[STIX]{x2202}_{\boldsymbol{A}_{e_{s}}},\end{eqnarray}$$

where the coefficient $\text{d}_{e}f=f_{v_{2}}-f_{v_{1}}$ corresponds to the oriented edge $e=[v_{1}v_{2}]$, and where the sums are taken over all temporal and spatial edges, respectively.

We have thus gathered the necessary data to complete the N2T construction of (2.6) for our DEC system. As in (3.9), we note

(3.29)$$\begin{eqnarray}\displaystyle \mathop{\sum }_{\unicode[STIX]{x1D6FC}}Q^{\unicode[STIX]{x1D6FC}}[f]\unicode[STIX]{x1D60C}_{u^{\unicode[STIX]{x1D6FC}}}(S) & = & \displaystyle \mathop{\sum }_{e_{t}}Q^{\unicode[STIX]{x1D719}_{e_{t}}}[f]\boldsymbol{\cdot }\unicode[STIX]{x1D60C}_{\unicode[STIX]{x1D719}_{e_{t}}}(S)+\mathop{\sum }_{e_{s}}Q^{\boldsymbol{A}_{e_{s}}}[f]\boldsymbol{\cdot }\unicode[STIX]{x1D60C}_{\boldsymbol{A}_{e_{s}}}(S)\nonumber\\ \displaystyle & = & \displaystyle \mathop{\sum }_{e_{t}}(\text{d}_{e_{t}}\,f)\boldsymbol{\cdot }(\unicode[STIX]{x1D6E5}^{e_{t}}\wedge \text{d}\star \,\text{d}A-\unicode[STIX]{x1D70C}^{e_{t}})\nonumber\\ \displaystyle & & \displaystyle +\,\mathop{\sum }_{e_{s}}(-\text{d}_{e_{s}}\,f)\boldsymbol{\cdot }(-\unicode[STIX]{x1D6E5}^{e_{s}}\wedge \text{d}\star \,\text{d}A+J^{e_{s}}).\end{eqnarray}$$

We now observe that

(3.30)$$\begin{eqnarray}(\text{d}_{e}f)\unicode[STIX]{x1D6E5}^{e}=\text{d}(f_{v}\unicode[STIX]{x1D6E5}^{v})=\text{d}f,\end{eqnarray}$$

so applying the Euler operator for $f_{v}$ at vertex  to (3.29) yields

(3.31)$$\begin{eqnarray}\displaystyle 0 & = & \displaystyle \unicode[STIX]{x1D60C}_{f_{v}}\left[\mathop{\sum }_{\unicode[STIX]{x1D6FC}}Q^{\unicode[STIX]{x1D6FC}}[f]\unicode[STIX]{x1D60C}_{u^{\unicode[STIX]{x1D6FC}}}(S)\right]\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D60C}_{f_{v}}\left[\mathop{\sum }_{V_{\unicode[STIX]{x1D70E}^{1}}}\text{d}f\wedge \text{d}\star \,\text{d}A-\mathop{\sum }_{e_{t}}\left(\text{d}_{e_{t}}\,f\right)\boldsymbol{\cdot }\unicode[STIX]{x1D70C}^{e_{t}}-\mathop{\sum }_{e_{s}}(\text{d}_{e_{s}}\,f)\boldsymbol{\cdot }J^{e_{s}}\right]\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D6E5}^{v}\wedge \star \unicode[STIX]{x1D6FF}\star \,\text{d}\star \,\text{d}A+(\unicode[STIX]{x1D70C}_{n+1/2}^{i}-\unicode[STIX]{x1D70C}_{n-1/2}^{i})+\mathop{\sum }_{j}J^{[ij]}\nonumber\\ \displaystyle & = & \displaystyle (\unicode[STIX]{x1D70C}_{n+1/2}^{i}-\unicode[STIX]{x1D70C}_{n-1/2}^{i})+\mathop{\sum }_{j}J^{[ij]}\end{eqnarray}$$

since up to a sign, $(\unicode[STIX]{x1D6FF}\star \,\text{d})=(\star \,\text{d}\star \star \,\text{d})=\star \text{d}^{2}=0$. The sum of $J^{[ij]}$ over $j$ captures all spatial edges that terminate on vertex . The last equality of (3.31) reveals the desired charge conservation law on $P_{C}$. By its very construction through N2T, this conservation law is guaranteed to be an off-shell differential identity, as was (3.10). We readily verify this fact as follows.

First, we restrict our sources $\unicode[STIX]{x1D70C}$ and $J$ to a particle of charge $q$ whose path over one time step remains within a single spatial tetrahedron; the general case follows without significant alteration. We then recall (Bossavit Reference Bossavit1988) that the Whitney 0-form $\unicode[STIX]{x1D711}^{i}$ interpolates from vertex $i$ via barycentric coordinates such that, over the tetrahedron $[ijk\ell ]$,

(3.32)$$\begin{eqnarray}\unicode[STIX]{x1D711}^{i}+\unicode[STIX]{x1D711}^{\,j}+\unicode[STIX]{x1D711}^{k}+\unicode[STIX]{x1D711}^{\ell }=1.\end{eqnarray}$$

In vector form, the 1-form $\unicode[STIX]{x1D711}^{ij}$ is then given by

(3.33)$$\begin{eqnarray}\unicode[STIX]{x1D711}^{ij}=\unicode[STIX]{x1D711}^{i}\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}^{\,j}-\unicode[STIX]{x1D711}^{\,j}\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}^{i}.\end{eqnarray}$$

Summing over the three spatial edges terminating on vertex $i$ of the tetrahedron containing the particle, therefore

(3.34)$$\begin{eqnarray}\displaystyle \mathop{\sum }_{j\neq i}\unicode[STIX]{x1D711}^{ij} & = & \displaystyle \unicode[STIX]{x1D711}^{i}\unicode[STIX]{x1D735}\left(\mathop{\sum }_{j\neq i}\unicode[STIX]{x1D711}^{\,j}\right)-\left(\mathop{\sum }_{j\neq i}\unicode[STIX]{x1D711}^{\,j}\right)\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}^{i}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D711}^{i}\unicode[STIX]{x1D735}(1-\unicode[STIX]{x1D711}^{i})-(1-\unicode[STIX]{x1D711}^{i})\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}^{i}\nonumber\\ \displaystyle & = & \displaystyle -\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}^{i}.\end{eqnarray}$$

It follows, then, that

(3.35)$$\begin{eqnarray}\displaystyle \mathop{\sum }_{j\neq i}J^{[ij]} & = & \displaystyle q\left(\frac{\boldsymbol{X}_{f}-\boldsymbol{X}_{i}}{h}\right)\boldsymbol{\cdot }\int _{t_{i}}^{t_{f}}\text{d}t\mathop{\sum }_{j\neq i}\unicode[STIX]{x1D711}^{ij}(\boldsymbol{X}(t))\nonumber\\ \displaystyle & = & \displaystyle -q\left(\frac{\boldsymbol{X}_{f}-\boldsymbol{X}_{i}}{h}\right)\boldsymbol{\cdot }\int _{t_{i}}^{t_{f}}\text{d}t\,\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}^{i}(\boldsymbol{X}(t))\nonumber\\ \displaystyle & = & \displaystyle -q\int _{i}^{f}\boldsymbol{v}\,\lrcorner \,\text{d}\unicode[STIX]{x1D711}^{i}\nonumber\\ \displaystyle & = & \displaystyle -q[\unicode[STIX]{x1D711}^{i}(\boldsymbol{X}_{f})-\unicode[STIX]{x1D711}^{i}(\boldsymbol{X}_{i})],\end{eqnarray}$$

where $\boldsymbol{v}\,\lrcorner \,\text{d}\unicode[STIX]{x1D711}^{i}$ is the interior product of the exact form $\text{d}\unicode[STIX]{x1D711}^{i}$ with respect to the velocity $\boldsymbol{v}:=(1/h)(\boldsymbol{X}_{f}-\boldsymbol{X}_{i})$, which is constant over a single time step of the particle. Upon comparison with the definition for $\unicode[STIX]{x1D70C}$ in (3.25a), it is clear that (3.31) holds off-shell, as desired. The N2T formalism of Hydon & Mansfield (Reference Hydon and Mansfield2011) has succeeded in deriving the off-shell, discrete conservation law.

Before we depart from the Lagrangian formalism, we note that an alternative approach to deriving the conservation laws of the continuous space–time and DEC Vlasov–Maxwell actions – equations (3.2) and (3.21) – entails gauge fixing these actions by setting $\unicode[STIX]{x1D719}(x)=0$ and $\unicode[STIX]{x1D719}_{n+1/2}^{i}=0$, respectively. Such a gauge fixing removes these systems’ degeneracy and uniquely determines solutions to their equations of motion. In such an approach, N1T is applied to the time-independent symmetry transformation $\boldsymbol{A}(t,\boldsymbol{x})\rightarrow \boldsymbol{A}(t,\boldsymbol{x})-\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}(\boldsymbol{x})$, thereby deriving a non-trivial conservation law in the form of a time evolution of Gauss’s law. In the ensuing sections, we pursue an analogous gauge-fixing approach for the Hamiltonian Vlasov–Maxwell system, employing the Hamiltonian formalism’s counterpart to N1T – the momentum map.

4.1 The Poisson structure of the Vlasov–Maxwell system

We first recall the Poisson bracket of Marsden & Weinstein (Reference Marsden and Weinstein1982) for the Vlasov–Maxwell system,

(4.1)$$\begin{eqnarray}\{\{F,G\}\}[f,\boldsymbol{A},\boldsymbol{Y}]=\int \text{d}\boldsymbol{x}\,\text{d}\boldsymbol{p}\,f\left\{\frac{\unicode[STIX]{x1D6FF}F}{\unicode[STIX]{x1D6FF}f},\frac{\unicode[STIX]{x1D6FF}G}{\unicode[STIX]{x1D6FF}f}\right\}_{\boldsymbol{x}\boldsymbol{p}}+\int \text{d}\boldsymbol{x}\left(\frac{\unicode[STIX]{x1D6FF}F}{\unicode[STIX]{x1D6FF}\boldsymbol{A}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x1D6FF}G}{\unicode[STIX]{x1D6FF}\boldsymbol{Y}}-\frac{\unicode[STIX]{x1D6FF}G}{\unicode[STIX]{x1D6FF}\boldsymbol{A}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x1D6FF}F}{\unicode[STIX]{x1D6FF}\boldsymbol{Y}}\right),\end{eqnarray}$$

with time evolution defined by the Hamiltonian

(4.2)$$\begin{eqnarray}H[f,\boldsymbol{A},\boldsymbol{Y}]=\frac{1}{2}\int f\boldsymbol{\cdot }|\boldsymbol{p}-\boldsymbol{A}|^{2}\,\text{d}\boldsymbol{x}\,\text{d}\boldsymbol{p}+\frac{1}{2}\int [|\boldsymbol{Y}|^{2}+|\unicode[STIX]{x1D735}\times \boldsymbol{A}|^{2}]\,\text{d}\boldsymbol{x}.\end{eqnarray}$$

Here, $F$ and $G$ represent arbitrary functionals of the distribution function $f(\boldsymbol{x},\boldsymbol{p})$, the 3-component vector potential $\boldsymbol{A}(\boldsymbol{x})$ and its conjugate momentum $\boldsymbol{Y}(\boldsymbol{x})$. As we shall see momentarily, $\boldsymbol{Y}$ can be readily identified as negative the electric field strength – (i.e. $\boldsymbol{Y}=-\boldsymbol{E}$). We note that our system is rendered in the temporal gauge, wherein the electric potential satisfies $\unicode[STIX]{x1D719}(\boldsymbol{x})=0$. As in Marsden & Weinstein (Reference Marsden and Weinstein1982), we denote the Poisson bracket in (4.1) by $\{\{\cdot ,\cdot \}\}$ merely to distinguish it from other Poisson structures.

The $\int \text{d}\boldsymbol{x}\,\text{d}\boldsymbol{p}\,f\{\unicode[STIX]{x1D6FF}_{f}\cdot ,{\unicode[STIX]{x1D6FF}_{f}\cdot \}}_{\boldsymbol{x}\boldsymbol{p}}$ operator in the first line of (4.1) is a Lie–Poisson bracket (Marsden & Ratiu Reference Marsden and Ratiu1999), which defines a Poisson structure for functions on a dual Lie algebra $\mathfrak{g}^{\ast }$. In general, the Lie–Poisson bracket on an arbitrary dual Lie algebra $\mathfrak{g}^{\ast }$ is defined to inherit the bracket $[\cdot ,\cdot ]$ of its underlying Lie algebra $\mathfrak{g}$ as follows:

(4.3)$$\begin{eqnarray}\{F,G\}(\unicode[STIX]{x1D6FC}):=-\left\langle \unicode[STIX]{x1D6FC},\left[\frac{\unicode[STIX]{x1D6FF}F}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FC}},\frac{\unicode[STIX]{x1D6FF}G}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FC}}\right]\right\rangle .\end{eqnarray}$$

The bracket of (4.3) is defined $\forall$$F,G\in C^{\infty }(\mathfrak{g}^{\ast })$ with respect to some fixed $\unicode[STIX]{x1D6FC}\in \mathfrak{g}^{\ast }$, where $\langle \cdot ,\cdot \rangle$ represents the linear pairing of elements of $\mathfrak{g}^{\ast }$ and $\mathfrak{g}$. The functional derivative $\unicode[STIX]{x1D6FF}F/\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FC}\in \mathfrak{g}^{\ast \ast }$ can be seen to define a linear function on $\mathfrak{g}^{\ast }$, in that it acts as a directional derivative on the functional $F$ at the ‘point’ $\unicode[STIX]{x1D6FC}\in \mathfrak{g}^{\ast }$. In particular, for arbitrary $\unicode[STIX]{x1D6FD}\in \mathfrak{g}^{\ast }$

(4.4)$$\begin{eqnarray}\left\langle \unicode[STIX]{x1D6FD},\frac{\unicode[STIX]{x1D6FF}F}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FC}}\right\rangle =D_{\unicode[STIX]{x1D6FC}}F\cdot \unicode[STIX]{x1D6FD}=\left.\frac{\text{d}}{\text{d}\unicode[STIX]{x1D716}}\right|_{\unicode[STIX]{x1D716}=0}F(\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D6FD}).\end{eqnarray}$$

Since $\mathfrak{g}^{\ast \ast }\cong \mathfrak{g}$, the functional derivative may be interpreted as an element of the Lie algebra.

In the present context, the Lie algebra $\mathfrak{g}$ corresponds to infinitesimal transformations of $(\boldsymbol{x},\boldsymbol{p})\cong \mathbb{R}^{6}$, the position–momentum phase space. Such transformations can be regarded as Hamiltonian vector fields on $\mathbb{R}^{6}$, which map via anti-homomorphism to their corresponding generating functions, i.e.

(4.5)$$\begin{eqnarray}[\boldsymbol{X}_{h},\boldsymbol{X}_{k}]=-\{h,k\}_{\boldsymbol{x}\boldsymbol{p}}.\end{eqnarray}$$

The bracket $\{\cdot ,\cdot \}_{\boldsymbol{x}\boldsymbol{p}}$ therefore serves as a Lie bracket, defined pointwise on $\mathbb{R}^{6}$

(4.6)$$\begin{eqnarray}\{h,k\}_{\boldsymbol{x}\boldsymbol{p}}:=(\unicode[STIX]{x2202}_{\boldsymbol{x}}h\boldsymbol{\cdot }\unicode[STIX]{x2202}_{\boldsymbol{p}}k-\unicode[STIX]{x2202}_{\boldsymbol{x}}k\boldsymbol{\cdot }\unicode[STIX]{x2202}_{\boldsymbol{p}}h).\end{eqnarray}$$

The dual Lie algebra $\mathfrak{g}^{\ast }$ is similarly identified by distribution densities on $\mathbb{R}^{6}$, which pair linearly to Hamiltonian functions via integration

(4.7)$$\begin{eqnarray}\langle f,h\rangle :=\int fh\,\text{d}\boldsymbol{x}\,\text{d}\boldsymbol{p}\end{eqnarray}$$

for $f\in \mathfrak{g}^{\ast },h\in \mathfrak{g}$.

In this way, the operator $\int \text{d}\boldsymbol{x}\,\text{d}\boldsymbol{p}~f\{\unicode[STIX]{x1D6FF}_{f}\cdot ,{\unicode[STIX]{x1D6FF}_{f}\cdot \}}_{\boldsymbol{x}\boldsymbol{p}}$ comprising the first term of (4.1) is seen to be a Lie–Poisson bracket of the form in (4.3). We note that the negative sign of (4.3) cancels with the negative sign of the anti-homomorphism of (4.5) to produce this operator.

The second term of (4.1) represents the electromagnetic ‘sector’ of our Poisson structure, and derives from the canonical symplectic structure on the cotangent space – $T^{\ast }Q=\{(\boldsymbol{A},\boldsymbol{Y})\}$ – of the configuration space $Q=\{\boldsymbol{A}\}$. Therefore, the complete setting of the Vlasov–Maxwell system is a Poisson manifold, given by

(4.8)$$\begin{eqnarray}M=\mathfrak{g}^{\ast }\times T^{\ast }Q,\end{eqnarray}$$

with its bracket defined in (4.1).

We now consider dynamics on this Poisson manifold $M$. To derive our Hamiltonian EOM, it is convenient to define functionals

(4.9)$$\begin{eqnarray}F(\boldsymbol{u}):=\int \text{d}\boldsymbol{u}^{\prime }\,F(\boldsymbol{u}^{\prime })\unicode[STIX]{x1D6FF}(\boldsymbol{u}-\boldsymbol{u}^{\prime })\end{eqnarray}$$

for $F\in \{\,f,\boldsymbol{A},\boldsymbol{Y}\}$ as in Kraus et al. (Reference Kraus, Kormann, Morrison and Sonnendrücker2017), where $\boldsymbol{u}=(\boldsymbol{x},\boldsymbol{p})$ or $\boldsymbol{u}=\boldsymbol{x}$, as appropriate. Plugging such functionals into (4.1)–(4.2), we find

(4.10)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\dot{f}}(\boldsymbol{x},\boldsymbol{p})=\{\{\,f,H\}\}=-[\unicode[STIX]{x2202}_{\boldsymbol{x}}\,f+\unicode[STIX]{x2202}_{\boldsymbol{p}}f\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\boldsymbol{A})]\boldsymbol{\cdot }(\boldsymbol{p}-\boldsymbol{A}),\\ \dot{\boldsymbol{A}}(\boldsymbol{x})=\{\{\boldsymbol{A},H\}\}=\boldsymbol{Y},\\ \displaystyle \dot{\boldsymbol{Y}}(\boldsymbol{x})=\{\{\boldsymbol{Y},H\}\}=\int f\boldsymbol{\cdot }(\boldsymbol{p}-\boldsymbol{A})\,\text{d}\boldsymbol{p}-\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D735}\times \boldsymbol{A}.\end{array}\right\}\end{eqnarray}$$

We observe that $\boldsymbol{Y}$ plays the role of $-\boldsymbol{E}$, as expected. For convenience, we note that the familiar form of the Vlasov equation may be recovered from the first line of (4.10) by defining a distribution density $\bar{f}$ on $(\boldsymbol{x},\boldsymbol{v})$ space where $\boldsymbol{v}=\boldsymbol{p}-\boldsymbol{A}$, i.e.

(4.11)$$\begin{eqnarray}f(\boldsymbol{x},\boldsymbol{p})=\bar{f}(\boldsymbol{x},\boldsymbol{p}-\boldsymbol{A})=\bar{f}(\boldsymbol{x},\boldsymbol{v}),\end{eqnarray}$$

such that $\unicode[STIX]{x2202}_{\boldsymbol{x}}\,f=\unicode[STIX]{x2202}_{\boldsymbol{x}}\,\bar{f}-(\unicode[STIX]{x1D735}\boldsymbol{A})\boldsymbol{\cdot }\unicode[STIX]{x2202}_{\boldsymbol{v}}\,\bar{f}$; $\unicode[STIX]{x2202}_{\boldsymbol{p}}f=\unicode[STIX]{x2202}_{\boldsymbol{v}}\,\bar{f}$; and ${\dot{f}}=\unicode[STIX]{x2202}_{t}\,\bar{f}-\dot{\boldsymbol{A}}\boldsymbol{\cdot }\unicode[STIX]{x2202}_{\boldsymbol{v}}\,\bar{f}$. Here, we use $\unicode[STIX]{x1D735}\equiv \unicode[STIX]{x2202}_{\boldsymbol{x}}$ interchangeably, and adopt the dyad convention

(4.12)$$\begin{eqnarray}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{A}\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{w}=v_{i}A_{i}B_{j}w_{j}\end{eqnarray}$$

in Einstein notation.

4.2 Gauge symmetry and the momentum map

With our Poisson and Hamiltonian structure in hand, we now examine the gauge symmetry of the Vlasov–Maxwell system. Continuing to follow Marsden & Weinstein (Reference Marsden and Weinstein1982), we define a group action $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D713}}:M\rightarrow M$ on our Poisson manifold $M=\mathfrak{g}^{\ast }\times T^{\ast }Q$ of the form

(4.13)$$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D713}}:(f,\boldsymbol{A},\boldsymbol{Y})\mapsto (f\circ \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D713}},\boldsymbol{A}-\unicode[STIX]{x1D735}\unicode[STIX]{x1D713},\boldsymbol{Y}),\end{eqnarray}$$

where

(4.14)$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D713}}(\boldsymbol{x},\boldsymbol{p}):=(\boldsymbol{x},\boldsymbol{p}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}).\end{eqnarray}$$

We emphasize that $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D713}}$ transforms $f$, and not $\boldsymbol{p}$ itself. It is straightforward to check that $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D713}}$ is a canonical group action, i.e. that the Poisson bracket is preserved by the pullback of $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D713}}$, namely $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D713}}^{\ast }\{\{F,G\}\}=\{\{\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D713}}^{\ast }F,\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D713}}^{\ast }G\}\}$.

We define such an arbitrary function $\unicode[STIX]{x1D713}\in {\mathcal{F}}$ as belonging to the abelian group ${\mathcal{F}}:=C^{\infty }(\mathbb{R}^{3})$ of smooth functions on $\mathbb{R}^{3}$, with the group composition law of addition. Its Lie algebra $\mathfrak{f}$ is also identifiable as the smooth functions on $\mathbb{R}^{3}$, while its dual $\mathfrak{f}^{\ast }$ is the set of densities over $\mathbb{R}^{3}$ that pair to elements of $\mathfrak{f}$ via integration over $\mathbb{R}^{3}$ – analogous to the $\mathbb{R}^{6}$ integration of (4.7).

Now let $\unicode[STIX]{x1D719}\in \mathfrak{f}$ denote an arbitrary Lie algebra element, such that $\exp (\unicode[STIX]{x1D716}\unicode[STIX]{x1D719})\in {\mathcal{F}}$$\forall$$\unicode[STIX]{x1D716}\in \mathbb{R}$. By differentiating the group action $\unicode[STIX]{x1D6F7}_{\exp (\unicode[STIX]{x1D716}\unicode[STIX]{x1D719})}$ on $M$, we may associate to any such $\unicode[STIX]{x1D719}\in \mathfrak{f}$ the corresponding vector field $\unicode[STIX]{x1D719}_{M}$ on $M$, namely

(4.15)$$\begin{eqnarray}\unicode[STIX]{x1D719}_{M}:=\left.\frac{\text{d}}{\text{d}\unicode[STIX]{x1D716}}\right|_{\unicode[STIX]{x1D716}=0}\unicode[STIX]{x1D6F7}_{\exp (\unicode[STIX]{x1D716}\unicode[STIX]{x1D719})}.\end{eqnarray}$$

The vector field $\unicode[STIX]{x1D719}_{M}$ is therefore the infinitesimal generator of the group action on $M$ corresponding to $\unicode[STIX]{x1D719}\in \mathfrak{f}$.

For any canonical group action on Poisson manifold $M$, we may seek a corresponding momentum map. The momentum map $\unicode[STIX]{x1D707}:M\rightarrow \mathfrak{f}^{\ast }$ of a group action is defined such that, $\forall$$\unicode[STIX]{x1D719}\in \mathfrak{f}$ and $m\in M$, the induced function

(4.16)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D707}^{\unicode[STIX]{x1D719}}:M\rightarrow \mathbb{R}\\ m\mapsto \langle \unicode[STIX]{x1D707}(m),\unicode[STIX]{x1D719}\rangle \end{array}\right\}\end{eqnarray}$$

satisfies

(4.17)$$\begin{eqnarray}\{\{F,\unicode[STIX]{x1D707}^{\unicode[STIX]{x1D719}}\}\}=\unicode[STIX]{x1D719}_{M}(F)\end{eqnarray}$$

for arbitrary $F\in C^{\infty }(M)$. Here, $\unicode[STIX]{x1D719}_{M}(F)$ is the Lie derivative of $F$ along the vector field $\unicode[STIX]{x1D719}_{M}$. In particular, the momentum map $\unicode[STIX]{x1D707}$ assigns a dual element of $\mathfrak{f}^{\ast }$ to each point of $M$ such that, when $\unicode[STIX]{x1D707}$ is everywhere paired with an element $\unicode[STIX]{x1D719}\in \mathfrak{f}$ of the Lie algebra, the resulting function $\unicode[STIX]{x1D707}^{\unicode[STIX]{x1D719}}$ on $M$ is a generating function of the associated vector field $\unicode[STIX]{x1D719}_{M}$.

The preceding definition of $\unicode[STIX]{x1D707}$ is general to arbitrary Poisson systems with canonical group actions, and we now apply it to find $\unicode[STIX]{x1D707}$ for the Vlasov–Maxwell system of interest. We first note that a single point $m\in M=\mathfrak{g}^{\ast }\times T^{\ast }Q$ specifies $(f,\boldsymbol{A},\boldsymbol{Y})$ over the entire $(\boldsymbol{x},\boldsymbol{p})$ phase space. Given the group action defined in (4.13)–(4.14), it is immediately seen that $\unicode[STIX]{x1D719}_{M}$ can be expressed as the following infinitesimal operator on $M$ corresponding to $\unicode[STIX]{x1D719}(\boldsymbol{x})\in \mathfrak{f}$:

(4.18)$$\begin{eqnarray}\{\{\cdot ,\unicode[STIX]{x1D707}^{\unicode[STIX]{x1D719}}\}\}=\int \text{d}\boldsymbol{x}\,\text{d}\boldsymbol{p}\,\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\boldsymbol{p}}\frac{\unicode[STIX]{x1D6FF}}{\unicode[STIX]{x1D6FF}f}-\int \text{d}\boldsymbol{x}\,\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}\boldsymbol{\cdot }\frac{\unicode[STIX]{x1D6FF}}{\unicode[STIX]{x1D6FF}\boldsymbol{A}}.\end{eqnarray}$$

Upon inspection, it is evident that to generate the operator of (4.18), the Poisson bracket of (4.1) requires that $\unicode[STIX]{x1D707}^{\unicode[STIX]{x1D719}}$ be given by

(4.19)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}^{\unicode[STIX]{x1D719}}(m) & := & \displaystyle \langle \unicode[STIX]{x1D707}(m),\unicode[STIX]{x1D719}\rangle \nonumber\\ \displaystyle & = & \displaystyle \int \text{d}\boldsymbol{x}\left[\int \text{d}\boldsymbol{p}\,f(\boldsymbol{x},\boldsymbol{p})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{Y}\right]\unicode[STIX]{x1D719}(\boldsymbol{x}),\end{eqnarray}$$

where $\langle \cdot ,\cdot \rangle =\int \text{d}\boldsymbol{x}$. Therefore, the momentum map must be

(4.20)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}(m) & = & \displaystyle \int \text{d}\boldsymbol{p}\,f(\boldsymbol{x},\boldsymbol{p})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{Y}\nonumber\\ \displaystyle & := & \displaystyle \unicode[STIX]{x1D70C}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{Y},\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}(\boldsymbol{x}):=\int \text{d}\boldsymbol{p}\,f(\boldsymbol{x},\boldsymbol{p})$. We note that $\unicode[STIX]{x1D707}^{\unicode[STIX]{x1D719}}:M\rightarrow \mathbb{R}$ is a real-valued function while $\unicode[STIX]{x1D707}(m)\in \mathfrak{f}^{\ast }$ is a density on $\mathbb{R}^{3}$, as desired.

For later use, we further observe that $\unicode[STIX]{x1D707}$ is group equivariant

(4.21)$$\begin{eqnarray}\unicode[STIX]{x1D707}\circ \unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D713}}=\text{Ad}_{\unicode[STIX]{x1D713}^{-1}}^{\ast }\circ \unicode[STIX]{x1D707},\end{eqnarray}$$

where $\text{Ad}_{\unicode[STIX]{x1D713}^{-1}}^{\ast }$ represents the coadjoint action (Marsden & Ratiu Reference Marsden and Ratiu1999) of $\unicode[STIX]{x1D713}\in {\mathcal{F}}$ on an element of $\mathfrak{f}^{\ast }$. In particular, it is clear by inspection of (4.20) that $\unicode[STIX]{x1D707}$ is invariant under ${\mathcal{F}}$ transformations, and since ${\mathcal{F}}$ is abelian, its coadjoint action on $\mathfrak{f}^{\ast }$ is simply the identity map.

4.3 Deriving the conservation law

The momentum map $\unicode[STIX]{x1D707}$ is defined as above – equations (4.15)–(4.17) – for any Poisson manifold $M$ with a canonical group action $\unicode[STIX]{x1D6F7}$. If it should happen that a Hamiltonian $H$ is furthermore defined on $M$ such that $H$ is invariant under $\unicode[STIX]{x1D6F7}$, then the momentum map $\unicode[STIX]{x1D707}$ so-constructed further guarantees a conservation law for the system.

Let us show this for our Vlasov–Maxwell system. We first note that $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D713}}$ leaves the Hamiltonian invariant, $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D713}}^{\ast }H=H$. By differentiating this expression with respect to $\unicode[STIX]{x1D713}$ as in (4.15), it is seen that, infinitesimally

(4.22)$$\begin{eqnarray}0=\unicode[STIX]{x1D719}_{M}(H)=\{\{H,\unicode[STIX]{x1D707}^{\unicode[STIX]{x1D719}}\}\}=-\{\{\unicode[STIX]{x1D707}^{\unicode[STIX]{x1D719}},H\}\}=-\dot{\unicode[STIX]{x1D707}}^{\unicode[STIX]{x1D719}}.\end{eqnarray}$$

Each linearly independent $\unicode[STIX]{x1D719}\in \mathfrak{f}$ therefore determines a unique first integral of the system – i.e. $\unicode[STIX]{x1D707}^{\unicode[STIX]{x1D719}}$.

We can make a stronger observation as well. Since $\dot{\unicode[STIX]{x1D707}}^{\unicode[STIX]{x1D719}}=0$ holds for arbitrary $\unicode[STIX]{x1D719}\in \mathfrak{f}$, the entire momentum map is invariant under the flow of $H$, that is

(4.23)$$\begin{eqnarray}\dot{\unicode[STIX]{x1D707}}=\{\{\unicode[STIX]{x1D707},H\}\}=0.\end{eqnarray}$$

This follows rigorously from the fundamental lemma of variational calculus applied to $\dot{\unicode[STIX]{x1D707}}^{\unicode[STIX]{x1D719}}$ via (4.19). As a result, we can apply the definition of (4.20) to derive

(4.24)$$\begin{eqnarray}0=\dot{\unicode[STIX]{x1D707}}=\dot{\unicode[STIX]{x1D70C}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\dot{\boldsymbol{Y}}.\end{eqnarray}$$

This completes the canonical derivation of the Vlasov–Maxwell local conservation law – $\dot{\unicode[STIX]{x1D707}}=0$ – in the continuous Hamiltonian formalism. We note that, setting $\boldsymbol{Y}=-\boldsymbol{E}$, equation (4.24) is the time evolution of Gauss’s law.

With an additional substitution to (4.24) from the EOM for $\dot{\boldsymbol{Y}}$ in (4.10), we may re-express this canonical conservation law in the form

(4.25)$$\begin{eqnarray}0=\dot{\unicode[STIX]{x1D70C}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{J},\end{eqnarray}$$

where $\boldsymbol{J}:=\int f\boldsymbol{\cdot }(\boldsymbol{p}-\boldsymbol{A})\,\text{d}\boldsymbol{p}$. Here $\unicode[STIX]{x1D70C}$ and $\boldsymbol{J}$ are (scalar and vector) densities over $\mathbb{R}^{3}$ and functionals in the sense of (4.9). This charge conservation law may be immediately checked by substituting the expression for ${\dot{f}}$ from (4.10). In the present Hamiltonian context, it is evident that (4.25) can no longer be regarded as an off-shell identity. (After all, time evolution itself is only ‘dynamically defined’, so to speak, by the Hamiltonian.)

4.4 Reduction of the Vlasov–Maxwell system

Finally, we undertake the Poisson reduction (Marsden & Weinstein Reference Marsden and Weinstein1974; Marsden & Ratiu Reference Marsden and Ratiu1986) of the Vlasov–Maxwell system. Given a Poisson manifold $(M,\{\cdot ,\cdot \}_{M})$ on which a Lie group $G$ acts by Poisson diffeomorphisms, the Poisson reduction of $M$ is the unique quotient map $\unicode[STIX]{x03C0}:(M,\{\cdot ,\cdot \}_{M})\rightarrow (M/G,\{\cdot ,\cdot \}_{M/G})$ satisfying $\unicode[STIX]{x03C0}^{\ast }\{\,f,g\}_{M/G}=\{\unicode[STIX]{x03C0}^{\ast }f,\unicode[STIX]{x03C0}^{\ast }g\}_{M}$. For a Poisson system equipped with a group-equivariant momentum map $\unicode[STIX]{x1D707}$ satisfying (4.21) – as in the Vlasov–Maxwell system of interest – such a quotient map may be defined on level sets of $\unicode[STIX]{x1D707}$, as we now describe.

Consider the preimage of an arbitrary $\unicode[STIX]{x1D6FC}\in \mathfrak{f}^{\ast }$ under $\unicode[STIX]{x1D707}:M\rightarrow \mathfrak{f}^{\ast }$, that is, the level set $\unicode[STIX]{x1D707}^{-1}(\unicode[STIX]{x1D6FC})\subset M$. We may take equivalence classes of this preimage under the full action of $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D713}}$$\forall$$\unicode[STIX]{x1D713}\in {\mathcal{F}}$. That is, we reduce $\unicode[STIX]{x1D707}^{-1}(\unicode[STIX]{x1D6FC})$ to the quotient manifold $\unicode[STIX]{x1D707}^{-1}(\unicode[STIX]{x1D6FC})/{\mathcal{F}}$, and thereby take a ‘slice’ of the orbit of $\unicode[STIX]{x1D707}^{-1}(\unicode[STIX]{x1D6FC})$ under the action of ${\mathcal{F}}$. Because $\unicode[STIX]{x1D707}$ is equivariant in the sense of (4.21), this quotient is well defined. The reduced manifold $\unicode[STIX]{x1D707}^{-1}(\unicode[STIX]{x1D6FC})/{\mathcal{F}}$ will again be a Poisson manifold, as we now show.

Let us consider the particular case $\unicode[STIX]{x1D6FC}=0$, and define $M_{0}:=\unicode[STIX]{x1D707}^{-1}(0)$. By (4.20), $M_{0}$ corresponds to the submanifold of $M$ on which $\unicode[STIX]{x1D70C}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{Y}$. We now take equivalence classes of $M_{0}$ under the orbit of ${\mathcal{F}}$ by defining new phase space coordinates that are invariant under the action of (4.13), namely

(4.26)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\bar{f}(\boldsymbol{x},\boldsymbol{v})=f(\boldsymbol{x},\boldsymbol{p}=\boldsymbol{v}+\boldsymbol{A}),\\ \boldsymbol{B}=\unicode[STIX]{x1D735}\times \boldsymbol{A},\\ \boldsymbol{E}=-\boldsymbol{Y}.\end{array}\right\}\end{eqnarray}$$

We therefore identify the manifold of equivalence classes $\tilde{M}_{0}:=M_{0}/{\mathcal{F}}$ with the manifold $(\bar{f},\boldsymbol{B},\boldsymbol{E})$ of densities $\bar{f}$ defined on $(\boldsymbol{x},\boldsymbol{v})$ space, vector fields $\boldsymbol{B}$ that satisfy $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}=0$, and vector fields $\boldsymbol{E}$ that satisfy $\bar{\unicode[STIX]{x1D70C}}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{E}$, where now $\bar{\unicode[STIX]{x1D70C}}:=\int \bar{f}\,\text{d}\boldsymbol{v}$. (We note that the choice to constrain $M$ to $\{m\in M\mid \unicode[STIX]{x1D707}(m)=0\}$ evidently corresponds to the physical case $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{E}-\bar{\unicode[STIX]{x1D70C}}=0$, in which no ‘external’ charges are present in the system. Such a choice must be made when determining the Vlasov–Maxwell system’s initial conditions, as we shall see.) Our reduction map is therefore summarized by

(4.27)$$\begin{eqnarray}\unicode[STIX]{x03C0}_{\text{red}}:\quad \begin{array}{@{}rcl@{}}\unicode[STIX]{x1D707}^{-1}(0)\subset M\hspace{10.00002pt} & \longrightarrow \hspace{10.00002pt} & \tilde{M}_{0}:=\unicode[STIX]{x1D707}^{-1}(0)/{\mathcal{F}}\\ (f(\boldsymbol{x},\boldsymbol{p}),\boldsymbol{A},\boldsymbol{Y})\hspace{10.00002pt} & \longmapsto \hspace{10.00002pt} & (\bar{f}(\boldsymbol{x},\boldsymbol{v}),\boldsymbol{B},\boldsymbol{E}).\end{array}\end{eqnarray}$$

As calculated in Marsden & Weinstein (Reference Marsden and Weinstein1982) § 7, the substitution of (4.26) into the bracket of (4.1) yields the following reduced Poisson bracket on $\tilde{M}_{0}$:

(4.28)$$\begin{eqnarray}\displaystyle & & \displaystyle \{\{F,G\}\}^{\text{red}}[\bar{f},\boldsymbol{B},\boldsymbol{E}]=\int \text{d}\boldsymbol{x}\,\text{d}\boldsymbol{v}\left[\bar{f}\left\{\frac{\unicode[STIX]{x1D6FF}F}{\unicode[STIX]{x1D6FF}\bar{f}},\frac{\unicode[STIX]{x1D6FF}G}{\unicode[STIX]{x1D6FF}\bar{f}}\right\}_{\boldsymbol{x}\boldsymbol{v}}\right.\nonumber\\ \displaystyle & & \displaystyle \quad +\left.\bar{f}\boldsymbol{B}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{v}}\frac{\unicode[STIX]{x1D6FF}F}{\unicode[STIX]{x1D6FF}\bar{f}}\times \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{v}}\frac{\unicode[STIX]{x1D6FF}G}{\unicode[STIX]{x1D6FF}\bar{f}}\right)+\left(\frac{\unicode[STIX]{x1D6FF}F}{\unicode[STIX]{x1D6FF}\boldsymbol{E}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\bar{f}}{\unicode[STIX]{x2202}\boldsymbol{v}}\frac{\unicode[STIX]{x1D6FF}G}{\unicode[STIX]{x1D6FF}\bar{f}}-\frac{\unicode[STIX]{x1D6FF}G}{\unicode[STIX]{x1D6FF}\boldsymbol{E}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\bar{f}}{\unicode[STIX]{x2202}\boldsymbol{v}}\frac{\unicode[STIX]{x1D6FF}F}{\unicode[STIX]{x1D6FF}\bar{f}}\right)\right]\nonumber\\ \displaystyle & & \displaystyle \quad +\,\int \text{d}\boldsymbol{x}\left(\frac{\unicode[STIX]{x1D6FF}F}{\unicode[STIX]{x1D6FF}\boldsymbol{E}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \frac{\unicode[STIX]{x1D6FF}G}{\unicode[STIX]{x1D6FF}\boldsymbol{B}}-\frac{\unicode[STIX]{x1D6FF}G}{\unicode[STIX]{x1D6FF}\boldsymbol{E}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \frac{\unicode[STIX]{x1D6FF}F}{\unicode[STIX]{x1D6FF}\boldsymbol{B}}\right).\end{eqnarray}$$

We note that this process of Poisson reduction preserves the $\dot{\unicode[STIX]{x1D707}}=0$ conservation law associated with the unreduced Poisson manifold $M$. After all, the image of $\unicode[STIX]{x03C0}_{\text{red}}$ restricts $M$ to (quotients of) a submanifold $M_{0}\subset M$ on which $\unicode[STIX]{x1D707}$ is constant – in particular, level sets of a single value of $\unicode[STIX]{x1D707}$. The conservation law of (4.24) is clearly respected by this reduction, and may simply be re-expressed in the phase space variables of the reduced manifold $\tilde{M}_{0}$, along with its form in (4.25), i.e.

(4.29)$$\begin{eqnarray}\displaystyle 0=\dot{\bar{\unicode[STIX]{x1D707}}} & = & \displaystyle \dot{\bar{\unicode[STIX]{x1D70C}}}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\dot{\boldsymbol{E}}\nonumber\\ \displaystyle & = & \displaystyle \dot{\bar{\unicode[STIX]{x1D70C}}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\bar{\boldsymbol{J}},\end{eqnarray}$$

where

(4.30)$$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D707}} & = & \displaystyle \int \text{d}\boldsymbol{v}\,\bar{f}(\boldsymbol{x},\boldsymbol{v})-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{E}\nonumber\\ \displaystyle & = & \displaystyle \bar{\unicode[STIX]{x1D70C}}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{E},\end{eqnarray}$$

and where $\bar{\unicode[STIX]{x1D70C}}:=\int \bar{f}\,\text{d}\boldsymbol{v}$ and $\bar{\boldsymbol{J}}:=\int \bar{f}\boldsymbol{v}\,\text{d}\boldsymbol{v}$.

We note that the reduced bracket of (4.28) is a well-defined Poisson bracket specifically on the quotient submanifold $\tilde{M}_{0}$. Some of the plasma physics literature (e.g. Morrison Reference Morrison1982, Reference Morrison2013; Kraus et al. Reference Kraus, Kormann, Morrison and Sonnendrücker2017) notes that (4.28) generally fails to satisfy a Jacobi identity, however, so we pause to elucidate the source of this contrasting point of view.

In particular, the aforementioned literature defines the Vlasov–Maxwell system on an augmented manifold that includes all unconstrained vector fields $\boldsymbol{E},\boldsymbol{B}\in \mathbb{R}^{3}$:

(4.31)$$\begin{eqnarray}\tilde{M}_{0}^{+}:=\tilde{M}_{0}\sqcup \{\boldsymbol{E},\boldsymbol{B}\mid \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{E}\neq \bar{\unicode[STIX]{x1D70C}},\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}\neq 0\}.\end{eqnarray}$$

When the bracket of (4.28) is defined on $\tilde{M}_{0}^{+}$ and not on $\tilde{M}_{0}$, it no longer everywhere obeys the Jacobi identity (Morrison Reference Morrison1982; Chandre et al. Reference Chandre, Guillebon, Back, Tassi and Morrison2013); in particular, the Jacobi identity is satisfied on $\tilde{M}_{0}^{+}$ only when $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}=0$. Indeed, the constraint $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}=0$ appears as an exogenous defect that must be satisfied for $(\tilde{M}_{0}^{+},\{\{\cdot ,\cdot \}\}^{\text{red}})$ to be considered a Poisson manifold. On $\tilde{M}_{0}^{+}$, the bracket of (4.28) also acquires additional Casimirs,

(4.32)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\{\{\cdot ,\bar{\unicode[STIX]{x1D70C}}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{E}\}\}=0\\ \{\{\cdot ,\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}\}\}=0,\end{array}\right\}\end{eqnarray}$$

in much the same way that a Poisson structure on $\mathbb{R}^{2}=\{(x,y)\}$ acquires a $z$ Casimir when the system is embedded in $\mathbb{R}^{3}$.

We adopt the point of view that it is more natural to regard the bracket of (4.28) as a Poisson bracket defined on the submanifold of physical interest – $\tilde{M}_{0}=\unicode[STIX]{x1D707}^{-1}(0)/{\mathcal{F}}$ – rather than a defected bracket defined on the larger manifold including arbitrary vector fields $\boldsymbol{E}$ and $\boldsymbol{B}$. In a sense, it is merely a lack of economical notation that leads us to coordinatize $\tilde{M}_{0}$ with vector symbols $\boldsymbol{E}$ and $\boldsymbol{B}$ that are more commonly defined over all of $\mathbb{R}^{3}$.

It is clear from this discussion, however, that care must be taken in any numerical implementation of the reduced Vlasov–Maxwell bracket to constrain one’s fields to the phase space of $\tilde{M}_{0}$; generic, unconstrained vector fields $\boldsymbol{E},\boldsymbol{B}\in \mathbb{R}^{3}$ are to be avoided.

6.1 An ‘unreduced’ Hamiltonian PIC method

With the formalism we have developed, we proceed to explore PIC methods in the Hamiltonian setting, by defining a PIC splitting method adapted from Xiao et al. (Reference Xiao, Qin, Liu, He, Zhang and Sun2015) and Qin et al. (Reference Qin, Liu, Xiao, Zhang, He, Wang, Sun, Burby, Ellison and Zhou2016). The latter of these references implements a symplectic-Euler integrator for the unreduced Poisson bracket of (4.1), while the former implements a splitting method for the reduced bracket of (4.28). We shall synthesize the two, defining a gauge-compatible splitting method for the unreduced bracket of (4.1) and, in so doing, demonstrate the merit of this new class of splitting methods. The result is an explicit-time-advance, canonical, locally charge-conserving PIC method, whose momentum map and conservation law we shall systematically derive.

In Qin et al. (Reference Qin, Liu, Xiao, Zhang, He, Wang, Sun, Burby, Ellison and Zhou2016), a Klimontovich–Maxwell PIC method is derived from the unreduced bracket of (4.1) by specifying the following form for the distribution function $f(\boldsymbol{x},\boldsymbol{p})$ of $L$ particles, analogous to (3.1)

(6.1)$$\begin{eqnarray}f(\boldsymbol{x},\boldsymbol{p})=\mathop{\sum }_{i=1}^{L}\unicode[STIX]{x1D6FF}^{(3)}(\boldsymbol{x}-\boldsymbol{X}_{i})\unicode[STIX]{x1D6FF}^{(3)}(\boldsymbol{p}-\boldsymbol{P}_{i}).\end{eqnarray}$$

Here, $(\boldsymbol{X}_{i},\boldsymbol{P}_{i})$ denotes the dynamical coordinates of particle $i$ in phase space. The fields $\boldsymbol{A}$ and $\boldsymbol{Y}$ are also discretized on a (three-dimensional) spatial lattice and are denoted $(\boldsymbol{A}_{n},\boldsymbol{Y}_{n})$ at lattice site $n$. We shall further require the interpolation of $\boldsymbol{A}$

(6.2)$$\begin{eqnarray}\boldsymbol{A}(\boldsymbol{x})=\mathop{\sum }_{n=1}^{N}\boldsymbol{A}_{n}W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{x}-\boldsymbol{x}_{n}),\end{eqnarray}$$

where $n$ is an index over all $N$ lattice sites and $W_{\unicode[STIX]{x1D70E}_{1}}$ is an (as yet unspecified) interpolation function for $\boldsymbol{A}$.

A Poisson bracket for this discrete system simply follows from the canonical symplectic structure of its variables. In particular, we define the symplectic manifold

(6.3)$$\begin{eqnarray}M_{d}=T^{\ast }X\times T^{\ast }Q,\end{eqnarray}$$

where $X=\mathbb{R}^{3L}$ is the space of particle position coordinates and $Q=\mathbb{R}^{3N}$ is the space of vector potentials on the lattice, such that $T^{\ast }X=\{(\boldsymbol{X}_{i},\boldsymbol{P}_{i})\}$ and $T^{\ast }Q=\{(\boldsymbol{A}_{n},\boldsymbol{Y}_{n})\}$. A point $m\in M_{d}$ correspondingly specifies $(\boldsymbol{X}_{i},\boldsymbol{P}_{i},\boldsymbol{A}_{n},\boldsymbol{Y}_{n})$$\forall$$i,n$ (where the subscript $d$ denotes discretization).

The Poisson bracket for this symplectic manifold therefore takes its usual Darboux-coordinate form

(6.4)$$\begin{eqnarray}\displaystyle \{\{F,G\}\}_{d}[\boldsymbol{X}_{i},\boldsymbol{P}_{i},\boldsymbol{A}_{n},\boldsymbol{Y}_{n}] & = & \displaystyle \mathop{\sum }_{i=1}^{L}\left(\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\boldsymbol{X}_{i}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}G}{\unicode[STIX]{x2202}\boldsymbol{P}_{i}}-\frac{\unicode[STIX]{x2202}G}{\unicode[STIX]{x2202}\boldsymbol{X}_{i}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\boldsymbol{P}_{i}}\right)\nonumber\\ \displaystyle & & \displaystyle +\,\mathop{\sum }_{n=1}^{N}\left(\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\boldsymbol{A}_{n}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}G}{\unicode[STIX]{x2202}\boldsymbol{Y}_{n}}-\frac{\unicode[STIX]{x2202}G}{\unicode[STIX]{x2202}\boldsymbol{A}_{n}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\boldsymbol{Y}_{n}}\right).\end{eqnarray}$$

We observe that, unlike its continuous counterpart in (4.1), the bracket of (6.4) is non-degenerate; it defines $M_{d}$ not only as a Poisson manifold, but as a symplectic manifold.

The discrete Hamiltonian of Qin et al. (Reference Qin, Liu, Xiao, Zhang, He, Wang, Sun, Burby, Ellison and Zhou2016) is derived from (4.2) by substituting the Klimontovich distribution of (6.1) and expanding terms of the form $|\boldsymbol{P}_{i}-\boldsymbol{A}(\boldsymbol{X}_{i})|^{2}$ using (6.2)

(6.5)$$\begin{eqnarray}\displaystyle & & \displaystyle H_{d}[\boldsymbol{X}_{i},\boldsymbol{P}_{i},\boldsymbol{A}_{n},\boldsymbol{Y}_{n}]=\frac{1}{2}\mathop{\sum }_{i=1}^{L}\left[\boldsymbol{P}_{i}^{2}-2\boldsymbol{P}_{i}\boldsymbol{\cdot }\mathop{\sum }_{n=1}^{N}\boldsymbol{A}_{n}W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{X}_{i}-\boldsymbol{x}_{n})\right.\nonumber\\ \displaystyle & & \displaystyle \quad +\left.\mathop{\sum }_{m,n=1}^{N}\boldsymbol{A}_{m}\boldsymbol{\cdot }\boldsymbol{A}_{n}W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{X}_{i}-\boldsymbol{x}_{m})W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{X}_{i}-\boldsymbol{x}_{n})\right]+\frac{1}{2}\mathop{\sum }_{n=1}^{N}[\boldsymbol{Y}_{n}^{2}+|\unicode[STIX]{x1D735}_{d}^{+}\times \boldsymbol{A}|_{n}^{2}].\end{eqnarray}$$

Here, the operator $(\unicode[STIX]{x1D735}_{d}^{\pm }\times )_{n}$ represents a discrete curl, defined by

(6.6)$$\begin{eqnarray}(\unicode[STIX]{x1D735}_{d}^{\pm }\times \boldsymbol{A})_{n}:=\pm \left(\begin{array}{@{}c@{}}\displaystyle \frac{A_{i,j\pm 1,k}^{3}-A_{i,j,k}^{3}}{\unicode[STIX]{x0394}y}-\frac{A_{i,j,k\pm 1}^{2}-A_{i,j,k}^{2}}{\unicode[STIX]{x0394}z}\\ \displaystyle \frac{A_{i,j,k\pm 1}^{1}-A_{i,j,k}^{1}}{\unicode[STIX]{x0394}z}-\frac{A_{i\pm 1,j,k}^{3}-A_{i,j,k}^{3}}{\unicode[STIX]{x0394}x}\\ \displaystyle \frac{A_{i\pm 1,j,k}^{2}-A_{i,j,k}^{2}}{\unicode[STIX]{x0394}x}-\frac{A_{i,j\pm 1,k}^{1}-A_{i,j,k}^{1}}{\unicode[STIX]{x0394}y}\end{array}\right)\end{eqnarray}$$

for $n=(i,j,k)$.

We now describe the gauge symmetry of this discrete Hamiltonian system. We define the group action $\unicode[STIX]{x1D6F7}_{f}$ on $M_{d}$ by analogy with (4.13)

(6.7)$$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{f}(\boldsymbol{X}_{i},\boldsymbol{P}_{i},\boldsymbol{A}_{n},\boldsymbol{Y}_{n})=(\boldsymbol{X}_{i},[\boldsymbol{P}_{i}-\unicode[STIX]{x1D735}_{d}^{+}f(\boldsymbol{X}_{i})],[\boldsymbol{A}_{n}-(\unicode[STIX]{x1D735}_{d}^{+}f)_{n}],\boldsymbol{Y}_{n}),\end{eqnarray}$$

where

(6.8)$$\begin{eqnarray}\unicode[STIX]{x1D735}_{d}^{+}f(\boldsymbol{x})=\mathop{\sum }_{n=1}^{N}(\unicode[STIX]{x1D735}_{d}^{+}f)_{n}W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{x}-\boldsymbol{x}_{n})\end{eqnarray}$$

and where $(\unicode[STIX]{x1D735}_{d}^{\pm })_{n}$ is a discrete gradient defined by

(6.9)$$\begin{eqnarray}(\unicode[STIX]{x1D735}_{d}^{\pm }f)_{n}:=\pm \left(\begin{array}{@{}c@{}}\displaystyle \frac{f_{i\pm 1,j,k}-f_{i,j,k}}{\unicode[STIX]{x0394}x}\\ \displaystyle \frac{f_{i,j\pm 1,k}-f_{i,j,k}}{\unicode[STIX]{x0394}y}\\ \displaystyle \frac{f_{i,j,k\pm 1}-f_{i,j,k}}{\unicode[STIX]{x0394}z}\end{array}\right).\end{eqnarray}$$

We note that $\unicode[STIX]{x1D735}_{d}^{\pm }\times \unicode[STIX]{x1D735}_{d}^{\pm }=0$ as an operator. (If the $\pm$ signs agree, this relation holds identically; if they disagree, it holds only after a summation over lattice points, $\sum _{n}$.) We also note that – in contrast with (4.13)–(4.14) – $\boldsymbol{P}_{i}$ and $\boldsymbol{A}_{n}$ are shifted in the same direction in (6.7), reflecting the fact that $\boldsymbol{p}$ and $\boldsymbol{P}_{i}$ have opposite signs in (6.1) when we reinterpret the transformation of (4.14) as a transformation of $\boldsymbol{P}_{i}$.

The function $f$ appearing in the group action of (6.7) is to be understood as a scalar function defined only at lattice points. In particular, $f\in {\mathcal{F}}_{d}$ is a group element of the set ${\mathcal{F}}_{d}$ of discrete scalar functions with an abelian composition law of addition. Its Lie algebra $\mathfrak{f}_{d}$ is also the set of discrete scalar functions on the lattice, while its dual $\mathfrak{f}_{d}^{\ast }$ is the set of densities, which pair to elements of $\mathfrak{f}_{d}$ by summing over pointwise products

(6.10)$$\begin{eqnarray}\langle \unicode[STIX]{x1D6FC},\unicode[STIX]{x1D719}\rangle :=\mathop{\sum }_{n=1}^{N}\unicode[STIX]{x1D6FC}_{n}\unicode[STIX]{x1D719}_{n}\quad \forall ~\unicode[STIX]{x1D6FC}\in \mathfrak{f}_{d}^{\ast },\unicode[STIX]{x1D719}\in \mathfrak{f}_{d}.\end{eqnarray}$$

We must verify that the group action is canonical, a task most easily approached infinitesimally. In particular, we ask whether the following infinitesimal form of $\{\{\unicode[STIX]{x1D6F7}_{f}^{\ast }F,\unicode[STIX]{x1D6F7}_{f}^{\ast }G\}\}_{d}=\unicode[STIX]{x1D6F7}_{f}^{\ast }\{\{F,G\}\}_{d}$ holds:

(6.11)$$\begin{eqnarray}\displaystyle & & \displaystyle \left\{\left\{-\unicode[STIX]{x1D735}_{d}^{+}\unicode[STIX]{x1D719}(\boldsymbol{X}_{i})\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\boldsymbol{P}_{i}}-\unicode[STIX]{x1D735}_{d}^{+}\unicode[STIX]{x1D719}_{n}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\boldsymbol{A}_{n}},G\right\}\right\}_{d}-(F\leftrightarrow G)\nonumber\\ \displaystyle & & \displaystyle \quad =-\unicode[STIX]{x1D735}_{d}^{+}\unicode[STIX]{x1D719}(\boldsymbol{X}_{i})\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\{\{F,G\}\}_{d}}{\unicode[STIX]{x2202}\boldsymbol{P}_{i}}-\unicode[STIX]{x1D735}_{d}^{+}\unicode[STIX]{x1D719}_{n}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\{\{F,G\}\}_{d}}{\unicode[STIX]{x2202}\boldsymbol{A}_{n}},\end{eqnarray}$$

where summation over repeated indices is implicit. After applying (6.4) to evaluate each bracket, equation (6.11) is seen to be true only when $\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D735}_{d}^{+}\unicode[STIX]{x1D719}(\boldsymbol{X}_{i})=0$. This requires the operator relation

(6.12)$$\begin{eqnarray}\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D735}_{d}^{+}=0.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D735}\equiv \unicode[STIX]{x2202}_{\boldsymbol{X}_{i}}$ is a continuous spatial gradient.

Equation (6.12) therefore necessitates the following condition on the interpolation function $W_{\unicode[STIX]{x1D70E}_{1}}$:

(6.13)$$\begin{eqnarray}\mathop{\sum }_{n=1}^{N}(\unicode[STIX]{x1D735}_{d}^{+}\unicode[STIX]{x1D719})_{n}W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{x}-\boldsymbol{x}_{n})=\unicode[STIX]{x1D735}\mathop{\sum }_{n=1}^{N}\unicode[STIX]{x1D719}_{n}W_{\unicode[STIX]{x1D70E}_{0}}(\boldsymbol{x}-\boldsymbol{x}_{n})\end{eqnarray}$$

for some interpolation function $W_{\unicode[STIX]{x1D70E}_{0}}$. This condition was already discovered in Xiao et al. (Reference Xiao, Qin, Liu, He, Zhang and Sun2015), and is analogous to a property of the simplicial Whitney forms described earlier (Bossavit Reference Bossavit1988). Our discussion of this condition merely contributes that, in a Hamiltonian context, the motivation for the constraint in (6.13) is the canonicality of the group action.

We now solve for $\unicode[STIX]{x1D707}_{d}$, the momentum map on $M_{d}$ associated with the group action of (6.7), using the symplectic structure of (6.4). First, we must find the infinitesimal generator $\unicode[STIX]{x1D719}_{M_{d}}$ of our group action on $M_{d}$, defined analogously to (4.15). Given the group action of (6.7) we expect $\unicode[STIX]{x1D719}_{M_{d}}$ to take the form (already implicitly used in (6.11))

(6.14)$$\begin{eqnarray}\{\{\cdot ,\unicode[STIX]{x1D707}_{d}^{\unicode[STIX]{x1D719}}\}\}_{d}=-\mathop{\sum }_{i=1}^{L}\unicode[STIX]{x1D735}_{d}^{+}\unicode[STIX]{x1D719}(\boldsymbol{X}_{i})\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{P}_{i}}-\mathop{\sum }_{n=1}^{N}(\unicode[STIX]{x1D735}_{d}^{+}\unicode[STIX]{x1D719})_{n}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{A}_{n}},\end{eqnarray}$$

where we denote the pairing of the momentum map with $\unicode[STIX]{x1D719}$ by $\langle \unicode[STIX]{x1D707}_{d},\unicode[STIX]{x1D719}\rangle =\unicode[STIX]{x1D707}_{d}^{\unicode[STIX]{x1D719}}$. The Poisson bracket of (6.4) therefore requires that $\unicode[STIX]{x1D707}_{d}^{\unicode[STIX]{x1D719}}$ be given by

(6.15)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}_{d}^{\unicode[STIX]{x1D719}}(m) & = & \displaystyle \mathop{\sum }_{n=1}^{N}(\unicode[STIX]{x1D735}_{d}^{+}\unicode[STIX]{x1D719})_{n}\boldsymbol{\cdot }\left[\mathop{\sum }_{i=1}^{L}\int _{-\infty }^{\boldsymbol{X}_{i}}\text{d}\boldsymbol{X}_{i}^{\prime }\,W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{X}_{i}^{\prime }-\boldsymbol{x}_{n})-\boldsymbol{Y}_{n}\right]\nonumber\\ \displaystyle & = & \displaystyle \mathop{\sum }_{n=1}^{N}\unicode[STIX]{x1D719}_{n}\unicode[STIX]{x1D735}_{d}^{-}\boldsymbol{\cdot }\left[-\mathop{\sum }_{i=1}^{L}\int _{-\infty }^{\boldsymbol{X}_{i}}\text{d}\boldsymbol{X}_{i}^{\prime }\,W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{X}_{i}^{\prime }-\boldsymbol{x}_{n})+\boldsymbol{Y}_{n}\right],\end{eqnarray}$$

where in the second line we have summed by parts (Hydon & Mansfield Reference Hydon and Mansfield2011) using the discrete divergence operator

(6.16)$$\begin{eqnarray}\unicode[STIX]{x1D735}_{d}^{\pm }\boldsymbol{\cdot }\boldsymbol{v}_{n}:=\pm \mathop{\sum }_{\unicode[STIX]{x1D6FC}=1}^{3}\frac{v_{n\pm \hat{\unicode[STIX]{x1D6FC}}}^{\unicode[STIX]{x1D6FC}}-v_{n}^{\unicode[STIX]{x1D6FC}}}{\unicode[STIX]{x0394}x^{\unicode[STIX]{x1D6FC}}}.\end{eqnarray}$$

Note that $\text{d}\boldsymbol{X}_{i}^{\prime }$ is treated in (6.15) and hereafter as a vector, with each component integrated individually. We observe that $\unicode[STIX]{x1D735}_{d}^{\pm }\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{d}^{\pm }\times =0$ as an operator (when $\pm$ signs agree).

Given the pairing defined in (6.10), the momentum map $\unicode[STIX]{x1D707}_{d}$ must therefore be

(6.17)$$\begin{eqnarray}(\unicode[STIX]{x1D707}_{d}(m))_{n}=-\unicode[STIX]{x1D735}_{d}^{-}\boldsymbol{\cdot }\mathop{\sum }_{i=1}^{L}\int _{-\infty }^{\boldsymbol{X}_{i}}\text{d}\boldsymbol{X}_{i}^{\prime }\,W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{X}_{i}^{\prime }-\boldsymbol{x}_{n})+\unicode[STIX]{x1D735}_{d}^{-}\boldsymbol{\cdot }\boldsymbol{Y}_{n}\end{eqnarray}$$

defined at each lattice site $n\in [1,N]$. Due to the gauge invariance of $H_{d}$ in (6.5) – that is, $\unicode[STIX]{x1D6F7}_{f}^{\ast }H_{d}=H_{d}$ – the full system evolved in continuous time by $H_{d}$ obeys the conservation law

(6.18)$$\begin{eqnarray}\dot{\unicode[STIX]{x1D707}}_{d}=0,\end{eqnarray}$$

as in the continuous Vlasov–Maxwell system of § 4. Equations (6.17)–(6.18) define the conservation law of our discrete Hamiltonian system in continuous time, systematically derived via the momentum map.

Following the analysis of (4.24)–(4.25), we may re-express this conservation law by deriving the continuous-time EOM of the full Hamiltonian $H_{d}$, as follows:

(6.19)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \dot{\boldsymbol{X}}_{i}=\{\{\boldsymbol{X}_{i},{H_{d}\}\}}_{d}=\boldsymbol{P}_{i}-\mathop{\sum }_{m=1}^{N}\boldsymbol{A}_{m}W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{X}_{i}-\boldsymbol{x}_{m}),\\ \displaystyle \dot{\boldsymbol{P}}_{i}=\{\{\boldsymbol{P}_{i},{H_{d}\}\}}_{d}=\mathop{\sum }_{m=1}^{N}(\dot{\boldsymbol{X}}_{i}\boldsymbol{\cdot }\boldsymbol{A}_{m})\unicode[STIX]{x1D735}W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{X}_{i}-\boldsymbol{x}_{m}),\\ \displaystyle \dot{\boldsymbol{A}}_{n}=\{\{\boldsymbol{A}_{n},{H_{d}\}\}}_{d}=\boldsymbol{Y}_{n},\\ \displaystyle \dot{\boldsymbol{Y}}_{n}=\{\{\boldsymbol{Y}_{n},{H_{d}\}\}}_{d}=\mathop{\sum }_{i=1}^{L}\dot{\boldsymbol{X}}_{i}W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{X}_{i}-\boldsymbol{x}_{n})-(\unicode[STIX]{x1D735}_{d}^{-}\times \unicode[STIX]{x1D735}_{d}^{+}\times \boldsymbol{A})_{n}.\end{array}\right\}\end{eqnarray}$$

Now substituting $\dot{\boldsymbol{Y}}_{n}$ into the charge conservation law equations (6.17)–(6.18), we note that

(6.20)$$\begin{eqnarray}0=\dot{\unicode[STIX]{x1D70C}}_{n}+\unicode[STIX]{x1D735}_{d}^{-}\boldsymbol{\cdot }\boldsymbol{J}_{n},\end{eqnarray}$$

where

(6.21)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D70C}_{n}:=-\unicode[STIX]{x1D735}_{d}^{-}\boldsymbol{\cdot }\mathop{\sum }_{i=1}^{L}\int _{-\infty }^{\boldsymbol{X}_{i}}\text{d}\boldsymbol{X}_{i}^{\prime }\,W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{X}_{i}^{\prime }-\boldsymbol{x}_{n}),\\ \displaystyle \boldsymbol{J}_{n}:=\mathop{\sum }_{i=1}^{L}\dot{\boldsymbol{X}}_{i}W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{X}_{i}-\boldsymbol{x}_{n}).\end{array}\right\}\end{eqnarray}$$

This is an alternative form of the charge conservation law equation (6.18) for the continuous-time evolution of the Hamiltonian system of Qin et al. (Reference Qin, Liu, Xiao, Zhang, He, Wang, Sun, Burby, Ellison and Zhou2016). (We observe that, unlike its counterpart in (4.25), it is an off-shell identity.)

The form of $\unicode[STIX]{x1D70C}_{n}$ in (6.21) can be justified by a schematic one-dimensional example in which $W_{\unicode[STIX]{x1D70E}_{1}}(x)=1$ on $0\leqslant x<\unicode[STIX]{x0394}x$ and 0 otherwise. For a single particle at $X_{i}=0.2$, we have

(6.22)$$\begin{eqnarray}\unicode[STIX]{x1D70C}_{n}=-\unicode[STIX]{x1D735}_{d}^{-}\boldsymbol{\cdot }\int _{-\infty }^{0.2}\text{d}X_{i}^{\prime }\,W_{\unicode[STIX]{x1D70E}_{1}}(X_{i}^{\prime }-x_{n})=\left\{\begin{array}{@{}ll@{}}0.8/\unicode[STIX]{x0394}x\quad & n=0,\\ 0.2/\unicode[STIX]{x0394}x\quad & n=1,\\ 0\quad & n\neq 0,1.\end{array}\right.\end{eqnarray}$$

This result demonstrates the appropriateness of the momentum map’s systematically derived charge density.

We now define an algorithmic solution of this Hamiltonian system via a splitting method, and examine the preservation of $\unicode[STIX]{x1D707}_{d}$. To algorithmically evolve this system in discrete time, we define a gauge-compatible splitting method adapted from He et al. (Reference He, Qin, Sun, Xiao, Zhang and Liu2015, Reference He, Sun, Qin and Liu2016). We define Hamiltonian subsystems

(6.23)$$\begin{eqnarray}H_{d}=\mathop{\sum }_{\unicode[STIX]{x1D6FC}=1}^{3}H_{\text{Klim}}^{\unicode[STIX]{x1D6FC}}+H_{\boldsymbol{ A}}+H_{\boldsymbol{Y}},\end{eqnarray}$$

where

(6.24)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle H_{\text{Klim}}^{\unicode[STIX]{x1D6FC}}:=\frac{1}{2}\mathop{\sum }_{i=1}^{L}(P_{i}^{\unicode[STIX]{x1D6FC}}-A^{\unicode[STIX]{x1D6FC}}(\boldsymbol{X}_{i}))^{2},\\ \displaystyle H_{\boldsymbol{A}}:=\frac{1}{2}\mathop{\sum }_{n=1}^{N}|\unicode[STIX]{x1D735}_{d}^{+}\times \boldsymbol{A}|_{n}^{2},\\ \displaystyle H_{\boldsymbol{Y}}:=\frac{1}{2}\mathop{\sum }_{n=1}^{N}\boldsymbol{Y}_{n}^{2}.\end{array}\right\}\end{eqnarray}$$

We immediately note that these subsystems are all gauge invariant for the group action of (6.7) – $\unicode[STIX]{x1D6F7}_{f}^{\ast }H_{i}=H_{i}$$\forall f\in {\mathcal{F}}_{d}$ and $i$ – and will therefore comprise a gauge-compatible splitting – and preserve $\unicode[STIX]{x1D707}_{d}$ – if they can be exactly solved.

Let us examine the EOM for each subsystem $H_{i}$ in turn

(6.25)

(6.26)

(6.27)

where $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D6FD}}\equiv \unicode[STIX]{x2202}/\unicode[STIX]{x2202}X_{i}^{\unicode[STIX]{x1D6FD}}$. (We emphasize that $\unicode[STIX]{x1D6FC}$ is fixed, and is not summed over in the expressions for $H_{\text{Klim}}^{\unicode[STIX]{x1D6FC}}$.) $H_{\boldsymbol{A}}$ and $H_{\boldsymbol{Y}}$ are exactly solvable at a glance. Furthermore, $H_{\text{Klim}}^{\unicode[STIX]{x1D6FC}}$ is seen to be exactly solvable by noting that $\ddot{X}_{i}^{\unicode[STIX]{x1D6FD}}=0$; ${\dot{X}}_{i}^{\unicode[STIX]{x1D6FD}}$ is therefore a constant determined by a time step’s initial conditions. The evolutions of $\dot{\boldsymbol{P}}_{i}$ and $\dot{\boldsymbol{Y}}_{n}$ in $H_{\text{Klim}}^{\unicode[STIX]{x1D6FC}}$ follow immediately from this analysis.

The exact time evolutions of $H_{\boldsymbol{A}}$, $H_{\boldsymbol{Y}}$ and $H_{\text{Klim}}^{\unicode[STIX]{x1D6FC}}$ are therefore explicitly solved, defining by construction an explicit-time-advance gauge-compatible splitting method that exactly preserves the momentum map, $\dot{\unicode[STIX]{x1D707}}_{d}=0$, as desired. The alternative form of the charge conservation law given in (6.20) – that is, $\dot{\unicode[STIX]{x1D70C}}_{n}+\unicode[STIX]{x1D735}_{d}^{-}\boldsymbol{\cdot }\boldsymbol{J}_{n}=0$ – is also exactly preserved in this algorithm, because the substitution that led from (6.18) to (6.20) – that is, $\unicode[STIX]{x1D735}_{d}^{-}\boldsymbol{\cdot }\dot{\boldsymbol{Y}}_{n}=\unicode[STIX]{x1D735}_{d}^{-}\boldsymbol{\cdot }\boldsymbol{J}_{n}$ – holds for each Hamiltonian subsystem above.

Finally, we note that the momentum map $\unicode[STIX]{x1D707}_{d}$ has significant ramifications for the appropriate initial conditions of the preceding algorithm. We refer the reader to a brief but important discussion of these initial conditions in the text following (6.30) below.

6.2 A ‘reduced’ Hamiltonian PIC method

We now examine the PIC method of Xiao et al. (Reference Xiao, Qin, Liu, He, Zhang and Sun2015), which employs a splitting method equivalent to that of the preceding section for the reduced Vlasov–Maxwell bracket of (4.28).

We will mirror Xiao et al. (Reference Xiao, Qin, Liu, He, Zhang and Sun2015) and derive this PIC scheme by undertaking the symplectic reduction (Marsden & Weinstein Reference Marsden and Weinstein1974) of the discrete canonical bracket defined in (6.4). As in § 4.4, we define a mapping to the reduced symplectic manifold $\tilde{M}_{d_{0}}=\unicode[STIX]{x1D707}_{d}^{-1}(0)/{\mathcal{F}}_{d}$, with coordinates given by

(6.28)$$\begin{eqnarray}\unicode[STIX]{x03C0}_{d,\text{red}}:\quad \begin{array}{@{}rcl@{}}\unicode[STIX]{x1D707}_{d}^{-1}(0)\subset M_{d}\hspace{10.00002pt} & \longrightarrow \hspace{10.00002pt} & \tilde{M}_{d_{0}}=\unicode[STIX]{x1D707}_{d}^{-1}(0)/{\mathcal{F}}_{d}\\ (\boldsymbol{X}_{i},\boldsymbol{P}_{i},\boldsymbol{A}_{n},\boldsymbol{Y}_{n})\hspace{10.00002pt} & \longmapsto \hspace{10.00002pt} & (\boldsymbol{X}_{i},\boldsymbol{V}_{i},\boldsymbol{B}_{n},\boldsymbol{E}_{n}),\end{array}\end{eqnarray}$$

where

(6.29)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\boldsymbol{X}_{i}=\boldsymbol{X}_{i},\\ \boldsymbol{V}_{i}=\boldsymbol{P}_{i}-\boldsymbol{A}(\boldsymbol{X}_{i}),\\ \boldsymbol{B}_{n}=(\unicode[STIX]{x1D735}_{d}^{+}\times \boldsymbol{A})_{n},\\ \boldsymbol{E}_{n}=-\boldsymbol{Y}_{n}.\end{array}\right\}\end{eqnarray}$$

As discussed earlier, care must be taken to ensure that the discrete fields $\boldsymbol{B}_{n}$ and $\boldsymbol{E}_{n}$ of $\tilde{M}_{d_{0}}$ obey the reduced manifold constraints

(6.30)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle (\unicode[STIX]{x1D735}_{d}^{+}\boldsymbol{\cdot }\boldsymbol{B})_{n}=0,\\ \displaystyle -\unicode[STIX]{x1D735}_{d}^{-}\boldsymbol{\cdot }\mathop{\sum }_{i=1}^{L}\int _{-\infty }^{\boldsymbol{X}_{i}}\text{d}\boldsymbol{X}_{i}^{\prime }W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{X}_{i}^{\prime }-\boldsymbol{x}_{n})-(\unicode[STIX]{x1D735}_{d}^{-}\boldsymbol{\cdot }\boldsymbol{E})_{n}=0.\end{array}\right\}\end{eqnarray}$$

We note that these constraints must also be satisfied by any initial condition of the algorithm. The former condition is necessary to enforce a physically valid magnetic field. If the latter condition (Gauss’s law) is not satisfied initially, it will have the effect of adding fixed ‘external’ charges at the corresponding vertex $n$. In particular, a non-zero initial Gauss’s law condition will evolve the system along some other reduced manifold $\tilde{M}_{d_{\unicode[STIX]{x1D6FC}}}=\unicode[STIX]{x1D707}_{d}^{-1}(\unicode[STIX]{x1D6FC})/{\mathcal{F}}_{d}$ with fixed external charge density $\unicode[STIX]{x1D6FC}$.

A similar initial condition must be determined for the unreduced algorithm of § 6.1 as well. The unreduced algorithm enforces the constraint $(\unicode[STIX]{x1D735}_{d}^{+}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{d}^{+}\times \boldsymbol{A})_{n}=0$ automatically. However, for simulations without external charges, care should be taken so that the value of $(\unicode[STIX]{x1D707}_{d}(m))_{n}$ in (6.17) is everywhere initialized to zero. (Alternatively, equation (6.17) can be used to properly initialize a simulation with external charges that remain fixed for all time.) We note that our derivation of the momentum map is essential to this correct specification of initial conditions.

We therefore proceed with the reduction of our discrete system and substitute equation (6.29) into the bracket of (6.4) to find

(6.31)$$\begin{eqnarray}\displaystyle & & \displaystyle \{\{F,G\}\}_{d}^{\text{red}}[\boldsymbol{X}_{i},\boldsymbol{V}_{i},\boldsymbol{B}_{n},\boldsymbol{E}_{n}]\nonumber\\ \displaystyle & & \displaystyle \quad =\mathop{\sum }_{i=1}^{L}\left(\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\boldsymbol{X}_{i}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}G}{\unicode[STIX]{x2202}\boldsymbol{V}_{i}}-\frac{\unicode[STIX]{x2202}G}{\unicode[STIX]{x2202}\boldsymbol{X}_{i}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\boldsymbol{V}_{i}}+\left[\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\boldsymbol{V}_{i}}\times \frac{\unicode[STIX]{x2202}G}{\unicode[STIX]{x2202}\boldsymbol{V}_{i}}\right]\boldsymbol{\cdot }\mathop{\sum }_{n=1}^{N}\boldsymbol{B}_{n}W_{\unicode[STIX]{x1D70E}_{2}}(\boldsymbol{X}_{i}-\boldsymbol{x}_{n})\right)\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\mathop{\sum }_{n=1}^{N}\left(\left[\mathop{\sum }_{i=1}^{L}\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\boldsymbol{V}_{i}}W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{X}_{i}-\boldsymbol{x}_{n})-\left(\unicode[STIX]{x1D735}_{d}^{-}\times \frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\boldsymbol{B}}\right)_{n}\right]\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}G}{\unicode[STIX]{x2202}\boldsymbol{E}_{n}}\right.\nonumber\\ \displaystyle & & \displaystyle \qquad -\left.\left[\mathop{\sum }_{i=1}^{L}\frac{\unicode[STIX]{x2202}G}{\unicode[STIX]{x2202}\boldsymbol{V}_{i}}W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{X}_{i}-\boldsymbol{x}_{n})-\left(\unicode[STIX]{x1D735}_{d}^{-}\times \frac{\unicode[STIX]{x2202}G}{\unicode[STIX]{x2202}\boldsymbol{B}}\right)_{n}\right]\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\boldsymbol{E}_{n}}\right).\end{eqnarray}$$

To derive the $\unicode[STIX]{x2202}_{\boldsymbol{V}_{i}}F\times \unicode[STIX]{x2202}_{\boldsymbol{V}_{i}}G\boldsymbol{\cdot }\boldsymbol{B}(\boldsymbol{X}_{i})$ term in the bracket above, our interpolation functions were required to satisfy an additional constraint

(6.32)$$\begin{eqnarray}\unicode[STIX]{x1D735}\times \mathop{\sum }_{n=1}^{N}\boldsymbol{A}_{n}W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{x}-\boldsymbol{x}_{n})=\mathop{\sum }_{n=1}^{N}(\unicode[STIX]{x1D735}_{d}^{+}\times \boldsymbol{A})_{n}W_{\unicode[STIX]{x1D70E}_{2}}(\boldsymbol{x}-\boldsymbol{x}_{n})\end{eqnarray}$$

for some interpolation function $W_{\unicode[STIX]{x1D70E}_{2}}$. As in (6.13), this is a generalized, higher-dimensional analogue of the simplicial Whitney interpolant constraint.

Lastly, we re-express the Hamiltonian in the reduced coordinates of $\tilde{M}_{d_{0}}$ as

(6.33)$$\begin{eqnarray}H_{d}^{\text{red}}[\boldsymbol{X}_{i},\boldsymbol{V}_{i},\boldsymbol{B}_{n},\boldsymbol{E}_{n}]=\frac{1}{2}\mathop{\sum }_{i=1}^{L}\boldsymbol{V}_{i}^{2}+\frac{1}{2}\mathop{\sum }_{n=1}^{N}(\boldsymbol{E}_{n}^{2}+\boldsymbol{B}_{n}^{2}).\end{eqnarray}$$

We have thus recovered the reduced Hamiltonian system of Xiao et al. (Reference Xiao, Qin, Liu, He, Zhang and Sun2015).

As discussed in § 4.4 for spatially continuous systems, this reduced Hamiltonian system is automatically guaranteed to preserve the momentum map of its parent, as long as it evolution is constrained to $\tilde{M}_{d_{0}}$. To see that this is the case, we may compute its evolution equations under the splitting scheme analogous to the unreduced case (He et al. Reference He, Qin, Sun, Xiao, Zhang and Liu2015, Reference He, Sun, Qin and Liu2016; Xiao et al. Reference Xiao, Qin, Liu, He, Zhang and Sun2015)

(6.34)$$\begin{eqnarray}H_{d}^{\text{red}}=\mathop{\sum }_{\unicode[STIX]{x1D6FC}=1}^{3}H_{\boldsymbol{ V}}^{\unicode[STIX]{x1D6FC}}+H_{\boldsymbol{ B}}+H_{\boldsymbol{E}},\end{eqnarray}$$

where

(6.35)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle H_{\boldsymbol{V}}^{\unicode[STIX]{x1D6FC}}:=\frac{1}{2}\mathop{\sum }_{i=1}^{L}(V_{i}^{\unicode[STIX]{x1D6FC}})^{2},\\ \displaystyle H_{\boldsymbol{B}}:=\frac{1}{2}\mathop{\sum }_{n=1}^{N}\boldsymbol{B}_{n}^{2},\\ \displaystyle H_{\boldsymbol{E}}:=\frac{1}{2}\mathop{\sum }_{n=1}^{N}\boldsymbol{E}_{n}^{2}.\end{array}\right\}\end{eqnarray}$$

These subsystems generate the following EOM:

(6.36)

(6.37)

(6.38)

We note again that $\unicode[STIX]{x1D6FC}$ is not summed over in the expressions for subsystem $H_{\boldsymbol{V}}^{\unicode[STIX]{x1D6FC}}$.

Upon inspection, it is evident that the $\tilde{M}_{d_{0}}$ constraints of (6.30) are obeyed in each subsystem when they are exactly solved. (As in the unreduced case, the above subsystems are readily exactly solved. In particular, note that $\dot{V}_{i}^{\unicode[STIX]{x1D6FC}}=0$ in $H_{\boldsymbol{V}}^{\unicode[STIX]{x1D6FC}}$.) Consequently, the exact conservation law of the reduced system is systematically derived by simply expressing the unreduced momentum map of (6.17) in $\tilde{M}_{d_{0}}$ coordinates

(6.39)$$\begin{eqnarray}\displaystyle (\bar{\unicode[STIX]{x1D707}}_{d})_{n} & = & \displaystyle -\unicode[STIX]{x1D735}_{d}^{-}\boldsymbol{\cdot }\mathop{\sum }_{i=1}^{L}\int _{-\infty }^{\boldsymbol{X}_{i}}\text{d}\boldsymbol{X}_{i}^{\prime }\,W_{\unicode[STIX]{x1D70E}_{1}}(\boldsymbol{X}_{i}^{\prime }-\boldsymbol{x}_{n})-\unicode[STIX]{x1D735}_{d}^{-}\boldsymbol{\cdot }\boldsymbol{E}_{n}\nonumber\\ \displaystyle & := & \displaystyle \unicode[STIX]{x1D70C}_{n}-\unicode[STIX]{x1D735}_{d}^{-}\boldsymbol{\cdot }\boldsymbol{E}_{n}\nonumber\\ \displaystyle & = & \displaystyle 0,\end{eqnarray}$$

where in the final line we have noted that $\bar{\unicode[STIX]{x1D707}}_{d}$ vanishes by our previous choice of reduction to the preimage submanifold ${\unicode[STIX]{x1D707}_{d}}^{-1}(0)$. Equation (6.39) is Gauss’s law, for which we are by construction guaranteed

(6.40)$$\begin{eqnarray}\dot{\bar{\unicode[STIX]{x1D707}}}_{d}=0,\end{eqnarray}$$

as desired. (An analogous observation was made for the reduced bracket in Kraus et al. (Reference Kraus, Kormann, Morrison and Sonnendrücker2017), where the momentum map was treated as a Casimir.) The local charge conservation law of (6.20) – whose expression is unmodified in the reduced submanifold – is furthermore satisfied, since $\unicode[STIX]{x1D735}_{d}^{-}\boldsymbol{\cdot }\dot{\boldsymbol{E}}_{n}=-\unicode[STIX]{x1D735}_{d}^{-}\boldsymbol{\cdot }\boldsymbol{J}_{n}$ holds in each subsystem of $H_{d}^{\text{red}}$.