Energetic bounds on gyrokinetic instabilities. Part 1. Fundamentals

Abstract

Upper bounds on the growth of free energy in gyrokinetics are derived. These bounds apply to all local gyrokinetic instabilities in the geometry of a flux tube, i.e. a slender volume of plasma aligned with the magnetic field, regardless of the geometry of field, the number of particle species or collisions. The results apply both to linear instabilities and to the nonlinear growth of finite-amplitude fluctuations.

1. Introduction

For the last six and a half decades, an enormous effort has been devoted to the study of microinstabilities in magnetically confined plasmas. Mathematically, such instabilities can be described by the Boltzmann equation for the plasma particles coupled to Maxwell's equations for the electric and magnetic fields, but it is often sufficient to consider the somewhat simpler gyrokinetic system of equations (Rutherford & Frieman 1968; Taylor & Hastie 1968; Catto 1978; Antonsen & Lane 1980; Catto, Tang & Baldwin 1981; Frieman & Chen 1982; Brizard & Hahm 2007; Krommes 2012; Catto 2019). These equations apply if the instability wavelength perpendicular to the magnetic field is comparable to the ion or electron gyroradius, but the wavelength is much longer in the direction along the field, which is normally the case for the most important microinstabilities and turbulence afflicting magnetised plasmas in the laboratory. Gyrokinetics also finds fruitful application in other parts of plasma physics, such as astrophysics (Schekochihin et al. 2009), and has been the subject of thousands of publications. Several millions of lines of computer code has been written for the purpose of numerically simulating gyrokinetic instabilities and turbulence (Kotschenreuther, Rewoldt & Tang 1995; Garbet et al. 2010).

As a result of this effort, a great deal of knowledge about various microinstabilities has been accumulated. Ion-temperature-gradient- (ITG-) and electron-temperature-gradient- driven modes, trapped-electron modes, kinetic ballooning modes and microtearing modes have, for instance, been found to be unstable and cause turbulence in tokamaks, stellarators and other fusion devices. However, a basic problem is that these and other instabilities tend to be sensitive to assumptions made about plasma parameters and the magnetic-field geometry. A cylindrical plasma does not have the same stability properties as a plasma slab, toroidal plasmas are different from cylindrical plasmas, and tokamaks and stellarators are also substantially different. As a result, little is known in general about gyrokinetic microinstabilities, despite the great effort devoted to their study.

In a recent publication (Helander & Plunk 2021), universal upper bounds on the growth rates of local gyrokinetic instabilities could nevertheless be derived in such a way that the results hold in any low-beta plasma, regardless of the magnetic geometry, number of particle species and collisions. The reason why these bounds are so general is they result from thermodynamic considerations. It is the budget of the Helmholtz free energy that constrains all instability growth rates to lie below the bounds in question. In the present paper, we provide more mathematical details of this calculation and extend it by showing how the bounds can be sharpened. In particular, we calculate the lowest possible bound on the growth rate that can be obtained from the free-energy budget of a plasma with ‘adiabatic’ electrons and a single kinetic ion species. In subsequent publications, such rates of ‘optimal growth’ will be derived in more complex cases that include both electrostatic and magnetic fluctuations. We will also show how the bounds can be lowered by simultaneously considering the budget of free energy and electrostatic energy, and compare them with gyrokinetic simulations. The present paper serves as an introduction to this series of publications.

2. Gyrokinetic system of equations

The mathematical setting of our considerations is that of local gyrokinetics. The distribution function of each species $a$ is written as (Catto 1978)

(2.1)\begin{equation} f_a({\boldsymbol r}, E_a, \mu_a, t) = F_{a0} ({\psi}, E_a) \left( 1 - \frac{e_a \delta \phi({\boldsymbol r})}{T_a} \right) + g_a({\boldsymbol R}, E_a, \mu_a, t), \end{equation}

where ${\boldsymbol r}$ denotes the particle position and ${\boldsymbol R} = {\boldsymbol r} - {\boldsymbol b} \times {\boldsymbol v} / \varOmega _a$ the gyrocentre position. Here, the magnetic field has been written as ${\boldsymbol B} = B {\boldsymbol b} = \boldsymbol {\nabla } \psi \times \boldsymbol {\nabla } \alpha$ in terms of Clebsch coordinates $(\psi,\alpha )$. If the magnetic field lines trace out toroidal surfaces, as in tokamaks and stellarators, a ballooning transform is necessary unless all field lines close on themselves. The gyrofrequency is $\varOmega _a = e_a B / m_a$, where $m_a$ denotes mass and $e_a$ charge. The equilibrium distribution function is taken to be Maxwellian, with density $n_a(\psi )$ and temperature $T_a(\psi )$ constant on magnetic surfaces, and no mean flow velocity. The particle velocity is denoted ${\boldsymbol v} = v_\| {\boldsymbol b} + {\boldsymbol v}_\perp$, the unperturbed energy by $E_a = m_a v^2 / 2 + e_a \varPhi (\psi )$, and the magnetic moment $\mu _a = m_a v_\perp ^2 / (2B)$ is a lowest-order constant of the motion. The geometry is taken to be that of a ‘flux tube’, i.e. a slender volume of plasma aligned with the magnetic field, with a rectangular cross-section in the $(\psi,\alpha )$-plane. Periodic boundary conditions on the fluctuations will be applied in this plane, so that all perturbations can be Fourier decomposed. For instance, the electrostatic potential fluctuations $\delta \phi$ are

(2.2)\begin{equation} \delta \phi(\psi,\alpha,l) = \sum_{\boldsymbol k} \delta \phi_{\boldsymbol k} (l) \exp [{{\rm i}(k_\psi \psi + k_\alpha \alpha)}], \end{equation}

where ${\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. The Fourier coefficients must satisfy $\delta \phi ^\ast _{\boldsymbol k} = \delta \phi _{-\boldsymbol k}$ in order that the potential be real.

The ‘non-adiabatic’ part of the distribution function $g_a$ evolves according to the nonlinear gyrokinetic equation (Frieman & Chen 1982)

(2.3)\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}') \bar{\chi}_{a, {\boldsymbol k}'} g_{a, {\boldsymbol k} - {\boldsymbol k}'}\nonumber\\ = \sum_b \left[ C_{ab}(g_{a, \boldsymbol k},F_{b0}) + C_{ab}(F_{a0},g_{b, \boldsymbol k}) \right] + \frac{e_a F_{a0}}{T_a} \left( \frac{\partial}{\partial t} + {\rm i} \omega_{{\ast} a}^T \right) \bar{\chi}_{a, \boldsymbol k} , \end{gather}

where $\omega _d = {\boldsymbol k} \boldsymbol {\cdot } {\boldsymbol v}_d$ denotes the drift frequency (with ${\boldsymbol v}_d$ being the unperturbed drift velocity),

(2.4)\begin{gather} \omega_{{\ast} a} = \frac{k_\alpha T_a}{e_a} \frac{{\rm d} \ln n_a}{{\rm d} \psi}, \end{gather}

(2.5)\begin{gather} \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{gather}

(2.6)\begin{gather} \bar{\chi}_{a \boldsymbol k} = {\rm J}_0 \left( \frac{k_\perp v_\perp}{\varOmega_a} \right)\left( \delta \phi_{{\boldsymbol k}} - v_\| \delta A_{\| {\boldsymbol k}} \right) + {\rm J}_1 \left( \frac{k_\perp v_\perp}{\varOmega_a} \right) \frac{v_\perp}{k_\perp} \delta B_{\| \boldsymbol k}, \end{gather}

and ${\rm J}_0$ and ${\rm J}_1$ are Bessel functions. The gyro-averaged and linearised collision operator between species $a$ and $b$ is denoted by $C_{ab}$, and the field perturbations are given by

(2.7)\begin{gather} \sum_a \lambda_a \delta \phi_{\boldsymbol k} = \sum_a e_a \int g_{a, {\boldsymbol k}} {\rm J}_{0a}\, {\rm d}^3v, \end{gather}

(2.8)\begin{gather} \delta A_{\| {\boldsymbol k}} = \frac{\mu_0}{k_\perp^2} \sum_a e_a \int v_\| g_{a, {\boldsymbol k}} {\rm J}_{0a}\,{\rm d}^3v, \end{gather}

(2.9)\begin{gather} \delta B_{\| {\boldsymbol k}} ={-} \frac{\mu_0}{k_\perp} \sum_a e_a \int v_\perp g_{a, {\boldsymbol k}} {\rm J}_{1a} \,{\rm d}^3v. \end{gather}

Here and in the following, we write $\lambda _a = {n_a e_a^2}/{T_a}$ and ${\rm J}_{na} = {\rm J}_n(k_\perp v_\perp / \varOmega _a)$. Equation (2.7) expresses quasineutrality, (2.8) Ampère's law and (2.9) the condition that the sum of the thermal pressure and the magnetic pressure should be constant on the short length scale of the fluctuations. The volume element in velocity space is

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

where the sum is taken over both values of $\sigma = v_\| / | v_\| | = \pm 1$.

Note that we restrict our attention to the original gyrokinetic equation (2.3) of Frieman & Chen (1982), which does not include equilibrium flows. We thus only consider instabilities caused by density and temperature gradients, but not those associated with velocity-space anisotropy or non-Maxwellian distribution functions, such as fast-ion-driven instabilities (Chen & Zonca 2016). Moreover, stabilisation or destabilisation associated with flow-velocity shear is not included in the analysis although it can be quite important in practice (see, e.g., Barnes et al. 2011). It should be possible to include such effects by adding an appropriate term to (2.3), at least in the case that the equilibrium is axisymmetric (Artun & Tang 1994; Parra, Barnes & Peeters 2011). In non-symmetric equilibria, the situation is fundamentally more complicated because any equilibrium flow must be small (Helander 2014).

As we show in the following, it is advantageous to introduce the function

(2.11)\begin{equation} \delta F_{a, \boldsymbol k} = g_{a, \boldsymbol k} - \frac{e_a {\rm J}_{0a} \delta \phi_{\boldsymbol k}}{T_a} F_{a0}, \end{equation}

where all quantities are evaluated at the gyrocentre position $\boldsymbol R$. The quasineutrality condition then becomes

(2.12)\begin{equation} \sum_a \lambda_a \left[1 - \varGamma_0(b_a) \right] \delta \phi_{\boldsymbol k} = \sum_a e_a \int \delta F_{a, {\boldsymbol k}} {\rm J}_{0a} \,{\rm d}^3v, \end{equation}

where $\varGamma _0(x) = {\rm I}_0(x) {\rm e}^{-x}$, $b_a = k_\perp ^2 \rho _a^2 = k_\perp ^2 T_a / (m_a \varOmega _a^2)$ and we have used an integral given in appendix A. In the following, we sometimes write $\varGamma _{0a}$ instead of $\varGamma _0(b_a)$.

3. Helmholtz free energy

The budget of Helmholtz free energy has been considered by several authors, e.g. Krommes & Hu (1993), Brizard (1994), Sugama et al. (1996), Garbet et al. (2005), Schekochihin et al. (2009), Banon Navarro et al. (2011), Hatch et al. (2016) and Stoltzfus-Dueck & Scott (2017), and is obtained by multiplying the gyrokinetic equation (2.3) by $T_a g_a^\ast / F_{a0}$, taking the real part, summing over all species and wavenumbers, integrating over velocity space and, finally, taking an average over the volume of the flux tube, which we denote by angular brackets,

(3.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}

We note that the average could also be defined keeping $L$ finite, e.g. for periodic systems, without affecting what follows. In order for the integral to converge, we require that the functions $\bar {\chi }_{\boldsymbol k}(l)$ should be bounded. On the left-hand side of (2.3), this operation,

(3.2)\begin{equation} {\rm Re} \, \sum_{a,{\boldsymbol k}} T_a \left\langle \int \left( \cdots \right) \frac{g_{a, \boldsymbol k}^\ast}{F_{a0}} \,{\rm d}^3v \right\rangle, \end{equation}

annihilates the second term because

(3.3)\begin{equation} {\rm Re} \, \left\langle \int v_\| \frac{g_{a, {\boldsymbol k}}^\ast}{F_{a0}} \frac{\partial g_{a, {\boldsymbol k}}}{\partial l} \,{\rm d}^3v \right\rangle = \left.\lim_{L\rightarrow \infty} \sum_\sigma \frac{{\rm \pi} \sigma}{m_a^2} \int_{{-}L}^L \,{\rm d}l \int \frac{{\rm d}E_a}{F_{a0} } \int \frac{\partial |g_{a, {\boldsymbol k}}|^2}{\partial l} \,{\rm d}\mu_a \right/ \int_{{-}L}^L \frac{{\rm d}l}{B} = 0, \end{equation}

where we have used (2.10) and assumed that $|g_{a, {\boldsymbol k}}|^2$ remains bounded as $l \rightarrow \infty$.¹ The operation also eliminates the third term because $\omega _{da}$ is real, and the fourth term because

(3.4)\begin{equation} {\rm Re} ({\boldsymbol k} \times {\boldsymbol k}') g^\ast_{a,\boldsymbol k} \bar{\chi}_{a, {\boldsymbol k}'} g_{a,{\boldsymbol k} - {\boldsymbol k}'} = {\rm Re} ({\boldsymbol k} \times {\boldsymbol q}) g_{a,-\boldsymbol k} \bar{\chi}_{a, {\boldsymbol k} + {\boldsymbol q}} g_{a, - {\boldsymbol q}}, \end{equation}

where ${\boldsymbol q} = {\boldsymbol k}' - {\boldsymbol k}$ and we have used $g^\ast _{a,\boldsymbol k} = g_{a,-\boldsymbol k}$. Because the right-hand side changes sign if ${\boldsymbol k}$ and $\boldsymbol q$ are interchanged, the result vanishes upon summation over $\boldsymbol k$ and ${\boldsymbol q}$. The remainder of the equation thus becomes

(3.5)\begin{equation} \frac{{\rm d}}{{\rm d}t} \sum_{a,{\boldsymbol k}} T_a \left\langle \int \frac{|g_{a \boldsymbol k}|^2}{2 F_{a0}} \,{\rm d}^3v \right\rangle = \sum_{\boldsymbol k} C({\boldsymbol k},t) + {\rm Re} \sum_{a,{\boldsymbol k}} e_a\left\langle \int g_{a,{\boldsymbol k}}^\ast \left( \frac{\partial}{\partial t} + {\rm i} \omega_{{\ast} a}^T \right) \bar{\chi}_{a \boldsymbol k} \,{\rm d}^3v \right\rangle, \end{equation}

where

(3.6)\begin{equation} C({\boldsymbol k},t) = {\rm Re} \, \sum_{a,b} T_a \left\langle \int \frac{g_{a, \boldsymbol k}^\ast}{ F_{a0}} \left[C_{ab}(g_{a, \boldsymbol k},F_{b0}) + C_{ab}(F_{a0},g_{b, \boldsymbol k}) \right] \,{\rm d}^3v \right\rangle \le 0 \end{equation}

is negative or vanishes by Boltzmann's $H$-theorem. By using the field equations (2.7)–(2.9), we find

(3.7)\begin{equation} \sum_{a}e_a \int g_{a,{\boldsymbol k}}^\ast \frac{\partial \bar{\chi}_{a \boldsymbol k}}{\partial t} \,{\rm d}^3v = \frac{1}{2} \frac{{\rm d}}{{\rm d}t} \left( \sum_{a} \lambda_a |\delta \phi_{\boldsymbol k}|^2 - \frac{| \delta {\boldsymbol B}_{\boldsymbol k} |^2}{\mu_0}\right), \end{equation}

where $| \delta {\boldsymbol B}_{\boldsymbol k} |^2 = | k_\perp \delta A_{\|\boldsymbol k} |^2 + | \delta B_{\|\boldsymbol k} |^2$ and, thus, we obtain our key equation:

(3.8)\begin{equation} \frac{{\rm d}}{{\rm d}t} \sum_{\boldsymbol k} H({\boldsymbol k},t) = 2 \sum_{\boldsymbol k} \left[C({\boldsymbol k},t) + D({\boldsymbol k},t)) \right], \end{equation}

where we have written

(3.9)\begin{gather} D({\boldsymbol k}, t) = {\rm Im} \, \sum_a e_a \left\langle \int g_{a, \boldsymbol k} \omega_{{\ast} a}^T \bar{\chi}^\ast_{a,\boldsymbol k} \,{\rm d}^3v \right\rangle, \end{gather}

(3.10)\begin{gather} H({\boldsymbol k},t) = \sum_a \left\langle T_a \int \frac{|g_{a, \boldsymbol k}|^2}{F_{a0}} \,{\rm d}^3v - \lambda_a |\delta \phi_{\boldsymbol k}|^2 \right\rangle + \left\langle \frac{| \delta {\boldsymbol B}_{\boldsymbol k} |^2}{\mu_0} \right\rangle. \end{gather}

It is helpful to write $H$ in terms of $\delta F_a$, defined in (2.11), instead of $g_a$:

(3.11)\begin{equation} H({\boldsymbol k},t) = \sum_a \left\langle T_a \int \frac{|\delta F_{a, \boldsymbol k}|^2}{F_{a0}} \,{\rm d}^3v + \lambda_a (1 - \varGamma_{0a} ) | \delta \phi_{\boldsymbol k}|^2 \right\rangle + \left\langle \frac{| \delta {\boldsymbol B}_{\boldsymbol k} |^2}{\mu_0} \right\rangle, \end{equation}

which makes it clear that $H$ can never be negative and only vanishes if all distribution-function perturbations $\delta F_a$ vanish everywhere in phase space. The first term in $H$ is recognised from the Gibbs entropy formula: if $F = F_0 + \delta F$, then to second order in $\delta F$,

(3.12)\begin{equation} - \int F \ln F \,{\rm d}^3v ={-} \int \left[ F_0 \ln F_0 + \left( 1 + \ln F_0 \right) \delta F + \frac{\delta F^2}{2 F_0} \right] \,{\rm d}^3v, \end{equation}

which motivates us to define

(3.13)\begin{equation} S_a({\boldsymbol k},t) ={-} \left\langle \int \frac{|\delta F_{a, \boldsymbol k}|^2}{F_{a0}} \,{\rm d}^3v \right\rangle. \end{equation}

Furthermore, we write

(3.14)\begin{equation} U({\boldsymbol k},t) = \left\langle \sum_{a} \lambda_a (1 - \varGamma_{0a} ) | \delta \phi_{\boldsymbol k}|^2 + \frac{| \delta {\boldsymbol B} |^2}{\mu_0} \right\rangle, \end{equation}

and note that, in the short-wavelength limit, $b_a = (k_\perp \rho _a)^2 \ll 1$, $\varGamma _0(b_a) = 1 - b_a + O(b_a^2)$, so that

(3.15)\begin{equation} U({\boldsymbol k},t) = \left\langle \sum_{a} \frac{m_a n_a k^2 | \delta \phi_{\boldsymbol k}|^2}{B^2} + \frac{| \delta {\boldsymbol B} |^2}{\mu_0} \right\rangle, \end{equation}

where the first term represents the kinetic energy of ${\boldsymbol E} \times {\boldsymbol B}$ motion and the second term magnetic energy. We thus arrive at the formula

(3.16)\begin{equation} H({\boldsymbol k},t) = U ({\boldsymbol k},t) - \sum_a T_a S_a({\boldsymbol k},t), \end{equation}

with $U$ denoting the energy of the fluctuations and $S_a$ their entropy, suggesting that $H$ describes the Helmholtz free energy of the fluctuations and (3.8) the budget of this energy. Indeed, on the right-hand side of this equation $C$ reflects the increase in entropy due to collisions, and $D$ can be written as

(3.17)\begin{align} D({\boldsymbol k},t) & = {\rm Re} \, \sum_a T_a \left\langle \int g_a \delta \dot {\boldsymbol R}_{a, \boldsymbol k}^\ast \boldsymbol{\cdot} \boldsymbol{\nabla} F_{a0} \,{\rm d}^3v \right\rangle\nonumber\\ & ={-} \sum_a \left ( T_a\varGamma_a \frac{{\rm d} \ln p_a}{{\rm d} \psi} + q_a \frac{{\rm d} \ln T_a}{{\rm d} \psi} \right). \end{align}

Here

(3.18)\begin{equation} \delta \dot {\boldsymbol R}_{a, \boldsymbol k} = \frac{{\rm i} \bar{\chi}_{a, \boldsymbol k}{\boldsymbol b} \times {\boldsymbol k}}{B} \end{equation}

describes the gyrocentre velocity perturbation due to the fluctuations, and the radial particle and heat fluxes are

(3.19)\begin{equation} \left. \begin{gathered} \varGamma_a({\boldsymbol k},t) = {\rm Re}\left\langle \int \delta F_{a, \boldsymbol k} (\delta \dot {\boldsymbol R}_{a, \boldsymbol k}^\ast \boldsymbol{\cdot} \boldsymbol{\nabla} \psi)\,{\rm d}^3v \right\rangle,\\ q_a({\boldsymbol k},t) = {\rm Re}\left\langle \int \delta F_{a, \boldsymbol k} \left( \frac{m_a v^2}{2} - \frac{5 T_a}{2} \right) (\delta \dot {\boldsymbol R}_{a, \boldsymbol k}^\ast \boldsymbol{\cdot} \boldsymbol{\nabla} \psi)\,{\rm d}^3v \right\rangle. \end{gathered} \right\} \end{equation}

The term in (3.17) involving $\varGamma _a$ is thus suggestive of the thermodynamic work performed by the particle flux against the pressure gradient, and the term involving $q_a$ relates to entropy production due to a heat flux down the temperature gradient.

Thanks to the nonlinear term in the gyrokinetic equation, free energy can be transferred between different wavenumbers and be ‘cascaded’ to small scales, where it is dissipated by collisions, much like kinetic energy in Navier–Stokes turbulence. The way in which this occurs and gives rise to a turbulent spectrum of fluctuations has been studied extensively in the literature (Schekochihin et al. 2009; Tatsuno et al. 2009; Banon Navarro et al. 2011; Stoltzfus-Dueck & Scott 2017). We shall use the free-energy budget (3.8) for a different purpose, namely, to derive rigorous upper bounds on linear and nonlinear growth rates. Outside the realm of gyrokinetics, this has earlier been accomplished for linear instabilities by Fowler and co-workers (Fowler 1964,  1968; Brizard et al. 1991).

4. Cauchy–Schwarz inequalities

For simplicity, we restrict our considerations to low-beta plasmas, where fluctuations in the magnetic-field strength can be neglected, $\delta B_\| = 0$. This approximation is common in the literature but will be removed in the next publication in this series of papers.

Our basic mathematical tools are the triangle and Cauchy–Schwarz inequalities, which limit the amplitude of field fluctuations that are possible given a certain entropy budget. For instance, it follows from the field equation (2.12) that the electrostatic potential is bounded by

(4.1)\begin{equation} \sum_a \lambda_a \left(1 - \varGamma_{0a} \right) |\delta \phi_{\boldsymbol k} | \le \sum_a | e_a | \left( \int \frac{ |\delta F_{a,{\boldsymbol k}} |^2}{F_{a0}} \,{\rm d}^3v \int F_{a0} {\rm J}_{0a}^2 \,{\rm d}^3v \right)^{1/2}. \end{equation}

Thus, if we measure the relative entropy perturbation at the scale $\boldsymbol k$ of each species $a$ by the dimensionless quantity

(4.2)\begin{equation} s_a({\boldsymbol k},t) = \frac{1}{n_a} \int \frac{|\delta F_{a \boldsymbol k} |^2}{F_{a0}} \,{\rm d}^3v, \end{equation}

then it follows that the electrostatic potential is subject to the bound

(4.3)\begin{equation} \sum_a \lambda_a \left(1 - \varGamma_{0a} \right) |\delta \phi_{\boldsymbol k} | \le \sum_a n_a |e_a| \sqrt{ \varGamma_{0a} s_a}. \end{equation}

Analogously, it follows from Ampère's law (2.8) that the magnetic potential is limited by

(4.4)\begin{equation} | \delta A_{\| \boldsymbol k} | \le \sum_a \frac{\mu_0 | e_a | }{k_\perp^2}\left( \int \frac{ |\delta F_{a,{\boldsymbol k}} |^2}{F_{a0}} \,{\rm d}^3v \int v_\|^2 F_{a0} {\rm J}_{0a}^2 \,{\rm d}^3v \right)^{1/2}, \end{equation}

i.e.

(4.5)\begin{equation} \frac{k_\perp | \delta A_{\| \boldsymbol k} |}{B} \le \sum_a \frac{\beta_a}{2 k_\perp \rho_a} \sqrt{\varGamma_{0a} s_a} \simeq \frac{\beta_e}{2 k_\perp \rho_e} \sqrt{\varGamma_{0e} s_e}, \end{equation}

where $\beta _a(l) = 2 \mu _0 n_a T_a / B^2$. In the last, approximate equality, we have recognised the fact that the sum is usually dominated by the contribution from the electrons thanks to their small gyroradius. Because $k_\perp \rho _e$ is small for most instabilities of interest, the inequality (4.5) is not very restrictive, but it is nevertheless valuable as it implies that gyrokinetic instabilities are electrostatic in the limit $\beta _e \rightarrow 0$. Moreover, it is used in the following to demonstrate that the growth rate remains bounded in the limit $k_\perp \rho _e \rightarrow \infty$.

We can also apply the triangle and Cauchy–Schwarz inequalities to the free-energy production rate (3.9):

(4.6)\begin{gather} D({\boldsymbol k},t) \le \sum_a |e_a| |n_a s_a|^{1/2} \left\langle \int F_{a0} (\omega_{{\ast} a}^T)^2 {\rm J}_0^2 \left(|\delta \phi_{\boldsymbol k} |^2 + v_\|^2 |\delta A_{\| \boldsymbol k}|^2 \right) \,{\rm d}^3v \right\rangle^{1/2}\nonumber\\ = \sum_a n_a |e_a \omega_{{\ast} a}| |s_a|^{1/2} \left\langle M(\eta_a, b_a) |\delta \phi_{\boldsymbol k} |^2 + N(\eta_a, b_a) \frac{T_a |\delta A_{\| \boldsymbol k}|^2}{m_a} \right\rangle^{1/2}, \end{gather}

where the functions

(4.7)\begin{gather} M(\eta_a ,b_a) = \frac{1}{n_a} \int \left[1 + \eta_a \left( \frac{m_a v^2}{2 T_a} - \frac{3}{2} \right) \right]^2 F_{a0} {\rm J}_{0a}^2 \,{\rm d}^3v, \end{gather}

(4.8)\begin{gather} N(\eta_a ,b_a) = \frac{1}{n_a} \int \frac{m_a v_\|^2}{T_a} \left[1 + \eta_a \left( \frac{m_a v^2}{2 T_a} - \frac{3}{2} \right) \right]^2 F_{a0} {\rm J}_{0a}^2 \,{\rm d}^3v \end{gather}

can be expressed in terms of modified Bessel functions as

(4.9)\begin{gather} M(\eta,b) = \left( 1 + \frac{3 \eta^2}{2} - 2\eta(1+ \eta) b + 2 \eta^2 b^2 \right) \varGamma_0(b) + \eta b\left( 2 + \eta - 2 \eta b \right) \varGamma_1(b), \end{gather}

(4.10)\begin{gather} N(\eta,b) = \left( 1 + 2 \eta + \frac{7 \eta^2}{2} - 2\eta(1+ 2\eta) b + 2 \eta^2 b^2 \right) \varGamma_0(b) + \eta b\left( 2 + 3 \eta - 2 \eta b \right) \varGamma_1(b), \end{gather}

using integrals given in appendix A. In the limits of very small and very large wavelength, respectively, the asymptotic forms of these functions are

(4.11)\begin{gather} M(\eta,b) \simeq \left\{ \begin{array}{@{}ll} 1 + \dfrac{3 \eta^2}{2}, & b \rightarrow 0 \\ \dfrac{1-\eta + \dfrac{5 \eta^2}{4}}{\sqrt{2 {\rm \pi}b}}, & b \rightarrow \infty, \end{array} \right. \end{gather}

(4.12)\begin{gather}N(\eta,b) \simeq \left\{ \begin{array}{@{}ll} 1 + 2 \eta + \dfrac{7 \eta^2}{2}, & b \rightarrow 0 \\ \dfrac{1+\eta + \dfrac{9 \eta^2}{4}}{\sqrt{2 {\rm \pi}b}}, & b \rightarrow \infty. \end{array} \right. \end{gather}

5. Upper bounds on linear growth rates

In this section, we temporarily consider linear instabilities and thus focus on a single pair of wavenumbers $(k_\psi,k_\alpha )$. Thanks to Boltzmann's $H$-theorem, the quantity $C({\boldsymbol k},t)$ is always negative and the relation (3.8) thus implies an upper bound on the linear growth rate

(5.1)\begin{equation} \gamma({\boldsymbol k}) \le \frac{D({\boldsymbol k},t)}{H({\boldsymbol k},t)}. \end{equation}

As we have already bounded $D$ from above, we merely need to find a suitable bound on

(5.2)\begin{equation} H({\boldsymbol k},t) = \sum_a \left\langle n_a T_a s_a + \lambda_a \left( 1 - \varGamma_{0a} \right) |\delta \phi_{\boldsymbol k}|^2 \right\rangle + \left\langle \frac{| k_\perp \delta A_{\| \boldsymbol k} |^2}{\mu_0} \right\rangle \end{equation}

from below to derive an upper bound on $\gamma ({\boldsymbol k})$. Some care is needed to construct reasonably tight bounds, but all results are largely independent of the geometry of the magnetic field because the second and third terms from (2.3) do not contribute to the free-energy balance equation (3.8). The bound (5.1) therefore only depends on the magnetic geometry through the two quantities $B(l)$ and $k_\perp (l) = | k_\psi \boldsymbol {\nabla } \psi + k_\alpha \boldsymbol {\nabla } \alpha |$.

5.1. Adiabatic electrons

We begin by considering the simplest case of a hydrogen plasma with a Boltzmann-distributed, or so-called ‘adiabatic’, electron response, where $g_e$ is taken to vanish. This is the traditionally simplest gyrokinetic model of ITG and trapped-ion instabilities, which account for a substantial fraction of the turbulence and transport in tokamaks and stellarators, and therefore has been the subject of hundreds, if not thousands, of publications. As $g_e$ vanishes and there are no magnetic fluctuations [in the approximation used in (4.5)], the free energy becomes

(5.3)\begin{equation} {H} = {n T_i}\left\langle s_i + \left(1 +\tau - \varGamma_{0i} \right) \left| \frac{e \delta \phi_{\boldsymbol k}}{T_i} \right|^2 \right\rangle. \end{equation}

where $n = n_i = n_e$ and $\tau = T_i / T_e$. Furthermore, the quasineutrality condition (2.12) reduces to

(5.4)\begin{equation} \left(1 +\tau - \varGamma_{0i} \right) \frac{e \delta \phi_{\boldsymbol k}}{T_i} = \frac{1}{n} \int \delta F_i {\rm J}_{0i} \,{\rm d}^3v, \end{equation}

and the bound (4.3) is thus replaced by the more stringent condition

(5.5)\begin{equation} \left(1 +\tau - \varGamma_{0i} \right) \frac{e |\delta \phi_{\boldsymbol k}|}{T_i} \le \sqrt{ \varGamma_{0i} s_{i}}. \end{equation}

Thanks to this inequality, the free energy satisfies

(5.6)\begin{equation} H \ge \left\langle \frac{1+\tau}{\varGamma_{0i}} \left(1 +\tau - \varGamma_{0i} \right) \left| \frac{e \delta \phi_{\boldsymbol k}}{T_i} \right|^2 \right\rangle. \end{equation}

The free-energy production term can be simplified somewhat because the quasineutrality condition (5.4) in the case of adiabatic electrons implies that there is no particle flux. Indeed, the flux from (3.19),

(5.7)\begin{equation} \varGamma_i({\boldsymbol k},t) ={-} {\rm Re}\left\langle \frac{{\rm i} \delta {\phi}_{\boldsymbol k}^\ast ({\boldsymbol b} \times {\boldsymbol k}) \boldsymbol{\cdot} \boldsymbol{\nabla} \psi}{B} \int \delta F_{i, \boldsymbol k} {\rm J}_{0i} \,{\rm d}^3v \right\rangle, \end{equation}

vanishes because of (2.12), and $D$ thus becomes

(5.8)\begin{equation} D({\boldsymbol k},t) = {\rm Im} \, \eta_i \omega_{{\ast} i} \left\langle e \delta \phi_{\boldsymbol k}^\ast\int g_{i\boldsymbol k} \left( \frac{m_i v^2}{2 T_i} - \frac{3}{2} \right) {\rm J}_{0i} \,{\rm d}^3v \right\rangle. \end{equation}

As a result, in the inequality (4.6), the function $M(\eta,b)$ can be replaced by

(5.9)\begin{equation} \tilde M(\eta,b) = \eta^2 \left[ \left(\frac{3 }{2} - 2b + 2 b^2 \right) \varGamma_0(b) + b\left(1 - 2 b \right) \varGamma_1(b) \right], \end{equation}

and the bound (5.1) becomes

(5.10)\begin{equation} \frac{\gamma}{\omega_{{\ast} i}} \le \frac{\left\langle \tilde M(\eta_i, b_i) |\delta \phi_{\boldsymbol k} |^2 \right\rangle^{1/2}}{ \left\langle (1+\tau) [(1+\tau) \varGamma_{0i}^{{-}1} - 1]|\delta \phi_{\boldsymbol k} |^2 \right\rangle^{1/2}}. \end{equation}

Here, $\tilde M (\eta _i, b_i)$ is a decreasing function of $b_i$, and the denominator is an increasing function of the same quantity. The right-hand side is thus maximised by choosing $|\delta \phi _{\boldsymbol k}(l)|^2 = \delta (l-l_0)$, where $l_0$ is the position along the field line where $b_i(l) = k_\perp ^2 \rho _i^2 \propto (k_\perp /B)^2$ is minimised. We thus obtain

(5.11)\begin{equation} \frac{\gamma}{\omega_{{\ast} i}} \le \sqrt{ \frac{ \tilde M(\eta_i, b_{\rm min})} {(1+\tau) \left[ (1 + \tau) \varGamma_0^{{-}1} (b_{\rm min}) - 1 \right]}}, \end{equation}

where $b_\textrm {min} = b_i(l_0)$. The result is plotted in figure 1. Note that all dependence on the geometry of the magnetic field has disappeared: our limit on the growth rate is spatially local in nature and only depends on the minimum value of $k_\perp \rho _i$.

Figure 1. (a) Upper bound (5.11) on the growth rate normalised to $\eta _i \omega _{\ast i} / (k_\perp \rho _i)$ of gyrokinetic instabilities for $k_\psi = 0$ and three different values of $\tau = T_i/T_e$ in a hydrogen plasma with adiabatic electrons as a function of the smallest value of $k_\perp \rho _i$ along the magnetic field. (b) The best possible bound (6.20) for free-energy growth, which is about a factor of two lower.

This bound, which applies to all local gyrokinetic instabilities in a plasma with adiabatic electrons, is not optimal and can be improved by a factor of approximately two, as we show in the next section. Nevertheless, it displays scalings that have been seen in many publications and numerical simulations over the years. For long wavelengths, $b_i \rightarrow 0$, it reduces to

(5.12)\begin{equation} \gamma \le |\eta_i \omega_{{\ast} i}| \sqrt{ \frac{3}{2\tau(1+\tau)}}. \end{equation}

Note that all dependence on the magnetic geometry has disappeared, and because $\omega _{\ast i} \propto k_\alpha$ the growth rate is proportional to $k_\alpha$ in this limit. For short wavelengths, $k_\perp \rho _i \gg 1$, the bound remains finite,

(5.13)\begin{equation} \gamma \le \frac{| \eta_i \omega_{{\ast} i} |}{1+\tau} \sqrt{ \frac{ 5}{8 {\rm \pi}b_{\rm min}} }, \end{equation}

because

(5.14)\begin{equation} b_{\rm min} = \min_l \left[ \left( k_\psi^2 |\boldsymbol{\nabla} \psi |^2 + 2 k_\psi k_\alpha \boldsymbol{\nabla} \psi \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha + k_\alpha^2 |\boldsymbol{\nabla} \alpha |^2 \right) \frac{T_i}{m_i \varOmega_i^2} \right] \end{equation}

is a positive-definite quadratic form in $k_\psi$ and $k_\alpha$. Indeed, $\gamma (k_\psi,k_\alpha )$ approaches a finite constant in the limit $k_\alpha \rightarrow \infty$ and vanishes if $k_\psi \rightarrow \infty$ at fixed $k_\alpha$. Moreover, at constant ion temperature, the bound (5.11) increases with the electron temperature through the scaling with $\tau$, which is a well-known feature of numerical simulations and analytical dispersion relations in explicitly tractable limits (Biglari, Diamond & Rosenbluth 1989; Romanelli 1989; Plunk et al. 2014; Zocco et al. 2018). This unfortunate scaling is thought to degrade energy confinement in electron-heated tokamaks and stellarators.

5.2. Electromagnetic instabilities

We now turn to the more general case of an arbitrary number of kinetic species, but still restrict our attention to instabilities with $\delta B_\| = 0$. No attempt will be made to make the bound as low as possible. Our main concern is to show that an upper bound exists and that it is itself bounded as a function of $\boldsymbol k$, so that there is a universal upper bound on the growth rate at any wavelength. This result will be of crucial importance when we consider nonlinear growth in a subsequent section. In the next publication of this series, we show how to extend the calculation to include fluctuations of the magnetic field strength and how to compute the lowest possible bounds in this context.

We begin by seeking lower bounds on $H$ under the constraints (4.3) and (4.5), which lead us to a simple quadratic minimisation problem treated in appendix B, where the minimum

(5.15)\begin{equation} \min_{x_1, x_2, \ldots} f(x_1, x_2, \ldots) = \sum_a q_a x_a^2 \end{equation}

subject to the constraint

(5.16)\begin{equation} \sum_a p_a x_a \ge c, \end{equation}

is found for the case that $q_a$ and $p_a$ are positive real numbers. In terms of this notation, we first choose $x_a = \sqrt {s_a}$, $p_a = n_a |e_a| \sqrt {\varGamma _{0a}}$, $q_a = n_a T_a$ and

(5.17)\begin{equation} c = \sum_a \lambda_a \left( 1 - \varGamma_{0a} \right) |\delta \phi_{\boldsymbol k}|, \end{equation}

and then obtain

(5.18)\begin{equation} \left.\sum_a n_a T_a s_a \ge \left[ \sum_a \lambda_a \left( 1 - \varGamma_{0a} \right) |\delta \phi_{\boldsymbol k}|\right ]^2 \right/ \sum_c \lambda_c \varGamma_{0c}. \end{equation}

As a result of this inequality, we conclude from (5.2) that $H \ge \left \langle L |\delta \phi _{\boldsymbol k} |^2 \right \rangle$ with

(5.19)\begin{equation} L(l) = \left.\left( \sum_a \lambda_a \right) \left( \sum_b \lambda_b (1-\varGamma_{0b}) \right) \right/ \left( \sum_c \lambda_c \varGamma_{0c} \right). \end{equation}

Similarly, by instead choosing $c = |k_\perp \delta A_{\|\boldsymbol k} |/\mu _0$ and

(5.20)\begin{equation} p_a = \frac{n_a |e_a|}{k_\perp} \sqrt{ \frac{T_a \varGamma_{0a}}{m_a}}, \end{equation}

we find

(5.21)\begin{equation} \left. \sum_a n_a T_a s_a \ge \frac{|k_\perp \delta A_{\|\boldsymbol k} |^2}{\mu_0} \right/ \sum_a \frac{\beta_a \varGamma_{0a}}{2 b_a}, \end{equation}

where $\beta _a = 2 \mu _0 n_a T_a / B^2$. Because the gyroradius of the electrons is usually much smaller than that of any ion species and $\varGamma _{a0} = \varGamma _0(b_a)$ is a decreasing function of particle mass, only the electrons need to be kept in the sum over species, and we conclude that $H$ is bounded from below by

(5.22)\begin{equation} H({\boldsymbol k},t) \ge \left\langle \frac{|k_\perp \delta A_{\|\boldsymbol k} |^2}{\mu_0} \left( 1 + \frac{2 b_e}{\beta_e \varGamma_{0e}} \right) \right\rangle = \frac{n{\rm e}^2}{m_e} \left\langle K |\delta A_{\|\boldsymbol k} |^2 \right\rangle, \end{equation}

with

(5.23)\begin{equation} K(l) = \frac{2 b_e}{\beta_e} \left( 1 + \frac{2 b_e}{\beta_e \varGamma_{0e}} \right). \end{equation}

We are now ready to apply our basic upper bound (5.1), where we use (4.6) and

(5.24)\begin{gather} H \ge \left\langle n_a T_a s_a \right\rangle^{1/2} \left\langle L |\delta \phi_{\boldsymbol k} |^2 \right\rangle^{1/2}, \end{gather}

(5.25)\begin{gather} H \ge \left\langle n_a T_a s_a \right\rangle^{1/2} \left\langle \frac{n{\rm e}^2}{m_e} K |\delta A_{\|\boldsymbol k} |^2 \right\rangle^{1/2}, \end{gather}

to conclude that

(5.26)\begin{equation} \gamma \le \sum_a |\omega_{{\ast} a}| \sqrt{\frac{\left\langle \lambda_a M(\eta_a, b_a) |\delta \phi_{\boldsymbol k}|^2 \right\rangle}{\left\langle L |\delta \phi_{\boldsymbol k}|^2 \right\rangle}} + |\omega_{{\ast} e}| \sqrt{\frac{\left\langle N(\eta_e, b_e) |\delta A_{\| \boldsymbol k} |^2 \right\rangle}{\left\langle K |\delta A_{\| \boldsymbol k} |^2\right\rangle}} \end{equation}

where the contribution from ions to the electromagnetic term in $D$ has been neglected, being a factor of order $m_e/m_i$ smaller than the electron contribution. As $L$ is an increasing function of the quantities $b_a$, which are all proportional to $(k_\perp / B)^2$, the first term on the right is maximised if $|\delta \phi _{\boldsymbol k}(l)|^2$ is chosen to be a delta function in the point $l_0$ where the function $k_\perp (l) / B(l)$ attains its minimum. Similarly, the second term is maximised by choosing $| \delta A_{\| \boldsymbol k}(l) |^2 \propto \delta (l-l_1)$ where $l_1$ is the point where $K(l)/N(l)$ is minimised. We thus arrive at the result

(5.27)\begin{equation} \gamma({\boldsymbol k}) \le \gamma_{\rm bound} ({\boldsymbol k}) = \sum_a |\omega_{{\ast} a}| \sqrt{\frac{\lambda_a M(\eta_a, b_a(l_0))}{L(l_0)}} + |\omega_{{\ast} e}| \sqrt{\frac{ N(\eta_e, b_e(l_1))}{K(l_1) }}. \end{equation}

Apart from the neglect of terms of order $m_e/m_i$ and fluctuations in the magnetic-field strength, $\delta B_\|$, this upper bound on the growth rate is completely general and applies to any local gyrokinetic instability. It applies to ITG and electron-temperature-gradient modes, kinetic and resistive ballooning modes, trapped-ion and trapped-electron modes and microtearing modes, as well as to the so-called universal and ubiquitous instabilities.

A particularly simple and important case is that of a hydrogen plasma without other ions and $k_\perp \rho _e \ll 1$. Noting that $\omega _{\ast i} = - \tau \omega _{\ast e}$ and using the asymptotic forms (4.11) and (4.12), we find

(5.28)\begin{equation} \frac{\gamma}{|\omega_{{\ast} e}|} \le \sqrt{\frac{\tau(\varGamma_{0i} + \tau)}{(1 + \tau) (1 - \varGamma_{0i})} } \left( \sqrt{\tau M(\eta_i, b_i)} + \sqrt{1 + \frac{3 \eta_e^2}{2}} \right) + \beta_e \sqrt{\frac{ 1 + 2 \eta_a + 7 \eta_e^2/2}{2 b_e \left(\beta_e + 2 b_e \right)}}, \end{equation}

where the first term on the right is evaluated at $l=l_0$ and the second term (which is proportional to $\beta _e$) at $l=l_1$. Both terms give an upper bound on $\gamma$ that remains finite in the long-wavelength limit because $\omega _{\ast e}$ is proportional to $k_\alpha$ and

(5.29)\begin{equation} 1 - \varGamma_{0i} \simeq b_i = (k_\perp \rho_i)^2, \end{equation}

in the limit $b_i \ll 1$. Furthermore, as long as $k_\perp \rho _e \ll 1$, the growth rate is subject to a bound equal to

(5.30)\begin{equation} \gamma < C_0 \left( 1 + \tau^{{-}1/2} \right) \frac{v_{Ti}}{L_\perp} + \frac{C_1 \beta_e}{\sqrt{\beta_e + 2 b_e}} \frac{v_{Te}}{L_\perp}, \end{equation}

where $C_0$ and $C_1$ are numbers of order unity, $v_{Ti}$ denotes the ion thermal speed and $L_\perp$ the length scale of the equilibrium density and temperature gradients. In the opposite limit, $k_\perp \rho _e \gg 1$, the term proportional to $\beta _e$ can be neglected and we instead obtain

(5.31)\begin{equation} {\gamma} \le \frac{\tau {|\omega_{{\ast} e}|}}{1 + \tau} \sqrt{\frac{1-\eta_e + 5 \eta_e^2/4}{2 {\rm \pi}b_e(l_0)}} = \frac{C_2 v_{Te}}{(1 + \tau^{{-}1}) L_\perp}, \end{equation}

where $v_{Te}$ denotes the electron thermal speed and $C_2$ is a number of order unity.

6. Optimal bounds

The bounds (5.11) and (5.27) are not optimal and can be improved. In this section, we derive the best possible bound, in a sense that will be made precise, for the simplest case of a hydrogen plasma with adiabatic electrons. If $\varphi = e \delta \phi _{\boldsymbol k} / T_i$ and $g = g_{i\boldsymbol k}$, we have

(6.1)\begin{gather} \varphi = \frac{1}{n (1+\tau)} \int g {\rm J}_{0} \,{\rm d}^3v, \end{gather}

(6.2)\begin{gather} H = nT_i \left\langle \frac{1}{n} \int \frac{|g|^2}{F_{i0}} \,{\rm d}^3v - (1+\tau) |\varphi|^2 \right\rangle, \end{gather}

(6.3)\begin{gather} D = \frac{\eta_i \omega_{{\ast} i} T_i}{2i} \left\langle \int \left(\varphi^\ast g - \varphi g^\ast\right) x^2 {\rm J}_{0i} \,{\rm d}^3v \right\rangle, \end{gather}

where $x^2 = m_i v^2 / 2 T_i$. Here $D$ and $H$ are thus quadratic functionals of $g$, and the challenge is to maximise the ratio $D[g]/H[g]$ over all such functions.

In order to do so, we first note that $D$ and $\varphi$ only depend on two moments of $g$, namely,

(6.4)\begin{equation} K_j[g] = \frac{1}{n} \int g x^{2j} {\rm J}_{0i} \,{\rm d}^3v, \end{equation}

where $j = 0$ or 1. We can therefore begin by minimising $H[g]$ over all functions with given values of these two moments. Using Lagrange multipliers, $c_0$ and $c_1$, we are thus led to minimise the functional

(6.5)\begin{equation} H[g] - 2c_0 K_0[g] - 2c_1 K_1[g], \end{equation}

which gives

(6.6)\begin{equation} g = \left( c_0 + c_1 x^2 \right) {\rm J}_{0i} F_{i0}. \end{equation}

We have thus reduced our problem to that of finding the maximum value of $D/H$ expressed as a ratio of two quadratic forms in the coefficients $c_j$. Note that conventional eigenmodes, i.e. functions satisfying the linearised version of (2.3) are, in general, not of the form (6.6). This equation describes modes of optimal free-energy growth, which are distinct from eigenmodes and will be studied in greater detail in Part 2 of this series of papers.

If we write

(6.7)\begin{equation} G_j(b_i) = \frac{1}{n} \int F_{i0} x^{2j} {\rm J}_{0i}^2 \,{\rm d}^3v, \end{equation}

so that

(6.8)\begin{gather} G_0(b_i) = \varGamma_0(b_i), \end{gather}

(6.9)\begin{gather} G_1(b_i) = \left( \frac{3}{2} - b_i \right) \varGamma_0(b_i) + b_i \varGamma_1(b_i), \end{gather}

(6.10)\begin{gather} G_2(b_i) = \left( \frac{15}{4} - 5b_i + 2 b_i^2 \right) \varGamma_0(b_i) + \left( 4 - 2 b_i \right) b_i \varGamma_1(b_i), \end{gather}

then

(6.11)\begin{equation} D = \frac{n T_i G(b_i) }{2i(1+\tau)} \left( c_0^\ast c_1 - c_0 c_1^\ast \right), \end{equation}

where

(6.12)\begin{equation} G(b) = G_0(b_i) G_2(b_i) - G_1^2(b_i) = \left( \frac{3}{2} - 2 b_i + b_i^2 \right) \varGamma_0^2(b_i) + b_i \varGamma_0 (b_i) \varGamma_1(b_i) - b_i^2 \varGamma_1^2 (b_i^2), \end{equation}

and

(6.13)\begin{align} H = n T_i \left[ G_0 \left( 1 - \frac{G_0}{1\!+\!\tau} \right) c_0 c_0^\ast{+} G_1 \left( 1 - \frac{G_0}{1+\tau} \right) \left( c_0^\ast c_1 + c_0 c_1^\ast \right) + \left( G_2 - \frac{G_1^2}{1+\tau} \right) c_1 c_1^\ast \right].\end{align}

In order to maximise the ratio and calculate

(6.14)\begin{equation} \hat \gamma = \max_{c_0, c_1} \left( \frac{D}{H} \right) , \end{equation}

we consider the variations

(6.15)\begin{gather} \delta D = \frac{n T_i G }{2i(1+\tau)} \left( c_1 \delta c_0^\ast{-} c_0 \delta c_1^\ast \right) + {\rm c.c.}, \end{gather}

(6.16)\begin{gather} \delta H = n T_i \left[ G_0 \left( 1 - \frac{G_0}{1+\tau} \right) c_0 \delta c_0^\ast{+} G_1 \left( 1 - \frac{G_0}{1+\tau} \right) \left( c_1 \delta c_0^\ast{+} c_0 \delta c_1^\ast \right) \right. \end{gather}

(6.17)\begin{gather}\left. + \left( G_2 - \frac{G_1^2}{1+\tau} \right) c_1 \delta c_1^\ast \right]+ {\rm c.c.}, \end{gather}

where $\textrm {c.c.}$ stands for the complex conjugate, and we note that the maximum is reached when

(6.18)\begin{equation} \delta D = \hat \gamma \delta H, \end{equation}

which gives a system of equations

(6.19)\begin{equation} \frac{2{\rm i}\hat \gamma}{\eta_i \omega_{{\ast} i}} \left[ \begin{array}{cc} G_0 \left( 1 + \tau - G_0 \right) & G_1 \left( 1 + \tau - G_0 \right)\\ G_1 \left( 1 + \tau - G_0 \right) & G_2(1 + \tau) - G_1^2 \end{array} \right] \left[ \begin{array}{c} c_0 \\ c_1 \end{array} \right] = G \left[ \begin{array}{c} -c_1 \\ c_0 \end{array} \right], \end{equation}

which has non-zero solutions if

(6.20)\begin{equation} \hat \gamma = \frac{|\eta_i \omega_{{\ast} i}|}{2} \sqrt{\frac{G(b_i)}{(1+\tau) [1 + \tau - G_0(b_i)]}}. \end{equation}

This is the ‘optimal’ bound on the growth rate that can be obtained within our formalism in the sense that no lower bound is possible. Indeed, growth of the free energy at this rate is realised if no collisions are present and the distribution function is chosen as dictated by (6.6) with $c_0$ and $c_1$ satisfying the eigenvalue problem (6.18). The bound (6.20) is shown in figure 1 and is lower than our previous result (5.11) by a factor of 2 and $\sqrt {5}$ in the limits of long and short wavelengths, respectively,

(6.21)\begin{equation} \hat \gamma \rightarrow \left\{ \begin{array}{@{}c l} \dfrac{|\eta_i \omega_{{\ast} i}|}{2} \sqrt{\dfrac{3}{2 \tau (1+\tau)}} , & b_i \ll 1\\ \dfrac{| \eta_i \omega_{{\ast} i}|}{(1+\tau) \sqrt{8 {\rm \pi}b_i}}, & b_i \gg 1. \end{array} \right. \end{equation}

7. Bounds on nonlinear growth

Our most general bound (5.27) is not optimal and will be improved substantially in our next publication, but its most important implication follows already from this crude form. The right-hand side is a bounded function of the mode numbers $(k_\psi, k_\alpha )$, and the linear growth rate can therefore never exceed the maximum

(7.1)\begin{equation} \gamma_{\rm max} = \sup_{\boldsymbol k} \; \gamma_{\rm bound} ({\boldsymbol k}). \end{equation}

As we now show, this conclusion also holds for nonlinear growth.

Consider the evolution of a set of fluctuations governed by the gyrokinetic system of equations starting from some arbitrary initial condition, specified by the distribution functions $\delta F_a$ of all species at $t=0$. According to (3.8) the instantaneous growth of the total free energy,

(7.2)\begin{equation} H_{\rm tot}(t) = \sum_{\boldsymbol k} H({\boldsymbol k},t) \end{equation}

is bounded by

(7.3)\begin{equation} \frac{{\rm d} H_{\rm tot}}{{\rm d}t} \le 2 \sum_{\boldsymbol k} D({\boldsymbol k},t), \end{equation}

where each term is subject to the bound

(7.4)\begin{equation} D({\boldsymbol k},t) \le \gamma_{\rm bound}({\boldsymbol k}) H({\boldsymbol k},t). \end{equation}

The growth rate of the total free energy is therefore limited by twice the maximum linear growth

(7.5)\begin{equation} \frac{{\rm d} \ln H_{\rm tot}}{{\rm d}t} \le 2 \gamma_{\rm max}. \end{equation}

This bound holds for fluctuations of arbitrary amplitude within the gyrokinetic formalism. In particular, it must hold in any gyrokinetic simulation of turbulence.

Moreover, if collisions are absent, then instantaneous growth of the free energy is possible at any positive rate up to the ‘optimal’ one, which for the particularly simple case of adiabatic electrons was derived in the previous subsection. To see this, suppose the bounds on the right-hand side of (7.1) are chosen optimally in the sense that

(7.6)\begin{equation} \gamma_{\rm bound} ({\boldsymbol k}) = \sup_{g} \, \frac{D[g,{\boldsymbol k}]}{H[g,{\boldsymbol k}]}, \end{equation}

where $D$ and $H$ are now considered to be quadratic functionals of the distribution functions $g = \{ g_a \}$ of all species. This means, then, that there is a choice of wavenumber and initial data such that the free energy grows at a rate arbitrarily close to $2 \gamma _\textrm {max}$. Conversely, there is a similar limit on the rate at which the free energy can decay in the absence of collisions,

(7.7)\begin{equation} \frac{{\rm d} \ln H_{\rm tot}}{{\rm d}t} \ge - 2 \gamma_{\rm max}, \end{equation}

as follows from the observation that $D[g,{\boldsymbol k}]$ is odd in the wavenumber $\boldsymbol k$ at fixed $g$ whereas $H$ is even. The transformation ${\boldsymbol k} \rightarrow - {\boldsymbol k}$ thus changes the sign of the ratio ${D[g,{\boldsymbol k}]}/{H[g,{\boldsymbol k}]}$ if $g$ is held constant.² Any upper bound on this ratio therefore automatically implies a similar lower bound when collisions are absent.

8. Conclusions

As we have shown, it is possible to derive rigorous upper bounds on the growth rate of linear instabilities and on the nonlinear growth of free energy in gyrokinetics. Unlike most other results in the field, these bounds are universal and hold in plasmas with any number of particle species regardless of collisionality and magnetic-field geometry. For simplicity, we have taken the plasma pressure (beta) to be sufficiently small that fluctuations in the magnetic-field strength can be neglected, $\delta B_\| = 0$, but this restriction will be removed in Part 2 in the present series of papers.

In the case of a plasma with a single kinetic ion species and ‘adiabatic’ electrons, the bound is given by (6.20) and is of order

(8.1)\begin{equation} \gamma_{\rm bound} \sim \frac{k_\perp \rho_i}{\sqrt{\tau (1+\tau)}} \cdot \frac{ v_{Ti}}{L_\perp } \end{equation}

for $k_\perp \rho _i < 1$ and

(8.2)\begin{equation} \gamma_{\rm bound} \sim \frac{v_{Ti}}{(1 + \tau) L_\perp} \end{equation}

for shorter wavelengths. The dependence on the parameter $\tau = T_i/T_e$ reflects a well-known unfavourable dependence of the ITG growth rate on electron temperature.

The bound (5.27) we found on instabilities with kinetic electrons is less restrictive and remains finite in the limit $k_\perp \rho _i \rightarrow 0$. It is a sum of two distinct contributions: an electrostatic term and an electromagnetic term that vanishes if $\beta _e \rightarrow 0$. As we shall show in the next publication of this series, this result is not qualitatively affected by the inclusion of parallel magnetic fluctuations.

Actual microinstability growth rates must lie below these bounds. For instance, toroidal ITG modes with adiabatic electrons and $k_\perp \rho _i \ll 1$ have growth rates

(8.3)\begin{equation} \gamma \sim \sqrt{\frac{\eta_i \omega_{{\ast} i} \omega_{di}}{\tau} } \sim \frac{k_\perp \rho_i}{\sqrt{\tau}} \cdot \frac{ v_{Ti}}{\sqrt{R L_\perp} } \end{equation}

in the strongly driven limit (Biglari et al. 1989; Romanelli 1989; Plunk et al. 2014; Zocco et al. 2018), and trapped-ion modes have a similar growth rate (Biglari et al. 1989). Here $R$ denotes the radius of curvature of the magnetic field, so that $\omega _{di} \sim (k_\perp \rho _i) v_{Ti} / R$. Due to the assumption $| \omega _{di} /\omega _{\ast i} |\sim L_\perp / R \ll 1$ (corresponding to strong instability drive) made in the derivation of this estimate, the growth rate is smaller than our upper bound. Similarly, in the theory of kinetic ballooning modes, the assumption $L_\perp / R \ll 1$ leads to growth rates of order (Tang, Connor & Hastie 1980; Aleynikova et al. 2018)

(8.4)\begin{equation} \gamma \sim \frac{\sqrt{\omega_{di} \left[ (1 + \eta_i) \omega_{{\ast} i} - (1 + \eta_e) \omega_{{\ast} e} \right]}}{k_\perp \rho_i}. \end{equation}

This growth rate never exceeds our bound (5.28) and scales as our estimate (5.30). In less strongly driven cases, the growth rate is lower.

Although all our results are quite general, they do not encompass all instabilities of interest. Kink modes and tearing modes sometimes need a gyrokinetic treatment in a thin layer around a resonant magnetic surface, where magnetic reconnection may occur, but take their energy from the exterior region and depend on the overall plasma current profile (Hazeltine, Dobrott & Wang 1975; Drake & Lee 1977). Such instabilities cannot adequately be described in the geometry of a magnetic flux tube (Connor et al. 2014,  2019) and are not subject to the bounds derived in the present paper. Mathematically, they are not covered by our treatment because the solution of the gyrokinetic equation involves matching to the exterior region, whose destabilising influence is usually described by a parameter $\varDelta '$, making these modes non-local in nature. However, microtearing modes which are driven by local gradients are subject to our bound (5.27) on electromagnetic instabilities.

As already remarked, instabilities driven by equilibrium flow shear would need an additional term in the gyrokinetic equation. The parallel-velocity-gradient instability, which could then be treated, is known to be capable of causing subcritical turbulence below the linear stability threshold due to transiently growing modes (Barnes et al. 2011). The latter would, however, be subject to bounds similar to those we have derived, because these also apply to nonlinear growth and thus limit the possibility of subcritical turbulence excitation.