# The zonal-flow residual does not tend to zero in the limit of small mirror
ratio

E. Rodríguez  G. G. Plunk  Max Planck Institute for Plasma Physics, 17491
Greifswald, Germany

###### Abstract

The intensity of the turbulence in tokamaks and stellarators depends on its
ability to excite and sustain zonal flows. Insight into this physics may be
gained by studying the “residual”, i.e. the late-time linear response of the
system to an initial perturbation. We investigate this zonal-flow residual in
the limit of a small magnetic mirror ratio, where we find that the typical
quadratic approximation to RH (Rosenbluth & Hinton, 1998) breaks down. Barely
passing particles are in this limit central in determining the resulting level
of the residual, which we estimate analytically. The role played by the
population with large orbit width provides valuable physical insight into the
response of the residual beyond this limit. Applying this result to tokamak,
quasi-symmetric and quasi-isodynamic equilibria, using a near-axis
approximation, we identify the effect to be more relevant (although small) in
the core of quasi-axisymmetric fields, where the residual is smallest. The
analysis in the paper also clarifies the relationship between the residual and
the geodesic acoustic mode, whose typical theoretical set-ups are similar.

##  1 Introduction

There exists a strong current interest in exploring the space of stellarators
(Spitzer Jr, 1958; Boozer, 1998; Helander, 2014), three-dimensional, toroidal
magnetic confinement fields. Optimising such fields in order to achieve plasma
confinement and ultimately controlled thermonuclear fusion requires of careful
design and shaping of the field for it to present desired physical properties.
In guiding this search, it is imperative to have a good understanding of the
key physics involved. Given the breadth of the stellarator concept, though,
this naturally requires stretching our understanding of physics that are
comparatively mature in the simpler case of the axisymmetric tokamak
(Mukhovatov & Shafranov, 1971; Wesson, 2011).

Amongst the critical elements that govern the behaviour of a stellarator,
turbulence is a particularly interesting and important one. Understanding the
neoclassical behaviour of stellarators has historically captivated much of the
focus of research, mainly because of its predominant role in the transport of
unoptimised stellarators through the so-called $1/\nu$ regime (Galeev et al.,
1969; Stringer, 1972; Ho & Kulsrud, 1987; Nemov et al., 1999; Mynick, 2006).
Progress over the last decades, and especially over the past years (Beidler et
al., 2021; Landreman & Paul, 2022; Goodman et al., 2023), has however brought
turbulence to the forefront, and it is now regarded as one of the key elements
determining the performance of stellarators.

Zonal flow dynamics are of particular interest in the study of turbulence
(Diamond et al., 2005), as they are understood to play a key role in
regulating turbulence by shearing eddies apart, lowering the overall intensity
of turbulent fluctuations. The description of full zonal-flow dynamics is
certainly complex, as an essentially non-linear response of the system.
However, one may learn some basic information about the ability for a given
magnetic equilibrium to sustain such flows by considering the behaviour of the
so-called zonal-flow residual (Rosenbluth & Hinton, 1998; Xiao & Catto, 2006;
Sugama & Watanabe, 2006; Monreal et al., 2016). The residual is the long-time
remnant of an initial radially varying perturbation of the electrostatic
potential. The prevalence of a large such remnant is, at least sometimes,
indicative of the system’s capacity to sustain zonal dynamics in a turbulent
state (Watanabe et al., 2008; Xanthopoulos et al., 2011). The calculation of
the residual thus serves as a reasonable starting point for the assessment of
zonal flows in a given magnetic equilibrium. The main theoretical
understanding of the residual behaviour was pioneered by Rosenbluth & Hinton
(1998), and subsequently refined and extended by others (Xiao & Catto, 2006;
Sugama & Watanabe, 2006; Monreal et al., 2016; Plunk & Helander, 2024),
including in the electromagnetic context (Catto et al., 2017).

The level of the residual depends strongly on the size of the orbit-width,
$\delta$, of the particles in the field, that is, the magnitude of the
particle deviation from flux surfaces as they move along field lines. The
dependence is so strong that, in a typical scenario (Rosenbluth & Hinton,
1998), it is the trapped particles (whose orbit widths are largest) that
contribute most to the residual. The larger the orbit widths, the lower the
residual levels, as the shielding from these becomes more effective
(Rosenbluth & Hinton, 1998; Xiao & Catto, 2006). In fact, it is conventionally
argued that in the limit of $B$ becoming flat (small mirror ratio), the large
trapped particle orbits cause the residual to vanish. Of course it is also in
this limit that there are also no trapped particles left in the problem,
somewhat complicating the asymptotic analysis.

In this paper we revisit the theoretical question of the zonal-flow residual
in this limit. An assessment is presented in Section 2, where we also draw
connections to the standard framework of geodesic-acoustic-modes (Conway et
al., 2021). We learn that barely passing particles play the dominant role in
determining the final finite value of the residual in the small mirror ratio
limit. This large-orbit-width part of the population behaves, we argue, as if
non-omnigeneous, as far as the residual is concerned. We find support for
these claims numerically through linear gyrokinetic simulations. We close the
discussion in Section 4 with an assessment of the relevance of this effect on
tokamaks and omnigeneous stellarators, which appears to be limited.

##  2 Residual calculation in the small mirror ratio limit

###  2.1 Brief derivation of the residual

Let us start our discussion on the zonal-flow residual by calculating it in
its most typical of set-ups. We follow closely the work of Rosenbluth & Hinton
(1998); Xiao & Catto (2006); Monreal et al. (2016); Plunk & Helander (2024),
but include a brief derivation for completeness and as a way of introduction
of notation.

By residual, which we denote $\phi(\infty)$, we mean the surface averaged
collisionless electrostatic potential in the long time limit. To describe it,
we take the linearised, electrostatic gyrokinetic equation as starting point
(Connor et al., 1978, 1980),

$\left(\frac{\partial}{\partial
t}+i\tilde{\omega}_{d}+v_{\parallel}\frac{\partial}{\partial\ell}\right)g=\frac{q}{T}F_{0}J_{0}\frac{\partial}{\partial
t}\phi,$ (1)

written in the ballooning formalism with the variation perpendicular to the
field line described by $\mathbf{k}_{\perp}=k_{\psi}\nabla\psi$. Here $\psi$
is the flux surface label (the toroidal flux over $2\pi$), so that the
electrostatic potential perturbation $\phi$ has a main strong off-surface
variation, which is the reason why there is no diamagnetic term in Eq. (1),
$\omega_{\star}=0$. Other symbols have their usual meaning: $F_{0}$ is the
background Maxwellian distribution, $J_{0}=J_{0}(x_{\perp}\sqrt{2b})$ the
Bessel function of the first kind representing Larmor radius effects and
$b=(k_{\psi}|\nabla\psi|\rho)^{2}/2$ the Larmor radius parameter, with
$\rho=v_{T}/\Omega$, $v_{T}=\sqrt{2T/m}$ and $\Omega=q\bar{B}/m$ (at this
point we are considering a general species of mass $m$, charge $q$ and
temperature $T$). The drift frequency
$\tilde{\omega}_{d}=\omega_{d}(v/v_{T})^{2}(1-\lambda B/2)$ and
$\omega_{d}=\mathbf{v}_{D}\cdot\mathbf{k}_{\perp}=v_{T}\rho\bar{B}k_{\psi}{\boldsymbol{\kappa}}\times\mathbf{B}\cdot\nabla\psi/B^{2}$,
with $\bar{B}$ a reference field, ${\boldsymbol{\kappa}}$ the curvature of the
field and the drift is considered in the low $\beta$ limit. The velocity space
variables are $\lambda=\mu/\mathcal{E}$ and particle velocity
$v=\sqrt{2\mathcal{E}/m}$, where $\mu$ is the first adiabatic invariant and
$\mathcal{E}$ the particle energy. The parallel velocity can then be written
as $v_{\parallel}=\sigma v\sqrt{1-\lambda B}$, where $\sigma$ is the sign of
$v_{\parallel}$.

Equation (1) is then a partial differential equation in time $t$ and the arc
length along the field line $\ell$, for the electrostatic potential $\phi$ and
the non-adiabatic part of the distribution function, $g$, with a dependence on
the velocity space variables $\\{\sigma,v,\lambda\\}$. Performing a Laplace
transform in time (Schiff, 2013, Theorem 2.7) yields

$\left(\omega-\tilde{\omega}_{d}+iv_{\parallel}\frac{\partial}{\partial\ell}\right)\hat{g}=\frac{q}{T}F_{0}J_{0}\omega\hat{\phi}+i\delta\\!F(0),$
(2)

where
$\delta\\!F(0)\stackrel{{\scriptstyle\cdot}}{{=}}g(0)-(q/T)J_{0}F_{0}\phi(0)$
can be interpreted as the initial perturbation of the system, and we are using
the hats to indicate the Laplace transform.

To eliminate the explicit $\ell$ dependence that the curvature,
$\tilde{\omega}_{d}$, brings into the equation, we shall define the orbit
width $\delta$,

$v_{\parallel}\frac{\partial}{\partial\ell}\delta=\tilde{\omega}_{d}-\overline{\tilde{\omega}_{d}}$
(3)

so that we may write,

$\left(iv_{\parallel}\frac{\partial}{\partial\ell}-\overline{\tilde{\omega}_{d}}+\omega\right)\hat{h}=\frac{q}{T}F_{0}\omega
J_{0}\hat{\phi}e^{i\delta}+ie^{i\delta}\delta\\!F(0),$ (4)

and $\hat{h}=\hat{g}e^{i\delta}$. The function $\delta$ describes the off-
surface displacement of particles (in $\psi$) as a function of $\ell$, for
each particle identified by its velocity space labels. The overline notation
indicates the bounce average,

$\overline{f}=\begin{cases}\begin{aligned}
&\frac{1}{\tau_{b}}\frac{1}{v}\int_{\mathrm{b}}\frac{1}{\sqrt{1-\lambda
B}}\sum_{\sigma}f\,\mathrm{d}\ell,\\\
&\lim_{L\rightarrow\infty}\frac{1}{\tau_{t}}\frac{1}{v}\int_{\mathrm{p}}\frac{f\,\mathrm{d}\ell}{\sqrt{1-\lambda
B}}.\end{aligned}\end{cases}$ (5)

The first expression applies to trapped particles, where the integral is taken
between the left and right bounce points and summed over both directions
($\sigma$) of the particle’s motion. The normalisation factor is the bounce
time, $\tau_{b}$, defined following $\overline{1}=1$. For passing particles,
the integral is taken over the whole flux surface (i.e. the infinite extent of
the field line explicitly indicated by the limit), and normalised by the
transit time, $\tau_{t}$.

When $\overline{\tilde{\omega}_{d}}=0$, Eq. (4) simplifies. This corresponds
to the physical interpretation of particles having no net off-surface drift.
This is the defining property of omnigeneity (Hall & McNamara, 1975a; Cary &
Shasharina, 1997; Helander & Nührenberg, 2009; Landreman & Catto, 2012), which
we shall assume to hold throughout this work. For a treatment of the non-
omnigeneous problem see Helander et al. (2011); Monreal et al. (2016).

Because we are interested in the behaviour at large time scales, we expand in
$\omega/\omega_{t}\sim\epsilon_{t}$, applying
$\hat{h}=\hat{h}^{(0)}+\hat{h}^{(1)}+\dots$ and
$\hat{\phi}=\hat{\phi}^{(0)}+\hat{\phi}^{(1)}+\dots$, and considering Eq. (4)
order by order,

$\displaystyle iv_{\parallel}\frac{\partial}{\partial\ell}\hat{h}^{(0)}$
$\displaystyle\approx 0,$ (6a) $\displaystyle
iv_{\parallel}\frac{\partial}{\partial\ell}\hat{h}^{(1)}+\omega\hat{h}^{(0)}$
$\displaystyle\approx\frac{q}{T}\omega
F_{0}J_{0}\hat{\phi}^{(0)}e^{i\delta}+ie^{i\delta}\delta\\!F(0),$ (6b)
$\displaystyle\vdots$

From Eq. (6a) it follows that,

$\hat{h}^{(0)}=\overline{\hat{h}^{(0)}}.$ (7)

Thus, bounce averaging Eq. (6b), and assuming that $\hat{\phi}^{(0)}$ is
$\ell$-independent, we may write down the leading order expression for
$\hat{g}^{(0)}$,

$\hat{g}^{(0)}=\frac{q}{T}F_{0}e^{-i\delta}\overline{J_{0}e^{i\delta}}\hat{\phi}^{(0)}+\frac{i}{\omega}e^{-i\delta}\overline{\delta\\!F(0)e^{i\delta}}.$
(8)

With this expression for $\hat{g}$, we may then apply the quasineutrality
condition (Connor et al., 1980) summing over ions and electrons. Explicitly,
and summing over electrons and ions (subscripts $e$ and $i$ respectively)

$\sum_{e,i}\int\mathrm{d}^{3}\mathbf{v}J_{0}\hat{g}=n\frac{q_{i}}{T_{i}}(1+\tau)\hat{\phi},$
(9)

where $\tau=T_{i}/ZT_{e}$ and $Z=-q_{i}/q_{e}$, then yields

$\hat{\phi}^{(0)}\approx\frac{1}{n}\left\langle\int\mathrm{d}^{3}\mathbf{v}J_{0}e^{-i\delta}\overline{J_{0}e^{i\delta}}F_{0}\right\rangle_{\psi}\hat{\phi}^{(0)}+\frac{i}{\omega}\frac{1}{n}\left\langle\int\mathrm{d}^{3}\mathbf{v}J_{0}e^{-i\delta}\overline{\delta\\!F(0)e^{i\delta}}\right\rangle_{\psi}.$
(10)

Here $\langle\dots\rangle_{\psi}$ denotes a flux surface average (Helander,
2014), and we have taken the limit of $m_{e}\ll m_{i}$, so that the limit of a
negligible electron Larmor radius and electron banana width may be taken; this
is equivalent to an adiabatic electron response
$\phi-\langle\phi\rangle_{\psi}$, making the final form of the residual
independent of electrons.

By inverse Laplace transforming this latest expression (Schiff, 2013, Theorem
2.36), we obtain,

$\phi(\infty)=\frac{\frac{1}{n}\left\langle\int\mathrm{d}^{3}\mathbf{v}J_{0}e^{-i\delta}\overline{\delta\\!F(0)e^{i\delta}}\right\rangle_{\psi}}{1-\frac{1}{n}\left\langle\int\mathrm{d}^{3}\mathbf{v}J_{0}e^{-i\delta}\overline{J_{0}e^{i\delta}}F_{0}\right\rangle_{\psi}}.$
(11)

To finalise the calculation of the residual, we must consider some initial
perturbation of the ion population. Following Rosenbluth & Hinton (1998);
Monreal et al. (2016), we perturb the density of the ions with
$\delta\\!F(0)=(\delta n/n)J_{0}F_{0}$, a perturbed Maxwellian, sidestepping
the issue of detailed initial-condition dependence of the residual, especially
important at shorter wavelengths (Monreal et al., 2016). Applying
quasineutrality at $t=0$, the density perturbation may be directly related to
the perturbed electrostatic potential $\phi(0)$. Assuming that $b$ is
independent of $\ell$ for simplicity, $\delta
n/n=\phi(0)(1-\Gamma_{0})/\Gamma_{0}$ where $\Gamma_{0}=e^{-b}I_{0}(b)$ and
$I_{0}$ is the Bessel function of the first kind. Therefore, the expression
for the residual at long times is,

$\frac{\phi(\infty)}{\phi(0)}\approx\frac{1-\Gamma_{0}}{\Gamma_{0}}\frac{\frac{1}{n}\left\langle\int\mathrm{d}^{3}\mathbf{v}J_{0}e^{-i\delta}\overline{J_{0}e^{i\delta}}F_{0}\right\rangle_{\psi}}{1-\frac{1}{n}\left\langle\int\mathrm{d}^{3}\mathbf{v}J_{0}e^{-i\delta}\overline{J_{0}e^{i\delta}}F_{0}\right\rangle_{\psi}}.$
(12)

###  2.2 Finite orbit width

In order to proceed with the evaluation of Eq. (12) we first need to study the
orbits of our particles, namely $\delta$. These will depend critically on both
$B(\ell)$ (which controls the time spent by particles along different segments
of the field-line), and the normal curvature $\omega_{d}$ (that determines the
off-surface velocity). Although in an actual equilibrium field these functions
are connected to each other, it is formally convenient to set this equilibrium
connection aside, and treat them as largely independent quantities in the
context of a single flux tube.

Despite this independence, it is important to respect some minimal properties.
First, for the choice of functions to appropriately represent the behaviour in
an omnigeneous field, they should prevent diverging particle orbits. We
prevent this ill-behaviour by ensuring that the critical points of $B(\ell)$
match points of zero radial drift; that is, $\omega_{d}(\ell)=0$ wherever
$\mathrm{d}B(\ell)/\mathrm{d}\ell=0$. This property is known as pseudosymmetry
(Mikhailov et al., 2002; Skovoroda, 2005), and is necessary to represent an
omnigeneous field. However, it is not sufficient. In addition, we must impose
that all the orbits $\delta$ are closed; that is, that they come back to the
same $\psi$ at bounce points, or for passing particles, after a period.

With this, we may write explicitly $\delta$ integrating Eq. (3), as

$\delta=\sigma\frac{v}{v_{T}}\int_{\bar{\ell}_{0}}^{\bar{\ell}}\frac{1-\lambda
B/2}{\sqrt{1-\lambda
B}}\frac{\omega_{d}(\bar{\ell}^{\prime})}{\omega_{t}}\mathrm{d}\bar{\ell}^{\prime}$
(13)

where we have introduced a normalised length scale $\bar{\ell}$ and an
associated transit frequency $\omega_{t}=v_{T}/L$, with $L$ some reference
length scale. The integral is defined so that $\delta(\bar{\ell}_{0})=0$,
where $\bar{\ell}_{0}$ corresponds to bounce points for trapped particles, and
the point $B=B_{\mathrm{max}}$ for passing ones to guarantee continuity across
the trapped-passing boundary.111Note that by virtue of omnigeneity it does not
matter which point of maximum $B$ or bounce point (left or right) along the
field line we choose, because $\delta=0$ at all of these by virtue of
omnigeneity. This property of omnigeneous fields is very important, and it
allows us to treat each well along the field line independently from every
other. This is so because there is no accumulation of radial displacement of
passing particles across maxima. Thus, the considerations that the paper
presents for a single well could be extended to multiple ommnigeneous wells,
treating each separately, and summing their contributions when considering
flux surface averages, as needed in Eq. (12).

The regularising role of pseudosymmetry at critical points of $B(\ell)$, where
it avoids diverging behaviour, can be seen directly from Eq. (13). This allows
us to rewrite $\delta$ in a form that avoids the explicit $1/\sqrt{\cdot}$
divergence using integration by parts,

$\delta=-\sigma\frac{v}{v_{T}}\left[\frac{B^{2}}{\partial_{\bar{\ell}}B}\frac{\omega_{d}(\bar{\ell})}{\omega_{t}}\frac{\sqrt{1-\lambda
B}}{B}\right]_{\bar{\ell}_{0}}^{\bar{\ell}}+\sigma\frac{v}{v_{T}}\int_{\bar{\ell}_{0}}^{\bar{\ell}}\frac{\sqrt{1-\lambda
B}}{B}\partial_{\bar{\ell}^{\prime}}\left(\frac{B^{2}}{\partial_{\bar{\ell}^{\prime}}B}\frac{\omega_{d}(\bar{\ell}^{\prime})}{\omega_{t}}\right)\mathrm{d}{\bar{\ell}}^{\prime}.$
(14)

This integrated form of the equation is also useful to numerically compute
$\delta$ near bounce points.

These expressions are so far quite general, and we shall now specialise to a
simple representative system. In particular, we assume to have a single unique
magnetic well along the field line222Along any fieldline of an omnigeneous
field, every time a maximum of $B$ is crossed, one falls into a new magnetic
well. In the case of a tokamak, all those wells are identical by virtue of
axisymmetry, and thus the consideration of a single unique well is sufficient.
Other optimised configurations, though, lack this exact symmetry, which
requires some additional interpretation. Some of this is discussed in Section
4. , described simply by $B=\bar{B}\left(1-\Delta\cos\pi\bar{\ell}\right)$ and
$\omega_{d}=\omega_{d}\sin\pi\bar{\ell}$, where the domain is taken to be
$\bar{\ell}\in[-1,1]$. Thus the scale $L$ can be interpreted as the connection
length in the problem, or the half-width of the well, $\Delta$ the mirror
ratio and $\omega_{d}$ the drift. This particular choice is convenient in two
ways: first, because the choice $\omega_{d}=c\partial_{\ell}B$, with $c$ some
proportioonality constant, simplifies Eq. (14) and conveniently guarantees the
closure of particle orbits; and second, because many of the integrals that
ensue may be carried out exactly for such simple analytic functions. Of
course, deforming these geometric functions away from these forms (in
particular, breaking the parity in $\bar{\ell}$) will directly affect the
orbit shape $\delta$ and ultimately the residual, but this model nonetheless
includes the essential ingredients.

Figure 1: Example of passing and trapped orbits. Numerical examples of trapped
and passing orbits for different values of $\lambda$ for the model field
considered in the paper. The plots were generated for $\Delta=0.05$. The
dotted line on top and bottom correspond to the $\delta(0)$ estimate in Eq.
(18) (grey line simply indicates the reference $\delta=0$ level). Critical
points are marked with solid points.

####  2.2.1 Passing particles

Let us start our description of the passing particle orbits by considering
their maximum deviation off the flux surface, i.e. their orbit widths
$\delta|_{\bar{\ell}=0}=\delta(0)$. By passing particles we refer to the
portion of velocity space with $\lambda\in(0,1/B_{\mathrm{max}})$, which we
may also label with the convenient shifted variable
$\hat{\lambda}=1/(1+\Delta)-\bar{B}\lambda$. In this case $\hat{\lambda}=0$
represents the trapped-passing boundary, and
$\hat{\lambda}=\bar{B}/B_{\mathrm{max}}$ is approached for the passing
particles far from the trapped-passing boundary, which we will refer to as
strongly passing. It is convenient to introduce yet an additional label for
passing particles, namely
$\kappa=2\lambda\bar{B}\Delta/[1-\lambda\bar{B}(1-\Delta)]$, which is bounded
$\kappa\in(0,1)$ and denotes barely passing particles by $\kappa=1$ and
strongly passing by $\kappa=0$.

For the model field considered, $\delta(0)$ may be evaluated exactly in terms
of $\lambda$, $\Delta$ and other parameters. However, it is more insightful to
consider some relevant asymptotic limits. In the limit of a small mirror ratio
$\Delta$, the passing population is naturally separated into three different
regimes, where we may write,

$\delta_{\mathrm{pass}}(0)\approx-\sigma\frac{v}{v_{T}}\frac{\omega_{d}}{\pi\omega_{t}}\times\begin{cases}\begin{aligned}
&\sqrt{\frac{2}{\Delta}}&\quad\text{if }\hat{\lambda}\ll\Delta,\\\
&\frac{1}{\sqrt{\hat{\lambda}}}&\quad\text{if }\Delta\ll\hat{\lambda}\ll 1,\\\
&\frac{1}{\sqrt{\hat{\lambda}}}+\sqrt{\hat{\lambda}}&\quad\text{if
}\hat{\lambda}\gg\Delta.\end{aligned}\end{cases}$ (15)

The orbits are widest within a layer of width $\Delta$ near the trapped-
passing boundary, where all barely passing particles have large, almost
identical orbits that scale like $\sim 1/\sqrt{\Delta}$. This is a consequence
of particles moving slowly along the field line by an amount
$v_{\parallel}\sim\sqrt{1-\lambda B}\sim\sqrt{\Delta}$. Thus, there always
exists a sufficiently small mirror ratio able to slow down barely passing
particles enough so as for them to have a sizeable orbit width; this is true
even for a small radial drift $\epsilon\equiv\omega_{d}/\pi\omega_{t}\ll 1$.

We estimate the size of the $v$-space layer that includes particles with a
sizeable orbit width (i.e. $|\delta_{\mathrm{pass}}(0)|>1$) in the limit of
$\epsilon\ll 1$ by taking the behaviour of a typical thermal particle
$v/v_{T}\sim 1$ as reference in Eq. (15), so that

$\hat{\lambda}<\epsilon^{2}=\left(\frac{\omega_{d}}{\pi\omega_{t}}\right)^{2}.$
(16)

Such a layer can only exist if the mirror ratio is sufficiently small,

$\Delta/\epsilon^{2}\ll 1.$ (17)

Not satisfying this mirror ratio ordering restores the standard view of
passing particles having small orbit widths (as in the quadratic approximation
of the residual in Rosenbluth & Hinton (1998)). The small mirror ratio
ordering alongside the $\epsilon\ll 1$ assumption are henceforth assumed.

####  2.2.2 Trapped particles

The procedure above may be repeated for trapped particles. Defining a trapped
particle label
$\bar{\kappa}=1/\kappa=[1/(\lambda\bar{B})-(1-\Delta)]/2\Delta$, deeply
trapped particles are denoted by $\bar{\kappa}=0$ and barely trapped ones by
$\bar{\kappa}=1$. The orbit width may then be written as,

$\delta_{\mathrm{trap}}(0)\approx-\sigma\frac{v}{v_{T}}\epsilon\sqrt{\frac{2\bar{\kappa}}{\Delta}},$
(18)

assuming $\Delta\ll 1$. Unlike passing particles, the majority of trapped
particles have a significant orbit width (in the $\sim 1/\sqrt{\Delta}$
sense), except for a minute fraction near the bottom of the well which barely
moves away from that point. This fraction may be estimated to be

$\bar{\kappa}<\frac{\Delta}{\epsilon^{2}},$ (19)

which we have already assumed small.

###  2.3 Evaluating the residual for small mirror ratio

Figure 2: Separation of particles into groups. The diagram depicts the
separation of the particle population into four different groups (I to IV).
Groups I and IV (light blue) represent the population with a small orbit
width, while II and III (light red) correspond to large ones. The diagram is a
schematic with the vertical representing $1/\lambda$, the horizontal
$\bar{\ell}$ and the black line representing the magnetic well
$B(\bar{\ell})$.

In the limit of a small mirror ratio, we have learned from the analysis of the
orbits that the particle population may be divided into four different groups.
Each of these groups is characterised by having a large or small $\delta$, and
thus a different contribution to Eq. (12). We refer to each of these groups by
Roman numerals I to IV, starting from strongly passing particles (see Figure
2).

To proceed with the residual integral, let us assume for simplicity the
finite-Larmor quantity $b$ to be small. This is compatible with $\epsilon$
being small (note that $\epsilon\propto k_{\perp}\rho_{i}$). With this, we may
write the integral in the denominator of the residual, Eq. (12),

$1-\frac{1}{n}\left\langle\int\mathrm{d}^{3}\mathbf{v}J_{0}e^{-i\delta}\overline{J_{0}e^{i\delta}}F_{0}\right\rangle_{\psi}\approx
b-\frac{1}{n}\left\langle\int\mathrm{d}^{3}\mathbf{v}\left(e^{-i\delta}\overline{e^{i\delta}}-1\right)F_{0}\right\rangle_{\psi},$
(20)

where we used,

$\frac{1}{n}\left\langle\int\mathrm{d}^{3}\mathbf{v}J_{0}^{2}F_{0}\right\rangle_{\psi}=\frac{1}{2}e^{-b}I_{0}(b),$
(21)

in the small $b$ limit and the velocity space integrals include all groups.
The integral remaining in Eq. (20) has been simplified by dropping finite-
Larmor radius corrections. For groups I and IV for which $\delta$ is small,
retaining $b$ would give an even smaller $O(\delta^{2}b)$ correction, which we
drop. For groups II and III, the correction would also be small in the sense
$O(b\sqrt{\Delta})$, under the assumption of small $\Delta$.

Now separating the integral left in Eq. (20) into the different group
contributions,

$I=\frac{1}{n}\left\langle\int\mathrm{d}^{3}\mathbf{v}\left(e^{-i\delta}\overline{e^{i\delta}}-1\right)F_{0}\right\rangle_{\psi}=\sum_{\mathrm{I,~{}IV}}\frac{1}{n}\left\langle\int\mathrm{d}^{3}\mathbf{v}\left(\overline{\delta}^{2}-\overline{\delta^{2}}\right)F_{0}\right\rangle_{\psi}+\\\
+\sum_{\mathrm{II,~{}III}}\frac{1}{n}\left\langle\int\mathrm{d}^{3}\mathbf{v}\left(e^{-i\delta}\overline{e^{i\delta}}-1\right)F_{0}\right\rangle_{\psi}.$
(22)

This separation enables us to exploit the smallness or largeness of $\delta$
accordingly. The smallness of the orbit width for groups I and IV has already
been exploited to write the leading order contribution in powers of $\delta$
in the first term of the right-hand side of Eq. (22). This contribution should
be familiar, as it has the quadratic form in which the Rosenbluth-Hinton
residual is customarily written (Rosenbluth & Hinton, 1998; Xiao & Catto,
2006; Plunk & Helander, 2024). We set this part of the calculation aside for
now, and focus on the new contributions by groups II and III.

####  2.3.1 Contribution from barely passing particles (group II)

Let us continue our analysis by looking at barely passing particles in group
II (see Fig. 2), and their contribution to Eq. (22),

$I_{\mathrm{II}}=\underbrace{\frac{1}{n}\left\langle\int_{\mathrm{II}}\mathrm{d}^{3}\mathbf{v}e^{-i\delta}\overline{e^{i\delta}}F_{0}\right\rangle_{\psi}}_{①}-\underbrace{\frac{1}{n}\left\langle\int_{\mathrm{II}}\mathrm{d}^{3}\mathbf{v}F_{0}\right\rangle_{\psi}}_{②}.$
(23)

First consider ①, and rewrite it following Xiao & Catto (2006) as,

$①=\frac{1}{n}\left\langle\int_{\mathrm{II}}\mathrm{d}^{3}\mathbf{v}\left(\overline{\cos\delta}^{2}+\overline{\sin\delta}^{2}\right)F_{0}\right\rangle_{\psi},$
(24)

where we have dropped terms odd in $v_{\parallel}$, annihilated by the
integral over velocity space. Note that, although tempting,
$\overline{\sin\delta}$ is generally nonzero according to our convention for
the bounce average in Eq. (5), where each direction of the passing particles
is treated separately.

To continue with the calculation, we need to evaluate $\overline{\cos\delta}$
explicitly, exploiting that within group II, the function $\delta$ has a large
amplitude. As a result, we expect the cosine of $\delta$ to oscillate quickly
along $\bar{\ell}$ resulting in an almost exact cancellation. The non-zero
contribution may be estimated through the well-known stationary phase
approximation (Bender & Orszag, 2013, Sec. 6.5),

$\overline{\cos\delta}=\frac{1}{\tau_{t}\omega_{t}}\frac{v_{T}}{v}\Re\left\\{\int_{-1}^{1}\frac{e^{i\delta}}{\sqrt{1-\lambda
B}}\mathrm{d}\bar{\ell}\right\\}\approx\frac{1}{\tau_{t}\omega_{t}}\frac{v_{T}}{v}\sum_{i}\sqrt{\frac{2\pi}{|\delta^{\prime\prime}(\ell_{i})|}}\frac{\cos[\delta(\ell_{i})-\pi/4]}{\sqrt{1-\lambda
B(\ell_{i})}},$ (25)

where the sum is over the turning points of $\delta$ in $\bar{\ell}\in[0,1]$.
Using the details of $\delta$ developed in Sec. 2.2.1 and Appendix A,

$\overline{\cos\delta}\approx\frac{2}{\tau_{t}\omega_{t}}\left(\frac{v_{T}}{v}\right)^{3/2}\sqrt{\frac{\omega_{t}}{\omega_{d}}}\left[(4\hat{\lambda})^{-1/4}+\frac{1}{(2\Delta+\hat{\lambda})^{1/4}}\cos\left(\frac{v}{v_{T}}\frac{\epsilon}{\sqrt{\Delta/2+\hat{\lambda}}}-\frac{\pi}{4}\right)\right].$
(26)

The first term inside the square brackets comes from the edge contribution,
and the second from the point of maximum excursion.

Now that we have $\overline{\cos\delta}$ we must integrate over velocity
space, Eq. (24). To do so we introduce the velocity space measure in the
$\\{v,\lambda,\sigma\\}$ coordinate system (already summed over $\sigma$ to
give a factor of 2) (Hazeltine & Meiss, 2003, Sec. 4.4),

$\mathrm{d}^{3}\mathbf{v}\rightarrow\frac{2\pi B}{\sqrt{1-\lambda
B}}v^{2}\mathrm{d}v\mathrm{d}\lambda.$ (27)

and noting that by definition any bounced averaged quantity is
$\bar{\ell}$-independent, write for any function $f$ in our single well,

$\left\langle\int_{\mathrm{II}}\mathrm{d}^{3}\mathbf{v}\bar{f}\right\rangle_{\psi}=\pi\bar{B}\int_{v=0}^{\infty}\int_{\mathrm{II}}v^{2}\frac{v}{v_{T}}\tau_{t}\omega_{t}\bar{f}\mathrm{d}v\mathrm{d}\lambda,$
(28)

correct to leading order in $\Delta$.

The simplifying assumption of a $v$-independent boundary layer in Eq. (16)
allows us to explicitly carry out the integral over $v$ first. Noting the that
with the ordering $\epsilon^{2}/\Delta\gg 1$ (large $A$),

$\displaystyle\int_{0}^{\infty}ve^{-v^{2}}\cos^{2}\left(Av-\frac{\pi}{4}\right)\mathrm{d}v\approx\frac{1}{4},$
(29a) $\displaystyle\int_{0}^{\infty}ve^{-v^{2}}\mathrm{d}v=\frac{1}{2},$
(29b)

we find using the explicit form of the Maxwellian $F_{0}$,

$①\approx\frac{2}{\sqrt{\pi}}\frac{1}{\omega_{d}}\int_{0}^{\epsilon^{2}}\frac{1}{\hat{\tau}_{t}}\left(\frac{1}{\sqrt{\hat{\lambda}}}+\frac{1}{\sqrt{2\Delta+\hat{\lambda}}}\right)\mathrm{d}\hat{\lambda},$
(30)

where $\hat{\tau}_{t}=\tau_{t}(v/v_{T})$ is a function of $\lambda$. In this
form of ① we have already included the contribution from
$\overline{\sin\delta}$, which can be easily shown to be equivalent to that of
the $\overline{\cos\delta}$. To carry out the integral over $\hat{\lambda}$ we
change variables to $\kappa$, defined in Sec. 2.2.1. The integration domain
becomes $\kappa\in[2\Delta/\epsilon^{2},1]$, with an integral measure

$\frac{\mathrm{d}\kappa}{\mathrm{d}\lambda}=2\Delta\bar{B}\left(1+\frac{1-\Delta}{2\Delta}\kappa\right)^{2}.$
(31)

The contribution from the edges of the orbit (the first term in Eq. (30)) can
be shown to be small upon integration over $\kappa$ in the limit of small
$\Delta$. All that is left is the contribution from the point of maximal
excursion, which can be approximated assuming $K(\kappa)\approx\pi/2$,

$①\approx\frac{\epsilon}{\pi^{3/2}}.$ (32)

This concludes the calculation of ①, but ② remains to be found. This
contribution corresponds to finding the fraction of phase space occupied by
the barely passing particles in group II. Using Eq. (28) and the definition of
region II, the integrals over $\kappa$ and $v$ yield,

$②\approx\epsilon.$ (33)

Altogether,

$I_{\mathrm{II}}\approx\frac{\epsilon}{\pi^{3/2}}(1-\pi^{3/2})\approx-0.26\frac{\omega_{d}}{\omega_{t}},$
(34)

yielding an overall negative contribution linear in $k_{\psi}\rho_{i}$.

####  2.3.2 Contribution from the bulk of trapped particles (group III)

A similar approach to that for the barely passing particles may be directly
applied to the trapped particles that constitute group III. Given the
similarities of the calculation we shall be less explicit here.

The evaluation of the integral starts once again by separating the integral
$I_{\mathrm{III}}$ into two parts, ① and ②, like in Eq. (23). In the
calculation of ①, and unlike for passing particles, we only need to consider
the $\overline{\cos\delta}$ term, as $\overline{\sin\delta}=0$ upon summing
over both particle directions, Eq. (5). The $\overline{\cos\delta}$ term may
be computed much like in the previous section, employing the stationary phase
approach. In this case, the only turning point of $\delta_{\mathrm{trap}}$ is
at the centre of the domain, $\bar{\ell}=0$. With that, using the expressions
for $\delta_{\mathrm{trap}}$ introduced in Sec. 2.2.2 and Appendix A, and
performing the integral over $v$ first,

$①\approx\frac{1}{\pi^{3/2}}\frac{\Delta}{\epsilon},$ (35)

which is a small contribution that vanishes in the limit of $\Delta\rightarrow
0$. The velocity space volume occupied by the bulk of trapped particles, ②, is
of course also small in the limit of a small mirror ratio,
$②\sim\sqrt{\Delta}$. Thus, the contribution to the residual from the trapped
population in group III is small in the limit of $\Delta\rightarrow 0$.

####  2.3.3 Final form of the residual

Gathering the pieces of the calculation above, the integral in Eq. (22)
evaluates to,

$I\approx-0.26\frac{\omega_{d}}{\omega_{t}},$ (36)

in the limit of $\Delta\ll\epsilon^{2}\ll 1$. The latter is particularly
important to argue that the contribution from the particles of groups I and IV
is subsidiary in this limit. We do not need to compute it explicitly to argue
that it scales like $\epsilon^{2}$, and thus is one order $\epsilon$ higher
than the contribution from barely passing particles. Therefore, we may drop
those contributions in writing the result in Eq. (36).

We may now write the expression for the residual itself, going back to Eq.
(12) using the definition of $I$ in Eq. (22),

$\frac{\phi(\infty)}{\phi(0)}\approx\frac{1}{1+0.26\frac{\omega_{d}}{b\omega_{t}}},$
(37)

which in the limit of $b\ll\omega_{d}$, say for very long radial wavelengths,
can be expressed as

$\frac{\phi(\infty)}{\phi(0)}\approx
1.92~{}k_{\perp}\rho_{i}\left(\frac{k_{\perp}\rho_{i}}{\omega_{d}/\omega_{t}}\right).$
(38)

##  3 Analysis of the residual in the small mirror ratio limit

The preceding analysis demonstrates that in the limit of a small mirror ratio
there remains a finite residual in the problem. Barely passing particles near
the passing-trapped boundary dominate the behaviour of the residual in this
limit. This is a result of a narrow $\lambda$-space layer of width
$\epsilon^{2}$ having sufficiently slow parallel velocities so that their
orbits are wide. The result is a partial shielding of the potential. Their
orbit width is so large, though, that their shielding is not as efficient as
it may be at smaller $\delta$, and thus the residual is larger than one would
a priori expect.

There are two important actors that determine the final value of the residual
in this limit: (i) the width of the layer, and (ii) the shape of the orbit.
Both of these may be identified directly in the derivation of the residual
above. The residual will be larger the smaller the layer is, as the shielding
population decreases. The shorter the time that the particles spend near the
point of maximal excursion, the larger the residual will also be; orbit shapes
that are flat near that point are detrimental to the residual.

Figure 3: Example of residual as a function of mirror ratio. The plots present
(a) the time evolution of the average electrostatic potential for different
mirror-ratios simulated with the gyrokinetic code stella, (b) comparison of
residual from the gyrokinetic code stella and numerical evaluation of Eq.
(12), and (c) relative contribution to the residual by passing/trapped
population, and by each $\lambda$. The simulation for (a) and (b) is based on
the cyclone-base-case with $|\mathbf{B}|$ modified, leaving the curvature
drift unchanged. The color code in (a) corresponds to the different mirror
ratios on the right plot, from lower (darker) to larger (brighter) values of
$\Delta$. The right plot (b) presents the residual values from stella as
scatter points (with errorbars indicating the variation of the potential in
the last 20% of the time trace), the triangle marker shows the simulation of
the flat-$B$ scenario, the solid line the numerical evaluation of Eq. (12),
the dotted black line the analytical estimate of Xiao-Catto (Xiao & Catto,
2006), and the red dotted line the asymptotic expression in Eq. (38). The
central bottom plot (c) shows the relative contribution to the residual by
trapped/passing particles. The plots left and right represent the relative
contribution to the residual by different parts of the population, where the
vertical coordinate represents $1/\lambda$, with the black line representing
$B$. The calculations are done at $k_{\perp}\rho_{i}\approx 0.048$
($k_{y}\rho_{i}=0.05$ in stella).

The behaviour of the residual at small mirror ratio can be checked against
both careful numerical integration of Eq. (12) and linear electrostatic
gyrokinetic simulations with the stella code (Barnes et al., 2019). We present
such a comparison in Figure 3. For that comparison, a local field along a flux
tube is constructed from a reference cyclone-base case (a simple Miller
geometry (Miller et al., 1998)) whose $B$ has been modified with varying
mirror ratios $\Delta$, while keeping all other elements of the geometry
unchanged. The numerical evaluation of Eq. (12) is done by careful treatment
of bounce integrals using double-exponential integration methods (Takahasi &
Mori, 1974) to appropriately deal with bounce points and logarithmic
divergences in $\lambda$-space (details on the python code may be found in the
Zenodo repository associated to this paper). The linear gyrokinetic
simulations are run with large velocity space resolution in an attempt to
resolve the boundary layer in velocity space to the best capacity within
reason. This means that they must also be run for long times, on the order of
the transit time of the smallest resolved velocity in order to reach the
residual. We take the residual from these simulations to be the value of the
potential at the latest time simulated.333We are running these simulations in
stella with $N_{v_{\parallel}}=2000$, $N_{\mu}=100$, $\Delta t=0.0125$ and
$N_{t}=64000$, considered high resolutions. The smallest mirror ratio cases
can be challenging to simulate and converge fully even under these extremely
resolved conditions. For the semi-quantitative considerations in this paper we
consider them to be sufficient, though. In addition to these numerical
niceties, the physical oscillations of the electrostatic potential also pose
an additional limitation, as these variations are not damped completely in the
time domain of consideration for the lowest mirror ratios. This can lead to an
inaccurate ‘measured’ residual, but is once again deemed sufficient in the
time domain considered for the semi-quantitative comparison here considered
(see error-bars in Figure 3). Having these two numerical forms of assessing
the residual provides us with additional forms to diagnose the results. In
particular, and given the good agreement between the simulations with the
numerical evaluation of the residual in Eq. (12), we can assess the
contribution from different regions of velocity space to the residual using
the latter (see Figure 3c).

In the small mirror ratio limit, as predicted, there is a dominant
contribution from a narrow boundary layer (group II). The analytic estimate of
the residual in the small mirror ratio, Eq. (38), agrees to a good degree
(within $\sim 5-10\%$) with the simulation and integration (see red line in
Figure 3b). As the mirror ratio increases the importance within velocity space
shifts (see Figure 3c) and the bulk of trapped particles becomes dominant (the
standard Rosenbluth & Hinton (1998) picture). In that limit the residual can
be estimated by Rosenbluth & Hinton (1998) (RH),

$\left.\frac{\phi(\infty)}{\phi(0)}\right|_{\mathrm{RH}}=\frac{1}{1+1.6\epsilon^{2}/(k_{\perp}\rho_{i})^{2}\sqrt{\Delta}},$
(39)

or more precisely by Xiao & Catto (2006), as explicitly shown in Fig. 3b
(black dotted line). The standard RH residual, Eq. (39), exhibits a stronger
dependence on the drift and transit time compared to the small mirror ratio
limit, although the physical mechanism behind the residual remains broadly
speaking the same. Namely, making the drift $\omega_{d}$ or the connection
length smaller, the orbit width becomes smaller, so does the finite orbit
polarisation and shielding power of the plasma, and thus the resulting
residual grows.

The preeminence of the RH or small mirror residual will change depending on
the parameters of both the field and perturbation. A clear example of the
latter is the dependence on $k_{\perp}\rho_{i}$. In fact, for any finite
$\Delta$, there always exists a perpendicular length-scale long enough for
which the RH scenario is recovered (formally, a value of $k_{\psi}$ below
which the ordering $\epsilon^{2}\gg\Delta$ is violated), leading to a finite
residual at small $k_{\perp}\rho_{i}$. Of course, the field parameters also
play a key role. Most clearly, the variation of the mirror ratio $\Delta$
explicitly involves a regime transition between the $\Delta$-independent
small-mirror residual, Eq. (38), and the RH residual (see Figure 3b). This
takes place when $\Delta\sim\epsilon^{2}$, which is approximately

$\Delta_{t}\approx
0.1(k_{\perp}\rho_{i})^{2}\left(\frac{\omega_{d}/\omega_{t}}{k_{\perp}\rho_{i}}\right)^{2}.$
(40)

If the orbit width of the bulk is made larger, then the small-mirror
contribution becomes relevant sooner. However, we must remain within the limit
$\epsilon^{2}\ll 1$, which we considered in the construction of our residual
calculation. Staying within that limit, the transition mirror ratio must obey
$\Delta_{t}<10^{-1}$, which implies that the transition occurs at small mirror
ratios of at most a few per-cent. Of course, the exact value of this
transition will generally not be as simple. We may compute it more accurately
by defining numerically $\Delta_{t}$ as the mirror ratio at which the low
$k_{\perp}\rho_{i}$ limit of the XC (Xiao & Catto, 2006) residual matches the
low-mirror ratio residual.

Before moving to an analysis of these effects on different equilibria, let us
turn to interpreting the time dependence of the residual observed in Figure
3a. There are clearly two oscillation time-scales in the problem set-up
considered: the faster damped geodesic-acoustic modes (GAMs) (Sugama &
Watanabe, 2006; Gao et al., 2006, 2008; Conway et al., 2021) and a slower
oscillation. The former appear rather invariant under $\Delta$ (as one would
expect from a passing ion dominated phenomenon), while the latter change
significantly. In fact, this slower time scale behaviour is reminiscent of the
slower oscillations attributed to the non-omnigeneous nature of stellarator
fields (Mishchenko et al., 2008; Mishchenko & Kleiber, 2012; Helander et al.,
2011; Monreal et al., 2017; Alonso et al., 2017). This provides us with an
additional way of interpreting the boundary layer contribution to the low-
mirror residual. Because of their long transit time compared to their radial
drift, these particles behave de facto as non-omnigeneous particles, at least
in a transient sense. The result are long time scale oscillations with a slow
damping rate. The damping and frequency of oscillations grow in their time
scale as $\Delta$ becomes smaller, which we attribute to the increasingly non-
omnigeneous behaviour of the particles in this limit. A more in-depth
investigation of this behaviour is left for future work.

###  3.1 Geodesic acoustic mode (GAM) connection

From the analysis of the time trace of our simulations, we observe that the
residual and GAMs are just different dynamical phases of the same system. One
then expects to see them both arise consistently in the same asymptotic limit.

GAMs are damped, oscillatory modes resulting from a balance between streaming
and off-surface drift, basic reigning elements in the residual as well. Thus,
these oscillatory modes are, like the residual, often studied as part of the
assessment of the field response to zonal flows. The basic theoretical set-up
for studying GAMs involves a flat-$|\mathbf{B}|$ field, where dynamics are
dominated by passing ions, and the only inhomogeneity along field-lines is
introduced by an oscillatory $\omega_{d}$. Under the assumption of a small
$\omega_{d}/\omega_{t}$ (equivalent to the small $\epsilon$ we have considered
in this paper), the behaviour of GAMs may be reduced to a simple dispersion
relation Sugama & Watanabe (2005, 2006); Gao et al. (2006, 2008). We reproduce
some of the details of this derivation and the dispersion relation in Appendix
B.

The key observation is that the limit $\omega\rightarrow 0$ of these
dispersion relations, which determine the long time behaviour of the
electrostatic potential (Schiff, 2013, Theorem. 2.36), yields no residual. But
we have shown just above that actually a finite residual remains in the limit
of vanishing mirror ratio. A natural question thus arises: where is this
residual hiding? It might be tempting to identify the slow GAM mode identified
by Gao et al. (2006) with the residual, due to its similar form. This purely
damped mode reads

$\frac{\phi(t\rightarrow\infty)}{\phi(0)}\approx\frac{1}{1+\frac{\epsilon^{2}}{4b}\left(1+\frac{\pi}{2(1+\tau)}\right)}e^{-\gamma
t},$ (41)

where,

$\frac{\gamma}{\omega_{t}}=\frac{\pi^{3/2}}{2}\left[\frac{2b}{\epsilon^{2}}+\left(\frac{1}{2}+\frac{\pi}{4(1+\tau)}\right)\right]^{-1}.$
(42)

The amplitude of the mode exhibits a quadratic finite orbit width dependence
much in the fashion of the RH residual. Although the damping of the mode can
be slow (with a characteristic decay time $\sim\epsilon^{2}/b\omega_{t}$), and
thus display an effective value of the residual (transiently), it does not
formally correspond to a collisionless, undamped residual.

In addition, it has a quadratic scaling rather than the linear one derived
above. To resolve this apparent inconsistency we must recognise the importance
of barely passing particles. For this subset of the population the transit
time is so long that the ordering $\omega_{t}\gg\omega_{d}$ is not accurate,
and thus the derivation of the usual GAM dispersion relation needs reworking.
We present the details of how to do this in Appendix B. Doing so, one can
recover a finite valued residual with the same scaling as derived above,
albeit with a different numerical factor. This difference is due to the
difference in the derivation, and gives a factor of 0.20 instead of a 0.26 in
Eq. (37). This reconciling of the residual and GAM calculations is a
theoretical relief.

##  4 Field survey

In the preceding analysis of the residual problem we learned that there are
two different regimes in which the behaviour of the residual is quite
different. One, the regime where the layer dynamics become dominating, which
occurs at small mirror ratios ($\Delta_{t}<10^{-1}$). And the more typical RH
residual one, occurring at moderate values of $\Delta$, in which the bulk of
the trapped particle population dominates the response of the system. We now
explore the question of which regime prevails under the conditions that arise
in different classes of magnetic equilibria.

Let us start with the simplest family of magnetic field configurations: the
circular tokamak. That is, an axisymmetric magnetic field configuration, with
circular cross-sections and thus a unique magnetic well, which is the closest
scenario to our idealised model-field. In such a scenario, we may reduce the
relevant field properties to a few parameters, namely the safety factor $q$,
the mirror ratio $\Delta$ and the radial wavenumber, $k_{\perp}\rho_{i}$. In
the context of the residual, one may think of the safety factor $q$ as
determining the ratio of the radial drift (in a tokamak $\omega_{d}\sim 1/R$)
to the connection length ($\omega_{t}^{-1}\sim qR$), explicitly
$q=\omega_{d}/(\pi k_{\perp}\rho_{i}\omega_{t})$. With that, the relevant
expressions for the residual read, following Eqs. (37) and (39),

$\left.\frac{\phi(\infty)}{\phi(0)}\right|_{\mathrm{lay}}\approx\frac{1}{1+1.63q/(k_{\perp}\rho_{i})},\quad\left.\frac{\phi(\infty)}{\phi(0)}\right|_{\mathrm{RH}}\approx\frac{1}{1+1.6q^{2}/\sqrt{\Delta}}.$
(43)

The larger the $q$, the larger the connection length, the larger the orbit
width $\delta$ and the the lower the residual. In terms of these tokamak
parameters, we may also rewrite the condition for the regime transition in Eq.
(40): the layer contribution becomes relevant for
$\Delta<\Delta_{t}\sim(qk_{\perp}\rho_{i})^{2}$. For a typical value of $q\sim
1$, and a wavenumber $k_{\perp}\rho_{i}\sim 0.1$, this implies mirror ratios
below a percent. This is a rather small mirror ratio, which will only be
reached sufficiently close to the magnetic axis (where $B$ is nearly constant
due to axisymmetry). For shorter wavelengths or larger safety factors (which
also reduce the residual) $\Delta_{t}$ will be larger. Because this occurs at
the expense of larger orbit width, taking this limit to its extreme will
ultimately lead to $\epsilon\sim 1$, implying $\delta>1$ for all particles,
corresponding to a completely different regime.444Large wavenumber behaviour
was explored by (Xiao & Catto, 2006; Monreal et al., 2016). Physically, as the
orbit sizes become large, they become less effective at shielding the original
potential perturbation, and the residual grows. Note however that this
large-$k_{\perp}\rho_{i}$ behaviour is more sensitive to initial conditions
(Monreal et al., 2016) and electron dynamics should be brought in for a
consistent treatment.

To extend the discussion beyond the rather simplified case of circularly
shaped tokamaks, we need some form in which to estimate the input parameters
to our residual calculation. We will focus on so-called optimised stellarator
configurations: namely, quasisymmetric (Boozer, 1983a; Nührenberg & Zille,
1988; Rodríguez et al., 2020) and quasi-isodynamic (Cary & Shasharina, 1997;
Helander & Nührenberg, 2009; Nührenberg, 2010) ones. The former can be seen as
the natural generalisation of the axisymmetric case, where the field has a
direction of symmetry on $|\mathbf{B}|$ instead of the whole vector
$\mathbf{B}$. The direction of symmetry can be toroidal (quasi-axisymmetry) or
helical (quasi-helical). This symmetry forces the magnetic wells along the
field line to be all nearly identical (same $B$ and $\omega_{d}$ (Boozer,
1983b), but different $k_{\perp}\rho_{i}$). In quasi-isodynamic fields, the
contours of $|\mathbf{B}|$ are closed poloidally, and carefully shaped to
grant omnigeneity (Bernardin et al., 1986; Cary & Shasharina, 1997; Hall &
McNamara, 1975b; Helander, 2014). As a result, wells are differently shaped,
but all share the feature of being omnigeneous; that is, the orbits described
by $\delta$ are closed as in Figure 1. The description will in that case have
to involve an average over wells.

Our approach now will be to construct effective model parameters for all of
these configuration types, that may be applied to the above familiar
expressions for the tokamak case, e.g Eqn. 43. These parameters will be
derived using the inverse-coordinate near-axis description of equilibria
(Garren & Boozer, 1991b; Landreman & Sengupta, 2019; Rodríguez et al., 2023;
Plunk et al., 2019), as detailed in Appendix C, and summarised in Table 1. We
have included the case of a shaped tokamak for comparison. Let us now discuss
the interpretation of these results.

| Tokamak | QS | QI
---|---|---|---
$q_{\mathrm{eff}}$ | $\displaystyle\frac{1}{\iota}\frac{\eta R_{\mathrm{ax}}}{\hat{\mathcal{G}}}$ | $\displaystyle\frac{1}{\iota-N}\frac{\eta R_{\mathrm{ax}}}{\hat{\mathcal{G}}}$ | $\displaystyle\frac{2}{\pi N_{\mathrm{nfp}}}\frac{\bar{d}R_{\mathrm{ax}}}{\hat{\mathcal{G}}}$
$\Delta$ | $r\eta$ | $r\eta$ | $\Delta$
Table 1: Characteristic near-axis residual-related parameters in optimised
stellarators. The table presents the value of the residual-relevant parameters
$q_{\mathrm{eff}}$ and $\Delta$ for tokamaks and different optimised
stellarator types, obtained using the near-axis description of the fields (see
Appendix C). The parameters are: $R_{\mathrm{ax}}$ the effective major radius
(the length of the magnetic axis divided by $2\pi$), $\iota$ the rotational
transform, $N$ the symmetry of the QS field, $N_{\mathrm{nfp}}$ number of
field periods, $\eta$ and $\bar{d}$ leading poloidal variation of
$|\mathbf{B}|$ over flux surfaces (roughly proportional to the axis curvature)
and $\hat{\mathcal{G}}$ geometric factor defined in Eq. (44).

The first important distinction between fields is with regards to the
behaviour of the mirror ratio. In tokamaks, as well as quasisymmetric
stellarators, the mirror ratio has a strong radial dependence. In particular,
because $|\mathbf{B}|$ has a direction of symmetry with a toroidal component,
$\Delta$ must decrease towards the axis and do so at a rate related to the
curvature of the field (within the near axis description it is proportional to
the distance form the axis and $\eta\sim\kappa$, see Appendix C). This implies
the appearance of a finite region near the magnetic axis where the low-mirror
residual becomes relevant. In practice, though, this region tends to be
narrow, and thus likely unimportant (see Figure 4s).

Figure 4: Residual and closeness to the residual transition as a function of
radius. The plot shows the residual (top) and the ratio of the mirror ratio
$\Delta$ to the residual regime tranition value $\Delta_{t}$ (bottom) for
DIII-D (equilibrium from Austin et al. (2019), shot 170680 at 2200ms)
(tokamak), precise QA (QA stellarator) and precise QH (QH stellarator)
configurations (Landreman & Paul, 2022). The residual is computed numerically
evaluating Eq. (12) using the global equilibria of the configurations to
estimate the simplified single-well parameters for the residual calculation.
The bottom plots are evaluated computing $\Delta_{t}$ as the mirror ratio
value at which the XC estimate of the residual equals the small mirror ratio
limit of the residual. It therefore is a measure of relevance of the low-
mirror residual regime. It is clear that the centre of the QA configuration is
where the low-mirror ratio is most relevant. The residual calculation was done
for $k_{\perp}\rho_{i}=0.1$ for these.

It is particularly narrow in tokamaks, where the safety factor decreases
towards the axis and can have a significant global shear, unlike
quasisymmetric stellarators (Landreman & Paul, 2022; Landreman, 2022;
Rodríguez et al., 2023; Giuliani, 2024). The consequence of this is also an
inversion of the behaviour of the residual with radius: it tends to be largest
in the core in a tokamak, but smallest for QS ones (see Figure 4). QI
stellarators are significantly different to both tokamaks and QS stellarators.
As a result of having poloidally closed contours, the on-axis $|\mathbf{B}|$
is not constant, and thus the mirror ratio tends to a non-zero constant on the
axis. This frees $\Delta$ from its strong radial dependence, preventing the
low-mirror residual region from manifesting.

In addition to the differences in $\Delta$, the changes in the magnitude of
the magnetic field gradient $\bnabla B$ (which affects $\omega_{d}$), the flux
surface shaping (which affects $k_{\perp}$) and the connection length (which
affects $\omega_{t}$) do also impact the residual. All of these physical
elements may be captured in a parameter $q_{\mathrm{eff}}=\omega_{d}/(\pi
k_{\perp}\rho_{i}\omega_{t})$, given in Table 1. We define such a parameter to
play the role that the safety factor takes in the circular-cross-section
scenario of the residual. In particular, one should interpret this
$q_{\mathrm{eff}}$ as a generalised form of $q$ in the residual expressed in
Eq. (43) and other places. As such, larger $q_{\mathrm{eff}}$ implies lower
residual and a higher relevance of the low-mirror residual regime. Let us
discuss what determines $q_{\mathrm{eff}}$ for each case in Table 1.

We start by analysing the role played by the perpendicular geometry (in
particular $\langle|\nabla\psi|^{2}\rangle$). This is captured by (see Eqns.
84 and 97),

$\hat{\mathcal{G}}^{2}=\frac{1}{2\pi}\int_{0}^{2\pi}\frac{\mathrm{d}\varphi}{\sin
2e},$ (44)

where we define the angle $e$ such that $\mathcal{E}=\tan e$ is the elongation
of the flux surfaces in the plane normal to the axis as a function of
$\varphi$ (Rodríguez, 2023) and we have considered the limit of small mirror
ratio ($\Delta\ll 1$). The angle $e\in(0,\pi/2)$ may be interpreted as the
angle subtended by a right-angle triangle with the major and minor axes as
catheti. Thus, a circular cross-section is represented by $e=\pi/4$, and the
corresponding $\hat{\mathcal{G}}=1$. Any elliptical shape will then have a
larger $\hat{\mathcal{G}}>1$ (as $\sin 2e<1$ for $e\neq\pi/4$ in the domain
considered). Increasing the elongation of flux surfaces increases the average
flux expansion, $\langle|\nabla\psi|^{2}\rangle$, leading to a decrease of
$q_{\mathrm{eff}}$, a larger residual and a decrease in the importance of the
low-mirror residual. This is consistent with Xiao et al. (2007). Physically,
increasing elongation brings flux surfaces closer together, and thus narrows
the orbit widths in real space. Any non-axisymmetric shape will necessarily
have $\hat{\mathcal{G}}>1$ (Landreman & Sengupta, 2019; Camacho Mata et al.,
2022; Rodríguez, 2023), but variations between optimised configurations will
be moderate given that limiting flux surface shaping is often an optimisation
criterion.

Let us now focus on the differences in the magnitude of the magnetic drifts.
The drift is controlled by the gradients of $|\mathbf{B}|$, which decrease the
residual the larger they become. The balance between magnetic gradients (and
thus magnetic pressure) and magnetic field line tension provides an important
observation: the more curved field lines are, the stronger the gradients. In
the near axis framework, this naturally leads to a picture in which the more
strongly shaped a magnetic axis is, the larger the gradients will be. This
behaviour is represented by parameters $\eta$ and $\bar{d}$ in Table 1 (see
Appendix C for a more precise description), which typically scale like
$\eta\sim\kappa$ (Rodríguez et al., 2023), where $\kappa$ is the axis
curvature. For similarly shaped cross-sections, $\eta$ (or $\bar{d}$) will be
larger for QH and QI stellarators compared to QA and tokamaks (Rodriguez et
al., 2022; Camacho Mata et al., 2022), and even more with the number of field
periods. The drift in the QI case deserves special consideration, because the
pointwise radial drift varies from field line to field line, vanishing on some
(Helander & Nührenberg, 2009; Landreman & Catto, 2012). Thus, on ‘average’,
the drift in these configurations is smaller (see Appendix C for the details),
which can enhance the residual. In brief, QH configurations are expected to
have the largest field gradients, followed by QIs in which the field-line
averaging reduces the effective gradients, and finally QAs and tokamaks.

The last element of consideration in $q_{\mathrm{eff}}$ is the connection
length, i.e. the length along the field line of a magnetic well. The
difference in the topology of the $|\mathbf{B}|$ contours (and their alignment
to magnetic field lines) leads to the following comparative scaling,
$R_{\mathrm{ax}}/\iota:R_{\mathrm{ax}}/(\iota-N):R_{\mathrm{ax}}/N_{\mathrm{nfp}}$.
Of course, this naturally leads to ordering the connection lengths to be
largest for QA and tokamaks, smaller for QHs and the smallest for QIs. This
follows from the observation that the number of field periods serves as an
upper bound of $\iota$ for QHs in practice.

The three elements discussed above compete with each other, but the
preeminence of the connection length on $q_{\mathrm{eff}}$ in practice leads
to the relative ordering,

$q_{\mathrm{eff,tok}}\sim q_{\mathrm{eff,QA}}>q_{\mathrm{eff,QH}}\gtrsim
q_{\mathrm{eff,QI}}.$ (45)

This should be regarded as a rough guide, not as a rigid rule; a similar
ordering for the overall size of the residual is argued by Plunk & Helander
(2024).

Figure 5: Parameter $q_{\mathrm{eff}}$ for QS and QI configurations.
Statistics of $q_{\mathrm{eff}}$ for QS and QI configurations. The left plots
represent the normalised (by total area) density of QH and QA configurations
by their value of $q_{\mathrm{eff}}$ in the QS near-axis database in Landreman
(2022), which serves as a representative population of optimised QS
configurations. The density for each number of field period (color) is stacked
vertically on top of one another, and represents the number of configurations
in the database satisfying those parameters. The rightmost plot shows the same
analysis through a QI near-axis database (Plunk, 2024). This shows the rough
relative ordering of $q_{\mathrm{eff}}$ between different omnigeneous fields,
as indicated in the text. Most QH configurations are $N=4$, and their
$q_{\mathrm{eff}}$ is the lowest for all $N$, while larger or smaller $N$ lead
roughly to larger $q_{\mathrm{eff}}$. This shows the complexity and detail of
The $N=2$ is the main QA.

To strengthen and illustrate this behaviour of $q_{\mathrm{eff}}$ across
different configurations, we use the large database of near-axis QS
configurations of Landreman (2022) and near-axis QI configurations of Plunk
(2024) to evaluate this parameter across configurations. This confirms that
one expects the residual to be smallest in tokamaks and QAs, with the small-
mirror regime barely becoming relevant near their core. We leave a more
complete analysis of these databases and the lessons to be learned from these
for the future. We also note that more complex field shaping beyond the simple
model used in this paper could change some of the exact quantitative behaviour
observed concerning especially the location of the residual transition, but we
also leave this to future investigations.

##  5 Conclusions

In this paper, we have carefully analysed the behaviour of the residual in the
limit of small mirror ratio. The contribution of barely passing particles
provides a finite residual in this limit, changing its usual scaling and
exchanging roles of the importance between trapped and passing particles. We
identify the role of such barely trapped particles and provide some analytical
estimates, that we compare to some gyrokinetic simulations. This limiting
behaviour, however, is shown to occur at very small mirror ratios
$\Delta<(\omega_{d}/\omega_{t})^{2}$, where $\omega_{d}$ is the radial drift
frequency and $\omega_{t}$ the transit frequency of a thermal particle to
travel a connection length. An analysis using near-axis theory of this effect
through tokamaks, quasisymmetric and quasi-isodynamic stellarators suggests
that although barely, the centre of quasi-axisymmetric stellarators is the
region in which some of these effects could manifest most clearly. This
analysis also shows (including a cross-check through a large database of
configurations) that the residual itself tends to be larger in quasi-
isodynamic stellarators, to be followed by quasi-helical and lastly quasi-
axisymmetric (and tokamak) ones.

## Data availability

The data that support the findings of this study are openly available at the
Zenodo repository with DOI/URL 10.5281/zenodo.12805697.

## Acknowledgements

We gratefully acknowledge fruitful discussion with R. Nies and W. Sengupta.

## Funding

E. R. was supported by a grant of the Alexander-von-Humboldt-Stiftung, Bonn,
Germany, through a postdoctoral research fellowship.

## Declaration of interest

The authors report no conflict of interest.

##  Appendix A Additional details on the orbit widths

In this appendix we complete the information about the finite orbit width
provided in Section 2.2, necessary to complete the residual calculation in
Section 2.3.

###  A.1 Passing particles

Let us consider the shape of the orbits described by the barely passing
particles living within the boundary layer defined in Section 2.2.1 (see
Figure 1). To evaluate the residual integrals in Eq. (12) we require
information about the turning points of $\delta$. In particular, besides the
location and value of $\delta$ extrema, the second derivative (Bender &
Orszag, 2013, Sec. 6.5). The second derivative at those points is,

$\delta^{\prime\prime}_{\mathrm{pass}}=\sigma\frac{v}{v_{T}}\frac{\epsilon\pi^{2}}{2}\times\begin{cases}\begin{aligned}
&\frac{1}{\sqrt{\hat{\lambda}}},&\quad(\bar{\ell}=\pm 1)\\\
&-\frac{1}{\sqrt{2\Delta+\hat{\lambda}}},&\quad(\bar{\ell}=0)\end{aligned}\end{cases}$
(46)

where we have used the definition of $\hat{\lambda}$ and $\Delta\ll 1$.

To complete the orbit description, we also need the transit time of passing
particles. In the simplified single well model, this is defined to be the time
taken by a particle to move from $\bar{\ell}=-1$ to $1$. The time can be
expressed (Helander & Sigmar, 2005, Eq. (7.27)) in terms of the elliptic
function $K$ (Olver et al., 2020, Sec. 19)(Abramowitz & Stegun, 1968, Eq.
(16.1.1)),

$\tau_{t}\omega_{t}=\frac{4}{\pi}\frac{v_{T}}{v}\frac{K(\kappa)}{\sqrt{1-\lambda\bar{B}(1-\Delta)}},$
(47)

where $\kappa=2\lambda\Delta/[1/\bar{B}-\lambda(1-\Delta)]$.

###  A.2 Trapped particles

The orbits described by trapped particles are ostensibly different. The
function $\delta(\bar{\ell})$ has a single turning point at the centre of the
orbit, point at which the second derivative is

$\delta_{\mathrm{trap}}^{\prime\prime}(0)\approx\sigma\frac{v}{v_{T}}\frac{\epsilon\pi^{2}}{\sqrt{2\bar{\kappa}\Delta}}.$
(48)

The orbits, unlike those of passing particles, are sharp at, in this case,
bounce points. This is a result of the particles spending longer at these
points, where the radial drift is non-zero. This difference in how particles
spend their time on different parts of their orbit also affects the expression
for the orbit time, here called bounce time (Connor et al., 1983)(Helander &
Sigmar, 2005, Eq. (7.28)),

$\tau_{b}\omega_{t}=\frac{2}{\pi}\frac{v_{T}}{v}\sqrt{\frac{2}{\lambda\bar{B}\Delta}}K(\bar{\kappa}).$
(49)

##  Appendix B Residual in a GAM scenario

In this Appendix we present how the description of geodesic acoustic modes
(GAMs) can be made to align with the finite residual result derived in the
main text. To that end, let us start by re-writing the linearised gyrokinetic
equation in Eq. (1) and dropping the initial condition,

$iv_{\parallel}\partial_{\ell}\hat{g}+(\omega-\tilde{\omega}_{d})\hat{g}-J_{0}F_{0}\omega\frac{q\hat{\phi}}{T}=0.$
(50)

As in the residual calculation, we have written the equation for
$k_{\alpha}=0$, which leads to vanishing of the diamagnetic drive.

Because we are here interested in the GAM dynamics, it is conventional to
specialise to an artificial flat-$B$ field, one in which the sole field
property that varies along the field-line is the curvature drift (i.e.
$k_{\perp}\rho_{i}$ is also constant). Modelling
$\omega_{d}(\ell)=\omega_{d}\cos(\pi\ell/L_{d})$, we may Fourier resolve Eq.
(50) writing $\hat{g}=\sum_{n=-\infty}^{\infty}\hat{g}_{n}e^{in\pi\ell/L_{d}}$
and $\hat{\phi}=\sum_{n=-\infty}^{\infty}\hat{\phi}_{n}e^{in\pi\ell/L_{d}}$.
Taking into account the coupling through $\omega_{d}$, and

$\hat{g}\cos\left(\frac{\pi\ell}{L_{d}}\right)=\frac{1}{2}\sum_{n=-\infty}^{\infty}(\hat{g}_{n+1}+\hat{g}_{n-1})e^{in\pi\ell/L_{d}},$
(51)

we may then write Eq. (50) as,

$\left(-nx_{\parallel}+\frac{\omega}{\omega_{t}}\right)\hat{g}_{n}-\frac{\tilde{\omega}_{d}}{2\omega_{t}}\left(\hat{g}_{n-1}+\hat{g}_{n+1}\right)=F_{0}J_{0}\frac{\omega}{\omega_{t}}\frac{q\hat{\phi}_{n}}{T},$
(52)

where $\omega_{t}=\pi v_{T}/L_{d}$ is the transit frequency over the
characteristic scale of the drift variation and
$x_{\parallel}=v_{\parallel}/v_{T}$.

The system has a sideband coupling through the drift, whose overlap is
controlled by $\omega_{d}/\omega_{t}$. Thus, ordering
$\epsilon=\omega_{d}/\omega_{t}\ll 1$ is particularly convenient to regularise
the problem and be able to truncate it. In fact, if we drive the system
uniformly, meaning we assume $\hat{\phi}_{0},~{}\hat{g}_{0}\sim O(1)$, we
expect to find small sidebands. That way, we may focus on the following
reduced system of equations,

$\displaystyle\left(x_{\parallel}+\frac{\omega}{\omega_{t}}\right)\hat{g}_{-1}-\frac{\tilde{\omega}_{d}}{2\omega_{t}}\hat{g}_{0}\approx
F_{0}J_{0}\frac{\omega}{\omega_{t}}\frac{q\hat{\phi}_{-1}}{T},$ (53a)
$\displaystyle\frac{\omega}{\omega_{t}}\hat{g}_{0}-\frac{\tilde{\omega}_{d}}{2\omega_{t}}(\hat{g}_{-1}+\hat{g}_{1})\approx
F_{0}J_{0}\frac{\omega}{\omega_{t}}\frac{q\hat{\phi}_{0}}{T},$ (53b)
$\displaystyle-\left(x_{\parallel}-\frac{\omega}{\omega_{t}}\right)\hat{g}_{1}-\frac{\tilde{\omega}_{d}}{2\omega_{t}}\hat{g}_{0}\approx
F_{0}J_{0}\frac{\omega}{\omega_{t}}\frac{q\hat{\phi}_{1}}{T}.$ (53c)

In addition to the gyrokinetic equation written in this form, we must complete
the eigenvalue problem with the quasineutrality condition. The condition, now
explicitly involving electrons ($e$) and ions ($i$), reads in this basis,

$\frac{T_{i}}{q_{i}}\sum_{s=e,i}\int
J_{0s}\hat{g}_{s,k}\mathrm{d}^{3}\mathbf{v}=n(1+\tau)\hat{\phi}_{k},$ (54)

where the sum is over both ions and electrons. To construct the final form of
the dispersion we shall eventually use $b_{e}/b_{i}\sim m_{e}/m_{i}\ll 1$,
$\zeta_{e}/\zeta_{i}\sim\sqrt{m_{e}/m_{i}}\ll 1$ and
$\epsilon_{e}/\epsilon_{i}\sim\sqrt{m_{i}/m_{e}}$.

###  B.1 GAM dispersion

The common form of the dispersion relation for GAMs is obtained by combining
the equations in Eqs. (53) to write $\hat{g}_{0}$ explicitly as function of
$\hat{\phi}_{0}$ to leading order in $O(\epsilon^{2})$ and performing the
appropriate velocity space integrals. The result (Gao et al., 2006, 2008;
Sugama & Watanabe, 2006),

$\mathcal{D}=1-\Gamma_{0}(b)+\frac{\epsilon^{2}}{2}\left[\mathcal{D}^{(2)}-\frac{(\mathcal{D}^{(1)})^{2}}{1+\tau+\mathcal{D}^{(0)}}\right],$
(55)

where,

$\displaystyle\mathcal{D}^{(2)}=~{}\frac{1}{\zeta}\left[\Gamma_{0}(b)\frac{\zeta}{2}\left(1+2\zeta^{2}(1+\zeta
Z(\zeta))\right)+F_{2}(b)\zeta(1+\zeta
Z(\zeta))+\frac{1}{4}F_{4}(b)Z(\zeta)\right],$ (56a)
$\displaystyle\mathcal{D}^{(1)}=~{}\Gamma_{0}(b)\zeta(1+\zeta
Z(\zeta))+\frac{1}{2}F_{2}(b)Z(\zeta),$ (56b)
$\displaystyle\mathcal{D}^{(0)}=\Gamma_{0}(b)\zeta Z(\zeta),$ (56c) $\int
F_{0}J_{0}^{2}\mathrm{d}^{3}\mathbf{v}=\Gamma_{0}(b),$ (56d)

and $\zeta=\omega/\omega_{t}$. The dispersion relation is consistent with
multiple modes, which have been explored in Gao et al. (2008). Note that in
those pieces of work (Gao et al., 2006, 2008; Sugama & Watanabe, 2006), the
problem is solved not using a Fourier resolution of the problem like we have
here, but instead using the integrating factor approach of Connor et al.
(1980).

The dispersion relation in Eq. (55) can be assessed near $\zeta\rightarrow 0$,
which is responsible for the long time response of the plasma (Schiff, 2013,
Theorem 2.36). It may be shown by expanding the dispersion function (Fried &
Conte, 2015), and taking for simplicity the small finite Larmor radius limit,

$\mathcal{D}\approx\frac{b}{\omega}\left[1+\frac{\epsilon^{2}}{4b}\left(1+\frac{\pi}{2(1+\tau)}\right)\right]\left(\omega-\omega_{0}\right),$
(57)

where

$\frac{\omega_{0}}{\omega_{t}}=-i\frac{\sqrt{\pi}}{2}\left[\frac{2b}{\epsilon^{2}}+\left(\frac{1}{2}+\frac{\pi}{4(1+\tau)}\right)\right]^{-1}.$
(58)

The system shows a purely damped mode, but no truly net residual.

###  B.2 Revival of the residual

This no residual conclusion is not consistent with the calculation in this
paper. So, where is the residual hiding? To see how the approach to the GAM
could have missed the residual contribution, let us go back to the truncated
system of equations where the $n=0,~{}\pm 1$ modes are retained, Eqs. (53),
and recombine them into

$\frac{T/q}{F_{0}J_{0}}g_{\pm}=\frac{1}{2}\frac{(\tilde{\omega}_{d}/\omega_{t})^{2}\mp
4\zeta(x_{\parallel}\pm\zeta)}{(\tilde{\omega}_{d}/\omega_{t})^{2}+2(x_{\parallel}^{2}-\zeta^{2})}\phi_{\pm}\mp\frac{\tilde{\omega}_{d}}{\omega_{t}}\frac{x_{\parallel}\pm\zeta}{(\tilde{\omega}_{d}/\omega_{t})^{2}+2(x_{\parallel}^{2}-\zeta^{2})}\phi_{0}\\\
-\frac{1}{2}\left(\frac{\tilde{\omega}_{d}}{\omega_{t}}\right)^{2}\frac{1}{(\tilde{\omega}_{d}/\omega_{t})^{2}+2(x_{\parallel}^{2}-\zeta^{2})}\phi_{\mp},$
(59a)
$\frac{T/q}{F_{0}J_{0}}g_{0}=\frac{2(x_{\parallel}^{2}-\zeta^{2})}{(\tilde{\omega}_{d}/\omega_{t})^{2}+2(x_{\parallel}^{2}-\zeta^{2})}\phi_{0}+\frac{\tilde{\omega}_{d}}{\omega_{t}}\frac{x_{\parallel}-\zeta}{(\tilde{\omega}_{d}/\omega_{t})^{2}+2(x_{\parallel}^{2}-\zeta^{2})}\phi_{-}\\\
-\frac{\tilde{\omega}_{d}}{\omega_{t}}\frac{x_{\parallel}+\zeta}{(\tilde{\omega}_{d}/\omega_{t})^{2}+2(x_{\parallel}^{2}-\zeta^{2})}\phi_{+},$
(59b)

where $\pm$ denote the $n=\pm 1$ sidebands. We did not use this full form of
the equations when deriving the dispersion relation for the GAMs, but instead
their limit when $\epsilon=\omega_{d}/\omega\ll 1$. Formally, this ordering
was used to expand the kinetic resonant denominators

$\mathcal{R}=\frac{1}{\tilde{\omega}_{d}^{2}/\omega_{t}^{2}+2(x_{\parallel}^{2}-\zeta^{2})},$
(60)

that are found ubiquitous in Eqs. (59b). For this expansion in the denominator
to be sound we must have, of course,
$x_{\parallel}^{2}-\zeta^{2}\gg\tilde{\omega}_{d}^{2}/\omega_{t}^{2}$, where
we shall not forget the velocity space dependence of
$\tilde{\omega}_{d}=\omega_{d}(x_{\parallel}^{2}+x_{\perp}^{2}/2)$. The GAM
dispersion relation thus fails to describe any physics where
$x_{\parallel}^{2}-\zeta^{2}\ll\epsilon^{2}x_{\perp}^{4}/4$. This is
especially problematic at long time scales (i.e. within a layer in
$\omega$-space where $\omega<\omega_{d}$) and for the part of the population
living within a narrow layer of order $x_{\parallel}\sim\epsilon$ in velocity
space near $x_{\parallel}=0$. I.e. the GAM description overlooks the
contribution from barely passing particles, whose transit time is
significantly longer than that of the bulk.

The question is then, how can one capture the behaviour from within this layer
properly in this GAM formalism? Can one recover a residual result like that in
Eq. (37)? To do so we must not expand in small $\tilde{\omega}_{d}$, but
instead do so in $\zeta\rightarrow 0^{+}$ (indicating approach from the
positive $\Im\\{\omega\\}$ direction). With this in mind, let us write the
quasineutrality condition applied to Eq. (59b) as

$(1+\tau-\mathcal{D}^{(2)})\hat{\phi}(0)\approx-\frac{\epsilon}{2}\left[\mathcal{D}^{(1)}_{-}\hat{\phi}(-1)-\mathcal{D}^{(1)}_{+}\hat{\phi}(1)\right],$
(61)

where

$\displaystyle\mathcal{D}^{(2)}=$ $\displaystyle~{}\frac{2}{\bar{n}}\int
F_{0}J_{0}^{2}(x_{\parallel}^{2}-\zeta^{2})\mathcal{R}\mathrm{d}^{3}\mathbf{v}$
(62) $\displaystyle\mathcal{D}^{(1)}_{\pm}=-\frac{2}{\bar{n}}\int
F_{0}J_{0}^{2}\left(x_{\parallel}^{2}+\frac{x_{\perp}^{2}}{2}\right)(x_{\parallel}\pm\zeta)\mathcal{R}\mathrm{d}^{3}\mathbf{v}.$
(63)

To evaluate these integrals, we rewrite $\mathcal{R}$ by separating it into a
sum over simple poles. To do so, we define,

$\Delta=~{}\sqrt{\frac{1}{\epsilon^{2}}+x_{\perp}^{2}+2\zeta^{2}},\quad\zeta_{\pm}=~{}\frac{\Delta}{\epsilon}\pm\left(\frac{1}{\epsilon^{2}}+\frac{x_{\perp}^{2}}{2}\right),$
(64)

so that

$\mathcal{R}=-\frac{1}{2\epsilon\Delta}\left[\frac{1}{x_{\parallel}^{2}+\zeta_{+}}-\frac{1}{x_{\parallel}^{2}-\zeta_{-}}\right].$
(65)

Choosing the negative branch of the square root for a correct continuation
from $\Im\\{\zeta\\}>0$ to the rest of the complex plane,

$\frac{1}{x_{\parallel}^{2}\pm\zeta_{\pm}}=\frac{1}{2\sqrt{\mp\zeta_{\pm}}}\left(\frac{1}{x_{\parallel}-\sqrt{\mp\zeta_{\pm}}}-\frac{1}{x_{\parallel}+\sqrt{\mp\zeta_{\pm}}}\right),$
(66)

in such a way that the integrals Eqs. (62)-(63) explicitly involve integrals
over $x_{\parallel}$. This form of $\mathcal{R}$ allows us to express
integrals in terms of plasma dispersion functions (Fried & Conte, 2015) upon
appropriate redefinition of the sign of $x_{\parallel}$ (which will annihilate
the contribution from odd $x_{\parallel}$ terms).555We shall here not be
extremely careful with the definition of branch cuts and the precise
deformation of the Laplace contour in $\zeta$-space. This would be needed for
a fuller description of the time response of the system (one that captures the
contribution from branch cuts for example), but here we content ourselves with
the $\zeta\rightarrow 0$ response. As a result, we may write the integrals as
a combination of

$\displaystyle I_{nm}=$ $\displaystyle~{}\frac{1}{\bar{n}}\int
x_{\parallel}^{2n}x_{\perp}^{2m}F_{0}J_{0}^{2}\mathcal{R}\mathrm{d}^{3}\mathbf{v}=-\frac{1}{\epsilon}\int_{0}^{\infty}x_{\perp}^{2m+1}J_{0}^{2}e^{-x_{\perp}^{2}}\frac{1}{\Delta}\left[\frac{Z_{n}(\sqrt{-\zeta_{+}})}{\sqrt{-\zeta_{+}}}-\frac{Z_{n}(\sqrt{\zeta_{-}})}{\sqrt{\zeta_{-}}}\right]\mathrm{d}x_{\perp},$
(67)

where we define,

$Z_{n}(x)=\frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty}\frac{x_{\parallel}^{2n}e^{-x_{\parallel}^{2}}}{x_{\parallel}-x}\mathrm{d}x_{\parallel},$
(68)

for $\Im\\{x\\}>0$, and analytically continued to the rest of the complex
plane. In particular, we may write

$\displaystyle\mathcal{D}_{\pm}^{(1)}=\mp
2\zeta\left(I_{10}+\frac{I_{01}}{2}\right),$ (69a)
$\displaystyle\mathcal{D}^{(2)}=2\left(I_{10}-\zeta^{2}I_{00}\right).$ (69b)

These integrals remain quite sophisticated, and simplifying them is paramount
to analytically proceed forward. A natural simplifying attempt is to use
asymptotic forms of the plasma dispersion function (Fried & Conte, 2015). The
argument $\zeta_{+}\approx
2/\epsilon^{2}+x_{\perp}^{2}-\epsilon^{2}x_{\perp}^{4}/8$, which is a large
and positive real part quantity owing to the largeness of $1/\epsilon^{2}$, we
may use the asymptotic form (Fried & Conte, 2015, Sec. IID)
$Z(x)\approx-\sum_{n=0}^{\infty}x^{-(2n+1)}(n-1/2)!/\sqrt{\pi}$ (the
exponential term is exponentially small). In the case of
$\zeta_{-}\approx\zeta^{2}-x_{\perp}^{4}\epsilon^{2}/8$ and we may consider an
expansion in this small argument. Namely, (Fried & Conte, 2015, Sec. IIC)
$Z(x)=i\sqrt{\pi}\exp(-x^{2})-x\sum_{n=0}^{\infty}(-x^{2})^{n}\sqrt{\pi}/(n+1/2)!$.
This introduces a leading order non-zero imaginary contribution.

With the above tools in place, we may proceed and compute the required
integrals to the necessary order.

####  B.2.1 Integrals for $\mathcal{D}^{(2)}$

Let us compute first the leading order $I_{00}$. Without having to go into the
complex details about the specific branch cuts and complex quadrant of $\zeta$
in the complex plane, one can show (Gradshteyn & Ryzhik, 2014, Eq. 3.387.7)

$I_{00}\approx~{}\int_{0}^{\infty}x_{\perp}J_{0}^{2}e^{-x_{\perp}^{2}}\frac{\sqrt{\pi}}{\sqrt{\frac{x_{\perp}^{4}\epsilon^{2}}{8}-\zeta^{2}}}\mathrm{d}x_{\perp}[1+O(\zeta,\epsilon^{2})]\propto\frac{1}{\epsilon}\ln\left(\frac{\epsilon}{\zeta\sqrt{2}}\right),$
(70)

where for this estimate we have assumed $b\ll 1$ to approximate $J_{0}\sim 1$
and we have kept the leading order term in $\zeta$ (in the limit of small
$\zeta$). So, in the limit of $\zeta\rightarrow 0$, this integral diverges
logarithmically, but its contribution to $\mathcal{D}^{(2)}$ vanishes, Eq.
(69b).

Computing then $I_{10}$, and using $Z_{1}(x)=x[1+xZ(x)]$,

$\displaystyle I_{10}\approx$
$\displaystyle~{}-\int_{0}^{\infty}x_{\perp}J_{0}^{2}e^{-x_{\perp}^{2}}\left[-1+\frac{\epsilon}{2}\sqrt{\frac{\pi}{2}}x_{\perp}^{2}+\frac{\epsilon^{2}}{4}(1+2x_{\perp}^{2}-x_{\perp}^{4})+O(\epsilon^{3})\right]\mathrm{d}x_{\perp}$
(71) $\displaystyle\approx$
$\displaystyle~{}\frac{1}{2}\left[\Gamma_{0}(b)-\frac{\epsilon}{2}\sqrt{\frac{\pi}{2}}F_{2}(b)-\frac{\epsilon^{2}}{4}(\Gamma_{0}(b)+2F_{2}(b)-F_{4}(b))\right]$
(72) $\displaystyle\approx$
$\displaystyle~{}-\frac{1}{2}\left(b-1+\frac{\epsilon}{2}\sqrt{\frac{\pi}{2}}+\frac{\epsilon^{2}}{4}\right)=\frac{1}{2}\mathcal{D}^{(2)},$
(73)

where we used the relevant Weber integrals (Gradshteyn & Ryzhik, 2014, Eq.
6.615) and the notation
$F_{n}=2\int_{0}^{\infty}x^{n+1}e^{-x^{2}}J_{0}^{2}(x\sqrt{2b})\mathrm{d}x$,
and in the last line considered the small $b$ limit. Importantly, there is a
term linear in $\epsilon$ which comes from the pole contribution to the plasma
dispersion function.

####  B.2.2 Integrals for $\mathcal{D}_{\pm}^{(1)}$

With $\mathcal{D}^{(2)}$ constructed, we may turn to $\mathcal{D}^{(1)}$, Eq.
(69a). The integral has an overall factor of $\zeta$, and thus to leading
order, it will vanish unless there is some $\zeta$-divergence. The term
$I_{10}$, which we have just computed, does not have such divergence, and thus
its contribution will vanish. So we only need to calculate $I_{01}$, which one
may show to be $I_{01}\approx\sqrt{2}/\epsilon$ to leading order. Thus,
$\mathcal{D}^{(1)}_{\pm}\sim O(\zeta)$, and thus it will vanish in the small
$\zeta$ limit. One may savely drop the coupling terms in Eq. (61) (the
sideband $\phi_{\pm}$ does not have any divergent behaviour neither).

####  B.2.3 Dispersion relation

Thus, the remaining dispersion function is,

$\mathcal{D}=1-\left[\Gamma_{0}(b)-\frac{\epsilon}{2}\sqrt{\frac{\pi}{2}}F_{2}(b)-\frac{\epsilon^{2}}{4}(\Gamma_{0}(b)+2F_{2}(b)-F_{4}(b))\right],$
(74)

where we have summed over species and taken the limit of $m_{e}/m_{i}\ll 1$,
and all quantities here should now be considered to represent ions. The value
of the residual can then be written666We are being loose here about initial
condition, but we may simply consider the RH initial condition of a uniformly
perturbed potential., assuming $b\ll 1$ for simplicity,

$\frac{\phi(\infty)}{\phi(0)}=\frac{1}{1+\frac{\epsilon}{2b}\sqrt{\frac{\pi}{2}}+\frac{\epsilon^{2}}{4b}}.$
(75)

It includes the leading order linear term in $\epsilon$, as the residual
expression in the main text does. The difference with the result in the main
text is the numerical factor in front of the linear term. As opposed to the
$0.26\omega_{d}/(b\omega_{t})$ obtained in the text, and realising that
$\omega_{t}$ as used in this appendix is $\pi$ times that in the main text,
the result here yields $(1/2\sqrt{2\pi})\omega_{d}/(b\omega_{t})\approx
0.20\omega_{d}/(b\omega_{t})$. This is a 30% discrepancy between both
estimates of the residual, but the same scaling nonetheless.

##  Appendix C Near-axis properties in optimised configurations

In this Appendix we present the near-axis calculations necessary to obtain the
expressions in Table 1 for the residual relevant parameters in different
omnigeneous magnetic fields. These should be taken as informed estimates for
the amplitudes of the simple model assumed in the main text. As we shall show,
this is a good fit for QS fields, but not so much for QI. We assume some basic
understanding of inverse-coordinate near-axis theory (Garren & Boozer, 1991b,
a), and shall not derive the basic building elements of it. We refer the
reader to the work by Landreman & Sengupta (2019) for the general equations
for magnetohydrostatic equilibrium and in particular in a quasisymmetric
configuration, and Plunk et al. (2019); Rodríguez & Plunk (2023) for quasi-
isodynamic ones. We shall here use, with further explicit reference to those
works, the elements needed for the evaluation of the appropriate quantities.

###  C.1 Quasisymmetric fields

Let us start by writing the magnetic field magnitude near the axis for a
quasisymmetric field (Garren & Boozer, 1991a, Eq. (A1)) (Landreman & Sengupta,
2019, Eq. (2.15)),

$B\approx B_{0}(1+r\eta\cos\chi),$ (76)

where $r=\sqrt{2\psi/\bar{B}}$ is a pseudo-radial coordinate normalised to a
reference $\bar{B}$, and $\chi=\theta-N\varphi$, where $N$ is the direction of
symmetry of the QS field and we are using Boozer coordinates. Because $B_{0}$
is a constant, it is clear from this form that the constant parameter $\eta$
measures the variation of the magnetic field within a surface (to leading
order). Thus, along a field line (at constant $\alpha$) the magnetic field
depends on $\chi=\alpha+\bar{\iota}\varphi$, and thus the mirror ratio is,

$\Delta=r\eta,$ (77)

as indicated in Table 1.

We now need to construct the other input important to the residual calculation
which is,

$q_{\mathrm{eff}}=\frac{1}{\pi}\frac{1}{k_{\perp}\rho_{i}}\frac{\omega_{d}}{\omega_{t}},$
(78)

whose definition is meant to take the place of $q$ in the RH residual. See the
main text, Section 4, for more details, including its connections to banana
widths (roughly $\sim\rho_{i}q_{\mathrm{eff}}/\sqrt{\Delta}$) and the
transition between the low-mirror and RH residual regimes.

Let us start by finding the amplitude of the drift frequency
$\omega_{d}(\chi)$. The curvature drift is by definition,

$\omega_{d}(\chi)=-v_{T}\frac{\mathbf{B}\times\nabla
B\cdot\nabla\psi}{B^{3}}\bar{B}k_{\psi}\rho_{i},$ (79)

where we have defined the ion Larmor radius $\rho_{i}=m_{i}v_{T}/q_{i}\bar{B}$
with respect to some reference field $\bar{B}$. The triple vector product may
be directly computed using the contravariant Boozer coordinate basis in the
near-axis framework (Jorge & Landreman, 2020, Eq. (45))777The expression in
Jorge & Landreman (2020) has an incorrect additional factor of $B_{0}$, as can
be checked dimensionally. This typo is unimportant., which yields

$\omega_{d}(\chi)=-v_{T}B_{0}r\eta k_{\psi}\rho_{i}\sin\chi+O(r^{2}).$ (80)

The coefficient $\omega_{d}$ may be directly read-off from the amplitude of
this expression. Note here that $\eta$ plays a primary role in controlling the
magnitude of the radial drift, as it controls the magnitude of the magnetic
field magnitude gradients.

To make sense of the typical magnitude of $\eta$, it is convenient to
introduce the description of flux surface shapes in the near-axis framework.
Flux surfaces are defined as a function of Boozer coordinates with respect to
the magnetic axis, $\mathbf{r}_{0}$, in the Frenet-Serret basis
$\\{\hat{b},\hat{\kappa},\hat{\tau}\\}$ (tangent, normal and binormal) of the
latter, so that
$\mathbf{r}(\psi,\theta,\varphi)-\mathbf{r}_{0}=X\hat{\kappa}+Y\hat{\tau}+Z\hat{b}$.
Thus $X$ is a function that gives the distance from flux surfaces to the axis
along the normal to the latter. To leading order this is proportional to
$X_{1}=r\eta/\kappa$, while along the binormal it scales like
$Y_{1}\sim\kappa/\eta$ (Landreman & Sengupta, 2019, Eq. (2.13)). Thus, in
order to avoid extreme shaping $\eta\sim\kappa$ (Rodríguez et al., 2023). As
$\kappa$ is generally a function of the toroidal angle and $\eta$ is not, the
shaping of flux surfaces will change toroidally, but one may take the
curvature as a scale for $\eta$. In the case of a circular cross section
tokamak one may show that $\eta=1/R$. This relation between the variation of
the magnetic field and the curvature of the axis (a field line after all) is a
physical consequence of the relation between the bending field lines and
magnetic pressure.

We now need to find an expression for the transit time
$\omega_{t}=v_{T}/L_{d}$, where $L_{d}$ is the connection length; the distance
from the trough to the top of the well. We thus need to compute $\ell$, the
distance along the field line. In quasisymmetry the length is simply a
rescaled form of the Boozer toroidal angle $\varphi$, so that (Landreman &
Sengupta, 2019, Eq. (A20))

$\frac{\mathrm{d}\chi}{\mathrm{d}\ell}\approx\frac{\bar{\iota}}{R_{\mathrm{ax}}},$
(81)

where $R_{\mathrm{ax}}=L_{\mathrm{ax}}/2\pi$ and $L_{\mathrm{ax}}$ is the
length of the magnetic axis, and $\bar{\iota}=\iota-N$. Given that in Eq. (76)
the magnetic field has a well of halfwidth $\pi$, then $L_{d}\approx\pi
R_{\mathrm{ax}}/\bar{\iota}$ and,

$\omega_{t}=\bar{\iota}\frac{v_{T}}{\pi R_{\mathrm{ax}}}.$ (82)

Finally, let us consider the normalized perpendicular wavenumber
$(k_{\perp}\rho_{i})^{2}=\langle|\nabla\psi|^{2}\rangle(k_{\psi}\rho_{i})^{2}$.
Note how we are using an averaged form of the flux expansion, which makes the
FLR parameter constant, as assumed in our model construction. The particular
form of $k_{\perp}\rho_{i}$ is motivated by the involvement of
$b=(k_{\perp}\rho_{i})^{2}/2$ in the residual, where it appears flux surface
averaged (Plunk & Helander, 2024) (including variation along the line would be
straightforward). We need $|\nabla\psi|^{2}$ from the near-axis description of
the field; using the contravariant basis once again (Jorge & Landreman, 2020,
Eq. (41)),

$|\nabla\psi|^{2}\approx
r^{2}\left(B_{0}\frac{\kappa}{\eta}\right)^{2}\left[\left(\frac{\eta}{\kappa}\right)^{4}\sin^{2}\chi+\left(\cos\chi-\sigma\sin\chi\right)^{2}\right],$
(83)

where $\sigma$ is a function of the toroidal angle $\varphi$, result of
solving a non-linear Riccati equation (Garren & Boozer, 1991a; Landreman &
Sengupta, 2019). The flux surface average of this expression can be carried
out straightforwardly, using to leading order
$\langle\dots\rangle\approx\int\mathrm{d}\chi\mathrm{d}\varphi\dots/(4\pi^{2})$,

$\left\langle|\nabla\psi|^{2}\right\rangle\approx(rB_{0}\hat{\mathcal{G}})^{2},$
(84)

where,

$\hat{\mathcal{G}}^{2}=\frac{1}{4\pi}\int_{0}^{2\pi}\left(\frac{\kappa}{\eta}\right)^{2}\left(1+\sigma^{2}+\frac{\eta^{4}}{\kappa^{4}}\right)\mathrm{d}\varphi.$
(85)

The involvement of $\sigma$ makes this geometric quantity rather obscure. In
fact $\sigma$ is directly related to the shaping of flux surfaces as
$Y_{1}=(\kappa/\eta)(\sin\chi+\sigma\cos\chi)$ (Landreman & Sengupta, 2019,
Eq. (2.13)), but its interpretation in simple terms is difficult (Rodríguez,
2023). Although it may be understood roughly as a measure of the rotation of
the elliptical cross-sections near the axis respect to the Frenet-Serret frame
(Rodríguez, 2023, Eq. (B4a)), it also affects the elongation of flux surfaces.
It would be beneficial in the discussion, thus, to provide a more direct
geometric interpretation to $\hat{\mathcal{G}}$. We do so using (Rodríguez,
2023, Eq. (3.2a)) to write,

$\hat{\mathcal{G}}^{2}=\frac{1}{2\pi}\int_{0}^{2\pi}\frac{1}{\sin
2e}\mathrm{d}\varphi$ (86)

where $\mathcal{E}=\tan e$ and $\mathcal{E}$ is the elongation of the flux
surfaces in the plane normal to the axis as a function of $\varphi$. The angle
$e\in(0,\pi/2)$ may be interpreted as the angle subtended by a right-angle
triangle with the major and minor axes of the ellipse as catheti. Thus the
geometric factor $\hat{\mathcal{G}}$ is a direct measure of the flux surface
elongation. A value of $\hat{\mathcal{G}}=1$ corresponds to all cross-section
being circular, any amount of shaping leading to $\hat{\mathcal{G}}>1$.

Putting everything together into $q_{\mathrm{eff}}$,

$q_{\mathrm{eff}}=\frac{1}{\iota-N}\frac{\eta
R_{\mathrm{ax}}}{\hat{\mathcal{G}}}.$ (87)

####  C.1.1 Tokamak limit

The case of the axisymmetric tokamak is a particularly simple limit of this.
Considering the limit of $\kappa\rightarrow 1/R$, where $R$ is the major
radius, then $R_{\mathrm{ax}}\rightarrow R$ and all quantities become
$\varphi$-independent. Then, we may write $q_{\mathrm{eff}}=q(\eta
R)/\hat{\mathcal{G}}_{\mathrm{tok}}$, where $q=1/\iota$ is the safety factor
and, $\hat{\mathcal{G}}_{\mathrm{tok}}^{2}=1/\sin 2e$. If we then consider a
circular cross-section tokamak (where $e=\pi/4$), then $\eta=1/R$,
$\hat{\mathcal{G}}=1$, and thus $q_{\mathrm{eff}}=q$. This is why we have
defined $q_{\mathrm{eff}}$ the way we have. As a reference
$\hat{\mathcal{G}}=2$ corresponds to $e=\pi/8$ and thus an elongation
$\mathcal{E}\approx 0.4$.

###  C.2 Quasi-isodynamic fields

Let us write the magnetic field of an exactly omnigeneous, QI, stellarator-
symmetric field near the axis (Plunk et al., 2019, Eq. (6.1)) (Rodríguez &
Plunk, 2023, Eqs. (8-9a)),

$B=B_{0}(\varphi)\left[1-rd(\varphi)\sin\alpha+O(r^{2})\right],$ (88)

where $B_{0}(\varphi)$ and $d(\varphi)$ are even and odd functions of
$\varphi$ respectively. The latter is required for the fulfilment of
omnigeneity. Note that $B$ is here an explicit function of $\alpha$, which
unless the rotational transform is integer, makes $B$ a non-periodic function.
This is the well-known impossibility of achieving omnigeneity exactly to
leading order near the axis with poloidal $|\mathbf{B}|$ contours (Plunk et
al., 2019). Acknowledging that in practice omnigeneity will have to be broken
in some buffer region near the tops (Plunk et al., 2019; Camacho Mata et al.,
2022), we shall consider Eq. (88) as given.

Let us now consider a simple model for the magnetic field on axis,

$B_{0}(\varphi)=\bar{B}\left(1-\Delta\cos N_{\mathrm{nfp}}\varphi\right),$
(89)

where $\Delta$ is the mirror ratio and $N_{\mathrm{nfp}}$ is the number of
field periods (the toroidal $N_{\mathrm{nfp}}$-fold symmetry). Unlike in the
QS scenario, the control of the on-axis magnetic field in a QI configuration
gives complete control of the mirror ratio.

The choice of this form of $B_{0}$ requires the curvature to have vanishing
points at $\varphi=n\pi/N_{\mathrm{nfp}}$ for $n\in\mathbb{Z}$, and non-
vanishing first derivative (often referred to as a first order zero). Not
doing so would lead to the loss of trapped particles as discussed in detail in
Rodríguez & Plunk (2023). As a result, the variation in the field $d(\varphi)$
must also share those zeroes with $\kappa$ to avoid extreme shaping (the
leading order shaping is analogous to the QS scenario). For now, let us keep
it general and construct the necessary coefficients as we did with the QS
case. Starting off the drift, and using (Jorge & Landreman, 2020, Eq. (37)),

$\omega_{d}(\theta)\approx-rv_{T}\bar{B}\kappa\left(X_{1c}\sin\theta-
X_{1s}\cos\theta\right)k_{\psi}\rho_{i},$ (90)

where $X_{1c}$ and $X_{1s}$ are the cosine and sine $\theta$-harmonics of
$X_{1}$ to leading order. Following their definition in terms of $B$
(Landreman & Sengupta, 2019, Eq. (A22)), and using the expression for $B$ in
Eq. (88), for an exactly omnigeneous field,

$\displaystyle X_{1c}=$ $\displaystyle\frac{d}{\kappa}\sin\iota\varphi,$ (91a)
$\displaystyle X_{1s}=$ $\displaystyle-\frac{d}{\kappa}\cos\iota\varphi,$
(91b)

so that Eq. (90) reduces to,

$\omega_{d}(\varphi)=-rv_{T}\bar{B}d(\varphi)k_{\psi}\rho_{i}\cos\alpha.$ (92)

We need the amplitude of this function to feed into $q_{\mathrm{eff}}$, Of
course, generally the shape of this function will not be that of a simple sine
as in the QS case. However, we may choose the simple form,

$d(\varphi)=\bar{d}\sin(N_{\mathrm{nfp}}\varphi),$ (93)

to give an amplitude $\omega_{d}\approx rv_{T}\bar{d}\bar{B}\cos\alpha$. Note
a significant difference with respect to the QS case, which is the explicit
$\alpha$ dependence. The amplitude of the field varies from field-line to
field-line. We have lost the field-line equivalence (Boozer, 1983b; Helander,
2014; Rodriguez et al., 2020) of quasisymmetry. To treat this difference
consistently within the residual treatment we would have to treat more
carefully the variation of the field over the surface. However, for a rough
estimate of the drift amplitude, let us keep it as is for now.

Let us now consider $|\nabla\psi|^{2}$ (Jorge & Landreman, 2020, Eq. (33)),

$|\nabla\psi|^{2}=r^{2}B_{0}^{2}\left[\left(X_{1c}\sin\theta-
X_{1s}\cos\theta\right)^{2}+\left(Y_{1c}\sin\theta-
Y_{1s}\cos\theta\right)^{2}\right],$ (94)

where for our ideal omnigeneneous field (Landreman & Sengupta, 2019, Eq.
(A25)),

$\displaystyle Y_{1c}=$
$\displaystyle\frac{\bar{B}}{B_{0}}\frac{\kappa}{d}\left(\cos\iota\varphi+\sigma\sin\iota\varphi\right),$
(95a) $\displaystyle Y_{1s}=$
$\displaystyle-\frac{\bar{B}}{B_{0}}\frac{\kappa}{d}\left(\sigma\cos\iota\varphi-\sin\iota\varphi\right).$
(95b)

Therefore,

$|\nabla\psi|^{2}\approx
r^{2}B_{0}^{2}\left[\left(\frac{d}{\kappa}\right)^{2}\cos^{2}\alpha+\left(\frac{\kappa}{d}\frac{\bar{B}}{B_{0}}\right)^{2}\left(\sin\alpha+\sigma\cos\alpha\right)^{2}\right].$
(96)

Assuming $\Delta\ll 1$ to simplify the flux surface averages and approximate
$B_{0}\approx\bar{B}$, integrating over $\alpha$ and $\varphi$,

$\left\langle|\nabla\psi|\right\rangle\approx\left(r\bar{B}\hat{\mathcal{G}}_{\mathrm{QI}}\right)^{2},$
(97)

where,

$\hat{\mathcal{G}}_{\mathrm{QI}}^{2}=\frac{1}{4\pi}\int_{0}^{2\pi}\left(\frac{\kappa}{d}\right)^{2}\left(1+\sigma^{2}+\frac{d^{4}}{\kappa^{4}}\right)\mathrm{d}\varphi.$
(98)

Note the similarity of this expression to the QS geometric factor Eq. (85). In
fact, Eq. (98) is exactly equivalent to Eq. (86), the expression in terms of
the elongation of flux surfaces in the plane normal to the magnetic axis.

Finally we compute the connection length, which under the approximation of
$\Delta\ll 1$ we may write as $L_{d}\approx\pi
r_{\mathrm{ax}}/N_{\mathrm{nfp}}$. Putting all together,

$q_{\mathrm{eff}}=\frac{1}{N_{\mathrm{nfp}}}\frac{\bar{d}R_{\mathrm{ax}}}{\hat{\mathcal{G}}_{\mathrm{QI}}}\cos\alpha.$
(99)

Note how this parameter changes from field line to field line. The
contribution to the total residual can be thought of as a sum over wells,
where each of these can be thought of separately, thanks to the condition of
omnigeneity. As we move along the field line then, we see different wells,
which assuming this to be the only element that changes from well to well, and
using

$\lim_{N\rightarrow\infty}\frac{1}{N}\sum_{n=0}^{N}|\cos(2\pi\iota
n)|=\frac{1}{2\pi}\int_{0}^{2\pi}|\cos\alpha|\mathrm{d}\alpha=\frac{2}{\pi},$
(100)

by application of Weyl’s lemma (Weyl, 1916, Eq. (2)) for irrational $\iota$,
we may construct an effective parameter $q_{\mathrm{eff}}$,

$q_{\mathrm{eff}}=\frac{1}{N_{\mathrm{nfp}}}\frac{2}{\pi}\frac{\bar{d}R_{\mathrm{ax}}}{\hat{\mathcal{G}}}.$
(101)

We shall not consider here any more sophisticated approach that deals with
these variations more carefully or takes additional differences between wells
into account.

## References

  * Abramowitz & Stegun (1968) Abramowitz, Milton & Stegun, Irene A 1968 Handbook of mathematical functions with formulas, graphs, and mathematical tables, , vol. 55. US Government printing office.
  * Alonso et al. (2017) Alonso, JA, Sánchez, E, Calvo, I, Velasco, JL, McCarthy, KJ, Chmyga, A, Eliseev, LG, Estrada, T, Kleiber, R, Krupnik, LI & others 2017 Observation of oscillatory radial electric field relaxation in a helical plasma. Physical Review Letters 118 (18), 185002.
  * Austin et al. (2019) Austin, Max E, Marinoni, A, Walker, ML, Brookman, MW, Degrassie, JS, Hyatt, AW, McKee, GR, Petty, CC, Rhodes, TL, Smith, SP & others 2019 Achievement of reactor-relevant performance in negative triangularity shape in the diii-d tokamak. Physical review letters 122 (11), 115001.
  * Barnes et al. (2019) Barnes, Michael, Parra, Felix I & Landreman, Matt 2019 stella: An operator-split, implicit–explicit $\delta$f-gyrokinetic code for general magnetic field configurations. Journal of Computational Physics 391, 365–380.
  * Beidler et al. (2021) Beidler, CD, Smith, HM, Alonso, A, Andreeva, T, Baldzuhn, J, Beurskens, MNA, Borchardt, Matthias, Bozhenkov, SA, Brunner, Kai Jakob, Damm, Hannes & others 2021 Demonstration of reduced neoclassical energy transport in wendelstein 7-x. Nature 596 (7871), 221–226.
  * Bender & Orszag (2013) Bender, Carl M & Orszag, Steven A 2013 Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory. Springer Science & Business Media.
  * Bernardin et al. (1986) Bernardin, M. P., Moses, R. W. & Tataronis, J. A. 1986 Isodynamical (omnigenous) equilibrium in symmetrically confined plasma configurations. The Physics of Fluids 29 (8), 2605–2611.
  * Boozer (1983a) Boozer, Allen H. 1983a Transport and isomorphic equilibria. The Physics of Fluids 26 (2), 496–499.
  * Boozer (1983b) Boozer, Allen H 1983b Transport and isomorphic equilibria. The Physics of Fluids 26 (2), 496–499.
  * Boozer (1998) Boozer, Allen H 1998 What is a stellarator? Physics of Plasmas 5 (5), 1647–1655.
  * Camacho Mata et al. (2022) Camacho Mata, K., Plunk, G. G. & Jorge, R. 2022 Direct construction of stellarator-symmetric quasi-isodynamic magnetic configurations. Journal of Plasma Physics 88 (5), 905880503.
  * Cary & Shasharina (1997) Cary, J. R. & Shasharina, S. G. 1997 Omnigenity and quasihelicity in helical plasma confinement systems. Physics of Plasmas 4 (9), 3323–3333, arXiv: https://pubs.aip.org/aip/pop/article-pdf/4/9/3323/12664528/3323_1_online.pdf.
  * Catto et al. (2017) Catto, Peter J, Parra, Felix I & Pusztai, István 2017 Electromagnetic zonal flow residual responses. Journal of Plasma Physics 83 (4), 905830402.
  * Connor et al. (1980) Connor, JW, Hastie, RJ & Taylor, JB 1980 Stability of general plasma equilibria. iii. Plasma Physics 22 (7), 757.
  * Connor et al. (1983) Connor, J. W., Hastie, R. J. & Martin, T. J. 1983 Effect of pressure gradients on the bounce-averaged particle drifts in a tokamak. Nuclear fusion 23 (12), 1702.
  * Connor et al. (1978) Connor, J W, Hastie, R J & Taylor, J B 1978 Phys. Rev. Lett. 40 (6), 396.
  * Conway et al. (2021) Conway, Garrard D, Smolyakov, Andrei I & Ido, Takeshi 2021 Geodesic acoustic modes in magnetic confinement devices. Nuclear Fusion 62 (1), 013001.
  * Diamond et al. (2005) Diamond, Patrick H, Itoh, SI, Itoh, K & Hahm, TS 2005 Zonal flows in plasma—a review. Plasma Physics and Controlled Fusion 47 (5), R35.
  * Fried & Conte (2015) Fried, Burton D & Conte, Samuel D 2015 The plasma dispersion function: the Hilbert transform of the Gaussian. Academic press.
  * Galeev et al. (1969) Galeev, Albert A, Sagdeev, RZ, Furth, HP & Rosenbluth, MN 1969 Plasma diffusion in a toroidal stellarator. Physical Review Letters 22 (11), 511.
  * Gao et al. (2006) Gao, Zhe, Itoh, K, Sanuki, H & Dong, JQ 2006 Multiple eigenmodes of geodesic acoustic mode in collisionless plasmas. Physics of plasmas 13 (10).
  * Gao et al. (2008) Gao, Zhe, Itoh, K, Sanuki, H & Dong, JQ 2008 Eigenmode analysis of geodesic acoustic modes. Physics of Plasmas 15 (7).
  * Garren & Boozer (1991a) Garren, D. A. & Boozer, A. H. 1991a Existence of quasihelically symmetric stellarators. Physics of Fluids B: Plasma Physics 3 (10), 2822–2834.
  * Garren & Boozer (1991b) Garren, D. A. & Boozer, A. H. 1991b Magnetic field strength of toroidal plasma equilibria. Physics of Fluids B: Plasma Physics 3 (10), 2805–2821.
  * Giuliani (2024) Giuliani, Andrew 2024 Direct stellarator coil design using global optimization: application to a comprehensive exploration of quasi-axisymmetric devices. Journal of Plasma Physics 90 (3), 905900303.
  * Goodman et al. (2023) Goodman, A.G., Camacho Mata, K., Henneberg, S.A., Jorge, R., Landreman, M., Plunk, G.G., Smith, H.M., Mackenbach, R.J.J., Beidler, C.D., Helander, P. & et al. 2023 Constructing precisely quasi-isodynamic magnetic fields. Journal of Plasma Physics 89 (5), 905890504.
  * Gradshteyn & Ryzhik (2014) Gradshteyn, Izrail Solomonovich & Ryzhik, Iosif Moiseevich 2014 Table of integrals, series, and products. Academic press.
  * Hall & McNamara (1975a) Hall, L. S. & McNamara, B. 1975a Three-dimensional equilibrium of the anisotropic, finite-pressure guiding-center plasma: Theory of the magnetic plasma. The Physics of Fluids 18 (5), 552–565, arXiv: https://pubs.aip.org/aip/pfl/article-pdf/18/5/552/12317924/552_1_online.pdf.
  * Hall & McNamara (1975b) Hall, Laurence S. & McNamara, Brendan 1975b Three-dimensional equilibrium of the anisotropic, finite-pressure guiding-center plasma: Theory of the magnetic plasma. The Physics of Fluids 18 (5), 552–565.
  * Hazeltine & Meiss (2003) Hazeltine, Richard D & Meiss, James D 2003 Plasma confinement. Courier Corporation.
  * Helander (2014) Helander, P. 2014 Theory of plasma confinement in non-axisymmetric magnetic fields. Reports on Progress in Physics 77 (8), 087001.
  * Helander et al. (2011) Helander, P, Mishchenko, A, Kleiber, R & Xanthopoulos, P 2011 Oscillations of zonal flows in stellarators. Plasma Physics and Controlled Fusion 53 (5), 054006.
  * Helander & Nührenberg (2009) Helander, P. & Nührenberg, J. 2009 Bootstrap current and neoclassical transport in quasi-isodynamic stellarators. Plasma Physics and Controlled Fusion 51 (5), 055004.
  * Helander & Sigmar (2005) Helander, Per & Sigmar, Dieter J 2005 Collisional transport in magnetized plasmas, , vol. 4. Cambridge University Press.
  * Ho & Kulsrud (1987) Ho, Darwin D.-M. & Kulsrud, Russell M. 1987 Neoclassical transport in stellarators. The Physics of Fluids 30 (2), 442–461.
  * Jorge & Landreman (2020) Jorge, Rogerio & Landreman, Matt 2020 The use of near-axis magnetic fields for stellarator turbulence simulations. Plasma Physics and Controlled Fusion 63 (1), 014001.
  * Landreman (2022) Landreman, Matt 2022 Mapping the space of quasisymmetric stellarators using optimized near-axis expansion. Journal of Plasma Physics 88 (6), 905880616.
  * Landreman & Catto (2012) Landreman, M. & Catto, P. J. 2012 Omnigenity as generalized quasisymmetry. Physics of Plasmas 19 (5), 056103.
  * Landreman & Paul (2022) Landreman, M. & Paul, E. 2022 Magnetic fields with precise quasisymmetry for plasma confinement. Physical Review Letters 128 (3), 035001.
  * Landreman & Sengupta (2019) Landreman, M. & Sengupta, W. 2019 Constructing stellarators with quasisymmetry to high order. Journal of Plasma Physics 85 (6), 815850601.
  * Mikhailov et al. (2002) Mikhailov, M. I., Shafranov, V. D., Subbotin, A. A., Isaev, M. Y., Nührenberg, J., Zille, R. & Cooper, W. A. 2002 42 (11), L23–L26.
  * Miller et al. (1998) Miller, RL, Chu, MS, Greene, JM, Lin-Liu, YR & Waltz, RE 1998 Noncircular, finite aspect ratio, local equilibrium model. Physics of Plasmas 5 (4), 973–978.
  * Mishchenko et al. (2008) Mishchenko, Alexey, Helander, Per & Könies, Axel 2008 Collisionless dynamics of zonal flows in stellarator geometry. Physics of Plasmas 15 (7).
  * Mishchenko & Kleiber (2012) Mishchenko, Alexey & Kleiber, Ralf 2012 Zonal flows in stellarators in an ambient radial electric field. Physics of Plasmas 19 (7).
  * Monreal et al. (2016) Monreal, Pedro, Calvo, Iván, Sánchez, Edilberto, Parra, Félix I, Bustos, Andrés, Könies, Axel, Kleiber, Ralf & Görler, Tobias 2016 Residual zonal flows in tokamaks and stellarators at arbitrary wavelengths. Plasma Physics and Controlled Fusion 58 (4), 045018.
  * Monreal et al. (2017) Monreal, Pedro, Sánchez, Edilberto, Calvo, Iván, Bustos, Andrés, Parra, Félix I, Mishchenko, Alexey, Könies, Axel & Kleiber, Ralf 2017 Semianalytical calculation of the zonal-flow oscillation frequency in stellarators. Plasma Physics and Controlled Fusion 59 (6), 065005.
  * Mukhovatov & Shafranov (1971) Mukhovatov, VS & Shafranov, VD 1971 Plasma equilibrium in a tokamak. Nuclear Fusion 11 (6), 605.
  * Mynick (2006) Mynick, H. E. 2006 Transport optimization in stellarators. Physics of Plasmas 13 (5), 058102.
  * Nemov et al. (1999) Nemov, V. V., Kasilov, S. V., Kernbichler, W. & Heyn, M. F. 1999 Evaluation of $1/\nu$ neoclassical transport in stellarators. Physics of Plasmas 6 (12), 4622–4632.
  * Nührenberg & Zille (1988) Nührenberg, J. & Zille, R. 1988 Quasi-helically symmetric toroidal stellarators. Physics Letters A 129 (2), 113 – 117.
  * Nührenberg (2010) Nührenberg, Jürgen 2010 Development of quasi-isodynamic stellarators. Plasma Physics and Controlled Fusion 52 (12), 124003.
  * Olver et al. (2020) Olver, F. W. J., Daalhuis, A. B. Olde, Lozier, D. W., Schneider, B. I., Boisvert, R. F., Clark, C. W., B. R. Mille and, B. V. Saunders, Cohl, H. S. & M. A. McClain, eds. 2020 Nist digital library of mathematical functions. http://dlmf.nist.gov/, Release 1.0.26 of 2020-03-15.
  * Plunk & Helander (2024) Plunk, GG & Helander, P 2024 The residual flow in well-optimized stellarators. Journal of Plasma Physics 90 (2), 905900205.
  * Plunk (2024) Plunk, G. G., et al 2024 A geometric approach to constructing quasi-isodynamic fields. In preparation.
  * Plunk et al. (2019) Plunk, G. G., Landreman, M. & Helander, P. 2019 Direct construction of optimized stellarator shapes. part 3. omnigenity near the magnetic axis. Journal of Plasma Physics 85 (6), 905850602.
  * Rodriguez et al. (2020) Rodriguez, E., Helander, P. & Bhattacharjee, A. 2020 Necessary and sufficient conditions for quasisymmetry. Physics of Plasmas 27 (6), 062501.
  * Rodriguez et al. (2022) Rodriguez, E., Sengupta, W. & Bhattacharjee, A. 2022 Phases and phase-transitions in quasisymmetric configuration space. Plasma Physics and Controlled Fusion 64 (10), 105006.
  * Rodríguez et al. (2023) Rodríguez, E, Sengupta, W & Bhattacharjee, A 2023 Constructing the space of quasisymmetric stellarators through near-axis expansion. Plasma Physics and Controlled Fusion 65 (9), 095004.
  * Rodríguez (2023) Rodríguez, E. 2023 Magnetohydrodynamic stability and the effects of shaping: a near-axis view for tokamaks and quasisymmetric stellarators. Journal of Plasma Physics 89 (2), 905890211.
  * Rodríguez et al. (2020) Rodríguez, E., Helander, P. & Bhattacharjee, A. 2020 Necessary and sufficient conditions for quasisymmetry. Physics of Plasmas 27 (6), 062501.
  * Rodríguez & Plunk (2023) Rodríguez, E. & Plunk, G. G. 2023 Higher order theory of quasi-isodynamicity near the magnetic axis of stellarators. Physics of Plasmas 30 (6), 062507.
  * Rosenbluth & Hinton (1998) Rosenbluth, MN & Hinton, FL 1998 Poloidal flow driven by ion-temperature-gradient turbulence in tokamaks. Physical review letters 80 (4), 724.
  * Schiff (2013) Schiff, Joel L 2013 The Laplace transform: theory and applications. Springer Science & Business Media.
  * Skovoroda (2005) Skovoroda, A. A. 2005 3d toroidal geometry of currentless magnetic configurations with improved confinement. Plasma Physics and Controlled Fusion 47 (11), 1911–1924.
  * Spitzer Jr (1958) Spitzer Jr, Lyman 1958 The stellarator concept. The Physics of Fluids 1 (4), 253–264.
  * Stringer (1972) Stringer, TE 1972 Effect of the magnetic field ripple on diffusion in tokamaks. Nuclear Fusion 12 (6), 689.
  * Sugama & Watanabe (2005) Sugama, H & Watanabe, T-H 2005 Dynamics of zonal flows in helical systems. Physical review letters 94 (11), 115001.
  * Sugama & Watanabe (2006) Sugama, Hideo & Watanabe, T-H 2006 Collisionless damping of zonal flows in helical systems. Physics of Plasmas 13 (1).
  * Takahasi & Mori (1974) Takahasi, Hidetosi & Mori, Masatake 1974 Double exponential formulas for numerical integration. Publications of the Research Institute for Mathematical Sciences 9 (3), 721–741.
  * Watanabe et al. (2008) Watanabe, T-H, Sugama, H & Ferrando-Margalet, S 2008 Reduction of turbulent transport with zonal flows enhanced in helical systems. Physical review letters 100 (19), 195002.
  * Wesson (2011) Wesson, John 2011 Tokamaks; 4th ed.. International series of monographs on physics . Oxford: Oxford Univ. Press.
  * Weyl (1916) Weyl, Hermann 1916 Über die gleichverteilung von zahlen mod. eins. Mathematische Annalen 77 (3), 313–352.
  * Xanthopoulos et al. (2011) Xanthopoulos, P, Mischchenko, A, Helander, P, Sugama, H & Watanabe, T-H 2011 Zonal flow dynamics and control of turbulent transport in stellarators. Physical review letters 107 (24), 245002.
  * Xiao & Catto (2006) Xiao, Yong & Catto, Peter J 2006 Short wavelength effects on the collisionless neoclassical polarization and residual zonal flow level. Physics of Plasmas 13 (10).
  * Xiao et al. (2007) Xiao, Yong, Catto, Peter J & Dorland, William 2007 Effects of finite poloidal gyroradius, shaping, and collisions on the zonal flow residual. Physics of plasmas 14 (5).