The residual flow in well-optimized stellarators

Abstract

The gyrokinetic theory of the residual flow, in the electrostatic limit, is revisited, with optimized stellarators in mind. We consider general initial conditions for the problem, and identify cases that lead to a non-zonal residual electrostatic potential, i.e. one having a significant component that varies within a flux surface. We investigate the behaviour of the ‘intermediate residual’ in stellarators, a measure of the flow that remains after geodesic acoustic modes have damped away, but before the action of the slower damping that is caused by unconfined particle orbits. The case of a quasi-isodynamic stellarator is identified as having a particularly large such residual, owing to the small orbit width achieved by optimization.

1 Introduction

The work of Rosenbluth & Hinton (Rosenbluth & Hinton 1998; Hinton & Rosenbluth 1999) established the idea of an undamped ‘residual’ potential in tokamaks. The idea is to initialize an electrostatic potential, varying only in the radial direction, and track its value over a time much longer than the transit period over which particles move along the magnetic field lines. Initially, the potential oscillates and diminishes in amplitude, due to geodesic acoustic mode activity, but eventually a steady residual signal emerges. Because all of this happens due to a collisionless dynamics, it was argued that full gyrokinetics is needed to properly model the turbulence, which remains the prevailing attitude to this day.

The residual proved very popular in part due to the fact that it admits exact predictions, commonly used to benchmark gyrokinetic codes. However, the question of what such calculations imply for the behaviour of fully developed turbulence, where a strongly nonlinear dynamics prevails, is a difficult one to grapple with. Part of the difficulty is that the dynamics of the driven potential cannot be easily separated from that of the turbulence that drives it. Indeed, modelling the turbulence as a steady source leads to unbounded growth of the residual (Rosenbluth & Hinton 1998). Other basic questions arise as to whether the residual is relevant in cases where the damped solutions called geodesic acoustic modes (GAMs) may be effectively driven by the turbulence (Waltz & Holland 2008).

There is a vast body of work concerning zonal flows in magnetic fusion (see for instance Diamond et al. (2005) for a start), with many simple limits having been considered theoretically, and a multitude of sometimes disparate models and explanations having been proposed. There is, however, broad agreement that zonal flows have a beneficial influence on turbulence, lowering fluctuation levels by shearing turbulent eddies (Hahm et al. 1999), promoting transport of energy to small scales and inducing coupling between unstable and stable eigenmodes (Makwana, Terry & Kim 2012). Zonal flows are also responsible for the Dimits shift, whereby turbulence is all but eliminated for a finite range above the threshold of the linear instability (Dimits et al. 2000; Rogers, Dorland & Kotschenreuther 2000; St-Onge 2017; Pueschel, Li & Terry 2021; Hallenbert & Plunk 2022). In the context of stellarators, it is clear that the strength of the geodesic curvature, related to both GAM damping and residual levels, has a strong effect on the overall turbulence levels (Xanthopoulos et al. 2011), and variation in the linear response of zonal flows is believed to underlie confinement differences between different configurations of the Large Helical Device (LHD) (Watanabe, Sugama & Ferrando-Margalet 2008). The neoclassical radial electric field in a stellarator can also stabilize turbulence via shearing etc., and is experimentally associated with enhanced confinement (see for instance Lore et al. 2010), but its origins are distinct from zonal flows, and there is not yet evidence to support such a role in the W-7X stellarator (Xanthopoulos et al. 2020). Overall, there is a strong motivation to deepen the understanding of the theoretical foundations of zonal flows in stellarators, of which the residual is a key part, especially to aid in the design of future devices.

On the most fundamental theoretical level, one observes that, according to gyrokinetic theory, the entire ‘$k_\alpha = 0$ component’ of the fluctuations is stable (${\boldsymbol {k}}_\perp = k_\psi \boldsymbol {\nabla } \psi + k_\alpha \boldsymbol {\nabla } \alpha$ is the wavenumber perpendicular to the magnetic field, where $\alpha$ and $\psi$ are defined such that ${\boldsymbol {B}} = \boldsymbol {\nabla }\psi \times \boldsymbol {\nabla }\alpha$ is the magnetic field), having no source of free energy. That is, in addition to the zonal potential, which by definition is constant within a flux surface, there are also components that vary within a flux surface by virtue of smooth dependence along the field line, including the electrostatic potential and other moments such as parallel ion flow or temperature perturbations. What is the fate of these components, and do they also have a role in regulating the turbulence?

In the present work we focus on the residual in the context of stellarators instead of tokamaks (although we note that our findings apply equally to the latter). The above issues also arise in a stellarator, and it is even more unclear for which cases the residual is a useful quantity as a predictor of the turbulence. Indeed, as was found by Mishchenko, Helander & Könies (2008) and Helander et al. (2011a), the residual in a stellarator is strongly affected by the presence of unconfined particle orbits, i.e. trapped particles bouncing back and forth in magnetic wells that drift radially between magnetic surfaces; see also Monreal et al. (2016). If this drift is non-zero, however small, it was found that the residual is sharply reduced, an apparent strike against the stellarator.

However, in cases where the radial drift is sufficiently small, such as optimized stellarators, the additional stellarator-specific damping must necessarily act on a time scale much longer that the drive and saturation of the turbulence. For such cases, we identify an ‘intermediate’ residual, namely, the solution that is found at times much longer than the transit time and much smaller than the time scale of the stellarator-specific damping. To calculate this residual, we revisit the initial value problem of gyrokinetics, this time without assuming the initial condition or the final state to be well approximated by a zonal potential. Indeed, we find that the residual is generally non-zonal, but we do identify conditions under which such assumptions are valid.

This calculation casts the optimized stellarator in a more positive light, showing that it can exhibit a large residual, especially in the case that the widths of particle orbits are small, which, as we demonstrate, is especially true in a particular class, so-called quasi-isodynamic stellarators. Indeed, in such stellarators the residual is found to be much larger than that in tokamaks, which could have significant consequences for the regulation of plasma turbulence.

2 Gyrokinetic solution of the initial value problem

We will follow some of the notation conventions of Helander et al. (2011a), 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. (2011a), 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$.¹ We also define the flux-surface average (Helander 2014)

(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. (2011a) 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. (2016), 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. (2011a) (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. 2008; Helander et al. 2011a; Monreal et al. 2017). 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. 2011a). 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. (2011a), 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 (1998), 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 (1998), 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 (1998), 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 (1998) 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 (1998) – 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 2002). 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 1983; Nührenberg & Zille 1988), 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 2006). 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,² 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 2014)

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

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

5 Conclusions

Although the long-time-asymptotic residual potential is expected to be small in a stellarator (Mishchenko et al. 2008; Helander et al. 2011a), we have argued that a well-optimized stellarator exhibits a larger effective residual on time scales important for turbulence. To assess this ‘intermediate’ residual, we have revisited the general initial value problem, allowing for arbitrary initial condition, and derive the resulting form of the residual, whether zonal or non-zonal. We identify two cases (in the limit of small orbit width and Larmor radius) depending on the charge induced by the double-orbit-averaged (bounce- and gyro-averaged) initial distribution functions, which can be described as follows.

If this charge is zonal (constant on a flux surface), we find that the residual potential is also zonal, and depends only on the initial potential, i.e. it is insensitive to other details of the initial distribution functions. In that case, we note that a large ‘intermediate’ residual is indeed possible in stellarators, even in cases when the ‘true’ time-asymptotic residual is negligibly small. It is argued that the intermediate residual will be largest in quasi-isodynamic stellarators, smaller in quasi-helically symmetric stellarators and smallest in tokamaks and quasi-axisymmetric stellarators. This may be counter-intuitive to some readers, as it is known that undamped equilibrium flows can be sustained on time scales exceeding the ion collision time in quasi-symmetric stellarators but not in quasi-isodynamic ones (Helander & Simakov 2008; Helander, Geiger & Maaßberg 2011b), but there is a distinction between those equilibrium flows and the small-scale zonal flows that arise spontaneously with micro-turbulence. For the latter flows, which regulate turbulence at low collisionality, it is the quasi-isodynamic stellarator that performs the best. These stellarators exhibit a much larger residual than tokamaks and quasi-axisymmetric stellarators, where the RH factor $1.6 q^2 \epsilon ^{-1/2}$ substantially exceeds unity. The collisional damping that occurs in quasi-isodynamic fields (but not in quasisymmetric ones) due to the lack of intrinsic ambipolarity only takes place on the longer time scale of ion collisions.

Formal complications arise in the calculation of the residual when the charge induced by the initial condition has a significant non-zonal component. We show in this case that the residual potential is non-zonal, i.e. varies in the flux surface, and generally depends on details of the initial distribution functions. We work out the general form of the complete solution, leaving its more detailed analysis for later, but note that its understanding may allow the consideration of a broader class of nonlinear drives, in other words a more general source for the $k_\alpha = 0$ component.

Although we derive the source ($S$) of the residual from the initial condition, i.e. a delta function in time, the actual source in the gyrokinetic equation is the nonlinear term, which provides free energy to the $k_\alpha = 0$ component. In the absence of a fully nonlinear theory describing the steady-state dynamics of the stable and unstable components of the turbulence, it is reasonable, in interpreting the result of the residual calculation, to consider what form a realistic turbulent source might take.

One possibility is that the source has a significant non-zonal temperature perturbation, as expected from secondary instability theory (Plunk & Navarro 2017), one mechanism by which zonal flows may be driven, which predicts that the temperature perturbation is both non-zonal $\delta T \neq \langle \delta T \rangle$ and acquires its size and spatial dependence from the instabilities (ion temperature gradient) that drive it. Other hints can be obtained directly from turbulence simulations. Although the perpendicular temperature of the $k_\alpha = 0$ components is generally observed to be small (Rogers et al. 2000), it is also observed to grow in relative amplitude at strong drive (Plunk, Navarro & Jenko 2015). It is therefore unclear whether the temperature perturbations of our source should be expected to be large enough to drive a strongly non-zonal potential ($\delta T/(e_i\phi ) \sim 1$), but it seems unlikely that they will always be so small that the non-zonal part of the residual can be assumed asymptotically small ($\delta T/(e_i\phi ) \sim b_i$). Similar questions also apply concerning other components of the source, such as parallel ion flow, and these will all have to be explored further in the future.

The work leaves ample opportunity for further studies, especially involving gyrokinetic simulations. Fully nonlinear simulations of the turbulence may help identify cases where the non-zonal solutions described by (3.6)–(3.7) may arise. On a more basic level, linear initial value simulations should also be conducted to verify the quantitative validity of these expressions, especially for recently found designs that satisfy the quasi-isodynamic condition to high precision (Goodman et al. 2023). It should be noted that the predictions of this work apply also for the much simpler context of tokamak geometry. In particular, (3.10) gives a prediction for the residual when the initial condition is non-zonal, which may already be tested for simple model tokamak geometries with spatially varying flux compression ($|\boldsymbol {\nabla } \psi |$). The inequality derived to bound the residual, (4.16), might also be further investigated in some limits, and it, along with related estimates, could prove useful in stellarator optimization for reduced turbulence.