2 Gyrokinetic solution of the initial value problem

We will follow some of the notation conventions of Helander et al. (Reference Helander, Mishchenko, Kleiber and Xanthopoulos2011a), but with some adaptations due to the fact that the derivation will use gyrokinetics instead of drift kinetics. We are interested in solving the gyrokinetic system of equations in the electrostatic limit for the time evolution of a linear mode with wavenumber ${\boldsymbol {k}}_\perp = k_r \boldsymbol {\nabla } r + k_\alpha \boldsymbol {\nabla } \alpha$, where $r(\psi )$ is an arbitrary radial coordinate that is constant on magnetic surfaces, and the magnetic field is expressed as ${\boldsymbol {B}} = \boldsymbol {\nabla }\psi \times \boldsymbol {\nabla }\alpha$ in terms of the toroidal flux function $\psi$ and Clebsch angle $\alpha$. Here, we are concerned only in the case $k_\alpha = 0$, where there is no source of free energy for the perturbation, so the linear mode is stable. The collisionless gyrokinetic equation in this limit is

(2.1)\begin{equation} \frac{\partial g_{a}}{\partial t} + v_{\|} \boldsymbol{\nabla}_{\|} g_a + {\rm i} \omega_{da} g_{a} = \frac{e_a F_{a0}}{T_a} \frac{\partial\phi}{\partial t} {\rm J}_{0a} , \end{equation}

with $\boldsymbol {\nabla }_{\|} = \hat {\boldsymbol {b}}\boldsymbol {\cdot }\boldsymbol {\nabla } = \partial /\partial l$, where $\hat {\boldsymbol {b}} = {\boldsymbol {B}}/B$ and $l$ is the arc length along field line. Here, $g_a$ is the gyrocentre distribution function for species ‘$a$’ and $F_{a0} = n_a (m_a/2{\rm \pi} T_a)^{3/2} \exp (-m_a v^2/2 T_a)$, with $n_a$ and $T_a$ the bulk density and temperature, $m_a$ the particle mass, and $\textbf{v}$ its velocity. The electrostatic potential $\phi$ is found from the quasi-neutrality constraint

(2.2)\begin{equation} \sum_a n_a\frac{e_a^2}{T_a} \phi = \sum_a e_a \int g_{a} {\rm J}_{0a} \,{\rm d}^3v, \end{equation}

where ${\rm J}_{na} = {\rm J}_n(k_\perp v_\perp / \varOmega _a)$, with $\varOmega _a = e_a B/m_a$, with $e_a$ the particle charge, and the velocity element is expressed in gyrokinetic phase space variables as

(2.3)\begin{equation} {\rm d}^3v = 2 {\rm \pi}v_\perp \,{\rm d}v_\perp v_{\|} = \sum_\sigma \frac{2 {\rm \pi}B \,{\rm d}E_a \,{\rm d}\mu_a}{m_a^2 | v_{\|} |} = \sum_\sigma \frac{{\rm \pi} B v^2 \,{\rm d}v \,{\rm d}{\lambda}}{\sqrt{1-{\lambda} B}}, \end{equation}

where we define $v_\perp = |\hat {\boldsymbol {b}}\times {\boldsymbol {v}}|$, $v_{\|} = \hat {\boldsymbol {b}}\boldsymbol {\cdot } {\boldsymbol {v}}$, $E_a = m_a v^2/2$, $\mu _a = m_a v_\perp ^2/(2 B)$, $v = |{\boldsymbol {v}}|$ and $\lambda = \mu _a/E_a$. In what follows, we mostly use phase space variables $v$ and $\lambda$, and express the parallel velocity as $v_{\|} = \sigma v \sqrt {1 - \lambda B}$ where $\sigma = v_{\|} / | v_{\|} | = \pm 1$. Following Helander et al. (Reference Helander, Mishchenko, Kleiber and Xanthopoulos2011a), the drift frequency is defined as $\omega _{da} = k_r {\boldsymbol {v}}_{da}\boldsymbol {\cdot }\boldsymbol {\nabla } r$ with

(2.4)\begin{equation} {\boldsymbol{v}}_{da}\boldsymbol{\cdot}\boldsymbol{\nabla} r = \bar{v}_{ra} + v_{{\parallel}} \boldsymbol{\nabla}_{\|}\delta_{ra}, \end{equation}

where $\bar {v}_{ra}$ denotes the transit-averaged radial drift, which is zero for passing particles and, in the case of tokamaks and omnigenous stellarators, also for trapped particles. The term $v_{\parallel } \boldsymbol {\nabla }_{\|}\delta _{ra}$ is zero under the orbit average, so $\delta _{ra}$ is the radial excursion, or ‘orbit width’, of particles, a periodic function on the torus. Equation (2.4) simply represents the splitting of the radial drift into mean and oscillatory behaviour with respect to the transit/bounce average. The quantities $\bar {v}_{ra}$ and $\delta _{ra}$ are defined by this equation, and may be obtained as solutions of it, with the appropriate boundary condition for $\delta _{ra}$ in $l$; for this, we take $\delta _{ra} = 0$ at bounce points of trapped-particle orbits, or at the field maximum for passing particles, implying that $\delta _{ra}$ is odd in $v_{\|}$. The transit average, which is designed to annihilate the operator $v_{\|}\boldsymbol {\nabla }_{\|}$, is defined as

(2.5)\begin{equation} \bar{f} = \left.\frac{1}{2}\sum_\sigma\int^{l_2}_{l_1} \frac{f}{\sqrt{1 - \lambda B(l)}}\,{\rm d}l\right/ \int^{l_2}_{l_1} \frac{1}{\sqrt{1 - \lambda B(l)}}\,{\rm d}l, \end{equation}

where, for trapped particles, $l_1$ and $l_2$ are the bounce points (where $v_{\|} = 0$) such that $B(l_1)=B(l_2) = 1/\lambda$, while for passing particles ($\lambda < 1/B_\mathrm {max}$) the average is understood in the limiting sense, as $l_1 \rightarrow -\infty$ and $l_2 \rightarrow \infty$.Footnote ¹ We also define the flux-surface average (Helander Reference Helander2014)

(2.6)\begin{equation} \left\langle \cdots \right\rangle = \left.\lim_{L\rightarrow \infty} \int_{-L}^L (\cdots ) \frac{{\rm d}l}{B} \right/ \int_{-L}^L \frac{{\rm d}l}{B}. \end{equation}

Following previous works, we take the Laplace transform of (2.1), defining $\hat {g}_a = \int _0^{\infty } \,{\rm d}t \exp (-pt) g_a(t)$, and introduce an integrating factor to absorb the orbit width term, defining $h_a = \exp ({\rm i} k_r \delta _{ra}) g_a$

(2.7)\begin{equation} (\,p + {\rm i} k_r \bar{v}_{ra})\hat{h}_a + v_{\|} \boldsymbol{\nabla}_{\|} \hat{h}_a = \left[\,p \frac{e_a\hat{\phi}}{T_a} {\rm J}_{0a} F_{a0} + \delta F_a(0)\right] {\rm e}^{{\rm i} k_r \delta_{ra}}, \end{equation}

where $\delta F_a(0)$ is the initial value of the gyrocentre distribution function

(2.8)\begin{equation} \delta F_{a} = g_{a} - \frac{e_a \phi}{T_a}{\rm J}_{0a} F_{a0}. \end{equation}

We are interested in times much longer that the transit/bounce time scale $\omega _b \hat {h} \sim v_{\|} \boldsymbol {\nabla }_{\|} \hat {h}$, which we order similar to the non-secular part of the drift frequency, i.e. $\omega _b \sim k_r v_{\|}\boldsymbol {\nabla }_{\|}\delta _r$, but $p \sim k_r \bar {v}_r \ll \omega _b$. This implies that, at dominant order, we have $\partial \hat {h}/\partial l = 0$. At next order, we transit average and use continuity of $g_a$ at bounce points, which, because $\delta _{ra} = 0$ at such points, implies $h_a|_{\sigma =1}=h_a|_{\sigma =-1}$ at all $l$, yielding

(2.9)\begin{equation} (\,p + {\rm i} k_r \bar{v}_{ra}) \hat{h}_a = \left(\,p\overline{\frac{e_a\hat{\phi}}{T_a} {\rm J}_{0a} {\rm e}^{{\rm i} k_r \delta_{ra}}}F_{a0} + \overline{\delta F_a(0) {\rm e}^{{\rm i} k_r \delta_{ra}}}\right)\!. \end{equation}

The solution for $\hat {g}_a$ is therefore

(2.10)\begin{equation} \hat{g}_a = \frac{1}{p + {\rm i} k_r \bar{v}_{ra}} \left(\,p\overline{\frac{e_a\hat{\phi}}{T_a} {\rm J}_{0a} {\rm e}^{{\rm i} k_r \delta_{ra}}}F_{a0} + \overline{\delta F_a(0) {\rm e}^{{\rm i} k_r \delta_{ra}}}\right) {\rm e}^{-{\rm i} k_r \delta_{ra}}.\end{equation}

This is the same as the result of Helander et al. (Reference Helander, Mishchenko, Kleiber and Xanthopoulos2011a) except that it retains the gyro-average (${\rm J}_{0a}$) and keeps $\phi$ and $\delta F_a$ under the orbit average, allowing them to vary within the flux surface. A similar comparison can be made with the result of Monreal et al. (Reference Monreal, Calvo, Sánchez, Parra, Bustos, Könies, Kleiber and Görler2016), which also retains the full gyro-average but not the mentioned flux-surface dependence. To obtain an equation simply for $\phi$, we substitute this expression into the quasi-neutrality condition (2.2), which gives

(2.11)\begin{gather} \sum_a \frac{e_a^2}{T_a}\left(n_a \hat{\phi} - \int \,{\rm d}^3v {\rm J}_{0a}F_{a0} \frac{p}{p + {\rm i} k_r \bar{v}_{ra}}\overline{\hat{\phi} {\rm J}_{0a} {\rm e}^{{\rm i} k_r \delta_{ra}}} {\rm e}^{-{\rm i} k_r \delta_{ra}}\right)\nonumber\\ = \sum_a e_a \int \,{\rm d}^3v {\rm J}_{0a}\frac{1}{p + {\rm i} k_r \bar{v}_{ra}}\overline{\delta F_a(0) {\rm e}^{{\rm i} k_r \delta_{ra}}} {\rm e}^{-{\rm i} k_r \delta_{ra}}. \end{gather}

Solving this equation, which can be compared with equation (6) of Helander et al. (Reference Helander, Mishchenko, Kleiber and Xanthopoulos2011a) (see also Appendix A), would give the fully general solution for the post-GAM zonal-flow dynamics, allowing for arbitrary orbit widths ($k \delta _r \sim k\rho \sim 1$). One would like to be able to solve (2.11) for $\hat {\phi }$ then invert the Laplace transform and obtain $\phi (t)$. Unfortunately, the situation is rather complicated, as $\hat {\phi }$ appears under an orbit average involving a resonant velocity integral, so we will have to take some limits to make further progress.

2.1 Limit of small orbit width and ion Larmor radius

Following previous works, we now consider the limit

(2.12)\begin{equation} k_r\delta_{ra} \sim k_\perp\rho_a \ll 1. \end{equation}

We will later take $k_\perp \rho _e = k_r \delta _{re} = 0$, since $\rho _e \ll \rho _i$ and $\delta _{re} \ll \delta _{ri}$. Some care is needed here as the polarization effects arising in the gyrokinetic equation enter at the same order as those which appear in the quasi-neutrality constraint, which itself is singular. This is made clear by recasting the quasi-neutrality condition (2.2) in terms of the gyrocentre distribution function,

(2.13)\begin{equation} \sum_a n_a \frac{e_a^2}{T_a} \left[1 - \varGamma_0(b_a) \right] \phi = \sum_a e_a \int \delta F_{a} {\rm J}_{0a} \,{\rm d}^3v = \sum_a e_a \delta n_a, \end{equation}

where $\varGamma _0(x) = {\rm I}_0(x) {\rm e}^{-x}$ and $b_a = k_\perp ^2 \rho _a^2 = k_\perp ^2 T_a / (m_a \varOmega _a^2)$. Note that, in the final expression, we introduce the ‘gyrocentre density’ $\delta n_a$. The point is that the gyrocentre density is small in the limit $b_i \ll 1$, i.e. taking $\varGamma _0(b) \approx 1 - b$, one sees that $\delta n/n \sim {O}(b_i e_i\phi /T_i)$

(2.14)\begin{equation} \sum_a e_a \delta n_a = b_i n_i \frac{e_i^2 \phi}{T_i}, \end{equation}

where $b_i = k_\perp ^2 \rho _i^2 = k_\perp ^2 m_i T_i/(e_i^2B^2)$, which means that it is necessary to include terms of order $b_i$ to solve for the electrostatic potential in this limit.

We now perform the expansion on (2.11), ordering $b_i e_i \phi /T_i \sim \delta n_a/n_a$. To put the final expression in a more convenient form, we rewrite the resonant term on the left-hand side using $p/(\,p + {\rm i} k_r \bar {v}_{ra}) = 1 - {\rm i} k_r \bar {v}_{ra}/(\,p + {\rm i} k_r \bar {v}_{ra})$, and thus obtain

(2.15)\begin{gather} \sum_a n_a\frac{e_a^2}{T_a}\left(\hat{\phi} - \frac{1}{n_a} \int \,{\rm d}^3v F_{a0} \bar{\hat{\phi}}\right) + \sum_a \frac{e_a^2}{T_a}\int \,{\rm d}^3 v \frac{{\rm i} k_r \bar{v}_{ra}}{p + {\rm i} k_r \bar{v}_{ra}}F_{a0}\bar{\hat{\phi}} \nonumber\\ + \sum_a \frac{e_a^2}{T_a}\int \,{\rm d}^3 v F_{a0}\left[\bar{\hat{\phi}} b_a x_\perp^2/2 + \overline{\hat{\phi} b_a x_\perp^2/2} - k_r^2 \overline{\hat{\phi}\delta_{ra}}\delta_{ra} + k_r^2 \overline{\hat{\phi}\delta_{ra}^2}/2+ k_r^2 \bar{\hat{\phi}}\delta_{ra}^2/2\right]\nonumber\\ = \frac{1}{p}\sum_a e_a \int \,{\rm d}^3v {\rm J}_{0a}\left(1 - \frac{{\rm i} k_r \bar{v}_{ra}}{p + {\rm i} k_r \bar{v}_{ra}}\right)\overline{\delta F_a(0) {\rm e}^{{\rm i} k_r \delta_{ra}}} {\rm e}^{-{\rm i} k_r \delta_{ra}}, \end{gather}

where we have written $x_\perp ^2 = m_a v_\perp ^2 / (2 T_a)$, expanded the Bessel function as ${\rm J}_{0a} (k_\perp v_\perp / \varOmega _a) \simeq 1 - b_a x_\perp ^2 / 2$ and recognized that terms that are linear in $\delta _{ra}$ are odd in $v_{\|}$ and vanish upon integration. This equation is still fairly complicated, but the terms can be identified. On the left-hand side, the first term shows contributions from the non-zonal part of the potential: it is zero under zonal average and also for the case $\phi = \langle \phi \rangle$; we return to this later. The second term contains the resonance that yields zonal-flow oscillations and damping in stellarators with $\bar {v}_{ra} \neq 0$ (Mishchenko et al. Reference Mishchenko, Helander and Könies2008; Helander et al. Reference Helander, Mishchenko, Kleiber and Xanthopoulos2011a; Monreal et al. Reference Monreal, Sánchez, Calvo, Bustos, Parra, Mishchenko, Könies and Kleiber2017). This is not the focus of the present work, although we discuss it briefly in Appendix A. Note that the finite orbit effects on this term are neglected, which is justified in the limit of small $k_r\bar {v}_{ra}/p$. On the second line we encounter all the finite orbit width (FOW) and finite Larmor radius (FLR) terms associated with the residual. These expressions will simplify significantly (and become more familiar) in limits when the potential is mostly zonal. Finally, on the right-hand side we have the contribution from the initial condition; note the separation into a resonant term and another which can be evaluated using quasi-neutrality in terms of the initial potential. We keep this term exact, for now, since we would like to discuss the consequences of several possible orderings for the initial condition itself, $\delta F_a(0)$ in the following section; see also Appendix C.

3 The residual potential

Let us consider the limit where the resonance can be neglected, i.e. let us take $p \gg k_r \bar {v}_{ra}$. Here, we neglect damping special to stellarators, which includes both exponential and algebraic decay (Helander et al. Reference Helander, Mishchenko, Kleiber and Xanthopoulos2011a). Although that damping process is slow, it can still manage to deplete most of the zonal amplitude. Indeed, as shown by equations (16)–(17) of Helander et al. (Reference Helander, Mishchenko, Kleiber and Xanthopoulos2011a), the final residual is independent of the size of the damping rate when that rate is small, i.e. any non-zero value of $\bar {v}_{ra}$ results in a strong correction to the residual that depends on the fraction of trapped particles. Therefore, the ‘residual’ that arises before this decay takes effect may be more relevant for understanding the interaction between zonal flows and turbulence, especially in optimized stellarators where $\bar {v}_{ra}$ is made to be as close to zero as possible. We therefore focus on the ‘intermediate residual’ defined to be the value of the potential long after the GAMs have decayed, $\tau _{G} \sim 1/\gamma _{\mathrm {GAM}}$, but long before the final residual is obtained, on the time scale of the stellarator-specific damping of Mishchenko and Helander, $\tau _{M} \sim 1/(k_r \bar {v}_{ra})$

(3.1)\begin{equation} \phi_\mathrm{res} \equiv \lim_{\frac{\tau_{M}}{t} \rightarrow \infty} \left(\lim_{\frac{t}{\tau_\mathrm{G}} \rightarrow \infty} \phi(t) \right)\!. \end{equation}

For tokamaks (and perfectly omnigenous stellarators) this quantity coincides with the conventional definition of the residual, as defined by Rosenbluth & Hinton (Reference Rosenbluth and Hinton1998), as shown in what follows.

To be slightly more formal, we note that the inner limit of (3.1) has already been taken much earlier in our calculation to obtain (2.9), and additional limits are to be considered subsidiary to that one. In particular, $(t \omega _b)^{-1} \ll k_r \bar {v}_{r} t \ll 1$, with the latter condition expressing the outer limit of (3.1). According to these orderings, $t$ is assumed to be both large and small, i.e. $t \gg \omega _b^{-1}$ but $t \ll (k_r \bar {v}_{r})^{-1}$ (implying $k_r \bar {v}_{r}/\omega _b \ll 1$), which, because $\bar {v}_{r} = 0$ is never exactly true in an actual stellarator, can only ever approximately be satisfied within a finite time interval, between $\omega _b^{-1}$ and $(k_r \bar {v}_{r})^{-1}$. The observation of $\phi _\mathrm {res}$ may therefore be a challenge in some cases, for instance in gyrokinetic simulations of stellarators for which these time scales are not well separated.

Obtaining the desired limit in our calculation is, however, a simpler matter, as we need only apply $k_r \bar {v}_{ra}/p \rightarrow 0$ to (2.15). The only remaining dependence on $p$ is the factor of $1/p$ on the source term, and the Laplace transform may be inverted to obtain

(3.2)\begin{gather} \sum_a \frac{e_a^2}{T_a}\int \,{\rm d}^3 v F_{a0}\left[\bar{\phi} b_a x_\perp^2/2 + \overline{\phi b_a x_\perp^2/2} - k_r^2 \overline{\phi\delta_{ra}}\delta_{ra} + k_r^2 \overline{\phi\delta_{ra}^2}/2 + k_r^2 \bar{\phi}\delta_{ra}^2/2\right]\nonumber\\ + \sum_a n_a\frac{e_a^2}{T_a}\left(\phi - \frac{1}{n_a} \int \,{\rm d}^3v F_{a0} \bar{\phi}\right) = S, \end{gather}

where the right-hand side denotes the source term (not yet expanded for small orbit width)

(3.3)\begin{equation} S = \sum_a e_a\int \,{\rm d}^3v {\rm J}_{0a}\overline{\delta F_a(0) {\rm e}^{{\rm i} k_r \delta_{ra}}} {\rm e}^{-{\rm i} k_r \delta_{ra}}. \end{equation}

Equation (2.15) can be compared with previous results in the small orbit width and small Larmor radius limits (for completeness, the result valid for arbitrary $k^2\delta _r^2$ and $b_i$ is given in Appendix B). We note differences coming from the fact that potential is kept under the bounce average, because we allow for $\phi \neq \langle \phi \rangle$, and the FLR terms take a similar form.

Let us consider what can be said about the general solution of (2.15). The final term on the left-hand side of (2.15) has the form of a linear operator on $\phi$, defined by

(3.4)\begin{equation} \sum_a n_a\frac{e_a^2}{T_a}\left(\phi - \frac{1}{n_a} \int \,{\rm d}^3v F_{a0} \bar{\phi}\right) = \sum_a n_a\frac{e_a^2}{T_a} {\mathcal{L}} \phi , \end{equation}

which is zero if and only if $\phi = \langle \phi \rangle$; see Appendix E. As a consequence, this operator is invertible on the non-zonal part of the potential, $\delta \phi$, defined by the following:

(3.5)\begin{equation} \phi = \delta\phi + \varPhi, \end{equation}

where $\varPhi = \langle \phi \rangle$. We may formally expand both $\varPhi = \varPhi ^{(0)} + \epsilon \varPhi ^{(1)} + \cdots$, $\delta \phi = \delta \phi ^{(0)} + \epsilon \delta \phi ^{(1)} + \cdots$ and $\delta F_a = \delta F_{a}^{(0)} + \cdots$ in our small parameter ($\epsilon \sim b_i \sim k_r^2\delta _{ri}^2$); it will not be necessary to keep these extra superscripts in what follows because we will only use zeroth quantities in our final expressions. With this expansion, the dominant contributions to (2.15) come from the right-hand side and the term $\sum _a n_a({e_a^2}/{T_a}) {\mathcal {L}} \phi$ on the left-hand side. The resulting equation can be formally solved for $\delta \phi ^{(0)}$ by use of the inverse ${\mathcal {L}}^{-1}$, yielding the non-zonal part of the residual potential

(3.6)\begin{equation} \delta\phi_\mathrm{res} = \left(\sum_a n_a\frac{e_a^2}{T_a}\right)^{-1} {\mathcal{L}}^{-1} \left[S_0 - \left\langle S_0\right\rangle\right]\!,\end{equation}

where $S_0 = \sum _a e_a\int \,{\rm d}^3v \overline {\delta F_{a}^{(0)}(0)}$. This equation implies that the residual potential essentially derives its non-zonal component from non-uniformity of initial charge distribution on the surface (or to be more precise, the charge density of the transit average of the gyro-average of the initial distribution functions). The main conclusion here, which may or may not be surprising, is that this component does not in fact decay away to zero.

At next order in our expansion, (2.15), the contributions from $\delta \phi ^{(1)}$ appearing under the operator ${\mathcal {L}}$ are eliminated by flux-surface average, leaving an equation for the zonal part of the potential, $\varPhi ^{(0)}$

(3.7)\begin{align} \varPhi_\mathrm{res} = \frac{\left\langle b_i \phi(0)\right\rangle +\tilde{\varPhi}_S - (2n_i)^{-1}\left\langle\int \,{\rm d}^3 v F_{i0}\left[\overline{\delta\phi_{{\rm res}}} b_i x_\perp^2 + \overline{\delta\phi_{{\rm res}} b_i x_\perp^2} + k_r^2 \overline{\delta\phi_{{\rm res}}\delta_{r}^2} + k_r^2 \overline{\delta\phi_{{\rm res}}}\delta_{r}^2\right]\right\rangle}{\left\langle b_i \right\rangle + n_i^{-1}\left\langle\int \,{\rm d}^3v F_{i0}k_r^2\delta_{r}^2 \right\rangle},\end{align}

where we use $\delta _{re} \ll \delta _{ri} \equiv \delta _{r}$ (dropping indices), $\rho _e \ll \rho _i$, the identity $\langle \bar {f}\rangle = \langle\, f \rangle$ (D2), $\overline {\delta _{ra}} = 0$ (due to oddness in $\sigma$ for trapped particles, and by choice of convention for passing particles) and quasi-neutrality for the initial condition. The latter can be written at each order in terms of the $n$th-order initial condition $\delta F_{a}^{(n)}(0)$, but these details are left for Appendix C. Additional contributions from the source at this order are included in the term

(3.8)\begin{align} \tilde{\varPhi}_S = \frac{T_i}{2n_ie_i} \left\langle\int \,{\rm d}^3v \left(2{\rm i} \overline{k_r \delta_{r}\delta F_i(0)} + \left(\delta F_i(0)-\overline{\delta F_i(0)}\right) b_i x_\perp^2 - k_r^2 \overline{\delta F_i(0)\delta_{r}^2} - k_r^2 \overline{\delta F_i(0)}\delta_{r}^2 \right)\right\rangle. \end{align}

Note the mixed orders of the terms: here, it is possible to consider, just for example, contributions to the initial condition $\delta F_i(0) \sim {O}(\delta _r k_r)$ that are odd in $v_{\|}$, as done by Rosenbluth & Hinton (Reference Rosenbluth and Hinton1998), or make even contributions at order $\delta F_i(0) \sim {O}(1)$, for instance due to pressure perturbations. However, we emphasize that such cases are generally inconsistent with the assumption that the residual is zonal, and the rather unwieldy expression of (3.7) must then be considered to determine the residual. The conditions under which this general solution prevails depend on details of the turbulence, and this will be discussed more later.

We note that a closed form expression for the zonal part of the residual ($\varPhi _\mathrm {res}$) in terms of the source may be obtained by substituting the solution for $\delta \phi _\mathrm {res}$, (3.6), into (3.7); we do not do this here as the resulting expression is not enlightening, especially without an explicit form of ${\mathcal {L}}^{-1}$.

3.1 Recovering the residual zonal flow

If the contribution to the charge from the initial gyro-centre distribution, i.e. the right-hand size of (2.15), is constant on a flux surface (to zeroth order), then the $\delta \phi$ term must balance with the small (FLR, FOW) terms, and we can conclude

(3.9)\begin{equation} \delta\phi \sim {O}(b_i \varPhi), \end{equation}

and $\delta \phi$ can be safely neglected in (3.7); those from $\tilde {\varPhi }_S$ can also be neglected if $\delta F_i(0) \sim {O}(b_i)$ is assumed. Solving the equation for $\varPhi$ we then obtain the residual

(3.10)\begin{equation} \phi_\mathrm{res} = \frac{\left\langle b_i \phi(0)\right\rangle}{\left\langle b_i \right\rangle + n_i^{-1}\left\langle\int \,{\rm d}^3v F_{i0}k_r^2\delta_{r}^2 \right\rangle}. \end{equation}

We see that, even in this limit, we do not exactly recover the result of Rosenbluth & Hinton (Reference Rosenbluth and Hinton1998), as we do not assume the initial potential to be zonal. The Rosenbluth–Hinton (RH) result can be written in our notation as follows:

(3.11)\begin{equation} \phi_\mathrm{res}^\mathrm{RH} = \varPhi(0)\frac{\left\langle b_i\right\rangle}{\left\langle b_i\right\rangle + n_{i}^{-1}\left\langle\int \,{\rm d}^3v F_{i0}k_r^2\delta_r^2 \right\rangle}, \end{equation}

where we have used that Rosenbluth & Hinton (Reference Rosenbluth and Hinton1998) assumed the initial potential to be zonal to write this result in terms of $\phi (0) = \langle \phi (0)\rangle = \varPhi (0)$. Evidently, there is one limit (for arbitrary $\phi (0)$) in which the RH result is obtained from (3.10), which is when $\langle b_i \rangle \approx b_i$, e.g. for the circular tokamak model.

3.2 Dependence of residual zonal flow on initial potential

More generally, we note that there is a class of initial conditions consistent with (3.9), including that traditionally assumed in calculations of the residual, $\phi (0) = \langle \phi (0)\rangle = \varPhi (0)$. Indeed, allowing non-zonal $\phi (0)$, (3.10) exhibits a certain variation in what can be obtained for the ratio $\phi _\mathrm {res}/\varPhi (0)$, as compared with the RH expression, i.e.

(3.12)\begin{equation} \frac{\phi_\mathrm{res}}{\phi_\mathrm{res}^\mathrm{RH}} = \frac{\left\langle b_i \phi(0)\right\rangle}{\left\langle b_i\right\rangle\left\langle \phi(0)\right\rangle}, \end{equation}

which arises only from the variation in $b_i$. An interesting and possibly useful case is that of initially zonal distribution functions (a convenient way to initialize a gyrokinetic simulation and therefore a good test case). In this case the initial charge is also zonal and from (C4) we have

(3.13)\begin{equation} \left\langle b_i \phi(0)\right\rangle = b_i \phi(0).\end{equation}

Dividing by $b_i$ and averaging we can solve for this for the zonal potential, $\varPhi (0) = \langle \phi (0)\rangle$, and obtain

(3.14)\begin{equation} \phi_\mathrm{res} = \varPhi(0) \frac{\left\langle b_i^{-1} \right\rangle^{-1}}{\left\langle b_i \right\rangle + n_i^{-1}\left\langle\int \,{\rm d}^3v F_{i0}k_r^2\delta_{r}^2 \right\rangle}. \end{equation}

Because of the inequality between the harmonic and arithmetic means, we find that the residual expressed by (3.14) is less than or equal to the RH expression, in particular $\phi _\mathrm {res}/\varPhi (0) \leq \phi _\mathrm {res}^\mathrm {RH}/\varPhi (0)$, with equality only in the case of uniform $b_i$.

4 The residual in stellarators and tokamaks

Having demonstrated how to calculate the intermediate residual for stellarators, the natural question arises about how different stellarators fare with respect to this measure, how they compare with tokamaks and in particular whether anything can be said about the different classes of optimized stellarators.

4.1 Tokamaks and quasi-symmetric stellarators

The residual is inversely proportional to a weighted average of $\delta _r^2$ and will thus be particularly large in a field where the radial width of most particle orbits is small. In a standard large-aspect-ratio tokamak with circular cross-section – the case considered by Rosenbluth & Hinton (Reference Rosenbluth and Hinton1998) – circulating ion orbits have radial excursions of order $q \rho _i$ whereas trapped ones have larger banana orbits of width

(4.1)\begin{equation} \delta_r \sim \frac{q \rho_i}{\epsilon^{1/2}}, \end{equation}

where $\epsilon \ll 1$ denotes the inverse aspect ratio and $q=\iota ^{-1}$ the inverse rotational transform (Helander & Sigmar Reference Helander and Sigmar2002). Although the latter only constitute a small fraction $f_t \sim \epsilon ^{1/2} \ll 1$ of the total number of particles, they dominate the average of $\delta _r^2$, which becomes of order

(4.2)\begin{equation} \frac{1}{n_i} \int \,{\rm d}^3v F_{i0} \delta_r^2 \sim (1-f_t) q^2 \rho_i^2 + f_t \left( \frac{q \rho_i}{\epsilon^{1/2}} \right)^2 \sim \frac{q^2 \rho_i^2 }{\epsilon^{1/2}}, \end{equation}

and qualitatively explains the RH result

(4.3)\begin{equation} \phi_\mathrm{res}^\mathrm{RH} = \frac{\varPhi(0)}{1 + 1.64 q^2\epsilon^{-1/2}}. \end{equation}

Since, in a typical tokamak, the term $1.64 q^2\epsilon ^{-1/2}$ is considerably larger than unity, this residual is relatively weak.

In quasisymmetric stellarators, particle trajectories are similar those in tokamaks (Boozer Reference Boozer1983; Nührenberg & Zille Reference Nührenberg and Zille1988), and the calculation is therefore mathematically identical if the symbols are suitably re-interpreted. In quasi-axisymmetric stellarators, the orbit width is equal to that in a tokamak, and the residual is therefore given by an expression like (4.3), except that the numerical factor 1.64 needs to be adjusted if the magnetic field strength does not vary sinusoidally along the field. A similar adjustment is required in shaped tokamaks (Xiao & Catto Reference Xiao and Catto2006). In quasihelically symmetric stellarators, the banana orbit width

(4.4)\begin{equation} \delta_r \sim \frac{\rho_i}{|N - \iota|}, \end{equation}

is smaller than that in a tokamak by a factor $|N/\iota - 1|>1$. Such stellarators thus have a larger residual than quasi-axisymmetric ones and tokamaks.

4.2 Quasi-isodynamic stellarators

The smallest orbit widths, and thus the largest residuals, are realized in so-called quasi-isodynamic stellarators,Footnote ² which are omnigenous stellarators with poloidally closed contours of constant field strength. Such stellarators usually do not carry any significant amount of net toroidal current, and the magnetic field can be written as (Helander Reference Helander2014)

(4.5)\begin{equation} {\boldsymbol{B}} = G(\psi) \boldsymbol{\nabla} \varphi + K(\psi, \alpha, \varphi) \boldsymbol{\nabla} \psi, \end{equation}

in Boozer coordinates, and the radial drift velocity becomes

(4.6)\begin{equation} v_r = \frac{v_{\|}^2 + v_\perp^2/2}{\varOmega} ({\boldsymbol{b}} \times \boldsymbol{\nabla} \ln B) \boldsymbol{\cdot} \boldsymbol{\nabla} r =- \frac{v^2 r'(\psi)}{\varOmega} \left(1 - \frac{\lambda B}{2} \right) \left( \frac{\partial B}{\partial \alpha} \right)_{\psi,\varphi}, \end{equation}

where $\varOmega = eB/m$ denotes the gyrofrequency. Equation (2.4) can be solved for the radial excursion

(4.7)\begin{equation} \delta_r = \frac{1}{v} \int v_r \,{\rm d}t, \end{equation}

where we use $\bar {v}_r = 0$ and define ${\rm d}t = {\rm d}l/\sqrt {1-\lambda B}$ with $t$ a time-like variable along the orbit. The lower limit of integration can be chosen so that $\overline {\delta _r} = 0$. For magnetically trapped orbits, i.e. for values of $\lambda$ less than $1/B_{\rm max}$, this is achieved by choosing the lower integration limit to correspond to a bounce point and for passing orbits to the point of maximum field strength. In the latter case

(4.8)\begin{equation} \overline{\delta_r} = \left.\lim_{L\rightarrow \infty} \int_0^L \frac{\delta_r \,{\rm d}l}{\sqrt{1-\lambda B}} \right/ \int_0^L \frac{{\rm d}l}{\sqrt{1-\lambda B}} =\left. \left\langle \frac{B \delta_r}{\sqrt{1-\lambda B}} \right\rangle \right/ \left\langle \frac{B}{\sqrt{1-\lambda B}} \right\rangle, \end{equation}

vanishes thanks to the $\alpha$-derivative in (4.6), because the flux-surface average can be written as (Helander & Nührenberg Reference Helander and Nührenberg2009)

(4.9)\begin{equation} \langle \cdots \rangle = \frac{1}{V'} \int_0^{2{\rm \pi}} \,{\rm d}\alpha \int_0^L\frac{{\rm d}l}{B} (\cdots ) , \end{equation}

where

(4.10)\begin{equation} V'(\psi) = \int_0^{2{\rm \pi}} \,{\rm d}\alpha \int_0^L \frac{{\rm d}l}{B}. \end{equation}

Here, an integral over the arc length $l$ is taken over one period of the device from one maximum of $B$ to the next, and the distance $L$ between them is independent of $\alpha$, which is true for perfectly quasi-isodynamic fields.

We shall not endeavour to calculate the zonal-flow response (3.10) explicitly but rather show that it is relatively large in quasi-isodynamic stellarators by finding an upper bound on the quantity

(4.11)\begin{equation} D = \frac{1}{n} \left\langle \int \,{\rm d}^3v F_{0} \delta_r^2 \right\rangle = \frac{1}{nV'} \int_0^{2{\rm \pi}} \,{\rm d}\alpha \int_0^L \frac{{\rm d}l}{B} \int_0^\infty F_0 2 {\rm \pi}v^2 \,{\rm d}v \int_0^{1/B} \frac{\delta_r^2 B \,{\rm d}\lambda}{\sqrt{1-\lambda B}}. \end{equation}

By interchanging the integrals over $\lambda$ and $l$, and replacing the latter by $t$, we find

(4.12)\begin{equation} D = \frac{1}{nV'} \int_0^\infty F_0 2 {\rm \pi}v^2 \,{\rm d}v \int_0^{2{\rm \pi}} \,{\rm d}\alpha \int_0^{1/B_{\rm min}} \,{\rm d}\lambda \int \delta_r^2 \,{\rm d}t, \end{equation}

where the $t$-integral is taken over the region where $1-\lambda B$ is positive, i.e. over the entire range for passing orbits and over the magnetic trapping well(s) for trapped ones. At the end points of each $t$-integral, the function $\delta _r$ vanishes. As a consequence of the Poincaré inequality

(4.13)\begin{equation} \int_0^{\tau_b} g^2(t) \,{\rm d}t \le \frac{\tau_b^2}{{\rm \pi}^2} \int_0^{\tau_b} \left( \frac{{\rm d}g}{{\rm d}t}\right)^2 \,{\rm d}t, \end{equation}

for functions such that $g(0) = g(\tau _b) = 0$, we thus conclude that $D$ is bounded from above by

(4.14)\begin{equation} D \le \frac{2}{{\rm \pi} nV'} \int_0^\infty F_0 \,{\rm d}v \int_0^{2{\rm \pi}} \,{\rm d}\alpha \int_0^{1/B_{\rm min}} \tau_b^2 \,{\rm d}\lambda \int v_r^2 \,{\rm d}t, \end{equation}

where

(4.15)\begin{equation} \tau_b(\lambda) = \int_{\lambda B(l) < 1} \frac{{\rm d}l}{\sqrt{1- \lambda B(l)}}. \end{equation}

Substituting (4.6) finally results in the rigorous inequality

(4.16)

The integrals in this expression depend on details in the spatial variation in the magnetic field strength, but we note that, generally, the $\lambda$-integral is of order $L^2/B$ and the $l$-integral of order $\epsilon ^2 L$, where $\epsilon$ denotes the relative poloidal variation of $B$ at constant $\varphi$. Since ${\rm d} \psi / {\rm d}r \sim rB$ and $V' \sim 2 {\rm \pi}L/B$, we thus obtain

(4.17)\begin{equation} D \lesssim \frac{3}{2 {\rm \pi}^2} \left( \frac{\epsilon \rho_i L}{r} \right)^2. \end{equation}

In a quasi-isodynamic stellarator, the level curves of constant magnetic field strength close poloidally, rather than toroidally, on each flux surface. The field strength varies along the magnetic axis, where it is a function only of $\varphi$, and in its vicinity, the quantity $\epsilon = \partial \ln B / \partial \alpha = \partial \ln B / \partial \theta$ appearing in (4.17) is thus small. In the typical large-aspect-ratio scenario, we can estimate $\epsilon \sim r\kappa$, where $\kappa$ is the curvature of the magnetic axis; see for example Plunk, Landreman & Helander (Reference Plunk, Landreman and Helander2019). A conservative bound for this curvature is $\kappa \lesssim 1/L$ (optimization may achieve somewhat lower values), yielding $\epsilon \lesssim r/L$. We therefore expect from (4.17) that $D \lesssim \rho _i^2$ for a quasi-isodynamic stellarator, and the residual (3.10) is therefore comparable to the initial perturbation, i.e. much larger than in a tokamak.