1 Introduction
This is the third paper in a series (Helander & Plunk Reference Helander and Plunk2022; Plunk & Helander Reference Plunk and Helander2022), in which we develop a linear and nonlinear stability theory based on gyrokinetic energy balance. The last two papers used Helmholtz free energy, and introduced the concept of optimal mode growth for fully electromagnetic gyrokinetics. The present paper proposes a generalized energetic measure of fluctuations, allowing the inclusion of additional instability mechanisms. We do this first for a simple case, namely the electrostatic limit (low plasma 
$\beta$) with only one kinetic species (ions), with the electrons being treated adiabatically. These simplifications limit the application to ion-temperature-gradient (ITG) driven turbulence, although the central result of the paper is capable of treating completely general details of the magnetic geometry.
Free energy is a useful concept for understanding nonlinear and linear aspects of plasma turbulence. At the level of linear instabilities it is common to speak of a source of free energy that drives modes. Indeed, without a source of free energy, provided by background plasma gradients (density, temperature, flows), there can be no linear instabilities (nor can there be subcritical turbulence Landreman, Plunk & Dorland Reference Landreman, Plunk and Dorland2015; Plunk & Helander Reference Plunk and Helander2022). However, there is usually another ingredient that arises in the detailed analysis of normal linear instabilities, namely the wave–particle resonance. In gyrokinetic theory, this involves parallel motion (along the magnetic field) and magnetic drift, and the resonance is physically linked to the work that the electrostatic field performs on gyrocentre motion. However, the terms needed to capture this do not contribute to free energy balance, and the influence of resonance therefore cannot be accounted for by the optimal modes that we introduced in our previous works.
In this work, we propose a new measure of gyrokinetic fluctuations, a generalization of the concept of free energy, that incorporates the resonance mechanism, and, via the magnetic drift, the full details of the background magnetic geometry. We demonstrate the existence of a class of quadratic measures closely related to the Helmholtz free energy that behave as positive–definite norms for fluctuations in the distribution function. The corresponding energy balance equation is then used to derive a theory of optimal modes that most efficiently extract this energy from its source. The growth rate of these optimal modes provides a rigorous upper bound on the growth rate of linear instabilities, and this bound is shown to be lower than that obtained previously from the Helmholtz free energy. By studying some simple limits, we show that we recover some expected behaviour of both the slab and toroidal branches of the ITG mode.
2 Definitions and gyrokinetic energy balance
We follow closely the conventions of the first paper in this series (Helander & Plunk Reference Helander and Plunk2022), focusing on local gyrokinetic theory in the geometry of a flux tube, whereby fluctuations in the distribution function may be considered small, and periodic in the coordinates perpendicular to the field line. Essential definitions are summarized in what follows, but more detail and background can be found in §§ 2 and 3 of Helander & Plunk (Reference Helander and Plunk2022) (henceforth also called ‘Part 1’).
The ion gyrokinetic equation in the electrostatic limit is written
\begin{equation} \frac{\partial g_{\boldsymbol{k}}}{\partial t} + v_{\|} \frac{\partial g_{\boldsymbol{k}}}{\partial l} + {\rm i} \tilde{\omega}_d g_{\boldsymbol{k}} + \frac{1}{B^2} \sum_{{\boldsymbol{k}}'} {\boldsymbol{B}} \boldsymbol{\cdot} ({\boldsymbol{k}} \times {\boldsymbol{k}}') \overline{\delta\phi}_{{\boldsymbol{k}}'} g_{{\boldsymbol{k}} - {\boldsymbol{k}}'} = \frac{e_i F_{0}}{T_i} \left( \frac{\partial}{\partial t} + {\rm i} \omega_{{\ast}}^T \right) \overline{\delta\phi}_{\boldsymbol{k}}, \end{equation}
where 
$g$ is the gyrocentre-dependent part of the perturbed ion distribution function, i.e. 
$f_i = ( 1 - e_i \delta \phi (\boldsymbol {r})/T_i )F_{i0} + g(\boldsymbol {R}, E_i, \mu _i, t).$ Its phase space variables are the energy 
$E_a = m_a v^2 / 2 + e_a \varPhi (\psi )$ and the magnetic moment 
$\mu _a = m_a v_\perp ^2 / (2B)$, and the perpendicular wavenumber is 
${\boldsymbol {k}} = {\boldsymbol {k}}_\perp = k_\psi \boldsymbol {\nabla } \psi + k_\alpha \boldsymbol {\nabla } \alpha$ with 
$k_\psi$ and 
$k_\alpha$ independent of the arc length 
$l$ along the magnetic field, and 
$\psi$ and 
$\alpha$ defined via 
${\boldsymbol {B}} = B {\boldsymbol {b}} = \boldsymbol {\nabla } \psi \times \boldsymbol {\nabla } \alpha$. We neglect collisions here,Footnote 1 and used the simplified notation 
$g_{\boldsymbol {k}} = g_{i, {\boldsymbol {k}}}$, and 
$\omega _{\ast } = \omega _{\ast i}$, etc. because the adiabatic approximation 
$g_{e, {\boldsymbol {k}}} = 0$ is assumed throughout; note that it is not necessary here to include the customary correction for the zonal component (i.e. for 
$k_\alpha = 0$; see Dorland & Hammett Reference Dorland and Hammett1993) because the growth of this component is zero due to the fact that it has no source of free energy (see below). We will also assume 
$k \rho _i \sim 1$, implying 
$k\rho _e \ll 1$.
As derived for the general case in § 3 of Part 1, the gyrokinetic free energy balance equation obtained in this limit reads
\begin{equation} \frac{{\rm d}}{{\rm d} t}\sum_{\boldsymbol{k}} H = 2 \sum_{\boldsymbol{k}} D,\end{equation}
where the drive term 
$D$ is
\begin{equation} D({\boldsymbol{k}}, t) = {\rm Im} \, e_i \left\langle \int g_{\boldsymbol{k}} \omega_{{\ast}}^T \overline{\delta\phi}^\ast_{\boldsymbol{k}}\, {\rm d}^3v \right\rangle, \end{equation}and the free energy, expressed in terms of the gyrocentre distribution function
\begin{equation} H({\boldsymbol{k}},t) = \left\langle T_i \int \frac{|g_{\boldsymbol{k}}|^2}{F_{i0}} {\rm d}^3v - \sum_a \frac{n_a e_a^2}{T_a} |\delta \phi_{\boldsymbol{k}}|^2 \right\rangle,\end{equation}where the space average is defined as (extensions are discussed in § 3 of Part 1)
\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}The diamagnetic frequencies are
\begin{equation}
\begin{gathered}
\omega_{{\ast} a} = \frac{k_\alpha T_a}{e_a} \frac{{\rm d}
\ln n_a}{{\rm d} \psi},\\ \omega_{{\ast} a}^T =
\omega_{{\ast} a} \left[1 + \eta_a \left( \frac{m_a v^2}{2
T_a} - \frac{3}{2} \right)\right],
\end{gathered}
\end{equation}
and the magnetic drift frequency is
\begin{equation} \tilde{\omega}_d = {\boldsymbol{k}} \boldsymbol{\cdot} \boldsymbol{v}_d, \end{equation}
where the magnetic drift velocity is 
$\boldsymbol {v}_d = \hat {\boldsymbol {b}}\times ((v_{\perp }^2/2)\boldsymbol {\nabla } \ln B + v_\|^2\boldsymbol {\kappa })/\varOmega _i$, 
$\boldsymbol {\kappa } = \hat {\boldsymbol {b}}\boldsymbol {\cdot }\boldsymbol {\nabla }\hat {\boldsymbol {b}}$ and 
$\varOmega _a = e_a B / m_a$ is the gyrofrequency. Noting 
$\hat {\boldsymbol {b}}\times (\boldsymbol {\nabla }\ln B -\boldsymbol {\kappa } - \mu _0/B^2 \boldsymbol {\nabla } p) = 0$, we can assume 
$\hat {\boldsymbol {b}}\times \boldsymbol {\nabla }\ln B \approx \hat {\boldsymbol {b}}\times \boldsymbol {\kappa }$ in the appropriate limit of low plasma 
$\beta$, which allows us to separate the drift frequency into velocity-dependent and space-dependent factors following Plunk et al. (Reference Plunk, Helander, Xanthopoulos and Connor2014):Footnote 2
\begin{equation} \tilde{\omega}_d = \omega_d(l) \left[ \frac{v_{{\perp}}^2}{2v_{{\rm th}}^2} + \frac{v_{{\parallel}}^2}{v_{{\rm th}}^2} \right]. \end{equation}The gyro-averaged electrostatic potential is denoted
\begin{equation} \overline{\delta\phi}_{\boldsymbol{k}} = {\rm J}_0 \left( \frac{k_\perp v_\perp}{\varOmega_i} \right) \delta \phi_{\boldsymbol{k}}, \end{equation}and the quasi-neutrality condition is
\begin{equation} \sum_a \frac{n_a e_a^2}{T_a} \,\delta \phi_{\boldsymbol{k}} = e_i \int g_{\boldsymbol{k}} {\rm J}_{0} {\rm d}^3v, \end{equation}
where 
$\textrm {J}_{n} = \textrm {J}_n(k_\perp v_\perp / \varOmega _i)$ is a Bessel function. Following our previous convention, we define the free energy as twice that which appears in some other publications. Henceforth, we suppress the 
$\boldsymbol {k}$-subscripts.
2.1 Electrostatic energy and positive–definiteness of free energy
It is useful to decompose the free energy into a part associated with a perturbed distribution function and a part associated with fluctuations in the electrostatic field, i.e.
\begin{equation} H = G + E,\end{equation}where
\begin{gather} G ={-}T_i S_i = \left\langle T_i \int \frac{|\delta F|^2}{F_{i0}} {\rm d}^3v \right\rangle \end{gather}
\begin{gather}E = \left\langle \left(\tau + 1 - \varGamma_{0} \right)\frac{n_i e_i^2}{T_i}|\delta \phi|^2 \right\rangle. \end{gather}
Recall the conventional definitions 
$\varGamma _n(b) = \exp (-b)\textrm {I}_n(b)$ and 
$b = k_\perp ^2 \rho _i^2 = k_\perp ^2 T_i / (m_i \varOmega _i^2)$, and 
$\tau = (e T_i)/(e_i T_e)$. Note that 
$\delta F = g - (e_i \overline {\delta \phi }/T_i) F_0$ is the gyro-averaged perturbed distribution function, and these two contributions to 
$H$ can be identified as the gyrokinetic perturbed entropy and the gyrokinetic field energy.
Although the general electromagnetic free energy admits a similar form as (2.11) ((3.11) Part 1), we note that the electrostatic limit is distinguished by the fact that the field contribution 
$E$ is itself a nonlinear invariant of the gyrokinetic system (Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009), by which we mean that it is conserved under the sole action of the nonlinearity, but not by all the linear terms; the same applies for the Helmholtz free energy. The conservation of 
$E$ may be viewed as an additional constraint on the nonlinear dynamics, with consequences e.g., for the cascade and production of large-scale 
$E\times B$ flows (Plunk et al. Reference Plunk, Cowley, Schekochihin and Tatsuno2010).
For what follows, we need the electrostatic energy balance equation. This is obtained by multiplying the ion gyrokinetic equation by 
$e_i \overline {\delta \phi }^\ast$, integrating over velocity, averaging over the parallel coordinate 
$l$ and summing over perpendicular wavenumber 
${\boldsymbol {k}}$, yielding (Helander, Proll & Plunk Reference Helander, Proll and Plunk2013)
\begin{equation} \frac{{\rm d}}{{\rm d}t}\sum_{\boldsymbol{k}} E = 2 \sum_{\boldsymbol{k}} K,\end{equation}
where the drive term 
$K$ is
\begin{equation} K ={-}{\rm Re} \, e_i \left\langle \int \overline{\delta\phi}^\ast\left(v_{\|}\frac{\partial }{\partial l} + {\rm i}\omega_{d}\right)g \,{\rm d}^3v \right\rangle. \end{equation}This is composed to two contributions, one coming from the parallel streaming term, and the other coming from the magnetic drift term. The first contribution has a simple physical interpretation, as the rate of energy exchanged between particles and the parallel electric field (i.e. the volume average of the parallel current multiplied by the parallel electric field), while the second term describes an analogous process in the perpendicular direction associated with the drift motion of gyrocentres.
Equation (2.11) is a physically transparent form that makes it clear that the free energy 
$H$ is a positive–definite norm for the distribution function 
$g$,Footnote 3 i.e.
\begin{equation} H \geqslant 0\quad\text{and}\quad H = 0 \text{ iff } g = 0, \end{equation}
over all of phase space, 
$\ell$ and 
$\boldsymbol {v}$. To see this, note that the quantities 
$G$ and 
$E$ are both positive, i.e. 
$G \geqslant 0$, obviously, and 
$E \geqslant 0$ because 
$\varGamma _0 \leqslant 1$. Therefore, if 
$H = 0$ then both 
$E = 0$ and 
$G = 0$. The first implies 
$\delta \phi = 0$ everywhere, while the second implies 
$\delta F = 0$ over all of phase space; 
$\delta \phi = 0$ and 
$\delta F = 0$ obviously implies 
$g = 0$.
We note that positive–definiteness is a desirable property of an energetic measure that can be useful for setting bounds on the growth rate of fluctuations; if a non-zero fluctuation (
$g\neq 0$) has zero measure 
$M$, while 
$\textrm {d}M/\textrm {d}t \neq 0$, then the rate 
$M^{-1}\,\textrm {d}M/\textrm {d}t$ is unbounded, so that the problem of optimal modes (i.e. determining the form of 
$g$ that maximizes growth) is ill posed. For an example of this, consider 
$M = E$. One may find functions 
$g_\epsilon$ such that 
$|\delta \phi | < \epsilon$ for arbitrarily small 
$\epsilon > 0$, while 
$g_\epsilon \sim 1$ is itself not small. In this case 
$\textrm {d}E/\textrm {d}t \sim \epsilon$ while 
$E \sim \epsilon ^2$ so the rate 
$E^{-1}\,\textrm {d}E/\textrm {d}t$ is divergent in the limit 
$\epsilon \rightarrow 0$.
Although we mainly consider a plasma with a single kinetic ion species and adiabatic electrons, the concepts and the formalism carry over to the more general case of a plasma with an arbitrary number of kinetic species, as shown in Appendix A. An important limitation, however, is that magnetic fluctuations and collisions are neglected.
2.2 Generalized free energy
The positive–definiteness of 
$H$ suggests a family of related quadratic energetic measures that are also positive–definite. In particular it is clear that something of the form
\begin{equation} \skew2\tilde{H} = H - \Delta E,\end{equation}
will be positive–definite, by the same arguments as the previous section, for particular values of the parameter 
$\varDelta$. For instance the choice 
$\varDelta < 1$ allows trivial generalization of the arguments, but we will see that the value can be extended beyond this.
To find a range of permissible values of 
$\varDelta$, and to help simplify subsequent derivations, we will consider a linear transformation on the distribution function, i.e. we define a new distribution function 
$\tilde {g}$ in terms of which the energy 
$\skew2\tilde{H}$ may be expressed using the Euclidean (or 
$L^2$) norm
\begin{equation} \skew2\tilde{H} = ||\tilde{g}||^2 = (\tilde{g}, \tilde{g}),\end{equation}where we have introduced the inner product
\begin{equation} (\tilde{g}_1, \tilde{g}_2) = \left\langle T_i \int \frac{\tilde{g}_1^{{\ast}}\tilde{g}_2^{}}{F_{0}} {\rm d}^3v \right\rangle. \end{equation}
We will refer to (2.18) as a ‘diagonal’ form of the norm, as it does not involve additional linear operations on the distribution function, compared with the form given by (2.11). To find the relationship between 
$\tilde {g}$ and 
$g$, we introduce the ansatz 
$\tilde {g} = g - \nu \textrm {J}_0 F_0 e_i\delta \phi /T_i$, and substitute this into (2.18). By evaluating velocity integrals using quasi-neutrality (2.10) and (2.13), one can find the free parameter 
$\nu$ that yields the form (2.17), that is
\begin{equation} \nu = \frac{1}{\varGamma_0}\left(1+\tau - \sqrt{(1+\tau - \varGamma_0)(1+\tau -\Delta \varGamma_0)}\right), \end{equation}
where we have taken the negative root for convenience. Observe that in order for 
$\nu$ to be real, we must have
\begin{equation} \varDelta \leqslant (1 + \tau)/\varGamma_0.\end{equation}
The parameter 
$\varDelta$ can of course be negative, in which case its magnitude is unbounded. Noting that 
$\varGamma _0$ generally depends on 
$k$, we may also assume the more restrictive 
$\varDelta \leqslant (1 + \tau )$ to ensure that 
$\skew2\tilde{H}$ remains a nonlinear invariant.
We pause to note that the choice 
$\varDelta = 0$ yields a novel form of the conventional (Helmholtz) free energy, immediately suggesting what can be considered as the phase-space density of free energy, namely the quantity 
$T_i |\tilde {g}|^2/F_0$, for which there is not yet an expression available.Footnote 4
It is useful now to write quasi-neutrality in terms of 
$\tilde {g}$
\begin{equation} \frac{e_i}{T_i}\delta\phi = \frac{\lambda}{n_i}\int \tilde{g} {\rm J}_{0} \,{\rm d}^3v,\end{equation}where
\begin{equation} \lambda = \frac{1}{\sqrt{(1+\tau -\varGamma_0)(1+\tau -\Delta\varGamma_0)}}.\end{equation} Finally, we can show that 
$\skew2\tilde{H}$ is positive–definite. First, positivity follows from (2.18), and it is obvious from (2.10) that if 
$g = 0$ then 
$\delta \phi = 0$ so that 
$E$ and 
$H$ both vanish, implying 
$\skew2\tilde{H} = 0$. On the other hand, if we assume that 
$\skew2\tilde{H} = 0$, then (2.18) implies that 
$\tilde {g} = 0$, and (2.22) implies that 
$\delta \phi = 0$, from which we conclude 
$g = 0$. In summary, 
$\skew2\tilde{H} \geqslant 0$ and 
$\skew2\tilde{H} = 0$ iff 
$g = 0$.
3 Modes of optimal growth
A key point in introducing the generalization of free energy 
$\skew2\tilde{H}$ is that this quantity introduces wave–particle effects (parallel resonance and drift resonance) that enter the electrostatic energy balance equation, (2.14).
We note that, for the choice 
$\varDelta = 0$, the energy 
$\skew2\tilde{H}$ reduces to the conventional Helmholtz free energy that we have studied in the previous papers. For this choice the modes of optimal growth correspond exactly to the modes introduced in § 6 of Part 1, and studied in fully electromagnetic limit in Plunk & Helander (Reference Plunk and Helander2022) (Part 2). Because those modes are included as a limit of our present analysis, we can be assured that the most stringent bound on growth obtained from the generalized free energy will be at least as good as the known bound obtained from the Helmholtz free energy.
Note that, as long as the parameter 
$\varDelta$ is independent of 
${\boldsymbol {k}}$, the quantity 
$\skew2\tilde{H}$ is conserved by the nonlinearity, i.e. under summation over 
${\boldsymbol {k}}$. This is because it is a linear combination of two nonlinear invariants. One simply combines (2.2) and (2.14) to obtain
\begin{equation} \frac{{\rm d}}{{\rm d} t}\sum_{\boldsymbol{k}} \skew2\tilde{H} = 2 \sum_{\boldsymbol{k}} (D - \Delta K),\quad\text{for $\varDelta$ independent of ${\boldsymbol{k}}$,}\end{equation}
i.e. the change of this measure is due to the drive terms of electrostatic and free energy, and is otherwise conserved by the turbulent interactions. It is potentially useful to also consider 
$\varDelta$ that does depend on 
${\boldsymbol {k}}$, for the purpose of obtaining bounds on linear growth, but the nonlinear implications will be less clear in that case.
In direct analogy to how modes of optimal free energy growth were defined, we introduce a rate 
$\varLambda$
\begin{equation} \varLambda = (D-\Delta K)/\skew2\tilde{H},\end{equation}
to be optimized over the space of ion distribution functions 
$g$. We note the bound on conventional gyrokinetic instability growth rates
\begin{equation} \gamma_L \leqslant \max_{g} \varLambda.\end{equation} Having already found a diagonal form of the generalized free energy, (2.18), we need not use a variational approach to find the states of extremal 
$\varLambda$. We simply identify the Hermitian linear operators associated with the input of free energy and electrostatic energy, i.e.
\begin{gather} D = (\tilde{g}, {\mathcal{D}} \tilde{g}), \end{gather}
\begin{gather}K = (\tilde{g}, {\mathcal{K}} \tilde{g}). \end{gather}
To obtain these forms note that, from (2.22), 
$\delta \phi$ can be regarded as the result of a linear operator on 
$\tilde {g}$. Then, using (2.22) and (2.3) and (2.15), and some straightforward algebra (see Appendix B), we obtain explicit forms for the operators. First we have
\begin{equation} {\mathcal{D}}\tilde{g} = \frac{{\rm i}\lambda}{2 n_i} {\rm J}_0 F_0 \eta\omega_\ast \int {\rm d}^3v^\prime \,{\rm J}_0^\prime \tilde{g}^\prime\left[ \left(\frac{v}{v_{{\rm th}}}\right)^2 - \left(\frac{v^\prime}{v_{{\rm th}}}\right)^2\right],\end{equation}
where primes denote evaluation at 
$\boldsymbol {v}^\prime$ and 
$v_{\textrm {th}} = \sqrt {2 T_i/m_i}$. For convenience, the operator 
${\mathcal {K}}$ can be split into its parallel and perpendicular components as 
${\mathcal {K}} = {\mathcal {K}}_{\|} + {\mathcal {K}}_d$, for which we obtain
\begin{equation} {\mathcal{K}}_d \tilde{g} = \frac{{\rm i}\lambda}{2 n_i}\omega_d(\ell) F_0 {\rm J}_0 \int {\rm d}^3v^\prime \,{\rm J}_0^\prime \tilde{g}^\prime \left[\left(\frac{v_{{\perp}}}{\sqrt{2}v_{{\rm th}}}\right)^2 + \left(\frac{v_{{\parallel}}}{v_{{\rm th}}}\right)^2 - \left(\frac{v_{{\perp}}^\prime}{\sqrt{2}v_{{\rm th}}}\right)^2 - \left(\frac{v_{{\parallel}}^\prime}{v_{{\rm th}}}\right)^2\right] ,\end{equation}and
\begin{align} {\mathcal{K}}_{\|} \tilde{g} &= \frac{\lambda}{2 n_i}F_0\left\{ {\rm J}_0\left[{-}B\frac{\partial}{\partial l}\left(\frac{1}{B}\int {\rm d}^3v^\prime \,v_{{\parallel}}^\prime {\rm J}_0^\prime \tilde{g}^\prime\right) + \int {\rm d}^3v^\prime \,v_{{\parallel}}^\prime \frac{\partial {\rm J}_0^\prime}{\partial l} \tilde{g}^\prime \right]\right.\nonumber\\ &\quad \left.+ \frac{v_{{\parallel}}}{\lambda}\frac{\partial}{\partial l}\left({\rm J}_0 \lambda\int {\rm d}^3v^\prime \,{\rm J}_0^\prime\tilde{g}^\prime \right) \right \}. \end{align}
In deriving (3.8), it is important to note that the parallel derivative is taken at fixed magnetic moment and particle energy, and that the velocity-space volume element 
$\textrm {d}^3v$ is proportional to 
$B/v_\|$ in these variables. More details are given in Appendix A. The kinetic eigenvalue problem can be stated now as
\begin{equation} \varLambda \tilde{g} = \left({\mathcal{D}} - \Delta {\mathcal{K}}\right)\tilde{g},\end{equation}
where solutions, i.e. pairs 
$\varLambda _n$ and 
$\tilde {g}_n(l, \boldsymbol {v})$, realize optimal growth of 
$\skew2\tilde{H}$. To see this, we can decompose the distribution function 
$\tilde {g}$ in terms of the orthogonal eigenmodes of (3.9), 
$\tilde {g} = \sum _n c_n \tilde {g}_n$ to obtain (choosing 
$||\tilde {g}_n|| = 1$)
\begin{equation} \frac{1}{2\skew2\tilde{H}}\frac{{\rm d} \skew2\tilde{H}}{{\rm d}t} = \frac{\sum_n |c_n|^2 \varLambda_n}{\sum_n |c_n|^2}, \end{equation}from which it is clear that the rate of energy growth of the system is maximized by setting the distribution function equal to the mode of largest growth rate, i.e.
\begin{equation} \max_{g} \frac{1}{2\skew2\tilde{H}}\frac{{\rm d} \skew2\tilde{H}}{{\rm d}t} = \max_{n} \varLambda_n = \varLambda_\text{max}, \end{equation}
which bounds the normal growth of the system, 
$\gamma _L$, in accordance with (3.3).
3.1 Moment form of eigenproblem
The analysis of (3.9) is greatly simplified by adopting a moment form. As found in the preceding papers, there are natural moments that appear in the energy input terms that can be identified to reduce the dimensionality of the problem substantially. Upon inspecting the energy balance equations one finds the following key dimensionless integrals:
\begin{gather} \kappa_1 = \int {\rm d}^3 v \,{\rm J}_0 \tilde{g}/n_i, \end{gather}
\begin{gather}\kappa_2 = \int {\rm d}^3 v \left(\frac{v^2}{v_{{\rm th}}^2}\right) {\rm J}_0 \tilde{g}/n_i, \end{gather}
\begin{gather}\kappa_3 = \int {\rm d}^3 v \left(\frac{v_{{\perp}}^2}{2v_{{\rm th}}^2} + \frac{v_{{\parallel}}^2}{v_{{\rm th}}^2}\right) {\rm J}_0 \tilde{g}/n_i, \end{gather}
\begin{gather}\kappa_4 = \int {\rm d}^3 v \left(\frac{v_{{\parallel}}}{v_{{\rm th}}}\right){\rm J}_0 \tilde{g}/n_i, \end{gather}
\begin{gather}\kappa_5 = \int {\rm d}^3 v \left(\frac{v_{{\parallel}}}{v_{{\rm th}}}\right) \frac{\partial {\rm J}_0}{\partial l} \tilde{g}/n_i, \end{gather}
where 
$\kappa _1$ is a density-like moment, 
$\kappa _2$ and 
$\kappa _3$ are pressure-like, 
$\kappa _4$ is parallel ion flow while 
$\kappa _5$ is more abstract.
It is easy to recognize these integrals on the right-hand side of (3.6), (3.7) and (3.8), and straightforward to rewrite those equations in moment form. The dimensional reduction is achieved by taking moments of the these equations to obtain a coupled set of five fluid equations. These, which are given in Appendix C, can be combined, leading, after lengthy algebra, to a relatively simple second-order ordinary differential equation, the main result of this paper
\begin{align} &\left(\frac{4\varLambda^2}{\lambda^2} + \left(\Delta \omega_d G_{3} - \eta\omega_\ast G_{1} \right)^2 - G_{0}\left[ (\eta\omega_\ast)^2 G_{2} - 2 \Delta \omega_d \eta\omega_\ast G_{4} + \varDelta^2\omega_d^2 G_{5}\right] \right) \varphi\nonumber\\ &\quad = \varDelta^2 v_{{\rm th}}^2 G_{0} B \left[ -\frac{\partial}{\partial l}\left(\frac{ G_{0,2}}{B}\frac{\partial \varphi}{\partial l}\right) + \varphi\frac{ G_{0,2}^{\prime\prime}}{B} - \varphi\frac{\partial}{\partial l}\left(\frac{ G_{0,2}^{\prime}}{B}\right) \right ], \end{align}
where 
$\varphi = e_i \delta \phi /T_i = \lambda \kappa _1$ is the normalized electrostatic potential. The functions 
$G_{m,n}$, 
$G^\prime _{m,n}$ and 
$G^{\prime \prime }_{m,n}$, which depend on arc length via 
$b(l)$ and 
$B(l)$, are defined in terms of integrals involving Bessel functions, and are evaluated in Appendix D; see, in particular, (D3)–(D4) and (D9)–(D12). The other 
$b$-dependent factors (
$G_0$–
$G_5$) can be expressed in terms of 
$G_{m,n}$, and are evaluated in terms of more elementary Bessel functions in Appendix D.1.
In (3.17) we see the eigenvalue 
$\varLambda$ entering quadratically, reflecting the fact that there will be two real roots, one positive and one negative, owing to Hermiticity and time-reversal symmetry of the full eigenproblem, (3.9). Note that the terms arising from the parallel drive of electrostatic energy are placed on the right-hand side. In the following section, we will consider some simple limits of this equation, and leave its more general solution for a future publication.
4 Simple limits
In this section we will consider some simple limits applied to (3.17), and draw some comparison with linear theory of the main instability targeted by limit of this paper, the ITG mode (see for instance Plunk et al. Reference Plunk, Helander, Xanthopoulos and Connor2014). To start, we note that taking 
$\varDelta = 0$, so that 
$\skew2\tilde{H}$ becomes the conventional Helmholtz free energy, yields
\begin{equation} \varLambda^2 = \frac{(\eta\omega_\ast)^2}{4(1+\tau -\varGamma_0)(1+\tau)} \left( G_{0}G_{2} - G_{1} ^2\right), \end{equation}which matches (6.20) of Part 1.
In considering other simplifications, we first should note that the adiabatic electron approximation already neglects the contribution from a trapped electron population, which requires either large electron collisionality, or uniform magnetic-field strength as measured along the field line. Let us assume the latter for simplicity
\begin{equation} \frac{\partial B}{\partial l} = 0. \end{equation}
Making this assumption simplifies somewhat (3.17), where all the explicit factors of 
$B$ drop out of the right-hand side. A more significant simplification is achieved by assuming unsheared and uniform magnetic geometry, in particular
\begin{gather} \frac{\partial b}{\partial l} = 0, \end{gather}
\begin{gather}\frac{\partial \omega_d}{\partial l} = 0. \end{gather}
In this limit, all of the coefficients of (3.17) are constants, and a simple dispersion relation is the obtained by taking 
$\partial \varphi /\partial l = \textrm {i} k_\| \varphi$. We find
\begin{equation} \frac{4\varLambda^2}{\lambda^2} + \left(\Delta \omega_d G_{3} - \eta\omega_\ast G_{1} \right)^2 - G_{0}\left[ (\eta\omega_\ast)^2 G_{2} - 2 \Delta \omega_d \eta\omega_\ast G_{4} + \varDelta^2\omega_d^2 G_{5}\right] = \varDelta^2 k_\|^2 v_{{\rm th}}^2 G_{0}^2/2,\end{equation}
were we have used 
$G_{0,2} = G_{0}/2$. As noted in § 2.2, the quantity 
$\varDelta$ is a free parameter, over which we can optimize 
$\varLambda$ to improve the bounds on the growth rate of fluctuations.
4.1 Slab ITG mode
Setting 
$\omega _d = 0$ leaves only the slab branch of the ITG mode, driven by the temperature gradient, and involving ion parallel resonance. Equation (4.5) reduces to
\begin{equation} \frac{4\varLambda^2}{(\eta\omega_\ast)^2\lambda^2} = G_{0}G_{2} - G_{1}^2 + \varDelta^2 \kappa_\|^{{-}2} G_{0}^2/2,\end{equation}
where 
$\kappa _\| = \eta \omega _\ast /(k_\|v_{\textrm {th}})$. Because 
$G_0G_2 - G_1^2 \geqslant 0$, the two contributions on the right-hand side are both positive but the solution for which 
$\varLambda$ is minimal is actually not obtained for 
$\varDelta = 0$, due to the implicit dependence of 
$\lambda$ on 
$\varDelta$ given by (2.23).
Although all values of the parameter 
$\varDelta$ satisfying (2.21) yield a valid bound on the growth rate of normal modes (
$\gamma _L \leqslant \varLambda$), the lowest value is the most stringent and serves as the closest approximation of 
$\gamma _L$. To obtain this ‘optimal bound’, we can consider the extrema of 
$\varLambda ^2/(\eta \omega _\ast )^2$, i.e.
\begin{equation} \frac{{\rm d}}{{\rm d}\varDelta}\left( \frac{G_{0}G_{2} - G_{1}^2 + \varDelta^2 \kappa_\|^{{-}2} G_{0}^2/2}{(1+\tau -\varGamma_0)(1+\tau -\Delta\varGamma_0)}\right) = 0.\end{equation}
This results in a quadratic equation for 
$\varDelta$ that is still rather complicated so we will consider the limit 
$b \rightarrow 0$; see Appendix D.2 for the relevant limits of 
$G_{m,n}$, etc. Applying the limit to (4.6) yields
\begin{equation} \frac{\varLambda^2}{(\eta\omega_\ast)^2} = \frac{3 + \varDelta^2/\kappa_\|^2}{8\tau(1 + \tau - \varDelta)}. \end{equation} This solution diverges as 
$\varDelta$ approaches 
$1+\tau$; recall that this is the upper limit allowed by (2.21). It also grows in an unbounded fashion as 
$\varDelta \rightarrow -\infty$. There is an optimal value giving minimal 
$|\varLambda |$, obtained by solving (4.7) in this limit. This solution, denoted as 
$\varDelta _\textrm {min}$, is
\begin{equation} \varDelta_\text{min} = 1 + \tau - \sqrt{(1+\tau)^2 + 3 \kappa_\|^2}, \end{equation}where the negative root has been selected to be consistent with (2.21). Substituting this solution into (4.8) gives
\begin{equation} \varLambda_{\min}^2 = \frac{(\eta\omega_\ast)^2}{4\bar{\kappa}_\|^2 \tau(1+\tau)}\left( \sqrt{1 + 3 \bar{\kappa}_\|^2} - 1\right),\end{equation}
where we define 
$\bar {\kappa }_\| = \kappa _\|/(1+\tau )$; see figure 1. This reaches its maximum value in the limit 
$\bar {\kappa }_\| \rightarrow 0$, and is a decreasing function of 
$|\bar {\kappa }_\||$, i.e.
\begin{equation} \varLambda_{\min}^2 = \begin{cases} \dfrac{3}{8 \tau (\tau +1)}(\eta\omega_\ast)^2, & \text{for $\bar{\kappa}_\| \rightarrow 0$,}\\[8pt] \dfrac{\sqrt{3}}{4\tau}|\eta\omega_\ast k_\parallel v_{{\rm th}}|, & \text{for $|\bar{\kappa}_\|| \gg 1$}. \end{cases} \end{equation}
Physically, the first result implies that, when drive (
$\eta \omega _\ast$) is much smaller than the parallel transit frequency (
$k_\parallel v _{\textrm {th}}$), the best bound is equal to that obtained by free energy (
$\varDelta = 0$). In this case, the bound is consistent from expectations of the growth rate of a resonant slab ITG mode, i.e. 
$\gamma _L \sim \eta \omega _\ast$.
Figure 1. Bound of the growth rate of the slab ITG mode (blue), obtained from the optimal growth of the generalized free energy, plotted vs the instability parameter 
$\bar {\kappa }_\| = \eta \omega _\ast /[k_\parallel v _{\textrm {th}} (1+\tau )]$. For comparison, the growth rate 
$\gamma _L$ is obtained by solving the linear dispersion relation (E1), and the quantity 
$\tau (1+\tau ) (\gamma _L/\eta \omega _*)^2$ is plotted for the case 
$\tau = 1$ (dashed black).
In the opposite limit (
$\bar {\kappa }_\| \gg 1$), however, when the drive large, i.e. in the so-called non-resonant or ‘fluid’ limit, we obtain a much lower bound, essentially the geometric mean of the drive and the parallel transit frequency 
$k_\|v_{\textrm {th}}$. We note that this bound is not as low as what is obtained from the non-resonant solution of the dispersion relation (without a density gradient), i.e. 
$\gamma _L \sim \eta \omega _\ast ^{1/3}(k_\parallel v _{\textrm {th}})^{2/3}$ (Plunk et al. Reference Plunk, Helander, Xanthopoulos and Connor2014), but nevertheless captures the expected weakening (relative to the resonant result) qualitatively. A plot of the optimal bound of (4.10) is provided in figure 2. Note that the resonant stabilization at high 
$k_\|$ (low 
$\bar {\kappa }_\|$) is not captured in this case.
Figure 2. Bound of the growth rate of the toroidal ITG mode, obtained from the optimal growth of the generalized free energy, plotted vs the instability parameter 
$\bar {\kappa }_d = \eta \omega _\ast /[\omega _d (1+\tau )]$. For comparison, the growth rate 
$\gamma _L$ is obtained by solving the linear dispersion relation (E1), and the quantity 
$\tau (1+\tau ) (\gamma _L/\eta \omega _*)^2$ is plotted for the case 
$\tau = 1$ (dashed black).
It is interesting to observe that this latter limit corresponds to 
$\varDelta _\textrm {min} \rightarrow -\infty$, making 
$\skew2\tilde{H}$ in some sense dominated by the electrostatic component.
4.2 Toroidal ITG mode
Now taking 
$k_\parallel v _{\textrm {th}}$ to be small, we can neglect the right-hand side of (4.5), leaving
\begin{equation} \frac{4\varLambda^2}{\lambda^2} = G_{0}\left[ (\eta\omega_\ast)^2 G_{2} - 2 \Delta \omega_d \eta\omega_\ast G_{4} + \varDelta^2\omega_d^2 G_{5}\right] - \left(\Delta \omega_d G_{3} - \eta\omega_\ast G_{1} \right)^2.\end{equation}
To derive the optimal choice of 
$\varDelta$, we again take the 
$b \rightarrow 0$ limit and obtain from (4.12)
\begin{equation} \frac{\varLambda^2}{(\eta\omega_\ast)^2} = \frac{3 \varDelta ^2-8 \Delta \kappa_d +6 \kappa_d^2}{16 \kappa_d^2 \tau (\tau + 1 - \varDelta)},\end{equation}
where we define 
$\kappa _d = \eta \omega _\ast /\omega _d$. Note the similar qualitative behaviour with 
$\varDelta$ as (4.6), namely its divergence at the 
$\varDelta \rightarrow 1 + \tau$, and unbounded growth as 
$\varDelta \rightarrow -\infty$. The key difference here arises in the linear term in the drive parameter 
$\kappa$; this is expected from the theory of the toroidal ITG mode since the sign of the drift frequency (associated with so-called ‘good’ and ‘bad’ magnetic curvature) is important for the resonance.
We now find the value of 
$\varDelta$ that minimizes 
$\varLambda$
\begin{equation} \varDelta_\text{min} = (1+\tau) \left(1 - \sqrt{2 \bar{\kappa}_d^2-8\bar{\kappa}_d/3 +1}\right), \end{equation}
where we define the parameter 
$\bar {\kappa }_d = \kappa _d/(1 + \tau )$. Substituting this into (4.13) yields
\begin{equation} \varLambda_\text{min}^2 = (\eta\omega_\ast)^2\left( \frac{8 \bar{\kappa}_d \left(\zeta \left(\bar{\kappa}_d\right)-1\right)+3 \left(\zeta \left(\bar{\kappa}_d\right)-1\right)^2+6 \bar{\kappa}_d^2}{16 \tau (\tau +1) \bar{\kappa}_d^2 \zeta \left(\bar{\kappa}_d\right)} \right),\end{equation}
with 
$\zeta = \sqrt {2 \bar {\kappa }_d^2-8\bar {\kappa }_d/3 + 1}$. This expression for 
$\varLambda _\textrm {min}$ is naturally separated into a factor that depends only on the instability parameter 
$\bar {\kappa }_d$, from which we can derive the asymptotic behaviour. To show the behaviour of this factor we plot the quantity 
$\tau (1 + \tau )\varLambda _\textrm {min}^2/(\eta \omega _\ast )^2$ in figure 2. The overall behaviour of 
$\varLambda _\textrm {min}$ is captured by the following limits:
\begin{equation} \varLambda_\text{min}^2 = \begin{cases} |\eta\omega_\ast\omega_d|\left(\dfrac{3 \sqrt{2} + 4 \sigma}{8 \tau}\right), & \text{for $|\bar{\kappa}_d| \gg 1$,}\\[10pt] \dfrac{(\eta\omega_\ast)^2}{24\tau(1+\tau)}, & \text{for $\bar{\kappa}_d \rightarrow 0$},\\[8pt] \dfrac{3(\eta\omega_\ast)^2}{8\tau(1+\tau)}, & \text{for $\bar{\kappa}_d \rightarrow 4/3$}, \end{cases} \end{equation}
where we denote 
$\sigma = \pm 1$ as the sign of 
$\kappa _d$. At large drive (
$|\bar {\kappa }_d| \gg 1$; 
$|\eta \omega _\ast | \gg |\omega _d| (1+\tau )$) we recover the expected non-resonant (‘fluid’) behaviour of the toroidal ITG mode with no density gradient, namely 
$\gamma _L \sim \sqrt {\eta \omega _\ast \omega _d}$. Note that this growth rate is much smaller than the bound found by merely considering the Helmholtz free energy (Helander & Plunk Reference Helander and Plunk2022). Although we do not see the complete stabilization (
$\varLambda = 0$) at negative values of 
$\bar {\kappa }_d$ (opposite sign of 
$\omega _d$ and 
$\eta \omega _\ast$) expected from theory, there is a strong asymmetry, with 
$|\varLambda |$ having its larger values at positive 
$\bar {\kappa }_d$ and being comparatively much smaller for negative 
$\bar {\kappa }_d$.
The value 
$\bar {\kappa }_d = 4/3$ (i.e. 
$\eta \omega _\ast = 4 (1+\tau )\omega _d/3$) achieves the maximal value of 
$\varLambda _\textrm {min}$ at fixed 
$\eta \omega _\ast$, and therefore is evocative of the resonance condition for the toroidal ITG modes 
$\eta \omega _\ast \sim \omega _d$ (Biglari, Diamond & Rosenbluth Reference Biglari, Diamond and Rosenbluth1989). This value of 
$\bar {\kappa }_d$ is obtained by solving for 
$\varDelta _\textrm {min} = 0$, explaining why it produces the worst bound, i.e. that given by optimal growth of Helmholtz free energy.
It is noteworthy that for the limit 
$\bar {\kappa }_d \rightarrow 0$ (
$\omega _d \gg \eta \omega _\ast /(1+\tau )$) our method yields a value of 
$|\varLambda |$ that is a factor of 
$1/3$ reduced as compared with the resonant case, again at least qualitatively reproducing the expected stabilization of the toroidal ITG mode in this limit.
5 Conclusion
We have demonstrated that the use of a generalized form of free energy 
$\skew2\tilde{H}$ introduces some of the physics of wave–particleresonance that is missing in the theory of optimal mode growth of the Helmholtz free energy (Helander & Plunk Reference Helander and Plunk2021, Reference Helander and Plunk2022; Plunk & Helander Reference Plunk and Helander2022). The growth rates of optimal modes of the generalized free energy provide a rigorous upper bound on the growth of conventional gyrokinetic instabilities (‘normal modes’), which, for appropriate choice of 
$\Delta$ is less than or equal to the Helmholtz bound. Moreover, optimal modes of generalized free energy depend on the magnetic-field geometry to a greater extent than those associated with Helmholtz free energy. The difference in growth rates can be very large. For instance, in the important case of a strongly driven toroidal ITG mode, the Helmholtz bound is larger by a factor of order 
$\eta \omega _\ast / \omega _d \gg 1$.
A single ordinary differential equation has been derived for optimal modes, allowing general magnetic geometry. We found solutions of this equation in some simple limits to demonstrate that it indeed recovers, at least qualitatively, some of the physical effects expected from the theory of linear ITG modes, including sensitivity to the ratio of the frequencies associated with drive and resonance, and transition of the instability when this ratio is near one. Density gradient dependence of the ITG mode is absent from both the electrostatic and free energy input terms, assuming adiabatic electrons, so its effect is not accounted for by the theory presented here.
The results of this work have possible implications for ‘turbulence optimization’, i.e. the endeavour to shape the equilibrium magnetic geometry of stellarators for low turbulence. The general result, (3.17), allows, in principle, for the inclusion of the complete geometric information that is needed to run gyrokinetic simulations. However, the solution of this equation should be far simpler and more efficient, due to the reduction of velocity space to a single moment. Stellarator optimization requires a large number of calculations to be performed, which so far has severely limited the use of direct gyrokinetic simulations, even linear ones. The simplifications obtained here, which also include the removal of the need for time integration, will need to be balanced against numerical costs associated with eigenvalue problems, and possible costs associated with determining the optimal value of the parameter 
$\varDelta$, if the absolute strongest bound is desired.
Our results found for the toroidal branch of the ITG mode hint at a possible optimization strategy based on optimal modes. Consider fixed plasma conditions, i.e. a given temperature gradient (
$\eta \omega _\ast$) and temperature ratio (
$\tau$): at high drive (
$\eta \omega _\ast / (1 + \tau ) > 4\omega _d/3$), minimization of the optimal growth rate 
$|\varLambda |$ is achieved by minimization of the magnetic drift 
$\omega _d$ (i.e. magnetic curvature), corresponding to minimization of the strongly driven(non-resonant) toroidal ITG mode. On the other hand, at low drive (
$\eta \omega _\ast /(1 + \tau ) < 4\omega _d/3$), the increase of 
$\omega _d$ is favoured, corresponding to a weakening of the marginally unstable ITG mode, i.e. an increase of the threshold of instability. The latter case corresponds to ‘critical gradient’ optimization, an idea which has recently been developed (Roberg-Clark, Plunk & Xanthopoulos Reference Roberg-Clark, Plunk and Xanthopoulos2022a; Roberg-Clark, Xanthopoulos & Plunk Reference Roberg-Clark, Xanthopoulos and Plunk2022b).
It is worth mentioning the larger context in which optimal modes are potentially interesting to study, and applications besides their use to bound or estimate the growth rate of normal instabilities. Although these modes (
$\tilde {g}_n$, 
$\varLambda _n$) do not arise as late-time-asymptotic solutions of initial-value gyrokinetic simulations, as normal modes do, they are still realizable in the sense that a gyrokinetic simulation initialized in one of these states 
$\tilde {g}_n$ will temporarily exhibit an energy growth rate exactly equal to 
$2\varLambda _n$. These modes are thus just as ‘real’ as conventional eigenmodes and may also be observed, to some extent, in nonlinear simulations, where fluctuations are continually driven away by nonlinear interactions from the form of normal linear instabilities. They are also potentially useful for studying systems without any unstable normal modes, i.e. ‘subcritical’ turbulence (indeed this is the context in which they were originally formulated). Furthermore, since optimal modes are based on quadratic norms that are nonlinearly conserved by the gyrokinetic equations, they can be used to bound the instantaneous growth rates observed in fully nonlinear solutions; this point was made in Parts 1 and 2, but applies equally well here.
More general solutions of the optimal mode equation, and the application to optimization, will be pursued in future works. Other special limits can also be explored including the limit of large electron-bounce frequency, appropriate for studying trapped-electron modes, or adiabatic-ion limits applied to universal instabilities. Electromagnetic generalizations are also possible: although it is not clear how to construct a positive–definite electromagnetic form of the generalized free energy 
$\skew2\tilde{H}$ that is a nonlinear invariant, it is certainly possible to consider related measures that focus on linear bounds.
Acknowledgements
Editor Alex Schekochihin thanks the referees for their advice in evaluating this article.
Funding
This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 – EUROfusion). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. This work was partly supported by a grant from the Simons Foundation (560651, PH).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Several kinetic species
For a plasma with an arbitrary number of particle species, we multiply each gyrokinetic equation
\begin{gather} \frac{\partial g_{a,{\boldsymbol{k}}}}{\partial t} + v_{\|} \frac{\partial g_{a,{\boldsymbol{k}}}}{\partial l} + {\rm i} \omega_{da} g_{a, {\boldsymbol{k}}} + \frac{1}{B^2} \sum_{{\boldsymbol{k}}'} \boldsymbol{B} \boldsymbol{\cdot} ({\boldsymbol{k}} \times {\boldsymbol{k}}') \overline{\delta\phi}_{{\boldsymbol{k}}'} g_{a, {\boldsymbol{k}} - {\boldsymbol{k}}'} \end{gather}
\begin{gather}= \frac{e_a F_{a0}}{T_a} \left( \frac{\partial}{\partial t} + {\rm i} \omega_{{\ast} a}^T \right) \overline{\delta\phi}_{\boldsymbol{k}} , \end{gather}
by 
$e_a \overline {\delta \phi }_{\boldsymbol {k}}^\ast$, integrate over velocity space, take the real part and the average 
$\langle \cdots \rangle$ over the flux tube and sum over all species 
$a$ and wave vectors 
$\boldsymbol {k}$. In other words, we apply the operator
\begin{equation} {\rm Re} \, \sum_{a,{\boldsymbol{k}}} e_a \left\langle \int \overline{\delta\phi}_{\boldsymbol{k}}^\ast \left( \cdots \right) {\rm d}^3v \right\rangle. \end{equation}Since the expression
\begin{gather} {\rm Re} ({\boldsymbol{k}} \times {\boldsymbol{k}}') \overline{\delta\phi}_{{\boldsymbol{k}}'}^\ast \overline{\delta\phi}_{\boldsymbol{k}} g_{a,{\boldsymbol{k}} - {\boldsymbol{k}}'} = \tfrac{1}{2} ({\boldsymbol{k}} \times {\boldsymbol{k}}') \left( \overline{\delta\phi}_{{\boldsymbol{k}}'}^\ast \overline{\delta\phi}_{\boldsymbol{k}} g_{a,{\boldsymbol{k}} - {\boldsymbol{k}}'} + \overline{\delta\phi}_{{\boldsymbol{k}}'} \overline{\delta\phi}_{\boldsymbol{k}}^\ast g_{a,{\boldsymbol{k}} - {\boldsymbol{k}}'}^\ast \right) \end{gather}
\begin{gather}= \tfrac{1}{2} ({\boldsymbol{k}} \times {\boldsymbol{k}}') \left( \overline{\delta\phi}_{-{\boldsymbol{k}}'} \overline{\delta\phi}_{\boldsymbol{k}} g_{a,{\boldsymbol{k}} - {\boldsymbol{k}}'} + \overline{\delta\phi}_{{\boldsymbol{k}}'} \overline{\delta\phi}_{-\boldsymbol{k}} g_{a,{-\boldsymbol{k}} + {\boldsymbol{k}}'} \right) \end{gather}
changes sign under an exchange of 
$\boldsymbol {k}$ and 
${\boldsymbol {k}}'$, the nonlinear terms cancel upon summation over 
$\boldsymbol {k}$ and 
${\boldsymbol {k}}'$, and we obtain
\begin{equation} {\rm Re} \, \sum_{a,{\boldsymbol{k}}} e_a \left\langle \int \overline{\delta\phi}_{\boldsymbol{k}}^\ast \left( \frac{\partial g_{a,{\boldsymbol{k}}}}{\partial t} + v_{\|} \frac{\partial g_{a,{\boldsymbol{k}}}}{\partial l} + {\rm i} \omega_{da} g_{a, {\boldsymbol{k}}} - \frac{e_a F_{a0}}{T_a} \frac{\partial \overline{\delta\phi}_{\boldsymbol{k}}}{\partial t} \right) {\rm d}^3v \right\rangle = 0. \end{equation}The quasi-neutrality equation (2.10) can be used to write the first term as
\begin{equation} {\rm Re} \, \sum_{a,{\boldsymbol{k}}} e_a \int \overline{\delta\phi}_{\boldsymbol{k}}^\ast \frac{\partial g_{a,{\boldsymbol{k}}}}{\partial t} {\rm d}^3v = \frac{{\rm d}}{{\rm d}t} \sum_{\boldsymbol{k}} \frac{n_a e_a^2}{2 T_a} \left| \delta \phi_{\boldsymbol{k}} \right|^2. \end{equation}We thus arrive at the electrostatic energy balance equation
\begin{gather} \frac{{\rm d}}{{\rm d}t} \sum_{\boldsymbol{k}} \frac{n_a e_a^2}{2 T_a} \left\langle \left[1 - \varGamma_0(b_{a{\boldsymbol{k}}}) \right] \left| \delta \phi_{\boldsymbol{k}} \right|^2 \right\rangle \end{gather}
\begin{gather}={-} {\rm Re} \, \sum_{a,\boldsymbol{k}} \left\langle e_a \int \overline{\delta\phi}_{\boldsymbol{k}}^\ast \left( v_\| \frac{\partial g_{a,{\boldsymbol{k}}}}{\partial l} + {\rm i} \omega_{da} g_{a,{\boldsymbol{k}}} \right) {\rm d}^3v \right\rangle, \end{gather}which is the generalization of (2.14) to several species. The right-hand side can be interpreted as minus the work done by the electric field on the various particle species. The first term in this expression contains
\begin{align} \int {\rm J}_{0a} v_\| \frac{\partial g_{a,{\boldsymbol{k}}}}{\partial l} {\rm d}^3v &= \sum_\sigma \frac{2 {\rm \pi}\sigma B}{m_a^2} \int_0^\infty {\rm d}E_a \int_0^{E_a/B} {\rm J}_{0a} \frac{\partial g_{a,{\boldsymbol{k}}}}{\partial l} {\rm d}\mu_a \end{align}
\begin{align} &= B \frac{\partial}{\partial l} \left( \frac{1}{B} \int {\rm J}_{0a} v_\| g_{a,{\boldsymbol{k}}} \,{\rm d}^3v \right) - \int \frac{\partial {\rm J}_{0a}}{\partial l} v_\| g_{a,{\boldsymbol{k}}} \,{\rm d}^3v, \end{align}
which we used in (3.8), with 
$\sigma = v_\| / | v_\| |$, 
$E_a = m_av^2/2$ and 
$\mu _a = m_a v_\perp ^2/(2B)$.
Appendix B. Derivation of operators 
${\mathcal {D}}$ and ${\mathcal {K}}$
We first write forms of 
$D$ and 
$K$ explicit in 
$\tilde {g}$, noting that the contribution to 
$D$ proportional to the density gradient is zero by use of quasi-neutrality with the adiabatic electron approximation (the factor 
$\eta \omega _\ast \propto \textrm {d}T_i/\textrm {d}\psi$ appears in what follows, but never 
$\omega _\ast \propto \textrm {d}n_i/\textrm {d}\psi$ alone). The terms proportional to 
$\nu$, involved in the transformation from 
$g$ to 
$\tilde {g}$, are also zero, due to oddness in 
$v_\|$, so expressing energy input in terms of 
$\tilde {g}$ merely has the consequence of introducing the overall factor 
$\lambda$. The expressions are
\begin{gather} D ={-} \left\langle T_i\frac{{\rm i}\lambda }{2 n_i} \, \int \tilde{g}(\boldsymbol{v})\tilde{g}^*(\boldsymbol{v}^\prime)\, \eta\omega_\ast \left(\frac{v^2}{v_{{\rm th}}^2}\right) {\rm J}_0 {\rm J}_0^\prime \,{\rm d}^3v \,{\rm d}^3v^{\prime}\right\rangle + \text{c.c.}, \end{gather}
\begin{gather}K_\| ={-} \left\langle T_i\frac{\lambda}{2n_i} \, \int \tilde{g}^*(\boldsymbol{v}^\prime)\left(v_{\|}\frac{\partial \tilde{g}(\boldsymbol{v})}{\partial l}\right) {\rm J}_0 {\rm J}_0^\prime \,{\rm d}^3v \,{\rm d}^3v^\prime \right\rangle + \text{c.c.}, \end{gather}
\begin{gather}K_d ={-} \left\langle T_i\frac{i\lambda}{2n_i} \, \int \tilde{g}^*(\boldsymbol{v}^\prime)\omega_{d} \tilde{g}(\boldsymbol{v}) {\rm J}_0 {\rm J}_0^\prime \,{\rm d}^3v \,{\rm d}^3v^\prime \right\rangle + \text{c.c.}, \end{gather}
where we separate 
$K = K_\| + K_d$ and ‘c.c.’ denotes complex conjugate. These can be re-expressed in terms of linear operators by writing them in the form
\begin{gather} D = (\tilde{g}, {\mathcal{D}}\tilde{g}) = \left\langle T_i\int {\rm d}^3v\frac{\tilde{g}^*}{F_0} {\mathcal{D}}\tilde{g} \right\rangle, \end{gather}
\begin{gather}K_\| = (\tilde{g}, {\mathcal{K}}_\|\tilde{g}) = \left\langle T_i\int {\rm d}^3v\frac{\tilde{g}^*}{F_0} {\mathcal{K}}_\|\tilde{g} \right\rangle, \end{gather}
\begin{gather}K_d = (\tilde{g}, {\mathcal{K}}_d\tilde{g}) = \left\langle T_i\int {\rm d}^3v\frac{\tilde{g}^*}{F_0} {\mathcal{K}}_d\tilde{g} \right\rangle. \end{gather}
Identifying 
${\mathcal {D}}$ and 
${\mathcal {K}}_d$ is simply a matter of exchanging labels of dummy variables of integration (
$v$ for 
$v^\prime$, etc.). Manipulating the expression for 
$K_\|$ to reveal 
${\mathcal {K}}_\|$ is more involved. We also need to integrate by parts in 
$l$ and will need to use the fact that 
$\partial /\partial l$ is performed at fixed phase-space variables 
$E_i$ and 
$\mu = m_i v_{\perp }^2/(2B)$. The velocity-space volume element contains an important factor of 
$1/v_\|$, which generally depends on 
$l$ and does not itself commute with 
$\partial /\partial l$
\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_i \,{\rm d}\mu_i}{m_i^2 | v_\| |}, \end{equation}
where 
$\sigma$ denotes the sign of 
$v_\|$.
Appendix C. Moment form of eigenproblem
The three terms of (3.9) can be rewritten in terms of the moments of 
$\tilde {g}$ (3.12)
\begin{gather} {\mathcal{D}} \tilde{g} = \frac{{\rm i}\lambda}{2}\eta\omega_\ast {\rm J}_0F_0\left[ \frac{v^2}{v_{{\rm th}}^2}\kappa_1 - \kappa_2\right], \end{gather}
\begin{gather}{\mathcal{K}}_\| \tilde{g} = \frac{\lambda}{2}F_0\left[ v_{{\rm th}} {\rm J}_0\left({-}B\frac{\partial}{\partial l} \left(\frac{\kappa_4}{B}\right) + \kappa_5 \right) + \frac{v_\|}{\lambda} \frac{\partial}{\partial l}\left({\rm J}_0 \lambda \kappa_1\right)\right], \end{gather}
\begin{gather}{\mathcal{K}}_d \tilde{g} = \frac{{\rm i}\lambda}{2} \omega_d {\rm J}_0F_0\left[ \left(\frac{v_\perp^2}{2v_{{\rm th}}^2} + \frac{v_\|^2}{v_{{\rm th}}^2}\right)\kappa_1 - \kappa_3\right]. \end{gather}Then, taking moments of (3.9) yields the following five equations:
\begin{gather} \frac{2\varLambda}{\lambda} \kappa_1 ={\rm i}\eta\omega_\ast\left(G_1\kappa_1 - G_0\kappa_2\right) - {\rm i}\Delta\omega_d\left(G_3\kappa_1 - G_0\kappa_3\right) -\Delta G_0v_{{\rm th}} \left(\kappa_5 - B\frac{\partial}{\partial l}\left(\frac{\kappa_4}{B}\right)\right), \end{gather}
\begin{gather}\frac{2\varLambda}{\lambda} \kappa_2 ={\rm i}\eta\omega_\ast\left(G_2\kappa_1 - G_1\kappa_2\right) - {\rm i}\Delta\omega_d\left( G_2 \kappa_1 - G_1 \kappa_3 \right) -\Delta G_1v_{{\rm th}} \left(\kappa_5 - B\frac{\partial}{\partial l}\left(\frac{\kappa_4}{B}\right)\right), \end{gather}
\begin{gather}\frac{2\varLambda}{\lambda} \kappa_3 ={\rm i}\eta\omega_\ast\left(G_4\kappa_1 - G_3\kappa_2\right) - {\rm i}\Delta\omega_d\left( G_5 \kappa_1 - G_3\kappa_3 \right) -\Delta G_3v_{{\rm th}} \left(\kappa_5 - B\frac{\partial}{\partial l}\left(\frac{\kappa_4}{B}\right)\right), \end{gather}
\begin{gather}\frac{2\varLambda}{\lambda} \kappa_4 ={-}\Delta v_{{\rm th}}\left( G_{0,2}^\prime\kappa_1 + \frac{G_{0,2}}{\lambda}\frac{\partial }{\partial l}(\lambda \kappa_1)\right), \end{gather}
\begin{gather}\frac{2\varLambda}{\lambda} \kappa_5 ={-}\Delta v_{{\rm th}}\left( G_{0,2}^{\prime\prime}\kappa_1 + \frac{G_{0,2}^\prime}{\lambda}\frac{\partial}{\partial l}(\lambda \kappa_1)\right). \end{gather}
See the next section where the integrals 
$G_{m,n}$, etc., are defined and evaluated. Note that the final two equations can be immediately used to eliminate 
$\kappa _4$ and 
$\kappa _5$, leaving a system of three equations. The second and third equations are used together to find forms for 
$\kappa _2$ and 
$\kappa _3$ in terms of 
$\kappa _1$, and these forms are substituted into the first equation to obtain the final form, in terms of 
$\kappa _1$ only, given by (3.17).
Appendix D. Bessel-type integrals
The following definitions, mostly copied from Plunk & Helander (Reference Plunk and Helander2022), are needed to perform the various integrals that appear in the moment equations for our eigenproblem. First, we need a general form of Weber's integral (Iyanaga & Kawada Reference Iyanaga and Kawada1980)
\begin{align} {\mathcal{I}}_{n}(p,a_1,a_2) &= \int_0^\infty \exp({-}p t^2) {\rm J}_{n}(a_1 t){\rm J}_{n}(a_2 t) t \,{\rm d}t\nonumber\\ &= \frac{1}{2p} \exp\left( \frac{-a_1^2-a_2^2}{4p}\right) {\rm I}_{n}\left(\frac{a_1 a_2}{2p}\right), \end{align}
where 
$\textrm {I}_{n}$ is the modified Bessel function of order 
$n$. The integrals we need to evaluate can be conveniently found in terms of 
${\mathcal {I}}_{n}$. We define
\begin{gather} G_{{\perp} m}(b) = 2 \int_0^{\infty} x_\perp^{m+1} \exp({-}x_\perp^2) {\rm J}_0^2(\sqrt{2 b} x_\perp) \,{\rm d}x_\perp, \end{gather}
\begin{gather}G_{{\perp} m}^{(1)}(b) = 2 \int_0^{\infty} x_\perp^{m+2} \exp({-}x_\perp^2) {\rm J}_0(\sqrt{2 b} x_\perp){\rm J}_1 (\sqrt{2 b} x_\perp) \,{\rm d}x_\perp, \end{gather}
\begin{gather}G_{{\perp} m}^{(2)}(b) = 2 \int_0^{\infty} x_\perp^{m+3} \exp({-}x_\perp^2) {\rm J}_1^2(\sqrt{2 b} x_\perp) \,{\rm d}x_\perp, \end{gather}
where 
$m$ is assumed to be even. Now we note that these integrals can be evaluated in terms of Weber's integral
\begin{gather} G_{{\perp} m}(b) = 2\left[ \left(-\frac{{\rm d}}{{\rm d}p}\right)^{m/2} {\mathcal{I}}_0(p,\sqrt{2b},\sqrt{2b})\right]_{p=1}, \end{gather}
\begin{gather}G_{{\perp} m}^{(1)}(b) = 2\left[ \left(-\frac{{\rm d}}{{\rm d}p}\right)^{m/2} \left(-\frac{{\rm d}}{{\rm d}\lambda}\right) { \mathcal{I}}_0(p,\lambda,\sqrt{2b})\right]_{p=1, \lambda = \sqrt{2b}}, \end{gather}
\begin{gather}G_{{\perp} m}^{(2)}(b) = 2\left[ \left(-\frac{{\rm d}}{{\rm d}p}\right)^{m/2}\left(-\frac{{\rm d}}{{\rm d}\lambda_1}\right)\left(-\frac{d}{d\lambda_2}\right) {\mathcal{I}}_0(p,\lambda_1,\lambda_2)\right]_{p=1, \lambda_1 = \lambda_2 = \sqrt{2b}}. \end{gather}The above relations allow us to evaluate the functions
\begin{gather} G_{m,n}(b) = G_{{\perp} m}(b) G_{{\parallel} n}, \end{gather}
\begin{gather}G_{m,n}^{(1)}(b) = G_{{\perp} m}^{(1)}(b) G_{{\parallel} n}, \end{gather}
\begin{gather}G_{m,n}^{(2)}(b) = G_{{\perp} m}^{(2)}(b) G_{{\parallel} n}, \end{gather}where
\begin{equation} G_{{\parallel} n} = \frac{1}{\sqrt{\rm \pi}}\int_{-\infty}^{\infty} \exp({-}x_\parallel^2) x_\parallel^n \,{\rm d}x_\parallel{=} \frac{1+({-}1)^n}{2\sqrt{\rm \pi}} \varGamma_E\left(\frac{1+n}{2}\right), \end{equation}
and 
$\varGamma _E$ is the gamma function. Finally, we can evaluate the integrals 
$G_{m,n}^{\prime }$ and 
$G_{m,n}^{\prime }$. We define
\begin{gather} G_{{\perp} m}^{\prime}(b) = 2 \int_0^{\infty} x_\perp^{m+1} \exp({-}x_\perp^2) {\rm J}_0 \frac{\partial {\rm J}_0}{\partial l} {\rm d}x_\perp, \end{gather}
\begin{gather}G_{{\perp} m}^{\prime\prime}(b) = 2 \int_0^{\infty} x_\perp^{m+1} \exp({-}x_\perp^2) \left(\frac{\partial {\rm J}_0}{\partial l}\right)^2 \,{\rm d}x_\perp. \end{gather}
Relating 
$x_\perp$ to the proper gyrokinetic phase-space variable 
$\mu$, that is 
$x_\perp = (\mu B/T_i)^{1/2}$, allows the derivatives to be evaluated
\begin{equation} \frac{\partial {\rm J}_0}{\partial l} ={-}\frac{1}{2} x_\perp\sqrt{2b} {\rm J}_1\frac{\partial}{\partial l}\ln(b B), \end{equation}so that we can write
\begin{gather} G_{{\perp} m}^{\prime}(b) ={-}\sqrt{b/2} \left(\frac{\partial}{\partial l}\ln(b B) \right) G_{{\perp} m}^{(1)}(b), \end{gather}
\begin{gather}G_{{\perp} m}^{\prime\prime}(b) = b/2 \left(\frac{\partial}{\partial l}\ln(b B)\right)^2 G_{{\perp} m}^{(2)}(b). \end{gather}These expressions allow us to evaluate
\begin{gather} G_{m,n}^\prime = G_{{\perp} m}^{\prime} G_{{\parallel} n}, \end{gather}
\begin{gather}G_{m,n}^{\prime\prime} = G_{{\perp} m}^{\prime\prime} G_{{\parallel} n}. \end{gather}D.1 Explicit expressions for some Bessel integrals
The 
$b(l)$-dependent factors in (3.17) can be written as
\begin{gather} G_{0} = G_{0,0}, \end{gather}
\begin{gather}G_{1} = G_{2,0} + G_{0,2}, \end{gather}
\begin{gather}G_{2} = G_{4,0} + 2 G_{2,2} + G_{0,4}, \end{gather}
\begin{gather}G_{3} = \tfrac{1}{2}G_{2,0} + G_{0,2}, \end{gather}
\begin{gather}G_{4} = \tfrac{1}{2}G_{4,0} + \tfrac{3}{2} G_{2,2} + G_{0,4}, \end{gather}
\begin{gather}G_{5} = \tfrac{1}{4}G_{4,0} + G_{2,2} + G_{0,4}, \end{gather}
which can be evaluated using the identities of the previous section in terms of the familiar 
$\varGamma _n(b)$ of gyrokinetic theory (suppressing its argument for succinctness)
\begin{gather} G_{0} = \varGamma_0, \end{gather}
\begin{gather}G_{1} = \left(\tfrac{3}{2}-b\right) \varGamma_0+b \varGamma_1, \end{gather}
\begin{gather}G_{2} = \tfrac{1}{4} \left(\left(6 b^2-20 b+15\right) \varGamma_0 + 2 b\left((10-4 b) \varGamma_1 + b \varGamma_2\right)\right), \end{gather}
\begin{gather}G_{3} = \tfrac{1}{2} \left(b \varGamma_1-(b-2) \varGamma_0\right), \end{gather}
\begin{gather}G_{4} = \tfrac{1}{4} \left(\left(3 b^2-11 b+10\right) \varGamma_0+b \left((11-4 b) \varGamma_1+b \varGamma_2\right)\right), \end{gather}
\begin{gather}G_{5} = \tfrac{1}{8} \left(\left(3 b^2-12 b+14\right) \varGamma_0+b \left(b \varGamma_2-4 (b-3) \varGamma_1\right)\right), \end{gather}where we recall
\begin{equation} \varGamma_n(b) = \exp({-}b){\rm I}_n(b). \end{equation}For completeness, we evaluate the few remaining factors that enter (3.17)
\begin{gather} G_{0,2} = \frac{\varGamma _0}{2}, \end{gather}
\begin{gather}G_{0,2}^{(1)} = \frac{\sqrt{b} \left(\varGamma _0-\varGamma _1\right)}{2 \sqrt{2}}, \end{gather}
\begin{gather}G_{0,2}^{(2)} = \tfrac{1}{8} \left(3 b \varGamma _0+(2-4 b) \varGamma _1+b \varGamma _2\right), \end{gather}and using (D11)
\begin{gather} G_{0,2}^\prime ={-}\sqrt{b/2} \frac{\partial}{\partial l}\left(\ln(b B)\right) \frac{\sqrt{b}\left(\varGamma_0 - \varGamma_1\right)}{2\sqrt{2}}, \end{gather}
\begin{gather}G_{0,2}^{\prime\prime} = b/2 \left(\frac{\partial}{\partial l}\ln(b B)\right)^2 {\left[3b\varGamma_0 + (2-4b)\varGamma_1 +b\varGamma_2 \right]}/8. \end{gather}D.2 Limit $b \rightarrow 0$
In the limit 
$b \rightarrow 0$ we obtain
\begin{gather} G_{0} = 1, \end{gather}
\begin{gather}G_{1} = \tfrac{3}{2}, \end{gather}
\begin{gather}G_{2} = \tfrac{15}{4}, \end{gather}
\begin{gather}G_{3} = 1, \end{gather}
\begin{gather}G_{4} = \tfrac{5}{2}, \end{gather}
\begin{gather}G_{5} = \tfrac{7}{4}, \end{gather}and
\begin{gather} G_{0,2} = \tfrac{1}{2}, \end{gather}
\begin{gather}G_{0,2}^{(1)} = 0, \end{gather}
\begin{gather}G_{0,2}^{(2)} = 0. \end{gather}Appendix E. Linear dispersion relation
For comparison, the following local linear dispersion relation can be solved numerically (here, we neglect finite Larmor radius effects 
$\textrm {J}_0 = 1$ and assume zero density gradient)
\begin{equation} 1 + \tau = \frac{2}{\sqrt{\rm \pi}}\int_0^{\infty} x_\perp \,{\rm d} x_\perp\int_{-\infty}^{\infty} {\rm d}x_\| \left[\frac{\varOmega - (x_\perp^2 + x_\|^2 - 3/2) }{\varOmega - \kappa_d^{{-}1}(x_\|^2 + x_\perp^2/2) - x_\| \kappa_\|^{{-}1}} \right]\exp({-}x_\|^2-x_\perp^2),\end{equation}
where 
$\varOmega = \omega /(\omega _* \eta )$. The velocity integrals can be evaluated separately for the 
$k_\| = 0$ and 
$\omega _d = 0$ cases in terms of the plasma dispersion function (Kadomtsev & Pogutse Reference Kadomtsev and Pogutse1970; Biglari et al. Reference Biglari, Diamond and Rosenbluth1989).
Open access
