Collisional alpha particle transport in a quasisymmetric stellarator with a single helicity imperfection

Abstract

The tangential drift of the trapped alpha particles in bounce or transit averaged kinetic treatments of stellarators reverses direction on each flux surface at a particular value of pitch angle. The vanishing of the tangential drift corresponds to a resonance that allows a narrow collisional boundary layer to form due to the presence of pitch angle scattering by the background ions. The alphas in and adjacent to this drift reversal layer are particularly sensitive to collisions because they are in or very close to resonance. As a result, enhanced collisional transport occurs due to the existence of this drift reversal resonance in a nearly quasisymmetric stellarator with a single helicity imperfection. Moreover, the value of the resonant pitch angle for drift reversal on neighbouring flux surfaces varies continuously, with the inner flux surfaces having a larger resonant pitch angle than the outer ones. This pitch angle dependence means phase space ‘tubes’ or ‘pods’ exist that connect the inner flux surfaces to the outer ones. These pods allow collisional radial transport of the alphas to extend over the entire radial cross section. When collisions are finite, but weak, and the single helicity departure from quasisymmetry large enough, the collisionless alpha particle motion remains constrained by collisions as they complete their drift trajectories in phase. In particular, the small radial scales introduced by the radial extent or width of the phase space pods require the retention of the nonlinear radial drift term in the kinetic equation. The associated collisional radial transport is evaluated and found to be significant, but is shown to preferentially remove slower speed alphas without substantially affecting birth alphas.

1. Introduction

Resonant plateau transport (Park, Boozer & Menard 2009) occurs whenever there is a resonance in the particle motion or singularity that must be resolved by collisions, with plateau added since it leads to transport independent of collision frequency. The presence of a collisional resonance normally implies the singular behaviour is removed by a boundary layer due to the diffusive nature of collisions (Calvo et al. 2017; Catto 2019; Catto & Tolman 2021; Catto, Tolman & Parra 2023). In a nearly quasisymmetric (QS) stellarator these references consider the resonance that occurs when the transit averaged tangential drift in the flux surface vanishes because the drift reverses. The tangential drift in a flux surface $\psi$ traced out by the magnetic field $\boldsymbol{B} = \boldsymbol{\nabla }\alpha \times \boldsymbol{\nabla }\psi$ is in the $\boldsymbol{\nabla }\alpha$ direction. In a small inverse aspect ratio $\textrm{(}\varepsilon \ll \textrm{1)}$ stellarator the location of this drift resonance depends on two phase space variables: the poloidal flux $\psi \propto \varepsilon$ and an adiabatic invariant pitch angle variable $\lambda$ in velocity space. A resonant plateau treatment collisionally resolves the resonance by solving a linearized kinetic equation. The resulting solution introduces small scale radial variation, $\Delta \varepsilon$, associated with the resonant interaction in the presence of a drive term. As a result, the neglect of the radial derivative of the perturbed distribution function compared with that of the unperturbed distribution limits the validity of such resonant plateau treatments once the phase space structure becomes larger than the collisional boundary layer. This behaviour suggests the existence of a transport regime in which the small radial scales introduced by pitch angle scattering collisions can give rise to well defined phase space regions that have an island structure when projected onto minor radius for a fixed resonant pitch angle and when projected onto pitch angle for a fixed minor radius. This double island phase space structure in normalized minor radius $(\varepsilon )$ and pitch angle $(\lambda )$ variables is tubular or pod shaped. Once this nonlinear behaviour has to be retained the collisional transport is reduced below the predicted resonant plateau level.

This seldom studied phase space pod, weak collisionality limit is the focus of the material presented, and is found to lead to radial transport of resonant alpha particles that is proportional to the collision frequency and the island width associated with a flux surface at fixed pitch angle. It is distinct from the weak collisional regime investigated by d'Herbemont et al. (2022) who considered non-resonant background ion transport in the more complicated magnetic wells associated with larger departures from quasisymmetry. The improved estimate they arrived at was also linear in the collision frequency, $\bar{\nu }$, as found by Mynick (1983) at low collisionality for non-resonant background ions in non-optimized stellarators. Unlike these previous evaluations, the electric field drift plays no role for the alphas as magnetic drift dominates. Interestingly the diffusivity found here is proportional to the island width as in Mynick (1983), who used a radial averaging procedure as in the bumpy torus evaluation of Hazeltine & Catto (1981). This same bumpy torus reference also introduced the truncated Taylor expansion method used by Shaing (2015) and Shaing & Hsu (2014) to obtain superbanana results similar to those of Mynick (1983). Unlike these earlier treatments, the one presented here does not require introducing the second adiabatic invariant ($J = \oint {\textrm{d}\ell {v_{||}}}$, with ${v_{||}}$ the parallel velocity and $\textrm{d}\ell$ the incremental length along $\boldsymbol{B}$). Instead, a reduced constant of the motion is employed as in Hamilton et al. (2023) to avoid the questionable Taylor expansion procedure. To avoid confusion, the imprecise terminology superbanana transport is avoided.

In the sections that follow, the linear in $\bar{\nu }$ regime is treated in a systematic and quantitative manner for trapped alpha particle energy transport in an imperfectly optimized, but nearly QS stellarator. The goal is to determine the behaviour of the radial energy transport when a single helicity departure from quasisymmetry is large enough that a quasilinear treatment fails. Before doing so, some background estimates are given in § 2. Then, the reduced kinetic equation for trapped alpha particle transport in a nearly QS stellarator is reduced to its fundamental form in § 3 where the nonlinear radial drift term is retained to allow the departure from quasisymmetry to form phase space ‘pods’ or ‘tubes’ in minor radius and pitch angle. In § 4 a detailed solution to the nonlinear kinetic equation is obtained in a limit in which the QS alpha particle motion is altered by the phase space pods and weakly disrupted by collisions. Then, in § 5 the energy flux is evaluated to confirm that the diffusivity is proportional to collision frequency and radial island width at fixed pitch angle, due to the departure of the magnetic field from quasisymmetry. A brief discussion follows in § 6.

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 1979) 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 1979) 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. 2022; White 2022).

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 1979; Shaing 2015; Catto 2019). 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. (2023) (see also Catto 2019) that assumes the QS field depends on the angle variable $\eta = M\vartheta - N\zeta$, with $\vartheta$ and $\zeta$ the Boozer (1981) poloidal and toroidal angle variables, respectively. Normally $M = 1$ (Cary & Shasharina 1997; Landreman & Paul 2022), 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 1981). 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 2012). 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 2022) 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 1981; Shaing & Hsu 2014; Shaing 2015) 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. (2023) 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. (2023) is of interest. The Hamilton et al. (2023) work was preceded by an incomplete solution by Petviachvili (1999).

To recover the precise Hamilton et al. (2023) 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. (2023):

(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. (2023).)

By keeping a time derivative, Hamilton et al. (2023) 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. (2023) 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. (2023). Moreover, based on the numerical solution of Hamilton et al. (2023), 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. (2023). Integrating from the separatrix, where ${g_1}(h = 1) = 0$ to be skew symmetric, the piecewise continuous solution is (Hamilton 2024)

(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. (2023). 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. (2023) 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.

6. Discussion

In an imperfect QS stellarator the departure from quasisymmetry can be large enough that the usual resonant plateau or superbanana plateau treatment of the collisional transport of alpha particles fails because the presence of narrow phase space pods introduces small radial scale lengths that require retaining nonlinear effects. When the normalized error magnetic field $\delta$ is large enough and the alpha collisions weak enough to satisfy $\varDelta \ll 1$ for ${v_\lambda }\sim {v_c}\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }v \le {v_0}$, that is, when

(6.1)\begin{equation}1 \gg \frac{{{B_{||}}}}{{\varepsilon {B_0}}} = \frac{\delta }{\varepsilon } \gg {\left( {\frac{{{R_0}{\Omega_0}a^{\prime}}}{{qv_c^2{\tau_s}\varepsilon }}} \right)^{2/3}},\end{equation}

the transport enters a regime not previously considered for alpha particle transport in which the pod structure plays a key role. This pod structure is kinetic in nature as the transport depends on both the minor radius and the pitch angle velocity space variable, with each flux surface having a slightly different resonant pitch angle defining the location of its pod centre. In the weakly collisional regime considered here, the lowest order motion of the alphas is altered by the phase space pods. The transit or bounce averaged kinetic calculation is performed by ignoring collisions in lowest order to find a next order constraint that flattens the perturbed alpha distribution function in the pod for the bound alphas. The barely unbound or circulating alphas are also pod perturbed. The collisional transport is evaluated by making use of the next order reduced kinetic equation that retains collisions. Only the unbound contribute to transport.

The single helicity diffusivity suggests that the collisional resonant plateau is of most concern for birth alphas in imperfect QS stellarators, as pods reduce collisional transport.

The kinetic procedure illustrated here for alphas is also expected to be relevant in other problems, including large amplitude toroidal Alfvén-eigenmode and neoclassical tearing mode driven alpha transport, and intense radio frequency heating and current drive. All that is required for similar physics to hold is the presence of a wave–particle resonance that is sensitive to diffusive collisions.