2. Background and estimates for trapped alpha particle transport in a nearly QS field

To illustrate why the new regime arises, consider the model trapped alpha, bounce or transit averaged drift kinetic equation

(2.1)\begin{equation}{\bar{\omega }_\alpha }\frac{{\partial \tilde{f}}}{{\partial \varphi }} - \bar{V}\sin \varphi \frac{{\partial (\bar{f} + \tilde{f})}}{{\partial r}} = \bar{\nu }\varepsilon \frac{{{\partial ^2}\tilde{f}}}{{\partial {\lambda ^2}}}\end{equation}

for a nearly QS, large aspect ratio stellarator. Here, $\bar{f}$ and $\tilde{f}$ are the unperturbed and perturbed alpha distribution functions; $\bar{V}\sin \varphi \partial \bar{f}/\partial r$ is the drive term associated with the departure from quasisymmetry with $\bar{V}$ the radial drift due to $\varphi$ the symmetry breaking, non-QS angular variation drive; ${\bar{\omega }_\alpha }$ is the transit averaged tangential drift of the alphas in a flux surface that reverses direction at some pitch angle (Galeev & Sagdeev Reference Galeev and Sagdeev1979) to be defined in detail shortly; $\bar{\nu }$ is the pitch angle scattering frequency of the trapped alphas by the background ions; $\lambda = 2\mu {B_0}/{v^2}$ is the adiabatic invariant pitch angle variable, with $\mu$ the magnetic moment, $v = |\boldsymbol{v}|$ the speed of the alphas and ${B_0}$ a normalizing magnetic field; and r is the minor radius flux surface label. The inverse aspect ratio $\varepsilon = r/{R_0} \ll 1$ on the right side of (2.1) accounts for the width in pitch angle, ${\varepsilon ^{1/2}}$, of the trapped region of velocity space, with ${R_0}$ the major radius. The pitch angle scattering frequency times the slowing down time, ${\tau _s}$, is small for the alphas with $\bar{\nu }{\tau _s}\sim v_\lambda ^3/v_0^3 \ll 1$, where ${v_0}$ the alpha birth speed and ${v_\lambda }$ the speed at which pitch angle scattering enters for a non-resonant birth alpha (the critical speed ${v_c}$ is where electron and ion drag are comparable and ${v_\lambda }\sim {v_c}$).

At large aspect ratio ${\bar{\omega }_\alpha } = {\bar{\omega }_\alpha }({\kappa ^2})$ vanishes on a flux surface when $\kappa _0^2 = 0.83$ (Galeev & Sagdeev Reference Galeev and Sagdeev1979) assuming magnetic shear is negligible, where ${\kappa ^2} = [1 - (1 - \varepsilon )\lambda ]/2\varepsilon \lambda$. However, while $\lambda$ is an adiabatic invariant, ${\kappa ^2}$ is not because of its $\varepsilon$ dependence, which appears because the trapped–passing boundary depends on inverse aspect ratio $\varepsilon \ll 1$, with $1/(1 + \varepsilon ) < \lambda < 1/(1 - \varepsilon )$ for the trapped alphas. As a result, the resonant pitch angle ${\lambda _0}$ depends on the inverse aspect ratio $\varepsilon$ of the flux surface of interest according to ${\lambda _0} = 1/[1 + (2\kappa _0^2 - 1)\varepsilon ] = 1/(1 + 0.66\varepsilon )$. In the vicinity of the resonance

(2.2)\begin{equation}{\bar{\omega }_\alpha } ={-} 2({\kappa ^2} - \kappa _0^2){\bar{\omega }^{\prime}_\alpha },\end{equation}

with ${\bar{\omega }^{\prime}_\alpha }\sim {\bar{\omega }_\alpha }\sim q{v^2}/{\varOmega _0}R_0^2\varepsilon$. The on axis gyrofrequency of the alphas is denoted by ${\varOmega _0}$, and q is the safety factor. Then, ${\kappa ^2} - \kappa _0^2$ depends on $\varepsilon$ as well as $\lambda$, as is seen from

(2.3)\begin{equation}{\kappa ^2} - \kappa _0^2 = \frac{{[1 - (1 - \varepsilon )\lambda - 2\kappa _0^2\varepsilon \lambda ]}}{{2\varepsilon \lambda }} \approx{-} \frac{{\lambda - [1 - (2\kappa _0^2 - 1)\varepsilon ]}}{{2\varepsilon }},\end{equation}

leading to a more insightful form of the model equation with the nonlinear term retained:

(2.4)\begin{equation}[\lambda - (1 + \varepsilon - 2\kappa _0^2\varepsilon )]\frac{{{{\bar{\omega }^{\prime}}_\alpha }}}{\varepsilon }\frac{{\partial \tilde{f}}}{{\partial \varphi }} - \frac{{\bar{V}}}{{{R_0}}}\sin \varphi \; \frac{{\partial (\bar{f} + \tilde{f})}}{{\partial \varepsilon }} = \bar{\nu }\varepsilon \frac{{{\partial ^2}\tilde{f}}}{{\partial {\lambda ^2}}}.\end{equation}

Notice that if ${\kappa ^2} - \kappa _0^2 = \Delta {\kappa ^2}\sim 1$, then $\lambda - [1 - (2\kappa _0^2 - 1)\varepsilon ]\sim \varepsilon \ll 1$ to account for the $\varepsilon$ factor in the denominator of (2.3). Moreover, each flux surface $\varepsilon$ has a slightly different resonant pitch angle ${\lambda _0}(\varepsilon )$, implying each island structure in normalized minor radius $\varepsilon$ is the cross-section of a tube or pod in phase space since the resonance also involves pitch angle ‘islands’. Neighbouring flux surfaces have slightly different resonant pitch angles so the phase space ‘pods’ extend from the inner (larger ${\lambda _0}$) to outer (smaller ${\lambda _0}$) flux surfaces.

For each and every flux surface $0 \le \varepsilon \ll 1$, there is a resonant pitch angle ${\lambda _0}(\varepsilon )$. Therefore, the phase space pods extend over the entire minor radius for $\varepsilon \ll 1$. They are referred to as pods since they have an island width in the pitch angle variable $\lambda$ as well as in the normalized radial or flux surface variable $\varepsilon$ in the presence of a departure from quasisymmetry due to the error field angle dependence denoted by $\varphi$. Consequently, the phase space pods considered herein are not quite the usual islands observed in Poincaré plots of collisionless particle motion (Paul et al. Reference Paul, Bhattacharjee, Landreman, Alex, Velasco and Nies2022; White Reference White2022).

Defining $\Delta \varepsilon = \varepsilon - (1 - \lambda )/(2\kappa _0^2 - 1) \ll 1$ at fixed $\lambda$ and balancing drift and nonlinear terms gives ${\bar{\omega }^{\prime}_\alpha }\Delta \varepsilon /\varepsilon \sim \bar{V}/{R_0}\Delta \varepsilon$, leading to a narrow radial island width estimate of

(2.5)\begin{equation}\Delta \varepsilon \sim {(\bar{V}\varepsilon /{\bar{\omega }_\alpha }{R_0})^{1/2}}\sim {(\varepsilon \delta )^{1/2}} \ll {\varepsilon ^{1/2}},\end{equation}

when the radial drift $\bar{V}\sim {\bar{\omega }_\alpha }{R_0}\delta \;$ is due to a very small normalized departure from quasisymmetry, $\delta = {B_{||}}/{B_0} \ll \varepsilon$.

Resonant plateau behaviour is found by solving the simpler linearized equation

(2.6)\begin{equation}(\lambda - {\lambda _0})\frac{{{{\bar{\omega }^{\prime}}_\alpha }}}{\varepsilon }\frac{{\partial \tilde{f}}}{{\partial \varphi }} - \frac{{\bar{V}}}{{{R_0}}}\sin \varphi \; \frac{{\partial \bar{f}}}{{\partial \varepsilon }} = \bar{\nu }\varepsilon \frac{{{\partial ^2}\tilde{f}}}{{\partial {\lambda ^2}}},\end{equation}

with ${\lambda _0} = 1 - (2\kappa _0^2 - 1)\varepsilon \approx 1 - 0.66\varepsilon$. In this limit, the nonlinear term is small because the fine-scale radial variation of $\tilde{f}$ is removed by the collisional boundary layer so no phase space pod formation occurs. The solution leads to resonant plateau transport or what is often referred to as superbanana plateau transport and is expected to be a dominant collisional loss mechanism for birth alphas (Galeev & Sagdeev Reference Galeev and Sagdeev1979; Shaing Reference Shaing2015; Catto Reference Catto2019). The width of the resonance, $\Delta \lambda = \lambda - {\lambda _0}$ at fixed $\varepsilon$, is estimated by balancing drift and collision terms, ${\bar{\omega }^{\prime}_\alpha }\Delta \lambda /\varepsilon \sim \bar{\nu }\varepsilon /{(\Delta \lambda )^2}$, thereby giving a resonant plateau boundary layer width:

(2.7)\begin{equation}\Delta \lambda \sim {({\varepsilon ^2}\bar{\nu }/{\bar{\omega }_\alpha })^{1/3}},\end{equation}

and an effective resonant plateau collision frequency:

(2.8)\begin{equation}{\bar{\nu }_{\textrm{eff}}}\sim \bar{\nu }\varepsilon /{(\Delta \lambda )^2}\sim \bar{\nu }\varepsilon {({\bar{\omega }_\alpha }/{\varepsilon ^2}\bar{\nu })^{2/3}}.\end{equation}

The resonant plateau (or superbanana plateau) diffusivity that results for birth alphas is then

(2.9)\begin{equation}{D_{\textrm{rp}}}\sim (\Delta \lambda /{\varepsilon ^{1/2}}){(\bar{V}/{\bar{\nu }_{\textrm{eff}}})^2}{\bar{\nu }_{\textrm{eff}}}\sim {\varepsilon ^{1/2}}{\bar{V}^2}/{\bar{\omega }_\alpha }\sim qv_0^2{\delta ^2}/{\varOmega _0}{\varepsilon ^{1/2}},\end{equation}

where $\bar{V}/{\bar{\nu }_{\textrm{eff}}}$ is the step size with $\bar{V}\sim {\bar{\omega }_\alpha }{R_0}\delta$, and $\Delta \lambda /{\varepsilon ^{1/2}}$ is the effective resonant trapped fraction of the collisional boundary layer.

Phase space pods do not form in the resonant plateau regime because the departure from quasisymmetry is so small that $\partial \tilde{f}/\partial \varepsilon \sim \tilde{f}/\Delta \varepsilon \ll \partial \bar{f}/\partial \varepsilon \sim \bar{f}/\varepsilon$. In this $\Delta \lambda \gg \Delta \varepsilon$ limit, collisions prevent phase space pod formation. Therefore, resonant plateau behaviour requires

(2.10)\begin{equation}{({\varepsilon ^2}\bar{\nu }/{\bar{\omega }_\alpha })^{1/3}} \gg {(\varepsilon \delta )^{1/2}},\end{equation}

giving the restriction

(2.11)\begin{equation}\tilde{f}/\bar{f} \ll \; \Delta \varepsilon /\varepsilon \sim {(\delta /\varepsilon )^{1/2}} \ll \Delta \lambda /\varepsilon \sim {(\bar{\nu }/{\bar{\omega }_\alpha }\varepsilon )^{1/3}},\end{equation}

where for pitch angle scatter, $\bar{\nu }{\tau _s}\sim v_\lambda ^3/v_0^3 \ll 1$, with ${\tau _s}$ the alpha slowing down time due to electron drag. For larger departures from quasisymmetry the preceding can be difficult to satisfy since $\bar{\nu }/{\bar{\omega }_\alpha }\varepsilon \sim {10^{ - 4}}$ for birth alphas with ${R_0}/{v_0}{\tau _s}\sim {10^{ - 5}}$ and ${\varOmega _0}{R_0}/q{v_0}\sim {10^2}$. A less well optimized stellarator might only have $1 \gg \delta /\varepsilon \mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel> \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }{10^{ - 2}}$.

The preceding estimates indicate that the maximal ordering of interest for a well optimized QS stellarator is to allow $\Delta \varepsilon \sim \Delta \lambda$ or

(2.12)\begin{equation}{(\varepsilon \delta )^{1/2}}\sim {({\varepsilon ^2}\bar{\nu }/{\bar{\omega }_\alpha })^{1/3}}.\end{equation}

Unfortunately, finding a completely general analytic solution in this limit is impractical. However, by assuming collisions are sufficiently weak and the departure from quasisymmetry sufficiently large, the nearly collisionless, reasonably well confined motion of an approximately QS field is perturbed by the phase space pod structure in a way that can be evaluated using a reduced constant of the motion. In this limit the phase space pods remain well defined because the very narrow collisional boundary layer is only at the pod boundary (or separatrix). The bound and barely circulating alphas are aware of the phase space tube or double island structure in their motion, but the weak collisions result in small, speed independent spatial steps of ${R_0}{(\varepsilon \delta )^{1/2}}$ due to a local flattening of the perturbed distribution function.

To avoid confusion, the alpha motion inside a phase space pod or tube is herein referred to as bound or librating, rather than trapped (since the alphas are already trapped in the nearly QS magnetic field of the stellarator). Due to collisions, the bound motion results in a local flattening of the perturbed distribution function in the phase space pod with a radial tube or pod width of ${R_0}{(\varepsilon \delta )^{1/2}}$. Only a constraint from the lowest order collision operator need be evaluated to determine this distribution function to lowest order, by matching to the pitch angle region away from the separatrices of the phase space pod. Alpha motion outside the phase space tube or pod is referred to as unbound or circulating so passing is reserved to refer to alpha motion in the nearly QS stellarator field.

In the finite pod width limit, the radial diffusivity is evaluated in detail in the following sections in the weak collisionality limit. A simple estimate that is consistent with this evaluation and the presence of a very narrow collisional boundary layer about the separatrices is obtained by taking the radial step size as ${R_0}{(\varepsilon \delta )^{1/2}}$, and the effective fraction as the bound and barely circulating fraction of the trapped fraction ${(\varepsilon \delta /\varepsilon )^{1/2}}$. In addition, the lowest order alpha motion is essentially collisionless with the alphas slowing to spend more time in the vicinity of the phase space pods. As a result, these alphas are acted on by electron drag as well as the strong pitch angle scattering associated with the separatrix boundary layer. Consequently, in the presence of finite pods, birth alphas are not depleted immediately. As a result, lower speed alphas as well as birth alphas contribute to the radial energy transport, thereby resulting in the replacement $\bar{\nu }{\partial ^2}/\partial {\lambda ^2}\sim v_\lambda ^3/{\tau _s}v_0^2{v_c}\varepsilon \delta$, giving ${\bar{\nu }_{\textrm{eff}}}\sim v_\lambda ^3/{\tau _s}v_0^2{v_c}\delta$ for the phase space pod limit considered herein. Not surprisingly, the bound resonant orbits become fully depleted, flattening the alpha distribution function in the pod, and resulting in only the barely circulating alphas in the vicinity of the separatrix contributing to transport. Based on these estimates, the energy diffusivity in the weak collisionality, large departure from quasisymmetry limit is of order

(2.13)\begin{equation}{D_\nu }\sim {\delta ^{1/2}}{[{R_0}{(\varepsilon \delta )^{1/2}}]^2}(v_\lambda ^3/{v_c}v_0^2{\tau _s}\delta )\sim {\delta ^{1/2}}\varepsilon R_0^2v_\lambda ^3/v_0^2{v_c}{\tau _s}.\end{equation}

The preceding estimate is verified by the detailed calculations performed in the following sections.

In general, the ratio of these two diffusivities,

(2.14)\begin{equation}\frac{{{D_\nu }}}{{{D_{\textrm{rp}}}}}\sim {\left( {\frac{\varepsilon }{\delta }} \right)^{3/2}}\frac{{{\varOmega _0}R_o^2v_\lambda ^3}}{{qv_0^4{v_c}{\tau _s}}},\end{equation}

is allowed to be comparable to or less than one. For example, they are comparable and give large transport $({\tau _s}D/{r^2}\sim 1)$ for an optimized stellarator with $\varepsilon \sim 0.1$ and $\delta \sim {10^{ - 2}}$, and for birth alphas with ${R_0}/{v_0}{\tau _s}\sim {10^{ - 5}}$ and ${\varOmega _0}{R_0}/q{v_0}\sim {10^2}$ (for which ${(\varepsilon \delta )^{1/2}}\sim {({\varepsilon ^2}\bar{\nu }/{\bar{\omega }_\alpha })^{1/3}}$ is satisfied within a factor of ${v_0}/{v_c}\sim 3$). As highly optimized QS stellarator coil sets will be difficult to fabricate and the diffusivity ${D_\nu }$ associated with the departure from quasisymmetry acts to remove lower speed alphas, the analytically tractable limit ${(\varepsilon \delta )^{1/2}} \gg {({\varepsilon ^2}\bar{\nu }/{\bar{\omega }_\alpha })^{1/3}}$ is the focus of the remaining sections.

3. Reduced kinetic equation for trapped alpha particles in a nearly QS field

The trapped alpha particles satisfy the reduced drift kinetic equation in Catto et al. (Reference Catto, Tolman and Parra2023) (see also Catto Reference Catto2019) that assumes the QS field depends on the angle variable $\eta = M\vartheta - N\zeta$, with $\vartheta$ and $\zeta$ the Boozer (Reference Boozer1981) poloidal and toroidal angle variables, respectively. Normally $M = 1$ (Cary & Shasharina Reference Cary and Shasharina1997; Landreman & Paul Reference Landreman and Paul2022), since only the region near the magnetic axis is of interest, and N is an integer ($N = 0$ for quasiaxisymmetry). For well defined flux surfaces, the magnetic field is

(3.1)\begin{equation}\boldsymbol{B} = \boldsymbol{\nabla }{\psi _t} \times \boldsymbol{\nabla }\vartheta + \boldsymbol{\nabla }\zeta \times \boldsymbol{\nabla }{\psi _p} = K(\boldsymbol{r})\boldsymbol{\nabla }{\psi _p} + G({\psi _p})\boldsymbol{\nabla }\vartheta + I({\psi _p})\boldsymbol{\nabla }\zeta .\end{equation}

In the preceding, ${\psi _t}$ and ${\psi _p}$ are the toroidal and poloidal flux functions, respectively, and are related by $\partial {\psi _t}/\partial {\psi _p} = q$, and the flux functions $I({\psi _p})$ and $G({\psi _p})$ are related to the poloidal and toroidal currents (Boozer Reference Boozer1981). Constant QS B contours or curves close on themselves when the magnetic field line label $\alpha = (\zeta - q\vartheta )$ changes by $2\mathrm{\pi }(M - qN)$ at fixed $\eta$ as $\vartheta$ and $\zeta$ change by $2\mathrm{\pi }N$ and $2\mathrm{\pi }M$, respectively. Denoting the two forms of $\boldsymbol{B}$ together gives ${B^2}/(qI + G) = \boldsymbol{\nabla }\vartheta \times \boldsymbol{\nabla }\zeta \boldsymbol{\cdot }\boldsymbol{\nabla }{\psi _p} = \boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla }\eta /(M - qN)$ with $(qI + G)/B \approx q{R_0}$ at large aspect ratio.

At large aspect ratio $\partial /\partial {\psi _p} \approx (2\mathrm{\pi }q/{B_0}a^{\prime})\partial /\partial r$, with $a(r)$ the area of the elliptical flux surfaces nearest the magnetic axis, where the magnetic field is assumed to be of the form

(3.2)\begin{equation}B = {B_0}[1 - \varepsilon \cos \eta + \delta \cos (m\vartheta - n\zeta )] = B = {B_0}(1 - \varepsilon \cos \eta ) + {B_{||}}(\eta ,\alpha ),\end{equation}

with ${B_{||}} = {B_0}\delta \cos (m\vartheta - n\zeta )$ a small departure from quasisymmetry ($\delta \ll \varepsilon$, with $m \ne M$ and $n \ne N$), and $\chi = m\vartheta - n\zeta = [(mN - nM)\alpha + \eta (m - qn)]/(M - qN)$.

The transit or bounce averaged kinetic equation for the trapped alphas assumes $\partial f/\partial \eta = 0$ lowest order. Then bounce averaging the nonlinear kinetic equation and using the periodicity of the trapped alpha motion gives

(3.3)\begin{equation}({\oint_\alpha {\textrm{d}\tau {\omega_\alpha }} } )\frac{{\partial f}}{{\partial \alpha }} + \left( {\int_\alpha {\textrm{d}\tau {\boldsymbol{v}_d}\boldsymbol{\cdot }\boldsymbol{\nabla }{\psi_p}} } \right)\; \frac{{\partial f}}{{\partial {\psi _p}}} = \oint_\alpha {\textrm{d}\tau } \left[ {C\{ f\} + \frac{{S({\psi_p})\delta (v - {v_0})}}{{4\mathrm{\pi }{v^2}}}} \right].\end{equation}

The transit average is performed at fixed $\alpha$ using $\textrm{d}\eta /\textrm{d}\tau = {v_{||}}\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla }\eta$ with ${v_{||}}$ the parallel alpha velocity, and $\textrm{d}\tau > 0$ the incremental change in the time along the trapped trajectory.

In the radial magnetic drift,

(3.4)\begin{equation}{\boldsymbol{v}_d}\boldsymbol{\cdot }\boldsymbol{\nabla }{\psi _p} = \left( {\frac{{MI + NG}}{{M - qN}}} \right){v_{||}}\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla }\left( {\frac{{{v_{||}}}}{\boldsymbol{\varOmega }}} \right) - {v_{||}}B\frac{\partial }{{\partial \alpha }}\left( {\frac{{{v_{||}}}}{\boldsymbol{\varOmega }}} \right),\end{equation}

the first term is the neoclassical transport drive term (Landreman & Catto Reference Landreman and Catto2012). It is unaffected by a small departure from quasisymmetry as no resonance occurs, and it is assumed to be unimportant compared with the second drive term due to the departure from quasisymmetry. The collision operator is denoted by $C\{ f\}$, $\varOmega = {Z_\alpha }eB/{M_\alpha }c$ is the alpha gyrofrequency (with ${Z_\alpha }$ and ${M_\alpha }$ the charge number and mass, e the charge on a proton and c the speed of light) and S is the isotropic birth rate of alphas born at speed ${v_0}$ as indicated by the delta function.

Shear is normally weak in optimized stellarators (Landreman & Paul Reference Landreman and Paul2022) and is ignored here. In the absence of magnetic shear, the trapped $(\lambda \approx 1)$ alpha drift on a QS flux surface for a large aspect ratio stellarator is adequately approximated by

(3.5)\begin{equation}{\boldsymbol{v}_d}\boldsymbol{\cdot }\boldsymbol{\nabla }\alpha = {\omega _\alpha } \approx {v_{||}}B\frac{\partial }{{\partial {\psi _p}}}\left( {\frac{{{v_{||}}}}{\varOmega }} \right) \approx{-} \frac{{\mathrm{\pi }q{v^2}\cos \eta }}{{{\varOmega _0}{R_0}a^{\prime}}},\end{equation}

resulting in

(3.6)\begin{equation}\oint_\alpha {\textrm{d}\tau } = \frac{{8q{R_0}K(\kappa )}}{{(M - qN)v\sqrt {2\varepsilon } }}\sim \frac{{q{R_0}}}{{v\sqrt \varepsilon }}\end{equation}

and

(3.7)\begin{equation}\oint_\alpha {\textrm{d}\tau {\omega _\alpha }} = \frac{{8\mathrm{\pi }{q^2}v[2E(\kappa ) - K(\kappa )]}}{{(M - qN){\varOmega _0}a^{\prime}\sqrt {2\varepsilon } }},\end{equation}

where ${\kappa ^2} = [1 - (1 - \varepsilon )\lambda ]/2\varepsilon \lambda$ and $\textrm{d}\tau \approx q{R_0}\textrm{d}\eta /(M - qN){v_{||}}$. Drift reversal occurs at $2E(\kappa ) = K(\kappa )$, when $\kappa _0^2 = 0.83$. The resonant $\lambda$ depends on $\varepsilon$ and satisfies $1/(1 + \varepsilon ) < {\lambda _0} = 1/[1 + (2\kappa _0^2 - 1)\varepsilon ] = 1/(1 + 0.66\varepsilon ) < 1 < 1/(1 - \varepsilon )$.

Next, using $v_{||}^2 = \; {\xi ^2}{v^2} = {v^2}(1 - \lambda B/{B_0})$ and $\varepsilon \ll 1$ gives

(3.8)\begin{align}\oint_\alpha {\textrm{d}\tau {\boldsymbol{v}_d}\boldsymbol{\cdot }\boldsymbol{\nabla }{\psi _p}} \to - \oint_\alpha {\textrm{d}\tau {v_{||}}B\frac{\partial }{{\partial \alpha }}\left( {\frac{{{v_{||}}}}{\varOmega }} \right)} = {v^2}\oint_\alpha {\frac{{\textrm{d}\tau }}{\varOmega }\left( {1 - \frac{{\lambda B}}{{2{B_0}}}} \right)\frac{{\partial {B_{||}}}}{{\partial \alpha }}} \approx \frac{{{v^2}}}{{2{\varOmega _0}}}\frac{\partial }{{\partial \alpha }}\oint_\alpha {\textrm{d}\tau {B_{||}}} ,\end{align}

where $\lambda B/{B_0} \approx 1$ is used for the trapped alphas. Defining $\eta = 0$ as the bottom of the magnetic well, and using $\int_\alpha {\textrm{d}\tau \sin [\eta (m - qn)/(M - qN)]} = 0$, only the cosine terms matter in

(3.9)\begin{align}\begin{aligned} {B_{||}} & = {B_0}\delta \{ \cos [\eta (m - qn)/(M - qN)\left] {\cos } \right[\alpha (mN - nM)/(M - qN)]\\ & \quad - \sin [\eta (m - qn)/(M - qN)\left] {\sin } \right[\alpha (mN - nM)/(M - qN)]\} . \end{aligned}\end{align}

Therefore, it is convenient to define $\varphi = p\alpha$ with $p = |(mN - nM)/(M - qN)|$ and

(3.10)\begin{equation}\varTheta = {{\oint_\alpha {\textrm{d}\tau \cos [\eta |(m - qn)/(M - qN)|]} } / {\oint_\alpha {\textrm{d}\tau } }},\end{equation}

along with ${\varOmega _0} = {Z_\alpha }e{B_0}/{M_\alpha }c$, to obtain the radial drift term:

(3.11)\begin{equation}\oint_\alpha {\textrm{d}\tau {\boldsymbol{v}_d}\boldsymbol{\cdot }\boldsymbol{\nabla }{\psi _p}} = \frac{{ - p{B_0}{v^2}\delta }}{{2{\varOmega _0}}}\sin \varphi \oint_\alpha {\textrm{d}\tau \cos \left( {\eta \left|{\frac{{m - qn}}{{M - qN}}} \right|} \right)} = \frac{{ - p{B_0}{v^2}\varTheta \delta }}{{2{\varOmega _0}}}\sin \varphi \oint_\alpha {\textrm{d}\tau } .\end{equation}

The full transit averaged form of the nonlinear reduced drift kinetic equation for the trapped alphas is then obtained by defining the transit average

(3.12)\begin{equation}\overline {( \cdots )} = {{\oint_\alpha {\textrm{d}\tau ( \cdots )} } / {\oint_\alpha {\textrm{d}\tau } }}\end{equation}

to obtain the form

(3.13)\begin{equation}{\bar{\omega }_\alpha }\frac{{\partial \tilde{f}}}{{\partial \varphi }} - \frac{{\bar{V}}}{{{R_0}}}\sin \varphi \; \frac{{\partial (\bar{f} + \tilde{f})}}{{\partial \varepsilon }} = \frac{{\bar{C}\{ \,\tilde{f}\} }}{p},\end{equation}

with

(3.14)\begin{align}{\bar{\omega }_\alpha } = \frac{{\mathrm{\pi }q{v^2}[2E(\kappa ) - K(\kappa )]}}{{{\varOmega _0}{R_o}a^{\prime}K(\kappa )}} \approx{-} \frac{{\mathrm{\pi }q{v^2}({\kappa ^2} - \kappa _0^2)}}{{4\kappa _0^2(1 - \kappa _0^2){\varOmega _0}{R_o}a^{\prime}}} = \frac{{\mathrm{\pi }q{v^2}[\lambda - (1 + \varepsilon - 2\kappa _0^2\varepsilon )]}}{{8\varepsilon \kappa _0^2(1 - \kappa _0^2){\varOmega _0}{R_o}a^{\prime}}}\end{align}

and

(3.15)\begin{equation}\bar{V} = \frac{{\mathrm{\pi }q{v^2}\varTheta \delta }}{{{\varOmega _0}a^{\prime}}}.\end{equation}

In ${\bar{\omega }_\alpha }$ an expansion around the drift reversal pitch angle $\kappa _0^2$ is performed and ${\lambda _0} \approx 1 - (2\kappa _0^2 - 1)\varepsilon$ inserted. Only the final form of (3.14) is required here. This resonance (3.14) depends on $\lambda$ and $\varepsilon$, and not just flux as assumed in earlier treatments (Hazeltine & Catto Reference Hazeltine and Catto1981; Shaing & Hsu Reference Shaing and Hsu2014; Shaing Reference Shaing2015) using truncated Taylor expansions about a reference flux surface.

The full alpha distribution function is now written as $f = \bar{f} + \tilde{f}$, with $\bar{f} \gg \tilde{f}$, and $\bar{f}$ the isotropic slowing down tail distribution function (for which $\partial \bar{f}/\partial \varphi = 0 = \partial \bar{f}/\partial \lambda$):

(3.16)\begin{equation}\bar{f} = \frac{{S{\tau _s}H({v_0} - v)}}{{4\mathrm{\pi }({v^3} + v_c^3)}},\end{equation}

where ${\tau _s} = 3{M_\alpha }T_e^{3/2}/4{(2\mathrm{\pi }{m_e})^{1/2}}Z_\alpha ^2{e^4}{n_e}\ell n\varLambda$ is the slowing down time, ${v_c}$ is the critical speed defined by $v_c^3 = 3{\mathrm{\pi }^{1/2}}{(2{m_e})^{ - 1/2}}n_e^{ - 1}T_e^{3/2}{\varSigma _i}(Z_i^2{n_i}/{M_i})$ and H is a step function.

For the narrow boundary layers of interest here only the diffusive terms due to pitch angle scattering of the alphas by the background ions are required (drag is negligible since $v_\lambda ^3/v_0^3 \gg {(\Delta \lambda )^2}\sim {({\varepsilon ^2}\bar{\nu }/{\bar{\omega }_\alpha })^{2/3}}$), giving

(3.17)\begin{equation}\bar{C}\{ \,\tilde{f}\} \approx \frac{{2v_\lambda ^3}}{{{\tau _s}{v^3}\oint_\alpha {\textrm{d}\tau } }}\frac{\partial }{{\partial \lambda }}\left[ {\lambda \left( {\oint_\alpha {\textrm{d}\tau } \frac{\xi }{B}} \right)\frac{{\partial \tilde{f}}}{{\partial \lambda }}} \right] = \frac{{4\varepsilon v_\lambda ^3}}{{{\tau _s}{v^3}K(\kappa )}}\frac{\partial }{{\partial \lambda }}\left\{ {[E(\kappa ) - (1 - {\kappa^2})K(\kappa )]\frac{{\partial \tilde{f}}}{{\partial \lambda }}} \right\} = p\bar{\nu }\varepsilon \frac{{{\partial ^2}\tilde{f}}}{{\partial {\lambda ^2}}}.\end{equation}

Here $v_\lambda ^3 = 3{\mathrm{\pi }^{1/2}}T_e^{3/2}{\varSigma _i}Z_i^2{n_i}/{(2{m_e})^{1/2}}{M_\alpha }{n_e}\sim v_c^3$ and

(3.18)\begin{equation}\bar{\nu } = \frac{{4v_\lambda ^3(\kappa _0^2 - 1/2)}}{{p{\tau _s}{v^3}}}\end{equation}

is defined after approximating the elliptic integrals in the collision operator by their drift reversal values by using $[E(\kappa ) - (1 - {\kappa ^2})K(\kappa )]/K(\kappa ) \to \kappa _0^2 - 1/2$.

The preceding definitions yield the desired nonlinear form of the reduced drift kinetic equation for the trapped alphas to be

(3.19)\begin{equation}{\bar{\omega }_\alpha }\frac{{\partial \tilde{f}}}{{\partial \varphi }} - \frac{{\bar{V}}}{{{R_0}}}\sin \varphi \; \frac{{\partial (\bar{f} + \tilde{f})}}{{\partial \varepsilon }} = \bar{\nu }\varepsilon \frac{{{\partial ^2}\tilde{f}}}{{\partial {\lambda ^2}}}.\end{equation}

When $\partial \tilde{f}/\partial \varepsilon \ll \partial \bar{f}/\partial \varepsilon$ the preceding equation has the usual resonant plateau or superbanana plateau solution (as in quasilinear theory) as noted earlier.

The problem of interest here is the alternative limit in which $\partial \tilde{f}/\partial \varepsilon \sim \partial \bar{f}/\partial \varepsilon = \bar{f}^{\prime}$ with $\bar{\nu }$ very small, but finite, and $\bar{f}^{\prime}$ effectively a constant since it varies slowly with $\varepsilon$ and is independent of $\eta ,\varphi = p\alpha$ and $\lambda$. Before finding a solution with the nonlinearity retained, it is convenient to cast the equation into a convenient form. Dividing by $\mathrm{\pi }q{v^2}/{\varOmega _0}a^{\prime}{R_0}$ leads to

(3.20)\begin{equation}\frac{{[\lambda - (1 + \varepsilon - 2\kappa _0^2\varepsilon )]}}{{8\varepsilon \kappa _0^2(1 - \kappa _0^2)}}\frac{{\partial \tilde{f}}}{{\partial \varphi }} - \varTheta \delta \sin \varphi \; \frac{{\partial (\bar{f} + \tilde{f})}}{{\partial \varepsilon }} = \frac{{\bar{\nu }{R_0}\varepsilon {\varOmega _0}a^{\prime}}}{{\mathrm{\pi }q{v^2}}}\frac{{{\partial ^2}\tilde{f}}}{{\partial {\lambda ^2}}}.\end{equation}

Only the $\varepsilon$ dependence of $\partial \tilde{f}/\partial \varepsilon$ and the resonance of the coefficient of $\partial \tilde{f}/\partial \varphi$ matter. All other $\varepsilon$ dependence is unimportant. In the next section, the preceding equation will be solved in a weak collisionality limit that allows the nonlinear kinetic pod structure to be retained.

4. Nonlinear kinetic equation solution in the presence of a pod

Defining the new radial variable x and pitch angle variable $\varLambda$ by letting

(4.1)\begin{equation}\varepsilon = xL\end{equation}

and

(4.2)\begin{equation}(1 - \lambda )/(2\kappa _0^2 - 1) = \; \varLambda L,\end{equation}

with $\bar{f}^{\prime} = \partial \bar{f}/\partial x = L\partial \bar{f}/\partial \varepsilon$ and

(4.3)\begin{equation}L = {\left[ {\frac{{8\varepsilon \kappa_0^2(1 - \kappa_0^2)\varTheta \delta }}{{2\kappa_0^2 - 1}}} \right]^{1/2}},\end{equation}

leads to the more compact form

(4.4)\begin{equation}(x - \varLambda )\frac{{\partial \tilde{f}}}{{\partial \varphi }} - \sin \varphi \; \left( {\bar{f}^{\prime} + \frac{{\partial \tilde{f}}}{{\partial x}}} \right) = \Delta \frac{{{\partial ^2}\tilde{f}}}{{\partial {\varLambda ^2}}},\end{equation}

where

(4.5)\begin{equation}\varDelta = \frac{{\bar{\nu }\sqrt {2\varepsilon } {R_0}{\varOmega _0}a^{\prime}}}{{4\mathrm{\pi }{{[\kappa _0^2(1 - \kappa _0^2)]}^{1/2}}]{{[(2\kappa _0^2 - 1)\varTheta \delta ]}^{3/2}}q{v^2}}}.\end{equation}

Notice that $\delta \to 0$ is not compatible with $\varDelta \ll 1$.

Ignoring a term that is a constant multiplying $\bar{f}^{\prime}$ for the moment by letting

(4.6)\begin{equation}\tilde{f} = g - (x - \varLambda )\bar{f}^{\prime},\end{equation}

and recalling $\partial \bar{f}^{\prime}/\partial \varphi = 0 = \partial \bar{f}^{\prime}/\partial \varLambda$, then $\bar{f}^{\prime} + \partial \tilde{f}/\partial x = \partial g/\partial x$ gives a convenient form similar to that considered by Hamilton et al. (Reference Hamilton, Tolman, Arzamasskiy and Duarte2023) in an astrophysical context, namely

(4.7)\begin{equation}(x - \varLambda )\frac{{\partial g}}{{\partial \varphi }} - \sin \varphi \; \frac{{\partial g}}{{\partial x}} = \Delta \frac{{{\partial ^2}g}}{{\partial {\varLambda ^2}}}.\end{equation}

In this compact form, the steady state, fully phase mixed solution in the presence of weak collisions $\textrm{(}\varDelta \ll 1)$ found by Hamilton et al. (Reference Hamilton, Tolman, Arzamasskiy and Duarte2023) is of interest. The Hamilton et al. (Reference Hamilton, Tolman, Arzamasskiy and Duarte2023) work was preceded by an incomplete solution by Petviachvili (Reference Petviachvili1999).

To recover the precise Hamilton et al. (Reference Hamilton, Tolman, Arzamasskiy and Duarte2023) form, it is convenient to introduce the reduced Hamiltonian or constant of the motion:

(4.8)\begin{equation}h = {\textstyle{1 \over 2}}{(x - \varLambda )^2} - \cos \varphi ,\end{equation}

with h = 1 the location of the separatrix between the bound $( - 1 < h < 1)$ and unbound $(h > 1)$ motion. The preceding allows the kinetic equation to be written as

(4.9)\begin{equation}\frac{{\partial h}}{{\partial x}}\frac{{\partial g}}{{\partial \phi }} - \; \frac{{\partial h}}{{\partial \phi }}\; \frac{{\partial g}}{{\partial x}} = \Delta \frac{{{\partial ^2}g}}{{\partial {\varLambda ^2}}}.\end{equation}

More interestingly, as x and $\varLambda$ only enter in the combination $x - \varLambda$, it is convenient to define a mixed kinetic variable j depending on minor radius x and pitch angle $\varLambda$ that can also be written in terms of the reduced constant of the motion h and the angular dependent departure from quasisymmetry $\varphi$ as

(4.10)\begin{equation}j(x,\varLambda ) = x - \varLambda = \sigma \sqrt {2(h + \cos \varphi )} = j(h,\varphi ),\end{equation}

with $\sigma ={\pm} 1$. This kinetic variable j simplifies the reduced kinetic equation to the steady state form of Hamilton et al. (Reference Hamilton, Tolman, Arzamasskiy and Duarte2023):

(4.11)\begin{equation}j\frac{{\partial g}}{{\partial \varphi }} - \sin \varphi \; \frac{{\partial g}}{{\partial j}} = \Delta \frac{{{\partial ^2}g}}{{\partial {j^2}}}.\end{equation}

The steady state solution for $\varDelta = 0.001$ shown in their figure 2(a) is reproduced here with their kind permission as figure 1. However, the kinetic variable j used here differs from that of Hamilton et al. as they are not interested in the distinction between x and $\varLambda$. This distinction is important here because the presence of a resonance requires drift reversal and thereby leads to pitch angle $\varLambda$ and radial x variation, with the radial flattening occurring in g and not $\bar{f}$. The distinction also matters when the energy flux is evaluated in the next section. Pods are defined in $j,\varphi$ space with configuration space x and pitch angle $\varLambda$ islands appearing for fixed $\varLambda$ and x, respectively. Discontinuities are resolved by the narrow collisional boundary layers about the $h = 1$ separatrices or pod boundaries.

Figure 1.

Contours of constant $g(j,\varphi )$, with the flattened bound region inside the separatrix at h = 1 and the two unbound regions above (in red and yellow) and below (in dark and light blue). Very narrow collisional boundary layers surround the separatrix. (Reprinted with permission from Hamilton et al. (Reference Hamilton, Tolman, Arzamasskiy and Duarte2023).)

By keeping a time derivative, Hamilton et al. (Reference Hamilton, Tolman, Arzamasskiy and Duarte2023) numerically solved the preceding equation for various values of $\varDelta$ to find skew symmetric steady state solutions satisfying

(4.12)\begin{equation}g(j,\; \varphi ) ={-} g( - j, - \; \varphi )\end{equation}

that relax to the steady state panels shown in their figure 2. The skew symmetric form of the solution means the oscillating function $\varTheta$ can be replaced by $|\varTheta |$ in L and $\varDelta$ since the same solution procedure is valid even if the sign of the $\sin \varphi$ term changes in (4.11).

Continuing to follow Hamilton et al. (Reference Hamilton, Tolman, Arzamasskiy and Duarte2023) by changing from $j,\varphi$ variables to $h,\varphi$ variables using

(4.13)\begin{equation}{\left. {\frac{{\partial g}}{{\partial \varphi }}} \right|_j} = {\left. {\frac{{\partial g}}{{\partial \varphi }}} \right|_h} + {\left. {\frac{{\partial h}}{{\partial \varphi }}} \right|_j}{\left. {\frac{{\partial g}}{{\partial h}}} \right|_\varphi } = {\left. {\frac{{\partial g}}{{\partial \varphi }}} \right|_h} + \sin \varphi {\left. {\frac{{\partial g}}{{\partial h}}} \right|_\varphi }\end{equation}

and

(4.14)\begin{equation}{\left. {\frac{{\partial g}}{{\partial j}}} \right|_\varphi } = {\left. {\frac{{\partial h}}{{\partial j}}} \right|_\varphi }{\left. {\frac{{\partial g}}{{\partial h}}} \right|_\varphi } = j{\left. {\frac{{\partial g}}{{\partial h}}} \right|_\varphi },\end{equation}

the kinetic equation simplifies further to become

(4.15)\begin{equation}j{\left. {\frac{{\partial g}}{{\partial \varphi }}} \right|_h} = \Delta j{\left. {\frac{\partial }{{\partial h}}} \right|_\varphi }\left( {j{{\left. {\frac{{\partial g}}{{\partial h}}} \right|}_\varphi }} \right).\end{equation}

In the $\varDelta \ll 1$ limit, the lowest order motion is collisionless, but the solution form is collisionally constrained. The solution can only depend on h for the isotropically born alphas as they are confined in a nearly QS magnetic field. Furthermore, to cancel the $- (x - \varLambda )\bar{f}^{\prime} ={-} j\bar{f}^{\prime}$ term in $\tilde{f}$ in the freely circulating limit (that is, at large pitch angle for each flux surface) far from the separatrix, a skew symmetric solution is required to lowest order, as in Hamilton et al. (Reference Hamilton, Tolman, Arzamasskiy and Duarte2023). Moreover, based on the numerical solution of Hamilton et al. (Reference Hamilton, Tolman, Arzamasskiy and Duarte2023), any fine structure associated with the bound or librating motion is collisionally phase mixed away to find the flattened solution $g = 0$ since it cannot depend on $\varphi$ and must vanish to satisfy skew symmetry. This solution will be verified shortly. Consequently, a solution of the form

(4.16)\begin{equation}g = {g_1}(h) + {g_2}(h,\varphi ) + \cdots \end{equation}

is assumed to satisfy

(4.17)\begin{equation}{\left. {\frac{{\partial {g_1}}}{{\partial \varphi }}} \right|_h} = 0,\end{equation}

to lowest order. To next order, the equation is simply

(4.18)\begin{equation}j{\left. {\frac{{\partial {g_2}}}{{\partial \varphi }}} \right|_h} = \Delta j{\left. {\frac{\partial }{{\partial h}}} \right|_\varphi }\left( {j{{\left. {\frac{{\partial {g_1}}}{{\partial h}}} \right|}_\varphi }} \right).\end{equation}

Solutions must be found for both the bound (or librating) alpha motion in the pods or tubes and the unbound (or circulating) alpha particle motion outside them. The outer limit of the unbound solution far from the separatrix must smoothly match onto the inner limit of the solution far from the separatrix defining the pod boundary.

Dividing by j and eliminating ${g_2}$ by integration over a full bound periodic orbit or from $- \mathrm{\pi }$ to $\mathrm{\pi }$ for an unbound (at fixed h) orbit, leads to the solubility constraint

(4.19)\begin{equation}{\left. {\frac{\partial }{{\partial h}}} \right|_\varphi }\left[ {({\oint_h {\textrm{d}\varphi j} } ){{\left. {\frac{{\partial {g_1}}}{{\partial h}}} \right|}_\varphi }} \right] = 0.\end{equation}

For the unbound using $j = \sigma \sqrt {2(h + 1) - 4{{\sin }^2}(\phi /2)}$ with $k = \sqrt {2/(h + 1)}$:

(4.20)\begin{equation}\int_{ - \mathrm{\pi }}^\mathrm{\pi } {\textrm{d}\varphi j} = \sigma 4\sqrt {2(h + 1)} \int_0^{\mathrm{\pi }/2} {\textrm{d}t\sqrt {1 - [2/(h + 1)]{{\sin }^2}(t)} } = \sigma 8{k^{ - 1}}E(k),\end{equation}

where k = 0 (or $h \to \infty$) are freely circulating and the separatrix is at k = 1 (or $h = 1$). Far from the separatrix $({j^2} \gg 2\cos \varphi )$ the alpha distribution function must go to the slowing down tail result $\bar{f}$. Consequently, ${g_1} \to (x - \varLambda )\bar{f}^{\prime} = j\bar{f}^{\prime} \approx \sigma \bar{f}^{\prime}\sqrt {2h}$ is required. Therefore, to match the freely circulating limit requires $\partial {g_1}/\partial h{|_\varphi } \to \sigma \bar{f}^{\prime}/\sqrt {2h} \approx \sigma k\bar{f}^{\prime}/2$, and gives the unbound or circulating solution

(4.21)\begin{equation}{\left. {\frac{{\partial {g_1}}}{{\partial h}}} \right|_\varphi } = \frac{{\sigma \mathrm{\pi }k\bar{f}^{\prime}}}{{4E(k)}}.\end{equation}

This form is equivalent to (A10) of Hamilton et al. (Reference Hamilton, Tolman, Arzamasskiy and Duarte2023). Integrating from the separatrix, where ${g_1}(h = 1) = 0$ to be skew symmetric, the piecewise continuous solution is (Hamilton Reference Hamilton2024)

(4.22)\begin{equation}{g_1} = \frac{{\mathrm{\pi }\sigma \bar{f}^{\prime}}}{4}\int_1^h {\frac{{\textrm{d}hk}}{{E(k)}}} = \mathrm{\pi }\sigma \bar{f}^{\prime}\int_k^1 {\frac{{\textrm{d}t}}{{{t^2}E(t)}}} ,\end{equation}

with $k = \sqrt {2/(h + 1)}$. As $k \to 0$,

(4.23)\begin{equation}\frac{\mathrm{\pi }}{2}\int_k^1 {\frac{{\textrm{d}t}}{{{t^2}E(t)}}} = \int_k^1 {\frac{{\textrm{d}t}}{{{t^2}}}} + \int_k^1 {\frac{{\textrm{d}t}}{{{t^2}}}\left[ {\frac{\mathrm{\pi }}{{2E(t)}} - 1} \right]} \approx \frac{1}{k} - 0.6894 - \frac{k}{4} + O({k^3}),\end{equation}

giving the desired freely circulating result

(4.24)\begin{equation}{g_1} \to \sigma \bar{f}^{\prime}(\sqrt {2h} - 1.379 + \cdots ) \approx (j - 1.379\sigma )\bar{f}^{\prime}.\end{equation}

Consequently, the lowest order unbound solution is

(4.25)\begin{equation}\tilde{f} = {g_1} - (j - 1.379\sigma )\bar{f}^{\prime} = \bar{f}^{\prime}\left[ {\sigma \mathrm{\pi }\int_k^1 {\frac{{\textrm{d}t}}{{{t^2}E(t)}}} - j + 1.379\sigma } \right].\end{equation}

To account for the constant from the integral, $\tilde{f} = g - (x - \varLambda - 1.379\sigma )\bar{f}^{\prime}$ generalizes (4.6) to get the solution for $\tilde{f}$. The skew symmetric factor $1.379\sigma \bar{f}^{\prime}$ is needed to properly account for the three phase space regions (two unbound with $\sigma ={\pm} 1$ and one bound with $\sigma = 0$) separated by collisional boundary layers at the separatrices as in figure 1.

To demonstrate that the bound solution ${g_1} = 0$ is valid, let $b\sin t = \sin (\varphi /2)$. Then

(4.26)\begin{equation}\int_h {\textrm{d}\varphi j/8(1 + h)} = \int_0^{\mathrm{\pi }/2} {\textrm{d}t{{(1 - {b^2}{{\sin }^2}t)}^{ - 1/2}}{{\cos }^2}t} = (1 - {b^{ - 2}})K(b) + {b^{ - 2}}E(b),\end{equation}

with ${b^2} = (1 + h)/2$. Integrating the bound solubility constraint up to the separatrix gives

(4.27)\begin{equation}[E(b) - (1 - {b^2})K(b )]{\left. {\frac{{\partial {g_1}}}{{\partial h}}} \right|_\varphi } = {\left. {\frac{{\partial {g_1}}}{{\partial h}}} \right|_{\varphi ,h = 1}}.\end{equation}

As $b \to 0$, $E(b) - (1 - {b^2})K(b) \to \mathrm{\pi }{b^2}/4$, which would require $\partial {g_1}/\partial h{|_\varphi } \to \infty$ to keep $\partial {g_1}/\partial h{|_{\varphi ,h = 1}}$ finite. Therefore, ${g_1} = 0$ is the only well behaved, skew-symmetric, fully phase mixed, bound solution in the presence of weak collisions that is consistent with the solubility constraint from the collision operator as shown numerically by Hamilton et al. (Reference Hamilton, Tolman, Arzamasskiy and Duarte2023). Consequently, the desired piecewise continuous skew symmetric solution for the bound is

(4.28)\begin{equation}\tilde{f} ={-} j\bar{f}^{\prime}.\end{equation}

The preceding unbound and bound solutions for g are in agreement with figure 2 of Hamilton et al. (Reference Hamilton, Tolman, Arzamasskiy and Duarte2023) for $\varDelta = 0.001$ in the steady state. The $1.379\sigma$ step in the solution ${f_1}$ at the separatrix, as seen by comparing (4.25) and (4.28), is smoothed by a narrow collisional boundary layer at the separatrix that need not be resolved by the solution procedure here. The behaviour in the boundary layer does not play a role in the results that follow. The piecewise continuous behaviour is removed in the numerical solution if boundary layers are resolved. In the next section these analytic solutions are used to form the alpha energy flux.

5. Energy flux

The radial energy flux at large aspect ratio and fixed x is

(5.1)\begin{equation}Q = \frac{{\mathrm{\pi }q{M_\alpha }}}{{{B_0}a^{\prime}}}\left\langle {\int {{\textrm{d}^3}vf{v^2}{\boldsymbol{v}_d}\boldsymbol{\cdot }\boldsymbol{\nabla }{\psi_p}} } \right\rangle .\end{equation}

It must be evaluated to determine the alpha particle loss due to the departure from quasisymmetry, where the flux surface average is

(5.2)\begin{equation}\langle \ldots \rangle = {{\oint {\frac{{\textrm{d}\vartheta \,\textrm{d}\zeta }}{{\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla }\vartheta }}( \ldots )} } / {\oint {\frac{{\textrm{d}\vartheta \,\textrm{d}\zeta }}{{\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla }\vartheta }}} }} = {{\oint {\frac{{\textrm{d}\eta \,\textrm{d}\varphi }}{{\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla }\eta }}( \ldots )} } / {\oint {\frac{{\textrm{d}\eta \,\textrm{d}\varphi }}{{\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla }\eta }}} }},\end{equation}

and ${\textrm{d}^3}v = 2\mathrm{\pi }{v^3}\,\textrm{d}v\,\textrm{d}\lambda B/|{v_{||}}|{B_0}$, where both signs of ${v_{||}}$ are summed over by $\int_\alpha {\textrm{d}\tau }$ in the numerator as only trapped alphas contribute.

Keeping only the radial drift due to a departure from quasisymmetry, recalling $\partial f/\partial \eta = 0$, and using the results from § 3 yields

(5.3)\begin{equation}\frac{{2\mathrm{\pi }q}}{{{B_0}a^{\prime}}}\oint_\alpha {\textrm{d}\tau {\boldsymbol{v}_d}\boldsymbol{\cdot }\boldsymbol{\nabla }{\psi _p}} ={-} p\bar{V}\sin \varphi \oint_\alpha {\textrm{d}\tau } .\end{equation}

As a result. the energy flux becomes

(5.4)\begin{align}Q = \frac{{2{\mathrm{\pi }^2}q{M_\alpha }\int {\textrm{d}v} \int {\textrm{d}\lambda {v^5}} \oint {\textrm{d}\varphi \tilde{f}} \oint_\alpha {\textrm{d}\tau {\boldsymbol{v}_d}\boldsymbol{\cdot }\boldsymbol{\nabla }{\psi _p}} }}{{B_o^2a^{\prime}\oint {\textrm{d}\eta \,\textrm{d}\varphi /\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla }\eta } }} ={-} \frac{{\mathrm{\pi }p{M_\alpha }\int {\textrm{d}v{v^5}} \int {\textrm{d}\lambda \bar{V}} \oint {\textrm{d}\varphi g\sin \varphi } \oint_\alpha {\textrm{d}\tau } }}{{{B_0}\oint {\textrm{d}\eta \,\textrm{d}\varphi /\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla }\eta } }},\end{align}

where only g from $\tilde{f} = g - (x - \varLambda - 1.379\sigma )\bar{f}^{\prime}$ contributes because of the $\varphi$ integral.

The reduced constant of the motion h means that collisionless orbits are confined and no collisional transport can occur due to ${g_1}(h)$. Recalling $\textrm{d}\lambda = (2\kappa _0^2 - 1)L\,\textrm{d}h/j$, accounting for both $\sigma$ signs by multiplying by 2 when integrating over h, and using (4.21)–(4.25), gives

(5.5)\begin{align}\begin{aligned} \int_{ - \infty }^\infty {\textrm{d}j} \oint {\dfrac{{\textrm{d}\varphi g\sin \varphi }}{{\oint {\textrm{d}\varphi } }}} & = 2\int_1^\infty {\textrm{d}h} \oint {\dfrac{{\textrm{d}\varphi g\sin \varphi }}{{j\oint {\textrm{d}\varphi } }}} ={-} 2\int_1^\infty {\textrm{d}h} \oint {\textrm{d}\varphi } \dfrac{g}{{\oint {\textrm{d}\varphi } }}{\left. {\dfrac{{\partial j}}{{\partial \varphi }}} \right|_h}\\ & = 2\int_1^\infty {\textrm{d}h} \oint {\textrm{d}\varphi } \dfrac{j}{{\oint {\textrm{d}\varphi } }}{\left. {\dfrac{{\partial {g_2}}}{{\partial \varphi }}} \right|_h} = \dfrac{{2\Delta }}{{\oint {\textrm{d}\varphi } }}\oint {\textrm{d}\varphi } \int_1^\infty {\textrm{d}hj} {\left. {\dfrac{\partial }{{\partial h}}} \right|_\varphi }\left( {j{{\left. {\dfrac{{\partial {g_1}}}{{\partial h}}} \right|}_\varphi }} \right)\\ & = 2\Delta \left[ {\left( {\dfrac{{\oint {\textrm{d}\varphi {j^2}} }}{{\oint {\textrm{d}\varphi } }}} \right)\left. {\dfrac{{\partial {g_1}}}{{\partial h}}} \right|_{\varphi ,h = 1}^{h \to \infty } - {{g_1}} |_{h = 1}^{h \to \infty }} \right] = 2\Delta \bar{f}^{\prime}\left[ {\left. {\dfrac{{\mathrm{\pi }\sqrt 2 h}}{{2E(k)\sqrt {h + 1} }}} \right|_{\varphi ,h = 1}^{h \to \infty } - (\sqrt {2h} - 1.379)} \right]\\ & ={-} (\mathrm{\pi } - 2.758)\Delta \bar{f}^{\prime} ={-} 0.122\mathrm{\pi }\Delta \bar{f}^{\prime}. \end{aligned}\end{align}

The lower limit on the h integral is a reminder that only the unbound $(h \ge 1)$ contribute to the radial transport. The unbound require $1 \le h = {j^2}/2 - \cos \varphi$ indicating ${j^2} \ge 2(1 + \cos \varphi ) > 0$, with $- \mathrm{\pi } < \varphi < \mathrm{\pi }$, so no singularity occurs at $j = 0$ as the unbound never make it to j = 0 (the centreline of the phase space pod).

Inserting (5.5) in energy flux and noticing $\bar{f}^{\prime} = L\partial \bar{f}/\partial \varepsilon$, Q becomes

(5.6)\begin{equation}Q ={-} \frac{{0.122{\mathrm{\pi }^2}p(2\kappa _0^2 - 1){M_\alpha }({v\oint_\alpha {\textrm{d}\tau } } )\int_0^{{v_0}} {\textrm{d}v{v^4}{L^2}\bar{V}\Delta \partial \bar{f}} /\partial \varepsilon }}{{{B_0}\oint {\textrm{d}\eta /\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla }\eta } }},\end{equation}

with $({v\oint_\alpha {\textrm{d}\tau } } )/{B_0}\oint {\textrm{d}\eta /\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla }\eta } \approx 4K({\kappa _0})/\mathrm{\pi }\sqrt {2\varepsilon }$, $L\sim {(\varepsilon \varTheta \delta )^{1/2}}$ and

(5.7)\begin{align}\begin{aligned} (2\kappa _0^2 - 1){L^2}\bar{V}\varDelta & = \dfrac{{{{[\kappa _0^2(1 - \kappa _0^2)]}^{1/2}}}}{{{{(2\kappa _0^2 - 1)}^{3/2}}}}{(2\varepsilon )^{3/2}}\bar{\nu }{R_0}{(|\varTheta |\delta )^{1/2}}\\ & = \dfrac{{2{{[\kappa _0^2(1 - \kappa _0^2)]}^{1/2}}}}{{p{{(2\kappa _0^2 - 1)}^{1/2}}}}\dfrac{{{{(2\varepsilon )}^{3/2}}{R_0}v_\lambda ^3}}{{{\tau _s}{v^3}}}{(|\varTheta |\delta )^{1/2}}. \end{aligned}\end{align}

Using the slowing down tail form for $\bar{f}$ and assuming ${v_0}/{v_c} \gg 1$ gives a result dominated by the contributions from the $v\sim {v_c}$ alphas, namely

(5.8)\begin{equation}\int_0^{{v_0}} {\textrm{d}vv} \frac{{\partial \bar{f}}}{{\partial r}} = \frac{1}{{4\mathrm{\pi }}}\frac{\partial }{{\partial r}}\left( {\frac{{S{\tau_s}}}{{{v_c}}}\int_0^{{v_0}/{v_c}} {\frac{{\textrm{d}\kern0.7pt xx}}{{{x^3} + 1}}} } \right) \approx \frac{\partial }{{\partial r}}\left( {\frac{{S{\tau_s}}}{{6\sqrt 3 {v_c}}}} \right).\end{equation}

As a result, the energy flux of the alphas is

(5.9)\begin{align}Q ={-} \frac{{0.976K({\kappa _0}){{[\kappa _0^2(1 - \kappa _0^2)]}^{1/2}}}}{{3\sqrt 3 {{(2\kappa _0^2 - 1)}^{1/2}}}}\frac{{\int {R_0^2{M_\alpha }v_\lambda ^3{{(|\varTheta |\delta )}^{1/2}}} }}{{{\tau _s}}}\frac{\partial }{{\partial r}}\left( {\frac{{S{\tau_s}}}{{{v_c}}}} \right) ={-} 0.203\frac{{\int {R_0^2{M_\alpha }v_\lambda ^3{{(|\varTheta |\delta )}^{1/2}}} }}{{{\tau _s}}}\frac{\partial }{{\partial r}}\left( {\frac{{S{\tau_s}}}{{{v_c}}}} \right),\end{align}

provided $\varDelta \ll 1$ or ${\delta ^{3/2}}qv/\sqrt {2\varepsilon } {\varOmega _0}a^{\prime} \gg \bar{\nu }{R_0}/v$. The ${\delta ^{1/2}}$ dependence of Q implies the energy flux is proportional to the radial pod width. In this limit, the normalized pod width is finite and satisfies ${(\varepsilon \delta )^{1/2}} \gg {(\bar{\nu }{R_0}/v)^{2/3}}{(\varepsilon {\varOmega _0}a^{\prime}/qv)^{2/3}}\sim {({R_0}/{v_c})^{2/3}}{(\varepsilon {\varOmega _0}a^{\prime}/q{v_\lambda })^{2/3}}$ for speeds as small as ${v_\lambda }\sim {v_c}$. At these speeds, removing a birth energy factor ${M_\alpha }v_0^2/2$ and assuming $\varTheta \approx 1$ to obtain an upper bound estimate of the transport, the diffusivity is of order

(5.10)\begin{equation}{D_\nu }\sim \varepsilon R_0^2v_\lambda ^3{\delta ^{1/2}}/{v_c}v_0^2{\tau _s},\end{equation}

as in (2.13), where the alpha density is ${n_\alpha } \approx S{\tau _s}/\ell n({v_0}/{v_c})\sim S{\tau _s}$. (When $|\varTheta |\ll 1\;$ such that $\varDelta \gg 1$, only resonant plateau transport occurs.) As indicated in § 2, this diffusivity can be comparable to or less than the diffusivity associated with resonant plateau transport. But unlike resonant plateau transport, the transport in the presence of the phase space pods acts to remove the alphas with speeds of the order of ${v_\lambda }\sim {v_c}$, rather than those just being born. Consequently, it might prove useful for ash removal, since the particle diffusivity will be of the order of $v_0^2/v_c^2$ larger than (5.10).

Large diffusive alpha energy loss in a slowing down time is avoided if ${\tau _s}{D_\nu }/{r^2} \ll 1$ or $\delta /\varepsilon \ll \varepsilon v_0^4v_c^2/v_\lambda ^6\sim 1\unicode{x2013}10$. Moreover, the normalized radial phase space pod width ${(\varepsilon \delta )^{1/2}}$ must be larger than the collisional boundary layer width ${({\varepsilon ^2}\bar{\nu }/{\bar{\omega }_\alpha })^{1/3}}$ to keep $\varDelta \ll 1$ requiring $\delta /\varepsilon \gg {(\bar{\nu }/\varepsilon {\bar{\omega }_\alpha })^{2/3}}\sim {10^{ - 2}}$ for birth alphas with ${R_0}/{v_0}{\tau _s}\sim {10^{ - 5}}$ and ${\varOmega _0}{R_0}/q{v_0}\sim {10^2}$. Ultimately, a comprehensive investigation will be required. However, it is encouraging that resonant plateau transport overestimates the transport level when the optimization is less than perfect as it is ${(\varepsilon \delta )^{3/2}}/({\varepsilon ^2}\bar{\nu }/{\bar{\omega }_\alpha })$ larger. Indeed, the resonant plateau estimate of (2.9) suggests that $\delta /\varepsilon \sim {10^{ - 2}}$ is required for good collisional confinement of birth alphas.