2.1. Slow manifold reduction

We will apply a dimension-reduction technique for dynamical systems with a time scale separation known as slow manifold reduction. An overview of this tool, tailored for an audience of plasma physicists, can be found in Burby (Reference Burby2020a), while a more mathematical review is given by MacKay (Reference MacKay2004). Readers can also consult Van Kampen (Reference Van Kampen1985) for additional applications and context. The reduction to a slow manifold has been used historically by Lorenz (Reference Lorenz1986) in the hopes of understanding how fast internal gravity dynamics can be decoupled from the slow Rossby waves consistent with meteorological observations to explain why the latter dominate. This was immediately followed by proving the non-existence of the aforementioned manifold in Lorenz & Krishnamurthy (Reference Lorenz and Krishnamurthy1987) and the clarification in Lorenz (Reference Lorenz1992). We note that when slow manifold existence and stability cannot be established rigorously for a given mathematical model, one may consider a ‘backwards’ theory alternative, e.g. Roberts (Reference Roberts2015), which states that there exists a system close to the original one with an exact slow manifold.

Let the system be describable by a set of ‘slow’ variables $x$ and ‘fast’ variables $y$, and suppose that the dynamics of the pair $(x,y)$ is prescribed by the system of (possibly infinite-dimensional) ordinary differential equations

(2.1)\begin{equation} \left.\begin{gathered}\dot x = g_\epsilon(x,y)\\ \epsilon\dot y = f_\epsilon(x,y), \end{gathered}\right\} \end{equation}

where $0 < \epsilon \ll 1$ is a small parameter. When $\epsilon = 0$, we see that the second differential equation degenerates to an algebraic equation that imposes the constraint $y = y_{0}^\star (x)$, where $y_{0}^\star$ is defined implicitly by

(2.2)\begin{equation} f_{0}(x, y_{0}^\star(x)) = 0. \end{equation}

It was shown by Fenichel (Reference Fenichel1979) in finite dimensions that if the manifold

(2.3)\begin{equation} \mathcal{M} = \{ (x,y): y = y_{0}^\star \} \end{equation}

is normally hyperbolic (all eigenvalues of $D_y f_\epsilon (x,y^\star (x))$ have non-zero real parts) then, under some regularity assumptions, the dynamics of $(x,y)$ with $1\gg \epsilon >0$ contains an invariant manifold that is close to $\mathcal {M}$ and that converges to $\mathcal {M}$ as $\epsilon \rightarrow 0$.

When $\mathcal {M}$ is not normally hyperbolic, the invariant manifold need not exist. However, almost invariant manifolds can still be constructed by introducing the formal power series

(2.4)\begin{equation} y^\star(x) = y^\star_{0} + \epsilon y^\star_{1} + \epsilon^2 y^\star_{2} + \dots, \end{equation}

and demanding that the ‘graph’ of $y^\star$ is invariant order by order in $\epsilon$. For non-negative integer $n$, the manifold

(2.5)\begin{equation} \mathcal{M}_\epsilon^{(n)} = \{(x,y)\, |\, y = (y_0^\star{+} \dots + \epsilon^n y_n^*)(x)\} \end{equation}

is almost invariant in the sense that the normal component of the vector field $(\dot {x},\dot {y})$ is small along $\mathcal {M}_\epsilon ^{(n)}$. Because the projection of $(\dot {x},\dot {y})$ along $\mathcal {M}_\epsilon ^{(n)}$ is $O(1)$, $\mathcal {M}_\epsilon ^{(n)}$ is known as a slow manifold of order $n$, while the formal power series $y^\star (x)$ is known as the formal slow manifold or the slaving function.

The various terms of the formal slow manifold are obtained from the appropriate ordering of the original equations of motion (2.1). One introduces Fréchet derivatives (for background see e.g. Courant & Hilbert Reference Courant and Hilbert1953; Lang Reference Lang1995), which are commonly computed using the directional derivative formula

(2.6)\begin{equation} DF(\psi)\left[\delta\psi \right] := \left.\frac{\textrm{d} }{\textrm{d} \epsilon}\right|_0 F(\psi + \epsilon\,\delta \psi), \end{equation}

where F is an arbitrary functional. The algorithm for obtaining these terms amounts to solving the partial differential equation (functional partial differential equation in infinite dimensions)

(2.7)\begin{equation} \epsilon D_x y^\star_\epsilon(x) \left[g_\epsilon(x,y^\star_\epsilon(x)) \right] = f_\epsilon(x, y_\epsilon^\star(x)), \end{equation}

order by order in $\epsilon$ and subscript x indicates Fréchet derivative with respect to x. For instance, the first-order terms are

(2.8)\begin{equation} D_x y^\star_{0}(x) \left[g_{0}(x,y^\star_{0}(x)) \right] = \,f_{1}(x, y_{0}^\star(x)) + D_y f_{0}(x,y^\star_{0}(x)) \left[y_{1}^\star(x) \right]. \end{equation}

Likewise, the second-order terms, albeit somewhat cumbersome, are written

(2.9)\begin{align} \varepsilon^{2} &\Rightarrow D y_{1}^{*}(x)\left[g_{0}\left(x, y_{0}^{*}(x)\right)\right]+D y_{0}^{*}(x)\left[g_{1}\left(x, y_{0}^*\right)(x)\right]\nonumber\\ &\quad +D y_{0}^{*}(x)\left[D_y g_{0}\left(x, y_{0}^{*}(x)\right)\left[y_{1}^{*}(x)\right]\right]=\,f_{2}\left(x, y_{0}^{*}(x)\right)\nonumber\\ &\quad +D_{y} f_{1}\left(x, y_{0}^{*}(x)\right)\left[y_{1}^{*}(x)\right]+\tfrac{1}{2} D_{y}^{2} f_{0}\left(x, y_{0}^{*}(x)\right)\left[y_{1}^{*}(x), y_{1}^{*}(x)\right]\nonumber\\ &\quad +D_{y} f_{0}\left(x, y_{0}^{*}(x)\left[y_{2}^{*}(x)\right]\right.. \end{align}

It is natural to wonder about the normal stability of the slow manifolds $\mathcal {M}_\epsilon ^{(n)}$. In other words, if a trajectory $(x(t),y(t))$ begins near $\mathcal {M}_\epsilon ^{(n)}$ then how long does it remain near $\mathcal {M}_\epsilon ^{(n)}$? When $\mathcal {M}$ is attractive, one expects stability over arbitrarily large time intervals. When $\mathcal {M}$ is normally elliptic (imaginary eigenvalues), as will be the case for the slow manifold studied in this article, the dynamics near a slow manifold is neutrally stable to leading order in perturbation theory, and therefore exhibits normal stability on at least an $O(1)$ time scale (note that the shortest time scale contained in the dynamical system (2.1) is $O(\epsilon )$). However, higher-order effects can lead to resonances that destabilize the slow manifold on $O(1/\epsilon )$ time scales. For a large class of Hamiltonian fast–slow systems with finitely many slow variables and elliptic slow manifolds, Kristiansen & Wulff (Reference Kristiansen and Wulff2016) showed that there is an optimal truncation of the slow manifold that remains normally stable on the $O(1/\epsilon )$ time scale. Establishing normal stability in the elliptic case over even larger time intervals may sometimes be achieved by identifying an adiabatic invariant whose critical manifold coincides with the formal slow manifold, as in Cotter & Reich (Reference Cotter and Reich2004). Long-term slow manifold stability is also consistent with the long-term equilibrium stability results obtained for the Vlasov–Maxwell system with non-relativistic scaling in Han-Kwan et al. (Reference Han-Kwan, Nguyen and Rousset2018). In the spirit of these previous studies, we conjecture that the slow manifold identified in this article exhibits normal stability on an $O(1/\epsilon )$ time scale, at least when optimally truncated. In § 2.4 we outline how the optimal truncation strategy developed in Kristiansen & Wulff (Reference Kristiansen and Wulff2016) applies formally to the Vlasov–Maxwell system, and therefore provide additional supporting evidence for our conjecture. In future work, we plan to investigate the possibility of our slow manifold arising as the zero level set of a wave action adiabatic invariant, which would provide further insight into normal stability.

2.2. The Vlasov–Maxwell equations as a fast–slow system

Various plasma models, including magnetohydrodynamics (MHD) and Hall MHD, rely on the assumption that one can ignore certain large-scale high-frequency modes, such as light waves. This drastically reduces the multi-scale difficulty of modelling a plasma. However, the approximation is often made in an ad hoc fashion. Here, we work with the collisionless Vlasov–Maxwell system of equations and systematically perform an asymptotic expansion such that formally the small parameter corresponds to the inverse normalized speed of light $c^{-1}$. In particular, in order to form a dimensionless small parameter one can consider a ratio between a typical speed of a particle and the speed of light. This expansion procedure bears a resemblance to the well-known Chapman–Enskog method, but retains the connection with the (infinite-dimensional) phase-space geometry underlying the closure. This is achieved by computing the slaving function $y^\star (x)$ instead of power-series expansions of solutions of the Vlasov–Maxwell equations. This more geometric perspective on the closure problem will be instrumental in uncovering the closure's Hamiltonian structure.

We assume that each plasma species $\sigma$ (typically two species are considered, but any number of species may be included without difficulty) obeys a collisionless relativistic Vlasov equation

(2.10)\begin{equation} \frac{\partial f_\sigma}{\partial t} + \boldsymbol{v}\boldsymbol{\cdot} \frac{\partial f_\sigma}{\partial \boldsymbol{r}} + e_\sigma\left( \boldsymbol{E} + \boldsymbol{v}\times \boldsymbol{B} \right) \boldsymbol{\cdot} \frac{\partial f_\sigma}{\partial \boldsymbol{p}} = 0, \end{equation}

coupled, as usual, with a self-consistent electromagnetic field defined on the triply periodic spatial domain $\mathbb {Q} = \mathbb {S}^1\times \mathbb {S}^1\times \mathbb {S}^1$. (Our slow manifold analysis is straightforward to adapt to other domain topologies and boundary conditions.) Here, as usual, B stands for the magnetic field, E electric and $\boldsymbol{v}$ and $\boldsymbol{p}$ are the velocities and momenta of the corresponding species. However, our analysis of the Hamiltonian structure on the slow manifold requires either an unbounded or periodic domain. The reason for this limitation is that the Hamiltonian structure underlying the Vlasov–Maxwell system is only understood in unbounded or periodic domains; identifying the Hamiltonian structure in a bounded domain with even simple boundary conditions like specular reflection remains an open research problem.) The field evolves according to Maxwell's equations in SI units

(2.11)\begin{gather} -\frac{1}{c^{2}} \frac{\partial \boldsymbol{E}}{\partial t}+\boldsymbol{\nabla} \times \boldsymbol{B}=\mu_{0} \boldsymbol{J} - \mu_0\boldsymbol{J}_H, \end{gather}

(2.12)\begin{gather}\frac{\partial \boldsymbol{B}}{\partial t}+\boldsymbol{\nabla} \times \boldsymbol{E}=0, \end{gather}

(2.13a,b)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{E}=\frac{\rho - \rho_H}{\varepsilon_{0}}\quad\text{and}\quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{B} = 0, \end{gather}

where μ₀ magnetic permeability of vacuum, $\mu _{0}\varepsilon _0=1/{c^{2}}$ with $c$ a dimensionful speed of light, $\rho _H = \int _{\mathbb {Q}} \rho \,\textrm {d}^3r / \int _{\mathbb {Q}}\,\textrm {d}^3r$, ρ charge density and $\boldsymbol {J}_H = \int _{\mathbb {Q}}\,\boldsymbol {J}\,\textrm {d}^3r/\int _{\mathbb {Q}}\,\textrm {d}^3r$ denotes the harmonic component of the current density (for more details see Appendix B). Note that the integral of the usual Gauss’ law $\epsilon _0\,\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {E} = \rho$ over the periodic domain $Q = (\mathbb {S}^1)^3$ imposes the restriction $\int _Q\rho \,\textrm {d}^3x = 0$, i.e. that the total charge vanishes. Therefore the $\rho _H$ in (2.13a,b) may be interpreted as the constant neutralizing background necessary to enforce this constraint when the plasma itself has a net charge. In order to find the small parameter, we pass to new dimensionless variables via a transformation that identifies typical length $L_0$ and time $t_0$ scales and introducing species mass m _σ

(2.14)\begin{equation} \left.\begin{gathered} t_0 = \sqrt{\frac{4{\rm \pi}\,\varepsilon_0m_0}{e_0\,\rho_0}}\\ m_\sigma = m_0 m_\sigma^\prime,\quad e_\sigma = e_0\,e_\sigma^\prime,\quad t = t_0 t^\prime,\\ f_\sigma(\boldsymbol{r},\boldsymbol{p}) = \frac{\rho_0}{e_0}\left(\frac{t_0}{m_0\,L_0}\right)^3\,f_\sigma^\prime\left(\frac{\boldsymbol{r}}{L_0},\frac{t_0\boldsymbol{p}}{m_0\,L_0}\right) ,\\ \boldsymbol{E}(\boldsymbol{r}) = \frac{L_0\,\rho_0}{4{\rm \pi}\varepsilon_0}\boldsymbol{E}^\prime\left( \frac{\boldsymbol{r}}{L_0} \right),\quad \boldsymbol{B}(\boldsymbol{r}) = \frac{L_0\,\rho_0}{4{\rm \pi}\varepsilon_0\,c}\,\boldsymbol{B}^\prime\left(\frac{\boldsymbol{r}}{L_0}\right). \end{gathered}\right\}\end{equation}

After the transformation is applied, primes are suppressed and a small parameter is identified as $\epsilon := L_0/(c t_0)$. In what follows, we will use the symbol $1/c$ as a convenient placeholder for $\epsilon$. The reader may observe that in this way the equations appear to be written in CGS units, which are better suited for the particular non-relativistic limit we study. After application of these steps the Vlasov equation (2.10) is re-written

(2.15)\begin{equation} \dot{f}_\sigma + \boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla} f_\sigma + e_\sigma\left( \boldsymbol{E} + \frac{1}{c}\boldsymbol{v}\times \boldsymbol{B} \right) \boldsymbol{\cdot} \frac{\partial f_\sigma}{\partial \boldsymbol{p}} = 0, \end{equation}

where the electric field $\boldsymbol {E} \in \varOmega ^{1}_T\oplus \varOmega ^{1}_L$ is a vector field on $\mathbb {Q}$ and can be decomposed into transverse and longitudinal components $\boldsymbol {E} = \boldsymbol {E}_T + \boldsymbol {E}_L =: \varPi _L \boldsymbol {E} + \varPi _T \boldsymbol {E}$, such that $\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {E}_T = 0$ and $\boldsymbol {\nabla }\times \boldsymbol {E}_L = 0$. See Appendix B for a review of Helmholtz–Hodge decomposition on compact manifolds without boundaries. (The most general decomposition of a smooth vector field on the triply periodic domain $\mathbb {Q}$ must allow also for a harmonic component $\boldsymbol {E}_H$ such that $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {E}_H = 0$ and $\boldsymbol {\nabla }\times \boldsymbol {E}_H = 0$. However, this harmonic component may be dropped self-consistently provided the harmonic component of the current density is subtracted from the Ampére–Maxwell equation, as in (2.11).) Note that the projection operators may be expressed as

(2.16)\begin{gather} \varPi_T = 1 - \varPi_L - \varPi_H ={-}\varDelta^{{-}1}\boldsymbol{\nabla}\times\boldsymbol{\nabla}\times \end{gather}

(2.17)\begin{gather}\varPi_L = \boldsymbol{\nabla}\varDelta^{{-}1}\boldsymbol{\nabla}\boldsymbol{\cdot} . \end{gather}

For details see (B17) and (B18).

Single-particle phase-space distribution functions $f_\sigma \in L^1(T^\star \mathbb {Q})\cap C^\infty (T^*\mathbb {Q})$ are integrable and smooth on $T^\star \mathbb {Q}= \mathbb {Q}\times \mathbb {R}^3$. Moreover, $\boldsymbol {p} = m_\sigma \gamma _\sigma \boldsymbol {v}$, where

(2.18)\begin{equation} \gamma_\sigma = \sqrt{1+\frac{p^2}{m^2_\sigma c^2}}. \end{equation}

The magnetic field $\boldsymbol {B}\in \varOmega _T^2$ is purely transverse $\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {B} = 0$. Likewise, in the new units we have (2.13a,b) a Gauss law constraint

(2.19)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{E}_L=4{\rm \pi} \sum_\sigma e_\sigma\int \,\textrm{d}^3 p\, f_\sigma - 4{\rm \pi}\rho_H. \end{equation}

For the transverse and longitudinal components of the electric field we have the evolution equations

(2.20)\begin{gather} \frac{1}{c} \frac{\partial \boldsymbol{E}_L}{\partial t} ={-}\frac{4{\rm \pi}}{c} \boldsymbol{J}_L, \end{gather}

(2.21)\begin{gather}\frac{1}{c} \frac{\partial \boldsymbol{E}_T}{\partial t} = \boldsymbol{\nabla}\times\boldsymbol{B}-\frac{4{\rm \pi}}{c} \boldsymbol{J}_T, \end{gather}

where

(2.22)\begin{equation} \boldsymbol{J} := \sum_\sigma e_\sigma\int \,\textrm{d}^3 p \, \boldsymbol{v}\, f_\sigma. \end{equation}

Using the gauge where the electrostatic potential vanishes, the relationship between $\boldsymbol {E}_L$ and the longitudinal part of the vector potential $\boldsymbol {A}_L$ is

(2.23)\begin{equation} \boldsymbol{E}_{L}={-}\frac{1}{c} \frac{\partial \boldsymbol{A}_{L}}{\partial t}. \end{equation}

Finally, the system is closed by Faraday's law

(2.24)\begin{equation} \frac{1}{c} \frac{\partial \boldsymbol{B}_T}{\partial t} ={-} \boldsymbol{\nabla}\times \boldsymbol{E}_T. \end{equation}

It has been proved by Degond & Raviart (Reference Degond and Raviart1992) that the Darwin model in three-dimensional cases approximates Maxwell's equations for appropriate initial data up to the second order of the dimensionless parameter $v/c$ for magnetic field $B$ and to the third order for electric field $E$, where $v$ is the characteristic velocity. (The somewhat mysterious notion of “appropriate” initial data is laid bare in the context of slow manifold reduction theory; well-prepared initial data are merely initial conditions chosen to lie on the slow manifold.) In what follows we will re-derive this result and go an order beyond Darwin to establish the next approximation.

2.3. The dark slow manifold: piezoelectric correction

The candidates for the slow variables in our case are $X = (L^{1}(T^\star \mathbb {Q})\cap C^\infty (T^*\mathbb {Q}))\times \varOmega _L^1 \ni (\,f, \boldsymbol {E}_L)$ and, following the recipe (2.4), we treat transverse fields as fast $Y = \varOmega _T^1\times \varOmega _L^2 \ni (\boldsymbol {E}_T,\boldsymbol {B})$ (which physically carry light waves). Our slaving functions are therefore

(2.25)\begin{equation} \left.\begin{gathered} \boldsymbol{E}^\star_T(\,f_\sigma,\boldsymbol{E}_L) = \boldsymbol{E}^\star_{T0}(\,f_\sigma,\boldsymbol{E}_L) + \frac{1}{c}\boldsymbol{E}^\star_{T1}(\,f_\sigma,\boldsymbol{E}_L)+\dots\\ \boldsymbol{B}^\star(\,f_\sigma,\boldsymbol{E}_L) = \boldsymbol{B}^\star_{0}(\,f_\sigma,\boldsymbol{E}_L) + \frac{1}{c}\boldsymbol{B}^\star_{1}(\,f_\sigma,\boldsymbol{E}_L)+\dots \end{gathered}\right\} .\end{equation}

Using (2.25) we can cast (2.21) and (2.24) as

(2.26)\begin{equation} \boldsymbol{\nabla}\times\boldsymbol{B}^\star{=} \frac{4{\rm \pi}}{c}\, \boldsymbol{J}_T +\frac{1}{c} D_f \boldsymbol{E}^{{\star}}_{T}\left[\dot{f} \right] + \frac{1}{c} D_{{\boldsymbol{E}}_{L}} \boldsymbol{E}^{{\star}}_{T}\left[\dot{\boldsymbol{E}}_L \right], \end{equation}

and

(2.27)\begin{equation} -\boldsymbol{\nabla}\times\boldsymbol{E}^\star_T =\frac{1}{c} D_f \boldsymbol{B}^\star\left[\dot{f} \right] + \frac{1}{c} D_{\boldsymbol{E}_{L}} \boldsymbol{B}^{{\star}}\left[\dot{\boldsymbol{E}}_{L}\right]. \end{equation}

We find trivially that in the zeroth-order expansion with respect to $1/c$

(2.28a,b)\begin{equation} \boldsymbol{\nabla} \times \boldsymbol{E}_{T0}^{{\star}}=0,\quad \boldsymbol{\nabla} \times \boldsymbol{B}_{0}^{{\star}}=0, \end{equation}

which implies

(2.29a,b)\begin{equation} \boldsymbol{E}^\star_{T0}= 0,\quad \boldsymbol{B}^\star_{0} =0, \end{equation}

because $\boldsymbol {E}_T$ and $\boldsymbol {B}$ are solenoidal. Notice that the first three non-zero terms in the expansion of $\boldsymbol {J}$ are

(2.30)\begin{equation} \left.\begin{gathered} \boldsymbol{J}_{0}=\sum_\sigma\frac{e_\sigma}{m_\sigma} \int \,\textrm{d}^{3} p\, f_\sigma \boldsymbol{p}, \quad \boldsymbol{J}_{2}={-}\sum_\sigma\frac{e_\sigma}{2 m^{3}_\sigma} \int \,\textrm{d}^{3} p\, f_\sigma p^{2} \boldsymbol{p}, \\ \boldsymbol{J}_{4}=\sum_\sigma\frac{3 \,e_\sigma}{8 m^{5}_\sigma} \int \,\textrm{d}^{3} p\, f_\sigma p^{4} \boldsymbol{p},\quad \dots \end{gathered}\right\}.\end{equation}

The absence of transverse fields to this order suggests absence of light waves. Below, we will show that, if we prepare plasma in this special state, light waves are going to be suppressed in future. In the first order

(2.31a,b)\begin{equation} \boldsymbol{\nabla} \times \boldsymbol{E}_{T1}^{*}=0,\quad \frac{1}{c} \boldsymbol{\nabla} \times \boldsymbol{B}_{1}^{*}={-}\frac{4 {\rm \pi}}{c} \boldsymbol{\nabla} \times \varDelta^{{-}1} \boldsymbol{\nabla}\times \boldsymbol{J}_{0}. \end{equation}

Thus the following (Darwin) magnetic field is recovered

(2.32)\begin{equation} \boldsymbol{B}^\star_{1} ={-} 4{\rm \pi} \,\varDelta^{{-}1} \boldsymbol{\nabla}\times \boldsymbol{J}_{0}:={-} \varDelta^{{-}1}\boldsymbol{\nabla}\times\sum_\sigma \frac{4{\rm \pi} e_\sigma}{m_\sigma}\int \,\textrm{d}^3 p \, f_\sigma \boldsymbol{p}. \end{equation}

In many circumstances we prefer to work with the transverse vector potential instead of the magnetic field

(2.33)\begin{equation} \boldsymbol{A}^\star_{T1} ={-}4{\rm \pi}\, \varDelta^{{-}1} \varPi_T\, \boldsymbol{J}_{0} , \end{equation}

where special care is taken to ensure that the Green's function acts on an object that is already in the transverse space. Meanwhile, it turns out that the electric field is absent at this order

(2.34)\begin{equation} \boldsymbol{E}^\star_{1} = 0. \end{equation}

This limit constitutes the well-known Darwin approximation. As stated earlier, it is typically applied in quasistatic situations. However, Eremin et al. (Reference Eremin, Hemke, Brinkmann and Mussenbrock2013) show its utility in the intermediate range for high-frequency capacitively coupled discharges. The Darwin approximation offers a simple and efficient way of carrying out electromagnetic simulations as it removes the Courant condition plaguing explicit electromagnetic algorithms and can be implemented as a straightforward modification of electrostatic algorithms.

From (2.32) and the $1/c^2$ expansion of (2.21) it is evident that, this time

(2.35)\begin{equation} \boldsymbol{B}^\star_{2} = 0. \end{equation}

On the other hand, from $1/c^2$ expansion of (2.27) we see

(2.36)\begin{equation} \boldsymbol{\nabla}\times \boldsymbol{E}^\star_{T2} ={-} D_f \boldsymbol{B}^\star_{1} \left[\dot{f}^{0} \right] \equiv{-}\sum_\sigma D_{f_\sigma} \boldsymbol{B}^\star_{1} \left[\dot{f}^{0}_\sigma \right], \end{equation}

where $D$ is a Fréchet derivative and we only need to keep zeroth-order terms (with respect to $1/c$ expansion) in the expression for $\dot {f}$ (see (2.15)). Because of (2.18) it is clear that the effects of relativity are not felt in this limit and thus we have $\boldsymbol {v} \approx \boldsymbol {p}/m$. Thus, to this order,

(2.37)\begin{equation} \dot{f}_\sigma^{0}+\frac{\boldsymbol{p}}{m_\sigma} \boldsymbol{\cdot} \boldsymbol{\nabla} f_\sigma+e_\sigma \boldsymbol{E}_{L} \boldsymbol{\cdot} \frac{\partial f_\sigma}{\partial \boldsymbol{p}}=0. \end{equation}

From (2.36) one obtains, after integration by parts,

(2.38)\begin{equation} \boldsymbol{E}^\star_{T2} = {-} \sum_\sigma\varDelta^{{-}1} \varPi_T\left(4{\rm \pi} e_\sigma\, \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{T}_{0\sigma} - \omega_{p\sigma}^2 \boldsymbol{E}_L \right), \end{equation}

where we have introduced the zeroth-order stress tensor

(2.39)\begin{equation} \boldsymbol{T}_{0\sigma} := \int \,\textrm{d}^3 p\, \frac{{\boldsymbol{p} \, \boldsymbol{p}}}{m^2_\sigma} f_\sigma,\end{equation}

and plasma frequency

(2.40a,b)\begin{equation} \omega_{p\sigma}^2 := \frac{4{\rm \pi} e^2_\sigma}{m_\sigma}\int \,\textrm{d}^3 p \,f_\sigma.\quad {\omega_p^2 := \sum_\sigma \omega_{p\sigma}}. \end{equation}

Note that, due to the spatial dependence of $\omega _{p\sigma }^2$, the last term in (2.38) does not vanish despite the fact that it is acted upon by the transverse projection operator.

According to (2.38), stress generates transverse electric fields. Obviously, in a quasineutral plasma if the species are moving with the same velocity and the distribution function, the term would vanish. Thus we conclude that it has a two-fluid nature. The origin of this term can be, in fact, traced back to the Ohm's law found in the Braginskii multifluid model with quasineutrality. Similar reformulation leading to this term has been proposed for Maxwell–Euler system by Degond, Deluzet & Doyen (Reference Degond, Deluzet and Doyen2017). A similar effect may be found in certain solids. In solids piezoelectric effect is the induction of an electric charge in response to an applied mechanical strain, see Zhang & Hoshino (Reference Zhang and Hoshino2014), for example. It has numerous applications, most notably in microphones.

We can show that the stress-tensor term can be simplified for a non-uniform drifting Maxwellian (since only zeroth-order corrections are required we can use the non-relativistic expression)

(2.41)\begin{equation} f(\boldsymbol{r},\boldsymbol{v}) = n(\boldsymbol{r}) \left( \frac{m}{{\rm \pi} T} \right)^{{3}/{2}} \exp\left(-\frac{m(\boldsymbol{v} - \boldsymbol{V}(\boldsymbol{r}))^2}{2 T}\right). \end{equation}

We have

(2.42)\begin{equation} \int \,\textrm{d}^3 v\, v_j v_k f = \delta_{jk} \frac{n T}{2 m} + n V_j V_k. \end{equation}

Next, because of the longitudinal form of the first term it vanishes under the transverse projection. The stress term survives, in general. Of course, the second term in (2.38) vanishes in case of a uniform density.

The second term in (2.39) demonstrates the coupling between the longitudinal mode, originating in Langmuir oscillations, and the transverse mode. This effect is akin to the generation of transverse waves from the longitudinal ones. However, in this approximation the field that is generated is of the static nature, rather than radiation.

To see what happens further we go to the next order. From the application of (2.35) and $\dot {f}_{1} = 0$ to (2.26) it is easy to see that

(2.43)\begin{equation} \boldsymbol{E}_{T3}^\star{=} 0. \end{equation}

On the other hand, we have a rather more complicated correction to the magnetic field determined by

(2.44)\begin{equation} \boldsymbol{\nabla}\times\boldsymbol{B}_{3}^\star{=} 4{\rm \pi} \boldsymbol{J}_{T2}\left[f \right] + D_f \boldsymbol{E}^\star_{T2}\left[\dot{f}^{0} \right] + D_{\boldsymbol{E}_L} \boldsymbol{E}^\star_{T2}\left[\dot{\boldsymbol{E}}_L^{0} \right], \end{equation}

which, using (2.20) and (2.30), leads to

(2.45)\begin{align} \boldsymbol{\nabla} \times \boldsymbol{B}_{3}^{*}&={-}\sum_\sigma\frac{4 {\rm \pi}e_\sigma}{2 m^{3}_\sigma} \varPi_{T} \int \,\textrm{d}^{3}p\, f_\sigma p^{2} \boldsymbol{p}{-}\sum_\sigma\frac{4 {\rm \pi}e_\sigma}{m_\sigma} \varDelta^{{-}1}{\varPi}_{T} \nonumber\\ &\quad \cdot\int \,\textrm{d}^{3} {p} \left(4{\rm \pi} e_\sigma f_\sigma {\varPi}_{L} \boldsymbol{J}_{0} +e_\sigma \boldsymbol{E}_{L}\left[\frac{\boldsymbol{p}}{m_\sigma} \boldsymbol{\cdot} \boldsymbol{{\nabla}} f_\sigma+\boldsymbol{E}_L \boldsymbol{\cdot}\frac{\partial f_\sigma}{\partial \boldsymbol{p}}\right]\right.\nonumber\\ &\quad \left. -\boldsymbol{\nabla} \boldsymbol{\cdot}\left[\frac{\boldsymbol{p}{\boldsymbol{p}}}{m_\sigma} \left(\frac{\boldsymbol{p}}{m_\sigma} \boldsymbol{\cdot} {\boldsymbol{\nabla}} f_\sigma +e_\sigma \boldsymbol{E}_{L} \boldsymbol{\cdot} \frac{\partial f_\sigma}{\partial \boldsymbol{p}}\right)\right]\right). \end{align}

Thus we collect these terms and use integration by parts to obtain

(2.46)\begin{align} \boldsymbol{A}_{T3}^\star &={-}\varDelta^{{-}2}\varPi_T\left[ \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{Q}_0 - \boldsymbol{\nabla}\boldsymbol{\cdot} \left(\boldsymbol{E}_L\, \boldsymbol{W} + \boldsymbol{W} \,\boldsymbol{E}_L \right) -\boldsymbol{E}_L \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{W} \right]\nonumber\\ &\quad - 4{\rm \pi}\varDelta^{{-}1} \varPi_T\boldsymbol{J}_{2} {+}\sum_\sigma 4{\rm \pi} \varDelta^{{-}2}\varPi_T\omega_{p\sigma}^2\varPi_L \boldsymbol{J}_{0} , \end{align}

where we have defined the charge-weighted heat flux tensor

(2.47)\begin{equation} \boldsymbol{Q}_0 := \sum_\sigma 4{\rm \pi} e_\sigma \int \,\textrm{d}^3p \, \frac{{\boldsymbol{p} \boldsymbol{p} \boldsymbol{p}}}{m^3_\sigma} f_\sigma, \end{equation}

and introduced the velocity weighted by the plasma frequency

(2.48)\begin{equation} \boldsymbol{W} := \sum_\sigma \frac{4{\rm \pi} e_\sigma^2 }{m_\sigma}\int \frac{\boldsymbol{p}}{m_\sigma}f_\sigma\,\textrm{d}^3p. \end{equation}

It is somewhat remarkable that a charge-density-weighted heat flux produces a magnetic field in a dark plasma. To the authors’ knowledge, this has not been reported in the literature yet.

2.4. Stability of the dark slow manifold

Our analysis of the slow manifold so far has been entirely based on asymptotic expansions. Therefore, it is unclear how long an initially dark Vlasov–Maxwell solution will remain dark. Any thorough analysis of this question will require a detailed functional analysis. However, we will now show that arguments from Kristiansen & Wulff (Reference Kristiansen and Wulff2016) may be applied formally to the Vlasov–Maxwell system in order to obtain a heuristic that suggests the time scale is at least $O(c)$ as $c\rightarrow \infty$ for initial conditions prepared ‘optimally close’ to the slow manifold. Thus, we believe there is compelling evidence to conjecture that darkness persists over large time intervals in non-relativistic plasmas. Kristiansen & Wulff (Reference Kristiansen and Wulff2016) study slow manifolds that arise in a certain class of Hamiltonian fast–slow systems. The system phase space is assumed to be of the form $\mathcal {P} = \mathcal {X}\times \mathcal {Y}\ni (x,y)$, where $x$ is the slow variable, $y$ is the fast variable and both $\mathcal {X}$ and $\mathcal {Y}$ are real Hilbert spaces. The fast space $\mathcal {Y}$ is allowed to be infinite-dimensional, but the dimension of the slow space $\mathcal {X}$ is assumed to be finite. While the rigorous arguments in this work use finite-dimensionality of $\mathcal {X}$ in an essential way, at a formal level the arguments apply to infinite-dimensional $\mathcal {X}$ as well. The Poisson tensor on $\mathcal {P}$ is given by $\mathcal {J} = \epsilon \,\mathcal {J}_{\mathcal {X}} + \mathcal {J}_{\mathcal {Y}}$, where $\epsilon \ll 1$ and $\mathcal {J}_{\mathcal {X}},\mathcal {J}_{\mathcal {Y}}$ are Poisson tensors on $\mathcal {X}$ and $\mathcal {Y}$, respectively. The system Hamiltonian is required to be of the form

(2.49)\begin{equation} H_\epsilon(x,y) = h_\epsilon(x) + \langle R_\epsilon(x),y \rangle+\tfrac{1}{2}\langle y,(L+a_\epsilon(x))y \rangle + \gamma_\epsilon(x,y),\end{equation}

where all $\epsilon$ dependence is continuous, $L,a_\epsilon (x)$ are self-adjoint, $L+a_\epsilon (x)$ is invertible and $\gamma _\epsilon (x,y) = O(||y||^3)$. It is also assumed that $h_\epsilon (x),R_\epsilon (x),a_\epsilon (x)$ and $\gamma _\epsilon (x,y)$ are real analytic. In this context, $y=0$ defines a limiting slow manifold.

The first major result from Kristiansen & Wulff (Reference Kristiansen and Wulff2016), Theorem 2.1, states that there is a near-identity canonical transformation $(x,y)\mapsto (\bar {x},\bar {y})$ that transforms the system Hamiltonian (2.49) into

(2.50)\begin{equation} \bar{H}_\epsilon(\bar{x},\bar{y}) = \bar{h}_\epsilon(\bar{x}) + \langle \bar{R}_\epsilon(\bar{x}),\bar{y} \rangle+\tfrac{1}{2}\langle \bar{y},(L+\bar{a}_\epsilon({\bar{x}}))\bar{y} \rangle + \bar{\gamma}_\epsilon(\bar{x},\bar{y}), \end{equation}

where $\bar {a}_\epsilon (\bar {x})$ is self-adjoint, all barred quantities are close to their un-barred counterparts and, crucially, $\bar {R}_\epsilon = O(\exp (-\kappa /\epsilon ))$ for some positive constant $\kappa$. The evolution equation for $\bar {y}$ is therefore $\dot {\bar {y}} = \mathcal {J}_{\mathcal {Y}}\,\boldsymbol {\nabla }_{\bar {y}}\bar {H}_\epsilon = \mathcal {J}_{\mathcal {Y}}(L+\bar {a}_\epsilon (x))\bar {y} + \mathcal {J}_{\mathcal {Y}}\,\boldsymbol {\nabla }_{\bar {y}}\bar {\gamma }_\epsilon + \mathcal {J}_{\mathcal {Y}}\bar {R}_\epsilon (\bar {x})$. Since the right-hand side of this equation is exponentially small at $\bar {y}=0$, the set $\{(\bar {x},\bar {y})\mid \bar {y}=0\}$ defines an exponentially accurate slow manifold.

The barred coordinate system is constructed by applying a sequence of $N$ near-identity canonical transformations such that, after the $n{\text {th}}$ ($n\leqslant N$) step in the iteration, the transformed $R_\epsilon$ is smaller by a factor of $\epsilon /\xi _n$ relative to the previous iteration. Here, $\xi _n$ is a parameter that controls the loss of regularity introduced by the $n{\text {th}}$ transformation; smaller $\xi$ corresponds to a smaller loss of regularity, but a less dramatic reduction in the size of the transformed $R_\epsilon$. In order to produce an exponential reduction in the size of the transformed $R_\epsilon$ while avoiding a catastrophic loss of regularity, the authors choose $\xi _n \sim 2\epsilon$ for large $n$ (so $\epsilon /\xi _n \sim 1/2$) and $N\sim 1/\epsilon$. These choices ensure that size of the transformed $R_\epsilon$ is halved after a single iteration, and reduced by a factor of $2^{-1/\epsilon } = \exp (-[\ln 2]/\epsilon )$ after $N$ iterations, which is the desired exponential effect.

The next major result from Kristiansen & Wulff (Reference Kristiansen and Wulff2016), Corollary 2.2, uses the special coordinates described in the previous paragraph to estimate the stability time scale of the $\bar {y}=0$ slow manifold. Assuming $L+a_\epsilon (x)$ is positive definite, the result states that a solution of the fast–slow system with $\bar {y}(\tau =0)=0$ will satisfy $\bar {y}(\tau ) = O(\exp (-\kappa _1/\epsilon ))$ for times $\tau \in [0,\kappa _2/\epsilon ^2]$. In other words, if a solution starts on the exponentially accurate slow manifold then it will remain exponentially close to that manifold on a time scale that is at least $O(\epsilon ^{-2})$ in $\tau$. The proof is a Lyapunov-type argument that exploits the near constancy of the function

(2.51)\begin{equation} I = \tfrac{1}{2}\langle \bar{y},(L+\bar{a}_\epsilon(\bar{x}))\bar{y}\rangle+\bar{\gamma}_\epsilon(\bar{x},\bar{y}) \end{equation}

along solutions of the fast–slow system. Near invariance of $I$ follows from the following simple calculation:

(2.52)\begin{align} \dot{I} &= \left\langle \mathcal{J}_{\mathcal{Y}}(L+\bar{a}_\epsilon(\bar{x}))\bar{y} + \mathcal{J}_{\mathcal{Y}}\,\boldsymbol{\nabla}_{\bar{y}}\bar{\gamma}_\epsilon + \mathcal{J}_{\mathcal{Y}}\bar{R}_\epsilon(\bar{x}) ,(L+\bar{a}_\epsilon(\bar{x}))\bar{y}+\boldsymbol{\nabla}_{\bar{y}}\bar{\gamma}_\epsilon\right\rangle+O(\epsilon) \nonumber\\ & = \left\langle \mathcal{J}_{\mathcal{Y}}\bar{R}_\epsilon(\bar{x}) ,(L+\bar{a}_\epsilon(\bar{x}))\bar{y}+ \boldsymbol{\nabla}_{\bar{y}}\bar{\gamma}_\epsilon \right\rangle+O(\epsilon)\nonumber\\ & = O(\exp(-\kappa/\epsilon)) + O(\epsilon). \end{align}

The arguments from Kristiansen & Wulff (Reference Kristiansen and Wulff2016) apply formally to the Vlasov–Maxwell system because of the following. First, we recall that the Vlasov–Maxwell system may be written as a Hamiltonian system on the space $\mathcal {P}=\mathcal {X}\times \mathcal {Y}$, where $\mathcal {X}\ni f$ is the space of canonical momentum distribution functions $f(\boldsymbol {x},\boldsymbol {{\rm \pi} })$, and $\mathcal {Y}\ni (\boldsymbol {E}_T,\boldsymbol {A}_T)$ comprises pairs of transverse vector potentials $\boldsymbol {A}_T$ and transverse electric fields $\boldsymbol {E}_T$. The Poisson bracket of functionals $F(\,f,\boldsymbol {E}_T,\boldsymbol {A}_T),G(\,f,\boldsymbol {E}_T,\boldsymbol {A}_T)$ is given by $\{F,G\} = \epsilon \,\{F,G\}_{\mathcal {X}}+\{F,G\}_{\mathcal {Y}}$, where $\epsilon = 1/c$ and

(2.53)\begin{equation} \left.\begin{gathered} \{F,G\}_{\mathcal{X}} = \int \left[\frac{\delta F}{\delta f},\frac{\delta G}{\delta f}\right]f\, \textrm{d}^3\boldsymbol{r}\,\textrm{d}^3\boldsymbol{\rm \pi}\\ \{F,G\}_{\mathcal{Y}} = 4{\rm \pi}\int \left(\frac{\delta F}{\delta\boldsymbol{E}_T} \boldsymbol{\cdot}\frac{\delta G}{\delta\boldsymbol{A}_T} - \frac{\delta G}{\delta\boldsymbol{E}_T}\boldsymbol{\cdot}\frac{\delta F}{\delta\boldsymbol{A}_T} \right) \textrm{d}^3\boldsymbol{r}. \end{gathered}\right\} \end{equation}

Here, $[\boldsymbol {\cdot },\boldsymbol {\cdot }]$ denotes the usual canonical Poisson bracket on $(\boldsymbol {x},\boldsymbol {{\rm \pi} })$-space, as appropriate when working with canonical momenta $\boldsymbol {{\rm \pi} }$. Note that this is the only section where we work in coordinates which depend on canonical momentum for consistency with Kristiansen & Wulff (Reference Kristiansen and Wulff2016). This bracket is obtained by quotioning the ‘canonical’ Vlasov–Maxwell bracket discussed in Marsden & Weinstein (Reference Marsden and Weinstein1982) by gauge transformations. The system Hamiltonian is

(2.54)\begin{align} H_\epsilon(x,y) &= \frac{1}{8{\rm \pi}}\int \boldsymbol{E}_T\boldsymbol{\cdot}\boldsymbol{E}_T\,\textrm{d}^3\boldsymbol{r} + \frac{1}{8{\rm \pi}}\int \boldsymbol{\nabla}\times\boldsymbol{A}_T\boldsymbol{\cdot}\boldsymbol{\nabla} \times\boldsymbol{A}_T\,\textrm{d}^3\boldsymbol{r} \nonumber\\ &\quad + \frac{1}{8{\rm \pi}}\int \boldsymbol{E}_L(\,f)\boldsymbol{\cdot}\boldsymbol{E}_L(\,f)\,\textrm{d}^3\boldsymbol{r} + \int \frac{1}{2m} \left(\boldsymbol{\rm \pi}-\epsilon\,e\boldsymbol{A}_T\right)^2f\,\textrm{d}^3\boldsymbol{r}\,\textrm{d}^3\boldsymbol{\rm \pi}, \end{align}

where the longitudinal electric field $\boldsymbol {E}_L(\,f) = -\boldsymbol {\nabla }\varphi (\,f)$ is the unique solution of the elliptic partial differential equation $-\Delta \varphi (\,f) =4{\rm \pi} e\,n(\,f) - 4{\rm \pi} \rho _0$, where $\rho _0$ denotes a constant neutralizing background charge and $n(\,f) = \int f(\boldsymbol {r},\boldsymbol {{\rm \pi} })\,\textrm {d}^3\boldsymbol {{\rm \pi} }$ is the number density. Note that the time variable $\tau$ for this Hamiltonian formulation of Vlasov–Maxwell is related to the time variable $t$ used elsewhere in this article by $\tau = c \,t$. As required by Kristiansen & Wulff (Reference Kristiansen and Wulff2016), the Poisson bracket has a product structure that is compatible with the fast–slow split, and the Hamiltonian has the form (2.49), with

(2.55)\begin{gather} L(\boldsymbol{E}_T,\boldsymbol{A}_T) = \frac{1}{4{\rm \pi}}(\boldsymbol{E}_T,-\Delta \boldsymbol{A}_T) \end{gather}

(2.56)\begin{gather}a_\epsilon(\,f)(\boldsymbol{E}_T,\boldsymbol{A}_T) = \epsilon^2\left(0,\frac{e^2\,n(\,f)}{m}\,\boldsymbol{A}_T\right) \end{gather}

(2.57)\begin{gather}h_\epsilon(\,f) = \int \frac{1}{2m} |\boldsymbol{\rm \pi}|^2\,f\,\textrm{d}^3\boldsymbol{r}\,\textrm{d}^3\boldsymbol{\rm \pi}+ \frac{1}{8{\rm \pi}}\int \boldsymbol{E}_L(\,f)\boldsymbol{\cdot}\boldsymbol{E}_L(\,f)\,\textrm{d}^3\boldsymbol{r} \end{gather}

(2.58)\begin{gather}R_\epsilon(\,f) ={-}\epsilon\,e\int \frac{\boldsymbol{\rm \pi}}{m}\,f\,\textrm{d}^3\boldsymbol{\rm \pi} \end{gather}

(2.59)\begin{gather}\gamma_\epsilon(\,f,\boldsymbol{E}_T,\boldsymbol{A}_T) = 0. \end{gather}

It follows that the formal transformation $(x,y)\mapsto (\bar {x},\bar {y})$, mentioned earlier, that ‘flattens’ the slow manifold can be derived for this system much as in Kristiansen & Wulff (Reference Kristiansen and Wulff2016). Analyticity of the Vlasov–Maxwell system in $(\,f,\boldsymbol {E}_T,\boldsymbol {A}_T)$ then suggests that optimal truncation of the formal transformation should produce an exponentially accurate slow manifold, as in Theorem 2.1 from Kristiansen & Wulff (Reference Kristiansen and Wulff2016). If this is indeed the case, then positive definiteness of $L+a_\epsilon$ (easily seen to be true) together with the proof of Corollary 2.2 in Kristiansen & Wulff (Reference Kristiansen and Wulff2016) implies normal stability of the exponentially accurate dark slow manifold for $\tau \in [0,c_2/\epsilon ^2]$, or equivalently $t\in [0,c_2/\epsilon ]$. In other words, there exist solutions of the Vlasov–Maxwell system with initially exponentially small light-wave activity that remain exponentially close to the dark slow manifold on $O(1/\epsilon )=O(c)$ time intervals.

While this argument suggests a concrete path toward proving persistence of dark plasma states in the Vlasov–Maxwell phase space on $O(1/\epsilon )$ time scales, it is important to understand the argument's shortcomings. First and foremost, in Kristiansen & Wulff (Reference Kristiansen and Wulff2016), it was crucial for the authors to understand precisely the loss of regularity introduced by each step in their sequence of canonical transformations. Due to specifics of the function-analytic setting considered by Kristiansen & Wulff (Reference Kristiansen and Wulff2016), the details of this regularity loss may very well differ for the Vlasov–Maxwell system. Any modifications that might arise need to be studied and accounted for appropriately in the optimal truncation procedure. Finally, this argument establishes long-term stability of the optimally truncated slow manifold, but not of lower-order truncations. Stability of these lower-order truncations must be assessed using other methods.

On the other hand, the issues inherent to this argument disappear if the Vlasov–Maxwell system is replaced with a continuous-time structure-preserving discretization as in Krauss et al. (Reference Krauss, Kormann, Morrison and Sonnendrücker2017) or Burby (Reference Burby2017a). It would be interesting to use this argument to study non-relativistic slow manifolds in discrete Vlasov–Maxwell systems in the future.

3.1. Mathematical preliminaries

First, we will review some well-known results from finite-dimensional Lagrangian mechanics. For details see Morrison (Reference Morrison1998), Cary & Brizard (Reference Cary and Brizard2009), José & Saletan (Reference José and Saletan1998) and Marsden & Ratiu (Reference Marsden and Ratiu1999), for instance, although we follow slightly different conventions. The Lagrangian for an individual charged particle moving in an electromagnetic field may be written in the ‘phase-space’ form as

(3.1)\begin{equation} L = \left(\boldsymbol{p}+\frac{e}{c} \boldsymbol{A}\right)\boldsymbol{\cdot} \dot{\boldsymbol{r}} - \left(e \,\varPhi(\boldsymbol{r}) + \frac{m v^2}{2} \right) =: {\theta}_B \boldsymbol{\cdot} {\dot{\boldsymbol{z}}} - H({\boldsymbol{z}}), \end{equation}

where we have introduced the Lagrange $1$-form $\theta _B = \left (\boldsymbol {p}+({e}/{c})\boldsymbol {A}\right )\boldsymbol {\cdot } \textrm {d}\boldsymbol {r}$ and the phase-space coordinate $\boldsymbol {z} = (\boldsymbol {r},\boldsymbol {p})$. For a review of differential forms for plasma physicist see MacKay (Reference MacKay2020). The symplectic $2$-form on the momentum phase space associated with $L$ can be obtained via ${\omega }_B = -\textrm {d}{\theta }_B$ and has the form

(3.2)\begin{equation} {\omega_B = \textrm{d}\boldsymbol{r}\wedge \textrm{d}\boldsymbol{p} + \frac{e}{c} \,\textrm{d}r^i \wedge \textrm{d}r^j \, \frac{\partial A_i}{\partial r^j}}, \end{equation}

which leads to the symplectic version of Hamilton's equations $\dot {q} = \partial H/\partial p$ and $\dot {p}= - \partial H/\partial q$

(3.3)\begin{equation} {\omega_B({\dot{\boldsymbol{z}}}, \boldsymbol{\cdot})} \equiv i_{\dot{\boldsymbol{z}}} {\omega}_B = \textrm{d} H. \end{equation}

The symplectic form can also be inverted to obtain the Poisson tensor

(3.4)\begin{equation} \mathcal{J}_B = \begin{pmatrix} 0 & \delta^{ij} \\ -\delta^{ij} & \epsilon_{ijk} \dfrac{e}{c} B^k \end{pmatrix},\end{equation}

and thus the Poisson bracket is

(3.5)\begin{equation} \{f,g\} :=\frac{\partial f}{\partial z^\alpha} J^{\alpha\beta} \frac{\partial g}{\partial z^\beta}= \frac{\partial f}{\partial \boldsymbol{r}} \boldsymbol{\cdot}\frac{\partial g}{\partial \boldsymbol{p}} - \frac{\partial g}{\partial \boldsymbol{r}} \boldsymbol{\cdot}\frac{\partial f}{\partial \boldsymbol{p}} + \frac{e\boldsymbol{B}}{c}\boldsymbol{\cdot} \frac{\partial f}{\partial \boldsymbol{p}}\times \frac{\partial g}{\partial \boldsymbol{p}}. \end{equation}

For example, in coordinates one recovers $\dot {z}^\alpha = \{z, H\} = J^{\alpha \beta } \partial H/\partial z^\beta$. When $\boldsymbol {B}=0$, we write $\theta _B = \theta _0 =\boldsymbol {p}\boldsymbol {\cdot } \textrm {d}{\boldsymbol {r}}$ and we recover the canonical symplectic form $\omega _0=-\textrm {d}\theta _0 = \textrm {d}\boldsymbol {r}\wedge \textrm {d}\boldsymbol {p}$ and Poisson tensor $\mathcal {J}_0$.

Next, we generalize to treat a distribution of particles. Neglecting the electromagnetic field for now, the analogue of $\theta _B$ for an ensemble of charged particles, which we will denote $\varTheta$, is a $1$-form on the infinite-dimensional space of Lagrangian configuration maps

(3.6)\begin{equation} g: \overset{\circ}{\mathcal{P}} \rightarrow \mathcal{P},\end{equation}

where $\overset {\circ }{\mathcal {P}}\ni (\boldsymbol {r}_0,\boldsymbol {p}_0)$ denotes the space of particle labels and $\mathcal {P}\ni (\boldsymbol {r},\boldsymbol {p})$ denotes the Eulerian phase space. (We remark that, as sets, $\overset {\circ }{\mathcal {P}}$ and $\mathcal {P}$ are the same. The notational distinction reflects the different physical interpretations of particle locations and particle labels.) In particular, using the notation introduced in Appendix A, we have

(3.7)\begin{equation} \iota_{\dot{g}}\varTheta_g = \int \iota_V\theta_0\,f\,\textrm{d}^3r\,\textrm{d}^3p, \end{equation}

where the Eulerian phase-space velocity is given by $V = \dot {g}\circ g^{-1}$ and $f$ is defined in terms of $g$ according to

(3.8)\begin{equation} f\,\textrm{d}^3r\,\textrm{d}^3p = g_{*}(\overset{\circ}{f}\,\textrm{d}^3r_0\,\textrm{d}^3p_0).\end{equation}

Here, $\overset {\circ }{f}$ is some fixed reference distribution on label space.

To find the symplectic form associated with $\varTheta$, we must compute the exterior derivative $\varOmega = -\textrm {d}\varTheta$. We proceed by applying the formula (A2) from Appendix A with $\alpha = \varTheta$. To begin, we introduce the functional $I$ on the space of curves $g(t)$ with $g(-1)$ and $g(1)$ fixed whose value at $g(t)$ is $I(g) = \int _{-1}^{1}\int \iota _V\theta _0\,f\,\textrm {d}^3r\,\textrm {d}^3p\,\textrm {d}t$. The first variation of $I$ is given by

(3.9)\begin{equation} \delta I(g)[\delta g] = \int_{{-}1}^{1}\int \iota_V\iota_\xi \,\textrm{d}\theta_0\,f\,\textrm{d}^3r\,\textrm{d}^3p\,\textrm{d}t, \end{equation}

where $\xi = \delta g\circ g^{-1}$ and we have used the identities $\delta V = \dot {\xi } + \mathcal {L}_{V}\xi$ and $\delta (\,f\,\textrm {d}^3r\,\textrm {d}^3p) = -\mathcal {L}_{\xi }(\,f\,\textrm {d}^3r\,\textrm {d}^3p)$, as in the Euler–Poincaré theory developed by Holm, Marsden & Ratiu (Reference Holm, Marsden and Ratiu1998). Now, introducing $g(t)$ and $\delta g(t)$ as in Appendix A, the formula (A2) implies

(3.10)\begin{align} \iota_{\delta g_2}\iota_{\delta g_1}\varOmega_g & ={-}\lim_{a\rightarrow 0}\frac{1}{2a} \delta I(g)[\mathbb{I}_{[{-}a,a]}\delta g]\nonumber\\ & ={-}\lim_{a\rightarrow 0}\frac{1}{2a}\int_{{-}a}^{a}\int \iota_{V}\iota_{\xi} \, \textrm{d}\theta_0\,f\,\textrm{d}^3r\,\textrm{d}^3p\,\textrm{d}t\nonumber\\ & ={-}\int \iota_{\xi_2}\iota_{\xi_1} \,\textrm{d}\theta_0\,f\,\textrm{d}^3r\,\textrm{d}^3p. \end{align}

To invert the $2$-form $\varOmega$ and obtain the Poisson bracket $\{\boldsymbol {\cdot },\boldsymbol {\cdot }\}$ on the space of Lagrangian configuration maps, we first compute the general Hamiltonian vector field $X_F = \dot {g}_F = V_F\circ g$ associated with a real-valued Hamiltonian functional $F(g)$. By definition, $X_F$ satisfies

(3.11)\begin{equation} \iota_{X_F}\varOmega= \textrm{d}F.\end{equation}

See Appendix A for a definition of exterior derivative of a scalar functional. In particular, if $W = \delta g = \xi \circ g$ is any vector on $g$-space then we must have $\iota _W\iota _{X_F}\varOmega = \iota _W \,\textrm {d}F$. After defining the functional derivative $\delta F/\delta g$ according to

(3.12)\begin{equation} \iota_W\,\textrm{d}F = \int \frac{\delta F}{\delta g}\boldsymbol{\cdot} \delta g\,\overset{\circ}{f}\,\textrm{d}^3r_0\,\textrm{d}^3p_0, \end{equation}

such that $\delta F/\delta g$ is a $1$-form at $g(\boldsymbol {r}_0,\boldsymbol {p}_0)$, we therefore infer that Hamilton's equation for $X_F$ is equivalent to

(3.13)\begin{equation} \int \iota_{\xi}\iota_{V_F}\omega_0\,f\,\textrm{d}^3r\,\textrm{d}^3p = \int \iota_{\xi}\left(\frac{\delta F}{\delta g}\circ g^{{-}1}\right)\,f\,\textrm{d}^3r\,\textrm{d}^3p, \end{equation}

for each vector field $\xi$ on the Eulerian phase space. Since $f$ is positive everywhere, we conclude $X_F$ is determined by

(3.14)\begin{equation} V_F = \mathcal{J}_0\boldsymbol{\cdot} \left(\frac{\delta F}{\delta g}\circ g^{{-}1}\right), \end{equation}

where $\mathcal {J}_0$ is the Poisson tensor associated with $\omega _0$. The following manipulation now immediately gives the Poisson bracket:

(3.15)\begin{align} \{G,F\} &= \mathcal{L}_{X_F}G\nonumber\\ &= \int \frac{\delta G}{\delta g}\boldsymbol{\cdot} \dot{g}_F\,\overset{\circ}{f}\,\textrm{d}^3r_0\,\textrm{d}^3p_0\nonumber\\ &= \int \left(\frac{\delta G}{\delta g}\circ g^{{-}1}\right)\boldsymbol{\cdot}\mathcal{J}_0\boldsymbol{\cdot}\left(\frac{\delta F}{\delta g}\circ g^{{-}1}\right)\,f\,\textrm{d}^3r\,\textrm{d}^3p. \end{align}

Notice that if $F(g) = \mathcal {F}(g_*\overset {\circ }{f}\,\textrm {d}^3r_0\,\textrm {d}^3p_0)$ and $G(g) = \mathcal {G}(g_*\overset {\circ }{f}\,\textrm {d}^3r_0\,\textrm {d}^3p_0)$ for some real-valued functionals $\mathcal {F},\mathcal {G}$ on the space of volume forms, then their Poisson bracket is given by

(3.16)\begin{equation} \{G,F\}(g) = \int \,\textrm{d}\,\frac{\delta \mathcal{G}}{\delta f}\boldsymbol{\cdot}\mathcal{J}_0\boldsymbol{\cdot} \textrm{d}\frac{\delta \mathcal{F}}{\delta f}\,f\,\textrm{d}^3r\,\textrm{d}^3p = \{\mathcal{F},\mathcal{G}\}_{LP}(g_* \overset{\circ}{f}\,\textrm{d}^3r_0\,\textrm{d}^3p_0),\end{equation}

where $\{\boldsymbol {\cdot },\boldsymbol {\cdot }\}_{LP}$ is the well-known Lie–Poisson bracket on the space of volume forms. Equation (3.16) says that the mapping $g\mapsto g_*(\overset {\circ }{f}\,d^3r_0\,d^3p_0)$ is a Poisson map between the space of Lagrangian configuration maps with the bracket $\{\boldsymbol {\cdot },\boldsymbol {\cdot }\}$ and the space of volume forms equipped with the Lie–Poisson bracket. We remark, however, that the image of this Poisson map is not equal to the space of volume forms, since pushforward preserves the total number of particles of $\overset {\circ }{f}$.

We can repeat the above procedure for the full multi-species Vlasov–Maxwell system. The infinite-dimensional phase space now comprises tuples of the form $z=(g_\sigma , \boldsymbol {A},\boldsymbol {E})$, where $g_\sigma$ is the Lagrangian configuration map for species $\sigma$, $\boldsymbol {A}$ is the vector potential and $\boldsymbol {E}$ is the electric field. The space of all such $z$ is denoted $Z$. In this case the Lagrange 1-form paired with the tangent vector $\dot z := (\dot {g}_\sigma , \dot {\boldsymbol {A}}, \dot {\boldsymbol {E}} )$ is

(3.17)\begin{equation} \iota_{\dot{z}}{\theta}_{\text{MV}\,z} = \sum_\sigma\int \iota_{V_\sigma}\theta_{B\sigma} { f}_\sigma\,\textrm{d}^3r\,\textrm{d}^3p - \frac{1}{4{\rm \pi} c}\int \boldsymbol{E}\boldsymbol{\cdot}\dot{\boldsymbol{A}}\,\textrm{d}^3r. \end{equation}

The formula (A2) with $\alpha = \theta _{\text {MV}}$ then implies that the symplectic form $\varOmega _{\text {MV}} = -\textrm {d}\theta _{\text {MV}}$ is given by

(3.18)\begin{align} \iota_{\delta z_2}\iota_{\delta z_1}\varOmega_{\text{MV}\,z}& = \sum_\sigma\int \iota_{\xi_{\sigma2}}\iota_{\xi_{\sigma1}}\omega_{B\sigma}\,f_\sigma\,\textrm{d}^3r\,\textrm{d}^3p\nonumber\\ & \quad +\frac{1}{4{\rm \pi} c}\int \left( \delta\boldsymbol{E}_1\boldsymbol{\cdot}\delta{\boldsymbol{A}}_2-\delta{\boldsymbol{E}}_2 \boldsymbol{\cdot}\delta\boldsymbol{A}_1\right)\textrm{d}^3r\nonumber\\ &\quad -\sum_\sigma \int \frac{e_\sigma}{c}\left(\delta\boldsymbol{A}_1\boldsymbol{\cdot} \xi_{\sigma2}^{\boldsymbol{r}} - \delta{\boldsymbol{A}}_2\boldsymbol{\cdot} \xi_{\sigma1}^{\boldsymbol{r}}\right)\,f_\sigma\,\textrm{d}^3r\,\textrm{d}^3p. \end{align}

To find the general Hamiltonian vector field $X_F$ associated with a Hamiltonian $F(g_\sigma ,\boldsymbol {E},\boldsymbol {A})$, we first record the definition of the functional derivatives $\delta F/\delta g_\sigma$, $\delta F/\delta \boldsymbol {E}$ and $\delta F/\delta \boldsymbol {A}$:

(3.19)\begin{equation} \iota_{\delta z}\,\textrm{d}F = \sum_\sigma \int \frac{\delta F}{\delta g_\sigma}\boldsymbol{\cdot} \delta g_\sigma\,\overset{\circ}{f}_\sigma\,\textrm{d}^3r_0\,\textrm{d}^3p_0 + \int \frac{\delta F}{\delta \boldsymbol{E}}\boldsymbol{\cdot} \delta\boldsymbol{E}\,\textrm{d}^3r + \int \frac{\delta F}{\delta \boldsymbol{A}}\boldsymbol{\cdot} \delta\boldsymbol{A}\,\textrm{d}^3r. \end{equation}

Then, we impose the condition $\iota _{\delta z}\iota _{X_F}\varOmega _{\text {MV}} = \iota _{\delta z}\,\textrm {d}F$ for all vectors $\delta z$, which leads to the following expression for $X_F=(V_{F\sigma }\circ g_\sigma ,\dot {\boldsymbol {E}}_F,\dot {\boldsymbol {A}}_F)$:

(3.20)\begin{gather} V_{F\sigma} = \mathcal{J}_{B\sigma}\boldsymbol{\cdot}\left(\frac{\delta F}{\delta g_\sigma}\circ g_\sigma^{{-}1} - 4{\rm \pi} e_\sigma \frac{\delta F}{\delta\boldsymbol{E}}\boldsymbol{\cdot} \textrm{d}\boldsymbol{r}\right) \end{gather}

(3.21)\begin{gather}\dot{\boldsymbol{E}}_F ={-}4{\rm \pi}\sum_\sigma e_\sigma \int \left[\mathcal{J}_{B\sigma}\boldsymbol{\cdot}\left(\frac{\delta F}{\delta g_\sigma}\circ g_\sigma^{{-}1}\right)\right]^{\boldsymbol{r}}\,f_\sigma\,\textrm{d}^3p + 4{\rm \pi} c \frac{\delta F}{\delta\boldsymbol{A}} \end{gather}

(3.22)\begin{gather}\dot{\boldsymbol{A}}_F ={-}4{\rm \pi} c\frac{\delta F}{\delta \boldsymbol{E}}. \end{gather}

Finally, we again use the identity $\{G,F\}_{\text {MV}} = \mathcal {L}_{X_F}G = \iota _{X_F}\,\textrm {d}G$ to find that the Poisson bracket for the Vlasov–Maxwell system is given by

(3.23)\begin{align} &\{G,F\}_{\text{MV}} \nonumber\\ &\quad = \sum_\sigma\int \left(\frac{\delta G}{\delta g_\sigma}\circ g_{\sigma}^{{-}1} - 4{\rm \pi} e_\sigma \frac{\delta G}{\delta\boldsymbol{E}}\boldsymbol{\cdot} \textrm{d}\boldsymbol{r}\right) \boldsymbol{\cdot}\mathcal{J}_{B\sigma}\boldsymbol{\cdot}\left(\frac{\delta F}{\delta g_\sigma} \circ g_{\sigma}^{{-}1} - 4{\rm \pi} e_\sigma \frac{\delta F}{\delta\boldsymbol{E}}\boldsymbol{\cdot} \textrm{d}\boldsymbol{r}\right)\,f_\sigma \,\textrm{d}^3r\,\textrm{d}^3p \nonumber\\ &\qquad + 4{\rm \pi} c\int \left(\frac{\delta G}{\delta\boldsymbol{E}}\boldsymbol{\cdot}\frac{\delta F}{\delta\boldsymbol{A}} - \frac{\delta G}{\delta\boldsymbol{A}}\boldsymbol{\cdot}\frac{\delta F}{\delta \boldsymbol{E}}\right)\textrm{d}^3r. \end{align}

This bracket is a special case of equation (3.252) in Burby (Reference Burby2015). It is related to the well-known Morrison–Marsden–Weinstein bracket (see Morrison Reference Morrison1980; Marsden & Weinstein Reference Marsden and Weinstein1982) on the space of tuples $(\,f_\sigma ,\boldsymbol {E},\boldsymbol {B})$ by the Poisson mapping

(3.24)\begin{equation} (g_\sigma,\boldsymbol{E},\boldsymbol{A})\mapsto (g_{\sigma*}(\overset{\circ}{f}_\sigma\,\textrm{d}^3r_0\,\textrm{d}^3p_0),\boldsymbol{E},\boldsymbol{\nabla}\times\boldsymbol{A}). \end{equation}

As is well known, the Hamiltonian for the Vlasov–Maxwell system is the sum of particle and field energies,

(3.25)\begin{equation} H_{\text{MV}}(g_\sigma,\boldsymbol{E},\boldsymbol{A}) = \sum_\sigma \int m_\sigma c^2\,\gamma_\sigma\,f_\sigma\,\textrm{d}^3r\,\textrm{d}^3p + \frac{1}{8{\rm \pi}}\int (|\boldsymbol{E}|^2 + |\boldsymbol{\nabla}\times\boldsymbol{A}|^2)\,\textrm{d}^3r,\end{equation}

where $\gamma _\sigma = \sqrt {1 + |\boldsymbol {p}|^2/(m_\sigma c)^2}$ is the Lorentz factor for species $\sigma$.

3.2. Why is dark dynamics Hamiltonian?

We will now give a precise explanation for why plasma dynamics on the dark slow manifold possesses a Hamiltonian structure. Further discussion of Hamiltonian structure on slow manifolds more generally is found in the review articles of MacKay (Reference MacKay2004) and Burby (Reference Burby2020a).

Suppose that $(Z,\varOmega )$ is a symplectic manifold with symplectic form $\varOmega$. If $\varLambda \subset Z$ is any submanifold of $Z$ with inclusion map $I_\varLambda :\varLambda \rightarrow Z$, then $\varLambda$ inherits a $2$-form $\varOmega _\varLambda = I_\varLambda ^*\varOmega$, where $I_\varLambda ^*$ denotes the pullback to $\varLambda$. This $2$-form is closed since $\textrm {d}\varOmega _\varLambda = \textrm {d}I_\varLambda ^*\varOmega = I_\varLambda ^*\,\textrm {d}\varOmega = 0$, which means that $\varLambda$ is intrinsically a presymplectic manifold.Footnote ¹ If $\varOmega _\varLambda$ happens to be non-degenerate everywhere on $\varLambda$, then $\varLambda$ is intrinsically a symplectic manifold. As explained in Sniatycki (Reference Sniatycki1974), the Poisson bracket induced by $\varOmega _\varLambda$ in the non-degenerate case coincides with the well-known Dirac bracket for constrained mechanical systems; the submanifold $\varLambda$ represents the constraint.

If $X_H$ is a Hamiltonian vector field on $(Z,\varOmega )$ and $\varLambda$ is an invariant submanifold, then $X_H$ is tangent to $\varLambda$ and the induced vector field $X_{H\varLambda }$ on $\varLambda$ is Hamiltonian with respect to the presymplectic structure $\varOmega _\varLambda$ discussed above. To see this, first note that since $X_H$ is Hamiltonian we have $\iota _{X_H}\varOmega = \textrm {d}H$ in $Z$. Pulling back this equation to $\varLambda$ along $I_\lambda$ gives $\iota _{X_{H\varLambda }}\varOmega _\varLambda = \textrm {d}H_\varLambda$ where $H_\varLambda = I_\varLambda ^*H$, which says that $X_{H\varLambda }$ is a Hamiltonian vector field on $\varLambda$ with Hamiltonian $H_\varLambda$. In this sense, the dynamics on invariant manifolds contained in symplectic manifolds always inherits its own intrinsic Hamiltonian structure.

As we explained in § 2.2, the Eulerian, gauge-invariant form of the Vlasov–Maxwell system contains a (formal) invariant manifold $S$ equal to the dark slow manifold. We constructed $S$ as the constraint set $\boldsymbol {E}_T = \boldsymbol {E}_T^\star (\,f_\sigma ,\boldsymbol {E}_L)$, $\boldsymbol {B} = \boldsymbol {B}_T^\star (\,f_\sigma ,\boldsymbol {E}_L)$ inside of $(\,f_\sigma ,\boldsymbol {E},\boldsymbol {B})$-space. Since the Eulerian, gauge-independent phase space $(\,f_\sigma ,\boldsymbol {E},\boldsymbol {B})$ is related to the Lagrangian, gauge-dependent phase space $(g_\sigma ,\boldsymbol {E},\boldsymbol {A})$ by elimination of the gauge and relabelling degrees of freedom, we may ‘undo’ the gauge and relabelling symmetries to construct an invariant manifold $\tilde {S}$ in $z=(g_\sigma ,\boldsymbol {E},\boldsymbol {A})$-space that projects onto the dark slow manifold $S$. To wit, $\tilde {S}\subset Z$ is defined by the constraints $\boldsymbol {A}_T = \boldsymbol {A}_T^{\star \star }(g_\sigma ,\boldsymbol {E}_L)$ and $\boldsymbol {E}_T =\boldsymbol {E}_T^{\star \star }(g_\sigma ,\boldsymbol {E}_L)$, where the slaving functions $\boldsymbol {A}_T^{\star \star }$ and $\boldsymbol {E}_T^{\star \star }$ are given by

(3.26)\begin{gather} \boldsymbol{A}_T^{{\star}{\star}}(g_\sigma,\boldsymbol{E}_L) = \boldsymbol{A}_T^{{\star}}(\,f_\sigma,\boldsymbol{E}_L) \end{gather}

(3.27)\begin{gather}\boldsymbol{E}_T^{{\star}{\star}}(g_\sigma,\boldsymbol{E}_L) = \boldsymbol{E}_T^{{\star}}(\,f_\sigma,\boldsymbol{E}_L) \end{gather}

(3.28)\begin{gather}g_{\sigma *}(\overset{\circ}{f}\,\textrm{d}^3r_0\,\textrm{d}^3p_0) =f_\sigma \,\textrm{d}^3r\,\textrm{d}^3p, \end{gather}

and $\boldsymbol {A}_T^\star (\,f_\sigma ,\boldsymbol {E}_L)$ is the unique transverse vector field on $Q$ whose curl is $\boldsymbol {B}^*(\,f_\sigma ,\boldsymbol {E}_L)$.

Since we have already shown that the gauge-dependent, Lagrangian form of the Vlasov–Maxwell system is a Hamiltonian system on the symplectic manifold $(Z,\varOmega _{\text {MV}})$, the preceding remarks with $\varLambda = \tilde {S}$ imply that the dark plasma dynamics (in gauge-dependent, Lagrangian form) must possess an intrinsic Hamiltonian structure. As we will see, the $2$-form on $\tilde {S}$ happens to be non-degenerate, which implies that the (all-orders) Poisson bracket on $\tilde {S}$ is a Dirac bracket.

To compute this structure, we must find both the Hamiltonian $H_{\tilde {S}}$ and the $2$-form $\varOmega _{\tilde {S}}$ on $\tilde {S}$. To that end, we introduce the inclusion map $I_{\tilde {S}}:(g_\sigma ,\boldsymbol {E}_L,\boldsymbol {A}_L)\mapsto (g_\sigma ,\boldsymbol {E},\boldsymbol {A})$ defined by

(3.29)\begin{gather} \boldsymbol{E} = \boldsymbol{E}_L + \boldsymbol{E}_T^\star(\,f_\sigma,\boldsymbol{E}_L) \end{gather}

(3.30)\begin{gather}\boldsymbol{A} = \boldsymbol{A}_L + \boldsymbol{A}_T^*(\,f_\sigma,\boldsymbol{E}_L). \end{gather}

Then, we pull back the primitive $1$-form $\theta _{\text {MV}}$ for $\varOmega _{\text {MV}}$ along $I_{\tilde {S}}$ to obtain the Lagrange $1$-form $\theta ^\star$ for dark plasma dynamics

(3.31)\begin{align} \iota_{\dot{x}}\theta^\star_x& = \sum_\sigma \int \left(\boldsymbol{p}+\frac{e_\sigma}{c}\boldsymbol{A}_L + \frac{e_\sigma}{c}\boldsymbol{A}_T^\star\right)\boldsymbol{\cdot} V_\sigma^{\boldsymbol{r}}\,f_\sigma\,\textrm{d}^3r\,\textrm{d}^3p\nonumber\\ &\quad - \frac{1}{4{\rm \pi} c}\int (\boldsymbol{E}_L \boldsymbol{\cdot} \dot{\boldsymbol{A}}_L + \boldsymbol{E}_{T}^\star \boldsymbol{\cdot} D_{\boldsymbol{E}_L}\boldsymbol{A}_T^\star[\dot{\boldsymbol{E}}_L])\,\textrm{d}^3r\nonumber\\ &\quad + \frac{1}{4{\rm \pi} c}\sum_\sigma\int ( \boldsymbol{E}_{T}^\star\boldsymbol{\cdot} D_f\boldsymbol{A}_T^\star[\text{div}(V_\sigma f_\sigma)])\,\textrm{d}^3r. \end{align}

Here, $x=(g_\sigma ,\boldsymbol {E}_L)$ parameterizes the dark slow manifold and $V_\sigma ,f_\sigma$ are defined in terms of $g_\sigma$ as in § 3.1. We also pull back the Vlasov–Maxwell Hamiltonian along $I_{\tilde {S}}$ to obtain the Hamilton function $H^\star$ governing the dark plasma dynamics

(3.32)\begin{equation} H^\star(x) = \sum_\sigma \int m_\sigma c^2\,\gamma_\sigma\,f_\sigma\,\textrm{d}^3r\,\textrm{d}^3p + \frac{1}{8{\rm \pi}}\int (|\boldsymbol{E}_L|^2 + |\boldsymbol{E}_T^\star|^2 + |\boldsymbol{\nabla}\times\boldsymbol{A}_T^\star|^2)\,\textrm{d}^3r. \end{equation}

The all-orders evolution equations for a dark plasma are then $\dot {x} = X_{H^\star }(x)$, where

(3.33)\begin{equation} \iota_{X_{H^\star}}\,\textrm{d}\theta^\star{=}{-}\textrm{d}H^\star.\end{equation}

We will refer to (3.33) as the dark plasma Hamilton equations. We remark that the preceding discussion implies that the formula $\dot {x} = X_{H^\star }(x)$ for the dynamics on the slow manifold must be equivalent to the general formula for the slow dynamics $\dot {x} = g_\epsilon (x,y^\star _\epsilon (x))$ to all orders.

At this point, we could in principle compute the Poisson tensor associated with the dark $2$-form $\varOmega ^\star = -\textrm {d}\theta ^\star$ in order to express the dark plasma dynamics in Poisson bracket form. However, such a computation would not be compatible with our desire to identify truncations of the formal power series $X_{H^\star }$ that possess a Hamiltonian structure. The essential issue is that the Poisson tensor associated with $\varOmega ^*$ is an infinite formal power series in $1/c$, naive truncations of which will fail to satisfy the Jacobi identity. This makes finding computable approximations of $X_{H^\star }$ with a true Hamiltonian structure extremely challenging. In the following subsection, we will therefore take a more nuanced approach to expressing the dark plasma dynamics as a Hamiltonian system that makes structure-preserving truncation much simpler.

Our approach will follow the example set in Burby (Reference Burby2017b), where the same issue was addressed in the context of the slow manifold underlying MHD. In particular, we will apply a near-identity non-canonical transformation $\mathcal {T}:x\mapsto \bar {x}$ on the dark slow manifold that causes the power-series expansion of $\bar {\varOmega }^\star = \mathcal {T}_*\varOmega ^\star$ to truncate at finite order. (In fact we will achieve $\bar {\varOmega }^\star = \bar {\varOmega }^\star _0 + ({1}/{c})\bar {\varOmega }^\star _1$.) We say that the transformation $\mathcal {T}$ rectifies the symplectic structure. This will allow us to identify a closed-form expression for the dark Poisson bracket at the cost of introducing some additional complexity into the dark Hamiltonian $\bar {H}^\star = \mathcal {T}_*H^\star$. Then, we will compute the power-series expansion of $\bar {H}^\star$ to the first post-Darwin order. By replacing $\bar {H}^\star$ with its post-Darwin approximation in the dark Hamilton equation $\iota _{\bar {X}_{\bar {H}^\star }}\bar {\varOmega }^\star = \textrm {d}\bar {H}^\star$, while retaining the full (exact) form of $\bar {\varOmega }^\star$, we will obtain a computable post-Darwin Hamiltonian approximation of the dark plasma dynamics.

3.3. Derivation of the rectifying transformation

The Lagrange $1$-form restricted to the dark slow manifold is given in (3.31). Note that the ‘coordinates’ we use on the dark slow manifold are $x = (g_\sigma ,\boldsymbol {E}_L)$, where $g_\sigma$ is the species-$\sigma$ Lagrangian configuration map and $\boldsymbol {E}_L$ is the longitudinal electric field. Since the constraint functions $\boldsymbol {A}_T^\star$ and $\boldsymbol {E}_T^\star$ that define the slow manifold are infinite formal power series in $1/c$, it is clear that $\theta ^\star$ is also an infinite formal power series in $1/c$. We wish to change variables from $x$ to $\bar {x}$ using a near-identity transformation $\mathcal {T}:x\mapsto \bar {x}$ such that the power series defining the dark Lagrange $1$-form truncates modulo exact $1$-forms at finite order when expressed in terms of $\bar {x}=(\bar {g}_\sigma ,\bar {\boldsymbol {E}}_L)$. We will call the transformation $\mathcal {T}$ a rectifying transformation. We work modulo addition of exact $1$-forms to $\theta ^\star$ because it is only the dark $2$-form $\varOmega ^\star = -\textrm {d}\theta ^\star$ that is physically important.

It may not be immediately clear why it should be possible to find a rectifying transformation. The explanation dwells in the fact that symplectic manifolds have no local invariants. This should be contrasted with Riemannian geometry, where the curvature tensor is a local invariant; curvature cannot be eliminated by a coordinate transformation. Thus, symplectic manifolds satisfy the Darboux theorem, which says that, locally, all symplectic forms are related by coordinate transformations. While the existence of a rectifying transformation, which we require to be globally defined, is not implied directly by the Darboux theorem, one can prove existence of such a transformation (as a formal power series) using the idea underlying Moser's celebrated proof of the Darboux theorem. See Burby (Reference Burby2017b) for details.

In order to streamline our derivation, we represent $\mathcal {T}$ as a composition of formal Lie transforms, i.e.

(3.34)\begin{equation} \mathcal{T} = \dots\circ\exp(G_3)\circ \exp(G_2)\circ \exp(G_1),\end{equation}

where the $G_k$ are vector fields on $x$-space (not to be confused with the $\boldsymbol {x}$-space $\mathbb {Q}$!) that we allow to be formal power series in $1/c$. To specify $\mathcal {T}$, we will derive formulas for the $G_k$. Our derivation will be facilitated by the well-known expression for the pushforward of a differential form $\alpha$ along a Lie transform $\exp (G)$

(3.35)\begin{equation} \exp(G)_*\alpha =\alpha -\mathcal{L}_{G}\alpha + \tfrac{1}{2}\mathcal{L}_{G}^2\alpha - \tfrac{1}{6}\mathcal{L}_G^3\alpha + \dots . \end{equation}

In order to work with these infinite-dimensional Lie transforms explicitly, we will always compute Lie derivatives by first applying Cartan's identity $\mathcal {L}_{G} = \iota _{G}\,\textrm {d} + \textrm {d}\iota _G$ and then using the formalism described in Appendix A for computing exterior derivatives on infinite-dimensional spaces.

The first step in our derivation is to list the leading five terms in the power-series expansion of the (pre-transformed) dark Lagrange $1$-form $\theta ^\star$. We have $\theta ^\star =\theta ^\star _0 + c^{-1}\,\theta ^\star _1 + c^{-2}\theta ^\star _2 + \dots$ with

(3.36)\begin{gather} \iota_{\dot{x}}\theta^\star_{0\,x} = \sum_{\sigma} \int \boldsymbol{p}\boldsymbol{\cdot} V_\sigma^{\boldsymbol{r}}\,f_\sigma\,\textrm{d}^3r\,\textrm{d}^3p \end{gather}

(3.37)\begin{gather}\iota_{\dot{x}}\theta^\star_{1\,x} = \sum_\sigma \int e_\sigma\boldsymbol{A}_L\boldsymbol{\cdot} V_\sigma^{\boldsymbol{r}}\,f_\sigma\,\textrm{d}^3r\,\textrm{d}^3p - \frac{1}{4{\rm \pi}}\int \boldsymbol{E}_L\boldsymbol{\cdot}\dot{\boldsymbol{A}}_L\,\textrm{d}^3r \end{gather}

(3.38)\begin{gather}\iota_{\dot{x}}\theta^\star_{2\,x} = \sum_\sigma \int e_\sigma \boldsymbol{A}_{T1}^\star\boldsymbol{\cdot} V_\sigma^{\boldsymbol{r}}\,f_\sigma\,\textrm{d}^3r\,\textrm{d}^3p \end{gather}

(3.39)\begin{gather}\iota_{\dot{x}}\theta^\star_{3\,x} =0 \end{gather}

(3.40)\begin{gather}\iota_{\dot{x}}\theta^\star_{4\,x} = \sum_\sigma \int e_\sigma \boldsymbol{A}_{T3}^\star\boldsymbol{\cdot}V_\sigma^{\boldsymbol{r}}\,f_\sigma\,\textrm{d}^3r\,\textrm{d}^3p + \frac{1}{4{\rm \pi}}\sum_\sigma\int \boldsymbol{E}_{T2}^\star\boldsymbol{\cdot} D_{f_\sigma}\boldsymbol{A}_{T1}^\star[\text{div}(V_\sigma f_\sigma)]\,\textrm{d}^3r, \end{gather}

where we have used $\boldsymbol {A}_{T0}^\star =\boldsymbol {A}_{T2}^\star = \boldsymbol {E}_{T0}^\star = \boldsymbol {E}_{T1}^\star =0$ and $D_{\boldsymbol {E}_L}\boldsymbol {A}_{T1}^\star = 0$. Higher-order contributions to $\theta ^\star$ will not be necessary in our analysis. With the help of (2.33) the second term of (3.40) can be written as

(3.41)\begin{equation} {\sum_\sigma \frac{4{\rm \pi} e_\sigma}{m_\sigma}\int f_\sigma\,\iota_{V_\sigma}\,\textrm{d}\left(\boldsymbol{p}\boldsymbol{\cdot}\varDelta^{{-}1} \boldsymbol{E}_{T2}^\star\right) }, \end{equation}

where we have re-written six-dimensional divergence as a Lie derivative and used the Cartan identity. Next, we record formulas for the exterior derivatives of $\theta ^\star _k$ for $k=0,\dots ,3$. (We will not need $\textrm {d}\theta ^\star _4$.) In deriving these formulas, we reuse the procedure from Appendix A for computing exterior derivatives in infinite dimensions

(3.42)\begin{align} \iota_{\delta x_2}\iota_{\delta x_1}\,\textrm{d}\theta^\star_{0\,x} &= \sum_\sigma\int \iota_{\xi_{\sigma 2}}\iota_{\xi_{\sigma 1}}\,\textrm{d}\theta_0\,f_\sigma\,\textrm{d}^3r\,\textrm{d}^3p\end{align}

(3.43)\begin{align} \iota_{\delta x_2}\iota_{\delta x_1}\,\textrm{d}\theta^\star_{1\,x} &=\sum_\sigma\int e_\sigma\left(\delta\boldsymbol{A}_{L1}\boldsymbol{\cdot}\xi_{\sigma 2}^{\boldsymbol{r}} - \delta\boldsymbol{A}_{L2}\boldsymbol{\cdot}\xi_{\sigma 1}^{\boldsymbol{r}}\right)\,f_\sigma\,\textrm{d}^3r\,\textrm{d}^3p\nonumber\\ &\quad - \frac{1}{4{\rm \pi}}\int\left(\delta\boldsymbol{E}_{L1}\boldsymbol{\cdot}\delta\boldsymbol{A}_{L2} - \delta\boldsymbol{E}_{L2}\boldsymbol{\cdot}\delta\boldsymbol{A}_{L1}\right)\textrm{d}^3r \end{align}

(3.44)\begin{align} \iota_{\delta x_2}\iota_{\delta x_1}\,\textrm{d}\theta^\star_{2\,x} &=\sum_\sigma \int e_\sigma \boldsymbol{B}_1^\star\boldsymbol{\cdot} \xi_{\sigma 1}^{\boldsymbol{r}}\times\xi_{\sigma 2}^{\boldsymbol{r}}\,f_\sigma\,\textrm{d}^3r\,\textrm{d}^3p\nonumber\\ &\qquad +4{\rm \pi}\int \left(\sum_\sigma\frac{e_\sigma}{m_\sigma}\boldsymbol{\nabla}\boldsymbol{\cdot}\int \xi_{\sigma 1}^{\boldsymbol{r}}\boldsymbol{p}\,f_\sigma\,\textrm{d}^3p\right)\boldsymbol{\cdot}\varDelta^{{-}1} \varPi_T\left(\sum_\sigma e_\sigma\int \xi_{\sigma 2}^{\boldsymbol{r}}\,f_\sigma\,\textrm{d}^3p\right)\textrm{d}^3r \nonumber\\ &\qquad -4{\rm \pi}\int \left(\sum_\sigma\frac{e_\sigma}{m_\sigma}\boldsymbol{\nabla}\boldsymbol{\cdot}\int \xi_{\sigma 2}^{\boldsymbol{r}}\boldsymbol{p}\,f_\sigma\,\textrm{d}^3p\right)\boldsymbol{\cdot}\varDelta^{{-}1} \varPi_T\left(\sum_\sigma e_\sigma\int \xi_{\sigma 1}^{\boldsymbol{r}}\,f_\sigma\,\textrm{d}^3p\right)\textrm{d}^3r \nonumber\\ &\qquad -4{\rm \pi}\int \left(\sum_\sigma\frac{e_\sigma}{m_\sigma}\int \xi_{\sigma 1}^{\boldsymbol{p}}\,f_\sigma\,\textrm{d}^3p\right)\boldsymbol{\cdot}\varDelta^{{-}1}\varPi_T\left(\sum_\sigma e_\sigma\int \xi_{\sigma 2}^{\boldsymbol{r}}\,f_\sigma\,\textrm{d}^3p\right)\textrm{d}^3r\nonumber\\ &\qquad +4{\rm \pi}\int \left(\sum_\sigma\frac{e_\sigma}{m_\sigma}\int \xi_{\sigma 2}^{\boldsymbol{p}}\,f_\sigma\,\textrm{d}^3p\right)\boldsymbol{\cdot}\varDelta^{{-}1}\varPi_T\left(\sum_\sigma e_\sigma\int \xi_{\sigma 1}^{\boldsymbol{r}}\,f_\sigma\,\textrm{d}^3p\right)\textrm{d}^3r \end{align}

(3.45)\begin{align} &\iota_{\delta x_2}\iota_{\delta x_1}\,\textrm{d}\theta^\star_{3\,x} =0. \end{align}

We now determine $G_1$ by requiring it solves the linear equation

(3.46)\begin{equation} \iota_{G_1}\,\textrm{d}\theta^\star_0 +c^{{-}1}\iota_{G_1}\,\textrm{d}\theta^\star_1 = c^{{-}2}\theta^\star_2.\end{equation}

The steps necessary to invert this equation have been outlined in (3.14), (3.20), (3.21) and (3.22). (The reason for this choice of $G_1$ will become clear soon.) Using the formulas (3.42), (3.43) and (3.38), we find that $G_1 = (Y_{\sigma 1}\circ g_\sigma , G_1^{\boldsymbol {E}_L},G_1^{\boldsymbol {A}_L})$ is given by

(3.47)\begin{gather} Y_{\sigma1} = \frac{e_\sigma}{c^2}\boldsymbol{A}_{T1}^\star\boldsymbol{\cdot}\partial_{\boldsymbol{p}} \end{gather}

(3.48)\begin{gather}G_1^{\boldsymbol{E}_L} = 0 \end{gather}

(3.49)\begin{gather}G_1^{\boldsymbol{A}_L} = 0. \end{gather}

Since $G_1$ is $O(c^{-2})$ and satisfies (3.46), the pushforward of $\theta ^\star$ along the Lie transform $\exp (G_1)$ has the power-series expansion

(3.50)\begin{equation} \exp(G_1)_*\theta^\star = \theta^\star_0 + c^{{-}1}\theta^\star_1 + c^{{-}4}\theta^\star_4 - \tfrac{1}{2}c^{{-}2}\,\iota_{G_1}\,\textrm{d}\theta^\star_2+ O(c^{{-}5}), \end{equation}

where we have used our freedom to add exact differentials to the Lagrange $1$-form in order to make the replacement $\exp (G_1)_* = \exp (-\mathcal {L}_{G_1})\rightarrow \exp (-\iota _{G_1}d)$. Thus, if we were to stop here and set $G_k = 0$ for $k>1$, we would have succeeded in making the difference between the transformed Lagrange $1$-form $\bar {\theta }^\star$ and $\theta ^\star _0 + c^{-1}\theta ^\star _1$ one order higher than prior to applying the transformation.

Finally, we turn to determining the $G_k$ for $k>1$. Since the first Lie transform, generated by $G_1$, brought the Lagrange $1$-form closer to $\theta ^\star _0 + c^{-1}\,\theta ^\star _1$, we will select the higher $G_k$ so that the transformed Lagrange form is equal to $\theta ^\star _0 + c^{-1}\,\theta ^\star _1$ to all orders in $1/c$. We give an inductive proof in Appendix C that this can indeed be done. For the purposes of this article, it is only necessary for us to explicitly perform the algebra that determines $G_2$. However, all subsequent $G_k$ may be determined in a similar manner. (The tedium of the required algebraic manipulations quickly becomes unreasonable for $k>2$.)

In keeping with the above strategy, we demand that $G_2$ satisfies the linear equation

(3.51)\begin{equation} \iota_{G_2}\,\textrm{d}\theta^\star_0+c^{{-}1}\iota_{G_2}\,\textrm{d}\theta^\star_1 = c^{{-}4}\theta^\star_4 - \tfrac{1}{2}c^{{-}2}\iota_{G_1}\,\textrm{d}\theta_2^\star.\end{equation}

What remains is to estimate the last term. From the last four formidable-looking lines of (3.44), only the third one will contribute to the final result because of the form of $G_1$, as we see from equation (3.51). Using formulas (3.42), (3.43), (3.40), (3.41) and (3.44), we find the solution $G_2 = (Y_{\sigma 2}\circ g_\sigma ,G_2^{\boldsymbol {E}_L},G_2^{\boldsymbol {A}_L})$ is given by

(3.52)\begin{gather} Y_{\sigma 2} =\frac{e_\sigma}{c^4}\left[\boldsymbol{A}_{T3}^\star{+} \boldsymbol{\nabla}\left(\frac{\boldsymbol{p}}{m_\sigma}\boldsymbol{\cdot} \varDelta^{{-}1}\boldsymbol{E}_{T2}^\star\right)+\frac{1}{2}\varDelta^{{-}1}\varPi_T\omega_p^2 \boldsymbol{A}_{T1}^\star\right]\boldsymbol{\cdot}\frac{\partial}{\partial\boldsymbol{p}} \nonumber\\ \quad -\frac{1}{c^4}\left[\frac{e_\sigma}{m_\sigma}\varDelta^{{-}1}\boldsymbol{E}_{T2}^\star\right] \boldsymbol{\cdot}\frac{\partial}{\partial\boldsymbol{r}} \end{gather}

(3.53)\begin{gather} G_2^{\boldsymbol{E}_L} = \frac{1}{c^4}\varPi_L\omega_p^2\varDelta^{{-}1}\boldsymbol{E}_{T2}^\star \end{gather}

(3.54)\begin{gather}G_2^{\boldsymbol{A}_L} = 0. \end{gather}

These formulas show in particular that $G_2 = O(c^{-4})$. Since it can be shown that $G_k = O(c^{-5})$ for $k> 2$, this, together with (3.51), implies

(3.55)\begin{equation} \bar{\theta}^\star = \theta^\star_0 + c^{{-}1}\,\theta^\star_1 + O(c^{{-}5}), \end{equation}

consistent with our claim that $\bar {\theta }^\star = \theta ^\star _0 + c^{-1}\,\theta ^\star _1$ to all orders in $1/c$.

To summarize, we have identified a near-identity non-canonical transformation on the dark slow manifold that rectifies the dark Poisson bracket. Using the variables $(g_\sigma , \boldsymbol {E}_L,\boldsymbol {A}_L)$ to parameterize the slow manifold, the transformation is given explicitly to $O(c^{-4})$ by

(3.56)\begin{align} {\boldsymbol{A}}_L &=\bar{\boldsymbol{A}}_L \end{align}

(3.57)\begin{align}{\boldsymbol{E}}_L &=\bar{\boldsymbol{E}}_L - \sum_\sigma\frac{1}{c^4} \varPi_L \bar{\omega}_{p\sigma}^2\varDelta^{{-}1} \overline{\boldsymbol{E}}_{T2}^\star{+} O(c^{{-}5}) \end{align}

(3.58)\begin{align} {{f}}_\sigma & =\bar{{f}}_\sigma + \frac{e_\sigma}{c^2}\bar{\boldsymbol{A}}_{T1}^\star\boldsymbol{\cdot} \partial_{\boldsymbol{p}}\bar{{f}}_\sigma+\frac{1}{2}\frac{e_\sigma^2}{c^4}\bar{\boldsymbol{A}}_{T1}^*\bar{ \boldsymbol{A}}_{T1}^*:\partial_{\boldsymbol{p}}^2\bar{{f}}_\sigma\nonumber\\ &\quad +\frac{e_\sigma}{c^4}\left[\bar{\boldsymbol{A}}_{T3}^\star{+} \frac{1}{2} \varDelta^{{-}1}{\varPi_T}\bar{\omega}_{p}^2 \bar{\boldsymbol{A}}_{T1}^\star{+} \varDelta^{{-}1} \boldsymbol{\nabla} \frac{\boldsymbol{p}\boldsymbol{\cdot} \bar{\boldsymbol{E}}_{T2}^\star}{m_\sigma} \right]\boldsymbol{\cdot} \partial_{\boldsymbol{p}}\bar{{f}}_\sigma\nonumber\\ &\quad - \frac{e_\sigma}{m_\sigma c^4} \varDelta^{{-}1} \bar{\boldsymbol{E}}^\star_{T2}\boldsymbol{\cdot}\boldsymbol{\nabla} \bar{f}_\sigma+O(1/c^5). \end{align}

These formulas follow from the definition of a Lie transform; if $G = (Y_\sigma \circ g_\sigma ,G^{\boldsymbol {E}_L},G^{\boldsymbol {A}_L})$ is the Lie generator then $\exp (G)(g_\sigma ,\boldsymbol {E}_L,\boldsymbol {A}_L) = (\bar {g}_\sigma ,\bar {\boldsymbol {E}}_L,\bar {\boldsymbol {A}}_L)$ is the $\lambda = 1$ solution of the system of evolution equations

(3.59a–c)\begin{equation} \frac{\textrm{d}g_\sigma}{\textrm{d}\lambda} = Y_\sigma\circ g_\sigma,\quad \frac{\textrm{d}\boldsymbol{E}_L}{\textrm{d}\lambda} = G^{\boldsymbol{E}_L},\quad \frac{\textrm{d}\boldsymbol{A}_L}{\textrm{d}\lambda} = G^{\boldsymbol{A}_L}. \end{equation}

After applying this rectifying transformation, the dark Lagrange $1$-form is given by $\bar {\theta }^\star = \theta ^\star _0 + c^{-1}\,\theta ^\star _1$ to all orders in $1/c$. Therefore, the (rectified) dark symplectic form $-\textrm {d}\bar {\theta }^\star =\bar {\varOmega }^\star$ is given by

(3.60)\begin{align} {\iota_{\delta x_2}\iota_{\delta x_1}}\bar{\varOmega}^\star & = \sum_\sigma\int \omega_0[\xi_{\sigma 1},\xi_{\sigma 2}]\,f_\sigma\,\textrm{d}^3r\,\textrm{d}^3p + \frac{1}{4{\rm \pi} c}\int (\delta\boldsymbol{E}_{L1}\boldsymbol{\cdot}\delta\boldsymbol{A}_{L2} - \delta\boldsymbol{E}_{L2}\boldsymbol{\cdot}\delta\boldsymbol{A}_{L1})\,\textrm{d}^3r\nonumber\\ &\quad -\sum_\sigma\int \frac{e_\sigma}{c}(\delta \boldsymbol{A}_{L1}\boldsymbol{\cdot}\xi_{\sigma 2}^{\boldsymbol{r}} - \delta\boldsymbol{A}_{L2}\boldsymbol{\cdot}\xi_{\sigma 1}^{\boldsymbol{r}})\,f_\sigma\,\textrm{d}^3r\,\textrm{d}^3p, \end{align}

and the associated dark Poisson bracket is

(3.61)\begin{align} \{\mathcal{F},\mathcal{G}\}^\star &=\sum_\sigma\int \,\textrm{d}^6 z \, f_\sigma \left[\frac{\delta \mathcal{F}}{\delta g_\sigma}\circ g^{{-}1}_\sigma -4{\rm \pi} e \frac{\delta \mathcal{F}}{\delta \boldsymbol{E}_L}\boldsymbol{\cdot} \textrm{d}\boldsymbol{r}\right]\boldsymbol{\cdot} \mathcal{J}_0 \boldsymbol{\cdot} \left[\frac{\delta \mathcal{G}}{\delta g_\sigma}\circ g^{{-}1}_\sigma -4{\rm \pi} e \frac{\delta \mathcal{G}}{\delta \boldsymbol{E}_L}\boldsymbol{\cdot} \textrm{d}\boldsymbol{r}\right] \nonumber\\ &\quad + 4{\rm \pi} c\int \,\textrm{d}^3r \left(\frac{\delta \mathcal{F}}{\delta \boldsymbol{E}_L}\boldsymbol{\cdot} \frac{\delta \mathcal{G}}{\delta \boldsymbol{A}_L} - \frac{\delta \mathcal{G}}{\delta \boldsymbol{E}_L}\boldsymbol{\cdot} \frac{\delta \mathcal{F}}{\delta \boldsymbol{A}_L} \right). \end{align}

We remark that the functional derivative with respect to a longitudinal vector field is defined (consistently and in accordance with the standard mathematical notion of a functional derivative) to be a longitudinal vector field. Therefore in particular if $\boldsymbol {W}$ is some vector field then $\delta /\delta \boldsymbol {E}_L\int \boldsymbol {E}_L\boldsymbol {\cdot } \boldsymbol {W}\,\textrm {d}^3 r = \varPi _L\boldsymbol {W}$, not $\boldsymbol {W}$. We also remark that the bracket $\{\boldsymbol {\cdot },\boldsymbol {\cdot }\}^\star$ is written in terms of Lagrangian configuration maps; the corresponding Eulerian bracket is

(3.62)\begin{align} \{\mathcal{F},\mathcal{G}\}_{\text{Dark}} &= \sum_\sigma\int \,\textrm{d}^6 z \, f_\sigma \left[\textrm{d}\frac{\delta \mathcal{F}}{\delta f_\sigma} -4{\rm \pi} e \frac{\delta \mathcal{F}}{\delta \boldsymbol{E}_L}\boldsymbol{\cdot} \textrm{d}\boldsymbol{r}\right]\boldsymbol{\cdot} \mathcal{J}_0 \boldsymbol{\cdot} \left[\textrm{d}\frac{\delta \mathcal{G}}{\delta f_\sigma} - 4{\rm \pi} e \frac{\delta \mathcal{G}}{\delta \boldsymbol{E}_L}\boldsymbol{\cdot} \textrm{d}\boldsymbol{r}\right] \nonumber\\ &\quad + 4{\rm \pi} c\int \,\textrm{d}^3r \left(\frac{\delta \mathcal{F}}{\delta \boldsymbol{E}_L}\boldsymbol{\cdot} \frac{\delta \mathcal{G}}{\delta \boldsymbol{A}_L} - \frac{\delta \mathcal{G}}{\delta \boldsymbol{E}_L}\boldsymbol{\cdot} \frac{\delta \mathcal{F}}{\delta \boldsymbol{A}_L} \right). \end{align}

3.4. Derivation of the post-Darwin Hamiltonian

In the previous subsection, we succeeded in identifying a rectifying transformation that allowed us to find a closed-form expression for the Poisson bracket on the slow manifold. We are now in position to derive approximations of the dark plasma evolution equations that possess a Hamiltonian structure. The remaining step is to derive a formula for the slow manifold Hamiltonian expressed in terms of the rectified variables $(\bar {g}_\sigma ,\bar {\boldsymbol {E}}_L,\bar {\boldsymbol {A}}_L)$. The purpose of this section is to derive that Hamiltonian $\bar {H}^\star$ to the first post-Darwin order. Our post-Darwin approximation to the dark plasma dynamics is then defined by the Hamilton equation $\dot {x} = \{x,\bar {H}^\star _{\textrm {PD}}\}_{\text {Dark}}$, where $\bar {H}^\star _{\textrm {PD}}$ is the post-Darwin Hamiltonian and the rectified dark Poisson bracket is given in (3.61).

The original Hamiltonian (3.25) reads

(3.63)\begin{equation} H^\star(\,f_\sigma,\boldsymbol{E}_L) = \sum_\sigma\int \,\textrm{d}^6 z\, f_\sigma \gamma_\sigma m_\sigma c^2 + \frac{1}{8{\rm \pi}} \int \,\textrm{d}^3r \left( E^2_L + E^{{\star} 2}_T + B^{{\star} 2} \right) \end{equation}

and can be expended in powers of $1/c$ as $H^\star (\,f,E_L) = c^2 H_{-2}^\star + H_0^\star + c^{-2} H_2^\star + c^{-4} H_2$ with the two lowest-order terms given as

(3.64a,b)\begin{equation} H^\star_{{-}2} = \sum_\sigma \int \,\textrm{d}^6 z\, f_\sigma m_\sigma c^2,\quad H^\star_{0} =\sum_\sigma \int \,\textrm{d}^6 z\, \frac{f_\sigma p^2}{2m_\sigma} + \frac{1}{8{\rm \pi}}\int \,\textrm{d}^3r\, E_L^2, \end{equation}

the (quadratic in $f_\sigma$) Darwin terms given as

(3.65)\begin{equation} H^\star_2 ={-} \sum_\sigma\int \,\textrm{d}^6 z \, \frac{f_\sigma p^4}{8 m^3_\sigma} + \frac{1}{8{\rm \pi}}\int \,\textrm{d}^3r \, {B}^{{\star} 2}_{1}, \end{equation}

and the ‘piezoelectric terms’ given as

(3.66)\begin{equation} H^\star_4 = \sum_\sigma\int \,\textrm{d}^6 z \frac{f_\sigma p^6}{16 m^5_\sigma} + \frac{1}{8{\rm \pi}} \int \,\textrm{d}^3 r\, |\boldsymbol{E}_{T2}^\star|^{ 2} + \frac{1}{4{\rm \pi}} \int \,\textrm{d}^3 r\, \boldsymbol{B}_{1}^\star\boldsymbol{\cdot} \boldsymbol{B}_{3}^\star. \end{equation}

Note that the second term in $H^\star _4$ is quadratic in $\,f_{\sigma }$, while the last term is cubic. The cubic term, which involves $\boldsymbol {B}^\star _3$, appears particularly formidable (cf. (2.46)). As we will soon see, much of the complexity introduced by this term is eliminated by the rectifying transformation (3.58). Also note that we are free to remove the rest-energy term $H^\star _{-2}$ since it is a Casimir.

By defining $\boldsymbol { {\rm \pi}}_\sigma := \boldsymbol {p} - ({e_\sigma }/{c^2}) \boldsymbol {A}_{T1}^\star$ and substituting the inverse transformation formulas (3.56)–(3.58) into the Hamiltonian on the slow manifold, we obtain the following expression for the post-Darwin Hamiltonian in terms of the rectified variables $(\bar {\boldsymbol {E}}_L,\bar {\boldsymbol {A}}_L,\bar {f}_\sigma )$:

(3.67)\begin{align} {\bar{H}^\star_{\text{PD}}}(\bar{\boldsymbol{E}}_L,\bar{f}_\sigma) &=\sum_\sigma \int \left(\frac{|\bar{\boldsymbol{\rm \pi}}_\sigma|^2}{2m_\sigma} - \frac{|\bar{\boldsymbol{\rm \pi}}_\sigma|^4}{8 m_\sigma^3 c^2}+ \frac{|\bar{\boldsymbol{\rm \pi}}_\sigma|^6}{16 m_\sigma^5c^4}\right)\,\bar{f}_\sigma\,\textrm{d}^3r\,\textrm{d}^3p\nonumber\\ &\quad -\frac{2{\rm \pi}}{c^2}\int \left[\varPi_T\bar{\boldsymbol{J}}_0 + \varPi_T\frac{\bar{\omega}_p^2}{c^2}\varDelta^{{-}1}\varPi_T\bar{\boldsymbol{J}}_0\right] \varDelta^{{-}1}\left[\varPi_T\bar{\boldsymbol{J}}_0+ \varPi_T\frac{\bar{\omega}_p^2}{c^2}\varDelta^{{-}1}\varPi_T\bar{\boldsymbol{J}}_0\right]\textrm{d}^3r\nonumber\\ &\quad-\frac{{1}}{8{\rm \pi} c^4}\int |\bar{\boldsymbol{E}}_{T2}^\star|^2\,\textrm{d}^3r +\frac{1}{8{\rm \pi}}\int |\bar{\boldsymbol{E}}_L|^2\,\textrm{d}^3r +\frac{2{\rm \pi}}{ c^4}\int \bar{\omega}_p^2\,|\varDelta^{{-}1}\varPi_T\bar{\boldsymbol{J}}_0|^2\,\textrm{d}^3r, \end{align}

where special care was taken to write the expression so that it is manifestly positive definite. The integrand on the first line comprises the first three non-constant terms of the Taylor expansion of $\sqrt {1+|\bar {\boldsymbol {{\rm \pi} }}_\sigma |^2/m_\sigma ^2 c^2}$, which is readily seen to be positive. The second line involves a $O(c^{-6})$ term which will be dropped below to obtain the $1/c$ expansion of the equations of motion, but otherwise it is quartic in $f$. Modulo $O(c^{-6})$ terms, the second line and the last term in the third line can be condensed into

(3.68)\begin{equation} -\frac{1}{8{\rm \pi} c^2} \int \,\textrm{d}^3 r\,\bar{\boldsymbol{A}}_{T1}^\star \boldsymbol{\cdot} \left(\varDelta + \frac{\bar{\omega}_p^2}{c^2} \right) \bar{\boldsymbol{A}}_{T1}^\star. \end{equation}

If we combine everything together and take variational derivatives we get

(3.69)\begin{align} \frac{\delta \bar{H}^\star_{\textrm{PD}}}{\delta \bar{f}_\sigma} &= {\frac{p^2}{2m_\sigma}-\frac{p^4}{8m_\sigma^3 c^2}+\frac{p^6}{16m_\sigma^5c^4}} {-\frac{e_\sigma \boldsymbol{p}\boldsymbol{\cdot} \boldsymbol{A}^\star_{T1}}{m_\sigma c^2}\left(1-\frac{p^2}{2m^2_\sigma c^2} \right)}\nonumber\\ & \quad + \frac{4{\rm \pi} e_\sigma}{m_\sigma c^4} \boldsymbol{p}\boldsymbol{\cdot} \varDelta^{{-}1}\varPi_T J_2 -\frac{{ e_\sigma} }{c^4}\left( \frac{\boldsymbol{p p}}{{m_\sigma}^2} : \varDelta^{{-}1}\boldsymbol{\nabla} \,\bar{\boldsymbol{E}}_{T2}^\star{+} \frac{e_\sigma}{m_\sigma} \bar{\boldsymbol{E}}_L\boldsymbol{\cdot} \varDelta^{{-}1} \bar{\boldsymbol{E}}_{T2}^\star \right), \end{align}

where the square roots indicate the expansion introduced in (3.67) and we have used equation (2.30), while the variational derivative with respect to the other remaining slow variable reads

(3.70)\begin{equation} \frac{\delta \bar{H}^\star_{\textrm{PD}}}{\delta \bar{\boldsymbol{E}}_L} = \frac{\bar{ \boldsymbol{E}}_L}{4{\rm \pi}} - \frac{{1}}{4{\rm \pi} c^4} {\varPi_L}\bar{\omega}_p^2 \varDelta^{{-}1} \bar{\boldsymbol{E}}_{T2}^\star. \end{equation}

This completes the derivation of the post-Darwin Hamiltonian.

Using the Poisson bracket (3.62), the corresponding equations of motion on the slow manifold are

(3.71)\begin{gather} \partial_t \bar{f}_\sigma + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\left[\partial_{\boldsymbol{p}}\frac{\delta \bar{H}^\star_{\textrm{PD}}}{\delta f_\sigma}\right]\bar{f}_\sigma\right) + \partial_{\boldsymbol{p}}\boldsymbol{\cdot} \left(\left[4{\rm \pi} e_\sigma \frac{\delta \bar{H}^\star_{\textrm{PD}}}{\delta\boldsymbol{E}_L} - \boldsymbol{\nabla}\frac{\delta \bar{H}^\star_{\textrm{PD}}}{\delta f_\sigma}\right]\bar{f}_\sigma\right) = 0 \end{gather}

(3.72)\begin{gather}\partial_t \bar{\boldsymbol{E}}_L +4{\rm \pi}\varPi_L\sum_\sigma e_\sigma \int \left(\partial_{\boldsymbol{p}}\frac{\delta \bar{H}^\star_{\textrm{PD}}}{\delta \bar{f}_\sigma}\right)\,\bar{f}_\sigma\,\textrm{d}^3p = 4 {\rm \pi}c\frac{\delta \bar{H}^\star_{\textrm{PD}}}{\delta \bar{{\boldsymbol{A}}_L}} \end{gather}

(3.73)\begin{gather}\partial_t\bar{{\boldsymbol{A}}}_L ={-}4{\rm \pi} c\frac{\delta \bar{H}^\star_{\textrm{PD}}}{\delta \bar{\boldsymbol{E}}_L}. \end{gather}

Note that, since $\bar {H}^\star _{\textrm {PD}}$ is independent of $\boldsymbol {A}_L$, (3.72) reproduces the longitudinal Ampère equation, which is equivalent to the preservation of Gauss's law. This result will persist to all orders in $c^{-1}$ as a result of the gauge invariance of the rectifying transformation. For the same reason, (3.71) and (3.72) comprise a closed system of evolution equations for the distribution function $\bar {f}_\sigma$ and the longitudinal electric field $\bar {\boldsymbol {E}}_L$. As such, the evolution equations for the dark manifold dynamics may be written in terms of the same dependent variables used in the Vlasov–Poisson system.

3.5. On the Hamiltonian nature of Darwin's approximation

We can collect all terms up to $1/c^2$ obtained in the previous section to get the transformed Darwin's Hamiltonian

(3.74)\begin{align} \bar{H}_D^\star(\bar{f}_\sigma,\bar{\boldsymbol{E}}_L) &= \sum_\sigma \int \,\textrm{d}^6 z\, \frac{\bar{f}_\sigma p^2}{2m_\sigma} + \frac{1}{8{\rm \pi}}\int \,\textrm{d}^3r\, \bar{E}_L^2 \nonumber\\ &\quad -\sum_\sigma\int \,\textrm{d}^6 z \, \frac{\bar{f}_\sigma p^4}{8 m^3_\sigma c^2} +\frac{2{\rm \pi}}{c^2}\int \,\textrm{d}^3r\, \bar{\boldsymbol{J}}_{0} \boldsymbol{\cdot} \varDelta^{{-}1} \varPi_T\, \bar{\boldsymbol{J}}_{0}. \end{align}

We note here that, to this order, the transformations (3.58) amount to a transformation from kinetic to canonical momentum.

While the Hamiltonian approach to the post-Darwin extension appears to be new, there is an instance of a Hamiltonian study of Darwin's approximation found in the literature, namely Krause, Apte & Morrison (Reference Krause, Apte and Morrison2007). In particular, in this work, action principles for the Darwin approximation in the Vlasov context were presented. However, the phase space only spans the field of density functions $f$, so that the Lie–Poisson structure is that of a Vlasov–Poisson bracket:

(3.75)\begin{equation} \{F, G\}_{D}=\int \,\textrm{d}^{6} \zeta \,f_{D}\left[\frac{\delta F}{\delta \,f_{D}}, \frac{\delta G}{\delta \,f_{D}}\right], \end{equation}

which naturally produces Vlasov equation ${\partial \,f_{D}}/{\partial t}=\{\,f_{D}, H\}$ when equipped with the Hamiltonian

(3.76)\begin{align} H\left[\,f_{D}\right]&= \int \,\textrm{d}^{6} \zeta \,f_{D}(\zeta, t) \left[\frac{{\rm \pi}^{2}}{2 m}-\frac{{\rm \pi}^{4}}{8 m^{3} c^{2}}\right] \nonumber\\ &\quad +\frac{e^{2}}{2} \int \,\textrm{d}^{6} \zeta \int \,\textrm{d}^{6} \zeta^{\prime} \,f_{D}(\zeta, t) \,f_{D}\left(\zeta^{\prime}, t\right) K\left(\boldsymbol{r} | \boldsymbol{r}^{\prime}\right) \nonumber\\ &\quad -\frac{e^{2}}{2 m^{2} c^{2}} \int \,\textrm{d}^{6} \zeta \int \,\textrm{d}^{6} \zeta^{\prime} \,f_{D}(\zeta, t) \,f_{D}\left(\zeta^{\prime}, t\right) K_{i j}\left(\boldsymbol{r} | \boldsymbol{r}^{\prime}\right) {\rm \pi}_{i} {\rm \pi}_{j}^{\prime}, \end{align}

where ${\rm \pi} _i$ correspond to Eulerianized full canonical momenta. To make a further comparison, we recognize from (15), (16) and (23) of Krause et al. (Reference Krause, Apte and Morrison2007) that $\int \,\textrm {d}^3 r^\prime \, K_{ij}(\boldsymbol {r} | \boldsymbol {r}^{\prime }) \equiv -4{\rm \pi} \varDelta ^{-1}\varPi _T$ and that they have single species. In addition, in our case (3.67), $\boldsymbol { {\rm \pi}}$ is a kinetic momentum up to $O(1/c^2)$, since the active transformations we performed in (3.58) can be viewed as a passive relabelling transformation in $\boldsymbol {p}$. So the ${\rm \pi} _i$ in Krause et al. (Reference Krause, Apte and Morrison2007) correspond to our new $\boldsymbol {p}$ and not $\boldsymbol {{\rm \pi} }$.

It is straightforward to recover this Hamiltonian formulation of the Darwin equation from ours by restricting to a level set of the momentum map associated with gauge symmetry (i.e. the residual of Gauss's law) and then quotienting by translations in $\boldsymbol {A}_L$.