2 Equilibrium and NTM perturbation

In this section, we first describe the tokamak magnetic equilibrium in which the alpha particle transport occurs. Then, we describe the NTM perturbation that modifies this equilibrium. Finally, we present the parameters that describe alpha particles in this equilibrium.

2.1 Description of equilibrium

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

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

Here, $\zeta$ is the toroidal angle, with $| \boldsymbol {\nabla } \zeta | = R^{-1}$, where $R$ is the major radius coordinate and $q(\psi )$ is the safety factor. We study the behaviour of alpha particles confined by a magnetic field which is stationary and axisymmetric except for the NTM perturbation. Such a field is stated in the Clebsch representation (Grad & Rubin Reference Grad and Rubin1958; Kruskal & Kulsrud Reference Kruskal and Kulsrud1958; D'haeseleer et al. Reference D'haeseleer, Hitchon, Callen and Shohet2012) 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 )$ characterising the strength of the toroidal magnetic field by $B_\zeta = I(\psi )/R$. We do not consider background electric fields because they do not affect alpha particle trajectories as strongly as the magnetic fields do (a discussion is found in § 2 of Tolman & Catto Reference Tolman and Catto2021). The unit vector corresponding to the magnetic field is

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

The poloidal angle is chosen such thatFootnote ¹

(2.4)\begin{equation} \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} \vartheta = \left|I \right|/\left[q\left(\psi \right)R^2 \right] = q \left(\psi \right)^{-1} \left|\boldsymbol{B} \boldsymbol{\cdot}\boldsymbol{\nabla} \zeta \right|. \end{equation}

At this point, $\vartheta$ is a fully general straightened field line poloidal coordinate which must only be chosen such that (2.4) holds. This allows for shaping and finite aspect ratio. Later, a circular-cross-section, high-aspect-ratio approximation is used.

The poloidal component of the field is denoted by $B_{\rm p}$, and can be found from $B_{\rm p} \approx \epsilon B/q$, with the inverse aspect ratio

(2.5)\begin{equation} \epsilon \approx r /R, \end{equation}

where $r$ is the local minor radius. The flux coordinate and the minor radius are related by

(2.6)\begin{equation} \partial/\partial \psi = \left(1/ RB_{\rm p} \right) \partial /\partial r. \end{equation}

The shear of the field is given by

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

The strength of the on-axis field is given by $B_0$.

The total magnetic field affecting the alpha particles includes both the equilibrium magnetic field and the magnetic field resulting from the NTM, which is discussed in § 2.2. The total magnetic field from all of these sources is denoted $\boldsymbol {B}_{\rm tot}$ (with unit vector $\hat {b}_{\rm 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.

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 being constructed (Creely et al. Reference Creely, Greenwald, Ballinger, Brunner, Canik, Doody, Fülöp, Garnier, Granetz and Gray2020; 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 Creely, Greenwald, Ballinger, Brunner, Canik, Doody, Fülöp, Garnier, Granetz and Gray2020; 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, we include the alpha particle birth speed, $v_0$, even though this parameter is the same in any tokamak.

2.2 Description of NTM perturbation

We consider situations where the equilibrium discussed previously is perturbed by an NTM, which we describe in this section. At rational surfaces $\psi = \psi _{\rm s}$ in the tokamak, the safety factor is given by

(2.8) \begin{equation} q\left( \psi = \psi_{\rm s} \right) = \frac{m}{n}. \end{equation}

An NTM is a resistively-driven instability at this surface that introduces a perturbation to the vector potential (Hirvijoki et al. Reference Hirvijoki, Snicker, Korpilo, Lauber, Poli, Schneller and Kurki-Suonio2012, Reference Hirvijoki, Asunta, Koskela, Kurki-Suonio, Miettunen, Sipilä, Snicker and Äkäslompolo2014) given approximately by

(2.9)\begin{equation} \boldsymbol{A_1} = A_\parallel \left( \psi , \vartheta, \zeta\right) \hat{b}, \end{equation}

where

(2.10)\begin{equation} A_\parallel \left( \psi , \vartheta, \zeta\right) =A_{{\parallel}} \left(\psi\right) \cos{\left(n\zeta - m \vartheta \right)} = {\rm Re}\left[ A_{{\parallel}} \left(\psi \right) \exp({{\rm i} n \alpha + {\rm i} (q n - m) \vartheta})\right]. \end{equation}

We consider only cases where a single NTM is present. This mode will be resonant with the magnetic field at the rational surface. At this surface, the perturbation will not vary at all along the direction of the magnetic field; its wavevector will be perpendicular to the magnetic field. The values of $n$ and $m$ are typically low, as higher mode number modes bend field lines more and, thus, are more stable (Buttery et al. Reference Buttery, Günter, Giruzzi, Hender, Howell, Huysmans, La Haye, Maraschek, Reimerdes and Sauter2000). A brief explanation of how this mode arises is given in Appendix A.

The perturbed magnetic field from (2.9) is given by

(2.11)\begin{equation} \boldsymbol{B}_1 = \boldsymbol{\nabla} \times \boldsymbol{A}_1 . \end{equation}

The overall magnetic field unit vector $\hat {b}_{\rm tot}$ is related to the perturbation and the background field by

(2.12)\begin{equation} \hat{b}_{\rm tot} \approx \hat{b} + B^{-1} \boldsymbol{\nabla} \times \left( B^{-1} A_\parallel \boldsymbol{B} \right) \approx \hat{b} + B^{-1} \boldsymbol{\nabla} \left(I A_\parallel{/}B \right) \times \boldsymbol{\nabla} \zeta, \end{equation}

where the last form preserves axisymmetry and perturbs the flux by

(2.13)\begin{equation} \psi_1 =-I A_\parallel{/}B, \end{equation}

assuming $B_{\rm p} \ll B$.

Later, it will be helpful to have expressions for some quantities involving the perturbed magnetic field. For example, we can state the convenient result, obtained using (2.2), that

(2.14)\begin{equation} \boldsymbol{B}_1 \boldsymbol{\cdot} \boldsymbol{\nabla} \psi = \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \frac{A_\parallel}{B} \boldsymbol{B} \times \boldsymbol{\nabla} \psi \right) = \boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} \left(I A_\parallel{/}B \right) - B \partial A_\parallel{/}\partial\zeta. \end{equation}

Then, we also have that

(2.15)\begin{equation} \boldsymbol{B}_1 \boldsymbol{\cdot}\boldsymbol{\nabla} \psi + \boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi_1 = \boldsymbol{B} \boldsymbol{\cdot}\boldsymbol{\nabla} \left( \psi_1 + I A_\parallel{/}B \right) - B \partial_\zeta A_\parallel=-B \partial_\zeta A_\parallel, \end{equation}

which vanishes upon toroidal averaging, implying that we can define perturbed, axisymmetric flux surfaces.

2.3 Description of particles

We study the behaviour of a population of alpha particles in tokamak magnetic fields perturbed by an NTM of the form (2.9). In this section, we present the system of quasilinear equations used to study this population. Alpha particles are characterised by their speed, $v$ (or, equivalently, energy normalised to mass, $\mathcal {E} \equiv v^2/2$), the sign of $v_\parallel$, the component of their velocity parallel to the equilibrium magnetic field, and their magnetic moment. The sign of $v_\parallel$ is given by the variable $\sigma$, defined as follows:

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

The magnetic moment can be defined relative to the total magnetic field,

(2.17)\begin{equation} \mu_{\rm tot} \equiv \frac{v_{{\perp},{\rm tot}}^2}{2 B_{\rm tot}}, \end{equation}

or relative to the unperturbed field,

(2.18)\begin{equation} \mu \equiv \frac{v_\perp^2}{2B}. \end{equation}

(Here, $v_{\perp,{\rm tot}}$ is the component of their velocity perpendicular to total field and $v_\perp$ is the component of their velocity perpendicular to the unperturbed field.) Note that these definitions of magnetic moment exclude the alpha particle mass. We can also define

(2.19)\begin{equation} \lambda \equiv 2 \mu B_0/v^2. \end{equation}

The alpha particle mass and charge number are given by $M_\alpha$ and $Z_\alpha$, respectively; $\varOmega _{\rm tot} \equiv Z_\alpha eB_{\rm 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 (in the unperturbed field) is $\varOmega _{\rm p} \equiv Z_\alpha eB_{\rm p}/M_\alpha c$. The alpha particle gyroradius is $\rho _\alpha \equiv v_\perp / \varOmega$ and the poloidal gyroradius $\rho _{{\rm p} \alpha } \equiv v_\perp / \varOmega _{\rm p}$.

3 Development of quasilinear formulation

Outside the narrow tearing layer, the NTM perturbation (2.9) has a radial scale length which is large compared with the alpha gyroradius. Therefore, we can begin our study from the drift-kinetic equation (Hazeltine Reference Hazeltine1973):

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

Note that this expression neglects the finite spread of alpha particle birth speeds. Here, the collision operator appropriate for alpha particles with speed much faster than the ion thermal speed is given by (Cordey Reference Cordey1976; Catto Reference Catto2018)

(3.2) \begin{equation} C\left\{\,f \right\} =\frac{1}{\tau_{{\rm s}} v^2} \frac{\partial}{\partial v} \left[ \left(v^3 + v_{{\rm c}}^3 \right) f \right] +\frac{2 v_\lambda^3 B_{0}}{\tau_{{\rm 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}

The first term represents electron and ion drag whereas the second represents pitch angle scattering off of bulk ions. Here, the alpha slowing down time is given by

(3.3) \begin{equation} \tau_{{\rm s}}\left(\psi \right) =\frac{3M_\alpha T_{{\rm e}}^{3/2}\left(\psi \right) }{4 \left( 2 {\rm \pi}m_{{\rm e}} \right)^{1/2} Z_\alpha ^2 e^4 n_{{\rm e}}\left(\psi \right) \ln \varLambda_{{\rm c}}}, \end{equation}

with $\ln \varLambda _{\rm c}$, $T_{\rm e}$, $n_{\rm e}$, and $m_{\rm e}$ the Coulomb logarithm, the electron temperature, the electron density, and the electron mass, respectively. The critical speed at which alpha particles switch from mainly losing energy to electrons to mainly losing energy to ions is found by summing over background ions,

(3.4) \begin{equation} v_{{\rm c}}^3\left(\psi \right) =\frac{3 {\rm \pi}^{1/2} T_{{\rm e}}^{3/2} \left(\psi \right)}{\left(2m_{{\rm e}}\right)^{1/2} n_{{\rm e}} \left(\psi \right)} \sum_{ i} \frac{Z_{\rm i}^2 n_{\rm i}\left(\psi \right)}{M_{\rm i}}, \end{equation}

with $Z_{\rm i}$, $n_{\rm i}$, and $M_{\rm i}$ 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:

(3.5) \begin{equation} v_\lambda^3\left(\psi \right) \equiv\frac{3 {\rm \pi}^{1/2} T_{{\rm e}}^{3/2} \left(\psi \right)}{\left(2m_{{\rm e}}\right)^{1/2} n_{{\rm e}} \left(\psi \right) M_\alpha } \sum_i Z_{{\rm i}}^2 n_{\rm i}\left(\psi \right). \end{equation}

The velocity coordinate contains a component along the magnetic field $\hat {b}_{\rm tot}$, $v_{\parallel,{\rm tot}}$, and a component perpendicular to it, $\boldsymbol {v}_{{\rm d},{\rm tot}}$. The drift velocity is given by

(3.6)\begin{align} \boldsymbol{v}_{{\rm d},{\rm tot}} = \frac{\lambda_{\rm tot} v^2}{2 B_0\varOmega_{\rm tot}}\hat{b}_{\rm tot} \times \boldsymbol{\nabla} B_{\rm tot} + \frac{v_{{\parallel}, {\rm tot}}^2}{\varOmega_{\rm tot}} \hat{b}_{\rm tot} \times \left( \hat{b}_{\rm tot} \boldsymbol{\cdot}\boldsymbol{\nabla} \hat{b}_{\rm tot} \right) \approx \varOmega_{\rm tot}^{-1} v_{{\parallel},{\rm tot}} \boldsymbol{\nabla}_{{\perp}, {\rm tot}} \times \left(v_{{\parallel},{\rm tot}} \hat{b}_{\rm tot} \right), \end{align}

where $\boldsymbol {\nabla }_{\perp,{\rm tot}} = - \hat {b}_{\rm tot} \times (\hat {b}_{\rm tot} \times \boldsymbol {\nabla })$. Later, we use $v_\parallel$ and $\boldsymbol {v}_{\rm d}$ to refer to the parallel velocity and drift in the presence of only the unperturbed field.

For a mode of form (2.9), there is no significant perturbed electric field as any time variation of the perturbed magnetic field is unimportant. In addition, the effect of the tokamak's unperturbed electric field is unimportant relative to the effect of the magnetic drift.Footnote ² Thus, the coefficient of $\partial f/\partial \mathcal{E}$ can be set to zero. In addition, we neglect the time derivative as we are seeking steady-state solutions in a slowly time-varying NTM field. Then, our kinetic equation becomes

(3.7)\begin{equation} \left(v_{{\parallel} ,{\rm tot}} \hat{b}_{\rm tot} + \boldsymbol{v}_{{\rm d},{\rm tot}} \right) \boldsymbol{\cdot} \boldsymbol{\nabla} f= C\left\{\,f\right\} +\frac{S_{\rm fus} \delta \left(v-v_0\right)}{4 {\rm \pi}v^2}. \end{equation}

Our goal is to develop approximate expressions from and solutions of (3.7) that allow us to evaluate the radial energy flux of the passing alpha particles under the influence of the NTM. We begin by recognising that

(3.8)\begin{equation} \left(v_{{\parallel} ,{\rm tot}} \hat{b}_{\rm tot} + \boldsymbol{v}_{{\rm d},{\rm tot}} \right) \boldsymbol{\cdot}\boldsymbol{\nabla} f = \frac{v_{{\parallel}, {\rm tot}}}{B_{\rm tot}} \boldsymbol{\nabla}\boldsymbol{\cdot} \left[ f \left(\boldsymbol{B}_{\rm tot} + \frac{B_{\rm tot}}{v_{{\parallel},{\rm tot}}} \boldsymbol{v}_{{\rm d},{\rm tot}} \right) \right]. \end{equation}

We can average (3.7) over $\zeta$ and over a transit using ${\rm d}\tau _{\rm tot} = {\rm d} \vartheta /(v_{\parallel,{\rm tot}} \hat {b}_{\rm tot} \boldsymbol {\cdot } \boldsymbol {\nabla }\vartheta )$:

(3.9)\begin{align} & \frac{1 }{\displaystyle 2 {\rm \pi}\oint {\rm d}\tau_{\rm tot} }\oint \frac{{\rm d}\tau_{\rm tot} \,{\rm d} \zeta v_{{\parallel},{\rm tot}}}{B_{\rm tot}} \boldsymbol{\nabla} \boldsymbol{\cdot} \left[ f \left(\boldsymbol{B}_{\rm tot} + \frac{B_{\rm tot}}{v_{{\parallel}, {\rm tot}}} \boldsymbol{v}_{{\rm d},{\rm tot}}\right) \right]\nonumber\\ & \quad = \frac{1 }{\displaystyle 2 {\rm \pi}\oint {\rm d}\tau_{\rm tot} } \oint {\rm d}\tau_{\rm tot} \,{\rm d} \zeta \left[C\left\{\,f\right\} +\frac{S_{\rm fus} \delta \left(v-v_0\right)}{4 {\rm \pi}v^2}\right]. \end{align}

Using the divergence in $\psi _{\rm tot} = \psi + \psi _1, \, \vartheta, \, \zeta$ variables, we can simplify the left-hand side of the previous equation:

(3.10)\begin{align} & \frac{1 }{2 \displaystyle{\rm \pi}\oint {\rm d}\tau_{\rm tot} }\oint \frac{{\rm d}\tau_{\rm tot} \,{\rm d} \zeta v_{{\parallel}, {\rm tot}}}{B_{\rm tot}} \boldsymbol{\nabla} \boldsymbol{\cdot} \left[ f \left(\boldsymbol{B}_{\rm tot} + \frac{B_{\rm tot}}{v_{{\parallel}, {\rm tot}}} \boldsymbol{v}_{{\rm d},{\rm tot}} \right) \right]\nonumber\\ & \quad= \frac{1 }{2 \displaystyle{\rm \pi}\oint {\rm d}\tau_{\rm tot} } \frac{\partial}{\partial \psi_{\rm tot}} \oint \frac{\sigma \,{\rm d}\vartheta \,{\rm d} \zeta }{\boldsymbol{B}_{\rm tot} \boldsymbol{\cdot} \boldsymbol{\nabla} \vartheta } \left[ f \left(\boldsymbol{B}_{\rm tot} + \frac{B_{\rm tot}}{v_{{\parallel}, {\rm tot}}} \boldsymbol{v}_{{\rm d},{\rm tot}} \right) \boldsymbol{\cdot} \boldsymbol{\nabla} \psi_{\rm tot} \right]. \end{align}

In this expression, the term $\boldsymbol {v}_{{\rm d},{\rm tot}} \boldsymbol {\cdot }\boldsymbol {\nabla } \psi _{\rm tot}$ drives neoclassical transport, which we do not evaluate here.Footnote ³ Let us now consider the term $f \boldsymbol {B}_{\rm tot} \boldsymbol {\cdot } \boldsymbol {\nabla } \psi _{\rm tot}$. First, consider that in the absence of an NTM, the distribution function must be the axisymmetric

(3.11)\begin{equation} f\left(\psi,\vartheta, v, \sigma, \lambda \right) = f_0 \left(\psi, \vartheta, v, \sigma, \lambda\right). \end{equation}

The introduction of the NTM causes two effects in the distribution function. First, the equilibrium distribution function acquires the structure of the perturbed flux surfaces, i.e. $f_0 (\psi, \vartheta, v, \sigma, \lambda ) \rightarrow f_0 (\psi + \psi _1, \vartheta, v, \sigma, \lambda )$. Second, the distribution function itself may be changed due to the motions resulting from the NTM, i.e. $f_0 (\psi + \psi _1, \vartheta, v, \sigma, \lambda ) \rightarrow f_0 (\psi + \psi _1, \vartheta, v, \sigma, \lambda ) + f_1 (\psi + \psi _1, \vartheta, \zeta, v, \sigma, \lambda )$. That is, we have

(3.12)\begin{equation} f =f_0 \left(\psi + \psi_1, \vartheta, v, \sigma, \lambda \right) + f_1 \left(\psi + \psi_1, \vartheta, \zeta, v, \sigma, \lambda \right). \end{equation}

With this expression, we can evaluate $f \boldsymbol {B}_{\rm tot} \boldsymbol {\cdot } \boldsymbol {\nabla } \psi _{\rm tot}$. Neglecting the very small $\zeta$ dependence in $\boldsymbol {B}_{\rm tot} \boldsymbol {\cdot } \boldsymbol {\nabla } \psi _{\rm tot}$, we note that all first-order terms go to zero when averaged over $\zeta$, and we can write

(3.13)\begin{align} & \frac{1 }{2 \displaystyle{\rm \pi}\oint {\rm d}\tau_{\rm tot} } \frac{\partial}{\partial \psi_{\rm tot}} \oint \frac{\sigma \,{\rm d}\vartheta \,{\rm d} \zeta }{\boldsymbol{B}_{\rm tot} \boldsymbol{\cdot} \boldsymbol{\nabla} \vartheta } f \boldsymbol{B}_{\rm tot} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi_{\rm tot} \nonumber\\ & \quad=\frac{1 }{\displaystyle 2 {\rm \pi}\oint {\rm d}\tau_{\rm tot} } \frac{\partial}{\partial \psi_{\rm tot}} \oint \frac{\sigma \,{\rm d}\vartheta \,{\rm d} \zeta }{\boldsymbol{B}_{\rm tot} \boldsymbol{\cdot} \boldsymbol{\nabla} \vartheta } f \left( \boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi_1 + \boldsymbol{B}_1 \boldsymbol{\cdot} \boldsymbol{\nabla} \psi \right). \end{align}

Making use of (2.15) we find that

(3.14)\begin{equation} \frac{\displaystyle S_{\rm fus} \delta \left(v -v_0\right)\oint_\zeta \,{\rm d}\tau }{4 {\rm \pi}v^2} + \oint_\zeta {\rm d}\tau \,C\left\{\,f_0 \right\}=- \frac{\partial}{\partial \psi} \oint \frac{ \sigma \,{\rm d}\vartheta \,{\rm d} \zeta }{\hat{b} \boldsymbol{\cdot}\boldsymbol{\nabla} \vartheta } f_1 \partial_\zeta A_\parallel. \end{equation}

The subscript $\zeta$ indicates that the integral is performed at a fixed value of $\zeta$. When the right-hand side of (3.14) is negligible, the value of $f_0$ is found to be the slowing down distribution,

(3.15) \begin{equation} f_0 \left(\psi, v\right) =\frac{S_{\rm fus}\left(\psi \right) \tau_{{\rm s}}\left(\psi \right) H\left(v_0-v\right)}{4 {\rm \pi} \left[v^3 + v_{{\rm c}}^3\left(\psi\right)\right]}, \end{equation}

with $H$ a unit step function. To find the energy flux $\varGamma _q$ that results from the right-hand side of (3.14), we can multiply by $M_\alpha v^2/2$ and $v_\parallel d^3 v/B$, where (summing over both signs of $v_\parallel$)

(3.16)\begin{equation} {\rm d}^3 v = \frac{2 {\rm \pi}B v^3 \,{\rm d}v\,{\rm d}\lambda}{ B_0 \left|v_\parallel \right| }, \end{equation}

and divide by $\oint _\zeta {\rm d}\vartheta /(\boldsymbol {B} \boldsymbol {\cdot }\boldsymbol {\nabla } \vartheta )$. We find, assuming large aspect ratio,

(3.17)\begin{align} & \frac{M_\alpha v_0^2}{2}S_{\rm fus} + \frac{M_\alpha }{\displaystyle 2 \oint_\zeta \,{\rm d} \vartheta /\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} \vartheta} \oint_{\zeta} \frac{{\rm d} \vartheta }{\boldsymbol{B} \boldsymbol{\cdot}\boldsymbol{\nabla} \vartheta } \int {\rm d}^3v\, v^2 C \left\{\,f_0 \right\}\nonumber\\ & \quad = \frac{-M_\alpha /2 }{\displaystyle 2 {\rm \pi}\oint {\rm d}\vartheta / \left(\boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{\nabla} \vartheta \right)} \partial_\psi \left[ \oint \frac{{\rm d}\vartheta \,{\rm d} \zeta }{\boldsymbol{B} \boldsymbol{\cdot}\boldsymbol{\nabla} \vartheta } \int {\rm d}^3v \,v^2 v_\parallel f_1 \partial_\zeta A_\parallel\right]\equiv \frac{1}{r} \frac{\partial}{\partial r } \left( r \varGamma_q \right). \end{align}

Then, we can write that the alpha particle energy flux is given by

(3.18) \begin{equation} \varGamma_q \approx-\frac{M_\alpha }{4 {\rm \pi}R B_{{\rm p}} \displaystyle\oint {\rm d}\vartheta / \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}\vartheta} \oint \frac{{\rm d} \vartheta \,{\rm d} \zeta }{\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} \vartheta} \int {\rm d}^3v\,v^2 v_\parallel f_1 \partial_\zeta A_\parallel. \end{equation}

In the next section we find a solution for $f_1$. If unperturbed trajectories are sufficient to calculate the distribution response, it is possible to neglect the nonlinear terms in the kinetic equation.Footnote ⁴ Then, we need only solve the linear equation

(3.19) \begin{equation} \left( v_\parallel \hat{b} + \boldsymbol{v}_{{\rm d}}\right)\boldsymbol{\cdot}\boldsymbol{\nabla} f_1 - v_\parallel \left( \partial_\zeta A_\parallel \right) \left( \partial f_0/\partial \psi \right) = C \left\{ f_1 \right\}. \end{equation}

4 Plasma response to perturbations

In this section, we calculate the plasma response to the NTM perturbation described in the previous section. Our expression (2.9) for the vector potential and the use of unperturbed trajectories in (3.19) suggests using drift kinetic angular momentum $\psi _\star = \psi - I v_\parallel /\varOmega$, $\vartheta$, and $\alpha _\star = \zeta - q_\star \vartheta$ as the variables, with $q_\star = q (\psi _\star )$, since $n\zeta - m\vartheta = n\alpha -(m - nq ) \vartheta = n\alpha _\star - (m -q_\star n )\vartheta$. (Note that these variables mean that to the requisite order $(v_\parallel \hat {b} + \boldsymbol {v}_{\rm d} ) \boldsymbol {\cdot }\boldsymbol {\nabla } \psi _\star = 0$.) Then, we can posit that

(4.1)\begin{equation} f_1 = {\rm Re} {\left[-{\rm i} F_1 \left(\psi_\star, \vartheta, v , \mu, \sigma \right) {\rm e}^{{\rm i} n \alpha_\star }\right]}. \end{equation}

Using (2.10) and (3.19), we find that the equation governing $F_1$ is

(4.2)\begin{align} & {\rm i} v_\parallel \hat{b}\boldsymbol{\cdot}\boldsymbol{\nabla} \vartheta \left. \partial_\vartheta F_1 \right|_{\alpha_\star} + {\rm i} n \omega_\alpha F_1 - \exp({{\rm i}nq_\star \vartheta}) C \left\{F_1 \exp({-{\rm i}n q_\star\vartheta })\right\} \nonumber\\ & \quad=-n v_\parallel \left(\frac{\partial f_0}{\partial \psi} \right) A_{{\parallel}} \exp({-{\rm i} \left(m - q_\star n \right) \vartheta}). \end{align}

Here, we have defined

(4.3) \begin{equation} \omega_{\alpha} \equiv \left(v_\parallel \hat{b} + \boldsymbol{v}_{{\rm d}} \right) \boldsymbol{\cdot}\boldsymbol{\nabla} \alpha_\star{\approx} v_\parallel \hat{b} \boldsymbol{\cdot}\boldsymbol{\nabla} \vartheta \frac{\partial}{\partial \psi} \left( \frac{q I v_\parallel}{\varOmega} \right) . \end{equation}

(Owing to the dominance of parallel streaming over drifts, we have that $v_\parallel \hat {b} \boldsymbol {\cdot }\boldsymbol {\nabla } \vartheta \gg \boldsymbol {v}_{\rm d} \boldsymbol {\cdot }\boldsymbol {\nabla } \vartheta$.) Let us integrate (4.2) over the particle trajectory, assuming $A_\parallel (\psi ) \approx A_\parallel ( \psi _\star )$ because $\psi _\star -\psi = - ( I v_\parallel /\varOmega ) (\partial q/\partial \psi )$ is assumed small, i.e.

(4.4)\begin{align} & \oint {\rm d} \tau \left[{\rm i} v_\parallel \hat{b}\boldsymbol{\cdot}\boldsymbol{\nabla} \vartheta \left. \partial_\vartheta F_1 \right|_{\alpha_\star} + {\rm i} n \omega_\alpha F_1 - \exp({{\rm i}nq_\star \vartheta}) C \left\{F_1 \exp({-{\rm i}n q_\star\vartheta })\right\} \right] \nonumber\\ & \quad =-n A_{{\parallel}} \left(\psi_\star\right) \left(\frac{\partial f_0}{\partial \psi} \right) \oint v_\parallel \exp({-{\rm i} \left(m - q_\star n \right) \vartheta})\,{\rm d}\tau . \end{align}

We note the radial ponderomotive force departure during a poloidal transit due to the nonlinear radial drift from the NTM drive is small compared with a poloidal gyroradius for our small NTM amplitude ordering as is seen by expanding $A_\parallel (\psi +\psi _1) \approx A_\parallel -(I/B)\,{\rm d}A_\parallel ^2/{\rm d}\psi$.

For trapped particles, whose trajectories traverse a closed path in $\vartheta$, the right-hand side of (4.4) vanishes, which shows that the trapped particle response vanishes. The drift caused by the perturbation is proportional to the particle parallel velocity, so particles on closed trajectories along which parallel velocity reverses sign experience no significant drift.

For passing particles there is a significant response when there is a resonance. For the passing particles, we can also Fourier transform in poloidal angle, writing that

(4.5)\begin{equation} f_1 = {\rm Re} {\left[-{\rm i} F_1 \left(\psi_\star, v , \mu, \sigma \right) \exp({{\rm i} n \alpha_\star- {\rm i}\left( m - q_\star n \right)\vartheta })\right]}. \end{equation}

With this new expression, we can again integrate over the passing trajectory, and recognise that $q \approx q_{\star }$ for the collision term and that the magnetic shear term in the drift term will cancel one in the streaming term. We can define a drift without shear,

(4.6)\begin{equation} q v_\parallel \hat{b} \boldsymbol{\cdot}\boldsymbol{\nabla} \vartheta \frac{\partial}{\partial \psi} \left( \frac{I v_\parallel}{\varOmega}\right) \equiv \omega_q, \end{equation}

and write that

(4.7) \begin{equation} 2 {\rm \pi}\sigma {\rm i} \left(nq -m \right) F_1 + {\rm i}n \bar{\omega}_{q} \tau_f F_{1} - \oint {\rm d} \tau C \left\{F_{1}\right\} =-\frac{2{\rm \pi} q R \sigma V \left(v \right)}{v \sqrt{\epsilon}} \frac{\partial f_0}{\partial r}, \end{equation}

where

(4.8) \begin{gather} \tau_f= \oint {\rm d} \tau \approx \frac{4qR \sqrt{\left(1-\epsilon\right) k^2 + 2 \epsilon} K \left(k \right)}{v \sqrt{2\epsilon}} \xrightarrow[k\rightarrow 1]{} - \left( \frac{2qR}{v \sqrt{2\epsilon}}\right) \log{\left(1-k \right)}, \end{gather}

(4.9) \begin{gather} \bar{\omega}_q \tau_f\equiv \oint {\rm d}\tau \,\omega_q, \end{gather}

and we have used quantities defined in table 2. We have adopted the following approximate expression for the strength of the axisymmetric magnetic field:

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

with $\epsilon \approx r /R \ll 1$ (as introduced in § 2). In this field, the parallel velocity is given by $v_\parallel = \pm v \sqrt {[1 - (1-\epsilon ) \lambda ] - 2 \epsilon \lambda \sin ^2{\vartheta /2}}$. We also introduce the parameter $k$, which is defined in terms of $\lambda$,

(4.11)\begin{equation} k^{-2} \equiv \frac{1-\left(1-\epsilon\right)\lambda}{2 \epsilon \lambda}. \end{equation}

A value of $k= 1$ corresponds to a barely passing particle and a value of $k= 0$ corresponds to a freely passing particle. [A pedagogical introduction to this parameter can be found in Helander & Sigmar (Reference Helander and Sigmar2005).]

Table 2. Values defined in order to simplify expressions in the paper.

At a sharp resonance, the pitch angle scattering part of the collision operator (3.2) is most important. The drift resonance is strongest for the barely passing ($k\rightarrow 1$) alpha particles because the transit time $\tau _{\rm f}$ in $\bar {\omega }_q\tau _{\rm f}$ diverges logarithmically. We can set $\lambda \approx 1$ in the pitch angle scattering part of the collision operator. (This simplification reflects that the most important $\lambda$-dependence in the collision operator comes from the second derivative with respect to $\lambda$ of $F_1$.) Then we can write, using quantities in tables 2 and 3,

(4.12)\begin{equation} \oint {\rm d} \tau \,C \left\{F_{1}\right\} \approx \frac{8\sqrt{2\epsilon} \epsilon \nu_{\rm pas}\left(v \right) qR }{v} \frac{\partial^2 F_1}{\partial \lambda^2} \approx \frac{\nu_{\rm pas}\left(v \right) qR}{\sqrt{2 \epsilon} v } \frac{\partial^2 F_1 }{\partial k^2}. \end{equation}

Here, we have also modified the derivative in the collision operator with the approximate expression $\partial / \partial \lambda \approx [ 1/( 4 \epsilon ) ] \partial / \partial k$.

Table 3. Quantities appearing in the kinetic equation. The complete elliptic integral of the first kind is denoted by $K(k )$; the complete integral of the second kind is $E(k)$. (This table uses quantities defined in table 2.)

With the definitions of $\lambda$, (2.18), and $k$, (4.11), we can evaluate the bounce time, the transit averaged drift, and an average present in the collision operator, which are given in table 3. Then, our kinetic equation (4.7) becomes

(4.13) \begin{equation} \omega_{{\rm d}} \left(v \right) \log{\left[\gamma \left(1-k \right) \right]} F_1 + {\rm i} \nu_{\rm pas}\left(v \right) \frac{\partial^2 F_1}{\partial k^2} = {\rm i} 2{\rm \pi} \sqrt{2} \sigma V \left(v \right) \frac{\partial f_0}{\partial r}, \end{equation}

where we have defined the normally large factor

(4.14) \begin{equation} \gamma \left(v \right) \equiv \exp\left(\frac{ 2 {\rm \pi}\sqrt{2\epsilon} \sigma v \left(qn-m \right)}{\omega_{{\rm d}} \left(v \right)q R}\right), \end{equation}

which increases with distance in flux space from the centre of the NTM for $\sigma (qn-m)>0$.

In order to simplify the kinetic equation, let us make an additional set of definitions and approximations. First, we define a variable $x$ such that

(4.15)\begin{equation} x \equiv \gamma \left(1 - k \right), \end{equation}

which represents distance from the trapped–passing boundary and accounts for the distance from the rational surface of interest. We note that $\gamma \gg 1$ and that expanding around $x=1$ gives

(4.16)\begin{equation} \log{x} \approx x -1. \end{equation}

We also define a small quantity

(4.17)\begin{equation} u \left(v \right) \equiv \frac{ \gamma^{2/3} \left( v \right)\nu_{\rm pas}^{1/3} \left(v \right)}{\omega_{d}^{1/3} \left(v \right)}, \end{equation}

which represents the width of the resonance, allowing the definition of an order-unity quantity

(4.18)\begin{equation} z \equiv \frac{x-1}{u}. \end{equation}

Furthermore, we can define a normalised version of the distribution function $\varUpsilon (z )$ such that

(4.19)\begin{equation} F_1 \left(z \right) \equiv-\left[\frac{2 {\rm \pi}\sqrt{2} \sigma u^2 V \left(v \right)}{\nu_{\rm pas} \gamma^2} \frac{\partial f_0}{\partial r }\right] \varUpsilon\left(z \right). \end{equation}

Then, our kinetic equation reduces to

(4.20)\begin{equation} \partial^2 \varUpsilon / \partial z^2 - {\rm i} z \varUpsilon =-1. \end{equation}

This has the particular solution of Su & Oberman (Reference Su and Oberman1968)

(4.21)\begin{equation} \varUpsilon_{{\rm SO}} \left(z \right) = \int_0^{\infty} {\rm d} \tau \exp({-{\rm i}z \tau -\tau^3/3}). \end{equation}

The homogeneous solution, subject to the constraint that it must vanish at $z \rightarrow \infty$, is

(4.22)\begin{equation} \varUpsilon_{h} \left(z \right) = {\rm Ai}\left(z{\rm e}^{{\rm i}{\rm \pi}/6}\right), \end{equation}

where Ai is the Airy function. To find the overall solution, we apply the constraint that the plasma response must vanish at the trapped–passing boundary $z = -1/u$, giving

(4.23)\begin{equation} \varUpsilon = \varUpsilon_{{\rm SO}} \left(z \right) + \varUpsilon_{\sqrt{\nu}} \left(z \right), \end{equation}

where we have defined

(4.24)\begin{equation} \varUpsilon_{\sqrt{\nu}} \left(z \right) \equiv- \varUpsilon_{{\rm SO}} \left(-1/u \right) {\rm Ai}\left(z {\rm e}^{{\rm i} {\rm \pi}/6}\right) / {\rm Ai}\left(-{\rm e}^{{\rm i}{\rm \pi}/6}/u \right). \end{equation}

This expression is plotted in figure 1, showing that $\gamma$ controls the proximity of the response to the freely passing location $k=0$ and $u$ controls the width of the resonance response.

Figure 1. Plasma response function $\varUpsilon$, (4.23), showing that $\gamma$ controls the proximity of the resonance to the trapped–passing boundary and $u$ controls the width of the resonance response.

5 Evaluation of flux

In this section, we show how to use the alpha distribution perturbation $f_1$ derived in the previous section to obtain the resulting alpha energy flux. This energy flux is the quantity of practical use in understanding and predicting tokamak discharges. In addition, we discuss key parts of the expression for this flux.

5.1 Setup of flux expression

The alpha particle energy flux of (3.18) is evaluated by recalling (3.16) to obtain

(5.1) \begin{equation} \varGamma_q \approx-\frac{M_\alpha }{2 R B_{{\rm p}} \displaystyle\oint {\rm d} \vartheta/ \boldsymbol{B} \boldsymbol{\cdot}\boldsymbol{\nabla} \vartheta} \int {\rm d} \vartheta \,{\rm d} \zeta \,{\rm d}v\,{\rm d} \lambda \frac{v^5 v_\parallel}{ \boldsymbol{B} \boldsymbol{\cdot}\boldsymbol{\nabla} \vartheta \left| v_\parallel \right|} f_1 \partial_\zeta A_\parallel. \end{equation}

Using (2.10), (4.5), and (4.19), and the approximate expression $d\lambda \approx 4 \epsilon \,{\rm d}k$, we can simplify this expression to

(5.2) \begin{equation} \varGamma_q =-\frac{4\sqrt{2} {\rm \pi}^2 \epsilon M_\alpha n A_\parallel }{R B_{{\rm p}} \nu_{\rm pas}} \int_0^{v_0} {\rm d}v \int_0^1 {\rm d}k \frac{v^5 u^2 V\left(v \right)}{\gamma^2} \frac{\partial f_0}{\partial r} {\rm Re}\left[ \varUpsilon \left(z \right) \right]. \end{equation}

Using the expression ${\rm d}k = - (u/\gamma )\,{\rm d}z$, the pitch angle integral is performed for $u \ll 1$ by using the asymptotic forms in Abramowitz & Stegun (Reference Abramowitz and Stegun1968) to find

(5.3)\begin{equation} \int_0^1 {\rm d}k\,{\rm Re}{\left[\varUpsilon\left(z \right) \right]} \approx \frac{u{\rm \pi}}{\gamma} \left( 1 + \frac{u^{3/2}}{\sqrt{2} {\rm \pi}} \right). \end{equation}

The first term in (5.3) is the resonant plateau contribution to the energy transport, whereas the $u^{3/2} \propto \nu _{\rm pas}^{1/2}$ term is the small correction due to the presence of the trapped–passing boundary. As important contributions to the speed integral are from the $v$ near the birth speed, we can use the approximate expression

(5.4) \begin{equation} \frac{\partial f_0}{\partial r} \approx\frac{\partial n_\alpha /\partial r }{4 {\rm \pi}v^3 \ln\left(v_0/v_{{\rm c}}\right)} , \end{equation}

along with the speed-dependent (4.17) and the expressions in table 2, to write our energy flux as

(5.5) \begin{equation} \varGamma_q =-\frac{\sqrt{2} {\rm \pi}^2 n A_\parallel^2 M_\alpha q^2 \varOmega_{{\rm p}}}{\sqrt{\epsilon} B_0^2 \ln{\left(v_0/v_{{\rm c}}\right)}}\frac{\partial n_\alpha}{\partial r} \int_0^{v_0} {\rm d}v \frac{v}{\gamma}. \end{equation}

Then, let us write

(5.6)\begin{equation} \gamma = \exp({2\varPi_0 v_0/v}), \end{equation}

where we have defined

(5.7) \begin{equation} \varPi_0 \equiv\frac{\sqrt{2 \epsilon} \varOmega_{{\rm p}} R {\rm \pi} \sigma \left(nq-m \right) }{nv_0q} , \end{equation}

such that we can write

(5.8)\begin{equation} \int_0^{v_0 } \frac{v}{\gamma} {\rm d}v \approx \frac{v_0^2 {\rm e}^{-2\varPi_0} }{2 + 2\varPi_0} \equiv \frac{v_0^2}{2} G \left(\varPi_0 \right), \end{equation}

where $2\varPi _0 \rightarrow 2 \varPi _0 +2$ is inserted to avoid singular behaviour at $qn = m$ and recover the appropriate value as $\varPi _0 \rightarrow 0$. Then, ignoring the small $u^{3/2}$ correction from the trapped–passing boundary, our energy flux is

(5.9)\begin{equation} \varGamma_q =-\frac{\sqrt{2} {\rm \pi}^2 n A_\parallel^2 M_\alpha q^2 \varOmega_{\rm p} v_0^2 G\left(\varPi_0\right)}{2 \sqrt{\epsilon} B_0^2 \ln{\left(v_0/v_{\rm c}\right)}} \frac{\partial n_\alpha}{\partial r}. \end{equation}

From here, we can also define a diffusion coefficient:

(5.10) \begin{equation} D_{\rm rp} =-\frac{2 \varGamma_q }{ \left(\partial n_\alpha/\partial r\right) M_\alpha v_0^2} =\frac{\sqrt{2} {\rm \pi}^2 n A_\parallel^2 q^2 \varOmega_{{\rm p}} G\left(\varPi_0 \right)}{\sqrt{\epsilon} B_0^2 \ln{\left(v_0/v_{{\rm c}}\right)}}. \end{equation}

We next show how a phenomenological estimate of this diffusivity can be made for $\varPi _0 \sim 1$.

5.2 Comparison to phenomenological estimates

Let us compare (5.10) with an expression that can be arrived at via phenomenological methods. In particular, the resonant plateau diffusivity can be estimated as

(5.11)\begin{equation} D_{\rm rp} = F_{\rm trap} V^2/\nu_{\rm eff}, \end{equation}

where $F_{\rm trap}$ is the fraction of particles in resonance, $V$ is the radial drift, and $\nu _{\rm eff}$ is the effective collision frequency. The width of the resonance in $k$ is given by

(5.12) \begin{equation} u_{{\rm k}} = \frac{u}{\gamma} = \left(\frac{\nu_{\rm pas}}{\gamma \omega_{{\rm d}}}\right)^{1/3},\end{equation}

where $u$ is defined in (4.17). (The additional $\gamma$ relative to (4.17) reflects the factor of $\gamma$ in (4.15).) Then, the fraction of passing particles in resonance is

(5.13)\begin{equation} F_{\rm trap} = \sqrt{\epsilon} u_k. \end{equation}

The radial drift is $V(v)$, defined in table 2. The effective collision frequency is

(5.14) \begin{equation} \nu_{\rm eff} = \frac{\nu_{\rm pas}}{u_{{\rm k}}^2} = \left( \nu_{\rm pas}^{1/3} \gamma^{2/3} \omega_{{\rm d}}^{2/3}\right). \end{equation}

Then, the diffusivity can be estimated as

(5.15) \begin{equation} D_{\rm rp} =\frac{\sqrt{\epsilon} V^2}{\gamma \omega_{{\rm d}}} =\frac{ n q^2 \varOmega_{{\rm p}} A_\parallel^2}{\gamma \sqrt{\epsilon} B_0^2}, \end{equation}

in qualitative agreement with (5.10).

5.3 Constraints on transport from quasilinear calculation

To avoid significant transport, we require

(5.16)\begin{align} 1 \gg\frac{D_{\rm rp} \tau_{\rm s}}{a^2} & =\frac{\sqrt{2} n A_\parallel^2 {\rm \pi}^2 q^2 \varOmega_{\rm p} \tau_{\rm s} G\left(\varPi_0 \right)} {\tau_{\rm s} \sqrt{\epsilon} B_0^2\ln{\left(v_0/v_{\rm c}\right)} a^2}\nonumber\\ & =\left( \frac{B_\perp}{B} \right)^2 \frac{ \sqrt{2} \epsilon^{3/2} R^2 {\rm \pi}^2 \varOmega_{\rm p} \tau_{\rm s} G\left(\varPi_0 \right)}{a^2 n \ln\left(v_0/v_{\rm c}\right)}, \end{align}

where in the final equality we have used $B_\perp /B_0 = nA_\parallel /(R B_{\rm p})$. Using values in table 1, we calculate that $\tau _{\rm s} = 0.2 \ {\rm s}$, $\varOmega _{\rm p} = 10^8 \ {\rm s}^{-1}$, $v_0/v_{\rm c} \approx 2.3$, $\epsilon \approx 1/5$, and our expression becomes

(5.17)\begin{equation} \left(B_\perp{/}B_0\right)^2 \ll 10^{-9} / G\left(\varPi_0 \right). \end{equation}

These results show that concern for NTM-driven alpha transport loss is well-founded as (5.17) implies even a very small NTM amplitude can lead to substantial alpha particle energy loss. In the next section, we estimate when we expect our treatment begins to fail because the unperturbed axisymmetric tokamak alpha particle trajectories are perturbed by the NTM magnetic field.