2. Action-angle variables

We first restrict ourselves to the collisionless Vlasov equation. Since the inclusion of collisions does not affect the fundamental approach, we examine them later in Appendix C. The Vlasov equation is

(2.1)\begin{equation} \frac{\partial f}{\partial t} + \{f, H\} = 0, \end{equation}

where $f$ is the distribution function, $H$ is the single-particle Hamiltonian and $\{\cdot , \cdot \}$ denotes the Poisson bracket. Using phase space coordinates, this can be written as

(2.2)\begin{equation} \frac{\partial f}{\partial t} + \dot{\boldsymbol{q}} \boldsymbol{\cdot} \frac{\partial f}{\partial \boldsymbol{q}} + \dot{\boldsymbol{p}} \boldsymbol{\cdot} \frac{\partial f}{\partial \boldsymbol{p}} = 0, \end{equation}

where for a single particle $\boldsymbol {q}$ is the position, $\boldsymbol {p}$ is the canonical momentum and the time derivatives are given by Hamilton's equations of motion. For electromagnetic fields relevant to a tokamak, the Hamiltonian of a single charged particle is non-trivial. Although this form of the Vlasov equation and others like it offer an intuitive physical picture, these coordinates can make solving the system quite cumbersome. QuaLiKiz instead employs an action-angle formalism to simplify the perturbative analysis. Such a formalism in the context of tokamak physics was first elaborated by Kaufman (Reference Kaufman1972) and expanded upon by Mahajan & Chen (Reference Mahajan and Chen1985). The core principle is to define a canonical transformation,

(2.3)\begin{equation} (\boldsymbol{q}, \boldsymbol{p}) \to (\boldsymbol{\alpha}, \boldsymbol{J}), \end{equation}

for which Hamilton's equations of motion simplify in the new phase space ($\boldsymbol {\alpha }, \boldsymbol {J}$). By restricting ourselves to a canonical transformation, we preserve the form of Vlasov's equation. The coordinates $\boldsymbol {\alpha }$ and $\boldsymbol {J}$ respectively correspond to the action angles and adiabatic invariants of our system. It is well known (Goldstein, Poole & Safko Reference Goldstein, Poole and Safko2001) that Hamilton's equations of motion then reduce to

(2.4)\begin{gather} \frac{\partial H}{\partial \boldsymbol{\alpha}} ={-} \dot{\boldsymbol{J}} = \boldsymbol{0}, \end{gather}

(2.5)\begin{gather}\frac{\partial H}{\partial \boldsymbol{J}} = \dot{\boldsymbol{\alpha}} = \boldsymbol{\varOmega}, \end{gather}

where $\boldsymbol {\varOmega }$ are the constant angular frequencies associated with each adiabatic invariant and the time derivatives are written in Newton's dot notation. At first glance, it may seem that we have simply shifted the difficulty of the problem to calculating this new canonical transformation itself. The power of this method comes from analysing the unperturbed system and then including electromagnetic fluctuations in the Hamiltonian.

We define the Hamiltonian to be

(2.6)\begin{equation} H = H_0 + \delta h, \end{equation}

where the unperturbed Hamiltonian is simply

(2.7)\begin{equation} H_0 = \frac{1}{2m} (\boldsymbol{p}^2 - e \boldsymbol{A}_0^2) + e \varPhi. \end{equation}

Here, $m$ and $e$ are respectively the mass and charge of the particle, $\boldsymbol {A}_0$ is the equilibrium vector potential, and $\varPhi$ is the equilibrium electrostatic potential. Since QuaLiKiz operates in the electrostatic limit, we therefore apply a perturbation $\delta h$ such that

(2.8)\begin{equation} \delta h = e \phi, \end{equation}

where $\phi$ is the electrostatic perturbation. We then define the action-angle coordinates in reference to the unperturbed Hamiltonian,

(2.9)\begin{gather} \frac{\partial H_0}{\partial \boldsymbol{\alpha}} = \boldsymbol{0}, \end{gather}

(2.10)\begin{gather}\frac{\partial H_0}{\partial \boldsymbol{J}} = \boldsymbol{\varOmega}. \end{gather}

Hamilton's equations of motion then become

(2.11)\begin{gather} \frac{\partial H}{\partial \boldsymbol{\alpha}} ={-} \dot{\boldsymbol{J}} = e \frac{\partial \phi}{\partial \boldsymbol{\alpha}}, \end{gather}

(2.12)\begin{gather}\frac{\partial H}{\partial \boldsymbol{J}} = \dot{\boldsymbol{\alpha}} = \boldsymbol{\varOmega} + e \frac{\partial \phi}{\partial \boldsymbol{J}}. \end{gather}

We note that because the unperturbed Hamiltonian is a function of $\boldsymbol {J}$ and not $\boldsymbol {\alpha }$, all equilibrium quantities are also only functions of $\boldsymbol {J}$. Furthermore, any function of $\boldsymbol {\alpha }$ is periodic with respect to $\boldsymbol {\alpha }$; thus, the perturbed quantities in our system admit a Fourier series expansion. Moreover, it can be shown that $\boldsymbol {\alpha } \text { and } \boldsymbol {J}$ are canonical coordinates even after introducing a perturbation (Mahajan & Chen Reference Mahajan and Chen1985). These features will simplify the derivation greatly.

The next task is to define the canonical transformation by specifying the action-angle variables in terms of the position $\boldsymbol {r}$ and the velocity $\boldsymbol {v}$ of the particle. The three adiabatic invariants in a tokamak correspond to the magnetic moment, the longitudinal invariant (also known as the bounce-transit action) and the poloidal flux. They are defined as follows:

(2.13)\begin{gather} J_1 = \frac{m}{e} \mu, \end{gather}

(2.14)\begin{gather}J_2 = \frac{1}{2 {\rm \pi}} \oint m p_\parallel\, \textrm{d} l = \frac{1}{2 {\rm \pi}} \oint (m v_\parallel{+} e A_\parallel)\, \textrm{d} l, \end{gather}

(2.15)\begin{gather}J_3 = P_{\varphi} = \frac{m v_\parallel R_0 B_\varphi^0}{B} - e \psi. \end{gather}

Here, $\mu = W_\perp /B$ is the magnetic moment, where $W_\perp = \frac {1}{2} m v_\perp ^2$ is the kinetic energy associated with the velocity perpendicular to the magnetic field $\boldsymbol {B}$. Meanwhile, $v_\parallel$ and $A_\parallel$ are the components of the velocity and vector potential parallel to the magnetic field, respectively, with $dl$ being the signed differential length along the particle orbit. We also define $\psi$ to be minus the poloidal magnetic flux normalized to $2 {\rm \pi}$, which is calculated by integrating the flux of the magnetic field through a disk tangent to the flux surface everywhere:

(2.16)\begin{equation} \psi ={-} \frac{1}{2{\rm \pi}} \int_S \boldsymbol{B} \boldsymbol{\cdot} \, \textrm{d} \boldsymbol{S}. \end{equation}

Subsections 2.1–2.3 discuss each of the three adiabatic invariants and define their associated action angles and angular frequencies. For the remainder of the derivation, we also use the spatial coordinates $\boldsymbol {r} = (r, \theta , \varphi )$, where $r$ is the minor radial position, $\theta$ is the geometric poloidal angle and $\varphi$ is the geometric toroidal angle. We use a right-handed coordinate system such that $\hat {\boldsymbol {r}} \times \hat {\boldsymbol {\theta }} = \hat {\boldsymbol {\varphi }}$. For further references characterizing the action angles $\boldsymbol {\alpha }$, we refer the reader to the works of Garbet (Reference Garbet2001) and Garbet et al. (Reference Garbet, Laurent, Mourgues, Roubin and Samain1990).

2.1. Magnetic moment

In the presence of a magnetic field, charged particles gyrate about the field line at the cyclotron frequency $\varOmega _1 = eB/m$. With a strong enough magnetic field, the cyclotron frequency is much larger than any other characteristic frequency in the plasma. Under such conditions, the magnetic moment $\mu$ is adiabatically conserved (Brizard & Hahm Reference Brizard and Hahm2007; Stephens et al. Reference Stephens, Brzozowski III and Jenko2017), and the gyromotion can be decoupled from the guiding centre motion:

(2.17)\begin{gather} r = r_G + \rho \cos(\alpha_1), \end{gather}

(2.18)\begin{gather}\theta = \theta_G + \frac{\rho}{r} \sin(\alpha_1), \end{gather}

(2.19)\begin{gather}\varphi = \varphi_G, \end{gather}

where $\alpha _1$ is equivalent to the gyrophase, $\rho$ is the gyroradius and the subscript ‘$G$’ refers to the location of the particle's guiding centre. These guiding centre variables obey the guiding centre equations of motion. Ordinarily, the exact invariant associated with the gyromotion depends on the electrostatic potential. For QuaLiKiz, we assume that the electrostatic field is small compared with the kinetic energy. Thus, we simply take $J_1$ to be the ordinary magnetic moment $\mu = W_\perp /B$.

Later in the derivation, we will need to average various functions over the gyrophase $\alpha _1$ by integrating over $\alpha _1$. We therefore consider the general integral

(2.20)\begin{equation} g_{n_1} = \int_{-{\rm \pi}}^{\rm \pi} \frac{\textrm{d} \alpha_1}{2 {\rm \pi}} g(\boldsymbol{r}) \, \textrm{e}^{- {\rm i} n_1 \alpha_1}, \end{equation}

where $n_1$ is an integer. It will be shown later that factors of $\textrm {e}^{- \textrm {i} n_1 \alpha _1}$ arise from taking Fourier expansions in terms of $\alpha _1$. We define the Fourier transform of $g$ to be

(2.21)\begin{equation} \tilde{g}(\boldsymbol{k}) = \int \, \textrm{d}^3r g(\boldsymbol{r}) \, \textrm{e}^{{\rm i} \boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{r}}, \end{equation}

with the corresponding inverse Fourier transform

(2.22)\begin{equation} g(\boldsymbol{r}) = \int \frac{\textrm{d}^3k}{(2{\rm \pi})^3} \tilde{g}(\boldsymbol{k}) \, \textrm{e}^{-{\rm i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{r}}. \end{equation}

We use the Fourier transform to obtain

(2.23)\begin{align} g_{n_1} = \int_{-{\rm \pi}}^{\rm \pi} \frac{\textrm{d} \alpha_1}{2 {\rm \pi}} \int \frac{\textrm{d}^3k}{(2{\rm \pi})^3}\tilde{g}(\boldsymbol{k})\, \textrm{e}^{-{\rm i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{r}- {\rm i} n_1 \alpha_1} = \int_{-{\rm \pi}}^{\rm \pi} \frac{\textrm{d} \alpha_1}{2 {\rm \pi}} \int \frac{\textrm{d}^3k}{(2{\rm \pi})^3}\tilde{g}(\boldsymbol{k})\, \textrm{e}^{-{\rm i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{\rho} - {\rm i} n_1 \alpha_1} \, \textrm{e}^{-{\rm i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{R}_G}. \end{align}

Here, we have decoupled the gyromotion from the guiding centre motion via $\boldsymbol {r} =\boldsymbol {R}_G + \boldsymbol {\rho }$. We then write

(2.24)\begin{equation} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{\rho} = k_\perp \rho \cos(\alpha_1), \end{equation}

where

(2.25)\begin{equation} k_\perp{=} |\boldsymbol{k} - \boldsymbol{k} \boldsymbol{\cdot} \hat{\boldsymbol{b}}| \approx \sqrt{k_r^2 + k_\theta^2}. \end{equation}

Because the toroidal magnetic field is much stronger than the poloidal field, we neglect $k_{\varphi }$ terms when computing $k_\perp$ and approximate it as above. Note that according to our definition of the Fourier transform, $k_r$ and $k_\theta$ are operators in real space such that

(2.26)\begin{gather} k_r \to {\rm i} \frac{\partial}{\partial r}, \end{gather}

(2.27)\begin{gather}k_\theta \to \frac{{\rm i}}{r} \frac{\partial}{\partial \theta}. \end{gather}

We may then integrate over $\alpha _1$ independently, which leads to

(2.28)\begin{equation} \int^{\rm \pi} _{-{\rm \pi}} \frac{\textrm{d} \alpha_1}{2 {\rm \pi}} \, \textrm{e}^{- {\rm i} k_\perp \rho \cos(\alpha_1) - {\rm i} n_1 \alpha_1} = ({-}{\rm i})^{n_1} J_{n_1}(k_\perp \rho), \end{equation}

where $J_n$ is the $n$th Bessel function of the first kind. Therefore, we finally have that

(2.29)\begin{equation} g_{n_1} = \int \frac{\textrm{d}^3 k}{(2 {\rm \pi})^3} ({-}{\rm i})^{n_1} J_{n_1}(k_\perp \rho) \tilde{g}(\boldsymbol{k}) \, \textrm{e}^{- {\rm i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{R}_G} = ({-}{\rm i})^{n_1} (J_{n_1}(k_\perp \rho) \boldsymbol{\cdot} g)(\boldsymbol{R}_G).\end{equation}

As a shorthand, we treat the Bessel function in real space as a differential operator that acts on $g$, after which we evaluate the resulting function at the guiding centre. The Bessel function is simply a scalar function in Fourier space instead of a differential operator. The case of $n_1 = 0$ corresponds to the well-known gyro-average. After completing the gyro-average, all functions are evaluated at the guiding centre. Thus, we drop the subscript ‘G’ for convenience and treat all spatial variables as those corresponding to the guiding centre. The adiabatic invariants $J_2$ and $J_3$ are explicitly calculated within the guiding centre framework where we hold $\mu$ constant and ignore the cyclotron motion.

2.2. Longitudinal invariant

To calculate $J_2$, we consider the guiding centre particle motion along a magnetic field line; such a particle completes bounce-transit orbits with frequency $\varOmega _2$. This is the bounce frequency for trapped particles and the transit frequency for passing particles. Here, we neglect excursions from the field line due to various guiding centre drifts by holding $r$ constant. For an extended treatment of bounce-transit motion, see the works of Brizard (Reference Brizard2011) and Stephens, Garbet & Jenko (Reference Stephens, Garbet and Jenko2020).

Assuming that the equilibrium electrostatic potential is small, the guiding centre velocity parallel to the magnetic field is

(2.30)\begin{equation} |v_\parallel| = \sqrt{\frac{2}{m}(E - \mu B)}, \end{equation}

where $E$ is the total kinetic energy of the particle. As an approximation, we take the typical equilibrium magnetic field to be of the form

(2.31)\begin{equation} \boldsymbol{B} = B_\varphi(r, \theta) {\hat{\boldsymbol{\varphi}}} + B_\theta(r, \theta) {\hat{\boldsymbol{\theta}}} = \frac{1}{1 + r/R_0 \cos(\theta)}( B^0_\varphi(r) {\hat{\boldsymbol{\varphi}}} + B^0_\theta(r) {\hat{\boldsymbol{\theta}}}), \end{equation}

where $R_0$ is the major radius. This corresponds to the magnetic field in a circular cross-section tokamak without any Shafranov shift. Defining the inverse aspect ratio $\epsilon = r/R_0$, we recognize that this circular equilibrium is the small $\epsilon$ limit of a more general axisymmetric equilibrium. QuaLiKiz is thus well suited to machines where the aspect ratio of the device is $\sim 3$ or larger. Devices with smaller aspect ratios such as spherical tokamaks, however, cannot be reliably simulated in QuaLiKiz.

A particle is considered trapped if it reflects at a bounce angle $\theta _b$, which requires

(2.32)\begin{equation} \frac{\mu B^0 (r)}{E} \ge 1 - \epsilon. \end{equation}

Otherwise, the particle is considered passing since it will simply continue to travel along the magnetic field line without reflecting. We rewrite $v_\parallel$ to be

(2.33)\begin{equation} v_\parallel{=} \sqrt{\frac{2 T}{m}} \epsilon_\parallel \sqrt{\xi} \sqrt{1 - \lambda b(r, \theta)}. \end{equation}

Here, $\xi = E/T$ where $T$ is the temperature and $\epsilon _\parallel = \pm 1$ determines the sign of the parallel velocity. We also define

(2.34)\begin{gather} \lambda = \frac{\mu B(r, \theta = 0)}{E} = \frac{v_\perp^2}{v^2 b(r, \theta)}, \end{gather}

(2.35)\begin{gather}b(r, \theta) = \frac{B(r, \theta)}{B(r, \theta = 0)}. \end{gather}

It is clear then that $\lambda$ is a pitch angle parameter and determines whether the particle is trapped or passing.

The bounce-transit frequency is defined as

(2.36)\begin{equation} |\varOmega_2| = \frac{2 {\rm \pi}}{T_2}, \end{equation}

where

(2.37)\begin{equation} T_2 = \oint \frac{\textrm{d} \theta}{\displaystyle \left|\frac{\textrm{d} \theta}{\textrm{d} t}\right|}. \end{equation}

We note that one full poloidal orbit for trapped particles includes both the forward motion, where $\theta$ goes from $-\theta _b$ to $\theta _b$, and the backward motion, where $\theta$ goes from $\theta _b$ to $-\theta _b$. For passing particles, the poloidal orbit only includes one full pass where $\theta$ goes from $-{\rm \pi}$ to ${\rm \pi}$. The sign of the transit frequency for passing particles is aligned with that of the parallel velocity and is thus determined by $\epsilon _\parallel$, while the bounce frequency is always positive for trapped particles. Assuming that $B_\varphi \gg B_\theta$, then $\hat {\boldsymbol {b}}$, the direction of the magnetic field, is approximately ${\hat {\boldsymbol {\varphi }}}$. We again emphasize that this approximation breaks down for devices such as spherical tokamaks. Therefore, we write

(2.38)\begin{equation} v_\parallel{=} \boldsymbol{v}\boldsymbol{\cdot}\hat{\boldsymbol{b}} \approx \dot{\varphi} (R_0 + r \cos(\theta)) \approx q R_0 \dot{\theta}, \end{equation}

where we have taken the toroidal field to lie mostly in the toroidal direction to make the above approximation. Here, we have defined the safety factor as

(2.39)\begin{equation} q(r) = \frac{\textrm{d} \varphi}{\textrm{d} \theta} \approx \frac{r B_\varphi}{R_0 B_\theta}. \end{equation}

The safety factor describes how many times a magnetic field line wraps around toroidally per poloidal turn. The magnitude of the bounce-transit frequency is then

(2.40)\begin{equation} |\varOmega_2| = \sqrt{\frac{2T}{m}} \frac{\sqrt{\xi}}{qR_0} \bar{\varOmega}_2 (r, \lambda), \end{equation}

where we define

(2.41)\begin{equation} \bar{\varOmega}_2 (r, \lambda) = \frac{2 {\rm \pi}}{\displaystyle \oint \, \textrm{d} \theta \frac{1}{\sqrt{1 - \lambda b(r, \theta)}}}. \end{equation}

Note that for passing particles, we take the sign of the transit frequency to be the sign of the parallel velocity and multiply by $\epsilon _\parallel$ accordingly. We then calculate $\bar {\varOmega }_2$ in the small $\epsilon$ limit to find

(2.42)\begin{equation} \bar{\varOmega}_2 \approx \left\{\begin{array}{@{}ll} \displaystyle \dfrac{{\rm \pi} \sqrt{\epsilon}}{2 \sqrt{2} {K}(\kappa)} & \text{if } 0\le \kappa < 1 \text{ (trapped)},\\ \displaystyle \dfrac{{\rm \pi} \kappa \sqrt{\epsilon}}{\sqrt{2} {K}(\kappa^{{-}1})} & \text{if } 1 < \kappa < \infty \text{ (passing)}. \end{array}\right. \end{equation}

Here, ${K}$ is the complete elliptic integral of the first kind and $\kappa$ is a trapped parameter defined such that

(2.43)\begin{equation} \lambda = 1 - 2 \epsilon \kappa^2. \end{equation}

In the small $\epsilon$ limit, $0 \leq \kappa < 1$ for trapped particles and $1 < \kappa < \infty$ for passing particles. We also calculate the bounce-transit action to be

(2.44)\begin{equation} J_2 \approx \left\{\begin{array}{@{}ll} \displaystyle \dfrac{8 m q R_0 \sqrt{E/m} \sqrt{\epsilon}}{\rm \pi} ({E}(\kappa) - (1-\kappa^2) {K}(\kappa)) & \text{if } 0\le \kappa < 1 \text{ (trapped)},\\ \displaystyle \dfrac{4 m q R_0 \sqrt{E/m} \sqrt{\epsilon}}{\rm \pi} \kappa {E}(\kappa^{{-}1}) + e \varPhi_t & \text{if } 1 < \kappa < \infty \text{ (passing)}, \end{array}\right. \end{equation}

where ${E}$ is the complete elliptic integral of the second kind and $\varPhi _t$ is the toroidal flux normalized by $2 {\rm \pi}$. The flux term is absent for trapped particles since the closed line integral of $A_\parallel$ is zero for trapped orbits.

Calculating the angular variable $\alpha _2$ requires the explicit equation of motion

(2.45)\begin{equation} \frac{\textrm{d} \alpha_2}{\textrm{d} t} = \varOmega_2. \end{equation}

This is of course the definition of $\alpha _2$ such that it is conjugate to the action variable $J_2$. To find an explicit expression for $\alpha _2$ in terms of the poloidal angle $\theta$, we make use of the chain rule,

(2.46)\begin{equation} \frac{\textrm{d} \alpha_2}{\textrm{d} t} = \frac{\textrm{d} \alpha_2}{\textrm{d} \theta} \dot{\theta} = \frac{\textrm{d} \alpha_2}{\textrm{d} \theta} \frac{v_\parallel}{q R_0} = \varOmega_2. \end{equation}

We emphasize that $\varOmega _2$ is not dependent on $\alpha _2$ or $\theta$. Thus, this differential equation can be integrated using elliptic functions, which leads to an expression of $\alpha _2$ in terms of $\theta$. We use the convention that $\alpha _2(\theta = 0) = 0$, which leads to

(2.47)\begin{equation} \alpha_2 = \int^{\theta}_{0} \, \textrm{d} \theta' \frac{q R_0 \varOmega_2}{v_{{\parallel}}} = \int^{\theta} \, \textrm{d} \theta' \frac{q R_0 \varOmega_2}{\sqrt{2T / m} \epsilon_{{\parallel}} \sqrt{\xi} \sqrt{1 - \lambda b}}. \end{equation}

For trapped particles, we must keep in mind that $\epsilon _{\parallel }$ switches sign after the particle bounces. The integral can then be simplified in the small $\epsilon$ limit, which leads to

(2.48)\begin{equation} \alpha_2 \approx \left\{ \begin{array}{@{}ll} \dfrac{\rm \pi}{2} \dfrac{{F}\left(\sin^{{-}1}\left[\kappa \sin\left(\dfrac{\theta}{2}\right)\right], \kappa\right)}{{K}(\kappa)} & \text{if } 0\le \kappa < 1 \text{ (trapped)},\\ {\rm \pi}\dfrac{{F}\left(\dfrac{\theta}{2}, \kappa^{{-}1}\right)}{{K}(\kappa^{{-}1})} & \text{if } 1 < \kappa < \infty \text{ (passing)}, \end{array} \right. \end{equation}

where ${F}$ is the incomplete elliptic integral of the first kind. Essentially, the integral takes the same form as when calculating $\varOmega _2$, the primary difference being that we integrate up to arbitrary $\theta$ rather than up to the bounce angle $\theta _b$ for trapped particles or up to ${\rm \pi}$ for passing particles.

Finally, let $G (\epsilon _\parallel , \theta )$ be a quantity that varies over the bounce-transit orbit along the field line. It is often of interest to time average $G$ over the orbit; we define the bounce-transit average $\langle G (\epsilon _\parallel , \theta )\rangle$ to be

(2.49)\begin{equation} \langle G (\epsilon_\parallel, \theta )\rangle = \frac{1}{T_2} \oint \frac{\textrm{d} \theta}{\displaystyle \left|\frac{\textrm{d} \theta}{\textrm{d} t}\right|} G (\epsilon_\parallel, \theta ) = \frac{\displaystyle \oint \, \textrm{d} \theta \frac{G (\epsilon_\parallel, \theta )}{\sqrt{1 - \lambda b}}}{\displaystyle \oint \, \textrm{d} \theta \frac{1}{\sqrt{1-\lambda b}}}. \end{equation}

For passing particles, the average is explicitly

(2.50)\begin{equation} \oint \, \textrm{d} \theta \frac{G (\epsilon_\parallel, \theta )}{\sqrt{1-\lambda b}} = \int_{-{\rm \pi}} ^{\rm \pi} \, \textrm{d} \theta \frac{G (\epsilon_\parallel, \theta )}{\sqrt{1-\lambda b}}, \end{equation}

while for trapped particles the average is instead

(2.51)\begin{equation} \oint \, \textrm{d} \theta \frac{G (\epsilon_\parallel, \theta )}{\sqrt{1-\lambda b}} = \int_{-\theta_b} ^{\theta_b} \, \textrm{d} \theta \frac{G(\epsilon_\parallel, \theta) + G(-\epsilon_\parallel, \theta)}{\sqrt{1-\lambda b}}. \end{equation}

Note that because the line integral must be closed, a sum over $\epsilon _\parallel$ must be performed for trapped particles so that quantities such as $v_\parallel$ average to 0.

In this discussion so far, we have neglected any magnetic drifts and excursions from the field line. We include such effects in § 2.3, as they characterize the third adiabatic invariant – the poloidal flux.

2.3. Poloidal flux

In an axisymmetric equilibrium, the canonical toroidal momentum, $P_\varphi$, is conserved since no external quantities depend explicitly on the toroidal angle $\varphi$. From guiding centre theory, we can write the canonical toroidal momentum as

(2.52)\begin{equation} P_\varphi = \frac{m v_\parallel R_0 B_\varphi^0}{B} - e \psi. \end{equation}

This is an exact invariant of the system. We construct $J_3$ such that it approximates $P_\varphi$ provided that the poloidal flux term dominates. For typical parameters in a tokamak plasma this is indeed the case, since $(P_\varphi + e \psi ) / (e \psi ) \sim \sqrt {m T}/(eB R_0)$. Inputting Joint European Torus (JET)-like parameters, $T = 5$ keV, $m = m_D, B = 3$ T, $R_0 = 3$ m, then $\sqrt {m T}/(eB R_0) \sim 10^{-3}$, making this a very reasonable approximation. We therefore write

(2.53)\begin{equation} J_3 ={-} e \psi. \end{equation}

To calculate the poloidal flux, we utilize Stoke's theorem; the surface integral of $\boldsymbol {B}$ simply becomes a closed line integral of $\boldsymbol {A}$ to find

(2.54)\begin{equation} A_\varphi(r, \theta) = \frac{1}{1 + \epsilon \cos(\theta)} \int^r B_\theta(r^\prime, \theta) (1 + r^\prime/R_0 \cos(\theta)) \, \textrm{d} r^\prime = \frac{\displaystyle \int^r B^0_\theta(r^\prime) \, \textrm{d} r^\prime}{1 + \epsilon \cos(\theta)}. \end{equation}

Thus, we obtain

(2.55)\begin{equation} J_3 ={-} e \psi ={-} \frac{1}{2{\rm \pi}} \int_0 ^{2{\rm \pi}} e(1 + \epsilon \cos(\theta)) A_\varphi \, \textrm{d} \varphi = e \int^r R_0 B_\theta^0(r^\prime) \, \textrm{d} r^\prime. \end{equation}

We see then that $J_3$ is purely a function of $r$ such that

(2.56)\begin{equation} \frac{\textrm{d} J_3}{\textrm{d} r} \approx{-}\frac{e r B}{q}. \end{equation}

We next calculate $\varOmega _3$, which is the toroidal precession frequency for trapped particles and the toroidal rotation frequency for passing particles. Along the bounce-transit orbit, guiding centre drifts cause radial excursions from the magnetic field line. In addition, passing particles wind around the magnetic field line toroidally due to the lack of any bounce point. To calculate this frequency, we need to first calculate deviations from the field line orbit, noting that radial excursions from the field line are of the order of the gyroradius. To aid in the calculation, we exploit the exact conservation of the canonical toroidal momentum:

(2.57)\begin{equation} \psi = \bar{\psi} + \frac{m B^0_\varphi R_0 v_\parallel} {e B} = \bar{\psi} + \psi_1. \end{equation}

Here, $\bar {\psi }$ corresponds to the reference magnetic flux surface defined to be

(2.58)\begin{equation} \bar{\psi} ={-} \frac{P_\varphi}{e}, \end{equation}

and $\psi _1$ is the deviation from that flux surface. Recall that the field line orbit assumed that

(2.59)\begin{equation} \frac{\textrm{d} \theta}{\textrm{d} t} \approx \frac{v_\parallel}{q R_0}, \end{equation}

where we hold $r$ and thus $\psi$ fixed. Since the exact field-line-following orbit breaks $P_\varphi$ conservation, we need to include deviations from the field line caused by conservation of $P_\varphi$ along with guiding centre drifts in order to consistently expand the guiding centre equation of motion with respect to the gyroradius.

The guiding centre equation of motion is

(2.60)\begin{equation} \frac{\textrm{d} x}{\textrm{d} t} = (v_\parallel \hat{\boldsymbol{b}} + \boldsymbol{v}_D ) \boldsymbol{\cdot} \boldsymbol{\nabla} x, \end{equation}

where $x$ is any spatial coordinate and $\boldsymbol {v}_D$ is the guiding centre drift velocity. We then expand the guiding centre equation of motion for variables $\psi$, $\theta$ and $\varphi$ and find that

(2.61)\begin{gather} \frac{\textrm{d} \psi}{\textrm{d} t} \approx (\boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla} \psi)|_{\bar{\psi}}, \end{gather}

(2.62)\begin{gather}\frac{\textrm{d} \theta}{\textrm{d} t} \approx \left. \left(\frac{v_\parallel}{q R_0} \right) \right|_{\bar{\psi}} + \frac{\textrm{d}}{\textrm{d} \psi} \left. \left(\frac{v_\parallel}{q R_0} \right) \right|_{\bar{\psi}} \psi_1 + (\boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla} \theta )|_{\bar{\psi}}, \end{gather}

(2.63)\begin{gather}\frac{\textrm{d} \varphi}{\textrm{d} t} \approx \left(\frac{v_\parallel}{q R_0} \frac{\textrm{d} \varphi}{\textrm{d} \theta} \right) + \frac{\textrm{d}}{\textrm{d} \psi} \left. \left(\frac{v_\parallel}{q R_0} \frac{\textrm{d} \varphi}{\textrm{d} \theta} \right) \right|_{\bar{\psi}} \psi_1 + (\boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla} \varphi ) |_{\bar{\psi}}. \end{gather}

We take the $\boldsymbol {v}_D$ to be the the sum of the classical curvature, grad-$B$, and $E$-cross-$B$ drifts:

(2.64)\begin{equation} \boldsymbol{v}_D = \frac{\hat{\boldsymbol{b}}}{e B} (m v_\parallel^2 \boldsymbol{k} + \mu \boldsymbol{\nabla} B + e \boldsymbol{\nabla} \varPhi), \end{equation}

where $\hat {\boldsymbol {k}}$ is the curvature vector defined such that

(2.65)\begin{equation} \hat{\boldsymbol{b}} \times \boldsymbol{k} = \boldsymbol{\nabla} \times \hat{\boldsymbol{b}} - (\hat{\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} \times \hat{\boldsymbol{b}}) \hat{\boldsymbol{b}}. \end{equation}

We can simplify the equation of motion in the toroidal direction by noting that

(2.66)\begin{align} \frac{\textrm{d}}{\textrm{d} \psi} \left(\frac{v_\parallel}{q R_0} \frac{\textrm{d} \varphi}{\textrm{d} \theta} \right) \psi_1 &= \frac{\textrm{d}}{\textrm{d} \psi} \left(\frac{v_\parallel}{q R_0} \right) \frac{\textrm{d} \varphi}{\textrm{d} \theta} \psi_1 + \frac{\textrm{d}}{\textrm{d} \psi} \left( \frac{\textrm{d} \varphi}{\textrm{d} \theta} \right) \frac{v_\parallel}{q R_0} \psi_1 \nonumber\\ &= \frac{\textrm{d} \varphi}{\textrm{d} \theta} \frac{\textrm{d} \theta}{\textrm{d} t} + \frac{\textrm{d}}{\textrm{d} \psi} \left( \frac{\textrm{d} \varphi}{\textrm{d} \theta} \right) \frac{v_\parallel}{q R_0} \psi_1 - \frac{\textrm{d} \varphi}{\textrm{d} \theta} \frac{v_\parallel}{q R_0} - \frac{\textrm{d} \varphi}{\textrm{d} \theta} (\boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla} \theta ), \end{align}

where we used the equation of motion in the poloidal direction. Substituting this in and evaluating ${\textrm {d} \varphi }/{\textrm {d} \theta }$, we obtain

(2.67)\begin{align} \dot{\varphi} = \frac{\textrm{d}}{\textrm{d} \psi} \left( \frac{q}{1 + \epsilon \cos(\theta)}\right) \psi_1 \frac{v_\parallel}{q R_0} + \boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla} \varphi - \frac{q}{1 + \epsilon \cos(\theta)} \boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla} \theta + \frac{q}{1 + \epsilon \cos(\theta)} \frac{\textrm{d} \theta}{\textrm{d} t}, \end{align}

where we evaluate all radial coordinates at $\bar {r}$ such that

(2.68)\begin{equation} \psi(\bar{r}) = \bar{\psi}. \end{equation}

Finally, we take the bounce-transit average of $\dot {\varphi }$ and find that

(2.69)\begin{equation} \varOmega_3 = \langle \dot{\varphi}\rangle = \langle \omega_d\rangle + \bar{\epsilon} q(\bar{r}) \varOmega_2 = \varOmega_d + \bar{\epsilon} q(\bar{r}) \varOmega_2. \end{equation}

Here, $\bar {\epsilon }$ is 0 for trapped particles and 1 for passing particles, $\varOmega _d$ is the frequency purely due to the guiding centre drifts, and $\omega _d$ is associated with the instantaneous deviation from the magnetic field line. The extra term for passing particles is due to the toroidal rotation from following the field line in a complete poloidal turn. This parallel velocity dependent term is absent for trapped particles since their average toroidal position does not change as a result of a complete field-line-following bounce. We approximate $\omega _d$ as

(2.70)\begin{equation} \omega_d \approx \frac{\textrm{d} q}{\textrm{d} \psi} \frac{\psi_1 v_\parallel}{q} + \boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla} \varphi - q \boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla} \theta. \end{equation}

The poloidal component of the magnetic drift dominates, thus we ignore the toroidal component. Using the $s\text {--}\alpha$ equilibrium, we calculate the guiding centre drift in Appendix B. Thus, $\varOmega _3$ is computed with a finite Shafranov shift. The poloidal component of the guiding centre drift is

(2.71)\begin{equation} \boldsymbol{v}_{D} \boldsymbol{\cdot} \boldsymbol{\nabla} \theta \approx{-}\frac{m}{e B r R_0} \left(v_\parallel^2 + \frac{v_\perp^2}{2} \right) ( \cos(\theta) - \alpha_M \sin^2 (\theta)) - \frac{E_r }{ r B}, \end{equation}

where $E_r$ is the radial electric field. We also define $\alpha _M$ such that

(2.72)\begin{equation} \alpha_M ={-} q^2 \beta \frac{R_0}{P} \frac{\textrm{d} P}{\textrm{d} r}. \end{equation}

Here, $\beta = 2 \mu _0 P/B^2$ with $\mu _0$ being the vacuum permeability. The $E$-cross-$B$ drift can be separated from the magnetic drifts, so that we obtain

(2.73)\begin{equation} \boldsymbol{v}_{D,B} \boldsymbol{\cdot} \boldsymbol{\nabla} \theta \approx{-} \frac{v_{D,B}}{r} (\cos(\theta)- \alpha_M \sin^2(\theta)), \end{equation}

where $v_{D,B}$ characterizes the magnetic drifts and can be written as

(2.74)\begin{equation} v_{D,B} = \frac{m}{eBR_0} \left(v_\parallel^2 + \frac{v_\perp^2}{2}\right) = \frac{\xi T}{e B R_0} (2 - \lambda b), \end{equation}

where the magnetic shear is

(2.75)\begin{equation} s = \frac{r}{q} \frac{\textrm{d} q}{\textrm{d} r}. \end{equation}

We can therefore rewrite $\varOmega _d$ as

(2.76)\begin{align} \varOmega_d = \frac{q}{r} \frac{T \xi}{e B R_0} \left\langle (2 - \lambda b) (\cos(\theta)- \alpha_M \sin^2(\theta)) + \frac{2 s}{\epsilon} (1 - \lambda b )\right\rangle + \left\langle \frac{q E_r}{r B}\right\rangle = \omega_{d0} \xi F_d(\kappa) + \omega_E, \end{align}

where $F = F_d(\kappa )$ is the bounce-transit averaged term, $\omega _{d0}$ is characteristic of the magnetic precession frequency and defined to be

(2.77)\begin{equation} \omega_{d0} = \frac{q T}{r e B R_0}, \end{equation}

and $\omega _E$ corresponds to the $E$-cross-$B$ velocity. We explicitly carry out the bounce-transit average by rewriting the $\lambda b$ terms in terms of $\kappa$ and $\theta$ and find that

(2.78)\begin{align} F_d(\kappa) \approx \left\{ \begin{array}{@{}ll} -1 + \dfrac{2 {E}(\kappa)}{{K}(\kappa)} + 4 s \left(\kappa^2 -1 + \dfrac{{E}(\kappa)}{{K}(\kappa)} \right) & \\ \quad - \dfrac{4 \alpha_M}{3} \left(1 - \kappa^2 - (1 - 2 \kappa^2) \dfrac{{E}(\kappa)}{{K}(\kappa)}\right) & \text{if } 0\le \kappa < 1 \text{ (trapped)},\\ - (2 \kappa^2 -1 )+ 2 \kappa^2 \dfrac{{E}(\kappa^{{-}1})}{{K}(\kappa^{{-}1})} + 4 s \kappa^2 \dfrac{{E}(\kappa^{{-}1})}{{K}(\kappa^{{-}1})} & \\ \quad - \dfrac{4 \alpha_M}{3} \kappa^2 \left((2 \kappa^2 -1 ) \dfrac{{E}(\kappa^{{-}1})}{{K}(\kappa^{{-}1})} - 2 \kappa^2 + 2 \right) & \text{if } 1 < \kappa < \infty \text{ (passing)}. \end{array}\right. \end{align}

Rather than calculating $\alpha _3$ explicitly, it suffices to write its generic integral form. We also include the various oscillating quantities associated with the precession motion. In general, we have

(2.79)\begin{gather} \psi = \bar{\psi} + \tilde{\psi}, \end{gather}

(2.80)\begin{gather}\varphi = \alpha_3 + q(\bar{r}) \tilde{\theta} + \tilde{\varphi}, \end{gather}

(2.81)\begin{gather}\theta = \bar{\epsilon} \alpha_2 + \tilde{\theta}. \end{gather}

Here, $\tilde {\psi }$ represents the excursion from the reference flux surface $\bar {\psi }$ during the poloidal orbit. Meanwhile, $\tilde {\varphi }$ is the difference in toroidal precession between a circular geometry and a more general equilibrium magnetic field. Later, we will use $\tilde {r}$ instead of $\tilde {\psi }$ with the understanding that $\psi (\tilde {r}) = \tilde {\psi }$. Meanwhile, $\tilde {\theta }$ is associated with the oscillatory poloidal motion. We define these quantities as

(2.82)\begin{gather} \tilde{\psi} = \int^{\alpha_2} \frac{\textrm{d} \alpha_2'}{\varOmega_2} \boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla} \psi = \psi_1 \end{gather}

(2.83)\begin{gather}\tilde{\theta} = \int^{\alpha_2} \, \textrm{d} \alpha_2' \left(\frac{1}{\varOmega_2} \frac{\textrm{d} \theta}{\textrm{d} t} - \bar{\epsilon}\right), \end{gather}

(2.84)\begin{gather}\tilde{\varphi} = \frac{\textrm{d} q}{\textrm{d} \psi} \tilde{\theta} \tilde{\psi} + \int^{\alpha_2} \frac{\textrm{d} \alpha_2'}{\varOmega_2} \left( \bar{\epsilon} \frac{\textrm{d} q}{\textrm{d} \psi} \tilde{\psi} \varOmega_2 + \boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla}\varphi - \frac{\textrm{d} q}{\textrm{d} \psi} \tilde{\theta} \boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla} \psi - q \boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla} \theta - \varOmega_d\right). \end{gather}

This guarantees that

(2.85)\begin{equation} \frac{\textrm{d} \alpha_3}{\textrm{d} t} = \varOmega_d + \bar{\epsilon} q \varOmega_2. \end{equation}

Note that in the above, $q$ and $\textrm {d} q/\textrm {d} \psi$ are evaluated at $\bar {r}$ and thus are time independent. Having characterized action-angle coordinates, we can proceed to solving the Vlasov equation using these coordinates.

3. The Vlasov equation

To begin, we write the Vlasov equation in action-angle variables:

(3.1)\begin{equation} \frac{\textrm{d} f}{\textrm{d} t} = \frac{\partial f}{\partial t} + \dot{\boldsymbol{\alpha}}\boldsymbol{\cdot} \frac{\partial f}{\partial \boldsymbol{\alpha}} + \dot{\boldsymbol{J}} \boldsymbol{\cdot} \frac{\partial f}{\partial \boldsymbol{J}} = 0. \end{equation}

We remind ourselves that Hamilton's equations of motion in these coordinates are

(3.2)\begin{gather} \dot{\boldsymbol{J}} ={-} \frac{\partial H}{\partial \boldsymbol{\alpha}} ={-} e \frac{\partial \phi}{\partial \boldsymbol{\alpha}}, \end{gather}

(3.3)\begin{gather}\dot{\boldsymbol{\alpha}} = \frac{\partial H}{\partial \boldsymbol{J}} = e \frac{\partial \phi}{\partial \boldsymbol{J}} + \boldsymbol{\varOmega}. \end{gather}

We later generalize the above equation with a Krook-style operator to add collisions for trapped electrons in Appendix C, but for now we work in the collisionless limit. The next step is to linearize the system by assuming the distribution function is composed of an equilibrium part $f_0 = f_0(\boldsymbol {J})$ and a perturbed part $\delta f = \delta f(\boldsymbol {\alpha }, \boldsymbol {J}, t).$ Dropping any quadratic perturbative terms, we obtain

(3.4)\begin{equation} \frac{\partial \delta f}{\partial t} + \boldsymbol{\varOmega} \boldsymbol{\cdot} \frac{\partial \delta f}{\partial \boldsymbol{\alpha}} - e \frac{\partial \phi}{\partial \boldsymbol{\alpha}} \boldsymbol{\cdot} \frac{\partial f_0}{\partial \boldsymbol{J}} = 0. \end{equation}

As stated earlier, any perturbative functions we consider must be periodic in the angular variables $\boldsymbol {\alpha }$. Therefore, we utilize a discrete Fourier transform in $\delta f$ and $\phi$:

(3.5)\begin{gather} \delta f = \sum_{\boldsymbol{n}} f_{\boldsymbol{n}}(\boldsymbol{J}) \, \textrm{e}^{{\rm i} (\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol {\alpha} - \omega t)}, \end{gather}

(3.6)\begin{gather}\phi = \sum_{\boldsymbol{n}} \phi_{\boldsymbol{n}}(\boldsymbol{J}) \, \textrm{e}^{{\rm i} (\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol {\alpha} - \omega t)}. \end{gather}

To extract the physical quantity, we take the real part of the Fourier series. Here, $\boldsymbol {n}$ corresponds to the mode number of the Fourier term and $\omega$ is the complex frequency of oscillation. We decompose the complex frequency as $\omega = \omega _r + {\rm i} \gamma$, where $\omega _r$ is the real frequency and $\gamma$ is the growth rate. Note that in QuaLiKiz, we only consider unstable modes with $\gamma > 0$ and ignore stable modes; although this does not change the fundamental approach, it does afford us some slight computational simplicity since we do not have to search for solutions in the entire complex plane. As an ansatz, we treat $\omega = \omega _{n_3}$ to be dependent on $n_3$ only, not $n_1$ and $n_2$. To consistently solve the dispersion relation, we will eventually need to sum over $n_1$ and $n_2$. The individual Fourier components can be calculated from the physical quantity via

(3.7)\begin{gather} f_{\boldsymbol{n}} = \int \frac{\textrm{d}^3\alpha}{(2 {\rm \pi})^3} \delta f(t=0) \, \textrm{e}^{-{\rm i}\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol {\alpha}}, \end{gather}

(3.8)\begin{gather}\phi_{\boldsymbol{n}} = \int \frac{\textrm{d}^3\alpha}{(2 {\rm \pi})^3} \phi(t=0) \, \textrm{e}^{-{\rm i}\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol {\alpha}}, \end{gather}

where we integrate each angular variable from $0$ to $2 {\rm \pi}$.

To proceed, we assume the equilibrium distribution is a shifted Maxwellian:

(3.9)\begin{align} f_0(\boldsymbol{J}) &= n_0 \left(\frac{m}{2 {\rm \pi}T}\right)^{3/2} \exp\left(- \frac{m (\boldsymbol{v} - \boldsymbol{U})^2} {2T}\right)\nonumber\\ &= n_0 \left(\frac{m}{2 {\rm \pi}T}\right)^{3/2} \exp\left( \frac{-H_0 + e \varPhi + m \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{U} - \dfrac{m U^2}{2}}{T} \right). \end{align}

Here, $n_0$ is the equilibrium number density, $T$ is the temperature and $\boldsymbol {U}$ is the equilibrium plasma rotation velocity. In general, $n_0$, $T$ and $\boldsymbol {U}$ will vary with position and therefore depend on $\boldsymbol {J}$. By only considering toroidal rotation, we make the approximation

(3.10)\begin{equation} U_\varphi \hat{\boldsymbol{\varphi}} \approx U_\parallel \hat{\boldsymbol{b}}, \end{equation}

which allows us to write

(3.11)\begin{gather} U^2 \approx U_\parallel^2, \end{gather}

(3.12)\begin{gather}\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{U} \approx v_\parallel U_\parallel. \end{gather}

We also take into account gradients of the parallel rotation velocity. Due to the presence of rotation, we also include the radial electric field as well as its gradient. We use the natural frequency parameter for the electric field shear $\gamma _E$ defined as

(3.13)\begin{equation} \gamma_E ={-}\frac{1}{B} \frac{\textrm{d} E_r}{\textrm{d} r}. \end{equation}

This will allow us to Taylor expand the characteristic $E$-cross-$B$ frequency such that

(3.14)\begin{equation} \omega_E \approx \omega_{E0} + \omega_E' x, \end{equation}

where we expand about the radial distance $x = 0$ and $\partial _r \omega _E = \omega _E'$ is related to the radial electric shear.

Additionally, in QuaLiKiz we ignore terms that go as the square of fluctuating quantities. Since we assume a small Mach number as well as a small derivative in the parallel velocity, we thus assume that

(3.15)\begin{equation} m U_\parallel \left|\frac{\partial v_\parallel}{\partial \boldsymbol{J}}\right| \ll \boldsymbol{\omega}_\ast. \end{equation}

This term is responsible for turbulent acceleration and arises from the presence of rotation in the equilibrium distribution function. We expect this term to be negligible for non-impurities in the low-Mach-number limit and will thus neglect it as an approximation (Garbet et al. Reference Garbet, Esteve, Sarazin, Abiteboul, Bourdelle, Dif-Pradalier, Ghendrih, Grandgirard, Latu and Smolyakov2013).

Substituting the above expressions as well as the Fourier series into the linearized Vlasov equation, we isolate each term mode by mode due to completeness and orthogonality of the Fourier series to find $f_{\boldsymbol {n}}$ in terms of $\phi _{\boldsymbol {n}}$. The result is

(3.16)\begin{equation} f_{\boldsymbol{n}} = \frac{e \phi_{\boldsymbol{n}} \boldsymbol{n} \boldsymbol{\cdot} \dfrac{\partial f_0}{\partial \boldsymbol{J}}}{\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\varOmega} - \omega} = \frac{f_0}{T}\frac{e \phi_{\boldsymbol{n}} \boldsymbol{n} \boldsymbol{\cdot} (\boldsymbol{\omega}_\ast{+} \boldsymbol{\omega}_E - \boldsymbol{\varOmega})}{\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\varOmega} - \omega}, \end{equation}

where the diamagnetic frequency $\boldsymbol {\omega }_\ast$ is

(3.17)\begin{equation} \boldsymbol{\omega}_\ast{=} T \left( \frac{1}{n_0} \frac{\partial n_0}{\partial \boldsymbol{J}} + \left(\xi - \frac{3}{2} - \frac{U_\parallel}{v_T^2} (2 v_\parallel{-} U_\parallel) \right) \frac{1}{T} \frac{\partial T}{\partial \boldsymbol{J}}\right) + \frac{2 (v_\parallel{-} U_\parallel )}{v_T^2} \frac{\partial U_\parallel}{\partial \boldsymbol{J}}, \end{equation}

where the thermal velocity is $v_T = \sqrt {2 T /m }$ and the frequency associated with the $E$-cross-$B$ drift is

(3.18)\begin{equation} \boldsymbol{\omega}_E = \frac{e}{T} \frac{\textrm{d} \varPhi}{\textrm{d} \boldsymbol{J}}. \end{equation}

We then rewrite the equation to be

(3.19)\begin{equation} f_{\boldsymbol{n}} ={-}\frac{e \phi_{\boldsymbol{n}}}{T} f_0 \left(1 - \frac{\omega - \boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\omega}_\ast{-} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_E}{\omega - \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\varOmega}} \right), \end{equation}

where it is now clear that there is an adiabatic part and a frequency dependent part of the equation.

The next step to solving the dispersion relation is to use Poisson's equation,

(3.20)\begin{equation} \nabla^2 \phi = \sum_s - \frac{e_s n_s}{\epsilon_0}, \end{equation}

where the $s$ subscript labels the particle species and $\epsilon _0$ is the vacuum permittivity. In the earlier parts of the derivation, we had suppressed the subscript for various quantities (e.g. $m, T, n_0, f_0, \delta f, \ldots$); we include the subscript for the time being. The total number density $n_s$ is

(3.21)\begin{equation} n_s = \int \, \textrm{d}^3 v f_s. \end{equation}

The perturbed electrostatic potential is calculated using the perturbed number density

(3.22)\begin{equation} \delta n_s = \int \, \textrm{d}^3 v \delta f_s. \end{equation}

We obtain the perturbed charge density by then multiplying the perturbed number density by the charge of the species. To enforce quasineutrality, we take the sum of the total charge density to be null and require that

(3.23)\begin{equation} \lambda_D \ll \left|\frac{\phi}{\boldsymbol{\nabla} \phi}\right|, \end{equation}

where $\lambda _D$ is the Debye length. Because we are interested in length scales much longer than the Debye length, the Laplacian term in Poisson's equation is negligible. The quasineutrality constraint then becomes

(3.24)\begin{equation} \sum_s \int \, \textrm{d}^3 v e_s \delta f_s = 0. \end{equation}

Since $\phi = \phi (\boldsymbol {r})$ is independent of velocity, if we multiply both sides of the above equation by $\phi ^\ast$, the complex conjugate of $\phi$, we can simply move it inside the integral. We then integrate over space, which results in

(3.25)\begin{equation} \sum_s \int \, \textrm{d}^3r \, \textrm{d}^3 v e_s \phi^\ast \delta f_s = 0. \end{equation}

By multiplying by the electrostatic potential and integrating, we have recast the differential equation via a weak formulation using the variational method (Samain Reference Samain1970; Garbet et al. Reference Garbet, Laurent, Mourgues, Roubin and Samain1990; Garbet Reference Garbet2001; Nguyen, Garbet & Smolyakov Reference Nguyen, Garbet and Smolyakov2008). Instead of solving for the exact functions $\phi$ or $\delta f$ that satisfy Poisson's equation, we can simply approximate $\phi$ and $\delta f$ with a suitable function and focus on the dispersion relation itself. Typically, when the Laplacian is kept, the differential equation is put into the weak formulation by integrating the Laplacian term by parts; this technique is well established in other fields such as finite element analysis (Johnson Reference Johnson1991).

We next substitute in the Fourier expansions and the expression relating $\phi _{\boldsymbol {n}}$ and $f_{\boldsymbol {n}}$. The result is

(3.26)\begin{align} \sum_s \sum_{\boldsymbol{n}, \boldsymbol{n}^\prime} \int \, \textrm{d}^3r \, \textrm{d}^3v \frac{e_s^2 \phi_{\boldsymbol{n}} \phi_{\boldsymbol{n}^\prime}^\ast}{T_s} f_{0,s} \left(1 - \frac{\omega - \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_\ast{-} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_E}{\omega - \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\varOmega}}\right) \, \textrm{e}^{{\rm i} (\boldsymbol{n} - \boldsymbol{n}^\prime) \boldsymbol{\cdot} \boldsymbol{\alpha}} \, \textrm{e}^{-{\rm i} (\omega - (\omega^\prime)^\ast) t} = 0. \end{align}

To simplify this integral, we first perform the change of variables $(\boldsymbol {r}, \boldsymbol {v}) \to (\boldsymbol {r}, \boldsymbol {p})$; the Jacobian of this transformation is simply $m_s^{-1}$. We then perform the change of variables $(\boldsymbol {r}, \boldsymbol {p}) \to (\boldsymbol {\alpha }, \boldsymbol {J})$; the Jacobian of this particular transformation is $1$ because this is guaranteed to be a canonical transformation. We therefore obtain

(3.27)\begin{align} \sum_s \sum_{\boldsymbol{n}, \boldsymbol{n}^\prime} \int \, \textrm{d}^3\alpha \, \textrm{d}^3 J \frac{e_s^2 \phi_{\boldsymbol{n}} \phi_{\boldsymbol{n}^\prime}^\ast}{m_s T_s} f_{0,s} \left(1 - \frac{\omega - \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_\ast{-} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_E}{\omega - \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\varOmega}}\right) \, \textrm{e}^{{\rm i} (\boldsymbol{n} - \boldsymbol{n}^\prime) \boldsymbol{\cdot} \boldsymbol{\alpha}} \, \textrm{e}^{-{\rm i} (\omega - (\omega^\prime)^\ast) t} = 0. \end{align}

We note that the exponential terms are $\boldsymbol {\alpha }$ dependent. We then use orthogonality of the Fourier series to find that

(3.28)\begin{equation} \sum_s \sum_{n_1, n_2} \int \, \textrm{d}^3 J \frac{e_s^2 |\phi_{\boldsymbol{n}}|^2 }{m_s T_s} f_{0,s} \left(1 - \frac{\omega - \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_\ast{-} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_E}{\omega - \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\varOmega}}\right) = 0. \end{equation}

Because we developed the Fourier series such that $\omega$ only depends on the mode number $n_3$, we can solve for each value of $n_3$ individually while summing over $n_1$ and $n_2$. While the summation arising from this convention seems to make the problem more difficult at first glance, we shall see later it allows for a variety of simplifications. Moreover, the integrand is now completely independent of $\boldsymbol {\alpha }$. As such, we integrate over the action angles again and transform back to conventional variables, which leads to

(3.29)\begin{equation} \sum_s \sum_{n_1, n_2} \int \, \textrm{d}^3r \, \textrm{d}^3v \frac{e_s^2 |\phi_{\boldsymbol{n}}|^2 }{T_s} f_{0,s} \left(1 - \frac{\omega - \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_\ast{-} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_E}{\omega - \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\varOmega}}\right) = 0. \end{equation}

Even though the integrand is a function of only $\boldsymbol {J}$, the parameters in the integrand are more naturally expressed in terms of other coordinates such as the minor radius and the pitch angle parameter. Thus, further coordinate transformations to simplify this expression are inevitable. As such, they are most easily carried out when starting from the typical configuration space variables $(\boldsymbol {r}, \boldsymbol {v}).$

For ease of notation, we split up the dispersion relation as follows:

(3.30)\begin{equation} \sum_s \mathcal{L}_{0, s} - \mathcal{L}_{\text{passing}, s} - \mathcal{L}_{\text{trapped}, s} = 0. \end{equation}

Here, $\mathcal {L}_0$ is the portion of the integral that is simply multiplied by $1$, which we call the adiabatic part. Then, $\mathcal {L}_{\text {passing}}$ is the portion of the integral that is frequency dependent and integrated over the part of velocity space that encompasses passing particles, while $\mathcal {L}_{\text {trapped}}$ consists of the trapped particles instead.

To proceed with solving the dispersion relation, we must first calculate $|\phi _{\boldsymbol {n}}|^2$. This requires a three-dimensional (3-D) integral over $\textrm{d}^3 \alpha$ as shown in (3.8). Once that is done, we then proceed to calculate the integral in the dispersion relation itself for the adiabatic part, trapped part and passing part separately. Although our expression appears to be a six-dimensional integral, we can utilize a number of symmetries, transformations and approximations to simplify the form down to at most 2-D integrals. Although integrals of higher dimension can be in principle calculated numerically, the curse of dimensionality renders such integrals computationally expensive. Thus, a reduction to two dimensions affords us a great deal of speed at the cost of some amount of accuracy.

4. Ballooning representation

Before integrating $|\phi _{\boldsymbol {n}}|^2$ with respect to the action angles, we review key results regarding the ballooning representation. Because $\phi (r,\theta ,\varphi )$ must be periodic in $\theta$ and $\varphi$, we may expand $\phi$ as a Fourier series,

(4.1)\begin{equation} \phi(r, \theta, \varphi) = \sum_{m, n} \phi_{m,n}(r - r_0) \, \textrm{e}^{{\rm i} (m \theta + n \varphi)}. \end{equation}

Here, $r_0$ is the location of the resonant flux surface for each given $m$ and $n$; in other words, $q(r_0) = -m/n$. We take these modes to be localized around the resonant flux surface. These modes are often radially localized such that the distance between any two adjacent resonant rational flux surfaces is much longer than the characteristic length scale of the plasma equilibrium. If that condition holds, then all modes $\phi _{m,n}$ have nearly identical radial envelopes where each radial profile is centred on their corresponding reference flux surface (Connor, Hastie & Taylor Reference Connor, Hastie and Taylor1979). These flux surfaces are all rational flux surfaces since $m$ and $n$ are integers. Meanwhile, the general ballooning representation of $\phi$ is

(4.2)\begin{equation} \phi(r, \theta, \varphi) = \sum_p \sum_n \hat{\phi}_{n} (\theta + 2 p {\rm \pi}, \theta_0 ) \, \textrm{e}^{{\rm i} n (\varphi - q(r) (\theta - \theta_0 + 2 p {\rm \pi}))}, \end{equation}

where $\theta _0$ is the ballooning angle and $p$ denotes the various harmonics. Here, we have approximated the potential by separating it into a quickly varying eikonal and a slowly varying envelope.This representation ultimately comes from the fact that the instabilities in question are strongly anisotropic and flute-like where $k_\parallel \gg k_{\perp }$. In the absence of toroidal rotation, the ballooning angle is typically taken to be zero since the most unstable modes are localized around $\theta _b = 0$. In the presence of finite toroidal rotation and an equilibrium electrostatic potential, the ballooning angle is shifted away from zero (Bourdelle et al. Reference Bourdelle, Citrin, Baiocchi, Casati, Cottier, Garbet and Imbeaux2016). However, it is known that the shift due to rotation is typically of the order of $10^{-1}$ in the tokamak core (Candy, Waltz & Rosenbluth Reference Candy, Waltz and Rosenbluth2004). Although up-down asymmetries and global effects can lead to an effective non-zero ballooning angle not necessarily close to zero, QuaLiKiz does not take into account global effects nor up-down asymmetries. For the rest of the derivation, we therefore take the ballooning angle to be zero as an approximation. This is equivalent to assuming that the envelope is radially independent. Moreover, if the profile is heavily localized around $\theta = 0$, we can use the strong ballooning approximation and ignore all harmonics except for $p = 0$, which leads to

(4.3)\begin{equation} \phi(r, \theta, \varphi) = \sum_n \hat{\phi}_n(\theta) \, \textrm{e}^{{\rm i} n (\varphi - q(r) \theta)}. \end{equation}

It is important to note that the decomposition in terms of $\phi _{m,n}$ describes how the same radial profile is localized about adjacent flux surfaces. Meanwhile, the decomposition in terms of $\hat {\phi }_n$ describes how the linear eigenmode balloons along the field line. This can be seen more explicitly if one considers that

(4.4)\begin{gather} \boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi = 0, \end{gather}

(4.5)\begin{gather}\boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} (\varphi - q \theta) = 0, \end{gather}

(4.6)\begin{gather}\boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} \theta \neq 0. \end{gather}

Thus, with the above set of variables, $\theta$ indicates the location on any given field line. Because the magnetic curvature is unfavourable on the low-field side of the tokamak when one considers the interplay between the curvature vector and the pressure gradient for normal tokamak profiles, we expect fluctuations to peak about $\theta = 0$. We can demonstrate a direct link between $\phi _{m,n}$ and $\hat {\phi }_n$ by calculating the Fourier components of $\phi$, which leads to

(4.7)\begin{equation} \int_{-{\rm \pi}}^{\rm \pi} \frac{\textrm{d} \varphi}{2{\rm \pi}} \int_{-{\rm \pi}}^{\rm \pi} \frac{\textrm{d} \theta}{2{\rm \pi}} \phi(r, \theta, \varphi) \, \textrm{e}^{- {\rm i} n \varphi - i m \theta} = \phi_{m,n}(r-r_0) = \int_{-{\rm \pi}}^{\rm \pi} \frac{\textrm{d} \theta}{2 {\rm \pi}} \hat{\phi}_n (\theta) \, \textrm{e}^{- {\rm i} \theta (n q(r) + m)}. \end{equation}

We then make two approximations. First, we Taylor expand the term in the eikonal around the reference flux surface,

(4.8)\begin{equation} n q(r) + m = n q_0 + \frac{r - r_0}{d} + m = \frac{r - r_0}{d} = \frac{x}{d} + \frac{\tilde{r}}{d}, \end{equation}

where $q_0 = q(r_0)$, the radial difference between different rational flux surfaces is defined as

(4.9)\begin{equation} \frac{1}{nd} = \left.\frac{\textrm{d} q}{\textrm{d} r}\right|_{r = r_0}, \end{equation}

and $x$ is defined as

(4.10)\begin{equation} x = \bar{r} - r_0, \end{equation}

where we ignore second derivatives of the safety factor. After doing so, we find that

(4.11)\begin{equation} \phi_{m,n}(r-r_0) \approx \int_{-{\rm \pi}}^{\rm \pi} \frac{\textrm{d} \theta}{2 {\rm \pi}} \hat{\phi}_n (\theta) \, \textrm{e}^{- {\rm i} ({\theta (r - r_0)}/{d})}. \end{equation}

Second, we invoke the strong ballooning approximation by treating $\hat {\phi }_n$ as heavily localized around $\theta = 0$; this allows us to integrate from $-\infty$ to $\infty$ instead of from $-{\rm \pi}$ to ${\rm \pi}$. The result is

(4.12)\begin{equation} \phi_{m,n}(r-r_0) \approx \int_{-\infty}^{\infty} \frac{\textrm{d} \theta}{2 {\rm \pi}} \hat{\phi}_n (\theta) \, \textrm{e}^{- {\rm i} ({\theta (r - r_0 )}/{d})}. \end{equation}

Comparing it with our previous definition of the Fourier transform, we find that $\hat {\phi }_n(\theta )$ is simply the Fourier transform of $\phi _{m,n}(r)$, with $k_r = \theta /d$. The transformation is given by

(4.13)\begin{equation} \hat{\phi}_n(\theta) = \int^{\infty}_{-\infty} \frac{\textrm{d} r}{|d|} \phi_{m,n}(r-r_0) \, \textrm{e}^{{\rm i} ({\theta (r-r_0 )}/{d})}. \end{equation}

We are now in a position to integrate over the action angles to fully calculate $\phi _{\boldsymbol {n}}$. The procedure to integrate over $\alpha _1$ has already been discussed in § 2, with the relevant result being (2.29). We therefore only need to discuss in detail the integrations over $\alpha _2$ and $\alpha _3$ while treating all variables within the guiding centre framework. Trapped particle motion and passing particle motion differ such that the two cases must be handled separately.

4.1. Trapped

For deeply trapped particles, the equations for the action variables simplify to

(4.14)\begin{gather} r = \bar{r} + \delta_b \cos(\alpha_2), \end{gather}

(4.15)\begin{gather}\theta = \theta_b \sin(\alpha_2), \end{gather}

(4.16)\begin{gather}\varphi = \alpha_3 + q(\bar{r}) \theta_b \sin(\alpha_2) + \tilde{\varphi}. \end{gather}

Here, we define the banana width $\delta _b$ as

(4.17)\begin{equation} \delta_b = \frac{q \rho}{\sqrt{\epsilon}}. \end{equation}

While more exact expressions for the bounce motion can be given using Jacobi elliptic functions, we use the above equations for all trapped particles as an approximation. We first integrate over $\alpha _2$, once again utilizing the Fourier transform,

(4.18)\begin{equation} \int^{\rm \pi}_{-{\rm \pi}} \frac{\textrm{d} \alpha_2}{2 {\rm \pi}} \phi(\boldsymbol{r}) \, \textrm{e}^{- {\rm i} n_2 \alpha_2} = \int^{\rm \pi}_{-{\rm \pi}} \frac{\textrm{d} \alpha_2}{2 {\rm \pi}} \int \frac{\textrm{d}^3 k}{(2 {\rm \pi})^3}\tilde{\phi}(\boldsymbol{k}) \, \textrm{e}^{-{\rm i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{r} - {\rm i} n_2 \alpha_2}. \end{equation}

We then proceed in a fashion similar to the gyro-average derivation in § 2 by noting that $\boldsymbol {k} \boldsymbol {\cdot } {\boldsymbol {r}} = \boldsymbol {k} \boldsymbol {\cdot } \bar {\boldsymbol {r}} + k_r \delta _b \cos (\alpha _2)$. The result is

(4.19)\begin{equation} \int^{\rm \pi}_{-{\rm \pi}} \frac{\textrm{d} \alpha_2}{2 {\rm \pi}} \phi(\boldsymbol{r}) \, \textrm{e}^{- {\rm i} n_2 \alpha_2} = \int \frac{\textrm{d}^3 k}{(2 {\rm \pi})^3} ({-}{\rm i})^{n_2} J_{n_2}(k_r \delta_b) \tilde{\phi}(\boldsymbol{k}) \, \textrm{e}^{-{\rm i} \boldsymbol{k} \boldsymbol{\cdot} \bar{\boldsymbol{r}}}. \end{equation}

In essence, we obtain a bounce average over the banana width (Depret et al. Reference Depret, Garbet, Bertrand and Ghizzo2000). We note that this in particular is a rather crude approximation. The particularities of the bounce motion, such as the bounce angle and the radial excursion, technically depend on the pitch angle of the particle; we are in essence smearing this out by taking a representative trapped particle such that the banana width is constant. The averaging procedure is also approximate as we only take into account the radial deviation. As seen in the work of Biglari & Chen (Reference Biglari and Chen1986), we would normally obtain a $\theta$ dependence in the argument of the Bessel function; the $k_r$ term manifests as the Fourier link established earlier. As a result, trapped particles have two Bessel operators acting on the potential corresponding to the gyromotion and the banana orbit respectively.

We now proceed to the integration over $\alpha _3$, for now leaving the Bessel functions aside and evaluating the position at $\boldsymbol {r} = \bar {\boldsymbol {r}}$. In doing so, we must be aware that for trapped particles $\bar {\theta } = 0$; that is, the variation of $\theta$ only comes from the bounce orbit which we averaged over. We also ignore $\tilde {\varphi }$ for the same reason. Moreover, because we assume the modes to have an identical radial envelope, we are free to keep only one of the poloidal harmonics as there is an effective decay in the other poloidal harmonics. Making the strong assumption that the actual radial envelope can be approximated in this way, we pick $m_0 = - n_3 q_0$, as this forces any $\theta$ dependence we approximately neglected in the eikonal to vanish. Thus, we obtain

(4.20)\begin{equation} \int_{-{\rm \pi}}^{\rm \pi} \frac{\textrm{d} \alpha_3}{2{\rm \pi}} \phi(\bar{\boldsymbol{r}}) \, \textrm{e}^{- {\rm i} n_3 \alpha_3} = \sum_{n} \int_{-{\rm \pi}}^{\rm \pi} \frac{\textrm{d} \alpha_3}{2{\rm \pi}} \phi_{m_0,n}(\bar{r}- r_0) \, \textrm{e}^{{\rm i} (n \alpha_3 - n_3 \alpha_3)} = \phi_{m_0,n_{3}}(\bar{r}- r_0). \end{equation}

To compensate for choosing only one poloidal harmonic, when we later perform the integrals required for the dispersion relation, we extend the radial limits of integration to $-\infty < r < \infty$. Essentially, because the mode structure corresponding to different flux surfaces does not change, integrating over a finite domain while retaining the sum over all harmonics is equivalent to taking one harmonic and integrating over an infinite domain (Garbet et al. Reference Garbet, Laurent, Mourgues, Roubin and Samain1990). Aside from the Bessel functions, nothing in the trapped part of the dispersion relation is dependent on $\theta = k_r d$. Therefore, we are free to take the amplitude squared of the averaged potential to obtain

(4.21)\begin{equation} |\phi_{\boldsymbol{n}}|^2 = |J_{n_1}(k_\perp \rho)J_{n_2}(\delta_b k_r) \boldsymbol{\cdot} \phi_{m_0, n_3}|^2(\bar{r} - r_0) = |J_{n_1}(k_\perp \rho)J_{n_2}(\delta_b k_r) \boldsymbol{\cdot} \phi_{m_0, n_3}|^2(x), \end{equation}

where we evaluate the function at $x = \bar {r} - r_0$.

4.2. Passing

We now calculate $\phi _{\boldsymbol {n}}$ for passing particles. Instead of utilizing the poloidal harmonics, it is more useful to use the ballooning representation directly. Substituting in the expression for $\alpha _3$ and then integrating over $\alpha _3$ leads to

(4.22)\begin{equation} \int_{-{\rm \pi}}^{\rm \pi}\frac{\textrm{d} \alpha_3}{2 {\rm \pi}} \phi(\boldsymbol{r}) \, \textrm{e}^{-{\rm i} n_3 \alpha_3} = \sum_n \int_{-{\rm \pi}}^{\rm \pi}\frac{\textrm{d} \alpha_3}{2 {\rm \pi}} \hat{\phi}_n(\theta(\alpha_2))\, \textrm{e}^{{\rm i} n(\alpha_3 + \tilde{\varphi} - q(r) \alpha_2 +(q(\bar{r}) - q(r)) \tilde{\theta}) - {\rm i} n_3 \alpha_3}. \end{equation}

Here, $\theta$ is taken to be a function of $\alpha _2$. It is crucial that we recognize not all the safety factors in the eikonal are evaluated at the same point. We have both $q(r) = q(\bar {r} + \tilde {r})$ and $q(\bar {r})$. The term $q(r) - q(\bar {r})$ can be Taylor expanded about $r_0$:

(4.23)\begin{equation} q(r) - q(\bar{r}) \approx q_0 + \frac{r - r_0}{nd} - q_0 - \frac{\bar{r}-r_0}{nd} = \frac{r - \bar{r}}{nd} = \frac{\tilde{r}}{nd} . \end{equation}

Carrying out the integral then gives us

(4.24)\begin{equation} \int_{-{\rm \pi}}^{\rm \pi}\frac{\textrm{d} \alpha_3}{2 {\rm \pi}} \phi(\boldsymbol{r}) \, \textrm{e}^{-{\rm i} n_3 \alpha_3} = \hat{\phi}_{n_3}(\theta(\alpha_2)) \, \textrm{e}^{{\rm i} n_3 (\tilde{\varphi} - q(r) \alpha_2 - ({\tilde{r}}/{n_3 d}) \tilde{\theta})}. \end{equation}

We now multiply by $\textrm {e}^{- \textrm {i} n_2 \alpha _2}$ and integrate with respect to $\alpha _2$. The eikonal can be simplified if we only keep $n_2 = m_0 = - n q_0$,

(4.25)\begin{align} {\rm i} \left(n_3 \tilde{\varphi} - n_3 q(r) \alpha_2 - \frac{\tilde{r}}{d} \tilde{\theta} - m_0 \alpha_2\right) &= {\rm i} \left(n_3 \tilde{\varphi} - \frac{x}{d} \alpha_2 - \frac{\tilde{r}}{d}(\alpha_2 + \tilde{\theta}))\right)\nonumber\\ &= {\rm i} \left(n_3 \tilde{\varphi} - \frac{x}{d} \alpha_2 - \frac{\tilde{r}}{d} \theta(\alpha_2)\right), \end{align}

where we have used

(4.26)\begin{align} -n_3 q(r) \alpha_2 - m_0 \alpha_2 \approx{-}n_3 q_0 \alpha_2 - \frac{r - r_0}{d} \alpha_2 - m_0 \alpha_2 ={-}\frac{r - r_0}{d} \alpha_2 ={-}\frac{x}{d} \alpha_2 - \frac{\tilde{r}}{d} \alpha_2 \end{align}

and the expression $\alpha _2 = \theta - \tilde {\theta }$. We then obtain

(4.27)\begin{align} \int_{-{\rm \pi}}^{\rm \pi} \frac{\textrm{d} \alpha_2}{2 {\rm \pi}} \hat{\phi}_{n_3} (\theta(\alpha_2))\, \textrm{e}^{{\rm i} (n_3 \tilde{\varphi} - ({x}/{d}) \alpha_2 - ({\tilde{r}}/{d}) \theta(\alpha_2))} \approx \int_{-\infty}^{\infty} \frac{\textrm{d} \alpha_2}{2 {\rm \pi}} \hat{\phi}_{n_3} (\theta(\alpha_2))\, \textrm{e}^{{\rm i} (n_3 \tilde{\varphi} - ({\tilde{r}}/{d}) \theta(\alpha_2))} \, \textrm{e}^{-{\rm i}({x}/{d})\alpha_2}, \end{align}

where we have invoked the strong ballooning approximation. We can see that this is simply an inverse Fourier transform going from $\alpha _2$ to $x$. Thus, we write that

(4.28)\begin{equation} \phi_{\boldsymbol{n}} = \phi_{n_1, m_0, n_3} = (J_{n_1}(k_\perp \rho) \boldsymbol{\cdot} \mathcal{F}^{{-}1}(\hat{\phi}_{n_3} (\theta(\alpha_2))\, \textrm{e}^{{\rm i} (n_3 \tilde{\varphi} - ({\tilde{r}}/{d}) \theta(\alpha_2))})) (x), \end{equation}

where $\mathcal {F}^{-1}$ inverts the Fourier transform as described above with respect to $\alpha _2$. While the $\tilde {\varphi }$ dependence can be approximately ignored in a circular geometry, the $\tilde {r} \theta (\alpha _2)$ dependence in the eikonal must be kept.

We shall see that the mode numbers $n_1$ and $n_2$ do not appear explicitly in the final expression. For convenience, we thus write $n_3 = n$ and identify it as the toroidal mode number.

4.3. Gaussian eigenfunction

We now introduce the functional form of the potential. We use the ansatz that the poloidal harmonic structure is a shifted Gaussian:

(4.29)\begin{equation} \phi_{m_0, n}(x) \sim \phi_0 \, \textrm{e}^{-{(x - x_0)^2}/{2 w^2}}. \end{equation}

This Gaussian corresponds to the lowest-order eigenfunction from the corresponding ballooning equation; the higher-order terms utilize the Hermite polynomials and are neglected here. This Gaussian has a complex width $w$ and shift $x_0$. In the limit of no rotation, $x_0 = 0$ and the Gaussian is centred about $x = 0$. An asymmetry in the eigenfunction about the origin is required to produce quasilinear angular momentum flux. We also note that we restrict ourselves to up-down symmetric geometries, thereby neglecting intrinsic momentum transport (Ball et al. Reference Ball, Parra, Barnes, Dorland, Hammett, Rodrigues and Loureiro2014). Moreover, even in the absence of rotation, a Gaussian eigenfunction cannot describe modes with odd parity in the electrostatic eigenfunction. Although ETG, ITG and TEM instabilities often have even parity in the electrostatic potential in the absence of rotation, these modes can have odd parity eigenfunctions at strongly driving gradients such as in the pedestal (Pueschel et al. Reference Pueschel, Hatch, Ernst, Guttenfelder, Terry, Citrin and Connor2019). In addition, electromagnetic modes such as the microtearing mode do not have even parity electrostatic potentials (Smith et al. Reference Smith, Guttenfelder, LeBlanc and Mikkelsen2011), not to mention that the parallel vector potential would need to be described. Thus, full incorporation of electromagnetic modes and odd parity electrostatic modes would require a generalization of the Gaussian eigenfunction assumption and thus the use of the aforementioned higher-order terms using the Hermite polynomials. Allowing more general geometry, relaxing the eigenfunction ansatz and including electromagnetic modes would allow for additional mechanisms related to the generation of angular momentum flux. Although the amplitude $\phi _0$ cannot be obtained from quasilinear theory, it factors out of the dispersion relation and does not affect the linear mode frequency calculation. Setting the amplitude will be necessary to calculate the quasilinear fluxes and requires the use of a saturation rule, which is detailed in § 9.

To obtain expressions for $w$ and $x_0$, we move into the high-frequency fluid limit. The original derivation can be found in the paper by Cottier et al. (Reference Cottier, Bourdelle, Camenen, Gürcan, Casson, Garbet, Hennequin and Tala2014) and an extensive, revised derivation can be found in that of Citrin (Reference Citrin2017); here, we shall only discuss the basic principle. We consider the dispersion relation

(4.30)\begin{equation} D(\omega) = \sum_s \int \, \textrm{d}^3 v \frac{f_0 e_s^2}{T_s} \left(1 - J_{0,s}^2 \frac{\omega -n \omega_E - n \omega_{{\ast}, s}}{\omega - n \omega_E - k_\parallel v_\parallel{-} n \omega_{d,s}}\right) \phi_{m_0, n}(x). \end{equation}

This is the local dispersion relation obtained if we consider the strong form of Poisson's equation rather than the weak form; we do not multiply by $\phi ^\ast$ and integrate over space. The Bessel function is such that $J_{0,s} = J_0(k_\perp \rho _s) J_0(k_r \delta _{b,s})$ for trapped particles and $J_{0,s} = J_0(k_\perp \rho _s)$ for passing particles. Meanwhile, we define the parallel wavenumber as $k_\parallel = (k_\theta s x)/(q R_0)$. The local drift frequency $\omega _{d,s}$ is

(4.31)\begin{equation} \omega_{d,s} = \left\{ \begin{array}{@{}ll} \displaystyle\omega_{d0, s} \xi \left({-}1 + \dfrac{2 {E}(\kappa)}{{K}(\kappa)} + 4 s \left(\kappa^2 -1 + \dfrac{{E}(\kappa)}{{K}(\kappa)} \right) \right. & \\ \quad \displaystyle-\left. \dfrac{4 \alpha_M}{3} \left(1 - \kappa^2 - (1 - 2 \kappa^2) \dfrac{{E}(\kappa)}{{K}(\kappa)}\right)\right) & \text{if trapped},\\ \omega_{d0, s} \xi (2 - \lambda b ) ( \cos(\theta) + (s \theta - \alpha_M \sin(\theta)) \sin(\theta) ) & \text{if passing}. \end{array}\right. \end{equation}

The passing form of the drift frequency is due to the radial structure of the eigenfunction as covered in § 7. We also note the Fourier link in the passing drift frequency that $\theta ^2 \to k_r^2 d^2$. To proceed, we take $\bar {\omega } = \omega - n \omega _E$ to be larger than $k_\parallel v_\parallel$ and $\omega _{d0,s}$, and for trapped particles we approximate $k_\parallel v_\parallel \approx 0$. We also take $\delta _{b,e} \ll \delta _{b,i}$ and $\rho _e \ll \rho _i$, where the ‘e’ subscript is for electrons and the ‘i’ subscript is for ions, to obtain

(4.32)\begin{align} D(\omega) &= \left[\frac{n_e}{T_e} \left(1 - \left\langle \left(1 - \frac{n \omega_{{\ast},e}}{\bar{\omega}}\right) \left(1 + \frac{n \omega_{d,e}}{\bar{\omega}} + \frac{n^2 \omega_{d,e}^2}{\bar{\omega}^2} \right) \right\rangle _t \right)\right. \nonumber\\ &\quad + \sum_i \frac{n_i Z_i}{T_i} \left(1 - \left\langle\!\left(1 - \frac{n \omega_{{\ast},i}}{\bar{\omega}} \right) \left(1 + \frac{n \omega_{d,i}}{\bar{\omega}} + \frac{n^2 \omega_{d,i}^2}{\bar{\omega}^2} \right) \left(1 - \frac{k_r^2 \delta_{b,i}^2}{2} \right) \left(1 - \frac{k_\theta^2 \rho_i^2}{2}\right)\!\right\rangle _t \right.\nonumber\\ & \quad - \left\langle \left(1 - \frac{n \omega_{{\ast},i}}{\bar{\omega}} \right) \left(1 + \frac{n \omega_{d,i}}{\bar{\omega}} + \frac{k_\parallel v_\parallel}{\bar{\omega}} + \left(\frac{n \omega_{d,i}}{\bar{\omega}} + \frac{k_\parallel v_\parallel}{\bar{\omega}} \right)^2 \right) \right. \nonumber\\ &\quad \left.\left.\left. \times\left(1 - \frac{k_\theta^2 \rho_i^2}{2} - \frac{k_r^2 \rho_i^2}{2} \right) \right\rangle _p \right) \right] \phi_{m_0, n}(x), \end{align}

where $Z_i$ is the proton number of the ion species. We define the averages over velocity space as

(4.33)\begin{gather} \langle g(\boldsymbol{v})\rangle _t = \int_{\text{trapped}} \, \textrm{d}^3v f_0 g(\boldsymbol{v}) , \end{gather}

(4.34)\begin{gather}\langle g(\boldsymbol{v})\rangle _p = \int_{\text{passing}} \, \textrm{d}^3v f_0 g(\boldsymbol{v}). \end{gather}

To do the $\theta$-dependent terms, we note that this is a differential equation. We approximate the differential operators on $\phi$ in the limit of small mode shift $x_0$, which leads to

(4.35)\begin{equation} \frac{\partial^2 \phi}{\partial x^2} = \left(\frac{x^2}{w^4} - \frac{2 x_0 x}{w^4} - \frac{1}{w^2} \right) \phi. \end{equation}

Next, we carry out the integrals both analytically and numerically as appropriate, and multiply the dispersion relation by $\omega ^3$ to obtain a modified dispersion relation,

(4.36)\begin{equation} \omega^3 D(\omega) = D_0 (\omega) + D_1 (\omega) x + D_2 (\omega) x^2 = 0. \end{equation}

Here, we separate terms proportional to $x^0, x^1$ and $x^2$. With three equations we can solve for the three unknowns $(\omega _0, w, x_0)$. We then find the solution $\omega _0$ such that

(4.37)\begin{equation} D_0 (\omega_0 ) = 0. \end{equation}

Having found this zeroth-order solution, we then find $x$ and $w$ such that

(4.38)\begin{gather} D_1 (\omega_0 ) = 0, \end{gather}

(4.39)\begin{gather}D_2 (\omega_0 ) = 0. \end{gather}

We do not cite the full solution here and direct the reader to the work of Citrin (Reference Citrin2017) for a complete derivation. Now that we have characterized $\phi _{\boldsymbol {n}}$ by calculating $x_0$ and $w$, we move to the dispersion relation itself, beginning with the adiabatic term.

5. Adiabatic functional

We first examine the adiabatic part of the functional, as it is the simplest to treat. It takes the form

(5.1)\begin{equation} \mathcal{L}_0 = \sum_{n_1, n_2} \int \, \textrm{d}^3r \, \textrm{d}^3v \frac{e^2 |\phi_{\boldsymbol{n}}|^2}{ T} f_0. \end{equation}

Here, we have suppressed the subscript $s$ as we will be working with each species independently. We first define a new function $\phi _{n} = \phi _{n}(\alpha _1, \alpha _2, \boldsymbol {J})$ such that

(5.2)\begin{equation} \phi_{n} = \int_{-{\rm \pi}}^{\rm \pi} \frac{\textrm{d} \alpha_3}{2{\rm \pi}} \phi(\boldsymbol{r}) \, \textrm{e}^{-{\rm i} n \alpha_3} = \sum_{n_1, n_2} \phi_{\boldsymbol{n}} \, \textrm{e}^{{\rm i} (n_1 \alpha_1 +n_2 \alpha_2 )}. \end{equation}

We then note due to the orthogonality of the Fourier series that

(5.3)\begin{equation} \int_{-{\rm \pi}}^{\rm \pi} \frac{\textrm{d} \alpha_1}{2 {\rm \pi}} \int_{-{\rm \pi}}^{\rm \pi} \frac{\textrm{d} \alpha_2}{2{\rm \pi}} |\phi_{n}|^2 = \sum_{n_1, n_2} |\phi_{\boldsymbol{n}}|^2. \end{equation}

Thus, it is more convenient to switch back to action-angle coordinates for an intermediate calculation:

(5.4)\begin{equation} \mathcal{L}_0 = \sum_{n_1, n_2} \int \, \textrm{d}^3\alpha \, \textrm{d}^3J \frac{e^2 |\phi_{\boldsymbol{n}}|^2}{m T} f_0 = \sum_{n_1, n_2} \int 4 {\rm \pi}^2 \textrm{d} \alpha_3 \, \textrm{d}^3J \frac{e^2 |\phi_{\boldsymbol{n}}|^2}{m T} f_0 = \int 4 {\rm \pi}^2 \textrm{d} \alpha_3 \, \textrm{d}^3J \frac{e^2 |\phi_n|^2}{m T} f_0. \end{equation}

This then simplifies to

(5.5)\begin{equation} \mathcal{L}_0 = \int 4 {\rm \pi}^2 \textrm{d} \alpha_3 \, \textrm{d}^3J \frac{e^2}{m T} f_0 \int_{-{\rm \pi}}^{\rm \pi} \frac{\textrm{d} \alpha_1}{2 {\rm \pi}} \int_{-{\rm \pi}}^{\rm \pi} \frac{\textrm{d} \alpha_2}{2{\rm \pi}} |\phi_{n}|^2 = \int \, \textrm{d}^3r \, \textrm{d}^3v \frac{e^2 |\phi_{n}|^2}{T} f_0. \end{equation}

The velocity space integration is straightforward,

(5.6)\begin{equation} \mathcal{L}_0 = \int \, \textrm{d}^3r \, \textrm{d}^3v \frac{e^2 |\phi_{n}|^2}{T} f_0 = \int \, \textrm{d}^3r \frac{e^2 n_0 |\phi_{n}|^2}{T}, \end{equation}

so all that is left is the spatial integration. Because we use toroidal coordinates, the differential volume element is

(5.7)\begin{equation} \textrm{d}^3r = r R_0 (1 + \epsilon \cos(\theta)) \, \textrm{d} r \, \textrm{d} \theta \, \textrm{d} \varphi. \end{equation}

We proceed to calculating $\phi _{n}$ using the poloidal harmonic expansion as detailed in § 4,

(5.8)\begin{equation} \phi(r, \theta, \varphi) = \sum_{m, n} \phi_{m,n}(r - r_0) \, \textrm{e}^{{\rm i} (m \theta + n \varphi)}. \end{equation}

When we examined the trapped Fourier modes, we already calculated $\phi _{n}$. We simply need to generalize it for passing particles as well, which results in

(5.9)\begin{equation} \phi_{n} = \sum_m \phi_{m, n}(r - r_0) \, \textrm{e}^{{\rm i} (m \theta + n q(\bar{r}) \tilde{\theta} + n \tilde{\varphi})}. \end{equation}

As before, we only keep the poloidal harmonic corresponding to $m_0 = -n q_0$ and expand the limits of integration for $r$ to compensate. The result is

(5.10)\begin{equation} |\phi_{n}|^2 = |\phi_{m_0, n}(r-r_0)|^2. \end{equation}

Because the integrand in the adiabatic functional now only depends on $r$, the integral simplifies to

(5.11)\begin{equation} \mathcal{L}_0 = \int_{-\infty}^{\infty} (2 {\rm \pi}R_0) 2{\rm \pi} r \, \textrm{d} r \frac{e^2 n_0}{T} |\phi_{m_0, n}(r-r_0)|^2 \approx \int {\textrm{d} x} R_0 r_0 (2 {\rm \pi})^2 \frac{e^2 n_0}{T} |\phi_{m_0, n}(x)|^2. \end{equation}

Here, we make use of the localization approximation which transforms the factor of $r$ in the integrand into $r_0$. Due to the Gaussian structure of $\phi _{m_0, n}$, this integral is easily performed and we find that

(5.12)\begin{equation} \mathcal{L}_0 = 4 {\rm \pi}^2 R_0 r_0 |\phi_0|^2 |w|^2 \exp\left(\frac{\textrm{Im}(x_0)^2} {\textrm{Re}(w^2)}\right) \sqrt{\frac{\rm \pi}{\textrm{Re}(w^2)}}. \end{equation}

Now that we have calculated the adiabatic functional, we next calculate the trapped functional.

6. Trapped functional

The trapped part of the dispersion relation reads

(6.1)\begin{equation} \mathcal{L}_{\text{trapped}} = \sum_{n_1, n_2} \int \, \textrm{d}^3r \, \textrm{d}^3v \frac{e^2}{T} f_0 \left(\frac{\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_\ast{+} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_E - \omega}{\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\varOmega} - \omega}\right) |\phi_{\boldsymbol{n}}|^2. \end{equation}

We emphasize that although this aspect of the derivation is collisionless, QuaLiKiz includes collisions for trapped electrons. Strictly speaking, this section concerns trapped ions. The majority of the derivation remains the same for trapped electrons, the key difference being that the eventual integral over the particle energy cannot be analytically simplified.

The first step is to determine the appropriate variables to integrate over. For the spatial variables, we once again use toroidal coordinates,

(6.2)\begin{equation} \textrm{d}^3r = R_0 \bar{r} (1 + \epsilon \cos(\theta)) \, \textrm{d} \bar{r} \, \textrm{d} \theta \, \textrm{d} \varphi. \end{equation}

For velocity space, we use the variables $(v, \lambda , v_\phi )$ which correspond to the speed $v$, pitch angle parameter $\lambda$ and cylindrical velocity phase $v_\phi$. The result is

(6.3)\begin{equation} \textrm{d}^3 v = \sum_{\epsilon_\parallel} v^2 \frac{b}{2\sqrt{1 - \lambda b}} \, \textrm{d} v \, \textrm{d} \lambda \, \textrm{d} v_\phi, \end{equation}

where the sum over $\epsilon _\parallel$ accounts for both possible signs of the parallel velocity. Because the integrand is independent of $\varphi$ or $v_\phi$, we obtain

(6.4)\begin{equation} \textrm{d}^3r \, \textrm{d}^3v = \sum_{\epsilon_\parallel} 2 {\rm \pi}^2 R_0 \bar{r} (1 + \epsilon \cos(\theta)) \, \textrm{d} \bar{r} \, \textrm{d} \theta v^2 \frac{b}{\sqrt{1 - \lambda b}} \, \textrm{d} v \, \textrm{d} \lambda. \end{equation}

It is important to note that the limits of integration depend on the order of integration. For a given $\theta$, the pitch angle parameter $\lambda$ for a trapped particle is bounded by

(6.5)\begin{equation} \frac{1 - \epsilon}{1 + \epsilon} \le \lambda \le \frac{1 + \epsilon \cos(\theta)}{1 + \epsilon}. \end{equation}

The lower bound corresponds to the trapped–passing boundary, while the upper bound corresponds to a particle that has $v_\parallel = 0$ at a given angle $\theta$. We can, however, exchange the order of integration as follows:

(6.6)\begin{equation} \int_{0}^{2 {\rm \pi}} \, \textrm{d} \theta \int_{({1-\epsilon})/({1+\epsilon})}^{({1+ \epsilon \cos(\theta)})/({1 + \epsilon})} \, \textrm{d} \lambda f(\theta, \lambda) = \int_{({1-\epsilon})/({1+\epsilon})}^{1} \, \textrm{d} \lambda \int_{-\theta_b}^{\theta_b} \, \textrm{d} \theta f(\theta, \lambda). \end{equation}

We recall that the definition of a bounce average is

(6.7)\begin{equation} \langle G(\epsilon_\parallel, \theta)\rangle = \frac{\displaystyle \int_{-\theta_b} ^{\theta_b} \, \textrm{d} \theta \frac{G(\epsilon_\parallel, \theta) + G(-\epsilon_\parallel, \theta)}{\sqrt{1-\lambda b}} }{\displaystyle \int_{-\theta_b} ^{\theta_b}\frac{2 \, \textrm{d} \theta}{\sqrt{1-\lambda b}} } = \sum_{\epsilon_\parallel} \frac{\bar{\varOmega}_2}{2 {\rm \pi}} \int_{-\theta_b}^{\theta_b} \, \textrm{d} \theta \frac{G(\epsilon_\parallel, \theta)}{\sqrt{1-\lambda b}}. \end{equation}

By exchanging our limits of integration and integrating over $\theta$ first, part of the trapped functional simplifies to become a bounce average.

Next, we approximate the equilibrium distribution function assuming the Mach number $U_\parallel / c_s$ is small, where $c_s = \sqrt {T/m}$ is the sound speed. We thus note that high-rotation tokamaks will break the assumed ordering in QuaLiKiz. In addition, even if the electron and ion species each have low Mach number, care must be taken with regards to heavy impurities, as the larger mass can lead to a higher Mach number for the same rotation velocity and temperature (Citrin Reference Citrin2017). Using the low-Mach-number ordering, the distribution function then simplifies to

(6.8)\begin{equation} f_0 \approx n_0 \left(\frac{m}{2 {\rm \pi}T}\right)^{3/2} \, \textrm{e}^{-\xi} \left( 1 + \frac{2 v_\parallel U_\parallel}{v_T^2} + \frac{U_\parallel^2}{v_T^2} \left(\frac{2 v_\parallel^2}{v_T^2} - 1\right) \right). \end{equation}

Moreover, because $|\varOmega _1|, |\varOmega _2| \gg |\omega |$, we can approximate this integral by truncating the sum at $n_1 = n_2 = 0$. We also perform a change of variables from $v$ to $\xi$ to obtain

(6.9)\begin{align} \mathcal{L}_{\text{trapped}} &= \int \, \textrm{d} \bar{r} \, \textrm{d} \lambda\, \textrm{d} \xi (2 {\rm \pi})^2 \frac{n_0 e^2}{T} \frac{R_0 \bar{r}}{\sqrt{\rm \pi}} \frac{\sqrt{\xi} \, \textrm{e}^{-\xi}}{\bar{\varOmega}_2} \nonumber\\ &\quad \times \frac{\displaystyle \left\langle (\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_\ast{+} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_E - \omega ) \left( 1 + \frac{2 v_\parallel U_\parallel}{v_T^2} + \frac{U_\parallel^2}{v_T^2} \left(\frac{2 v_\parallel^2}{v_T^2} - 1\right) \right) \right\rangle }{n \varOmega_3 - \omega} |\phi_{0,0,n}|^2. \end{align}

The Bessel functions from the gyromotion and the banana motion are implicit in $\phi _{0,0,n}$. We next simplify the partial derivatives with respect to $\boldsymbol {J}$. Because $n_1 = n_2 = 0$, we only keep the partial derivative with respect to $J_3$. Knowing that $J_3 = J_3(\bar {r})$, we perform a change in variables from $J_3$ to $\bar {r}$ and find that

(6.10)\begin{equation} \frac{\partial g}{\partial J_3} = \frac{\partial g}{\partial \bar{r}} \frac{\textrm{d} \bar{r}}{d J_3} \approx \frac{-R_0 \omega_{d0}}{T} \frac{\partial g}{\partial \bar{r}}, \end{equation}

where $g$ is a generic scalar function. We then define the following normalized gradients:

(6.11)\begin{gather} A_n ={-} \frac{R_0}{n} \frac{\textrm{d} n}{\textrm{d}\bar{r}}, \end{gather}

(6.12)\begin{gather}A_T ={-} \frac{R_0}{T} \frac{\textrm{d} T}{\textrm{d} \bar{r}}, \end{gather}

(6.13)\begin{gather}A_U ={-} \frac{R_0}{v_T} \frac{\textrm{d} U_{{\parallel}}}{\textrm{d} \bar{r}}. \end{gather}

To perform the bounce average, we note that only $v_\parallel$ is dependent on $\theta$. We perform the calculation explicitly to find that

(6.14)\begin{gather} \langle v_\parallel\rangle = \langle v_\parallel^3\rangle =0 , \end{gather}

(6.15)\begin{gather}\langle v_\parallel^2\rangle = \frac{4 E \epsilon}{m} \frac{({E}(\kappa) - (1 - \kappa^2) {K}(\kappa))}{{K}(\kappa)} = v_T^2 \xi H(\kappa), \end{gather}

where we define

(6.16)\begin{equation} H(\kappa)= \frac{2 \epsilon ({E}(\kappa) - (1 - \kappa^2) {K}(\kappa))}{{K}(\kappa)}. \end{equation}

To simplify our expressions, we also fold $\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {\omega }_E$ into the mode frequency such that

(6.17)\begin{equation} \bar{\omega} = \omega - \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_E \approx \omega - n \omega_E \approx \omega - n \omega_{E0} - n \omega_E' x, \end{equation}

where, as discussed earlier, we Taylor expand $\omega _E$ about $x = 0$ and $\partial _r \omega _E = \omega _E'$ is related to the radial electric shear. Rather than including $x$ fully, we instead approximate the term by averaging it over the Gaussian eigenfunctions:

(6.18)\begin{equation} \langle x\rangle _r = \frac{\displaystyle \int_{-\infty}^{\infty} {\textrm{d} x} |\phi_{m_0, n}(x)|^2 x }{\displaystyle \int_{-\infty}^{\infty} {\textrm{d} x} |\phi_{m_0, n}(x)|^2} = \textrm{Re}(x_0) + \frac{\textrm{Im}(x_0) \textrm{Im}(w^2)}{\textrm{Re}(w^2)}. \end{equation}

We then obtain

(6.19)\begin{equation} \bar{\omega} \approx \omega - n \omega_{E0} - n \omega_E' \langle x\rangle _r. \end{equation}

Ignoring all terms that are order cubic or higher with the Mach number, we then find that

(6.20)\begin{equation} \left\langle ( \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_\ast{+} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_E - \omega ) \left(1 + \frac{2 v_\parallel U_\parallel}{v_T^2} + \frac{U_\parallel^2}{v_T^2} \left( \frac{2 v_\parallel^2}{v_T^2} - 1\right)\right) \right\rangle = n\omega_{d0}( \mathcal{A}_t + \mathcal{B}_t \xi + \mathcal{C}_t \xi^2 ), \end{equation}

where

(6.21)\begin{align} \mathcal{A}_t & = \left(1 - \frac{U_\parallel^2}{v_T^2}\right) \left(A_n - \frac{3}{2} A_T - z^2 F_d (\kappa) \right) - \frac{U_\parallel}{v_T} \left(2 A_U - \frac{U_\parallel}{v_T} A_T\right), \end{align}

(6.22)\begin{align} \mathcal{B}_t &= \left(1 - \frac{U_\parallel^2}{v_T^2}\right) A_T + 4 A_U \frac{U_\parallel}{v_T} H(\kappa) + \frac{U_\parallel^2}{v_T^2} H(\kappa) (2 A_n - 7 A_T - 2 z^2 F_d (\kappa) ) \nonumber\\ & \approx \left(1 -\frac{U_\parallel^2}{v_T^2}\right) A_T + 4 A_U \frac{U_\parallel}{v_T} H(\kappa), \end{align}

(6.23)\begin{align} \mathcal{C}_t &= 2 A_T \frac{U_\parallel^2}{v_T^2} H(\kappa) \approx 0. \end{align}

Here, we have defined

(6.24)\begin{equation} z^2 = \frac{\bar{\omega}}{n \omega_{d0} F_d (\kappa)}. \end{equation}

Moreover, we take note that $H(\kappa ) \sim {O}(\epsilon )$; since the inverse aspect ratio $\epsilon$ is small, we can safely ignore all terms proportional to $U_\parallel ^2 H(\kappa ) / v_T^2$.

Substituting the above into the integrand, we obtain

(6.25)\begin{equation} \mathcal{L}_\text{trapped} = \int \, \textrm{d} \bar{r} \, \textrm{d} \lambda\, \textrm{d} \xi (2 {\rm \pi})^2 \frac{n_0 e^2}{T} \frac{R_0 \bar{r}}{\sqrt{\rm \pi}} \frac{\sqrt{\xi} \, \textrm{e}^{-\xi}}{\bar{\varOmega}_2} |\phi_{0,0,n}|^2 \frac{\mathcal{A}_t + \mathcal{B}_t \xi}{F_d (\kappa)(\xi - z^2)}. \end{equation}

Due to the localization of the mode, we evaluate any functions of $\bar {r}$ at $r_0$ in the above expression aside from the electrostatic potential. We then rewrite the trapped functional as

(6.26)\begin{align} \mathcal{L}_\text{trapped} = \int \, \textrm{d} \bar{r} \, \textrm{d} \lambda\, \textrm{d} \xi (2 {\rm \pi})^2 \frac{n_0 e^2}{T} \frac{R_0 r_0}{\sqrt{\rm \pi}} \frac{\sqrt{\xi} \, \textrm{e}^{-\xi}}{\bar{\varOmega}_2} |J_0(k_\perp \rho) J_0 (k_r \delta_b) \boldsymbol{\cdot} \phi_{m_0, n}|^2 \frac{\mathcal{A}_t + \mathcal{B}_t \xi}{F_d (\kappa)(\xi - z^2)}. \end{align}

The gyromotion and bounce motion appear in two separate Bessel functions. Since the only explicit radial dependence is contained in the electrostatic potential, we can change variables using Plancherel's theorem to integrate over $k_r$,

(6.27)\begin{equation} \int_{-\infty}^{\infty} {\textrm{d} x} f(x) g(x)^\ast{=} \int_{-\infty}^{\infty} \frac{\textrm{d} k_r}{2{\rm \pi}} \hat{f}(k_r) \hat{g}(k_r)^\ast. \end{equation}

After transforming to Fourier space, we treat the Bessel functions as normal scalar functions instead of differential operators. We next note that the Bessel functions are dependent on velocity through the gyroradius and banana width,

(6.28)\begin{gather} \rho = \frac{v_\perp}{\varOmega_1}, \end{gather}

(6.29)\begin{gather}\delta_b \approx \frac{q}{\sqrt{\epsilon}} \rho. \end{gather}

We approximate this energy dependence by averaging each Bessel function separately over velocity space using a Maxwellian distribution. Doing so allows us to retain finite Larmor radius and finite banana width effects while also making the energy and pitch angle integration tractable. We find that

(6.30)\begin{equation} \frac{\displaystyle \int \, \textrm{d}^3v J_0(k_\perp \rho)^2 f_0}{\displaystyle \int \, \textrm{d}^3v f_0} = \textrm{e}^{-{k_\perp^2 \rho_{\text{th}}^2}/{2}} I_0\left(\frac{k_\perp^2 \rho_{\text{th}}^2}{2}\right) = \varGamma_0 (k_\perp \rho_{\text{th}}), \end{equation}

where $I_0$ is a modified Bessel function of the first kind and the characteristic thermal gyroradius $\rho _\text {th}$ is defined as

(6.31)\begin{equation} \rho_{\text{th}} = \frac{\sqrt{2T/m}}{\varOmega_1}. \end{equation}

Similarly, for the average over the banana orbit we obtain

(6.32)\begin{equation} \frac{\displaystyle \int \, \textrm{d}^3v J_0(k_\perp \rho)^2 f_0}{\displaystyle \int \, \textrm{d}^3v f_0} = \textrm{e}^{-{k_\perp^2 \delta_{b,\text{th}}^2}/{2}} I_0\left(\frac{k_\perp^2 \delta_{b,\text{th}}^2}{2}\right) = \varGamma_0 (k_\perp \delta_{b,\text{th}}), \end{equation}

where the thermal banana width is

(6.33)\begin{equation} \delta_{b, \text{th}} = \frac{q}{\sqrt{\epsilon}} \rho_\text{th}. \end{equation}

Note that $k_\perp ^2$ is written as

(6.34)\begin{equation} k_\perp^2 = k_r^2 + k_\theta^2 = k_r^2 + \frac{n^2 q_0^2}{r_0^2}, \end{equation}

where we evaluate $k_\theta$ at $r_0$. That $k_\theta ^2 = n^2q^2/r^2$ comes from differentiating with respect to $\theta$ in the ballooning expansion due to the eikonal term. Because the $k_r$ dependence is now completely separable from the $\kappa$ and $\xi$ dependence, we write the trapped functional as

(6.35)\begin{align} \mathcal{L}_\text{trapped} &= \int \, \textrm{d} \lambda\, \textrm{d} \xi (2 {\rm \pi})^2 \frac{n_0 e^2}{T} \frac{R_0 r_0}{\sqrt{\rm \pi}} \frac{\sqrt{\xi} \, \textrm{e}^{-\xi}}{\bar{\varOmega}_2} \frac{\mathcal{A}_t + \mathcal{B}_t \xi} {F_d (\kappa)(\xi - z^2)} \nonumber\\ & \quad \times \int_{-\infty}^{\infty} \frac{\textrm{d} k_r}{2 {\rm \pi}} \varGamma_0 (k_\perp \rho_\text{th}) \varGamma_0 (k_r \delta_{b,\text{th}})| \textrm{d}\hat{\phi}_n (k_r d) |^2, \end{align}

where $\hat {\phi }_n$ is computed using a Fourier transform:

(6.36)\begin{equation} \hat{\phi}_n (k_r d) = \sqrt{2 {\rm \pi}} w \phi_0 \, \textrm{e}^{-({k_r^2 w^2}/{2})} \, \textrm{e}^{i k_r x_0}. \end{equation}

We next simplify the integral over $\xi$, which is of the form

(6.37)\begin{equation} \int_{0}^{\infty}\, \textrm{d} \xi \frac{\sqrt{\xi}}{\sqrt{\rm \pi}} \frac{\mathcal{A}_t + \mathcal{B} \xi}{\xi - z^2} \, \textrm{e}^{-\xi} = \int_{-\infty}^{\infty} \, \textrm{d} u \frac{u^2}{\sqrt{\rm \pi}} \frac{\mathcal{A}_t + \mathcal{B}_t u^2}{(u + z)(u - z)} \, \textrm{e}^{{-}u^2}, \end{equation}

where we performed the change of variables $\xi = u^2$. Using the plasma dispersion function detailed in Appendix A, this simplifies to

(6.38)\begin{align} \int_{-\infty}^{\infty} \, \textrm{d} u \frac{u^2}{\sqrt{\rm \pi}} \frac{\mathcal{A}_t + \mathcal{B}_t u^2}{(u+z)(u-z)} \, \textrm{e}^{{-}u^2} = \mathcal{A}_t G_2(z, -z) + \mathcal{B}_t G_4(z, -z) = \frac{\mathcal{A}_t Z_2 (z) + \mathcal{B}_t Z_4 (z)}{z}, \end{align}

where the final simplification is made using the fact that $Z_{2n}$ is an even function for $n \ge 0$. Meanwhile, we rewrite the integration over $\lambda$ with a change in variables,

(6.39)\begin{equation} \frac{\textrm{d} \lambda}{\bar{\varOmega}_2} = 4 f_t {K}(\kappa) \kappa \, \textrm{d} \kappa, \end{equation}

where we utilize the transformation

(6.40)\begin{equation} \lambda \approx 1 - 2 \epsilon \kappa^2, \end{equation}

and define the flux surface averaged trapped particle fraction

(6.41)\begin{equation} f_t = \frac{2 \sqrt{2 \epsilon}}{\rm \pi}. \end{equation}

Thus, the trapped functional simplifies to

(6.42)\begin{equation} \mathcal{L}_\text{trapped} = (2 {\rm \pi})^3 \frac{e^2 n_0}{T} r_0 R_0 f_t \langle \mathcal{I}_t\rangle _{\xi, \kappa} \langle \varGamma_0 (k_\perp \rho_{\text{th}}) \varGamma_0 (k_r \delta_{b, \text{th}}) |\, \textrm{d} \hat{\phi}(k_r d)|^2\rangle_{k_r}, \end{equation}

where

(6.43)\begin{equation} \langle \mathcal{I}_t\rangle _{\xi, \kappa} = \frac{2}{\rm \pi} \int_0^1 \, \textrm{d} \kappa \frac{{K}(\kappa) \kappa}{z F_d (\kappa)} (\mathcal{A}_t Z_2(z) + \mathcal{B}_t Z_4 (z)), \end{equation}

and

(6.44)\begin{equation} \langle \varGamma_0 (k_\perp \rho_{\text{th}}) \varGamma_0 (k_r \delta_{b, \text{th}}) | \textrm{d} \hat{\phi}(k_r d)|^2\rangle _{k_r} = \int_{-\infty}^{\infty} \frac{\textrm{d} k_r}{2 {\rm \pi}} \varGamma_0 (k_\perp \rho_{\text{th}}) \varGamma_0 (k_r \delta_{b, \text{th}}) | \textrm{d} \hat{\phi}(k_r d)|^2. \end{equation}

The remaining integrals are to be calculated numerically, where we note that $z$ is a function of both $\bar {\omega }$ and $\kappa$. Thus, the trapped functional is the product of two separate 1-D integrals, one of which is $\omega$ independent; we therefore characterize the trapped functional as a 1-D integral that must be calculated numerically. Now that we have simplified the expression for the trapped functional, we turn to calculating the passing functional.

7. Passing functional

The passing part of the dispersion relation reads

(7.1)\begin{equation} \mathcal{L}_{\text{passing}} = \sum_{n_1, n_2} \int \, \textrm{d}^3r \, \textrm{d}^3v \frac{e^2}{T} f_0 \left(\frac{\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_\ast{+} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_E - \omega}{\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\varOmega} - \omega}\right) |\phi_{\boldsymbol{n}}|^2. \end{equation}

We reuse many of the same arguments in § 6 regarding changes in variables and approximating the equilibrium distribution function. One key difference is that instead of the bounce average, we use the transit average

(7.2)\begin{equation} \langle G(\epsilon_\parallel, \theta) \rangle = \frac{\displaystyle \int_{-{\rm \pi}}^{\rm \pi} \, \textrm{d} \theta \frac{G(\epsilon_\parallel, \theta)}{\sqrt{1 - \lambda b}}}{\displaystyle \int_{-{\rm \pi}}^{\rm \pi} \frac{ \textrm{d} \theta}{\sqrt{1 - \lambda b}}} = \frac{\bar{\varOmega}_2}{2 {\rm \pi}} \int_{-{\rm \pi}}^{\rm \pi} \, \textrm{d} \theta \frac{G(\epsilon_\parallel, \theta)}{\sqrt{1 - \lambda b}}. \end{equation}

We note here that the bounce angle $\theta _b$ is set to ${\rm \pi}$ and that we do not perform a sum over $\epsilon _\parallel$ to compute the transit average. Moreover, the integration bounds for $\lambda$ are such that

(7.3)\begin{equation} 0 \le \lambda \le \frac{1 - \epsilon}{1+\epsilon}. \end{equation}

These bounds hold regardless of whether we integrate over $\theta$ before or after integrating over $\lambda$. Since they are independent of $\theta$, the order of integration of the two variables can be freely interchanged. As in the trapped case, we only keep $n_1 = 0$ since $|\varOmega _1| \gg |\omega |$. As discussed in § 4, $n_2$ refers to the poloidal harmonic. We keep only $n_2 = m_0$ and use the approximation that

(7.4)\begin{equation} m_0 + n q(\bar{r}) \approx \frac{x}{d}. \end{equation}

In the resonant denominator we then obtain

(7.5)\begin{equation} n \varOmega_3 + m_0 \varOmega_2 - \omega \approx n \omega_{d0} \xi F + \frac{x}{d} \varOmega_2 + n \omega_{E0} + n \omega_{E}' x - \omega, \end{equation}

where we also expand $\omega _E$ about $x = 0$. The passing functional is then

(7.6)\begin{align} \mathcal{L}_{\text{passing}} &= \sum_{\epsilon_\parallel} \int \, \textrm{d} \bar{r} \, \textrm{d} \lambda\, \textrm{d} \xi (2 {\rm \pi})^2 \frac{n_0 e^2}{T} \frac{R_0 \bar{r}}{\sqrt{\rm \pi}} \frac{\sqrt{\xi} \, \textrm{e}^{-\xi}}{\bar{\varOmega}_2} \nonumber\\ & \quad \times \frac{\displaystyle \left\langle (\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_\ast{+} n \omega_{E0} + n \omega_E' x - \omega ) \left( 1 + \frac{2 v_\parallel U_\parallel}{v_T^2} + \frac{U_\parallel^2}{v_T^2} \left(\frac{2 v_\parallel^2}{v_T^2} - 1\right) \right) \right\rangle }{\displaystyle n \omega_{d0} \xi F + \frac{x}{d} \varOmega_2 + n \omega_{E0} + n \omega_{E}' x - \omega } |\phi_{0,m_0,n}|^2. \end{align}

Here, we have evaluated all functions at $r = r_0$ except for the terms proportional to $x$ in the resonant denominator and numerator. These terms must be kept if we wish to take into account the effects of the poloidal motion as well as the radial electric field shear. We now evaluate the integration over $\bar {r}$ while leaving aside the term proportional to $x$ in the numerator.

To proceed, we use Plancherel's theorem to integrate over $\alpha _2$ instead of $\bar {r}$,

(7.7)\begin{equation} \int_{-\infty}^{\infty} \, \textrm{d} \bar{r} f(\bar{r}) g(\bar{r})^\ast{=} \int_{-\infty}^{\infty} \frac{\textrm{d} \alpha_2}{2{\rm \pi}} |d| \hat{f}(\alpha_2) \hat{g}(\alpha_2)^\ast. \end{equation}

For convenience, we compute the radial integral in isolation and relabel variables,

(7.8)\begin{equation} \int_{-\infty}^{\infty} \, \textrm{d} \bar{r} \frac{1}{a \dfrac{x}{d} - b} |\phi_{0, m_0, n}|^2 = \int_{-\infty}^{\infty} \frac{\textrm{d} \alpha_2}{2 {\rm \pi}} |d| \mathcal{F}\left(\frac{\phi_{0, m_0, n} }{a \dfrac{x}{d} - b} \right) \mathcal{F}( \phi_{0, m_0, n})^\ast. \end{equation}

We calculated in § 4 that

(7.9)\begin{equation} \mathcal{F}(\phi_{0, m_0, n}) = J_0(k_\perp(\alpha_2) \rho) \hat{\phi}_n (\theta(\alpha_2)) \, \textrm{e}^{{\rm i} n \tilde{\varphi}(\alpha_2) - {\rm i} ({\tilde{r}(\alpha_2)}/{d}) \theta(\alpha_2)}. \end{equation}

We note that $k_\perp$ is defined such that

(7.10)\begin{equation} k_\perp(\alpha_2)^2 = \frac{\theta(\alpha_2)^2}{d^2} + \frac{n^2 q_0^2}{r_0^2}. \end{equation}

We next use the convolution theorem to calculate the other Fourier transform,

(7.11)\begin{equation} \mathcal{F} \left( \frac{\phi_{0, m_0, n} }{a \dfrac{x}{d} - b}\right) = \frac{1}{2 {\rm \pi}} \mathcal{F}( \phi_{0, m_0, n} ) \ast \mathcal{F}\left( \frac{1 }{a \dfrac{x}{d} - b}\right). \end{equation}

Computing the Fourier transform of both functions and performing the convolution, we find that

(7.12)\begin{align} \mathcal{F} \left( \phi_{0, m_0, n} \frac{1}{1 + a \dfrac{x}{d} - b}\right) &= \int_{-\infty}^{\infty} \, \textrm{d} \alpha_2' \frac{{\rm i}}{|a|} \varTheta\left(\frac{\alpha_2 - \alpha_2'}{a}\right) J_0(k_\perp(\alpha_2') \rho) \hat{\phi}_n (\theta(\alpha_2')) \nonumber\\ &\quad \times \textrm{e}^{({-{\rm i}b}/{a}) (\alpha_2 ' - \alpha_2)} \, \textrm{e}^{{\rm i} n \tilde{\varphi}(\alpha_2') - {\rm i} ({\tilde{r}(\alpha_2)}/{d}) \theta(\alpha_2')}, \end{align}

where $\varTheta$ is the Heaviside step function. Here, we have assumed that $\textrm {Im}(b) > 0$. This is justified since $b \sim \omega$ and we are only interested in positive growth rates. Combining the results, we obtain for the passing integral $I_{p,r}$ that

(7.13)\begin{align} I_{p,r} &= \int \, \textrm{d}\bar{r} \frac{1}{a \dfrac{x}{d} - b } |\phi_{0,m_0,n}|^2 = \int \frac{\textrm{d} \alpha_2 \, \textrm{d} \alpha_2'}{2{\rm \pi}} \frac{{\rm i} |d|}{|a|} \varTheta\left(\frac{\alpha_2 - \alpha_2'}{a}\right) J_0^\ast (k_\perp(\alpha_2) \rho) J_0(k_\perp(\alpha_2') \rho)\nonumber\\ &\quad \times \hat{\phi}_n(\theta(\alpha_2))^\ast \hat{\phi}_n(\theta(\alpha_2')) \, \textrm{e}^{\varLambda(\alpha_2) - \varLambda(\alpha_2')}, \end{align}

where

(7.14)\begin{equation} \varLambda(\alpha_2) = {\rm i} \left( \frac{b}{a} \alpha_2 - n \tilde{\varphi}(\alpha_2) + \frac{\tilde{r}(\alpha_2)}{d} \theta(\alpha_2) \right). \end{equation}

We then substitute in

(7.15)\begin{gather} a = \varOmega_2 + n d \omega_E', \end{gather}

(7.16)\begin{gather}b = \omega - n \omega_{d0} \xi F - n \omega_{E0} \end{gather}

and rewrite the eikonal term to obtain

(7.17)\begin{equation} \varLambda(\alpha_2) = {\rm i} \left( \frac{\omega - n \omega_{d0} \xi F - n \omega_{E0}}{\varOmega_2 + n d \omega_E'} \alpha_2 - n \tilde{\varphi}(\alpha_2) + \frac{\tilde{r}(\alpha_2)}{d} \theta(\alpha_2) \right). \end{equation}

It is important to recognize the physical importance of $\varLambda$. In the ballooning representation, we encoded a certain particle trajectory in the eikonal that differs from the magnetic drift trajectory. The function $\varLambda$ encapsulates the phase difference between these two trajectories.

Before proceeding, we must recognize that integrating over $\alpha _2$ and $\alpha _2'$ is inconvenient. The function $\hat {\phi }_n$ has Gaussian structure in $\theta$, but not in $\alpha _2$. Thus, the next goal is to write the integrand in terms of $\theta$ and $\theta '$. First, we introduce new variables,

(7.18)\begin{gather} \theta_+ = \frac{\theta + \theta'}{2}, \end{gather}

(7.19)\begin{gather}\theta_- = \theta - \theta'. \end{gather}

We next Taylor expand the terms in the exponential around $\theta _+$ to find that

(7.20)\begin{equation} \varLambda(\alpha_2) - \varLambda(\alpha_2') \approx \theta_- \frac{\textrm{d} \alpha_2}{\textrm{d} \theta}(\theta = \theta_+ ) \varLambda^\prime (\theta = \theta_+ ), \end{equation}

where $\varLambda '$ denotes the derivative of $\varLambda$ with respect to $\alpha _2$. Due to the rapidly varying phase in the exponential, the factor of $i \omega \sim -\gamma$ in the exponential and the Gaussian integrand, we ignore higher-order terms to obtain the dominant contribution. From the equations listed in § 2, we find

(7.21)\begin{equation} \frac{\textrm{d}}{\textrm{d} \alpha_2} \left(n \tilde{\varphi} - \frac{\tilde{r}}{d} \theta\right) = \frac{n}{\varOmega_2} \left(\boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla} \varphi - q \boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla} \theta - \theta \frac{\textrm{d} q}{\textrm{d} \psi} \boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla} \psi - \varOmega_d \right). \end{equation}

The leading terms can be computed explicitly in much the same manner as when calculating the magnetic precession frequency,

(7.22)\begin{align} \bar{\varOmega}_d &= \boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla} \varphi - q \boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla} \theta - \theta \frac{\textrm{d} q}{\textrm{d} \psi} \boldsymbol{v}_D \boldsymbol{\cdot} \boldsymbol{\nabla} \psi \nonumber\\ &= \omega_{d0} \xi (2 - \lambda b ) (\cos(\theta) + (s \theta - \alpha_M \sin(\theta) \sin^2(\theta) )) + \frac{q E_r}{r B}. \end{align}

Although somewhat similar to the magnetic drift frequency proper, there are two key differences. First, the magnetic shear term is different and proportional to $\theta \sin (\theta )$. Second, this frequency is explicitly $\theta$-dependent and no bounce-transit average is performed. In carrying out the calculation the bounce-averaged magnetic drift terms partially cancel; for sufficiently small radial electric field shear we obtain

(7.23)\begin{equation} \varLambda(\alpha_2) - \varLambda(\alpha_2') \approx{-}i \theta_- \left(\frac{n \bar{\varOmega}_d (\theta_+ ) - \omega}{\varOmega_2 + n d \omega_E'}\right), \end{equation}

where

(7.24)\begin{equation} n \bar{\varOmega}_d \approx n \omega_{d0} \xi (2 - \lambda b ) (\cos(\theta) + (s \theta - \alpha_M \sin(\theta) \sin^2(\theta) )) + n \omega_{E0}. \end{equation}

We next change the variables of integration from $\alpha _2, \alpha _2'$ to $\theta , \theta '$ using

(7.25)\begin{equation} \frac{\textrm{d} \theta}{\textrm{d} \alpha_2} = \frac{\sqrt{1 - \lambda b}}{\bar{\varOmega}_2}. \end{equation}

The integral then becomes

(7.26)\begin{align} I_{p,r} &= \int \frac{\textrm{d} \theta \, \textrm{d} \theta'}{2{\rm \pi}} \frac{{\rm i} |d|}{|\varOmega_2 + n \,d \omega_E' |} \varTheta\left(\frac{\theta_-}{\varOmega_2 + n \,d \omega_E' }\right) J_0 \hat{\phi}_n(\theta)^\ast J_0' \hat{\phi}_n(\theta') \nonumber\\ &\quad \times \frac{\bar{\varOmega}_2}{\sqrt{1 - \lambda b(\theta)}} \frac{\bar{\varOmega}_2}{\sqrt{1 - \lambda b(\theta')}} \, \textrm{e}^{-{\rm i} \theta_- ({\bar{\varOmega}_2}/{\sqrt{1 - \lambda b(\theta_+)}}) (({n \bar{\varOmega}_d (\theta_+) - \omega})/({\varOmega_2 + n d \omega_E'}))}, \end{align}

where the Bessel functions are evaluated in terms of $\theta$ and $\theta '$. We now substitute in an expression for $\hat {\phi }_n$ in terms of a Fourier transform to obtain

(7.27)\begin{align} I_{p,r}& = \int \frac{\textrm{d} \theta \, \textrm{d} \theta' {\textrm{d} x} {\textrm{d} x}'}{2{\rm \pi} } \frac{{\rm i} \bar{\varOmega}_2^2}{|\varOmega_2 + n d \omega_E' | |d|} \varTheta\left(\frac{\theta_-}{\varOmega_2 + n \, \textrm{d} \omega_E'}\right) J_0 \phi^\ast_{m_0,n}(x) J_0' \phi_{m_0,n}(x') \nonumber\\ &\quad \times \frac{\displaystyle \textrm{e}^{-{\rm i} \theta_- ({\bar{\varOmega}_2}/{\sqrt{1 - \lambda b(\theta_+)}}) ({n \bar{\varOmega}_d(\theta_+) - \omega}/({\varOmega_2 + n d \omega_E'}))} \, \textrm{e}^{- {\rm i} ({\theta x}/{d})} \, \textrm{e}^{{\rm i} ({\theta' x'}/{d})}}{\sqrt{1 - \lambda b(\theta)} \sqrt{1 - \lambda b(\theta')}}. \end{align}

We then make the following substitutions:

(7.28)\begin{gather} x_+ = \frac{x + x'}{2}, \end{gather}

(7.29)\begin{gather}x_- = x - x', \end{gather}

(7.30)\begin{gather}k_+ = \frac{\theta_+}{|d|}, \end{gather}

(7.31)\begin{gather}k_- = \frac{\theta_-}{|d|}, \end{gather}

(7.32)\begin{gather}\textrm{d} \theta \, \textrm{d} \theta' \,{\textrm{d} x} \,{\textrm{d} x}' = \textrm{d} k_+ \textrm{d} k_- \textrm{d} x_+ \textrm{d} x_- |d|^2, \end{gather}

(7.33)\begin{gather}\frac{\theta'x'}{d} - \frac{\theta x}{d} = (- k_- x_+{-} k_+ x_- ) \frac{|d|}{d}, \end{gather}

to obtain

(7.34)\begin{align} I_{p,r} &= \int \frac{\textrm{d} k_+ \textrm{d} k_- \textrm{d} x_+ \textrm{d} x_-}{2{\rm \pi}} \frac{{\rm i} |d| (\bar{\varOmega}_2)^2}{|\varOmega_2 + n d \omega_E'|} \varTheta\left(\frac{k_-}{\varOmega_2 + n d \omega_E'}\right) J_0 \phi^\ast_{m_0,n}\left(x_+{+} \frac{x_-}{2}\right) J_0' \phi_{m_0,n}\left(x_+{-} \frac{x_-}{2}\right)\nonumber\\ &\quad \times \frac{\displaystyle \exp \left( \frac{-{\rm i} k_- |d|\bar{\varOmega}_2}{\sqrt{1 - \lambda b(k_+ |d|)}} \left(\frac{n \bar{\varOmega}_d(k_+ |d|) - \omega}{\varOmega_2 + n d \omega_E'}\right)\right) \, \textrm{e}^{-{\rm i} (k_- x_+{+} k_+ x_- ) ({|d|}/{d})}}{\displaystyle \sqrt{1 - \lambda b\left(k_+ |d| + \frac{k_- |d| }{2} \right)}\sqrt{1 - \lambda b\left(k_+ |d| - \frac{k_- |d| }{2} \right)}} . \end{align}

At first glance, it seems like we have only made the derivation more difficult. We are now performing a four-dimensional integration over variables which do not have a convenient Gaussian structure. Fortunately, this simplifies. First, we notice that the integration over $k_-$ can be done via an integration by parts procedure. In general, for a complex parameter $c$ we obtain

(7.35)\begin{align} \int_{a}^{b} \, \textrm{d} s g (s) \, \textrm{e}^{{\rm i} c s} &= \int_{a}^{b} \, \textrm{d} s \frac{1}{ik} \left(\frac{\textrm{d}}{\textrm{d} s} (g(s) \, \textrm{e}^{{\rm i} c s}) - \frac{\textrm{d} g}{\textrm{d} s} \, \textrm{e}^{{\rm i} c s}\right) \nonumber\\ & = \left. \frac{g(s) \, \textrm{e}^{{\rm i} c s}}{{\rm i} c} \right|_{a}^{b} - \int_a^b \, \textrm{d} s \frac{1}{({\rm i} c )^2} \left( \frac{\textrm{d}}{\textrm{d} s} \left(\frac{\textrm{d} g}{\textrm{d} s} \, \textrm{e}^{{\rm i} c s}\right) - \frac{d^2 g}{\textrm{d} s^2} \, \textrm{e}^{{\rm i} cs}\right).\nonumber\\ & \approx \sum_{m=0}^{N} \frac{({-}1)^m}{({\rm i} c)^{m+1}}\left[ \frac{\textrm{d}^m g}{\textrm{d} s^m} (b) - \frac{\textrm{d}^m g}{\textrm{d} s^m} (a)\right]. \end{align}

This is the asymptotic expansion for sufficiently large $c$. We apply a similar expansion to the integral over $k_-$ and keep only the first term. Because $\textrm {Im}(\omega ) > 0$ and the integrand contains a Heaviside step function, the first term is guaranteed to converge. Note that we would normally need to apply the method of steepest descent to properly approximate the integral; however, this requires that the term in the exponential have a saddle point somewhere in the complex plane. Due to our previous approximation, the term in the exponential is monotonic in $k_-$, thus the method of steepest descent is not necessary for our purposes. We find then that

(7.36)\begin{align} &\int_{-\infty}^{\infty} \, \textrm{d} k_- \frac{{\rm i} |d|}{ |\varOmega_2 + n d \omega_E'|} g( k_- ) \, \textrm{e}^{- {\rm i} k_- ({|d|}/{d})} \exp\left( \frac{-{\rm i} k_- |d|\bar{\varOmega}_2}{\sqrt{1 - \lambda b(k_+ |d|)}} \left(\frac{n \bar{\varOmega}_d(k_+ |d|) - \omega}{\varOmega_2 + n d \omega_E'}\right)\right) \nonumber\\ &\quad \approx \frac{g(0)}{\dfrac{\bar{\varOmega}_2}{\sqrt{1 - \lambda b (k_+ |d|)}} (n \bar{\varOmega}_d (k_+ |d| ) - \omega) + \varOmega_2 \dfrac{x_+}{d} + n \omega_E' x_+}. \end{align}

For convenience, we next replace all instances of $k_+ | d |$ with $k_+ d$; this is allowed since $\bar {\varOmega }_d$ and $b$ are even functions and the bounds of integration are symmetric. We obtain

(7.37)\begin{equation} I_{p,r} \approx \int \frac{\textrm{d} k_+ \textrm{d} x_+ \textrm{d} x_-}{ 2 {\rm \pi}} \frac{\bar{\varOmega}_2 \, \textrm{e}^{- i k_+ x_-}}{\sqrt{1 - \lambda b (k_+ d)}} \frac{\displaystyle J_0(\rho k_\perp)^2 \phi_{m0,n}\left(x_+{+} \frac{x_-}{2}\right)^\ast \phi_{m0,n}\left(x_+{-} \frac{x_-}{2}\right) }{\displaystyle n \bar{\varOmega}_d (k_+ d) \,{+} \left(\frac{\varOmega_2}{\bar{\varOmega}_2} \frac{x_+}{d} + n \frac{\omega_E'}{\bar{\varOmega}_2} x_+\right)\sqrt{1 \,{-}\, \lambda b(k_+ d)} \,{-}\, \omega}, \end{equation}

where

(7.38)\begin{equation} k_\perp^2 = k_+^2 + \frac{n^2 q_0^2}{r_0^2}. \end{equation}

As with the trapped functional, we separately average over the Bessel functions,

(7.39)\begin{equation} \frac{\displaystyle \int \, \textrm{d}^3v J_0 (k_\perp \rho)^2 f_0}{\displaystyle \int \, \textrm{d}^3v f_0} = \varGamma_0 (k_\perp \rho_\text{th}). \end{equation}

We next carry out the integral over $x_-$ by identifying it as the inverse Fourier transform of the product of two Gaussians, which leads to

(7.40)\begin{equation} I_{p,r} = \int \frac{\textrm{d} x_+ \textrm{d} k_+}{\sqrt{\rm \pi}} \frac{\bar{\varOmega}_2}{\sqrt{1 - \lambda b (k_+ d)}}\frac{\displaystyle \varGamma_0 (k_\perp \rho_{\text{th}}) \sqrt{\frac{\textrm{Im}(w^2)^2}{\textrm{Re}(w^2)} + \textrm{Re}(w^2)} \, \textrm{e}^{-\rho_\ast^2} \, \textrm{e}^{{-}k_\ast^2} \exp\left(\frac{\textrm{Im}(x_0)^2}{\textrm{Re}(w^2)}\right) }{\displaystyle n \bar{\varOmega}_d (k_+ d) + \left(\frac{\varOmega_2}{\bar{\varOmega}_2} \frac{x_+}{d} + n \frac{\omega_E'}{\bar{\varOmega}_2} x_+\right)\sqrt{1 - \lambda b(k_+ d)} - \omega}, \end{equation}

where

(7.41)\begin{gather} \rho_\ast{=} \frac{x_+{+} k_+ \textrm{Im}(w^2) - \textrm{Re}(x_0)}{\sqrt{\textrm{Re}(w^2)}}, \end{gather}

(7.42)\begin{gather}k_\ast{=} \frac{k_+ \textrm{Re}(w^2) + \textrm{Im} (x_0)}{{\sqrt{\textrm{Re}(w^2)}}}. \end{gather}

It is more convenient to numerically integrate this over $k_\ast$ and $\rho _\ast$ to take advantage of the explicit Gaussian structure. Because the Jacobian of this variable transformation is $1$, the change of variables is easily carried out. In addition, we approximate the $\lambda$-dependent terms by averaging over the pitch angle parameter. We also use the extremely passing particle limit, where $\theta \approx \alpha _2$. We then obtain

(7.43)\begin{equation} I_{p,r} \approx \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{\textrm{d} \rho_\ast \, \textrm{d} k_\ast}{\sqrt{\rm \pi}} \frac{\displaystyle \varGamma_0 (k_\perp \rho_{\text{th}}) \sqrt{\frac{\textrm{Im}(w^2)^2}{\textrm{Re}(w^2)} + \textrm{Re}(w^2)} \, \textrm{e}^{-\rho_\ast^2} \, \textrm{e}^{{-}k_\ast^2} \exp\left(\frac{\textrm{Im}(x_0)^2}{\textrm{Re}(w^2)}\right) }{\displaystyle n \omega_{d0} \xi F_p(k_+ d) + \epsilon_\parallel \sqrt{\xi} \frac{x_+}{d} \frac{\sqrt{2T/m}}{q R_0} - \bar{\omega}}, \end{equation}

where

(7.44)\begin{equation} \bar{\omega} = \omega - n \omega_{E0} - n \omega_E' x_+, \end{equation}

and

(7.45)\begin{equation} F_p(k_+ d) = \tfrac{4}{3} (\cos(k_+ d) + (s k_+ d - \alpha_M \sin(k_+ d))\sin(k_+ d)). \end{equation}

We note here that the factor of $4/3$ comes from taking the pitch angle average of $2 - \lambda b$ in the small $\epsilon$ limit. This approximation can be improved by considering higher-order $\epsilon$ terms, although this is not done in the current formulation of QuaLiKiz (QuaLiKiz version 2.8.1).

We now address terms in the numerator of the original integrand that are proportional to $x$; these terms arise from the radial electric field shear. In principle, their inclusion can be treated fully consistently by using the appropriate Fourier transforms as well as the convolution theorem in much the same way as we did before. However, as a crude approximation, we simply map $x \to x_+$ in the numerator as is effectively done in the denominator.

Next, we address the integration over $\lambda$ and $\xi$ in the full passing functional. Once these integrals are calculated, we fold them into the integration over $\rho _\ast$ and $k_\ast$. We wish to compute

(7.46)\begin{align} I_{p,E} &= \sum_{\epsilon_\parallel} \int \, \textrm{d} \lambda\, \textrm{d} \xi (2 {\rm \pi})^2 \frac{n_0 e^2}{T} \frac{R_0 \bar{r}}{\sqrt{\rm \pi}} \frac{\sqrt{\xi} \, \textrm{e}^{-\xi}}{\bar{\varOmega}_2}\nonumber\\ &\quad \times\frac{\displaystyle \left\langle (\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_\ast{-} \bar{\omega} ) \left( 1 + \frac{2 v_\parallel U_\parallel}{v_T^2} + \frac{U_\parallel^2}{v_T^2} \left(\frac{2 v_\parallel^2}{v_T^2} - 1\right) \right) \right\rangle }{\displaystyle n \omega_{d0} \xi F_p(k_+ d) + \epsilon_\parallel \sqrt{\xi} \frac{x_+}{d} \frac{\sqrt{2T/m}}{q R_0} - \bar{\omega}}. \end{align}

Because we averaged out the pitch angle dependence in the denominator of the integrand, the pitch angle integration in the numerator is completely separable and only dependent on the inverse aspect ratio $\epsilon$. This is perhaps the largest single approximation used in the passing part of the dispersion; it is necessary to ensure that the numerical integral is 2-D rather than 3-D. It is of potential interest to study the impact of this approximation. Calculating a more exact (albeit slower) integral to quantify the exact impact this approximation has on the resulting solutions and flux calculations is the subject of future work.

Since only $v_\parallel$ terms in the numerator are dependent on $\lambda$. We also use the fact that

(7.47)\begin{equation} \int_{0}^{({1-\epsilon})/({1+\epsilon})} \frac{\textrm{d} \lambda }{2 \bar{\varOmega}_2} = f_p, \end{equation}

where $f_p = 1 - f_t$ is the flux surface averaged passing particle fraction. We then compute

(7.48)\begin{equation} \int_{0}^{({1-\epsilon})/({1+\epsilon})} \frac{\textrm{d} \lambda}{\bar{\varOmega}_2} \langle v_\parallel^m\rangle = 2 f_p v_T^{m} \epsilon_\parallel^m \xi^{m/2} \lambda_m, \end{equation}

where we define $\lambda _m$ as

(7.49)\begin{equation} \lambda_m = \frac{\displaystyle \int_{0}^{({1-\epsilon})/({1+\epsilon})}\, \textrm{d} \lambda \int_{-{\rm \pi}}^{\rm \pi} \frac{\textrm{d} \theta \left(\sqrt{1 - \lambda b}\right)^m}{\sqrt{1-\lambda b}}}{\displaystyle \int_{0}^{({1-\epsilon})/({1+\epsilon})}\, \textrm{d} \lambda \int_{-{\rm \pi}}^{\rm \pi} \frac{\textrm{d} \theta}{\sqrt{1-\lambda b}}}. \end{equation}

We numerically calculate $\lambda _m$ separately from the rest of the dispersion relation since $\lambda _m$ is only dependent on $\epsilon$. Once again ignoring terms order cubic or higher with the Mach number, we find that

(7.50)\begin{align} &\int_0^{({1-\epsilon})/({1+\epsilon})} \frac{\textrm{d} \lambda}{\bar{\varOmega}_2} \left\langle ( \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_\ast{-} \bar{\omega} ) \left(1 + \frac{2 v_\parallel U_\parallel}{v_T^2} + \frac{U_\parallel^2}{v_T^2} \left( \frac{2 v_\parallel^2}{v_T^2} - 1\right)\right) \right\rangle \nonumber\\ & \quad = 2 f_p n\omega_{d0}( \mathcal{A}_p + \mathcal{B}_p \epsilon_\parallel \xi^{1/2} + \mathcal{C}_p \xi + \mathcal{D}_p \epsilon_\parallel \xi^{3/2} + \mathcal{E}_p \xi^2), \end{align}

where we define the terms

(7.51)\begin{gather} \mathcal{A}_p = \left(1 - \frac{U_\parallel^2}{v_T^2}\right) \left(A_n - \frac{3}{2} A_T - z^2 F_p \right) - \frac{U_\parallel}{v_T} \left(2 A_U - \frac{U_\parallel}{v_T} A_T\right), \end{gather}

(7.52)\begin{gather}\mathcal{B}_p = \left(A_U \left(2 - 6 \frac{U_\parallel^2}{v_T^2} \right) + U_\parallel (2 A_n - 5 A_T - 2 z^2 F_p )\right) \lambda_1, \end{gather}

(7.53)\begin{gather}\mathcal{C}_p = \left(1 - \frac{U_\parallel^2}{v_T^2}\right) A_T + 4 A_U \frac{U_\parallel}{v_T} \lambda_2 + \frac{U_\parallel^2}{v_T^2} \lambda_2 (2 A_n - 7 A_T - 2 z^2 F_p ), \end{gather}

(7.54)\begin{gather}\mathcal{D}_p = 2 A_T \frac{U_\parallel}{v_T} \lambda_1 + 4 A_U \frac{U_\parallel^2}{v_T^2} \lambda_3, \end{gather}

(7.55)\begin{gather}\mathcal{E}_p = 2 A_T \frac{U_\parallel^2}{v_T^2} \lambda_2, \end{gather}

and where

(7.56)\begin{equation} z^2 = \frac{\bar{\omega}}{n \omega_{d0} F_p(k_+ d)}. \end{equation}

Thus, the integral simplifies to

(7.57)\begin{align} I_{p,E} = \sum_{\epsilon_\parallel} 2 f_p \int_0^\infty \, \textrm{d}\xi (2 {\rm \pi})^2 \frac{n_0 e^2}{T} \frac{R_0 r_0}{\sqrt{\rm \pi}} \sqrt{\xi} \frac{\mathcal{A}_p + \mathcal{B}_p \epsilon_\parallel \xi^{1/2} + \mathcal{C}_p \xi + \mathcal{D}_p \epsilon_\parallel \xi^{3/2} + \mathcal{E}_p \xi^2}{\displaystyle F_p(k_+ d) \left( \xi + \epsilon_\parallel \sqrt{\xi} \frac{x_+}{d} \frac{v_T}{q R_0 F_p(k_+ d)} - z^2\right)}. \end{align}

We then perform a change in variables to $u = \sqrt {\xi }$ and note that

(7.58)\begin{equation} \sum_{\epsilon_\parallel} \int_0^{\infty}\, \textrm{d} \xi \sqrt{\xi} g(\epsilon_\parallel \sqrt{\xi}) = \sum_{\epsilon_\parallel} \int_0^{\infty}\, \textrm{d} u 2 u^2 g(\epsilon_\parallel u) = \int_{-\infty}^{\infty} \, \textrm{d} u 2 u^2 g(u). \end{equation}

The integral then becomes

(7.59)\begin{equation} I_{p,E} = 4 f_p \int_{-\infty}^\infty \, \textrm{d} u (2 {\rm \pi})^2 \frac{n_0 e^2}{T} \frac{R_0 r_0}{\sqrt{\rm \pi}} u^2 \frac{ \mathcal{A}_p + \mathcal{B}_p u + \mathcal{C}_p u^2 + \mathcal{D}_p u^3 + \mathcal{E}_p u^4}{F_p(k_+ d) \left( u^2 + u \dfrac{x_+}{d} \dfrac{v_T}{q R_0 F_p(k_+ d)} - z^2\right)}. \end{equation}

To simplify this integral further, we rewrite the denominator as

(7.60)\begin{equation} u^2 + u \frac{x_+}{d} \frac{v_T}{q R_0 F_p(k_+ d)} - z^2 = (u - z_+) (u - z_-), \end{equation}

where

(7.61)\begin{equation} z_{{\pm}} ={-}\frac{1}{2} \frac{x_+}{d} \frac{v_T}{q R_0 F_p(k_+ d)} \pm \sqrt{\left(\frac{1}{2} \frac{x_+}{d} \frac{v_T}{q R_0 F_p(k_+ d)}\right)^2 + z^2}. \end{equation}

This allows us to simplify the integral using the plasma dispersion functions defined in Appendix A, which allows us to obtain

(7.62)\begin{equation} I_{p,E} = 4 f_p (2 {\rm \pi})^2 \frac{n_0 e^2}{T} \frac{R_0 r_0}{F_p(k_+ d)} ( \mathcal{A}_p G_2 + \mathcal{B}_p G_3 + \mathcal{C}_p G_4 + \mathcal{D}_p G_5 + \mathcal{E}_p G_6 ), \end{equation}

where the associated Fried and Conte integrals $G_n = G_n(z_+, z_-)$ are evaluated at $z_+$ and $z_-$. Thus, the passing functional simplifies to

(7.63)\begin{align} \mathcal{L}_\text{passing} &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{\textrm{d} \rho_\ast \, \textrm{d} k_\ast }{\sqrt{\rm \pi}} (2 {\rm \pi})^3 \frac{e^2 n_0}{T} r_0 R_0 f_p \langle \mathcal{I}_p\rangle _{\xi, \lambda} \varGamma_0 (k_\perp \rho_{\text{th}}) \nonumber\\ &\quad \times \sqrt{\frac{\textrm{Im}(w^2)^2}{\textrm{Re}(w^2)} + \textrm{Re}(w^2)} \, \textrm{e}^{-\rho_\ast^2} \, \textrm{e}^{{-}k_\ast^2} \exp\left(\frac{\textrm{Im}(x_0)^2}{\textrm{Re}(w^2)}\right), \end{align}

where

(7.64)\begin{equation} \langle \mathcal{I}_p\rangle _{\xi, \lambda} = \frac{2}{{\rm \pi} F_p (k_+)}(\mathcal{A}_p G_2 + \mathcal{B}_p G_3 + \mathcal{C}_p G_4 + \mathcal{D}_p G_5 + \mathcal{E}_p G_6 ). \end{equation}

We have now reduced all parts of the dispersion relation to a numerically tractable form. The adiabatic piece can be calculated analytically, whereas the trapped and passing functionals require one- and two-dimensional integrals, respectively. With the dispersion relation in hand, we can proceed to apply quasilinear theory.

8. Quasilinear approximation

The core principle of quasilinear theory is to consider the slow time variation of the total distribution function $f$ and the resultant fluxes that attempt to drive the distribution function back to equilibrium. The validity of the quasilinear approximation depends on the decorrelation time of the potential being shorter than the eddy turnover time. The ratio of these two quantities is known as the Kubo number (Kubo Reference Kubo1963; Krommes Reference Krommes2002). The single particle analogue to this is that the individual particle must not be trapped in the field; this allows the dynamics to be characterized as a random walk process, which leads to a justification for the quasilinear approach. These characteristic times have been calculated and compared for both ETG and ITG–TEM turbulence in various parameter regimes (Lin et al. Reference Lin, Rice, Wukitch, Greenwald, Hubbard, Ince-Cushman, Lin, Porkolab, Reinke and Tsujii2008; Casati et al. Reference Casati, Bourdelle, Garbet, Imbeaux, Candy, Clairet, Dif-Pradalier, Falchetto, Gerbaud, Grandgirard, Gürcan, Hennequin, Kinsey, Ottaviani, Sabot, Sarazin, Vermare and Waltz2009; Citrin et al. Reference Citrin, Bourdelle, Cottier, Escande, Gürcan, Hatch, Hogeweij, Jenko and Pueschel2012). For all cases reported in the aforementioned references, the Kubo number is less than unity, even for those with normalized gradients of ${\sim }20$. Well-developed turbulence for these tokamak plasma parameters thus manifest random walk processes, and thus driven anomalous transport in the core is successfully modelled by quasilinear methods.

Moreover, it has been found that quasilinear models are successful in reproducing experimental results such as temperature profiles within $15\,\%$ root-mean-square error (Kinsey, Staebler & Waltz Reference Kinsey, Staebler and Waltz2008). We acknowledge, however, that successful quasilinear modelling of the total fluxes is not strictly indicative of the validity of the quasilinear approximation. For instance, the saturation rule for a quasilinear model can in principle correct for errors in the unsaturated quasilinear flux calculations. Thus, validation of quasilinear flux ratios before saturation is an important piece in confirming the use of the quasilinear approximation. We note that investigating the full range of validity for the quasilinear approximation is an active area of research, and that cases whereby the quasilinear approximation breaks down have been found (Laval & Pesme Reference Laval and Pesme1983; Besse et al. Reference Besse, Elskens, Escande and Bertrand2011). For QuaLiKiz, we restrict ourselves to the aforementioned parameter regimes whereby the quasilinear approach has been thoroughly verified.

To proceed, we first recall the Vlasov equation for a given species (again omitting the species label):

(8.1)\begin{equation} \frac{\partial f}{\partial t} + \dot{\boldsymbol{\alpha}} \boldsymbol{\cdot} \frac{\partial f}{\partial \boldsymbol{\alpha}} + \dot{\boldsymbol{J}} \boldsymbol{\cdot} \frac{\partial f}{\partial \boldsymbol{J}} = 0. \end{equation}

When we obtained the dispersion relation, we considered the linear response and neglected terms that are quadratic in the fluctuations. Moreover, we also assumed $f_0$ was time-independent. To proceed with the quasilinear approximation, we now suppose that $f_0$ varies slowly in time on a time scale longer than that of the linear modes. We may then perform a time average over the Vlasov equation such that $\langle f\rangle _t = f_0$ and the linear response averages to zero. We define the time average as

(8.2)\begin{equation} \langle g(t)\rangle _t = \frac{1}{T} \int_{{-}T/2}^{T/2} g(t + t') \, \textrm{d} t', \end{equation}

where $T$ is the time scale associated with the linear modes. The time-averaged Vlasov equation then reads

(8.3)\begin{equation} \left\langle \frac{\partial f}{\partial t} + \dot{\boldsymbol{\alpha}} \boldsymbol{\cdot} \frac{\partial f}{\partial \boldsymbol{\alpha}} + \dot{\boldsymbol{J}} \boldsymbol{\cdot} \frac{\partial f}{\partial \boldsymbol{J}}\right\rangle _t \approx \frac{\partial f_0}{\partial t} + \langle \{\textrm{Re}(\delta f), \textrm{Re}(e \phi) \} \rangle _t = 0. \end{equation}

Here, we take the real part of $\delta f$ or $\phi$ to obtain the physical quantity in accordance with our convention. To proceed, we rewrite the Poisson bracket as

(8.4)\begin{equation} \{\textrm{Re}(\delta f), \textrm{Re}(e \phi) \} = \frac{\partial}{\partial \boldsymbol{\alpha}} \boldsymbol{\cdot} \left(\textrm{Re}(\delta f) \frac{\partial \textrm{Re}(e \phi)}{\partial \boldsymbol{J}}\right) - \frac{\partial}{\partial \boldsymbol{J}} \boldsymbol{\cdot} \left(\textrm{Re}(\delta f) \frac{\partial \textrm{Re}(e \phi)}{\partial \boldsymbol{\alpha}}\right). \end{equation}

The time average can be simplified by noting that for any two general vectors $\boldsymbol {A}$ and $\boldsymbol {B}$ we have

(8.5)\begin{equation} \langle \textrm{Re}(\boldsymbol{A} \, \textrm{e}^{- {\rm i} \omega t}) \boldsymbol{\cdot} \textrm{Re}(\boldsymbol{B} \, \textrm{e}^{- {\rm i} \omega t})\rangle _t = \tfrac{1}{2} \textrm{Re}(\boldsymbol{A} \boldsymbol{\cdot} \boldsymbol{B}^\ast). \end{equation}

Due to the Fourier structure of $\delta f$ and $\phi$, we also note that

(8.6)\begin{equation} \frac{\partial}{\partial \boldsymbol{\alpha}} \langle \textrm{Re}(\delta f) \textrm{Re}(e \phi) \rangle _t = \boldsymbol{0}. \end{equation}

Essentially, the $\boldsymbol {\alpha }$ dependence disappears after performing the time averaging. Moreover, taking the real part of $\delta f$ and $\phi$ commutes with taking derivatives of real variables. We therefore obtain

(8.7)\begin{equation} \frac{\partial f_0}{\partial t} + \frac{\partial}{\partial \boldsymbol{J}} \boldsymbol{\cdot} \boldsymbol{\varGamma}_Q = 0, \end{equation}

where we define the quasilinear flux $\boldsymbol {\varGamma }_Q$ as

(8.8)\begin{align} \boldsymbol{\varGamma}_Q = \frac{1}{2} \textrm{Re}\left(\sum_{\boldsymbol{n}} {\rm i} \boldsymbol{n} f_{\boldsymbol{n}} e \phi_{\boldsymbol{n}}^\ast \right) ={-}\frac{1}{2} \textrm{Im} \left(\sum_{\boldsymbol{n}} \boldsymbol{n} \frac{e^2 |\phi_{\boldsymbol{n}}|^2}{T} f_0 \left(1 - \frac{\omega - \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_\ast{-} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\omega}_E}{\omega - \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\varOmega}}\right) \right). \end{align}

Here, $f_{\boldsymbol {n}}$ and $\phi _{\boldsymbol {n}}$ are related via the dispersion relation in the linearized problem. Thus, the quasilinear flux is computed by substituting in the solution of the dispersion relation including the found eigenvalues $\omega$, again only considering unstable modes. Modes that lack unstable solutions do not contribute to the quasilinear flux.

We are now in a position to calculate the flux surface averaged particle, toroidal angular momentum and energy fluxes by averaging the Vlasov equation over velocity and space. This is analogous to calculating the fluid equations by taking moments of the Vlasov equation. The radial fluxes can be calculated via a change in variables from $J_3$ to $r$. We find that

(8.9)\begin{gather} \frac{\partial \langle n\rangle }{\partial t} + \frac{\textrm{d} \varGamma}{\textrm{d} r} = 0, \end{gather}

(8.10)\begin{gather}\frac{\partial \langle m n R U_\parallel\rangle }{\partial t} + \frac{\textrm{d} \varPi}{\textrm{d} r} = 0, \end{gather}

(8.11)\begin{gather}\frac{3}{2} \frac{\partial \langle p \rangle }{\partial t} + \frac{\textrm{d} Q}{\textrm{d} r} = 0, \end{gather}

where $\varGamma$, $\varPi$ and $Q_E$ are the particle, toroidal momentum and energy fluxes defined as

(8.12)\begin{gather} \varGamma = \frac{1}{4 {\rm \pi}^2 d} \int \, \textrm{d}^3v \, \textrm{d}^3r \frac{1}{2} \textrm{Im}\left(\sum_{\boldsymbol{n}} \frac{n q}{r B} f_{\boldsymbol{n}} \phi_{\boldsymbol{n}}^\ast \right), \end{gather}

(8.13)\begin{gather}\varPi = \frac{1}{4 {\rm \pi}^2 d} \int \, \textrm{d}^3v \, \textrm{d}^3r \frac{m R v_\parallel}{2} \textrm{Im}\left(\sum_{\boldsymbol{n}} \frac{n q}{r B} f_{\boldsymbol{n}} \phi_{\boldsymbol{n}}^\ast \right), \end{gather}

(8.14)\begin{gather}Q_E = \frac{1}{4 {\rm \pi}^2 d} \int \, \textrm{d}^3v \, \textrm{d}^3r \frac{m(v^2 - U_\parallel^2)}{4} \textrm{Im}\left(\sum_{\boldsymbol{n}} \frac{n q}{r B} f_{\boldsymbol{n}} \phi_{\boldsymbol{n}}^\ast \right). \end{gather}

Here, we can see that the integrations to calculate the particle, toroidal momentum and energy fluxes are of the same form to solve the dispersion relation. The particle flux calculation is identical. Meanwhile, we must take into account an extra factor of $v_\parallel$ and $v^2$ for the angular momentum flux and energy flux integrations, respectively. These changes can be easily accommodated for without affecting the fundamental approach. For instance, the inclusion of $v^2$ simply changes the associated Fried and Conte integral. The physical significance of these fluxes can be further solidified by examining the perturbed $E$-cross-$B$ velocity. We find that

(8.15)\begin{equation} \delta \boldsymbol{v}_{E\times B} \boldsymbol{\cdot} \hat{\boldsymbol{r}} = \hat{\boldsymbol{r}} \boldsymbol{\cdot} -\frac{\boldsymbol{\nabla} \delta \phi \times \boldsymbol{B}}{B^2} \approx \sum_{\boldsymbol{n}} \frac{{\rm i} k_\theta \phi_{ \boldsymbol{n}}}{B} = \sum_{\boldsymbol{n}} \frac{{\rm i} n q}{r B} \phi_{\boldsymbol{n}}, \end{equation}

where we have again used the convention that $k_\theta \to (i/r) \partial _\theta$. We then find that

(8.16)\begin{equation} \langle \textrm{Re}(\delta \boldsymbol{v}_{E\times B} \boldsymbol{\cdot} \hat{\boldsymbol{r}}) \textrm{Re}(\delta f) \rangle _t = \frac{1}{2} \textrm{Re}\left( \sum_{\boldsymbol{n}} f_{\boldsymbol{n}} \left(\frac{{\rm i} n q}{r B} \phi_{\boldsymbol{n}} \right)^\ast \right) = \frac{1}{2} \textrm{Im}\left( \sum_{\boldsymbol{n}} f_{\boldsymbol{n}} \frac{n q}{rB} \phi_{\boldsymbol{n}}^\ast \right). \end{equation}

This lets us write the fluxes as

(8.17)\begin{gather} \varGamma = \langle \delta (n) \delta \boldsymbol{v}_{E\times B} \boldsymbol{\cdot} \hat{\boldsymbol{r}} \rangle _{t, r}, \end{gather}

(8.18)\begin{gather}\varPi = \langle \delta (m n R U_\parallel) \delta \boldsymbol{v}_{E\times B} \boldsymbol{\cdot} \hat{\boldsymbol{r}} \rangle _{t, r}, \end{gather}

(8.19)\begin{gather}Q_E = \langle \delta (P) \delta \boldsymbol{v}_{E\times B} \boldsymbol{\cdot} \hat{\boldsymbol{r}} \rangle _{t, r}. \end{gather}

Therefore, the particle, angular momentum and energy fluxes are simply related to moments of the perturbed distribution function integrated against the perturbed $E$-cross-$B$ velocity, where $\langle \dots \rangle _{t,r}$ denotes a time and spatial average. We also define

(8.20)\begin{gather} \delta n = \int \, \textrm{d}^3v \delta f, \end{gather}

(8.21)\begin{gather}\delta (n m R U_\parallel ) = \int \, \textrm{d}^3v m v_\parallel \delta f, \end{gather}

(8.22)\begin{gather}\delta P = \int \, \textrm{d}^3v \tfrac{1}{2} m (v^2 - U_\parallel^2)\delta f, \end{gather}

where we calculate the particle, angular momentum and energy fluxes for every species. We note that the toroidal angular momentum flux is only non-zero in the presence of rotations. The energy flux calculation can be approximated by noting in the small-Mach-number limit that

(8.23)\begin{equation} \delta P = \delta (n T) ={-} \tfrac{1}{2} m U_\parallel^2 \delta n + \int \, \textrm{d}^3v \tfrac{1}{2} m v^2 \approx \int \, \textrm{d}^3v \tfrac{1}{2} m v^2. \end{equation}

We also note that often we are concerned with the heat flux $Q$ relative to the convective energy flux $\frac {3}{2} T \varGamma$ (Horton Reference Horton1984). The heat flux is simply

(8.24)\begin{equation} Q = Q_E - \frac{3T}{2} \varGamma. \end{equation}

It is important to note that while we may obtain quasilinear flux ratios from the above procedure, we cannot with linear physics alone obtain the physical fluxes. Throughout the derivation, we have kept the amplitude of the fluctuating potential $\delta \phi$ arbitrary. The amplitude $\phi _0$ can only be obtained through the use of nonlinear physics by saturating the amplitude. Thus, the complete calculation of these fluxes must be obtained via a saturation rule obtained from a nonlinear computational code, in this case from the Gyrokinetic Electromagnetic Numerical Experiment (GENE) (Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000). This saturation rule is the topic of § 9.

9. Saturation rule

To formulate a saturation rule, we introduce the well-known mixing length estimate with an effective diffusivity $D$:

(9.1)\begin{equation} D = \left. \frac{\gamma_n}{\langle k_\perp^2\rangle } \right|_{\max}, \end{equation}

where we compute the value of $\gamma _n$ such that the quantity $\gamma _n / \langle k_\perp ^2\rangle$ is at its maximum over the linear spectrum for a given mode. Meanwhile, we average $k_\perp ^2$ over the electrostatic mode. We enforce this mixing length estimate for our various flux calculations by approximating the underlying process as a random walk (Bourdelle et al. Reference Bourdelle, Garbet, Imbeaux, Casati, Dubuit, Guirlet and Parisot2007). For instance, we mandate that the particle flux for a given species must be

(9.2)\begin{equation} \varGamma_s = \sum_{n} C_{\text{NL}} \frac{S_n}{R_0 n_{0s}} \left. \frac{\gamma_n}{\langle k_\perp^2\rangle } \right|_{\max} \frac{k_\theta}{k_{\theta, \text{max}}} L_{s, n, 0}, \end{equation}

where $C_{\text {NL}}$ is a dimensionless constant from nonlinear physics, the form factor $S_n$ is a mode-dependent form factor, $k_{\theta , \text {max}}$ corresponds to the mode that maximizes $\gamma _n / k_\perp ^2$ and $L_{s,n,0}$ is the dimensionless integral that actually computes the flux terms. The above expression is only valid when there is only one mode present in the linear spectrum. We can generalize the expression to account for the existence of multiple types of linear modes by introducing another form factor $S_{n'}$ into the expression and summing over both $n$ and $n'$, while we compute the maximum $\gamma _n / k_\perp ^2$ for a given $n'$.

We model $C_{\text {NL}}$ with the use of nonlinear gyrokinetic simulations. We distinguish between ITG scales, which we define as $k_\theta \rho _s < 2$, and ETG scales, which we define as $k_\theta \rho _s > 2$. Here, $\rho _s$ is the gyroradius of the main ion species such that $\rho _s = \sqrt {T_s / m_s} / \varOmega _{1,s}$ (note that lack of $\sqrt {2}$). The ITG scales are tuned to the GA-Standard nonlinear ion heat flux computed by GENE, whereas the ETG scales are tuned to a single-scale nonlinear GENE simulation based on JET parameters (Citrin Reference Citrin2017). These parameters are current as of QuaLiKiz version 2.8.1 and are subject to future change depending on updates to the nonlinear physics. The result is

(9.3)\begin{equation} C_{\text{NL}} = \left\{ \begin{array}{@{}ll} 271 / s_{\text{fac}} & \text{if } k_\theta \rho_s < 2 \text{ (ITG)},\\ 122 f_{\text{multi-scale}}/ s_{\text{fac}} & \text{if } k_\theta \rho_s > 2 \text{ (ETG)}. \end{array}\right. \end{equation}

Here, we have also introduced an ad hoc factor $s_{\text {fac}}$ for the case of low magnetic shear (Citrin et al. Reference Citrin, Bourdelle, Cottier, Escande, Gürcan, Hatch, Hogeweij, Jenko and Pueschel2012),

(9.4)\begin{equation} s_{\text{fac}} = \left\{ \begin{array}{@{}ll} 2.5 (1 - |s|) & \text{if } |s| < 0.6,\\ 1 & \text{if } |s| > 0.6, \end{array}\right. \end{equation}

as well as a multi-scale rule determined from the maximum of the respective spectra,

(9.5)\begin{equation} f_{\text{multi-scale}} = \frac{1}{\displaystyle 1 + \exp\left(-\frac{1}{5} \left( \frac{\gamma_{\text{ETG}, \max}}{\gamma_{\text{ITG}, \max}} - \sqrt{\frac{m_i}{m_e}}\right)\right)}, \end{equation}

where $m_e$ and $m_i$ are the masses of the electron and main ion, respectively. Here, the sigmoid guarantees a smooth transition from a strongly driven ion-scale mode regime and a strongly driven electron-scale mode regime, since it has been observed that ETG turbulence is suppressed when the ion-scale instability dominates (Maeyama et al. Reference Maeyama, Idomura, Watanabe, Nakata, Yagi, Miyato, Ishizawa and Nunami2015; Howard et al. Reference Howard, Holland, White, Greenwald and Candy2016).

Lastly, we provide an explicit expression for $k_\perp ^2$. In the ITG regime, we need to take into account contributions to $k_r^2$ that arise from the magnetic shear, the mode structure of the electrostatic perturbation and nonlinear effects. Meanwhile, in the ETG regime we assume full isotropization of the mode such that $k_r^2 = k_\theta ^2$. The result is

(9.6)\begin{equation} \langle k_\perp^2\rangle = \left\{ \begin{array}{@{}ll} k_\theta^2 + (k_{r-\text{NL}} + k_{r-\text{shear}})^2 & \text{if } k_\theta \rho_s < 2 \text{ (ITG)},\\ 2 k_\theta^2 & \text{if } k_\theta \rho_s > 2 \text{ (ETG)}. \end{array} \right. \end{equation}

The shear contribution can be calculated analytically as

(9.7)\begin{equation} k_{r-\text{shear}} = k_\theta |s| \sqrt{\langle \theta^2\rangle } = \frac{k_\theta s d}{\sqrt{2} \textrm{Re}(w^2)} \sqrt{\textrm{Re}(w^2) + 2 \textrm{Im}(x_0)^2}, \end{equation}

where we use

(9.8)\begin{equation} \langle \theta^2\rangle = \frac{\displaystyle \int_{-\infty}^{\infty} \theta^2 |\hat{\phi}(\theta)|^2 \, \textrm{d} \theta }{\displaystyle \int_{-\infty}^{\infty} |\hat{\phi}(\theta)|^2 \, \textrm{d} \theta} = \frac{d^2}{2 \textrm{Re}(w^2)} + \left(d \frac{\textrm{Im}(x_0)}{\textrm{Re}(w^2)}\right)^2. \end{equation}

Meanwhile, the nonlinear contribution has been tuned (Citrin et al. Reference Citrin, Bourdelle, Cottier, Escande, Gürcan, Hatch, Hogeweij, Jenko and Pueschel2012) such that

(9.9)\begin{equation} k_{r-\text{NL}} \rho_s = 0.4 \, \textrm{e}^{{-}2 |s|} q^{{-}0.5} + 1.5 \max\{k_\theta \rho_s - 0.2, 0\}. \end{equation}

Having now fully derived analytic expressions for the dispersion relation and quasilinear fluxes, we now discuss the numerical implementation of QuaLiKiz.

10. Numerical implementation

Recall that the dispersion relation is written as

(10.1)\begin{equation} \sum_s \mathcal{L}_{0, s} - \mathcal{L}_{\text{passing},s} - \mathcal{L}_{\text{trapped},s} = 0. \end{equation}

The trapped and passing functionals discussed in §§ 6 and 7 are both functions of the complex frequency $\omega$. Solving the dispersion relation is therefore a matter of finding the zeros of the complex analytic function $D(\omega )$, where

(10.2)\begin{equation} D(\omega) = \sum_s \mathcal{L}_{0, s} - \mathcal{L}_{\text{passing}, s} - \mathcal{L}_{\text{trapped}, s}. \end{equation}

To solve this, we use the Davies method, a numerical technique developed by Davies (Reference Davies1986) to find the zeros of an analytic function within the complex plane. The strategy takes advantage of the argument principle in complex analysis, which states that given a meromorphic function $f(z)$ that

(10.3)\begin{equation} \frac{1}{2 {\rm \pi}i} \oint_C \frac{f'(z)}{f(z)} \, \textrm{d} z = N - P, \end{equation}

where $N$ and $P$ are respectively the number of zeros and poles of $f(z)$ contained within the simple counter-clockwise contour $C$. Here, zero multiplicity and pole order are taken into account. For our purposes, we assume that $f(z)$ has no poles, which leads to

(10.4)\begin{equation} \frac{1}{2 {\rm \pi}i} \oint_C \frac{f'(z)}{f(z)} \, \textrm{d} z = N. \end{equation}

The key of the method is to recognize from Cauchy's residue theorem that, for integer $n$ such that $1 \le n \le N$, we can calculate the integral $S_n$ such that

(10.5)\begin{equation} S_n = \frac{1}{2 {\rm \pi}i} \oint_C z^n \frac{f'(z)}{f(z)} \, \textrm{d} z = \sum_{j=1}^{N} z_{0j}^n, \end{equation}

where $z_{0j}$ is the $j$th root of $f(z)$ (counting repeated roots as separate). We then construct the polynomial

(10.6)\begin{equation} P_N(z) = \prod_{j=1}^{N} (z - z_{0j}) = \sum_{j=0}^N A_j z^{N-j}, \end{equation}

where the coefficients $A_j$ can be computed from the relations

(10.7)\begin{gather} A_0 = 1, \end{gather}

(10.8)\begin{gather}S_1 + A_1 = 0, \end{gather}

(10.9)\begin{gather}S_2 + A_1 S_1 + 2 A_2 = 0, \end{gather}

(10.10)\begin{gather}S_{n} + A_1 S_{n-1} + A_2 S_{n-2} + \cdots + n A_n = 0, \quad n = 1,2,\ldots,N. \end{gather}

Excluding the trivial $A_0$ term, this is a linear system of $N$ equations. After solving this system, we can then construct the polynomial $P_N$ which has zeros that are precisely the solutions of the dispersion relation. We then extract a zero from the polynomial $P_N$ using a Newton solver and then define a new set of coefficients such that

(10.11)\begin{equation} S^{(1)}_n = S_n - z_{01}^n, \quad n = 0, 1, \ldots, N-1 , \end{equation}

where $z_{01}$ is the first zero found. With this new set of coefficients, we may then construct a new polynomial $P_{N-1}(z)$ and extract another zero. This process is repeated until all zeros are found. If the contour of integration is a unit circle, then a clever integration by parts results in

(10.12)\begin{equation} S_n ={-}\frac{n}{2 {\rm \pi}{\rm i}} \int_{0}^{2{\rm \pi}} \, \textrm{d} \theta \, \textrm{e}^{{\rm i} n \theta} \ln (\textrm{e}^{- {\rm i} N \theta} f(\textrm{e}^{{\rm i} \theta})), \quad n > 0. \end{equation}

The inclusion of $\textrm {e}^{- \textrm {i} N \theta }$ inside the logarithm is to handle the branch cut of the logarithm, and can be obtained by using $z^{-N} f(z)$ instead of $f(z)$ in the preceding formulae; such a substitution does not affect the value of $S_n$ for $n > 0$. For the $n= 0$ case, we simply compute the total change of the argument of $f(\textrm {e}^{\textrm {i} \theta })$ for $0 \le \theta \le 2 {\rm \pi}$ while keeping track of any jumps in the argument that would indicate a full winding. Thus, $S_n$ can be computed via standard quadrature methods for one-dimensional integration.

To apply this to the dispersion relation, we make use of a bijective mapping $\omega = \omega (z)$ (to be determined momentarily). This will allow us to retain the simplifications that come from integrating around a unit circle. The first step is to define $f(z)$ such that

(10.13)\begin{equation} f(z) = D(\omega(z)). \end{equation}

Then, we compute $S_n$ via numerical quadrature, which leads to roots $z_{0n}$ such that

(10.14)\begin{equation} f(z_{0n}) = D(\omega(z_{0n})) = 0. \end{equation}

Because the mapping is bijective, we may then simply apply the mapping onto the roots $z_{n0}$ to obtain

(10.15)\begin{equation} \omega_{0n} = \omega(z_{0n}), \end{equation}

where $\omega _{0n}$ are all the roots within the contour $C$ in the complex $\omega$-plane such that

(10.16)\begin{equation} D(\omega_{0n}) = 0. \end{equation}

The only task remaining is to define a suitable bijective mapping $\omega (z)$. Because QuaLiKiz only considers unstable modes, we demand that $\textrm {Im}(\omega ) > 0$ along the entirety of the contour in the $z$-plane. We first define the bijective mapping $(u, v) \to (x, y)$ as

(10.17)\begin{gather} x(u, v) = \frac{\text{sgn} (uv)}{v \sqrt{2}} \sqrt{u^2 + v^2 - \sqrt{(u^2 + v^2) (u^2 + v^2 - 4 u^2 v^2 )} }, \end{gather}

(10.18)\begin{gather}y(u, v) = \frac{\text{sgn} (uv) }{u \sqrt{2}} \sqrt{u^2 + v^2 - \sqrt{(u^2 + v^2) (u^2 + v^2 - 4 u^2 v^2 )} }. \end{gather}

The inverse mapping is given by

(10.19)\begin{gather} u(x, y) = \frac{x \sqrt{x^2 + y^2 - x^2 y^2}}{\sqrt{x^2 + y^2}}, \end{gather}

(10.20)\begin{gather}v(x, y) = \frac{y \sqrt{x^2 + y^2 - x^2 y^2}}{\sqrt{x^2 + y^2}}. \end{gather}

Since this mapping does not satisfy the Cauchy–Riemann equations, it is merely bijective, not conformal. This is known as a squircle mapping since it appears to be a square with rounded edges, and this specific kind was first formulated by Guasti (Reference Guasti1992). Denoting $\omega = x' + {\rm i} y'$ and $z = u + {\rm i} v$, we modify this mapping such that

(10.21)\begin{gather} x' = R_{x} + \frac{r_x}{a} x(a u, a v), \end{gather}

(10.22)\begin{gather}y' = R_{y} +\frac{1}{a} (R_{y} - \epsilon_y ) y (a u, a v). \end{gather}

With $C$ being the unit circle in the complex $z$-plane, let $C'$ be the mapped curve in the complex $\omega$-plane. Here, $(R_x, R_y)$ determines the approximate centre $C'$, $r_x$ and $a$ are scaling factors chosen to manipulate $C'$ into a rectangular shape, and $\epsilon _y$ is chosen to guarantee that $C'$ lies above the real axis. While the mapping is not conformal, it is sufficient for our method, since not only is it bijective but points interior to $C$ are mapped to the interior of $C'$. Thus, if we make the interior area of $C'$ sufficiently large and place it slightly above the real axis in the complex $\omega$-plane, then we will determine all eigenmodes of interest to us. After the solution frequencies are found, they are then refined using a standard Newton root-finding method.

While the contour integral and the Newton root-finding are done when QuaLiKiz is used on its own, when coupled to an integrated modelling suite a slight modification is made to the algorithm. We assume that the quasilinear transport changes slowly compared with the time scale of evolution of the plasma equilibrium. A typical transport solver iterates on a time step that is of the order of $\lesssim 10^{-2} \tau _E$, where $\tau _E$ is the energy confinement time. To speed up the code, QuaLiKiz will often only use the previous solution as an initial guess for the Newton solver rather than perform the full contour integral. Since codes like QuaLiKiz are often the bottleneck for the whole integrated modelling suite, such a speedup is necessary to make the simulation tractably feasible. In practice, QuaLiKiz will only perform the full contour integral once every ${\sim }10$ iterations.

Lastly, we discuss the numerical integration scheme currently in use by QuaLiKiz to calculate the trapped and passing functionals, which require 2-D integrations. Although QuaLiKiz used to rely on integration routines provided by the Numerical Algorithms Group (NAG), it now uses open source routines based on the Genz and Malik algorithm, dubbed ‘hcubature’. This algorithm was originally developed by Genz & Malik (Reference Genz and Malik1980); the current implementation is based on the C++ implementation (Johnson Reference Johnson2017). The version of the algorithm in QuaLiKiz has been ported to Fortran and is slightly modified as a result.

The goal of hcubature is to estimate

(10.23)\begin{equation} \boldsymbol{I} = \int_{a_1}^{b_1} \int_{a_2}^{b_2} \cdots \int_{a_n}^{b_n} \boldsymbol{f}(\boldsymbol{x}) \, \text{d}^n x. \end{equation}

Here, $\boldsymbol {I}$ is the estimate of the integral, while $a_i$ and $b_i$ are respectively the individual components of the lower and upper bounds of the integral $\boldsymbol {a}$ and $\boldsymbol {b}$, which are both constant vectors with dimension $n$. Meanwhile, $\boldsymbol {f}$ is a vector function of arbitrary dimension, and $\boldsymbol {x}$ is the argument of the function $\boldsymbol {f}$ and is of dimension $n$. The vectors $\boldsymbol {I}$ and $\boldsymbol {f}$ are of the same dimension. Thus, hcubature approximately integrates a vector integrand over a hyper-rectangle (or equivalently a scaled hypercube, hence the name ‘cubature’). The routine terminates when the global estimate of the absolute or relative error of the integral reach the desired tolerance and also calculates an error vector $\boldsymbol {e}$ with the same dimensionality as the integrand. While calculating the error vector is straightforward, incorporating it into the convergence criterion is non-trivial. In general, to estimate the error, we make a higher-order estimate $\boldsymbol {I}_0$ and a lower-order, less accurate estimate $\boldsymbol {I}_1$ and set the $i$th component of $\boldsymbol {\epsilon }$ to be

(10.24)\begin{equation} \epsilon_i = |(I_0)_i - (I_1)_i|. \end{equation}

For simplicity, we first consider a scalar function that we integrate over a hypercube,

(10.25)\begin{equation} I = \int_{{-}1}^1 \int_{{-}1}^1 \cdots \int_{{-}1}^1 f(x_1, x_2, \ldots, x_n) \, \text{d}^n x. \end{equation}

We estimate the integral using the following rule:

(10.26)\begin{align} I \approx I_0 &= w_1 f(0, 0, \ldots, 0) + w_2 \sum f(\lambda_2, 0, 0,\ldots, 0) + w_3\sum f(\lambda_3, 0, 0,\ldots,0) \nonumber\\ &\quad + w_4 \sum f(\lambda_4, \lambda_4, 0, 0,\ldots,0) + w_5 \sum f(\lambda_5, \lambda_5,\ldots,\lambda_5). \end{align}

Here, we sum over all possible permutations of coordinates while also allowing for sign changes. For example, if $f$ takes three arguments, then

(10.27)\begin{align} &\sum f(\lambda_1, 0,0) = f(\lambda_1,0,0) + f(0,\lambda_1,0) + f(0,0,\lambda_1) \nonumber\\ &\quad + f(-\lambda_1,0, 0) + f(0,-\lambda_1,0) + f(0,0,-\lambda_1). \end{align}

Genz and Malik constrain the parameters $w_i$ and $\lambda _i$ by requiring that the integration be exact for the following functions:

(10.28)\begin{equation} \left.\begin{array}{c@{}} f_1(x_1, x_2,\ldots,x_n) = 1,\\ f_2(x_1, x_2,\ldots,x_n) = x_1^2,\\ f_3(x_1, x_2,\ldots,x_n) = x_1^4,\\ f_4(x_1, x_2,\ldots,x_n) = x_1^6,\\ f_5(x_1, x_2,\ldots,x_n) = x_1^2 x_2^2,\\ f_6(x_1, x_2,\ldots,x_n) = x_1^4 x_2^2,\\ f_7(x_1, x_2,\ldots,x_n) = x_1^2 x_2^2 x_3^2. \end{array}\right\} \end{equation}

In addition, they also fix the parameters $\lambda _3 = \lambda _4$ to be a specific number, and solve the resulting nonlinear system of equations. The result can be found in the paper by Genz & Malik (Reference Genz and Malik1980). To estimate the error, we reuse $\lambda _i$ but calculate different weights $w'_i$ to make a lower-order estimate,

(10.29)\begin{align} I \approx I_1 &= w'_1 f(0, 0,\ldots,0) + w'_2 \sum f(\lambda_2, 0, 0,\ldots,0) \nonumber\\ &\quad + w'_3 \sum f(\lambda_3,0,0,\ldots,0) + w'_4 \sum f(\lambda_4,\lambda_4,0,0,\ldots,0). \end{align}

We calculate the weights with the same method as previously discussed and require the integration be exact for the functions $f_1, f_2, f_3$ and $f_5$. By keeping $\lambda _i$ the same, we can estimate the error with the reuse of function evaluations. The error is taken to be

(10.30)\begin{equation} \epsilon = |I_0 - I_1|. \end{equation}

The estimate procedure easily generalizes to that of a hyper-rectangle by using linear transformations. The calculation of $I_0$, $I_1$ and $\boldsymbol {\epsilon }$ can also be extended to the case of vector integrands by integrating every component simultaneously.

In the case that $n=1$, the above rule no longer applies. Instead, hcubature uses a 15-point Kronrod extension of a 7-point Gaussian quadrature rule. For $n$-point Gaussian quadrature, we estimate the integral

(10.31)\begin{equation} I = \int_{{-}1}^1 f(x) \, \text{d}x \approx \sum_{i=1}^n w'_i f(x_i). \end{equation}

To calculate the weights $w_i$ and the abscissa $x_i$, we require that the integration be exact for all polynomials up to degree $2 n - 1$. It can be shown using Lagrange interpolating polynomials and the theory of orthogonal polynomials that the abscissae $x_i$ correspond to the roots of the Legendre polynomial $P_n$ and that the weights are

(10.32)\begin{equation} w'_i = \frac{2}{(1 - x_i^2)\dfrac{\textrm{d} P_n}{\textrm{d} x}(x_i)^2}, \end{equation}

where the Legendre polynomials are normalized such that $P_n(1) = 1$.

One downside to this method is that the abscissa will in general be completely different for different order rules. Thus, naively comparing an $n$-point rule with an $n+1$-point rule to estimate the error is inefficient. Kronrod discovered that for an $n$-point Gaussian quadrature rule, one could add $n+1$ abscissae to exactly integrate polynomials up to order $3n + 1$, reusing the previous abscissa and computing new weights $w_i$. These new nodes correspond to the zeros of Legendre–Stieltjes polynomials, and their derivation will not be covered here. Thus, the 15-point rule corresponds to

(10.33)\begin{equation} I \approx I_0 = \sum_{i=1}^{15} w_i f(x_i), \end{equation}

the 7-point rule to

(10.34)\begin{equation} I \approx I_1 = \sum_{i=1}^7 w'_i f(x_i), \end{equation}

and the estimated error

(10.35)\begin{equation} \epsilon = |I_0 - I_1|. \end{equation}

Extending this to more general limits of integration simply requires a linear transformation.

Now that we have our integration schemes and error estimation rules for arbitrary $n$, we may proceed to describe the general algorithm.

Algorithm 1:

hcubature

Here, $\boldsymbol {f}$ is the vector integrand, $\boldsymbol {a}$ and $\boldsymbol {b}$ are respectively the lower and upper bounds of the integrand, $\epsilon _a$ and $\epsilon _r$ are respectively the requested absolute and relative error tolerances, maxEval is the maximum number of function evaluations to be allowed by the routine, and norm determines the convergence criterion (in conjunction with the requested error tolerances). The integer eval keeps track of the total number of function evaluations, the vectors $\boldsymbol {I}_0$ and $\boldsymbol {I}_1$ correspond to the integration estimates for a given hyper-rectangle, $\boldsymbol {\epsilon }$ is the error estimate for the hyper-rectangle, and $s$ is the suggested dimension of splitting. As for the output, $\boldsymbol {I}$ is the total integration estimate, $\boldsymbol {e}$ is the total error and ifail is an integer denoting whether any errors occurred while carrying out the procedure or whether the eval reached maxEval before convergence. Convergence is determined using the global error vector $\boldsymbol {e}$.

The algorithm splits the initial hyper-rectangle into pieces and stores them in a binary heap. The heap is sorted according to the largest component of the local error vector, where the root of the heap corresponds to the region with the largest error. Until the integral converges, we pop a hyper-rectangle from the root of the heap, split it into two regions, evaluate both regions accordingly, update the global integration and error estimates, and push both regions into the heap. This guarantees that the split region contributes the greatest to the global error. To determine along which direction to split the hyper-rectangle, we calculate a fourth divided difference using the same evaluation points,

(10.36)\begin{align} D_i &= \sum_j \left|f_j(0,0,\ldots,0,-\lambda_2,0,0,\ldots,0) \vphantom{\frac{\lambda_1^2}{\lambda_2^2}}\right.\nonumber\\ &\quad - 2 f_j(0,0,\ldots,0) + f_j(0,0,\ldots,0,\lambda_2,0,0\ldots,0) \nonumber\\ &\quad -\frac{\lambda_1^2}{\lambda_2^2} [f_j(0,0,\ldots,0,-\lambda_1,0,0,\ldots,0) - 2 f_j(0,0,\ldots,0)\nonumber\\ & \quad +f_j(0,0,\ldots,0,\lambda_1,0,0\ldots,0)] \left.\vphantom{\frac{\lambda_1^2}{\lambda_2^2}}\right|. \end{align}

Here, $i$ corresponds to the dimension at which we evaluate the functions. For example, if $i=2$, then

(10.37)\begin{align} D_2 &= \sum_j \left| \vphantom{\frac{\lambda_1^2}{\lambda_2^2}} f_j(0,-\lambda_2,0,0,\ldots,0) - 2 f_j(0,0,\ldots,0) + f_j(0,\lambda_2,0,0\ldots,0) \right. \nonumber\\ &\quad -\left.\frac{\lambda_1^2}{\lambda_2^2} [f_j(0,-\lambda_1,0,0,\ldots,0) - 2 f_j(0,0,\ldots,0) + f_j(0,\lambda_1,0,0\ldots,0)]\right|. \end{align}

Note that here we take the difference along each component of $\boldsymbol {f}$ and sum the absolute value of each difference. We determine $s$, the dimension along which we split the hyper-rectangle, by calculating the maximum component of $\boldsymbol {D}$. The coordinate corresponding to the maximum of $\boldsymbol {D}$ is the one in which we split the hyper-rectangle in half. For the 1-D case using the Gauss–Kronrod rule, no such calculation is required. We continually split the whole hyper-rectangle into smaller and smaller pieces until convergence is achieved.