2 Gyrokinetic system

To begin, we consider a local, electrostatic gyrokinetic system, where the domain is a narrow flux tube centred around a magnetic field line in the plasma. In this local approximation, a Fourier decomposition may be applied in the directions perpendicular to the magnetic field assuming periodicity in these coordinates. The electrostatic gyrokinetic equation for a species $a$, and the Fourier component $\boldsymbol {k}$ is

(2.1)\begin{align} & \frac{\partial g_{a, \boldsymbol{k}}}{\partial t} + v_{\|}\frac{\partial g_{a, \boldsymbol{k}}}{\partial l} + {\rm i}\omega_{{\rm da}} g_{a, \boldsymbol{k}}+ \frac{1}{B^2}\sum_{\boldsymbol{k}'}\boldsymbol{B}\boldsymbol{\cdot} (\boldsymbol{k}\times \boldsymbol{k}')\delta\phi_{\boldsymbol{k'}}{\rm J}_{0a} \,g_{a,\boldsymbol{k}-\boldsymbol{k}'} \nonumber\\ & \quad = \frac{e_a F_{a0}}{T_a}\left(\frac{\partial}{\partial t} + {\rm i}\omega_{*a}^T\right)\delta \phi_{\boldsymbol{k}} {\rm J}_{0a}, \end{align}

where $g_{a,\boldsymbol {k}}$ is the non-adiabatic part of the perturbed distribution function $g_a(\boldsymbol {R}, E_a, \mu _a, t) = \delta f(\boldsymbol {r}, \boldsymbol {v}, t) + e_a\delta \phi (\boldsymbol {r}, t) F_{a0}/T_a$, $F_{a0}$ is the background Maxwellian distribution and $\delta \phi$ is the electrostatic potential. The phase-space coordinates for the gyrocentre distribution function $g_a$ are the magnetic moment, $\mu _a = m_a v_\perp ^2/(2B)$ and the kinetic energy, $E_a = m_a v^2/2$. We also use Clebsch magnetic coordinates, in which the magnetic field is of the form $\boldsymbol {B} = \boldsymbol {\nabla } \psi \times \boldsymbol {\nabla } \alpha$, where $\psi$ is the toroidal magnetic flux (acting as a radial coordinate in toroidal geometry), and $\alpha$ is the binormal coordinate, which labels different field lines on surfaces of constant $\psi$. In these coordinates, the field-line following coordinate is $l$, the unit vector along the magnetic field is $\boldsymbol {b} = \boldsymbol {B}/B$ and the perpendicular wavenumber is ${\boldsymbol{k}} ={\boldsymbol{k}}_\perp = k_\alpha \nabla \alpha + k_\psi \nabla \psi$. The argument of the Bessel function is ${\rm J}_0 = {\rm J}_0(k_\perp v_\perp /\varOmega _a)$ with the gyrofrequency defined as $\varOmega _a = e_a B/m_a$. The magnetic drift frequency is $\omega _{da}= \boldsymbol {k}_\perp \boldsymbol {\cdot } \boldsymbol {v}_d$, where $\boldsymbol {v}_d = {\boldsymbol {b}}\times (v_\perp ^2/2 \boldsymbol {\nabla } \ln B + v_\|^2 \boldsymbol {\kappa })/{\varOmega _a}$. In the low-plasma pressure limit (assuming $\boldsymbol {\nabla } \ln B \approx \boldsymbol {\kappa }$), the drift frequency can be approximated by

(2.2)\begin{equation} \omega_{{\rm da}} \approx \hat{\omega}_{{\rm da}}(l)\left( \frac{v_\perp^2}{2 v_{{\rm Ta}}^2} + \frac{v_\|^2}{v_{{\rm Ta}}^2}\right), \end{equation}

where $v_{{\rm Ta}} = \sqrt {2T_a/m_a}$ is the thermal speed and $\hat {\omega }_{{\rm da}}(l)$ is a geometry dependent factor. The energy-dependent diamagnetic drift frequency is

(2.3)\begin{equation} \omega^T_a = \omega_{*a}\left[1 + \eta_a\left(\frac{v^2}{v_{{\rm Ta}}^2} - \frac{3}{2}\right)\right], \end{equation}

where the plasma gradients enter via $\omega _{*a} = (k_\alpha T_a/e_a)\,\mathrm {d} \ln n_a /\mathrm {d} \psi$ and the gradient ratio $\eta _a = \mathrm {d}\ln T_a/ \mathrm {d} \ln n_a$. In the electrostatic limit, the system is closed by the quasineutrality condition

(2.4)\begin{equation} \sum_a \frac{e_a^2 n_a}{T_a}\delta \phi_{\boldsymbol{k}} = \sum_a e_a\int g_{a,\boldsymbol{k}} {\rm J}_{0a} \, \mathrm{d}^3v. \end{equation}

In Clebsch coordinates, we may write the nonlinear term of (2.1) in the more convenient form

(2.5)\begin{equation} \frac{1}{B^2}\sum_{\boldsymbol{k}'}\boldsymbol{B} \boldsymbol{\cdot} (\boldsymbol{k}\times \boldsymbol{k}')\delta\phi_{\boldsymbol{k'}}{\rm J}_{0a} \ g_{a,\boldsymbol{k}-\boldsymbol{k}'} = \sum_{\boldsymbol{k}'} \left(k_\psi k_\alpha' -k_\alpha k_\psi' \right)\delta\phi_{\boldsymbol{k'}}{\rm J}_{0a} \, g_{a,\boldsymbol{k}-\boldsymbol{k}'}. \end{equation}

3 Bounce-averaged electrons

To simplify the electron gyrokinetic equation, we employ a standard ordering. We assume that the time scale of the dynamics of interest in the gyrokinetic system $\tau _D$ (i.e. the growth time of an instability) is much longer than the time scale on which thermal electrons transit the typical parallel length scale of the system, such that $\tau _D \gg L/v_{{\rm Te}}$. Expanding the electron distribution in orders of $L/(v_{{\rm Te}}\tau _D)$ as $g_{e,\boldsymbol {k}} = g_{e,\boldsymbol {k}}^0 + g_{e,\boldsymbol {k}}^1 \ldots$ we find that, to leading order,

(3.1)\begin{equation} v_\parallel \frac{\partial g_{e,\boldsymbol{k}}^0}{\partial l} = 0, \end{equation}

such that the lowest-order electron distribution function is constant along the magnetic field line. This, it turns out, is a highly consequential result of this ordering, and how it is interpreted is the key to tailoring the forthcoming optimal mode analysis to capture the behaviour of the electron response for ion-scale instabilities that satisfy our ordering in toroidal geometry.

The constancy of $g_{e}^0 = g_e^0(E_e,\mu _e, \psi, \alpha )$ in $l$ implies that its value is entirely determined by the parallel boundary conditions. For the passing particles in toroidal geometry, and for $k_\alpha \neq 0$, the relevant boundary conditions are the so-called ‘incoming’ boundary conditions of ballooning space where

(3.2)\begin{gather} g_{e,\boldsymbol{k}}(v_\parallel{>} 0, l ={-}\infty) = 0, \end{gather}

(3.3)\begin{gather}g_{e,\boldsymbol{k}}(v_\parallel{<} 0, l = \infty) = 0. \end{gather}

As explained in an appendix by Plunk et al. (Reference Plunk, Helander, Xanthopoulos and Connor2014), these boundary conditions are required to preserve causality (Connor, Hastie & Taylor Reference Connor, Hastie and Taylor1980), such that information from the boundary, which is infinitely far away, cannot be communicated in a finite time. Therefore, it must be the case that $g_{e}^0(k_\alpha \neq 0) = 0$ in the passing region of velocity space to satisfy the boundary conditions. For the $k_\alpha = 0$ component, ballooning space boundary conditions are not appropriate because there is no mechanism that should cause the mode to decay along the flux tube. Thus, for $k_\alpha = 0$ the electron distribution function has a zeroth-order contribution from the passing electrons.

The trapped particles on the other hand (in the absence of collisions), are confined to a limited region of ballooning space. As such, the appropriate boundary conditions for these particles are the ‘turning-point’ boundary conditions (Connor et al. Reference Connor, Hastie and Taylor1980),

(3.4)\begin{gather} g_{e,\boldsymbol{k}}(v_\parallel{>} 0, l = l_1) = g_{e,\boldsymbol{k}}(v_\parallel{<} 0, l = l_1), \end{gather}

(3.5)\begin{gather}g_{e,\boldsymbol{k}}(v_\parallel{>} 0, l = l_2) = g_{e,\boldsymbol{k}}(v_\parallel{<} 0, l = l_2) \end{gather}

where $l_{1,2}$ are the bounce points for a particle defined by $\lambda B(l_{1,2}) =1$, where $\lambda$ is the particle pitch angle, $\lambda = v_\perp ^2/(v^2 B)$. These boundary conditions can be satisfied with $g_{e,\boldsymbol {k}}^0 \neq 0$, and so, for $k_\alpha \neq 0$, only the trapped electrons contribute to the non-adiabatic electron response in this limit. Note that this is not the case in closed-field-line geometries where the parallel boundary conditions can be periodic, allowing $g_{e,\boldsymbol {k}}^0$ to be non-zero in the passing region of velocity space to zeroth order.

At next order in $L/(v_{\rm Te}\tau _D)$, we find

(3.6)\begin{align} & \frac{\partial g_{e,\boldsymbol{k}}^0}{\partial t} + v_\parallel \frac{\partial g_{e,\boldsymbol{k}}^1}{\partial l} + {\rm i} \omega_{{\rm de}}g_{e,\boldsymbol{k}}^0 + \sum_{\boldsymbol{k}'} \left(k_\psi k_\alpha' -k_\alpha k_\psi' \right)\delta\phi_{\boldsymbol{k}'}{\rm J}_{0e} g_{e,\boldsymbol{k}-\boldsymbol{k}'}^0 \nonumber\\ & \quad={-}\frac{e F_{e0}}{T_e}\left(\frac{\partial}{\partial t} +{ \rm i}\omega_{*e}^T\right)\delta \phi_{\boldsymbol{k}} {\rm J}_{0e}. \end{align}

We now define the bounce average as

(3.7)\begin{equation} \overline{(\ldots)} = \frac{\displaystyle\int\nolimits_{l_1}^{l_2} (\ldots)\,{\mathrm{d}l}/{v_\|}}{\displaystyle\int\nolimits_{l_1}^{l_2}\,{\mathrm{d}l}/{v_\|}} = \frac{1}{\tau_B}\int_{l_1}^{l_2} \frac{(\ldots) \,\mathrm{d}l}{\sqrt{1 -\lambda B}}, \end{equation}

where we have introduced the bounce time $\tau _B = \int _{l_1}^{l_2}({\mathrm {d}l}/{\sqrt {1 -\lambda B}})$.Footnote ¹ In pitch-angle coordinates, the velocity space element is defined as

(3.8)\begin{equation} \mathrm{d}^3 v = \sum_\sigma \frac{{\rm \pi} v^2 B \, \mathrm{d} v \, \mathrm{d}\lambda}{\sqrt{1 -\lambda B}}, \end{equation}

where $\sigma = v_\parallel /|v_\parallel |$. Applying the bounce average to (3.6) annihilates the $g_{e, \boldsymbol {k}}^1$ contribution, leaving an equation for $g^0_{e,\boldsymbol {k}}$. Additionally ordering $k_\perp \rho _e\ll 1$, where $\rho _a = v_{{\rm Ta}}/|\varOmega _a|$, for simplicity and omitting the superscript ‘$0$’ for brevity gives

(3.9)\begin{equation} \frac{\partial g_{e,\boldsymbol{k}}}{\partial t} + {\rm i} \overline{\omega}_{{\rm de}}g_{e,\boldsymbol{k}}+ \sum_{\boldsymbol{k}'} \left(k_\psi k_\alpha' -k_\alpha k_\psi' \right)\overline{\delta \phi}_{\boldsymbol{k'}} \, g_{e,\boldsymbol{k}-\boldsymbol{k}'} ={-}\frac{e F_{e0}}{T_e}\left(\frac{\partial}{\partial t} + {\rm i}\omega_{*e}^T\right)\overline{\delta \phi}_{\boldsymbol{k}}. \end{equation}

We now consider (3.9) as the governing equation for the bounce-averaged electrons in our gyrokinetic system, which is accurate to lowest order in $L/(v_{{\rm Te}}\tau _D)$. Here, we restrict our view to toroidal geometry such that, for $k_\alpha \neq 0$, $g_e$ is only non-zero in the trapped region of velocity space.

4 Energy balance

In the limit of $\tau _D \gg L/v_{\rm Te}$, we have a gyrokinetic system comprised of ions described by (2.1), with finite ion-Larmor-radius effects, and bounce-averaged electrons described by (3.9). As described in Parts 1, 2 and 3, we can construct energy norms for this two-species system that are conserved by the nonlinear terms of the gyrokinetic equations. These nonlinear invariants are the Helmholtz free energy, detailed in Parts 1 and 2, and the electrostatic energy, described in Part 3. We now derive these invariants in our system with bounce-averaged electrons.

4.1 Helmholtz free energy

For simplicity, we will consider a two-species hydrogen plasma with $n_i = n_e = n$, but the following can be easily generalised to multiple ion species. Constructing the Helmholtz free energy of this system proceeds in the same manner as in the general multispecies case, shown in Part 2, except for our assumption that $k_\perp \rho _e\ll 1$. This is done by multiplying the gyrokinetic equation for each species by the factor $T_a g_{a,\boldsymbol {k}}^*/F_{a0}$, integrating over velocity space, applying the flux-tube average defined by

(4.1)\begin{equation} \langle \ldots \rangle = \left.\lim_{L\to \infty}{\int_{{-}L}^{L}(\ldots)\frac{\mathrm{d}l}{B}} \right/ {\int_{{-}L}^{L}\frac{\mathrm{d}l}{B}} = \lim_{L\to \infty} \frac{1}{V} {\int_{{-}L}^{L}(\ldots)\frac{ \mathrm{d}l}{B}}, \end{equation}

summing the equations for both species, and taking the real part. Applying the flux-tube average has no impact on the bounce-averaged electron equation, but we leverage the identity

(4.2)\begin{equation} \left\langle \int \overline{f(\boldsymbol{v},l)} \, \mathrm{d}^3 v \right\rangle = \left\langle \int {f}(\boldsymbol{v},l)\, \mathrm{d}^3 v \right\rangle \end{equation}

to allow us to combine the electron and ion contributions to the Helmholtz free energy. After using quasineutrality, we arrive at the Helmholtz free energy, which looks similar to that of Parts 1 and 2,

(4.3)\begin{equation} {H}(\boldsymbol{k}, t) = \sum_{a} \left\langle T_a\int \frac{|g_{a,\boldsymbol{k}}|^2}{F_{a0}}\mathrm{d}^3 v - \frac{e_a^2 n}{T_a}| \delta \phi_{\boldsymbol{k}}|^2 \right\rangle. \end{equation}

The Helmholtz free energy balance reads, as before,

(4.4)\begin{align} \sum_{\boldsymbol{k}} \frac{ \mathrm{d}}{ \mathrm{d} t} H(\boldsymbol{k}, t) & = 2 \operatorname{Re} \sum_{a, \boldsymbol{k}} \left\langle e_a \int {\rm i}\omega_{*a}^T \delta \phi_{\boldsymbol{k}} {\rm J}_{0a} g_{a,\boldsymbol{k}}^* \, \mathrm{d}^3v \right\rangle \nonumber\\ & = 2 \sum_{\boldsymbol{k}} \left(D_i(\boldsymbol{k},t) + D_e^{\mathrm{tr}}(\boldsymbol{k},t)\right); \end{align}

the only difference here is that we have considered the limit ${\rm J}_{0e} \approx 1$ and $g_e = g_e(E_e,\mu _e, \psi, \alpha )$ with $g_e = 0$ for $\lambda < 1/B_{\mathrm {max}}$. We have denoted the electron contribution with a ‘tr’ superscript to signify that only the trapped region of velocity space contributes. The total free energy drive for a given $\boldsymbol {k}$ is then given by $D(\boldsymbol {k}, t) = D_i(\boldsymbol {k},t) + D_e^{\mathrm {tr}}(\boldsymbol {k},t)$, which determines the rate at which the system can grow given the presence of gradients in the plasma.

4.2 Electrostatic energy

To construct the generalised free energy, as demonstrated for a single kinetic species in Part 3, we first need to derive the electrostatic energy balance of the system with bounce-averaged electrons. As shown in Appendix A of Part 3, the electrostatic energy balance equation of a multispecies system may be constructed by applying the operator

(4.5)\begin{equation} \operatorname{Re}\sum_{\boldsymbol{k}} \left\langle \int e_a \delta \phi_{\boldsymbol{k}}^* {\rm J}_{0a}\,(\ldots) \,\mathrm{d^3}v\right\rangle \end{equation}

to the gyrokinetic equations for each species and summing over species. When this operator is applied to the electron equation, the nonlinear term is annihilated. This can be shown by noticing that its contribution,

(4.6)\begin{equation} \operatorname{Re} \sum_{\boldsymbol{k},\boldsymbol{k}'} \left \langle \int \left(k_\psi k_\alpha' -k_\alpha k_\psi' \right)\delta \phi_{\boldsymbol{k}}^*\overline{\delta \phi}_{\boldsymbol{k}'} \, g_{e,\boldsymbol{k}-\boldsymbol{k}'} \, \mathrm{d}^3 v\right\rangle, \end{equation}

when expressed in pitch-angle coordinates (3.8) and using the identity (Helander et al. Reference Helander, Proll and Plunk2013)

(4.7)\begin{equation} \int_{-\infty}^{+\infty} \, \mathrm{d} l \int_{1/B_{\mathrm{max}}}^{1/B(l)} \, \mathrm{d} \lambda \to \int_{1/B_{\mathrm{max}}}^{1/B_{\mathrm{min}}} \, \mathrm{d} \lambda \sum_j \int_{l_{1,j}}^{l_{2,j}} \, \mathrm{d}l , \end{equation}

where the subscript $j$ denotes the different magnetic trapping wells along the flux tube, becomes

(4.8)\begin{equation} \operatorname{Re} \sum_{\boldsymbol{k},\boldsymbol{k}', \sigma, j} \frac{\rm \pi}{V} \int_{1/B_{\mathrm{max}}}^{1/B_{\mathrm{min}}} \mathrm{d} \lambda \,\tau_{B,j} \int_{0}^{\infty} \mathrm{d} v \, v^2 \left(k_\psi k_\alpha' - k_\alpha k_\psi' \right)\overline{\delta\phi}_{\boldsymbol{k}}^*\overline{\delta \phi}_{\boldsymbol{k}'} \,g_{e,\boldsymbol{k}-\boldsymbol{k}'}. \end{equation}

Since $\delta \phi ^*_{\boldsymbol {k}} = \delta \phi _{-\boldsymbol {k}}$ and $g_{e,\boldsymbol {k}}^* = g_{e,-\boldsymbol {k}}$, upon taking the real part, this contribution is proportional to

(4.9)\begin{equation} \left(k_\psi k_\alpha' - k_\alpha k_\psi' \right)\left(\overline{\delta\phi}_{-\boldsymbol{k}}\overline{\delta \phi}_{\boldsymbol{k'}} \,g_{e,\boldsymbol{k}-\boldsymbol{k}'} + \overline{\delta\phi}_{\boldsymbol{k}}\overline{\delta \phi}_{-\boldsymbol{k'}} \, g_{e,-\boldsymbol{k}+\boldsymbol{k}'}\right) \end{equation}

which changes sign if $\boldsymbol {k}$ and $\boldsymbol {k}'$ are swapped, such that this term vanishes upon summation over $\boldsymbol {k}$ and $\boldsymbol {k}'$.

What remains of the electron equation when this operator is applied is

(4.10)\begin{align} & \operatorname{Re}\sum_{\boldsymbol{k}} \left\langle \int -e \delta \phi_{\boldsymbol{k}}^*\left(\frac{\partial g_{e,\boldsymbol{k}}}{\partial t} + {\rm i} \overline{\omega}_{{\rm de}} g_{e,\boldsymbol{k}}\right) \mathrm{d^3}v\right\rangle\nonumber\\ & \quad= \operatorname{Re}\sum_{\boldsymbol{k}}\left\langle \frac{e^2}{T_e}\int F_{e0} \delta\phi^*_{\boldsymbol{k}}\left(\frac{\partial}{\partial t} + i \omega_{*e}^T \right)\overline{\delta \phi}_{\boldsymbol{k}} \,\mathrm{d^3}v \right\rangle. \end{align}

We may also use the identity (4.7) to rewrite the right-hand side of this expression. Upon taking the real part of the right-hand side, the term proportional to $\omega _{*e}^T$ vanishes leaving

(4.11)\begin{equation} \operatorname{Re}\sum_{\boldsymbol{k}} \left\langle \int -e \delta \phi_{\boldsymbol{k}}^*\left(\frac{\partial g_{e,\boldsymbol{k}}}{\partial t} + {\rm i} \overline{\omega}_{{\rm de}} g_{e, \boldsymbol{k}}\right) \mathrm{d^3}v\right\rangle = \frac{1}{4V}\frac{\mathrm{d}}{\mathrm{d} t} \sum_{\boldsymbol{k}} \frac{e^2 n}{T_e}\sum_j \int_{1/B_{\mathrm{min}}}^{1/B_{\mathrm{max}}} \tau_{B,j} |\overline{\delta \phi}_{\boldsymbol{k},j} |^2\, \mathrm{d} \lambda. \end{equation}

Applying the operator (4.5) to the ion gyrokinetic equation proceeds similarly, albeit without bounce averages, yielding

(4.12)\begin{equation} \operatorname{Re}\sum_{\boldsymbol{k}} \left\langle \int e \delta \phi_{\boldsymbol{k}}^*{\rm J}_{0i}\left(\frac{\partial g_{i,\boldsymbol{k}}}{\partial t}+ v_\parallel \frac{\partial g_{i,\boldsymbol{k}}}{\partial l} + {\rm i} {\omega}_{{\rm di}} g_{i,\boldsymbol{k}}\right) \mathrm{d^3}v\right\rangle = \frac{\mathrm{d}}{\mathrm{d} t } \sum_{\boldsymbol{k}}\frac{e^2 n}{T_i}\left\langle \varGamma_{0i}(b)|{\delta \phi}_{\boldsymbol{k}}|^2 \right\rangle, \end{equation}

where $\varGamma _{0i}(b) = {\rm I}_0(b)\exp (-b)$ and $b = k_\perp ^2\rho _i^2$. If we now sum the contributions from the ions and the electrons, given by (4.12) and (4.11), after using quasineutrality, we arrive at electrostatic energy balance,

(4.13)\begin{equation} \frac{\mathrm{d}}{\mathrm{d} t}\sum_{\boldsymbol{k}} E(\boldsymbol{k}, t) = 2 \sum_{\boldsymbol{k}}\left({\rm K}_i(\boldsymbol{k},t) + {\rm K}_e^{\mathrm{tr}}(\boldsymbol{k},t) \right). \end{equation}

Here, $E(\boldsymbol {k}, t)$ is electrostatic energy contained in the potential fluctuations in this system given by

(4.14)\begin{equation} E(\boldsymbol{k},t) = \frac{e^2 n}{T_i}\left[\left\langle \left( 1 + \tau - \varGamma_{0i}(b)\right)|\delta \phi_{\boldsymbol{k}}|^2 \right\rangle- \frac{\tau}{2V}\sum_j \int_{1/B_{\mathrm{min}}}^{1/B_{\mathrm{max}}} \tau_{B,j} |\overline{\delta \phi}_{j,\boldsymbol{k}} |^2\, \mathrm{d} \lambda \right ], \end{equation}

where $\tau = T_i/T_e$. The drive terms ${\rm K}_i(\boldsymbol {k}, t)$ and ${\rm K}_e^{\mathrm {tr}}(\boldsymbol {k},t)$ are the work performed by the ions and trapped electrons on the electric field, respectively, and are given by

(4.15)\begin{equation} {\rm K}_i(\boldsymbol{k},t) ={-}\operatorname{Re} \left\langle e \int \delta \phi_{\boldsymbol{k}}^*{\rm J}_{0i}\left(v_\parallel \frac{\partial g_{i,\boldsymbol{k}}}{\partial l} + {\rm i} {\omega}_{{\rm di}} g_{i,\boldsymbol{k}}\right) \mathrm{d^3}v\right\rangle, \end{equation}

and

(4.16)\begin{equation} {\rm K}_e^{\mathrm{tr}}(\boldsymbol{k},t) = \operatorname{Re} \left\langle e \int {\rm i}\overline{\omega}_{{\rm de}} \delta \phi_{\boldsymbol{k}}^* g_{e,\boldsymbol{k}} \,\mathrm{d^3}v\right\rangle. \end{equation}

The total electrostatic energy drive is then $K(\boldsymbol {k}, t) = {\rm K}_i(\boldsymbol {k},t) + {\rm K}_e^{\mathrm {tr}}(\boldsymbol {k},t)$. The drive, $K$, determines the rate at which the electrostatic energy of the system can grow, given the energy transfer possible by wave–particle interactions. This form of electrostatic energy balance is similar to forms derived by Helander et al. (Reference Helander, Proll and Plunk2013) and Plunk et al. (Reference Plunk, Connor and Helander2017), where they were concerned with normal modes.

4.3 Generalised free energy

Now that we have found forms for both the electrostatic energy and the Helmholtz free energy, we can construct the generalised free energy, as shown in Part 3, by considering the linear combination

(4.17)\begin{equation} \tilde{H} = H - \Delta E, \end{equation}

where $\varDelta$ is a real constant. Here $\tilde {H}$ can be guaranteed to be a positive-definite energy norm if $\varDelta$ is chosen correctly; for example, any negative real choice of $\varDelta$ gives a positive definite energy measure. However, positive-definiteness will not be fulfilled beyond a certain magnitude of positive $\varDelta$. Note that $\tilde {H}$ does not define a single energetic norm, but rather an infinite set of norms for each allowable value of $\varDelta$, all of which are invariant under nonlinear interactions between different wavenumbers if $\varDelta$ is chosen to be independent of $k$ (Plunk & Helander Reference Plunk and Helander2023). The Helmholtz free energy is a special case in this set of energy norms for which $\varDelta = 0$.

The energy balance equation for $\tilde {H}$ is

(4.18)\begin{equation} \frac{\mathrm{d} }{\mathrm{d} t } \sum_{\boldsymbol{k}} \tilde{H}(\boldsymbol{k},t) = 2 \sum_{\boldsymbol{k}}\left(D(\boldsymbol{k},t)- \Delta K(\boldsymbol{k},t) \right). \end{equation}

This energy balance is now ‘aware’ of both the source of free energy growth (the gradients), and the mechanism by which it is extracted (wave–particle interactions). As defined in Part 3, the instantaneous growth rate of the generalised free energy is

(4.19)\begin{equation} \varLambda = \left(D - \Delta K \right)/\tilde{H}. \end{equation}

Note that because $\tilde {H}$ defines a positive definite square norm, the linear, normal mode growth $\gamma _L$ is bounded by the maximum possible value of $\varLambda$,

(4.20)\begin{equation} \gamma_L \leqslant \max_{g_e, g_i} \varLambda. \end{equation}

5 Optimal modes

We now seek the distribution functions $g_i$ and $g_e$ which extremise the value of $\varLambda$ – the optimal modes. These optimal modes are found by varying (4.19) with respect to the distribution function for each species and seeking an extremum ($\delta \varLambda = 0$). This optimisation procedure yields two coupled variational problems, one for each species, given by

(5.1)\begin{equation} \varLambda \frac{\delta \tilde{H}}{\delta g_a} = \frac{\delta D}{\delta g_a} -\varDelta\frac{\delta K}{\delta g_a}. \end{equation}

However, when performing the variation, we must take into consideration that the electron distribution function is independent of the field-line-following coordinate $l$. To enforce this, when computing the variation of a functional $F[g_e]$ as

(5.2)\begin{equation} \frac{\delta F}{\delta g_e} = \frac{{\rm d}}{{\rm d}\varepsilon}\bigg\rvert_{\varepsilon= 0}F[g_e + \varepsilon h_e], \end{equation}

the perturbation $h_e$ is taken to be independent of $l$. To make the impact of $l$-independence clearer, we define the inner product on the space of trapped-electron distribution functions as

(5.3)\begin{equation} \left(g_e, g_e\right)_{{\rm tr}} = \frac{2{\rm \pi}}{V}\sum_{j}\int_{1/B_{\mathrm{max}}}^{1/B_{\mathrm{min}}} \, \mathrm{d} \lambda \int_{0}^{\infty} \mathrm{d} v \,v^2 \tau_{B,j} T_e \frac{|g_e|^2}{F_{e0}}, \end{equation}

which can be derived from the inner product for ion distributions (which we also use for the ion distributions here) defined in Part 3,

(5.4)\begin{equation} (g_i, g_i) = \left\langle T_i \int \, \frac{|g_i|^2}{F_{i0}} \, \mathrm{d}^3v \right\rangle \end{equation}

by expressing the velocity space integration in pitch-angle coordinates, under the assumptions of $l$-independence, and $g_e$ being zero in the passing region of velocity space.

Evaluating the two variational problems (5.1) and expressing the result in terms of the inner products defined above, for the respective species over which the variation is performed, yields two coupled eigenvalue problems for the optimal modes and the growth rate $\varLambda$,

(5.5)\begin{equation} \varLambda \sum_b \tilde{\mathcal{H}}_{{\rm ab}} g_b = \sum_b \left( \mathcal{D}_{{\rm ab}}g_b - \Delta \mathcal{K}_{{\rm ab}}g_b\right), \end{equation}

where the operators $\tilde {\mathcal {H}}_{{\rm ab}}$, ${\mathcal {D}}_{{\rm ab}}$ and $\mathcal {K}_{{\rm ab}}$ are given in Appendix A. Here, and from now on, we will omit the subscript ‘$\boldsymbol {k}$’ for simplicity. These operators are expressed in terms of a finite set of velocity-space moments of the distributions,

(5.6)\begin{gather} \kappa_{1a} = \frac{1}{n}\int g_a {\rm J}_{0a} \, \mathrm{d}^3 v, \end{gather}

(5.7)\begin{gather}\kappa_{2a} = \frac{1}{n}\int \left( \frac{v^2}{v_{Ta}^2}\right) g_a {\rm J}_{0a} \, \mathrm{d}^3 v, \end{gather}

(5.8)\begin{gather}\kappa_{3e} = \frac{1}{n}\int\overline{f(l)(1 - \lambda B /2)}\left( \frac{v^2}{v_{\rm Te}^2}\right) g_e \, \mathrm{d}^3 v, \end{gather}

(5.9)\begin{gather}\kappa_{3i} = \frac{1}{n}\int\left( \frac{v_\perp^2}{2 v_{{\rm Ti}}^2} + \frac{v_\parallel^2}{v_{{\rm Ti}}^2}\right) g_i {\rm J}_{0i} \, \mathrm{d}^3 v, \end{gather}

(5.10)\begin{gather}\kappa_{4i} = \frac{1}{n}\int \left( \frac{v_\parallel}{v_{{\rm Ti}}}\right) g_i {\rm J}_{0i} \, \mathrm{d}^3 v, \end{gather}

(5.11)\begin{gather}\kappa_{5i} =\frac{1}{n}\int \left( \frac{v_\parallel}{v_{{\rm Ti}}}\right) g_i \frac{\partial {\rm J}_{0i}}{\partial l} \, \mathrm{d}^3 v, \end{gather}

where we have written $\hat {\omega }_{{\rm de}}(l) = \tilde {\omega }_{{\rm de}} f(l)$, $f(l)$ is a dimensionless function which captures the $l$-dependence of $\hat {\omega }_{{\rm de}}$, and $\tilde {\omega }_{{\rm de}}$ is a constant with units of inverse time. Note that, due to the appearance of the bounce average in each of the electron operators in Appendix A, the equation which results from the variation of the electron distribution function gives an optimal $g_e$ which is independent of $l$ as required.

This eigenvalue problem for the optimal mode growth has ‘full’ dependence on the magnetic geometry of the flux tube, i.e. the magnetic field strength, the curvature and the metric coefficients, all as a function of $l$. The system can be closed exactly in terms of the moments $\kappa _{n a}$ by taking moments of the two equations given by (5.5). This results in an $8\times 8$ system of integrodifferential equations which must be solved for the eigenvalue $\varLambda$.

This system of equations can be solved for any choice of $\varDelta$, for which $\tilde {H}$ is positive definite, allowing us to compute the optimal mode growth rate for any element in the set of energy norms defined by $\tilde {H}$. Thus, for a given scenario (e.g. plasma parameters and geometry), we can seek the value of $\varDelta$ which corresponds to the most restrictive norm on the dynamics of the system within this set. This will give the tightest possible bound on linear instability growth in accordance with (4.20).

In practical terms, this means that for each scenario we can optimise over $\varDelta$ to find the tightest upper bound on instability growth in the system. Here, this optimisation will be carried out numerically. As explained in Part 3, this tightest upper bound is guaranteed to be less than (or equal to) that of the Helmholtz free energy at $\varDelta = 0$ and the magnitude of the optimal $\varDelta$ can be thought of as the relative importance of wave–particle effects in restricting the optimal growth for a given scenario.

Note that, due to the complexity of the system in general magnetic geometry, we are unable to derive generally the limiting value of $\varDelta$ beyond which positive definiteness is no longer guaranteed. However, this does not pose a problem when solving the system numerically as positive definiteness can be checked at each step when optimising over $\varDelta$.

6 Solving the system in the slow-ion-transit limit

To make finding the solution to the eigenvalue problem more straightforward, we consider the limit $L/v_{{\rm Ti}} \gg \tau _D$, where the ion transit time is much longer than the typical instability time scale. In this limit, we can neglect the derivatives along the magnetic field line that arise due to the parallel ion motion. This reduces the dimensionality of the system by removing $\kappa _{4i}$ and $\kappa _{5i}$ from the problem, leaving a $6 \times 6$ system of integroalgebraic equations. Note that this ordering is the typical ‘toroidal’ ITG and TEM ordering often applied in linear, normal-mode analyses (Biglari, Diamond & Rosenbluth Reference Biglari, Diamond and Rosenbluth1989; Helander et al. Reference Helander, Proll and Plunk2013; Plunk et al. Reference Plunk, Connor and Helander2017). The system of moment equations in this limit is given in Appendix B.

6.1 A square magnetic well

We cannot make much progress towards an analytical solution of this system for a general magnetic geometry. Thus, to gain some insight into the behaviour of the system, we consider a ‘square’ magnetic field of the form

(6.1)\begin{equation} B(l) = \begin{cases} B_{\mathrm{min}}, & -L/2 < l < L/2, \\ B_{\mathrm{max}}, & |l| \geqslant L/2, \end{cases} \end{equation}

with $\omega _{{\rm da}}(l) = \mathrm {const.}$ and $k_\perp = \mathrm {const}$. In this simplified geometry, the $\kappa _{{\rm na}}$ eigenfunctions are constant inside the magnetic well between $-L/2 < l < L/2$ and the $\kappa _{ne}$ moments must vanish for $|l|\geqslant L/2$. In this case, there are two possibilities, either the ion moments must also vanish in this region, or the electron moments are zero everywhere, corresponding to the system with adiabatic electrons. Both of these solutions are realisable, with the adiabatic electron solution being described in Part 3. Here, we concern ourselves with the eigenfunctions that vanish at $B_\mathrm {max}$, whilst keeping in mind that the adiabatic electron solution is also present in the system outside the well.

The solution to the eigenvalue problem in this square magnetic geometry is given in appendix C and is denoted as $\varLambda _{\mathrm {SW}}$. We find that the value of $\varDelta$ beyond which the system is not guaranteed to be Hermitian (i.e. implying a violation of positive-definiteness) is

(6.2)\begin{equation} \varDelta < \frac{1+ \tau}{\varGamma_0 + \tau \varepsilon}, \end{equation}

where $\varepsilon = \sqrt {1 - B_{\textrm {min}}/B_{\textrm {max}}}$. While we will mostly employ this solution as a benchmark for comparison with numerical solutions, it can provide some insight into the behaviour of the optimal mode system in this limit. Firstly, in the Helmholtz limit with $\varDelta =0$, we see that $\varLambda _{\mathrm {SW}}$ reduces to an expression of the same form as the fully kinetic electron Helmholtz bound derived in Appendix D of Part 2, with the only difference being the lack of finite-electron-Larmor-radius effects in the solution here, and the reduction of the size of the electron response by the trapped-particle fraction $\varepsilon$. Later, we will see that this reduction of the population of electrons which contributes to instability growth has a profound impact on the qualitative behaviour of the bound. Moreover, in the limit of $\varepsilon \to 0$, we recover the adiabatic electron Helmholtz bound of Part 1. Thus, the system studied here behaves like an intermediate limit between the fully kinetic and adiabatic electron limits.

For the generalised free energy case, with $\varDelta \neq 0$, we also note that, in the limit of $\varepsilon \to 0$, $\varLambda _{\textrm {SW}}$ reduces to the adiabatic electron result of Part 3 in the toroidal ITG limit.

7 Helmholtz free energy ($\varDelta = 0$) bounds

For the case of $\varDelta = 0$, the generalised free energy reduces to the Helmholtz free energy, and the $\kappa _{3a}$ moments do not contribute to the system of equations. This reduces the dimensionality of the system to $4\times 4$. In this limit, the optimal modes do not depend on magnetic curvature, but retain dependency on $B(l)$ and the variation of $k_\perp$ with $l$. We will refer to the optimal mode growth rate of this system with $\varDelta = 0$ as the Helmholtz bound. To get a sense of how the optimal modes behave in response to the presence of different plasma gradients, we consider three cases: the case of a pure ion-temperature gradient in the plasma ($\boldsymbol {\nabla } n = 0$, $\boldsymbol {\nabla } T_e = 0$) that we refer to as the ‘pure ITG case’; a pure density gradient ($\boldsymbol {\nabla } T_i = 0$, $\boldsymbol {\nabla } T_e = 0$) case that we refer to as the ‘pure TEM case’; and a case with only an electron temperature gradient ($\boldsymbol {\nabla } n = 0$, $\boldsymbol {\nabla } T_i = 0$ which we refer to as the ‘pure ETG-TEM case’. While we have only considered single gradients here for simplicity, we note that the optimal modes can equally be applied to mixed gradient cases.

In figure 1, we show the Helmholtz bound of the system as a function of $b_i = k_\perp ^2\rho _i^2$. Here, for the numerical solution, we consider a magnetic field of the form $B(l) = B_0 - B_1 \cos (l)$ with $B_0 = 1$ and $B_1= 0.2$, and use $k_\perp = \mathrm {const.}$ for simplicity. In a sinusoidal field such as this, the magnetic trapping wells are identical and are separated by identical maxima (like in an omnigenous field), thus we need only solve the system in a single well with $-{\rm \pi} \leqslant l \leqslant {\rm \pi}$. Note that numerically solving the integroalgebraic system results in a spectrum of eigenvalues (an example of which can be seen in figure 2), which are all real due to the Hermiticy of the system, and come in pairs of positive and negative values of equal magnitude due to the time reversibility of the system (Plunk & Helander Reference Plunk and Helander2022). In figure 1 we only show the fastest-growing optimal mode growth rate (the largest positive eigenvalue) for clarity.

Figure 1.

Optimal mode growth rates $\varLambda$ of the Helmholtz free energy with $\tau = 1$ as a function of $k_\perp ^2\rho _i^2$ for the various electron models. Panel (a) shows the pure ITG case, (b) shows the pure TEM case and (c) shows the pure ETG-TEM case.

Figure 2.

Spectrum of numerical solutions to the kinetic eigenvalue problem (5.5) for a sinusoidal magnetic field strength with $\tau = 1$ for the pure-ITG case alongside the adiabatic electron optimal growth rate. Shown are the six largest eigenvalues found numerically.

We also show the analytical Helmholtz bound of the system in a square magnetic trapping well of the same depth, $\varLambda _{\mathrm {SW}}$, given by (C18) with $\varDelta = 0$. We compare these Helmholtz bounds with the bounce-averaged trapped electron response to the Helmholtz bounds with adiabatic electrons, given by (6.20) in Helander & Plunk (Reference Helander and Plunk2022), and the Helmholtz bound with fully kinetic electrons, taken from Appendix C of Plunk & Helander (Reference Plunk and Helander2022). For all gradients considered, we find that the numerical solution in the sinusoidal magnetic field closely follows the analytical solution for the square magnetic trapping well.

For the pure ITG case (figure 1a), we find that the Helmholtz bound of the system with the bounce-averaged trapped electron response lies between the bounds for the adiabatic and fully kinetic electron systems. Moreover, in figure 2 we see that this is true of the entire spectrum of solutions, which appear to be bounded from below in magnitude by the adiabatic electron bound. We note that, for all gradients, the optimal mode growth rate tends to zero as $k_\perp \to 0$ in the bounce-averaged system. This signifies the removal of the ‘magnetohydrodynamic-like’ instability mechanism present in the fully kinetic electron Helmholtz bound, which causes the bound to approach a finite value in this limit. Mathematically, the cause of this difference is illustrated by the square well solution $\varLambda _{\mathrm {SW}}$ where as $k_\perp \to 0$, the numerator tends to zero and the denominator is proportional to $(1 - \varepsilon )$, which is finite unless the trapped particle fraction is unity. Thus, we have $\varLambda _{\mathrm {SW}} \to$ 0 as $k_\perp \to 0$. This is in contrast to the fully kinetic electron bound of Part 2, whose denominator goes to zero in the limit $k_\perp \to 0$ in such a way that the bound remains finite. Physically, the Helmholtz bound with fully kinetic electrons overestimates the degree of instability because it contains a contribution from non-adiabatic passing electrons which, in toroidal geometry, does not contribute on time scales longer than the electron transit time.

The bounce-averaged electron system also yields a Helmholtz bound for the pure TEM case and the pure ETG-TEM case; see figures 1(b) and 1(c). In this bounce-averaged electron system, this provides a bound on the growth of $\boldsymbol {\nabla } n$ and $\boldsymbol {\nabla } T_e$-driven TEMs (this does not extend to ‘true’ ETGs, which exist at larger values of $k_\perp$ thus violating our assumption of $k_\perp \rho _e \ll 1$ and require a non-adiabatic passing electron response).

In figure 3, an example of two of the fastest-growing eigenmodes is shown for the pure TEM case. The eigenfunctions are similar (in $l$-variation) to the expectation from the linear theory, with the fastest growing mode having a maximum where the largest fraction of trapped electrons reside at the minimum of the magnetic field strength. The fastest-growing eigenmode has a significant contribution from the $\kappa _{{\rm ne}}$ moments, indicating their importance for extracting free energy at these parameters.

Figure 3.

Eigenfunctions of the kinetic eigenvalue problem (5.5), obtained numerically, for $\varDelta = 0$ (Helmholtz limit of the generalised free energy) with $\tau = 1$, $b_i = 1.5$, in the pure TEM case. Here, the absolute values of the $\kappa _{{\rm na}}$ moments are shown as a function of $l$, alongside the magnetic field strength $B(l)$, for the eigenfunction with the largest eigenvalue (a) and the second largest (b).

8 Generalised free energy bounds

We now turn our attention to the general problem in which $\varDelta$ is a free parameter. For this case where $\varDelta \neq 0$, the species moments $\kappa _{3a}$, are retained. Thus, the optimal modes of this system depend on the magnetic curvature drift, $\omega _{{\rm da}}$, for each species. We seek values of $\varDelta$ which minimise the maximum growth rate $\varLambda$ of the generalised free energy, giving the lowest possible upper bound. The largest resulting optimal growth rate, at this minimising value of $\varDelta$, gives an upper bound on instabilities, in the considered limit, that depends on the magnetic curvature.

Because of the complexity of the general system, we perform the optimisation over $\varDelta$ numerically and exclude values for which positive definiteness is violated. However, we have found that the optimisation space of $\varLambda ( \varDelta )$ is very smooth and that a global minimum is easily found.

As a first test for the numerical solution, to facilitate comparison with the analytical solution for square magnetic trapping well with constant curvature, we solve the system numerically for the case of a magnetic field of the form $B(l) = B_0 - B_1 \cos (l)$ with $B_0 = 1$ and $B_1= 0.1$ with constant curvature, $\hat {\omega }_{{\rm da}}(l) = \mathrm {const}$. Following the analysis of the generalised free energy with adiabatic electrons of Part 3, we consider the instability drive parameter $\kappa _d = \omega _{*i}/\hat {\omega }_{{\rm di}}$ or ($\kappa _d= \omega _{*i}\eta _i/\hat {\omega }_{{\rm di}}$ for the pure ITG case and $\kappa _d= |\omega _{*e}\eta _e|/\hat {\omega }_{{\rm di}}$ for the pure ETG-TEM case). We thus have $\hat {\omega }_{{\rm da}}/\omega _{*i} = \mathrm {sign}(e_a)\kappa _d^{-1}T_a/{T_i}$ for each species. Under this normalisation, if $\kappa _d > 0$, then $\hat {\omega }_{{\rm da}}\omega _{*a} > 0$ for both species, corresponding to ‘bad’ curvature where linear instability is expected. The case of $\kappa _d < 0$ corresponds to ‘good’ curvature, where the system is expected to be linearly stable.

In figure 4, we show the optimal mode growth for the pure TEM case, the pure ETG-TEM case, and for the pure ITG case, as the curvature-drive parameter $\kappa _d$ is varied. We see that the square and the cosine trapping wells follow similar trends with $\kappa _d$, attaining a maximum for a finite positive value of $\kappa _d$, where the optimal bound is obtained for $\varDelta = 0$ and therefore coincides with the Helmholtz bound. This feature of the generalised free energy bound overlapping with the Helmholtz bound at a particular value of the curvature drive was also found for the adiabatic electron analysis in Part 3. The optimal bound decreases monotonically away from this maximum towards the ‘strongly driven’ limit of large positive $\kappa _d$ and in the direction of small $\kappa _d$ in the ‘resonant limit’. The bound also decreases for $\kappa _d \leqslant 0$ where the curvature is favourable (‘good’).

Figure 4.

The optimal mode growth rate of generalised free energy, at the minimising value of $\varDelta$, versus the drive parameter $\kappa _d$ at $\tau = 1$. The analytical square well with constant curvature solution is shown alongside the numerical solutions of the eigenvalue problem in a cosine magnetic well with constant curvature. These are computed for a density-gradient-driven case (pure TEM case), an ion-temperature-gradient-driven case (pure ITG case) and an electron-temperature-gradient-driven case (pure ETG-TEM case).

In figure 5, we compare the optimal mode growth in the square well, $\varLambda _{\mathrm {SW}}$, with constant curvature to the solution to the linear dispersion relation with bounce-averaged electrons (details given in appendix D) in the same geometry, in the pure ITG case. Once again, we vary the drive parameter $\kappa _d$ to compare how each solution depends on the curvature inside the well. We confirm that the linear result lies below the upper bound, and we find that, much like the behaviour of the upper bound, the presence of bounce-averaged electrons in the linear dispersion relation ‘lifts up’ the linear growth rate of the ITG. In regions where the normal modes are stable, the optimal modes have a finite but reduced growth rate, as was found in Part 3 for the adiabatic electron system.

Figure 5.

Comparison between the generalised upper bound $\varLambda _{\mathrm {SW}}$, and the linear growth rate, for a pure ITG case in a square magnetic well with $B_{\mathrm {min}}/B_\mathrm {max} = 0.1$, $\tau = 1$, and constant curvature (see Appendix D). Here, the drift-kinetic limit has been considered for simplicity (${\rm J}_{0i} \approx 1$).

8.1 Impact of the maximum-$J$ property

A class of magnetic configurations in which the gyrokinetic system with bounce-averaged electrons is particularly interesting is provided by so-called ‘maximum-$J$’ geometries, where the radial profile of the parallel adiabatic invariant,

(8.1)\begin{equation} J(v,\lambda, \psi, \alpha) = m_e v\int_{l_1}^{l_2}\sqrt{1 - \lambda B} \,\mathrm{d} l, \end{equation}

has a maximum at the magnetic axis, and decreases radially (Rosenbluth Reference Rosenbluth1968). In such geometries, the electrons experience good curvature everywhere because

(8.2)\begin{equation} \omega_{*e}\overline{\omega}_{{\rm de}} ={-}\frac{T_e k_\alpha^2}{e^2\tau_B}\frac{\mathrm{d}\ln n}{\mathrm{d}\psi}\frac{\partial J}{\partial \psi} \end{equation}

is negative for all trapped electrons if $\partial J / \partial \psi < 0$ and ${\mathrm {d}\ln n}/{\mathrm {d}\psi } < 0$, as is usually the case in fusion plasmas (Proll et al. Reference Proll, Helander, Connor and Plunk2012; Helander et al. Reference Helander, Proll and Plunk2013). Quasi-isodynamic stellarators are a particular type of optimised stellarator (Helander & Nührenberg Reference Helander and Nührenberg2009; Nührenberg Reference Nührenberg2010) which can be designed to have this property (Sánchez et al. Reference Sánchez, Velasco, Calvo and Mulas2023; Rodríguez, Helander & Goodman Reference Rodríguez, Helander and Goodman2024). In such devices, the local curvature (as seen by the ions) may be unfavourable, but the bounce-averaged curvature for trapped electrons is always favourable. Many results from normal-mode theory, as well as gyrokinetic simulations, have found that the trapped electron population in these geometries have a stabilising influence if a plasma density gradient is present (Proll et al. Reference Proll, Helander, Connor and Plunk2012; Plunk et al. Reference Plunk, Helander, Xanthopoulos and Connor2014; Alcusón et al. Reference Alcusón, Xanthopoulos, Plunk, Helander, Wilms, Turkin, von Stechow and Grulke2020; Proll et al. Reference Proll, Plunk, Faber, Görler, Helander, McKinney, Pueschel, Smith and Xanthopoulos2022). The physical reason for this enhanced stability can be traced back to the fact that relatively little of the thermal energy is ‘available’ for conversion to instabilities (Helander Reference Helander2017; Mackenbach et al. Reference Mackenbach, Proll, Wakelkamp and Helander2023b).

To investigate if this property is present in the optimal mode system, we require a test geometry where the ions may experience bad curvature locally, while the bounce-averaged curvature of the electrons can be independently varied. A simple model magnetic geometry can be seen in figure 6. In this geometry, we choose a curvature of the form $\hat {\omega }_{{\rm da}}(l)/\omega _{*i} = \mathrm {sign}(e_a)\sigma _d \cos (l) T_a/{T_i}$ and $B(l) = B_0 + B_1\cos (l)$, where once again, we take $B_0 = 1.0$ and $B_1 = 0.1$. In this toy geometry, $\sigma _d$ is a parameter which determines the relative magnitude of the gradient length scale and the maximum radius of curvature, similar to the familiar $R/L_n$ or $R/L_{T_a}$ normalisation, where $R$ is the major radius of the torus and $L_{n,T_i}$ is the scale length of the density or temperature gradients, respectively. With this model geometry, by changing the sign of $\sigma _d$ we can choose between a mostly maximum-$J$ geometry for $\sigma _d < 0$, where most trapped electrons, particularly those which are deeply trapped, experience bounce-averaged good curvature, and a mostly minimum-$J$ geometry for $\sigma _d > 0$. Importantly, we can make this change between maximum-$J$ and minimum-$J$ geometry without changing the field-line-average curvature (zero in this case), whilst always having a region of unfavourable curvature along the field line.

Figure 6.

An example of an approximately maximum-$J$ toy geometry, where most trapped electrons, particularly those that are deeply trapped, experience bounce-averaged good curvature. (a) The curvature is shown alongside the magnetic field strength, where negative values indicate bad curvature. (b) The $\lambda$-dependence of the bounce-averaged electron drift is shown.

In figure 7, we show the optimal growth rate $\varLambda$ as $\sigma _d$ is varied, the pure ITG case, pure TEM case and pure ETG-TEM case. The upper bound in the pure TEM case is lower in the maximum-$J$ geometries (negative $\sigma _d$) than in the equivalent minimum-$J$ geometries (positive $\sigma _d$). To explore the physical mechanism at play, we also test the dependence of the optimal mode growth on the curvature drift for each species separately. For $\sigma _d<0$ (maximum-$J$), removing the electron drift frequency (the instability mechanism that remains is analogous to the ITEM of Plunk et al. (Reference Plunk, Connor and Helander2017)) increases growth rates if $|\sigma _d|$ is small, and has little effect for $\sigma _d$ large and negative. Thus, $\overline {\omega }_{{\rm de}}$ has a stabilising effect in the maximum-$J$ geometries, in agreement with the expectation linear theory. Moreover, in maximum-$J$ geometries at large negative values of $\sigma _d$, the instability mechanism is largely reliant on the ion curvature drift, as is expected from the linear ITEM or ubiquitous mode (Coppi & Pegoraro Reference Coppi and Pegoraro1977).

Figure 7.

The optimal mode growth rate of generalised free energy, for the minimising $\varDelta$, versus the drive parameter $\sigma _d$. Negative and positive values of this parameter correspond to the ‘maximum-$J$’ and ‘minimum-$J$’ cases of the model geometry shown in figure 6. (a) The pure TEM case is shown at $k_\perp \rho _i = 1.5$, (b) the pure ITG case in the drift kinetic limit (${\rm J}_{0i} \approx 1$) is shown and (c) shows the pure ETG-TEM case at $k_\perp \rho _i = 1.5$. The impact of species-dependent curvature drifts on each case is explored by zeroing out $\hat {\omega }_{{\rm di}}$ or $\overline {\omega }_{{\rm de}}$ independently.

For $\sigma _d > 0$ (minimum-$J$), the removal of the curvature drift for either species in the pure TEM case serves to move the maximum value of the optimal growth rate towards lower values of $\sigma _d$ and greatly reduces the growth rate for large $\sigma _d$. This suggests that the drifts of both species are important for the optimal extraction of free energy from the gradients.

In the pure ETG-TEM case, we see a reduction of the optimal mode growth for $\sigma _d < 0$. This is in keeping with the expectation from the dispersion relation of Plunk et al. (Reference Plunk, Connor and Helander2017), which shows that the maximum-$J$ property, associated with the relative direction of the bounce-averaged drift to the electron diamagnetic frequency, also exerts a stabilising influence in the strongly driven limit when $\eta _e$ is finite. The optimal mode analysis here suggests that this benefit extends into the resonant limit (small $\sigma _d$).

We observe that maximum-$J$ configurations also exhibit lower optimal growth in the ITG scenario ($\boldsymbol {\nabla } n = 0$, $\boldsymbol {\nabla } T_e = 0$). This is despite the absence of any electron free energy source, something which has not been predicted previously from linear normal-mode theory. It is found that beyond a negative value of $\sigma _d$ of sufficient magnitude, the ITG upper bound with bounce-averaged electrons coincides with that of the adiabatic electron upper bound of Part 3 computed in this geometry (because the bound with adiabatic electrons is local in this limit, here we evaluate the bound at each point in $l$ in the toy geometry and plot the largest value found). This is also reflected in the eigenfunction (not shown), which tends towards a delta function at the location of maximum magnetic field strength, where the trapped electron contribution is zero, suggesting that the adiabatic ITG solution to (5.5) becomes the solution with the largest $\varLambda$ for these values of $\sigma _d$. As shown in figure 7, this mechanism of stabilisation in maximum-$J$ geometry is present even if $\overline{\omega }_{\rm de} = 0$. Thus, this stabilising mechanism in maximum-$J$ geometries is distinct from the known mechanism that follows from the difference in sign between $\overline {\omega }_{\rm de}$ and $\omega _{*e}^T$ (Proll et al. Reference Proll, Plunk, Faber, Görler, Helander, McKinney, Pueschel, Smith and Xanthopoulos2022).

8.1.1 Interpreting with linear theory

To interpret the observed reduction of the ITG-driven optimal mode growth in maximum-$J$ geometries without a source of electron free energy, we turn to expressions from linear normal-mode theory (i.e solutions to the gyokinetic equation which evolve in time like $\exp ({-{\rm i} \omega t})$, where $\omega$ is the complex frequency). The comparable linear dispersion relation is given by (3.5) in Plunk et al. (Reference Plunk, Connor and Helander2017), for the case of a pure ITG, which is valid if the mode frequency satisfies $\omega /(\omega _{*i}\eta _i) \sim (\omega _{{\rm di}}/\omega _{*i}\eta _i)^{1/2} \ll 1$ (i.e. in the regime of large $\kappa _d$). Taking the drift-kinetic limit (${\rm J}_{0i} \approx 1)$, the normal mode frequency is given by

(8.3)\begin{align} \omega & = \left.\pm \sqrt{-\omega_{*i}\eta_i \int_{-\infty}^{+\infty} \hat{\omega}_{{\rm di}}(l) | \delta \phi|^2 \frac{\mathrm{d}l}{B}} \right/\nonumber\\ & \quad {\left(\tau \int_{-\infty}^{+\infty}| \delta \phi|^2 \frac{\mathrm{d}l}{B} - \frac{\tau}{2}\sum_j \int_{1/B_{\mathrm{min}}}^{1/B_{\mathrm{max}}} \tau_{B,j} |\overline{\delta \phi}_{j} |^2\, \mathrm{d} \lambda\right)^{{1}/{2}}}, \end{align}

where instability growth is given by a positive imaginary part of $\omega$ which arises due to regions of bad curvature along the magnetic field ($\omega _{*i}\eta _i \hat {\omega }_{{\rm di}} > 0$). We see that when the ITG is unstable, the trapped electron contribution in the denominator acts to further destabilise the mode by making the denominator (which is always positive (Helander et al. Reference Helander, Proll and Plunk2013)) smaller. This destabilising effect is due to the depletion of the Boltzmann response by the population of non-adiabatic trapped electrons.

The effect of the maximum-$J$ property on this dispersion relation is most easily illustrated when the case of the ‘square well’ magnetic field with constant curvature inside the well is considered. Inside the well, the ITG frequency becomes $\omega = \pm \sqrt {-\omega _{*i}\eta _i \hat {\omega }_{{\rm di}}/ (\tau (1- \varepsilon ))}$, where $\varepsilon$ is the trapped particle fraction, and outside the well it is given by $\omega = \pm \sqrt {-\omega _{*i}\eta _i \hat {\omega }_{{\rm di}}/\tau }$. If we consider a series of these square magnetic trapping wells with minima of $B = B_{\mathrm {min}}$ separated by regions of $B = B_{\mathrm {max}}$, with a magnetic curvature that alternates sign (from ‘good’ curvature to ‘bad’ curvature) between regions of $B = B_{\mathrm {max}}$ and $B = B_{\mathrm {min}}$, then there are two possibilities: either the bad curvature is located at the minimum of magnetic field strength, corresponding to a minimum-$J$ geometry; or the bad curvature is located outside the well, corresponding to a maximum-$J$ geometry.

In the minimum-$J$ case, the ITG is unstable inside the well, with a growth rate given by $\omega = \pm \sqrt {-\omega _{*i}\eta _i \hat {\omega }_{{\rm di}}/ (\tau (1- \varepsilon ))}$. In this scenario, the non-adiabatic population of trapped electrons can raise the ITG growth rate by depleting the Boltzmann response. On the other hand, in the maximum-$J$ case, the ITG is stable inside the well because the curvature is favourable there ($\hat {\omega }_{{\rm di}} \omega _{*i}\eta _i> 0$). Instead, the unstable ITG is outside of the well at $B= B_{\mathrm {max}}$ and the growth rate is $\omega = \pm \sqrt {-\omega _{*i}\eta _i \hat {\omega }_{{\rm di}}/\tau }$, corresponding to the ITG growth rate if the electrons are treated adiabatically.

This behaviour can also be inferred from the ITG growth rate in general geometry (8.3), in which it can be seen that for instability to occur, the average of $\omega _{*i}\eta _i\hat {\omega }_{{\rm di}}| \delta \phi |^2$ along the magnetic field line must be positive, such that $\delta \phi$ must have large amplitude in the regions of bad curvature. In maximum-$J$ geometry, the regions of bad curvature are at different locations in $l$ than the minima of magnetic field strength. In these geometries, for the instability criterion to be fulfilled, the non-adiabatic electron response, which is proportional to the average potential seen by the trapped electrons as they bounce in the magnetic well, will be necessarily small. Otherwise, if the electron contribution were large, the potential would have a large amplitude in a region of favourable curvature, thus impacting the instability criterion. As a result, the destabilising impact of trapped electrons on the ITG is diminished in such geometries.

This is not the case in minimum-$J$ geometries, where the regions of bad curvature are located at the minima of $B$ along the field line. In these geometries, the instability criterion and the magnitude of the trapped electron contribution are not in conflict, thus the ITG growth can be increased by the trapped electron population without impacting the instability mechanism.

In the optimal mode analysis, this effect is slightly obscured due to the optimisation over the $\varDelta$ parameter, but the same basic behaviour can be observed in the expression for $E(\boldsymbol {k}, t)$ given by (4.14), which is of the same form as the denominator in (8.3). Thus, the magnitude of $E(\boldsymbol {k},t)$ affects the magnitude of $\varLambda$ via the denominator of (4.19), for non-zero values of $\varDelta$. This stabilising effect of maximum-$J$ can be seen in figure 7, where for the ITG scenario, the upper bound with bounce averaged electrons approaches the adiabatic electron result for large negative values of $\sigma _d$, corresponding to an approximately maximum-$J$ geometry.

To summarise this section, the impact of the trapped electron response on ITG instability in the absence of electron free energy is always destabilising, but the amount of destabilisation is impacted by the location of the most deeply trapped particles relative to the regions of bad curvature. To use a term coined in the literature for zonal-flows (Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005), the location of the trapped particles changes the inertia of the ITG mode, with ITGs possessing lower inertia in minimum-$J$ configurations than in those which are maximum-$J$. This may go some way towards explaining why simulations performed in newly optimised, highly maximum-$J$ stellarators (Goodman et al. Reference Goodman, Xanthopoulos, Plunk, Smith, Nührenberg, Beidler, Henneberg, Roberg-Clark, Drevlak and Helander2024) show a smaller increase in ITG-driven heat fluxes with the addition of kinetic electrons in comparison with configurations which do not satisfy the maximum-$J$ property to the same degree.