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 (Reference Rosenbluth and Hinton1998), Xiao & Catto (Reference Xiao and Catto2006), Monreal et al. (Reference Monreal, Calvo, Sánchez, Parra, Bustos, Könies, Kleiber and Görler2016) and Plunk & Helander (Reference Plunk and Helander2024), 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. Reference Connor, Hastie and Taylor1978, Reference Connor, Hastie and Taylor1980)

(2.1) \begin{equation} \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 , \end{equation}

written in the ballooning formalism with the variation perpendicular to the field line described by the wave-vector $\boldsymbol{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 (2.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=\boldsymbol{v}_D\mathbf{\mathbf{\cdot}} \boldsymbol{k}_\perp =v_T \rho \bar {B} k_\psi {\boldsymbol \kappa }\times \boldsymbol{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 (2.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 Reference Schiff2013, Theorem 2.7) yields

(2.2) \begin{equation} \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), \end{equation}

where $\delta \!F(0)\stackrel {\mathbf{\cdot} }{=}g(0)-(q/T)J_0F_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$ ,

(2.3) \begin{equation} v_\parallel \frac {\partial }{\partial \ell }\delta =\tilde {\omega }_d - \overline {\tilde {\omega }_d} ,\end{equation}

so that we may write

(2.4) \begin{equation} \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), \end{equation}

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

(2.5) \begin{equation} \overline {f}=\begin{cases} \begin{aligned} &\frac {1}{\tau _b}\frac {1}{v}\int _{\textit{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} \end{equation}

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$ , (2.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 Reference Hall and McNamara1975a ; Cary & Shasharina Reference Cary and Shasharina1997; Helander & Nührenberg Reference Helander and Nührenberg2009; Landreman & Catto Reference Landreman and Catto2012), which we shall assume to hold throughout this work. For a treatment of the non-omnigenous problem see Helander et al. (Reference Helander, Mishchenko, Kleiber and Xanthopoulos2011) and Monreal et al. (Reference Monreal, Calvo, Sánchez, Parra, Bustos, Könies, Kleiber and Görler2016).

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)}+\cdots$ and $\hat {\phi }=\hat {\phi }^{(0)}+\hat {\phi }^{(1)}+\cdots$ , and considering (2.4) order by order

(2.6a) \begin{align} iv_\parallel \frac {\partial }{\partial \ell }\hat {h}^{(0)} &\approx 0, \end{align}

(2.6b) \begin{align} iv_\parallel \frac {\partial }{\partial \ell }\hat {h}^{(1)}+\omega \hat {h}^{(0)} &\approx \frac {q}{T}\omega F_0J_0\hat {\phi }^{(0)}e^{i\delta }+ie^{i\delta }\delta \!F(0), \end{align}

(2.6c) \begin{align} & \vdots \end{align}

From (2.6a ) it follows that

(2.7) \begin{equation} \hat {h}^{(0)}=\overline {\hat {h}^{(0)}}. \end{equation}

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

(2.8) \begin{equation} \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 }}. \end{equation}

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

(2.9) \begin{equation} \sum _{e,i}\int \mathrm{d}^3\boldsymbol{v}J_0 \hat {g}=n\frac {q_i}{T_i}(1+\tau )\hat {\phi }, \end{equation}

where $\tau =T_i/ZT_e$ and $Z = -q_i/q_e$ , then yields

(2.10) \begin{equation} \hat {\phi }^{(0)}\approx \frac {1}{n}\left \langle \int \mathrm{d}^3\boldsymbol{v}J_0e^{-i\delta }\overline {J_0e^{i\delta }}F_0 \right \rangle _\psi \hat {\phi }^{(0)}+\frac {i}{\omega }\frac {1}{n}\left \langle \int \mathrm{d}^3\boldsymbol{v}J_0e^{-i\delta }\overline {\delta \!F(0)e^{i\delta }} \right \rangle _\psi . \end{equation}

Here $\langle \cdots \rangle _\psi$ denotes a flux surface average (Helander Reference Helander2014), 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 Reference Schiff2013, Theorem 2.36), we obtain

(2.11) \begin{equation} \phi (\infty )=\frac {({1}/{n})\left \langle \int \mathrm{d}^3\boldsymbol{v}J_0e^{-i\delta }\overline {\delta \!F(0)e^{i\delta }} \right \rangle _\psi }{1- ({1}/{n})\left \langle \int \mathrm{d}^3\boldsymbol{v}J_0e^{-i\delta }\overline {J_0e^{i\delta }}F_0 \right \rangle _\psi }. \end{equation}

To finalise the calculation of the residual, we must consider some initial perturbation of the ion population. Following Rosenbluth & Hinton (Reference Rosenbluth and Hinton1998) and Monreal et al. (Reference Monreal, Calvo, Sánchez, Parra, Bustos, Könies, Kleiber and Görler2016), 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. Reference Monreal, Calvo, Sánchez, Parra, Bustos, Könies, Kleiber and Görler2016). 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

(2.12) \begin{equation} \frac {\phi (\infty )}{\phi (0)}\approx \frac {1-\Gamma _0}{\Gamma _0}\frac {({1}/{n})\left \langle \int \mathrm{d}^3\boldsymbol{v}J_0e^{-i\delta }\overline {J_0 e^{i\delta }} F_0 \right \rangle _\psi }{1-({1}/{n})\left \langle \int \mathrm{d}^3\boldsymbol{v}J_0e^{-i\delta }\overline {J_0 e^{i\delta }} F_0 \right \rangle _\psi }. \end{equation}

2.2. Finite orbit width

In order to proceed with the evaluation of (2.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. This treatment of $B$ and $\omega _d$ as independent quantities is particularly important to capture the behaviour of relevant stellarator fields such as quasi-isodynamic ones. In such fields, the variation of $|\boldsymbol{B}|$ along the field line is dominated by a strong toroidal variation (a variation that may be thought as a generalisation of the concept of linked magnetic mirrors), while its poloidal variation controls the drift and thus remains rather independent of the former. Thus, the small mirror limit (at fixed drift) as formulated in this paper is a sensible one. Of course, in other scenarios such as that of a tokamak, the connection is stronger, but we may always relate $B$ and $\omega _d$ at a later stage.

Despite this independence, it is important to respect some minimal properties. First, for the choice of functions to appropriately represent the behaviour in an omnigenous 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. Reference Mikhailov, Shafranov, Subbotin, Isaev, Nührenberg, Zille and Cooper2002; Skovoroda Reference Skovoroda2005), and is necessary to represent an omnigenous 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 (2.3), as

(2.13) \begin{equation} \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 }')}{\omega _t}\mathrm{d}\bar {\ell }' ,\end{equation}

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.Footnote ¹

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

(2.14) \begin{equation} \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 }'}\left (\frac {B^2 }{\partial _{\bar {\ell }'} B}\frac {\omega _d(\bar {\ell }')}{\omega _t}\right )\mathrm{d}{\bar {\ell }}'. \end{equation}

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 line,Footnote ² 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, $\varDelta$ 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 proportionality constant, simplifies (2.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.

The shape of the orbits may also undergo significant changes when the scale of variation of the background field is significant compared with the orbit $\delta$ . Noting that the orbit width in real space is $\delta /k_\perp$ , we can expect those non-local effects to become relevant when $L/\rho _i\sim \delta /(k_\perp \rho _i)$ , where $L$ is the scale of variation of the background. A breakdown is thus expected sufficiently close to the magnetic axis of the field, when $L\sim r$ , the distance from the axis. Physically, it is also true that studying a zonal perturbation with a radial scale that is longer that the distance to the magnetic axis is incompatible.

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 $\varDelta =0.05$ . The dotted line on top and bottom correspond to the $\delta (0)$ estimate in (2.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

(2.15) \begin{equation} \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 }}}+\sqrt {\hat {\lambda }} &\quad \text{if } \hat {\lambda }\gg \Delta . \end{aligned} \end{cases} \end{equation}

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)|\gt 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 (2.15), so that

(2.16) \begin{equation} \hat {\lambda }\lt \epsilon ^2=\left (\frac {\omega _d}{\pi \omega _t}\right )^2. \end{equation}

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

(2.17) \begin{equation} \Delta /\epsilon ^2\ll 1. \end{equation}

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 (Reference Rosenbluth and Hinton1998)). 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

(2.18) \begin{equation} \delta _{\mathrm{trap}}(0) \approx -\sigma \frac {v}{v_T}\epsilon \sqrt {\frac {2\bar {\kappa }}{\Delta }}, \end{equation}

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

(2.19) \begin{equation} \bar {\kappa }\lt \frac {\Delta }{\epsilon ^2}, \end{equation}

which we have already assumed small.

2.3. Evaluating the residual for small mirror ratio

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 (2.12). We refer to each of these groups by Roman numerals I to IV, starting from strongly passing particles (see figure 2).

Figure 2.

Separation of particles into groups. The diagram depicts the separation of the particle population into four different groups (I–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 })$ .

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 (2.12)

(2.20) \begin{equation} 1-\frac {1}{n}\left \langle \int \mathrm{d}^3\boldsymbol{v}J_0e^{-i\delta }\overline {J_0 e^{i\delta }} F_0 \right \rangle _\psi \approx b-\frac {1}{n}\left \langle \int \mathrm{d}^3\boldsymbol{v}\left (e^{-i\delta }\overline {e^{i\delta }}-1\right )F_0 \right \rangle _\psi , \end{equation}

where we used

(2.21) \begin{equation} \frac {1}{n}\left \langle \int \mathrm{d}^3\boldsymbol{v}J_0^2 F_0 \right \rangle _\psi =\frac {1}{2}e^{-b}I_0(b), \end{equation}

in the small $b$ limit and the velocity space integrals include all groups. The integral remaining in (2.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 ^2b)$ 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 (2.20) into the different group contributions

(2.22) \begin{align} I=\frac {1}{n}\left \langle \int \mathrm{d}^3\boldsymbol{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\boldsymbol{v}\left (\overline {\delta }^2-\overline {\delta ^2}\right )F_0 \right \rangle _\psi \nonumber \\ +\sum _{\mathrm{II,\,III}}\frac {1}{n}\left \langle \int \mathrm{d}^3\boldsymbol{v}\left (e^{-i\delta }\overline {e^{i\delta }}-1\right )F_0 \right \rangle _\psi . \end{align}

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 (2.22). This contribution should be familiar, as it has the quadratic form in which the Rosenbluth–Hinton residual is customarily written (Rosenbluth & Hinton Reference Rosenbluth and Hinton1998; Xiao & Catto Reference Xiao and Catto2006; Plunk & Helander Reference Plunk and Helander2024). 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 figure 2), and their contribution to (2.22)

(2.23) \begin{equation} I_{\mathrm{II}}=\underbrace {\frac {1}{n}\left \langle \int _{\mathrm{II}} \mathrm{d}^3\boldsymbol{v}e^{-i\delta }\overline {e^{i\delta }} F_0 \right \rangle _\psi }_{\unicode {x2460}} - \underbrace {\frac {1}{n}\left \langle \int _{\mathrm{II}} \mathrm{d}^3\boldsymbol{v}F_0 \right \rangle _\psi }_{\unicode {x2461}}. \end{equation}

First consider $\unicode {x2460}$ , and rewrite it following Xiao & Catto (Reference Xiao and Catto2006) as

(2.24) \begin{equation} \unicode {x2460} = \frac {1}{n}\left \langle \int _{\mathrm{II}} \mathrm{d}^3\boldsymbol{v}\left (\overline {\cos \delta }^2+\overline {\sin \delta }^2\right ) F_0 \right \rangle _\psi , \end{equation}

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 non-zero according to our convention for the bounce average in (2.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 Reference Bender and Orszag2013, § 6.5)

(2.25) \begin{equation} \overline {\cos \delta }=\frac {1}{\tau _t\omega _t}\frac {v_T}{v}\textit{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 ''(\ell _i)|}}\frac {\cos [\delta (\ell _i)-\pi /4]}{\sqrt {1-\lambda B(\ell _i)}}, \end{equation}

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

(2.26) \begin{equation} \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 ]. \end{equation}

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, (2.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 Reference Hazeltine and Meiss2003, § 4.4)

(2.27) \begin{equation} \mathrm{d}^3\boldsymbol{v}\rightarrow \frac {2\pi B}{\sqrt {1-\lambda B}}v^2\mathrm{d}v\mathrm{d}\lambda , \end{equation}

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

(2.28) \begin{equation} \left \langle \int _{\mathrm{II}}\mathrm{d}^3\boldsymbol{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 , \end{equation}

correct to leading order in $\Delta$ .

The simplifying assumption of a $v$ -independent boundary layer in (2.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$ )

(2.29a) \begin{align} &\int _0^\infty ve^{-v^2}\cos ^2\left (Av-\frac {\pi }{4}\right )\mathrm{d}v\approx \frac {1}{4}, \end{align}

(2.29b) \begin{align} &\int _0^\infty ve^{-v^2}\mathrm{d} v=\frac {1}{2}, \end{align}

we find using the explicit form of the Maxwellian $F_0$

(2.30) \begin{equation} \unicode {x2460}\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 }, \end{equation}

where $\hat {\tau }_t=\tau _t(v/v_T)$ is a function of $\lambda$ . In this form of $\unicode {x2460}$ 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 § 2.2.1. The integration domain becomes $\kappa \in [2\Delta /\epsilon ^2,1]$ , with an integral measure

(2.31) \begin{equation} \frac {\mathrm{d}\kappa }{\mathrm{d}\lambda }=2\Delta \bar {B}\left (1+\frac {1-\Delta }{2\Delta }\kappa \right )^2. \end{equation}

The contribution from the edges of the orbit (the first term in (2.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$

(2.32) \begin{equation} \unicode {x2460}\approx \frac {\epsilon }{\pi ^{3/2}}. \end{equation}

This concludes the calculation of $\unicode {x2460}$ , but $\unicode {x2461}$ remains to be found. This contribution corresponds to finding the fraction of phase space occupied by the barely passing particles in group II. Using (2.28) and the definition of region II, the integrals over $\kappa$ and $v$ yield

(2.33) \begin{equation} \unicode {x2461}\approx \epsilon . \end{equation}

Altogether

(2.34) \begin{equation} I_{\mathrm{II}}\approx \frac {\epsilon }{\pi ^{3/2}}(1-\pi ^{3/2})\approx -0.26\frac {\omega _d}{\omega _t}, \end{equation}

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, $\unicode {x2460}$ and $\unicode {x2461}$ , as in (2.23). In the calculation of $\unicode {x2460}$ , 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, (2.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 § 2.2.2 and Appendix A, and performing the integral over $v$ first

(2.35) \begin{equation} \unicode {x2460}\approx \frac {1}{\pi ^{3/2}}\frac {\Delta }{\epsilon }, \end{equation}

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, $\unicode {x2461}$ , is of course also small in the limit of a small mirror ratio, $\unicode {x2461}\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 (2.22) evaluates to

(2.36) \begin{equation} I\approx -0.26\frac {\omega _d}{\omega _t}, \end{equation}

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 (2.36).

We may now write the expression for the residual itself, going back to (2.12) using the definition of $I$ in (2.22)

(2.37) \begin{equation} \frac {\phi (\infty )}{\phi (0)}\approx \frac {1}{1+0.26({\omega _d}/{b\omega _t})}, \end{equation}

which in the limit of $b\ll \omega _d$ , say for very long radial wavelengths, can be expressed as

(2.38) \begin{equation} \frac {\phi (\infty )}{\phi (0)}\approx 1.92\,k_\perp \rho _i\left (\frac {k_\perp \rho _i}{\omega _d/\omega _t}\right ). \end{equation}

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, however, 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 (2.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 $|\boldsymbol{B}|$ modified, leaving the curvature drift unchanged. The colour 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 error bars 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 (2.12), the dotted black line the analytical estimate of Xiao & Catto (Reference Xiao and Catto2006) and the red dotted line the asymptotic expression in (2.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 (2.12) and linear electrostatic gyrokinetic simulations with the stella code (Barnes et al. Reference Barnes, Parra and Landreman2019). 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. Reference Miller, Chu, Greene, Lin-Liu and Waltz1998)) whose $B$ has been modified with varying mirror ratios $\Delta$ , while keeping all other elements of the geometry unchanged. The numerical evaluation of (2.12) is done by careful treatment of bounce integrals using double-exponential integration methods (Takahasi & Mori Reference Takahasi and Mori1974) 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 with 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, of 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.Footnote ³ 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 (2.12), we can assess the contribution from different regions of velocity space to the residual using the latter (see figure 3 c).

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, (2.38), agrees to a good degree (within ${\sim}5\,\%{-}10\,\%$ ) with the simulation and integration (see red line in figure 3 b). As the mirror ratio increases the importance within velocity space shifts (see figure 3 c) and the bulk of trapped particles becomes dominant (the standard Rosenbluth & Hinton (Reference Rosenbluth and Hinton1998) picture). In that limit the residual can be estimated by Rosenbluth & Hinton (Reference Rosenbluth and Hinton1998) (RH)

(3.1) \begin{equation} \left .\frac {\phi (\infty )}{\phi (0)}\right |_{\mathrm{RH}}=\frac {1}{1+1.6 \epsilon ^2/(k_\perp \rho _i)^2\sqrt {\Delta }}, \end{equation}

or more precisely by Xiao & Catto (Reference Xiao and Catto2006), as explicitly shown in figure 3(b) (black dotted line). The standard RH residual, (3.1), exhibits a stronger dependence on the drift and transit time compared with 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, (2.38), and the RH residual (see figure 3 b). This takes place when $\Delta \sim \epsilon ^2$ , which is approximately

(3.2) \begin{equation} \Delta _t\approx 0.1(k_\perp \rho _i)^2\left (\frac {\omega _d/\omega _t}{k_\perp \rho _i}\right )^2. \end{equation}

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\lt 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 Xiao & Catto (Reference Xiao and Catto2006) residual (XC) 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 3(a). There are clearly two oscillation time scales in the problem set-up considered: the faster damped geodesic acoustic modes (GAMs) (Sugama & Watanabe Reference Sugama and Watanabe2006; Gao et al. Reference Gao, Itoh, Sanuki and Dong2006, Reference Gao, Itoh, Sanuki and Dong2008; Conway et al. Reference Conway, Smolyakov and Ido2021) 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-omnigenous nature of stellarator fields (Mishchenko, Helander & Könies Reference Mishchenko, Helander and Könies2008; Helander et al. Reference Helander, Mishchenko, Kleiber and Xanthopoulos2011; Mishchenko & Kleiber Reference Mishchenko and Kleiber2012; Alonso et al. Reference Alonso2017; Monreal 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 with their radial drift, these particles behave de facto as non-omnigenous 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-omnigenous 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 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.

The 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- $|\boldsymbol{B}|$ field, where the dynamics is 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 Reference Sugama and Watanabe2005, Reference Sugama and Watanabe2006; Gao et al. Reference Gao, Itoh, Sanuki and Dong2006, Reference Gao, Itoh, Sanuki and Dong2008). 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 Reference Schiff2013, 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. (Reference Gao, Itoh, Sanuki and Dong2006) with the residual, due to its similar form. This purely damped mode reads

(3.3) \begin{equation} \frac {\phi (t\rightarrow \infty )}{\phi (0)}\approx \frac {1}{1+({\epsilon ^2}/{4b})\left [1+({\pi }/({2(1+\tau )}))\right ]}e^{-\gamma t}, \end{equation}

where

(3.4) \begin{equation} \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}. \end{equation}

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 (2.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 becomes dominating, which occurs at small mirror ratios ( $\Delta _t\lt 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. This analysis is particularly important in order to understand when the formal consideration of a small mirror ratio at a fixed magnetic drift is physically relevant.

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 (2.37) and (3.1)

(4.1) \begin{equation} \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.6 q^2/\sqrt {\Delta }}. \end{equation}

The larger the $q$ , the larger the connection length, the larger the orbit width $\delta$ and the lower the residual. In terms of these tokamak parameters, we may also rewrite the condition for the regime transition in (3.2): the layer contribution becomes relevant for $\Delta \lt \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 per cent. 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 \gt 1$ for all particles, corresponding to a completely different regime.Footnote ⁴

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 (QS) (Boozer Reference Boozer1983a ; Nührenberg & Zille Reference Nührenberg and Zille1988; Rodriguez, Helander & Bhattacharjee Reference Rodriguez, Helander and Bhattacharjee2020) and quasi-isodynamic (QI) (Cary & Shasharina Reference Cary and Shasharina1997; Helander & Nührenberg Reference Helander and Nührenberg2009; Nührenberg Reference Nührenberg2010) ones. The former can be seen as the natural generalisation of the axisymmetric case, where the field has a direction of symmetry on $|\boldsymbol{B}|$ instead of the whole vector $\boldsymbol{B}$ . The direction of symmetry can be toroidal quasi-axisymmetry (QA) or helical quasi-helical (QH). This symmetry forces the magnetic wells along the field line to be all nearly identical (same $B$ and $\omega _d$ (Boozer Reference Boozer1983b ), but different $k_\perp \rho _i$ ). In quasi-isodynamic fields, the contours of $|\boldsymbol{B}|$ are closed poloidally, and carefully shaped to grant omnigeneity (Hall & McNamara Reference Hall and McNamara1975; Bernardin, Moses & Tataronis Reference Bernardin, Moses and Tataronis1986; Cary & Shasharina Reference Cary and Shasharina1997; Helander Reference Helander2014). As a result, wells are differently shaped, but all share the feature of being omnigenous; 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. (4.1). These parameters will be derived using the inverse-coordinate near-axis description of equilibria (Garren & Boozer Reference Garren and Boozer1991b ; Landreman & Sengupta Reference Landreman and Sengupta2019; Plunk, Landreman & Helander Reference Plunk, Landreman and Helander2019; Rodríguez, Sengupta & Bhattacharjee Reference Rodríguez, Sengupta and Bhattacharjee2023), as detailed in Appendix C, and summarised in table 1. We have included the case of a shaped tokamak for comparison. The near-axis description of optimised stellarators has proven a powerful tool in understanding the behaviour of these fields. Although the residual calculation in the paper is bound to fail close to the axis (in the gyroradius sense), note that there remains a scale separation between this and the region of validity of the near-axis description, that is concerned with equilibrium scales, and thus insight may be still gained. Let us now discuss the interpretation of these results.

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 $|\boldsymbol{B}|$ over flux surfaces (roughly proportional to the axis curvature) and $\hat {\mathcal{G}}$ geometric factor defined in (4.2).

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 transition value $\Delta _t$ (bottom) for DIII-D (equilibrium from Austin et al. (Reference Austin2019), shot 170 680 at 2200 ms) (tokamak), precise QA (QA stellarator) and precise QH (QH stellarator) configurations (Landreman & Paul Reference Landreman and Paul2022). The black broken line indicates the region close to the axis where non-local effects on the residual start becoming relevant, namely $r\sim \rho _i/(k_\perp \rho _i)$ (we took the DIIID shot as a reference with $B\sim 2\mathrm{T}$ and $T_i\sim 4\,\mathrm{keV}$ for all cases). The residual is computed numerically evaluating (2.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.

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 $|\boldsymbol{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, however, this region tends to be narrow, and thus likely unimportant (see figure 4). We note that the proximity to the axis limits the local description of the residual (see e.g. the broken black line in figure 4), although this will generally depend on the particular realisation of the equilibrium field, most importantly, its physical size and strength of magnetic field, which can aid in the scale separation. 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 Reference Landreman and Paul2022; Landreman Reference Landreman2022; Rodríguez, Sengupta, & Bhattacharjee Reference Rodríguez, Sengupta and Bhattacharjee2023; Giuliani Reference Giuliani2024). 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 $|\boldsymbol{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, negating the presence of a region where $\Delta \rightarrow 0$ . In return, this provides additional control on $\Delta$ , which can be imposed rather independently from $\omega _d$ .

In addition to the differences in $\Delta$ , the changes in the magnitude of the magnetic field gradient $\boldsymbol \nabla 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 (4.1) 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 (C.9) and (C.22))

(4.2) \begin{equation} \hat {\mathcal{G}}^2=\frac {1}{2\pi }\int _0^{2\pi } \frac {\mathrm{d}\varphi }{\sin 2e}, \end{equation}

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 Reference Rodríguez2023) 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}}\gt 1$ (as $\sin 2e\lt 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. (Reference Xiao, Catto and Dorland2007). 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}}\gt 1$ (Landreman & Sengupta Reference Landreman and Sengupta2019; Camacho Mata, Plunk & Jorge Reference Camacho Mata, Plunk and Jorge2022; Rodríguez Reference Rodríguez2023), 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 $|\boldsymbol{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, Sengupta & Bhattacharjee Reference Rodríguez, Sengupta and Bhattacharjee2023), where $\kappa$ is the axis curvature. For similarly shaped cross-sections, $\eta$ (or $\bar {d}$ ) will be larger for QH and QI stellarators compared with QA and tokamaks (Rodriguez, Sengupta & Bhattacharjee Reference Rodriguez, Sengupta and Bhattacharjee2022; Camacho Mata et al. Reference Camacho Mata, Plunk and Jorge2022), 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 Reference Helander and Nührenberg2009; Landreman & Catto Reference Landreman and Catto2012). 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 $|\boldsymbol{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

(4.3) \begin{equation} q_{\mathrm{eff, tok}}\sim q_{\mathrm{eff, QA}}\gt q_{\mathrm{eff, QH}}\gtrsim q_{\mathrm{eff, QI}}. \end{equation}

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 (Reference Plunk and Helander2024). 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 (Reference Landreman2022) 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.

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 (Reference Landreman2022), which serves as a representative population of optimised QS configurations. The densities for each number of field period (colour) are stacked vertically on top of one another, and represent the number of configurations in the database satisfying those parameters. The rightmost plot shows the same analysis through a QI near-axis database (Plunk et al. Reference Plunk2024). This shows the rough relative ordering of $q_{\mathrm{eff}}$ between different omnigenous 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.