% - \Mu \nabla (B_0 + \mathfrak{B}_{1,\shortparallel})
- \Mu \nabla B_0
\\
\dot{\gy{U}}_\shortparallel &=
\frac{1}{m_s} \vc{b}_s^\effective \bcdot \roundpar{
q_s \vc{\mathcal{E}}_1^\reduced
% - \Mu \nabla (B_0 + \mathfrak{B}_{1,\shortparallel})
- \Mu \nabla B_0
\end{align}

where the reduced FLR corrected electric field is defined as (cf. <ref>)

\begin{equation}\label{eq:reduced_flc_electricfield}
\vc{\mathcal{E}}_1^\reduced \defeq -\nabla \gav{\evalperturbedgc{\phi}_1} - \pard{\gav{\evalperturbedgc{A}_{1,\shortparallel}}}{t} \hb.
\end{equation}

Note that the FLR corrected electromagnetic fields now coincide exactly with the appropriately gyro-averaged electromagnetic fields, since the reduced model neglects exactly the $\mathcal{O}(\varepsilon_B)$ terms still present in <ref>.

For the derivation of the field equations we again give Low's action explicitly resulting in (cf. <ref>)

\begin{multline}\label{eq:low_split_reduced_explicit}
\mathfrak{A}^\reduced(\gy{\vc{Z}}, \phi_1, \Apa)
= \sum_s \int \gy{f}_s(\gyzzero, \tzero) \biggl[
q_s \roundpar{
\vc{A}_{0,s}^\effective
+ \gav{\evalperturbedgc{A}_{1,\shortparallel}} \hb} \bcdot \dot{\gy{\vc{R}}}
+ \frac{m_s \gy{\Mu}}{q_s} \dot{\gy{\Theta}}
% - \gy{\Mu} (B_0 + \mathfrak{B}_{1,\shortparallel})
- \gy{\Mu} B_0
- \frac{m_s}{2} \gy{U}_\shortparallel^2
\\
- q_s \gav{\evalperturbedgc{\phi}_1}
\biggr] \jacobiangy(\gyzzero) \diff^6 \gyzzerod \diff t
+ \sum_s \int \gyfzero \biggl(
\frac{m_s}{2 B_0^2} \abs{\nabla_\perp (\phi_1 - \gy{u}_\shortparallel \Apa)}^2
- \frac{\gy{\mu}}{2 B_0} \abs{\nabla_\perp \Apa}^2
\biggr) \jacobiangcslab \diff^6 \gy{z} \diff t\\
- \frac{1}{2 \mu_0} \int \abs{\vc{B}_0 + \nabla_\perp \Apa \times \hb}^2 \diff^3 x \diff t.
\end{multline}

The quasi-neutrality equation (cf. <ref>) and Ampère's law (cf. <ref>) are given by

\begin{multline}
= \sum_s \int \gy{f}_s
q_s \gav{\cp} \dot{\gy{R}}_\shortparallel
\jacobiangy \diff^6 \gy{z} \diff t\\
+ \sum_s \int \gyfzero \biggl(
- \frac{m_s \gy{u}_\shortparallel}{B_0^2} \nabla_\perp (\phi_1 - \gy{u}_\shortparallel \Apa) \bcdot \nabla_\perp \cp
- \frac{\gy{\mu}}{B_0} \nabla_\perp \Apa \bcdot \nabla_\perp \cp
\biggr) \jacobiangcslab \diff^6 \gy{z} \diff t
- \frac{1}{\mu_0} \int \nabla_\perp \Apa \bcdot \nabla_\perp \cp \diff^3 x \diff t.
\end{multline}

\begin{align}
- \sum_s \int \vc{\mathcal{P}}_{1,s}^\reduced \bcdot \nabla_\perp \cp \diff^3 \gy{r}
\sum_s q_s \int \gy{f}_s
% \roundpar{
\gav{\evalperturbedgc{\cp}}
% - \frac{\gy{\mu}\mathfrak{C}(1)}{1 + \mathfrak{C}(\gy{\mu}/q)} \cp
% }
\jacobiangy \diff^6 \gy{z}
\frac{1}{\mu_0} \int \nabla_\perp \Apa \bcdot \nabla_\perp \cp \diff^3 x
\sum_s \int \vc{\mathcal{M}}_{1,s}^\reduced \bcdot (\nabla_\perp \cp \times \hb) \diff^3 \gy{r}
+ \sum_s q_s \int \gy{f}_s
% \roundpar{
\dot{\gy{R}}_\shortparallel \gav{\evalperturbedgc{\cp}}
% - \frac{\gy{\mu} \mathfrak{C}(\gy{u}_\shortparallel)}{1 + \mathfrak{C}(\gy{\mu}/q)} \cp
% }
\jacobiangy \diff^6 \gy{z},\label{eq:reduced_ampere}
\end{align}

where the reduced polarization (cf. <ref>) and magnetization (cf. <ref>) are given by

\begin{align}
\vc{\mathcal{P}}_{1,s}^\reduced &\defeq \int \gyfzero \vc{P}_{1,s}^\reduced \jacobiangcslab \diff^3 \gy{u}
, \quad
\vc{P}_{1,s}^\reduced \defeq \frac{m_s}{q_s B_0^2} \vc{F}_{1,\perp}^\reduced
, \quad
\vc{F}_{1,\perp}^\reduced \defeq -q_s \nabla_\perp (\phi_1 - \gy{u}_\shortparallel \Apa),\label{eq:reduced_polarization}\\
\vc{\mathcal{M}}_{1,s}^\reduced &\defeq \int \gyfzero \vc{M}_{1,s}^\reduced \jacobiangcslab \diff^3 \gy{u}
, \quad
\vc{M}_{1,s}^\reduced \defeq
- \hb \times \roundpar{\gy{u}_\shortparallel \vc{P}_{1,s}^\reduced
- \frac{\gy{\mu}}{B_0} \nabla_\perp \Apa}
\end{align}

When considering the EOMs (<ref>) as well as field equations (<ref>) of the reduced model, we find that this model coincides with reduced parallel models from the literature.
We note that this is only due to the choice $(\gparamr, \gparamtheta) = (1, 0)$, which yields the appropriate gyro-averages on both the scalar and vector potential.

§.§.§ Strong formulation

In order to be able to interpret the equations more easily, we present the strong formulation of the field equations as follows

\begin{align}
-\frac{1}{\mu_0} \nabla \bcdot \nabla_\perp \Apa + \nabla \bcdot \squarepar{\mathcal{C}(\gy{u}_\shortparallel) \nabla_\perp \phi_1 + \frac{p_{0,\perp} - p_{0,\shortparallel}}{B_0^2} \nabla_\perp \Apa} &=  \evalperturbedgcadjoint{\mathcal{J}}{}^{\free,\reduced}_\shortparallel,\\
- \nabla \bcdot \squarepar{\mathcal{C}(1) \nabla_\perp \phi_1 - \mathcal{C}(\gy{u}_\shortparallel) \nabla_\perp \Apa} &= \evalperturbedgcadjoint{\mathcal{R}}{}^\free,
\end{align}

where we have defined the spatially varying functions $\mathcal{C}(\zeta)$ as

\begin{equation}\label{eq:coeffC}
\mathcal{C}(\zeta) \defeq \sum_s \frac{m_s}{B_0^2} \int \gyfzero \zeta \jacobiangcslab \diff^3 \gy{u}
\end{equation}

The gyro-average adjoint of the parallel component of the free current $\evalperturbedgcadjoint{{\mathcal{J}}}{}^{\free,\reduced}_\shortparallel$ is different from the one discussed in <ref>, and it is defined in a weak sense as (cf. <ref>)

\begin{equation}
\int \evalperturbedgcadjoint{{\mathcal{J}}}{}^{\free,\reduced}_\shortparallel \cp \diff^3 \gy{r}
\defeq
\sum_s q_s \int \gy{f}_s \dot{\gy{{R}}}_\shortparallel \gav{\evalperturbedgc{\cp}}  \jacobiangy \diff^6 \gy{z}.
\end{equation}

The field equations (<ref>) become decoupled when a centred Maxwellian background particle distribution function is considered, which is defined as

\begin{equation}\label{eq:maxwellian}
\gyfzerocm \defeq \frac{\gy{\mathcal{N}}_{00,s}}{\sqrt{\pi}^{3} \uths^3}
\exp{\left(-\frac{m_s \gy{u}_\shortparallel^2 + 2 B_0 \gy{\mu}}{m_s \uths^2} \right)}.
\end{equation}

The following identities, which hold for a centred Maxwellian background distribution, result in significant simplifications

\begin{equation}\label{eq:centeredMax_integrate0}
\int \gyfzero \gy{u}_\shortparallel \jacobiangcslab \diff^3 \gy{u} = 0
, \quad
\int \gyfzero \roundpar{m_s \gy{u}_\shortparallel^2 - \gy{\mu} B_0} \jacobiangcslab \diff^3 \gy{u} = 0.
\end{equation}

Note that <ref>, in physical terms, coincides with the absence of a parallel background current as well as the isotropy of the temperature.
In particular, this results in symmetry $\mathcal{C}(\gy{u}_\shortparallel) = 0$ and isotropy $p_{0,\perp} = p_{0,\shortparallel}$, thereby fully decoupling the field equations (<ref>).

§.§ Symplectic gyrokinetic model from Brizard–Hahm

Besides the gauge invariant model described in [Burby & Brizard, 2019], there are also gauge variant gyrokinetic models which include $\Adperp$.
In particular, we consider the symplectic model from <cit.>, which we hereafter refer to as the Brizard–Hahm (BH) model.
In our terminology, this model coincides with the choice of gyrocentre coordinates given by <ref>.

§.§.§ The Hamiltonian

The first- and second-order Hamiltonian are derived from <cit.> resulting in (cf. <ref>)

\begin{equation}\label{eq:H1_bh}
\gy{H}_{1,s}^\bh \defeq q_s \gav{\evalperturbedgc{\phi}_1} + \Mu \dgav{\evalperturbedgcburby{B}_{1,\shortparallel}}
\end{equation}

and  (cf. <ref>)

\begin{equation}\label{eq:H2_bh}
\gy{H}_{2,s}^{\bh} \defeq
- \frac{m_s}{2 B_0^2} \abs{\nabla_\perp (\phi_1 - \gy{U}_\shortparallel A_{1,\shortparallel})}^2
+ \frac{\Mu}{2 B_0} \abs{\nabla_\perp \Apa}^2,
\end{equation}

The derivation, which neglects terms of $\mathcal{O}(\varepsilon_\perp^3)$ in $\gy{H}_{2,s}^{\bh}$, can be found in <ref>.
We note that this result can alternatively be derived from <ref> with $\gparamr = 1$, when making use of the approximation $\vc{B}_1 \approx \nabla_\perp \Apa \times \hb + (\curl \Adperp)_\shortparallel \hb$ and neglecting ${\partial \vc{A}_{1,\perp}} / {\partial t}$.

Note that the second-order Hamiltonian coincides exactly with that of the parallel-only model: $\gy{H}_{2,s}^{\bh} = \gy{H}_{2,s}^\reduced$.
This is remarkable, as it implies an absence of the polarization current as well as an absence of a contribution to the magnetization from the perpendicular part of the vector potential.

Both our proposed gauge invariant model and the Brizard–Hahm model [Brizard & Hahm, 2007] reduce to the parallel-only model when neglecting the perpendicular part of the vector potential.
Note that the same cannot be said about the gauge invariant model from [Burby & Brizard, 2019], which coincides with $(\gparamr, \gparamtheta) = (0, 0)$.
Therein, we find that, e.g. $\vc{\mathcal{A}}_1 = \Ad$ such that the vector potential without gyro-averaging appears in the EOMs.

§.§.§ Principle of least action

The EOMs are found by substituting the expression for the particle Hamiltonian $\gy{H}^{\mathrm{part},\bh} = \gy{H}_0 + \gy{H}_1^\bh$ into the general form of the EOMs given by <ref>, where we now have $\vc{\mathcal{A}}_1 = \gav{\evalperturbedgc{\vc{A}}_1}$, which therefore also changes the definition of the effective gyrocentre magnetic field $\vc{B}_s^\effective$ given by <ref>.
This results in (cf. <ref>)

\begin{align}
\dot{\gy{\vc{R}}} &= \gy{U}_\shortparallel \vc{b}_s^\effective - \frac{1}{q_s B_{s,\shortparallel}^\effective} \hb \times \squarepar{q_s \vc{\mathcal{E}}_1^\bh - \Mu \nabla \roundpar{B_0 + \dgav{\evalperturbedgcburby{B}_{1,\shortparallel}}}},\label{eq:gy_eoms_bh_R}\\
\dot{\gy{U}}_\shortparallel &= \frac{1}{m_s} \vc{b}_s^\effective \bcdot \squarepar{q_s \vc{\mathcal{E}}_1^\bh - \Mu \nabla \roundpar{B_0 + \dgav{\evalperturbedgcburby{B}_{1,\shortparallel}}}},\label{eq:gy_eoms_bh_Vpa}
\end{align}

where the reduced FLR corrected electric field is defined as (cf. <ref>)

\begin{equation}\label{eq:bh_flc_electricfield}
\vc{\mathcal{E}}_1^\bh \defeq -\nabla \gav{\evalperturbedgc{\phi}_1} - \pard{\gav{\evalperturbedgc{\vc{A}}_{1}}}{t}.
\end{equation}

Note that

the EOMs given by <ref> coincide with the ones found in <ref> when $\varepsilon_B = 0$.

The action is explicitly given by

\begin{multline}\label{eq:low_split_bh_explicit}
\mathfrak{A}^\bh(\gy{\vc{Z}}, \phi_1, \Ad, \lambda)
= \sum_s \int \gy{f}_s(\gyzzero, \tzero) \biggl[
q_s \roundpar{
\vc{A}_{0,s}^\effective
+ \gav{\evalperturbedgc{\vc{A}}_{1}}} \bcdot \dot{\gy{\vc{R}}}
+ \frac{m_s \gy{\Mu}}{q_s} \dot{\gy{\Theta}}
% - \gy{\Mu} (B_0 + \mathfrak{B}_{1,\shortparallel})
- \gy{\Mu} B_0
- \frac{m_s}{2} \gy{U}_\shortparallel^2 \\
- q_s \gav{\evalperturbedgc{\phi}_1} - \Mu \dgav{\evalperturbedgcburby{B}_{1,\shortparallel}}
\biggr] \jacobiangy(\gyzzero) \diff^6 \gyzzerod \diff t
+ \sum_s \int \gyfzero \biggl(
\frac{m_s}{2 B_0^2} \abs{\nabla_\perp (\phi_1 - \gy{u}_\shortparallel \Apa)}^2\\
- \frac{\gy{\mu}}{2 B_0} \abs{\nabla_\perp \Apa}^2
\biggr) \jacobiangcslab \diff^6 \gy{z} \diff t
- \frac{1}{2 \mu_0} \int \abs{\vc{B}_0 + \vc{B}_1}^2 \diff^3 x \diff t
+ \int C^\bh(\gy{\vc{Z}}, \phi_1, \Ad, \lambda) \diff^3 x \diff t,
\end{multline}

where a Lagrange multiplier $\lambda$ is introduced, with an associated constraint $C^\bh$, to ensure that a well-posed system of equations is found, as will be discussed in more detail in <ref>.
Therein, we will show that the constraint necessarily depends on all unknowns of our problem, including the characteristics $\gy{\vc{Z}}$.
Hence, when using a variational formulation of the Brizard–Hahm model [Brizard & Hahm, 2007] we find that well-posedness of the model implies that the Lagrange multiplier $\lambda$ affects the EOMs as well as each of the field equations.
We do not show this dependence here explicitly, and for now we ignore the contribution due to the constraint.

The action (<ref>) results in the same quasi-neutrality equation as found for the reduced parallel-only model, as given by <ref>.
We find that Ampère's law is now given by

\begin{align}
\frac{1}{\mu_0} \int (\vc{B}_0 + \vc{B}_1) \bcdot \curl \ca \diff^3 x
= {} & \sum_s \int \nabla \bcdot (\vc{\mathcal{M}}_{1,s}^\reduced \times \hb) \capa \diff^3 \gy{r} \nonumber \\ &
+ \sum_s \int \gy{f}_s \squarepar{q_s \gav{\evalperturbedgc{\ca}} \bcdot \dot{\gy{\vc{R}}} - \gy{\mu} \dgav{(\evalperturbedgcburby{\nabla} \times \ca_\perp)_\shortparallel}} \jacobiangy \diff^6 \gy{z}
% + \int G^\bh(\lambda, \ca) \diff^3 x
\end{align}

§.§.§ Strong formulation

As usual, we consider the strong formulation of the field equations, this results in the following Ampère's law and quasi-neutrality equation (cf. <ref>)

\begin{align}
\frac{1}{\mu_0} \curl \vc{B}
&= \overbrace{\sum_s \squarepar{\hb \nabla \bcdot (\vc{\mathcal{M}}_{1,s}^\reduced \times \hb)}}^{\vc{\mathcal{J}}^{\mathrm{b},\bh} = \vc{\mathcal{J}}^{\mathrm{m},\bh}}
+ \evalperturbedgcadjoint{\vc{\mathcal{J}}}{}^{\free,\bh}
% + \nabla_\perp \lambda
, \label{eq:bh_ampere}\\
\sum_s \nabla \bcdot \vc{\mathcal{P}}_{1,s}^\reduced &= \evalperturbedgcadjoint{\mathcal{R}}{}^\free.\label{eq:bh_gauss}
\end{align}

The meaning of the gyro-average adjoint of the free current $\evalperturbedgcadjoint{\vc{\mathcal{J}}}{}^{\free,\bh}$ is different from the one discussed in <ref>, and it is defined in a weak sense as (cf. <ref>)

\begin{equation}
\int \evalperturbedgcadjoint{\vc{\mathcal{J}}}{}^{\free,\bh} \bcdot \ca \diff^3 \gy{r}
\defeq
\sum_s q_s \int \gy{f}_s \squarepar{\dot{\gy{\vc{R}}} \bcdot \gav{\evalperturbedgc{\ca}} - \gy{\mu} \dgav{(\evalperturbedgcburby{\nabla} \times \ca)_\shortparallel}} \jacobiangy \diff^6 \gy{z}.
\end{equation}

This free current adjoint coincides with the free current adjoint from our proposed model, as given by <ref>, whenever $\varepsilon_B = 0$ thanks to <ref>.

Contrary to our proposed gauge invariant model discussed in strong formulation in <ref>, we now fail to recognize the structure of the macroscopic Maxwell's equations in <ref>.
This is primarily due to the absence of a polarization current, but also due to the fact that the magnetization current $\vc{\mathcal{J}}^{\mathrm{m},\bh}$ cannot be written as the curl of a magnetization.
Even when neglecting $\mathcal{O}(\varepsilon_B)$ contributions this is not possible, in which case we find

\begin{equation}
\vc{\mathcal{J}}{^\mathrm{m,\bh}} = (\curl \vc{\mathcal{M}}_{1,s}^\reduced)_\shortparallel + \mathcal{O}(\varepsilon_B).
\end{equation}

In order to further simplify the field equations we assume that our background particle distribution function is given by a centred Maxwellian (<ref>).
The resulting field equations are given by

\begin{align}
\frac{1}{\mu_0} \curl (\vc{B}_0 + \curl \Ad)
&= \evalperturbedgcadjoint{\vc{\mathcal{J}}}{}^{\free,\bh}
, \label{eq:bh_cm_ampere}\\
-\nabla \bcdot (\mathcal{C}(1) \nabla_\perp \phi_1) &= \evalperturbedgcadjoint{\mathcal{R}}{}^\free.\label{eq:bh_cm_gauss}
\end{align}

§.§.§ Well-posedness of the field equations

For the symplectic Brizard–Hahm model [Brizard & Hahm, 2007] we note that the free current density is defined differently from the free current density of our proposed model (<ref>).
This implies, in particular, that the gyro-average of the free charge continuity equation is no longer satisfied

\begin{equation}\label{eq:gy_continuity_gavadjoint_bh}
% \int \roundpar{\pard{\mathcal{R}^\free}{t} + \nabla \bcdot \vc{\mathcal{J}}^\free} \gav{\evalperturbedgc{\cp}} \diff^3 \gy{r} = 0
% \quad \iff \quad
\int \roundpar{\pard{\evalperturbedgcadjoint{\mathcal{R}}{}^\free}{t} + \nabla \bcdot \evalperturbedgcadjoint{\vc{\mathcal{J}}}{}^{\free,\bh}} \cp \diff^3 \gy{r}
= \int \vc{\mathcal{J}}^\free \bcdot \roundpar{\nabla \gav{\evalperturbedgc{\cp}} - \gav{\evalperturbedgc{\nabla} \cp}} \diff^3 \gy{r}
\end{equation}

Since the right-hand side does not vanish in general we find that a result analogous to <ref> is not obtained.

Upon inspecting <ref>, we find that the bound current is given by

\begin{equation}
\vc{\mathcal{J}}^{\bound,\bh}
= \vc{\mathcal{J}}^{\mathrm{m},\bh}
= \sum_s \hb \nabla \bcdot (\vc{\mathcal{M}}_{1,s}^\reduced \times \hb),
\end{equation}

which is not divergence free.
Computation of the divergence of Ampère's law (<ref>), with the contribution from the Lagrange multiplier $\lambda$, results in the constraint

\begin{equation}\label{eq:total_continuity_bh}
= \pard{}{t} (\evalperturbedgcadjoint{\mathcal{R}}{}^\free + \mathcal{R}^{\bound,\bh})
+ \nabla \bcdot (\evalperturbedgcadjoint{\vc{\mathcal{J}}}{}^{\free,\bh} + \vc{\mathcal{J}}^{\bound,\bh} + \mathcal{C}(1) \nabla_\perp \lambda)
\end{equation}

upon addition of the quasi-neutrality equation (<ref>) where we let $\mathcal{R}^{\bound,\bh} \defeq - \sum_s \nabla \bcdot \vc{\mathcal{P}}_{1,s}^\reduced$.
In this case, the gyro-average adjoint of the free charge continuity equation does not vanish, and we are therefore left with a total continuity equation which is enforced by means of a Lagrange multiplier.
Due to <ref> we find that the constraint <ref> also depends on the characteristics $\gy{\vc{Z}}$, which is very undesirable as it implies that the Lagrange multiplier $\lambda$ also affects the EOMs.

Hence, for the symplectic Brizard–Hahm model [Brizard & Hahm, 2007] we find that $\lambda \neq 0$ for three independent reasons: the polarization current is missing, the magnetization current is not the curl of a magnetization, and the gyro-average adjoint of the free charge continuity equation does not hold.
Note that a polarization current can be included by considering higher-order approximations of the first-order generating function, see e.g. [Qin et al., 1999].

§.§ Gyrokinetic Darwin approximation

The proposed gauge-invariant gyrokinetic model keeps more physics than the popular reduced parallel-only model discussed in <ref> and has more structure than the fully-symplectic Brizard–Hahm model [Brizard & Hahm, 2007].
In particular, we find that Ampère's law (<ref>) now contains a polarization current, which is not present in either of the two other models.
Part of this polarization current, however, gives rise to a fast compressional Alfvén wave as will be demonstrated in <ref>, and in most situations such waves are undesired due to their high frequency.

One way to deal with such high frequency components is to damp them numerically using implicit time integration methods.
Another option is to remove such waves from the underlying model, i.e. from the Lagrangian.
For instance, to remove light-speed waves from Maxwell's equations, rather than considering the limit of quasi-neutrality, one can consider the Darwin approximation, wherein only the contribution due to the scalar potential is kept in the $\epsilon_0 \partial \vc{E}_1 / \partial t$ term in Ampère's law.
Thus, the solenoidal, i.e. divergence free, part (when assuming that the Coulomb gauge is used) of the displacement current is neglected, resulting in a model which yields a second- and third-order accurate electric and magnetic field, respectively, in the small parameter $\varepsilon_c = v / c$ [Degond & Raviart, 1992].

§.§.§ The Hamiltonian

We follow an approach similar to the Darwin approximation to remove the fast compressional Alfvén waves.
Recall that the second-order gyrocentre Hamiltonian is given by (cf. <ref>, with $(\gparamr, \gparamtheta) = (1, 0)$)

\begin{equation}
\gy{H}_{2,s} =
- \frac{m_s}{2 B_0^2} \abs{\vc{E}_{1,\perp}}^2
+ \frac{m_s \gy{U}_\shortparallel}{B_0^2} \hb \bcdot (\vc{E}_{1,\perp} \times \vc{B}_1)
+ \frac{
\gy{\Mu} B_0
- m_s \gy{U}_\shortparallel^2
}{2 B_0^2} {\abs{\vc{B}_{1,\perp}}^2},
\end{equation}

where we have substituted the definition of the Lorentz force (<ref>), and this result thereby agrees exactly with <cit.>.
We note that the compressional Alfvén wave comes from the $\abs{\partial \Adperp / \partial t}^2$ term (which itself comes from the $\abs{\vc{E}_{1,\perp}}^2$ term), and it should therefore be sufficient to remove exactly this term from $\gy{H}_{2,s}$.

However, when well-posedness of the resulting Darwin model is considered we find that the $\hb \bcdot (\vc{E}_1 \times \vc{B}_1)$ term from $\gy{H}_{2,s}$ must also be omitted.
The contribution to the field Lagrangian due to the second-order part of the Hamiltonian involves integration over the background particle distribution function and therefore omitting this term is justified whenever

\begin{equation}
\sum_s m_s \int \gyfzero \gy{u}_\shortparallel \jacobiangcslab \diff^3 \gy{u} \approx 0.
\end{equation}

Note that this constraint is satisfied whenever we consider a nearly symmetric background distribution as given by <ref> with $\varepsilon_{U,s} \ll 1$.
As we require $\varepsilon_{U,s} \ll 1$ already for justifying $\varepsilon_\delta \ll 1$, we find that the term can be omitted from $\gy{H}_{2,s}$ without introducing any additional assumptions.

When removing both the compressional Alfvén wave and the $\hb \bcdot (\vc{E}_1 \times \vc{B}_1)$ term from $\gy{H}_{2,s}$ we find

\begin{equation}\label{eq:gy_darwin_H2}
\gy{H}_{2,s}^\darwin \defeq
- \frac{m_s}{2 B_0^2} \roundpar{\abs{\nabla_\perp \phi_1}^2 + 2 \nabla_\perp \phi_1 \bcdot \pard{\Adperp}{t}}
% + \frac{m_s \gy{U}_\shortparallel}{B_0^2} \hb \bcdot (\vc{E}_{1,\perp} \times \vc{B}_1)
+ \frac{
\gy{\Mu} B_0
- m_s \gy{U}_\shortparallel^2
}{2 B_0^2} {\abs{\vc{B}_{1,\perp}}^2},
\end{equation}

where we note that this an approximation breaks gauge invariance.

§.§.§ Principle of least action

As the particle part of the Hamiltonian is still given by <ref>, and the symplectic part of the Lagrangian is unchanged when compared to the gauge invariant model, we find that the EOMs are still given by <ref>.

The explicit expression of Low's action, while using <ref> and omitting the contribution due to the gauge condition (cf. <ref>), is given by

\begin{multline}\label{eq:low_split_explicit_darwin}
\mathfrak{A}^\darwin(\gy{\vc{Z}}, \phi_1, \Ad, \lambda)
= \sum_s \int \gyfzero(\gyzzero) \biggl[
q_s \roundpar{
\vc{A}_{0,s}^\effective
+ \vc{A}_1
+ \radgav{\evalperturbedgcburby{\vc{B}}_1\times \vc\rho}} \bcdot \dot{\gy{\vc{R}}}
+ \frac{m_s \gy{\Mu}}{q_s} \dot{\gy{\Theta}} \\
- \gy{\Mu} (B_0 + \dgav{\evalperturbedgcburby{B}_{1,\shortparallel}})
- \frac{m_s}{2} \gy{U}_\shortparallel^2
- q_s \phi_1
+ q_s \rho \radgav{\evalperturbedgcburby{E}_{1,\rho}}
\biggr] \jacobiangy(\gyzzero, \tzero) \diff^6 \gyzzerod \diff t\\
+ \sum_s \int \gyfzero \biggl[
\frac{m_s}{2 B_0^2} \roundpar{\abs{\nabla_\perp \phi_1}^2 + 2 \nabla_\perp \phi_1 \bcdot \pard{\Adperp}{t}}
% - \frac{m_s \gy{u}_\shortparallel}{B_0^2} \hb \bcdot (\vc{E}_{1,\perp} \times \vc{B}_{1,\perp})\\
- \frac{\gy{\mu} B_0 - m_s \gy{u}_\shortparallel^2}{2 B_0^2} \abs{\vc{B}_{1,\perp}}^2
\biggr] \jacobiangcslab \diff^6 \gy{z} \diff t\\
- \frac{1}{2 \mu_0} \int \abs{\vc{B}_0 + \vc{B}_1}^2 \diff^3 x \diff t
+ \int \mathcal{C}(1) \nabla_\perp \lambda \bcdot \Ad \diff^3 x \diff t
\end{multline}

Again, a Lagrange multiplier $\lambda$ is introduced which ensures well-posedness of the resulting equations.
In this case, we find that the constraint depends linearly on $\lambda$ and $\Ad$, and therefore only affects Ampère's law.
This constraint is motivated in <ref>.

The action (<ref>) results in the same quasi-neutrality equation as found in <ref>.

The Darwin approximation of Ampère's law is given by (cf. <ref>)

\begin{align}
\frac{1}{\mu_0} \int (\vc{B}_0 + \vc{B}_1) \bcdot \curl \ca \diff^3 x
% - \int G(\chi, \ca) \diff^3 x
= {} &
\sum_s \int \roundpar{\pard{\vc{\mathcal{P}}_{1,s}^\darwin}{t} + \curl \vc{\mathcal{M}}_{1,s}} \bcdot \ca \diff^3 \gy{r}
+ \int \evalperturbedgcadjointalt{\vc{\mathcal{J}}}{}^\free \bcdot \ca \diff^3 \gy{r} \nonumber\\ &
+ \int \mathcal{C}(1) \nabla_\perp \lambda \bcdot \ca \diff^3 x
\end{align}

where the Darwin polarization is given by (cf. <ref>)

\begin{equation}\label{eq:darwin_polarization}
\vc{\mathcal{P}}_{1,s}^\darwin \defeq \int \gyfzero \vc{P}_{1,s}^\darwin \jacobiangcslab \diff^3 \gy{u}
, \quad
\vc{P}_{1,s}^\darwin \defeq \frac{m_s}{q_s B_0^2} \vc{F}_{1,\perp}^\darwin
, \quad
% \vc{F}_{1,\perp}^\darwin \defeq q_s \roundpar{- \nabla_\perp \phi_1 + \gy{u}_\shortparallel \hb \times \vc{B}_1}.
\vc{F}_{1,\perp}^\darwin \defeq - q_s \nabla_\perp \phi_1.
\end{equation}

Setting the variation w.r.t. the Lagrange multiplier $\lambda$ to zero results in the following constraint

\begin{equation}\label{eq:darwin_weak_constraint}
\int \mathcal{C}(1) \nabla_\perp \cp \bcdot \Ad \diff^3 x = 0.
\end{equation}

§.§.§ Well-posedness of the macroscopic Darwin equations

We find that the resulting model can be written in strong formulation as

\begin{align}
\curl \vc{\mathcal{H}} &= \pard{\vc{\mathcal{D}}^\darwin}{t} + \evalperturbedgcadjointalt{\vc{\mathcal{J}}}{}^\free + \mathcal{C}(1) \nabla_\perp \lambda,\label{eq:darwin_ampere_strong}\\
% \pard{\vc{B}_1}{t} &= -\curl \vc{E}_1,\label{eq:darwin_faraday}\\
\nabla \bcdot \roundpar{\vc{\mathcal{D}}^\darwin - \mathcal{C}(1) \pard{\vc{A}_{1,\perp}}{t}} &= \evalperturbedgcadjoint{\mathcal{R}}{}^\free,\label{eq:darwin_gauss}
% \nabla \bcdot \roundpar{\Adperp \sum_s \int \gyfzero \frac{m_s}{B_0^2} \jacobiangcslab \diff^3 \gy{u}} &= 0,
% \nabla \bcdot \vc{B}_1 &= 0,\label{eq:darwin_maggauss}
\end{align}

where the Darwin displacement field is defined as

\begin{equation}\label{eq:darwin_displacement}
\vc{\mathcal{D}}^\darwin \defeq \sum_s \vc{\mathcal{P}}_{1,s}^\darwin
% = -\sum_s \frac{m_s}{B_0^2} \int \gyfzero \nabla_\perp \phi_1\jacobiangcslab \diff^3 \gy{u}
= -\mathcal{C}(1) \nabla_\perp \phi_1
% , \quad
% \vc{\mathcal{H}} = \frac{1}{\mu_0} (\vc{B}_0 + \vc{B}_1) - \sum_s \vc{\mathcal{M}}_{1,s}.
\end{equation}

Since the Darwin displacement field explicitly depends on the potential $\phi_1$ it follows that the Darwin model cannot directly be expressed in terms of the electromagnetic fields.

Regarding well-posedness of the resulting field equations it is important to note that the displacement field in Ampère's law (<ref>) (resulting in the displacement current) does not coincide with the displacement field from the quasi-neutrality equation (<ref>).
However, for compatibility of the field equations, which is a necessary condition for well-posedness, we have found in <ref> that the bound continuity equation must be satisfied.
For the Darwin model the bound continuity equation is given by (cf. <ref>)

\begin{equation}
\pard{\mathcal{R}^{\bound}}{t} + \nabla \bcdot (\vc{\mathcal{J}}^{\bound,\darwin} + \mathcal{C}(1) \nabla_\perp \lambda) = 0
\quad \iff \quad
\pardd{}{t} \nabla \bcdot (\mathcal{C}(1) \vc{A}_{1,\perp}) = - \nabla \bcdot (\mathcal{C}(1) \nabla_\perp \lambda),
\end{equation}

where the bound current density is given by $\vc{\mathcal{J}}^{\bound,\darwin} \defeq  {\partial \vc{\mathcal{D}}^\darwin} / {\partial t} + \sum_s \curl \vc{\mathcal{M}}_{1,s}$, and we have substituted <ref>.
This shows that the constraint $\nabla \bcdot (\mathcal{C}(1) \vc{A}_{1,\perp}) = 0$, which in weak form is given by <ref>, is sufficient to ensure that the field equations are compatible, provided that we can show that the Lagrange multiplier $\lambda$ vanishes.

When substituting the constraint (<ref>) into the quasi-neutrality equation (<ref>) we find

\begin{equation}
- \nabla \bcdot (\mathcal{C}(1) \nabla_\perp \phid) = \evalperturbedgcadjoint{\mathcal{R}}{}^\free
\end{equation}

which is decoupled from Ampère's law and thereby reduces to the well-known (and well-posed) perpendicular Laplace equation for $\phi_1$.

For the vector potential $\vc{A}_1$ and Lagrange multiplier $\lambda$ we have to solve the following system of equations

\begin{align}
\curl \squarepar{\frac{1}{\mu_0} \curl \Ad + \frac{p_{0,\perp} - p_{0,\shortparallel}}{B_0^2} (\curl \Ad)_\perp} - \mathcal{C}(1) \nabla_\perp \lambda = {} & \evalperturbedgcadjointalt{\vc{\mathcal{J}}}{}^\free - \mathcal{C}(1) \pard{}{t} \nabla_\perp \phid \nonumber\\ & - \frac{1}{\mu_0} \curl \vc{B}_0,\label{eq:darwin_ampere_saddle}\\
\nabla \bcdot \roundpar{\mathcal{C}(1) \Adperp} = & 0,\label{eq:darwin_constraint_saddle}
\end{align}

where we have substituted the definition of the magnetizing field <ref> and magnetization <ref> (where the $\hb \times \vc{E}_1$ term is no longer present due to the removal of the $\hb \bcdot (\vc{E}_{1,\perp} \times \vc{B}_1)$ term from the second-order gyrocentre Hamiltonian) into <ref> and the known terms are moved to the right-hand side.
The key property of <ref> is that it is invariant under $\vc{A}_1 \mapsto \vc{A}_1 + \nabla \eta$ and a gauge condition is therefore needed to fix this freedom.
This gauge condition is exactly provided by the constraint (<ref>).
The well-posedness of <ref> for an isotropic pressure is discussed in <ref>.

Note that including the $\hb \bcdot (\vc{E}_{1,\perp} \times \vc{B}_1)$ term in the second-order gyrocentre Hamiltonian would have resulted in an overdetermined system of equations for $\Ad$, and would therefore have required a non-zero $\lambda$.
It is exactly for this reason that we have omitted this term from $\gy{H}_{2,s}$.

§.§.§ Energy conservation

\begin{multline}
\pard{}{t} (\vc{\mathcal{D}}^\darwin \bcdot \vc{E}) =
2\pard{\vc{\mathcal{D}}^\darwin}{t} \bcdot \vc{E}
+ \sum_s \frac{m_s}{B_0^2} \int \gyfzero \gy{u}_\shortparallel \squarepar{(\hb \times \vc{B}_1) \bcdot \pard{\vc{E}}{t} - \roundpar{\hb \times \pard{\vc{B}_1}{t}} \bcdot \vc{E}} \jacobiangcslab \diff^3 \gy{u}\\
+ \nabla \bcdot \squarepar{\sum_s \frac{m_s}{B_0^2} \int \gyfzero \jacobiangcslab \diff^3 \gy{u} \roundpar{\pardd{\Adperp}{t} \phi_1 - \pard{\Adperp}{t} \pard{\phi_1}{t}}}
\end{multline}

For the Darwin model we find that evolution equation for the kinetic energy density (<ref>) is unchanged.
The equivalent to Poynting's theorem can also be shown, except that the potential energy density now evolves according to (cf. <ref>)

\begin{equation}
\pard{\mathcal{U}^\darwin}{t} + \nabla \bcdot \squarepar{\vc{E} \times \vc{\mathcal{H}} - \frac{\mathcal{C}(1)}{2} \roundpar{\pardd{\Adperp}{t} \phi_1 - \pard{\Adperp}{t} \pard{\phi_1}{t}}} = - \evalperturbedgcadjointalt{\vc{\mathcal{J}}}{}^\free \bcdot \vc{E}.
\end{equation}

This results in the following conserved energy

\begin{equation}
\mathfrak{E}^\darwin \defeq \int (\mathcal{T} + \mathcal{U}^\darwin) \diff^3 \gy{r},
\end{equation}

where the kinetic energy is still defined by <ref> and the potential energy density is given by

\begin{equation}\label{eq:dar_energy_pot}
\mathcal{U}^\darwin \defeq
\frac{1}{2} (
\vc{\mathcal{D}}^\darwin \bcdot \vc{E}
\vc{\mathcal{H}} \bcdot \vc{B}
% + \frac{\mathcal{C}(1)}{2} \pard{\Adperp}{t} \bcdot \nabla_\perp \phi_1
\end{equation}

§.§ Linearized models in a slab

In order to study the dispersive properties of the models under consideration, we consider the case of slab geometry, wherein the background magnetic field $\vc{B}_0$ is assumed to be constant.
In this case the FLR corrected electromagnetic fields exactly coincide with what one would expect (cf. <ref>)

\begin{equation}
\vc{\mathcal{B}}_1 = \gav{\evalperturbedgc{\vc{B}}_1}
, \quad
\vc{\mathcal{E}}_1 = \gav{\evalperturbedgc{\vc{E}}_1}
\end{equation}

since $\varepsilon_B = 0$.
Furthermore, the term that multiplies $\dot{\gy{\vc{R}}}$ in Ampère's law now exactly coincides with $\gav{\evalperturbedgc{\ca}}$ (cf. <ref>).
We reiterate that each of these three identities hold only because of the specific choice of our gyrocentre coordinate transformation $(\gparamr, \gparamtheta) = (1, 0)$.

§.§.§ Proposed gauge invariant model

We make use of the so-called susceptibility tensor in order to study the dispersive properties of our model.
The susceptibility tensor represents the linearized model with a Fourier ansatz.
More precisely, we find that computing the time derivative of Ampère's law (<ref>) results in

\begin{equation}\label{eq:ampere_E}
-\frac{1}{\mu_0} \curl (\curl \vc{E}_1)
\pardd{}{t} \roundpar{\epsilon_0 \vc{E}_1 + \sum_s \vc{\mathcal{P}}_{1,s}} + \sum_s \curl \pard{\vc{\mathcal{M}}_{1,s}}{t} + \pard{\evalperturbedgcadjointalt{\vc{\mathcal{J}}}{}^\free}{t}
\end{equation}

where we have substituted <ref>, and note that we include a possibly non-zero vacuum permittivity $\epsilon_0$ as this is customary when discussing the susceptibility tensor.
By substituting the expressions for the polarization (<ref>), magnetization (<ref>), and the gyrocentre current, and by repeatedly using Faraday's law (<ref>) we can obtain an equation which is expressed entirely in terms of the perturbed electric field $\vc{E}_1$.
We then linearize this equation and substitute the Fourier ansatz

\begin{equation}\label{eq:reduced_fourierE}
\vc{E}_1 = \fourier{\vc{E}}_1 \euler^{\ii (\vc{k} \bcdot \gy{\vc{r}} - \omega t)},
\end{equation}

which results in an equation of the form <cit.>

\begin{equation}\label{eq:wave_dispersionE}
\omega^2 (\mat{I}_3 + \gyFIXED{\mat{X}}) \fourier{\vc{E}}_1 + c^2 \vc{k} \times (\vc{k} \times \fourier{\vc{E}}_1) = \vc{0}_3,
\end{equation}

where $\gyFIXED{\mat{X}}$ is referred to as the gyrokinetic susceptibility tensor, $\mat{I}_3$ denotes the $3 \times 3$ identity matrix and $c = \sqrt{1 / (\mu_0 \epsilon_0)}$ denotes the speed of light in vacuum.
Note that $\mat{I}_3 + \gyFIXED{\mat{X}}$ is often referred to as the dielectric tensor.
The susceptibility tensor contains the (linearized) contributions to Ampère's law from the polarization, magnetization as well as the gyrocentre current and reduces the linearized gyrokinetic model to a material property: the permittivity of the `material' is given by $\epsilon_0 (\mat{I}_3 + \gyFIXED{\mat{X}})$.

We follow the discussion from [Zonta et al., 2021], wherein the gyrokinetic susceptibility tensor $\gyFIXED{\mat{X}}^\zlr$ is derived for the drift kinetic model from [Burby & Brizard, 2019] (i.e. our proposed model with $(\gparamr, \gparamtheta) = (0, 0)$ in the ZLR limit $\rho \rightarrow 0$) and subsequently compared to the ZLR limit <cit.> of the Vlasov–Maxwell susceptibility tensor <cit.>.

The derivation of the susceptibility tensor can be found in <ref>, and the expression for the susceptibility tensor in its most general form can be found in <ref>.
We now make two simplifications: firstly, we consider the ZLR limit of the susceptibility tensor, resulting in the drift kinetic susceptibility tensor given by <ref>, which coincides exactly with the low-frequency and ZLR limit of the Vlasov–Maxwell susceptibility tensor found in <cit.>.
And secondly, we consider the use of a centred Maxwellian background particle distribution as defined in <ref>.
Note that a constant background magnetic field $\vc{B}_0$ combined with a centred Maxwellian background distribution with constant density trivially satisfies the equilibrium conditions discussed in <ref>.
When using the identities given by <ref> we find that the resulting susceptibility tensor is given by

\begin{equation}\label{eq:gy_susceptibility_zlr_cm}
\begin{aligned}
\frac{\omega^2}{c^2} \gyFIXED{\mat{X}}= {} &
\sum_s \frac{\beta_{0,s}}{\uths^2}
\begin{pmatrix}
\omega^2 & \ii \omega \omega_{\mathrm{c},s} & 0\\
-\ii \omega \omega_{\mathrm{c},s} & \omega^2 - \kperp^2 \uths^2 & 0\\
0 & 0 & 0\\
\end{pmatrix}\\
& - \sum_s \frac{\beta_{0,s}}{\gy{\mathcal{N}}_{00,s} \uths^2} \int \frac{
\kpa^2 \gyfzerocm
{(\kpa \gy{u}_\shortparallel - \omega)^2}
\begin{pmatrix}
0 & 0 & 0\\
0 & \gy{u}_\tau^4 \kperp^2 \frac{1}{4} & \ii \frac{\kperp}{\kpa} \gy{u}_\tau^2 \frac{\omega \omega_{\mathrm{c},s}}{2}\\
0 & -\ii \gy{u}_\tau^2 \frac{\kperp}{\kpa} \frac{\omega \omega_{\mathrm{c},s}}{2} & \frac{\omega^2 \omega_{\mathrm{c},s}^2}{\kpa^2}
\end{pmatrix} \jacobiangcslab \diff^3 \gy{u},
\end{aligned}
\end{equation}

where we recall that the plasma-$\beta$ is defined in <ref>.

The remaining integrals can be expressed in terms of the plasma dispersion function $\mathcal{Z}$

\begin{equation}\label{eq:plasma_disp_fun}
\mathcal{Z}(\zeta) = \frac{1}{\sqrt{\pi}} \int_{-\infty}^{\infty} \frac{\euler^{-u^2}}{u - \zeta} \diff u.
\end{equation}

Substitution of <ref> in <ref> results in

\begin{equation}\label{eq:wave_dispersionE_mat}
\gyFIXED{\mat{D}}\fourier{\vc{E}}_1 = \vc{0}_3
, \quad
\gyFIXED{\mat{D}}
\frac{\omega^2}{c^2} \mat{I}_3
+ \frac{\omega^2}{c^2} \gyFIXED{\mat{X}}
+ \mat{K}^2,
\end{equation}

where the matrix $\mat{K}$ is such that

\begin{equation}\label{eq:Kmatrix}
% \mat{I}_\perp \vc{S} = (\mat{I}_3 - \hb \hb^\transpose) \vc{S} = \vc{S}_\perp
% , \quad
\mat{K} \vc{S} = \vc{S} \times \vc{k},
\end{equation}

and the Hermitian matrix $\gyFIXED{\mat{D}}$ can explicitly be written as

Note that

\begin{equation}
\mathcal{Z}'(\zeta)
= \frac{1}{\sqrt{\pi}} \int_{-\infty}^{\infty} \frac{\euler^{-u^2}}{(u - \zeta)^2} \diff u
% = -2 (1 + \zeta \mathcal{Z}(\zeta))
\end{equation}

such that

\begin{equation}
\int \frac{\kpa^2 \gy{\mathcal{F}}_{0,s}^\cm}{(\kpa \gy{u}_\shortparallel - \omega)^2} \gy{u}_\tau^{2p} \jacobiangcslab \diff^3 \gy{u}
\mathfrak{I}_\shortparallel
\mathfrak{I}_\mu(p)
\end{equation}

where the integrals are defined as

\begin{equation}
\mathfrak{I}_\shortparallel
\frac{2}{\uths^4}
\mathcal{Z}'(\omega / (\kpa \uths))
, \quad
\mathfrak{I}_\mu(p)
\int_0^\infty
\exp{\left(-\frac{2 B_0 \gy{\mu}}{m_s \uths^2} \right)}
\roundpar{\frac{2 B_0 \gy{\mu}}{m_s}}^p \frac{B_0}{m_s}
\diff \gy{\mu}
\end{equation}

for which

\begin{equation}
\mathfrak{I}_\mu(0) = \frac{\uths^2}{2}
, \quad
\mathfrak{I}_\mu(1) = \frac{\uths^4}{2}
, \quad
\mathfrak{I}_\mu(2) = \uths^6
\end{equation}

This results in

\begin{multline}
\omega^2 \gyFIXED{\mat{X}}=
\sum_s \beta_{0,s} \frac{c^2}{\uths^2}
\begin{pmatrix}
\omega^2 & \ii \omega \omega_{\mathrm{c},s} & 0\\
-\ii \omega \omega_{\mathrm{c},s} & \omega^2 - \kperp^2 \uths^2 & 0\\
0 & 0 & 0\\
\end{pmatrix}\\
- \sum_s \beta_{0,s} \frac{c^2}{\uths^2} \mathcal{Z}'({\omega} / (\kpa \uths))
\begin{pmatrix}
0 & 0 & 0\\
0 & \frac{\kperp^2 \uths^2}{2} & \ii \frac{\kperp \omega \omega_{\mathrm{c},s}}{2\kpa}\\
0 & -\ii \frac{\kperp \omega \omega_{\mathrm{c},s}}{2\kpa} & \frac{\omega^2 \omega_{\mathrm{c},s}^2}{\kpa^2 \uths^2}
\end{pmatrix}.
\end{multline}

\begin{equation}\label{eq:gy_dispersionmat_zlr_cm}
\begin{aligned}
\gyFIXED{\mat{D}}
= {} &
\begin{pmatrix}
\frac{\omega^2}{c^2} - \kpa^2 & 0 & \kpa \kperp\\
0 & \frac{\omega^2}{c^2} - k^2 & 0\\
\kpa \kperp & 0 & \frac{\omega^2}{c^2} - \kperp^2
\end{pmatrix}
+ \sum_s \frac{\beta_{0,s}}{\uths^2}
\begin{pmatrix}
\omega^2 & \ii \omega \omega_{\mathrm{c},s} & 0\\
-\ii \omega \omega_{\mathrm{c},s} & \boxed{\omega^2} - \kperp^2 \uths^2 & 0\\
0 & 0 & 0\\
\end{pmatrix}\\
& - \sum_s \frac{\beta_{0,s}}{\uths^2} \mathcal{Z}'({\omega} / (\kpa \uths))
\begin{pmatrix}
0 & 0 & 0\\
0 & \frac{\kperp^2 \uths^2}{2} & \ii \frac{\kperp \omega \omega_{\mathrm{c},s}}{2\kpa}\\
0 & -\ii \frac{\kperp \omega \omega_{\mathrm{c},s}}{2\kpa} & \frac{\omega^2 \omega_{\mathrm{c},s}^2}{\kpa^2 \uths^2}
\end{pmatrix}.
\end{aligned}
\end{equation}

We highlight the boxed ${\omega^2}$ term for the discussion that will follow in <ref>.

Non-trivial solutions to <ref> exist if and only if the dispersion matrix given by <ref> is singular and to this end we implicitly define the dispersion relation $\omega(\vc{k})$ via

\begin{equation}
\det \gyFIXED{\mat{D}}(\omega(\vc{k}), \vc{k}) = 0.
\end{equation}

We nondimensionalize the dispersion matrix using

\begin{equation}
\massratios = \frac{m_s}{m_\ii}
, \quad
\chargeratios = \frac{q_s}{q_\ii}
, \quad
\nondim{\omega} = \frac{\omega}{\uthi \kpa}
, \quad
\nondim{\vc{k}} = \varrho_\ii \vc{k}
, \quad
\nondim{c} = \frac{c}{\uthi}
\end{equation}

and consequently

\begin{equation}
\uths = \frac{\uthi}{\sqrt{\massratios}}
, \quad
\beta_{0,s} = \frac{m_s \mathcal{N}_{0} \mu_0 \uths^2}{B_0^2} = \beta_\ii
, \quad
\omega_{\mathrm{c},s} = \frac{q_s B_0}{m_s} = \frac{\chargeratios}{\massratios} \omega_{\mathrm{c},\ii}
, \quad
\varrho_s
= \frac{m_s \uths}{q_s B_0}
= \frac{\sqrt{\massratios}}{\chargeratios}  \varrho_\ii
\end{equation}

while assuming $\mathcal{N}_{s} = \mathcal{N}_{0}$.
This results in

\begin{multline}
\nondim{\gyFIXED{\mat{D}}}
\begin{pmatrix}
\frac{\nondim{\omega}^2}{\nondim{c}^2} \nondim{k}_\shortparallel^2 - \nondim{k}_\shortparallel^2 & 0 & \nondim{k}_\shortparallel \nondim{k}_\perp\\
0 & \frac{\nondim{\omega}^2}{\nondim{c}^2} \nondim{k}_\shortparallel^2 - \nondim{k}^2 & 0\\
\nondim{k}_\shortparallel \nondim{k}_\perp & 0 & \frac{\nondim{\omega}^2}{\nondim{c}^2} \nondim{k}_\shortparallel^2 - \nondim{k}_\perp^2
\end{pmatrix}
+ \beta_\ii \sum_s
\begin{pmatrix}
\massratios \nondim{k}_\shortparallel^2 \nondim{\omega}^2 & \ii \nondim{k}_\shortparallel \nondim{\omega} \chargeratios & 0\\
-\ii \nondim{k}_\shortparallel \nondim{\omega} \chargeratios & \massratios \nondim{k}_\shortparallel^2 \nondim{\omega}^2 - \nondim{k}_\perp^2 & 0\\
0 & 0 & 0\\
\end{pmatrix}\\
- \frac{\beta_\ii}{2} \sum_s \mathcal{Z}'(\sqrt{\massratios} \nondim\omega)
\begin{pmatrix}
0 & 0 & 0\\
0 & \nondim{k}_\perp^2 & \ii \nondim{k}_\perp \nondim{\omega} \chargeratios\\
0 & -\ii \nondim{k}_\perp \nondim{\omega} \chargeratios & 2\nondim{\omega}^2 \chargeratios^2
\end{pmatrix},
\end{multline}

after having multiplied out the common $\varrho_\ii^2$ factor.

§.§.§ Other models

It is a priori not clear what effect the Darwin approximation has on the dispersive properties of our model and to this end we also consider the gyrokinetic Darwin susceptibility tensor

for which a derivation can be found in <ref>.
When using a centred Maxwellian, as we did in <ref>, we find that the susceptibility tensor is given by <ref> where the boxed $\omega^2$ term is removed.

As discussed in <ref>, the EOMs of the symplectic Brizard–Hahm (BH) model [Brizard & Hahm, 2007] coincide with the ones from our proposed gauge invariant model if $\varepsilon_B = 0$, which is what we assume here.
Similarly, we find that the magnetization term vanishes, as it does in our proposed model, for the centred Maxwellian background particle distribution function that we consider here.
Hence, the only difference between the proposed model and the Brizard–Hahm model, under the current simplifying assumptions, is the missing polarization current.

However, when reconsidering <ref> under the current simplifying assumptions, we find that the compatibility constraint reduces to

\begin{equation}
\nabla \bcdot \roundpar{\mathcal{C}(1) \nabla_\perp \lambda}
= \pard{}{t} \nabla \bcdot (\mathcal{C}(1) \nabla_\perp \phi_1)
\end{equation}

upon substitution of <ref>, and it therefore follows that in this case $\lambda = \partial \phi_1 / \partial t$.
This implies that the Lagrange multiplier for the Brizard–Hahm model restores the polarization current and thereby coincides with the Darwin approximation of our model.
Note that this only works under the simplifying assumptions that we use here: the model is linearized, we assume $\varepsilon_B = 0$, use a centred Maxwellian background distribution, and we consider the ZLR limit.

For the parallel-only model we make use of the dispersion relation provided in <cit.>.

§.§.§ Comparison of the dispersion relations

The models are now compared by numerically evaluating the dispersion relations (the code can be found in <cit.>).
We consider electron and ion (deuteron) species, with $T_\ii = T_\euler$, $q_\ii = - q_\euler$ and an electron to ion mass ratio given by $m_\euler / m_\ii = 1 / 3671$.
Furthermore, we nondimensionalize the frequency and wave vector as

\begin{equation}
\nondim{\omega} \defeq \frac{\omega}{\uthi \kpa}
, \quad
\nondim{\vc{k}} \defeq \varrho_\ii \vc{k}
\end{equation}

In this case, we expect to find a shear and compressional Alfvén wave, which have the following frequencies (in the quasi-neutral limit $c \rightarrow \infty$)

\begin{equation}\label{eq:aflven_frequencies}
\nondim{\omega}_\mathrm{As} \defeq \frac{1}{\sqrt{{\beta_{0}}}}
, \quad
\nondim{\omega}_\mathrm{Ac} \defeq \frac{k}{\kpa} \nondim{\omega}_\mathrm{As}
\end{equation}

We note that the shear Alfvén frequency is constant w.r.t. $\nondim{k}_\perp$ whereas the compressional Alfvén frequency increases linearly with increasing wavenumber.
Hence, the compressional Alfvén wave is a fast wave.
This is especially true when small perpendicular length scales are considered (turbulence), since then the compressional Alfvén frequency becomes comparable to the cyclotron frequency

\begin{equation}
\varepsilon_\perp \sim 1, \quad \omega = \omega_{\mathrm{Ac}}
\quad \implies \quad
\varepsilon_\omega = \frac{\omega_{\mathrm{Ac}}}{\omega_{\mathrm{c},\ii}} = \frac{\nondim{k}}{\sqrt{\beta_{0}}} \sim 1.
\end{equation}

The presence of compressional Alfvén waves in our gyrokinetic model is therefore incompatible with the low-frequency approximation of the first-order generating function (<ref>), which relies on $\varepsilon_\omega \ll 1$.
This incompatibility seems to suggest that we have made a mistake in the derivation of our gyrokinetic model.
The origin of this issue lies in the low-frequency approximation of the first-order generating function as discussed in <ref>, where we have intentionally kept the gauge invariant parts of the right-hand side of <ref>.
While this ensures gauge invariance, it neglects the fact that the part of the electric field that comes from the vector potential has a time derivative, and this term would therefore be neglected if we had neglected all $\mathcal{O}(\varepsilon_\omega)$ terms.

The real part.

The negative imaginary part.

The dispersion relations for a fixed value of $\nondim{k}_\shortparallel = 2 \times 10^{-3}$ and $\beta_{0} = 10 \%$.
The black dotted line corresponds to $\nondim\omega = \nondim\omega_\mathrm{As}$, whereas the black dashed line corresponds to $\nondim\omega = \nondim\omega_\mathrm{Ac}$.

In <ref> we show the dispersion relations for a fixed value of $\nondim{k}_\shortparallel = 2 \times 10^{-3}$ and $\beta_{0} = 10 \%$.
The gauge invariant model with $(\gparamr, \gparamtheta) = (1, 0)$ results in two waves, which have the correct real part of the frequency, whereas the imaginary part has the correct sign which results in damping.
The dispersion relation of the other models result only in the shear Alfvén wave, as expected.

$R_0$     $a$     $\varrho_\mathrm{i}$     $\nondim{k}_\shortparallel$     $\nondim{k}_\perp$

AUG     1.6 m     0.8 m     $3.4 \times 10^{-3}$ m     $2.1\times 10^{-3}$     $7.2\times 10^{-3}$

ITER     6.2 m     2 m     $2.4 \times 10^{-3}$ m     $3.9\times 10^{-4}$     $2.0 \times 10^{-3}$

W7-X     5.5 m     0.53 m     $4.4 \times 10^{-3}$ m     $7.9\times 10^{-4}$     $1.4 \times 10^{-2}$

The length scales used for determining the wave vector $\nondim{\vc{k}}$, as obtained from [Zoni & Possanner, 2021].
The nondimensional wavenumbers $\nondim{k}_\shortparallel, \nondim{k}_\perp$ are computed according to <ref>.

In order to further explore the parameter space, we now use the geometrical parameters of the tokamak fusion devices ASDEX Upgrade (AUG) and ITER as well as the stellarator Wendelstein 7-X (W7-X), to determine the value of the wave vector $\nondim{\vc k}$.
For both the parallel and perpendicular direction we consider the lowest non-trivial wavenumber

\begin{equation}\label{eqn:alfven:eval:wavenumber_comp}
\nondim{k}_\shortparallel = \frac{\varrho_\ii}{R_0}, \quad \nondim{k}_\perp = \frac{2\pi \varrho_\ii}{a\sqrt{4 + \pi^2}} ,
\end{equation}

where $R_0$ and $a$ denote the major and minor radii of the tokamak, respectively.
The values of $R_0, a, \varrho_\ii$ are shown in <ref> and are taken from [Zoni & Possanner, 2021] for the tokamak fusion devices and from [Grieger et al., 1992, Klinger et al., 2019] for the stellarator W7-X.
We moreover show the resulting wavenumbers according to <ref>.

The real part.

The negative imaginary part.

The dispersion relations for fixed values of $\nondim{k}_\perp$ and $\nondim{k}_\shortparallel$, as determined from <ref>.
Only the shear Alfvén wave is shown.

For each of the three models and each of the three machines, we compute the dispersion curve where we now vary $\beta_{0}$ and keep the wave vectors fixed.
The results are shown in <ref>.
We find that the real part of the dispersion curves all overlap and agrees with <ref>, which suggests that the shear Alfvén frequency depends only on $\beta_{0}$.
When considering the imaginary part of the dispersion curve, we find differences not only between the machines, but also between the models, in particular for larger values of $\beta_{0}$.
The gauge invariant and Darwin model agree well, which shows that removal of the compressional Alfvén wave has not altered the shear Alfvén wave.
Moreover, we find that the parallel-only model yields a different imaginary part of the shear Alfvén frequency.

§.§ Summary of comparison

$\gyFIXED{\vc{\gamma}}_{1,\vc{R}}$     Reduces to
parallel-only     Gauge
invariant     Compressional
Alfvén wave     $\lambda = 0$

Parallel-only     $q \gav{\evalperturbedgc{A}_{1,\shortparallel}} \hb$     N/A     N     N     N/A

Brizard–Hahm     $q \gav{\evalperturbedgc{\vc{A}}_1}$     Y     N     $^\dagger$N$^\dagger$     N

Burby–Brizard     $q \vc{A}_1$     N     Y     N     N

Darwin approx.     $q \Ad + q \radgav{\evalperturbedgcburby{\vc{B}}_1 \times \vc{\rho}}$     Y     N     N     Y

Proposed     $q \Ad + q \radgav{\evalperturbedgcburby{\vc{B}}_1 \times \vc{\rho}}$     Y     Y     Y     Y

Properties of the different gyrokinetic models under consideration.
The models proposed in this paper are in boldface: the gauge invariant model with $(\gparamr, \gparamtheta) = (1, 0)$ and its corresponding Darwin approximation.
$^\dagger$ This can be a `yes (Y)' if the approach from [Qin et al., 1999] is followed.

We have compared five different gyrokinetic models, all of which are derived from an action principle.
An overview of the most essential properties is given in <ref>, wherein we have also included the symplectic part of the gyrocentre single-particle phase-space Lagrangian $\gyFIXED{\vc{\gamma}}_{1,\vc{R}}$ – with all other components equal to zero – which defines the gyrocentre coordinate transformation.
We now summarize the comparison of the models.

The parallel-only model [Kleiber et al., 2016], as described in <ref>, neglects the perpendicular part of the vector potential and thereby results in two scalar well-posed perpendicular Poisson-type field equations which decouple for a centred Maxwellian background particle distribution function.
The model is not gauge invariant.
Due to the absence of the perpendicular part of the vector potential, we find that the EOMs of the parallel model do not contain a grad-$B_{1,\shortparallel}$ drift.
This model is still the `working horse' and is widely used for the global simulation of turbulence in fusion devices [Garbet et al., 2010, Mishchenko et al., 2023].

The symplectic Brizard–Hahm model [Brizard & Hahm, 2007] (see <ref>) contains the full vector potential and results in a curl-curl type Ampère's law which is coupled to the same quasi-neutrality equation as found in the parallel-only model.
This model is also not gauge invariant and reduces to the parallel-only model when neglecting the perpendicular part of the vector potential.
The grad-$B_{1,\shortparallel}$ drift is present in the EOMs; this is the case for each of the gyrokinetic models which have the full vector potential.
Interestingly, the second-order Hamiltonian of the Brizard–Hahm model coincides with that of the parallel-only model.
This implies that this model has no polarization current, and the magnetization effects do not depend on the perpendicular part of the vector potential.
Moreover, the magnetization current is not the curl of a magnetization, and we find that the continuity equation is not satisfied for the gyro-average adjoint of the free charge and current density.
For those reasons the total continuity equation, consisting of the free and bound charge, is not satisfied, which is a necessary condition for well-posedness of the field equations.
It follows that well-posedness of the Brizard–Hahm model [Brizard & Hahm, 2007] requires the introduction of a Lagrange multiplier, which affects both the EOMs and the field equations in an unphysical way.

The gyrokinetic model from [Burby & Brizard, 2019] has played an essential role in this paper, as it has guided us to the development of a family of gauge invariant gyrokinetic models.
The model from [Burby & Brizard, 2019] coincides with $(\gparamr, \gparamtheta) = (0, 0)$, except that the authors propose to neglect the polarization current altogether, as it is a higher-order contribution in $\varepsilon_\omega$.
This decision however does imply that a Lagrange multiplier will be needed for this model to make it well-posed.
We have not discussed this model in any detail beyond the computation of the first- and second-order Hamiltonian, as found in <ref>.
The reason that we have not considered this parameter combination any further is because $\gparamr = 0$ results in a model where the symplectic part of the model is essentially drift-kinetic, does not reduce to the parallel-only model if the perpendicular part of the vector potential is neglected, and the Hamiltonian part of the model contains terms which close themselves off a physical interpretation.

Instead, we propose to use $(\gparamr, \gparamtheta) = (1, 0)$, which results in a gauge invariant gyrokinetic model for which all terms in both the EOMs and the field equations can be interpreted from a physical point of view.
Moreover, it yields the smallest coordinate transformation while resulting in a gauge invariant model wherein the gyrocentre magnetic moment is an invariant.
This model does reduce to the parallel-only model if the perpendicular part of the vector potential is neglected.
The single-particle Lagrangian and resulting Vlasov–Maxwell action principle of this model are derived and described in detail in <ref>, respectively.
We find that the continuity equation is satisfied for the gyro-average adjoint of the free charge and current density, and, since the model is gauge invariant and derived from an action principle, it follows that the bound charge continuity equation is satisfied as well.
These are necessary conditions for well-posedness of the field equations, which are satisfied without the need of a Lagrange multiplier.

The presence of the polarization current in our proposed gauge invariant model results in a compressional Alfvén wave, which is often undesired due to its relatively high frequency.
To this end, we have proposed a Darwin approximation to our gauge invariant gyrokinetic model, which is discussed in <ref>.
Therein, we deliberately remove the part of the polarization current which is responsible for the compressional Alfvén wave, but we thereby break gauge invariance of the model.
However, we have shown that the field equations are well-posed and decouple without the need of a Lagrange multiplier.

Finally, we have derived a dispersion relation for each of the models, which has confirmed the expected properties of the models: each model has the shear Alfvén wave and only the proposed gauge invariant model includes the compressional Alfvén wave.
In both cases the real part of the frequency agrees well with the theory, whereas the imaginary part is negative and thereby results in damping of the wave.

§ CONCLUSIONS

Motivated by the need for a more complete gyrokinetic model, wherein the perpendicular part of the vector potential is kept, we have discussed the gyrocentre coordinate transformation in detail.
The purpose of the gyrocentre coordinate transformation is to transform the perturbed single-particle guiding-centre phase-space Lagrange in such a way that it becomes independent of the gyro-phase.
This results in a reduction of phase-space dimension from six to five (the gyrocentre position, the parallel velocity and the invariant magnetic moment), and when moreover considering the limit of quasi-neutrality, removes the fastest time scales from the large range of length and time scales present in the Vlasov–Maxwell model.

When gyrokinetic modelling is considered at the level of the Lagrangian, which thereby utilizes a variational principle, it is ensured that the gyrokinetic model is structure preserving.
In particular, we find that energy is conserved regardless of the modelling error introduced by the truncated gyrocentre coordinate transformation.
However, an aspect that is often overlooked in gyrokinetics is the property of gauge invariance [Brizard & Hahm, 2007], which ensures that the model is invariant under the gauge transformation, as is the case of the Vlasov–Maxwell model.
In particular, the traditionally used parallel-only models are not gauge invariant.
To this end, we have generalized the approach proposed in [Burby & Brizard, 2019], wherein a gauge invariant gyrokinetic model is introduced.
We have derived sufficient conditions on the gyrocentre coordinate transformation which ensure that the resulting single-particle gyrocentre phase-space Lagrangian is gauge invariant.
Despite this additional restriction, the gyrocentre coordinate transformation is by no means uniquely defined, and this approach therefore results in a family of gauge invariant gyrokinetic models.
The family of models is parametrized by the two parameters $\gparamr, \gparamtheta$, where the model from [Burby & Brizard, 2019] coincides with $(\gparamr, \gparamtheta) = (0, 0)$.
Our derivation is presented without making use of the formalism of differential geometry.

In an effort to obtain a gyrokinetic model for which each of the equations of motion as well as the field equations can be interpreted from a physical point of view, we have chosen $(\gparamr, \gparamtheta) = (1, 0)$.
This choice leads to the smallest coordinate transformation which results in a gauge invariant model wherein the gyrocentre magnetic moment is an invariant.

We find that our model reduces to the parallel-only model when the perpendicular component of the vector potential is neglected, which does not hold for the gauge invariant model from [Burby & Brizard, 2019].
The resulting model has been derived to second-order in the perturbation parameter $\varepsilon_\delta$, and the second-order part of the Lagrangian contains polarization and magnetization effects which have a clear physical meaning.
Due to gauge invariance the model can be expressed directly in terms of the electromagnetic fields, rather than the potentials.
The gyrokinetic model thereby results in the macroscopic Maxwell's equations.
Moreover, we have derived equilibrium conditions on the background distribution function and magnetic field which justify the smallness of the perturbation parameter $\varepsilon_\delta$.

We find that our gauge invariant model possesses a magnetization current which is the curl of a magnetization.
In addition, it has a polarization current, and we find that our free charge continuity equation holds for the gyro-average adjoints of the free charge and current.
Each of those three properties is essential for showing that the bound charge continuity equation is naturally satisfied, which is necessary for well-posedness of the field equations.
Hence, unlike the Brizard–Hahm model [Brizard & Hahm, 2007], we find that the field equations of our gauge invariant model are compatible without the need of a Lagrange multiplier.
A brief summary of the comparison between each of the models under consideration is found in <ref>.

In addition, we have derived the gyrokinetic susceptibility tensor, which covers the material properties of the linearized gyrokinetic model in the macroscopic Maxwell's equations, for each of the models under consideration.
In the ZLR limit, we find that the gyrokinetic susceptibility tensor of our gauge invariant model agrees with that of the Vlasov–Maxwell system.
Moreover, the resulting dispersion relation shows that, besides the usual shear Alfvén wave, a fast compressional Alfvén wave is present in our gauge invariant model.
In fact, the compressional Alfvén wave is necessarily present when gauge invariance is sought after, due to the need for a gauge invariant approximation of the first-order generating function.
Due to the potentially high frequency of this wave we have proposed a Darwin approximation to our proposed gauge invariant model which successfully removes the fast compressional Alfvén wave without breaking the structure of the field equations.
We find that our Darwin model is still well-posed without the need of a Lagrange multiplier.

In future work, we plan to implement the two proposed models in the gyrokinetic particle-in-cell code EUTERPE [Jost et al., 2001, Kleiber et al., 2024] and compare them to the well-established and traditionally used parallel-only model.

§ ACKNOWLEDGMENTS

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 – EUROfusion).
Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission.
Neither the European Union nor the European Commission can be held responsible for them.

We thank
A. Bottino,
M. Campos Pinto,
R. Kleiber,
A. Könies,
P. Lauber,
O. Maj,
A. Mishchenko,
S. Possanner,
and B. Scott
for helpful discussions.

§ INVERSION OF THE LAGRANGE MATRIX

We consider the inversion of the Lagrange matrix

\begin{equation}
\mat{W} = \begin{pmatrix}
\mat{W}_{11} & \mat{W}_{12}\\
-\mat{W}_{12}^\transpose & \mat{W}_{22}
\end{pmatrix},
\end{equation}

\begin{equation}
\mat{W}_{11} = \begin{pmatrix}
q \mat{B} & -m \hb\\
m \hb^\transpose & 0
\end{pmatrix},\quad
\mat{W}_{12} = \begin{pmatrix}
\frac{m}{q} \vc{w} & \vc{0}_3\\
0 & 0\\
\end{pmatrix}, \quad
\mat{W}_{22} = \begin{pmatrix}
0 & \frac{m}{q} \\
-\frac{m}{q} & 0\\
\end{pmatrix}.
\end{equation}

Here $\mat{B}$ is such that

\begin{equation}
\mat{B} \vc{S} = \vc{S} \times \vc{B}
\end{equation}

for some magnetic field $\vc{B}$.

The inverse of $\mat{W}_{22}$ is readily given by

\begin{equation}
\mat{W}_{22}^{-1} = \begin{pmatrix}
0 & - \frac{q}{m}\\
\frac{q}{m} & 0\\
\end{pmatrix}.
\end{equation}

Using the Schur complement

\begin{equation}
\mat{S}
% = \mat{W} / \mat{W}_{22}
= \mat{W}_{11} + \mat{W}_{12} \mat{W}_{22}^{-1} \mat{W}_{12}^\transpose,
\end{equation}

we find that the inverse of the $2\times 2$ block matrix $\mat{W}$ is given by

\begin{equation}
\mat{W}^{-1} = \begin{pmatrix}
\mat{S}^{-1} & - \mat{S}^{-1} \mat{W}_{12} \mat{W}_{22}^{-1}\\
\mat{W}_{22}^{-1} \mat{W}_{12}^\transpose \mat{S}^{-1} & \mat{W}_{22}^{-1} - \mat{W}_{22}^{-1} \mat{W}_{12}^\transpose \mat{S}^{-1} \mat{W}_{12} \mat{W}_{22}^{-1}
\end{pmatrix}.
\end{equation}

The expressions for the block matrices result in

\begin{align}
\mat{W}_{12} \mat{W}_{22}^{-1} \mat{W}_{12}^\transpose
&= \begin{pmatrix}
\frac{m}{q} \vc{w} & \vc{0}_3 & \vc{0}_3 & \vc{d}\\
0 & 0 & 0 & 0\\
\end{pmatrix} \begin{pmatrix}
0 & - \frac{q}{m} & 0 & 0\\
\frac{q}{m} & 0 & 0 & 0\\
0 & 0 & 0 & 1\\
0 & 0 & -1 & 0\\
\end{pmatrix} \begin{pmatrix}
\frac{m}{q} \vc{w}^\transpose & 0\\
\vc{0}_3^\transpose & 0\\
\vc{0}_3^\transpose & 0\\
\vc{d}^\transpose & 0\\
\end{pmatrix} \nonumber\\
&= \begin{pmatrix}
\frac{m}{q} \vc{w} & \vc{0}_3 & \vc{0}_3 & \vc{d}\\
0 & 0 & 0 & 0\\
\end{pmatrix} \begin{pmatrix}
\vc{0}_3^\transpose & 0\\
\vc{w}^\transpose & 0\\
\vc{d}^\transpose & 0\\
\vc{0}_3^\transpose & 0\\
\end{pmatrix} \nonumber\\
&= \mat{0}_4
\end{align}

\begin{equation}\label{eq:matinverse_S_W11}
\mat{S} = \mat{W}_{11}
\end{equation}

for which the inverse is given by

\begin{equation}
\mat{S}^{-1} = \begin{pmatrix}
-\frac{\mat{B}_0}{q B_0 B_{\shortparallel}} & \frac{\vc{B}}{m B_{\shortparallel}}\\
-\frac{\vc{B}^\transpose}{m B_{\shortparallel}} & 0
\end{pmatrix}
\end{equation}

as can be verified by computing the product

\begin{equation}
\mat{S}\mat{S}^{-1} = \begin{pmatrix}
q \mat{B} & -m \hb\\
m \hb^\transpose & 0
\end{pmatrix} \begin{pmatrix}
-\frac{\mat{B}_0}{q B_0 B_{\shortparallel}} & \frac{\vc{B}}{m B_{\shortparallel}}\\
-\frac{\vc{B}^\transpose}{m B_{\shortparallel}} & 0
\end{pmatrix}
= \begin{pmatrix}
\frac{-\mat{B} \mat{B}_0 + \vc{B}_0 \vc{B}^\transpose}{B_0 B_{\shortparallel}} & \frac{q \mat{B}\vc{B}}{m B_{\shortparallel}} \\
-\frac{m \hb^\transpose \mat{B}_0}{q B_0 B_{\shortparallel}} & 1
\end{pmatrix}
\end{equation}

and observing that using <ref>

\begin{equation}
\mat{B} \vc{B} = \vc{B} \times \vc{B} = \vc{0}, \quad
\mat{B}_0^\transpose \hb = \vc{B}_{0} \times \hb = \vc{0}
\end{equation}

as well as

\begin{align}
\roundpar{-\mat{B} \mat{B}_0 + \vc{B}_0 \vc{B}^\transpose} \vc{S}
&= -\mat{B} (\vc{S} \times \vc{B}_0) + \vc{B}_0 (\vc{B} \bcdot \vc{S})\nonumber\\
&= -(\vc{S} \times \vc{B}_0) \times \vc{B}  + \vc{B}_0 (\vc{B} \bcdot \vc{S})\nonumber\\
&= B_{\shortparallel} B_0 \vc{S}
\end{align}

for any vector $\vc{S}$, which implies that $-\mat{B} \mat{B}_0 + \vc{B}_0 \vc{B}^\transpose = B_{\shortparallel} B_0 \mat{I}_3$.

Using these intermediate results we are ready to evaluate the blocks of the inverse of $\mat{W}$

\begin{align}
&= - \mat{S}^{-1} \mat{W}_{12} \mat{W}_{22}^{-1} \nonumber\\
&= -\begin{pmatrix}
-\frac{\mat{B}_0}{q B_0 B_{\shortparallel}} & \frac{\vc{B}}{m B_{\shortparallel}}\\
-\frac{\vc{B}^\transpose}{m B_{\shortparallel}} & 0
\end{pmatrix}
\begin{pmatrix}
\frac{m}{q} \vc{w} & \vc{0}_3 & \vc{0}_3 & \vc{d}\\
0 & 0 & 0 & 0\\
\end{pmatrix}
\begin{pmatrix}
0 & - \frac{q}{m} & 0 & 0\\
\frac{q}{m} & 0 & 0 & 0\\
0 & 0 & 0 & 1\\
0 & 0 & -1 & 0\\
\end{pmatrix} \nonumber\\
&= -\begin{pmatrix}
-\frac{\mat{B}_0}{q B_0 B_{\shortparallel}} & \frac{\vc{B}}{m B_{\shortparallel}}\\
-\frac{\vc{B}^\transpose}{m B_{\shortparallel}} & 0
\end{pmatrix}
\begin{pmatrix}
\vc{0}_3 & - \vc{w} & -\vc{d} & \vc{0}_3\\
0 & 0 & 0 & 0\\
\end{pmatrix}\nonumber\\
&= \begin{pmatrix}
\vc{0}_3 & - \frac{\vc{w} \times \hb}{q B_{\shortparallel}}  & -\frac{\vc{d} \times \hb}{q B_{\shortparallel}} & \vc{0}_3\\
0 & -\frac{\vc{B} \bcdot \vc{w}}{m B_{\shortparallel}} & -\frac{\vc{B} \bcdot \vc{d}}{m B_{\shortparallel}}  & 0\\
\end{pmatrix}
\end{align}

\begin{equation}
= - \mat{S}^{-1} \mat{W}_{12} \mat{W}_{22}^{-1}
= \begin{pmatrix}
\vc{0}_3 & - \frac{\vc{w} \times \hb}{q B_{\shortparallel}} \\
0 & -\frac{\vc{B} \bcdot \vc{w}}{m B_{\shortparallel}}\\
\end{pmatrix}
\end{equation}

\begin{align}
-\mat{W}_{22}^{-1} \mat{W}_{12}^\transpose \mat{S}^{-1} \mat{W}_{12} \mat{W}_{22}^{-1}
&= \mat{W}_{22}^{-1} \mat{W}_{12}^\transpose (\mat{W}^{-1})_{12} \nonumber\\
\\ &= \begin{pmatrix}
\vc{0}_3^\transpose & 0\\
\vc{w}^\transpose & 0\\
\vc{d}^\transpose & 0\\
\vc{0}_3^\transpose & 0\\
\end{pmatrix} \begin{pmatrix}
\vc{0}_3 & - \frac{\vc{w} \times \hb}{q B_{\shortparallel}}  & -\frac{\vc{d} \times \hb}{q B_{\shortparallel}} & \vc{0}_3\\
0 & -\frac{\vc{B} \bcdot \vc{w}}{m B_{\shortparallel}} & -\frac{\vc{B} \bcdot \vc{d}}{m B_{\shortparallel}}  & 0\\
\end{pmatrix}\nonumber\\
&= \begin{pmatrix}
0 & 0 & 0 & 0\\
0 & 0 & -\frac{1}{q B_{\shortparallel}} \hb \bcdot\vc{w} \times \vc{d}& 0\\
0 & \frac{1}{q B_{\shortparallel}} \hb \bcdot \vc{w} \times \vc{d} & 0 & 0\\
0 & 0 & 0 & 0\\
\end{pmatrix}
\end{align}

\begin{equation}
= \mat{W}_{22}^{-1} - \mat{W}_{22}^{-1} \mat{W}_{12}^\transpose \mat{S}^{-1} \mat{W}_{12} \mat{W}_{22}^{-1}
= \begin{pmatrix}
0 & - \frac{q}{m} \\
\frac{q}{m} & 0 \\
\end{pmatrix},
\end{equation}

such that the inverse of the Lagrange matrix is given by

\begin{equation}
\mat{W}^{-1} = \begin{pmatrix}
-\frac{\mat{B}_0}{q B_0 B_{\shortparallel}} & \frac{\vc{B}}{m B_{\shortparallel}} & \vc{0}_3 & - \frac{\vc{w} \times \hb}{q B_{\shortparallel}}\\
-\frac{\vc{B}^\transpose}{m B_{\shortparallel}} & 0 & 0 & -\frac{\vc{B} \bcdot \vc{w}}{m B_{\shortparallel}}\\
\vc{0}_3^\transpose & 0 & 0 & -\frac{q}{m}\\
\frac{(\vc{w} \times \hb)^\transpose }{q B_{\shortparallel}} & \frac{\vc{B} \bcdot \vc{w}}{m B_{\shortparallel}} & \frac{q}{m}& 0\\
\end{pmatrix}.
\end{equation}

§ GYROCENTRE COORDINATE TRANSFORMATION

The gyrocentre Lagrangian $\gy{L}$ is defined according to <ref>.
The transformation rules for the Hamiltonian and symplectic part follow directly by substituting <ref> into <ref> and by subsequently making use of a Taylor series expansion centred around the gyrocentre coordinate $\gy{\vc{Z}}$.
The contribution due to the generating function $\gy{S}$ is omitted here.

For the guiding-centre Hamiltonian this results in

\begin{equation}
\gc{H}(\gc{\vc{Z}})
= \gc{H}
+ \pard{\gc{H}}{\gy{\vc{Z}}} \bcdot \roundpar{- \gy{\vc{G}}_1 + \frac{1}{2} \pard{\gy{\vc{G}}_1}{\gy{\vc{Z}}} \gy{\vc{G}}_1 - \gy{\vc{G}}_2}
+ \frac{1}{2} \pardd{\gc{H}}{\gy{\vc{Z}}} \bcddot \gy{\vc{G}}_1 \otimes \gy{\vc{G}}_1
+ \mathcal{O}(\varepsilon_\delta^3)
\end{equation}

upon substitution of <ref> and using a Taylor-series expansion of $\gc{H}$ centred around the gyrocentre coordinate $\gy{\vc{Z}}$.
We use the notational convention that an absence of arguments implies evaluation at the gyrocentre coordinate, and we let $\otimes$ denote the outer product of two vectors

\begin{equation}
(\vc{a} \otimes \vc{b})_{ij} = a_i b_j,
\end{equation}

such that

\begin{equation}
\mat{M} \bcddot \vc{a} \otimes \vc{b} = \sum_{i,j = 1}^{3} \matcmp{M}_{ij} a_i b_j.
\end{equation}

This expression can be simplified to

\begin{equation}\label{eq:gc_gy_hamtransform}
\gc{H}(\gc{\vc{Z}})
= \gc{H}
- \pard{\gc{H}}{\gy{\vc{Z}}} \bcdot \roundpar{\gy{\vc{G}}_1 + \gy{\vc{G}}_2}
+ \frac{1}{2} \pard{}{\gy{\vc{Z}}}\roundpar{\pard{\gc{H}}{\gy{\vc{Z}}} \bcdot \gy{\vc{G}}_1} \bcdot \gy{\vc{G}}_1
+ \mathcal{O}(\varepsilon_\delta^3)
\end{equation}

The same approach is followed for the symplectic part of the guiding-centre Lagrangian.
Let $\gc{\Gamma}(\gc{\vc{Z}}) = \gc{\vc{\gamma}}(\gc{\vc{Z}}) \bcdot \dot{\gc{\vc{Z}}}$.
Substitution of <ref> results in

\begin{align}
\gc{\Gamma}(\gc{\vc{Z}})
= {} & \squarepar{
\gc{\vc{\gamma}}
- \pard{\gc{\vc{\gamma}}}{\gy{\vc{Z}}} \roundpar{\gy{\vc{G}} - \frac{1}{2} \pard{\gy{\vc{G}}_1}{\gy{\vc{Z}}} \gy{\vc{G}}_1}
+ \frac{1}{2} \pardd{\gc{\vc{\gamma}}}{\gy{\vc{Z}}} \bcddot \gy{\vc{G}}_1 \otimes \gy{\vc{G}}_1
} \nonumber\\ & \bcdot \squarepar{\dot{\gy{\vc{Z}}} - \totd{}{t} \roundpar{\gy{\vc{G}} - \frac{1}{2} \pard{\gy{\vc{G}}_1}{\gy{\vc{Z}}} \gy{\vc{G}}_1}}
+ \mathcal{O}(\varepsilon_\delta^3)
\end{align}

where we write $\gy{\vc{G}} \defeq \gy{\vc{G}}_1 + \gy{\vc{G}}_2$ for brevity, and we use the notational convention that an absence of arguments implies evaluation at the gyrocentre coordinate $\gy{\vc{Z}}$.
We are permitted to add a total derivative, by virtue of considering an action principle; we add the total derivative of

\begin{equation}
\roundpar{\gc{\vc{\gamma}} - \frac{1}{2} \pard{\gc{\vc{\gamma}}}{\gy{\vc{Z}}} \gy{\vc{G}}} \bcdot \roundpar{\gy{\vc{G}} - \frac{1}{2} \pard{\gy{\vc{G}}_1}{\gy{\vc{Z}}} \gy{\vc{G}}_1},
\end{equation}

resulting in the following symplectic part of the gyrocentre Lagrangian

\begin{align}
\gc{\Gamma}(\gc{\vc{Z}})
= {} &
\squarepar{
\gc{\vc{\gamma}}
- \pard{\gc{\vc{\gamma}}}{\gy{\vc{Z}}} \roundpar{\gy{\vc{G}} - \frac{1}{2} \pard{\gy{\vc{G}}_1}{\gy{\vc{Z}}} \gy{\vc{G}}_1}
+ \frac{1}{2} \pardd{\gc{\vc{\gamma}}}{\gy{\vc{Z}}} \bcddot \gy{\vc{G}}_1 \otimes \gy{\vc{G}}_1
} \bcdot \dot{\gy{\vc{Z}}}
+ \frac{1}{2} \roundpar{
\pard{\gc{\vc{\gamma}}}{\gy{\vc{Z}}} \gy{\vc{G}}
} \bcdot \roundpar{\pard{\gy{\vc{G}}}{t} + \pard{\gy{\vc{G}}}{\gy{\vc{Z}}} \dot{\gy{\vc{Z}}}} \\ &
+ \roundpar{\pard{\gc{\vc{\gamma}}}{t} + \pard{\gc{\vc{\gamma}}}{\gy{\vc{Z}}} \dot{\gy{\vc{Z}}}} \bcdot \roundpar{\gy{\vc{G}} - \frac{1}{2} \pard{\gy{\vc{G}}_1}{\gy{\vc{Z}}} \gy{\vc{G}}_1}
- \frac{1}{2} \squarepar{\pard{\gc{\vc{\gamma}}_0}{\gy{\vc{Z}}} \pard{\gy{\vc{G}}_1}{t} + \pard{}{\gy{\vc{Z}}} \roundpar{\pard{\gc{\vc{\gamma}}_0}{\gy{\vc{Z}}} \gy{\vc{G}}_1} \dot{\gy{\vc{Z}}}} \bcdot \gy{\vc{G}}_1
+ \mathcal{O}(\varepsilon_\delta^3).
\end{align}

\begin{align}
\gy{\vc{\gamma}}
= {} &
\gc{\vc{\gamma}}
+ \gc{\mat{W}} \roundpar{\gy{\vc{G}} - \frac{1}{2} \pard{\gy{\vc{G}}_1}{\gy{\vc{Z}}} \gy{\vc{G}}_1}
+ \frac{1}{2} \pardd{\gc{\vc{\gamma}}}{\gy{\vc{Z}}} \bcddot \gy{\vc{G}}_1 \otimes \gy{\vc{G}}_1 \nonumber \\ &
+ \frac{1}{2} \roundpar{\squarepar{\pard{\gy{\vc{G}}_1}{\gy{\vc{Z}}}}^\transpose \pard{\gc{\vc{\gamma}}}{\gy{\vc{Z}}}
- \squarepar{\pard{}{\gy{\vc{Z}}} \roundpar{\pard{\gc{\vc{\gamma}}_0}{\gy{\vc{Z}}} \gy{\vc{G}}_1}}^\transpose} \gy{\vc{G}}_1
\end{align}

where the following piece of the symplectic part of the transformed guiding-centre Lagrangian moves to the Hamiltonian part of the gyrocentre Lagrangian (while switching sign)

\begin{equation}\label{eq:gc_gy_symtransform_hampart1}
\pard{}{t} \roundpar{\gc{\vc{\gamma}}_1
- \frac{1}{2} \pard{\gc{\vc{\gamma}}_0}{\gy{\vc{Z}}} \gy{\vc{G}}_1} \bcdot \gy{\vc{G}}_1,
\end{equation}

where we have made use of $\partial \gc{\vc{\gamma}}_0 / \partial t = \vc{0}_6$.

To further simplify <ref> we now write it component wise as

\begin{align}
\gy{{\gamma}}_i
= {} &
\gc{{\gamma}}_i
+ \gc{\matcmp{W}}_{ij} \roundpar{\gy{{G}}_j - \frac{1}{2} \pard{\gy{{G}}_{1,j}}{\gy{{Z}}_k} \gy{{G}}_{1,k}}
+ \frac{1}{2} \frac{\partial \gc{{\gamma}}_i}{\partial \gy{{Z}}_j \partial \gy{{Z}}_k} \gy{{G}}_{1,j} \gy{{G}}_{1,k}
+ \frac{1}{2} \pard{\gy{{G}}_{1,j}}{\gy{{Z}}_i} \pard{\gc{{\gamma}}_j}{\gy{{Z}}_k} \gy{{G}}_{1,k} \nonumber\\ &
- \frac{1}{2} \pard{}{\gy{{Z}}_i} \roundpar{\pard{\gc{{\gamma}}_{0,j}}{\gy{{Z}}_k} \gy{{G}}_{1,k}} \gy{{G}}_{1,j}
\end{align}

This can be expanded to

\begin{align}
\gy{{\gamma}}_i
= {} &
\gc{{\gamma}}_i
+ \gc{\matcmp{W}}_{ij} \roundpar{\gy{{G}}_j - \frac{1}{2} \pard{\gy{{G}}_{1,j}}{\gy{{Z}}_k} \gy{{G}}_{1,k}}
+ \frac{1}{2} \roundpar{\frac{\partial^2 \gc{{\gamma}}_i}{\partial \gy{{Z}}_j \partial \gy{{Z}}_k} - \frac{\partial^2\gc{{\gamma}}_{0,j}}{\partial \gy{{Z}}_i \partial \gy{{Z}}_k}} \gy{{G}}_{1,j} \gy{{G}}_{1,k} \nonumber \\ &
+ \frac{1}{2} \gc{\matcmp{W}}_{kj} \pard{\gy{{G}}_{1,j}}{\gy{{Z}}_i} \gy{{G}}_{1,k}
\end{align}

where we have simplified the expression, and we note that a repeated index implies summation.
Note that

\begin{align}
\pard{}{\gy{Z}_i} \roundpar{\gc{\matcmp{W}}_{kj} \gy{{G}}_{1,j}} - \pard{}{\gy{Z}_k} \roundpar{\gc{\matcmp{W}}_{ij} \gy{{G}}_{1,j}}
= {} &
\roundpar{\frac{\partial^2 \gc{\gamma}_{i}}{\partial \gy{Z}_j \partial \gy{Z}_k} - \frac{\partial^2 \gc{\gamma}_{k}}{\partial \gy{Z}_j \partial \gy{Z}_i}}  \gy{{G}}_{1,j}
\nonumber \\ & + \gc{\matcmp{W}}_{kj} \pard{\gy{{G}}_{1,j}}{\gy{Z}_i} - \gc{\matcmp{W}}_{ij} \pard{\gy{{G}}_{1,j}}{\gy{Z}_k},
\end{align}

which allows <ref> to be written as

\begin{align}
\gy{{\gamma}}_i
= {} &
\gc{{\gamma}}_i
+ \gc{\matcmp{W}}_{ij}\gy{{G}}_j
+ \squarepar{\pard{}{\gy{Z}_i} \roundpar{\gc{\matcmp{W}}_{kj} \gy{{G}}_{1,j}} - \pard{}{\gy{Z}_k} \roundpar{\gc{\matcmp{W}}_{ij} \gy{{G}}_{1,j}}} \gy{{G}}_{1,k}
\end{align}

Finally, the gyrocentre Hamiltonian follows from subtracting <ref> from <ref> resulting in

\begin{equation}
\gy{H}
= \gc{H}
- \pard{\gc{H}}{\gy{\vc{Z}}} \bcdot \roundpar{\gy{\vc{G}}_1 + \gy{\vc{G}}_2}
+ \squarepar{
- \pard{}{t} \roundpar{\gc{\vc{\gamma}}_1 - \frac{1}{2} \pard{\gc{\vc{\gamma}}_0}{\gy{\vc{Z}}} \gy{\vc{G}}_1}
+ \frac{1}{2} \pard{}{\gy{\vc{Z}}}\roundpar{\pard{\gc{H}}{\gy{\vc{Z}}} \bcdot \gy{\vc{G}}_1}
} \bcdot \gy{\vc{G}}_1
\end{equation}

§ SUFFICIENT CONDITION FOR GAUGE INVARIANCE

We assume that $\gy{\vc{\gamma}}_1 - \gc{\vc{\gamma}}_1$ is gauge invariant.
Note that the perturbed guiding-centre Lagrangian (<ref>) is gauge invariant (up to a total derivative)

\begin{equation}
\gc{L}_1
%= q (\vc{A}_1 \bcdot \dot{\gc{\vc{R}}} - \phi_1)
\overset{\text{\cref{eq:gauge_transform}}}{\mapsto}
\gc{L}_1 + q \roundpar{\nabla \eta \bcdot \dot{\gc{\vc{R}}} + \pard{\eta}{t}}
\gc{L}_1 + q \totd{\eta}{t}.
\end{equation}

From <ref> it follows that the first-order gyrocentre Lagrangian can be written as

\begin{equation}
\gy{L}_1 = \gc{L}_1 + \roundpar{\gc{\mat{W}}_0 \gy{\vc{G}}_1 + \pard{\gy{S}_1}{\gy{\vc{Z}}}} \bcdot \dot{\gy{\vc{Z}}} + \pard{\gc{H}_0}{\gy{\vc{Z}}} \bcdot \gy{\vc{G}}_1 + \pard{\gy{S}_1}{t},
\end{equation}

from which it follows that $\gy{L}_1$ is gauge invariant if the first-order generating vector $\gy{\vc{G}}_1$ and first-order generating function $\gy{S}_1$ are gauge invariant.

Using <ref> we find that

\begin{equation}\label{eq:gy_S1_pde}
\pard{\gy{S}_1}{t} + \bracketgc{\gy{S}_1}{\gy{H}_0}
= q \wt{\psi}_1
= \wt{\gc{H}_1} + \dot{\gc{\vc{Z}}} \bcdot \wt{\roundpar{\gy{\vc{\gamma}}_1 - \gc{\vc{\gamma}}_1}}
= \wt{\gc{H}_1^\flr} + \dot{\gc{\vc{Z}}} \bcdot \wt{\roundpar{\gy{\vc{\gamma}}_1 - \gc{\vc{\gamma}}_1}},
\end{equation}

where we note that $\gc{H}_1^\flr$ is gauge invariant as it is expressed in terms of the electric field, and therefore $\gy{S}_1$ is gauge invariant since $\gy{\vc{\gamma}}_1 - \gc{\vc{\gamma}}_1$ is assumed to be gauge invariant.
It follows that the generating vector $\gy{\vc{G}}_1$, by <ref>, and thereby also the first-order gyrocentre Lagrangian are gauge invariant.

We additionally assume that $\gy{\vc{\gamma}}_2$ is gauge invariant.
The second-order gyrocentre Lagrangian is given by

\begin{equation}\label{eq:gy_L2}
\gy{L}_2 = \gy{\vc{\gamma}}_2 \bcdot \dot{\gy{\vc{Z}}} - \gy{H}_2,
\end{equation}

where the second-order gyrocentre Hamiltonian is given by <ref>.
In order for the second-order gyrocentre Lagrangian to be gauge invariant, we therefore need to show that the corresponding Hamiltonian is gauge invariant, which, upon inspection of <ref>, requires the vector $\vc{T}_1$

to be gauge invariant.

The vector $\vc{T}_1$, as given by <ref>, can be written as

\begin{equation}
\vc{T}_1
= \underbrace{\squarepar{\gc{\mat{W}}_1 \dot{\gc{\vc{Z}}} - \pard{\gc{H}_1}{\gy{\vc{Z}}} - \pard{\gc{\vc{\gamma}}_1}{t}}}_{\vc{T}_1^\dagger}
+ \frac{1}{2}\squarepar{\pard{\gc{\vc{\gamma}}_0}{\gy{\vc{Z}}} \pard{\gy{\vc{G}}_1}{t} + (\gy{\mat{W}}_1 - \gc{\mat{W}}_1) \dot{\gc{\vc{Z}}} + \pard{}{\gy{\vc{Z}}}\roundpar{\pard{\gc{H}_0}{\gy{\vc{Z}}} \bcdot \gy{\vc{G}}_1}},
\end{equation}

where we note that $\gy{\mat{W}}_1 - \gc{\mat{W}}_1$ is gauge invariant because $\gy{\vc{G}}_1$ is gauge invariant (see <ref>).
It therefore remains to be shown that $\vc{T}_1^\dagger$ is gauge invariant, which is expressed in terms of the first-order symplectic and Hamiltonian part of guiding-centre Lagrangian.
We now note that the FLR part of the first-order guiding-centre Lagrangian is expressed in terms of the electromagnetic fields (see <ref>) and is therefore gauge invariant.
Hence, we need only to show that the ZLR part $\vc{T}_1^{\dagger,\zlr}$ is gauge invariant, which is given by

\begin{align}
\vc{T}_1^{\dagger,\zlr}
= \gc{\mat{W}}_1^\zlr \dot{\gc{\vc{Z}}} - \pard{\gc{H}_1^\zlr}{\gy{\vc{Z}}} - \pard{\gc{\vc{\gamma}}_1^\zlr}{t}
% &= \begin{pmatrix}
%   q \mat{B}_1 & \vc{0}_3 & \vc{0}_3 & \vc{0}_3\\
%   \vc{0}_3^\transpose & 0 & 0 & 0\\
%   \vc{0}_3^\transpose & 0 & 0 & 0\\
%   \vc{0}_3^\transpose & 0 & 0 & 0\\
% \end{pmatrix}
% \begin{pmatrix}
%   \dot{\gc{\vc{R}}}\\
%   0\\
%   0\\
%   \omega_\mathrm{c}\\
% \end{pmatrix}
%  +
= \begin{pmatrix}
q (\vc{E}_1 + \dot{\gc{\vc{R}}} \times \vc{B}_1)\\
\end{pmatrix}
\end{align}

and is therefore gauge invariant.

§ GYRO-AVERAGING IDENTITIES

§.§ Taylor-series expansions

Let $Q = Q(\vc{r}, t)$ be a smooth function.
We denote by $\mat{H} Q$ the Hessian matrix of second-order partial derivatives of $Q$

\begin{equation}
(\mat{H} Q)_{ij} = \frac{\partial^2 Q}{\partial r_i \partial r_j},
\end{equation}

which allows us to write the Taylor series expansion of $Q$ centred around $\vc{r}$ as

\begin{equation}\label{eq:series_expansion}
\evalperturbedgc{Q} = Q(\vc{r} + \rho\hrho) = Q + \rho\hrho \bcdot\nabla Q + \frac{\rho^2}{2} \mat{H} Q \bcddot \hrho \otimes \hrho + \mathcal{O}(\rho^3).
\end{equation}

Where we recall <ref>, and note that an absence of arguments implies evaluation at $(\vc{r}, t)$.

Nondimensionalization of <ref> results in

\begin{equation}\label{eq:series_expansion_nondim}
\nondim{Q}(\vc{r} + \rho\hrho) = \nondim{Q} + \varepsilon_\perp \hrho \bcdot\nabla \nondim{Q} + \frac{\varepsilon_\perp^2}{2} \nondim{\mat{H}} \nondim{Q} \bcddot \hrho \otimes \hrho + \mathcal{O}(\varepsilon_\perp^3),
\end{equation}

where $Q = \units{Q} \nondim{Q}$, and we recall that the nondimensional perpendicular wavenumber is defined in <ref>.

We note that the $\mathcal{O}(\varepsilon_\perp^3)$ term in <ref> contains a $\hrho \otimes \hrho \otimes \hrho$ term, which yields zero upon computing the gyro-average.
It follows that the gyro-average of $\evalperturbedgc{Q}$ has the following Taylor series expansion

\begin{equation}\label{eq:gyro_identity_gav_lwa}
\gav{\evalperturbedgc{Q}} =
% + \frac{\rho^2}{4} \nabla_\perp^2 Q
+ \frac{\rho^2}{4} \mat{H} Q \bcddot (\vc{\hat{e}}_1 \otimes \vc{\hat{e}}_1 + \vc{\hat{e}}_2 \otimes \vc{\hat{e}}_2)
+ \mathcal{O}(\varepsilon_\perp^4),
\end{equation}

where we have multiplied <ref> by $\units{Q}$ and have made use of

\begin{align}
\gav{\hrho} &= \vc{0}_3,\\
\gav{\hrho \otimes \hrho} &= \frac{\vc{\hat{e}}_1 \otimes \vc{\hat{e}}_1 + \vc{\hat{e}}_2 \otimes \vc{\hat{e}}_2}{2}.\label{eq:gavrhorho}
\end{align}

§.§ Stokes' theorem

Using the definition of the gyro-average, as found in <ref>, we find that

\begin{equation}
\gav{\evalperturbedgc{S}_\tau}
= \frac{1}{2\pi} \int_0^{2\pi} \vc{S}(\gy{\vc{R}} + \vc{\rho}) \bcdot \htau(\gy{\vc{R}}, \gy{\theta}) \diff \gy{\theta}
= -\frac{1}{2\pi \rho} \int_{\partial D_\rho} \vc{S} \bcdot \vc{\hat{t}} \diff{l}
\end{equation}

where $\vc{\hat{t}} \perp \hb$ is the counter-clockwise tangent to the boundary of the disk $D_\rho$ centred at $\gy{\vc{R}}$ (i.e. the shaded disk shown in <ref>), which results in the minus sign.
Using Stokes' theorem and <ref>, we then find that

\begin{equation}\label{eq:gyro_identity_gav_tau}
\gav{\evalperturbedgc{S}_\tau}
= -\frac{1}{2\pi \rho} \int_{D_\rho} (\nabla \times \vc{S})_\shortparallel \diff^2 x
= - \frac{\rho}{2} \frac{1}{\pi\rho^2} \int_{D_\rho} (\nabla \times \vc{S})_\shortparallel \diff^2 x
= -\frac{\rho}{2} \dgav{\evalperturbedgcburby{\nabla} \times \vc{S}}_\shortparallel,
\end{equation}

where $\diff^2 x$ is an infinitesimal area element on the disk.

§.§ The gradient theorem

Terms containing a radial average (i.e. an integral over $\burby$) can often be interpreted in some alternative way.
For instance, the gradient theorem (or fundamental theorem of calculus for line integrals) results in

\begin{equation}\label{eq:gradient_theorem_scalar}
\rho \radgav{\hrho \bcdot \evalperturbedgcburby{\nabla} Q}
= \frac{1}{2\pi} \int_0^{2\pi}\int_0^1 \vc\rho \bcdot (\nabla Q)(\gy{\vc{R}} + \burby \vc\rho) \diff \burby \diff\gy{\theta}
= \gav{\evalperturbedgc{Q}} - Q.
\end{equation}

The same identity can be derived for vector fields

\begin{equation}\label{eq:gradient_theorem_vector}
\rho \radgav{\evalperturbedgcburby{\nabla} \vc{S} \hrho}
% = \frac{1}{2\pi} \int_0^{2\pi}\int_0^1 \vc\rho \bcdot (\nabla Q)(\gy{\vc{R}} + \burby \vc\rho) \diff \burby \diff\gy{\theta}
= \gav{\evalperturbedgc{\vc{S}}} - \vc{S}.
\end{equation}

By making use of $\nabla (\vc{a} \bcdot \vc{b}) = \vc{a} \times (\curl \vc{b}) + \vc{b} \times (\curl \vc{a}) + (\nabla \vc{b}) \vc{a} + (\nabla \vc{a}) \vc{b}$, while neglecting $\mathcal{O}(\varepsilon_B)$ contributions, we then find that

\begin{equation}\label{eq:gradient_theorem_vector_approx}
\radgav{\evalperturbedgcburby{\nabla} (\rho S_\rho)}
+ \rho \radgav{(\evalperturbedgcburby{\nabla} \times \vc{S}) \times \hrho}
+ \mathcal{O}(\varepsilon_B)
\gav{\evalperturbedgc{\vc{S}}} - \vc{S}.
\end{equation}

Similarly, we find that application of <ref> to the curl of a vector field results in

\begin{equation}
\rho \radgav{\evalperturbedgcburby{\nabla} (\curl \vc{S}) \hrho}
= \gav{\evalperturbedgc{\nabla} \times \vc{S}} - \curl \vc{S}.
\end{equation}

By making use of $\curl(\vc{a} \times \vc{b}) = \vc{a} \nabla \bcdot \vc{b} - \vc{b} \nabla \bcdot \vc{a} + (\nabla \vc{a}) \vc{b} - (\nabla \vc{b}) \vc{a}$, while neglecting $\mathcal{O}(\varepsilon_B)$ contributions, we then find that

\begin{equation}\label{eq:gradient_theorem_curlvector_approx}
\rho \radgav{\evalperturbedgcburby{\nabla} \times [(\curl \vc{S}) \times \hrho]}
+ \mathcal{O}(\varepsilon_B)
= \gav{\evalperturbedgc{\nabla} \times \vc{S}} - \curl \vc{S}.
\end{equation}

§ APPROXIMATION OF THE SECOND-ORDER HAMILTONIAN

We consider <ref>, with the first-order generating vector approximated by <ref> and $\vc{T}_{1}$ is as defined by <ref>, which we repeat here for convenience

\begin{equation}\label{eq:gy_Tvector_split}
\vc{T}_1 =
\underbrace{\frac{1}{2} (\gc{\mat{W}}_1 + \gy{\mat{W}}_1) \dot{\gc{\vc{Z}}}}_{\circled{A}}
- \underbrace{\pard{}{t} \roundpar{\gc{\vc{\gamma}}_1
- \frac{1}{2} \pard{\gc{\vc{\gamma}}_0}{\gy{\vc{Z}}} \gy{\vc{G}}_1}}_{\circled{B}}
- \underbrace{\pard{}{\gy{\vc{Z}}}\roundpar{\gc{H}_1 - \frac{1}{2} \pard{\gc{H}_0}{\gy{\vc{Z}}} \bcdot \gy{\vc{G}}_1}}_{\circled{C}}.
\end{equation}

In this section of the appendix, we use the following shorthand notation

\begin{equation}
\paramsum{Q}{R} \defeq \frac{1}{2} (Q + \gparamr \gav{Q}).
\end{equation}

For the approximation of $\gy{H}_2$ we omit all derivatives of the background magnetic field as well as of the perturbed electromagnetic fields.

Using <ref> we find that

\begin{equation}
\frac{1}{2} (\gc{\vc{\gamma}}_{1} + \gy{\vc{\gamma}}_{1}) =
\begin{pmatrix}
q\vc{A}_1 + q \rho \int_0^1 \paramsum{(\evalperturbedgcburby{\vc{B}}_1 \times \hrho)}{R} \diff \burby \\
- \frac{q \rho^2}{2} \int_0^1 \burby (\evalperturbedgcburby{B}_{1,\shortparallel}) \diff \burby\\
\end{pmatrix},
\end{equation}

from which it follows that

\begin{equation}
\frac{1}{2} \pard{}{\gy{\vc{Z}}} (\gc{\vc{\gamma}}_{1} + \gy{\vc{\gamma}}_{1}) =
\begin{pmatrix}
q\nabla \vc{A}_1 & \vc{0}_3 & \frac{q \rho}{2 \gy{\Mu}} \int_0^1 \paramsum{(\evalperturbedgcburby{\vc{B}}_1 \times \hrho)}{R} \diff \burby & \frac{q \rho}{2} \int_0^1  \evalperturbedgcburby{\vc{B}}_1 \times \htau \diff \burby\\
\vc{0}_3^\transpose & 0 & 0 & 0\\
\vc{0}_3^\transpose & 0 & 0 & 0\\
\vc{0}_3^\transpose & 0 & - \frac{q \rho^2}{2\gy{\Mu}} \int_0^1 \burby (\evalperturbedgcburby{B}_{1,\shortparallel}) \diff \burby & 0\\
\end{pmatrix}
\end{equation}

and therefore

\begin{equation}
\frac{1}{2} (\gc{\mat{W}}_{1} + \gy{\mat{W}}_{1}) =
\int_0^1 \begin{pmatrix}
q \mat{B}_1 & \vc{0}_3 & -\frac{q \rho}{2 \gy{\Mu}} \paramsum{(\evalperturbedgcburby{\vc{B}}_1 \times \hrho)}{R} & -\frac{q \rho}{2} \evalperturbedgcburby{\vc{B}}_1 \times \htau\\
\vc{0}_3^\transpose & 0 & 0 & 0\\
\frac{q \rho}{2 \gy{\Mu}} \squarepar{\paramsum{(\evalperturbedgcburby{\vc{B}}_1 \times \hrho)}{R}}^\transpose & 0 & 0 & -\frac{q \rho^2}{2\gy{\Mu}} \burby (\evalperturbedgcburby{B}_{1,\shortparallel})\\
\frac{q \rho}{2} \roundpar{ \evalperturbedgcburby{\vc{B}}_1 \times \htau}^\transpose & 0 & \frac{q \rho^2}{2\gy{\Mu}} \burby (\evalperturbedgcburby{B}_{1,\shortparallel}) & 0\\
\end{pmatrix} \diff\burby,
\end{equation}

where the matrix $\mat{B}_1$ is defined analogously to <ref>.
We have omitted all derivatives of the background magnetic field as well as of the perturbed electromagnetic fields.
Subsequent evaluation of the matrix-vector product with $\dot{\gc{\vc{Z}}}$ (while neglecting $\mathcal{O}(\varepsilon_B)$ contributions), where the guiding-centre EOMs are given by <ref>, results in

\begin{equation}
\begin{pmatrix}
q \mat{B}_1 & \vc{0}_3 & -\frac{q \rho}{2 \gy{\Mu}} \paramsum{(\evalperturbedgcburby{\vc{B}}_1 \times \hrho)}{R} & -\frac{q \rho}{2} \evalperturbedgcburby{\vc{B}}_1 \times \htau & \vc{0}_3 & -q \pard{\vc{A}_1}{t}\\
\vc{0}_3^\transpose & 0 & 0 & 0 & 0 & 0\\
\frac{q \rho}{2 \gy{\Mu}} \squarepar{\paramsum{(\evalperturbedgcburby{\vc{B}}_1 \times \hrho)}{R}}^\transpose & 0 & 0 & -\frac{q \rho^2}{\gy{\Mu}} \burby \paramsum{(\evalperturbedgcburby{B}_{1,\shortparallel})}{2} & 0 & 0\\
\frac{q \rho}{2} \roundpar{ \evalperturbedgcburby{\vc{B}}_1 \times \htau}^\transpose & 0 & \frac{q \rho^2}{2\gy{\Mu}} \burby (\evalperturbedgcburby{B}_{1,\shortparallel}) & 0 & 0 & 0\\
\vc{0}_3^\transpose & 0 & 0 & 0 & 0 & 0\\
q \roundpar{\pard{\vc{A}_1}{t}}^\transpose & 0 & 0 & 0 & 0 & 0\\
\end{pmatrix} \begin{pmatrix}
\gc{U}_\shortparallel \hb\\
\omega_\mathrm{c}\\
\end{pmatrix}
\end{equation}

\begin{equation}\label{eq:gy_H2_W1_times_Zdot}
\circled{A}
= \frac{1}{2} (\gc{\mat{W}}_{1} + \gy{\mat{W}}_{1}) \dot{\gc{\vc{Z}}}
= \int_0^1 \begin{pmatrix}
q \gc{U}_\shortparallel \hb \times \vc{B}_1 - \frac{q \rho \omega_\mathrm{c}}{2} \evalperturbedgcburby{\vc{B}}_1 \times \htau\\
\frac{q \rho \gc{U}_\shortparallel}{2 \gy{\Mu}} \paramsum{(\evalperturbedgcburby{\vc{B}}_1 \times \hrho)}{R} \bcdot \hb -\frac{\omega_\mathrm{c} q \rho^2}{2\gy{\Mu}} \burby (\evalperturbedgcburby{B}_{1,\shortparallel})\\
\frac{q \rho \gc{U}_\shortparallel}{2} (\evalperturbedgcburby{\vc{B}}_1 \times \htau) \bcdot \hb\\
\end{pmatrix} \diff \burby.
\end{equation}

The second term in $\vc{T}_1$ is given by

\begin{equation}
-\circled{B} = -
\begin{pmatrix}
q\pard{\vc{A}_1}{t}\\ % + q \int_0^1  \evalperturbedgcburby{\vc{B}}_1 \diff \burby \times \vc\rho\\
0\\%- q \rho^2 \int_0^1 \burby \evalperturbedgcburby{\vc{B}}_1 \diff \burby \bcdot \hb
\end{pmatrix}
+ \frac{1}{2}
\begin{pmatrix}
q \pard{\vc{A}_0}{\gy{\vc{R}}} & \vc{0}_3 & \vc{0}_3 & \vc{0}_3\\
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0\\
0 & 0 & \frac{m}{q} & 0
\end{pmatrix}
\begin{pmatrix}
\vc{0}_3\\
- \frac{\rho}{B_0} q \pard{E_{1,\rho}}{t}\\
-\frac{\rho}{2 B_0 \gy{\Mu}} q \pard{E_{1,\tau}}{t}
\end{pmatrix}
\end{equation}

\begin{equation}\label{eq:gy_T1_B}
\circled{B} =
\begin{pmatrix}
\end{pmatrix}
\end{equation}

by making use of <ref> as well as Faraday's equation.
The third term in $\vc{T}_1$ involves the average of the first-order guiding-centre and gyrocentre Hamiltonians

\begin{equation}
\frac{1}{2} (\gc{H}_1 + \gy{H}_1) =
q \phi_1
-  \int_0^1 \squarepar{\frac{q \rho}{2} (\evalperturbedgcburby{E}_{1,\rho} + \gav{\evalperturbedgcburby{E}_{1,\rho}})
+ (1 - \gparamr) \frac{q \rho}{2} \gy{U}_\shortparallel  \gav{\evalperturbedgcburby{B}_{1,\tau}}
- \gy{\Mu}  \burby \gav{\evalperturbedgcburby{B}_{1,\shortparallel}}} \diff\burby,
\end{equation}

which follows from substituting <ref>.
This results in

\begin{equation}\label{eq:gy_H2_partialderivativesH1}
\circled{C}
% = \frac{1}{2} \pard{}{\gy{\vc{Z}}} (\gc{H}_1 + \gy{H}_1)
= \int_0^1 \begin{pmatrix}
q \nabla \phi_1\\
-\frac{q \rho (1 - \gparamr)}{2} \gav{\evalperturbedgcburby{B}_{1,\tau}}\\
- \frac{q \rho}{4 \gy{\Mu}} (\evalperturbedgcburby{E}_{1,\rho} + \gav{\evalperturbedgcburby{E}_{1,\rho}}) - (1 - \gparamr) \frac{q \rho \gc{U}_\shortparallel}{4 \gy{\Mu}} \gav{\evalperturbedgcburby{B}_{1,\tau}}
+ \burby \gav{\evalperturbedgcburby{B}_{1,\shortparallel}}\\
- \frac{q \rho}{2} \evalperturbedgcburby{E}_{1,\tau} \\
\end{pmatrix} \diff \burby,
\end{equation}

where we have again neglected all derivatives of the background and perturbed electromagnetic fields.

When combining <ref> we find that $\vc{T}_1$ is given by

\begin{equation}\label{eq:gy_H2_T1}
\vc{T}_1 =
\int_0^1 \begin{pmatrix}
\vc{F}_1 - \frac{q \rho \omega_\mathrm{c}}{2} \evalperturbedgcburby{\vc{B}}_1 \times \htau\\
\frac{q \rho (1 - \gparamr)}{2} \gav{\evalperturbedgcburby{B}_{1,\tau}}\\
\frac{\rho}{4 \gy{\Mu}} (\evalperturbedgcburby{F}_{1,\rho} + \gav{\evalperturbedgcburby{F}_{1,\rho}})
- \burby (\evalperturbedgcburby{B}_{1,\shortparallel} + \gav{\evalperturbedgcburby{B}_{1,\shortparallel}})\\
\frac{\rho}{2} \evalperturbedgcburby{F}_{1,\tau} \\
\end{pmatrix} \diff \burby,
\end{equation}

where we have made use of the definition of the Lorentz force $\vc{F}_1$ as given by <ref>.
We are now ready to evaluate the second-order gyrocentre Hamiltonian, as expressed in <ref>, by substituting <ref>

\begin{align}
\gy{H}_{2} = &
\sgav{\squarepar{\vc{F}_1 - \frac{q \rho \omega_\mathrm{c}}{2} \int_0^1 \evalperturbedgcburby{\vc{B}}_1 \times \htau \diff \burby} \bcdot \squarepar{\frac{\hb \times}{B_{0,\shortparallel}^\effective} \int_0^1 \roundpar{\gparamr \gav{\evalperturbedgcburby{\vc{B}}_1 \times \vc\rho} - \evalperturbedgcburby{\vc{B}}_1 \times \vc\rho} \diff\burby + \frac{\rho B_{1,\rho}}{B_0} \hb}}\nonumber\\ &
+ \sgav{\squarepar{\frac{q \rho (1 - \gparamr)}{2} \int_0^1 \gav{\evalperturbedgcburby{B}_{1,\tau}} \diff \burby} \squarepar{- \frac{q \rho}{m} \int_0^1 \roundpar{\gparamr\gav{\evalperturbedgcburby{B}_{1,\tau}} - \evalperturbedgcburby{B}_{1,\tau}} \diff\burby}}\nonumber\\ &
- \sgav{\int_0^1 \squarepar{\frac{\rho}{4 \gy{\Mu}} (\evalperturbedgcburby{F}_{1,\rho} + \gav{\evalperturbedgcburby{F}_{1,\rho}}) - \burby (\evalperturbedgcburby{B}_{1,\shortparallel} + \gav{\evalperturbedgcburby{B}_{1,\shortparallel}}) }\diff \burby \squarepar{\frac{q^2 \rho^2}{m} \int_0^1 \burby \evalperturbedgcburby{B}_{1,\shortparallel} \diff\burby - \frac{\rho}{B_0} F_{1,\rho}}}\nonumber\\ &
+ \sgav{\squarepar{\frac{\rho}{2} \int_0^1 \evalperturbedgcburby{F}_{1,\tau} \diff \burby} \squarepar{-\frac{\rho}{2 B_0 \gy{\Mu}} F_{1,\tau}}}.
\end{align}

We consider the ZLR limit of all the terms, except for the term which moves to the particle part of the Hamiltonian, as discussed in <ref>.
That is, we make use of the approximation $\evalperturbedgcburby{Q} \approx Q$ resulting in

\begin{equation}
\begin{aligned}
\gy{H}_{2} = {} &
\frac{\rho}{B_0}\sgav{\squarepar{\vc{F}_1 - \frac{q \rho \omega_\mathrm{c}}{2} {\vc{B}}_1 \times \htau} \bcdot \squarepar{\hb \times \roundpar{\gparamr \gav{{\vc{B}}_1 \times \hrho} - {\vc{B}}_1 \times \hrho} + B_{1,\rho} \hb}}\\
& + \frac{1}{2} \sgav{\squarepar{\frac{\rho}{2 \gy{\Mu}} ({F}_{1,\rho} + \gav{{F}_{1,\rho}}) - {B}_{1,\shortparallel}} \squarepar{-\frac{q^2 \rho^2}{2 m} {B}_{1,\shortparallel} - \frac{\rho}{B_0} F_{1,\rho}}}\\
& - \frac{\rho^2}{4 B_0 \gy{\Mu}} \sgav{F_{1,\tau}^2}
+ (1 - \gparamr)^2 \frac{q^2 \rho^2}{2 m} \radgav{\evalperturbedgcburby{B}_{1,\tau}}^2,
\end{aligned}
\end{equation}

where we have moreover approximated $B_{0,\shortparallel}^\effective \approx B_0$.
When subsequently making use of $\gav{\htau} = \gav{\hrho} = \vc{0}_3$ as well as

\begin{equation}
\gav{S_\rho^2}
= \gav{S_\tau^2}
= \vc{S} \otimes \vc{S} \bcddot \gav{\hrho \otimes \hrho}
= \frac{1}{2} \vc{S} \otimes \vc{S} \bcddot (\vc{\hat{e}}_1 \otimes \vc{\hat{e}}_1 + \vc{\hat{e}}_2 \otimes \vc{\hat{e}}_2)
= \frac{1}{2} \abs{\vc{S}_\perp}^2,
\end{equation}

we find that

\begin{equation}\label{eq:gy_H2_gcF10}
\gy{H}_{2} =
\frac{\gy{\Mu}}{2 B_0} \abs{\vc{B}_{1,\perp}}^2
- \frac{m}{2 q^2 B_0^2} \abs{\vc{F}_{1,\perp}}^2
+ (1 - \gparamr)^2 \frac{q^2 \rho^2}{2 m} \radgav{\evalperturbedgcburby{B}_{1,\tau}}^2.
\end{equation}

§ COMPUTATION OF THE JACOBIAN

The aim is to compute the Jacobian of the transformation from the physical coordinates $\tilde{\vc{Z}} = (\vc{R}, \vc{U})$ to the gyrocentre coordinates $\gy{\vc{Z}}$.
To this end, we write the Lagrangian in physical coordinates $\tilde{\vc{Z}}$ as

\begin{equation}\label{eq:tilde_lagrangian}
\tilde{L} = \tilde{\vc{\gamma}} \bcdot \dot{\tilde{\vc{Z}}} - \tilde{H}.
\end{equation}

Here the Hamiltonian $\tilde{H}$ is defined such that, when the physical Lagrangian is transformed to gyrocentre coordinates, we find

\begin{equation}\label{eq:tilde_eq_gy}
\tilde{L}(\tilde{\vc{Z}}) = \gy{L}(\gy{\vc{Z}}).
\end{equation}

Letting $\tilde{\vc{Z}} = \tilde{\vc{Z}}(\gy{\vc{Z}})$, we find that

\begin{equation}
\dot{\tilde{\vc{Z}}} = \pard{\tilde{\vc{Z}}}{\gy{\vc{Z}}} \dot{\gy{\vc{Z}}},
\end{equation}

such that <ref> imply

\begin{equation}
\tilde{\vc{\gamma}} \bcdot \roundpar{\pard{\tilde{\vc{Z}}}{\gy{\vc{Z}}} \dot{\gy{\vc{Z}}}} = \gy{\vc{\gamma}} \bcdot \dot{\gy{\vc{Z}}}
\quad \implies \quad
\gy{\vc{\gamma}} = \roundpar{\pard{\tilde{\vc{Z}}}{\gy{\vc{Z}}}}^\transpose \tilde{\vc{\gamma}}.
\end{equation}

This results in the following expression for the gyrocentre Lagrange matrix

\begin{equation}
\gy{\mat{W}}
= \roundpar{\pard{\gy{\vc{\gamma}}}{\gy{\vc{Z}}}}^\transpose - \pard{\gy{\vc{\gamma}}}{\gy{\vc{Z}}}
= \roundpar{\pard{\tilde{\vc{\gamma}}}{\gy{\vc{Z}}}}^\transpose \pard{\tilde{\vc{Z}}}{\gy{\vc{Z}}}  - \roundpar{\pard{\tilde{\vc{Z}}}{\gy{\vc{Z}}}}^\transpose \pard{\tilde{\vc{\gamma}}}{\gy{\vc{Z}}}
= \roundpar{\pard{\tilde{\vc{Z}}}{\gy{\vc{Z}}}}^\transpose \tilde{\mat{W}} \pard{\tilde{\vc{Z}}}{\gy{\vc{Z}}},
\end{equation}

from which it follows that

\begin{equation}
\det \pard{\tilde{\vc{Z}}}{\gy{\vc{Z}}} = \sqrt{\frac{\det \gy{\mat{W}}}{\det \tilde{\mat{W}}}}.
\end{equation}

We note that

\begin{equation}
\det \tilde{\mat{W}}
= \det \begin{pmatrix}
q \mat{B}_0 & -m\mat{I}_3\\
m\mat{I}_3 & \mat{0}_3 \\
\end{pmatrix}
= m^6
\end{equation}

\begin{equation}
\det \gy{\mat{W}}
= \det \gy{\mat{S}} \det \gy{\mat{W}}_{22}
= \det \gy{\mat{W}}_{11} \frac{m^2}{q^2}
\end{equation}

by the Schur complement determinant formula, where we use the block structure of the Lagrange matrix as discussed in <ref>.
Finally, we note that direct computation shows that

\begin{equation}
\det \gy{\mat{W}}_{11}
= \det \begin{pmatrix}
q \mat{B}^\effective & -m \hb\\
m \hb^\transpose & 0
\end{pmatrix}
= (q m B_\shortparallel^\effective)^2,
\end{equation}

such that

\begin{equation}
\det \pard{\tilde{\vc{Z}}}{\gy{\vc{Z}}} = \frac{B_\shortparallel^\effective}{m}.
\end{equation}

§ PROOF OF LIOUVILLE'S THEOREM

The gyrocentre EOMs (<ref>) imply that

\begin{align}
\pard{B_{s,\shortparallel}^\effective}{t} &= \pard{\mathcal{B}_{1,\shortparallel}}{t},\\
\nabla \bcdot (B_{s,\shortparallel}^\effective \dot{\gy{\vc{R}}}) &= \frac{\vc{B}_s^\effective}{m_s} \bcdot \nabla \pard{\gy{H}_s}{\gy{U}_\shortparallel} + \frac{1}{q_s}\roundpar{\nabla \gy{H}_s + q_s \pard{\vc{\mathcal{A}}_1}{t}} \bcdot (\curl \hb) - \pard{\mathcal{B}_{1,\shortparallel}}{t},\\
\pard{}{\gy{U}_\shortparallel} (B_{s,\shortparallel}^\effective \dot{\gy{U}}_\shortparallel) &= -\frac{1}{q_s} (\curl \hb) \bcdot \roundpar{\nabla \gy{H}_s + q_s \pard{\vc{\mathcal{A}}_1}{t}} - \frac{\vc{B}_s^\effective}{m_s} \bcdot \pard{}{\gy{U}_\shortparallel} \nabla \gy{H}_s,
\end{align}

from which <ref> follows by simply adding the contributions.

We consider the coordinate transformation $\gy{\vc{z}} = \gy{\vc{Z}}(t; \gyzzero, \tzero) \mapsto \gyzzero$ applied to the left-hand side of <ref> resulting in

\begin{multline}\label{eq:liouville_intermediate}
\int \gy{f}_s(\gy{\vc{z}}, t) \mathcal{F}(\gy{\vc{z}}, t) \jacobiangy(\gy{\vc{z}}, t) \diff^6 \gy{z}
= \\
\int \gyfzero(\gyzzero) \mathcal{F}(\gy{\vc{Z}}(t; \gyzzero, \tzero)) \jacobiangy(\gy{\vc{Z}}(t; \gyzzero, \tzero), t) \mathfrak{H}(\gy{\vc{Z}}(t; \gyzzero, \tzero), t) \diff^6 \gyzzerod
\end{multline}

by defining the Jacobian

\begin{equation}
\mathfrak{H} = \det \pard{\gy{\vc{Z}}}{\gyzzero}
\end{equation}

and by making use of <ref>.
Therefore, we must show that

\begin{equation}
\totd{}{t} (\jacobiangy \mathfrak{H}) = 0,
\end{equation}

where we note that the product $\jacobiangy \mathfrak{H}$ is the Jacobian of the coordinate transformation from initial gyrocentre coordinates to physical coordinates.

We note that <ref> implies

\begin{equation}\label{eq_gy_liouville_res0}
\dotjacobiangy = - \jacobiangy \pard{}{\gy{\vc{Z}}} \bcdot \dot{\gy{\vc{Z}}}.
\end{equation}

In order to proceed, we make use of (for a proof we refer to <cit.>)

\begin{equation}
\dot{\mathfrak{H}}
% = \mathfrak{H} \tr \pard{\dot{\gy{\vc{Z}}}}{\gy{\vc{Z}}}
= \mathfrak{H} \pard{}{\gy{\vc{Z}}} \bcdot \dot{\gy{\vc{Z}}}.
\end{equation}

When combining this result with <ref> we find

\begin{equation}
\totd{}{t} (\jacobiangy \mathfrak{H})
= \dotjacobiangy \mathfrak{H} + \jacobiangy \dot{\mathfrak{H}}
% = \jacobiangy \roundpar{1 - \det \pard{\gy{\vc{Z}}}{\gyzzero}} \pard{}{\gy{\vc{Z}}} \bcdot \dot{\gy{\vc{Z}}}
= \dotjacobiangy \mathfrak{H} + \jacobiangy \mathfrak{H} \pard{}{\gy{\vc{Z}}} \bcdot \dot{\gy{\vc{Z}}}
= \dotjacobiangy \mathfrak{H} - \mathfrak{H} \dotjacobiangy
= 0
\end{equation}

It follows that

\begin{equation}
\jacobiangy \mathfrak{H}
= \jacobiangy(\gyzzero, \tzero) \mathfrak{H}(\gyzzero, \tzero)
= \jacobiangy(\gyzzero, \tzero),
\end{equation}

which, upon substitution in <ref>, results in <ref>.

§ PROOF OF LOCAL ENERGY CONSERVATION

A direct computation of the partial derivative of the kinetic energy density (<ref>) per species results in

\begin{align}
\pard{\gy{\mathcal{K}}_s}{t}
={} & \int \pard{}{t} (\gy{f}_s \jacobiangy) \gy{K}_s \diff^3 \gy{u}
- \int \gy{f}_s \gy{\mu} \dgav{(\evalperturbedgcburby{\nabla} \times \vc{E}_{1})_\shortparallel} \jacobiangy \diff^3 \gy{u} \nonumber\\
={} & - \nabla \bcdot \int \gy{f}_s \dot{\gy{\vc{R}}} \gy{K}_s \jacobiangy \diff^3 \gy{u}
+ \int \gy{f}_s \gy{\mu} \dot{\gy{\vc{R}}} \bcdot \nabla (B_0 + \dgav{\evalperturbedgcburby{B}_{1,\shortparallel}}) \jacobiangy \diff^3 \gy{u} \nonumber\\ &
+ m_s \int \gy{f}_s \dot{\gy{U}}_\shortparallel \gy{u}_\shortparallel \jacobiangy \diff^3 \gy{u}
- \int \gy{f}_s \gy{\mu} \dgav{(\evalperturbedgcburby{\nabla} \times \vc{E}_{1})_\shortparallel} \jacobiangy \diff^3 \gy{u} \nonumber\\
={} & - \nabla \bcdot \int \gy{f}_s \dot{\gy{\vc{R}}} \gy{K}_s \jacobiangy \diff^3 \gy{u}
+ \frac{1}{m_s} \int \gy{f}_s \gy{\mu} \roundpar{\gy{u}_\shortparallel \vc{B}_s^\effective
- \hb \times \vc{\mathcal{E}}_1} \bcdot \nabla (B_0 + \dgav{\evalperturbedgcburby{B}_{1,\shortparallel}}) \diff^3 \gy{u} \nonumber\\ &
+ \frac{1}{m_s} \int \gy{f}_s  \vc{B}_s^\effective \bcdot \squarepar{
q_s \vc{\mathcal{E}}_1
- \gy{\mu} \nabla (B_0 + \dgav{\evalperturbedgcburby{B}_{1,\shortparallel}})} \gy{u}_\shortparallel \diff^3 \gy{u}
- \int \gy{f}_s \gy{\mu} \dgav{(\evalperturbedgcburby{\nabla} \times \vc{E}_{1})_\shortparallel} \jacobiangy \diff^3 \gy{u} \nonumber\\
={} & - \nabla \bcdot \int \gy{f}_s \dot{\gy{\vc{R}}} \gy{K}_s \jacobiangy \diff^3 \gy{u}
+ q_s \int \gy{f}_s \dot{\gy{\vc{R}}} \bcdot \vc{\mathcal{E}}_1 \jacobiangy \diff^3 \gy{u}
- \int \gy{f}_s \gy{\mu} \dgav{(\evalperturbedgcburby{\nabla} \times \vc{E}_{1})_\shortparallel} \jacobiangy \diff^3 \gy{u},
\end{align}

upon substitution of Faraday's equation (<ref>), the conservative form of the Vlasov equation (<ref>) and the gyrocentre EOMs (<ref>).
It follows that the total kinetic energy evolves as <ref> upon substitution of the definition of the free current density (<ref>) and summation over the species $s$.

For computing the partial time derivative of the potential energy density (<ref>) we note that

\begin{equation}\label{eq:pot_evo_part1}
\pard{}{t} (\vc{\mathcal{D}} \bcdot \vc{E}) =
2\pard{\vc{\mathcal{D}}}{t} \bcdot \vc{E}
% + \sum_s \frac{m_s}{B_0^2} \int \gyfzero \pard{\vc{E}_{1,\perp}}{t} \jacobiangcslab \diff^3 \gy{u} \bcdot \vc{E}
+ \sum_s \frac{m_s}{B_0^2} \int \gyfzero \gy{u}_\shortparallel \squarepar{(\hb \times \vc{B}_1) \bcdot \pard{\vc{E}}{t} - \roundpar{\hb \times \pard{\vc{B}_1}{t}} \bcdot \vc{E}} \jacobiangcslab \diff^3 \gy{u}
\end{equation}

and similarly

\begin{equation}\label{eq:pot_evo_part2}
\pard{}{t} (\vc{\mathcal{H}} \bcdot \vc{B}) =
2\vc{\mathcal{H}} \bcdot \pard{\vc{B}}{t}
+ \sum_s \roundpar{\pard{\vc{B}}{t} \bcdot \vc{\mathcal{M}}_{1,s} - \pard{\vc{\mathcal{M}}_{1,s}}{t} \bcdot \vc{B}}.
\end{equation}

Upon substitution of the strong form of Ampère's law (<ref>) as well as Faraday's equation (<ref>) we find that

\begin{equation}
\pard{\vc{\mathcal{D}}}{t} \bcdot \vc{E} + \vc{\mathcal{H}} \bcdot \pard{\vc{B}}{t}
= (\curl \vc{\mathcal{H}} - \evalperturbedgcadjointalt{\vc{\mathcal{J}}}{}^\free) \bcdot \vc{E} - \vc{\mathcal{H}} \bcdot (\curl \vc{E})
= - \nabla \bcdot (\vc{E} \times \vc{\mathcal{H}}) - \evalperturbedgcadjointalt{\vc{\mathcal{J}}}{}^\free \bcdot \vc{E}.
\end{equation}

It remains to be shown that the two remaining terms from <ref> vanish.
The sum of the two terms is written as

\begin{align}
\sum_s \frac{m_s}{B_0^2} \int \gyfzero \Biggl( &
\gy{u}_\shortparallel \squarepar{(\hb \times \vc{B}_1) \bcdot \pard{\vc{E}_1}{t} - \roundpar{\hb \times \pard{\vc{B}_1}{t}} \bcdot \vc{E}_1} \nonumber\\ &
- \pard{\vc{B}_1}{t} \bcdot \squarepar{\gy{u}_\shortparallel \hb \times \vc{E}_{1,\perp} + \roundpar{\frac{\gy{\mu} B_0}{m_s} - \gy{u}_\shortparallel^2} \vc{B}_{1,\perp}} \nonumber \\ &
+ \pard{}{t} \squarepar{\gy{u}_\shortparallel \hb \times \vc{E}_{1,\perp} + \roundpar{\frac{\gy{\mu} B_0}{m_s} - \gy{u}_\shortparallel^2} \vc{B}_{1,\perp}} \bcdot \vc{B}_1
\biggr) \jacobiangcslab \diff^3 \gy{u}
\end{align}

upon substitution of the definition of the magnetization and polarization as defined in Eqs. (<ref>) and (<ref>), respectively.
Indeed, the sum vanishes and therefore <ref> also holds.

§ DERIVATION OF THE BRIZARD–HAHM HAMILTONIAN

The aim is to compute the ZLR approximation of the second-order Hamiltonian as found in <cit.>.
Their first- and second-order Hamiltonian are given by

\begin{equation}
\gy{H}_1^\bh = q \gav{{\psi}_1^\bh}
, \quad
{\psi}_1^\bh \defeq \evalperturbedgc{\phi}_1 - \gy{u}_\shortparallel \wt{A}_{1,\shortparallel} - \gy{u}_\tau \evalperturbedgc{A}_{1,\tau}
\end{equation}

\begin{equation}
\gy{H}_2^\bh = - \frac{q}{2} \gav{\bracketgc{{\gy{S}}_1^\bh}{\wt{\psi}{}_1^\bh}} + \frac{q^2}{2m} \squarepar{\gav{\abs{\evalperturbedgc{\vc{A}}_{1,\perp}}^2} + \gav{(\wt{A}_{1,\shortparallel})^2}} + \frac{q}{B_0} \gav{\evalperturbedgc{\vc{A}}_{1,\perp}} \bcdot (\hb \times \nabla \gav{{\psi}_1^\bh}),
\end{equation}

The first-order generating function $\gy{S}_1^\bh$ again satisfies <ref> to zeroth-order in $\varepsilon_\omega$.

§.§ ZLR approximation of the Poisson bracket

We use the following shorthand notation for the components of a vector field $\vc{S}$ in terms of the local coordinates $\vc{\hat{e}}_1, \vc{\hat{e}}_2$

\begin{equation}
S_{\hat{i}} \defeq \vc{S} \bcdot \vc{\hat{e}}_i
\end{equation}

as well as the corresponding directional derivative of a scalar function

\begin{equation}
\hatdd{i} Q \defeq \vc{\hat{e}}_i \bcdot \nabla Q.
\end{equation}

Provided with $\gy{S}_1^\bh$, we want to approximate the leading order term in $\varepsilon_\perp$ of the gyro-averaged guiding-centre Poisson bracket (<ref>) while neglecting $\mathcal{O}(\varepsilon_\omega)$ and $\mathcal{O}(\varepsilon_B)$ terms.
This results in

\begin{equation}\label{eq:H2_poissonbracket_term_tmp}
\gav{\bracketgc{{\gy{S}}_1^\bh}{\wt{\psi}{}_1^\bh}} =
\underbrace{\frac{q}{B_0} \sgav{\pard{}{\gy{\mu}} (\wt{\psi}{}_1^\bh)^2}}_{\circled{A}}
- \underbrace{\frac{\hb}{q B_0}\bcdot \sgav{\nabla_\perp {\gy{S}}_1^\bh \times \nabla_\perp \wt{\psi}{}_1^\bh}}_{\circled{B}}
+ \mathcal{O}(\varepsilon_\perp^3).
\end{equation}

Application of <ref> to $\wt{\psi}{}_1^\bh$ results in

\begin{align}\label{eq:psit_smallLR}
\wt{\psi}{}_1^\bh
= {} &
- \frac{B_0^2}{m} \vc{\rho} \bcdot \vc{Q}_{1}
- \rho \vtau \hrho \bcdot \nabla_\perp {A_{1,\tau}}
- \frac{\rho \vtau}{2} (\curl \Adperp)_\shortparallel
- \frac{\rho^2 \vtau}{2} \hrho \otimes \hrho \bcddot \mat{H} {A_{1,\tau}}
\nonumber\\ &
+ \frac{\rho^2}{2} \hrho \otimes \hrho \bcddot \mat{H} (\phid - \gy{u}_\shortparallel \Apa)
- \frac{\rho^2}{2} \nabla_\perp^2 (\phid - \gy{u}_\shortparallel \Apa)
+ \mathcal{O}(\varepsilon_\perp^3),
\end{align}

which shows that the first term can be written as

\begin{align}
\circled{A}
% &= \frac{q}{B_0} \pard{}{\gy{\mu}} \frac{1}{2\pi}\int_0^{2\pi}
%   \squarepar{
%     \frac{B_0^2}{m} \vc{\rho} \bcdot \vc{Q}_{1}
%     + \rho \vtau \hrho \bcdot \nabla_\perp \roundpar{\htau \bcdot \Adperp}
%     + \frac{\rho \vtau}{2} (\curl \Adperp)_\shortparallel
% }^2
% \diff\theta\nonumber\\
= {} & \frac{q}{B_0} \pard{}{\gy{\mu}}
\Biggl\langle
\underbrace{\roundpar{\frac{B_0^2}{m} \vc{\rho} \bcdot \vc{Q}_{1}}^2}_{\circled{A$_1$}}
+ \underbrace{\frac{\rho^2 \vtausqr}{4} [(\curl \Adperp)_\shortparallel]^2}_{\circled{A$_2$}}
+ \underbrace{\rho^2 \vtausqr \hrho \bcdot \nabla_\perp \roundpar{\htau \bcdot \Adperp} (\curl \Adperp)_\shortparallel}_{\circled{A$_3$}}\nonumber\\
& + \underbrace{\squarepar{\rho \vtau \hrho \bcdot \nabla_\perp \roundpar{\htau \bcdot \Adperp}}^2}_{\circled{A$_4$}}
+ \underbrace{\rho^2 \vtausqr A_{1,\tau} \hrho \otimes \hrho \bcddot \mat{H} {A_{1,\tau}}}_{\circled{A$_5$}}
\Biggr\rangle,
\end{align}

where we made use of the fact that integration over the interval $[0, 2\pi]$ of odd powers of cosine and/or sine functions yields zero, and we have defined the vector $\vc{Q}_{1}$ as

\begin{equation}\label{eq:polarization_bh}
\vc{Q}_{1} \defeq -\frac{m}{B_0^2} \squarepar{\nabla_\perp (\phi_1 - \gy{u}_\shortparallel A_{1,\shortparallel}) + \omega_{\mathrm{c}} \Adperp \times \hb}.
\end{equation}

Moreover, we used the fact that the $\Adperp$ term in $\vc{Q}_{1}$ is the only $\mathcal{O}(\varepsilon_\perp^0)$ term of $\wt{\psi}_1$, and therefore the only term that remains upon multiplication by the $\rho^2 \mat{H} A_{1,\tau}$ term, when neglecting $\mathcal{O}(\varepsilon_\perp^3)$ terms and after integration over $\theta$.
Recall that the outer product is defined in <ref>.

The first contribution can be evaluated as follows

\begin{equation}
\circled{A$_1$} = \frac{B_0^2}{q m} \abs{\vc{Q}_{1}}^2,
\end{equation}

where we made use of <ref>.
The second contribution can trivially be evaluated as

\begin{equation}
\circled{A$_2$} = \frac{2 \gy{\mu}}{q B_0} [(\curl \Adperp)_\shortparallel]^2,
\end{equation}

whereas the third contribution yields

\begin{equation}
\circled{A$_3$}
= \frac{4 \gy{\mu}}{q B_0} (\curl \Adperp)_\shortparallel \roundpar{\hatdd{2} A_{1,\hat{1}} - \hatdd{1} A_{1,\hat{2}}}
= -\frac{4 \gy{\mu}}{q B_0} (\curl \Adperp)_\shortparallel^2
\end{equation}

by making use of

\begin{equation}
\gav{\hrho \otimes \htau}
% = \frac{1}{2 \pi}\int_0^{2\pi} (\vc{\hat{e}}_1 \cos\theta - \vc{\hat{e}}_2 \sin\theta)(- \vc{\hat{e}}_1 \sin\theta - \vc{\hat{e}}_2 \cos\theta)^\transpose \diff\theta
= \frac{1}{2} (\vc{\hat{e}}_2 \otimes \vc{\hat{e}}_1 - \vc{\hat{e}}_1 \otimes \vc{\hat{e}}_2).
\end{equation}

For the fourth contribution we expand the gyro-phase dependent part as follows

\begin{equation}
\circled{A$_4$} = \frac{q}{B_0} \pard{}{\gy{\mu}} \sgav{\rho^2\vtausqr
\squarepar{\cos\theta \hatdd{1} \roundpar{\sin\theta A_{1,\hat{1}} + \cos\theta A_{1,\hat{2}}} - \sin\theta \hatdd{2} \roundpar{\sin\theta A_{1,\hat{1}} + \cos\theta A_{1,\hat{2}}}}^2},
\end{equation}

which can further be expanded as follows

\begin{multline}
\circled{A$_4$} = \frac{q}{B_0} \pard{}{\gy{\mu}} \rho^2\vtausqr
\biggl\langle
\sin^4\theta (\hatdd{2} A_{1,\hat{1}})^2
+ \cos^4\theta (\hatdd{1} A_{1,\hat{2}})^2
+ \sin^2\theta \cos^2\theta \roundpar{
(\hatdd{1} A_{1,\hat{1}})^2 + (\hatdd{2} A_{1,\hat{2}})^2
\\
-2 \sin^2\theta \cos^2\theta \roundpar{
(\hatdd{2} A_{1,\hat{1}}) (\hatdd{1} A_{1,\hat{2}})
+ (\hatdd{2} A_{1,\hat{2}}) (\hatdd{1} A_{1,\hat{1}})
\biggr\rangle ,
\end{multline}

The fourth contribution can be written as

\begin{equation}
\circled{A$_4$} = \frac{\gy{\mu}}{q B_0}
\roundpar{
\abs{\nabla_\perp \Adperp}^2
+ 2 [(\curl \Adperp)_\shortparallel]^2
- 2 \det \mat{G}_\perp
\end{equation}

where we make use of

\begin{equation}\label{eq:sin4_cos4}
\gav{\sin^4\theta} = \gav{\cos^4\theta} = \frac{3}{8}
, \quad
\gav{\sin^2\theta \cos^2\theta} = \frac{1}{8}
\end{equation}

and define

\begin{equation}
\mat{G}_\perp \defeq \begin{pmatrix}
\hatdd{1} A_{1,\hat{1}} & \hatdd{2} A_{1,\hat{1}}\\
\hatdd{1} A_{1,\hat{2}} & \hatdd{2} A_{1,\hat{2}}
\end{pmatrix}.
\end{equation}

Finally, the fifth contribution can be written as

\begin{multline}
\circled{A$_5$} = \frac{q}{B_0} \pard{}{\gy{\mu}} \biggl\langle \rho^2 \vtausqr
\squarepar{
\cos^2\theta \vc{\hat{e}}_1 \vc{\hat{e}}_1^\transpose
+ \sin^2\theta \vc{\hat{e}}_2 \vc{\hat{e}}_2^\transpose
- \sin\theta \cos\theta (\vc{\hat{e}}_2  \vc{\hat{e}}_1^\transpose
+ \vc{\hat{e}}_1 \vc{\hat{e}}_2^\transpose)
\bcddot \squarepar{
\sin^2\theta A_{1,\hat{1}} \mat{H} A_{1,\hat{1}}
+ \cos^2\theta A_{1,\hat{2}} \mat{H} A_{1,\hat{2}}
+ \sin\theta \cos\theta (A_{1,\hat{1}} \mat{H} A_{1,\hat{2}}
+ A_{1,\hat{2}} \mat{H} A_{1,\hat{1}})
\biggr\rangle.
\end{multline}

By making use of the fact that integration over the interval $[0, 2\pi]$ of odd powers of cosine and/or sine functions yields zero, we find that this can be further simplified to

\begin{multline}
\circled{A$_5$} = \frac{q}{B_0} \pard{}{\gy{\mu}} \frac{1}{2\pi}\int_0^{2\pi} \rho^2 \vtausqr\\
\biggl[\roundpar{
\cos^2\theta \vc{\hat{e}}_1 \vc{\hat{e}}_1^\transpose
+ \sin^2\theta \vc{\hat{e}}_2 \vc{\hat{e}}_2^\transpose
\bcddot \roundpar{
\sin^2\theta A_{1,\hat{1}} \mat{H} A_{1,\hat{1}}
+ \cos^2\theta A_{1,\hat{2}} \mat{H} A_{1,\hat{2}}
\sin^2\theta \cos^2\theta \roundpar{
\vc{\hat{e}}_2  \vc{\hat{e}}_1^\transpose
+ \vc{\hat{e}}_1 \vc{\hat{e}}_2^\transpose
\bcddot \roundpar{A_{1,\hat{1}} \mat{H} A_{1,\hat{2}}
+ A_{1,\hat{2}} \mat{H} A_{1,\hat{1}}
\diff\theta,
\end{multline}

which can be written as

\begin{multline}
\circled{A$_5$} = \frac{\gy{\mu}}{q B_0}
\biggl[
3 A_{1,\hat{1}} \hatdd{2}^2 A_{1,\hat{1}}
+ 3 A_{1,\hat{2}} \hatdd{1}^2 A_{1,\hat{2}}
+   A_{1,\hat{1}} \hatdd{1}^2 A_{1,\hat{1}}
+   A_{1,\hat{2}} \hatdd{2}^2 A_{1,\hat{2}}\\
- 2 A_{1,\hat{1}} \hatdd{1} \hatdd{2} A_{1,\hat{2}}
- 2 A_{1,\hat{2}} \hatdd{1} \hatdd{2} A_{1,\hat{1}}
\biggr],
\end{multline}

by making use of <ref> and $\rho^2 \vtausqr = 4\gy{\mu}^2 / q^2$.
This expression can be written more compactly as

\begin{equation}
\circled{A$_5$} = \frac{\gy{\mu}}{q B_0} \squarepar{
\Adperp \bcdot \nabla_\perp^2 \Adperp
- 2 (\Adperp \times \nabla_\perp [(\curl \Adperp)_\shortparallel])_\shortparallel
% - 2 A_{1,\hat{1}} \hatdd{2} (\curl \Adperp)_\shortparallel
% + 2 A_{1,\hat{2}} \hatdd{1} (\curl \Adperp)_\shortparallel
\end{equation}

Since the $\circled{B}$ term in <ref> already has a factor $\varepsilon_\perp^2$, we need to only keep the zeroth order term of $\wt{\psi}_1$ (see <ref>) for obtaining an $\mathcal{O}(\varepsilon_\perp^3)$ approximation

\begin{equation}
\wt{\psi}_1 = \vtau A_{1,\tau} + \mathcal{O}(\varepsilon_\perp)
\end{equation}

and similarly for $S_1$ for the computation of the second term in <ref> resulting in

\begin{equation}\label{eq:bracket_other_contribution}
\circled{B}
% = \frac{\vtausqr}{B_0 \omega_{\mathrm{c}}} \frac{1}{2\pi} \int_0^{2\pi} \hb \bcdot\nabla_\perp A_{1,\rho} \times \nabla_\perp A_{1,\tau} \diff\theta + \mathcal{O}(\varepsilon_\perp^3)
= -\frac{\vtausqr}{B_0 \omega_{\mathrm{c}}}
\hb \bcdot \nabla_\perp A_{1,\hat{1}} \times \nabla_\perp A_{1,\hat{2}} + \mathcal{O}(\varepsilon_\perp^3)
= -\frac{2 \gy{\mu}}{q B_0} \det \mat{G}_\perp + \mathcal{O}(\varepsilon_\perp^3)
% &\approx -\frac{2 \gy{\mu}}{q B_0} \det \mat{G}_\perp + \mathcal{O}(\varepsilon_\perp^3).
\end{equation}

By making use of

\begin{equation}
\frac{1}{2} \nabla_\perp^2 \abs{\Adperp}^2 = \abs{\nabla_\perp \Adperp}^2 + \Adperp \bcdot \nabla_\perp^2 \Adperp,
\end{equation}

we then find that the contribution to $\gy{H}_2^\bh$ due to the first-order generating function can be approximated by (neglecting $\mathcal{O}(\varepsilon_\perp^3)$ terms)

\begin{equation}
\gav{\bracketgc{{\gy{S}}_1^\bh}{\wt{\psi}_1^\bh}} =
\frac{B_0^2}{q m} \abs{\vc{Q}_{1}}^2
+ \frac{\gy{\mu}}{q B_0} \roundpar{
\frac{1}{2} \nabla_\perp^2 \abs{\Adperp}^2
- 2 \hb \bcdot \Adperp \times \nabla \squarepar{(\curl \Adperp)_\shortparallel}
\end{equation}

Note that in principle a high-frequency approximation of the first-order generating function $\gy{S}_1^\bh$ can also be considered, such as in the work of [Qin et al., 1999].

§.§ ZLR approximation of the Hamiltonian

The remaining gyro-averaged terms in $\gy{H}_2^\bh$ can be approximated in the ZLR limit as follows

\begin{align}
\frac{q^2}{2m} \gav{(\wt{A}_{1,\shortparallel})^2}
&= \frac{\gy{\mu}}{2 B_0} \abs{\nabla_\perp \Apa}^2 + \mathcal{O}(\varepsilon_\perp^4),\label{eq:fluc_mag_parallel_approx}\\
\frac{q^2}{2m} \gav{\abs{\evalperturbedgc{\vc{A}}_{1,\perp}}^2} &= \frac{q^2}{2m} \abs{\Adperp}^2 + \frac{\gy{\mu}}{4 B_0} \nabla_\perp^2 \abs{\Adperp}^2 + \mathcal{O}(\varepsilon_\perp^4),\label{eq:fluc_mag_perp_approx}\\
\gav{\evalperturbedgc{\vc{A}}_{1,\perp}} &= \vc{A}_{1,\perp} + \frac{\rho^2}{4} \nabla_\perp^2 \vc{A}_{1,\perp} + \mathcal{O}(\varepsilon_\perp^4)
\end{align}

by making use of <ref>.
Moreover, we note that

\begin{equation}
\gav{\evalperturbedgc{\psi}_1^\bh} =
\phi_1 - \gy{u}_\shortparallel A_{1,\shortparallel} + \frac{\rho^2}{4} \nabla_\perp^2 (\phi_1 - \gy{u}_\shortparallel A_{1,\shortparallel}) + \frac{\gy{\mu}}{q} (\nabla \times \vc{A}_1)_\shortparallel + \mathcal{O}(\varepsilon_\perp^3),
\end{equation}

by making use of <ref>.
Combining the ZLR approximations results in the following ZLR approximation of the second-order Hamiltonian

\begin{equation}
\gy{H}_2^{\bh} = - \frac{B_0^2}{2 m} \abs{\vc{Q}_{1}}^2
+ \frac{\gy{\mu}}{2 B_0} \abs{\nabla_\perp \Apa}^2
+ \frac{q^2}{2m} \abs{\Adperp}^2
+ \frac{q}{B_0} \vc{A}_{1,\perp} \bcdot [\hb \times \nabla_\perp (\phi_1 - \gy{u}_\shortparallel A_{1,\shortparallel})]
% + \frac{q}{B_0} \frac{\rho^2}{4} \roundpar{\vc{A}_{1,\perp} \bcdot \hb \times \nabla \nabla_\perp^2 \phi_1
%   +  \nabla_\perp^2 \vc{A}_{1,\perp} \bcdot \hb \times \nabla \phi_1},
\end{equation}

which can be simplified to <ref>.

§ WELL-POSEDNESS OF THE SADDLE-POINT PROBLEM

Here we briefly, and rather informally, discuss the well-posedness  of the saddle-point problem given by <ref> under the simplifying assumptions of a constant background magnetic field, an isotropic pressure $p_{0,\perp} = p_{0,\shortparallel}$ as well as a spatially constant coefficient $\mathcal{C}(1)$.
Therefore, the following system of equations are considered

\begin{align}
\curl (\curl \vc{A}) - \nabla_\perp \lambda &= \vc{\mathcal{J}},\label{eq:saddle_point_curl}\\
\nabla_\perp \bcdot \vc{A} &= 0,\label{eq:saddle_point_div}
\end{align}

where $\nabla \bcdot \vc{\mathcal{J}} = 0$ and with appropriate boundary conditions.
Assume that $\lambda = 0$ and that $\vc{A}$ satisfies

\begin{align}
-\nabla_\perp^2 A_\shortparallel &= \mathcal{J}_\shortparallel,\label{eq:saddle_point_parallel}\\
-\nabla^2 \vc{A}_\perp &= \vc{\mathcal{J}}_\perp - \nabla_\perp (\hb \bcdot \nabla A_\shortparallel),\label{eq:saddle_point_perp}
\end{align}

then we can show that $(\lambda, \vc{A})$ also satisfy <ref> provided that $\varepsilon_B = 0$.
Since each of the PDEs in <ref> is well-posed provided with appropriate boundary conditions, we find that <ref> is well-posed as well.

The equivalence can be shown as follows.
Assume that $\vc{A}$ solves <ref>.
Computing the divergence of <ref> results in

\begin{equation}
\nabla^2 (\nabla_\perp \bcdot \vc{A})
= \nabla_\perp \bcdot \vc{\mathcal{J}} - \hb \bcdot \nabla (\nabla_\perp^2 A_\shortparallel)
= \nabla_\perp \bcdot \vc{\mathcal{J}} + \hb \bcdot \nabla (\mathcal{J}_\shortparallel)
= \nabla \bcdot \vc{\mathcal{J}}
= 0
\end{equation}

upon substitution of <ref>.
If homogenous Dirichlet boundary conditions are `included' in the Laplace operator, then this implies that <ref> holds.
Finally, we add $\hb$ times <ref> to <ref> resulting in

\begin{equation}
\hb (\hb \bcdot \nabla)^2 A_\shortparallel - \nabla^2 \vc{A}
= \vc{\mathcal{J}} - \nabla_\perp (\hb \bcdot \nabla A_\shortparallel).
\end{equation}

When rearranging the terms and by making use of <ref>, we find that

\begin{equation}
\nabla (\nabla \bcdot \vc{A}) - \nabla^2 \vc{A}
= \vc{\mathcal{J}},
\end{equation}

which shows that $\vc{A}$ satisfies <ref>.

§ THE SUSCEPTIBILITY TENSOR IN A SLAB

For the computation of the susceptibility tensor we require a linearization of our model under a Fourier ansatz.
In particular, we require the Fourier component of the linearized distribution function, the free current, the polarization current, and finally the magnetization current.
Once each of these Fourier components is known, we combine the results into the gyrokinetic susceptibility tensor, and we compute the drift kinetic limit.

§.§ The Fourier component of the linearized distribution function

We need a closed-form expression for $\gyfone$, which can be obtained upon substitution of the Fourier ansatz

\begin{equation}\label{eq:fourier_f1}
\gyfone = \fourier{\gyfone}(\gy{u}_\shortparallel, \gy{\mu}) \euler^{\ii (\vc{k} \bcdot \gy{\vc{r}} - \omega t)}
\end{equation}

into the linearized Vlasov equation (cf. <ref>)

\begin{equation}
B_{0,\shortparallel}^\effective \pard{\gyfone}{t}
+ B_{0,\shortparallel}^\effective \gy{u}_\shortparallel \hb \bcdot \nabla \gyfone
% + \squarepar{
%       \gy{u}_\shortparallel \gav{\evalperturbedgc{\vc{B}}_{1,\perp}}
%     + \roundpar{
%         \gav{\evalperturbedgc{\vc{E}}_1}
%       - \frac{\mu}{q_s} \nabla \dgav{\evalperturbedgcburby{B}_{1,\shortparallel}}
%     } \times \hb
% } \bcdot \nabla \gyfzero\\
+ (B_{s,\shortparallel}^\effective \dot{\gy{U}}_\shortparallel)_0 \pard{\gyfone}{\gy{U}_\shortparallel}
+ (B_{s,\shortparallel}^\effective \dot{\gy{U}}_\shortparallel)_1 \pard{\gyfzero}{\gy{U}_\shortparallel}
= 0.
\end{equation}

\begin{equation}
\pard{\gyfone}{t}
+ \gy{u}_\shortparallel \hb \bcdot \nabla \gyfone
% + \frac{1}{B_{0,\shortparallel}^\effective}\squarepar{
%       \gy{u}_\shortparallel \gav{\evalperturbedgc{\vc{B}}_{1,\perp}}
%     + \roundpar{
%         \gav{\evalperturbedgc{\vc{E}}_1}
%       - \frac{\mu}{q_s} \nabla \dgav{\evalperturbedgcburby{B}_{1,\shortparallel}}
%     } \times \hb
% } \bcdot \nabla \gyfzero\\
- \frac{1}{m_s} \hb \bcdot \roundpar{\mu \nabla \dgav{\evalperturbedgcburby{B}_{1,\shortparallel}}
- q_s \gav{\evalperturbedgc{\vc{E}}_1}} \pard{\gyfzero}{\gy{u}_\shortparallel}
= 0,
\end{equation}

where we have made use of

\begin{equation}\label{eq:gy_eomR_linpart}
(B_{s,\shortparallel}^\effective \dot{\gy{\vc{R}}})_0 = \vc{B}_0 \gy{u}_\shortparallel
, \quad
(B_{s,\shortparallel}^\effective \dot{\gy{\vc{R}}})_1 =
\gy{u}_\shortparallel \gav{\evalperturbedgc{\vc{B}}_{1}}
+ \roundpar{
\gav{\evalperturbedgc{\vc{E}}_1}
- \frac{\mu}{q_s} \nabla \dgav{\evalperturbedgcburby{B}_{1,\shortparallel}}
} \times \hb
\end{equation}

using the gyrocentre EOM (<ref>) and <ref>

as well as

\begin{equation}
(B_{s,\shortparallel}^\effective \dot{\gy{U}}_\shortparallel)_0 = 0
, \quad
(B_{s,\shortparallel}^\effective \dot{\gy{U}}_\shortparallel)_1 = -\frac{1}{m}
% \roundpar{
\vc{B}_0
% + \gav{\evalperturbedgc{\vc{B}}_{1,\perp}}}
\bcdot \roundpar{\mu \nabla \dgav{\evalperturbedgcburby{B}_{1,\shortparallel}}
- q_s \gav{\evalperturbedgc{\vc{E}}_1}},
\end{equation}

which follows from substituting <ref> into the gyrocentre EOM (<ref>) and having assumed $\nabla \gyfzero = \vc{0}_3$.
Subsequent substitution of the Fourier ansatzes (<ref>) and Eq. (<ref>) (and similarly for the magnetic field) results in

\begin{equation}\label{eq:lin_vlasov_fourier_gav}
- \ii \omega \fourier{\gyfone}
+ \ii \gy{u}_\shortparallel \hb \bcdot \vc{k} \fourier{\gyfone}
% + \frac{1}{B_{0,\shortparallel}^\effective} \squarepar{
%       \gy{u}_\shortparallel \fourier{\gav{\evalperturbedgc{\vc{B}}_{1,\perp}}}
%     + \roundpar{
%         \fourier{\gav{\evalperturbedgc{\vc{E}}_1}}
%       - \ii \frac{\mu}{q_s} \vc{k} \fourier{\dgav{\evalperturbedgcburby{B}_{1,\shortparallel}}}
%     } \times \hb
% } \bcdot \nabla \gyfzero\\
- \frac{1}{m_s} \hb \bcdot \roundpar{\ii \mu \vc{k} \fourier{\dgav{\evalperturbedgcburby{B}_{1,\shortparallel}}}
- q_s \fourier{\gav{\evalperturbedgc{\vc{E}}_1}}} \pard{\gyfzero}{\gy{u}_\shortparallel}
= 0.
\end{equation}

In order to proceed, we must express the Fourier component of the gyro-average of a function in terms of the Fourier component of that function itself.
We write the wavevector $\vc k$ in terms of its parallel and perpendicular components as

\begin{equation}\label{eq:wavevector}
\vc k = \kpa \hb + \kperp(\vc{\hat{e}}_1 \cos \alpha + \vc{\hat{e}}_2 \sin \alpha)
\end{equation}

for some angle $\alpha$.
It follows that the gyro-average of a function $Q$

\begin{equation}
Q(\vc{r}) = \int \fourier{Q}(\vc{k}) \euler^{\ii\vc{k} \bcdot \vc{r}} \diff^3 k
\end{equation}

can be expressed in terms of the Fourier component in the following way

\begin{equation}
\gav{\evalperturbedgc{Q}}(\vc{r})
% = \frac{1}{2\pi} \int_0^{2\pi} \int \fourier{Q} \euler^{\ii\vc{k} \bcdot (\vc{r} + \vc{\rho}(\theta))} \diff^3 k \diff \theta
% = \frac{1}{2\pi} \int_0^{2\pi} \int \fourier{Q} \euler^{\ii\vc{k} \bcdot \vc{r}} \euler^{\ii \kperp \rho (\cos\theta\cos\alpha - \sin\theta\sin\alpha)} \diff^3 k \diff \theta
% = \frac{1}{2\pi} \int \fourier{Q} \euler^{\ii\vc{k} \bcdot \vc{r}} \int_0^{2\pi} \euler^{\ii \kperp \rho \cos (\theta + \alpha)} \diff \theta \diff^3 k
= \frac{1}{2\pi} \int \fourier{Q} \euler^{\ii\vc{k} \bcdot \vc{r}} \int_0^{2\pi} \euler^{\ii \kperp \rho \cos \theta} \diff \theta \diff^3 k
= \int \fourier{Q} \euler^{\ii\vc{k} \bcdot \vc{r}} J_0 \diff^3 k\label{eq:gav_bessel0}
\end{equation}

and therefore

\begin{equation}\label{eq:fourier_gav}
\fourier{\gav{\evalperturbedgc{Q}}} = J_0 \fourier{Q},
\end{equation}

where $J_n$ denotes the $n$-th order Bessel function of the first kind evaluated at $\kperp \rho$

\begin{equation}\label{eq:bessel_n}
J_n = \frac{1}{2\pi \ii^n} \int_0^{2\pi} \euler^{\ii \kperp \rho \cos \theta} \euler^{\ii \theta n} \diff \theta.
\end{equation}

A similar computation for the disc average, <ref>, yields

\begin{equation}\label{eq:fourier_dgav}
\fourier{\dgav{\evalperturbedgcburby{Q}}}
= 2 \int_0^1 \burby J_0(\burby \kperp \rho) \diff \burby \, \fourier{Q}
= \frac{2 J_1}{\kperp \rho} \fourier{Q}.
\end{equation}

Substitution of <ref> in <ref> yields

\begin{equation}
(- \ii \omega + \ii \kpa \gy{u}_\shortparallel) \fourier{\gyfone}
% + \frac{1}{B_0} \squarepar{
%       \gy{u}_\shortparallel J_0 \fourier{\vc{B}}_{1,\perp}
%     + \roundpar{
%         J_0 \fourier{\vc{E}}_1
%       - \ii \frac{\mu}{q_s} \frac{2 J_1}{\kperp \rho} \fourier{B}_{1,\shortparallel} \vc{k}
%     } \times \hb
% } \bcdot \nabla \gyfzero\\
+ \frac{q_s}{m_s} \hb \bcdot \roundpar{
J_0 \fourier{\vc{E}}_1
- \ii \frac{\mu}{q_s} \frac{2 J_1}{\kperp \rho} \fourier{B}_{1,\shortparallel} \vc{k}
} \pard{\gyfzero}{\gy{u}_\shortparallel}
= 0,
\end{equation}

\begin{multline}
\fourier{\gyfone}
\frac{\gy{u}_\shortparallel J_0 \vc{k} \times \fourier{\vc{E}}_1 \bcdot \nabla_\perp \gyfzero}
{\ii \omega B_0 (\omega - \kpa \gy{u}_\shortparallel)}\\
+ \frac{\squarepar{
\omega J_0 \fourier{\vc{E}}_1
- \ii \frac{\mu}{q_s} \frac{2 J_1}{\kperp \rho} (\vc{k} \times \fourier{\vc{E}}_1 \bcdot \hb) \vc{k}
} \bcdot \roundpar{\hb \times \nabla \gyfzero + \omega_{\mathrm{c},s} \pard{\gyfzero}{\gy{u}_\shortparallel} \hb}}
{\ii \omega B_0 (\omega - \kpa \gy{u}_\shortparallel)}
\end{multline}

\begin{equation}\label{eq:lin_vlasov_fourier}
\fourier{\gyfone} = \pertfvec \bcdot \fourier{\vc{E}}_1
, \quad
\pertfvec
\defeq
\roundpar{
J_0 \hb
- \ii \frac{2 \gy{\mu} \kpa J_1}{\omega q_s \kperp \rho} \hb \times \vc{k}}
\frac{
\omega_{\mathrm{c},s} \pard{\gyfzero}{\gy{u}_\shortparallel}
{\ii B_0 (\omega - \kpa \gy{u}_\shortparallel)}
\end{equation}

where we have made use of Faraday's equation in Fourier space

\begin{equation}\label{eq:Faraday_fourier}
\fourier{\vc{B}}_1 = \frac{1}{\omega} \vc{k} \times \fourier{\vc{E}}_1.
\end{equation}

§.§ The Fourier component of the linearized gyrocentre free current

We write the gyrocentre current (<ref>) as

\begin{equation}
\evalperturbedgcadjointalt{\vc{\mathcal{J}}}{}^\free = \sum_s \evalperturbedgcadjoint{\vc{\mathcal{J}}}_{s}
, \quad
\evalperturbedgcadjoint{\vc{\mathcal{J}}}{}^\free_s \defeq \int \roundpar{
q_s \gav{\evalperturbedgcadjoint{\gy{f}_s \jacobiangy  \dot{\gy{\vc{R}}}}}
- \gy{\mu} \dgav{\evalperturbedgcadjoint{\curl (\gy{f}_s \jacobiangy \hb)}}
} \diff^3 \gy{u},
\end{equation}

where we have made use of <ref> with $\varepsilon_B = 0$ and define the gyro-average adjoint $\gav{\evalperturbedgcadjoint{Q}}$ as (and similarly for the disc average)

\begin{equation}
\int \gav{\evalperturbedgcadjoint{Q}} \cp \diff^3 \gy{r}
\defeq
\int Q \gav{\evalperturbedgc{\cp}} \diff^6 \gy{z}
\end{equation}

for all suitable test functions $\cp$.

Linearization of the gyrocentre free current (<ref>), as is required for computing the susceptibility tensor, results in

\begin{equation}\label{eq:djdt_tmp}
\evalperturbedgcadjoint{\vc{\mathcal{J}}}{}^\free_{1,s} =
\frac{q_s}{m_s} \int \biggl(
\gyfzero \gav{\evalperturbedgcadjoint{(B_{s,\shortparallel}^\effective \dot{\gy{\vc{R}}})_1}}
+ \gav{\evalperturbedgcadjoint{\gyfone}} \vc{B}_0 \gy{u}_\shortparallel
- \gyfzero \frac{\gy{\mu}}{q_s} \dgav{\evalperturbedgcadjoint{\curl \vc{B}_{1,s,\shortparallel}^\effective}}