1. Introduction
An accelerated charged particle emits electromagnetic radiation and loses energy and momentum in reaction, in accordance with the radiation–reaction force (Abraham Reference Abraham1905; Lorentz Reference Lorentz1936; Dirac Reference Dirac1938; Pauli Reference Pauli1958; Landau & Lifshitz Reference Landau and Lifshitz1975). In a magnetized plasma, the radiation resulting from the gyration around the field lines is often referred to as the synchrotron emission. The radiation–reaction force (RR-force) increases with the particle energy and the accelerating force, which itself depends on the velocity in the case of the Lorentz force. The effect of the radiation on the particle motion can be significant for very energetic particles. In particular, it contributes to limit the energy reached by runaway electrons in tokamak plasmas (Bakhtiari et al. Reference Bakhtiari, Kramer, Takechi, Tamai, Miura, Kusama and Kamada2005), and explains the observation of an elevated critical electric field for runaway-electron generation (Stahl et al. Reference Stahl, Hirvijoki, Decker, Embréus and Fülöp2015).
One of the consequences of magnetic field non-uniformity in axisymmetric configurations such as dipole or tokamak fields is the superposition of three periodic motions in the particle trajectories, namely the gyromotion, bounce or transit motion, and drift precession. The acceleration associated with each periodic motion would in turn contribute to the radiation losses. Another consequence of the magnetic non-uniformity is that the RR-force could induce transport of particles across the magnetic flux surfaces. As the RR-force increases with increasing particle energy while the collisional force decreases with velocity, radial transport associated with the RR-force could overcome collisional transport at high relativistic energies.
In the presence of a weak magnetic non-uniformity such that the gyroradius of the charged particle is much smaller than the magnetic field non-uniformity length scale, the magnetic moment,
${\it\mu}$
, is an adiabatic invariant of the Hamiltonian particle motion. In this case, it is often useful to separate the gyromotion from the rest of the particle motion, and study longer time scales. Making use of the adiabatic invariance, Lie-transform perturbation methods can be used to eliminate the fast gyromotion and to derive the underlying guiding-centre dynamics. Cumulative effects of the gyromotion in a non-uniform magnetic field, such as the mirror force and magnetic drifts, are entirely retained in the guiding-centre dynamics, which is typically much easier to compute than the full particle dynamics. This approach is one of the classical results in modern plasma physics (Littlejohn Reference Littlejohn1983) and has been summarized in a review paper by Cary & Brizard (Reference Cary and Brizard2009).
Previous attempts to include the effect of magnetic field non-uniformity in the RR-force were made without going through a proper guiding-centre transformation. In the first paper, by Andersson, Helander & Eriksson (Reference Andersson, Helander and Eriksson2001), a contribution from the magnetic field-line curvature was argued in addition to the uniform-field formulation. This corrective term, which is of second order in magnetic-field non-uniformity, is interesting since it does not vanish for
${\it\mu}\rightarrow 0$
. Their approach, however, neglects a number of terms that are of first order in magnetic-field non-uniformity. Furthermore, some of the first-order terms and all second-order terms in the Hamiltonian motion are neglected which leads to inconsistent treatment of the dissipative and Hamiltonian guiding-centre dynamics. One can understand how keeping only the leading-order terms is necessary for conducting analytical studies of complex phenomena, as done in Andersson et al. (Reference Andersson, Helander and Eriksson2001), but in particle-following applications where the trajectory is calculated numerically one should be as consistent as possible to prevent e.g. numerical drift of the energy due to improper equations. In papers by Guan, Qin & Fisch (Reference Guan, Qin and Fisch2010) and Liu et al. (Reference Liu, Qin, Fisch, Teng and Wang2014), a rather different approach is adopted by treating the RR-force as an ‘effective electric field’ which is then added into the guiding-centre Lagrangian as a time-depending perturbation in the vector potential. This is not correct: the RR-force is of dissipative nature and no practical Lagrangian formulation exists within the framework of classical electrodynamics (here we do not discuss the quantum mechanical treatments). Moreover, the ‘effective field’ in these papers (Guan et al.
Reference Guan, Qin and Fisch2010; Liu et al.
Reference Liu, Qin, Fisch, Teng and Wang2014) is given only for a simplified toroidal geometry and if one calculates the corresponding equations of motion using the Euler–Lagrange equation, the result does not give the ‘effective electric field’ that they start with.
As the radiatiative momentum losses, however, are important for the dynamics of relativistic charged particles, we see that a guiding-centre description that is consistent with the Hamiltonian formalism is necessary. In the present paper, we derive the guiding-centre RR-force in a weakly non-uniform magnetic field using Lie-transform perturbation methods. In § 2, we first introduce the particle phase-space RR-force and its expression in a non-uniform magnetic field. A general method for including non-Hamiltonian forces into guiding-centre formalism is described in § 3. The guiding-centre transformation is carried out explicitly in § 4, where corrections to the guiding-centre equations of motion arising from the radiation losses are derived consistently with the first-order guiding-centre theory. Applications of the guiding-centre RR-force are discussed in the conclusion.
2. Radiation–reaction force in a non-uniform magnetic field
The radiation–reaction force was first described for a classical non-relativistic point charge by Lorenz (see Lorentz Reference Lorentz1936). Later, Abraham (Reference Abraham1905) and Dirac (Reference Dirac1938) generalized it to relativistic energies, obtaining the Lorentz–Abraham–Dirac (LAD) force (see e.g. Pauli Reference Pauli1958)
$$\begin{eqnarray}\boldsymbol{K}=\frac{e^{2}{\it\gamma}^{2}}{6{\rm\pi}{\it\varepsilon}_{0}c^{3}}\left[\ddot{\boldsymbol{v}}+\frac{3{\it\gamma}^{2}}{c^{2}}\left(\boldsymbol{v}\boldsymbol{\cdot }\dot{\boldsymbol{v}}\right)\dot{\boldsymbol{v}}+\frac{{\it\gamma}^{2}}{c^{2}}\left(\boldsymbol{v}\boldsymbol{\cdot }\ddot{\boldsymbol{v}}+\frac{3{\it\gamma}^{2}}{c^{2}}\left(\boldsymbol{v}\boldsymbol{\cdot }\dot{\boldsymbol{v}}\right)^{2}\right)\!\boldsymbol{v}\right],\end{eqnarray}$$
where
$e$
is the particle charge,
${\it\gamma}=1/\sqrt{1-v^{2}/c^{2}}=\sqrt{1+p^{2}/(mc)^{2}}$
is the relativistic factor,
$\boldsymbol{p}={\it\gamma}m\boldsymbol{v}$
is the particle momentum,
$\boldsymbol{v}$
and
$m$
are the particles velocity and mass respectively, and
$c$
is the speed of light.
The LAD-force, however, contains third-order time derivatives with respect to particle position and cannot be uniquely solved given initial values for the particle position and velocity, therefore violating causality. Another well-known problem with the LAD-force is that in the absence of external forces it allows the existence of so-called runaway solutions (see e.g. Rohrlich Reference Rohrlich2007) that lead to exponential growth of the particle velocity. These issues have generated discussion regarding what expression to use for the RR-force (see e.g. Griffiths, Proctor & Schroeter Reference Griffiths, Proctor and Schroeter2010). Landau & Lifshitz (Reference Landau and Lifshitz1975), for example, suggest a perturbative approach where the velocity derivatives in (2.1) are expressed in terms of the external force. Ford & O’Connell (Reference Ford and O’Connell1993) note that this approach is in fact the correct one. In Spohn (Reference Spohn2000), it is shown that the LAD-force should be limited to a so-called critical surface to avoid unphysical solutions, and that this actually corresponds to the perturbative approach. In this paper, we have chosen to use the formula from Landau & Lifshitz (Reference Landau and Lifshitz1975), as it does not suffer from the unphysical behaviour associated with the LAD-force.
In magnetized plasmas, the particle motion is typically dominated by the magnetic force
$\boldsymbol{F}_{m}=e\,\boldsymbol{v}\times \boldsymbol{B}$
. If we neglect the electric field, and replace the velocity derivatives in (2.1) with the Lorentz force, the expression for the RR-force simplifies to (Landau & Lifshitz Reference Landau and Lifshitz1975)
$$\begin{eqnarray}\boldsymbol{K}=-{\it\nu}_{r}\left(\boldsymbol{p}_{\bot }+\frac{p_{\bot }^{2}}{(mc)^{2}}\boldsymbol{p}\right)-{\it\epsilon}\,{\it\nu}_{r}{\it\Omega}^{-1}B^{-1}\dot{\boldsymbol{B}}\times \boldsymbol{p},\end{eqnarray}$$
where the relativistic Larmor frequency is
${\it\Omega}=eB/({\it\gamma}m)$
. The perpendicular momentum in (2.2) is
$\boldsymbol{p}_{\bot }=(\unicode[STIX]{x1D644}-\hat{\boldsymbol{b}}\hat{\boldsymbol{b}})\boldsymbol{\cdot }\boldsymbol{p}$
, with the magnetic field unit vector being
$\hat{\boldsymbol{b}}=\boldsymbol{B}/B$
. The characteristic time for the RR-force is
$$\begin{eqnarray}{\it\nu}_{r}^{-1}=\frac{6{\rm\pi}{\it\varepsilon}_{0}{\it\gamma}(mc)^{3}}{e^{4}B^{2}}=\frac{3c}{2{\it\gamma}r_{e}{\it\Omega}^{2}},\end{eqnarray}$$
where
$r_{e}=e^{2}/(4{\rm\pi}{\it\varepsilon}_{0}mc^{2})$
is the classical electron radius. With this notation
${\it\nu}_{r}/{\it\Omega}=(2r_{e}{\it\gamma}{\it\Omega})/(3c)$
, which is typically much smaller than one (for electrons
${\it\nu}_{r}/{\it\Omega}\sim 10^{-12}B/[T]$
). It is thus obvious that the magnetic force
$|\boldsymbol{F}_{m}|\sim {\it\Omega}p$
dominates the RR-force
$|\boldsymbol{K}|\sim {\it\nu}_{r}p$
, and that the RR-force therefore acts as a dissipation to the Hamiltonian motion in the external magnetic field.
In (2.2), we have also introduced the dimensionless parameter
${\it\epsilon}$
as the dimensionless guiding-centre ordering parameter (renormalized charge
$e\rightarrow {\it\epsilon}^{-1}e$
) which will be used throughout the guiding-centre transformation. Physical results are obtained by setting
${\it\epsilon}=1$
.
3. Dissipative forces in the guiding-centre formalism
To transform general dissipative forces into the guiding-centre formalism, we cannot apply the Lie-transformation exactly as is done for deriving the Hamiltonian guiding-centre equations of motion: general dissipative forces do not necessarily have a phase-space Lagrangian formulation. Instead, we proceed via detour. Transforming the particle phase-space continuity equation we will identify the components of the guiding-centre force for any desired combination of phase-space coordinates and simultaneously guarantee the density conservation. This approach can be further used to include the RR-force into the guiding-centre kinetic equation.
Starting from the particle phase-space continuity equation
$$\begin{eqnarray}\frac{\partial f}{\partial t}+\frac{\partial }{\partial \boldsymbol{z}}\boldsymbol{\cdot }(\dot{\boldsymbol{z}}\,f)+\frac{\partial }{\partial \boldsymbol{p}}\boldsymbol{\cdot }(\boldsymbol{K}f)=0,\end{eqnarray}$$
where
$\boldsymbol{z}=(\boldsymbol{x},\boldsymbol{p})$
are the phase-space coordinates and
$\boldsymbol{K}$
now denotes a general dissipative force, we first express the continuity equation in terms of the charged-particle Hamiltonian
$$\begin{eqnarray}H={\it\gamma}mc^{2},\end{eqnarray}$$
and the non-canonical charged-particle Poisson bracket
$$\begin{eqnarray}\{f,g\}=\frac{\partial f}{\partial \boldsymbol{x}}\boldsymbol{\cdot }\frac{\partial g}{\partial \boldsymbol{p}}-\frac{\partial f}{\partial \boldsymbol{p}}\boldsymbol{\cdot }\frac{\partial g}{\partial \boldsymbol{x}}+e\boldsymbol{B}\boldsymbol{\cdot }\frac{\partial f}{\partial \boldsymbol{p}}\times \frac{\partial g}{\partial \boldsymbol{p}},\end{eqnarray}$$
to obtain the Poisson-bracket formulation of the particle phase-space continuity equation
$$\begin{eqnarray}\frac{\partial f}{\partial t}+\{f,H\}+\{x^{i},K_{i}\,f\}=0.\end{eqnarray}$$
Here
$x^{i}$
is the Cartesian component of the position
$\boldsymbol{x}$
, and summation over repeated indices is assumed. We have also made use of the identity
$$\begin{eqnarray}{\dot{z}}^{{\it\alpha}}\equiv \{z^{{\it\alpha}},H\},\end{eqnarray}$$
that describes the Hamiltonian contribution to the equations of motion, and of the properties of the Poisson bracket
$$\begin{eqnarray}\{f,g\}\equiv \frac{\partial f}{\partial z^{{\it\alpha}}}\{z^{{\it\alpha}},z^{{\it\beta}}\}\frac{\partial g}{\partial z^{{\it\beta}}}\equiv \frac{1}{{\mathcal{J}}}\frac{\partial }{\partial z^{{\it\alpha}}}({\mathcal{J}}\{f,z^{{\it\alpha}}\}\,g),\end{eqnarray}$$
where
${\mathcal{J}}$
is the phase-space Jacobian satisfying
$\text{d}\boldsymbol{z}={\mathcal{J}}\,\text{d}^{6}z$
.
The guiding-centre transformation of the particle phase-space continuity equation is then given by applying the guiding-centre push-forward
${\mathcal{T}}_{gc}^{-1}$
according to
$$\begin{eqnarray}\displaystyle & & \displaystyle {\mathcal{T}}_{gc}^{-1}\left(\frac{\partial f}{\partial t}+\{f,H\}+\{x^{i},K_{i}\,f\}\right)=0\nonumber\\ \displaystyle & & \displaystyle \quad \Rightarrow \quad \frac{\partial F}{\partial t}+\{F,H_{gc}\}_{gc}+\{{\mathcal{T}}_{gc}^{-1}x^{i},({\mathcal{T}}_{gc}^{-1}K_{i})\,F\}_{gc}=0,\end{eqnarray}$$
where
$F\equiv {\mathcal{T}}_{gc}^{-1}f$
and
$\{\cdot ,\cdot \}_{gc}\equiv {\mathcal{T}}_{gc}^{-1}\{{\mathcal{T}}_{gc}\,\cdot ,{\mathcal{T}}_{gc}\,\cdot \}$
are now the guiding-centre distribution function and Poisson bracket, respectively. For an explanation of the transformation rules, see Brizard (Reference Brizard2004). Expressing the guiding-centre Poisson bracket in a divergence form we obtain
$$\begin{eqnarray}\frac{\partial F}{\partial t}+\frac{1}{{\mathcal{J}}_{gc}}\frac{\partial }{\partial Z^{{\it\alpha}}}[{\mathcal{J}}_{gc}({\dot{Z}}^{{\it\alpha}}+\{{\mathcal{T}}_{gc}^{-1}x^{i},Z^{{\it\alpha}}\}_{gc}({\mathcal{T}}_{gc}^{-1}K_{i}))F]=0,\end{eqnarray}$$
where
$Z^{{\it\alpha}}$
are the guiding-centre phase-space coordinates and
${\mathcal{J}}_{gc}$
is the guiding-centre phase-space Jacobian,
$\text{d}\boldsymbol{Z}={\mathcal{J}}_{gc}\,\text{d}^{6}Z$
.
If the characteristic time scale,
${\it\nu}_{k}^{-1}$
, of the dissipative force,
$|\boldsymbol{K}|\sim {\it\nu}_{k}p$
, is much longer than the time scale related to the gyromotion, i.e.
${\it\nu}_{k}/{\it\Omega}\ll 1$
, a closure scheme to obtain an equation for a gyro-averaged distribution function,
$\langle F\rangle$
, proceeds in a similar manner as was presented for the collision operator in Brizard (Reference Brizard2004). To lowest order in
${\it\epsilon}_{{\it\nu}}={\it\nu}_{k}/{\it\Omega}$
, the gyro-averaged guiding-centre continuity equation thus becomes
$$\begin{eqnarray}\frac{\partial \langle F\rangle }{\partial t}+\frac{1}{{\mathcal{J}}_{gc}}\frac{\partial }{\partial Z^{{\it\alpha}}}[{\mathcal{J}}_{gc}({\dot{Z}}^{{\it\alpha}}+\langle \{{\mathcal{T}}_{gc}^{-1}x^{i},Z^{{\it\alpha}}\}_{gc}({\mathcal{T}}_{gc}^{-1}K_{i})\rangle )\langle F\rangle ]=0,\end{eqnarray}$$
where we have used the fact that the Hamiltonian guiding-centre equations of motion
${\dot{Z}}^{{\it\alpha}}$
are, by construction, independent of the gyro-angle.
From the gyro-averaged guiding-centre continuity equation, we finally obtain an expression for a dissipative gyro-averaged guiding-centre force
$$\begin{eqnarray}{\mathcal{K}}_{gc}^{{\it\alpha}}=\langle \{{\mathcal{T}}_{gc}^{-1}x^{i},Z^{{\it\alpha}}\}_{gc}({\mathcal{T}}_{gc}^{-1}K_{i})\rangle \equiv \langle {\it\Delta}^{i{\it\alpha}}\cdot {\mathcal{T}}_{gc}^{-1}K_{i}\rangle ,\end{eqnarray}$$
where
$\langle \cdots \rangle =1/(2{\rm\pi})\int _{0}^{2{\rm\pi}}\,\text{d}{\it\theta}\,\dots$
is the gyro-angle average, and
$$\begin{eqnarray}{\it\bf\Delta}^{{\it\alpha}}\equiv \hat{\boldsymbol{e}}_{i}{\it\Delta}^{i{\it\alpha}}\equiv \hat{\boldsymbol{e}}_{i}\{{\mathcal{T}}_{gc}^{-1}x^{i},Z^{{\it\alpha}}\}_{gc}\end{eqnarray}$$
is the so-called guiding-centre projection coefficient.
4. First-order transformation in
$(\boldsymbol{X}\!,p_{\Vert },{\it\mu})$
phase-space
As the time scale related to the radiation reaction satisfies the condition
${\it\nu}_{r}/{\it\Omega}\ll 1$
, we can calculate the guiding-centre transformation of the RR-force according to the method described in § 3. In order to proceed, explicit expressions for the projection coefficients
${\it\Delta}^{i{\it\alpha}}$
are derived in appendix A.
Considering the phase-space
$(\boldsymbol{X},p_{\Vert },{\it\mu})$
, often used in particle tracing, we obtain
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\Delta}^{ij}=-{\it\epsilon}\frac{\hat{\boldsymbol{b}}}{eB_{\Vert }^{\star }}\times (\hat{\boldsymbol{e}}^{i}+\boldsymbol{{\rm\nabla}}^{\star }{\it\rho}_{{\it\epsilon}}^{i})\boldsymbol{\cdot }\hat{\boldsymbol{e}}^{\,j}-\frac{B^{\star j}}{B_{\Vert }^{\star }}\frac{\partial {\it\rho}_{{\it\epsilon}}^{i}}{\partial p_{\Vert }}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\Delta}^{ip_{\Vert }}=\frac{\boldsymbol{B}^{\star }}{B_{\Vert }^{\star }}\boldsymbol{\cdot }(\hat{\boldsymbol{e}}^{i}+\boldsymbol{{\rm\nabla}}^{\star }{\it\rho}_{{\it\epsilon}}^{i}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\Delta}^{i{\it\mu}}={\it\epsilon}^{-1}\,\frac{e}{m}\frac{\partial {\it\rho}_{{\it\epsilon}}^{i}}{\partial {\it\theta}}, & \displaystyle\end{eqnarray}$$
${\it\epsilon}$
is the dimensionless guiding-centre ordering parameter and the transformed gyroradius vector is defined as 
$$\begin{eqnarray}{\bf\rho}_{{\it\epsilon}}\equiv {\mathcal{T}}_{gc}^{-1}\boldsymbol{x}-\boldsymbol{X}\equiv {\it\epsilon}{\bf\rho}_{0}+{\it\epsilon}^{2}{\bf\rho}_{1}+\cdots \,,\end{eqnarray}$$
with the sub-indices referring to the order with respect to magnetic field non-uniformity. Thus, keeping the gyroradius up to the term
${\bf\rho}_{1}$
gives projection coefficients that are valid up to first order in the magnetic field non-uniformity.
The expressions for the so-called symplectic or effective magnetic field
$\boldsymbol{B}^{\star }$
and the modified gradient operator
$\boldsymbol{{\rm\nabla}}^{\star }$
as well as for the zeroth- and first-order gyroradius vectors,
${\bf\rho}_{0}$
and
${\bf\rho}_{1}$
, are given in appendix A.
Similarly as for
${\bf\rho}_{{\it\epsilon}}$
, we have for the push-forward of the particle phase-space RR-force
$$\begin{eqnarray}{\mathcal{T}}_{gc}^{-1}\boldsymbol{K}\equiv \boldsymbol{K}_{{\it\epsilon}}=\boldsymbol{K}_{0}+{\it\epsilon}\boldsymbol{K}_{1}+\cdots \,,\end{eqnarray}$$
where the sub-indices again refer to the order with respect to the magnetic field non-uniformity. An explicit expression for
$\boldsymbol{K}_{{\it\epsilon}}$
is given in appendix B. Using the expressions for the projection coefficients we then find the components for the guiding-centre radiation–reaction force
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{K}}^{\boldsymbol{X}}=-{\it\epsilon}\frac{\hat{\boldsymbol{b}}}{eB_{\Vert }^{\star }}\times \langle \boldsymbol{K}_{{\it\epsilon}}+\boldsymbol{{\rm\nabla}}^{\star }{\bf\rho}_{{\it\epsilon}}\boldsymbol{\cdot }\boldsymbol{K}_{{\it\epsilon}}\rangle -\frac{\boldsymbol{B}^{\star }}{B_{\Vert }^{\star }}\left\langle \frac{\partial {\bf\rho}_{{\it\epsilon}}}{\partial p_{\Vert }}\boldsymbol{\cdot }\boldsymbol{K}_{{\it\epsilon}}\right\rangle , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{K}}^{p_{\Vert }}=\frac{\boldsymbol{B}^{\star }}{B_{\Vert }^{\star }}\boldsymbol{\cdot }\langle \boldsymbol{K}_{{\it\epsilon}}+\boldsymbol{{\rm\nabla}}^{\star }{\bf\rho}_{{\it\epsilon}}\boldsymbol{\cdot }\boldsymbol{K}_{{\it\epsilon}}\rangle , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{K}}^{{\it\mu}}={\it\epsilon}^{-1}\,\frac{e}{m}\left\langle \frac{\partial {\bf\rho}_{{\it\epsilon}}}{\partial {\it\theta}}\boldsymbol{\cdot }\boldsymbol{K}_{{\it\epsilon}}\right\rangle . & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{K}}^{\boldsymbol{X}}=-{\it\epsilon}\frac{\hat{\boldsymbol{b}}}{eB_{\Vert }^{\star }}\times \langle \boldsymbol{K}_{0}+{\it\epsilon}\boldsymbol{K}_{1}+{\it\epsilon}\boldsymbol{{\rm\nabla}}^{\star }{\bf\rho}_{0}\boldsymbol{\cdot }\boldsymbol{K}_{0}\rangle -{\it\epsilon}^{2}\hat{\boldsymbol{b}}\left\langle \frac{\partial {\bf\rho}_{1}}{\partial p_{\Vert }}\boldsymbol{\cdot }\boldsymbol{K}_{0}\right\rangle , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{K}}^{p_{\Vert }}=\frac{\boldsymbol{B}^{\star }}{B_{\Vert }^{\star }}\boldsymbol{\cdot }\langle \boldsymbol{K}_{0}\rangle +{\it\epsilon}\hat{\boldsymbol{b}}\boldsymbol{\cdot }\langle \boldsymbol{K}_{1}+\boldsymbol{{\rm\nabla}}^{\star }{\bf\rho}_{0}\boldsymbol{\cdot }\boldsymbol{K}_{0}\rangle , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{K}}^{{\it\mu}}=\frac{e}{m}\left\langle \frac{\partial {\bf\rho}_{0}}{\partial {\it\theta}}\boldsymbol{\cdot }(\boldsymbol{K}_{0}+{\it\epsilon}\boldsymbol{K}_{1})\right\rangle +{\it\epsilon}\frac{e}{m}\left\langle \frac{\partial {\bf\rho}_{1}}{\partial {\it\theta}}\boldsymbol{\cdot }\boldsymbol{K}_{0}\right\rangle , & \displaystyle\end{eqnarray}$$
${\bf\rho}_{0}$
is independent of
$p_{\Vert }$
and that
$\boldsymbol{{\rm\nabla}}^{\star }{\bf\rho}_{1}$
is of second order in magnetic field non-uniformity.
Gyro-averages of the expressions in (4.9)–(4.11) are carried out in appendix C. The resulting components of the guiding-centre RR-force, acting as dissipative terms in the equations of motion for the corresponding coordinate, are for the guiding-centre position:
$$\begin{eqnarray}{\mathcal{K}}^{\boldsymbol{X}}=-{\it\epsilon}_{B}\,\frac{{\it\nu}_{r}}{{\it\Omega}_{\Vert }^{\star }}\frac{2{\it\mu}B}{mc^{2}}(\hat{\boldsymbol{b}}\times \dot{\boldsymbol{X}}+3v_{\Vert }\,{\it\varrho}_{\Vert }\,{\bf\kappa}),\end{eqnarray}$$
for the parallel momentum:
$$\begin{eqnarray}{\mathcal{K}}^{p_{\Vert }}=-{\it\nu}_{r}\,p_{\Vert }\frac{{\it\mu}B}{mc^{2}}(2+{\it\epsilon}_{B}{\it\varrho}_{\Vert }{\it\tau}_{B})-{\it\epsilon}_{B}\,{\it\nu}_{r}\frac{m{\it\mu}{\it\gamma}^{2}}{e}{\it\tau}_{B},\end{eqnarray}$$
and for the magnetic moment:
$$\begin{eqnarray}{\mathcal{K}}^{{\it\mu}}=-{\it\nu}_{r}\,{\it\mu}\left(1+\frac{2{\it\mu}B}{mc^{2}}\right)(2+{\it\epsilon}_{B}{\it\varrho}_{\Vert }{\it\tau}_{B}),\end{eqnarray}$$
where the parameter
${\it\epsilon}_{B}$
is introduced to explicitly point out the contribution from magnetic field non-uniformity. The parallel gyroradius and the modified gyro-frequency are defined as
$$\begin{eqnarray}{\it\varrho}_{\Vert }=\frac{p_{\Vert }}{eB},\quad {\it\Omega}_{\Vert }^{\star }={\it\Omega}(1+{\it\epsilon}_{B}{\it\varrho}_{\Vert }{\it\tau}_{B}).\end{eqnarray}$$
The first-order corrections in (4.13) and (4.14) relate to the magnetic field-line twist parameter,
${\it\tau}_{B}=\hat{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\times \hat{\boldsymbol{b}}$
, which also appears in Hamiltonian guiding-centre motion as the phase-space Jacobian is
${\mathcal{J}}_{gc}\equiv B_{\Vert }^{\star }\equiv B\,(1+{\it\epsilon}_{B}{\it\varrho}_{\Vert }{\it\tau}_{B})$
.
While the phase-space
$(\boldsymbol{X},p_{\Vert },{\it\mu})$
was used to carry out the explicit guiding-centre transformation, properties of the Poisson bracket (3.3) provide general rules for transforming between different phase-spaces according to
$$\begin{eqnarray}{\mathcal{K}}_{gc}^{{\it\alpha}}={\mathcal{K}}^{\boldsymbol{X}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\alpha}+{\mathcal{K}}^{p_{\Vert }}\frac{\partial {\it\alpha}}{\partial p_{\Vert }}+{\mathcal{K}}^{{\it\mu}}\frac{\partial {\it\alpha}}{\partial {\it\mu}}.\end{eqnarray}$$
For example, using a set of common phase-space coordinates
$(\boldsymbol{X},p,{\it\xi})$
, where
$p$
is the total guiding-centre momentum and
${\it\xi}=p_{\Vert }/p$
, one can calculate the partial derivatives for the total momentum:
$$\begin{eqnarray}\frac{\partial p}{\partial p_{\Vert }}={\it\xi},\quad \frac{\partial p}{\partial {\it\mu}}=\frac{mB}{p},\quad \boldsymbol{{\rm\nabla}}p=0,\end{eqnarray}$$
and for the pitch-angle cosine:
$$\begin{eqnarray}\frac{\partial {\it\xi}}{\partial p_{\Vert }}=\frac{1-{\it\xi}^{2}}{p},\quad \frac{\partial {\it\xi}}{\partial {\it\mu}}=-\frac{mB}{p^{2}}{\it\xi},\quad \boldsymbol{{\rm\nabla}}{\it\xi}=-\frac{1-{\it\xi}^{2}}{2{\it\xi}}\boldsymbol{{\rm\nabla}}\ln B,\end{eqnarray}$$
and evaluate the guiding-centre radiation–reaction coefficients for the momentum component:
$$\begin{eqnarray}{\mathcal{K}}^{p}=-{\it\nu}_{r}{\it\gamma}^{2}p(1-{\it\xi}^{2})(1+{\it\epsilon}_{B}{\it\varrho}_{\Vert }{\it\tau}_{B}),\end{eqnarray}$$
as well for the pitch-angle component:
$$\begin{eqnarray}{\mathcal{K}}^{{\it\xi}}={\it\nu}_{r}\,{\it\xi}(1-{\it\xi}^{2})(1+{\it\epsilon}_{B}{\it\varrho}_{\Vert }{\it\tau}_{B}/2)-{\it\epsilon}_{B}{\it\nu}_{r}{\it\gamma}^{2}(1-{\it\xi}^{2})^{2}{\it\varrho}{\it\tau}_{B}/2.\end{eqnarray}$$
Here we have introduced
${\it\varrho}=p/(eB)$
. In the limit of uniform magnetic field, these coincide with the ones used e.g. in Stahl et al. (Reference Stahl, Hirvijoki, Decker, Embréus and Fülöp2015).
5. Conclusions
For the first time, a consistent guiding-centre transformation of the radiation–reaction force for a particle travelling in an external magnetic field is presented up to the first order in the magnetic-field non-uniformity. As a result, we observe corrections that are proportional to the magnetic field-line twist, which itself plays an important role in the Hamiltonian equations of motion. As the magnetic moment is an exact invariant of Hamiltonian guiding-centre motion, the first-order correction, especially to the dissipative evolution of the magnetic moment, could be important: numerical simulations of tokamak first wall power loads from fast particles are very sensitive to the details of the magnetic field non-uniformities and to the guiding-centre phase-space trajectories. Due to the presence of the three periodic motions (gyro, bounce, and precession), small deviations in the guiding-centre trajectory may cumulate and change the wall power loads. As magnetic perturbations are considered as an option to mitigate the formation of dangerously large runaway-electron beams in tokamaks, equations to model the dissipative guiding-centre motion must be accurate and consistent with the rest of the tools.
The equations derived in this paper are applicable also to particle dynamics in astrophysical plasmas as well as to any magnetically confined laboratory plasmas as long as the guiding-centre formalism itself is valid. The paper also provides a method for transforming general dissipative forces to guiding-centre phase-space. The procedure is valid as long as the relative momentum loss over a gyroperiod is sufficiently small.
Carrying the transformation up to second order in magnetic field non-uniformity would provide an additional term in
${\mathcal{K}}^{p_{\Vert }}$
that does not vanish for
${\it\mu}\rightarrow 0$
. This component would further contribute to the guiding-centre motion along field lines in a curved magnetic field, and would dominate in the limit
${\it\mu}\rightarrow 0$
. The corresponding second-order guiding-centre transformation of the radiation–reaction force will be presented in a future contribution.
As the momentum loss due to the RR-force is also typically small over a particle bounce or transit time, performing an orbit averaging operation for axisymmetric configurations as prescribed in Brizard et al. (Reference Brizard, Decker, Peysson and Duthoit2009) would yield a reduced orbit-averaged guiding-centre RR-force operator in a three-dimensional phase-space. Radial transport coefficients including neoclassical effects could then be explicitly derived. Such an operator could then be readily implemented in a 3-D orbit-averaged guiding-centre kinetic code (Decker et al. Reference Decker, Peysson, Brizard and Duthoit2010).
Acknowledgements
We would like to thank Dr I. Pusztai, Mr A. Stahl, and Professor T. Fülöp for fruitful discussions on improving the manuscript. Work by A.J.B. was supported by a US DoE grant under contract no. DE-SC0006721.
Appendix A. Relativistic guiding-centre transformation
The relativistic guiding-centre Lagrangian one-form for the guiding-centre phase-space coordinates
$Z^{{\it\alpha}}=(\boldsymbol{X},p_{\Vert },{\it\mu},{\it\theta})$
is
$$\begin{eqnarray}{\it\Gamma}_{gc}\equiv {\it\Gamma}_{{\it\alpha}}\,\text{d}Z^{{\it\alpha}}-H_{gc}\,\text{d}t=({\it\epsilon}^{-1}e\boldsymbol{A}+p_{\Vert }\hat{\boldsymbol{b}})\boldsymbol{\cdot }\text{d}\boldsymbol{X}+{\it\epsilon}\frac{m{\it\mu}}{e}(\text{d}{\it\theta}-\boldsymbol{R}^{\star }\boldsymbol{\cdot }\text{d}\boldsymbol{X})-{\it\gamma}mc^{2}\,\text{d}t,\end{eqnarray}$$
where
${\it\gamma}=\sqrt{1+p_{\Vert }^{2}/(mc)^{2}+2{\it\mu}B/(mc^{2})}$
and
${\it\epsilon}$
is the guiding-centre ordering parameter, the modified gyrogauge field is
$\boldsymbol{R}^{\star }=\boldsymbol{R}+({\it\tau}_{B}/2)\hat{\boldsymbol{b}}$
with
$\boldsymbol{R}=\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{\bot }}\boldsymbol{\cdot }\hat{{\bf\rho}}\equiv \boldsymbol{{\rm\nabla}}\hat{\mathbf{1}}\boldsymbol{\cdot }\hat{\mathbf{2}}$
the Littlejohn’s gyrogauge vector, and
${\it\tau}_{B}=\hat{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\times \hat{\boldsymbol{b}}$
the magnetic field line torsion. The two right-handed orthogonal unit vector sets,
$(\hat{\boldsymbol{b}}(\boldsymbol{X}),\hat{\boldsymbol{\bot }}(\boldsymbol{X},{\it\theta}),\hat{{\bf\rho}}(\boldsymbol{X},{\it\theta}))$
and
$(\hat{\boldsymbol{b}}(\boldsymbol{X}),\hat{\mathbf{1}}(\boldsymbol{X}),\hat{\mathbf{2}}(\boldsymbol{X}))$
are
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{{\bf\rho}}=\cos {\it\theta}\,\hat{\mathbf{1}}-\sin {\it\theta}\,\hat{\mathbf{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{\boldsymbol{\bot }}=-\text{sin}\,{\it\theta}\,\hat{\mathbf{1}}-\cos {\it\theta}\,\hat{\mathbf{2}}. & \displaystyle\end{eqnarray}$$
${\it\Gamma}_{{\it\alpha}}\,\text{d}Z^{{\it\alpha}}$
is 
$$\begin{eqnarray}\displaystyle \{F,G\}_{gc} & = & \displaystyle {\it\epsilon}^{-1}\frac{e}{m}\left(\frac{\partial F}{\partial {\it\theta}}\frac{\partial G}{\partial {\it\mu}}-\frac{\partial F}{\partial {\it\mu}}\frac{\partial G}{\partial {\it\theta}}\right)\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\boldsymbol{B}^{\star }}{B_{\Vert }^{\star }}\boldsymbol{\cdot }\left(\boldsymbol{{\rm\nabla}}^{\star }F\frac{\partial G}{\partial p_{\Vert }}-\frac{\partial F}{\partial p_{\Vert }}\boldsymbol{{\rm\nabla}}^{\star }G\right)-{\it\epsilon}\,\frac{\hat{\boldsymbol{b}}}{eB_{\Vert }^{\star }}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}^{\star }F\times \boldsymbol{{\rm\nabla}}^{\star }G,\end{eqnarray}$$
where the modified gradient operator is
$\boldsymbol{{\rm\nabla}}^{\star }=\boldsymbol{{\rm\nabla}}+\boldsymbol{R}^{\star }\partial /\partial {\it\theta}$
, the phase-space Jacobian is
${\mathcal{J}}_{gc}=\boldsymbol{B}^{\star }\boldsymbol{\cdot }\hat{\boldsymbol{b}}\equiv B_{\Vert }^{\star }$
, and the effective magnetic field is
$$\begin{eqnarray}\boldsymbol{B}^{\star }=\boldsymbol{{\rm\nabla}}\times \left(\boldsymbol{A}+{\it\epsilon}\frac{p_{\Vert }}{e}\hat{\boldsymbol{b}}\right)+O({\it\epsilon}^{2})=\boldsymbol{B}+{\it\epsilon}\frac{p_{\Vert }}{e}\boldsymbol{{\rm\nabla}}\times \hat{\boldsymbol{b}}+O({\it\epsilon}^{2}).\end{eqnarray}$$
We will also find useful the expression
$$\begin{eqnarray}\frac{\boldsymbol{B}^{\star }}{B_{\Vert }^{\star }}\equiv \hat{\boldsymbol{b}}+{\it\epsilon}\frac{p_{\Vert }}{eB_{\Vert }^{\star }}\hat{\boldsymbol{b}}\times {\bf\kappa}+O({\it\epsilon}^{2}).\end{eqnarray}$$
The generating functions
$G_{n}^{{\it\alpha}}$
that define the coordinate transformations between the guiding-centre coordinates
$Z^{{\it\alpha}}$
and particle coordinates
$z^{{\it\alpha}}$
according to
$$\begin{eqnarray}\displaystyle & \displaystyle Z^{{\it\alpha}}=z^{{\it\alpha}}+{\it\epsilon}\,G_{1}^{{\it\alpha}}+{\it\epsilon}^{2}\left(G_{2}^{{\it\alpha}}+\frac{1}{2}G_{1}^{{\it\beta}}\frac{\partial G_{1}^{{\it\alpha}}}{\partial z^{{\it\beta}}}\right)+O({\it\epsilon}^{2}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle z^{{\it\alpha}}=Z^{{\it\alpha}}-{\it\epsilon}\,G_{1}^{{\it\alpha}}-{\it\epsilon}^{2}\left(G_{2}^{{\it\alpha}}-\frac{1}{2}G_{1}^{{\it\beta}}\frac{\partial G_{1}^{{\it\alpha}}}{\partial Z^{{\it\beta}}}\right)+O({\it\epsilon}^{2}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle G_{1}^{\boldsymbol{X}}=-{\bf\rho}_{0}\equiv -\sqrt{\frac{2m{\it\mu}}{e^{2}B}}\hat{{\bf\rho}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle G_{1}^{p_{\Vert }}=-p_{\Vert }\,{\bf\rho}_{0}\boldsymbol{\cdot }{\bf\kappa}+\frac{m{\it\mu}}{e}({\it\tau}_{B}+\unicode[STIX]{x1D622}_{1}:\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{b}}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle G_{1}^{{\it\mu}}={\bf\rho}_{0}\boldsymbol{\cdot }\left({\it\mu}\boldsymbol{{\rm\nabla}}\ln B+\frac{p_{\Vert }^{2}}{mB}{\bf\kappa}\right)-\frac{{\it\mu}\,p_{\Vert }}{eB}({\it\tau}_{B}+\unicode[STIX]{x1D622}_{1}:\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{b}}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle G_{1}^{{\it\theta}}=-{\bf\rho}_{0}\boldsymbol{\cdot }\boldsymbol{R}+\frac{p_{\Vert }}{eB}\left(\unicode[STIX]{x1D622}_{2}:\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{b}}\right)+\frac{\partial {\bf\rho}_{0}}{\partial {\it\theta}}\boldsymbol{\cdot }\left(\boldsymbol{{\rm\nabla}}\ln B+\frac{p_{\Vert }^{2}}{2m\,{\it\mu}B}{\bf\kappa}\right), & \displaystyle\end{eqnarray}$$
${\bf\kappa}=\hat{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{b}}$
. We also need the spatial component of the second-order generating function 
$$\begin{eqnarray}\displaystyle G_{2}^{\boldsymbol{X}} & = & \displaystyle \left[\frac{2p_{\Vert }}{eB}\left(\frac{\partial {\bf\rho}_{0}}{\partial {\it\theta}}\boldsymbol{\cdot }{\bf\kappa}\right)+\frac{m{\it\mu}}{e^{2}B}(\unicode[STIX]{x1D622}_{2}:\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{b}})\right]\hat{\boldsymbol{b}}+\frac{p_{\Vert }}{eB}\,{\it\tau}_{B}\,{\bf\rho}_{0}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{2}(G_{1}^{{\it\mu}}-{\it\mu}{\bf\rho}_{0}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\ln B)\frac{\partial {\bf\rho}_{0}}{\partial {\it\mu}}+\frac{1}{2}(G_{1}^{{\it\theta}}+{\bf\rho}_{0}\boldsymbol{\cdot }\boldsymbol{R})\frac{\partial {\bf\rho}_{0}}{\partial {\it\theta}}.\end{eqnarray}$$
The dyads
$\unicode[STIX]{x1D622}_{1}$
and
$\unicode[STIX]{x1D622}_{2}$
are
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D622}_{1}\equiv -\frac{1}{2}(\hat{{\bf\rho}}\hat{\boldsymbol{\bot }}+\hat{\boldsymbol{\bot }}\hat{{\bf\rho}})=\frac{\partial \unicode[STIX]{x1D622}_{2}}{\partial {\it\theta}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D622}_{2}\equiv \frac{1}{4}(\hat{\boldsymbol{\bot }}\hat{\boldsymbol{\bot }}-\hat{{\bf\rho}}\hat{{\bf\rho}})=-\frac{1}{4}\frac{\partial \unicode[STIX]{x1D622}_{1}}{\partial {\it\theta}}. & \displaystyle\end{eqnarray}$$
With the guiding-centre Poisson bracket and Hamiltonian given, obtaining the Hamiltonian equations of motion for each phase-space coordinate is then straightforward. For the phase-space
$Z^{{\it\alpha}}=(\boldsymbol{X},p_{\Vert },{\it\mu},{\it\theta})$
we find
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\boldsymbol{X}}=\{\boldsymbol{X},H_{gc}\}=\frac{p_{\Vert }}{{\it\gamma}m}\frac{\boldsymbol{B}^{\star }}{B_{\Vert }^{\star }}+{\it\epsilon}\frac{\hat{\boldsymbol{b}}}{eB_{\Vert }^{\star }}\times \frac{{\it\mu}}{{\it\gamma}}\boldsymbol{{\rm\nabla}}B, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{p}}_{\Vert }=\{p_{\Vert },H_{gc}\}=-\frac{\boldsymbol{B}^{\star }}{B_{\Vert }^{\star }}\boldsymbol{\cdot }\frac{{\it\mu}}{{\it\gamma}}\boldsymbol{{\rm\nabla}}B, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{{\it\mu}}=\{{\it\mu},H_{gc}\}={\it\epsilon}^{-1}\frac{e}{m}\frac{\partial H_{gc}}{\partial {\it\theta}}\equiv 0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{{\it\theta}}={\it\epsilon}^{-1}\frac{eB}{{\it\gamma}m}+\dot{\boldsymbol{X}}\boldsymbol{\cdot }\boldsymbol{R}^{\star }\equiv {\it\epsilon}^{-1}{\it\Omega}+\dot{\boldsymbol{X}}\boldsymbol{\cdot }\boldsymbol{R}^{\star }. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\boldsymbol{{\rm\nabla}}^{\star }{\bf\rho}_{0}=-\frac{1}{2}\boldsymbol{{\rm\nabla}}\ln B\,{\bf\rho}_{0}+\frac{1}{2}{\it\tau}_{B}\hat{\boldsymbol{b}}\frac{\partial {\bf\rho}_{0}}{\partial {\it\theta}}-(\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{b}}\boldsymbol{\cdot }{\bf\rho}_{0})\hat{\boldsymbol{b}},\end{eqnarray}$$
where
$\boldsymbol{{\rm\nabla}}\times \hat{\boldsymbol{b}}={\it\tau}_{B}\hat{\boldsymbol{b}}+\hat{\boldsymbol{b}}\times {\bf\kappa}$
, and also
$$\begin{eqnarray}\hat{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}^{\star }{\bf\rho}_{0}=\frac{{\it\tau}_{B}}{2}\frac{\partial {\bf\rho}_{0}}{\partial {\it\theta}}-{\it\mu}\,\boldsymbol{{\rm\nabla}}_{\Vert }\ln B\frac{\partial {\bf\rho}_{0}}{\partial {\it\mu}}-({\bf\rho}_{0}\boldsymbol{\cdot }{\bf\kappa})\hat{\boldsymbol{b}}.\end{eqnarray}$$
Appendix B. Push-forward of the radiation–reaction force
Noting that
${\mathcal{T}}_{gc}^{-1}{\it\gamma}\equiv {\it\gamma}$
, the push-forwards of particle momentum and magnetic field time derivative become
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{T}}_{gc}^{-1}\boldsymbol{p}\equiv {\mathcal{T}}_{gc}^{-1}({\it\gamma}m\boldsymbol{v})={\it\gamma}m\left({\mathcal{T}}_{gc}^{-1}\frac{\text{d}}{\text{d}t}{\mathcal{T}}_{gc}\right)({\mathcal{T}}_{gc}^{-1}\boldsymbol{x}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{T}}_{gc}^{-1}\dot{\boldsymbol{B}}\equiv \left({\mathcal{T}}_{gc}^{-1}\frac{\text{d}}{\text{d}t}{\mathcal{T}}_{gc}\right)({\mathcal{T}}_{gc}^{-1}\boldsymbol{B}). & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\left({\mathcal{T}}_{gc}^{-1}\frac{\text{d}}{\text{d}t}{\mathcal{T}}_{gc}\right)\equiv \frac{\partial }{\partial t}+\dot{\boldsymbol{X}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}+{\dot{p}}_{\Vert }\frac{\partial }{\partial p_{\Vert }}+\dot{{\it\theta}}\frac{\partial }{\partial {\it\theta}}\equiv \frac{\partial }{\partial t}+\dot{\boldsymbol{X}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}^{\star }+{\dot{p}}_{\Vert }\frac{\partial }{\partial p_{\Vert }}+{\it\epsilon}^{-1}{\it\Omega}\frac{\partial }{\partial {\it\theta}},\end{eqnarray}$$
and because it involves a term of order
${\it\epsilon}^{-1}$
, the push-forward of particle position
$$\begin{eqnarray}{\mathcal{T}}_{gc}^{-1}\boldsymbol{x}\equiv \boldsymbol{X}-{\it\epsilon}G_{1}^{\boldsymbol{X}}-{\it\epsilon}^{2}\left(G_{2}^{\boldsymbol{X}}-\frac{1}{2}G_{1}^{{\it\beta}}\frac{\partial G_{1}^{\boldsymbol{X}}}{\partial Z^{{\it\beta}}}\right)+O({\it\epsilon}^{3})\equiv \boldsymbol{X}+{\it\epsilon}{\bf\rho}_{0}+{\it\epsilon}^{2}{\bf\rho}_{1}+O({\it\epsilon}^{3}),\end{eqnarray}$$
and the push-forward of the magnetic field
$$\begin{eqnarray}{\mathcal{T}}_{gc}^{-1}\boldsymbol{B}\equiv \boldsymbol{B}-{\it\epsilon}G_{1}^{\boldsymbol{X}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{B}-{\it\epsilon}^{2}\left[G_{2}^{\boldsymbol{X}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{B}-\frac{1}{2}G_{1}^{{\it\beta}}\frac{\partial }{\partial Z^{{\it\beta}}}(G_{1}^{\boldsymbol{X}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{B})\right]+O({\it\epsilon}^{3}),\end{eqnarray}$$
have to be evaluated up to second order in
${\it\epsilon}^{2}$
. The explicit expression for the first-order Larmor radius vector is
$$\begin{eqnarray}\displaystyle {\bf\rho}_{1} & = & \displaystyle -\!\left[2\frac{p_{\Vert }}{eB}\left({\bf\kappa}\boldsymbol{\cdot }\frac{\partial {\bf\rho}_{0}}{\partial {\it\theta}}\right)-\frac{m{\it\mu}}{e^{2}B}\left(\unicode[STIX]{x1D622}_{2}:\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{b}}-\frac{1}{2}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\hat{\boldsymbol{b}}\right)\right]\hat{\boldsymbol{b}}\nonumber\\ \displaystyle & & \displaystyle -\!\left[\frac{1}{2}\frac{p_{\Vert }}{eB}({\it\tau}_{B}-\unicode[STIX]{x1D622}_{1}:\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{b}})+\frac{\partial {\bf\rho}_{0}}{\partial {\it\mu}}\boldsymbol{\cdot }\left({\it\mu}\boldsymbol{{\rm\nabla}}\ln B+\frac{p_{\Vert }^{2}}{mB}{\bf\kappa}\right)\right]{\bf\rho}_{0}\nonumber\\ \displaystyle & & \displaystyle -\!\left[\frac{p_{\Vert }}{eB}(\unicode[STIX]{x1D622}_{2}:\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{b}})+\frac{\partial {\bf\rho}_{0}}{\partial {\it\theta}}\boldsymbol{\cdot }\left(\boldsymbol{{\rm\nabla}}\ln B+\frac{p_{\Vert }^{2}}{2m{\it\mu}B}{\bf\kappa}\right)\right]\frac{\partial {\bf\rho}_{0}}{\partial {\it\theta}}.\end{eqnarray}$$
Now, up to first order in
${\it\epsilon}$
, we find the push-forward of the particle momentum
$$\begin{eqnarray}{\mathcal{T}}_{gc}^{-1}\boldsymbol{p}={\it\gamma}m\dot{\boldsymbol{X}}+eB\frac{\partial {\bf\rho}_{0}}{\partial {\it\theta}}+{\it\epsilon}\left(p_{\Vert }\hat{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}^{\star }{\bf\rho}_{0}+eB\frac{\partial {\bf\rho}_{1}}{\partial {\it\theta}}\right)+O({\it\epsilon}^{2})\equiv \boldsymbol{p}_{0}+{\it\epsilon}\,\boldsymbol{p}_{1}+O({\it\epsilon}^{2}),\end{eqnarray}$$
where the zeroth and first components are given by
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{p}_{0}=p_{\Vert }\hat{\boldsymbol{b}}+eB\frac{\partial {\bf\rho}_{0}}{\partial {\it\theta}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{p}_{1}=\frac{\hat{\boldsymbol{b}}}{eB_{\Vert }^{\star }}\times (m{\it\mu}\boldsymbol{{\rm\nabla}}B+p_{\Vert }^{2}{\bf\kappa})+p_{\Vert }\hat{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}^{\star }{\bf\rho}_{0}+eB\frac{\partial {\bf\rho}_{1}}{\partial {\it\theta}}. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}{\mathcal{T}}_{gc}^{-1}{\it\nu}_{r}={\it\nu}_{r}(1+2\,{\it\epsilon}\,{\bf\rho}_{0}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\ln B)+O({\it\epsilon}^{2}),\end{eqnarray}$$
the push-forward of the perpendicular momentum
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{T}}_{gc}^{-1}\boldsymbol{p}_{\bot }=eB\frac{\partial {\bf\rho}_{0}}{\partial {\it\theta}}+{\it\epsilon}(\boldsymbol{p}_{1}+G_{1}^{p_{\Vert }}\hat{\boldsymbol{b}}-p_{\Vert }{\bf\rho}_{0}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{b}})+O({\it\epsilon}^{2}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{T}}_{gc}^{-1}p_{\bot }^{2}=2m{\it\mu}B\left(1+{\it\epsilon}{\bf\rho}_{0}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\ln B-{\it\epsilon}\frac{G_{1}^{{\it\mu}}}{{\it\mu}}\right)+O({\it\epsilon}^{2}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\mathcal{T}}_{gc}^{-1}\dot{\boldsymbol{B}} & = & \displaystyle \frac{p_{\Vert }}{{\it\gamma}m}\hat{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{B}+{\it\Omega}\frac{\partial {\bf\rho}_{0}}{\partial {\it\theta}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{B}+{\it\epsilon}\frac{p_{\Vert }}{{\it\gamma}m}\hat{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}^{\star }({\bf\rho}_{0}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{B})\nonumber\\ \displaystyle & & \displaystyle -\,{\it\epsilon}{\it\Omega}\left[\frac{\partial G_{2}^{\boldsymbol{X}}}{\partial {\it\theta}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{B}+\frac{1}{2}\frac{\partial }{\partial {\it\theta}}\left(G_{1}^{{\it\beta}}\frac{\partial }{\partial Z^{{\it\beta}}}({\bf\rho}_{0}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{B})\right)\right]\nonumber\\ \displaystyle & & \displaystyle +\,{\it\epsilon}\left[\frac{\hat{\boldsymbol{b}}}{eB_{\Vert }^{\star }}\times \left(\frac{{\it\mu}}{{\it\gamma}}\boldsymbol{{\rm\nabla}}B+\frac{p_{\Vert }^{2}}{{\it\gamma}m}{\bf\kappa}\right)\right]\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{B}+O({\it\epsilon}^{2}).\end{eqnarray}$$
We also note that
${\mathcal{T}}_{gc}^{-1}({\it\nu}_{r}{\it\Omega}^{-1}B^{-1})\equiv {\it\nu}_{r}{\it\Omega}^{-1}B^{-1}$
.
Combining the above expressions, we calculate the push-forward of the particle radiation–reaction force
$$\begin{eqnarray}\displaystyle {\mathcal{T}}_{gc}^{-1}\boldsymbol{K} & = & \displaystyle -({\mathcal{T}}_{gc}^{-1}{\it\nu}_{r})\left[({\mathcal{T}}_{gc}^{-1}\boldsymbol{p}_{\bot })+\frac{({\mathcal{T}}_{gc}^{-1}p_{\bot }^{2})}{(mc)^{2}}({\mathcal{T}}_{gc}^{-1}\boldsymbol{p})\right]\nonumber\\ \displaystyle & & \displaystyle -\,{\it\epsilon}{\it\nu}_{r}{\it\Omega}^{-1}B^{-1}({\mathcal{T}}_{gc}^{-1}\dot{\boldsymbol{B}})\times ({\mathcal{T}}_{gc}^{-1}\boldsymbol{p}).\end{eqnarray}$$
Expressed as
${\mathcal{T}}_{gc}^{-1}\boldsymbol{K}=\boldsymbol{K}_{0}+{\it\epsilon}\boldsymbol{K}_{1}+O({\it\epsilon}^{2})$
, the zeroth-order term is given by
$$\begin{eqnarray}\boldsymbol{K}_{0}=-{\it\nu}_{r}\left(eB\frac{\partial {\bf\rho}_{0}}{\partial {\it\theta}}+\frac{2{\it\mu}B}{mc^{2}}\boldsymbol{p}_{0}\right),\end{eqnarray}$$
and for the first-order term we have
$$\begin{eqnarray}\displaystyle \boldsymbol{K}_{1} & = & \displaystyle 2({\bf\rho}_{0}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\ln B)\boldsymbol{K}_{0}-{\it\nu}_{r}{\it\Omega}^{-1}B^{-1}\left(\frac{\boldsymbol{p}_{0}}{{\it\gamma}m}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{B}\right)\times \boldsymbol{p}_{0}\nonumber\\ \displaystyle & & \displaystyle -\,{\it\nu}_{r}\left[\left(1+\frac{2{\it\mu}B}{mc^{2}}\right)\boldsymbol{p}_{1}+G_{1}^{p_{\Vert }}\hat{\boldsymbol{b}}-p_{\Vert }{\bf\rho}_{0}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{b}}+\frac{2{\it\mu}B}{mc^{2}}\left({\bf\rho}_{0}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\ln B-\frac{G_{1}^{{\it\mu}}}{{\it\mu}}\right)\boldsymbol{p}_{0}\right]\!.\qquad \quad\end{eqnarray}$$
Appendix C. Gyro-averages of the guiding-centre radiation–reaction force
With the push-forward of the particle phase-space RR-force given in appendix B, we can evaluate the necessary gyro-averages with the help of a useful identity
$$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{b}}:\langle \hat{\boldsymbol{\bot }}\hat{{\bf\rho}}\rangle =-\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{b}}:\langle \hat{{\bf\rho}}\hat{\boldsymbol{\bot }}\rangle =\frac{{\it\tau}_{B}}{2}.\end{eqnarray}$$
This helps us evaluate the gyro-averages of the push-forwarded RR-force:
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \boldsymbol{K}_{0}\rangle =-{\it\nu}_{r}p_{\Vert }\,\frac{2{\it\mu}B}{mc^{2}}\,\hat{\boldsymbol{b}},\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \boldsymbol{K}_{1}\rangle ={\it\nu}_{r}\frac{2{\it\mu}B}{mc^{2}}\frac{\hat{\boldsymbol{b}}}{eB_{\Vert }^{\star }}\times (m{\it\mu}\boldsymbol{{\rm\nabla}}B+p_{\Vert }^{2}{\bf\kappa})-3{\it\nu}_{r}\frac{2{\it\mu}B}{mc^{2}}p_{\Vert }\hat{\boldsymbol{b}}\times ({\it\varrho}_{\Vert }{\bf\kappa})-{\it\nu}_{r}\frac{2{\it\mu}B}{mc^{2}}p_{\Vert }{\it\varrho}_{\Vert }{\it\tau}_{B}\hat{\boldsymbol{b}},\qquad & \displaystyle\end{eqnarray}$$
${\it\varrho}_{\Vert }=p_{\Vert }/(eB)$
, as well as the rest of the gyro-averages that are needed in the expressions for
${\mathcal{K}}^{\boldsymbol{X}}$
,
${\mathcal{K}}^{p_{\Vert }}$
, and
${\mathcal{K}}^{{\it\mu}}$
: 
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \boldsymbol{{\rm\nabla}}^{\star }{\bf\rho}_{0}\boldsymbol{\cdot }\boldsymbol{K}_{0}\rangle =-{\it\nu}_{r}{\it\tau}_{B}\frac{m{\it\mu}}{e}\left(1+\frac{2{\it\mu}B}{mc^{2}}\right)\hat{\boldsymbol{b}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left\langle \frac{\partial {\bf\rho}_{1}}{\partial p_{\Vert }}\boldsymbol{\cdot }\boldsymbol{K}_{0}\right\rangle =0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left\langle \frac{\partial {\bf\rho}_{0}}{\partial {\it\theta}}\boldsymbol{\cdot }\boldsymbol{K}_{0}\right\rangle =-2{\it\nu}_{r}\frac{m{\it\mu}}{e}\left(1+\frac{2{\it\mu}B}{mc^{2}}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left\langle \frac{\partial {\bf\rho}_{0}}{\partial {\it\theta}}\boldsymbol{\cdot }\boldsymbol{K}_{1}\right\rangle =-2{\it\nu}_{r}\frac{m{\it\mu}}{e}\left(1+\frac{2{\it\mu}B}{mc^{2}}\right)\,{\it\varrho}_{\Vert }{\it\tau}_{B}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left\langle \frac{\partial {\bf\rho}_{1}}{\partial {\it\theta}}\boldsymbol{\cdot }\boldsymbol{K}_{0}\right\rangle ={\it\nu}_{r}\frac{m{\it\mu}}{e}\left(1+\frac{2{\it\mu}B}{mc^{2}}\right){\it\varrho}_{\Vert }{\it\tau}_{B}. & \displaystyle\end{eqnarray}$$
Now, the spatial component of the guiding-centre radiation–reaction force becomes
$$\begin{eqnarray}\displaystyle {\mathcal{K}}^{\boldsymbol{X}} & = & \displaystyle -{\it\epsilon}\frac{\hat{\boldsymbol{b}}}{eB_{\Vert }^{\star }}\times \langle \boldsymbol{K}_{0}+{\it\epsilon}\boldsymbol{K}_{1}+{\it\epsilon}\boldsymbol{{\rm\nabla}}^{\star }{\bf\rho}_{0}\boldsymbol{\cdot }\boldsymbol{K}_{0}\rangle -{\it\epsilon}^{2}\hat{\boldsymbol{b}}\left\langle \frac{\partial {\bf\rho}_{1}}{\partial p_{\Vert }}\boldsymbol{\cdot }\boldsymbol{K}_{0}\right\rangle ,\nonumber\\ \displaystyle & = & \displaystyle -{\it\epsilon}^{2}\,\frac{{\it\nu}_{r}}{{\it\Omega}_{\Vert }^{\star }}\frac{2{\it\mu}B}{mc^{2}}(\hat{\boldsymbol{b}}\times \dot{\boldsymbol{X}}+3v_{\Vert }\,{\it\varrho}_{\Vert }\,{\bf\kappa}),\end{eqnarray}$$
where we have introduced the modified gyro frequency
$$\begin{eqnarray}{\it\Omega}_{\Vert }^{\star }=(eB_{\Vert }^{\star })/({\it\gamma}m)={\it\Omega}\,(1+{\it\epsilon}{\it\varrho}_{\Vert }{\it\tau}_{B}).\end{eqnarray}$$
For the parallel momentum component we find
$$\begin{eqnarray}\displaystyle {\mathcal{K}}^{p_{\Vert }} & = & \displaystyle \frac{\boldsymbol{B}^{\star }}{B_{\Vert }^{\star }}\boldsymbol{\cdot }\langle \boldsymbol{K}_{0}\rangle +{\it\epsilon}\hat{\boldsymbol{b}}\boldsymbol{\cdot }\langle \boldsymbol{K}_{1}+\boldsymbol{{\rm\nabla}}^{\star }{\bf\rho}_{0}\boldsymbol{\cdot }\boldsymbol{K}_{0}\rangle ,\nonumber\\ \displaystyle & = & \displaystyle -{\it\nu}_{r}p_{\Vert }\frac{{\it\mu}B}{mc^{2}}(2+{\it\epsilon}{\it\varrho}_{\Vert }{\it\tau}_{B})-{\it\epsilon}{\it\nu}_{r}\frac{m{\it\mu}{\it\gamma}^{2}}{e}{\it\tau}_{B},\end{eqnarray}$$
and for the magnetic moment
${\it\mu}$
the force becomes
$$\begin{eqnarray}\displaystyle {\mathcal{K}}^{{\it\mu}} & = & \displaystyle \frac{e}{m}\left\langle \frac{\partial {\bf\rho}_{0}}{\partial {\it\theta}}\boldsymbol{\cdot }(\boldsymbol{K}_{0}+{\it\epsilon}\boldsymbol{K}_{1})\right\rangle +{\it\epsilon}\frac{e}{m}\left\langle \frac{\partial {\bf\rho}_{1}}{\partial {\it\theta}}\boldsymbol{\cdot }\boldsymbol{K}_{0}\right\rangle ,\nonumber\\ \displaystyle & = & \displaystyle -{\it\nu}_{r}{\it\mu}\left(1+\frac{2{\it\mu}B}{mc^{2}}\right)(2+{\it\epsilon}{\it\varrho}_{\Vert }{\it\tau}_{B}).\end{eqnarray}$$
