# Quasilinear gyrokinetic theory:
A derivation of QuaLiKiz

C. D. Stephens1,2<EMAIL_ADDRESS>X. Garbet3  J. Citrin4  C. Bourdelle3
K. L. van de Plassche4  and F. Jenko2  1 University of California, Los
Angeles, 475 Portola Plaza, Los Angles, California 90095, USA 2 Max Planck
Institute for Plasma Physics, Boltzmannstr. 2, 85748 Garching, Germany 3 CEA,
IRFM, F-13108 Saint Paul-lez-Durance, France 4 DIFFER—Dutch Institute for
Fundamental Energy Research, De Zaale 20, 5612 AJ Eindhoven,
The Netherlands

###### Abstract

In order to predict and analyze turbulent transport in tokamaks, it is
important to model transport that arises from microinstabilities. For this
task, quasilinear codes have been developed that seek to calculate particle,
angular momentum, and heat fluxes both quickly and accurately. In this
tutorial, we present a derivation of one such code known as QuaLiKiz, a
quasilinear gyrokinetic transport code. The goal of this derivation is to
provide a self-contained and complete description of the underlying physics
and mathematics of QuaLiKiz from first principles. This work serves both as a
comprehensive overview of QuaLiKiz specifically as well as an illustration for
deriving quasilinear models in general.

##  1 Introduction

The development of tractable transport models is crucial to further the study
and operation of tokamaks. Accurately characterizing the particle, angular
momentum, and heat transport in the tokamak core requires the understanding of
turbulence driven by microinstabilities, as these instabilities drive much of
the particle, momentum, and heat transport in the core. Integrated modeling
codes seek to predict and simulate tokamak discharges via the inclusion of
various different physics and sources, including from microinstabilities.
Nonlinear simulations of the kinetic equations are the most accurate way to
compute the transport from microinstabilities. For reference, the cost of such
a nonlinear simulation is on the order of $10^{4}$ CPUh to $10^{5}$ CPUh at a
single radial point, while integrated modeling frameworks require thousands
flux calculations for every second of a plasma discharge in a large tokamak
device (Citrin et al., 2017). Multi-scale simulations that take into account
the interplay of instabilities across wide ranges of time scales are even more
expensive (Waltz et al., 2007; Görler & Jenko, 2008; Howard et al., 2016).
Even linear kinetic simulations can prove to be intractable for integrated
modeling if not reduced enough. Thus, it is imperative to develop and refine
kinetic models that are both accurate enough to account for transport from
microinstabilities and fast enough to be coupled to an integrated modeling
framework.

QuaLiKiz is a quasilinear gyrokinetic transport model originally based on the
linear eigenvalue code Kinezero. Pieces of the derivation have been published
throughout the years including in Bourdelle (2000), Bourdelle et al. (2002),
and Bourdelle et al. (2007). Th underlying principles of the code regarding
the variational and action-angle approaches can be traced to Garbet et al.
(1990), and upgrades to the physics including angular momentum transport
(Cottier et al., 2014) and numerical improvements (Citrin et al., 2017) have
been made since its original development. The goal of QuaLiKiz is to calculate
the quasilinear transport that arise from microinstabilities. The core
principle is to develop linearize the kinetic equations and solve the
dispersion relation to find the complex frequencies for microinstabilities,
namely the ion temperature gradient (ITG), electron temperature gradient
(ETG), and trapped electron mode (TEM) instabilities. Upon solving the linear
problem, we then incorporate nonlinear physics to compute particle, angular
momentum, and heat fluxes. We do so via a quasilinear approach by coupling the
linear characteristics of the problem together and using previously performed
nonlinear kinetic simulations to saturate the perturbed state. Thus, while the
amplitudes of the modes are set by nonlinear physics, the key transport
features can be constructed from the linear regime. Quasilinear methods have
been shown to be valid in the tokamak core. Moreover, the quasilinear codes
are much faster than fully nonlinear kinetic codes. QuaLiKiz in particular can
perform a full computation in $\sim 1$ CPUs per wavenumber (Citrin et al.,
2017).

As a gyrokinetic code, QuaLiKiz is well suited to model the core of tokamak
devices which are strongly magnetized. Gyrokinetics is a popular approach to
investigate turbulent phenomena in magnetized plasmas such as those of fusion
devices (Brizard & Hahm, 2007; Cary & Brizard, 2009). Gyrokinetics is well
suited in scenarios where the microscopic dynamics are subject to the
gyrokinetic ordering. Essentially, we apply gyrokinetics to situations where
we can decouple the fast gyromotion of the charged particle from the slow
drift motion; this can be done when the time scale of the gyromotion is
significantly faster than all other time scales in the system and when the
gyroradius is smaller than almost all other length scales in the system. In
such a scenario, the magnetic moment is conserved, leading to a significant
reduction in the complexity of the dynamics (Stephens et al., 2017). Moreover,
gyrokinetics incorporates an ordering where the modes are anisotropic and
flute-like, meaning that the characteristic parallel wavelength of the mode is
large but perpendicular wavelengths can be comparable to the gyroradius. Thus,
gyrokinetics is well suited for theoretical and quantitative investigations of
magnetized plasma microturbulence. As a result, gyrokinetics has been used and
applied in a wide variety of systems (Wan et al., 2005; Rogers et al., 2007;
Wang et al., 2008; Pueschel et al., 2011, 2014; Howes et al., 2008, 2011; Told
et al., 2015; Navarro et al., 2016; Told et al., 2016). Even beyond tokamaks,
progress is being made in simulating stellarator plasmas in gyrokinetic codes
(Jenko & Kendl, 2002; Xanthopoulos & Jenko, 2007; Mynick et al., 2010; Nunami
et al., 2010; Baumgaertel et al., 2011). QuaLiKiz in particular, however,
assumes an axisymmetric geometry to simplify the dynamics, meaning QuaLiKiz is
only suitable for tokamaks and not stellarators.

Aside from the well established gyrokinetic approach, the key assumption
behind QuaLiKiz is the quasilinear approximation. In nonlinear simulations,
turbulent fluctuations eventually saturate due to coupling mechanisms between
different modes. However, it has been found that the nonlinear mode structure
can resemble the underlying linear mode structure; in particular, the cross
phases between fluctuating quantities in nonlinear simulations are identical
to that of linear simulations (Dannert & Jenko, 2005; Jenko et al., 2005). In
such situations, one also finds that ratios of the particle and heat fluxes
calculated in the linear regime match those calculated in the nonlinear regime
and that the real part of the nonlinear mode frequency resembles that of the
linear mode (Merz & Jenko, 2008; Görler & Jenko, 2008). Moreover, it has been
found that when different instabilities are found in the linear regime, their
interplay can manifest in the nonlinear regime (Merz & Jenko, 2010). This
motivates a quasilinear approach where the equilibrium distribution function
slowly evolves in comparison to the time scale of the instability, essentially
taking a mean field theory approach. Then, the linear response is acquired and
used to inform the the first order nonlinear behavior of the system.
Quasilinear flux ratios are then calculated and each flux is appropriately
saturated to the correct magnitude using a nonlinear saturation rule informed
by nonlinear physics. The approach allows us to exploit the fact that the
nonlinear state resembles the linear state to perform flux calculations
without needing to carry out a full nonlinear simulation (Citrin et al.,
2012).

However, constructing a quasilinear code instead of a nonlinear code is alone
not enough to increase the speed of calculations. Rather, a litany of
approximations and reductions are necessary. Aside from other typical
approximations for gyrokinetic tokamak codes (e.g. nonrelativistic particles,
quasineutrality), QuaLiKiz makes use of the following assumptions:

  * •

Adiabatic invariance. By exploiting the adiabatic invariants of the system, we
can formulate the Vlasov equation with action-angle variables. This requires
that the single-particle Hamiltonian be slowly varying in time in comparison
to the characteristic frequencies of motion. These frequencies correspond to
the cyclotron motion, the bounce-transit motion, and the toroidal drift and
precession.

  * •

Shifted Maxwellian with low Mach number and the $\delta f$ approximation.
QuaLiKiz linearizes the Vlasov equation by assuming a small perturbation from
the shifted Maxwellian. Although we include the effect of bulk plasma
rotation, we operate in the limit that the Mach number associated with the
rotation is small.

  * •

Electrostatic fluctuations. The code allows for electrostatic perturbations
and an equilibrium electric field. The absence of magnetic perturbations
allows for the exclusive use of Poisson’s equation while neglecting Ampere’s
law, thus simplifying the linear problem. To simplify the guiding center
motion, we require that the equilibrium electrostatic potential is small
compared to the characteristic thermal energy.

  * •

Trapped electron collisions. As an approximation, we utilize a Krook collision
operator for trapped electrons and neglect collisions entirely for passing
electrons and all ions.

  * •

Shifted circle geometry with small inverse aspect ratio. This simplified
geometry is used to calculate the magnetic drifts and perform integrals over
the pitch angle with ease. The $s-\alpha$ model gives rise to a radial shift
in the concentric flux surfaces called the Shafranov shift. The effect of this
shift is included when calculating the magnetic drifts, but ignored when
considering the bounce-transit motion. Thus, the treatment of guiding center
motion with respect to the geometry is inherently inconsistent. Moreover, the
$s-\alpha$ model is ad-hoc and does not solve the Grad-Shafranov equation.

  * •

Gaussian eigenfunctions. Instead of using a self-consistent eigenfunction for
the electrostatic modes, QuaLiKiz assumes the modes take the form of a
Gaussian. The shift and width of the Gaussian are calculated in the high mode
frequency limit as functions of the mode frequency, and substituted back into
the dispersion relation.

  * •

Strong ballooning. The electrostatic modes are assumed to be heavily localized
around their rational flux surface. This allows for a Fourier link between the
minor radius $r$ and the poloidal angle $\theta$, thus simplifying the
calculation. The localization also creates a separation of scales, thus
allowing the integrals to be more easily approximated.

  * •

Strongly passing and strongly trapped particles. Trapped and passing particles
are considered to be respectively strongly trapped and strongly passing. For
trapped particles, this greatly simplifies the relation between the physical
toroidal and poloidal angles and the action angles and leads to a kinetic
bounce average that is similar to the gyro-average. For passing particles, the
strongly passing assumption simplifies the integrals over the pitch angle due
to the dominating parallel velocity.

The goal of this work is to derive the analytic equations for QuaLiKiz step by
step. Although various overviews of the QuaLiKiz and Kinezero framework have
already been published (Bourdelle, 2015; Bourdelle et al., 2016; Citrin et
al., 2017), no combination of currently published works derive the entirety of
the model from first principles. We seek to fill this gap by offering a
comprehensive and complete formulation of QuaLiKiz. This work will then as a
result serve as a guide for improving upon QuaLiKiz and attaining physical and
mathematical intuition as to its key principles, approximations, and
computational methods. In addition, we also outline the new computational
method used to numerically calculate 1-dimensional and 2-dimensional
integrals. Moreover, this sort of work serves as a tutorial for those seeking
to understand the fundamental considerations in the formulation of any
quasilinear tokamak code. While many codes offer comprehensive manuals and
describe the key principles at play, the process of creating such a code from
scratch can often appear opaque and unintuitive. Thus, this derivation also
serves as a tutorial for those who seek to understand the physical,
mathematical, and computational aspects of quasilinear modeling in all their
gory details.

The paper is organized as follows: Section 2 reviews the action-angle
formalism and derives explicit expressions for the action-angle variables from
physical variables. In Section 3, we linearize the Vlasov equation and expand
the perturbed distribution function and electrostatic potential using a
Fourier series to derive the dispersion relation. To solve the dispersion
relation, we must integrate over all of phase space, resulting in a functional
that depends on the complex frequency of the mode. Section 4 examines the
ballooning transform and its role in simplifying the dispersion relation as
well as the characteristics of the electrostatic perturbation. Sections 5 \- 7
calculate the adiabatic, trapped, and passing parts of the functional,
respectively, resulting in a reduced expression for the dispersion relation.
Section 8 applies these results to the quasilinear problem to derive
expressions for the particle, toroidal angular momentum, and heat fluxes.
Section 9 connects the quasilinear results with nonlinear physics with the use
of a saturation rule. Section 10 explains the method of contour integration
used in QuaLiKiz to find the eigenmodes and the newly implemented numerical
integration method based on the Genz and Malik algorithm (Genz & Malik, 1980).
Finally, we summarize our work Section 11. We include Appendix A to serve as a
brief explanation of Fried and Conte integrals. In addition, we derive the
magnetic drift velocity in an $s-\alpha$ equilibrium in Appendix B and briefly
discuss the inclusion trapped electron collisions in Appendix C. The
derivation is performed in SI units, and we set the Boltzmann constant
$k_{B}=1$ such that our temperatures are in units of energy.

##  2 Action Angle Variables

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

$\frac{\partial f}{\partial t}+\left\\{f,H\right\\}=0,$ (1)

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

$\frac{\partial f}{\partial t}+\dot{\bm{q}}\cdot\frac{\partial
f}{\partial\bm{q}}+\dot{\bm{p}}\cdot\frac{\partial f}{\partial\bm{p}}=0,$ (2)

where, for a single particle $\bm{q}$ is the position, $\bm{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
(1972) and expanded upon by Mahajan & Chen (1985). The core principle is to
define a canonical transformation,

$(\bm{q},\bm{p})\to(\bm{\alpha},\bm{J}),$ (3)

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

$\displaystyle\frac{\partial H}{\partial\bm{\alpha}}$
$\displaystyle=-\dot{\bm{J}}=\bm{0},$ (4) $\displaystyle\frac{\partial
H}{\partial\bm{J}}$ $\displaystyle=\dot{\bm{\alpha}}=\bm{\Omega},$ (5)

where $\bm{\Omega}$ are the constant angular frequencies associated with each
adiabatic invariant. 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 analyzing the unperturbed system
and then including electromagnetic fluctuations in the Hamiltonian.

We define the Hamiltonian to be

$H=H_{0}+\delta h,$ (6)

where the unperturbed Hamiltonian is simply

$H_{0}=\frac{1}{2m}\left(\bm{p}^{2}-e\bm{A}_{0}^{2}\right)+e\Phi.$ (7)

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

$\delta h=e\phi,$ (8)

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

$\displaystyle\frac{\partial H_{0}}{\partial\bm{\alpha}}$
$\displaystyle=\bm{0},$ (9) $\displaystyle\frac{\partial
H_{0}}{\partial\bm{J}}$ $\displaystyle=\bm{\Omega}.$ (10)

Hamilton’s equations of motion then become

$\displaystyle\frac{\partial H}{\partial\bm{\alpha}}$
$\displaystyle=-\dot{\bm{J}}=e\frac{\partial\phi}{\partial\bm{\alpha}},$ (11)
$\displaystyle\frac{\partial H}{\partial\bm{J}}$
$\displaystyle=\dot{\bm{\alpha}}=\bm{\Omega}+e\frac{\partial\phi}{\partial\bm{J}}.$
(12)

We note that because the unperturbed Hamiltonian is a function of $\bm{J}$ and
not $\bm{\alpha}$, all equilibrium quantities are also only functions of
$\bm{J}$. Furthermore, any function of $\bm{\alpha}$ is periodic with respect
to $\bm{\alpha}$; thus, the perturbed quantities in our system admit a Fourier
series expansion. Moreover, it can be shown that $\bm{\alpha},\bm{J}$ are
canonical coordinates even after introducing a perturbation (Mahajan & Chen,
1985). 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 $\bm{r}$ and the velocity
$\bm{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:

$\displaystyle J_{1}=\frac{m}{e}\mu,$ (13) $\displaystyle
J_{2}=\frac{1}{2\pi}\oint
mp_{\parallel}dl=\frac{1}{2\pi}\oint\left(mv_{\parallel}+eA_{\parallel}\right)dl,$
(14) $\displaystyle
J_{3}=P_{\varphi}=\frac{mv_{\parallel}R_{0}B_{\varphi}^{0}}{B}-e\psi.$ (15)

Here, $\mu=W_{\perp}/B$ is the magnetic moment, where
$W_{\perp}=\frac{1}{2}mv_{\perp}^{2}$ is the kinetic energy associated with
the velocity perpendicular to the magnetic field $\bm{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\pi$, which is
calculated by integrating the flux of the magnetic field through a disk
tangent to the flux surface everywhere:

$\psi=-\frac{1}{2\pi}\int_{S}\bm{B}\cdot d\bm{S}.$ (16)

The following subsections 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
$\bm{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{\bm{r}}\times\hat{\bm{\theta}}=\hat{\bm{\varphi}}$. For further
references characterizing the action angles $\bm{\alpha}$, we refer the reader
to Garbet’s work in Garbet (2001) and Garbet et al. (1990).

###  2.1 Magnetic Moment

In the presence of a magnetic field, charged particles gyrate about the field
line at the cyclotron frequency $\Omega_{1}=eB/m$. With a strong enough
magnetic field, the cylcotron frequency is much larger than any other
characteristic frequency in the plasma. Under such conditions, the magnetic
moment $\mu$ is adiabatically conserved (Kruskal, 1962; Baños, 1967; Brizard &
Hahm, 2007; Cary & Brizard, 2009; Stephens et al., 2017), and the gyromotion
can be decoupled from the guiding center motion:

$\displaystyle r=r_{G}+\rho\cos(\alpha_{1}),$ (17)
$\displaystyle\theta=\theta_{G}+\frac{\rho}{r}\sin(\alpha_{1}),$ (18)
$\displaystyle\varphi=\varphi_{G},$ (19)

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 center.
These guiding center variables obey the guiding center 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 to 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 take average various functions over
the gyrophase $\alpha_{1}$ by integrating over $\alpha_{1}$. We therefore
consider the general integral

$g_{n_{1}}=\int_{-\pi}^{\pi}\frac{d\alpha_{1}}{2\pi}g(\bm{r})e^{-in_{1}\alpha_{1}},$
(20)

where $n_{1}$ is an integer. It will be shown that later that factors of
$e^{-in_{1}\alpha_{1}}$ arise from taking Fourier expansions in terms of
$\alpha_{1}$. We define the Fourier transform of $g$ to be

$\tilde{g}(\bm{k})=\int d^{3}rg(\bm{r})e^{i\bm{k}\cdot\bm{r}},$ (21)

with the corresponding inverse Fourier transform

$g(\bm{r})=\int\frac{d^{3}k}{(2\pi)^{3}}\tilde{g}(\bm{k})e^{-i\bm{k}\cdot\bm{r}}.$
(22)

We use the Fourier transform to obtain

$g_{n_{1}}=\int_{-\pi}^{\pi}\frac{d\alpha_{1}}{2\pi}\int\frac{d^{3}k}{(2\pi)^{3}}\tilde{g}(\bm{k})e^{-i\bm{k}\cdot\bm{r}-in_{1}\alpha_{1}}=\int_{-\pi}^{\pi}\frac{d\alpha_{1}}{2\pi}\int\frac{d^{3}k}{(2\pi)^{3}}\tilde{g}(\bm{k})e^{-i\bm{k}\cdot\bm{\rho}-in_{1}\alpha_{1}}e^{-i\bm{k}\cdot\bm{R}_{G}}.$
(23)

Here, we have decoupled the gyromotion from the guiding center motion via
$\bm{r}=\bm{R}_{G}+\bm{\rho}$. We then write

$\bm{k}\cdot\bm{\rho}=k_{\perp}\rho\cos(\alpha_{1}),$ (24)

where

$k_{\perp}=\left|\bm{k}-\bm{k}\cdot\hat{\bm{b}}\right|\approx\sqrt{k_{r}^{2}+k_{\theta}^{2}}.$
(25)

Note that according to our definition of the Fourier transform, $k_{r}$ and
$k_{\theta}$ are operators in real space such that

$\displaystyle k_{r}$ $\displaystyle\to i\frac{\partial}{\partial r},$ (26)
$\displaystyle k_{\theta}$
$\displaystyle\to\frac{i}{r}\frac{\partial}{\partial\theta}.$ (27)

We may then integrate over $\alpha_{1}$ independently, leading to

$\int^{\pi}_{-\pi}\frac{d\alpha_{1}}{2\pi}e^{-ik_{\perp}\rho\cos(\alpha_{1})-in_{1}\alpha_{1}}=(-i)^{n_{1}}J_{n_{1}}(k_{\perp}\rho),$
(28)

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

$g_{n_{1}}=\int\frac{d^{3}k}{(2\pi)^{3}}(-i)^{n_{1}}J_{n_{1}}(k_{\perp}\rho)\tilde{g}(\bm{k})e^{-i\bm{k}\cdot\bm{R}_{G}}=(-i)^{n_{1}}\left(J_{n_{1}}(k_{\perp}\rho)\cdot
g\right)(\bm{R}_{G}).$ (29)

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 center. 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 center. Thus, we drop the subscript “G" for
convenience and treat all spatial variables as those corresponding to the
guiding center. The adiabatic invariants $J_{2}$ and $J_{3}$ are explicitly
calculated within the guiding center framework where we hold $\mu$ constant
and ignore the cyclotron motion.

###  2.2 Longitudinal Invariant

To calculate $J_{2}$, we consider the guiding center particle motion along a
magnetic field line; such a particle completes bounce-transit orbits with
frequency $\Omega_{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 center drifts by holding $r$ constant.
For an extended treatment of bounce-transit motion, see Brizard (2011) and
Stephens et al. (2020).

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

$\left|v_{\parallel}\right|=\sqrt{\frac{2}{m}(E-\mu B)},$ (30)

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

$\bm{B}=B_{\varphi}(r,\theta)\hat{\bm{\bm{\varphi}}}+B_{\theta}(r,\theta)\hat{\bm{\bm{\theta}}}=\frac{1}{1+r/R_{0}\cos(\theta)}\left(B^{0}_{\varphi}(r)\hat{\bm{\bm{\varphi}}}+B^{0}_{\theta}(r)\hat{\bm{\bm{\theta}}}\right),$
(31)

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 which 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

$\frac{\mu B^{0}(r)}{E}\geq 1-\epsilon.$ (32)

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

$v_{\parallel}=\sqrt{\frac{2T}{m}}\epsilon_{\parallel}\sqrt{\xi}\sqrt{1-\lambda
b(r,\theta)}.$ (33)

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

$\displaystyle\lambda$ $\displaystyle=\frac{\mu
B\left(r,\theta=0\right)}{E}=\frac{v_{\perp}^{2}}{v^{2}b(r,\theta)},$ (34)
$\displaystyle b(r,\theta)$ $\displaystyle=\frac{B(r,\theta)}{B(r,\theta=0)}.$
(35)

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

$\left|\Omega_{2}\right|=\frac{2\pi}{T_{2}},$ (36)

where

$T_{2}=\oint\frac{d\theta}{\left|\frac{d\theta}{dt}\right|}.$ (37)

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 $-\pi$ to $\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{\bm{b}}$, the direction of the magnetic field, is approximately
$\hat{\bm{\bm{\varphi}}}$. We again emphasize that this approximation breaks
down for devices such as spherical tokamaks. Therefore, we write

$v_{\parallel}=\bm{v}\cdot\hat{\bm{b}}\approx\dot{\varphi}\left(R_{0}+r\cos(\theta)\right)\approx
qR_{0}\dot{\theta}.$ (38)

Here, we have defined the safety factor

$q(r)=\frac{d\varphi}{d\theta}\approx\frac{rB_{\varphi}}{R_{0}B_{\theta}}.$
(39)

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

$\left|\Omega_{2}\right|=\sqrt{\frac{2T}{m}}\frac{\sqrt{\xi}}{qR_{0}}\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu_{2}(r,\lambda),$
(40)

where we define

$\mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}(r,\lambda)=\frac{2\pi}{\oint d\theta\frac{1}{\sqrt{1-\lambda
b(r,\theta)}}}.$ (41)

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 $\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu_{2}$ in the small
$\epsilon$ limit to

$\mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}\approx\left\\{\begin{aligned}
&\frac{\pi\sqrt{\epsilon}}{2\sqrt{2}K(\kappa)}&&\text{ if }0\leq\kappa<1\text{
(trapped)},\\\
&\frac{\pi\kappa\sqrt{\epsilon}}{\sqrt{2}K\left(\kappa^{-1}\right)}&&\text{ if
}1<\kappa<\infty\text{ (passing)}.\end{aligned}\right.$ (42)

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

$\lambda=1-2\epsilon\kappa^{2}.$ (43)

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

$J_{2}\approx\left\\{\begin{aligned}
&\frac{8mqR_{0}\sqrt{E/m}\sqrt{\epsilon}}{\pi}\left(E\left(\kappa\right)-\left(1-\kappa^{2}\right)K\left(\kappa\right)\right)&&\text{
if }0\leq\kappa<1\text{ (trapped)},\\\
&\frac{4mqR_{0}\sqrt{E/m}\sqrt{\epsilon}}{\pi}\kappa
E\left(\kappa^{-1}\right)+e\Phi_{t}&&\text{ if }1<\kappa<\infty\text{
(passing)},\end{aligned}\right.$ (44)

where $E$ is the complete elliptic integral of the second kind and $\Phi_{t}$
is the toroidal flux normalized by $2\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

$\frac{d\alpha_{2}}{dt}=\Omega_{2}.$ (45)

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,

$\frac{d\alpha_{2}}{dt}=\frac{d\alpha_{2}}{d\theta}\dot{\theta}=\frac{d\alpha_{2}}{d\theta}\frac{v_{\parallel}}{qR_{0}}=\Omega_{2}.$
(46)

We emphasize that $\Omega_{2}$ is not dependent on $\alpha_{2}$ or $\theta$.
Thus, this differential equation can be integrated using elliptic functions,
leading to an expression of $\alpha_{2}$ in terms of $\theta$. We use the
convention that $\alpha_{2}(\theta=0)=0$, leading to

$\alpha_{2}=\int^{\theta}_{0}d\theta^{\prime}\frac{qR_{0}\Omega_{2}}{v_{\parallel}}=\int^{\theta}d\theta^{\prime}\frac{qR_{0}\Omega_{2}}{\sqrt{2T/m}\epsilon_{\parallel}\sqrt{\xi}\sqrt{1-\lambda
b}}.$ (47)

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, leading to

$\alpha_{2}\approx\left\\{\begin{aligned}
&\frac{\pi}{2}\frac{F\left(\sin^{-1}\left[\kappa\sin\left(\frac{\theta}{2}\right)\right],\kappa\right)}{K\left(\kappa\right)}&&\text{
if }0\leq\kappa<1\text{ (trapped)},\\\
&\pi\frac{F\left(\frac{\theta}{2},\kappa^{-1}\right)}{K\left(\kappa^{-1}\right)}&&\text{
if }1<\kappa<\infty\text{ (passing)},\end{aligned}\right.$ (48)

where $F$ is the incomplete integral of the first kind. Essentially, the
integral takes the same form as when calculating $\Omega_{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 $\pi$ for passing
particles.

Finally, let $G\left(\epsilon_{\parallel},\theta\right)$ 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 $\left<G\left(\epsilon_{\parallel},\theta\right)\right>$ to be

$\left<G\left(\epsilon_{\parallel},\theta\right)\right>=\frac{1}{T_{2}}\oint\frac{d\theta}{\left|\frac{d\theta}{dt}\right|}G\left(\epsilon_{\parallel},\theta\right)=\frac{\oint
d\theta\frac{G\left(\epsilon_{\parallel},\theta\right)}{\sqrt{1-\lambda
b}}}{\oint d\theta\frac{1}{\sqrt{1-\lambda b}}}.$ (49)

For passing particles, the average is explicitly

$\oint d\theta\frac{G\left(\epsilon_{\parallel},\theta\right)}{\sqrt{1-\lambda
b}}=\int_{-\pi}^{\pi}d\theta\frac{G\left(\epsilon_{\parallel},\theta\right)}{\sqrt{1-\lambda
b}},$ (50)

while for trapped particles the average is instead

$\oint d\theta\frac{G\left(\epsilon_{\parallel},\theta\right)}{\sqrt{1-\lambda
b}}=\int_{-\theta_{b}}^{\theta_{b}}d\theta\frac{G\left(\epsilon_{\parallel},\theta\right)+G\left(-\epsilon_{\parallel},\theta\right)}{\sqrt{1-\lambda
b}}.$ (51)

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 the next section,
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 center theory, we can write the
canonical toroidal momentum as

$P_{\varphi}=\frac{mv_{\parallel}R_{0}B_{\varphi}^{0}}{B}-e\psi.$ (52)

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
$\left(P_{\varphi}+e\psi\right)/\left(e\psi\right)\sim\sqrt{mT}/(eBR_{0})$.
Inputting JET-like parameters, $T=5$ keV, $m=m_{D},B=3$ T, $R_{0}=3$ m, then
$\sqrt{mT}/(eBR_{0})\sim 10^{-3}$, making this a very reasonable
approximation. We therefore write

$J_{3}=-e\psi.$ (53)

To calculate the poloidal flux, we utilize Stoke’s theorem; the surface
integral of
$\bm{B$}simplybecomesaclosedlineintegralof$A$tofind\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))dr^{\prime}=\frac{\int^{r}B^{0}_{\theta}(r^{\prime})dr^{\prime}}{1+\epsilon\cos\left(\theta\right)}.\end{equation}Thus,weobtain\begin{equation}J_{3}=-e\psi=-\frac{1}{2\pi}\int_{0}^{2\pi}e(1+\epsilon\cos(\theta))A_{\varphi}d\varphi=e\int^{r}R_{0}B_{\theta}^{0}(r^{\prime})dr^{\prime}.\end{equation}Weseethenthat$J_3$ispurelyafunctionof$r$suchthat\begin{equation}\frac{dJ_{3}}{dr}\approx-\frac{erB}{q}.\end{equation}\par
Wenextcalculate$Ω_3$,whichisthetoroidalprecessionfrequencyfortrappedparticlesandthetoroidalrotationfrequencyforpassingparticles.Alongthebounce-
transitorbit,guidingcenterdriftscauseradialexcursionsfromthemagneticfieldline.Inaddition,passingparticleswindaroundthemagneticfieldlinetoroidallyduetothelackofanybouncepoint.Tocalculatethisfrequency,weneedtofirstcalculatedeviationsfromthefieldlineorbit,notingthatradialexcursionsfromthefieldlineareoftheorderofthegyroradius.Toaidinthecalculation,weexploittheexactconservationofthecanonicaltoroidalmomentum:\begin{equation}\psi=\bar{\psi}+\frac{mB^{0}_{\varphi}R_{0}v_{\parallel}}{eB}=\bar{\psi}+\psi_{1}.\end{equation}Here,$¯ψ$correspondstothereferencemagneticfluxsurfacedefinedtobe\begin{equation}\bar{\psi}=-\frac{P_{\varphi}}{e}\end{equation}and$ψ_1$isthedeviationfromthatfluxsurface.Recallthatthefieldlineorbitassumedthat\begin{equation}\frac{d\theta}{dt}\approx\frac{v_{\parallel}}{qR_{0}},\end{equation}wherewehold$r$andthus$ψ$fixed.Sincetheexactfieldline-
followingorbitbreaks$P_φ$conservation,weneedtoincludedeviationsfromthefieldlinecausedbyconservationof$P_φ$alongwithguidingcenterdriftsinordertoconsistentlyexpandtheguidingcenterequationofmotionwithrespecttothegyroradius.\par
Theguidingcenterequationofmotionis\begin{equation}\frac{dx}{dt}=\left(v_{\parallel}\hat{\bm{b}}+\bm{v}_{D}\right)\cdot\nabla
x,\end{equation}where$x$isanyspatialcoordinateand$v_D$aretheguidingcenterdrifts.Wethenexpandtheguidingcenterequationofmotionforvariables$ψ,θ$,and$φ$andfindthat\begin{aligned}
\frac{d\psi}{dt}&\approx\left.\left(\bm{v}_{D}\cdot\nabla\psi\right)\right|_{\bar{\psi}},\\\
\frac{d\theta}{dt}&\approx\left.\left(\frac{v_{\parallel}}{qR_{0}}\right)\right|_{\bar{\psi}}+\frac{d}{d\psi}\left.\left(\frac{v_{\parallel}}{qR_{0}}\right)\right|_{\bar{\psi}}\psi_{1}+\left.\left(\bm{v}_{D}\cdot\nabla\theta\right)\right|_{\bar{\psi}},\\\
\frac{d\varphi}{dt}&\approx\left.\left(\frac{v_{\parallel}}{qR_{0}}\frac{d\varphi}{d\theta}\right)\right.+\frac{d}{d\psi}\left.\left(\frac{v_{\parallel}}{qR_{0}}\frac{d\varphi}{d\theta}\right)\right|_{\bar{\psi}}\psi_{1}+\left.\left(\bm{v}_{D}\cdot\nabla\varphi\right)\right|_{\bar{\psi}}.\end{aligned}Wetakethe$v_D$tobethethesumoftheclassicalcurvature,grad-$B$,and$E$-cross-$B$drifts:\begin{equation}\bm{v}_{D}=\frac{\hat{\bm{b}}}{eB}\left(mv_{\parallel}^{2}\bm{k}+\mu\nabla
B+e\nabla\Phi\right),\end{equation}where$^k$isthecurvaturevectordefinedsuchthat\begin{equation}\hat{\bm{b}}\times\bm{k}=\nabla\times\hat{\bm{b}}-\left(\hat{\bm{b}}\cdot\nabla\times\hat{\bm{b}}\right)\hat{\bm{b}}.\end{equation}Wecansimplifytheequationofmotioninthetoroidaldirectionbynotingthat\begin{equation}\begin{split}\frac{d}{d\psi}\left(\frac{v_{\parallel}}{qR_{0}}\frac{d\varphi}{d\theta}\right)\psi_{1}&=\frac{d}{d\psi}\left(\frac{v_{\parallel}}{qR_{0}}\right)\frac{d\varphi}{d\theta}\psi_{1}+\frac{d}{d\psi}\left(\frac{d\varphi}{d\theta}\right)\frac{v_{\parallel}}{qR_{0}}\psi_{1}\\\
&=\frac{d\varphi}{d\theta}\frac{d\theta}{dt}+\frac{d}{d\psi}\left(\frac{d\varphi}{d\theta}\right)\frac{v_{\parallel}}{qR_{0}}\psi_{1}-\frac{d\varphi}{d\theta}\frac{v_{\parallel}}{qR_{0}}-\frac{d\varphi}{d\theta}\left(\bm{v}_{D}\cdot\nabla\theta\right),\end{split}\end{equation}whereweusedtheequationofmotioninthepoloidaldirection.Substitutingthisinandevaluating$d
φdθ$,weobtain\begin{equation}\dot{\varphi}=\frac{d}{d\psi}\left(\frac{q}{1+\epsilon\cos\left(\theta\right)}\right)\psi_{1}\frac{v_{\parallel}}{qR_{0}}+\bm{v}_{D}\cdot\nabla\varphi-\frac{q}{1+\epsilon\cos\left(\theta\right)}\bm{v}_{D}\cdot\nabla\theta+\frac{q}{1+\epsilon\cos\left(\theta\right)}\frac{d\theta}{dt},\end{equation}whereweevaluateallradialcoordinatesat$¯r$suchthat\begin{equation}\psi\left(\bar{r}\right)=\bar{\psi}.\end{equation}Finally,wetakethebounce-
transitaverageof$˙φ$andfindthat\begin{equation}\Omega_{3}=\left<\dot{\varphi}\right>=\left<\omega_{d}\right>+\bar{\epsilon}q(\bar{r})\Omega_{2}=\Omega_{d}+\bar{\epsilon}q(\bar{r})\Omega_{2}.\end{equation}Here,$¯ϵ$is0fortrappedparticlesand1forpassingparticles,$Ω_d$isthefrequencypurelyduetotheguidingcenterdrifts,and$ω_d$isassociatedwiththeinstantaneousdeviationfromthemagneticfieldline.Theextratermforpassingparticlesisduetothetoroidalrotationfromfollowingthefieldlineinacompletepoloidalturn.Thisparallelvelocitydependenttermisabsentfortrappedparticlessincetheiraveragetoroidalpositiondoesnotchangeasaresultofacompletefieldline-
followingbounce.Weapproximate$ω_d$as\begin{equation}\omega_{d}\approx\frac{dq}{d\psi}\frac{\psi_{1}v_{\parallel}}{q}+\bm{v}_{D}\cdot\nabla\varphi-q\bm{v}_{D}\cdot\nabla\theta.\end{equation}Thepoloidalcomponentofthemagneticdriftdominates,thusweignorethetoroidalcomponent.Usingthe$s-α$equilibrium,wecalculatetheguidingcenterdriftinAppendix~{}\ref{Magnetic}.Thus,$Ω_3$iscomputedwithafiniteShafranovshift.Thepoloidalcomponentoftheguidingcenterdriftis\begin{equation}\bm{v}_{D}\cdot\nabla\theta\approx-\frac{m}{eBrR_{0}}\left(v_{\parallel}^{2}+\frac{v_{\perp}^{2}}{2}\right)\left(\cos\left(\theta\right)-\alpha\sin^{2}\left(\theta\right)\right)-\frac{E_{r}}{rB},\end{equation}where$E_r$istheradialelectricfield.Wealsodefine$α$suchthat\begin{equation}\alpha=-q^{2}\beta\frac{R_{0}}{P}\frac{dP}{dr}.\end{equation}Here,$
β= 2 μ_0
P/B^2$with$μ_0$beingthevacuumpermeability.The$E$-cross-$B$driftcanbeseparatedfromthemagneticdrifts,sothatweobtain\begin{equation}\bm{v}_{D,B}\cdot\nabla\theta\approx-\frac{v_{D,B}}{r}\left(\cos\left(\theta\right)-\alpha\sin^{2}\left(\theta\right)\right),\end{equation}where$v_D,B$characterizesthemagneticdriftsandcanbewrittenas\begin{equation}v_{D,B}=\frac{m}{eBR_{0}}\left(v_{\parallel}^{2}+\frac{v_{\perp}^{2}}{2}\right)=\frac{\xi
T}{eBR_{0}}\left(2-\lambda
b\right),\end{equation}wherethemagneticshearis\begin{equation}s=\frac{r}{q}\frac{dq}{dr}.\end{equation}Wecanthereforerewrite$
Ω_d$as\begin{equation}\Omega_{d}=\frac{q}{r}\frac{T\xi}{eBR_{0}}\left<\left(2-\lambda
b\right)\left(\cos(\theta)-\alpha\sin^{2}(\theta)\right)+\frac{2s}{\epsilon}\left(1-\lambda
b\right)\right>+\left<\frac{qE_{r}}{rB}\right>=\omega_{d0}\xi
F_{d}(\kappa)+\omega_{E},\end{equation}where$F = F_d(κ)$isthebounce-
transitaveragedterm,$ω_d0$ischaracteristicofthemagneticprecessionfrequencyanddefinedtobe\begin{equation}\omega_{d0}=\frac{qT}{reBR_{0}},\end{equation}and$ω_E$correspondstothe$E$-cross-$B$velocity.Weexplicitlycarryoutthebounce-
transitaveragebyrewritingthe$λb$termsintermsof$κ$and$θ$andfindthat\begin{equation}F_{d}(\kappa)\approx\left\\{\begin{aligned}
&-1+\frac{2E\left(\kappa\right)}{K\left(\kappa\right)}+4s\left(\kappa^{2}-1+\frac{E\left(\kappa\right)}{K\left(\kappa\right)}\right)&&\\\
&\qquad\qquad-\frac{4\alpha}{3}\left(1-\kappa^{2}-\left(1-2\kappa^{2}\right)\frac{E\left(\kappa\right)}{K\left(\kappa\right)}\right)&&\text{
if }0\leq\kappa<1\text{ (trapped)},\\\
&-\left(2\kappa^{2}-1\right)+2\kappa^{2}\frac{E\left(\kappa^{-1}\right)}{K\left(\kappa^{-1}\right)}+4s\kappa^{2}\frac{E\left(\kappa^{-1}\right)}{K\left(\kappa^{-1}\right)}&&\\\
&\qquad\qquad-\frac{4\alpha}{3}\kappa^{2}\left(\left(2\kappa^{2}-1\right)\frac{E\left(\kappa^{-1}\right)}{K\left(\kappa^{-1}\right)}-2\kappa^{2}+2\right)&&\text{
if }1<\kappa<\infty\text{ (passing)}.\end{aligned}\right.\end{equation}\par
Ratherthancalculating$α_3$explicitly,itsufficestowriteitsgenericintegralform.Wealsoincludethevariousoscillatingquantitiesassociatedwiththeprecessionmotion.Ingeneral,wehave\begin{aligned}
&\psi=\bar{\psi}+\tilde{\psi},\\\
&\varphi=\alpha_{3}+q(\bar{r})\tilde{\theta}+\tilde{\varphi},\\\
&\theta=\bar{\epsilon}\alpha_{2}+\tilde{\theta}.\end{aligned}Here,$~ψ$representstheexcursionfromthereferencefluxsurface$¯ψ$duringthepoloidalorbit.Meanwhile,$~φ$isthedifferenceintoroidalprecessionbetweenacirculargeometryandamoregeneralequilibriummagneticfield.Later,wewilluse$~r$insteadof$~ψ$withtheunderstandingthat$ψ(~r)
=
~ψ$.Meanwhile,$~θ$isassociatedwiththeoscillatorypoloidalmotion.Wedefinethesequantitiesas\begin{aligned}
\tilde{\psi}&=\int^{\alpha_{2}}\frac{d\alpha_{2}^{\prime}}{\Omega_{2}}\bm{v}_{D}\cdot\nabla\psi=\psi_{1}\\\
\tilde{\theta}&=\int^{\alpha_{2}}d\alpha_{2}^{\prime}\left(\frac{1}{\Omega_{2}}\frac{d\theta}{dt}-\bar{\epsilon}\right),\\\
\tilde{\varphi}&=\frac{dq}{d\psi}\tilde{\theta}\tilde{\psi}+\int^{\alpha_{2}}\frac{d\alpha_{2}^{\prime}}{\Omega_{2}}\left(\bar{\epsilon}\frac{dq}{d\psi}\tilde{\psi}\Omega_{2}+\bm{v}_{D}\cdot\nabla\varphi-\frac{dq}{d\psi}\tilde{\theta}\bm{v}_{D}\cdot\nabla\psi-q\bm{v}_{D}\cdot\nabla\theta-\Omega_{d}\right).\\\
\end{aligned}Thisguaranteesthat\begin{equation}\frac{d\alpha_{3}}{dt}=\Omega_{d}+\bar{\epsilon}q\Omega_{2}.\end{equation}Notethatintheabove,$q$and$dq/dψ$areevaluatedat$¯r$andthusaretimeindependent.Havingcharacterizedaction-
anglecoordinates,wecanproceedtosolvingtheVlasovequationusingthesecoordinates.\par$

##  3 The Vlasov Equation

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

$\frac{df}{dt}=\frac{\partial f}{\partial
t}+\dot{\bm{\alpha}}\cdot\frac{\partial
f}{\partial\bm{\alpha}}+\dot{\bm{J}}\cdot\frac{\partial f}{\partial\bm{J}}=0.$
(77)

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

$\displaystyle\dot{\bm{J}}=-\frac{\partial
H}{\partial\bm{\alpha}}=-e\frac{\partial\phi}{\partial\bm{\alpha}},$ (78)
$\displaystyle\dot{\bm{\alpha}}=\frac{\partial
H}{\partial\bm{J}}=e\frac{\partial\phi}{\partial\bm{J}}+\bm{\Omega},$ (79)

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}(\bm{J})$
and a perturbed part $\delta f=\delta f(\bm{\alpha},\bm{J},t).$ Dropping any
quadratic perturbative terms, we obtain

$\frac{\partial\delta f}{\partial t}+\bm{\Omega}\cdot\frac{\partial\delta
f}{\partial\bm{\alpha}}-e\frac{\partial\phi}{\partial\bm{\alpha}}\cdot\frac{\partial
f_{0}}{\partial\bm{J}}=0.$ (80)

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

$\displaystyle\delta
f=\sum_{\bm{n}}f_{\bm{n}}(\bm{J})e^{i(\bm{n}\cdot\bm{\alpha}-\omega t)},$ (81)
$\displaystyle\phi=\sum_{\bm{n}}\phi_{\bm{n}}(\bm{J})e^{i(\bm{n}\cdot\bm{\alpha}-\omega
t)}.$ (82)

To extract the physical quantity, we take the real part of the Fourier series.
Here, $\bm{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}+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

$\displaystyle f_{\bm{n}}=\int\frac{d^{3}\alpha}{(2\pi)^{3}}\delta
f(t=0)e^{-i\bm{n}\cdot\bm{\alpha}},$ (83)
$\displaystyle\phi_{\bm{n}}=\int\frac{d^{3}\alpha}{(2\pi)^{3}}\phi(t=0)e^{-i\bm{n}\cdot\bm{\alpha}},$
(84)

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

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

$f_{0}(\bm{J})=n_{0}\left(\frac{m}{2\pi
T}\right)^{3/2}\exp\left(-\frac{m\left(\bm{v}-\bm{U}\right)^{2}}{2T}\right)=n_{0}\left(\frac{m}{2\pi
T}\right)^{3/2}\exp\left(\frac{-H_{0}+e\Phi+m\bm{v}\cdot\bm{U}-\frac{mU^{2}}{2}}{T}\right).$
(85)

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

$U_{\varphi}\hat{\bm{\varphi}}\approx U_{\parallel}\hat{\bm{b}},$ (86)

which allows us to write

$\displaystyle U^{2}$ $\displaystyle\approx U_{\parallel}^{2},$ (87)
$\displaystyle\bm{v}\cdot\bm{U}$ $\displaystyle\approx
v_{\parallel}U_{\parallel}.$ (88)

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 natural frequency parameter for the electric
field shear $\gamma_{E}$ defined as

$\gamma_{E}=-\frac{1}{B}\frac{dE_{r}}{dr}.$ (89)

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

$\omega_{E}\approx\omega_{E0}+\omega_{E}^{\prime}x,$ (90)

where we expand about the radial distance $x=0$ and
$\partial_{r}\omega_{E}=\omega_{E}^{\prime}$ 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

$mU_{\parallel}\left|\frac{\partial
v_{\parallel}}{\partial\bm{J}}\right|\ll\bm{\omega}_{\ast}.$ (91)

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., 2013).

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_{\bm{n}}$ in
terms of $\phi_{\bm{n}}$. The result is

$f_{\bm{n}}=\frac{e\phi_{\bm{n}}\bm{n}\cdot\frac{\partial
f_{0}}{\partial\bm{J}}}{\bm{n}\cdot\bm{\Omega}-\omega}=\frac{f_{0}}{T}\frac{e\phi_{\bm{n}}\bm{n}\cdot\left(\bm{\omega}_{\ast}+\bm{\omega}_{E}-\bm{\Omega}\right)}{\bm{n}\cdot\bm{\Omega}-\omega},$
(92)

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

$\bm{\omega}_{\ast}=T\left(\frac{1}{n_{0}}\frac{\partial
n_{0}}{\partial\bm{J}}+\left(\xi-\frac{3}{2}-\frac{U_{\parallel}}{v_{T}^{2}}\left(2v_{\parallel}-U_{\parallel}\right)\right)\frac{1}{T}\frac{\partial
T}{\partial\bm{J}}\right)+\frac{2\left(v_{\parallel}-U_{\parallel}\right)}{v_{T}^{2}}\frac{\partial
U_{\parallel}}{\partial\bm{J}},$ (93)

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

$\bm{\omega}_{E}=\frac{e}{T}\frac{d\Phi}{d\bm{J}}.$ (94)

We then rewrite the equation to be

$f_{\bm{n}}=-\frac{e\phi_{\bm{n}}}{T}f_{0}\left(1-\frac{\omega-\bm{n}\cdot\bm{\omega}_{\ast}-\bm{n}\cdot\bm{\omega}_{E}}{\omega-\bm{n}\cdot\bm{\Omega}}\right),$
(95)

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,

$\nabla^{2}\phi=\sum_{s}-\frac{e_{s}n_{s}}{\epsilon_{0}},$ (96)

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

$n_{s}=\int d^{3}vf_{s}.$ (97)

The perturbed electrostatic potential is calculated using the perturbed charge
density

$\delta n_{s}=\int d^{3}v\delta f_{s}.$ (98)

To enforce quasineutrality, we take the sum of the total charge density to be
0 and require that

$\lambda_{D}\ll\left|\frac{\phi}{\nabla\phi}\right|,$ (99)

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. We thus obtain

$\sum_{s}\int d^{3}ve_{s}\delta f_{s}=0.$ (100)

Since $\phi=\phi(\bm{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,
resulting in

$\sum_{s}\int d^{3}rd^{3}ve_{s}\phi^{\ast}\delta f_{s}=0.$ (101)

By multiplying by the electrostatic potential and integrating, we have recast
the differential equation via a weak formulation using the variational
method(Samain, 1970; Garbet et al., 1990; Garbet, 2001; Nguyen et al., 2008).
Instead of solving for the exact function $\phi$ or $\delta f$ that satisfies
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, 1991).

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

$\sum_{s}\sum_{\bm{n},\bm{n}^{\prime}}\int
d^{3}rd^{3}v\frac{e_{s}^{2}\phi_{\bm{n}}\phi_{\bm{n}^{\prime}}^{\ast}}{T_{s}}f_{0,s}\left(1-\frac{\omega-\bm{n}\cdot\bm{\omega}_{\ast}-\bm{n}\cdot\bm{\omega}_{E}}{\omega-\bm{n}\cdot\bm{\Omega}}\right)e^{i(\bm{n}-\bm{n}^{\prime})\cdot\bm{\alpha}}e^{-i(\omega-(\omega^{\prime})^{\ast})t}=0.$
(102)

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

$\sum_{s}\sum_{\bm{n},\bm{n}^{\prime}}\int d^{3}\alpha
d^{3}J\frac{e_{s}^{2}\phi_{\bm{n}}\phi_{\bm{n}^{\prime}}^{\ast}}{m_{s}T_{s}}f_{0,s}\left(1-\frac{\omega-\bm{n}\cdot\bm{\omega}_{\ast}-\bm{n}\cdot\bm{\omega}_{E}}{\omega-\bm{n}\cdot\bm{\Omega}}\right)e^{i(\bm{n}-\bm{n}^{\prime})\cdot\bm{\alpha}}e^{-i(\omega-(\omega^{\prime})^{\ast})t}=0.$
(103)

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

$\sum_{s}\sum_{n_{1},n_{2}}\int
d^{3}J\frac{e_{s}^{2}\left|\phi_{\bm{n}}\right|^{2}}{m_{s}T_{s}}f_{0,s}\left(1-\frac{\omega-\bm{n}\cdot\bm{\omega}_{\ast}-\bm{n}\cdot\bm{\omega}_{E}}{\omega-\bm{n}\cdot\bm{\Omega}}\right)=0.$
(104)

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 $\bm{\alpha}$. As such, we integrate over the
action angles again and transform back to conventional variables, leading to

$\sum_{s}\sum_{n_{1},n_{2}}\int
d^{3}rd^{3}v\frac{e_{s}^{2}\left|\phi_{\bm{n}}\right|^{2}}{T_{s}}f_{0,s}\left(1-\frac{\omega-\bm{n}\cdot\bm{\omega}_{\ast}-\bm{n}\cdot\bm{\omega}_{E}}{\omega-\bm{n}\cdot\bm{\Omega}}\right)=0.$
(105)

Even though the integrand is a function of only $\bm{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 $\left(\bm{r},\bm{v}\right).$

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

$\sum_{s}\mathcal{L}_{0,s}-\mathcal{L}_{\text{passing},s}-\mathcal{L}_{\text{trapped},s}=0.$
(106)

Here, $\mathcal{L}_{0}$ is the portion of the integral that is simply
multiplied by $1$, which we call the adiabatic part.
$\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
$\left|\phi_{\bm{n}}\right|^{2}$. This requires a 3-dimensional integral over
$d^{3}\alpha$. 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 6-dimensional
integral, we can utilize a number of symmetries, transformations, and
approximations to simplify the form down to at most 2-dimensional 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 $\left|\phi_{\bm{n}}\right|^{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,

$\phi(r,\theta,\varphi)=\sum_{m,n}\phi_{m,n}(r-r_{0})e^{i(m\theta+n\varphi)}.$
(107)

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_{n,m}$ all have
nearly identical radial envelopes where each radial profile is centered on
their corresponding reference flux surface (Connor et al., 1979). These flux
surfaces are all rational flux surfaces since $m$ and $n$ are integers.
Meanwhile, the general ballooning representation of $\phi$ is

$\phi(r,\theta,\varphi)=\sum_{p}\sum_{n}\hat{\phi}_{n}(\theta+2p\pi,\theta_{0})e^{in\left(\varphi-q(r)\left(\theta-\theta_{0}+2p\pi\right)\right)},$
(108)

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 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. However, this shift is typically
on the order of $10^{-1}$ in relevant cases (Candy, 2009). Thus, for the rest
of the derivation we 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$,
leading to

$\phi(r,\theta,\varphi)=\sum_{n}\hat{\phi}_{n}(\theta)e^{in\left(\varphi-q(r)\theta\right)}.$
(109)

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

$\displaystyle\bm{B}\cdot\nabla\psi$ $\displaystyle=0,$ (110)
$\displaystyle\bm{B}\cdot\nabla\left(\varphi-q\theta\right)$
$\displaystyle=0,$ (111) $\displaystyle\bm{B}\cdot\nabla\theta$
$\displaystyle\neq 0,$ (112)

indicating that with the above set of variables that $\theta$ indicates the
location on any given field line. Because the magnetic curvature is
unfavorable 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$, leading to

$\int_{-\pi}^{\pi}\frac{d\varphi}{2\pi}\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\phi(r,\theta,\varphi)e^{-in\varphi-
im\theta}=\phi_{m,n}(r-r_{0})=\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\hat{\phi}_{n}(\theta)e^{-i\theta(nq(r)+m)}.$
(113)

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

$nq(r)+m=nq_{0}+\frac{r-r_{0}}{d}+m=\frac{r-r_{0}}{d}=\frac{x}{d}+\frac{\tilde{r}}{d},$
(114)

where $q_{0}=q(r_{0})$, the radial difference between different rational flux
surfaces is defined as,

$\frac{1}{nd}=\left.\frac{dq}{dr}\right|_{r=r_{0}},$ (115)

and $x$ is defined as

$x=\bar{r}-r_{0},$ (116)

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

$\phi_{m,n}(r-r_{0})\approx\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\hat{\phi}_{n}(\theta)e^{-i\frac{\theta\left(r-r_{0}\right)}{d}}.$
(117)

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 $-\pi$ to $\pi$. The
result is

$\phi_{m,n}(r-r_{0})\approx\int_{-\infty}^{\infty}\frac{d\theta}{2\pi}\hat{\phi}_{n}(\theta)e^{-i\frac{\theta\left(r-r_{0}\right)}{d}}.$
(118)

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

$\hat{\phi}_{n}(\theta)=\int^{\infty}_{-\infty}\frac{dr}{\left|d\right|}\phi_{m,n}(r-r_{0})e^{i\frac{\theta(r-r_{0})}{d}}.$
(119)

We are now in a position to integrate over the action angles to fully
calculate $\phi_{\bm{n}}$. The procedure to integrate over $\alpha_{1}$ has
already been discussed in Section 2, where we find that

$\int_{-\pi}^{\pi}\frac{d\alpha_{1}}{2\pi}g(\bm{r})e^{-in_{1}\alpha_{1}}=\left(-i\right)^{n_{1}}\left(J_{n_{1}}\left(k_{\perp}\rho\right)\cdot
g\right)\left(\bm{R}_{G}\right).$ (120)

We therefore only need to discuss in detail the integrations over $\alpha_{2}$
and $\alpha_{3}$ while treating all variables within the guiding center
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

$\displaystyle r$ $\displaystyle=\bar{r}+\delta_{b}\cos(\alpha_{2}),$ (121)
$\displaystyle\theta$ $\displaystyle=\theta_{b}\sin(\alpha_{2}),$ (122)
$\displaystyle\varphi$
$\displaystyle=\alpha_{3}+q(\bar{r})\theta_{b}\sin(\alpha_{2})+\tilde{\varphi}.$
(123)

Here, we define the banana width $\delta_{b}$ as

$\delta_{b}=\frac{q\rho}{\sqrt{\epsilon}}.$ (124)

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,

$\int^{\pi}_{-\pi}\frac{d\alpha_{2}}{2\pi}\phi(\bm{r})e^{-in_{2}\alpha_{2}}=\int^{\pi}_{-\pi}\frac{d\alpha_{2}}{2\pi}\int\frac{d^{3}k}{(2\pi)^{3}}\tilde{\phi}(\bm{k})e^{-i\bm{k}\cdot\bm{r}-in_{2}\alpha_{2}}.$
(125)

We then proceed in fashion similar to the gyro-average derivation in Section 2
by noting that
$\bm{k}\cdot{\bm{r}}=\bm{k}\cdot\bar{\bm{r}}+k_{r}\delta_{b}\cos(\alpha_{2})$.
The result is

$\int^{\pi}_{-\pi}\frac{d\alpha_{2}}{2\pi}\phi(\bm{r})e^{-in_{2}\alpha_{2}}=\int\frac{d^{3}k}{(2\pi)^{3}}(-i)^{n_{2}}J_{n_{2}}(k_{r}\delta_{b})\tilde{\phi}(\bm{k})e^{-i\bm{k}\cdot\bar{\bm{r}}}.$
(126)

In essence, we obtain a bounce average over the banana width (Depret et al.,
2000). 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 Biglari & Chen
(1986), 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 integrating over $\alpha_{3}$, for now leaving the Bessel
functions aside and evaluating the position at $\bm{r}=\bar{\bm{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 structure, we are free to keep
only one of the 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

$\int_{-\pi}^{\pi}\frac{d\alpha_{3}}{2\pi}\phi(\bar{\bm{r}})e^{-in_{3}\alpha_{3}}=\sum_{n}\int_{-\pi}^{\pi}\frac{d\alpha_{3}}{2\pi}\phi_{m_{0},n}(\bar{r}-r_{0})e^{i\left(n\alpha_{3}-n_{3}\alpha_{3}\right)}=\phi_{m_{0},n_{3}}(\bar{r}-r_{0}).$
(127)

To compensate for choosing only one poloidal harmonic, we must extend the
radial limits of integration to $-\infty<r<\infty$. 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

$\left|\phi_{\bm{n}}\right|^{2}=\left|J_{n_{1}}(k_{\perp}\rho)J_{n_{2}}(\delta_{b}k_{r})\cdot\phi_{m_{0},n_{3}}\right|^{2}(\bar{r}-r_{0})=\left|J_{n_{1}}(k_{\perp}\rho)J_{n_{2}}(\delta_{b}k_{r})\cdot\phi_{m_{0},n_{3}}\right|^{2}(x),$
(128)

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

###  4.2 Passing

We now calculate $\phi_{\bm{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

$\int_{-\pi}^{\pi}\frac{d\alpha_{3}}{2\pi}\phi(\bm{r})e^{-in_{3}\alpha_{3}}=\sum_{n}\int_{-\pi}^{\pi}\frac{d\alpha_{3}}{2\pi}\hat{\phi}_{n}(\theta(\alpha_{2}))e^{in(\alpha_{3}+\tilde{\varphi}-q(r)\alpha_{2}+(q(\bar{r})-q(r))\tilde{\theta})-in_{3}\alpha_{3}}.$
(129)

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}$:

$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}.$
(130)

Carrying out the integral then gives us

$\int_{-\pi}^{\pi}\frac{d\alpha_{3}}{2\pi}\phi(\bm{r})e^{-in_{3}\alpha_{3}}=\hat{\phi}_{n_{3}}(\theta(\alpha_{2}))e^{in_{3}(\tilde{\varphi}-q(r)\alpha_{2}-\frac{\tilde{r}}{n_{3}d}\tilde{\theta})}.$
(131)

We now multiply by $e^{-in_{2}\alpha_{2}}$ and integrate with respect to
$\alpha_{2}$. The eikonal can be simplified if we only keep
$n_{2}=m_{0}=-nq_{0}$,

$i\left(n_{3}\tilde{\varphi}-n_{3}q(r)\alpha_{2}-\frac{\tilde{r}}{d}\tilde{\theta}-m_{0}\alpha_{2}\right)=i\left(n_{3}\tilde{\varphi}-\frac{x}{d}\alpha_{2}-\frac{\tilde{r}}{d}\left(\alpha_{2}+\tilde{\theta})\right)\right)=i\left(n_{3}\tilde{\varphi}-\frac{x}{d}\alpha_{2}-\frac{\tilde{r}}{d}\theta(\alpha_{2})\right),$
(132)

where we have used

$-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}$
(133)

and the expression $\alpha_{2}=\theta-\tilde{\theta}$. We then obtain

$\int_{-\pi}^{\pi}\frac{d\alpha_{2}}{2\pi}\hat{\phi}_{n_{3}}(\theta(\alpha_{2}))e^{i\left(n_{3}\tilde{\varphi}-\frac{x}{d}\alpha_{2}-\frac{\tilde{r}}{d}\theta(\alpha_{2})\right)}\approx\int_{-\infty}^{\infty}\frac{d\alpha_{2}}{2\pi}\hat{\phi}_{n_{3}}(\theta(\alpha_{2}))e^{i\left(n_{3}\tilde{\varphi}-\frac{\tilde{r}}{d}\theta(\alpha_{2})\right)}e^{-i\frac{x}{d}\alpha_{2}},$
(134)

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

$\phi_{\bm{n}}=\phi_{n_{1},m_{0},n_{3}}=\left(J_{n_{1}}(k_{\perp}\rho)\cdot\mathcal{F}^{-1}\left(\hat{\phi}_{n_{3}}(\theta(\alpha_{2}))e^{i\left(n_{3}\tilde{\varphi}-\frac{\tilde{r}}{d}\theta(\alpha_{2})\right)}\right)\right)(x),$
(135)

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:

$\phi_{m_{0},n}(x)\sim\phi_{0}e^{-\frac{\left(x-x_{0}\right)^{2}}{2w^{2}}}.$
(136)

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
centered about $x=0$. 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 Section 9.

To obtain expressions for $w$ and $x_{0}$, we move into the high-frequency
fluid limit. The original derivation can be found in Cottier et al. (2014) and
an extensive, revised derivation can be found in Citrin et al. (2017); here,
we shall only discuss the basic principle. We consider the dispersion relation

$D\left(\omega\right)=\sum_{s}\int
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).$
(137)

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}\left(k_{\perp}\rho_{s}\right)J_{0}\left(k_{r}\delta_{b,s}\right)$
for trapped particles and $J_{0,s}=J_{0}\left(k_{\perp}\rho_{s}\right)$ for
passing particles. Meanwhile, we define the parallel wave number as
$k_{\parallel}=\left(k_{\theta}sx\right)/\left(qR_{0}\right)$. The local drift
frequency $\omega_{d,s}$ is

$\omega_{d,s}=\left\\{\begin{aligned}
&\omega_{d0,s}\xi\Biggl{(}-1+\frac{2E\left(\kappa\right)}{K\left(\kappa\right)}+4s\left(\kappa^{2}-1+\frac{E\left(\kappa\right)}{K\left(\kappa\right)}\right)&&\\\
&\qquad\qquad\quad-\frac{4\alpha}{3}\left(1-\kappa^{2}-\left(1-2\kappa^{2}\right)\frac{E\left(\kappa\right)}{K\left(\kappa\right)}\right)\Biggr{)}&&\text{
if trapped},\\\ &\omega_{d0,s}\xi\left(2-\lambda
b\right)\left(\cos\left(\theta\right)+\left(s\theta-\alpha\sin\left(\theta\right)\right)\sin\left(\theta\right)\right)&&\text{
if passing}.\end{aligned}\right.$ (138)

The passing form of the drift frequency is due to the radial structure of the
eigenfunction as covered in Section 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 $\mkern 1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu=\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

$\begin{split}D\left(\omega\right)=\Biggl{[}\frac{n_{e}}{T_{e}}&\left(1-\left<\left(1-\frac{n\omega_{\ast,e}}{\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu}\right)\left(1+\frac{n\omega_{d,e}}{\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu}+\frac{n^{2}\omega_{d,e}^{2}}{\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu^{2}}\right)\right>_{t}\right)\\\
&+\sum_{i}\frac{n_{i}Z_{i}}{T_{i}}\Biggl{(}1-\left<\left(1-\frac{n\omega_{\ast,i}}{\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu}\right)\left(1+\frac{n\omega_{d,i}}{\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu}+\frac{n^{2}\omega_{d,i}^{2}}{\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu^{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>_{t}\\\
&\qquad-\left<\left(1-\frac{n\omega_{\ast,i}}{\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu}\right)\left(1+\frac{n\omega_{d,i}}{\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu}+\frac{k_{\parallel}v_{\parallel}}{\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu}+\left(\frac{n\omega_{d,i}}{\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu}+\frac{k_{\parallel}v_{\parallel}}{\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu}\right)^{2}\right)\right.\\\
&\qquad\qquad\times\left.\left(1-\frac{k_{\theta}^{2}\rho_{i}^{2}}{2}-\frac{k_{r}^{2}\rho_{i}^{2}}{2}\right)\right>_{p}\Biggr{)}\Biggr{]}\phi_{m_{0},n}(x),\end{split}$
(139)

where $Z_{i}$ is the proton number of the ion species. We define the averages
over velocity space as

$\displaystyle\left<g\left(\bm{v}\right)\right>_{t}$
$\displaystyle=\int_{\text{trapped}}d^{3}vf_{0}g\left(\bm{v}\right),$ (140)
$\displaystyle\left<g\left(\bm{v}\right)\right>_{p}$
$\displaystyle=\int_{\text{passing}}d^{3}vf_{0}g\left(\bm{v}\right).$ (141)

Do to the $\theta$ dependent terms, this is a differential equation. We
approximate the differential operators on $\phi$ in the limit of small mode
shift $x_{0}$, leading to

$\frac{\partial^{2}\phi}{\partial
x^{2}}=\left(\frac{x^{2}}{w^{4}}-\frac{2x_{0}x}{w^{4}}-\frac{1}{w^{2}}\right)\phi.$
(142)

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,

$\omega^{3}D\left(\omega\right)=D_{0}\left(\omega\right)+D_{1}\left(\omega\right)x+D_{2}\left(\omega\right)x^{2}=0.$
(143)

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

$D_{0}\left(\omega_{0}\right)=0.$ (144)

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

$\displaystyle D_{1}\left(\omega_{0}\right)=0,$ (145) $\displaystyle
D_{2}\left(\omega_{0}\right)=0.$ (146)

We do not cite the full solution here and direct the reader to Citrin et al.
(2017) for a complete derivation. Now that we have characterized
$\phi_{\bm{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

$\mathcal{L}_{0}=\sum_{n_{1},n_{2}}\int
d^{3}rd^{3}v\frac{e^{2}\left|\phi_{\bm{n}}\right|^{2}}{T}f_{0}.$ (147)

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},\bm{J})$ such that

$\phi_{n}=\int_{-\pi}^{\pi}\frac{d\alpha_{3}}{2\pi}\phi(\bm{r})e^{-in\alpha_{3}}=\sum_{n_{1},n_{2}}\phi_{\bm{n}}e^{i(n_{1}\alpha_{1}+n_{2}\alpha_{2})}.$
(148)

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

$\int_{-\pi}^{\pi}\frac{d\alpha_{1}}{2\pi}\int_{-\pi}^{\pi}\frac{d\alpha_{2}}{2\pi}\left|\phi_{n}\right|^{2}=\sum_{n_{1},n_{2}}\left|\phi_{\bm{n}}\right|^{2}.$
(149)

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

$\mathcal{L}_{0}=\sum_{n_{1},n_{2}}\int d^{3}\alpha
d^{3}J\frac{e^{2}\left|\phi_{\bm{n}}\right|^{2}}{mT}f_{0}=\sum_{n_{1},n_{2}}\int
4\pi^{2}d\alpha_{3}d^{3}J\frac{e^{2}\left|\phi_{\bm{n}}\right|^{2}}{mT}f_{0}=\int
4\pi^{2}d\alpha_{3}d^{3}J\frac{e^{2}\left|\phi_{n}\right|^{2}}{mT}f_{0}.$
(150)

This then simplifies to

$\displaystyle\mathcal{L}_{0}=\int
4\pi^{2}d\alpha_{3}d^{3}J\frac{e^{2}}{mT}f_{0}\int_{-\pi}^{\pi}\frac{d\alpha_{1}}{2\pi}\int_{-\pi}^{\pi}\frac{d\alpha_{2}}{2\pi}\left|\phi_{n}\right|^{2}=\int
d^{3}rd^{3}v\frac{e^{2}\left|\phi_{n}\right|^{2}}{T}f_{0}.$ (151)

The velocity space integration is straightforward,

$\mathcal{L}_{0}=\int
d^{3}rd^{3}v\frac{e^{2}\left|\phi_{n}\right|^{2}}{T}f_{0}=\int
d^{3}r\frac{e^{2}n_{0}\left|\phi_{n}\right|^{2}}{T},$ (152)

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

$d^{3}r=rR_{0}\left(1+\epsilon\cos\left(\theta\right)\right)drd\theta
d\varphi.$ (153)

We proceed to calculating $\phi_{n}$ using the poloidal harmonic expansion as
detailed in Section 4,

$\phi(r,\theta,\varphi)=\sum_{m,n}\phi_{m,n}(r-r_{0})e^{i(m\theta+n\varphi)}.$
(154)

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

$\phi_{n}=\sum_{m}\phi_{m,n}(r-r_{0})e^{i(m\theta+nq(\bar{r})\tilde{\theta}+n\tilde{\varphi})}.$
(155)

As before, we only keep the poloidal harmonic corresponding to $m_{0}=-nq_{0}$
and expand the limits of integration for $r$ to compensate. The result is

$\left|\phi_{n}\right|^{2}=\left|\phi_{m_{0},n}(r-r_{0})\right|^{2}.$ (156)

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

$\mathcal{L}_{0}=\int_{-\infty}^{\infty}\left(2\pi R_{0}\right)2\pi
rdr\frac{e^{2}n_{0}}{T}\left|\phi_{m_{0},n}(r-r_{0})\right|^{2}\approx\int
dxR_{0}r_{0}\left(2\pi\right)^{2}\frac{e^{2}n_{0}}{T}\left|\phi_{m_{0},n}(x)\right|^{2}.$
(157)

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

$\mathcal{L}_{0}=4\pi^{2}R_{0}r_{0}\left|\phi_{0}\right|^{2}\left|w\right|^{2}\exp\left(\frac{\Imag\left(x_{0}\right)^{2}}{\Real\left(w^{2}\right)}\right)\sqrt{\frac{\pi}{\Real\left(w^{2}\right)}}.$
(158)

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

$\mathcal{L}_{\text{trapped}}=\sum_{n_{1},n_{2}}\int
d^{3}rd^{3}v\frac{e^{2}}{T}f_{0}\left(\frac{\bm{n}\cdot\bm{\omega}_{\ast}+\bm{n}\cdot\bm{\omega}_{E}-\omega}{\bm{n}\cdot\bm{\Omega}-\omega}\right)|\phi_{\bm{n}}|^{2}.$
(159)

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 use once again use toroidal coordinates,

$d^{3}r=R_{0}\bar{r}\left(1+\epsilon\cos\left(\theta\right)\right)d\bar{r}d\theta
d\varphi.$ (160)

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

$d^{3}v=\sum_{\epsilon_{\parallel}}v^{2}\frac{b}{2\sqrt{1-\lambda
b}}dvd\lambda dv_{\phi},$ (161)

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

$d^{3}rd^{3}v=\sum_{\epsilon_{\parallel}}2\pi^{2}R_{0}\bar{r}\left(1+\epsilon\cos\left(\theta\right)\right)d\bar{r}d\theta
v^{2}\frac{b}{\sqrt{1-\lambda b}}dvd\lambda.$ (162)

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

$\frac{1-\epsilon}{1+\epsilon}\leq\lambda\leq\frac{1+\epsilon\cos\left(\theta\right)}{1+\epsilon}.$
(163)

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:

$\int_{0}^{2\pi}d\theta\int_{\frac{1-\epsilon}{1+\epsilon}}^{\frac{1+\epsilon\cos\left(\theta\right)}{1+\epsilon}}d\lambda
f(\theta,\lambda)=\int_{\frac{1-\epsilon}{1+\epsilon}}^{1}d\lambda\int_{-\theta_{b}}^{\theta_{b}}d\theta
f(\theta,\lambda).$ (164)

We recall that the definition of a bounce average is

$\left<G(\epsilon_{\parallel},\theta)\right>=\frac{\int_{-\theta_{b}}^{\theta_{b}}d\theta\frac{G\left(\epsilon_{\parallel},\theta\right)+G\left(-\epsilon_{\parallel},\theta\right)}{\sqrt{1-\lambda
b}}}{\int_{-\theta_{b}}^{\theta_{b}}\frac{2d\theta}{\sqrt{1-\lambda
b}}}=\sum_{\epsilon_{\parallel}}\frac{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}{2\pi}\int_{-\theta_{b}}^{\theta_{b}}d\theta\frac{G\left(\epsilon_{\parallel},\theta\right)}{\sqrt{1-\lambda
b}}.$ (165)

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. Since the electron and ion rotation velocity is quite small
compared to the sound speed in tokamak plasmas, expanding to second-order in
the Mach number will be sufficient. The distribution function then simplifies
to

$f_{0}\approx n_{0}\left(\frac{m}{2\pi
T}\right)^{3/2}e^{-\xi}\left(1+\frac{2v_{\parallel}U_{\parallel}}{v_{T}^{2}}+\frac{U_{\parallel}^{2}}{v_{T}^{2}}\left(\frac{2v_{\parallel}^{2}}{v_{T}^{2}}-1\right)\right).$
(166)

Moreover, because
$\left|\Omega_{1}\right|,\left|\Omega_{2}\right|\gg\left|\omega\right|$, 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

$\begin{split}\mathcal{L}_{\text{trapped}}=\int d\bar{r}&d\lambda
d\xi\left(2\pi\right)^{2}\frac{n_{0}e^{2}}{T}\frac{R_{0}\bar{r}}{\sqrt{\pi}}\frac{\sqrt{\xi}e^{-\xi}}{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu_{2}}\\\
&\times\frac{\left<\left(\bm{n}\cdot\bm{\omega}_{\ast}+\bm{n}\cdot\bm{\omega}_{E}-\omega\right)\left(1+\frac{2v_{\parallel}U_{\parallel}}{v_{T}^{2}}+\frac{U_{\parallel}^{2}}{v_{T}^{2}}\left(\frac{2v_{\parallel}^{2}}{v_{T}^{2}}-1\right)\right)\right>}{n\Omega_{3}-\omega}\left|\phi_{0,0,n}\right|^{2}.\end{split}$
(167)

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
$\bm{J}$. Because $n_{1}=n_{2}=0$, we only keep the partial derivative with
respect to $J_{3}$. Knowing that $J_{3}=J_{3}\left(\bar{r}\right)$, we perform
a change in variables from $J_{3}$ to $\bar{r}$ and find that

$\frac{\partial g}{\partial J_{3}}=\frac{\partial
g}{\partial\bar{r}}\frac{d\bar{r}}{dJ_{3}}\approx\frac{-R_{0}\omega_{d0}}{T}\frac{\partial
g}{\partial\bar{r}},$ (168)

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

$\displaystyle A_{n}$ $\displaystyle=-\frac{R_{0}}{n}\frac{dn}{d\bar{r}},$
(169) $\displaystyle A_{T}$
$\displaystyle=-\frac{R_{0}}{T}\frac{dT}{d\bar{r}},$ (170) $\displaystyle
A_{U}$ $\displaystyle=-\frac{R_{0}}{v_{T}}\frac{dU_{\parallel}}{d\bar{r}}.$
(171)

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

$\displaystyle\left<v_{\parallel}\right>$
$\displaystyle=\left<v_{\parallel}^{3}\right>=0,$ (172)
$\displaystyle\left<v_{\parallel}^{2}\right>$
$\displaystyle=\frac{4E\epsilon}{m}\frac{\left(E\left(\kappa\right)-\left(1-\kappa^{2}\right)K\left(\kappa\right)\right)}{K\left(\kappa\right)}=v_{T}^{2}\xi
H\left(\kappa\right),$ (173)

where we define

$H\left(\kappa\right)=\frac{2\epsilon\left(E\left(\kappa\right)-\left(1-\kappa^{2}\right)K\left(\kappa\right)\right)}{K\left(\kappa\right)}.$
(174)

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

$\mkern 1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu=\omega-\bm{n}\cdot\bm{\omega}_{E}\approx\omega-n\omega_{E}\approx\omega-n\omega_{E0}-n\omega_{E}^{\prime}x,$
(175)

where, as discussed earlier, we Taylor expand $\omega_{E}$ about $x=0$ and
$\partial_{r}\omega_{E}=\omega_{E}^{\prime}$ is related to the radial electric
shear. Rather than including $x$ fully, we instead approximate the term by
averaging it over the Gaussian eigenfunctions:

$\left<x\right>_{r}=\frac{\int_{-\infty}^{\infty}dx\left|\phi_{m_{0},n}\left(x\right)\right|^{2}x}{\int_{-\infty}^{\infty}dx\left|\phi_{m_{0},n}\left(x\right)\right|^{2}}=\Real\left(x_{0}\right)+\frac{\Imag\left(x_{0}\right)\Imag\left(w^{2}\right)}{\Real\left(w^{2}\right)}.$
(176)

We then obtain

$\mkern 1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu\approx\omega-n\omega_{E0}-n\omega_{E}^{\prime}\left<x\right>_{r}.$ (177)

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

$\left<\left(\bm{n}\cdot\bm{\omega}_{\ast}+\bm{n}\cdot\bm{\omega}_{E}-\omega\right)\left(1+\frac{2v_{\parallel}U_{\parallel}}{v_{T}^{2}}+\frac{U_{\parallel}^{2}}{v_{T}^{2}}\left(\frac{2v_{\parallel}^{2}}{v_{T}^{2}}-1\right)\right)\right>=n\omega_{d0}\left(\mathcal{A}_{t}+\mathcal{B}_{t}\xi+\mathcal{C}_{t}\xi^{2}\right),$
(178)

where

$\displaystyle\mathcal{A}_{t}$
$\displaystyle=\left(1-\frac{U_{\parallel}^{2}}{v_{T}^{2}}\right)\left(A_{n}-\frac{3}{2}A_{T}-z^{2}F_{d}\left(\kappa\right)\right)-\frac{U_{\parallel}}{v_{T}}\left(2A_{U}-\frac{U_{\parallel}}{v_{T}}A_{T}\right),$
(179) $\displaystyle\mathcal{B}_{t}$
$\displaystyle=\left(1-\frac{U_{\parallel}^{2}}{v_{T}^{2}}\right)A_{T}+4A_{U}\frac{U_{\parallel}}{v_{T}}H\left(\kappa\right)+\frac{U_{\parallel}^{2}}{v_{T}^{2}}H\left(\kappa\right)\left(2A_{n}-7A_{T}-2z^{2}F_{d}\left(\kappa\right)\right)$
$\displaystyle\approx\left(1-\frac{U_{\parallel}^{2}}{v_{T}^{2}}\right)A_{T}+4A_{U}\frac{U_{\parallel}}{v_{T}}H\left(\kappa\right),$
(180) $\displaystyle\mathcal{C}_{t}$
$\displaystyle=2A_{T}\frac{U_{\parallel}^{2}}{v_{T}^{2}}H\left(\kappa\right)\approx
0.$ (181)

Here, we have defined

$z^{2}=\frac{\mkern 1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu}{n\omega_{d0}F_{d}\left(\kappa\right)}.$ (182)

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

Substituting the above into the integrand, we obtain

$\mathcal{L}_{\text{trapped}}=\int d\bar{r}d\lambda
d\xi\left(2\pi\right)^{2}\frac{n_{0}e^{2}}{T}\frac{R_{0}\bar{r}}{\sqrt{\pi}}\frac{\sqrt{\xi}e^{-\xi}}{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}\left|\phi_{0,0,n}\right|^{2}\frac{\mathcal{A}_{t}+\mathcal{B}_{t}\xi}{F_{d}\left(\kappa\right)\left(\xi-z^{2}\right)}.$
(183)

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

$\mathcal{L}_{\text{trapped}}=\int d\bar{r}d\lambda
d\xi\left(2\pi\right)^{2}\frac{n_{0}e^{2}}{T}\frac{R_{0}r_{0}}{\sqrt{\pi}}\frac{\sqrt{\xi}e^{-\xi}}{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}\left|J_{0}\left(k_{\perp}\rho\right)J_{0}\left(k_{r}\delta_{b}\right)\cdot\phi_{m_{0},n}\right|^{2}\frac{\mathcal{A}_{t}+\mathcal{B}_{t}\xi}{F_{d}\left(\kappa\right)\left(\xi-z^{2}\right)}.$
(184)

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 Parseval’s theorem to integrate over
$k_{r}$,

$\int_{-\infty}^{\infty}dxf(x)g(x)^{\ast}=\int_{-\infty}^{\infty}\frac{dk_{r}}{2\pi}\hat{f}(k_{r})\hat{g}(k_{r})^{\ast}.$
(185)

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,

$\displaystyle\rho$ $\displaystyle=\frac{v_{\perp}}{\Omega_{1}},$ (186)
$\displaystyle\delta_{b}$ $\displaystyle\approx\frac{q}{\sqrt{\epsilon}}\rho.$
(187)

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

$\frac{\int d^{3}vJ_{0}\left(k_{\perp}\rho\right)^{2}f_{0}}{\int
d^{3}vf_{0}}=e^{-\frac{k_{\perp}^{2}\rho_{\text{th}}^{2}}{2}}I_{0}\left(\frac{k_{\perp}^{2}\rho_{\text{th}}^{2}}{2}\right)=\Gamma_{0}\left(k_{\perp}\rho_{\text{th}}\right),$
(188)

where $I_{0}$ is a modified Bessel function of the first kind and the
characteristic thermal gyroradius $\rho_{\text{th}}$ is defined as

$\rho_{\text{th}}=\frac{\sqrt{2T/m}}{\Omega_{1}}.$ (189)

Similarly, for the average over the banana orbit we obtain

$\frac{\int d^{3}vJ_{0}\left(k_{\perp}\rho\right)^{2}f_{0}}{\int
d^{3}vf_{0}}=e^{-\frac{k_{\perp}^{2}\delta_{b,\text{th}}^{2}}{2}}I_{0}\left(\frac{k_{\perp}^{2}\delta_{b,\text{th}}^{2}}{2}\right)=\Gamma_{0}\left(k_{\perp}\delta_{b,\text{th}}\right),$
(190)

where the thermal banana width is

$\delta_{b,\text{th}}=\frac{q}{\sqrt{\epsilon}}\rho_{\text{th}}.$ (191)

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

$k_{\perp}^{2}=k_{r}^{2}+k_{\theta}^{2}=k_{r}^{2}+\frac{n^{2}q_{0}^{2}}{r_{0}^{2}},$
(192)

where we have evaluate $k_{\theta}$ at $r_{0}$. That
$k_{\theta}^{2}=n^{2}q^{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

$\begin{split}\mathcal{L}_{\text{trapped}}=\int d\lambda
d\xi\left(2\pi\right)^{2}&\frac{n_{0}e^{2}}{T}\frac{R_{0}r_{0}}{\sqrt{\pi}}\frac{\sqrt{\xi}e^{-\xi}}{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}\frac{\mathcal{A}_{t}+\mathcal{B}_{t}\xi}{F_{d}\left(\kappa\right)\left(\xi-z^{2}\right)}\\\
&\times\int_{-\infty}^{\infty}\frac{dk_{r}}{2\pi}\Gamma_{0}\left(k_{\perp}\rho_{\text{th}}\right)\Gamma_{0}\left(k_{r}\delta_{b,\text{th}}\right)\left|d\hat{\phi}_{n}\left(k_{r}d\right)\right|^{2},\end{split}$
(193)

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

$\hat{\phi}_{n}\left(k_{r}d\right)=\sqrt{2\pi}w\phi_{0}e^{-\frac{k_{r}^{2}w^{2}}{2}}e^{ik_{r}x_{0}}.$
(194)

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

$\int_{0}^{\infty}d\xi\frac{\sqrt{\xi}}{\sqrt{\pi}}\frac{\mathcal{A}_{t}+\mathcal{B}\xi}{\xi-z^{2}}e^{-\xi}=\int_{-\infty}^{\infty}du\frac{u^{2}}{\sqrt{\pi}}\frac{\mathcal{A}_{t}+\mathcal{B}_{t}u^{2}}{\left(u+z\right)\left(u-z\right)}e^{-u^{2}},$
(195)

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

$\int_{-\infty}^{\infty}du\frac{u^{2}}{\sqrt{\pi}}\frac{\mathcal{A}_{t}+\mathcal{B}_{t}u^{2}}{\left(u+z\right)\left(u-z\right)}e^{-u^{2}}=\mathcal{A}_{t}G_{2}\left(z,-z\right)+\mathcal{B}_{t}G_{4}\left(z,-z\right)=\frac{\mathcal{A}_{t}Z_{2}\left(z\right)+\mathcal{B}_{t}Z_{4}\left(z\right)}{z},$
(196)

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

$\frac{d\lambda}{\mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}=4f_{t}K\left(\kappa\right)\kappa d\kappa,$ (197)

where we utilize the transformation

$\lambda\approx 1-2\epsilon\kappa^{2}$ (198)

and define the flux surface averaged trapped particle fraction

$f_{t}=\frac{2\sqrt{2\epsilon}}{\pi}.$ (199)

Thus, the trapped functional simplifies to

$\mathcal{L}_{\text{trapped}}=\left(2\pi\right)^{3}\frac{e^{2}n_{0}}{T}r_{0}R_{0}f_{t}\left<\mathcal{I}_{t}\right>_{\xi,\kappa}\left<\Gamma_{0}\left(k_{\perp}\rho_{\text{th}}\right)\Gamma_{0}\left(k_{r}\delta_{b,\text{th}}\right)\left|d\hat{\phi}\left(k_{r}d\right)\right|^{2}\right>_{k_{r}},$
(200)

where

$\left<\mathcal{I}_{t}\right>_{\xi,\kappa}=\frac{2}{\pi}\int_{0}^{1}d\kappa\frac{K\left(\kappa\right)\kappa}{zF_{d}\left(\kappa\right)}\left(\mathcal{A}_{t}Z_{2}\left(z\right)+\mathcal{B}_{t}Z_{4}\left(z\right)\right)$
(201)

and

$\left<\Gamma_{0}\left(k_{\perp}\rho_{\text{th}}\right)\Gamma_{0}\left(k_{r}\delta_{b,\text{th}}\right)\left|d\hat{\phi}\left(k_{r}d\right)\right|^{2}\right>_{k_{r}}=\int_{-\infty}^{\infty}\frac{dk_{r}}{2\pi}\Gamma_{0}\left(k_{\perp}\rho_{\text{th}}\right)\Gamma_{0}\left(k_{r}\delta_{b,\text{th}}\right)\left|d\hat{\phi}\left(k_{r}d\right)\right|^{2}.$
(202)

The remaining integrals are to be calculated numerically, where we note that
$z$ is a function of both $\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern 1.5mu$ and $\kappa$.
Thus, the trapped functional is the product of two separate $1$-dimensional
integrals, one of which is $\omega$ independent; we therefore characterize the
trapped functional as a $1$-dimensional 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

$\mathcal{L}_{\text{passing}}=\sum_{n_{1},n_{2}}\int
d^{3}rd^{3}v\frac{e^{2}}{T}f_{0}\left(\frac{\bm{n}\cdot\bm{\omega}_{\ast}+\bm{n}\cdot\bm{\omega}_{E}-\omega}{\bm{n}\cdot\bm{\Omega}-\omega}\right)|\phi_{\bm{n}}|^{2}.$
(203)

We reuse many of the same arguments in Section 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

$\left<G\left(\epsilon_{\parallel},\theta\right)\right>=\frac{\int_{-\pi}^{\pi}d\theta\frac{G\left(\epsilon_{\parallel},\theta\right)}{\sqrt{1-\lambda
b}}}{\int_{-\pi}^{\pi}\frac{d\theta}{\sqrt{1-\lambda b}}}=\frac{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}{2\pi}\int_{-\pi}^{\pi}d\theta\frac{G\left(\epsilon_{\parallel},\theta\right)}{\sqrt{1-\lambda
b}}.$ (204)

We note here that the bounce angle $\theta_{b}$ is set to $\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

$0\leq\lambda\leq\frac{1-\epsilon}{1+\epsilon}.$ (205)

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
$\left|\Omega_{1}\right|\gg\left|\omega\right|$. As discussed in Section 4,
$n_{2}$ refers to the poloidal harmonic. We keep only $n_{2}=m_{0}$ and use
the approximation that

$m_{0}+nq\left(\bar{r}\right)\approx\frac{x}{d}.$ (206)

In the resonant denominator we then obtain

$n\Omega_{3}+m_{0}\Omega_{2}-\omega\approx n\omega_{d0}\xi
F+\frac{x}{d}\Omega_{2}+n\omega_{E0}+n\omega_{E}^{\prime}x-\omega,$ (207)

where we also expand $\omega_{E}$ about $x=0$. The passing functional is then

$\begin{split}\mathcal{L}_{\text{passing}}=\sum_{\epsilon_{\parallel}}\int
d\bar{r}&d\lambda
d\xi\left(2\pi\right)^{2}\frac{n_{0}e^{2}}{T}\frac{R_{0}\bar{r}}{\sqrt{\pi}}\frac{\sqrt{\xi}e^{-\xi}}{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu_{2}}\\\
&\times\frac{\left<\left(\bm{n}\cdot\bm{\omega}_{\ast}+n\omega_{E0}+n\omega_{E}^{\prime}x-\omega\right)\left(1+\frac{2v_{\parallel}U_{\parallel}}{v_{T}^{2}}+\frac{U_{\parallel}^{2}}{v_{T}^{2}}\left(\frac{2v_{\parallel}^{2}}{v_{T}^{2}}-1\right)\right)\right>}{n\omega_{d0}\xi
F+\frac{x}{d}\Omega_{2}+n\omega_{E0}+n\omega_{E}^{\prime}x-\omega}\left|\phi_{0,m_{0},n}\right|^{2}.\end{split}$
(208)

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 Parseval’s theorem to integrate over $\alpha_{2}$ instead
of $\bar{r}$,

$\int_{-\infty}^{\infty}d\bar{r}f(\bar{r})g(\bar{r})^{\ast}=\int_{-\infty}^{\infty}\frac{d\alpha_{2}}{2\pi}\left|d\right|\hat{f}(\alpha_{2})\hat{g}(\alpha_{2})^{\ast}.$
(209)

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

$\int_{-\infty}^{\infty}d\bar{r}\frac{1}{a\frac{x}{d}-b}\left|\phi_{0,m_{0},n}\right|^{2}=\int_{-\infty}^{\infty}\frac{d\alpha_{2}}{2\pi}\left|d\right|\mathcal{F}\left(\frac{\phi_{0,m_{0},n}}{a\frac{x}{d}-b}\right)\mathcal{F}\left(\phi_{0,m_{0},n}\right)^{\ast}.$
(210)

We calculated in Section 4 that

$\mathcal{F}\left(\phi_{0,m_{0},n}\right)=J_{0}(k_{\perp}(\alpha_{2})\rho)\hat{\phi}_{n}(\theta(\alpha_{2}))e^{in\tilde{\varphi}(\alpha_{2})-i\frac{\tilde{r}(\alpha_{2})}{d}\theta(\alpha_{2})}.$
(211)

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

$k_{\perp}\left(\alpha_{2}\right)^{2}=\frac{\theta\left(\alpha_{2}\right)^{2}}{d^{2}}+\frac{n^{2}q_{0}^{2}}{r_{0}^{2}}.$
(212)

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

$\mathcal{F}\left(\frac{\phi_{0,m_{0},n}}{a\frac{x}{d}-b}\right)=\frac{1}{2\pi}\mathcal{F}\left(\phi_{0,m_{0},n}\right)\ast\mathcal{F}\left(\frac{1}{a\frac{x}{d}-b}\right).$
(213)

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

$\begin{split}\mathcal{F}\left(\phi_{0,m_{0},n}\frac{1}{1+a\frac{x}{d}-b}\right)=\int_{-\infty}^{\infty}&d\alpha_{2}^{\prime}\frac{i}{\left|a\right|}\Theta\left(\frac{\alpha_{2}-\alpha_{2}^{\prime}}{a}\right)J_{0}(k_{\perp}(\alpha_{2}^{\prime})\rho)\hat{\phi}_{n}(\theta(\alpha_{2}^{\prime}))\\\
&\times
e^{\frac{-ib}{a}\left(\alpha_{2}^{\prime}-\alpha_{2}\right)}e^{in\tilde{\varphi}(\alpha_{2}^{\prime})-i\frac{\tilde{r}(\alpha_{2})}{d}\theta(\alpha_{2}^{\prime})},\end{split}$
(214)

where $\Theta$ is the Heaviside step function. Here, we have assumed that
$\Imag\left(b\right)>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

$\begin{split}I_{p,r}=\int
d\bar{r}\frac{1}{a\frac{x}{d}-b}\left|\phi_{0,m_{0},n}\right|^{2}=\int\frac{d\alpha_{2}d\alpha_{2}^{\prime}}{2\pi}\frac{i|d|}{\left|a\right|}&\Theta\left(\frac{\alpha_{2}-\alpha_{2}^{\prime}}{a}\right)J_{0}^{\ast}(k_{\perp}(\alpha_{2})\rho)J_{0}(k_{\perp}(\alpha_{2}^{\prime})\rho)\\\
&\times\hat{\phi}_{n}(\theta(\alpha_{2}))^{\ast}\hat{\phi}_{n}(\theta(\alpha_{2}^{\prime}))e^{\Lambda(\alpha_{2})-\Lambda(\alpha_{2}^{\prime})},\end{split}$
(215)

where

$\Lambda(\alpha_{2})=i\left(\frac{b}{a}\alpha_{2}-n\tilde{\varphi}\left(\alpha_{2}\right)+\frac{\tilde{r}\left(\alpha_{2}\right)}{d}\theta\left(\alpha_{2}\right)\right).$
(216)

We then substitute in

$\displaystyle a$ $\displaystyle=\Omega_{2}+nd\omega_{E}^{\prime},$ (217)
$\displaystyle b$ $\displaystyle=\omega-n\omega_{d0}\xi F-n\omega_{E0}$ (218)

and rewrite the eikonal term to obtain

$\Lambda(\alpha_{2})=i\left(\frac{\omega-n\omega_{d0}\xi
F-n\omega_{E0}}{\Omega_{2}+nd\omega_{E}^{\prime}}\alpha_{2}-n\tilde{\varphi}\left(\alpha_{2}\right)+\frac{\tilde{r}\left(\alpha_{2}\right)}{d}\theta\left(\alpha_{2}\right)\right).$
(219)

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

Before proceeding, we must recognize that integrating over $\alpha_{2}$ and
$\alpha_{2}^{\prime}$ 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^{\prime}$. First,
we introduce new variables,

$\displaystyle\theta_{+}$ $\displaystyle=\frac{\theta+\theta^{\prime}}{2},$
(220) $\displaystyle\theta_{-}$ $\displaystyle=\theta-\theta^{\prime}.$ (221)

We next Taylor expand the terms in the exponential around $\theta_{+}$ to find
that

$\Lambda(\alpha_{2})-\Lambda(\alpha_{2}^{\prime})\approx\theta_{-}\frac{d\alpha_{2}}{d\theta}\left(\theta=\theta_{+}\right)\Lambda^{\prime}\left(\theta=\theta_{+}\right),$
(222)

where $\Lambda^{\prime}$ denotes the derivative of $\Lambda$ 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 Section 2, we find

$\frac{d}{d\alpha_{2}}\left(n\tilde{\varphi}-\frac{\tilde{r}}{d}\theta\right)=\frac{n}{\Omega_{2}}\left(\bm{v}_{D}\cdot\nabla\varphi-q\bm{v}_{D}\cdot\nabla\theta-\theta\frac{dq}{d\psi}\bm{v}_{D}\cdot\nabla\psi-\Omega_{d}\right).$
(223)

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

$\begin{split}\mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{d}&=\bm{v}_{D}\cdot\nabla\varphi-q\bm{v}_{D}\cdot\nabla\theta-\theta\frac{dq}{d\psi}\bm{v}_{D}\cdot\nabla\psi\\\
&=\omega_{d0}\xi\left(2-\lambda
b\right)\left(\cos\left(\theta\right)+\left(s\theta-\alpha\sin\left(\theta\right)\sin^{2}\left(\theta\right)\right)\right)+\frac{qE_{r}}{rB}.\end{split}$
(224)

Although somewhat similar to the magnetic drift frequency proper, there are
two key differences. Firstly, the magnetic shear term is different and
proportional to $\theta\sin\left(\theta\right)$. Secondly, 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

$\Lambda\left(\alpha_{2}\right)-\Lambda\left(\alpha_{2}^{\prime}\right)\approx-i\theta_{-}\left(\frac{n\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{d}\left(\theta_{+}\right)-\omega}{\Omega_{2}+nd\omega_{E}^{\prime}}\right),$
(225)

where

$n\mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu_{d}\approx
n\omega_{d0}\xi\left(2-\lambda
b\right)\left(\cos\left(\theta\right)+\left(s\theta-\alpha\sin\left(\theta\right)\sin^{2}\left(\theta\right)\right)\right)+n\omega_{E0}.$
(226)

We next change the variables of integration from
$\alpha_{2},\alpha_{2}^{\prime}$ to $\theta,\theta^{\prime}$ using

$\frac{d\theta}{d\alpha_{2}}=\frac{\sqrt{1-\lambda b}}{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu_{2}}.$ (227)

The integral then becomes

$\begin{split}I_{p,r}=\int\frac{d\theta
d\theta^{\prime}}{2\pi}&\frac{i\left|d\right|}{\left|\Omega_{2}+nd\omega_{E}^{\prime}\right|}\Theta\left(\frac{\theta_{-}}{\Omega_{2}+nd\omega_{E}^{\prime}}\right)J_{0}\hat{\phi}_{n}(\theta)^{\ast}J_{0}^{\prime}\hat{\phi}_{n}(\theta^{\prime})\\\
&\times\frac{\mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}{\sqrt{1-\lambda b(\theta)}}\frac{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}{\sqrt{1-\lambda b(\theta^{\prime})}}e^{-i\theta_{-}\frac{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}{\sqrt{1-\lambda b(\theta_{+})}}\left(\frac{n\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{d}(\theta_{+})-\omega}{\Omega_{2}+nd\omega_{E}^{\prime}}\right)},\end{split}$
(228)

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

$\begin{split}I_{p,r}=\int\frac{d\theta
d\theta^{\prime}dxdx^{\prime}}{2\pi}&\frac{i\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}^{2}}{\left|\Omega_{2}+nd\omega_{E}^{\prime}\right|\left|d\right|}\Theta\left(\frac{\theta_{-}}{\Omega_{2}+nd\omega_{E}^{\prime}}\right)J_{0}\phi^{\ast}_{m_{0},n}(x)J_{0}^{\prime}\phi_{m_{0},n}(x^{\prime})\\\
&\times\frac{e^{-i\theta_{-}\frac{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}{\sqrt{1-\lambda b(\theta_{+})}}\left(\frac{n\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{d}(\theta_{+})-\omega}{\Omega_{2}+nd\omega_{E}^{\prime}}\right)}e^{-i\frac{\theta
x}{d}}e^{i\frac{\theta^{\prime}x^{\prime}}{d}}}{\sqrt{1-\lambda
b\left(\theta\right)}\sqrt{1-\lambda
b\left(\theta^{\prime}\right)}}.\end{split}$ (229)

We then make the following substitutions

$\displaystyle x_{+}$ $\displaystyle=\frac{x+x^{\prime}}{2},$ (230)
$\displaystyle x_{-}$ $\displaystyle=x-x^{\prime},$ (231) $\displaystyle
k_{+}$ $\displaystyle=\frac{\theta_{+}}{\left|d\right|},$ (232) $\displaystyle
k_{-}$ $\displaystyle=\frac{\theta_{-}}{\left|d\right|},$ (233) $\displaystyle
d\theta d\theta^{\prime}dxdx^{\prime}$
$\displaystyle=dk_{+}dk_{-}dx_{+}dx_{-}\left|d\right|^{2},$ (234)
$\displaystyle\frac{\theta^{\prime}x^{\prime}}{d}-\frac{\theta x}{d}$
$\displaystyle=\left(-k_{-}x_{+}-k_{+}x_{-}\right)\frac{\left|d\right|}{d},$
(235)

to obtain

$\begin{split}I_{p,r}=\int\frac{dk_{+}dk_{-}dx_{+}dx_{-}}{2\pi}&\frac{i\left|d\right|\left(\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}\right)^{2}}{|\Omega_{2}+nd\omega_{E}^{\prime}|}\Theta\left(\frac{k_{-}}{\Omega_{2}+nd\omega_{E}^{\prime}}\right)J_{0}\phi^{\ast}_{m_{0},n}\left(x_{+}+\frac{x_{-}}{2}\right)J_{0}^{\prime}\phi_{m_{0},n}\left(x_{+}-\frac{x_{-}}{2}\right)\\\
&\times\frac{\exp\left(\frac{-ik_{-}\left|d\right|\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}{\sqrt{1-\lambda b(k_{+}\left|d\right|)}}\left(\frac{n\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{d}(k_{+}\left|d\right|)-\omega}{\Omega_{2}+nd\omega_{E}^{\prime}}\right)\right)e^{-i\left(k_{-}x_{+}+k_{+}x_{-}\right)\frac{\left|d\right|}{d}}}{\sqrt{1-\lambda
b\left(k_{+}\left|d\right|+\frac{k_{-}\left|d\right|}{2}\right)}\sqrt{1-\lambda
b\left(k_{+}\left|d\right|-\frac{k_{-}\left|d\right|}{2}\right)}}.\end{split}$
(236)

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

$\begin{split}\int_{a}^{b}dsg\left(s\right)e^{ics}&=\int_{a}^{b}ds\frac{1}{ik}\left(\frac{d}{ds}\left(g\left(s\right)e^{ics}\right)-\frac{dg}{ds}e^{ics}\right)\\\
&=\left.\frac{g\left(s\right)e^{ics}}{ic}\right|_{a}^{b}-\int_{a}^{b}ds\frac{1}{\left(ic\right)^{2}}\left(\frac{d}{ds}\left(\frac{dg}{ds}e^{ics}\right)-\frac{d^{2}g}{ds^{2}}e^{ics}\right).\\\
&\approx\sum_{m=0}^{N}\frac{\left(-1\right)^{m}}{\left(ic\right)^{m+1}}\left[\frac{d^{m}g}{ds^{m}}\left(b\right)-\frac{d^{m}g}{ds^{m}}\left(a\right)\right].\end{split}$
(237)

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 $\Imag\left(\omega\right)>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

$\begin{split}\int_{-\infty}^{\infty}dk_{-}\frac{i\left|d\right|}{\left|\Omega_{2}+nd\omega_{E}^{\prime}\right|}&g\left(k_{-}\right)e^{-ik_{-}\frac{\left|d\right|}{d}}\exp\left(\frac{-ik_{-}\left|d\right|\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}{\sqrt{1-\lambda b(k_{+}\left|d\right|)}}\left(\frac{n\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{d}(k_{+}\left|d\right|)-\omega}{\Omega_{2}+nd\omega_{E}^{\prime}}\right)\right)\\\
&\approx\frac{g(0)}{\frac{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}{\sqrt{1-\lambda b\left(k_{+}\left|d\right|\right)}}\left(n\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{d}\left(k_{+}\left|d\right|\right)-\omega\right)+\Omega_{2}\frac{x_{+}}{d}+n\omega_{E}^{\prime}x_{+}}.\end{split}$
(238)

For convenience, we next replace all instances of $k_{+}\left|d\right|$ with
$k_{+}d$; this is allowed since $\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu_{d}$ and $b$ are
even functions and the bounds of integration are symmetric. We obtain

$I_{p,r}\approx\int\frac{dk_{+}dx_{+}dx_{-}}{2\pi}\frac{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}e^{-ik_{+}x_{-}}}{\sqrt{1-\lambda
b\left(k_{+}d\right)}}\frac{J_{0}\left(\rho
k_{\perp}\right)^{2}\phi_{m0,n}\left(x_{+}+\frac{x_{-}}{2}\right)^{\ast}\phi_{m0,n}\left(x_{+}-\frac{x_{-}}{2}\right)}{n\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{d}\left(k_{+}d\right)+\left(\frac{\Omega_{2}}{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}\frac{x_{+}}{d}+n\frac{\omega_{E}^{\prime}}{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}x_{+}\right)\sqrt{1-\lambda b\left(k_{+}d\right)}-\omega},$ (239)

where

$k_{\perp}^{2}=k_{+}^{2}+\frac{n^{2}q_{0}^{2}}{r_{0}^{2}}.$ (240)

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

$\frac{\int d^{3}vJ_{0}\left(k_{\perp}\rho\right)^{2}f_{0}}{\int
d^{3}vf_{0}}=\Gamma_{0}\left(k_{\perp}\rho_{\text{th}}\right).$ (241)

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

$\begin{split}I_{p,r}=\int\frac{dx_{+}dk_{+}}{\sqrt{\pi}}\frac{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}{\sqrt{1-\lambda
b\left(k_{+}d\right)}}\frac{\Gamma_{0}\left(k_{\perp}\rho_{\text{th}}\right)\sqrt{\frac{\Imag\left(w^{2}\right)^{2}}{\Real\left(w^{2}\right)}+\Real\left(w^{2}\right)}e^{-\rho_{\ast}^{2}}e^{-k_{\ast}^{2}}\exp\left(\frac{\Imag\left(x_{0}\right)^{2}}{\Real\left(w^{2}\right)}\right)}{n\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{d}\left(k_{+}d\right)+\left(\frac{\Omega_{2}}{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}\frac{x_{+}}{d}+n\frac{\omega_{E}^{\prime}}{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}x_{+}\right)\sqrt{1-\lambda
b\left(k_{+}d\right)}-\omega},\end{split}$ (242)

where

$\displaystyle\rho_{\ast}$
$\displaystyle=\frac{x_{+}+k_{+}\Imag\left(w^{2}\right)-\Real\left(x_{0}\right)}{\sqrt{\Real\left(w^{2}\right)}},$
(243) $\displaystyle k_{\ast}$
$\displaystyle=\frac{k_{+}\Real\left(w^{2}\right)+\Imag\left(x_{0}\right)}{{\sqrt{\Real\left(w^{2}\right)}}}.$
(244)

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

$\begin{split}I_{p,r}\approx\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{d\rho_{\ast}dk_{\ast}}{\sqrt{\pi}}\frac{\Gamma_{0}\left(k_{\perp}\rho_{\text{th}}\right)\sqrt{\frac{\Imag\left(w^{2}\right)^{2}}{\Real\left(w^{2}\right)}+\Real\left(w^{2}\right)}e^{-\rho_{\ast}^{2}}e^{-k_{\ast}^{2}}\exp\left(\frac{\Imag\left(x_{0}\right)^{2}}{\Real\left(w^{2}\right)}\right)}{n\omega_{d0}\xi
F_{p}\left(k_{+}d\right)+\epsilon_{\parallel}\sqrt{\xi}\frac{x_{+}}{d}\frac{\sqrt{2T/m}}{qR_{0}}-\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern 1.5mu},\end{split}$ (245)

where

$\mkern 1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu=\omega-n\omega_{E0}-n\omega_{E}^{\prime}x_{+},$ (246)

and

$F_{p}\left(k_{+}d\right)=\frac{4}{3}\left(\cos\left(k_{+}d\right)+\left(sk_{+}d-\alpha\sin\left(k_{+}d\right)\right)\sin\left(k_{+}d\right)\right).$
(247)

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.

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 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

$I_{p,E}=\sum_{\epsilon_{\parallel}}\int d\lambda
d\xi\left(2\pi\right)^{2}\frac{n_{0}e^{2}}{T}\frac{R_{0}\bar{r}}{\sqrt{\pi}}\frac{\sqrt{\xi}e^{-\xi}}{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}\frac{\left<\left(\bm{n}\cdot\bm{\omega}_{\ast}-\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu\right)\left(1+\frac{2v_{\parallel}U_{\parallel}}{v_{T}^{2}}+\frac{U_{\parallel}^{2}}{v_{T}^{2}}\left(\frac{2v_{\parallel}^{2}}{v_{T}^{2}}-1\right)\right)\right>}{n\omega_{d0}\xi
F_{p}\left(k_{+}d\right)+\epsilon_{\parallel}\sqrt{\xi}\frac{x_{+}}{d}\frac{\sqrt{2T/m}}{qR_{0}}-\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern 1.5mu}.$ (248)

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-dimensional rather than 3-dimensional. It is of potential interest to study
the impact this approximation has; one could calculate a more exact (albeit
slower) integral to quantify the exact impact this has on the resulting
solutions and flux calculations.

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

$\int_{0}^{\frac{1-\epsilon}{1+\epsilon}}\frac{d\lambda}{\mkern
1.5mu\overline{\mkern-1.5mu2\Omega\mkern-1.5mu}\mkern 1.5mu_{2}}=f_{p},$ (249)

where $f_{p}=1-f_{t}$ is the flux surface averaged passing particle fraction.
We then compute

$\int_{0}^{\frac{1-\epsilon}{1+\epsilon}}\frac{d\lambda}{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}\left<v_{\parallel}^{m}\right>=2f_{p}v_{T}^{m}\epsilon_{\parallel}^{m}\xi^{m/2}\lambda_{m},$
(250)

where we define $\lambda_{m}$ as

$\lambda_{m}=\frac{\int_{0}^{\frac{1-\epsilon}{1+\epsilon}}d\lambda\int_{-\pi}^{\pi}\frac{d\theta\left(\sqrt{1-\lambda
b}\right)^{m}}{\sqrt{1-\lambda
b}}}{\int_{0}^{\frac{1-\epsilon}{1+\epsilon}}d\lambda\int_{-\pi}^{\pi}\frac{d\theta}{\sqrt{1-\lambda
b}}}.$ (251)

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

$\begin{split}\int_{0}^{\frac{1-\epsilon}{1+\epsilon}}\frac{d\lambda}{\mkern
1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{2}}&\left<\left(\bm{n}\cdot\bm{\omega}_{\ast}-\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu\right)\left(1+\frac{2v_{\parallel}U_{\parallel}}{v_{T}^{2}}+\frac{U_{\parallel}^{2}}{v_{T}^{2}}\left(\frac{2v_{\parallel}^{2}}{v_{T}^{2}}-1\right)\right)\right>=\\\
&\qquad\qquad
2f_{p}n\omega_{d0}\left(\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}\right),\end{split}$
(252)

where we define the terms

$\displaystyle\mathcal{A}_{p}$
$\displaystyle=\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(2A_{U}-\frac{U_{\parallel}}{v_{T}}A_{T}\right),$
(253) $\displaystyle\mathcal{B}_{p}$
$\displaystyle=\left(A_{U}\left(2-6\frac{U_{\parallel}^{2}}{v_{T}^{2}}\right)+U_{\parallel}\left(2A_{n}-5A_{T}-2z^{2}F_{p}\right)\right)\lambda_{1},$
(254) $\displaystyle\mathcal{C}_{p}$
$\displaystyle=\left(1-\frac{U_{\parallel}^{2}}{v_{T}^{2}}\right)A_{T}+4A_{U}\frac{U_{\parallel}}{v_{T}}\lambda_{2}+\frac{U_{\parallel}^{2}}{v_{T}^{2}}\lambda_{2}\left(2A_{n}-7A_{T}-2z^{2}F_{p}\right),$
(255) $\displaystyle\mathcal{D}_{p}$
$\displaystyle=2A_{T}\frac{U_{\parallel}}{v_{T}}\lambda_{1}+4A_{U}\frac{U_{\parallel}^{2}}{v_{T}^{2}}\lambda_{3},$
(256) $\displaystyle\mathcal{E}_{p}$
$\displaystyle=2A_{T}\frac{U_{\parallel}^{2}}{v_{T}^{2}}\lambda_{2},$ (257)

and where

$z^{2}=\frac{\mkern 1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu}{n\omega_{d0}F_{p}\left(k_{+}d\right)}.$ (258)

Thus, the integral simplifies to

$I_{p,E}=\sum_{\epsilon_{\parallel}}2f_{p}\int_{0}^{\infty}d\xi\left(2\pi\right)^{2}\frac{n_{0}e^{2}}{T}\frac{R_{0}r_{0}}{\sqrt{\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}}{F_{p}\left(k_{+}d\right)\left(\xi+\epsilon_{\parallel}\sqrt{\xi}\frac{x_{+}}{d}\frac{v_{T}}{qR_{0}F_{p}\left(k_{+}d\right)}-z^{2}\right)}.$
(259)

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

$\sum_{\epsilon_{\parallel}}\int_{0}^{\infty}d\xi\sqrt{\xi}g\left(\epsilon_{\parallel}\sqrt{\xi}\right)=\sum_{\epsilon_{\parallel}}\int_{0}^{\infty}du2u^{2}g\left(\epsilon_{\parallel}u\right)=\int_{-\infty}^{\infty}du2u^{2}g\left(u\right).$
(260)

The integral then becomes

$\begin{split}I_{p,E}=4f_{p}\int_{-\infty}^{\infty}du\left(2\pi\right)^{2}\frac{n_{0}e^{2}}{T}\frac{R_{0}r_{0}}{\sqrt{\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}\left(k_{+}d\right)\left(u^{2}+u\frac{x_{+}}{d}\frac{v_{T}}{qR_{0}F_{p}\left(k_{+}d\right)}-z^{2}\right)}.\end{split}$
(261)

To simplify this integral further, we rewrite the denominator as

$u^{2}+u\frac{x_{+}}{d}\frac{v_{T}}{qR_{0}F_{p}\left(k_{+}d\right)}-z^{2}=\left(u-z_{+}\right)\left(u-z_{-}\right),$
(262)

where

$z_{\pm}=-\frac{1}{2}\frac{x_{+}}{d}\frac{v_{T}}{qR_{0}F_{p}\left(k_{+}d\right)}\pm\sqrt{\left(\frac{1}{2}\frac{x_{+}}{d}\frac{v_{T}}{qR_{0}F_{p}\left(k_{+}d\right)}\right)^{2}+z^{2}}.$
(263)

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

$\begin{split}I_{p,E}=4f_{p}\left(2\pi\right)^{2}\frac{n_{0}e^{2}}{T}\frac{R_{0}r_{0}}{F_{p}\left(k_{+}d\right)}\left(\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}\right),\end{split}$
(264)

where the associated Fried and Conte integrals
$G_{n}=G_{n}\left(z_{+},z_{-}\right)$ are evaluated at $z_{+}$ and $z_{-}$.
Thus, the passing functional simplifies to

$\begin{split}\mathcal{L}_{\text{passing}}=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}&\frac{d\rho_{\ast}dk_{\ast}}{\sqrt{\pi}}\left(2\pi\right)^{3}\frac{e^{2}n_{0}}{T}r_{0}R_{0}f_{p}\left<\mathcal{I}_{p}\right>_{\xi,\lambda}\Gamma_{0}\left(k_{\perp}\rho_{\text{th}}\right)\\\
&\times\sqrt{\frac{\Imag\left(w^{2}\right)^{2}}{\Real\left(w^{2}\right)}+\Real\left(w^{2}\right)}e^{-\rho_{\ast}^{2}}e^{-k_{\ast}^{2}}\exp\left(\frac{\Imag\left(x_{0}\right)^{2}}{\Real\left(w^{2}\right)}\right),\end{split}$
(265)

where

$\left<\mathcal{I}_{p}\right>_{\xi,\lambda}=\frac{2}{\pi
F_{p}\left(k_{+}\right)}\left(\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}\right).$
(266)

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 $1$\- and $2$-dimensional
integrals, respectively. With the dispersion relation in hand, we can proceed
to applying 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 turn-over time. The ratio of these two
quantities is known as the Kubo number (Kubo, 1963; Krommes, 2002). 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, leading to a justification for the quasilinear approach. These
characteristic times have been calculated and compared for both ETG and ITG-
TEM turbulence (Lin et al., 2008; Casati et al., 2009; Citrin et al., 2012).
For these general cases, the Kubo number is less than unity and well developed
turbulence for tokamak plasma parameters manifests random walk processes.
Moreover, it has been found that quasilinear models are successful in
reproducing experimental results such as temperature profiles within $15\%$
rms error (Kinsey et al., 2008).

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

$\frac{\partial f}{\partial t}+\dot{\bm{\alpha}}\cdot\frac{\partial
f}{\partial\bm{\alpha}}+\dot{\bm{J}}\cdot\frac{\partial f}{\partial\bm{J}}=0.$
(267)

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 $\left<f\right>_{t}=f_{0}$ and the linear
response averages to zero. We define the time average as

$\left<g\left(t\right)\right>_{t}=\frac{1}{T}\int_{-T/2}^{T/2}g\left(t+t^{\prime}\right)dt^{\prime},$
(268)

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

$\left<\frac{\partial f}{\partial t}+\dot{\bm{\alpha}}\cdot\frac{\partial
f}{\partial\bm{\alpha}}+\dot{\bm{J}}\cdot\frac{\partial
f}{\partial\bm{J}}\right>_{t}\approx\frac{\partial f_{0}}{\partial
t}+\left<\left\\{\Real(\delta f),\Real(e\phi)\right\\}\right>_{t}=0.$ (269)

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

$\left\\{\Real(\delta
f),\Real(e\phi)\right\\}=\frac{\partial}{\partial\bm{\alpha}}\cdot\left(\Real(\delta
f)\frac{\partial\Real(e\phi)}{\partial\bm{J}}\right)-\frac{\partial}{\partial\bm{J}}\cdot\left(\Real(\delta
f)\frac{\partial\Real(e\phi)}{\partial\bm{\alpha}}\right).$ (270)

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

$\left<\Real\left(\bm{A}e^{-i\omega t}\right)\cdot\Real\left(\bm{B}e^{-i\omega
t}\right)\right>_{t}=\frac{1}{2}\Real\left(\bm{A}\cdot\bm{B}^{\ast}\right).$
(271)

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

$\frac{\partial}{\partial\bm{\alpha}}\left<\Real\left(\delta
f\right)\Real\left(e\phi\right)\right>_{t}=\bm{0}.$ (272)

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

$\frac{\partial f_{0}}{\partial
t}+\frac{\partial}{\partial\bm{J}}\cdot\bm{\Gamma}_{Q}=0,$ (273)

where we define the quasilinear flux $\bm{\Gamma}_{Q}$ as

$\begin{split}\bm{\Gamma}_{Q}=\frac{1}{2}\Real\left(\sum_{\bm{n}}i\bm{n}f_{\bm{n}}e\phi_{\bm{n}}^{\ast}\right)=-\frac{1}{2}\Imag\left(\sum_{\bm{n}}\bm{n}\frac{e^{2}\left|\phi_{\bm{n}}\right|^{2}}{T}f_{0}\left(1-\frac{\omega-\bm{n}\cdot\bm{\omega}_{\ast}-\bm{n}\cdot\bm{\omega}_{E}}{\omega-\bm{n}\cdot\bm{\Omega}}\right)\right).\end{split}$
(274)

Here, $f_{\bm{n}}$ and $\phi_{\bm{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

$\displaystyle\frac{\partial\left<n\right>}{\partial t}+\frac{d\Gamma}{dr}$
$\displaystyle=0,$ (275)
$\displaystyle\frac{\partial\left<mnRU_{\parallel}\right>}{\partial
t}+\frac{d\Pi}{dr}$ $\displaystyle=0,$ (276)
$\displaystyle\frac{3}{2}\frac{\partial\left<p\right>}{\partial
t}+\frac{dQ}{dr}$ $\displaystyle=0,$ (277)

where $\Gamma$, $\Pi$, and $Q_{E}$ are the particle, toroidal momentum, and
energy fluxes defined as

$\displaystyle\Gamma$ $\displaystyle=\frac{1}{4\pi^{2}d}\int
d^{3}vd^{3}r\frac{1}{2}\Imag\left(\sum_{\bm{n}}\frac{nq}{rB}f_{\bm{n}}\phi_{\bm{n}}^{\ast}\right),$
(278) $\displaystyle\Pi$ $\displaystyle=\frac{1}{4\pi^{2}d}\int
d^{3}vd^{3}r\frac{mRv_{\parallel}}{2}\Imag\left(\sum_{\bm{n}}\frac{nq}{rB}f_{\bm{n}}\phi_{\bm{n}}^{\ast}\right),$
(279) $\displaystyle Q_{E}$ $\displaystyle=\frac{1}{4\pi^{2}d}\int
d^{3}vd^{3}r\frac{m\left(v^{2}-U_{\parallel}^{2}\right)}{4}\Imag\left(\sum_{\bm{n}}\frac{nq}{rB}f_{\bm{n}}\phi_{\bm{n}}^{\ast}\right).$
(280)

Here, we can see that the integrations to calculate the particle, toroidal