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

$\delta\bm{v}_{E\times
B}\cdot\hat{\bm{r}}=\hat{\bm{r}}\cdot-\frac{\nabla\delta\phi\times\bm{B}}{B^{2}}\approx\sum_{\bm{n}}\frac{ik_{\theta}\phi_{\bm{n}}}{B}=\sum_{\bm{n}}\frac{inq}{rB}\phi_{\bm{n}},$
(281)

where we have again used the convention that
$k_{\theta}\to\left(i/r\right)\partial_{\theta}$. We then find that

$\left<\Real\left(\delta\bm{v}_{E\times
B}\cdot\hat{\bm{r}}\right)\Real\left(\delta
f\right)\right>_{t}=\frac{1}{2}\Real\left(\sum_{\bm{n}}f_{\bm{n}}\left(\frac{inq}{rB}\phi_{\bm{n}}\right)^{\ast}\right)=\frac{1}{2}\Imag\left(\sum_{\bm{n}}f_{\bm{n}}\frac{nq}{rB}\phi_{\bm{n}}^{\ast}\right).$
(282)

This lets us write the fluxes as

$\displaystyle\Gamma$
$\displaystyle=\left<\delta\left(n\right)\delta\bm{v}_{E\times
B}\cdot\hat{\bm{r}}\right>_{t,r},$ (283) $\displaystyle\Pi$
$\displaystyle=\left<\delta\left(mnRU_{\parallel}\right)\delta\bm{v}_{E\times
B}\cdot\hat{\bm{r}}\right>_{t,r},$ (284) $\displaystyle Q_{E}$
$\displaystyle=\left<\delta\left(P\right)\delta\bm{v}_{E\times
B}\cdot\hat{\bm{r}}\right>_{t,r}.$ (285)

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

$\displaystyle\delta n$ $\displaystyle=\int d^{3}v\delta f,$ (286)
$\displaystyle\delta\left(nmRU_{\parallel}\right)$ $\displaystyle=\int
d^{3}vmv_{\parallel}\delta f,$ (287) $\displaystyle\delta P$
$\displaystyle=\int
d^{3}v\frac{1}{2}m\left(v^{2}-U_{\parallel}^{2}\right)\delta f,$ (288)

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

$\displaystyle\delta
P=\delta\left(nT\right)=-\frac{1}{2}mU_{\parallel}^{2}\delta n+\int
d^{3}v\frac{1}{2}mv^{2}\approx\int d^{3}v\frac{1}{2}mv^{2}.$ (289)

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

$Q=Q_{E}-\frac{3T}{2}\Gamma.$ (290)

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

##  9 Saturation Rule

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

$D=\left.\frac{\gamma_{n}}{\left<k_{\perp}^{2}\right>}\right|_{\text{max}},$
(291)

where we compute the value of $\gamma_{n}$ such that the quantity
$\gamma_{n}/\left<k_{\perp}^{2}\right>$ is at its maximum over the linear
spectrum for a given mode. Meanwhile, we average $k_{\perp}^{2}$ over the
electrostatic mode. We enforce this mixing length estimate for our various
flux calculations by approximating the underlying process as a random walk
(Bourdelle et al., 2007). For instance, we mandate that the particle flux for
a given species must be

$\Gamma_{s}=\sum_{n}C_{\text{NL}}\frac{S_{n}}{R_{0}n_{0s}}\left.\frac{\gamma_{n}}{\left<k_{\perp}^{2}\right>}\right|_{\text{max}}\frac{k_{\theta}}{k_{\theta,\text{max}}}L_{s,n,0},$
(292)

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

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

$C_{\text{NL}}=\left\\{\begin{aligned}
&271/s_{\text{fac}}\qquad\qquad\quad&&\text{ if }k_{\theta}\rho_{s}<2\text{
(ITG)},\\\ &122f_{\text{multi-scale}}/s_{\text{fac}}&&\text{ if
}k_{\theta}\rho_{s}>2\text{ (ETG)}.\end{aligned}\right.$ (293)

Here, we have also introduced an ad hoc factor $s_{\text{fac}}$ for the case
of low magnetic shear (Citrin et al., 2012),

$s_{\text{fac}}=\left\\{\begin{aligned}
&2.5\left(1-\left|s\right|\right)\qquad&&\text{ if }\left|s\right|<0.6,\\\
&1&&\text{ if }\left|s\right|>0.6,\end{aligned}\right.$ (294)

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

$f_{\text{multi-
scale}}=\frac{1}{1+\exp\left(-\frac{1}{5}\left(\frac{\gamma_{\text{ETG,
max}}}{\gamma_{\text{ITG, max}}}-\sqrt{\frac{m_{i}}{m_{e}}}\right)\right)},$
(295)

where $m_{e}$ and $m_{i}$ are the masses of the electron and main ion
respectively. Here, the sigmoid guarantees a smooth transition from a strongly
driven ion-scale mode regime and a strongly driven electron-scale mode regime,
since it has been observed that ETG turbulence is suppressed when the ion-
scale instability dominates.

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

$\left<k_{\perp}^{2}\right>=\left\\{\begin{aligned}
&k_{\theta}^{2}+\left(k_{r-\text{NL}}+k_{r-\text{shear}}\right)^{2}\qquad&&\text{
if }k_{\theta}\rho_{s}<2\text{ (ITG)},\\\ &2k_{\theta}^{2}&&\text{ if
}k_{\theta}\rho_{s}>2\text{ (ETG)}.\end{aligned}\right.$ (296)

The shear contribution can be calculated analytically as

$k_{r-\text{shear}}=k_{\theta}\left|s\right|\sqrt{\left<\theta^{2}\right>}=\frac{k_{\theta}sd}{\sqrt{2}\Real\left(w^{2}\right)}\sqrt{\Real\left(w^{2}\right)+2\Imag\left(x_{0}\right)^{2}},$
(297)

where we use

$\left<\theta^{2}\right>=\frac{\int_{-\infty}^{\infty}\theta^{2}\left|\hat{\phi}\left(\theta\right)\right|^{2}d\theta}{\int_{-\infty}^{\infty}\left|\hat{\phi}\left(\theta\right)\right|^{2}d\theta}=\frac{d^{2}}{2\Real\left(w^{2}\right)}+\left(d\frac{\Imag\left(x_{0}\right)}{\Real\left(w^{2}\right)}\right)^{2}.$
(298)

Meanwhile, the nonlinear contribution has been tuned (Citrin et al., 2012)
such that

$k_{r-\text{NL}}\rho_{s}=0.4e^{-2\left|s\right|}q^{-0.5}+1.5\max\left\\{k_{\theta}\rho_{s}-0.2,0\right\\}.$
(299)

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

##  10 Numerical Implementation

Recall that the dispersion relation is written as

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

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

$D\left(\omega\right)=\sum_{s}\mathcal{L}_{0,s}-\mathcal{L}_{\text{passing},s}-\mathcal{L}_{\text{trapped},s}.$
(301)

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

$\frac{1}{2\pi
i}\oint_{C}\frac{f^{\prime}\left(z\right)}{f\left(z\right)}dz=N-P,$ (302)

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

$\frac{1}{2\pi
i}\oint_{C}\frac{f^{\prime}\left(z\right)}{f\left(z\right)}dz=N.$ (303)

The key of the method is to recognize from Cauchy’s residue theorem that, for
integer $n$ such that $1\leq n\leq N$, we can calculate the integral $S_{n}$
such that

$S_{n}=\frac{1}{2\pi
i}\oint_{C}z^{n}\frac{f^{\prime}\left(z\right)}{f\left(z\right)}dz=\sum_{j=1}^{N}z_{0j}^{n},$
(304)

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

$P_{N}\left(z\right)=\prod_{j=1}^{N}\left(z-z_{0j}\right)=\sum_{j=0}^{N}A_{j}z^{N-j},$
(305)

where the coefficients $A_{j}$ can be computed from the relations

$\displaystyle A_{0}$ $\displaystyle=1,$ (306) $\displaystyle S_{1}+A_{1}$
$\displaystyle=0,$ (307) $\displaystyle S_{2}+A_{1}S_{1}+2A_{2}$
$\displaystyle=0,$ (308) $\displaystyle
S_{n}+A_{1}S_{n-1}+A_{2}S_{n-2}+\dots+nA_{n}$ $\displaystyle=0,\qquad
n=1,2,\dots,N,$ (309)

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

$S^{\left(1\right)}_{n}=S_{n}-z_{01}^{n},\qquad n=0,1,\dots,N-1\\\ $ (310)

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

$S_{n}=-\frac{n}{2\pi i}\int_{0}^{2\pi}d\theta
e^{in\theta}\ln\left(e^{-iN\theta}f\left(e^{i\theta}\right)\right),\qquad
n>0.$ (311)

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

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

$f\left(z\right)=D\left(\omega\left(z\right)\right).$ (312)

Then, we compute $S_{n}$ via numerical quadrature, leading to roots $z_{0n}$
such that

$f\left(z_{0n}\right)=D\left(\omega\left(z_{0n}\right)\right)=0$ (313)

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

$\omega_{0n}=\omega\left(z_{0n}\right),$ (314)

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

$D\left(\omega_{0n}\right)=0.$ (315)

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

$\displaystyle x\left(u,v\right)$
$\displaystyle=\frac{\text{sgn}\left(uv\right)}{v\sqrt{2}}\sqrt{u^{2}+v^{2}-\sqrt{\left(u^{2}+v^{2}\right)\left(u^{2}+v^{2}-4u^{2}v^{2}\right)}},$
(316) $\displaystyle y\left(u,v\right)$
$\displaystyle=\frac{\text{sgn}\left(uv\right)}{u\sqrt{2}}\sqrt{u^{2}+v^{2}-\sqrt{\left(u^{2}+v^{2}\right)\left(u^{2}+v^{2}-4u^{2}v^{2}\right)}}.$
(317)

The inverse mapping is given by

$\displaystyle u\left(x,y\right)$
$\displaystyle=\frac{x\sqrt{x^{2}+y^{2}-x^{2}y^{2}}}{\sqrt{x^{2}+y^{2}}},$
(318) $\displaystyle v\left(x,y\right)$
$\displaystyle=\frac{y\sqrt{x^{2}+y^{2}-x^{2}y^{2}}}{\sqrt{x^{2}+y^{2}}}.$
(319)

Since this mapping does not satisfy the Cauchy-Riemann equations, it is merely
bijective, not conformal. This is known as a squircle mapping since it appears
to be a square with rounded edges, and this specific kind was first formulated
in Guasti (1992). Denoting $\omega=x^{\prime}+iy^{\prime}$ and $z=u+iv$, we
modify this mapping such that

$\displaystyle x^{\prime}$
$\displaystyle=R_{x}+\frac{r_{x}}{a}x\left(au,av\right),$ (320) $\displaystyle
y^{\prime}$
$\displaystyle=R_{y}+\frac{1}{a}\left(R_{y}-\epsilon_{y}\right)y\left(au,av\right),$
(321)

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

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

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

The goal of hcubature is to estimate

$\bm{I}=\int_{a_{1}}^{b_{1}}\int_{a_{2}}^{b_{2}}\dots\int_{a_{n}}^{b_{n}}\bm{f}(\bm{x})\text{d}^{n}x.$
(322)

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

$\epsilon_{i}=\left|(I_{0})_{i}-(I_{1})_{i}\right|.$ (323)

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

$I=\int_{-1}^{1}\int_{-1}^{1}\dots\int_{-1}^{1}f(x_{1},x_{2},\dots,x_{n})\text{d}^{n}x.$
(324)

We estimate the integral using the following rule,

$\begin{split}I\approx I_{0}=&w_{1}f\left(0,0,\dots,0\right)+w_{2}\sum
f\left(\lambda_{2},0,0,\dots,0\right)+w_{3}\sum
f\left(\lambda_{3},0,0,\dots,0\right)\\\ &+w_{4}\sum
f\left(\lambda_{4},\lambda_{4},0,0,\dots,0\right)+w_{5}\sum
f\left(\lambda_{5},\lambda_{5},\dots,\lambda_{5}\right).\end{split}$ (325)

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

$\begin{split}\sum
f\left(\lambda_{1},0,0\right)=f\left(\lambda_{1},0,0\right)+&f\left(0,\lambda_{1},0\right)+f\left(0,0,\lambda_{1}\right)\\\
&+f\left(-\lambda_{1},0,0\right)+f\left(0,-\lambda_{1},0\right)+f\left(0,0,-\lambda_{1}\right).\end{split}$
(326)

Genz and Malik constrain the parameters $w_{i}$ and $\lambda_{i}$ by requiring
that the integration be exact for the functions

$\begin{split}&f_{1}\left(x_{1},x_{2},\dots,x_{n}\right)=1,\\\
&f_{2}\left(x_{1},x_{2},\dots,x_{n}\right)=x_{1}^{2},\\\
&f_{3}\left(x_{1},x_{2},\dots,x_{n}\right)=x_{1}^{4},\\\
&f_{4}\left(x_{1},x_{2},\dots,x_{n}\right)=x_{1}^{6},\\\
&f_{5}\left(x_{1},x_{2},\dots,x_{n}\right)=x_{1}^{2}x_{2}^{2},\\\
&f_{6}\left(x_{1},x_{2},\dots,x_{n}\right)=x_{1}^{4}x_{2}^{2},\\\
&f_{7}\left(x_{1},x_{2},\dots,x_{n}\right)=x_{1}^{2}x_{2}^{2}x_{3}^{2}.\end{split}$
(327)

In addition, they also fix the parameters $\lambda_{3}=\lambda_{4}$ to be a
specific number, and solve the resulting nonlinear system of equations. The
result can be found in Genz & Malik (1980). To estimate the error, we reuse
$\lambda_{i}$ but calculate different weights $w^{\prime}_{i}$ to make a
lower-order estimate,

$\begin{split}I\approx
I_{1}=w^{\prime}_{1}f\left(0,0,\dots,0\right)+&w^{\prime}_{2}\sum
f\left(\lambda_{2},0,0,\dots,0\right)\\\ &+w^{\prime}_{3}\sum
f\left(\lambda_{3},0,0,\dots,0\right)+w^{\prime}_{4}\sum
f\left(\lambda_{4},\lambda_{4},0,0,\dots,0\right).\end{split}$ (328)

We calculate the weights with the same method as previously discussed and
require the integration be exact for the functions $f_{1},f_{2},f_{3},$ and
$f_{5}$. By keeping $\lambda_{i}$ the same, we can estimate the error without
reusing function evaluations. The error is taken to be

$\epsilon=\left|I_{0}-I_{1}\right|.$ (329)

The estimate procedure easily generalizes to that of a hyperrectangle by using
linear transformations. The calculation of $I_{0}$, $I_{1}$, and
$\bm{\epsilon}$ can also be extended to the case of vector integrands by
integrating every component simultaneously.

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

$I=\int_{-1}^{1}f(x)\text{d}x\approx\sum_{i=1}^{n}w^{\prime}_{i}f(x_{i}).$
(330)

To calculate the weights $w_{i}$ and the abscissa $x_{i}$, we require that the
integration be exact for all polynomials up to degree $2n-1$. It can be shown
using Lagrange interpolating polynomials and the theory of orthogonal
polynomials that the abscissa $x_{i}$ correspond to the roots of the Legendre
polynomial $P_{n}$ and that the weights are

$w^{\prime}_{i}=\frac{2}{\left(1-x_{i}^{2}\right)\frac{dP_{n}}{dx}\left(x_{i}\right)^{2}},$
(331)

where the Legendre polynomials are normalized such that $P_{n}(1)=1$.

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

$I\approx I_{0}=\sum_{i=1}^{15}w_{i}f(x_{i}),$ (332)

the 7-point rule to

$I\approx I_{1}=\sum_{i=1}^{7}w^{\prime}_{i}f(x_{i}),$ (333)

and the estimated error

$\epsilon=\left|I_{0}-I_{1}\right|.$ (334)

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

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

Input: $\bm{f}$, $\bm{a}$, $\bm{b}$, $\epsilon_{a}$, $\epsilon_{r}$, maxEval,
norm

Output: $\bm{I}$, $\bm{e}$, ifail

1 Initialize $\text{eval}=0$

2 Create a hyperrectangle from $\bm{a}$, $\bm{b}$

3 Calculate $\bm{I}_{0}$, $\bm{I}_{1}$, $\bm{\epsilon}$ in the hyperrectangle

4 Calculate $s$, the suggested dimension along which to further discretize, in
the hyperrectangle

5 Update eval to be the number of points evaluated so far

6 $\bm{I}=\bm{I}_{0},\bm{e}=\bm{\epsilon}$

7 Push hyperrectangle into the binary heap with associated value
$\max_{i}\left|\epsilon_{i}\right|$ and with stored values
$\bm{I}_{0},\bm{\epsilon},s$

8 if  _converged_ then exit;

9 else

10 repeat

11 Pop a hyperrectangle from the binary heap

12 Update $\bm{I}=\bm{I}-\bm{I}_{0}$, $\bm{e}=\bm{e}-\bm{\epsilon}$ from the
popped hyperrectangle

13 Split the hyperrectangle in half along the suggested dimension $s$

14 Calculate $\bm{I}_{0}$, $\bm{I}_{1}$, $\bm{\epsilon}$, and $s$ for each
hyperrectangle

15 Update eval to be the number of points evaluated so far

16 Update $\bm{I}=\bm{I}+\sum\bm{I_{0}},\bm{e}=\bm{e}+\sum\bm{\epsilon}$ from
the two hyperrectangles

17 Push each hyperrectangle into the binary heap with associated value
$\max_{i}\left|\epsilon_{i}\right|$ and with stored values
$\bm{I}_{0},\bm{\epsilon},s$

18

19 until  _converged or $\text{eval}\geq\text{maxEval}$_;

20

Algorithm 1 hcubature

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

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

$\begin{split}D_{i}=\sum_{j}&\bigl{|}f_{j}(0,0,\dots,0,-\lambda_{2},0,0,\dots,0)\\\
&-2f_{j}(0,0,\dots,0)+f_{j}(0,0,\dots,0,\lambda_{2},0,0\dots,0)\bigr{.}\\\
&-\bigl{.}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}\left[f_{j}(0,0,\dots,0,-\lambda_{1},0,0,\dots,0)-2f_{j}(0,0,\dots,0)\right.\\\
&\left.\qquad+f_{j}(0,0,\dots,0,\lambda_{1},0,0\dots,0)\right]\bigr{|}.\end{split}$
(335)

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

$\begin{split}D_{2}=\sum_{j}&\bigl{|}f_{j}(0,-\lambda_{2},0,0,\dots,0)-2f_{j}(0,0,\dots,0)+f_{j}(0,\lambda_{2},0,0\dots,0)\bigr{.}\\\
&-\bigl{.}\frac{\lambda_{1}^{2}}{\lambda_{2}^{2}}\left[f_{j}(0,-\lambda_{1},0,0,\dots,0)-2f_{j}(0,0,\dots,0)+f_{j}(0,\lambda_{1},0,0\dots,0)\right]\bigr{|}.\end{split}$
(336)

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

##  11 Conclusions and Outlook

In this work, we derived the linear dispersion relation of quasilinear
gyrokinetic transport code QuaLiKiz from first principles. With the aid of
nonlinear simulations, we also extended the linear physics to a quasilinear
regime to calculate particle, toroidal angular momentum, and heat fluxes. The
formulation of QuaLiKiz relies upon multiple theoretical principles in fusion
plasma physics. First, we examined single particle motion in a circular
magnetic geometry and identified the adiabatic invariants of motion within a
guiding center framework. This allowed us to characterize electrostatic
perturbations to the system with the aid of action-angle variables. We used
this formulation to analyze the linearized Vlasov equation and Poisson’s
equation. We then simplified the resulting dispersion relation using the
ballooning representation, an eigenfunction ansatz, and various
approximations. The solution of this dispersion relation is computed using the
Davies method and numerical cubature methods. Finally, upon finding the
eigenmodes of the system, we use the solutions to compute the quasilinear
fluxes with the aid of a saturation rule informed by nonlinear simulations.

This derivation serves not only to help explain the inner workings of the
model, but also to guide potential improvement in QuaLiKiz. With the
formulation finally laid out, it is now clear where each individual
approximation enters the derivation. This will ease future QuaLiKiz
development that aims to extend the underlying physics or relax the various
approximations. Examples of such work includes introducing electromagnetic
perturbations, incorporating a more general magnetic geometry, and a more
accurate pitch angle integration for passing particles. Improvements made to
QuaLiKiz will allow for more accurate integrated modeling as well as further
optimization of the code.

An additional goal of this work is to provide an extensive, line-by-line
derivation for the sake of demonstrating how such a model can be formulated in
principle. Explicitly drawing upon multiple theoretical principles, such as
the action-angle variable formalism, helps to illustrate the utility of these
principles and their physical motivation. It is also useful to lay out the
various mathematical and numerical techniques necessary in a model such as
this, since many such tricks, methods, or approximations are often crucial in
making a problem computationally tractable. We hope that this work will
function not just as a tutorial for understanding and improving QuaLiKiz, but
also further development in quasilinear fusion codes in general.

##  Appendix A Fried and Conte Integrals

The Fried and Conte integral, also known as the plasma dispersion function, is
utilized frequently in kinetic plasma physics. It is defined as

$Z(x)=\begin{cases}\frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty}dv\frac{e^{-v^{2}}}{v-x},&\text{if}\
\Imag(x)>0,\\\
\mathcal{P}\frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty}dv\frac{e^{-v^{2}}}{v-x}+\sqrt{\pi}ie^{-x^{2}},&\text{if}\
\Imag(x)=0,\\\
\frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty}dv\frac{e^{-v^{2}}}{v-x}+2\sqrt{\pi}ie^{-x^{2}},&\text{if}\
\Imag(x)<0,\\\ \end{cases}$ (337)

where the case $\Imag(x)\leq 0$ is calculated by analytically continuing the
integral defined for $x>0$. When solving the Vlasov equation as an initial
value problem in time, a Laplace transform is implied when obtaining this
integral. To apply the Laplace transform correctly for the case of stable
modes, we must analytically continue the function. Luckily, since we only
consider unstable modes, we are free to restrict ourselves instead to the
related function

$Z_{0}(x)=\frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty}dv\frac{e^{-v^{2}}}{v-x}.$
(338)

If $\Imag(x)=0$, we take the Cauchy principle value of $Z_{0}$.

In carrying out the calculation, we define a generalization of the plasma
dispersion function defined as

$Z_{m}(x)=\frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty}dv\frac{v^{m}e^{-v^{2}}}{x-v},$
(339)

where $m$ is a nonnegative integer. It can be shown that these associated
Fried and Conte integrals can be written in terms of $Z_{0}(x)$:

$Z_{m}(x)=\begin{cases}x^{m}Z_{0}(x)+\frac{1}{\sqrt{\pi}}\sum_{k=0}^{\frac{m-1}{2}}x^{2k}\Gamma(\frac{m}{2}-k),&\text{if
m odd},\\\
x^{m}Z_{0}(x)+\frac{1}{\sqrt{\pi}}\sum_{k=0}^{\frac{m}{2}-1}x^{2k+1}\Gamma(\frac{m-1}{2}-k),&\text{if
m even},\end{cases}$ (340)

where $\Gamma(x)$ is the gamma function. For integer $n$ we note that
$Z_{2n+1}\left(x\right)$ is an even function and $Z_{2n}\left(x\right)$ is
odd. The first few of these associated Fried and Conte integrals are

$\displaystyle Z_{1}(x)$ $\displaystyle=1+xZ_{0}(x),$ (341) $\displaystyle
Z_{2}(x)$ $\displaystyle=x+x^{2}Z_{0}(x),$ (342) $\displaystyle Z_{3}(x)$
$\displaystyle=\frac{1}{2}+x^{2}+x^{3}Z_{0}(x),$ (343) $\displaystyle
Z_{4}(x)$ $\displaystyle=\frac{x}{2}+x^{3}+x^{4}Z_{0}(x).$ (344)

We also define a further generalization of the Fried and Conte integral as
described in Gürcan (2014):

$G_{m}(x_{1},x_{2})=\frac{1}{\sqrt{\pi}}\int^{\infty}_{-\infty}dv\frac{v^{m}e^{-v^{2}}}{\left(v-x_{1}\right)\left(v-x_{2}\right)}$
(345)

Through partial fraction decomposition, we can rewrite this as

$G_{m}(x_{1},x_{2})=\frac{1}{\sqrt{\pi}}\int^{\infty}_{-\infty}dv\left(\frac{1}{v-x_{1}}+\frac{x_{2}}{\left(v-x_{1}\right)\left(v-x_{2}\right)}\right)v^{m-1}e^{-v^{2}}=Z_{m-1}(x_{1})+x_{2}G_{m-1}(x_{1},x_{2}).$
(346)

Because $G_{m}(x_{1},x_{2})=G_{m}(x_{2},x_{1})$, we obtain

$G_{m}(x_{1},x_{2})=Z_{m-1}(x_{1})+x_{2}G_{m-1}(x_{1},x_{2})=G_{m}(x_{2},x_{1})=Z_{m-1}(x_{2})+x_{1}G_{m-1}(x_{1},x_{2}),$
(347)

which allows us to write

$G_{m}(x_{1},x_{2})=\frac{Z_{m}(x_{1})-Z_{m}(x_{2})}{x_{1}-x_{2}}.$ (348)

Note that $G_{m}(x_{1},x_{2})=G_{m}(-x_{1},-x_{2})$.

##  Appendix B Derivation of the Magnetic Drift Velocity

The goal of this section is to calculate the magnetic drift velocity
$\bm{v}_{D,B}$ in the $s-\alpha$ equilibrium by including a finite Shafranov
shift. We define the right-handed coordinate system
$\left(r,\theta,\varphi\right)$ using Cartesian coordinates and include the
Shafranov shift explicitly:

$\displaystyle x$
$\displaystyle=\left(R_{0}+r\cos(\theta)+\Delta(r)\right)\cos(\varphi),$ (349)
$\displaystyle y$
$\displaystyle=\left(R_{0}+r\cos(\theta)+\Delta(r)\right)\sin(\varphi),$ (350)
$\displaystyle z$ $\displaystyle=r\sin(\theta).$ (351)

Here, $\Delta$ is the outward radial shift of the circular flux surface’s
center. The coordinate system $\left(r,\theta,\varphi\right)$ is right-handed
but not orthogonal, so we must specify the metric coefficients. They are

$\displaystyle g_{rr}$
$\displaystyle=1+\left(\Delta^{\prime}\right)^{2}+2\Delta^{\prime}\cos\left(\theta\right),$
(352) $\displaystyle g_{r\theta}$ $\displaystyle=g_{\theta
r}=-\Delta^{\prime}r\sin\left(\theta\right),$ (353) $\displaystyle
g_{\theta\theta}$ $\displaystyle=r^{2},$ (354) $\displaystyle
g_{\varphi\varphi}$
$\displaystyle=\left(R_{0}+r\cos\left(\theta\right)+\Delta\right)^{2},$ (355)
$\displaystyle g_{r\varphi}$ $\displaystyle=g_{\varphi
r}=g_{\theta\varphi}=g_{\varphi\theta}=0,$ (356)

where $\Delta^{\prime}=\partial_{r}\Delta$. This leads to the Jacobian

$J=\sqrt{g}=\frac{1}{\nabla
r\cdot\left(\nabla\theta\times\nabla\varphi\right)}=r\left(R_{0}+r\cos\left(\theta\right)+\Delta\right)\left(1+\Delta^{\prime}\cos\left(\theta\right)\right).$
(357)

We next define a magnetic field for the $s-\alpha$ equilibrium. As an
approximation, we ignore the poloidal magnetic field and only consider the
toroidal magnetic. Thus, the magnetic field is

$\bm{B}\approx B_{0}R_{0}\nabla\varphi.$ (358)

This guarantees that the magnetic field strength is

$B=\frac{B_{0}R_{0}}{R_{0}+r\cos\left(\theta\right)+\Delta}=\frac{B_{0}R_{0}}{R},$
(359)

where $R=R\left(r,\theta\right)$. It is well known that one can obtain an
approximate expression for $\Delta^{\prime}$ from the Grad-Shafranov equation
to lowest-order in $\epsilon$ (Connor et al., 1983; Candy, 2009; Linder,
2016). The expression is

$\Delta^{\prime}\approx-\alpha=q^{2}\beta\frac{R_{0}}{P}\frac{dP}{dr}.$ (360)

The next step is to calculate the magnetic drift velocity,

$\bm{v}_{D,B}=\frac{m}{eB}\left(v_{\parallel}^{2}+\frac{v_{\perp}^{2}}{2}\right)\frac{\bm{B}\times\nabla
B}{B^{2}}+\frac{mv_{\parallel}^{2}}{eB}\frac{\beta}{2p}\frac{\bm{B}\times\nabla
p}{B}.$ (361)

The first term is the sum of the grad-$B$ drift as well as the dominant
component of the curvature drift. The second term is the portion of the
curvature drift that arises from considering the lowest-order MHD equilibrium.
Since QuaLiKiz is applied in the regime where $\alpha$ is small, we ignore the
second term entirely; this is equivalent to assuming that the magnetic field
is approximately curl-free. Taking note that we are not using an orthogonal
coordinate system, we find that the relevant cross product is

$\left(\bm{B}\times\nabla B\right)=R_{0}B_{0}\left(\frac{\partial B}{\partial
r}\nabla\varphi\times\nabla r+\frac{\partial
B}{\partial\theta}\nabla\varphi\times\nabla\theta\right).$ (362)

We can evaluate each component of the expression to obtain

$\displaystyle\left(\bm{B}\times\nabla B\right)\cdot\nabla r$
$\displaystyle=R_{0}B_{0}\frac{\partial B}{\partial\theta}\nabla
r\cdot\left(\nabla\varphi\times\nabla\theta\right)=-\frac{R_{0}B_{0}}{J}\frac{\partial
B}{\partial\theta}=-\frac{B^{2}\sin\left(\theta\right)}{R\left(1+\Delta^{\prime}\cos\left(\theta\right)\right)},$
(363) $\displaystyle\left(\bm{B}\times\nabla B\right)\cdot\nabla\theta$
$\displaystyle=R_{0}B_{0}\frac{\partial B}{\partial
r}\nabla\theta\cdot\left(\nabla\varphi\times\nabla
r\right)=\frac{R_{0}B_{0}}{J}\frac{\partial B}{\partial
r}=-\frac{B^{2}\left(\cos\left(\theta\right)+\Delta^{\prime}\right)}{Rr\left(1+\Delta^{\prime}\cos\left(\theta\right)\right)},$
(364) $\displaystyle\left(\bm{B}\times\nabla B\right)\cdot\nabla\varphi$
$\displaystyle=0.$ (365)

We then use the approximation

$\displaystyle\frac{1}{R_{0}\left(1+\Delta^{\prime}\cos\left(\theta\right)\right)}\approx\frac{1}{R_{0}}\left(1-\Delta^{\prime}\cos\left(\theta\right)\right)$
(366)

and substitute in $\Delta^{\prime}=-\alpha$ to obtain to lowest-order

$\displaystyle\bm{v}_{D}\cdot\nabla r$ $\displaystyle\approx-
v_{D,B}\sin\left(\theta\right),$ (367)
$\displaystyle\bm{v}_{D}\cdot\nabla\theta$
$\displaystyle\approx-\frac{v_{D,B}}{r}\left(\cos\left(\theta\right)-\alpha\sin^{2}\left(\theta\right)\right),$
(368) $\displaystyle\bm{v}_{D}\cdot\nabla\varphi$ $\displaystyle\approx 0,$
(369)

where we define the characteristic magnetic drift speed to be

$v_{D,B}=\frac{m}{eBR_{0}}\left(v_{\parallel}^{2}+\frac{v_{\perp}^{2}}{2}\right).$
(370)

##  Appendix C Collisions

The main sections of this work only consider the collisionless Vlasov
equation. In actuality, QuaLiKiz implements a Krook-type collision operator
for trapped electrons. To account for its inclusion, we modify the Vlasov
equation to

$\frac{\partial\delta f_{s}}{\partial t}+\bm{\Omega}\cdot\frac{\partial\delta
f_{s}}{\partial\bm{\alpha}}-e_{s}\frac{\partial\phi}{\partial\bm{\alpha}}\cdot\frac{\partial
f_{0s}}{\partial\bm{J}}=-\nu\left(\delta
f_{s}+\frac{e\phi}{T_{s}}f_{0s}\right),$ (371)

where $\nu$ is the collision frequency. Note that the $e_{s}\phi f_{0s}/T_{s}$
term accounts for the adiabatic response from the electrostatic perturbation.
We only include this term for electron-ion collisions as ion-ion collisions
and electron-electron collisions would produce only a small correction. Thus,
we drop the “s” in favor of “e” and take $e_{s}\to-e$. Substituting in our
Fourier expressions for $\delta f$ and $\phi$, we find that

$f_{\bm{n}}=\frac{f_{0e}}{T_{e}}\frac{-e\phi_{\bm{n}}\left(\bm{n}\cdot\bm{\omega}_{\ast}-\bm{n}\cdot\bm{\Omega}-\nu\right)}{\bm{n}\cdot\bm{\Omega}-\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu-i\nu}=\frac{e\phi_{\bm{n}}}{T_{e}}f_{0e}\left(1-\frac{\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu-\bm{n}\cdot\bm{\omega}_{\ast}}{\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu+i\nu-\bm{n}\cdot\bm{\Omega}}\right).$ (372)

Therefore, we can simply substitute $\omega\to\omega+i\nu$ in the denominator
of the resonant term to capture the effect of this collision operator. The
drawback is that we lose the ability to simplify the functional. In QuaLiKiz,
we take the collisional frequency to be

$\nu_{e}(\xi,\lambda,\epsilon)=\nu_{ei}\left(\xi\right)^{-3/2}Z_{\text{eff}}\frac{\epsilon}{(1-2\epsilon-\lambda)^{2}}\frac{0.111\delta+1.31}{11.79\delta+1},$
(373)

$\nu_{ei}$ is the electron-ion Coulomb collision frequency, $Z_{\text{eff}}$
is the effective charge of the ion species interacting with the electrons, and
the parameter $\delta$ is defined as

$\delta=12.0\left(\frac{\left|\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu\right|\epsilon}{\nu_{ei}Z_{eff}}\right)^{3/2}.$ (374)

The explicit definition of $\nu_{ei}$ is

$\nu_{ei}=\frac{e^{4}\lambda_{e}}{4\pi\epsilon_{0}^{2}\left(2T_{e}\right)^{3/2}m_{e}^{1/2}},$
(375)

where $\lambda_{e}$ is the Coulomb logarithm relevant for electron collisions.
Details for this collision operator can be found in Romanelli et al. (2007).
The numerical values as well as the derivation of $\delta$ were originally
calculated in Kotschenreuther et al. (1995) and then modified for QuaLiKiz’s
purposes. Because $\nu$ is a function non-trivial of $\xi$, we cannot simplify
the functional using this collision operator using Fried and Conte integrals,
and the integration over the energy must be done numerically. The inability to
simplify the $\xi$ integration results in $2$-dimensional integral. That
aside, all other aspects of the trapped functional derivation remain intact.

We note this specific form of the collision operator is modified in comparison
to the one found in Romanelli et al. (2007). It was found that the previous
form of the collision operator led to incorrect predictions for density
profiles when used in QuaLiKiz (coupled to integrated modeling suites) in
highly collisional regimes. In response, numerical parameters in the Krook
operator were tuned to linear simulations in GENE. In doing so, we keep
unchanged the generic dependence and numerical parameters derived from
fundamental principles unchanged (Stephens et al., 2021).

## References

  * Baños (1967) Baños, A. 1967 The guiding centre approximation in lowest order. Journal of Plasma Physics 1 (3), 305–316.
  * Baumgaertel et al. (2011) Baumgaertel, J. A., Belli, E. A., Dorland, W., Guttenfelder, W., Hammett, G. W., Mikkelsen, D. R., Rewoldt, G., Tang, W. M. & Xanthopoulos, P. 2011 Simulating gyrokinetic microinstabilities in stellarator geometry with GS2. Physics of Plasmas 18 (12), 122301.
  * Biglari & Chen (1986) Biglari, H. & Chen, L. 1986 Theory of energetic trapped particle-induced resistive interchange-ballooning modes. The Physics of Fluids 29 (9), 2960–2974.
  * Bourdelle (2000) Bourdelle, C. 2000 Analyse de stabilité de plasmas de tokamak. PhD thesis, Université Joseph Fourier-Grenoble I, France.
  * Bourdelle (2015) Bourdelle, C. 2015 Turbulent Transport in Tokamak Plasmas: briding theory and experiment. Habilitation thesis, Université de Provence, Aix-Marseille I, France.
  * Bourdelle et al. (2016) Bourdelle, C., Citrin, J., Baiocchi, B., Casati, A., Cottier, P., Garbet, X. & and, F. I. 2016 Core turbulent transport in tokamak plasmas: bridging theory and experiment with QuaLiKiz. Plasma Physics and Controlled Fusion 58 (1), 014036\.
  * Bourdelle et al. (2002) Bourdelle, C., Garbet, X., Hoang, G., Ongena, J. & Budny, R. 2002 Stability analysis of improved confinement discharges: internal transport barriers in tore supra and radiative improved mode in TEXTOR. Nuclear Fusion 42 (7), 892–902.
  * Bourdelle et al. (2007) Bourdelle, C., Garbet, X., Imbeaux, F., Casati, A., Dubuit, N., Guirlet, R. & Parisot, T. 2007 A new gyrokinetic quasilinear transport model applied to particle transport in tokamak plasmas. Physics of Plasmas 14 (11), 112501.
  * Brizard (2011) Brizard, A. J. 2011 Compact formulas for guiding-center orbits in axisymmetric tokamak geometry. Physics of Plasmas 18 (2), 022508.
  * Brizard & Hahm (2007) Brizard, A. J. & Hahm, T. S. 2007 Foundations of nonlinear gyrokinetic theory. Reviews of Modern Physics 79, 421–468.
  * Candy (2009) Candy, J. 2009 A unified method for operator evaluation in local grad–shafranov plasma equilibria. Plasma Physics and Controlled Fusion 51 (10), 105009.
  * Cary & Brizard (2009) Cary, J. R. & Brizard, A. J. 2009 Hamiltonian theory of guiding-center motion. Reviews of Modern Physics 81, 693–738.
  * Casati et al. (2009) Casati, A., Bourdelle, C., Garbet, X., Imbeaux, F., Candy, J., Clairet, F., Dif-Pradalier, G., Falchetto, G., Gerbaud, T., Grandgirard, V., Gürcan, Ö. D., Hennequin, P., Kinsey, J., Ottaviani, M., Sabot, R., Sarazin, Y., Vermare, L. & Waltz, R. 2009 Validating a quasi-linear transport model versus nonlinear simulations. Nuclear Fusion 49 (8), 085012.
  * Citrin et al. (2012) Citrin, J., Bourdelle, C., Cottier, P., Escande, D. F., Gürcan, . D., Hatch, D. R., Hogeweij, G. M. D., Jenko, F. & Pueschel, M. J. 2012 Quasilinear transport modelling at low magnetic shear. Physics of Plasmas 19 (6), 062305\.
  * Citrin et al. (2017) Citrin, J. & others 2017 Tractable flux-driven temperature, density, and rotation profile evolution with the quasilinear gyrokinetic transport model qualikiz. Plasma Physics and Controlled Fusion 59 (12), 124005.
  * Connor et al. (1983) Connor, J., Hastie, R. & Martin, T. 1983 Effect of pressure gradients on the bounce-averaged particle drifts in a tokamak. Nuclear Fusion 23 (12), 1702–1704.
  * Connor et al. (1979) Connor, J. W., Hastie, R. J. & Taylor, J. B. 1979 High mode number stability of an axisymmetric toroidal plasma. Proceedings of the Royal Society A 365, 1–17.
  * Cottier et al. (2014) Cottier, P., Bourdelle, C., Camenen, Y., Gürcan, Ö. D., Casson, F. J., Garbet, X., Hennequin, P. & Tala, T. 2014 Angular momentum transport modeling: achievements of a gyrokinetic quasi-linear approach. Plasma Physics and Controlled Fusion 56 (1), 015011.
  * Dannert & Jenko (2005) Dannert, T. & Jenko, F. 2005 Gyrokinetic simulation of collisionless trapped-electron mode turbulence. Physics of Plasmas 12 (7), 072309.
  * Davies (1986) Davies, B. 1986 Locating the zeros of an analytic function. Journal of Computational Physics 66 (1), 36–49.
  * Depret et al. (2000) Depret, G., Garbet, X., Bertrand, P. & Ghizzo, A. 2000 Trapped-ion driven turbulence in tokamak plasmas. Plasma Physics and Controlled Fusion 42 (9), 949–971.
  * Garbet (2001) Garbet, X. 2001 Instabilities, turbulence and transport in a magnetized plasma. Habilitation thesis, Université de Provence, Aix-Marseille I, France.
  * Garbet et al. (2013) Garbet, X., Esteve, D., Sarazin, Y., Abiteboul, J., Bourdelle, C., Dif-Pradalier, G., Ghendrih, P., Grandgirard, V., Latu, G. & Smolyakov, A. 2013 Turbulent acceleration and heating in toroidal magnetized plasmas. Physics of Plasmas 20 (7), 072502.
  * Garbet et al. (1990) Garbet, X., Laurent, L., Mourgues, F., Roubin, J. & Samain, A. 1990 Variational calculation of electromagnetic instabilities in tokamaks. Journal of Computational Physics 87 (2), 249–269.
  * Genz & Malik (1980) Genz, A. & Malik, A. 1980 Remarks on algorithm 006: An adaptive algorithm for numerical integration over an N-dimensional rectangular region. Journal of Computational and Applied Mathematics 6 (4), 295–302.
  * Goldstein et al. (2001) Goldstein, H., Poole, C. & Safko, J. 2001 Classical Mechanics. Addison-Wesley, San Francisco.
  * Görler & Jenko (2008) Görler, T. & Jenko, F. 2008 Scale separation between electron and ion thermal transport. Phys. Rev. Lett. 100, 185002.
  * Guasti (1992) Guasti, M. 1992 Analytic geometry of some rectilinear figures. International Journal of Mathematical Education in Science and Technology 23, 895–902.
  * Gürcan (2014) Gürcan, Ö. D. 2014 Numerical computation of the modified plasma dispersion function with curvature. Journal of Computational Physics 269, 156–167.
  * Horton (1984) Horton, W. 1984 Drift wave turbulence and anomalous transport. In Handbook of Plasma Physics, Basic Plasma Physics, Volume II (ed. A A Galeev & R N Sudan). Elsevier, New York.
  * Howard et al. (2016) Howard, N., Holland, C., White, A., Greenwald, M. & Candy, J. 2016 Multi-scale gyrokinetic simulation of tokamak plasmas: enhanced heat loss due to cross-scale coupling of plasma turbulence. Nuclear Fusion 56 (1), 014004.
  * Howes et al. (2008) Howes, G. G., Dorland, W., Cowley, S. C., Hammett, G. W., Quataert, E., Schekochihin, A. A. & Tatsuno, T. 2008 Kinetic simulations of magnetized turbulence in astrophysical plasmas. Phys. Rev. Lett. 100, 065004.
  * Howes et al. (2011) Howes, G. G., TenBarge, J. M., Dorland, W., Quataert, E., Schekochihin, A. A., Numata, R. & Tatsuno, T. 2011 Gyrokinetic simulations of solar wind turbulence from ion to electron scales. Phys. Rev. Lett. 107, 035004.
  * Jenko et al. (2005) Jenko, F., Dannert, T. & Angioni, C. 2005 Heat and particle transport in a tokamak: advances in nonlinear gyrokinetics. Plasma Physics and Controlled Fusion 47 (12B), B195–B206.
  * Jenko et al. (2000) Jenko, F., Dorland, W., Kotschenreuther, M. & Rogers, B. N. 2000 Electron temperature gradient driven turbulence. Physics of Plasmas 7 (5), 1904–1910.
  * Jenko & Kendl (2002) Jenko, F. & Kendl, A. 2002 Stellarator turbulence at electron gyroradius scales. New Journal of Physics 4, 35–35.
  * Johnson (1991) Johnson, C. 1991 Numerical Solution of Partial Differential Equations by the Finite Element Method. Cambridge University Press, Cambridge.
  * Johnson (2017) Johnson, S. G. 2017 Cubature.
  * Kaufman (1972) Kaufman, A. N. 1972 Quasilinear diffusion of an axisymmetric toroidal plasma. The Physics of Fluids 15 (6), 1063–1069.
  * Kinsey et al. (2008) Kinsey, J. E., Staebler, G. M. & Waltz, R. E. 2008 The first transport code simulations using the trapped gyro-landau-fluid model. Physics of Plasmas 15 (5), 055908.
  * Kotschenreuther et al. (1995) Kotschenreuther, M., Rewoldt, G. & Tang, W. 1995 Comparison of initial value and eigenvalue codes for kinetic toroidal plasma instabilities. Computer Physics Communications 88 (2), 128–140.
  * Krommes (2002) Krommes, J. A. 2002 Fundamental statistical descriptions of plasma turbulence in magnetic fields. Physics Reports 360 (1), 1–352.
  * Kruskal (1962) Kruskal, M. 1962 Asymptotic theory of hamiltonian and other systems with all solutions nearly periodic. Journal of Mathematical Physics 3 (4), 806–828.
  * Kubo (1963) Kubo, R. 1963 Stochastic liouville equations. Journal of Mathematical Physics 4 (2), 174–183.
  * Lin et al. (2008) Lin, Y., Rice, J. E., Wukitch, S. J., Greenwald, M. J., Hubbard, A. E., Ince-Cushman, A., Lin, L., Porkolab, M., Reinke, M. L. & Tsujii, N. 2008 Observation of ion-cyclotron-frequency mode-conversion flow drive in tokamak plasmas. Phys. Rev. Lett. 101, 235002.
  * Linder (2016) Linder, O. 2016 Comparison of tokamak linear microstability calculations between the gyrokinetic codes QuaLiKiz and GENE. Intership report, Eindhoven University of Technology.
  * Mahajan & Chen (1985) Mahajan, S. M. & Chen, C. Y. 1985 Plasma kinetic theory in action-angle variables. The Physics of Fluids 28 (12), 3538–3545.
  * Merz & Jenko (2008) Merz, F. & Jenko, F. 2008 Nonlinear saturation of trapped electron modes via perpendicular particle diffusion. Phys. Rev. Lett. 100, 035005.
  * Merz & Jenko (2010) Merz, F. & Jenko, F. 2010 Nonlinear interplay of TEM and ITG turbulence and its effect on transport. Nuclear Fusion 50 (5), 054005.
  * Mynick et al. (2010) Mynick, H. E., Pomphrey, N. & Xanthopoulos, P. 2010 Optimizing stellarators for turbulent transport. Phys. Rev. Lett. 105, 095004.
  * Navarro et al. (2016) Navarro, A. B., Teaca, B., Told, D., Groselj, D., Crandall, P. & Jenko, F. 2016 Structure of plasma heating in gyrokinetic alfvénic turbulence. Phys. Rev. Lett. 117, 245101.
  * Nguyen et al. (2008) Nguyen, C., Garbet, X. & Smolyakov, A. I. 2008 Variational derivation of the dispersion relation of kinetic coherent modes in the acoustic frequency range in tokamaks. Physics of Plasmas 15 (11), 112502.
  * Nunami et al. (2010) Nunami, M., Watanabe, T.-H. & Sugama, H. 2010 Gyrokinetic vlasov code including full three-dimensional geometry of experiments. Plasma and Fusion Research 5, 016–016.
  * Pueschel et al. (2011) Pueschel, M. J., Jenko, F., Told, D. & Büchner, J. 2011 Gyrokinetic simulations of magnetic reconnection. Physics of Plasmas 18 (11), 112102.
  * Pueschel et al. (2014) Pueschel, M. J., Told, D., Terry, P. W., Jenko, F., Zweibel, E. G., Zhdankin, V. & Lesch, H. 2014 Magnetic reconnection turbulence in strong guide fields: Basic properties and application to coronal heating. The Astrophysical Journal Supplement Series 213 (2), 30.
  * Rogers et al. (2007) Rogers, B. N., Kobayashi, S., Ricci, P., Dorland, W., Drake, J. & Tatsuno, T. 2007 Gyrokinetic simulations of collisionless magnetic reconnection. Physics of Plasmas 14 (9), 092110.
  * Romanelli et al. (2007) Romanelli, M., Regnoli, G. & Bourdelle, C. 2007 Numerical study of linear dissipative drift electrostatic modes in tokamaks. Physics of Plasmas 14 (8), 082305, arXiv: https://doi.org/10.1063/1.2755981.
  * Samain (1970) Samain, A. 1970 Dynamic stabilization of a confined plasma. Nuclear Fusion 10 (3), 325–335.
  * Stephens et al. (2017) Stephens, C. D., Brzozowski III, R. W. & Jenko, F. 2017 On the limitations of gyrokinetics: Magnetic moment conservation. Physics of Plasmas 24 (10), 102517.
  * Stephens et al. (2021) Stephens, C. D., Citrin, J., van de Plassche, K., Bourdelle, C., Tala, T., Salmi, A. & Jenko, F. 2021 Quasilinear modeling of collisional trapped electron modes. to be submitted to Nuclear Fusion .
  * Stephens et al. (2020) Stephens, C. D., Garbet, X. & Jenko, F. 2020 Analytic guiding center formulas for bounce-transit motion in a concentric circular, finite inverse aspect ratio tokamak geometry. Physics of Plasmas 27 (5), 052504.
  * Told et al. (2016) Told, D., Cookmeyer, J., Muller, F., Astfalk, P. & Jenko, F. 2016 Comparative study of gyrokinetic, hybrid-kinetic and fully kinetic wave physics for space plasmas. New Journal of Physics 18 (6), 065011.
  * Told et al. (2015) Told, D., Jenko, F., TenBarge, J. M., Howes, G. G. & Hammett, G. W. 2015 Multiscale nature of the dissipation range in gyrokinetic simulations of alfvénic turbulence. Phys. Rev. Lett. 115, 025003.
  * Waltz et al. (2007) Waltz, R. E., Candy, J. & Fahey, M. 2007 Coupled ion temperature gradient and trapped electron mode to electron temperature gradient mode gyrokinetic simulations. Physics of Plasmas 14 (5), 056116.
  * Wan et al. (2005) Wan, W., Chen, Y. & Parker, S. E. 2005 Gyrokinetic delta-f simulation of the collisionless and semicollisional tearing mode instability. Physics of Plasmas 12 (1), 012311\.
  * Wang et al. (2008) Wang, X. Y., Lin, Y., Chen, L. & Lin, Z. 2008 A particle simulation of current sheet instabilities under finite guide field. Physics of Plasmas 15 (7), 072103.
  * Xanthopoulos & Jenko (2007) Xanthopoulos, P. & Jenko, F. 2007 Gyrokinetic analysis of linear microinstabilities for the stellarator Wendelstein 7-X. Physics of Plasmas 14 (4), 042501.