Polarization effects in higher-order guiding-centre Lagrangian dynamics

Abstract

The extended guiding-centre Lagrangian equations of motion are derived by the Lie-transform perturbation method under the assumption of time-dependent and inhomogeneous electric and magnetic fields that satisfy the standard guiding-centre space–time orderings. Polarization effects are introduced into the Lagrangian dynamics by the inclusion of the polarization drift velocity in the guiding-centre velocity and the appearance of finite-Larmor-radius corrections in the guiding-centre Hamiltonian and guiding-centre Poisson bracket.

1. Introduction

Polarization effects have a rich history in plasma physics (Pfirsch 1984; Pfirsch & Morrison 1985; Cary & Brizard 2009). Their importance stems from the assumption of quasineutrality in a strongly magnetized plasma and the dielectric properties of a guiding-centre plasma (Hinton & Robertson 1984). While these effects are traditionally associated with the presence of an electric field in a magnetized plasma (Itoh & Itoh 1996; Hazeltine & Meiss 2003; Wang & Hahm 2009; Joseph 2021; Brizard 2023a), they are also associated with magnetic drifts (Kaufman 1986; Brizard 2013; Tronko & Brizard 2015).

Recently, second-order terms in guiding-centre Hamiltonian theory (in the absence of an electric field) were shown to be crucial (Brizard & Hodgeman 2023) in assessing the validity of the guiding-centre representation in determining whether guiding-centre orbits were numerically faithful to the particle orbits in axisymmetric magnetic geometries, which partially confirmed earlier numerical studies in axisymmetric tokamak plasmas (Belova, Gorelenkov & Cheng 2003). In particular, it was shown that a second-order correction associated with guiding-centre polarization (Kaufman 1986; Brizard 2013; Tronko & Brizard 2015) was needed in order to obtain faithful guiding-centre orbits.

Indeed, without the inclusion of second-order effects, it was shown that, within a few bounce periods after leaving the same physical point in particle phase space, a first-order guiding-centre orbit deviated noticeably from its associated particle orbit, while a second-order guiding-centre orbit followed the particle orbit to a high degree of precision (Brizard & Hodgeman 2023). In addition, as initially reported by Belova et al. (2003), the guiding-centre Hamiltonian formulation is a faithful representation of the particle toroidal angular momentum (Tronko & Brizard 2015; Brizard & Hodgeman 2023), which is an exact particle constant of motion in an axisymmetric magnetic field, only if second-order effects are included.

1.1. Lagrangian dynamics in extended phase space

In the present work, we consider the motion of a charged particle (of mass m and charge e) in time-dependent and inhomogeneous electric (E) and magnetic (B = B $\hat{\textsf {b}}$) fields (which still satisfy the guiding-centre space–time orderings $|\boldsymbol {\nabla }|^{-1} \gg \rho _{\rm th} = v_{\rm th}/\varOmega \def\mathsfbit#1{\boldsymbol{\mathsf{#1}}}$, i.e., characteristic spatial scales are long compared to the thermal gyroradius $\rho _{\rm th}$, and $\partial /\partial t \ll \varOmega = eB/mc$, i.e., time scales are long compared to the gyroperiod, where c denotes the speed of light) and we assume that the $E\times B$ velocity $\boldsymbol {u}_{E} = \boldsymbol {E}\boldsymbol {\times } c\hat {\textsf {b}}/B$ is comparable to the particle's thermal velocity $v_{\rm th}$. Because of the explicit time dependence of the electromagnetic fields, the Lagrangian charged-particle dynamics takes place in an odd-dimensional space $(\boldsymbol {q},\boldsymbol {p},t)$, where the non-autonomous Hamiltonian $H(\boldsymbol {q},\boldsymbol {p},t)$ is a function of the canonical coordinates $(\boldsymbol {q},\boldsymbol {p})$, from which the canonical Hamilton equations ${\rm d}\boldsymbol {q}/{\rm d} t = \partial H/\partial \boldsymbol {p}$ and ${\rm d}\boldsymbol {p}/{\rm d} t = -\partial H/\partial \boldsymbol {q}$ are derived, as well as time $t$, from which we obtain the energy equation ${\rm d} H/{\rm d} t = \partial H/\partial t$ (i.e. energy is not conserved).

The use of an extended phase space is a well-known method in classical mechanics (Lanczos 1970) used to deal with a time-dependent Hamiltonian system by transforming it into an autonomous Hamiltonian system in an even-dimensional symplectic setting. Here, the canonical time–energy coordinates $(t,w)$ are included in the extended phase-space coordinates $(\boldsymbol {q},t;\boldsymbol {p},w)$, where the space–time coordinates $(\boldsymbol {q},t)$ are canonically conjugate to the momentum–energy coordinates $(\boldsymbol {p},w)$, with the extended Hamilton equations ${\rm d} w/{\rm d} s = \partial \mathcal {H}/\partial t = \partial H/\partial t$ and ${\rm d} t/{\rm d} s = -\partial \mathcal {H}/\partial w = 1$, where $\mathcal {H} \equiv H(\boldsymbol {q},\boldsymbol {p},t) - w$ is the extended Hamiltonian and a particle orbit in extended phase space (parametrized by $s$) takes place on the energy surface $\mathcal {H} = 0$, i.e. $w = H(\boldsymbol {q},\boldsymbol {p},t)$.

Using the dimensional ordering parameter $\epsilon$ associated with the renormalized particle mass $m \rightarrow \epsilon m$ (Brizard 1995), instead of the standard ordering $e \rightarrow e/\epsilon$ (Kulsrud 1983; Littlejohn 1983), we begin with the extended phase-space particle Lagrangian one-form

(1.1)\begin{equation} \gamma = \left(\frac{e}{c}\boldsymbol{A} + \epsilon \boldsymbol{p}_{0}\right) \boldsymbol{\cdot}\textsf{d}\boldsymbol{x} - w\,\textsf{d} t \equiv \gamma_{0} + \epsilon \gamma_{1}, \end{equation}

and the extended particle Hamiltonian

(1.2)\begin{equation} \mathcal{H} = e \varPhi - w + \epsilon |\boldsymbol{p}_{0}|^{2}/2m \equiv \mathcal{H}_{0} + \epsilon \mathcal{H}_{1}, \end{equation}

where $\boldsymbol {p}_{0}$ denotes the local particle kinetic momentum at position $\boldsymbol {x}$. In the present work, we consider the standard ordering (Kulsrud 1983) for the parallel electric field: $\boldsymbol {E} = \boldsymbol {E}_{\bot } + \epsilon E_{\|}\hat {\textsf {b}}$. In contrast to Madsen (2010) and Frei, Jorge & Ricci (2020), who used the same mass ordering ($m \rightarrow \epsilon m$), we use extended (eight-dimensional) phase space in (1.1)–(1.2), where the energy coordinate $w$ is canonically conjugate to time $t$ (Littlejohn 1981; Cary & Brizard 2009). This extended phase-space formulation yields a simple form for the extended Poisson bracket (see (4.5)), also adopted (without derivation) by Madsen (2010), which plays an important role in the variational formulation of the guiding-centre Vlasov–Maxwell equations (Brizard 2023b).

Here, the electric-field ordering implies that the local particle momentum

(1.3)\begin{equation} \boldsymbol{p}_{0} \equiv p_{\|0}\hat{\textsf{b}}(\boldsymbol{x},t) + \boldsymbol{P}_{E}(\boldsymbol{x},t) + \boldsymbol{q}_{\bot0}(J_{0},\theta_{0};\boldsymbol{x},t) \end{equation}

is decomposed into the gyroangle-independent parallel component $p_{\|0} \equiv \boldsymbol {p}_{0}\boldsymbol {\cdot }\hat {\textsf {b}}$ and the $E\times B$ momentum

(1.4)\begin{equation} \boldsymbol{P}_{E} \equiv \boldsymbol{E}\boldsymbol{\times}\frac{e\hat{\textsf{b}}}{\varOmega} = m\boldsymbol{u}_{E}, \end{equation}

and the gyroangle-dependent perpendicular momentum $\boldsymbol {q}_{\bot 0} \equiv |\boldsymbol {q}_{\bot 0}|\hat {\bot }$, respectively, where $J_{0} \equiv |\boldsymbol {q}_{\bot 0}|^{2}/(2m\varOmega )$ represents the lowest-order gyroaction and the gyroangle derivative $\partial \boldsymbol {q}_{\bot 0}/\partial \theta _{0} \equiv \boldsymbol {q}_{\bot 0}\boldsymbol {\times }\hat {\textsf {b}} = -|\boldsymbol {q}_{\bot 0}|\hat {\rho }$ introduces the rotating orthogonal unit-vector fields $(\hat {\textsf {b}},\hat {\bot },\hat {\rho })$.

1.2. Purpose of the present work

The purpose of the present work is motivated by the need to derive higher-order guiding- centre equations that can accurately describe the magnetic confinement of charged particles in regions with steep gradients (e.g. the pedestal region of advanced tokamak plasmas). For many situations of practical interest, the presence of a strong electric field is associated with strong plasma flows with steep sheared rotation profiles for which second-order effects (including finite-Larmor-radius effects) must be included in a self-consistent guiding-centre theory (Hahm 1996; Chang, Ku & Weitzner 2004; Lanthaler et al. 2019; Frei et al. 2020).

Guiding-centre equations of motion with second-order corrections in the presence of time-independent electric and magnetic fields were derived using Lie-transform perturbation method by Brizard (1995) and Hahm (1996), following the earlier work of Littlejohn (1981). These perturbation methods were also used by Miyato et al. (2009), Madsen (2010) and Frei et al. (2020), who derived self-consistent guiding-centre Vlasov–Maxwell equations that included guiding-centre polarization and magnetization effects. Not all second-order effects were included in these models, however, and it is the purpose of the present work to derive a more complete higher-order guiding-centre Vlasov–Maxwell theory, with a full representation of guiding-centre polarization that can be directly derived by the guiding-centre push-forward method (Brizard 2013; Tronko & Brizard 2015).

2. Guiding-centre Lie-transform perturbation analysis

The derivation of the guiding-centre equations of motion by Lie-transform perturbation method is based on a phase-space transformation from the (local) particle extended (eight-dimensional) phase-space coordinates $z_{0}^{\alpha } = (\boldsymbol {x}, p_{\|0}; J_{0},\theta _{0};w_{0},t)$, where the energy–time canonical coordinates $(w_{0},t)$ are included, to the guiding-centre phase-space coordinates $Z^{\alpha } = (\boldsymbol {X}, P_{\|}; J, \theta ;W,t)$ generated by the vector fields $(\boldsymbol{\mathsf{G}}_{1}, \boldsymbol{\mathsf{G}}_{2}, \ldots )$

(2.1)\begin{equation} Z^{\alpha} = z_{0}^{\alpha} + \epsilon G_{1}^{\alpha} + \epsilon^{2}\left(G_{2}^{\alpha} + \tfrac{1}{2}\boldsymbol{\mathsf{G}}_{1} \boldsymbol{\cdot}\textsf{d} G_{1}^{\alpha}\right) + \cdots, \end{equation}

and its inverse

(2.2)\begin{equation} z_{0}^{\alpha} = Z^{\alpha} - \epsilon G_{1}^{\alpha} - \epsilon^{2}\left(G_{2}^{\alpha} - \tfrac{1}{2}\boldsymbol{\mathsf{G}}_{1} \boldsymbol{\cdot}\textsf{d} G_{1}^{\alpha}\right) + \cdots. \end{equation}

In order for the particle time to be invariant under the guiding-centre transformation (2.1), we require that $G_{n}^{t} \equiv 0$ to all orders $n \geq 1$, i.e. the guiding-centre time is identical to the particle time. From these generating vectors fields, the pull-back and push-forward Lie-transform operators $\boldsymbol{\mathsf{T}}_{\rm gc} = \exp (\epsilon \pounds _{1})\exp (\epsilon ^{2}\pounds _{2})\cdots$ and $\boldsymbol{\mathsf{T}}_{\rm gc}^{-1} = \cdots \exp (-\epsilon ^{2}\pounds _{2})\exp (-\epsilon \pounds _{1})$ are constructed in terms of Lie derivatives $(\pounds _{1}, \pounds _{2}, \ldots )$ generated by the vector fields $(\boldsymbol{\mathsf{G}}_{1}, \boldsymbol{\mathsf{G}}_{2}, \ldots )$. More details about the Lie-transform perturbation method used in guiding-centre theory can be found in Littlejohn (1982) (as well as the unpublished UCLA report Variational Principles for Guiding Center Motion (PPG-611) written by Littlejohn in 1982), while the notation used here is taken from Tronko & Brizard (2015) and Brizard & Tronko (2016).

We now wish to derive the extended guiding-centre phase-space Lagrangian one-form

(2.3)\begin{equation} \varGamma_{\rm gc} \equiv {\boldsymbol{\mathsf{T}}}_{\rm gc}^{-1}\gamma + \textsf{d} \sigma = \varGamma_{0{\rm gc}} + \epsilon \varGamma_{1{\rm gc}} + \epsilon^{2}\varGamma_{2{\rm gc}} + \cdots, \end{equation}

where the gauge scalar field $\sigma = \epsilon \sigma _{1} + \epsilon ^{2}\sigma _{2} + \cdots$ is chosen at each order in order to simplify the transformation, with

(2.4)\begin{gather} \varGamma_{1{\rm gc}} = \gamma_{1} - \iota_{1} \boldsymbol{\cdot}\boldsymbol{\omega}_{0} + \textsf{d} \sigma_{1}, \end{gather}

(2.5)\begin{gather}\varGamma_{2{\rm gc}} \equiv-\iota_{2}\boldsymbol{\cdot}\boldsymbol{\omega}_{0} - \tfrac{1}{2} \iota_{1}\boldsymbol{\cdot}(\boldsymbol{\omega}_{1} +\boldsymbol{\omega}_{{\rm gc}1}) + \textsf{d} \sigma_{2}, \end{gather}

and the extended guiding-centre Hamiltonian

(2.6)\begin{equation} \mathcal{H}_{\rm gc} \equiv {\boldsymbol{\mathsf{T}}}_{\rm gc}^{-1} \mathcal{H} = \mathcal{H}_{0{\rm gc}} + \epsilon\mathcal{H}_{1{\rm gc}} + \epsilon^{2}\mathcal{H}_{2{\rm gc}} + \cdots, \end{equation}

where

(2.7)\begin{gather} \mathcal{H}_{1{\rm gc}} = \mathcal{H}_{1} - \boldsymbol{\mathsf{G}}_{1}\boldsymbol{\cdot}\textsf{d}\mathcal{H}_{0}, \end{gather}

(2.8)\begin{gather}\mathcal{H}_{2{\rm gc}} =-\boldsymbol{\mathsf{G}}_{2}\boldsymbol{\cdot}\textsf{d}\mathcal{H}_{0} - \tfrac{1}{2}\boldsymbol{\mathsf{G}}_{1}\boldsymbol{\cdot}\textsf{d} (\mathcal{H}_{1} + \mathcal{H}_{1{\rm gc}}). \end{gather}

Here, we use the formulas $\iota _{n}\boldsymbol {\cdot }\boldsymbol {\omega } \equiv G_{n}^{\alpha }\omega _{\alpha \beta }\,\textsf {d} Z^{\beta }$ (for an arbitrary two-form $\boldsymbol {\omega }$) and $\boldsymbol{\mathsf{G}}_{n}\boldsymbol {\cdot }\textsf {d}\mathcal {K} \equiv G_{n}^{\alpha }\partial \mathcal {K}/\partial Z^{\alpha }$ (for an arbitrary scalar field $\mathcal {K}$), where the summation rule is used. We note that the guiding-centre phase-space transformation considered in the present work will contain all terms associated with first-order space–time derivatives of the electric and magnetic fields, which will require us to consider some terms at third order in $\epsilon$ in (2.3).

In order to construct the extended guiding-centre Lagrangian one-form (2.3), we will need to evaluate the contractions $\iota _{n}\boldsymbol {\cdot }\boldsymbol {\omega }_{0}$ generated by the vector fields $(\boldsymbol{\mathsf{G}}_{1}, \boldsymbol{\mathsf{G}}_{2}, \ldots )$ on the zeroth-order two-form

(2.9)\begin{equation} \boldsymbol{\omega}_{0} = \textsf{d}\gamma_{0}=\frac{e}{c} \left(\frac{\partial A_{j}}{\partial x^{i}}\textsf{d} x^{i} + \frac{\partial A_{j}}{\partial t}\textsf{d} t\right) \wedge\textsf{d} x^{j} - \textsf{d} w\wedge \textsf{d} t, \end{equation}

so that we obtain the $n$th-order expression

(2.10)\begin{equation} \iota_{n}\boldsymbol{\cdot}\boldsymbol{\omega}_{0} = \frac{e}{c} \boldsymbol{B}\boldsymbol{\times} G_{n}^{\boldsymbol{x}}\boldsymbol{\cdot}\textsf{d}\boldsymbol{X} - \left(\frac{e}{c}G_{n}^{\boldsymbol{x}}\boldsymbol{\cdot}\frac{\partial \boldsymbol{A}}{\partial t} + G_{n}^{w}\right)\textsf{d} t, \end{equation}

where $G_{n}^{\boldsymbol {x}}$ and $G_{n}^{w}$ denote the spatial and energy components of the $n$th-order generating vector field $\boldsymbol{\mathsf{G}}_{n}$. Similarly, we will need to evaluate the contractions $\iota _{n}\boldsymbol {\cdot }\boldsymbol {\omega }_{1}$ generated by the vector fields $(\boldsymbol{\mathsf{G}}_{1}, \boldsymbol{\mathsf{G}}_{2}, \ldots )$ on the first-order two-form $\boldsymbol {\omega }_{1} = \textsf {d}\gamma _{1}$, which yields the $(n + 1)$th-order expression

(2.11)\begin{equation} \iota_{n}\boldsymbol{\cdot}\boldsymbol{\omega}_{1} \equiv D_{n}(\,\boldsymbol{p}_{0}) \boldsymbol{\cdot}\textsf{d}\boldsymbol{X} - G_{n}^{\boldsymbol{x}}\boldsymbol{\cdot} \left(\frac{\partial \boldsymbol{p}_{0}}{\partial t}\textsf{d} t + \hat{\textsf{b}}\,\textsf{d} p_{\|0} + \frac{\partial \boldsymbol{q}_{\bot0}}{\partial J_{0}}\textsf{d}\kern0.06em J_{0} + \frac{\partial \boldsymbol{q}_{\bot0}}{\partial \theta_{0}}\textsf{d}\theta_{0}\right), \end{equation}

where $\boldsymbol {p}_{0}$ is defined in (1.3) and the spatial components are expressed in terms of the operator (Tronko & Brizard 2015)

(2.12)\begin{equation} D_{n}(\boldsymbol{C}) \equiv \left( G_{n}^{p_{\|}}\frac{\partial \boldsymbol{C}}{\partial p_{\|0}} + G_{n}^{J} \frac{\partial \boldsymbol{C}}{\partial J_{0}} + G_{n}^{\theta}\frac{\partial \boldsymbol{C}}{\partial \theta_{0}}\right) - G_{n}^{\boldsymbol{x}}\boldsymbol{\times}\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{C}, \end{equation}

where $\boldsymbol {C}$ is an arbitrary vector function on the guiding-centre phase space. In what follows, unless it is necessary, we will omit writing the subscript $0$ on local particle phase-space coordinates (i.e. $p_{\|0}$ is written as $p_{\|}$).

The purpose of the Lie-transform expressions (2.3) and (2.6) is to construct a gyroangle-independent extended guiding-centre phase-space Lagrangian one-form $\varGamma _{\rm gc}$ and extended guiding-centre Hamiltonian $\mathcal {H}_{\rm gc}$ in terms of which guiding-centre equations are derived. The generic forms considered for the extended guiding-centre phase-space Lagrangian one-form (2.3) and the extended guiding-centre Hamiltonian (2.6) in the present work are

(2.13)\begin{gather} \varGamma_{\rm gc} \equiv \frac{e}{c}\boldsymbol{A}\boldsymbol{\cdot}\textsf{d}\boldsymbol{X} - W\,\textsf{d} t + \epsilon\boldsymbol{P}_{\rm gc}\boldsymbol{\cdot}\textsf{d}\boldsymbol{X} + \epsilon^{2}J(\textsf{d}\theta - \boldsymbol{R}\boldsymbol{\cdot}\textsf{d}\boldsymbol{X} - \mathcal{S}\,\textsf{d} t), \end{gather}

(2.14)\begin{gather}\mathcal{H}_{\rm gc} \equiv e\varPhi + \epsilon K_{\rm gc} - W, \end{gather}

where the gyroangle-independent guiding-centre kinetic momentum $\boldsymbol {P}_{\rm gc}$ and guiding- centre kinetic energy $K_{\rm gc}$ are expressed as asymptotic series in powers of $\epsilon$, while the vector field $\boldsymbol {R} \equiv \boldsymbol {\nabla }\hat {\bot }\boldsymbol {\cdot }\hat {\rho }$ and the scalar field $\mathcal {S} \equiv (\partial \hat {\bot }/\partial t)\boldsymbol {\cdot }\hat {\rho }$ are required to preserve gyrogauge invariance. See appendix A of the recent paper by Brizard (2023b) for an updated discussion on gyrogauge invariance introduced by Littlejohn (1981, 1983, 1988). We note that the separation of the guiding-centre transformations of the extended guiding-centre phase-space Lagrangian one-form (2.3) and the extended guiding-centre Hamiltonian (2.6) might enable the application of computer algorithms previously used by Burby, Squire & Qin (2013), but this consideration falls outside the scope of the present work.

3. Symplectic polarization guiding-centre theory

In the present work, we will consider the symplectic polarization guiding-centre theory, where the guiding-centre kinetic momentum $\boldsymbol {P}_{\rm gc}$ retains the contribution from the $E\times B$ momentum (1.4), which will then introduce the polarization drift velocity explicitly in the guiding-centre equations of motion.

In the alternate Hamiltonian polarization guiding-centre theory, on the other hand, $E\times B$ momentum (1.4) is removed from the guiding-centre kinetic momentum and polarization effects enter solely through the guiding-centre Hamiltonian. The Hamiltonian case generally requires a different ordering for the electric field (because it produces a polarization displacement that must be compared with the characteristic lowest-order particle gyroradius) and will be considered in a future publication. This dual representation is analogous to the treatment of the perturbed magnetic field in nonlinear gyrokinetic theory (Brizard & Hahm 2007).

The reader interested in results of the guiding-centre Lie-transform perturbation analysis can skip § 3 and go to § 4, where we present the extended guiding-centre Hamiltonian structure as well as the (regular) guiding-centre Lagrangian, from which guiding-centre polarization and magnetization can be derived by functional derivatives with respect to the electric field $\boldsymbol {E}$ and magnetic field $\boldsymbol {B}$, respectively.

3.1. First-order perturbation analysis

3.1.1. First-order symplectic structure

We begin our perturbation analysis by considering the first-order guiding-centre symplectic one-form (2.4), which is now explicitly written in the symplectic polarization representation as

(3.1)\begin{align} \varGamma_{1{\rm gc}} & = (P_{\|}\hat{\textsf{b}} + \boldsymbol{P}_{E} + \boldsymbol{q}_{\bot}) \boldsymbol{\cdot}\textsf{d}\boldsymbol{X} - \frac{e}{c}\boldsymbol{B}\boldsymbol{\times} G_{1}^{\boldsymbol{x}} \boldsymbol{\cdot}\textsf{d}\boldsymbol{X} + \left(\frac{e}{c}G_{1}^{\boldsymbol{x}}\boldsymbol{\cdot} \frac{\partial \boldsymbol{A}}{\partial t} + G_{1}^{w}\right)\textsf{d} t \nonumber\\ & = (P_{\|}\hat{\textsf{b}} + \boldsymbol{P}_{E}) \boldsymbol{\cdot}\textsf{d}\boldsymbol{X} \equiv \boldsymbol{P}_{0}\boldsymbol{\cdot}\textsf{d}\boldsymbol{X}, \end{align}

where the first-order gauge scalar field $\sigma _{1}$ is not needed and the gyroangle-dependent terms on the right of (3.1) are eliminated by selecting

(3.2)\begin{equation} G_{1}^{\boldsymbol{x}} = \boldsymbol{q}_{\bot}\boldsymbol{\times}\frac{c\hat{\textsf{b}}}{eB} = \frac{1}{m\varOmega}\frac{\partial \boldsymbol{q}_{\bot}}{\partial \theta} \equiv-\boldsymbol{\rho}_{0}, \end{equation}

and

(3.3)\begin{equation} G_{1}^{w} =-\frac{e}{c}G_{1}^{\boldsymbol{x}}\boldsymbol{\cdot}\frac{\partial \boldsymbol{A}}{\partial t} = \frac{e}{c}\frac{\partial \boldsymbol{A}}{\partial t}\boldsymbol{\cdot}\boldsymbol{\rho}_{0}, \end{equation}

where the gyroangle-dependent vector $\boldsymbol {\rho }_{0}$ represents the lowest-order particle gyroradius.

With $G_{1}^{\boldsymbol {x}}$ defined by (3.2), the resulting first-order guiding-centre phase-space Lagrangian $\varGamma _{1{\rm gc}} = \boldsymbol {P}_{0}\boldsymbol {\cdot }\textsf {d}\boldsymbol {X}$ yields the $(n+1)$th-order contraction

(3.4)\begin{equation} \iota_{n}\boldsymbol{\cdot}\boldsymbol{\omega}_{{\rm gc}1} \equiv D_{n}(\boldsymbol{P}_{0}) \boldsymbol{\cdot}\textsf{d}\boldsymbol{X} - G_{n}^{\boldsymbol{x}}\boldsymbol{\cdot} \left(\hat{\textsf{b}}\,\textsf{d} P_{\|} + \frac{\partial \boldsymbol{P}_{0}}{\partial t}\textsf{d} t\right), \end{equation}

where, using the operator (2.12), the spatial components in (3.4) are

(3.5)\begin{equation} D_{n}(\boldsymbol{P}_{0})= G_{n}^{p_{\|}}\hat{\textsf{b}} - G_{n}^{\boldsymbol{x}} \boldsymbol{\times}\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{P}_{0}, \end{equation}

which contain gyroangle-dependent and -independent contributions.

3.1.2. First-order Hamiltonian

The first-order guiding-centre Hamiltonian is determined from (2.7) as

(3.6)\begin{align} \mathcal{H}_{1{\rm gc}} & = \mu B + \frac{|\boldsymbol{P}_{0}|^{2}}{2m} + \boldsymbol{P}_{E}\boldsymbol{\cdot}\frac{\boldsymbol{q}_{\bot}}{m} - eG_{1}^{\boldsymbol{x}}\boldsymbol{\cdot}\boldsymbol{\nabla}\varPhi + G_{1}^{w} \nonumber\\ & = \mu B + \frac{|\boldsymbol{P}_{0}|^{2}}{2m} + \boldsymbol{P}_{E}\boldsymbol{\cdot} \frac{\boldsymbol{q}_{\bot}}{m} - e\boldsymbol{\rho}_{0}\boldsymbol{\cdot}\boldsymbol{E}, \end{align}

where the components (3.2)–(3.3) were substituted, and the electric field is defined as $\boldsymbol {E} = -\boldsymbol {\nabla }\varPhi - c^{-1}\partial \boldsymbol {A}/\partial t$. By using the identity

(3.7)\begin{equation} \boldsymbol{P}_{E}\boldsymbol{\cdot}\boldsymbol{q}_{\bot}/m = \boldsymbol{E}\boldsymbol{\times} \frac{c\hat{\textsf{b}}}{B}\boldsymbol{\cdot}\boldsymbol{q}_{\bot} = e\boldsymbol{\rho}_{0}\boldsymbol{\cdot}\boldsymbol{E}, \end{equation}

the first-order guiding-centre Hamiltonian is automatically gyroangle independent

(3.8)\begin{equation} \mathcal{H}_{1{\rm gc}} = \mu B + |\boldsymbol{P}_{0}|^{2}/(2m) = \mu B + P_{\|}^{2}/(2m) + (m/2)|\boldsymbol{E}\boldsymbol{\times} c\hat{\textsf{b}}/B|^{2}, \end{equation}

which corresponds to the kinetic energy in the frame drifting with the $E\times B$ velocity.

3.2. Second-order perturbation analysis

3.2.1. Second-order symplectic structure

We now proceed with the second-order guiding-centre symplectic one-form (2.5), which is explicitly expressed as

(3.9)\begin{align} \varGamma_{2{\rm gc}} & =-\left[\frac{e}{c}\boldsymbol{B}\boldsymbol{\times} G_{2}^{\boldsymbol{x}} + D_{1}(\boldsymbol{P}_{2})\right] \boldsymbol{\cdot}\textsf{d}\boldsymbol{X} + \frac{1}{2}G_{1}^{\boldsymbol{x}}\boldsymbol{\cdot}\left(\frac{\partial \boldsymbol{q}_{\bot}}{\partial J} \textsf{d} J + \frac{\partial \boldsymbol{q}_{\bot}}{\partial \theta}\textsf{d}\theta \right) \nonumber\\ & \quad +\left(G_{1}^{\boldsymbol{x}}\boldsymbol{\cdot}\frac{\partial \boldsymbol{P}_{2}}{\partial t} + G_{2}^{w} + \frac{e}{c}G_{2}^{\boldsymbol{x}}\boldsymbol{\cdot}\frac{\partial \boldsymbol{A}}{\partial t}\right)\textsf{d} t \nonumber\\ & \equiv \boldsymbol{\varPi}_{1}\boldsymbol{\cdot}\textsf{d}\boldsymbol{X} + J\left(\textsf{d}\theta -\boldsymbol{R}\boldsymbol{\cdot}\textsf{d}\boldsymbol{X} - \mathcal{S}\,\textsf{d} t\right), \end{align}

where the second-order gauge scalar field $\sigma _{2}$ is not needed, and, using the definition (2.12) with $\boldsymbol {P}_{2} \equiv \boldsymbol {P}_{0} + \frac {1}{2}\boldsymbol {q}_{\bot }$, we find

(3.10)\begin{equation} D_{1}(\boldsymbol{P}_{2}) = D_{1}(\boldsymbol{P}_{0}) + \frac{1}{2} \left(g_{1}^{J}\frac{\partial \boldsymbol{q}_{\bot}}{\partial J} + g_{1}^{\theta} \frac{\partial \boldsymbol{q}_{\bot}}{\partial \theta} \right) + J \left[\boldsymbol{R} - \left(\frac{1}{2}\tau + \alpha_{1}\right)\hat{\textsf{b}} \right], \end{equation}

where $g_{1}^{J} \equiv G_{1}^{J} - J\boldsymbol {\rho }_{0}\boldsymbol {\cdot }\boldsymbol {\nabla }\ln B$ and $g_{1}^{\theta } \equiv G_{1}^{\theta } + \boldsymbol {\rho }_{0}\boldsymbol {\cdot }\boldsymbol {R}$, while $\tau \equiv \hat {\textsf {b}}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {\times }\hat {\textsf {b}}$ and $\alpha _{1} \equiv \mathsfbit {a}_{1}:\boldsymbol {\nabla }\hat {\textsf {b}}$ is defined in terms of the gyroangle-dependent dyadic tensor $\boldsymbol{\mathsf{a}}_{1} \equiv -\frac {1}{2}(\hat {\bot }\hat {\rho } + \hat {\rho }\hat {\bot })$ (Tronko & Brizard 2015). The first-order guiding-centre symplectic momentum $\boldsymbol {\varPi }_{1}$ in (3.9), which is assumed to be gyroangle independent, will be determined based on the consistency of the Lie-transform perturbation analysis at the third order (see § 3.3) as well as the guiding-centre push-forward derivation of the guiding-centre polarization in the absence of a background electric field (Brizard 2013; Tronko & Brizard 2015).

Substituting these expressions into (3.9), we obtain the vector equation

(3.11)\begin{equation} \boldsymbol{\varPi}_{1} =-\frac{e}{c}\boldsymbol{B}\boldsymbol{\times} G_{2}^{\boldsymbol{x}} - g_{1}^{p_{\|}}\hat{\textsf{b}} - \boldsymbol{\rho}_{0}\boldsymbol{\times}\boldsymbol{\nabla}\boldsymbol{\times} \boldsymbol{P}_{0} - \frac{1}{2} \left(g_{1}^{J}\frac{\partial \boldsymbol{q}_{\bot}}{\partial J} + g_{1}^{\theta} \frac{\partial \boldsymbol{q}_{\bot}}{\partial \theta}\right), \end{equation}

where $g_{1}^{p_{\|}} \equiv G_{1}^{p_{\|}} - J(\frac {1}{2}\tau + \alpha _{1})$, while we choose the second-order energy component

(3.12)\begin{equation} G_{2}^{w} = \boldsymbol{\rho}_{0}\boldsymbol{\cdot}\frac{\partial \boldsymbol{P}_{0}}{\partial t} - \frac{e}{c}\frac{\partial \boldsymbol{A}}{\partial t}\boldsymbol{\cdot} G_{2}^{\boldsymbol{x}}, \end{equation}

where we used

(3.13)\begin{equation} G_{1}^{\boldsymbol{x}}\boldsymbol{\cdot}\frac{\partial \boldsymbol{P}_{2}}{\partial t} =-\boldsymbol{\rho}_{0}\boldsymbol{\cdot}\frac{\partial \boldsymbol{P}_{2}}{\partial t} =-\boldsymbol{\rho}_{0}\boldsymbol{\cdot}\frac{\partial \boldsymbol{P}_{0}}{\partial t} - J\mathcal{S}. \end{equation}

Next, from the parallel component of (3.11), we obtain the first-order component

(3.14)\begin{equation} G_{1}^{p_{\|}} = \frac{\partial \boldsymbol{\rho}_{0}}{\partial \theta}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{\times}\boldsymbol{P}_{0} + J \left(\frac{1}{2}\tau + \alpha_{1}\right) - \varPi_{1\|}, \end{equation}

where $\varPi _{1\|} \equiv \hat {\textsf {b}}\boldsymbol {\cdot }\boldsymbol {\varPi }_{1}$ is the parallel component of the first-order symplectic momentum $\boldsymbol {\varPi }_{1}$. From the perpendicular components of (3.11), on the other hand, we find

(3.15)\begin{equation} G_{2}^{\boldsymbol{x}} = G_{2\|}^{\boldsymbol{x}}\hat{\textsf{b}} + \left(\frac{\hat{\textsf{b}}}{m\varOmega}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{P}_{0}\right) \boldsymbol{\rho}_{0} + \frac{1}{2} \left(g_{1}^{J}\frac{\partial \boldsymbol{\rho}_{0}}{\partial J} + g_{1}^{\theta} \frac{\partial \boldsymbol{\rho}_{0}}{\partial \theta}\right) - \boldsymbol{\varPi}_{1}\boldsymbol{\times}\frac{\hat{\textsf{b}}}{m\varOmega}, \end{equation}

where $G_{2\|}^{\boldsymbol {x}} \equiv \hat {\textsf {b}}\boldsymbol {\cdot } G_{2}^{\boldsymbol {x}}$. The interpretation of the first-order guiding-centre symplectic momentum $\boldsymbol {\varPi }_{1}$ will be given in § 3.3.2.

3.2.2. Second-order Hamiltonian

The second-order guiding-centre Hamiltonian is determined from (2.8) as

(3.16)\begin{align} \mathcal{H}_{2{\rm gc}} & = e\boldsymbol{E}\boldsymbol{\cdot} \left(G_{2}^{\boldsymbol{x}} - \frac{1}{2}\boldsymbol{\mathsf{G}}_{1} \boldsymbol{\cdot}\textsf{d}\boldsymbol{\rho}_{0}\right) + \frac{e}{2}\boldsymbol{\rho}_{0}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{E}\boldsymbol{\cdot}\boldsymbol{\rho}_{0} - \frac{P_{\|}}{m}G_{1}^{p_{\|}} - \varOmega G_{1}^{J} \nonumber\\ & \quad +\boldsymbol{\rho}_{0}\boldsymbol{\cdot}\left(\mu\boldsymbol{\nabla} B + \boldsymbol{\nabla}\boldsymbol{u}_{E}\boldsymbol{\cdot} m\boldsymbol{u}_{E} + \frac{\partial \boldsymbol{P}_{0}}{\partial t}\right), \end{align}

where $G_{1}^{p_{\|}}$ is given by (3.14).

Next, we introduce the particle gyroradius

(3.17)\begin{equation} \boldsymbol{\rho} \equiv \boldsymbol{x} - {\boldsymbol{\mathsf{T}}}_{\rm gc} \boldsymbol{X} = \epsilon\boldsymbol{\rho}_{0} - \epsilon^{2} \left(G_{2}^{\boldsymbol{x}} - \tfrac{1}{2}\boldsymbol{\mathsf{G}}_{1} \boldsymbol{\cdot}\textsf{d} \boldsymbol{\rho}_{0} \right) + \cdots = \epsilon \boldsymbol{\rho}_{0} + \epsilon^{2}\boldsymbol{\rho}_{1} + \cdots, \end{equation}

which is defined as the difference between the particle position $\boldsymbol {x}$ and the pull-back $\boldsymbol{\mathsf{T}}_{\rm gc}\boldsymbol {X}$ of the guiding-centre position $\boldsymbol {X}$, where $\boldsymbol {\rho }_{1} \equiv -G_{2}^{\boldsymbol {x}} + \frac {1}{2}\boldsymbol{\mathsf{G}}_{1}\boldsymbol {\cdot }\textsf {d}\boldsymbol {\rho }_{0}$ is the first-order particle gyroradius, where

(3.18)\begin{equation} \boldsymbol{\mathsf{G}}_{1}\boldsymbol{\cdot}\textsf{d}\boldsymbol{\rho}_{0} = g_{1}^{J}\frac{\partial \boldsymbol{\rho}_{0}}{\partial J} + g_{1}^{\theta}\frac{\partial \boldsymbol{\rho}_{0}}{\partial \theta} + \frac{J}{m\varOmega} \left(\boldsymbol{\nabla}\boldsymbol{\cdot}\hat{\textsf{b}} - 4\alpha_{2}\right)\hat{\textsf{b}}. \end{equation}

Hence, using (3.15) and (3.18), we obtain the first-order particle gyroradius vector

(3.19)\begin{align} \boldsymbol{\rho}_{1} & = \boldsymbol{\varPi}_{1}\boldsymbol{\times} \frac{\hat{\textsf{b}}}{m\varOmega} - \left[G_{2\|}^{\boldsymbol{x}} - \frac{J}{m\varOmega} \left(\frac{1}{2}\boldsymbol{\nabla} \boldsymbol{\cdot}\hat{\textsf{b}} - 2\alpha_{2}\right)\right]\hat{\textsf{b}} \nonumber\\ & \quad +\left(\frac{1}{2}\boldsymbol{\rho}_{0}\boldsymbol{\cdot}\boldsymbol{\nabla}\ln B - \frac{\hat{\textsf{b}}}{m\varOmega}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{\times}\boldsymbol{P}_{0}\right)\boldsymbol{\rho}_{0}, \end{align}

where $4\alpha _{2} \equiv -\partial \alpha _{1}/\partial \theta$.

We now note that the gyroangle-dependent part $\tilde {G}_{1}^{J} \equiv G_{1}^{J} - \langle G_{1}^{J}\rangle$ can be defined such that the right side of (3.16) only contains terms that are gyroangle independent. Hence, we find

(3.20)\begin{align} \varOmega\tilde{G}_{1}^{J} & =-\frac{p_{\|}}{m}\tilde{G}_{1}^{p_{\|}} - e\tilde{\boldsymbol{\rho}}_{1}\boldsymbol{\cdot}\boldsymbol{E} - \frac{2J}{m\varOmega}e {\boldsymbol{\mathsf{a}}}_{2}:\boldsymbol{\nabla}\boldsymbol{E} \nonumber\\ & \quad +\boldsymbol{\rho}_{0}\boldsymbol{\cdot}\left(\mu\boldsymbol{\nabla} B + \boldsymbol{\nabla}\boldsymbol{P}_{E} \boldsymbol{\cdot}\boldsymbol{u}_{E} + \frac{\partial \boldsymbol{P}_{0}}{\partial t}\right), \end{align}

where $\boldsymbol{\mathsf{a}}_{2} \equiv \frac {1}{4}(\hat {\bot }\hat {\bot } - \hat {\rho }\hat {\rho }) = -\frac {1}{4}\partial \boldsymbol{\mathsf{a}}_{1}/\partial \theta$ and the gyroangle-dependent part of (3.14) is

(3.21)\begin{equation} \tilde{G}_{1}^{p_{\|}} \equiv G_{1}^{p_{\|}} - \langle G_{1}^{p_{\|}}\rangle = \frac{\partial \boldsymbol{\rho}_{0}}{\partial \theta} \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{P}_{0} + J\alpha_{1}. \end{equation}

The second-order guiding-centre Hamiltonian is, thus, defined as

(3.22)\begin{align} \mathcal{H}_{2{\rm gc}} & = \frac{P_{\|}}{m}\left(\varPi_{1\|} - \frac{1}{2}J\tau \right) - \varOmega \langle G_{1}^{J}\rangle + \frac{Jc}{2B}(\mathbb{I} - \hat{\textsf{b}}\hat{\textsf{b}}):\boldsymbol{\nabla}\boldsymbol{E} - \langle\boldsymbol{\rho}_{1}\rangle\boldsymbol{\cdot} e\boldsymbol{E} \nonumber\\ & = \frac{P_{\|}}{m}\left(\varPi_{1\|} - \frac{1}{2}J\tau\right) - \varOmega\langle G_{1}^{J}\rangle + \boldsymbol{\nabla}\boldsymbol{\cdot} \left(\frac{e}{2}\langle\boldsymbol{\rho}_{0}\boldsymbol{\rho}_{0}\rangle \boldsymbol{\cdot}\boldsymbol{E}\right) + \left(\boldsymbol{\varPi}_{1} - \frac{1}{2}J\boldsymbol{\nabla}\boldsymbol{\times}\hat{\textsf{b}}\right)\boldsymbol{\cdot}\boldsymbol{u}_{E}, \end{align}

where $\langle G_{1}^{p_{\|}}\rangle = \frac {1}{2}J\tau - \varPi _{1\|}$ and the gyroangle-averaged first-order particle gyroradius is

(3.23)\begin{align} \langle\boldsymbol{\rho}_{1}\rangle & = \frac{J}{2m\varOmega} [(\mathbb{I} - \hat{\textsf{b}}\hat{\textsf{b}})\boldsymbol{\cdot}\boldsymbol{\nabla}\ln B + (\boldsymbol{\nabla}\boldsymbol{\cdot}\hat{\textsf{b}})\hat{\textsf{b}}] + \boldsymbol{\varPi}_{1}\boldsymbol{\times}\frac{\hat{\textsf{b}}}{m\varOmega} \nonumber\\ & \equiv-\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{1}{2}\langle\boldsymbol{\rho}_{0}\boldsymbol{\rho}_{0}\rangle\right) + \left(\boldsymbol{\varPi}_{1} - \frac{J}{2}\boldsymbol{\nabla}\boldsymbol{\times}\hat{\textsf{b}}\right)\boldsymbol{\times}\frac{\hat{\textsf{b}}}{m\varOmega}, \end{align}

where $\langle \boldsymbol {\rho }_{0}\boldsymbol {\rho }_{0}\rangle = (J/m\varOmega )(\mathbb {I} - \hat {\textsf {b}}\hat {\textsf {b}})$. Hence, we now need expressions for $\langle G_{1}^{J}\rangle$ and $\boldsymbol {\varPi }_{1}$ in order to obtain an explicit expression for the second-order guiding-centre Hamiltonian (3.22), which are determined at third order in our perturbation analysis.

3.3. Third-order perturbation analysis

The third-order term in the guiding-centre phase-space Lagrangian one-form (2.3) is expressed as

(3.24)\begin{equation} \varGamma_{3{\rm gc}} =-\iota_{3}\boldsymbol{\cdot}\boldsymbol{\omega}_{0} - \iota_{2}\boldsymbol{\cdot}\boldsymbol{\omega}_{{\rm gc}1} + \frac{\iota_{1}}{3} \cdot\textsf{d}\left(\iota_{1}\boldsymbol{\cdot}\boldsymbol{\omega}_{1} + \frac{\iota_{1}}{2}\boldsymbol{\cdot}\boldsymbol{\omega}_{{\rm gc}1}\right) + \textsf{d} \sigma_{3}. \end{equation}

In what follows, the gauge function $\sigma _{3}$ will play an important role in completing the guiding-centre phase-space transformation, while the third-order guiding-centre Hamiltonian will not be needed in the present guiding-centre formulation.

3.3.1. Third-order symplectic structure

The remaining components $(G_{2\|}^{\boldsymbol {x}}, \langle G_{1}^{J}\rangle, G_{1}^{\theta })$ and the first-order guiding-centre momentum $\boldsymbol {\varPi }_{1}$ will now be determined from the momentum components of the third-order guiding-centre symplectic one-form (3.24)

(3.25)\begin{align} \varGamma_{3\boldsymbol{p}} & \equiv \left[G_{2\|}^{\boldsymbol{x}} + \frac{\partial D_{1}(\boldsymbol{P}_{3})}{\partial p_{\|}} \boldsymbol{\cdot}\boldsymbol{\rho}_{0} + \frac{\partial \sigma_{3}}{\partial p_{\|}} \right]\textsf{d} p_{\|} + \left[\frac{2}{3}G_{1}^{\theta} + \frac{\partial D_{1}(\boldsymbol{P}_{3})}{\partial J} \boldsymbol{\cdot}\boldsymbol{\rho}_{0} + \frac{\partial \sigma_{3}}{\partial J}\right]\textsf{d} J \nonumber\\ & \quad +\left[-\frac{2}{3}G_{1}^{J} + \frac{\partial D_{1}(\boldsymbol{P}_{3})}{\partial \theta} \boldsymbol{\cdot}\boldsymbol{\rho}_{0} + \frac{\partial \sigma_{3}}{\partial \theta}\right]\textsf{d} \theta, \end{align}

where $\boldsymbol {P}_{3} \equiv \frac {1}{2}\boldsymbol {P}_{0} + \frac {1}{3}\boldsymbol {q}_{\bot }$, so that

(3.26)\begin{equation} D_{1}(\boldsymbol{P}_{3}) = \frac{1}{2}G_{1}^{p_{\|}}\hat{\textsf{b}} + \frac{1}{3} \left(G_{1}^{J}\frac{\partial \boldsymbol{q}_{\bot}}{\partial J} + G_{1}^{\theta} \frac{\partial \boldsymbol{q}_{\bot}}{\partial \theta}\right) + \boldsymbol{\rho}_{0}\boldsymbol{\times}\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{P}_{3}. \end{equation}

Since $\partial \boldsymbol {\rho }_{0}/\partial p_{\|} = 0$, the $p_{\|}$-component of (3.25) suggests that we define the new gauge function

(3.27)\begin{equation} \bar{\sigma}_{3} \equiv \sigma_{3} + D_{1}(\boldsymbol{P}_{3})\boldsymbol{\cdot}\boldsymbol{\rho}_{0} = \sigma_{3} - \frac{2}{3}JG_{1}^{\theta}, \end{equation}

where the last expression follows from (3.26). Using the new gauge function (3.27), the momentum components (3.25), therefore, become

(3.28)\begin{align} \varGamma_{3\boldsymbol{p}} & = \left(G_{2\|}^{\boldsymbol{x}} + \frac{\partial \bar{\sigma}_{3}}{\partial p_{\|}}\right)\textsf{d} p_{\|} + \left(G_{1}^{\theta} +\frac{\partial \bar{\sigma}_{3}}{\partial J}\right)\textsf{d} J \nonumber\\ & \quad +\left(\frac{\partial \bar{\bar{\sigma}}_{3}}{\partial \theta} - G_{1}^{J} + \frac{J\hat{\textsf{b}}}{m\varOmega}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{\times}\boldsymbol{P}_{0}\right)\textsf{d} \theta, \end{align}

where, using (3.26), we introduced the identities

(3.29)\begin{equation} \left.\begin{gathered} D_{1}(\boldsymbol{P}_{3})\boldsymbol{\cdot}\frac{\partial \boldsymbol{\rho}_{0}}{\partial J} \equiv-\frac{1}{3}G_{1}^{\theta}, \\ D_{1}(\boldsymbol{P}_{3})\boldsymbol{\cdot}\frac{\partial \boldsymbol{\rho}_{0}}{\partial \theta} \equiv \frac{1}{3}G_{1}^{J} + \frac{2J\hat{\textsf{b}}}{m\varOmega}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{P}_{3}, \end{gathered}\right\} \end{equation}

so that we can also introduce yet another gauge function

(3.30)\begin{equation} \bar{\bar{\sigma}}_{3} \equiv \bar{\sigma}_{3} - \frac{1}{3} \left(2J\boldsymbol{\rho}_{0}\boldsymbol{\cdot}\boldsymbol{R} + J \frac{\partial \boldsymbol{\rho}_{0}}{\partial \theta}\boldsymbol{\cdot}\boldsymbol{\nabla}\ln B\right) \end{equation}

in the $\theta$-component of (3.25). By requiring that the momentum components (3.28) vanish, we now obtain the definitions

(3.31)\begin{gather} G_{1}^{J} \equiv-\frac{J\hat{\textsf{b}}}{m\varOmega}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{P}_{0} + \frac{\partial \bar{\bar{\sigma}}_{3}}{\partial \theta}, \end{gather}

(3.32)\begin{gather}G_{2\|}^{\boldsymbol{x}} \equiv-\frac{\partial \bar{\sigma}_{3}}{\partial p_{\|}}, \end{gather}

(3.33)\begin{gather}G_{1}^{\theta} \equiv-\frac{\partial \bar{\sigma}_{3}}{\partial J}. \end{gather}

Hence, the components $G_{2\|}^{\boldsymbol {x}}$ and $G_{1}^{\theta }$ are determined from the third-order gauge function $\bar {\sigma }_{3}$, which is determined from (3.30), while $\bar {\bar {\sigma }}_{3}$ is determined from the gyroangle-dependent part $\tilde {G}_{1}^{J} \equiv \partial \bar {\bar {\sigma }}_{3}/\partial \theta$ obtained from (3.20). Since these gyroangle-dependent components are not needed in what follows, however, they will not be derived here.

3.3.2. Second-order guiding-centre Hamiltonian

From (3.31), we immediately conclude that $\langle G_{1}^{J}\rangle$ must be defined as

(3.34)\begin{equation} \langle G_{1}^{J}\rangle \equiv-\frac{J\hat{\textsf{b}}}{m\varOmega} \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{P}_{0}, \end{equation}

which was obtained in previous derivations (Brizard 1995; Madsen 2010; Frei et al. 2020), so that the second-order guiding-centre Hamiltonian (3.22) becomes

(3.35)\begin{equation} \mathcal{H}_{2{\rm gc}} =-\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{e}{2} \langle\boldsymbol{\rho}_{0}\boldsymbol{\rho}_{0}\rangle\boldsymbol{\cdot}\boldsymbol{E}\right) + \left(\boldsymbol{\varPi}_{1} + \frac{1}{2}J\boldsymbol{\nabla}\boldsymbol{\times}\hat{\textsf{b}}\right) \boldsymbol{\cdot}\frac{\boldsymbol{P}_{0}}{m}, \end{equation}

where we used the identity

(3.36)\begin{equation} \frac{J}{2m\varOmega}(\mathbb{I} - \hat{\textsf{b}}\hat{\textsf{b}})\boldsymbol{\cdot} e\boldsymbol{E} \equiv \frac{e}{2}\langle\boldsymbol{\rho}_{0}\boldsymbol{\rho}_{0}\rangle\boldsymbol{\cdot} \boldsymbol{E} \equiv \frac{1}{2}J\hat{\textsf{b}}\boldsymbol{\times}\boldsymbol{u}_{E}. \end{equation}

In previous works (Littlejohn 1981; Hahm 1996; Miyato et al. 2009; Madsen 2010; Frei et al. 2020), the choice for the first-order symplectic momentum $\boldsymbol {\varPi }_{1} = -\frac {1}{2}J\tau \hat {\textsf {b}}$ was used to eliminate the Baños drift (Baños 1967; Northrop & Rome 1978) from the guiding-centre velocity (i.e. $\partial \mathcal {H}_{2\textrm {gc}}/\partial P_{\|} = 0$), which is instead added to the definition of the guiding-centre parallel momentum (3.14). This choice, therefore, yields the second-order guiding-centre Hamiltonian $\mathcal {H}_{2\textrm {gc}} = (J\hat {\textsf {b}}/2)\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {\times }\boldsymbol {u}_{E}$.

A different choice adopted by Tronko & Brizard (2015) for the first-order symplectic momentum

(3.37)\begin{equation} \boldsymbol{\varPi}_{1{\rm pol}} \equiv-\frac{1}{2}J\boldsymbol{\nabla}\boldsymbol{\times}\hat{\textsf{b}}, \end{equation}

on the other hand, was previously shown to yield an exact Lie-transform perturbation derivation of the standard guiding-centre polarization derived by Kaufman (1986) in the absence of an electric field. While this choice still eliminates the Baños drift from the guiding-centre velocity, it also yields an expression for the second-order guiding-centre Hamiltonian (3.35) that exactly represents the guiding-centre finite-Larmor-radius (FLR) correction to the electrostatic potential energy $e\varPhi$ (Brizard 2023a):

(3.38)\begin{equation} \mathcal{H}_{2{\rm gc}} =-\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{e}{2} \langle\boldsymbol{\rho}_{0}\boldsymbol{\rho}_{0}\rangle\boldsymbol{\cdot}\boldsymbol{E}\right). \end{equation}

Finally, we note that recent numerical studies of particle and guiding-centre orbits in axisymmetric magnetic fields (Brizard & Hodgeman 2023) have shown that guiding-centre orbits are faithful (i.e. remain close) to particle orbits only if second-order effects, including the correction (3.37), are included, which confirms earlier results by Belova et al. (2003).

4. Guiding-centre Hamiltonian dynamics

In this section, we summarize the results of the Lie-transform perturbation analysis of the guiding-centre Lagrangian dynamics presented in § 3, and we remove the explicit $\epsilon$ scaling by restoring the physical mass $\epsilon m \rightarrow m$. Hence, we find the guiding-centre phase-space extended one-form

(4.1)\begin{equation} \varGamma_{\rm gc} = \left(\frac{e}{c}\boldsymbol{A} +\boldsymbol{\varPi}_{\rm gc}\right)\boldsymbol{\cdot}\textsf{d}\boldsymbol{X} + J\left(\textsf{d}\theta -\boldsymbol{R}\boldsymbol{\cdot}\textsf{d}\boldsymbol{X} - \mathcal{S}\,\textsf{d} t\right) - W\,\textsf{d} t, \end{equation}

where the guiding-centre symplectic momentum

(4.2)\begin{equation} \boldsymbol{\varPi}_{\rm gc} = P_{\|}\hat{\textsf{b}} + \boldsymbol{P}_{E} - \frac{J}{2}\boldsymbol{\nabla}\boldsymbol{\times}\hat{\textsf{b}} \end{equation}

includes the higher-order polarization correction (3.37). The extended guiding-centre Hamiltonian, on the other hand, is expressed as

(4.3)\begin{equation} \mathcal{H}_{\rm gc} = e\varPhi + K_{\rm gc} - W, \end{equation}

where the guiding-centre kinetic energy in the drifting frame is

(4.4)\begin{equation} K_{\rm gc} = \mu B + \frac{P_{\|}^{2}}{2m} +\frac{m}{2}|\boldsymbol{u}_{E}|^{2} - \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{J\hat{\textsf{b}}}{2}\boldsymbol{\times}\boldsymbol{u}_{E}\right), \end{equation}

which includes the FLR correction (3.38) to the electrostatic potential energy $e\varPhi$. We note that the presence of the gyrogauge fields $(\mathcal {S},\boldsymbol {R})$ in (4.1) guarantees gyrogauge invariance of the guiding-centre equations of motion derived from them.

4.1. Extended guiding-centre Poisson bracket

The extended guiding-centre Poisson bracket $\{\,,\,\}_\textrm {gc}$ is obtained by, first, constructing an 8$\times$8 matrix out of the components of the extended guiding-centre Lagrange two-form $\boldsymbol {\omega }_\textrm {gc} = \textsf {d}\varGamma _\textrm {gc}$ and, then, inverting this matrix to obtain the extended guiding-centre Poisson matrix, whose components are the fundamental brackets $\{Z^{\alpha }, Z^{\beta }\}_\textrm {gc}$. From these components, we obtain the extended guiding-centre Poisson bracket

(4.5)\begin{align} \{\mathcal{F}, \mathcal{G}\}_{\rm gc} & = \left(\frac{\partial \mathcal{F}}{\partial W}\frac{\partial^{*}\mathcal{G}}{\partial t} - \frac{\partial^{*}\mathcal{G}}{\partial t}\frac{\partial \mathcal{G}}{\partial W}\right) + \frac{\boldsymbol{B}^{*}}{B_{\|}^{*}}\boldsymbol{\cdot}\left(\boldsymbol{\nabla}^{*}\mathcal{F}\frac{\partial \mathcal{G}}{\partial P_{\|}} - \frac{\partial \mathcal{F}}{\partial P_{\|}}\boldsymbol{\nabla}^{*}\mathcal{G}\right) \nonumber\\ & \quad -\frac{c\hat{\textsf{b}}}{eB_{\|}^{*}}\boldsymbol{\cdot}\boldsymbol{\nabla}^{*}\mathcal{F}\boldsymbol{\times} \boldsymbol{\nabla}^{*}\mathcal{G} + \left(\frac{\partial \mathcal{F}}{\partial \theta}\frac{\partial \mathcal{G}}{\partial J} - \frac{\partial \mathcal{F}}{\partial J}\frac{\partial \mathcal{G}}{\partial \theta}\right), \end{align}

where

(4.6)\begin{equation} \frac{e}{c}\boldsymbol{B}^{*} = \frac{e}{c}\boldsymbol{B} + \boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{\varPi}_{\rm gc} - J\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{R}, \end{equation}

and the guiding-centre Jacobian is $\mathcal {J}_\textrm {gc} = (e/c)B_{\|}^{*} \equiv (e/c)\hat {\textsf {b}}\boldsymbol {\cdot }\boldsymbol {B}^{*}$. In addition, we have introduced the definitions

(4.7)\begin{gather} \frac{\partial^{*}}{\partial t} \equiv \frac{\partial }{\partial t} + \mathcal{S}\frac{\partial }{\partial \theta}, \end{gather}

(4.8)\begin{gather}\boldsymbol{\nabla}^{*} \equiv \boldsymbol{\nabla} + \boldsymbol{R}^{*}\frac{\partial }{\partial \theta} - \left(\frac{e}{c}\frac{\partial \boldsymbol{A}^{*}}{\partial t} + J\boldsymbol{\nabla}\mathcal{S}\right)\frac{\partial }{\partial W}, \end{gather}

where $\boldsymbol {R}^{*} \equiv \boldsymbol {R} + \frac {1}{2}\boldsymbol {\nabla }\boldsymbol {\times }\hat {\textsf {b}}$. We note that the Poisson bracket (4.5) can be expressed in divergence form as

(4.9)\begin{equation} \{\mathcal{F}, \mathcal{G}\}_{\rm gc} = \frac{1}{B_{\|}^{*}}\frac{\partial }{\partial Z^{\alpha}} (B_{\|}^{*}\mathcal{F}\{Z^{\alpha},\mathcal{G}\}_{\rm gc}), \end{equation}

and that it automatically satisfies the Jacobi identity

(4.10)\begin{equation} \{\mathcal{F},\{\mathcal{G},\mathcal{K}\}_{\rm gc}\}_{\rm gc} + \{\mathcal{G},\{\mathcal{K},\mathcal{F}\}_{\rm gc}\}_{\rm gc} + \{\mathcal{K},\{\mathcal{F},\mathcal{G}\}_{\rm gc}\}_{\rm gc} = 0, \end{equation}

since the extended guiding-centre Lagrange two-form $\boldsymbol {\omega }_\textrm {gc} = \textsf {d}\varGamma _\textrm {gc}$ is exact (i.e. $\textsf {d}\boldsymbol {\omega }_\textrm {gc} = \textsf {d}^{2}\varGamma _\textrm {gc} = 0$) provided $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {B}^{*} = 0$.

Next, we note that the operators (4.7) and (4.8) contain the gyrogauge-invariant combinations $\partial /\partial t + \mathcal {S}\partial /\partial \theta$ and $\boldsymbol {\nabla } + \boldsymbol {R}\partial /\partial \theta$, while (4.6) and (4.8) include the gyrogauge-independent vector fields (Ye & Kaufman 1992; Brizard 2023b)

(4.11)\begin{equation} \left(\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{R}, \boldsymbol{\nabla}\mathcal{S} - \frac{\partial \boldsymbol{R}}{\partial t}\right) = \left(-\frac{1}{2}\epsilon_{ijk}b^{i} \boldsymbol{\nabla} b^{j}\boldsymbol{\times}\boldsymbol{\nabla} b^{k}, -\boldsymbol{\nabla}\hat{\textsf{b}}\boldsymbol{\times}\hat{\textsf{b}}\boldsymbol{\cdot}\frac{\partial \hat{\textsf{b}}}{\partial t}\right), \end{equation}

where $\epsilon _{ijk}$ denotes the completely antisymmetric Levi-Civita tensor.

4.2. Guiding-centre Hamilton equations

The guiding-centre Hamilton equations $\dot {Z}^{\alpha } \equiv \{ Z^{\alpha }, \mathcal {H}_\textrm {gc}\}_\textrm {gc}$ include the guiding-centre velocity

(4.12)\begin{equation} \dot{\boldsymbol{X}} = \frac{P_{\|}}{m}\frac{\boldsymbol{B}^{*}}{B_{\|}^{*}} + \boldsymbol{E}^{*}\boldsymbol{\times}\frac{c\hat{\textsf{b}}}{B_{\|}^{*}}, \end{equation}

where $\hat {\textsf {b}}\boldsymbol {\cdot }\dot {\boldsymbol {X}} = P_{\|}/m$ defines the parallel guiding-centre velocity, and the guiding-centre parallel force

(4.13)\begin{equation} \dot{P}_{\|} = e\boldsymbol{E}^{*}\boldsymbol{\cdot}\frac{\boldsymbol{B}^{*}}{B_{\|}^{*}}, \end{equation}

where the modified electric field is represented as

(4.14)\begin{equation} e\boldsymbol{E}^{*} = e\boldsymbol{E} - \frac{\partial \boldsymbol{\varPi}_{\rm gc}}{\partial t} - \boldsymbol{\nabla} K_{\rm gc} + J \left(\frac{\partial \boldsymbol{R}}{\partial t} - \boldsymbol{\nabla}\mathcal{S}\right), \end{equation}

and the gyroangle angular velocity

(4.15)\begin{equation} \dot{\theta} \equiv \frac{\partial K_{\rm gc}}{\partial J} + \mathcal{S} + \dot{\boldsymbol{X}}\boldsymbol{\cdot}\boldsymbol{R}^{*} = \varOmega - \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{\hat{\textsf{b}}}{2}\boldsymbol{\times}\boldsymbol{u}_{E}\right) + \mathcal{S} + \dot{\boldsymbol{X}}\boldsymbol{\cdot}\boldsymbol{R}^{*}. \end{equation}

We note that the reduced guiding-centre equations of motion (4.12)–(4.13) satisfy the guiding-centre Liouville equation

(4.16)\begin{equation} \frac{\partial B_{\|}^{*}}{\partial t} =-\boldsymbol{\nabla}\boldsymbol{\cdot}(B_{\|}^{*}\dot{\boldsymbol{X}}) - \frac{\partial }{\partial P_{\|}}(B_{\|}^{*}\dot{P}_{\|}), \end{equation}

where

(4.17)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}(B_{\|}^{*}\dot{\boldsymbol{X}}) = \boldsymbol{\nabla}\boldsymbol{\times} \boldsymbol{E}^{*}\boldsymbol{\cdot} c\hat{\textsf{b}} - e\boldsymbol{E}^{*}\boldsymbol{\cdot} \frac{c}{e}\boldsymbol{\nabla}\boldsymbol{\times}\hat{\textsf{b}} =-\hat{\textsf{b}}\boldsymbol{\cdot}\frac{\partial \boldsymbol{B}^{*}}{\partial t} - e\boldsymbol{E}^{*}\boldsymbol{\cdot}\frac{\partial \boldsymbol{B}^{*}}{\partial P_{\|}} \end{equation}

and

(4.18)\begin{equation} \frac{\partial }{\partial P_{\|}}(B_{\|}^{*}\dot{P}_{\|}) = e\boldsymbol{E}^{*}\boldsymbol{\cdot} \frac{\partial \boldsymbol{B}^{*}}{\partial P_{\|}} + \boldsymbol{B}^{*} \boldsymbol{\cdot} e\frac{\partial \boldsymbol{E}^{*}}{\partial P_{\|}} = e\boldsymbol{E}^{*}\boldsymbol{\cdot} \frac{\partial \boldsymbol{B}^{*}}{\partial P_{\|}} - \boldsymbol{B}^{*}\boldsymbol{\cdot} \frac{\partial \hat{\textsf{b}}}{\partial t}, \end{equation}

where we made use of the modified Faraday's law $\partial \boldsymbol {B}^{*}/\partial t = -c\boldsymbol {\nabla }\boldsymbol {\times }\boldsymbol {E}^{*}$.

4.3. Eulerian variations of the guiding-centre Lagrangian

The results of the Lie-transform perturbation analysis carried out in § 3.3 can also be used to construct the following (regular) guiding-centre Lagrangian:

(4.19)\begin{align} L_{\rm gc} & = \left(\frac{e}{c}\boldsymbol{A} + \boldsymbol{\varPi}_{\rm gc} - J\boldsymbol{R}\right)\boldsymbol{\cdot}\dot{\boldsymbol{X}} + J\dot{\theta} - (e\varPhi +K_{\rm gc} + J\mathcal{S}) \nonumber\\ & \equiv (e/c)\boldsymbol{A}^{*}\boldsymbol{\cdot} \dot{\boldsymbol{X}} + J\dot{\theta} - H_{\rm gc}, \end{align}

where $\boldsymbol {\varPi }_\textrm {gc}$ and $K_\textrm {gc}$ are defined in (4.2) and (4.4), respectively. The guiding-centre Euler–Lagrange equations are derived from this Lagrangian as

(4.20)\begin{gather} \dot{P}_{\|}\hat{\textsf{b}} = e\boldsymbol{E}^{*} + (e/c)\dot{\boldsymbol{X}}\boldsymbol{\times}\boldsymbol{B}^{*}, \end{gather}

(4.21)\begin{gather}\hat{\textsf{b}}\boldsymbol{\cdot}\dot{\boldsymbol{X}} = P_{\|}/m, \end{gather}

which are associated with virtual displacements $\delta \boldsymbol {X}$ and $\delta P_{\|}$, respectively. From these equations, we easily recover the guiding-centre Hamilton equations (4.12) and (4.13). Likewise, the guiding-centre Euler–Lagrange equation associated with the virtual displacement $\delta J$ yields (4.15), while the virtual displacement $\delta \theta$ yields $\dot {J} = 0$ as a result of Noether's theorem.

In addition to variations with respect to guiding-centre phase-space coordinates, the guiding-centre Lagrangian (4.19) can also be varied with respect to the electric and magnetic fields $(\delta \boldsymbol {E}, \delta \boldsymbol {B})$, which respectively yield expressions for the guiding-centre polarization and magnetization (Brizard 2008). Here, the Eulerian variation of the guiding-centre Lagrangian (4.19) is expressed as (Brizard 2023b)

(4.22)\begin{align} \delta L_{\rm gc} & \equiv \left(\frac{e}{c}\delta\boldsymbol{A} + \delta\boldsymbol{\varPi}_{\rm gc} \right)\boldsymbol{\cdot}\dot{\boldsymbol{X}} - (e\delta\varPhi +\delta K_{\rm gc}) - J (\delta\mathcal{S}+ \dot{\boldsymbol{X}}\boldsymbol{\cdot}\delta\boldsymbol{R}) \nonumber\\ & = \left(\frac{e}{c}\delta\boldsymbol{A}\boldsymbol{\cdot}\dot{\boldsymbol{X}} - e\delta\varPhi\right) + \boldsymbol{\rm \pi}_{\rm gc}\boldsymbol{\cdot}\delta\boldsymbol{E} + \left(\boldsymbol{\mu}_{\rm gc} + \boldsymbol{\rm \pi}_{\rm gc}\boldsymbol{\times} \frac{\boldsymbol{P}_{0}}{mc}\right)\boldsymbol{\cdot}\delta\boldsymbol{B} + ({\rm FLR}), \end{align}

where the FLR corrections, which are ignored in (4.22), are calculated to first order in a recent paper (Brizard 2023b). Here, the guiding-centre electric-dipole moment $\boldsymbol {{\rm \pi} }_\textrm {gc}$ is defined as (Pfirsch 1984; Pfirsch & Morrison 1985)

(4.23)\begin{equation} \boldsymbol{\rm \pi}_{\rm gc} \equiv \frac{e\hat{\textsf{b}}}{\varOmega}\boldsymbol{\times} (\dot{\boldsymbol{X}}-\boldsymbol{u}_{E}), \end{equation}

while the guiding-centre magnetic dipole moment $\boldsymbol {\mu }_\textrm {gc} + \boldsymbol {{\rm \pi} }_\textrm {gc}\boldsymbol {\times }\boldsymbol {P}_{0}/(mc)$ is defined as the sum of the intrinsic magnetic dipole moment

(4.24)\begin{equation} \boldsymbol{\mu}_{\rm gc} \equiv \mu \left(-\hat{\textsf{b}} + \frac{1}{\varOmega}\frac{{\rm d}\hat{\textsf{b}}}{{\rm d} t}\boldsymbol{\times}\hat{\textsf{b}}\right), \end{equation}

which includes the gyrogauge correction associated with $\delta \mathcal {S} + \dot {\boldsymbol {X}}\boldsymbol {\cdot }\delta \boldsymbol {R}$ (Brizard 2023b), and the moving electric-dipole contribution (Jackson 1975), expressed in terms of the lowest-order guiding-centre momentum $\boldsymbol {P}_{0} = P_{\|}\hat {\textsf {b}} + \boldsymbol {P}_{E}$. By ignoring these FLR corrections, the guiding-centre polarization and magnetization are defined as moments of the guiding-centre electric and magnetic dipole moments

(4.25)\begin{gather} \boldsymbol{\mathcal{P}}_{\rm gc} = \int_{P} F_{\rm gc}\boldsymbol{\rm \pi}_{\rm gc}, \end{gather}

(4.26)\begin{gather}\boldsymbol{\mathcal{M}}_{\rm gc} = \int_{P} F_{\rm gc}\left(\boldsymbol{\mu}_{\rm gc} + \boldsymbol{\rm \pi}_{\rm gc}\boldsymbol{\times}\frac{\boldsymbol{P}_{0}}{mc}\right), \end{gather}

where the guiding-centre phase-space density $F_\textrm {gc} \equiv \mathcal {J}_\textrm {gc}F$ includes the guiding-centre Jacobian $\mathcal {J}_\textrm {gc}$ and the notation $\int _{P}$ includes an integral over guiding-centre momentum space as well as a sum over particle species.

Finally, we note that, in the absence of an electric field, the classical expression (Kaufman 1986) for the guiding-centre electric-dipole moment $\boldsymbol {{\rm \pi} }_\textrm {gc} = (e\hat {\textsf {b}}/\varOmega )\boldsymbol {\times }\dot {\boldsymbol {X}}$ is derived by Lie-transform perturbation method (Tronko & Brizard 2015) only if the first-order polarization correction (3.37) is used.

5. Summary

In the present work, a set of higher-order guiding-centre Hamilton equations was derived by Lie-transform perturbation method for the case of time-dependent electric and magnetic fields that satisfy the standard guiding-centre space–time orderings. The second-order corrections in the guiding-centre Hamiltonian represented FLR corrections of the lowest-order electrostatic potential energy $e\varPhi$. Additional second-order corrections in the guiding-centre Lagrangian (4.19) included gyrogauge-invariance contributions to the guiding-centre Hamiltonian and Poisson bracket as well as corrections leading to the correct guiding-centre polarization.

When we turn our attention to the self-consistent interactions of the charged-particle guiding centres and the electromagnetic fields associated with plasma confinement, we need to derive a set of higher-order guiding-centre Vlasov–Maxwell equations. Work presented elsewhere (Brizard 2023b) considered the variational formulation of the higher-order guiding-centre Vlasov–Maxwell equations derived from the guiding-centre Lagrangian (4.19) and the Maxwell Lagrangian density.

According to this variational principle, using the guiding-centre Liouville equation (4.16), the guiding-centre Vlasov equation for the guiding-centre phase-space density $F_\textrm {gc} \equiv \mathcal {J}_\textrm {gc}F$ is written in divergence form as

(5.1)\begin{equation} \frac{\partial F_{\rm gc}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\dot{\boldsymbol{X}}F_{\rm gc}) + \frac{\partial }{\partial P_{\|}} (\dot{P}_{\|}F_{\rm gc}) = 0, \end{equation}

while the Maxwell equations with particle sources are expressed as

(5.2)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{E} = 4{\rm \pi}\varrho \equiv 4{\rm \pi} (\varrho_{\rm gc}-\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\mathcal{P}}_{\rm gc}), \end{gather}

(5.3)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{B} - \frac{1}{c}\frac{\partial \boldsymbol{E}}{\partial t} = \frac{4{\rm \pi}}{c}\boldsymbol{J} \equiv \frac{4{\rm \pi}}{c}\left(\boldsymbol{J}_{\rm gc} + \frac{\partial \boldsymbol{\mathcal{P}}_{\rm gc}}{\partial t} + c\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{\mathcal{M}}_{\rm gc}\right), \end{gather}

where the guiding-centre charge and current densities are $\varrho _\textrm {gc} = \int _{P} eF_\textrm {gc}$ and $\boldsymbol {J}_\textrm {gc} = \int _{P} e\dot {\boldsymbol {X}}F_\textrm {gc}$. Here, the guiding-centre polarization charge density $\varrho _\textrm {pol} \equiv -\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {\mathcal {P}}_\textrm {gc}$ and current density $\boldsymbol {J}_\textrm {pol} \equiv \partial \boldsymbol {\mathcal {P}}_\textrm {gc}/\partial t$ are derived from the guiding-centre polarization (4.25), while the guiding-centre magnetization current density $\boldsymbol {J}_\textrm {mag} \equiv c\boldsymbol {\nabla }\boldsymbol {\times }\boldsymbol {\mathcal {M}}_\textrm {gc}$ is derived from the guiding-centre magnetization (4.26). The remaining source-free Maxwell equations are Faraday's law $\partial \boldsymbol {B}/\partial t = -c\boldsymbol {\nabla }\boldsymbol {\times }\boldsymbol {E}$ and Gauss's law $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {B} = 0$. We note that the guiding-centre Vlasov–Maxwell variational principle also guarantees the existence of exact energy-momentum conservation laws, derived by the Noether method (Brizard 2008). Our recent work (Brizard 2023b) has considered the set of higher-order guiding-centre Vlasov–Maxwell equations obtained by explicitly imposing the quasineutrality constraint $\varrho _\textrm {gc} = \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\mathcal {P}}_\textrm {gc}$.

Future work will explore the Hamiltonian structure of the guiding-centre Vlasov–Maxwell equations, when guiding-centre polarization and magnetization are included, which will generalize our previous work (Brizard 2021), and its development is motivated by the desire to construct structure-preserving numerical algorithms (Morrison 2017) using an important set of reduced Vlasov–Maxwell equations.