# Faithful guiding-center orbits in an axisymmetric magnetic field

Alain J. Brizard1,a and Brook C. Hodgeman1 1Department of Physics, Saint
Michael’s College, Colchester, VT 05439, USA
aAuthor to whom correspondence should be addressed<EMAIL_ADDRESS>
###### Abstract

The problem of the charged-particle motion in an axisymmetric magnetic
geometry is used to assess the validity of higher-order Hamiltonian guiding-
center theory, which includes higher-order corrections associated with
gyrogauge invariance as well as guiding-center polarization induced by
magnetic-field non-uniformity. Two axisymmetric magnetic geometries are
considered: a magnetic mirror geometry and a simple tokamak geometry. When a
magnetically-confined charged-particle orbit is regular (i.e., its guiding-
center magnetic moment is adiabatically invariant), the guiding-center
approximation, which conserves both energy and azimuthal canonical angular
momentum, is shown to be faithful to the particle orbit when higher-order
corrections are taken into account.

##  I Introduction

The guiding-center representation Littlejohn (1983); Brizard (1989); Cary and
Brizard (2009); Tronko and Brizard (2015) of charged-particle orbits is at the
foundation of most particle simulations of magnetized plasmas Qin and Guan
(2008); White (2014); Burby and Ellison (2017); Albert et al. (2020); Bierwage
et al. (2022). The faithfulness of this representation relies on the adiabatic
invariance of the magnetic moment $\mu$, which is expressed as an asymptotic
expansion based on the nonuniformity of the magnetic field. In addition, when
the magnetic field is axisymmetric, the connection between the exact particle
canonical azimuthal angular momentum and its guiding-center representation can
be used as a test for the faithfulness of the guiding-center approximation
Belova et al. (2003).

For each guiding-center orbit parametrized by the guiding-center magnetic
moment $\mu$ (with initial guiding-center position ${\bf X}_{0}$ and initial
parallel guiding-center momentum $P_{\|0}$), there corresponds an infinite set
of particle orbits (with local initial conditions ${\bf x}_{0}$ and $\dot{\bf
x}_{0}$) that are labeled by an initial gyroangle, measured on an initial
Kruskal ring Burby and Qin (2012); Burby (2020) (also parametrized by the
magnetic moment $\mu$) defined on the two-dimensional plane perpendicular to
the local magnetic field. (In a uniform magnetic field ${\bf
B}=B\,\widehat{\sf z}$, the Kruskal ring is a circle in the $(x,y)$-plane of
radius $\sqrt{2\mu\,B/m\Omega^{2}}$ centered at the guiding-center position
${\bf X}_{0}$, where $\Omega=eB/mc$ denotes the gyrofrequency of a charged
particle of mass $m$ and charge $e$ and $c$ denotes the speed of light.) By
adopting a guiding-center formulation that is gyrogauge invariant Littlejohn
(1983); Brizard (1989) (i.e., a formulation that is not only gyroangle
invariant but also independent how the gyroangle is measured), the guiding-
center orbit can be compared to an arbitrary particle orbit belonging to the
same initial Kruskal ring.

The purpose of our present work is to explore how higher-order Hamiltonian
guiding-center theory is faithful to charged-particle dynamics in an
axisymmetric magnetic field. In particular, compared to the exact particle
canonical angular momentum invariant, we will investigate the faithfulness of
the guiding-center canonical angular momentum derived either in the truncated
(lowest-order) guiding-center model White and Chance (1984); Cary and Brizard
(2009), in the standard work of Littlejohn Littlejohn (1983); Brizard (1989),
which includes gyrogauge corrections, or the extended work of Brizard Brizard
(2013) and Tronko and Brizard Tronko and Brizard (2015), which includes
effects due to guiding-center polarization Kaufman (1986). For this purpose,
we will consider particle and guiding-center orbits in axisymmetric mirror
geometry (Sec. III) and in axisymmetric (simple) tokamak geometry (Sec. IV).

##  II Particle and Guiding-center Orbits in General Axisymmetric Magnetic
Geometry

In the present Section, we consider the problem of charged particle motion in
a generic nonuniform magnetic field. Since the primary focus of our work
involves the effects of magnetic nonuniformity, we assume that the magnetic
field is stationary and an electric field is absent in our formulation.

###  II.1 Lagrangian particle dynamics

The orbits of a charged particle moving in a nonuniform magnetic field ${\bf
B}({\bf x})$ are represented as solutions of the Euler-Lagrange equations
obtained from the particle Lagrangian

$\displaystyle L$ $\displaystyle=$ $\displaystyle\left(\frac{e}{c}\,{\bf
A}({\bf x})\;+\;m\,{\bf v}\right)\,\mbox{\boldmath$\cdot$}\,\dot{\bf
x}\;-\;\frac{m}{2}\,|{\bf v}|^{2}$ (1) $\displaystyle=$
$\displaystyle\frac{e}{c}\,{\bf A}({\bf x})\,\mbox{\boldmath$\cdot$}\,\dot{\bf
x}\;+\;\frac{m}{2}\;\left|\dot{\bf x}\right|^{2},$

where the particle velocity is ${\bf v}=\dot{\bf x}$ and the magnetic field
${\bf B}\equiv\nabla\,\mbox{\boldmath$\times$}\,{\bf A}\equiv B\,\widehat{{\sf
b}}$ (which is decomposed in terms of its magnitude $B=|{\bf B}|$ and its
direction unit vector $\widehat{{\sf b}}={\bf B}/|{\bf B}|$) is represented in
terms of a vector potential ${\bf A}$. From this Lagrangian, we obtain the
Euler-Lagrange equations

$\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{\bf
x}}\right)\;=\;\frac{\partial L}{\partial\bf x},$ (2)

which yield the usual Lorentz force equation

$m\;\ddot{\bf x}\;=\;\frac{e}{c}\left(\nabla{\bf
A}\,\mbox{\boldmath$\cdot$}\,\dot{\bf x}\;-\frac{}{}\dot{\bf
x}\,\mbox{\boldmath$\cdot$}\,\nabla{\bf A}\right)\;=\;\frac{e}{c}\,\dot{\bf
x}\,\mbox{\boldmath$\times$}\,{\bf B},$ (3)

which is solved subject to the initial conditions $({\bf x}_{0},\dot{\bf
x}_{0}={\bf v}_{0})$. Since we are interested in orbital solutions of the
Lorentz force equation (3) over long time scales compared to the short
gyration period, which is inversely proportional to the gyrofrequency
$\Omega_{0}=eB_{0}/(mc)$, where $B_{0}$ denotes the characteristic strength of
the magnetic field, we introduce a dimensionless time
$t^{\prime}=\epsilon\,\Omega_{0}\,t$, where $\epsilon\ll 1$ denotes a small
ordering parameter, so that $\dot{\bf x}=\epsilon\,\Omega_{0}\;{\bf
x}^{\prime}$ (a prime denotes a derivative with respect to $t^{\prime}$).
Hence, the Lorentz force equation (3) becomes

$\epsilon\;\overline{\bf x}^{\prime\prime}\;=\;\overline{\bf
x}^{\prime}\,\mbox{\boldmath$\times$}\,\overline{\bf B}(\overline{\bf x}),$
(4)

where we have introduced a characteristic length scale $R_{0}$ associated with
the magnetic field, so that $\overline{\bf x}\equiv{\bf x}/R_{0}$ is
dimensionless, and the magnetic field ${\bf B}=B_{0}\,\overline{\bf
B}(\overline{\bf x})$ is expressed in terms of a dimensionless field
$\overline{\bf B}(\overline{\bf x})$. The solution for this equation of
motion, which can be carried out as an asymptotic expansion in powers of
$\epsilon$ Kruskal (1958), will be carried out numerically in this paper.

In the event the magnetic field is axisymmetric, the particle Lagrangian (1)
is independent of the particle azimuthal angle $\varphi$, and the azimuthal
canonical angular momentum

$P_{\varphi}\;\equiv\;\frac{\partial
L}{\partial\dot{\varphi}}\;=\;\left(\frac{e}{c}\,{\bf A}\;+\;m\,\dot{\bf
x}\right)\,\mbox{\boldmath$\cdot$}\,\frac{\partial\bf x}{\partial\varphi}$ (5)

is a constant of the motion for particle orbits. With the normalization
discussed above, we note that the azimuthal canonical angular momentum (5)
becomes

$\frac{P_{\varphi}}{mR_{0}^{2}\Omega_{0}}\;\equiv\;\overline{P}_{\varphi}\;=\;\left(\overline{\bf
A}\;+\frac{}{}\epsilon\,\overline{\bf
x}^{\prime}\right)\,\mbox{\boldmath$\cdot$}\,\frac{\partial\overline{\bf
x}}{\partial\varphi},$ (6)

with the dimensionless magnetic vector potential $\overline{\bf A}\equiv{\bf
A}/(B_{0}R_{0})$.

###  II.2 Lagrangian guiding-center dynamics

For most particle orbits that are solutions of the Lorentz force equation (3),
the lowest-order magnetic moment

$\mu_{0}\equiv\frac{m}{2B}\;|\widehat{{\sf
b}}\,\mbox{\boldmath$\times$}\,\dot{\bf
x}|^{2}\;=\;\frac{m}{2B}\;|\mbox{\boldmath$\rho$}_{0}|^{2}\,\Omega^{2}$ (7)

is an adiabatic invariant (where
$\mbox{\boldmath$\rho$}_{0}\equiv\widehat{{\sf
b}}\,\mbox{\boldmath$\times$}\,\dot{\bf x}/\Omega$ denotes the lowest-order
gyroradius), i.e., while the time derivative

$\dot{\mu}_{0}\;=\;-\,\mu_{0}\,{\bf v}\,\mbox{\boldmath$\cdot$}\,\nabla\ln
B-\frac{mv_{\|}}{B}\;{\bf v}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,{\bf v}_{\bot}\;\neq\;0$ (8)

does not vanish for a general magnetic field, its average over the fast
gyromotion time scale yields

$\langle\dot{\mu}_{0}\rangle=-\,\mu_{0}\,v_{\|}\left(\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\ln
B\;+\frac{}{}\nabla\,\mbox{\boldmath$\cdot$}\,\widehat{{\sf
b}}\right)=-\,\frac{\mu_{0}\,v_{\|}}{B}\;(\nabla\,\mbox{\boldmath$\cdot$}\,{\bf
B}),$

where $v_{\|}\equiv{\bf v}\,\mbox{\boldmath$\cdot$}\,\widehat{{\sf b}}$ is the
local parallel velocity, $\langle{\bf v}_{\bot}\rangle=0$ and $\langle m{\bf
v}_{\bot}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,{\bf
v}_{\bot}\rangle=\mu_{0}B\,(\nabla\,\mbox{\boldmath$\cdot$}\,\widehat{{\sf
b}})$. Since magnetic fields are divergenceless, we immediately find that
$\langle\dot{\mu}_{0}\rangle=0$, i.e., $\mu_{0}$ is an invariant over time
scales that are slow compared to the fast gyromotion time scale.

The purpose of the guiding-center transformation is to construct an expression
for the guiding-center moment

$\mu\;=\;\mu_{0}\;+\;\epsilon\,\mu_{1}\;+\;\cdots$ (9)

represented as an asymptotic series in powers of the dimensionless parameter
$\epsilon$, where the first-order correction

$\displaystyle\mu_{1}$ $\displaystyle=$ $\displaystyle\left(\mu_{0}\nabla\ln
B+\frac{p_{\|}^{2}\,\mbox{\boldmath$\kappa$}}{2\,mB}\right)\,\mbox{\boldmath$\cdot$}\,\mbox{\boldmath$\rho$}_{0}-\frac{3}{2}\,\mu_{0}\left(\frac{p_{\|}\,\tau}{m\Omega}\right)$
(10) $\displaystyle+\;\frac{p_{\|}}{2B}\;\frac{d\widehat{{\sf
b}}}{dt}\,\mbox{\boldmath$\cdot$}\,\mbox{\boldmath$\rho$}_{0}$

involves first-order magnetic-field nonuniformity associated with magnetic
curvature $\mbox{\boldmath$\kappa$}=\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf b}}$ and magnetic twist
$\tau=\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf
b}}$, with $d\widehat{{\sf b}}/dt\equiv\dot{\bf
x}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf b}}$ in a time-independent
nonuniform magnetic field. The new expression (10), which is derived from the
standard expression found in Refs. Littlejohn (1983); Brizard (1989); Cary and
Brizard (2009); Tronko and Brizard (2015) in App. A, is easily computed from
the particle dynamics. Hence, from the magnetic-moment analysis of the
particle orbit yields a relatively accurate value for the guiding-center
magnetic moment $\mu=\mu_{0}+\epsilon\,\mu_{1}$, which can then be used as a
label for the guiding-center orbit.

The reduced guiding-center representation of charged-particle dynamics in a
nonuniform magnetic field Cary and Brizard (2009) is obtained by an asymptotic
decoupling of the fast gyromotion from the slow magnetic-drift motion in a
reduced dynamical phase space with guiding-center coordinates
$Z^{\alpha}=({\bf X},P_{\|})$, while the fast gyromotion is represented by the
canonically-conjugate guiding-center action-angle coordinates $(J,\zeta)$,
where the gyroaction $J\equiv\mu B_{0}/\Omega_{0}$ (defined in terms of the
magnetic moment $\mu$) is canonically conjugate to the guiding-center
gyroangle $\zeta$.

The guiding-center Lagrangian is expressed up to first order in magnetic-field
nonuniformity as

$\displaystyle L_{\rm gc}$ $\displaystyle=$
$\displaystyle\left(\frac{e}{\epsilon\,c}\,{\bf A}({\bf
X})\;+\;P_{\|}\,\widehat{{\sf b}}({\bf
X})\;-\;\epsilon\,J\;\mbox{\boldmath$\cal R$}^{*}({\bf
X})\right)\,\mbox{\boldmath$\cdot$}\,\dot{\bf X}$ (11)
$\displaystyle+\;\epsilon\,J\;\dot{\zeta}\;-\;\left(\frac{P_{\|}^{2}}{2m}\;+\;\mu
B({\bf X})\right),$

where the $\epsilon$-ordering introduced in Eq. (11) is based on the standard
macroscopic ordering associated with the renormalization of the electric
charge $e\rightarrow e/\epsilon$. The selection of the vector field

$\mbox{\boldmath$\cal R$}^{*}({\bf
X})\;=\;\left\\{\begin{array}[]{lr}0&(\mbox{A})\\\ \mbox{\boldmath$\cal
R$}\;+\;\frac{1}{2}\,\tau\,\widehat{{\sf b}}&(\mbox{B})\\\
\mbox{\boldmath$\cal
R$}\;+\;\frac{1}{2}\,\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf
b}}&(\mbox{C})\end{array}\right.$ (12)

is based on whether the gyrogauge vector field $\mbox{\boldmath$\cal
R$}\equiv\nabla\widehat{\sf e}_{1}\,\mbox{\boldmath$\cdot$}\,\widehat{\sf
e}_{2}$ (defined Littlejohn (1983) in terms of the local orthogonal unit
vectors $\widehat{\sf e}_{1}$ and $\widehat{\sf e}_{2}\equiv\widehat{{\sf
b}}\,\mbox{\boldmath$\times$}\,\widehat{\sf e}_{1}$) is kept (B & C) or not
(A), and whether the guiding-center polarization correction Tronko and Brizard
(2015) is kept (C) or not (B). The truncated guiding-center model (A) is the
simplest guiding-center model that is used in several guiding-center orbit
codes (e.g., White and Chance (1984); White (2014)) and is reviewed in Ref.
Cary and Brizard (2009). The standard guiding-center model (B) was derived
Littlejohn (1983); Brizard (1989) to ensure that the guiding-center equations
of motion are not only independent of the gyroangle but also independent of
how this gyroangle is locally measured in the perpendicular plane spanned by
the unit vectors $(\widehat{\sf e}_{1},\widehat{\sf e}_{2})$. The extended
guiding-center model (C) was derived Tronko and Brizard (2015) to ensure that
the guiding-center transformation accurately represents the guiding-center
polarization Kaufman (1986); Brizard (2013).

The guiding-center Euler-Lagrange equations

$d(\partial L_{\rm gc}/\partial\dot{Z}^{\alpha})/dt\;=\;\partial L_{\rm
gc}/\partial Z^{\alpha}$

lead to the reduced guiding-center equations of motion

$\displaystyle\dot{\bf X}$ $\displaystyle=$
$\displaystyle\frac{P_{\|}}{m}\;\frac{{\bf
B}^{*}}{B_{\|}^{*}}\;+\;\frac{\epsilon\,c\widehat{{\sf
b}}}{eB_{\|}^{*}}\,\mbox{\boldmath$\times$}\,\mu\,\nabla B,$ (13)
$\displaystyle\dot{P}_{\|}$ $\displaystyle=$ $\displaystyle-\;\frac{{\bf
B}^{*}}{B_{\|}^{*}}\,\mbox{\boldmath$\cdot$}\,\mu\,\nabla B,$ (14)

where

${\bf B}^{*}\;=\;{\bf
B}\;+\;\frac{\epsilon\,c}{e}\left(P_{\|}\;\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf
b}}\;-\frac{}{}\epsilon\,J\;\nabla\,\mbox{\boldmath$\times$}\,\mbox{\boldmath$\cal
R$}^{*}\right),$ (15)

and $B_{\|}^{*}\equiv\widehat{{\sf b}}\,\mbox{\boldmath$\cdot$}\,{\bf B}^{*}$
can be used as the guiding-center Jacobian. These guiding-center equations are
solved subject to the initial conditions $({\bf X}_{0},P_{\|0})$, once again
labeled by the guiding-center magnetic moment $\mu$ obtained from the particle
orbit. Because the particle and guiding-center orbits share the same values of
energy $E$ and magnetic moment $\mu$, the initial guiding-center parallel
momentum can be chosen from the initial condition
$P^{2}_{\|0}/2m=E-\mu\,B({\bf X}_{0})$, where the initial guiding-center
position ${\bf X}_{0}$ is calculated from the initial particle position ${\bf
x}_{0}$ according to the guiding-center transformation Littlejohn (1983);
Brizard (1989)

${\bf X}\;=\;{\bf x}+\epsilon\,G_{1}^{{\bf x}}+\epsilon^{2}\,G_{2}^{{\bf
x}}+\frac{1}{2}\,\epsilon^{2}\,{\sf G}_{1}\cdot{\sf d}G_{1}^{{\bf x}}+\cdots,$
(16)

which implies that the initial guiding-center position ${\bf X}_{0}$ is
shifted from the initial particle position ${\bf x}_{0}$. Hence, the initial
guiding-center position ${\bf X}_{0}$ will depend on the guiding-center model
used in Eq. (12), which differs at second order through $G_{2}^{{\bf x}}$
Tronko and Brizard (2015).

Lastly, when the magnetic field is axisymmetric, the guiding-center Lagrangian
(11) is independent of the guiding-center azimuthal angle $\Phi$, and the
guiding-center azimuthal canonical angular momentum

$P_{{\rm gc}\Phi}=\left(\frac{e{\bf A}}{\epsilon\,c}+P_{\|}\,\widehat{{\sf
b}}-\epsilon\,J\;\mbox{\boldmath$\cal
R$}^{*}\right)\,\mbox{\boldmath$\cdot$}\,\frac{\partial\bf
X}{\partial\Phi}\equiv\frac{e{\bf A}^{*}}{\epsilon
c}\,\mbox{\boldmath$\cdot$}\,\frac{\partial\bf X}{\partial\Phi}$ (17)

is an exact guiding-center invariant. We note that the terms of third order in
$\epsilon^{3}$ in ${\bf A}^{*}$ (i.e., second order in magnetic-field
nonuniformity) and higher are ignored.

###  II.3 Validity of the guiding-center representation in general
axisymmetric magnetic geometry

While the azimuthal canonical angular momenta (5) and (17) are constants of
motion of their respective equations of motion, they can only be compared when
the guiding-center azimuthal canonical angular momentum (17) is pulled back
into particle phase space:

$\displaystyle{\sf T}_{\rm gc}P_{{\rm gc}\Phi}$ $\displaystyle=$
$\displaystyle P_{{\rm
gc}\Phi}+\epsilon\left(G_{1}^{\mu}\frac{\partial}{\partial\mu}+G_{1}^{p_{\|}}\frac{\partial}{\partial
p_{\|}}-\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\right)P_{{\rm
gc}\Phi}$ (18)
$\displaystyle+\;\epsilon^{2}\left(\frac{1}{2}\,\mbox{\boldmath$\rho$}_{0}\mbox{\boldmath$\rho$}_{0}:\nabla\nabla-\mbox{\boldmath$\rho$}_{1}\mbox{\boldmath$\cdot$}\nabla\right)P_{{\rm
gc}\Phi}+\cdots$ $\displaystyle\equiv$ $\displaystyle P_{\varphi},$

where $(G_{1}^{\mu},G_{1}^{p_{\|}})$ denote the first-order corrections to the
guiding-center magnetic moment and guiding-center parallel momentum,
respectively, and the first-order gyroradius is defined in particle phase
space as

$\mbox{\boldmath$\rho$}_{1}\;=\;-\;\frac{\mu_{0}B\,\mbox{\boldmath$\kappa$}}{2m\Omega^{2}}+\left(\frac{\mbox{\boldmath$\rho$}_{0}}{2}\,\mbox{\boldmath$\cdot$}\,\nabla\ln
B-\frac{p_{\|}\tau}{m\Omega}\right)\mbox{\boldmath$\rho$}_{0}+\rho_{1\|}\,\widehat{{\sf
b}},$ (19)

where, while an explicit expression for the parallel component
$\rho_{1\|}\equiv\mbox{\boldmath$\rho$}_{1}\,\mbox{\boldmath$\cdot$}\,\widehat{{\sf
b}}$ will not be needed in the present work, we note that its gyroangle
average is $\langle\rho_{1\|}\rangle=(\mu
B/2m\Omega^{2})\nabla\,\mbox{\boldmath$\cdot$}\,\widehat{{\sf b}}$. We note
that the first term in Eq. (19) appears as a result of the extension Tronko
and Brizard (2015) of the standard guiding-center transformation Littlejohn
(1983) that correctly calculates the guiding-center polarization (see App. B
for details).

The guiding-center representation is faithful to particle dynamics if the
identity (18) is satisfied up to an arbitrary order in $\epsilon$. This
identity guarantees that the guiding-center push-forward of the particle
conservation law $dP_{\varphi}/dt=0$ yields the guiding-center conservation
law

$0={\sf T}_{\rm gc}^{-1}\left(\frac{dP_{\varphi}}{dt}\right)=\left[{\sf
T}_{\rm gc}^{-1}\left(\frac{d}{dt}{\sf T}_{\rm gc}\right)\right]{\sf T}_{\rm
gc}^{-1}P_{\varphi}\equiv\frac{d_{\rm gc}P_{{\rm gc}\Phi}}{dt},$

where $d_{\rm gc}/dt$ is the time derivative generated by the guiding-center
Lagrangian dynamics and the guiding-center invariant $P_{{\rm
gc}\Phi}\equiv{\sf T}_{\rm gc}^{-1}P_{\varphi}$ may be truncated at an
arbitrary order in $\epsilon$.

In previous work, Belova et al. Belova et al. (2003) considered energetic-
particle orbits in the National Spherical Torus Experiment (NSTX) that
satisfied the adiabatic invariance of the higher-order guiding-center magnetic
moment (9), and numerically investigated the validity of the guiding-center
representation by verifying that the explicit expression for the guiding-
center pull-back of the guiding-center toroidal canonical angular momentum

${\sf T}_{\rm gc}P_{{\rm gc}\Phi}=P_{{\rm gc}\Phi}\left({\bf
x}-\epsilon\mbox{\boldmath$\rho$},p_{0\|}+\epsilon\,G_{1}^{p_{\|}},\mu_{0}+\epsilon
G_{1}^{\mu}\right)$ (20)

is nearly invariant, where $\mbox{\boldmath$\rho$}={\bf x}-{\sf T}_{\rm
gc}{\bf X}$ includes the first-order corrections (19) due to magnetic
nonuniformity. In Fig. 3 of Ref. Belova et al. (2003), Belova et al. Belova et
al. (2003) show that the standard guiding-center expression
$\psi^{*}=\psi-\epsilon\,(cP_{\|}/e)b_{\Phi}+\epsilon^{2}(cJ/e){\cal
R}_{\Phi}^{*}$, where ${\cal R}_{\Phi}^{*}=b_{z}+\frac{1}{2}\tau b_{\Phi}$,
yields an improved particle canonical angular momentum invariant (20) compared
with the truncated guiding-center expression
$\psi^{*}=\psi-\epsilon\,(cP_{\|}/e)b_{\Phi}$.

In the present paper, the validity of the guiding-center representation is
assessed on the basis of verifying that the expansion (18) is exactly valid at
each order in $\epsilon$ for two axisymmetric magnetic geometries: mirror
geometry (Sec. III) and simple tokamak geometry (Sec. IV). Here, the guiding-
center representation is shown to be faithful up to (and including) first
order in magnetic-field nonuniformity, which requires retaining all terms
appearing in Eq. (17).

##  III Magnetic Mirror Geometry

We begin by considering the problem of charged-particle motion in axisymmetric
magnetic mirror geometry, where the magnetic field is represented, using the
cylindrical coordinates $(r,\varphi,z)$, by the dimensionless expression

$\displaystyle{\bf B}$ $\displaystyle=$
$\displaystyle\nabla\psi(r,z)\,\mbox{\boldmath$\times$}\,\nabla\varphi\;=\;\frac{1}{r}\left(\frac{\partial\psi}{\partial
r}\;\widehat{\sf z}\;-\;\frac{\partial\psi}{\partial z}\;\widehat{\sf
r}\right)$ (21) $\displaystyle=$ $\displaystyle B_{r}(r,z)\;\widehat{\sf
r}\;+\;B_{z}(z)\;\widehat{\sf z},$

where the magnetic flux $\psi(r,z)$ is defined as

$\psi(r,z)\;=\;\frac{1}{2}\;r^{2}\,(1+z^{2}).$ (22)

Figure 1 shows the magnetic-mirror vector field (21) in the $(x,z)$-plane.
With a field line passing through the point $r=r_{0}$ on the equatorial plane
$z=0$, which is labeled by the magnetic flux $\psi(r_{0},0)=r_{0}^{2}/2$, the
radial coordinate $r(z)=r_{0}/\sqrt{1+z^{2}}$ of the field line can be
expressed as a function of $z$.

Figure 1: Plot of the magnetic-mirror vector field (21) in the $(x,z)$-plane.

The magnitude of the magnetic field (21) is

$B(r,z)\;=\;\sqrt{r^{2}\,z^{2}\;+\;(1+z^{2})^{2}},$ (23)

while the unit vector along the magnetic field is

$\widehat{{\sf b}}(r,\varphi,z)\;=\;b_{r}(r,z)\;\widehat{\sf
r}(\varphi)\;+\;b_{z}(r,z)\;\widehat{\sf z}$ (24)

where $b_{r}(r,z)=-\,r\,z/B(r,z)$, $b_{z}(r,z)=(1+z^{2})/B(r,z)$, and
$\widehat{\sf r}(\varphi)=\cos\varphi\;\widehat{\sf
x}+\sin\varphi\;\widehat{\sf y}$. Next, we can calculate

$\displaystyle\widehat{{\sf b}}\,\mbox{\boldmath$\times$}\,\nabla\ln B$
$\displaystyle=$
$\displaystyle-\frac{r\,b_{r}}{B}\left(3\,zb_{z}-rb_{r}\right)\nabla\varphi\equiv-\,K_{1}\nabla\varphi,$
(25) $\displaystyle\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf b}}$
$\displaystyle=$
$\displaystyle-\,\left(r^{2}/B+K_{1}\right)\nabla\varphi\equiv-\,K_{2}\nabla\varphi,$
(26)

with $\nabla\,\mbox{\boldmath$\times$}\,{\bf B}=-r^{2}\,\nabla\varphi$.

###  III.1 Particle dynamics

The normalized equations of motion (4) are expressed in cylindrical
coordinates as

$\epsilon\,{\bf
x}^{\prime\prime}=r\,\varphi^{\prime}\left(B_{z}\,\widehat{r}-B_{r}\,\widehat{\sf
z}\right)-\left(B_{z}\,r^{\prime}-B_{r}\,z^{\prime}\right)\widehat{\varphi},$
(27)

where ${\bf
x}^{\prime}=r^{\prime}\,\widehat{r}+r\,\varphi^{\prime}\,\widehat{\varphi}+z^{\prime}\,\widehat{\sf
z}$, with $\widehat{\varphi}=\partial\widehat{r}/\partial\varphi$. These
dimensionless equations are numerically solved for $\epsilon=1/20$, with the
initial conditions $(r_{0},\varphi_{0},z_{0})=(1,0,0)$ and
$(r_{0}^{\prime},\varphi_{0}^{\prime},z_{0}^{\prime})=(0,1/10,\sqrt{24}/10)$
associated with the dimensionless energy $E=1/4$.

We note that, because of the azimuthal symmetry of the magnetic field (21),
i.e., the components $(B_{r},B_{z})$ are independent of the azimuthal angle
$\varphi$, the azimuthal canonical angular momentum

$P_{\varphi}\;=\;\frac{1}{\epsilon}\,\psi(r,z)\;+\;r^{2}\varphi^{\prime}$ (28)

is a constant of the motion.

Figure 2 shows that, while the lowest-order normalized magnetic moment
(normalized to $\epsilon^{2}m\Omega_{0}^{2}/2B_{0}$)

$\mu_{0}\;=\;\left[\left(r\,\varphi^{\prime}\right)^{2}\;+\;\left(b_{z}\,r^{\prime}\;-\;b_{r}\,z^{\prime}\right)^{2}\right]/B$
(29)

is relatively well conserved when the particle is near the orbital bounce
points, its adiabatic invariance is compromised as the particle crosses the
equatorial plane $(z=0)$ between $t^{\prime}=28$ and 29. The addition of the
first-order correction

$\mu_{1}\;=\;\left(\mu_{0}\;K_{1}\;+\;\overline{p}_{\|}^{2}\;K_{2}/B\;+\;\overline{p}_{\|}z^{\prime}\;K_{3}\right)\varphi^{\prime}/B$
(30)

computed from Eq. (10) (with $\tau=\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf
b}}=0$), where $\overline{p}_{\|}=b_{r}\,r^{\prime}+b_{z}\,z^{\prime}$ and

$\begin{array}[]{rcl}\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\ln
B&=&\epsilon\,\varphi^{\prime}K_{1}/B,\\\
\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,(\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}})&=&\epsilon\,\varphi^{\prime}K_{2}/B,\\\
\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\widehat{{\sf
b}}^{\prime}&=&\epsilon\,z^{\prime}\varphi^{\prime}K_{3},\end{array}$ (31)

greatly improves the adiabatic invariance of the magnetic moment and validates
the guiding-center representation for the particle orbits in the magnetic-
mirror vector field (21). Here, $K_{1}$ and $K_{2}$ are defined in Eqs.
(25)-(26) and $K_{3}=r(3zb_{r}+rb_{z})/B^{2}$. Hence, the adiabatic invariance
of the magnetic moment justifies our use of the guiding-center approximation
in describing particle orbits in magnetic mirror geometry. The numerical value
$\mu\simeq 0.01257$ will be used as the normalized guiding-center magnetic
moment in the guiding-center equations of motion in magnetic mirror geometry.

Figure 2: Plots of the normalized lowest-order magnetic moment $\mu_{0}$
(gray) and the normalized magnetic moment $\mu=\mu_{0}+\epsilon\,\mu_{1}\simeq
0.01257$ (black), with first-order correction $\mu_{1}$ given by Eq. (30).
Here, $\epsilon=1/20$ and the normalized kinetic energy is $E_{0}=1/4$, with
the initial conditions $(r_{0},\varphi_{0},z_{0})=(1,0,0)$ and
$(r_{0}^{\prime},\varphi_{0}^{\prime},z_{0}^{\prime})=(0,1/10,\sqrt{24}/10)$.

###  III.2 Guiding-center dynamics

The guiding-center Lagrangian in magnetic mirror geometry is

$L_{\rm gc}\;=\;\frac{1}{\epsilon}\psi^{*}\,\Phi^{\prime}+P_{\|}\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,{\bf
X}^{\prime}\;-\;\left(\frac{1}{2}\,P_{\|}^{2}+J\,B\right),$ (32)

where we use the extended guiding-center model (C) in Eq. (12), ${\bf
X}=(R,\Phi,Z)$ denote the guiding-center position in cylindrical geometry,
$P_{\|}$ denotes the normalized guiding-center parallel momentum, and
$J\equiv\mu/2$ denotes the normalized guiding-center gyroaction. In addition,
in magnetic mirror geometry, we may choose the perpendicular unit vectors
$\widehat{\sf e}_{1}=\widehat{\Phi}$ and $\widehat{\sf e}_{2}=\widehat{{\sf
b}}\,\mbox{\boldmath$\times$}\,\widehat{\Phi}$, so that the gyrogauge vector
$\mbox{\boldmath$\cal R$}=\nabla\widehat{\sf
e}_{1}\,\mbox{\boldmath$\cdot$}\,\widehat{\sf e}_{2}=b_{z}\,\nabla\Phi$ and
the effective magnetic flux

$\psi^{*}\;\equiv{}\;\psi(R,Z)\;-\;\epsilon^{2}J\,\left(b_{z}-\frac{1}{2}\,K_{2}\right)$
(33)

is expressed in terms of the extended guiding-center model (C) in Eq. (12):

$\mbox{\boldmath$\cal
R$}^{*}\;=\;\left(b_{z}-\frac{1}{2}K_{2}\right)\nabla\Phi.$ (34)

From this Lagrangian, we obtain Euler-Lagrange equations that can be expressed
as Eqs. (13)-(14), where

$\displaystyle{\bf B}^{*}$ $\displaystyle=$
$\displaystyle\nabla\psi^{*}\,\mbox{\boldmath$\times$}\,\nabla\Phi\;-\;\epsilon\,P_{\|}\;K_{2}\,\nabla\Phi,$
$\displaystyle B_{\|}^{*}$ $\displaystyle=$ $\displaystyle\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\psi^{*}\,\mbox{\boldmath$\times$}\,\nabla\Phi,$

The guiding-center equations are thus expressed as

$\displaystyle(\dot{R},\;\dot{Z})$ $\displaystyle=$
$\displaystyle\frac{P_{\|}}{RB_{\|}^{*}}\left(-\;\frac{\partial\psi^{*}}{\partial
Z},\frac{\partial\psi^{*}}{\partial R}\right),$ (36) $\displaystyle\dot{\Phi}$
$\displaystyle=$
$\displaystyle-\;\frac{\epsilon}{R^{2}B_{\|}^{*}}\left(P_{\|}^{2}\;K_{2}\;+\;J\,B\;K_{1}\right),$
(37) $\displaystyle\dot{P}_{\|}$ $\displaystyle=$
$\displaystyle-\;\frac{J}{RB_{\|}^{*}}\left(\frac{\partial B}{\partial
Z}\,\frac{\partial\psi^{*}}{\partial R}-\frac{\partial B}{\partial
R}\,\frac{\partial\psi^{*}}{\partial Z}\right),$ (38)

which exactly conserve the (dimensionless) guiding-center energy ${\cal
E}=P_{\|}^{2}/2+JB$ (where $J=\mu/2$ is obtained from Fig. 2) and the guiding-
center canonical azimuthal angular momentum $P_{{\rm
gc}\Phi}=\epsilon^{-1}\psi^{*}$. These dimensionless guiding-center equations
are solved for $\epsilon=1/20$ with the initial conditions
$(R_{0},\Phi_{0},Z_{0})=(201/200,0,0)$, which takes into account the radial
shift (16) from the initial particle position
$(r_{0},\varphi_{0},z_{0})=(1,0,0)$, and $P_{\|0}=\sqrt{E-J\,B(R_{0},Z_{0})}$.

Figure 3 shows shows plots of the particle azimuthal angle $\varphi$ (gray)
and the guiding-center azimuthal angle $\Phi$ (black) for a particle orbit in
the magnetic-mirror vector field (21) during a full bounce period. These
orbital solutions are obtained by numerical integration of the normalized
equations of motion (27) and (36)-(38), while conserving energy and azimuthal
canonical angular momentum within machine precision. We note that the guiding-
center azimuthal angle changes very slowly when the particle is near a turning
point, while it changes rapidly as the particle crosses the equatorial plane
($z=0$).

Figure 4, on the other hand, shows plots of the particle orbit (gray) and the
guiding-center orbit (black) during a bounce period in the $(y,z)$-plane.
While the motion is periodic in $(r,z)$, there is a slow drift motion in the
azimuthal direction, which can be seen in Figs. 3 and 4 (the slow drift motion
is proceeds to the left on the $y$-axis). The rapid, small-amplitude
oscillations that are noticeable in Figs. 3 and 4 are due to the fast
gyromotion of a charged particle about a magnetic-field line.

Figure 3: Plots of the particle azimuthal angle $\varphi$ (gray) and the
guiding-center azimuthal angle $\Phi$ (black) during a bounce period. Figure
4: Plots of particle orbit (gray) and guiding-center orbit (black) in the
$(y,z)$ plane during a bounce period. Here, the particle and guiding-center
orbits begin at point A and, after one bounce period, the guiding-center orbit
has drifted (in the negative-y direction) to point B.

###  III.3 Validity of the guiding-center representation in mirror geometry

We now show that the guiding-center representation of particle motion in
magnetic mirror geometry is faithful to the exact particle motion, by showing
that the guiding-center pull-back ${\sf T}_{\rm gc}P_{{\rm
gc}\Phi}=P_{\varphi}$ of the guiding-center canonical azimuthal angular
momentum is equal to the particle canonical azimuthal angular momentum.

Up to second order in $\epsilon$, the guiding-center pull-back ${\sf T}_{\rm
gc}P_{{\rm gc}\Phi}$ is expressed as

$\displaystyle{\sf T}_{\rm gc}P_{{\rm gc}\Phi}$ $\displaystyle=$
$\displaystyle\frac{1}{\epsilon}\left(\psi-\epsilon\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\psi-\epsilon^{2}\mbox{\boldmath$\rho$}_{1}\,\mbox{\boldmath$\cdot$}\,\nabla\psi\right)$
(39)
$\displaystyle+\;\frac{\epsilon}{2}\;\mbox{\boldmath$\rho$}_{0}\mbox{\boldmath$\rho$}_{0}:\nabla\nabla\psi-\epsilon\,J\left(b_{z}-\frac{1}{2}\,K_{2}\right),$

where the first-order gyroradius correction is given by Eq. (19) and the
contribution associated with the magnetic twist $\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf
b}}=0$ vanishes for magnetic mirror geometry. Here, we find
$\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\psi=-\,r^{2}\varphi^{\prime}$,
$\mbox{\boldmath$\kappa$}\,\mbox{\boldmath$\cdot$}\,\nabla\psi=-\,K_{2}B$, and
$\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\ln
B=K_{1}\,\varphi^{\prime}/B$, so that

$-\,\mbox{\boldmath$\rho$}_{1}\,\mbox{\boldmath$\cdot$}\,\nabla\psi\;=\;-\,\frac{1}{2}\,J\,K_{2}\;+\;\frac{r^{2}\varphi^{\prime
2}}{2\;B}\,K_{1},$

while

$\frac{1}{2}\,\mbox{\boldmath$\rho$}_{0}\mbox{\boldmath$\rho$}_{0}:\nabla\nabla\psi\;=\;J\,b_{z}\;-\;\frac{r^{2}\varphi^{\prime
2}}{2\;B}\,K_{1}.$

Hence, up to second order in $\epsilon$, the guiding-center pull-back (39)
yields

${\sf T}_{\rm gc}P_{{\rm
gc}\Phi}\;=\;\frac{1}{\epsilon}\,\psi\;+\;r^{2}\varphi^{\prime}\;\equiv\;P_{\varphi},$
(40)

which confirms the validity of the guiding-center representation in magnetic
mirror geometry.

##  IV Simple Tokamak Magnetic Geometry

We now turn our attention to the problem of charged-particle motion in a
simple magnetic tokamak geometry, with circular concentric magnetic surfaces
without Shafranov shift. The magnetic field is represented, using the quasi-
cylindrical coordinates $(r,\vartheta,\varphi)$, by the dimensionless
expression

$\displaystyle{\bf
B}\;\equiv\;\frac{r^{2}}{qh}\,\nabla\vartheta\;+\;\nabla\varphi,$ (41)

where $q(r)=q_{0}+\sigma r^{2}/2$ is the safety factor (we will use $q_{0}=1$
and $\sigma=2$ in our numerical calculations) and $h=1+r\cos{\vartheta}$ is
the normalized distance from the magnetic axis $(r=0)$ to the particle
position in the poloidal plane (which is normalized by the major radius
$R_{0}$ of the magnetic axis). Since the magnetic field (41) is divergenceless

${\bf B}\equiv\nabla\,\mbox{\boldmath$\times$}\,{\bf
A}\;=\;\nabla\Psi\,\mbox{\boldmath$\times$}\,\nabla\vartheta+\nabla\varphi\,\mbox{\boldmath$\times$}\,\nabla\psi,$
(42)

it can be written in terms of the vector potential

${\bf A}\;=\;\Psi\,\nabla\vartheta\;-\;\psi\,\nabla\varphi.$ (43)

By comparing with Eq. (41), we find $\partial\Psi(r,\vartheta)/\partial
r=r/h(r,\vartheta)$ and $d\psi/dr=r/q(r)$, so that the toroidal and poloidal
magnetic fluxes $\Psi$ and $\psi$ are

$\displaystyle\Psi(r,\vartheta)$ $\displaystyle=$
$\displaystyle\int_{0}^{r}\frac{u\,du}{h(u,\vartheta)}=\frac{1}{\cos{\vartheta}}\left(r-\frac{\ln
h(r,\vartheta)}{\cos\vartheta}\right),$ (44) $\displaystyle\psi(r)$
$\displaystyle=$
$\displaystyle\int_{0}^{r}\frac{u\,du}{q(u)}=\frac{1}{\sigma}\,\ln\left(\frac{q(r)}{q_{0}}\right),$
(45)

where we chose $\Psi(0,\vartheta)=0=\psi(0)$. In what follows, we use the
quasi-cylindrical unit vectors $\widehat{\sf
r}=\cos\vartheta\,\widehat{\rho}+\sin\vartheta\,\widehat{\sf z}$,
$\widehat{\vartheta}=\partial\widehat{\sf
r}/\partial\vartheta=-\,\sin\vartheta\,\widehat{\rho}+\cos\vartheta\,\widehat{\sf
z}$, and $\widehat{\varphi}\equiv\widehat{\sf
r}\,\mbox{\boldmath$\times$}\,\widehat{\vartheta}=\widehat{\rho}\,\mbox{\boldmath$\times$}\,\widehat{\sf
z}=\partial\widehat{\rho}/\partial\varphi$, with the quasi-cylindrical
Jacobian ${\cal J}=(\nabla
r\,\mbox{\boldmath$\times$}\,\nabla\vartheta\,\mbox{\boldmath$\cdot$}\,\nabla\varphi)^{-1}=r\,h(r,\vartheta)$.

The magnitude of the magnetic field (41) is

$B(r,\vartheta)\;=\;\frac{1}{h(r,\vartheta)}\sqrt{1+r^{2}/q(r)^{2}}\;\equiv\;\frac{\beta(r)}{h(r,\vartheta)},$
(46)

and the unit vector along the magnetic field is

$\widehat{{\sf
b}}(r,\vartheta)=\frac{r}{q\beta}\;\widehat{\vartheta}\;+\;\frac{1}{\beta}\;\widehat{\varphi}\;\equiv\;b_{\vartheta}(r)\,\nabla\vartheta\;+\;b_{\varphi}(r,\vartheta)\,\nabla\varphi,$
(47)

with $\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf b}}=\nabla
b_{\vartheta}\,\mbox{\boldmath$\times$}\,\nabla\vartheta+\nabla
b_{\varphi}\,\mbox{\boldmath$\times$}\,\nabla\varphi$. Next, we can calculate

$\displaystyle\nabla\ln B$ $\displaystyle=$
$\displaystyle\left(\frac{r\,g}{q\beta^{2}}-\frac{\cos\vartheta}{h}\right)\;\widehat{\sf
r}+\frac{\sin\vartheta}{h}\;\widehat{\vartheta},$ (48) $\kappa$
$\displaystyle=$
$\displaystyle\frac{\sin\vartheta}{h\beta^{2}}\left(\widehat{\vartheta}-\frac{r}{q}\widehat{\varphi}\right)-\frac{\widehat{\sf
r}}{\beta^{2}}\left(\frac{\cos\vartheta}{h}+\frac{r}{q^{2}}\right),$ (49)
$\displaystyle\tau$ $\displaystyle=$ $\displaystyle\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf
b}}\;=\;\frac{1}{\beta^{2}}\left(g\;+\;\frac{1}{hq}\right),$ (50)
$\displaystyle\widehat{{\sf b}}^{\prime}$ $\displaystyle=$
$\displaystyle-\,\frac{\beta^{\prime}}{\beta}\,\widehat{{\sf
b}}\;-\;\left(\varphi^{\prime}\cos\vartheta+\frac{r}{q}\,\vartheta^{\prime}\right)\frac{\widehat{\sf
r}}{\beta}\;+\;\frac{gr^{\prime}}{\beta}\;\widehat{\vartheta}$ (51)
$\displaystyle+\;\frac{\varphi^{\prime}\sin\vartheta}{\beta}\left(\widehat{\vartheta}\;-\;\frac{r}{q}\,\widehat{\varphi}\right),$

where $g(r)\equiv d(r/q)/dr$ and $\beta^{\prime}=(rg/q\beta)\,r^{\prime}$.

###  IV.1 Particle dynamics

The dimensionless particle Lagrangian is

$L=\frac{1}{\epsilon}\left(\Psi\,\vartheta^{\prime}\>-\>\psi\,\varphi^{\prime}\right)\;+\;\frac{1}{2}\left(r^{\prime
2}+r^{2}\,\vartheta^{\prime 2}+h^{2}\,\varphi^{\prime 2}\right),$ (52)

from which we obtain the following equations of motion

$\displaystyle r^{\prime\prime}$ $\displaystyle=$
$\displaystyle\frac{r}{\epsilon}\left(\frac{\vartheta^{\prime}}{h}-\frac{\varphi^{\prime}}{q}\right)+r\,\vartheta^{\prime
2}+h\,\cos\vartheta\,\varphi^{\prime 2},$ (53)
$\displaystyle\vartheta^{\prime\prime}$ $\displaystyle=$
$\displaystyle\frac{-1}{r}\left(\frac{r^{\prime}}{\epsilon\,h}+h\,\sin\vartheta\,\varphi^{\prime
2}+2\,r^{\prime}\,\vartheta^{\prime}\right),$ (54)
$\displaystyle\varphi^{\prime\prime}$ $\displaystyle=$
$\displaystyle\frac{r\,r^{\prime}}{\epsilon\,h^{2}\,q}-\frac{2\,\varphi^{\prime}}{h}\,\left(r^{\prime}\,\cos\vartheta-r\,\vartheta^{\prime}\,\sin\vartheta\right).$
(55)

These dimensionless equations of motion are solved numerically for
$\epsilon=1/100$, with the initial conditions
$(r_{0},\vartheta_{0},\varphi_{0})=(1/2,0,0)$ and
$(r_{0}^{\prime},\vartheta_{0}^{\prime},\varphi_{0}^{\prime})=(0,8,\sqrt{22}/3)$
associated with a dimensionless energy $E=43/2$.

We note that, because of the azimuthal symmetry of the magnetic field (41),
i.e., the components $(B_{\vartheta},B_{\varphi})$ are independent of the
azimuthal angle $\varphi$, the azimuthal canonical angular momentum (5),
expressed in dimensionless form as

$P_{\varphi}\;=\;-\,\frac{1}{\epsilon}\,\psi(r)\>+\>h^{2}(r,\vartheta)\,\varphi^{\prime},$
(56)

is a constant of the motion.

Before moving on to guiding center theory, we calculate the normalized
magnetic moment $\mu=\mu_{0}+\epsilon\,\mu_{1}+\cdots$, where

$\mu_{0}=\left(r^{\prime 2}+r^{2}\omega^{2}\right)/B,$

with

$\omega\;\equiv\;\left(\vartheta^{\prime}\;-\;h\,\varphi^{\prime}/q\right)/\beta,$
(57)

and the dimensionless lowest-order gyroradius is

$\mbox{\boldmath$\rho$}_{0}=-\,\frac{r\omega}{B}\,\widehat{\sf
r}\;+\;\frac{r^{\prime}}{\beta
B}\,\left(\widehat{\vartheta}-\frac{r}{q}\widehat{\varphi}\right).$ (58)

Figure 5 shows that the lowest-order normalized magnetic moment is poorly
conserved, especially as the particle crosses the equatorial plane
$(\vartheta=0)$. We can greatly improve the adiabatic invariance of the
magnetic moment by calculating $\mu_{1}$ from Eq. (10):

$\mu_{1}=\mu_{0}\,G_{1}+\frac{p_{\parallel}^{2}}{B}\,G_{2}+\frac{p_{\parallel}}{B}\,G_{3}-\frac{3\,\mu_{0}\,p_{\parallel}}{2\,B}\,\tau,$
(59)

where

$\displaystyle G_{1}$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\ln
B=\frac{r^{\prime}}{\beta^{2}}\,\sin\vartheta-\frac{r\omega}{B}\left(\frac{r\,g}{q\beta^{2}}-\frac{\cos\vartheta}{h}\right),$
(60) $\displaystyle G_{2}$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\mbox{\boldmath$\kappa$}=\frac{r^{\prime}}{\beta^{2}}\,\sin\vartheta+\frac{r\omega}{\beta^{2}B}\left(\frac{\cos\vartheta}{h}+\frac{r}{q^{2}}\right),$
(61)

and

$\displaystyle
G_{3}\;=\;\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\widehat{{\sf
b}}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{r^{\prime}}{B}\left(\frac{g\,r^{\prime}}{\beta^{2}}+\varphi^{\prime}\sin\vartheta\right)$
(62) $\displaystyle+\;\frac{r\omega}{\beta
B}\left(\frac{r\vartheta^{\prime}}{q}+\varphi^{\prime}\cos\vartheta\right).$

Figure 5 shows lowest order magnetic moment $\mu_{0}$ and the improved
$\mu=\mu_{0}+\epsilon\,\mu_{1}$. This first-order correction causes $\mu\simeq
12.23$ to be a good adiabatic invariant, numerically validating the guiding-
center representation for the particle orbits in the simple tokamak magnetic
field (41). This adiabatic invariance of the magnetic moment justifies our use
of the guiding-center approximation in describing particle orbits in simple
magnetic tokamak geometry.

Figure 5: Plots of the normalized lowest-order magnetic moment $\mu_{0}$
(gray) and the normalized magnetic moment $\mu=\mu_{0}+\epsilon\,\mu_{1}\simeq
12.23$ (black) for one bounce period, with first-order correction $\mu_{1}$
given by Eq. (59). Here, $\epsilon=1/100$ and the normalized kinetic energy is
$E_{0}=21.5$, with the initial conditions
$(r_{0},\vartheta_{0},\varphi_{0})=(1/2,0,0)$ and
$(r_{0}^{\prime},\vartheta_{0}^{\prime},\varphi_{0}^{\prime})=(0,8,\sqrt{22}/3)$.

###  IV.2 Guiding-center dynamics

The guiding-center Lagrangian in simple tokamak geometry is

$L_{\rm
gc}=\frac{1}{\epsilon}\left(\Psi^{*}\,\Theta^{\prime}-\psi^{*}\,\Phi^{\prime}\right)-\epsilon\,J\,R^{\prime}\,{\cal
R}_{R}^{*}-\left(\frac{1}{2}\,P_{\|}^{2}+J\,B\right),$ (63)

where ${\bf X}=(R,\Theta,\Phi)$ denote the guiding-center position in quasi-
cylindrical geometry, $J=\mu/2$ denotes the dimensionless guiding-center
magnetic moment, and $P_{\|}$ denotes the dimensionless guiding-center
momentum parallel to the magnetic field. In addition, the effective poloidal
and toroidal magnetic fluxes

$\displaystyle\Psi^{*}$ $\displaystyle=$
$\displaystyle\Psi\;+\;\epsilon\,P_{\parallel}\,b_{\Theta}\;-\;\epsilon^{2}\,J\,{\cal
R}^{*}_{\Theta},$ (64) $\displaystyle\psi^{*}$ $\displaystyle=$
$\displaystyle\psi\;-\;\epsilon\,P_{\parallel}\,b_{\Phi}\;+\;\epsilon^{2}\,J\,{\cal
R}^{*}_{\Phi}$ (65)

are expressed in terms of the extended guiding-center model (C) in Eq. (12):
$\mbox{\boldmath$\cal R$}^{*}=\mbox{\boldmath$\cal
R$}+\frac{1}{2}\,\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf b}}$. Here,
we calculate the gyrogauge vector $\mbox{\boldmath$\cal R$}=\nabla\widehat{\sf
e}_{1}\,\mbox{\boldmath$\cdot$}\,\widehat{\sf e}_{2}$ by choosing
$\widehat{\sf e}_{1}=\widehat{\sf r}$ and $\widehat{\sf e}_{2}=\widehat{{\sf
b}}\,\mbox{\boldmath$\times$}\,\widehat{\sf r}$, so that we obtain
$\mbox{\boldmath$\cal R$}=\beta^{-1}\nabla\Theta-b_{z}\nabla\Phi$, where
$b_{z}\equiv\widehat{{\sf b}}\,\mbox{\boldmath$\cdot$}\,\widehat{\sf
z}=(R/q\beta)\,\cos\Theta$, and, according to the extended guiding-center
model (C) in Eq. (12), we find

$\displaystyle\mbox{\boldmath$\cal R$}^{*}$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\cal R$}\;+\;\frac{1}{2}\left(\nabla
b_{\Phi}\,\mbox{\boldmath$\times$}\,\nabla\Phi\;+\frac{}{}\nabla
b_{\Theta}\,\mbox{\boldmath$\times$}\,\nabla\Theta\right)$ (66)
$\displaystyle\equiv$ $\displaystyle{\cal R}_{R}^{*}\,\nabla
R\>+\>{\cal{R}}_{\Theta}^{*}\,\nabla\Theta\>+\>{\cal{R}}_{\Phi}^{*}\,\nabla\Phi.$

From this Lagrangian, we obtain Euler-Lagrange equations that can be expressed
as Eqs. (13)-(14), where

${\bf
B}^{*}=\nabla\Psi^{*}\,\mbox{\boldmath$\times$}\,\nabla\Theta-\nabla\psi^{*}\,\mbox{\boldmath$\times$}\,\nabla\Phi+\nabla\chi\,\mbox{\boldmath$\times$}\,\nabla
R,$ (67)

with $\chi\equiv-\,\epsilon^{2}\,J{\cal R}_{R}^{*}$ and
$B_{\|}^{*}=\widehat{{\sf b}}\,\mbox{\boldmath$\cdot$}\,{\bf B}^{*}$. The
guiding-center equations are thus expressed as

$\displaystyle R^{\prime}$ $\displaystyle=$
$\displaystyle-\;\frac{P_{\|}}{{\cal J}_{\rm
gc}}\;\frac{\partial\psi^{*}}{\partial\Theta}\;-\;\frac{\epsilon\,J}{{\cal
J}_{\rm gc}}\;b_{\Phi}\frac{\partial B}{\partial\Theta},$ (68)
$\displaystyle\Theta^{\prime}$ $\displaystyle=$
$\displaystyle\frac{P_{\|}}{{\cal J}_{\rm
gc}}\;\frac{\partial\psi^{*}}{\partial R}\;+\;\frac{\epsilon\,J}{{\cal J}_{\rm
gc}}\;b_{\Phi}\frac{\partial B}{\partial R},$ (69)
$\displaystyle\Phi^{\prime}$ $\displaystyle=$
$\displaystyle\frac{P_{\|}}{{\cal J}_{\rm
gc}}\left(\frac{\partial\Psi^{*}}{\partial
R}+\epsilon^{2}J\;\frac{\partial{\cal
R}_{R}^{*}}{\partial\Theta}\right)\;-\;\frac{\epsilon\,J}{{\cal J}_{\rm
gc}}\;b_{\Theta}\frac{\partial B}{\partial R},$ (70) $\displaystyle
P_{\|}^{\prime}$ $\displaystyle=$ $\displaystyle-\;\frac{J}{{\cal J}_{\rm
gc}}\left(\frac{\partial\psi^{*}}{\partial
R}\;\frac{\partial\overline{B}}{\partial\Theta}\;-\;\frac{\partial\psi^{*}}{\partial\Theta}\;\frac{\partial
B}{\partial R}\right),$ (71)

where ${\cal J}_{\rm gc}\equiv{\cal J}\,B_{\|}^{*}$ combines the quasi-
cylindrical Jacobian ${\cal J}$ and the guiding-center Jacobian $B_{\|}^{*}$.
We note that these equations exactly conserve the guiding-center energy ${\cal
E}=P_{\|}^{2}/2+JB$ and the guiding-center canonical azimuthal angular
momentum

$P_{{\rm
gc}\Phi}\;=\;-\;\frac{1}{\epsilon}\;\psi^{*}\;=\;-\,\frac{1}{\epsilon}\;\psi\;+\;P_{\|}\,b_{\Phi}\;-\;\epsilon\,J\,{\cal
R}_{\Phi}^{*}.$ (72)

The dimensionless guiding-center equations are solved numerically for
$\epsilon=1/100$ with the initial conditions
$(R_{0},\Theta_{0},\Phi_{0})=(0.5380,0,0)$, which takes into account the
radial shift (16) from the initial particle position
$(r_{0},\varphi_{0},z_{0})=(0.5,0,0)$, and
$P_{\|0}=\sqrt{E-J\,B(R_{0},\Theta_{0})}$.

Figure 6: Plots of the particle radial position (gray) and the guiding-center
radial position (black) during the first two bounce periods. Figure 7: Plots
of the particle position (gray) and the guiding-center position (black)
projected into the poloidal plane during the first two bounce periods. Figure
8: Plots of the three-dimensional particle orbit (gray) and the guiding-center
orbit (black) during the first two bounce periods.

Figures 6-8 show plots of the particle position (gray) and the guiding-center
position (black), obtained from the extended guiding-center model C in Eq.
(12), during the first bounce periods. Figure 6 shows plots of the particle
radial position (gray) and the guiding-center radial position (black) during
the first two bounce periods, Fig. 7 shows the classic “closed” guiding-center
banana orbit (black) projected onto the poloidal plane (at constant toroidal
angle), and Fig. 8 shows that the three-dimensional guiding-center orbit
(black) follows the three-dimensional particle orbit (gray) very well over two
bounce periods.

###  IV.3 Higher-order guiding-center orbits

Figure 9: Plots of the particle position (gray) and the guiding-center
positions projected into the poloidal plane near the upper turning point,
according to three different guiding-center models in Eq. (12): (black, solid)
Extended guiding-center model (C); (black, dashed) Standard guiding-center
model (B); and (black, dotted) Truncated guiding-center model (A).

We now make a few remarks on the three guiding-center models presented in Eq.
(12). In order to compare their effectiveness at approximating the particle
orbit, we need to ensure that the initial conditions for these guiding-center
orbits are consistent with the initial conditions for the particle orbit. This
consistency is achieved by connecting the initial conditions through the
guiding-center transformation (16).

In what follows, the guiding-center orbit A is generated from the initial
condition obtained from the lowest-order relation ${\bf X}={\bf
x}-\epsilon\,\mbox{\boldmath$\rho$}_{0}$ (i.e., the initial condition only
takes into account the lowest-order guiding-center transformation), which
yields the initial radial position $R_{0}|_{A}=0.5396$ from the particle
initial conditions. The initial conditions for the guiding-center orbits B and
C, on the other hand, are distinguished by the guiding-center polarization
correction $G_{2}^{\bf x}|_{\rm pol}=(J/2m\Omega)\,\mbox{\boldmath$\kappa$}$
for the extended guiding-center model C Tronko and Brizard (2015), while
$G_{2}^{\bf x}|_{\rm pol}=0$ for the standard guiding-center model B
Littlejohn (1983); Brizard (1989). Hence, we use the initial radial positions
$R_{0}|_{C}=0.5380$ and $R_{0}|_{B}=0.5384$, which are calculated when the
particle initial conditions are inserted in the guiding-center transformation
(16).

Figure 9 shows the particle position (gray) and the guiding-center positions
(labeled A, B, and C) projected into the poloidal plane near the upper turning
point. Because the guiding-center orbits have different guiding-center
toroidal canonical angular momenta (65), Fig. 9 shows that each guiding-center
turning point occurs on a different magnetic surface. In addition, while the
standard guiding-center orbit B slightly overshoots the particle center of
gyration at the turning point, the truncated guiding-center orbit A largely
undershoots the particle center of gyration. We note that the three guiding-
center orbits are nearly indistinguishable away from the turning-point regions
(see Fig. 7).

Figure 10 shows the three-dimensional particle orbit (gray) and three guiding-
center orbits (labeled A, B, and C) during the first two bounce periods. Here,
we see that the lowest-order truncated guiding-center orbit A clearly
separates from the particle orbit (i.e., the guiding-center orbit is located
well outside of the particle’s gyration radius), while the two higher-order
guiding-center orbits B and C are still nearly indistinguishable over the
first two bounce periods, except when the standard guiding-center orbit B
overshoots the turning point. We note that the standard guiding-center orbit B
noticeably separates from the particle orbit during the next two bounce
periods, while the extended guiding-center orbit C still remains at the center
of the particle orbit.

Figure 10: Plots of the three-dimensional particle orbit (gray) and the
guiding-center orbits (labeled A, B, and C) during two bounce periods,
according to three different guiding-center models in Eq. (12): (black, solid)
Extended guiding-center model (C); (black, dashed) Standard guiding-center
model (B); and (black, dotted) Truncated guiding-center model (A). Here, the
truncated guiding-center orbit A has clearly separated from the particle
orbit, while the higher-order guiding-center orbits B and C are nearly
indistinguishable, except when the standard guiding-center orbit B overshoots
the turning point.

###  IV.4 Validity of the guiding-center representation in simple tokamak
magnetic geometry

We now show that the guiding-center approximation is faithful to the particle
motion in simple tokamak magnetic geometry by showing that the guiding-center
pull-back ${\sf T}_{\rm gc}P_{{\rm gc}\Phi}=P_{\varphi}$ of the guiding-center
canonical azimuthal angular momentum is equal to the particle canonical
azimuthal angular momentum. Up to second order in $\epsilon$, the guiding-
center pull-back ${\sf T}_{\rm gc}P_{{\rm gc}\Phi}$ is expressed as

$\displaystyle{\sf T}_{\rm gc}P_{{\rm gc}\Phi}$ $\displaystyle=$
$\displaystyle-\frac{1}{\epsilon}\,\psi\;+\;\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\psi\;+\;p_{\|}\;b_{\varphi}$
$\displaystyle+\;\epsilon\;\mbox{\boldmath$\rho$}_{1}\,\mbox{\boldmath$\cdot$}\,\nabla\psi\;-\;\frac{\epsilon}{2}\,\mbox{\boldmath$\rho$}_{0}\mbox{\boldmath$\rho$}_{0}:\nabla\nabla\psi$
$\displaystyle+\;\epsilon\left(G_{1}^{p_{\|}}\,b_{\varphi}-p_{\|}\,\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla
b_{\varphi}\right)\;-\;\epsilon\,J\,{\cal R}_{\varphi}^{*},$

where the first-order gyroradius correction $\mbox{\boldmath$\rho$}_{1}$ is
given by Eq. (19) and the first-order correction to the guiding-center
parallel momentum is given by Eq. (83):

$G_{1}^{p_{\|}}\;=\;-\;\frac{p_{\|}}{2}\,G_{2}\;+\;\frac{3}{2}\,J\;\tau\;-\;\frac{1}{2}\,G_{3},$
(74)

where $(\tau,G_{2},G_{3})$ are defined in Eqs. (50) and (61)-(62).

First, using $\mbox{\boldmath$\rho$}_{0}=\widehat{{\sf
b}}\,\mbox{\boldmath$\times$}\,{\bf x}^{\prime}/B$ and the simple-tokamak
identity ${\bf B}\,\mbox{\boldmath$\times$}\,\partial{\bf
x}/\partial\varphi=\nabla\psi$, we find

$\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\psi\;=\;\frac{\widehat{{\sf
b}}}{B}\,\mbox{\boldmath$\times$}\,{\bf
x}^{\prime}\,\mbox{\boldmath$\cdot$}\,{\bf
B}\,\mbox{\boldmath$\times$}\,\frac{\partial\bf
x}{\partial\varphi}\;=\;h^{2}\varphi^{\prime}\;-\;p_{\|}\,b_{\varphi},$

so that, at first order (i.e., zeroth order in magnetic-field nonuniformity),
we find

$\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\psi\;+\;p_{\|}\;b_{\varphi}\;=\;h^{2}\,\varphi^{\prime}.$
(75)

Hence, we now need to show that, at second order (i.e., first order in
magnetic-field nonuniformity), we find the identity

$\displaystyle J{\cal R}_{\varphi}^{*}$ $\displaystyle\equiv$
$\displaystyle\mbox{\boldmath$\rho$}_{1}\,\mbox{\boldmath$\cdot$}\,\nabla\psi\;-\;\frac{1}{2}\mbox{\boldmath$\rho$}_{0}\mbox{\boldmath$\rho$}_{0}:\nabla\nabla\psi$
(76)
$\displaystyle+\;G_{1}^{p_{\|}}b_{\varphi}\;-\;p_{\|}\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla
b_{\varphi},$

so that Eq. (IV.4) becomes

${\sf T}_{\rm gc}P_{{\rm
gc}\Phi}\;=\;-\frac{1}{\epsilon}\,\psi\;+\;h^{2}\,\varphi^{\prime}\;\equiv\;P_{\varphi},$
(77)

which guarantees the validity of the guiding-center representation in simple
tokamak magnetic geometry. A complete proof of the identity (76) is given in
App. C.

##  V Summary

In previous work Brizard (2017), we showed that the guiding-center
approximation was valid in a straight magnetic field with constant
perpendicular magnetic gradient, even in the presence of strong gradients. In
addition, based on the existence of an exact analytical solution for the
particle orbits, this work also confirmed that the guiding-center polarization
corresponded exactly with an orbit-averaged particle displacement.

In the present work, we extended our investigation of the validity of the
guiding-center approximation in describing charged single-particle motion in a
nonuniform magnetic field. Here, we considered regular particle orbits in
azimuthally symmetric magnetic mirror geometry and simple tokamak magnetic
geometry, in which the azimuthal angular canonical momentum is conserved and
the guiding-center magnetic moment is an adiabatic invariant. We successfully
validated the guiding-center approximation in describing particle motion in an
azimuthally symmetric magnetic field provided higher-order guiding-center
corrections are taken into account, which had already been noted for the case
of an axisymmetric tokamak magnetic field Belova et al. (2003). In particular,
the guiding-center polarization correction in the guiding-center azimuthal
angular canonical momentum, not taken into account in the standard guiding-
center approximation Littlejohn (1983); Brizard (1989), proved crucial in
establishing a faithful guiding-center representation for regular particle
orbits in axisymmetric magnetic geometry.

Lastly, we note that the truncated guiding-center model White (2014) is used
extensively in guiding-center particle simulations, despite the fact that it
is not as faithful to particle orbits as higher-order guiding-center models.
In most applications, however, the truncated guiding-center model is either
used to analyze particle orbits in the presence of perturbed electric and/or
magnetic fields, or as the unperturbed component for the gyrocenter orbit used
in nonlinear gyrokinetic theory Brizard (1989). Future work may look into the
issue of faithfulness for these applications. In addition, the faithfulness of
the guiding-center representation for particle orbits in non-axisymmetric
magnetic geometries may be explored.

##  Appendix A Guiding-center Transformation

The standard expression for the first-order correction to the guiding-center
magnetic moment is Littlejohn (1983); Tronko and Brizard (2015)

$G_{1}^{\mu}=\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\left(\mu_{0}\,\nabla\ln
B+\frac{p_{\|}^{2}\,\mbox{\boldmath$\kappa$}}{m\,B}\right)-\mu_{0}\;\varrho_{\|}\,(\tau+\alpha_{1}),$
(78)

where $\varrho_{\|}\equiv p_{\|}/(m\Omega)$ and
$\alpha_{1}\equiv-\,\frac{1}{2}\,(\widehat{\bot}\widehat{\rho}+\widehat{\rho}\widehat{\bot}):\nabla\widehat{{\sf
b}}$ is constructed from the gyrorangle-dependent unit vectors
$\widehat{\bot}\equiv\partial\widehat{\rho}/\partial\zeta=\widehat{\rho}\,\mbox{\boldmath$\times$}\,\widehat{{\sf
b}}$. Using the identity

$\alpha_{1}\;=\;\frac{1}{2}\,\tau\;-\;\widehat{\bot}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\widehat{\rho},$ (79)

we obtain

$\tau+\alpha_{1}\;=\;\frac{3}{2}\,\tau\;-\;\widehat{\bot}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\widehat{\rho}.$ (80)

Next, we write

$\displaystyle\mu_{0}\;\varrho_{\|}\,\left(\widehat{\bot}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\widehat{\rho}\right)$ $\displaystyle=$
$\displaystyle\frac{p_{\|}}{2B}\;\left({\bf
v}_{\bot}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\mbox{\boldmath$\rho$}_{0}\right)$ (81)
$\displaystyle=$ $\displaystyle\frac{p_{\|}}{2B}\;\left(\frac{d\widehat{{\sf
b}}}{dt}\;-\;v_{\|}\,\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}\right)\,\mbox{\boldmath$\cdot$}\,\mbox{\boldmath$\rho$}_{0}$
$\displaystyle=$ $\displaystyle\left(\frac{p_{\|}}{2B}\;\frac{d\widehat{{\sf
b}}}{dt}\;-\;\frac{p_{\|}^{2}\mbox{\boldmath$\kappa$}}{2\,mB}\right)\,\mbox{\boldmath$\cdot$}\,\mbox{\boldmath$\rho$}_{0},$

which yields Eq. (10):

$\displaystyle G_{1}^{\mu}$ $\displaystyle=$
$\displaystyle\left(\mu_{0}\nabla\ln
B+\frac{p_{\|}^{2}\,\mbox{\boldmath$\kappa$}}{2\,mB}\right)\,\mbox{\boldmath$\cdot$}\,\mbox{\boldmath$\rho$}_{0}-\frac{3}{2}\,\mu_{0}\left(\frac{p_{\|}\,\tau}{m\Omega}\right)$
(82) $\displaystyle+\;\frac{p_{\|}}{2B}\;\frac{d\widehat{{\sf
b}}}{dt}\,\mbox{\boldmath$\cdot$}\,\mbox{\boldmath$\rho$}_{0}.$

Using the same identity (80), the standard expression for the first-order
correction to the guiding-center parallel momentum Littlejohn (1983); Tronko
and Brizard (2015) is replaced with the new expression

$\displaystyle G_{1}^{p_{\|}}$ $\displaystyle=$
$\displaystyle-\;p_{\|}\,\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\mbox{\boldmath$\kappa$}\;+\;\frac{\mu_{0}B}{\Omega}\;(\tau+\alpha_{1})$
(83) $\displaystyle=$
$\displaystyle-\;\frac{p_{\|}}{2}\,\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\mbox{\boldmath$\kappa$}\;+\;\frac{3}{2}\,\frac{\mu_{0}B}{\Omega}\;\tau\;-\;\frac{m}{2}\,\frac{d\widehat{{\sf
b}}}{dt}\,\mbox{\boldmath$\cdot$}\,\mbox{\boldmath$\rho$}_{0}.$

##  Appendix B Guiding-center Polarization

The guiding-center polarization was calculated directly from the guiding-
center transformation in our previous works Brizard (2013); Tronko and Brizard
(2015). It is formally defined by the multipole expansion

$\mbox{\boldmath$\pi$}_{\rm gc}\;\equiv\;e\,\langle\mbox{\boldmath$\rho$}_{\rm
gc}\rangle\;-\;\nabla\,\mbox{\boldmath$\cdot$}\,\left(\frac{e}{2}\;\langle\mbox{\boldmath$\rho$}_{\rm
gc}\mbox{\boldmath$\rho$}_{\rm gc}\rangle\right)\;+\;\cdots,$ (84)

where the dipole and quadrupole moments are shown here, while
$\mbox{\boldmath$\rho$}_{\rm gc}\equiv{\sf T}_{\rm gc}^{-1}{\bf x}-{\bf X}$ is
the guiding-center gyroradius. We note that the guiding-center gyroradius
$\mbox{\boldmath$\rho$}_{\rm gc}$ is related to the particle gyroradius
$\mbox{\boldmath$\rho$}\equiv{\bf x}-{\sf T}_{\rm gc}{\bf X}$ by the identity
$\mbox{\boldmath$\rho$}_{\rm gc}\equiv{\sf T}_{\rm
gc}^{-1}\mbox{\boldmath$\rho$}$. Using the guiding-center transformation
presented by Tronko and Brizard Tronko and Brizard (2015), we find the dipole
contribution

$\displaystyle\langle\mbox{\boldmath$\rho$}_{\rm gc}\rangle$ $\displaystyle=$
$\displaystyle\epsilon^{2}\left(\langle\mbox{\boldmath$\rho$}_{1}\rangle-\frac{p_{\|}^{2}\mbox{\boldmath$\kappa$}}{m^{2}\Omega^{2}}\right)$
(85) $\displaystyle-\;\frac{\epsilon^{2}\mu
B}{m\Omega^{2}}\left[2\,\nabla_{\bot}\ln
B\;+\frac{}{}\left(\nabla\,\mbox{\boldmath$\cdot$}\,\widehat{{\sf
b}}\right)\widehat{{\sf b}}\right],$

and the quadrupole contribution

$\displaystyle-\,\nabla\,\mbox{\boldmath$\cdot$}\,\left(\frac{1}{2}\;\langle\mbox{\boldmath$\rho$}_{\rm
gc}\mbox{\boldmath$\rho$}_{\rm gc}\rangle\right)$ $\displaystyle=$
$\displaystyle-\,\nabla\,\mbox{\boldmath$\cdot$}\,\left(\frac{\epsilon^{2}\mu
B}{2m\Omega^{2}}\;(\mathbb{I}-\widehat{{\sf b}}\widehat{{\sf b}})\right)$
$\displaystyle=$ $\displaystyle\frac{\epsilon^{2}\mu
B}{2m\Omega^{2}}\;\nabla_{\bot}\ln B$ $\displaystyle+\;\frac{\epsilon^{2}\mu
B}{2m\Omega^{2}}\left[\mbox{\boldmath$\kappa$}\;+\frac{}{}\left(\nabla\,\mbox{\boldmath$\cdot$}\,\widehat{{\sf
b}}\right)\widehat{{\sf b}}\right],$

which both appear at $\epsilon^{2}$ at their lowest orders (i.e., first order
in magnetic-field nonuniformity).

Hence, the guiding-center polarization (84) is expressed as

$\displaystyle\mbox{\boldmath$\pi$}_{\rm gc}$ $\displaystyle=$
$\displaystyle\epsilon^{2}e\left(\langle\mbox{\boldmath$\rho$}_{1}\rangle_{\rm
pol}\;+\;\frac{\mu B}{2\,m\Omega^{2}}\;\mbox{\boldmath$\kappa$}\right)$ (87)
$\displaystyle+\;\frac{\epsilon^{2}e\widehat{{\sf
b}}}{\Omega}\,\mbox{\boldmath$\times$}\,\left[\frac{\widehat{{\sf
b}}}{m\Omega}\,\mbox{\boldmath$\times$}\,\left(\mu\,\nabla
B+\frac{p_{\|}^{2}}{m}\,\mbox{\boldmath$\kappa$}\right)\right],$

where $\langle\mbox{\boldmath$\rho$}_{1}\rangle_{\rm pol}$ is the polarization
correction not included in the standard guiding-center transformation Brizard
(1989); Cary and Brizard (2009). We, therefore, recover the standard guiding-
center polarization $\mbox{\boldmath$\pi$}_{\rm
gc}\equiv\epsilon^{2}(e\widehat{{\sf
b}}/\Omega)\,\mbox{\boldmath$\times$}\,d{\bf X}/dt$ Kaufman (1986) only if we
choose

$\langle\mbox{\boldmath$\rho$}_{1}\rangle_{\rm pol}\;=\;-\,\frac{\mu
B}{2\,m\Omega^{2}}\;\mbox{\boldmath$\kappa$},$ (88)

which appears as the first term on the right side of Eq. (19). With this
choice, the magnetic vector potential ${\bf A}^{*}$ defined in Eq. (17)
becomes

$\displaystyle\frac{e{\bf A}^{*}}{\epsilon c}$ $\displaystyle=$
$\displaystyle\frac{e{\bf A}}{\epsilon c}+P_{\|}\,\widehat{{\sf
b}}\;-\;\epsilon J\left(\mbox{\boldmath$\cal R$}+\frac{1}{2}\tau\widehat{{\sf
b}}\right)\;+\;\frac{m\Omega}{\epsilon}\,\langle\mbox{\boldmath$\rho$}_{1}\rangle_{\rm
pol}$ (89) $\displaystyle=$ $\displaystyle\frac{e{\bf A}}{\epsilon
c}+P_{\|}\,\widehat{{\sf b}}\;-\;\epsilon J\left(\mbox{\boldmath$\cal
R$}+\frac{1}{2}\,\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf b}}\right),$

where the standard correction $\frac{1}{2}\tau\widehat{{\sf b}}$ Brizard
(1989); Cary and Brizard (2009) is replaced with the correction
$\frac{1}{2}\,\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf b}}$ Brizard
(2013); Tronko and Brizard (2015).

Lastly, it is important to keep in mind that the guiding-center polarization
discussed here occurs in the absence of an external electric field and is
simply due to magnetic-field non-uniformity. Since the guiding-center
polarization effect is inversely proportional the the gyrofrequency, it is an
important effect for ions. The electric field generated by the charge
separation associated with the guiding-center polarization, therefore,
requires a self-consistent treatment that must include an electric field as an
integral part of the guiding-center formulation Brizard (1995). This self-
consistent analysis, however, is outside the scope of this paper.

##  Appendix C Proof of Identity (76)

In this Appendix, we proceed with a proof of the identity (76). Using
$(\mbox{\boldmath$\kappa$}/B)\,\mbox{\boldmath$\cdot$}\,\nabla\psi=\tau\,b_{\varphi}-2\,\left(b_{z}+{\cal
R}_{\varphi}^{*}\right)$, we write

$\displaystyle\mbox{\boldmath$\rho$}_{1}\,\mbox{\boldmath$\cdot$}\,\nabla\psi$
$\displaystyle=$
$\displaystyle-\,\frac{1}{2}\,J\left[\tau\,b_{\varphi}\;-\frac{}{}2\,\left(b_{z}+{\cal
R}_{\varphi}^{*}\right)\right]$ (90)
$\displaystyle+\;\left(\frac{1}{2}\,G_{1}-\frac{p_{\|}}{B}\,\tau\right)\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\psi,$

Next, since $\nabla\psi=(r/q)\,\widehat{\sf r}$, we find

$\displaystyle\nabla\nabla\psi$ $\displaystyle=$ $\displaystyle
g\;\widehat{\sf r}\,\widehat{\sf
r}\;+\;\frac{1}{q}\,\widehat{\vartheta}\widehat{\vartheta}\;+\;B\,b_{z}\;\widehat{\varphi}\widehat{\varphi}$
(91) $\displaystyle=$ $\displaystyle B\,b_{z}\;{\sf
I}\;+\;\left(g\;-\;\frac{r\cos\vartheta}{qh}\right)\;\widehat{\sf
r}\,\widehat{\sf r}\;+\;\frac{\widehat{\vartheta}\widehat{\vartheta}}{qh},$

so that

$\frac{1}{2}\mbox{\boldmath$\rho$}_{0}\mbox{\boldmath$\rho$}_{0}:\nabla\nabla\psi=Jb_{z}+\frac{1}{2}\left(g-\frac{r\cos\vartheta}{qh}\right)\rho_{0r}^{2}+\frac{\rho_{0\vartheta}^{2}}{2\,qh},$
(92)

where $\rho_{0r}=-\,r\omega/B$, $\rho_{0\vartheta}=r^{\prime}/\beta B$, and we
used the lowest-order expression

$J\;=\;\frac{B}{2}\left(\rho_{0r}^{2}\;+\frac{}{}\beta^{2}\rho_{0\vartheta}^{2}\right)\;=\;\frac{1}{2B}\left(r^{\prime
2}\;+\;r^{2}\omega^{2}\right),$

where $\rho_{0\varphi}=-(r/q)\,\rho_{0\vartheta}$. With these expressions, and
using Eq. (74), the identity (76) becomes

$\displaystyle J\,{\cal R}_{\varphi}^{*}$ $\displaystyle=$
$\displaystyle-\,\frac{J}{2}\left[\tau\,b_{\varphi}-2\frac{}{}\left(b_{z}+{\cal
R}_{\varphi}^{*}\right)\right]\;-\;J\,b_{z}\;+\;\frac{3}{2}\,J\;\tau\,b_{\varphi}$
$\displaystyle+\;\frac{p_{\|}}{B}\,\left(G_{1}\;-\;\frac{1}{2}\,G_{2}\;-\;\tau\;\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\psi\right)\;-\;\frac{G_{3}}{2B}$
$\displaystyle-\;\frac{\rho_{0\vartheta}^{2}}{2\,hq}-\frac{1}{2}\left(g-\frac{r\cos\vartheta}{hq}\right)\rho_{0r}^{2}+\frac{1}{2}\,G_{1}\;\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\psi,$

which, after cancellations, yields an expression for $G_{3}$:

$\displaystyle G_{3}$ $\displaystyle=$ $\displaystyle
2\,J\,\tau\;+\;p_{\|}\left(2\,G_{1}\;-\;G_{2}\;-\frac{}{}2\,\tau\;\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\psi\right)$
$\displaystyle+\;G_{1}\,B\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\psi\;-\;\frac{B\,\rho_{0\vartheta}^{2}}{hq}\;-\;\left(\beta^{2}\,\tau-\frac{1}{q}\right)B\rho_{0r}^{2},$

where we used

$g-\frac{r\cos\vartheta}{hq}\;=\;g+\frac{1}{q}\left(\frac{1}{h}-1\right)\;=\;\beta^{2}\,\tau-\frac{1}{q},$

which follows from the definition (50) for $\tau$.

We now compare Eq. (C) with Eq. (62), which requires expressing
$(\vartheta^{\prime},\varphi^{\prime})$ in terms of $(p_{\|},\omega)$, where
$\omega$ is defined in Eq. (57) and
$p_{\|}=\varphi^{\prime}/B+r^{2}\vartheta^{\prime}/(q\beta)$ is the lowest-
order dimensionless particle parallel momentum. Hence, after substituting

$\displaystyle\vartheta^{\prime}$ $\displaystyle=$
$\displaystyle(\omega+p_{\|}/q)/\beta,$ $\displaystyle\varphi^{\prime}$
$\displaystyle=$ $\displaystyle(p_{\|}-r^{2}\omega/q)/(h\beta),$

into Eq. (62), we obtain a second expression for $G_{3}$:

$\displaystyle G_{3}$ $\displaystyle=$ $\displaystyle
2\,J\,\tau\;+\;p_{\|}\,G_{2}\;+\;G_{1}\,B\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\psi$
$\displaystyle-\;\frac{2J}{hq\,\beta^{2}}\;-\;\frac{r^{2}\omega^{2}}{B}\left[g\;-\;\frac{(h\beta^{2}+1-\beta^{2})}{hq\,\beta^{2}}\right],$

where we used $1+\beta^{2}r\cos\vartheta=h\beta^{2}+1-\beta^{2}$. By comparing
Eqs. (C) and (C), we obtain the following expression

$p_{\|}\left(G_{1}\;-\;G_{2}\;-\frac{}{}\tau\;\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\psi\right)\;=\;0,$
(96)

after carrying out several cancellations on the right side of Eq. (96).
Lastly, using Eqs. (60)-(61), we obtain

$\displaystyle G_{1}-G_{2}$ $\displaystyle=$
$\displaystyle-\,\frac{r\omega}{B}\left(\frac{rg}{q\beta^{2}}-\frac{\cos\vartheta}{h}\right)-\frac{r\omega}{B\beta^{2}}\left(\frac{\cos\vartheta}{h}-\frac{r}{q^{2}}\right)$
(97) $\displaystyle=$
$\displaystyle-\,\tau\;\frac{r^{2}\omega}{qB}\;\equiv\;\tau\,\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\psi,$

which confirms Eq. (96) and completes the proof of the identity (76).

###### Acknowledgements.

The present work was supported by the National Science Foundation grant
PHY-2206302.  Data Availability Statement The Mathematica code used to
generate the plots in the present manuscript is available upon request.

## References

  * Littlejohn (1983) R. Littlejohn, J. Plasma Phys. 29, 111 (1983).
  * Brizard (1989) A. Brizard, J. Plasma Phys. 41, 541 (1989).
  * Cary and Brizard (2009) J. R. Cary and A. J. Brizard, Rev. Mod. Phys. 81, 693 (2009).
  * Tronko and Brizard (2015) N. Tronko and A. J. Brizard, Phys. Plasmas 22, 112507 (2015).
  * Qin and Guan (2008) H. Qin and X. Guan, Phys. Rev. Lett. 100, 035006 (2008).
  * White (2014) R. B. White, _The Theory Of Toroidally Confined Plasmas, 3rd ed._ (Imperial College Press, 2014).
  * Burby and Ellison (2017) J. Burby and C. Ellison, Phys. Plasmas 24, 110703 (2017).
  * Albert et al. (2020) C. G. Albert, S. V. Kasilov, and W. Kernbichler, J. Comp. Phys. 403, 109065 (2020).
  * Bierwage et al. (2022) A. Bierwage, R. B. White, and A. Matsuyama, Phys. Plasmas 29, 113905 (2022).
  * Belova et al. (2003) E. V. Belova, N. N. Gorelenkov, and C. Z. Cheng, Phys. Plasmas 10, 3240 (2003).
  * Burby and Qin (2012) J. Burby and H. Qin, Phys. Plasmas 19, 052106 (2012).
  * Burby (2020) J. W. Burby, J. Math. Phys. 61, 012703 (2020).
  * White and Chance (1984) R. White and M. Chance, Phys. Fluids 27, 2455 (1984).
  * Brizard (2013) A. J. Brizard, Phys. Plasmas 20, 092309 (2013).
  * Kaufman (1986) A. N. Kaufman, Phys. Fluids 29, 1736 (1986).
  * Kruskal (1958) M. Kruskal, Princeton University (NJ) Project Mattehorn (1958).
  * Brizard (2017) A. J. Brizard, Phys. Plasmas 20, 042115 (2017).
  * Brizard (1995) A. J. Brizard, Phys. Plasmas 2, 459 (1995).