Drift kinetic theory of alpha transport by tokamak perturbations

The coordinates describing the tokamak equilibrium are $\psi$ (the poloidal flux function), $\vartheta$ (the poloidal angle)³³3At this point, $\vartheta$ is a fully general straightened field line poloidal coordinate which must only be chosen such that $q(\psi)=(\hat{b}\cdot\nabla\zeta)/(\hat{b}\cdot\nabla\vartheta)$. This allows for shaping and finite aspect ratio. Later, in section 5, a circular-cross-section, high-aspect-ratio approximation is used to obtain analytic expressions for the flux., and $\alpha$, defined by

$$\alpha\equiv\zeta-q\left(\psi\right)\vartheta.$$

(1)

Here, $\zeta$ is the toroidal angle. The safety factor is $q$, which characterizes the twist of the background magnetic field by $q\left(\psi\right)\equiv\hat{b}\cdot\nabla\zeta/\left(\hat{b}\cdot\nabla\vartheta\right)$ with $\hat{b}$ giving the direction of the equilibrium magnetic field.⁴⁴4Later in the paper, an approximate form of this definition, $q\left(\psi\right)\approx\left(R\hat{b}\cdot\nabla\vartheta\right)^{-1}$, with $R$ the major radius coordinate, will be used to simplify expressions. We study the behavior 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

$$\vec{B}=\nabla\alpha\times\nabla\psi=I\left(\psi\right)\nabla\zeta+\nabla\zeta\times\nabla\psi,$$

(2)

with $I\left(\psi\right)$ characterizing the strength of the toroidal magnetic field by $B_{\zeta}=I\left(\psi\right)/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 the magnetic fields do.⁵⁵5This follows from the expression for the particle drifts, given in (13), where we can find that for alphas the tangential component of the $\vec{E}\times\vec{B}$ drift is given by $\vec{v}_{\vec{E}\times\vec{B}}\cdot\nabla\alpha=-c\partial\Phi_{0}/\partial\psi$, with $\Phi_{0}$ the background potential, while the $\nabla B$ drift is given by $\vec{v}_{\nabla B}\cdot\nabla\alpha\sim-\rho_{\alpha}v_{0}/\left(aR\right)$, with $\rho_{\alpha}$ the alpha gyroradius and $v_{0}$ the alpha birth speed. Typically, $\partial\Phi_{0}/\partial\psi\sim-\left(T_{i}/en_{i}\right)\partial n_{i}/\partial\psi$, with $T_{i}$ and $n_{i}$ the background ion temperature and density, respectively. Taking the scale length of the ion population to be the device minor radius $a$ gives $\vec{v}_{\vec{E}\times\vec{B}}\cdot\nabla\alpha\sim\rho_{i}v_{i}/a^{2}$, with $\rho_{i}$ the ion gyroradius. The tangential component of the $\nabla B$ drift is then larger than the tangential component of the $\vec{E}\times\vec{B}$ drift by a factor of $\epsilon\rho_{\alpha}v_{0}/\left(\rho_{i}v_{i}\right)$, which is large. As noted previously, the unit vector corresponding to this field is

$$\frac{\vec{B}}{B}\equiv\hat{b}.$$

(3)

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=\left(1/RB_{p}\right)\partial/\partial r$.] The shear of the field is given by

$$s\equiv\left(r/q\right)\partial q/\partial r.$$

(4)

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 section 3. The total magnetic field from all of these sources is denoted $\vec{B}_{tot}$ (with unit vector $\hat{b}_{tot}$) and the total electric field is $\vec{E}_{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 velocity, $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,

$$\lambda_{tot}=\frac{B_{0}v_{\perp,tot}^{2}}{B_{tot}v^{2}},$$

(5)

with $B_{0}$ the equilibrium on-axis magnetic field strength, or relative to the unperturbed field,

$$\lambda=\frac{B_{0}v_{\perp}^{2}}{Bv^{2}}.$$

(6)

(Here, $v_{\perp,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; $\Omega_{tot}\equiv Z_{\alpha}eB_{tot}/M_{\alpha}c$ is the gyrofrequency in the total field and $\Omega\equiv Z_{\alpha}eB/M_{\alpha}c$ the gyrofrequency in the unperturbed field. The alpha particle poloidal gyrofrequency is $\Omega_{p}\equiv Z_{\alpha}eB_{p}/M_{\alpha}c$. The alpha particle gyroradius is $\rho_{\alpha}\equiv v_{\perp}/\Omega$, and the poloidal gyroradius $\rho_{p\alpha}\equiv v_{\perp}/\Omega_{p}$.

We study the effect of perturbations on the energetic alpha particle distribution function $f$ using Hazeltine’s drift kinetic equation (Hazeltine, 1973). This formalism is appropriate when $\rho_{\alpha}$ (which at the alpha particle birth speed is typically a couple of centimeters) 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,

$$a_{\alpha}\equiv-\frac{n_{\alpha}}{RB_{p}\partial n_{\alpha}/\partial\psi},$$

(7)

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

$$\frac{\partial f}{\partial t}+\left(v_{\parallel,tot}\hat{b}_{tot}+\vec{v}_{d,tot}\right)\cdot\nabla f+\left[\frac{Z_{\alpha}e}{M_{\alpha}}\left(v_{\parallel,tot}\hat{b}_{tot}+\vec{v}_{d,tot}\right)\cdot\vec{E}_{tot}+\frac{\lambda_{tot}v^{2}}{2B_{0}}\frac{\partial B_{tot}}{\partial t}\right]\frac{\partial f}{\partial\mathcal{E}}\\ =C\left\{f\right\}+\frac{S_{fus}\delta\left(v-v_{0}\right)}{4\pi v^{2}}.$$

(8)

Here, $S_{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 $\vec{E}_{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,tot}$. The energetic alpha collision operator represents collisions with the background plasma electrons and ions, and is given by (Cordey, 1976; Catto, 2018)

$$C\left\{f\right\}=\frac{1}{\tau_{s}v^{2}}\frac{\partial}{\partial v}\left[\left(v^{3}+v_{c}^{3}\right)f\right]+\frac{2v_{\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).$$

(9)

(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 of bulk ions. Here, the alpha slowing down time is given by

$$\tau_{s}\left(\psi\right)=\frac{3M_{\alpha}T_{e}^{3/2}\left(\psi\right)}{4\left(2\pi m_{e}\right)^{1/2}Z_{\alpha}^{2}e^{4}n_{e}\left(\psi\right)\ln\Lambda_{c}},$$

(10)

with $\ln\Lambda_{c}$, $T_{e}$, $n_{e}$, and $m_{e}$, the Coulomb logarithm, 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,

$$v_{c}^{3}\left(\psi\right)=\frac{3\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}},$$

(11)

with $Z_{i}$, $n_{i}$, and $M_{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 behavior of the equilibrium energetic alpha population, which does not have the narrow boundary layers encountered later:

$$v_{\lambda}^{3}\left(\psi\right)\equiv\frac{3\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).$$

(12)

The drift velocity is given by

$$\vec{v}_{d,tot}=\frac{c}{B_{tot}^{2}}\vec{E}_{tot}\times\vec{B}_{tot}+\frac{\lambda_{tot}v^{2}}{2B_{0}\Omega_{tot}}\hat{b}_{tot}\times\nabla B_{tot}+\frac{v_{\parallel}^{2}}{\Omega_{tot}}\hat{b}_{tot}\times\left(\hat{b}_{tot}\cdot\nabla\hat{b}_{tot}\right),$$

(13)

where the different terms represent, respectively, the $\vec{E}\times\vec{B}$ drift, the $\nabla B$ drift, and the curvature drift. We use $\vec{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 (8), $f_{0}$ can be calculated to be the familiar slowing down distribution:

$$f_{0}\left(\psi,v\right)=\frac{S_{fus}\left(\psi\right)\tau_{s}\left(\psi\right)H\left(v_{0}-v\right)}{4\pi\left[v^{3}+v_{c}^{3}\left(\psi\right)\right]}.$$

(14)

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

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., 2020; Rodriguez-Fernandez et al., 2020).

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,

$$\vec{A_{1}}=\vec{A}_{\perp}\left(\psi,\vartheta,\alpha,t\right)+A_{\parallel}\left(\psi,\vartheta,\alpha,t\right)\hat{b},$$

(15)

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

$$\hat{b}_{tot}\approx\hat{b}+B^{-1}\left(\nabla A_{\parallel}\times\hat{b}\right).$$

(16)

A perturbation can be characterized by a frequency $\omega$, a toroidal mode number $n$, a poloidal mode number $m$, and an effective radial wave number $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 sections 6 and 7.

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. (2009); information on ripple in SPARC in Scott et al. (2020).] 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.

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 magnetohydrodynamics (Van Zeeland et al., 2006). The TAE is destabilized by energetic particles (Fu & Van Dam, 1989). 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, 2013).⁶⁶6Inclusion of finite ion gyroradius effects introduces a parallel electric field, but these effects are small and are outside the scope of the MHD model often used to describe TAE structure. They are represented by a vector potential with $\vec{A}_{\perp}=0$ and an electric potential that obeys

$$-\frac{1}{c}\frac{\partial A_{\parallel}}{\partial t}-\hat{b}\cdot\nabla\Phi=E_{\parallel}=0.$$

(17)

[This representation of the modes is used frequently in orbit-following codes (Hirvijoki et al., 2012).] The TAE frequency is given by $\omega\approx v_{A}/\left(2qR\right)$, with $v_{A}$ the on-axis Alfvén speed, $v_{A}\equiv B_{0}/\sqrt{4\pi n_{i}m_{i}}$, with $m_{i}$ the average ion mass (Heidbrink, 2008). 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., 2015). A simple analytic estimate predicts that the most strongly driven TAEs are characterized by toroidal mode numbers $n$ ranging from $\sqrt{\epsilon}R\Omega_{p}/\left(qv_{A}\right)$ to $R\Omega_{p}/\left(qv_{A}\right)$, depending on if TAEs are driven primarily by trapped or passing particles (Fu & Cheng, 1992; Breizman & Sharapov, 1995; Fülöp et al., 1996; Rodrigues et al., 2015).⁷⁷7This estimate results from balancing the width of the mode $1/k_{\psi}\sim\epsilon R/\left(nq\right)$, with the width of the particle orbit of resonant alpha particles (which scales with $v_{A}\sqrt{\epsilon}/\Omega_{p}$ for trapped particles and with $v_{A}/\Omega_{p}$ for passing particles). Various variations of this balance exist in the literature. 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 behavior.) The poloidal mode number $m$ is given by $nq-1/2$ (Heidbrink, 2008) and the radial wave number $k_{\psi}$ is given by $k_{\psi}=nq/\left(\epsilon R\right)$ (Breizman & Sharapov, 1995; Rodrigues et al., 2016; Tolman et al., 2019). The condition (17) and the TAE parameters can be used to find $\Phi=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.

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

$$\partial_{t}f_{1}+\left(v_{\parallel}\hat{b}+\vec{v}_{d}\right)\cdot\nabla f_{1}+\left[\left(v_{\parallel}\hat{b}_{tot}+\vec{v}_{d,tot}\right)\cdot\nabla\psi\right]_{1}\frac{\partial f_{0}}{\partial\psi}+\left[\frac{Z_{\alpha}e}{M_{\alpha}}\left(v_{\parallel}\hat{b}+\vec{v}_{d}\right)\cdot\vec{E}_{1}\right]\frac{\partial f_{0}}{\partial\mathcal{E}}\\ =C\left\{f_{1}\right\},$$

(18)

where $f_{1}$ is the distribution function response to the perturbation and $\left[\left(v_{\parallel}\hat{b}_{tot}+\vec{v}_{d,tot}\right)\cdot\nabla\psi\right]_{1}$ is the component of $\left[\left(v_{\parallel}\hat{b}_{tot}+\vec{v}_{d,tot}\right)\cdot\nabla\psi\right]$ 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:⁸⁸8This drift kinetic angular momentum is a constant of the alpha particle’s motion because the left-hand side of (102) operating on it is zero: $\left(v_{\parallel}\hat{b}+\vec{v}_{d}\right)\cdot\nabla\psi_{\star}=\vec{v}_{d}\cdot\nabla\psi-v_{\parallel}\hat{b}\cdot\nabla\left(Iv_{\parallel}/\Omega\right)-\vec{v}_{d}\cdot\nabla\left(Iv_{\parallel}/\Omega\right)$. It can be shown (Helander & Sigmar, 2005; Parra & Catto, 2010) that $\vec{v}_{d}\cdot\nabla\psi=v_{\parallel}\hat{b}\cdot\nabla\left(Iv_{\parallel}/\Omega\right)$ and $\vec{v}_{d}\cdot\nabla\left(Iv_{\parallel}/\Omega\right)=0$, such that $\left(v_{\parallel}\hat{b}+\vec{v}_{d}\right)\cdot\nabla\psi_{\star}=0$.

$$\psi_{\star}\equiv\psi-\frac{Iv_{\parallel}}{\Omega}.$$

(19)

The coordinate $\alpha$ is modified from (1) in accordance:

$$\alpha_{\star}\equiv\zeta-q_{\star}\vartheta,$$

(20)

with

$$q_{\star}=q\left(\psi_{\star}\right)$$

(21)

the safety factor evaluated at $\psi_{\star}$ instead of $\psi$. Consistent with the approximation (111), 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+\left(Z_{\alpha}e\Phi/M_{\alpha}\right)\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 $\left[\left(v_{\parallel}\hat{b}_{tot}+\vec{v}_{d,tot}\right)\cdot\nabla\psi\right]_{1}$ gives that the equation governing $h$ is

$$\partial_{t}h+\left(v_{\parallel}\hat{b}+\vec{v}_{d}\right)\cdot\nabla h-C\left\{h+\frac{Z_{\alpha}e\Phi}{M_{\alpha}}\frac{\partial f_{0}}{\partial\mathcal{E}}\right\}=\\ -\frac{Z_{\alpha}e}{M_{\alpha}}\left(\frac{\partial\Phi}{\partial t}-\frac{v_{\parallel}}{c}\frac{\partial A_{\parallel}}{\partial t}+\frac{\lambda v^{2}}{2c\Omega}\frac{\partial B_{1\parallel}}{\partial t}\right)\frac{\partial f_{0}}{\partial\mathcal{E}}-c\left(\frac{\partial\Phi}{\partial\alpha_{\star}}-\frac{v_{\parallel}}{c}\frac{\partial A_{\parallel}}{\partial\alpha_{\star}}+\frac{\lambda v^{2}}{2c\Omega}\frac{\partial B_{1\parallel}}{\partial\alpha_{\star}}\right)\frac{\partial f_{0}}{\partial\psi_{\star}}.$$

(22)

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

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

$$\Phi\left(\psi,\vartheta,\alpha\right)=\sum_{m,n,\omega}\Phi_{mn\omega}\cos{\left[n\alpha-\omega t+\left(nq-m\right)\vartheta+\int^{\psi}d\psi^{\prime}\frac{k_{\psi}}{RB_{p}}\right]}\\ =\sum_{m,n,\omega}\Re{\left[\Phi_{mn\omega}e^{in\alpha-i\omega t+i\left(nq-m\right)\vartheta+i\int^{\psi}d\psi^{\prime}\frac{k_{\psi}}{RB_{p}}}\right]}.$$

(23)

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

$$\Phi\left(\psi,\vartheta,\alpha\right)=\sum_{m,n,\omega}\Re{\left[\Phi_{mn\omega}e^{in\alpha_{\star}-i\omega t+i\left(nq_{\star}-m\right)\vartheta+i\int^{\psi_{\star}}d\psi^{\prime}\frac{k_{\psi}}{RB_{p}}+\frac{ik_{\psi}Iv_{\parallel}}{RB_{p}\Omega}}\right]}.$$

(24)

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

$$\partial_{t}h+v_{\parallel}\hat{b}\cdot\nabla\vartheta\frac{\partial h}{\partial\vartheta}+\omega_{\alpha_{\star}}\frac{\partial h}{\partial\alpha_{\star}}\\ =C\left\{h+\frac{Z_{\alpha}e\Phi}{M_{\alpha}}\frac{\partial f_{0}}{\partial\mathcal{E}}\right\}-i\sum_{m,n}D_{n\omega}S_{mn\omega}e^{in\alpha_{\star}-i\omega t+i\left(nq_{\star}-m\right)\vartheta+i\int^{\psi_{\star}}d\psi\frac{k_{\psi}}{RB_{p}}+\frac{ik_{\psi}Iv_{\parallel}}{RB_{p}\Omega}},$$

(25)

where⁹⁹9The quantity $\omega_{\alpha_{\star}}$ is found as follows: $\left(v_{\parallel}\hat{b}+\vec{v}_{d}\right)\cdot\nabla\alpha_{\star}=\left(v_{\parallel}\hat{b}+\vec{v}_{d}\right)\cdot\left(\nabla\zeta-q_{\star}\nabla\vartheta\right)$ $=\left(v_{\parallel}\hat{b}+\vec{v}_{d}\right)\cdot\left[\nabla\zeta-q\nabla\vartheta+\left(Iv_{\parallel}/\Omega\right)\left(\partial q/\partial\psi\right)\nabla\vartheta\right]$ $=\left(v_{\parallel}\hat{b}+\vec{v}_{d}\right)\cdot\left[\nabla\alpha+\vartheta\nabla q+\left(Iv_{\parallel}/\Omega\right)\left(\partial q/\partial\psi\right)\nabla\vartheta\right]$ $\approx\vec{v}_{d}\cdot\nabla\alpha+v_{\parallel}\hat{b}\cdot\nabla\left[\vartheta\left(Iv_{\parallel}/\Omega\right)\left(\partial q/\partial\psi\right)\right]$. In the first equality, we have used that $\left(v_{\parallel}\hat{b}+\vec{v}_{d}\right)\cdot\nabla\psi_{\star}=0$; in the last equality, we have used that $v_{d}\cdot\nabla\psi=v_{\parallel}\hat{b}\cdot\left(Iv_{\parallel}/\Omega\right)$. More information about the evaluation of $\vec{v}_{d}\cdot\nabla\alpha$ can be found in Catto (2019b).

$$\omega_{\alpha_{\star}}\equiv\vec{v}_{d}\cdot\nabla\alpha+v_{\parallel}\hat{b}\cdot\nabla\left(\vartheta\frac{Iv_{\parallel}}{\Omega}\frac{\partial q}{\partial\psi}\right)=v_{\parallel}\hat{b}\cdot\nabla\vartheta\frac{\partial}{\partial\psi}\left(\frac{Bv_{\parallel}}{\Omega\hat{b}\cdot\nabla\vartheta}\right)$$

(26)

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

$$S_{mn\omega}\left(\vartheta,v,\lambda\right)=\Phi_{mn\omega}-\frac{v_{\parallel}}{c}A_{\parallel mn\omega}+\frac{\lambda v^{2}}{2c\Omega}B_{1\parallel mn\omega}.$$

(27)

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

$$D_{n\omega}\left(v,\psi_{\star}\right)\equiv cn\frac{\partial f_{0}}{\partial\psi_{\star}}\left(1-\frac{\omega}{n\omega_{\star}}\right)$$

(28)

with

$$\omega_{\star}\equiv\frac{cM_{\alpha}\partial f_{0}/\partial\psi_{\star}}{Z_{\alpha}e\partial f_{0}/\partial\mathcal{E}}.$$

(29)

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

$$h\left(\psi_{\star},\vartheta,\alpha_{\star},t,v,\lambda,\sigma\right)=\sum_{n,\omega}h_{n\omega}\left(\vartheta,v,\lambda,\sigma\right)e^{in\alpha_{\star}-i\omega t+i\int^{\psi_{\star}}d\psi\frac{k_{\psi}}{RB_{p}}}.$$

(30)

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

$$\sigma=\left\{\begin{array}[]{@{}l@{\thinspace}l}0,\,\text{trapped particles}\\ \frac{v_{\parallel}}{\left|v_{\parallel}\right|},\,\text{passing particles}\\ \end{array}\right..$$

(31)

At this point, a discussion of the role of collisionality in the transport reveals the appropriate treatment of the collision operator in (25). This role is directly analogous to that played by collisionality in superbanana plateau tokamak transport [see discussion in Shaing (2015); Calvo et al. (2017); Catto (2019a)]. It is also similar to the role played by collisionality in the plateau regime of neoclassical transport [see Helander & Sigmar (2005)], in the damping of plasma echos [see Su & Oberman (1968)], and in other wave-particle resonance processes (Duarte et al., 2019; Catto, 2020).

Our analysis of the transport equation (25) 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 side of (25) vanish. If there were no collisionality, the effect of the perturbation drift [represented by the second term on the left side of (25)], on such resonant particles would diverge. So collisionality is necessary to resolve the singularity of the resonance.¹⁰¹⁰10Here, we refer only to behavior existing within the equation we solve: a linearized drift kinetic equation with a perturbation of zero or real frequency. Physically, it is possible for other effects to resolve the resonance. For instance, resonant particles will dephase from the perturbation as they are moved by it, resulting in nonlinear resolution of the resonance through the formation of phase space islands. The addition of a small imaginary component to the perturbation frequency is also able to resolve the resonance. However, it is not clear that phase space island formation (in the absence of resonance overlap) or an imaginary component of the frequency provide the decorellation and dissipation necessary for transport. The relationship between these mechanisms and pitch angle collisions in resolving the resonance and in causing transport is an interesting avenue for future work. One possibility is that resonance resolution by these mechanisms is similar enough to the collisional resolution that fluxes calculated in this paper are valid in cases where collisional resolution is less important than other mechanisms. 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 (9), rather than drag, is the important collisional behavior 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 picture 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, 2005).

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

$$Q_{\alpha}\sim M_{\alpha}v_{res}^{2}\delta\lambda\left\{\left[\left(v_{\parallel}\hat{b}_{tot}+\vec{v}_{d,tot}\right)\cdot\nabla\psi\right]_{1}\Delta t\right\}^{2}\Delta t^{-1}\frac{\partial n_{\alpha}}{\partial\psi},$$

(32)

with $M_{\alpha}$ the alpha particle mass, $v_{res}$ a characteristic alpha resonant velocity, $\delta\lambda$ the fraction of particles that are resonant with the mode, $\left[\left(v_{\parallel}\hat{b}_{tot}+\vec{v}_{d,tot}\right)\cdot\nabla\psi\right]_{1}$ the radial velocity caused by the perturbation, $\Delta t$ the time that a particle moves before becoming decorrelated by a collision (such that $\left[\left(v_{\parallel}\hat{b}_{tot}+\vec{v}_{d,tot}\right)\cdot\nabla\psi\right]_{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 side of (25) (we use $\omega_{\alpha_{\star}}$ to include pitch angle dependence) with the pitch angle scattering part of the collision operator, (9),

$$n\omega_{\alpha_{\star}}\delta\lambda\sim\frac{\nu_{pas}}{\delta\lambda^{2}},$$

(33)

with $\nu_{pas}$ representing the coefficient of the second term of (9). [Note that the effect of the collision operator on the adiabatic component, $\left(Z_{\alpha}e\Phi/M_{\alpha}\right)\partial f_{0}/\partial\mathcal{E}$, is negligible in comparison to the sharp variation of $h$ in the resonant boundary layer.] Thus, we have

$$\delta\lambda\sim\left(\frac{\nu_{pas}}{n\omega_{\alpha_{\star}}}\right)^{1/3}.$$

(34)

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

$$Q_{\alpha}\sim\frac{M_{\alpha}v_{res}^{2}\left[\left(v_{\parallel}\hat{b}_{tot}+\vec{v}_{d,tot}\right)\cdot\nabla\psi\right]_{1}^{2}}{n\omega_{\alpha_{\star}}}\frac{\partial n_{\alpha}}{\partial\psi}.$$

(35)

Note that this flux is independent of the collision frequency $\nu_{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 et al., 1954; Mynick, 1986), that is, an operator for which $C\left\{h\right\}=-\nu h$. The Krook operator $\nu$ is related to $\nu_{pas}$ by $\nu\sim\nu_{pas}/\delta\lambda^{2}$ [recall that the pitch angle scattering operator (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 which are not qualitatively sensitive to the full form of the collision operator.

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

$$v_{\parallel}\hat{b}\cdot\nabla\vartheta\frac{\partial h_{n\omega}}{\partial\vartheta}-i\left(\omega-n\omega_{\alpha_{\star}}+i\nu\right)h_{n\omega}=-iD_{n\omega}\sum_{m}S_{mn\omega}e^{i\left(nq_{\star}-m\right)\vartheta+\frac{ik_{\psi}Iv_{\parallel}}{RB_{p}\Omega}}.$$

(36)

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

$$\Lambda\left(\tau\right)\equiv\omega-n\omega_{\alpha_{\star}}-v_{\parallel}\hat{b}\cdot\nabla\left[\left(nq_{\star}-m\right)\vartheta+\frac{k_{\psi}Iv_{\parallel}}{RB_{p}\Omega}\right],$$

(37)

(36) may be written

$$\frac{\partial}{\partial\tau}\left[h_{n\omega}e^{-i\int_{\tau_{0}}^{\tau}d\tau^{\prime}\left(\omega-n\omega_{\alpha_{\star}}+i\nu\right)}\right]=-iD_{n\omega}\sum_{m}S_{mn\omega}e^{-i\int_{\tau_{0}}^{\tau}d\tau^{\prime}\left[\Lambda\left(\tau^{\prime}\right)+i\nu\right]}e^{\frac{ik_{\psi}Iv_{\parallel}\left(\vartheta=0\right)}{RB_{p}\Omega}}.$$

(38)

Here, we have chosen $\tau_{0}$ such that $\vartheta\left(\tau_{0}\right)=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\left(\tau=0\right)=\vartheta$,

$$\left.h_{n\omega}e^{-i\int_{\tau_{0}}^{\tau}d\tau^{\prime}\left(\omega-n\omega_{\alpha_{\star}}+i\nu\right)}\right|_{\tau=0}=-iD_{n\omega}e^{\frac{ik_{\psi}Iv_{\parallel}\left(\vartheta=0\right)}{RB_{p}\Omega}}\sum_{m}\int_{-\infty}^{0}S_{mn\omega}e^{-i\int_{\tau_{0}}^{\tau}d\tau^{\prime}\left[\Lambda\left(\tau^{\prime}\right)+i\nu\right]}d\tau.$$

(39)

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

$$h_{n\omega}=-iD_{n\omega}\sum_{m}e^{i\left(nq_{\star}-m\right)\vartheta+\frac{ik_{\psi}Iv_{\parallel}}{RB_{p}\Omega}}\int_{-\infty}^{0}S_{mn\omega}e^{-i\int_{0}^{\tau}d\tau^{\prime}\left[\Lambda\left(\tau^{\prime}\right)+i\nu\right]}d\tau.$$

(40)

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 which can be more easily evaluated and understood. The expression (40) can be split into contributions from each complete bounce (for a trapped particle) or poloidal transit (for a passing particle) of duration given by

$$\tau_{b}\equiv\oint d\tau,$$

(41)

such that

$$h_{n\omega}=-iD_{n\omega}\sum_{m}e^{i\left(nq_{\star}-m\right)\vartheta+\frac{ik_{\psi}Iv_{\parallel}}{RB_{p}\Omega}}\sum_{j=0}^{\infty}\left[\int_{-\left(j+1\right)\tau_{b}}^{-j\tau_{b}}\right]S_{mn\omega}e^{-i\int_{0}^{\tau}d\tau^{\prime}\left[\Lambda\left(\tau^{\prime}\right)+i\nu\right]}d\tau.$$

(42)

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

$$h_{n\omega}=-iD_{n\omega}\sum_{m}e^{i\left(nq_{\star}-m\right)\vartheta+\frac{ik_{\psi}Iv_{\parallel}}{RB_{p}\Omega}}\sum_{j=0}^{\infty}e^{ij\oint d\tau\left[\Lambda\left(\tau\right)+i\nu\right]}\int_{-\tau_{b}}^{0}S_{mn\omega}e^{-i\int_{0}^{\tau}d\tau^{\prime}\left[\Lambda\left(\tau^{\prime}\right)+i\nu\right]}d\tau.$$

(43)

This transformation requires that $\oint d\tau\Lambda\left(\tau\right)$ be the same (up to a factor of $2\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 fulfill this condition, because $\oint d\tau\Lambda\left(\tau\right)$ 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, (40) 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 will reveal in  (46) 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 (43) is convergent and may be summed as

$$h_{n\omega}=\sum_{m}\frac{-iD_{n\omega}e^{i\left(nq_{\star}-m\right)\vartheta+\frac{ik_{\psi}Iv_{\parallel}}{RB_{p}\Omega}}\int_{-\tau_{b}}^{0}S_{mn\omega}e^{-i\int_{0}^{\tau}d\tau^{\prime}\left[\Lambda\left(\tau^{\prime}\right)+i\nu\right]}d\tau}{1-e^{i\oint d\tau\left[\Lambda\left(\tau\right)+i\nu\right]}}.$$

(44)

Expanding the denominator about its zeros using $e^{i2\pi l}=1$, with $l$ any integer, denoting the bounce or poloidal transit harmonic, gives

$$h_{n\omega}=\sum_{m,l}\frac{D_{n\omega}e^{i\left(nq_{\star}-m\right)\vartheta+\frac{ik_{\psi}Iv_{\parallel}}{RB_{p}\Omega}}\int_{-\tau_{b}}^{0}S_{mn\omega}e^{-i\int_{0}^{\tau}d\tau^{\prime}\left[\Lambda\left(\tau^{\prime}\right)+i\nu\right]}d\tau}{\oint d\tau\left[\Lambda\left(\tau\right)+i\nu\right]-2\pi l}.$$

(45)

Alpha particles are in general very unlikely to pitch angle scatter because of their high energy. Thus, for nearly all particles, $\nu\oint d\tau\ll\oint\Lambda d\tau-2\pi l$. In this limit, analogous to the limit that occurs in the plateau regime of neoclassical transport [see Helander & Sigmar (2005)], (45) approaches

$$h_{n\omega}=-i\pi D_{n\omega}\sum_{m,l}\delta\left(Q_{l}\right)e^{i\left(nq_{\star}-m\right)\vartheta+\frac{ik_{\psi}Iv_{\parallel}}{RB_{p}\Omega}-i\int_{0}^{\tau_{0}}d\tau^{\prime}\Lambda\left(\tau^{\prime}\right)}\oint S_{mn\omega}e^{-i\int_{\tau_{0}}^{\tau}d\tau^{\prime}\Lambda\left(\tau^{\prime}\right)}d\tau.$$

(46)

Here, we have also split the integral in the exponential into two parts which both have limits at $\tau_{0}$, where $\vartheta\left(\tau_{0}\right)=0$. In addition, we have introduced a delta function $\delta\left(Q_{l}\right)$ with a resonance function as its argument,

$$Q_{l}\left(v,\lambda\right)\equiv\oint\Lambda d\tau-2\pi l.$$

(47)

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.