2.1. Equilibrium and phase-space coordinates

The coordinates describing the tokamak equilibrium are $\psi$ (the poloidal flux function), $\vartheta$ (the poloidal angle)Footnote ² and $\alpha$, defined by

(2.1)\begin{equation} \alpha \equiv \zeta - q\left(\psi\right) \vartheta. \end{equation}

Here, $\zeta$ is the toroidal angle. The safety factor is $q$, which characterizes the twist of the background magnetic field by $q (\psi ) \equiv \hat {\boldsymbol{b}} \boldsymbol {\cdot } \boldsymbol {\nabla } \zeta / ( \hat {\boldsymbol{b}} \boldsymbol {\cdot } \boldsymbol {\nabla } \vartheta )$ with $\hat {\boldsymbol{b}}$ giving the direction of the equilibrium magnetic field.Footnote ³ We study the behaviour of alpha particles confined by a magnetic field which is axisymmetric and stationary except for the perturbation. Such a field is stated to zeroth order in the perturbation via the Clebsch representation as

(2.2)\begin{equation} \boldsymbol{B} = \boldsymbol{\nabla} \alpha \times \boldsymbol{\nabla} \psi =I \left( \psi \right) \boldsymbol{\nabla} \zeta + \boldsymbol{\nabla} \zeta \times \boldsymbol{\nabla} \psi, \end{equation}

with $I(\psi )$ characterizing the strength of the toroidal magnetic field by $B_\zeta = I(\psi )/R$, where $R$ is the major radius coordinate. We do not consider background electric fields because they do not affect alpha particle trajectories as strongly as do the magnetic fields.Footnote ⁴ As noted previously, the unit vector corresponding to this field is

(2.3)\begin{equation} \frac{\boldsymbol{B}}{B} \equiv \hat{\boldsymbol{b}}. \end{equation}

The poloidal component of the field is denoted $B_p$, and can be found from $B_p \approx \epsilon B/q$, with the inverse aspect ratio $\epsilon \approx r /R$, where $r$ is local the minor radius. (The flux coordinate and the minor radius are related by $\partial /\partial \psi = (1/ RB_p ) \partial /\partial r$.) The shear of the field is given by

(2.4)\begin{equation} s\equiv \left(r/q\right) \partial q /\partial r. \end{equation}

The total magnetic and electric fields affecting the alpha particles will include both the equilibrium magnetic field and magnetic and electric fields resulting from the perturbations, which are discussed in § 3. The total magnetic field from all of these sources is denoted $\boldsymbol {B}_{\textrm {tot}}$ (with unit vector $\hat {\boldsymbol{b}}_{\textrm {tot}}$) and the total electric field is $\boldsymbol {E}_{\textrm {tot}}$. Several parameters can be defined in terms of the total field or the unperturbed field; the total quantities are in general only used before the perturbation analysis is carried out.

Alpha particles in these fields are characterized by their speed $v$ (or, equivalently, energy normalized to mass, $\mathcal {E} \equiv v^2/2$), the sign of $v_\parallel$, their velocity parallel to the equilibrium magnetic field, and their pitch angle, which may be defined relative to the total magnetic field,

(2.5)\begin{equation} \lambda_{\textrm{tot}} = \frac{B_{0} v_{{\perp},\textrm{tot}}^2}{B_{\textrm{tot}} v^2}, \end{equation}

with $B_0$ the equilibrium on-axis magnetic field strength, or relative to the unperturbed field,

(2.6)\begin{equation} \lambda= \frac{B_{0} v_\perp^2}{B v^2}. \end{equation}

(Here, $v_{\perp ,\textrm {tot}}$ is the speed perpendicular to total field and $v_\perp$ is the speed perpendicular to the unperturbed field.)

The alpha particle mass and charge are given by $M_\alpha$ and $Z_\alpha$, respectively; $\varOmega _{\textrm {tot}} \equiv Z_\alpha eB_{\textrm {tot}}/M_\alpha c$ is the gyrofrequency in the total field and $\varOmega \equiv Z_\alpha eB/M_\alpha c$ the gyrofrequency in the unperturbed field. The alpha particle poloidal gyrofrequency is $\varOmega _p \equiv Z_\alpha eB_p/M_\alpha c$. The alpha particle gyroradius is $\rho _\alpha \equiv v_\perp / \varOmega$, and the poloidal gyroradius $\rho _{p \alpha } \equiv v_\perp / \varOmega _p$.

2.2. Drift kinetic equation and equilibrium alpha particle distribution

We study the effect of perturbations on the energetic alpha particle distribution function $f$ using Hazeltine's drift kinetic equation (Hazeltine Reference Hazeltine1973). This formalism is appropriate when $\rho _\alpha$ (which at the alpha particle birth speed is typically a couple of centimetres) is much less than the length scales relevant to the problem, including the perpendicular scale length of the perturbation and the scale length of the alpha particle density,

(2.7)\begin{equation} a_\alpha \equiv{-} \frac{n_\alpha}{R B_p \partial n_\alpha /\partial \psi}, \end{equation}

with $n_\alpha$ the alpha particle density. The drift kinetic equation reads

(2.8)\begin{align} &\frac{\partial f}{\partial t} + (v_{{\parallel} , \textrm{tot}} \hat {\boldsymbol{b}}_{\textrm{tot}} + \boldsymbol{v}_{d,\textrm{tot}}) \boldsymbol{\cdot} \boldsymbol{\nabla} f + \left[\frac{Z_\alpha e}{M_\alpha } (v_{{\parallel},\textrm{tot}} \hat {\boldsymbol{b}}_{\textrm{tot}} + \boldsymbol{v}_{d,\textrm{tot}})\boldsymbol{\cdot} \boldsymbol{E}_{\textrm{tot}} + \frac{\lambda_{\textrm{tot}} v^2}{2 B_0} \frac{\partial B_{\textrm{tot}}}{\partial t} \right] \frac{\partial f}{\partial \mathcal{E}}\nonumber\\ &\quad = C\left\{f\right\} +\frac{S_{\textrm{fus}} \delta \left(v-v_0\right)}{4 {\rm \pi} v^2}. \end{align}

Here, $S_{\textrm {fus}}$ is the source rate due to fusion, $v_0$ is the alpha particle birth speed (${1.3\times 10^{9}}$ cm s$^{-1}$) and $\boldsymbol {E}_{\textrm {tot}}$ is the total electric field (equal to the contribution from the perturbations, since we neglect the background electric field). The velocity parallel to the total field is $v_{\parallel ,\textrm {tot}}$. The energetic alpha collision operator represents collisions with the background plasma electrons and ions, and is given by (Cordey Reference Cordey1976; Catto Reference Catto2018)

(2.9)\begin{equation} C\left\{f \right\} = \frac{1}{\tau_s v^2} \frac{\partial}{\partial v} [(v^3 + v_c^3)f] + \frac{2 v_\lambda^3 B_{0}}{\tau_s v^3 B} \frac{v_\parallel}{v} \frac{\partial}{\partial \lambda} \left( \lambda \frac{v_\parallel}{v} \frac{\partial f}{\partial \lambda} \right). \end{equation}

(Note that the distinction between the total and the unperturbed field is unimportant in the collision operator.) The first term represents electron and ion drag while the second represents pitch-angle scattering off bulk ions. Here, the alpha slowing down time is given by

(2.10)\begin{equation} \tau_s\left(\psi \right) = \frac{3M_\alpha T_e^{3/2}\left(\psi \right) }{4 \left( 2 {\rm \pi} m_e \right)^{1/2} Z_\alpha ^2 e^4 n_e\left(\psi \right) \ln \varLambda_c}, \end{equation}

with $\ln \varLambda _c$, $T_e$, $n_e$ and $m_e$ being the Coulomb logarithm and the electron temperature, density and mass, respectively. The critical speed at which alpha particles switch from being mainly decelerated by electrons to being mainly decelerated by ions is found by summing over background ions:

(2.11)\begin{equation} v_c^3\left(\psi \right) = \frac{3 {\rm \pi}^{1/2} T_e^{3/2} \left(\psi \right)}{\left(2m_e\right)^{1/2} n_e \left(\psi \right)} \sum_i \frac{Z_i^2 n_i\left(\psi \right)}{M_i}, \end{equation}

with $Z_i$, $n_i$ and $M_i$ being the charge, density and mass of each of the background species. This is of similar size to $v_\lambda$, the speed at which pitch-angle scattering is important to the behaviour of the equilibrium energetic alpha population, which does not have the narrow boundary layers encountered later:

(2.12)\begin{equation} v_\lambda^3\left(\psi \right) \equiv \frac{3 {\rm \pi}^{1/2} T_e^{3/2} \left(\psi \right)}{\left(2m_e\right)^{1/2} n_e \left(\psi \right) M_\alpha } \sum_i Z_i^2 n_i\left(\psi \right). \end{equation}

The drift velocity is given by

(2.13)\begin{equation} \boldsymbol{v}_{d,\textrm{tot}} = \frac{c}{B_{\textrm{tot}}^2} \boldsymbol{E}_{\textrm{tot}}\times \boldsymbol{B}_{\textrm{tot}} + \frac{\lambda_{\textrm{tot}} v^2}{2 B_0\varOmega_{\textrm{tot}}}\hat {\boldsymbol{b}}_{\textrm{tot}} \times \boldsymbol{\nabla} B_{\textrm{tot}} + \frac{v_\parallel^2}{\varOmega_{\textrm{tot}}} \hat {\boldsymbol{b}}_{\textrm{tot}} \times (\hat {\boldsymbol{b}}_{\textrm{tot}} \boldsymbol{\cdot} \boldsymbol{\nabla} \hat {\boldsymbol{b}}_{\textrm{tot}}), \end{equation}

where the different terms represent, respectively, the $\boldsymbol {E} \times \boldsymbol {B}$ drift, the $\boldsymbol {\nabla } B$ drift and the curvature drift. We use $\boldsymbol {v}_d$ to refer to the zeroth-order drift in the presence of only the unperturbed field.

The alpha particle distribution, $f$, consists of an equilibrium component, $f_0$, and a response to any electromagnetic perturbations, $f_1$. Using (2.8), $f_0$ can be calculated to be the familiar slowing down distribution:

(2.14)\begin{equation} f_0\left(\psi, v \right)= \frac{S_{\textrm{fus}}\left(\psi \right) \tau_s\left(\psi \right) H\left(v_0-v\right)}{4 {\rm \pi} [v^3 + v_c^3(\psi)]}. \end{equation}

The derivation of this expression for suprathermal alphas is given in appendix A.

2.3. Example tokamak parameters

When numerical examples are needed in this paper, we use the equilibrium parameters given in table 1, which are similar to those planned for SPARC, a DT tokamak experiment currently under development (Creely et al. Reference Creely2020; Rodriguez-Fernandez et al. Reference Rodriguez-Fernandez, Howard, Greenwald, Creely, Hughes, Wright, Holland, Lin and Sciortino2020).

Table 1.

Example tokamak parameters used in this paper, similar to those planned for SPARC (Creely et al. Reference Creely2020; Rodriguez-Fernandez et al. Reference Rodriguez-Fernandez, Howard, Greenwald, Creely, Hughes, Wright, Holland, Lin and Sciortino2020). The bulk plasma is assumed to be an equal mix of deuterium and tritium. For convenience, this table includes the alpha particle birth speed, $v_0$, even though this parameter is the same in any tokamak.

3.1. General description of perturbations

As discussed in § 2.2, the drift kinetic treatment bounds the scale length of the perturbation to be much longer than the alpha particle gyroradius, but the form of the perturbation is otherwise flexible. In general, perturbations to the tokamak background field can be represented by a perturbed vector potential:

(3.1)\begin{equation} \boldsymbol{A_1} = \boldsymbol{A}_\perp \left( \psi, \vartheta, \zeta,t \right) + A_\parallel \left( \psi , \vartheta, \zeta, t\right) \hat {\boldsymbol{b}}, \end{equation}

and a perturbed electric potential $\varPhi (\psi , \vartheta , \zeta ,t )$. These correspond to a perturbed electric field $\boldsymbol {E}_1 = -\boldsymbol {\nabla } \varPhi - (1/c ) \partial \boldsymbol {A}_1/\partial t$ and a perturbed magnetic field $\boldsymbol {B}_1 = \boldsymbol {\nabla } \times \boldsymbol {A}_1$. We refer to the parallel perturbed magnetic field resulting from $\boldsymbol {A}_{\perp }$ as $B_{1 \parallel }$. The overall magnetic field unit vector $\hat {\boldsymbol{b}}_{\textrm {tot}}$ is related to the perturbation and the background field by

(3.2)\begin{equation} \hat {\boldsymbol{b}}_{\textrm{tot}} \approx \hat {\boldsymbol{b}} + B^{{-}1} (\boldsymbol{\nabla} A_\parallel{\times} \hat {\boldsymbol{b}}). \end{equation}

A perturbation can be characterized by a frequency $\omega$, a toroidal mode number $n$, a poloidal mode number $m$ and an effective radial wavenumber $k_\psi$, or by a sum of components each characterized by these quantities. These fully general forms are maintained in the paper until specific examples are computed in §§ 6 and 7.

3.2. Example perturbations: ripple and TAE

The first example perturbation considered in this paper is ripple, which results from the discrete nature of tokamak field coils. Ripple is composed of a perpendicular vector potential perturbation and corresponding parallel magnetic field $B_{1\parallel }$ and is time-independent. It is characterized by a high $n \approx 18$ and low $m$ in most devices; its radial scale length is roughly determined by the device size or some fraction thereof, so $k_\psi \sim 1/a$ (note that $k_\psi$ depends on several other parameters, including $n$, but we require only knowledge of the rough size of $k_\psi$ for ripple). These characteristics mean that the drift kinetic formalism is appropriate. Ripple is typically of strongest magnitude on the low-field side in regions of high safety factor $q$, high minor radius normalized to major radius $\epsilon$ and high shear $s$. (Information on ripple in ITER can be found in Kurki-Suonio et al. (Reference Kurki-Suonio, Asunta, Hellsten, Hynönen, Johnson, Koskela, Lönnroth, Parail, Roccella and Saibene2009) and information on ripple in SPARC in Scott et al. (Reference Scott, Kramer, Tolman, Snicker, Varje, Särkimäki, Wright and Rodriguez-Fernandez2020).) A summary of these typical values, as well as the specific example values of key ripple parameters used to demonstrate the evaluation of alpha heat flux, is given in table 2.

Table 2.

Defining characteristics of this paper's two demonstration perturbations and typical equilibrium parameters where these perturbations exist. The value $n$ is, for ripple, a common ripple periodicity and, for TAEs, the range of strongly driven toroidal mode numbers. However, other values of $n$ are also sometimes of interest for these perturbations. The specific values of these quantities used for numerical examples in this paper are also included. ‘N/A’ denotes parameters for which specific example values are not needed.

The second example perturbation is the TAE. The TAE is a special class of Alfvén wave that exists in a tokamak; its structure is described by MHD (Van Zeeland et al. Reference Van Zeeland, Kramer, Austin, Boivin, Heidbrink, Makowski, McKee, Nazikian, Solomon and Wang2006). The TAE is destabilized by energetic particles (Fu & Van Dam Reference Fu and Van Dam1989). The TAEs, and many other MHD modes, are time-dependent, with magnetic fields perpendicular to the background field, and result in no parallel electric field (White Reference White2013).Footnote ⁵ They are represented by a vector potential with $\boldsymbol {A}_\perp = 0$ and an electric potential that obeys

(3.3)\begin{equation} {-}\frac{1}{c}\frac{\partial A_\parallel }{\partial t}- \hat {\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} \varPhi = E_\parallel{=} 0. \end{equation}

(This representation of the modes is used frequently in orbit-following codes (Hirvijoki et al. Reference Hirvijoki, Snicker, Korpilo, Lauber, Poli, Schneller and Kurki-Suonio2012).) The TAE frequency is given by $\omega \approx v_A/ (2q R)$, with $v_A$ the on-axis Alfvén speed, $v_A \equiv B_0/\sqrt {4 {\rm \pi} n_i m_i}$, with $m_i$ the average ion mass (Heidbrink Reference Heidbrink2008). The TAEs mostly exist in the core and outer core of the tokamak, where $q$, $\epsilon$ and $s$ have low to moderate values (Rodrigues et al. Reference Rodrigues, Figueiredo, Ferreira, Coelho, Nabais, Borba, Loureiro, Oliver and Sharapov2015). A simple analytic estimate predicts that the most strongly driven TAEs are characterized by toroidal mode numbers $n$ ranging from $\sqrt {\epsilon } R \varOmega _p /(q v_A)$ to $R \varOmega _p /(q v_A)$, depending on whether TAEs are driven primarily by trapped or passing particles (Fu & Cheng Reference Fu and Cheng1992; Breizman & Sharapov Reference Breizman and Sharapov1995; Fülöp et al. Reference Fülöp, Lisak, Kolesnichenko and Anderson1996; Rodrigues et al. Reference Rodrigues, Figueiredo, Ferreira, Coelho, Nabais, Borba, Loureiro, Oliver and Sharapov2015).Footnote ⁶ Alpha-driven TAEs interact with both trapped and passing particle populations, such that the most destabilized values of $n$ will be intermediate between these values. However, for the specific value of $n$ we use to compute numerical examples in this paper, we select a value near the bottom of the range in order to strictly obey the drift kinetic ordering. (Higher $n$ modes are narrower and their width may approach the alpha gyroradius, such that drift kinetics is not a good description of their behaviour.) The poloidal mode number $m$ is given by $nq-1/2$ (Heidbrink Reference Heidbrink2008) and the radial wavenumber $k_\psi$ is given by $k_\psi = nq /(\epsilon R)$ (Breizman & Sharapov Reference Breizman and Sharapov1995; Rodrigues et al. Reference Rodrigues, Figueiredo, Borba, Coelho, Fazendeiro, Ferreira, Loureiro, Nabais, Pinches and Polevoi2016; Tolman et al. Reference Tolman, Loureiro, Rodrigues, Hughes and Marmar2019). The condition (3.3) and the TAE parameters can be used to find $\varPhi = v_A A_\parallel /c$. A summary of the typical TAE characteristics and the specific values used as examples in this paper are given in table 2.

4.1. Initial equation set-up

To first order in the perturbation, the drift kinetic equation (2.8) is

(4.1)\begin{align} &\frac{\partial f_1}{\partial t} +(v_\parallel \hat {\boldsymbol{b}} + \boldsymbol{v}_{d} )\boldsymbol{\cdot} \boldsymbol{\nabla} f_1+ [(v_\parallel \hat {\boldsymbol{b}}_{\textrm{tot}} + \boldsymbol{v}_{d,\textrm{tot}})\boldsymbol{\cdot} \boldsymbol{\nabla} \psi]_1\frac{\partial f_0}{ \partial \psi} \nonumber\\ &\quad + \left[ \frac{Z_\alpha e}{M_\alpha } \left(v_\parallel \hat {\boldsymbol{b}} +\boldsymbol{v}_d \right)\boldsymbol{\cdot} \boldsymbol{E}_1 \right] \frac{\partial f_0}{\partial \mathcal{E}}= C\left\{f_1\right\}, \end{align}

where $f_1$ is the distribution function response to the perturbation and $[(v_\parallel \hat {\boldsymbol{b}}_{\textrm {tot}} + \boldsymbol {v}_{d,\textrm {tot}} )\boldsymbol {\cdot } \boldsymbol {\nabla } \psi ]_1$ is the component of $[(v_\parallel \hat {\boldsymbol{b}}_{\textrm {tot}} + \boldsymbol {v}_{d,\textrm {tot}} )\boldsymbol {\cdot } \boldsymbol {\nabla } \psi ]$ resulting from the perturbation. Further analysis of this equation is simplified by working in a coordinate that follows the alpha particle orbit, which has a finite departure from a flux surface. This variable is the following constant of the alpha particle's motion:Footnote ⁷

(4.2)\begin{equation} \psi_\star{\equiv} \psi - \frac{I v_\parallel}{\varOmega}. \end{equation}

The coordinate $\alpha$ is modified from (2.1) accordingly:

(4.3)\begin{equation} \alpha_\star{\equiv} \zeta - q_\star \vartheta, \end{equation}

with

(4.4)\begin{equation} q_\star{=} q \left(\psi_\star\right) \end{equation}

being the safety factor evaluated at $\psi _\star$ instead of $\psi$. Consistent with the approximation (A 10), we can take that $\partial f_0/\partial \psi _\star \approx \partial f_0/ \partial \psi$.

The perturbed distribution $f_1$ has two parts, such that $f_1 = h + (Z_\alpha e\varPhi /M_\alpha ) \partial f_0/\partial \mathcal {E}$. The second term is an adiabatic response to the perturbation, which does not participate in transport (it is exactly out of phase with the perturbed velocity that causes transport). Evaluation of $[ (v_\parallel \hat {\boldsymbol{b}}_{\textrm {tot}} + \boldsymbol {v}_{d,\textrm {tot}} )\boldsymbol {\cdot } \boldsymbol {\nabla } \psi ]_1$ gives that the equation governing $h$ is

(4.5)\begin{align} &\frac{\partial h} {\partial t} + ( v_\parallel \hat {\boldsymbol{b}}+ \boldsymbol{v}_{d}) \boldsymbol{\cdot} \boldsymbol{\nabla} h - C\left\{ h + \frac{Z_\alpha e\varPhi}{M_\alpha } \frac{\partial f_0}{\partial \mathcal{E}} \right\}={-}\frac{Z_\alpha e}{M_\alpha } \left( \frac{\partial \varPhi}{\partial t} - \frac{v_\parallel}{c} \frac{\partial A_\parallel}{\partial t} + \frac{\lambda v^2}{ 2 c \varOmega} \frac{\partial B_{1\parallel}}{\partial t}\right)\frac{\partial f_0}{\partial \mathcal{E}}\nonumber\\ &\quad - c \left(\frac{\partial \varPhi}{\partial \alpha_\star} - \frac{v_\parallel}{c} \frac{\partial A_\parallel}{\partial \alpha_\star } + \frac{\lambda v^2}{ 2c \varOmega} \frac{\partial B_{1\parallel}}{\partial \alpha_\star} \right) \frac{\partial f_0}{\partial \psi_\star}. \end{align}

Note that the collision operator acting on the adiabatic component is not zero, as it sometimes is in similar calculations.

The perturbed fields $\varPhi$, $A_\parallel$ and $B_{1\parallel }$ are functions of $\psi$, $\zeta$ and $\vartheta$ and can be Fourier-analysed in those coordinates, then rewritten in terms of $\psi$, $\alpha$ and $\vartheta$. For example, the electric potential is given by

(4.6)\begin{align} &\varPhi \left(\psi, \vartheta, \zeta \right) = \sum_{m,n,\omega} \varPhi_{mn\omega} \cos{\left[n\alpha - \omega t + \left(nq -m \right) \vartheta + \int^\psi \,\textrm{d}\psi' \frac{k_\psi}{R B_p} \right]}\nonumber\\ &\quad = \sum_{m,n, \omega}\mathrm{Re}{\left[\varPhi_{mn\omega}\exp\left({\textrm{i} n \alpha - \textrm{i}\omega t +\textrm{i}\left(n q-m \right)\vartheta + \textrm{i} \int^{\psi}\,\textrm{d}\psi' \frac{k_\psi}{R B_p} }\right)\right]}. \end{align}

Moreover, for use in (4.5), the fields can also be written in terms of the starred coordinates, so that, for example,

(4.7)\begin{align} \varPhi \left(\psi, \vartheta, \zeta \right) =\sum_{m,n,\omega}\mathrm{Re}{\left[\varPhi_{mn\omega} \exp\left({\textrm{i} n \alpha_\star{-} \textrm{i}\omega t +\textrm{i}\left(n q_\star{-}m \right) \vartheta + \textrm{i} \int^{\psi_\star}\,\textrm{d}\psi' \frac{k_\psi}{R B_p} + \frac{ik_\psi I v_\parallel}{R B_p\varOmega}}\right)\right]}. \end{align}

Due to the dominance of parallel streaming over drifts, we may take that $\boldsymbol {v}_d \boldsymbol {\cdot } \boldsymbol {\nabla } \vartheta \ll v_\parallel \hat {\boldsymbol{b}}\boldsymbol {\cdot } \boldsymbol {\nabla } \vartheta$. Also, the form of the magnetic field introduced in (2.2) can be used to show $v_\parallel \hat {\boldsymbol{b}} \boldsymbol {\cdot } \boldsymbol {\nabla } h = v_\parallel \hat {\boldsymbol{b}} \boldsymbol {\cdot } \boldsymbol {\nabla } \vartheta \partial h / \partial \vartheta$. Then, (4.5) becomes

(4.8)\begin{align} &\frac{\partial h} {\partial t} + v_\parallel \hat {\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} \vartheta \frac{\partial h}{\partial \vartheta} + \omega_{\alpha_\star} \frac{\partial h}{\partial \alpha_\star}= C\left\{ h +\frac{Z_\alpha e\varPhi}{M_\alpha } \frac{\partial f_0}{\partial \mathcal{E}} \right\}\nonumber\\ &\quad -\textrm{i} \sum_{m,n}D_{n \omega}S_{mn\omega} \exp\left({\textrm{i} n \alpha_\star{-} \textrm{i} \omega t + \textrm{i} \left(nq_\star{-}m \right) \vartheta+ \textrm{i}\int^{\psi_\star}\,\textrm{d}\psi \frac{k_\psi}{RB_p} + \frac{\textrm{i} k_\psi I v_\parallel}{R B_p\varOmega}}\right) , \end{align}

whereFootnote ⁸

(4.9)\begin{equation} \omega_{\alpha_\star} \equiv \boldsymbol{v}_d \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha + v_\parallel \hat {\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} \left( \vartheta \frac{I v_\parallel}{\varOmega} \frac{\partial q}{\partial \psi} \right) = v_\parallel \hat {\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} \vartheta \frac{\partial}{\partial \psi} \left( \frac{B v_\parallel}{\varOmega \hat {\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} \vartheta} \right) \end{equation}

and $S_{mn\omega }$ is given by

(4.10)\begin{equation} S_{mn\omega}\left(\vartheta, v, \lambda \right) = \varPhi_{mn\omega} - \frac{v_\parallel}{c} A_{{\parallel} mn\omega} + \frac{\lambda v^2}{2 c \varOmega} B_{1 \parallel mn\omega}. \end{equation}

(The $\vartheta$ dependence displayed in $S_{mn\omega }$ is due to the dependence of $v_\parallel$ on this quantity.) Also, we define

(4.11)\begin{equation} D_{n \omega}\left(v, \psi_\star\right) \equiv c n \frac{\partial f_0}{\partial \psi_\star}\left(1 - \frac{\omega}{n \omega_\star} \right) \end{equation}

with

(4.12)\begin{equation} \omega_\star{\equiv} \frac{c M_\alpha \partial f_0/ \partial \psi_\star}{Z_\alpha e \partial f_0/\partial \mathcal{E}}. \end{equation}

Note that the drifts on the left-hand side of (4.8) are caused by the unperturbed fields. Thus, we can employ the following form of $h$:

(4.13)\begin{equation} h\left(\psi_\star, \vartheta, \alpha_\star,t, v, \lambda, \sigma \right) =\sum_{n,\omega} h_{n\omega}\left(\vartheta, v, \lambda, \sigma \right) \exp\left({\textrm{i} n \alpha_\star{-}\textrm{i} \omega t+ \textrm{i} \int^{\psi_\star}\,\textrm{d} \psi \frac{k_\psi}{RB_p}}\right). \end{equation}

Here, we have introduced the variable $\sigma$, defined as follows:

(4.14)\begin{equation} \sigma = \left\{ \begin{array}{@{}ll} 0, & \text{trapped particles},\\ \dfrac{v_\parallel}{\left|v_\parallel\right|}, & \text{passing particles}. \end{array} \right. \end{equation}

4.2. Role of collisionality, phenomenological description of transport and form of collision operator

At this point, a discussion of the role of collisionality in the transport reveals the appropriate treatment of the collision operator in (4.8). This role is directly analogous to that played by collisionality in superbanana plateau tokamak transport (see discussion in Shaing (Reference Shaing2015), Calvo et al. (Reference Calvo, Parra, Velasco and Alonso2017) and Catto (Reference Catto2019a)). It is also similar to the role played by collisionality in the plateau regime of neoclassical transport (see Helander & Sigmar Reference Helander and Sigmar2005), in the damping of plasma echoes (see Su & Oberman Reference Su and Oberman1968) and in other wave–particle resonance processes (Duarte et al. Reference Duarte, Gorelenkov, White and Berk2019; Catto Reference Catto2020).

Our analysis of the transport equation (4.8) will reveal that the effect of the perturbations on particle motion is highest for particles which are in resonance with the mode, that is, those for which all terms on the left-hand side of (4.8) vanish. If there were no collisionality, the effect of the perturbation drift (represented by the second term on the left-hand side of (4.8)) on such resonant particles would diverge. So collisionality is necessary to resolve the singularity of the resonance.Footnote ⁹ A small boundary layer of finite width in pitch-angle space and velocity space will collisionally disrupt resonance with the mode, and the perturbed distribution function will vary sharply near this layer. Pitch-angle scattering, represented by the second term in (2.9), rather than drag, is the important collisional behaviour because the second derivative of the pitch-angle scattering operator is stronger near the resonance boundary layer than the first derivative of the drag. A schematic representation of $h$ about a resonant location in pitch-angle space is given in figure 1. A possible physical interpretation of this resonance structure is that particles drift while resonant with the wave until they pitch-angle scatter, giving rise to gradual diffusion (Helander & Sigmar Reference Helander and Sigmar2005).

Figure 1.

A schematic representation of the boundary layer around the resonance. A specific value of pitch angle, $\lambda =\lambda _{\textrm {res}}$, is resonant with the mode, and the non-adiabatic alpha response $h_{n\omega }$ is largest for this value of $\lambda$. Pitch-angle scattering collisions broaden $h_{n\omega }$ about this resonant value. The sharp variation of $h$ near the resonant $\lambda _{\textrm {res}}$ enhances the importance of the pitch-angle scattering operator (given by the second term in (2.9)), which contains a second derivative with respect to $\lambda$.

This understanding yields a phenomenological estimate of the heat flux which should result from the transport being considered. The heat flux will be given by

(4.15)\begin{equation} Q_\alpha\sim M_\alpha v_{\textrm{res}}^2 \delta \lambda \{ [ (v_\parallel \hat {\boldsymbol{b}}_{\textrm{tot}} + \boldsymbol{v}_{d,\textrm{tot}} )\boldsymbol{\cdot} \boldsymbol{\nabla} \psi ]_1 \Delta t\}^2 \Delta t^{{-}1} \frac{\partial n_\alpha}{\partial \psi}, \end{equation}

with $M_\alpha$ the alpha particle mass, $v_{\textrm {res}}$ a characteristic alpha resonant velocity, $\delta \lambda$ the fraction of particles that are resonant with the mode, $[ (v_\parallel \hat {\boldsymbol{b}}_{\textrm {tot}} + \boldsymbol {v}_{d,\textrm {tot}} )\boldsymbol {\cdot } \boldsymbol {\nabla } \psi ]_1$ the radial velocity caused by the perturbation and $\Delta t$ the time that a particle moves before becoming decorrelated by a collision (such that $[ (v_\parallel \hat {\boldsymbol{b}}_{\textrm {tot}} + \boldsymbol {v}_{d,\textrm {tot}} )\boldsymbol {\cdot } \boldsymbol {\nabla } \psi ]_1 \Delta t$ is the step a particle takes before becoming decorrelated). The fraction of particles affected by the mode is estimated from the width of the boundary layer representing the broadened resonance, given by balancing the left-hand side of (4.8) (we use $\omega _{\alpha _\star }$ to include pitch-angle dependence) with the pitch-angle scattering part of the collision operator, (2.9):

(4.16)\begin{equation} n \omega_{\alpha_\star} \delta \lambda \sim \frac{\nu_{\textrm{pas}}}{\delta \lambda^2}, \end{equation}

with $\nu _{\textrm {pas}}$ representing the coefficient of the second term of (2.9). (Note that the effect of the collision operator on the adiabatic component, $(Z_\alpha e\varPhi /M_\alpha ) \partial f_0 / \partial \mathcal {E}$, is negligible in comparison with the sharp variation of $h$ in the resonant boundary layer.) Thus, we have

(4.17)\begin{equation} \delta \lambda \sim \left( \frac{\nu_{\textrm{pas}}}{n \omega_{\alpha_\star} }\right)^{1/3}. \end{equation}

The decorrelation time $\Delta t$ is given by $\Delta t \sim \delta \lambda ^2/\nu _{\textrm {pas}}$, such that estimated flux is

(4.18)\begin{equation} Q_\alpha \sim \frac{M_\alpha v_{\textrm{res}}^2[(v_\parallel \hat {\boldsymbol{b}}_{\textrm{tot}} + \boldsymbol{v}_{d,\textrm{tot}})\boldsymbol{\cdot} \boldsymbol{\nabla} \psi]_1^2}{n\omega_{\alpha_\star}} \frac{\partial n_\alpha}{\partial \psi}. \end{equation}

Note that this flux is independent of the collision frequency $\nu _{\textrm {pas}}$ because $\delta \lambda$ times $\Delta t$ is independent of the collision frequency. The role of the collisions is crucial to the transport, but collisionality does not appear in the final expression for the flux. Thus, in the following steps we use a Krook collision operator (Bhatnagar, Gross & Krook Reference Bhatnagar, Gross and Krook1954; Mynick Reference Mynick1986), that is, an operator for which $C \{ h\} = - \nu h$. The Krook operator $\nu$ is related to $\nu _{\textrm {pas}}$ by $\nu \sim \nu _{\textrm {pas}}/ \delta \lambda ^2$ (recall that the pitch-angle scattering operator (2.9) depends on the second derivative of $h$ with respect to $\lambda$). We also ignore the effect of the collision operator on the adiabatic response because the sharp variation of $h$ is more important. Whenever collisions play a role in the calculation, they allow the convergence of a sum or resolve a singularity, actions that are not qualitatively sensitive to the full form of the collision operator.

4.3. Solution for $h$ with the Krook operator

Inserting the Krook operator and the form (4.13) into (4.8) reveals the equation to be solved for each $h_{n\omega }$:

(4.19)\begin{align} v_\parallel \hat {\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} \vartheta \frac{\partial h_{n\omega}}{\partial \vartheta} - \textrm{i} (\omega - n\omega_{\alpha_\star} + \textrm{i}\nu)h_{n\omega} ={-}\textrm{i}D_{n \omega} \sum_m S_{mn\omega} \exp\left({\textrm{i} \left(nq_\star{-}m \right)\vartheta +\frac{\textrm{i} k_\psi I v_\parallel}{R B_p \varOmega}}\right). \end{align}

Now, recall the variable $\tau$, which characterizes the progression of particles along the magnetic field (see (A 6)). Then, defining

(4.20)\begin{equation} \varLambda\left(\tau\right) \equiv \omega - n \omega_{\alpha_\star}- v_\parallel \hat {\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} \left[\left(nq_\star{-}m \right)\vartheta + \frac{k_\psi I v_\parallel}{R B_p\varOmega} \right], \end{equation}

(4.19) may be written

(4.21)\begin{align} &\frac{\partial}{\partial \tau} \left[h_{n\omega} \exp\left({-\textrm{i}\int_{\tau_0}^\tau\,\textrm{d}\tau ' (\omega - n \omega_{\alpha_\star} + \textrm{i} \nu )}\right) \right] \nonumber\\ &\quad ={-}\textrm{i} D_{n\omega}\sum_mS_{mn\omega} \exp\left({-\textrm{i}\int_{\tau_0}^\tau\,\textrm{d}\tau ' [\varLambda (\tau' ) + \textrm{i} \nu]}\right) \exp\left({\frac{\textrm{i}k_\psi I v_\parallel \left(\vartheta = 0 \right)}{R B_p\varOmega}}\right) . \end{align}

Here, we have chosen $\tau _0$ such that $\vartheta (\tau _0) = 0$. The perturbed distribution function $h_{n\omega }$ can then be found by integrating from $\tau = -\infty$, where $h_{n\omega }= 0$, forward to the ‘present’ time $\tau =0$, with $\vartheta (\tau = 0) = \vartheta$:

(4.22)\begin{align} &\left. h_{n\omega} \exp\left({-\textrm{i}\int_{\tau_0}^\tau\, \textrm{d}\tau ' (\omega - n \omega_{\alpha_\star} + \textrm{i} \nu)}\right) \right|_{\tau = 0} ={-}\textrm{i}D_{n\omega} \exp\left({\frac{\textrm{i}k_\psi I v_\parallel \left(\vartheta = 0 \right)}{R B_p\varOmega}}\right)\nonumber\\ &\quad \times\sum_m \int_{-\infty}^0 S_{mn\omega} \exp\left({-\textrm{i}\int_{\tau_0}^\tau\, \textrm{d}\tau ' [\varLambda(\tau')+ \textrm{i} \nu ]}\right) \,\textrm{d} \tau . \end{align}

Moving the exponential on the left-hand side to the right-hand side gives

(4.23)\begin{equation} h_{n\omega} ={-}\textrm{i}D_{n \omega} \sum_m \exp\left({\textrm{i}\left(nq_\star{-}m\right)\vartheta + \frac{\textrm{i} k_\psi I v_\parallel}{R B_p \varOmega} }\right) \int_{-\infty}^0 S_{mn\omega} \exp\left({-\textrm{i}\int_{0}^\tau\,\textrm{d}\tau ' [\varLambda(\tau') + \textrm{i} \nu]}\right)\,\textrm{d} \tau. \end{equation}

The preceding expression represents the response of the plasma distribution to the presence of a perturbation of the tokamak fields. Next, we work to write this expression in a form that can be more easily evaluated and understood. The expression (4.23) can be split into contributions from each complete bounce (for a trapped particle) or poloidal transit (for a passing particle) of duration given by

(4.24)\begin{equation} \tau_b\equiv \oint \,\textrm{d}\tau, \end{equation}

such that

(4.25)\begin{align} h_{n\omega} &={-}\textrm{i}D_{n\omega} \sum_m \exp\left({\textrm{i}\left(nq_\star{-}m\right)\vartheta + \frac{\textrm{i} k_\psi I v_\parallel}{R B_p\varOmega}}\right) \sum_{j=0}^{\infty} \left[ \int_{-\left(j+1\right)\tau_b }^{{-}j\tau_b}\right] S_{mn\omega} \nonumber\\ &\quad \times \exp\left({-\textrm{i}\int_{0}^\tau \,\textrm{d}\tau' [\varLambda(\tau' ) + \textrm{i} \nu]}\right)\,\textrm{d} \tau. \end{align}

Note that at the trapped–passing boundary, the value of $\tau _b$ for a trapped particle is twice that for a passing particle. This complication will not affect the calculation that follows as velocity space integrations account for both signs of $v_\parallel$. Each portion of the integral may be modified with the change of variables $\tau \rightarrow \tau + j \tau _b$, such that the entire integral reduces to the geometric series

(4.26)\begin{align} h_{n\omega} &={-}\textrm{i}D_{n \omega} \sum_m \exp\left({\textrm{i}\left(nq_\star{-}m\right)\vartheta + \frac{\textrm{i} k_\psi I v_\parallel}{RB_p\varOmega}}\right) \sum_{j=0}^{\infty} \exp\left({\textrm{i}j\oint\, \textrm{d}\tau \left[\varLambda \left(\tau \right) + \textrm{i} \nu \right]}\right) \nonumber\\ &\quad \times \int_{-\tau_b}^0 S_{mn\omega}\exp\left({-\textrm{i}\int_{0}^\tau\,\textrm{d}\tau' [\varLambda(\tau') + \textrm{i} \nu ]}\right)\,\textrm{d} \tau. \end{align}

This transformation requires that $\oint \textrm {d} \tau \varLambda (\tau )$ be the same (up to a factor of $2{\rm \pi}$) for each bounce or transit, which corresponds to the existence of periodicity. Trapped particles, which traverse the same values of $\vartheta$ on each orbit, automatically fulfil this condition, because $\oint \,\textrm {d} \tau \varLambda (\tau )$ is trivially the same for each bounce. Passing particles do not reverse course and traverse values of $\vartheta$ ranging from negative infinity to infinity. For them, the split is only possible when the particle motion is resonant with the wave, such that the contributions to transport from each transit are the same. When the particle is not in phase with the wave, however, (4.23) will approach zero integrated over the entire passing particle trajectory because of the rapid oscillation of the phase factor, unmitigated by periodicity. The following analysis is appropriate for all trapped particles and for resonant passing particles. The analysis reveals in (4.29) onwards that transport comes from resonant particles only, such that the formulation is appropriate for all transport-causing particles, passing and trapped.

Because of the small imaginary part provided by the collisionality, the series (4.26) is convergent and may be summed as

(4.27)\begin{align} h_{n\omega} = \sum_m\frac{\displaystyle-\textrm{i}D_{n\omega} \exp\left({\textrm{i}\left(nq_\star{-}m\right)\vartheta + \dfrac{\textrm{i} k_\psi I v_\parallel}{RB_p\varOmega}}\right) \int_{-\tau_b}^0 S_{mn\omega} \exp\left({-\textrm{i}\int_{0}^\tau \,\textrm{d}\tau' [\varLambda(\tau')+ \textrm{i} \nu]}\right)\,\textrm{d} \tau}{\displaystyle 1- \exp\left({\textrm{i}\oint\,\textrm{d}\tau \left[\varLambda \left(\tau \right)+ \textrm{i} \nu\right]}\right) }. \end{align}

Expanding the denominator about its zeros using $\textrm {e}^{\textrm {i}2 {\rm \pi} l}=1$, with $l$ any integer, denoting the bounce or poloidal transit harmonic or resonant path index, gives

(4.28)\begin{align} h_{n\omega} =\sum_{m,l} \frac{\displaystyle D_{n\omega} \exp\left({\textrm{i}\left(nq_\star{-}m\right)\vartheta + \dfrac{\textrm{i} k_\psi I v_\parallel}{RB_p\varOmega}}\right) \int_{-\tau_b}^0 S_{mn\omega} \exp\left({-\textrm{i}\int_{0}^\tau\,\textrm{d}\tau ' [\varLambda(\tau')+\textrm{i} \nu ]}\right)\,\textrm{d} \tau}{\displaystyle\oint \,\textrm{d}\tau \left[\varLambda \left(\tau \right) + \textrm{i} \nu \right]- 2 {\rm \pi} l }. \end{align}

Alpha particles are in general very unlikely to pitch-angle scatter because of their high energy. Thus, for nearly all particles, $\nu \oint \,\textrm {d}\tau \ll \oint \varLambda \,\textrm {d}\tau - 2 {\rm \pi} l$. In this limit, analogous to the limit that occurs in the plateau regime of neoclassical transport (see Helander & Sigmar Reference Helander and Sigmar2005), (4.28) approaches

(4.29)\begin{align} h_{n\omega} &={-}\textrm{i} {\rm \pi} D_{n \omega} \sum_{m,l} \delta\left( Q_l \right) \exp\left({\textrm{i}\left(nq_\star{-}m\right)\vartheta + \frac{\textrm{i} k_\psi I v_\parallel}{RB_p\varOmega}- \textrm{i} \int_{0}^{\tau_0}\,\textrm{d}\tau ' \varLambda(\tau') }\right)\nonumber\\ &\quad \times \oint S_{mn\omega} \exp\left({- \textrm{i} \int_{\tau_0}^\tau\, \textrm{d}\tau' \varLambda(\tau')}\right)\,\textrm{d} \tau. \end{align}

Here, we have also split the integral in the exponential into two parts which both have limits at $\tau _0$, where $\vartheta (\tau _0)=0$. In addition, we have introduced a delta function $\delta (Q_l)$ with a resonance function as its argument:

(4.30)\begin{equation} Q_l\left(v,\lambda\right) \equiv\oint \varLambda \,\textrm{d}\tau -2 {\rm \pi} l. \end{equation}

Only particles for which $Q_l$ is very close to zero participate in transport, and the strength of pitch-angle scattering is accentuated for them from the normal weak level. The sharp variation in $h_{n\omega }$ between resonant particles and non-resonant particles (visible schematically in figure 1) leads to this accentuation.

These resonances, discussed at greater length in § 5.4, are bounce harmonic resonances, which occur along paths in velocity space when the mode frequency, the particle bounce frequency and the tangential precession frequency resonate. Such resonances are familiar from treatments of neoclassical toroidal viscosity of bulk tokamak plasmas (Linsker & Boozer Reference Linsker and Boozer1982; Mynick Reference Mynick1986; Park et al. Reference Park, Boozer and Menard2009; Kim et al. Reference Kim, Park and Boozer2013; Logan et al. Reference Logan, Park, Kim, Wang and Berkery2013) and also appear in treatments of Alfvén eigenmode stability (e.g. Borba & Kerner Reference Borba and Kerner1999) and in single-particle analyses of energetic particle transport by neoclassical tearing modes, another type of tokamak perturbation (Poli et al. Reference Poli, García-Muñoz, Fahrbach and Günter2008).

Finally, using (4.13) to obtain the full form of $h$ in terms of the unstarred coordinates gives the expression needed to form the heat flux:

(4.31)\begin{align} h &={-}\sum_{m,n,\omega, l} \left\{ i {\rm \pi} \left[ \exp\left({\textrm{i} n \zeta -\textrm{i} \omega t -\textrm{i}m\vartheta + \textrm{i} \int^\psi \,\textrm{d} \psi' \frac{k_\psi}{RB_p} - \textrm{i} \int_{0}^{\tau_0}\,\textrm{d}\tau ' \varLambda(\tau')}\right) \right]D_{n\omega}\delta \left(Q_l\right)\right.\nonumber\\ &\quad \left. \times \oint S_{mn\omega} \exp\left({- \textrm{i} \int_{\tau_0}^\tau\, \textrm{d}\tau' \varLambda(\tau')}\right)\, \textrm{d} \tau \right\}. \end{align}

5.1. Set-up of flux expression

To find the heat flux, the perturbation to the distribution function, $h$, found in (4.31), is multiplied by the perturbed radial velocity that also results from the electromagnetic perturbation, $[ (v_\parallel \hat {\boldsymbol{b}}_{\textrm {tot}} + \boldsymbol {v}_{d,\textrm {tot}} )\boldsymbol {\cdot } \boldsymbol {\nabla } \psi ]_1$.Footnote ¹⁰ The result is weighted by energy and integrated over the distribution function. Then, this quantity is averaged over the flux surface:

(5.1)\begin{equation} Q_\alpha = \left\langle \int\,\textrm{d}^3v (M_\alpha v^2/2) h [(v_\parallel \hat {\boldsymbol{b}}_{\textrm{tot}} + \boldsymbol{v}_{d,\textrm{tot}})\boldsymbol{\cdot} \boldsymbol{\nabla}\psi]_1\right\rangle. \end{equation}

Here, the flux surface average is defined by

(5.2)\begin{equation} \left\langle A \right\rangle \equiv \frac{\oint \dfrac{A \textrm{d}\vartheta\,\textrm{d} \zeta}{\boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{\nabla} \vartheta}}{ \oint \dfrac{\textrm{d}\vartheta\,\textrm{d} \zeta}{\boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{\nabla} \vartheta} } . \end{equation}

To enable analytic evaluation of the flux, we now adopt the following approximate expression for the strength of the axisymmetric magnetic field:

(5.3)\begin{equation} B = B_{0} \left[1- \epsilon\left(\psi \right) \cos \vartheta \right], \end{equation}

with $\epsilon \approx r /R$ (as introduced in § 2) small. Then,

(5.4)\begin{equation} \left\langle A \right\rangle \approx \frac{B_0}{4 {\rm \pi}^2 q R } \oint \frac{A \textrm{d}\vartheta\,\textrm{d} \zeta}{\boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{\nabla} \vartheta}. \end{equation}

The perturbation to the distribution function $h$ may be found in (4.31). The perturbed radial velocity $[ (v_\parallel \hat {\boldsymbol{b}}_{\textrm {tot}} + \boldsymbol {v}_{d,\textrm {tot}} )\boldsymbol {\cdot } \boldsymbol {\nabla } \psi ]_1$ may be read off from (4.5), giving

(5.5)\begin{align} [(v_\parallel \hat {\boldsymbol{b}}_{\textrm{tot}} + \boldsymbol{v}_{d,\textrm{tot}})\boldsymbol{\cdot} \boldsymbol{\nabla}\psi]_1 = \sum_{m,n,\omega} \textrm{i}nc S_{mn\omega} \exp\left({\textrm{i}n\zeta - \textrm{i} \omega t - \textrm{i} m\vartheta + \textrm{i} \int^\psi\, \textrm{d}\psi' \frac{k_\psi}{RB_p}}\right). \end{align}

The flux has contributions from both trapped and passing particles, and will be evaluated in slightly different ways for each of these populations. In general, we illustrate the trapped particle techniques first, then modify them for passing particles later. One of the key differences between the trapped and the passing calculations is the treatment of the particle pitch angle $\lambda$ (2.6). For trapped particles, integration in pitch angle is carried out in terms of the trapping variable $\kappa$, which is defined in terms of $\lambda$:

(5.6)\begin{equation} \kappa^2 = \frac{1-\left(1-\epsilon\right)\lambda}{2 \epsilon \lambda}. \end{equation}

This variable is defined such that, for a trapped particle, the angle of the turning point in the magnetic field (5.3) is given by $2 \sin ^{-1}(\kappa )$. Deeply trapped particles have $\kappa = 0$ and barely trapped particles have $\kappa = 1$. For passing particles, the analogous variable is

(5.7)\begin{equation} k \equiv \frac{1}{\kappa}, \end{equation}

where a value of $k= 1$ corresponds to a barely passing particle and a value of $k= 0$ corresponds to a fully passing particle. (A pedagogical introduction to the trapping parameter can be found in Helander & Sigmar (Reference Helander and Sigmar2005).) With the definitions of $\lambda$ (2.6) and $\kappa$ (5.6) we can write the differential in (5.1) as

(5.8a,b)\begin{equation} \textrm{d}^3v = \frac{{\rm \pi} B v^3\,\textrm{d}v \,\textrm{d} \lambda}{B_{0} v_\parallel}; \quad \textrm{d} \lambda = \frac{-4 \epsilon \kappa d \kappa }{(1 - \epsilon +2 \epsilon \kappa^2)^2}. \end{equation}

An equivalent expression in terms of $k$ holds for passing particles.

5.3. Assumption of $q_\star \approx q$ and final statement of flux

The flux $Q_{\alpha ,mn \omega }$ is a function of flux $\psi$; the perturbation is also a function of this variable. However, $\psi _\star$, defined in (4.2), not $\psi$, is a constant of the alpha particle's motion, such that the alpha particle response is a function of $\psi _\star$. In the flux expression (5.9), $\psi _\star$ appears only in the variable $q_\star$ through the dependence of $Q_l$ and $\varLambda$ on this quantity. Recalling (4.4) and (4.2), we have

(5.12)\begin{equation} q_\star{\approx} q\left(\psi\right) - \frac{\partial q}{\partial \psi} \frac{I v_\parallel}{\varOmega} = q - \frac{qs v_\parallel}{\epsilon R \varOmega_p} \end{equation}

(note that here we have provided $q_\star$ in terms of both $\partial q/\partial \psi$ and shear $s$, (2.4)). Because $v_\parallel$ is a function of $\vartheta$, this quantity varies poloidally in a complicated manner. Such complication introduces significant challenge in the evaluation of the flux, in particular making the term that enforces resonance, $\delta (Q_l)$, dependent on $\vartheta$.

In this work, we thus choose to neglect the difference between $q$ and $q_\star$, leaving a detailed study of the effect of the difference between these quantities to future work. This simplification places constraints on the values of shear $s$ for which our work is valid. Because $q_\star$ appears in the combination $(n q_\star -m )$, we can compare $q_\star$ with the other terms in this combination to develop the following condition for the validity of setting $q_\star = q$:

(5.13)\begin{equation} s \ll \frac{\epsilon R \varOmega_p \left(nq-m\right)}{nqv}. \end{equation}

The degree of restriction that this condition places on the magnitude of the shear depends on the characteristics of the perturbation in question. For ripple, the condition is

(5.14)\begin{equation} s \ll \frac{\epsilon R}{\rho_{p\alpha}}, \end{equation}

a generous condition which allows the approximation to be valid even in the high-shear regions near the tokamak edge where the ripple is important. For the TAE, the equivalent expression is

(5.15)\begin{equation} s \ll \frac{\epsilon R}{nq \rho_{p \alpha}}, \end{equation}

which, given the fairly high values of $n$ characterizing TAEs, is a more restrictive condition, though it should still hold in much of the low-shear core plasma where TAE activity is most common.

With the replacement $q_{\star }\rightarrow q$, the outer $\tau$ integral in (5.9) can be commuted past $\delta (Q_l)$. In addition, using (2.14), we can insert

(5.16)\begin{equation} \frac{\partial f_0}{\partial \psi} \approx \frac{\partial n_\alpha /\partial \psi }{4{\rm \pi} (v^3 + v_c^3) \ln (v_0/v_c) } \approx \frac{\partial n_\alpha /\partial \psi}{4{\rm \pi} v^3 \ln (v_0/v_c)} \end{equation}

into the definition of $D_{n\omega }$ (4.11). (Recall that we have also taken $\partial f_0 /\partial \psi _{\star } \approx \partial f_0 / \partial \psi$.) This gives, for trapped particles,

(5.17)\begin{align} Q_{\alpha,mn\omega ,t} ={-} \sum_l \frac{\epsilon M_\alpha c^2 n^2 \partial n_\alpha/ \partial \psi}{8qR \ln{\left(v_0/v_c\right)}} \int_0^{v_0}\,\textrm{d}v \int_0^1 \,\textrm{d}\kappa \frac{\kappa v^2 \delta\left(Q_l\right) \left(1-\dfrac{\omega}{n\omega_\star} \right)(\mathbb{C}_{mn\omega}^2+ \mathbb{S}_{mn\omega}^2)}{(1-\epsilon +2 \epsilon \kappa^2)^2} . \end{align}

The equivalent expression for passing particles is the following, adding a subscript $p$ to make clear that the expression refers to passing particles:

(5.18)\begin{align} Q_{\alpha,mn \omega ,p} ={-}\sum_l \frac{\epsilon M_\alpha c^2 n^2 \partial n_\alpha/ \partial \psi}{8qR \ln{\left(v_0/v_c\right)}} \int_0^{v_0} \,\textrm{d}v \int_0^1\,\textrm{d}k \frac{k v^2 \delta \left( Q_l \right) \left(1-\dfrac{\omega}{n \omega_\star} \right)(\mathbb{C}_{mn\omega}^2+ \mathbb{S}_{mn\omega}^2)}{[(1-\epsilon)k^2 +2 \epsilon]^2} . \end{align}

These expressions are central results for this paper and can be used to evaluate the alpha flux caused by a wide range of tokamak perturbations. In later sections, they are applied to the flux resulting from ripple and TAEs as examples. Already, the expressions display clear and expected trends. That the flux is driven by the alpha density gradient is shown by the dependence of the expressions on $\partial n_\alpha /\partial \psi$. That the flux increases with the perturbed amplitude squared is also clear. The following subsections examine the elements of these fluxes that remain opaque: the term $\delta (Q_l)$, which enforces the resonance condition, and the term $\mathbb {C}_{mn\omega }^2+ \mathbb {S}_{mn\omega }^2$, which we refer to as a ‘phase factor’.

5.4. Resonance condition

The flux expressions (5.17) and (5.18) include the delta function $\delta (Q_l)$. This delta function enforces the resonance given in (4.30) and determines which particles are able to be transported by the mode and to contribute to the alpha flux. The expression for the resonance (4.30) can be evaluated to read (setting $q_\star = q$, as discussed in § 5.3)

(5.19)\begin{equation} Q_l\left(v,\kappa\right) = \omega \tau_b - n \bar{\omega}_{\alpha_\star} \tau_b - 2 {\rm \pi} \sigma \left(nq-m\right) - 2 {\rm \pi} l, \end{equation}

where the bounce or transit time $\tau _b$ is defined in (4.24), $\sigma$ is defined in (4.14) and the transit average drift is defined by

(5.20)\begin{equation} \bar{\omega}_{\alpha_\star} \equiv \tau_b^{{-}1} \oint \omega_{\alpha_\star}\,\textrm{d} \tau. \end{equation}

The terms in (5.19) can be evaluated for trapped and passing particles to give the expressions in terms of $v$ and $\kappa$ or $k$ in table 3.Footnote ¹¹ The term $\delta (Q_l)$ in the flux expression thus relates the pitch angle and speed of particles that participate in transport and contribute to the flux.

Table 3.

Quantities involved in resonance condition. The complete elliptic integral of the first kind is denoted by $K(\kappa )$; the complete integral of the second kind is $E(\kappa )$.

5.5. Evaluation of phase factor

The term $\mathbb {C}_{mn\omega }^2 + \mathbb {S}_{mn\omega }^2$ in (5.17) and (5.18) depends on integrals defined in (5.10) and (5.11). These integrals are evaluated along particle trajectories, which advance with $\tau$ (A 6). This complex term, which we refer to as a ‘phase factor’, determines how effective a given resonance is at creating a perturbation to the alpha distribution function. In this subsection, we write the phase factor for trapped and passing particles in a form that can be integrated numerically or approximated.

We begin with trapped particles. Trapped particles follow banana orbits, represented schematically in figure 2, with two legs on which the particle velocity is in opposite directions. Our evaluation assumes that the particle orbit is in resonance with the mode at a particular harmonic index $l$ via (5.19). Only such resonant particles are able to contribute to the flux. The integrals in the phase factor are evaluated along each of these legs, which contribute in distinct ways that depend on $l$ and finite orbit effects to the overall phase factor.

Figure 2.

A banana orbit, with legs 1 and 2 defined. Note that banana orbits always move in the sense indicated (i.e. with positive velocity on the outer side of the banana) for plasma current in the direction of positive $\zeta$ because of the conservation of $\psi _\star$ (4.2).

We begin by introducing the variable $x$, which is the following function of $\vartheta$:

(5.21)\begin{equation} \kappa \sin x \equiv \sin \left(\vartheta/2 \right). \end{equation}

The definition of $\kappa$ (5.6) can be used to show that the banana bounce points occur at $x = \pm {\rm \pi} /2$, so that the use of $x$ standardizes the necessary integration. Now, to write $\mathbb {C}_{mn\omega }$ and $\mathbb {S}_{mn\omega }$ explicitly, we must write various quantities that appear in $\tau$ (A 6) or $\varLambda$ (4.20) in terms of $x$. The parallel velocity in the magnetic field (5.3) is given by

(5.22)\begin{equation} v_\parallel{=} \pm \kappa v \sqrt{2\epsilon \lambda} \cos{x} \approx{\pm} \kappa v \sqrt{2\epsilon} \cos{x} , \end{equation}

where the last expression exploits that for trapped particles $\lambda$ differs from $1$ by at most $\pm \epsilon$. The drift term appearing in $\varLambda$, (4.9), is given by

(5.23)\begin{equation} \omega_{\alpha_\star} = \frac{v^2}{2 \varOmega_p R^2} (1 -2 \kappa^2 \sin^2 x + 4 s \kappa^2 \cos^2 x). \end{equation}

We will also use $\hat {\boldsymbol{b}} \boldsymbol {\cdot } \boldsymbol {\nabla } \vartheta \approx 1/(q R)$ to simplify our expressions.

To make our presentation of the phase factor more compact, we define a set of expressions [$a_t(x)$, $b_t$, $c_t(x)$ and $d_t(x)$] given in table 4. Then, the contribution of the outer leg $1$, which has velocity in the direction of positive $\vartheta$, to $\mathbb {C}_{mn\omega }$ (5.10) can be shown to be given by (where we have again taken $\lambda \approx 1$)

(5.24)\begin{align} \mathbb{C}_{mn \omega, \textrm{outer}} &= \frac{2q R \left( \varPhi_{m n \omega} + \dfrac{v^2}{2c \varOmega}B_{1 \parallel m n\omega} \right)}{v\sqrt{2\epsilon}} \displaystyle\int_{-{\rm \pi}/2}^{{\rm \pi}/2}\,{\textrm{d}x} \frac{\cos{\left[ a_t\left(x\right) +b_t -c_t\left(x\right) -d_t\left(x\right) \right]}}{\sqrt{1-\kappa^2 \sin^2 x}}\nonumber\\ &\quad - \frac{ A_{{\parallel} mn \omega} q R}{c} \int_{-{\rm \pi}/2}^{{\rm \pi}/2}{\textrm{d}x} \frac{2\kappa \cos{x} \cos{\left[a_t\left(x\right) +b_t -c_t\left(x\right) -d_t\left(x\right) \right]} }{\sqrt{1-\kappa^2 \sin^2 x}}. \end{align}

Meanwhile, the inner leg is given by

(5.25)\begin{align} \mathbb{C}_{mn \omega, \textrm{inner}} &= \frac{2q R \left( \varPhi_{m n \omega} + \dfrac{v^2}{2c \varOmega}B_{1 \parallel m n\omega} \right)}{v\sqrt{2\epsilon}} \int_{{\rm \pi}/2}^{-{\rm \pi}/2} {\textrm{d}x} \frac{- \cos{\left[ {\rm \pi} l - a_t\left(x\right) +b_t +c_t\left(x\right) -d_t\left(x\right) \right]}}{\sqrt{1-\kappa^2 \sin^2 x}}\nonumber\\ &\quad - \frac{ A_{{\parallel} mn \omega} q R}{c} \int_{{\rm \pi}/2}^{-{\rm \pi}/2} \,{\textrm{d}x} \frac{ 2\kappa \cos{x} \cos{\left[ {\rm \pi} l - a_t\left(x\right) +b_t +c_t\left(x\right) -d_t\left(x\right) \right]}}{\sqrt{1-\kappa^2 \sin^2 x}}. \end{align}

Here, our use of $l$ follows from the recognition that only particles which are in resonance according to (5.19) contribute to the flux, and that the quantity $2{\rm \pi} l$ accumulates from $\omega$ and $n \omega _{\alpha _\star }$ integrated over the entire orbit starting from $x =0$. The inner leg value can be restated by applying the relationship $\cos (x+ {\rm \pi} l ) = (-1)^l \cos x$ and flipping the limits so they agree with those given for the outer leg:

(5.26)\begin{align} \mathbb{C}_{mn\omega l, \textrm{inner}} &= \frac{2q R \left({-}1\right)^l \left( \varPhi_{m n \omega} + \dfrac{ v^2}{2c \varOmega}B_{1 \parallel m n\omega} \right)}{v\sqrt{2\epsilon}} \int_{-{\rm \pi}/2}^{{\rm \pi}/2} {\textrm{d}x} \frac{ \cos{\left[ - a_t\left(x\right) +b_t +c_t\left(x\right) -d_t\left(x\right) \right]}}{\sqrt{1-\kappa^2 \sin^2 x}}\nonumber\\ &\quad + \frac{ A_{{\parallel} mn \omega} q R \left({-}1\right)^l}{c} \int_{-{\rm \pi}/2}^{{\rm \pi}/2} \,{\textrm{d}x} \frac{ 2\kappa \cos{x} \cos{\left[ - a_t\left(x\right) +b_t +c_t\left(x\right) -d_t\left(x\right) \right]}}{\sqrt{1-\kappa^2 \sin^2 x}}. \end{align}

Equivalent expressions hold for $\mathbb {S}_{mn\omega }$. Forming and evaluating the quantity $\mathbb {C}_{mn\omega }^2 + \mathbb {S}_{mn\omega }^2$ reveals that the dependence on $b_t$ vanishes, as is to be expected given that this parameter depends on the parallel velocity evaluated at $\vartheta (\tau _0)= 0$, an arbitrary location. Suppressing the $x$ dependence of $a_t$, $c_t$ and $d_t$ for compactness, the full phase factor reads

(5.27)\begin{align} \frac{\mathbb{C}_{mn\omega}^2 + \mathbb{S}_{mn\omega}^2}{4q^2 R^2} &= \left\{\frac{ \varPhi_{m n \omega} + \dfrac{ v^2}{2c \varOmega}B_{1 \parallel m n\omega}}{v\sqrt{2\epsilon}} \int_{-{\rm \pi}/2}^{{\rm \pi}/2}\,{\textrm{d}x} \frac{ \cos{\left( a_t -c_t -d_t \right)} + \left({-}1\right)^l \cos{\left(- a_t +c_t -d_t \right) }}{\sqrt{1-\kappa^2 \sin^2 x}}\right.\nonumber\\ &\qquad \left. - \frac{A_{{\parallel} mn\omega} \kappa}{c} \int_{-{\rm \pi}/2}^{{\rm \pi}/2}\,{\textrm{d}x} \frac{\cos{x} [\cos{\left( a_t -c_t -d_t \right)} - \left({-}1\right)^l \cos{\left(- a_t +c_t -d_t \right)}]}{\sqrt{1-\kappa^2\sin^2{x}}}\vphantom{\frac{ \varPhi_{m n \omega} + \dfrac{ v^2}{2c \varOmega}B_{1 \parallel m n\omega}}{v\sqrt{2\epsilon}}}\right\}^2\nonumber\\ &\quad + \left\{\frac{ \varPhi_{m n \omega} + \dfrac{ v^2}{2c \varOmega}B_{1 \parallel m n\omega}}{v\sqrt{2\epsilon}} \int_{-{\rm \pi}/2}^{{\rm \pi}/2}\,{\textrm{d}x} \frac{ \sin{\left( a_t -c_t -d_t \right)} + \left({-}1\right)^l \sin{\left(- a_t +c_t -d_t \right) }}{\sqrt{1-\kappa^2 \sin^2 x}}\right.\nonumber\\ &\qquad \left. - \frac{A_{{\parallel} mn\omega} \kappa}{c} \int_{-{\rm \pi}/2}^{{\rm \pi}/2}\,{\textrm{d}x} \frac{\cos{x} [\sin{\left( a_t -c_t -d_t \right)} - \left({-}1\right)^l \sin{\left(- a_t +c_t -d_t \right)}]}{\sqrt{1-\kappa^2\sin^2{x}}}\vphantom{\frac{ \varPhi_{m n \omega} + \dfrac{ v^2}{2c \varOmega}B_{1 \parallel m n\omega}}{v\sqrt{2\epsilon}}}\right\}^2. \end{align}

Though this expression is daunting, in practise many simplifications are possible. First, most modes have only a subset of $\varPhi _{mn\omega }$, $B_{1 \parallel m n \omega }$ and $A_{\parallel mn\omega }$, and these quantities may be proportional to each other. In addition, it is sometimes possible to neglect some of the terms in the arguments of the cosine. The integrals can also be approximated in a variety of ways.

Table 4.

Definitions of quantities used in the trapped particle phase factor.

Passing particles are simpler to treat than trapped particles. For passing particles the variable $x$ is defined as

(5.28)\begin{equation} x \equiv \frac{\vartheta}{2}. \end{equation}

Recalling the variable $k$ defined in (5.7), the parallel velocity in the magnetic field (5.3) for passing particles is given by

(5.29)\begin{equation} v_\parallel{=} \pm \frac{v \sqrt{2\epsilon \lambda}}{k} \sqrt{1-k^2 \sin^2{x}}. \end{equation}

Note that $\lambda$ can be very small for passing particles and thus cannot be set to 1 as it was for the trapped particles. For passing particles, the drift (4.9) is

(5.30)\begin{equation} \omega_{\alpha_\star} = \frac{\lambda v^2}{2 \varOmega_p R^2 k^2} [k^2 -2 k^2 \sin^2 x + 4 s(1-k^2\sin^2x)]. \end{equation}

We again define certain expressions to make the presentation of the phase factor more compact in table 5.

A passing particle orbit goes around the entire device, from $x = -{\rm \pi} /2$ to $x= {\rm \pi} /2$. Its $\mathbb {C}_{mn\omega }$ is given by

(5.31)\begin{align} \mathbb{C}_{mn\omega} &= \frac{2kq R \left(\varPhi_{m n \omega} + \dfrac{\lambda v^2}{2c \varOmega}B_{1 \parallel m n\omega} \right)}{v \sqrt{2\epsilon \lambda}}\int_{-{\rm \pi}/2}^{{\rm \pi}/2} \frac{\textrm{d}x}{\sqrt{1-k^2\sin^2{x}}} \cos{(a_p + b_p -c_p -d_p)}\nonumber\\ &\quad - \frac{2\sigma q R A_{{\parallel} m n \omega}}{c} \int_{-{\rm \pi}/2}^{{\rm \pi}/2} \,{\textrm{d}x} \cos{(a_p +b_p -c_p -d_p)}. \end{align}

An analogous expression holds for $\mathbb {S}_{mn\omega }$. Forming the sum of the two terms gives

(5.32)\begin{align} \frac{\mathbb{C}_{mn\omega}^2 + \mathbb{S}_{mn\omega}^2}{4q^2 R^2} &= \left\{\frac{ k\left( \varPhi_{m n \omega} + \dfrac{\lambda v^2}{2c \varOmega}B_{1 \parallel m n\omega}\right)}{v\sqrt{2\epsilon \lambda}} \int_{-{\rm \pi}/2}^{{\rm \pi}/2} \,{\textrm{d}x} \frac{ \cos{( a_p -c_p -d_p)} }{\sqrt{1-k^2 \sin^2 x}}\right. \nonumber\\ &\quad \left. - \frac{\sigma A_{{\parallel} mn\omega}}{c} \int_{-{\rm \pi}/2}^{{\rm \pi}/2}\,{\textrm{d} x} \cos{( a_p -c_p -d_p )} \vphantom{\frac{ k\left( \varPhi_{m n \omega} + \dfrac{\lambda v^2}{2c \varOmega}B_{1 \parallel m n\omega}\right)}{v\sqrt{2\epsilon \lambda}}}\right\}^2+ \left\{\frac{ k\left( \varPhi_{m n \omega} + \dfrac{\lambda v^2}{2c \varOmega}B_{1 \parallel m n\omega}\right)}{v\sqrt{2\epsilon \lambda}} \right.\nonumber\\ &\quad\left. \times \int_{-{\rm \pi}/2}^{{\rm \pi}/2} \,{\textrm{d}x} \frac{ \sin{(a_p -c_p -d_p)} }{\sqrt{1-k^2 \sin^2 x}} - \frac{\sigma A_{{\parallel} mn\omega} }{c} \int_{-{\rm \pi}/2}^{{\rm \pi}/2}\,{\textrm{d}x} \sin{( a_p -c_p -d_p)} \vphantom{\frac{ k\left( \varPhi_{m n \omega} + \dfrac{\lambda v^2}{2c \varOmega}B_{1 \parallel m n\omega}\right)}{v\sqrt{2\epsilon \lambda}}}\right\}^2. \end{align}

Again, in practise this complicated expression can be simplified in a variety of ways.

Table 5.

Definitions of quantities used in the passing particle phase factor. Note that $\sigma$ is the parameter introduced (4.14), which characterizes the direction of a particle's motion in $\vartheta$.

7.1. Evaluation of flux

The resonance structures for the TAE are shown in figures 7–9. Trapped particles are able to resonate with the TAE at low values of the harmonic index $l$, as seen in figure 7. Near the alpha particle birth speed (the right-hand edge of the plot) the mode resonates with a small subset of pitch angles nearer the trapped–passing boundary ($\kappa = 1$). As the alpha particles slow down (moving left along the $x$ axis), a wider range of pitch angles can resonate with the mode. For a given value of $\kappa$, lower values of $l$ are able to resonate with higher-speed particles.

Figure 7.

Trapped particle ($\sigma = 0$) resonance structure (i.e. where $Q_l(v, \kappa ) = 0$) for a TAE with parameters given in table 2 in a tokamak described by the values in table 1.

Figure 8.

Passing, $\sigma = 1$, particle resonance structure (i.e. where $Q_l(v, \kappa ) = 0$) for a TAE with parameters given in table 2 in a tokamak described by the values in table 1.

Figure 9.

Passing, $\sigma = -1$, particle resonance structure (i.e. where $Q_l(v, \kappa ) = 0$) for a TAE with parameters given in table 2 in a tokamak described by the values in table 1.

Passing particle resonances (with both $\sigma = 1$ and $\sigma = -1$) with the TAE are similar, with alpha particles near birth resonating with barely passing particles and those that have slowed down resonating at a wider variety of pitch angles. These resonance structures can be seen in figures 8 and 9. Again, lower values of $l$ tend to resonate with higher-velocity particles for a given value of $k$. As the drift $\overline {\omega _{\alpha _\star }}$ is negative, an $l=0$ resonance is only possible for $\sigma = 1$, but we will soon find the transport associated with this sign of $\sigma$ is very small. Note that fully passing particles (with $k=0$) resonate at the Alfvén speed and one-third of the Alfvén speed. This aligns with the resonance structure described in simplified treatments of TAE resonance which focus on freely passing particles (Betti & Freidberg Reference Betti and Freidberg1992; Heidbrink Reference Heidbrink2008).

Now, we move to the evaluation of the TAE flux, starting with trapped particles. Again, this procedure begins with consideration of the phase factor (5.27). Full evaluation of this quantity must be done numerically. However, a tolerable approximation of its value can be made for low values of $l$. At low $l$ for typical TAE parameters, the dominant term in the sinusoids of the phase factor is the finite orbit width term $c_t$. This can be verified numerically. After neglecting $a_t$ and $d_t$, the terms in the phase factor which are proportional to $A_\parallel$ can be neglected as higher order in $\epsilon$. Then, the terms $\cos {(c_t)} = \cos {(b_t \cos {x} )}$ (for even $l$) and $\sin {(c_t)} = \sin {(b_t \cos {x} )}$ (for odd $l$) can be written in terms of Bessel functions using the Bessel generating function $\textrm {e}^{\textrm {i} b_t \cos {x}} = \sum _{n = - \infty }^{n= \infty } i^n \textrm {e}^{\textrm {i}nx} \textrm {J}_n(b_t )$. Here, $\textrm {J}_n$ represents the Bessel function of the first kind of order $n$. Higher values of $n$ in the sum are more oscillatory and do not contribute much to the phase factor integral. So, only the lowest $n$ are maintained. This process gives

(7.1)\begin{equation} \frac{\mathbb{C}_{mn\omega}^2 + \mathbb{S}_{mn\omega}^2}{4q^2 R^2}\approx \begin{cases} \dfrac{8 \varPhi_{ m n \omega}^2 \textrm{J}_0\left(b_t\right)^2 K\left(\kappa\right)^2}{\epsilon v^2} , & l \text{ even},\\ \dfrac{32 \varPhi_{mn\omega}^2 \textrm{J}_1\left(b_t\right)^2 (\sin^{{-}1}{\kappa})^2}{\epsilon \kappa^2 v^2}, & l \text{ odd}. \end{cases} \end{equation}

The approximation in (7.1) is compared to a full numerical evaluation of (5.27) for $l=0$, $l=1$ and $l=2$ in figure 10. This comparison shows that, although the approximation is not exact, it captures key trends in the phase factor for low values of $l$. The approximation is less good for higher values of $l$, where $a_t$ becomes more important and cannot be neglected.

Figure 10.

Exact, numerically evaluated phase factor for trapped particles interacting with the TAE compared to the approximate form given in (7.1), for (a) $l=0$, (b) $l=1$ and (c) $l=2$. The approximation reproduces overall trends well for low values of $l$. For higher values of $l$, the neglect of $a_t$ in the approximation is inappropriate; however, higher values of $l$ have very small contributions to the flux. The normalizations in these plots are found in the same way as those in figure 6.

However, higher values of $l$ will have minimal contribution to the flux for three reasons. First, the resonance structures in figures 7–9 make clear that, for a given value of $\kappa$, higher values of $l$ correspond to lower values of $v$. These values of $l$ contribute less heat flux, which depends on $v^2$. Second, the second term in $1 - \omega /(n\omega _\star )$, which appears in the flux expression (5.17), scales with velocity like $1/v^2$ (as we will see later on in (7.3)), further reducing the contribution of lower velocities to the flux. Finally, similar arguments to those used in the consideration of the ripple phase factor show that the sinusoids in the phase factor will become highly oscillatory for high values of $l$, such that the value of the phase factor must shrink by the Riemann–Lebesgue lemma. This can be verified computationally. In figure 10, the phase factor for $l=2$ is already smaller than for $l=0$ and $l=1$, which is an example of this effect. Because higher $l$ indices have such small contributions to the flux, in this section we compute the contributions from $l=0,1,2$ only.

We now proceed with the evaluation of (5.17). To lowest order in $\epsilon$, the denominator of (5.17) is 1. The delta function enforcing the resonance is $\delta (Q_l)$, with $Q_l$ as evaluated in (5.19) with $\sigma = 0$ for trapped particles. This delta function is used to evaluate the velocity integral, recalling the general property that

(7.2)\begin{equation} \delta \left[f \left(x \right)\right] = \sum_i \frac{\delta \left(x - a_i \right)}{\left| \dfrac{\textrm{d}f}{\textrm{d}x}\left(a_i \right) \right|} \end{equation}

for any function $f(x )$ with roots $a_i$. Furthermore, we can write that

(7.3)\begin{equation} 1 - \frac{\omega}{n \omega_\star} = 1 -\frac{3 \varOmega_p \omega R a_\alpha}{n v^2}, \end{equation}

with $a_\alpha$ the alpha scale length defined in (2.7). For typical values of $a_\alpha$, $a_\alpha /R\sim \epsilon$; inserting this scaling and the resonant value of $v$ where $Q_l(v,\kappa ) = 0$ gives that

(7.4)\begin{equation} 1 - \frac{\omega}{n \omega_\star} = 1 - \textit{O}\left(\epsilon \right), \end{equation}

such that $1 - \omega /(n\omega _\star )$ may be set to $1$ to lowest order in $\epsilon$ (at higher values of $l$ than those considered here, where the resonant $v$ becomes very small for many values of $\kappa$, this replacement would not be appropriate). Then, the contribution to the flux from each $l$ reduces to a fairly compact integral over $\kappa$. This integral is different for $l$ even and $l$ odd:

(7.5)\begin{align} Q_{\alpha,mn\omega,t} ={-}\frac{ 4M_\alpha n^2 c^2 \varPhi_{ m n \omega}^2 qR}{\ln{\left(v_0/v_c\right)} } \frac{\partial n_\alpha }{\partial \psi } \sum_l \left\{ \begin{array}{@{}l@{\hspace{12pt}}l} \displaystyle\int_0^{\kappa_0}\,\textrm{d}\kappa \dfrac{\kappa \textrm{J}_0 \left(\dfrac{\sqrt{2\epsilon } k_\psi v_{\textrm{res}} \kappa}{ \varOmega_p}\right)^2 K\left(\kappa\right)^2}{\left|\dfrac{\partial Q_l}{\partial v} \left(v_{\textrm{res}}\right) \right| }, & l\text{ even},\\ \displaystyle\int_0^{\kappa_0}\,\textrm{d} \kappa \dfrac{4 \textrm{J}_1 \left(\dfrac{\sqrt{2\epsilon } k_\psi v_{\textrm{res}} \kappa}{ \varOmega_p}\right)^2 (\sin^{{-}1}{\kappa})^2}{\kappa\left|\dfrac{\partial Q_l}{\partial v} \left(v_{\textrm{res}}\right) \right|}, & l\text{ odd}. \end{array} \right. \end{align}

Here, $v_{\textrm {res}}$ is the resonant velocity at which $Q_l (v_{\textrm {res}}, \kappa ) = 0$. The value of $\kappa$ in resonance at the birth velocity $v_0$ is called $\kappa _0$, i.e. using (5.19) and table 3,

(7.6)\begin{equation} 0 = \frac{8qR K\left(\kappa_0 \right) \omega}{v_0 \sqrt{2\epsilon}} - \frac{4nqv_0 \left[2 E\left(\kappa_0\right) - K\left(\kappa_0\right) \right]}{\varOmega_p R \sqrt{2\epsilon}} - 2 {\rm \pi} l. \end{equation}

The integrals in (7.5) can easily be evaluated numerically for specific values of the relevant parameters.

To obtain a closed expression, we now make another set of approximations for the values of $\kappa _0$. For $l= 0$, $\kappa _0$ is low, and the resonance condition (7.6) is expanded about $\kappa _0 = 0$ to ${\textit{O}}( \kappa _0^2 )$ and then solved for $\kappa _0$. For higher values of $l$, with higher values of $\kappa _0$, different approximations are used. The quantity $2E(\kappa _0 ) - K(\kappa _0)$ vanishes at the point $\kappa _0 \approx 0.91$ and is small near that value. Resonances with $l>0$ have values of $\kappa _0$ in this region. For $l=1$, we can simply approximate $\kappa _0$ with this value. For resonances of higher $l$, we can neglect the term proportional to $2E(\kappa _0 ) - K(\kappa _0)$ in the resonance condition (7.6). Also, we make the replacement $K(\kappa _0)\rightarrow \ln {(4/\sqrt {1-\kappa _0^2})}$ in the remaining terms (this is the limit of $K(\kappa _0)$ as $\kappa _0 \rightarrow 1$). This process gives, where we have inserted $\omega \approx v_A/(2 q R)$ to simplify the expressions,Footnote ¹⁵

(7.7)\begin{equation} \kappa_0 \approx \left\{ \begin{array}{@{}ll} \sqrt{1 - \dfrac{v_A \varOmega_p R}{nqv_0^2}}, & l=0,\\ 0.91, & l=1,\\ \sqrt{1-16 \exp\left({-\dfrac{{\rm \pi} l v_0 \sqrt{2\epsilon}}{v_A}}\right)}, & l > 1. \\ \end{array} \right. \end{equation}

The integrands can be expanded about $\kappa = 0$ to lowest non-trivial order in $\kappa$ (${\textit{O}}(\kappa )$ for $l$ even and ${\textit{O}}(\kappa ^3)$ for $l$ odd). Prior to making the expansion for odd $l$, the velocity in the argument of $\textrm {J}_1$ is replaced with the birth velocity. Comparison of the resulting approximate integrands with the numerically evaluated integrands shows that the approximations are good throughout the integration interval. Then, we calculate

(7.8)\begin{equation} Q_{\alpha,mn\omega,t} \approx{-}\frac{\sqrt{\epsilon} M_\alpha n{\rm \pi} c^2 \varPhi_{m n \omega}^2 \varOmega_p R^2}{4\sqrt{2} \ln{\left(v_0/v_c\right)} }\frac{\partial n_\alpha }{\partial \psi } \sum_l C_{l,t}, \end{equation}

where $C_{l,t}$ is a parameter given in table 6 for the most important low $l$ indices.

Table 6.

Coefficients of key harmonic indices of contributions to the trapped flux, used in (7.8). Only values for $l=0,1,2$ are given because higher values of $l$ have negligible contributions to the flux. We have inserted $\omega \approx v_A/(2qR )$ to simplify the coefficients.

This expression, giving the trapped particle alpha heat flux from a TAE of a given magnitude, is a central result of this paper. For the example parameters in this paper, $\sum _l C_{l,t} \approx 1.07$.

Now, we consider passing particles, with flux given by the expression in (5.18). The first step is to approximate the phase factor (5.32). The dominant term in the phase factor sinusoids is again $c_p$. However, it is not possible to use the Bessel generating function for the passing particles, so instead we simply take $\cos {c_p} \approx \cos {b_p}$ and $\sin {c_p} \approx \sin {b_p}$. Then the phase factor becomes, for $l$ even or odd,

(7.9)\begin{equation} \frac{\mathbb{C}^2_{mn\omega} + \mathbb{S}_{mn\omega}^2 }{4q^2 R^2} = \frac{{\rm \pi}^2 k^2 \varPhi_{m n \omega}^2 }{2\epsilon \lambda v^2} \left(\frac{2K\left(k\right)}{{\rm \pi}} - \frac{\sigma v \sqrt{2\epsilon \lambda} }{kv_A} \right)^2. \end{equation}

This approximation is compared with the full numerical evaluation of the phase factor in table 7.

Table 7.

Exact, numerically evaluated phase factor for passing particles interacting with the TAE compared with the approximate form given in (7.9). The approximation reproduces overall trends well for low values of $l$. Note that the phase factor is much smaller for $\sigma =1$ than for $\sigma =-1$. The normalizations in the plots in this table are different from those used in the analogous plots for trapped particles, found in figures 6 and 10. In addition to being in terms of $k$ and $\lambda$, the normalization is a factor of $2$ higher, which reflects that passing phase factors are intrinsically larger.

Notably, the phase factor is very small for $\sigma = 1$ in both the numerically evaluated and the approximated forms. This results from the cancellation between the contribution to the phase factor from the electric potential and that from the parallel vector potential. Thus, moving forward we evaluate only the contribution from $\sigma = -1$.

The denominator in (5.18) cannot be set to $1$ as it was for the evaluation of the trapped flux (5.17); however, the delta function is applied in the same way as for the trapped particles, using (5.19), and $1- {\omega }/{n \omega _\star }$ is again $1$ to lowest order in $\epsilon$. We set $\sigma = -1$ because, as discussed previously, this is the only significant contribution to the flux. Then the integral to be evaluated is

(7.10)\begin{equation} Q_{\alpha,mn\omega, p} ={-}\frac{M_\alpha n^2 {\rm \pi}^2 c^2 \varPhi_{ m n \omega}^2 qR}{4 \ln{\left(v_0/v_c\right)}} \frac{\partial n_\alpha}{\partial \psi} \sum_l\int_0^{k_0}\,\textrm{d}k \frac{k \left(\dfrac{2 K\left(k\right)}{{\rm \pi}} + \dfrac{v_{\textrm{res}} \sqrt{2\epsilon \lambda}}{kv_A} \right)^2 }{[(1-\epsilon)k^2 + 2 \epsilon] \left| \dfrac{\partial Q_l}{\partial v} \left(v_{\textrm{res}}\right)\right|}, \end{equation}

where $v_{\textrm {res}}$ is the resonant velocity at which $Q_l (v_{\textrm {res}}, k ) = 0$. The upper limit of integration is the value of $k$ in resonance at the birth velocity $v_0$, which can be found for each value of $l$ using (5.19) and the passing quantities in table 3.

For $l=1$, a good approximation for $k_0$ is found from expanding the resonance condition approximately $k_0 = 0$ to ${\textit{O}} ( k_0^2 )$. In contrast, for $l=2$, a good approximation is found by expanding about $k_0 = 1$ to ${\textit{O}} [ \ln (1-k_0) ]$. This process gives, with $\omega \approx v_A/(2 qR)$ used to simplify the expressions,

(7.11)\begin{equation} k_0 = \left\{\begin{array}{@{}ll} \sqrt{2 \epsilon} \sqrt{\dfrac{v_0^2}{v_A^2} -1}, & l=1,\\ 1 - 8 \exp\left({\dfrac{-v_0 (3 \sqrt{2\epsilon} \varOmega_p {\rm \pi} R + 4 n q v_0)} {nqv_0^2 + \varOmega_p R v_A}}\right), & l > 1. \end{array} \right. \end{equation}

The integrand can again be expanded about $k=0$ to lowest non-trivial order in $k$. Then, we find, where we have inserted $\omega \approx v_A/(2qR)$,

(7.12)\begin{equation} Q_{\alpha,mn\omega, p} \approx{-}\frac{M_\alpha n^2 {\rm \pi} v_0 c^2\varPhi_{m n \omega}^2 qR }{ \ln{\left(v_0/v_c\right)} } \frac{\partial n_\alpha}{\partial \psi}\sum_l C_{l,p} , \end{equation}

where $C_{l,p}$ are parameters given in table 8. Again, we include only harmonic indices $l\leq 2$. Higher harmonic indices have small contribution to the flux. This expression, giving the passing particle alpha heat flux from a TAE of a given magnitude, is another central result of this paper. For the example parameters in this paper, $\sum _l C_{l,p} \approx 0.41$. The passing flux is larger than the trapped flux. This results from multiple factors. For instance, there are more passing particles than trapped particles in the isotropic alpha population. Also, the passing phase factor tends to be larger than the trapped phase factor.

Table 8.

Coefficients of key harmonic indices of contributions to the passing flux, used in (7.12).

The expressions for the trapped, (7.8), and passing, (7.12), fluxes can be compared with the phenomenological flux, (4.18). In particular, for the trapped flux, (7.8), setting $n\omega _{\alpha _\star }\sim nv_{\textrm {res}}^2/(\varOmega _p R^2 )$ and ${(\boldsymbol {v}_\parallel \hat {\boldsymbol{b}}_{\textrm {tot}} + \boldsymbol {v}_{d,\textrm {tot}} )_1\boldsymbol {\cdot }\boldsymbol {\nabla }\psi \sim n c \varPhi _{m n \omega }}$ reproduces the dependencies seen in (7.8). For the passing flux, setting $n\omega _{\alpha _\star }\sim \omega \approx v_A/(2qR)$, $v_{\textrm {res}} \sim v_A\sim v_0$ and making the same replacement as for the trapped flux for the radial velocity reproduces the dependencies seen in (7.12). As intuition suggests, the heat flux scales with the square of the perturbation amplitude, $\varPhi _{m n \omega }$, and increases with the periodicity of the mode, $n$.

7.2. Interpretation of flux and development of saturation condition

The magnitude of the alpha flux caused by the TAE is understood by comparing the rate of the diffusion it causes with the rate of slowing down. Slowing down (described in appendix A) removes alpha particles from the energetic population and transfers their energy to the local bulk plasma, leaving the alpha particles as helium ash. So, with $D$ the diffusion coefficient describing the TAE flux and $a$ the device minor radius, if

(7.13)\begin{equation} \frac{D\tau_s}{a^2} \ll 1, \end{equation}

TAE diffusion is weak and does not remove alpha particles before they can give their energy to the background plasma and transition to ash. When

(7.14)\begin{equation} \frac{D\tau_s}{a^2} \sim 1, \end{equation}

TAE diffusion is very strong. It would be capable of diffusing alpha particles across the whole poloidal tokamak cross-section (with an area that scales as $a^2$) in a slowing down time if mode activity of similar strength existed across that area. (Note that the width of an individual mode is significantly smaller than $a$.) At this level, alpha particles diffuse sufficiently fast to cause significant modification to the alpha heat deposition profile. (At this point the perturbative approach used in this paper becomes inappropriate. For losses at this level, a quasilinear treatment that retains the radial flux losses in (A 1) for $f_0$ is needed.) The diffusion coefficient is related to the flux by

(7.15)\begin{equation} D \equiv \frac{-2Q_\alpha}{R^2 B_p^2 \dfrac{\partial n_\alpha }{\partial \psi} M_\alpha v_0^2}, \end{equation}

such that the condition for avoiding strong TAE diffusion, (7.13), becomes

(7.16)\begin{equation} \frac{\sqrt{2\epsilon} \varOmega_p {\rm \pi} v_A^2 R^2 \sum_l C_{l,t} }{4 nv_0^2} \left(\frac{B_{1 m n \omega} }{B} \right)^2 \frac{\tau_s}{a^2} \ll 1 \end{equation}

for trapped particles and

(7.17)\begin{equation} \frac{2 {\rm \pi} q v_A^2 R \sum_l C_{l,p} }{v_0} \left(\frac{B_{1 m n \omega} }{B} \right)^2 \frac{\tau_s}{a^2}\ll 1 \end{equation}

for passing particles. Here, we have stated the flux in terms of the strength of the perturbed magnetic field, $B_{1mn\omega }/B \sim nc \varPhi _{mn\omega } / ( v_A R B_p )$, which is more commonly seen in the literature than $\varPhi _{mn\omega }$. In addition, we have used $\ln {(v_0/v_c) }\approx 1$. For fixed tokamak equilibrium parameters and fixed radial location $\epsilon$, these conditions constrain the maximum amplitude for the mode to avoid serious alpha depletion. We plot the value of $D \tau _s/a^2$ as a function of the normalized mode amplitude, $B_{1mn\omega }/B$, in figure 11.

Figure 11.

The strength of the diffusion caused by a TAE (i.e. the left-hand sides of (7.16) and (7.17)) as a function of amplitude for trapped and passing particles. This quantity must remain well below 1 to avoid significant alpha depletion. The saturated TAE amplitude suggested by a simple model presented in this paper, (7.19), is indicated by the purple dash-dotted line. However, experiments and other models and numerical treatments of saturation often find amplitudes ranging in magnitude from $1 \times 10^{-5}$ to $100 \times 10^{-5}$.

Understanding how this constraint compares to typical TAE amplitudes requires a rough estimate for mode saturation. We now develop such an estimate. Toroidal Alfvén eigenmode saturation has been estimated using a variety of methods. Some works balance the linear growth rate of the instability with the rate at which particles get trapped by the wave and flatten the gradient driving the instability, sometimes including collisions that restore the original gradient (Wu & White Reference Wu and White1994; Fu & Park Reference Fu and Park1995; Wang & Briguglio Reference Wang and Briguglio2016; Zhou & White Reference Zhou and White2016; Todo Reference Todo2019). Other works consider the effects of zonal flows and nonlinear mode couplings (Spong, Carreras & Hedrick Reference Spong, Carreras and Hedrick1994; Chen & Zonca Reference Chen and Zonca2012). Our technique for estimating saturation is similar to the first of these approaches in that it also considers TAE drive reduction due to the flattening of the alpha gradient and collisional refilling of this gradient. However, we consider the removal of the drive to be due to flattening, rather than due to particle trapping.

The mechanism of saturation is represented schematically in figure 12. As the TAE grows, the non-adiabatic particle response $h$ reduces the alpha gradient responsible for mode drive. However, collisions concurrently act to restore the original particle gradient by counteracting this flattening. The balance of these processes results in mode saturation. (If the original drive were not restored fast enough by collisions, $h$ and the mode amplitude would shrink, while if it were restored faster than it was reduced, $h$ and mode amplitude would grow.) Mathematically, this balance is stated as

(7.18)\begin{equation} cnS_{mn\omega}\frac{\partial h_{mn\omega}}{\partial \psi} \sim \nu_{\textrm{pas}} \frac{\partial^2 h_{mn\omega}} {\partial \lambda^2}. \end{equation}

The left-hand side is the rate at which the drive for $h_{mn\omega }$, $cn S_{mn\omega } \partial f_0/\partial \psi$, is reduced by $h_{mn\omega }$. The right-hand side, $\nu _{\textrm {pas}} \partial ^2 h_{mn\omega }/\partial \lambda ^2$, is the rate at which the original drive is restored by collisions.

Figure 12.

Schematic representation of TAE saturation. The non-adiabatic particle response $h$ reduces the alpha gradient $\partial f_0/\partial \psi$ and the drive to the TAE. Collisions (represented by $\nu$) restore the original particle gradient. The balance of these processes, as expressed in (7.18), results in mode saturation.

This saturation condition can be evaluated in two different limits. The first, low-collisionality, limit occurs when $\partial h_{mn\omega } /\partial \psi \sim ( \partial h_{mn\omega }/\partial \kappa ) ( \partial \kappa / \partial \psi )\sim h_{mn\omega }/( R^2 B_p \delta \lambda )$, with $\delta \lambda$ given by (4.17). (This expression is for trapped particles; the equivalent passing expression is in terms of $k$.) The second, high-collisionality, limit occurs when $\partial h_{mn\omega } /\partial \psi \sim k_\psi h_{mn\omega }/(R B_p )$, with $k_\psi$ given by the value for the TAE given in table 2. The transition between the two limits occurs when $k_\psi \sim 1/(R \delta \lambda )$. Inserting $\partial ^2 h_{mn\omega } /\partial \lambda ^2 \sim h_{mn\omega }/ (\delta \lambda )^2$, with $\delta \lambda$ given by (4.17) with $n\omega _{\alpha _\star } \sim \omega$, and $cn S_{mn\omega } \sim c n \varPhi _{m n \omega } \sim v_A R B_p B_1/ B$ reveals the saturated amplitude of the TAE:

(7.19)\begin{equation} \frac{B_1}{B} = \begin{cases} \dfrac{R \nu_{\textrm{pas}}^{2/3} \omega^{1/3} }{v_A}, & \dfrac{\nu_{\textrm{pas}}}{\omega} < \left(\frac{\epsilon}{nq} \right)^3,\\ \dfrac{\epsilon R \nu_{\textrm{pas}}^{1/3} \omega^{2/3} }{nq v_A}, & \dfrac{\nu_{\textrm{pas}}}{\omega} > \left(\frac{\epsilon}{nq} \right)^3. \end{cases} \end{equation}

Here, $\nu _{\textrm {pas}}\sim v_\lambda ^3/(v_0^3 \tau _s)$ (from the second term in the collision operator (2.9)). Similar scalings of saturated amplitude with collisionality are found elsewhere in the literature (Berk et al. Reference Berk, Breizman and Ye1992; Slaby et al. Reference Slaby, Könies, Kleiber and García-Regaña2018; White et al. Reference White, Duarte, Gorelenkov and Meng2019). For the parameters used in tables 1 and 2, the low-collisionality limit is appropriate, and gives an amplitude of roughly $B_1/B \approx 1.1 \times 10^{-5}$. This saturation estimate is comfortably below the stochastic threshold given in Berk et al. (Reference Berk, Breizman and Ye1992), making it consistent with our assumption throughout the paper that the TAE amplitude is below the threshold for stochasticity.

Saturated Alfvén eigenmode amplitudes observed experimentally (Fasoli et al. Reference Fasoli, Borba, Gormezano, Heeter, Jaun, Jacquinot, Kerner, King, Lister and Sharapov1997; Nazikian et al. Reference Nazikian, Fu, Batha, Bell, Bell, Budny, Bush, Chang, Chen and Cheng1997; Van Zeeland et al. Reference Van Zeeland, Kramer, Austin, Boivin, Heidbrink, Makowski, McKee, Nazikian, Solomon and Wang2006) and predicted computationally (Fitzgerald et al. Reference Fitzgerald, Sharapov, Rodrigues and Borba2016) or analytically (Zonca et al. Reference Zonca, Romanelli, Vlad and Kar1995) using other models have shown values ranging from $10^{-5}$ to $10^{-3}$ across a wide range of machine parameters and types of energetic particle populations; in extreme cases, amplitudes of $10^{-2}$ have even been found (Todo, Berk & Breizman Reference Todo, Berk and Breizman2003). Thus, our estimated amplitude is at the bottom of the range typically found in the literature.

The expected saturated amplitude for the example parameters used in this paper, $1.1 \times 10^{-5}$, is indicated on the plot of diffusion strength, figure 11, with a purple dash-dotted line, revealing that at this amplitude, $D\tau _s/a^2$ is small for both the trapped and the passing populations. At this amplitude, significant redistribution of alphas due to this TAE is not expected, though some localized alpha gradient flattening may occur. (In addition, synergistic action with overlapping modes of different $m$ and $n$ could change the transport; since TAEs consist of coupled poloidal harmonics, this change is likely to occur.) An increase in saturated amplitude from the level predicted in our simplified model, but still within the $10^{-5}\text {--}10^{-3}$ range commonly seen in the literature, could lead to significant redistribution. This significant activity is especially likely for the passing population.